/[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 4102 by gross, Wed Dec 12 06:17:03 2012 UTC revision 4122 by gross, Thu Dec 20 05:42:35 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, grad, inner, integrate, interpolate, kronecker, boundingBoxEdgeLengths, vol, sqrt, length  from esys.escript import 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 94  class Regularization(CostFunction): Line 94  class Regularization(CostFunction):
94                    
95    
96          """          """
         if numLevelSets>1:  
               raise ValueError("Currently only numLevelSets<=1 is supported.")  
           
97          if w0 == None and w1==None:          if w0 == None and w1==None:
98                raise ValueError("Values for w0 or for w1 must be given.")                raise ValueError("Values for w0 or for w1 must be given.")
99      if wc == None and  numLevelSets>1:      if wc == None and  numLevelSets>1:
# Line 251  class Regularization(CostFunction): Line 248  class Regularization(CostFunction):
248          :rtype: `LinearPDE`          :rtype: `LinearPDE`
249          """          """
250          return self.__pde          return self.__pde
251        
252      def getDualProduct(self, m, r):      def getDualProduct(self, m, r):
253          """          """
254          returns the dual product of a gradient represented by X=r[1] and Y=r[0]          returns the dual product of a gradient represented by X=r[1] and Y=r[0]
# Line 410  class Regularization(CostFunction): Line 407  class Regularization(CostFunction):
407          DIM=self.getDomain().getDim()          DIM=self.getDomain().getDim()
408          numLS=self.getNumLevelSets()          numLS=self.getNumLevelSets()
409    
         n=0  
   
410          if self.__w0 is not None:          if self.__w0 is not None:
411              Y = m * self.__w0 * mu              Y = m * self.__w0 * mu
412          else:          else:
413              Y = Data()          if numLS == 1:
414          n+=numLS            Y = Scalar(0,  grad_m.getFunctionSpace())
415            else:
416                  Y = Data(0, (numLS,) , grad_m.getFunctionSpace())
417    
418          if self.__w1 is not None:          if self.__w1 is not None:
419            X=grad_m*self.__w1
420              if numLS == 1:              if numLS == 1:
421                  X=grad_m* self.__w1*mu                  X=grad_m* self.__w1*mu
422              else:              else:
                 X=grad_m[k,:]*self.__w1[k,:]  
423                  for k in xrange(numLS):                  for k in xrange(numLS):
424                      X[k,:]*=mu[k]                      X[k,:]*=mu[k]
425          else:          else:
426              X = Data()              X = Data(0, grad_m.getShape(), grad_m.getFunctionSpace())
427          if numLS > 1:          
428              raise NotImplementedError          # cross gradient terms:
429            if numLS > 1:    
430          for  k in xrange(numLS):
431                 grad_m_k=grad_m[k,:]
432             l2_grad_m_k = length(grad_m_k)**2
433             for  l in xrange(k):
434               grad_m_l=grad_m[l,:]
435               l2_grad_m_l = length(grad_m_l)**2
436               grad_m_lk = inner(grad_m_l, grad_m_k)
437               f=  mu_c[l,k]* self.__wc[l,k]
438               X[l,:] += f * ( l2_grad_m_l *  grad_m_l - grad_m_lk * grad_m_k )
439               X[k,:] += f * ( l2_grad_m_k *  grad_m_k - grad_m_lk * grad_m_l )
440            
441          return ArithmeticTuple(Y, X)          return ArithmeticTuple(Y, X)
442    
443    
# Line 462  class Regularization(CostFunction): Line 471  class Regularization(CostFunction):
471                          for i in xrange(DIM): A[k,i,k,i]=self.__w1[k,i] * mu[k]                          for i in xrange(DIM): A[k,i,k,i]=self.__w1[k,i] * mu[k]
472                                                    
473              if numLS > 1:              if numLS > 1:
474                  raise NotImplementedError             for  k in xrange(numLS):
475                     grad_m_k=grad_m[k,:]
476                 l2_grad_m_k = length(grad_m_k)**2
477                 o_kk=outer(grad_m_k, grad_m_k)
478                 for  l in xrange(k):
479                    grad_m_l=grad_m[l,:]
480                    l2_grad_m_l = length(grad_m_l)**2
481                    i_lk = inner(grad_m_l, grad_m_k)
482                    o_lk = outer(grad_m_l, grad_m_k)
483                    o_kl = outer(grad_m_k, grad_m_l)
484                    o_ll=outer(grad_m_l, grad_m_l)
485                    f=  mu_c[l,k]* self.__wc[l,k]
486                    
487                    A[l,:,l:] += f * ( l2_grad_m_k * kronecker(DIM) - o_kk )
488                    A[l,:,k:] += f * ( 2 * o_lk -   o_kl - i_lk * kronecker(DIM) )
489                    A[k,:,l:] += f * ( 2 * o_kl -   o_lk - i_lk * kronecker(DIM) )
490                    A[k,:,k:] += f * ( l2_grad_m_l * kronecker(DIM) - o_ll )
491              self.getPDE().setValue(A=A)              self.getPDE().setValue(A=A)
492          #self.getPDE().resetRightHandSideCoefficients()          #self.getPDE().resetRightHandSideCoefficients()
493          #self.getPDE().setValue(X=r[1])          #self.getPDE().setValue(X=r[1])

Legend:
Removed from v.4102  
changed lines
  Added in v.4122

  ViewVC Help
Powered by ViewVC 1.1.26