/[escript]/trunk/downunder/py_src/regularizations.py
ViewVC logotype

Diff of /trunk/downunder/py_src/regularizations.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 4078 by gross, Thu Nov 15 03:45:24 2012 UTC revision 4079 by gross, Fri Nov 16 07:59:01 2012 UTC
# Line 27  from costfunctions import CostFunction Line 27  from costfunctions import CostFunction
27    
28  import numpy as np  import numpy as np
29  from esys.escript.linearPDEs import LinearPDE, IllegalCoefficientValue  from esys.escript.linearPDEs import LinearPDE, IllegalCoefficientValue
30  from esys.escript import Data, grad, inner, integrate, kronecker, boundingBoxEdgeLengths, interpolate  from esys.escript import Data, grad, inner, integrate, kronecker, boundingBoxEdgeLengths, interpolate, vol
31  from esys.escript.pdetools import ArithmeticTuple  from esys.escript.pdetools import ArithmeticTuple
32    
33  class Regularization(CostFunction):  class Regularization(CostFunction):
# Line 105  class Regularization(CostFunction): Line 105  class Regularization(CostFunction):
105                         # we count s0, then s1, then sc (k<l).                         # we count s0, then s1, then sc (k<l).
106      # THE S0 weighting factor      # THE S0 weighting factor
107      n=0      n=0
108        VV=vol(domain)
109          if not s0 is None:          if not s0 is None:
110          s0 = interpolate(s0,self.__pde.getFunctionSpaceForCoefficient('D'))              s0 = interpolate(s0,self.__pde.getFunctionSpaceForCoefficient('D'))    
111          s=s0.getShape()          s=s0.getShape()
# Line 113  class Regularization(CostFunction): Line 114  class Regularization(CostFunction):
114               V=integrate(s0)               V=integrate(s0)
115               if V > 0:               if V > 0:
116                 self.__weight_index.append(n)                 self.__weight_index.append(n)
117                 s0*=1./V                 s0*=VV/V
118                   else:                   else:
119                 s0=None                 s0=None
120           else:           else:
# Line 122  class Regularization(CostFunction): Line 123  class Regularization(CostFunction):
123               if s == (numLevelSets,):               if s == (numLevelSets,):
124               for k in xrange(numLevelSets):               for k in xrange(numLevelSets):
125                  V=integrate(s0[k])                  V=integrate(s0[k])
126                  if Lsup(V) > 0:                  if V > 0:
127                      self.__weight_index.append(n+k)                      self.__weight_index.append(n+k)
128                      s0[k]*=1/V                      s0[k]*=VV/V
129           else:           else:
130               raise ValueError("Unexpected shape %s for weight s0."%s)                 raise ValueError("Unexpected shape %s for weight s0."%s)  
131          self.__s0=s0          self.__s0=s0
# Line 137  class Regularization(CostFunction): Line 138  class Regularization(CostFunction):
138                    
139          if numLevelSets == 1 :          if numLevelSets == 1 :
140               if s == (DIM,) :               if s == (DIM,) :
141                 V=integrate(inner(s1, self.__L2))                 V=integrate(inner(s1, 1/self.__L2))
142                 if V > 0:                 if V > 0:
143                self.__weight_index.append(n)                self.__weight_index.append(n)
144                    s1*=1/V                    s1*=VV/V
145                   print "REG SCALE = ",s1
146           else:           else:
147               raise ValueError("Unexpected shape %s for weight s1."%s)               raise ValueError("Unexpected shape %s for weight s1."%s)
148              else:              else:
# Line 148  class Regularization(CostFunction): Line 150  class Regularization(CostFunction):
150               for k in xrange(numLevelSets):               for k in xrange(numLevelSets):
151                  for i in xrange(DIM):                  for i in xrange(DIM):
152                  ww=s1[k,:]                  ww=s1[k,:]
153                  V=integrate(inner(ww,self.__L2))                  V=integrate(inner(ww,1/self.__L2))
154                      if V > 0:                      if V > 0:
155                         self.__weight_index.append(n+i)                         self.__weight_index.append(n+i)
156                         s1[k,:]=ww*(1./V)                         s1[k,:]=ww*(VV/V)
157           else:           else:
158               raise ValueError("Unexpected shape %s for weight s1."%s)                     raise ValueError("Unexpected shape %s for weight s1."%s)      
159                        
# Line 173  class Regularization(CostFunction): Line 175  class Regularization(CostFunction):
175                  V=integrate(ww)                  V=integrate(ww)
176                      if V > 0:                      if V > 0:
177                         self.__weight_index.append(n+k*numLevelSets+l)                         self.__weight_index.append(n+k*numLevelSets+l)
178                         sc[l,k]=ww/V*self.__L4                         sc[l,k]=ww*VV/V*self.__L4
179                         sc[k,l]=0                         sc[k,l]=0
180                else:                else:
181               raise ValueError("Unexpected shape %s for weight s0."%s)                     raise ValueError("Unexpected shape %s for weight s0."%s)      
# Line 305  class Regularization(CostFunction): Line 307  class Regularization(CostFunction):
307              if self.__s0 is not None:              if self.__s0 is not None:
308              if numLS == 1:              if numLS == 1:
309                   D=self.__s0 * mu[n]                   D=self.__s0 * mu[n]
                  print "D =",D  
310                  else:                    else:  
311                       D=self.getPDE().createCoefficient("D")                       D=self.getPDE().createCoefficient("D")
312                   for k in xrange(numLS): D[k,k]=self.__s0[k] * mu[n+k]                   for k in xrange(numLS): D[k,k]=self.__s0[k] * mu[n+k]
# Line 316  class Regularization(CostFunction): Line 317  class Regularization(CostFunction):
317              if self.__s1 is not None:              if self.__s1 is not None:
318             if numLS == 1:             if numLS == 1:
319             for i in xrange(DIM): A[i,i]=self.__s1[i] * mu[n]             for i in xrange(DIM): A[i,i]=self.__s1[i] * mu[n]
            print "A =",A  
320             else:             else:
321                 for k in xrange(numLS):                 for k in xrange(numLS):
322                      for i in xrange(DIM): A[k,i,k,i]=self.__s1[k,i] * mu[n+k]                      for i in xrange(DIM): A[k,i,k,i]=self.__s1[k,i] * mu[n+k]
323              n+=numLS              n+=numLS
324              if self.__sc is not None:              if self.__sc is not None:
325                  raise NotImplementedError                  raise NotImplementedError
326                print A
327              self.getPDE().setValue(A=A)              self.getPDE().setValue(A=A)
328          self.getPDE().resetRightHandSideCoefficients()          self.getPDE().resetRightHandSideCoefficients()
329          print A          self.getPDE().setValue(X=r[1])
330          print r[1]          print "X only: ",self.getPDE().getSolution()
331          print r[0]          self.getPDE().resetRightHandSideCoefficients()
332            self.getPDE().setValue(Y=r[0])
333            print "Y only: ",self.getPDE().getSolution()
334            
335            self.getPDE().resetRightHandSideCoefficients()
336          self.getPDE().setValue(X=r[1], Y=r[0])          self.getPDE().setValue(X=r[1], Y=r[0])
337          return self.getPDE().getSolution()          return self.getPDE().getSolution()
338                    
Legend:
Removed from v.4078  
changed lines
  Added in v.4079
  ViewVC Help
Powered by ViewVC 1.1.26