/[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 4123 by gross, Thu Dec 20 05:42:35 2012 UTC revision 4124 by gross, Fri Dec 21 01:35:43 2012 UTC
# Line 25  __all__ = ['Regularization'] Line 25  __all__ = ['Regularization']
25  from costfunctions import CostFunction  from costfunctions import CostFunction
26    
27  import numpy as np  import numpy as np
28  from esys.escript import Data, Scalar, grad, inner, integrate, interpolate, kronecker, boundingBoxEdgeLengths, vol, sqrt, length  from esys.escript import ReducedFunction, outer, Data, Scalar, grad, inner, integrate, interpolate, kronecker, boundingBoxEdgeLengths, vol, sqrt, length
29  from esys.escript.linearPDEs import LinearPDE, IllegalCoefficientValue  from esys.escript.linearPDEs import LinearPDE, IllegalCoefficientValue
30  from esys.escript.pdetools import ArithmeticTuple  from esys.escript.pdetools import ArithmeticTuple
31    
# Line 165  class Regularization(CostFunction): Line 165  class Regularization(CostFunction):
165          else:          else:
166               wc = interpolate(wc,self.__pde.getFunctionSpaceForCoefficient('A'))               wc = interpolate(wc,self.__pde.getFunctionSpaceForCoefficient('A'))
167               sc=wc.getShape()               sc=wc.getShape()
168               raise ValueError("Unexpected shape %s for weight wc."%sc)               if not sc == (numLevelSets, numLevelSets):
169                    raise ValueError("Unexpected shape %s for weight wc."%(sc,))
170          # ============= now we rescale weights: ======================================          # ============= now we rescale weights: ======================================
171          L2s=np.asarray(boundingBoxEdgeLengths(domain))**2          L2s=np.asarray(boundingBoxEdgeLengths(domain))**2
172          L4=1/np.sum(1/L2s)**2          L4=1/np.sum(1/L2s)**2
173          if numLevelSets is 1 :          if numLevelSets == 1 :
174              A=0              A=0
175              if w0 is not None:              if w0 is not None:
176              A = integrate(w0)              A = integrate(w0)
# Line 184  class Regularization(CostFunction): Line 185  class Regularization(CostFunction):
185              else:              else:
186             raise ValueError("Non-positive weighting factor detected.")             raise ValueError("Non-positive weighting factor detected.")
187          else:          else:
188    
189           for k in xrange(numLevelSets):           for k in xrange(numLevelSets):
190               A=0               A=0
191                   if w0 is not None:                   if w0 is not None:
# Line 200  class Regularization(CostFunction): Line 202  class Regularization(CostFunction):
202                 raise ValueError("Non-positive weighting factor for level set component %d detected."%k)                 raise ValueError("Non-positive weighting factor for level set component %d detected."%k)
203                    
204               # and now the cross-gradient:               # and now the cross-gradient:
205               for l in xrange(k):                 if wc is not None:
206                   A = integrate(wc[l,k])/L4                   for l in xrange(k):  
207                   if A > 0:                      A = integrate(wc[l,k])/L4
208                      f = scale_c[l,k]/A                      if A > 0:
209                      wc[l,k]*=f                         f = scale_c[l,k]/A
210                       else:                         wc[l,k]*=f
211                      raise ValueError("Non-positive weighting factor for cross-gradient level set components %d and %d detected."%(l,k))  #                        else:
212    #                      raise ValueError("Non-positive weighting factor for cross-gradient level set components %d and %d detected."%(l,k))
213                                    
214          self.__w0=w0          self.__w0=w0
215          self.__w1=w1          self.__w1=w1
216          self.__wc=wc          self.__wc=wc
           
   
217    
218          self.__pde_is_set=False                  self.__pde_is_set=False        
219          if self.__numLevelSets > 1:          if self.__numLevelSets > 1:
# Line 348  class Regularization(CostFunction): Line 349  class Regularization(CostFunction):
349          numLS=self.getNumLevelSets()          numLS=self.getNumLevelSets()
350          if mu_c is None or numLS < 2:          if mu_c is None or numLS < 2:
351          self.__mu_c = np.ones((numLS,numLS))          self.__mu_c = np.ones((numLS,numLS))
352        if isinstance(mu_c, float) or isinstance(mu_c, int):
353            self.__mu_c = np.zeros((numLS,numLS))
354            self.__mu_c[:,:]=mu_c
355      else:      else:
356          mu_c=np.asarray(mu_c)          mu_c=np.asarray(mu_c)
357          if mu_c.shape == (numLS,numLS):          if mu_c.shape == (numLS,numLS):
# Line 407  class Regularization(CostFunction): Line 411  class Regularization(CostFunction):
411          DIM=self.getDomain().getDim()          DIM=self.getDomain().getDim()
412          numLS=self.getNumLevelSets()          numLS=self.getNumLevelSets()
413    
414            print "WARNING: WRONG FUNCTION SPACE"
415            
416            grad_m=grad(m, ReducedFunction(m.getDomain()))
417          if self.__w0 is not None:          if self.__w0 is not None:
418              Y = m * self.__w0 * mu              Y = m * self.__w0 * mu
419          else:          else:
# Line 484  class Regularization(CostFunction): Line 491  class Regularization(CostFunction):
491                  o_ll=outer(grad_m_l, grad_m_l)                  o_ll=outer(grad_m_l, grad_m_l)
492                  f=  mu_c[l,k]* self.__wc[l,k]                  f=  mu_c[l,k]* self.__wc[l,k]
493                                    
494                  A[l,:,l:] += f * ( l2_grad_m_k * kronecker(DIM) - o_kk )                  A[l,:,l,:] += f * ( l2_grad_m_k * kronecker(DIM) - o_kk )
495                  A[l,:,k:] += f * ( 2 * o_lk -   o_kl - i_lk * kronecker(DIM) )                  A[l,:,k,:] += f * ( 2 * o_lk -   o_kl - i_lk * kronecker(DIM) )
496                  A[k,:,l:] += f * ( 2 * o_kl -   o_lk - i_lk * kronecker(DIM) )                  A[k,:,l,:] += f * ( 2 * o_kl -   o_lk - i_lk * kronecker(DIM) )
497                  A[k,:,k:] += f * ( l2_grad_m_l * kronecker(DIM) - o_ll )                  A[k,:,k,:] += f * ( l2_grad_m_l * kronecker(DIM) - o_ll )
498              self.getPDE().setValue(A=A)              self.getPDE().setValue(A=A)
499          #self.getPDE().resetRightHandSideCoefficients()          #self.getPDE().resetRightHandSideCoefficients()
500          #self.getPDE().setValue(X=r[1])          #self.getPDE().setValue(X=r[1])

Legend:
Removed from v.4123  
changed lines
  Added in v.4124

  ViewVC Help
Powered by ViewVC 1.1.26