/[escript]/trunk/escript/py_src/linearPDEs.py
ViewVC logotype

Diff of /trunk/escript/py_src/linearPDEs.py

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

revision 104 by jgs, Fri Dec 17 07:43:12 2004 UTC revision 114 by jgs, Fri Mar 4 07:12:37 2005 UTC
# Line 10  import escript Line 10  import escript
10  import util  import util
11  import numarray  import numarray
12    
 def identifyDomain(domain=None,data={}):  
      """  
      @brief Return the Domain which is equal to the input domain (if not None)  
      and is the domain of all Data objects in the dictionary data.  
      An exception is raised if this is not possible  
   
      @param domain  
      @param data  
      """  
      # get the domain used by any Data object in the list data:  
      data_domain=None  
      for d in data.itervalues():  
           if isinstance(d,escript.Data):  
              if not d.isEmpty(): data_domain=d.getDomain()  
      # check if domain and data_domain are identical?  
      if domain == None:  
          if data_domain == None:  
               raise ValueError,"Undefined PDE domain. Specify a domain or use a Data class object as coefficient"  
      else:  
          if data_domain == None:  
               data_domain=domain  
          else:  
            if not data_domain == domain:  
                  raise ValueError,"Domain of coefficients doesnot match specified domain"  
      # now we check if all Data class object coefficients are defined on data_domain:  
      for i,d in data.iteritems():  
          if isinstance(d,escript.Data):  
             if not d.isEmpty():  
                if not data_domain==d.getDomain():  
                  raise ValueError,"Illegal domain for coefficient %s."%i  
      # done:  
      return data_domain  
   
 def identifyNumEquationsAndSolutions(dim,coef={}):  
      # get number of equations and number of unknowns:  
      numEquations=0  
      numSolutions=0  
      for i in coef.iterkeys():  
         if not coef[i].isEmpty():  
            res=_PDECoefficientTypes[i].estimateNumEquationsAndNumSolutions(coef[i].getShape(),dim)  
            if res==None:  
                raise ValueError,"Illegal shape %s of coefficient %s"%(coef[i].getShape().__str__(),i)  
            else:  
                numEquations=max(numEquations,res[0])  
                numSolutions=max(numSolutions,res[1])  
      return numEquations,numSolutions  
   
13    
14  def _CompTuple2(t1,t2):  def _CompTuple2(t1,t2):
15     """     """
# Line 70  def _CompTuple2(t1,t2): Line 23  def _CompTuple2(t1,t2):
23     elif dif>0: return -1     elif dif>0: return -1
24     else: return 0     else: return 0
25    
26  class PDECoefficientType:  class PDECoefficient:
27      """      """
28      @brief      @brief
29      """      """
30      # identifier for location of Data objects defining coefficients      # identifier for location of Data objects defining COEFFICIENTS
31      INTERIOR=0      INTERIOR=0
32      BOUNDARY=1      BOUNDARY=1
33      CONTACT=2      CONTACT=2
34      CONTINUOUS=3      CONTINUOUS=3
35      # identifier in the pattern of coefficients:      # identifier in the pattern of COEFFICIENTS:
36      # the pattern is a tuple of EQUATION,SOLUTION,DIM where DIM represents the spatial dimension, EQUATION the number of equations and SOLUTION the      # the pattern is a tuple of EQUATION,SOLUTION,DIM where DIM represents the spatial dimension, EQUATION the number of equations and SOLUTION the
37      # number of unknowns.      # number of unknowns.
38      EQUATION=3      EQUATION=3
# Line 96  class PDECoefficientType: Line 49  class PDECoefficientType:
49         self.what=where         self.what=where
50         self.pattern=pattern         self.pattern=pattern
51         self.altering=altering         self.altering=altering
52           self.resetValue()
53    
54        def resetValue(self):
55           """
56           @brief resets coefficient value to default
57           """
58           self.value=escript.Data()
59    
60      def getFunctionSpace(self,domain):      def getFunctionSpace(self,domain):
61         """         """
# Line 108  class PDECoefficientType: Line 68  class PDECoefficientType:
68         elif self.what==self.CONTACT: return escript.FunctionOnContactZero(domain)         elif self.what==self.CONTACT: return escript.FunctionOnContactZero(domain)
69         elif self.what==self.CONTINUOUS: return escript.ContinuousFunction(domain)         elif self.what==self.CONTINUOUS: return escript.ContinuousFunction(domain)
70    
71        def getValue(self):
72           """
73           @brief returns the value of the coefficient:
74           """
75           return self.value
76        
77        def setValue(self,newValue):
78           """
79           @brief set the value of the coefficient to new value
80           """
81           self.value=newValue
82        
83      def isAlteringOperator(self):      def isAlteringOperator(self):
84          """          """
85      @brief return true if the operator of the PDE is changed when the coefficient is changed      @brief return true if the operator of the PDE is changed when the coefficient is changed
# Line 168  class PDECoefficientType: Line 140  class PDECoefficientType:
140                  s=s+(dim,)                  s=s+(dim,)
141          return s          return s
142    
 _PDECoefficientTypes={  
 "A"         : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,PDECoefficientType.DIM,PDECoefficientType.SOLUTION,PDECoefficientType.DIM),PDECoefficientType.OPERATOR),  
 "B"         : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,PDECoefficientType.DIM,PDECoefficientType.SOLUTION),PDECoefficientType.OPERATOR),  
 "C"         : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,PDECoefficientType.SOLUTION,PDECoefficientType.DIM),PDECoefficientType.OPERATOR),  
 "D"         : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,PDECoefficientType.SOLUTION),PDECoefficientType.OPERATOR),  
 "X"         : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,PDECoefficientType.DIM),PDECoefficientType.RIGHTHANDSIDE),  
 "Y"         : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,),PDECoefficientType.RIGHTHANDSIDE),  
 "d"         : PDECoefficientType(PDECoefficientType.BOUNDARY,(PDECoefficientType.EQUATION,PDECoefficientType.SOLUTION),PDECoefficientType.OPERATOR),  
 "y"         : PDECoefficientType(PDECoefficientType.BOUNDARY,(PDECoefficientType.EQUATION,),PDECoefficientType.RIGHTHANDSIDE),  
 "d_contact" : PDECoefficientType(PDECoefficientType.CONTACT,(PDECoefficientType.EQUATION,PDECoefficientType.SOLUTION),PDECoefficientType.OPERATOR),  
 "y_contact" : PDECoefficientType(PDECoefficientType.CONTACT,(PDECoefficientType.EQUATION,),PDECoefficientType.RIGHTHANDSIDE),  
 "r"         : PDECoefficientType(PDECoefficientType.CONTINUOUS,(PDECoefficientType.EQUATION,),PDECoefficientType.RIGHTHANDSIDE),  
 "q"         : PDECoefficientType(PDECoefficientType.CONTINUOUS,(PDECoefficientType.SOLUTION,),PDECoefficientType.BOTH),  
 }  
   
143  class LinearPDE:  class LinearPDE:
144     """     """
145     @brief Class to define a linear PDE     @brief Class to handel a linear PDE
146        
147     class to define a linear PDE of the form     class to define a linear PDE of the form
148    
# Line 204  class LinearPDE: Line 161  class LinearPDE:
161           u_i=r_i where q_i>0           u_i=r_i where q_i>0
162    
163     """     """
164       TOL=1.e-13
165     DEFAULT_METHOD=util.DEFAULT_METHOD     DEFAULT_METHOD=util.DEFAULT_METHOD
166     DIRECT=util.DIRECT     DIRECT=util.DIRECT
167     CHOLEVSKY=util.CHOLEVSKY     CHOLEVSKY=util.CHOLEVSKY
# Line 214  class LinearPDE: Line 172  class LinearPDE:
172     SSOR=util.SSOR     SSOR=util.SSOR
173     GMRES=util.GMRES     GMRES=util.GMRES
174     PRES20=util.PRES20     PRES20=util.PRES20
175       ILU0=util.ILU0
176       JACOBI=util.JACOBI
177    
178     def __init__(self,domain,numEquations=0,numSolutions=0):     def __init__(self,domain,numEquations=0,numSolutions=0):
179       """       """
# Line 221  class LinearPDE: Line 181  class LinearPDE:
181    
182       @param args       @param args
183       """       """
184         # COEFFICIENTS can be overwritten by subclasses:
185         self.COEFFICIENTS={
186           "A"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM,PDECoefficient.SOLUTION,PDECoefficient.DIM),PDECoefficient.OPERATOR),
187           "B"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),
188           "C"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION,PDECoefficient.DIM),PDECoefficient.OPERATOR),
189           "D"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),
190           "X"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM),PDECoefficient.RIGHTHANDSIDE),
191           "Y"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
192           "d"         : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),
193           "y"         : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
194           "d_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),
195           "y_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
196           "r"         : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
197           "q"         : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.SOLUTION,),PDECoefficient.BOTH)}
198    
199       # initialize attributes       # initialize attributes
200       self.__debug=None       self.__debug=None
# Line 246  class LinearPDE: Line 220  class LinearPDE:
220       self.__sym=False       self.__sym=False
221       self.__lumping=False       self.__lumping=False
222    
223       def createCoefficient(self, name):
224         """
225         @brief create a data object corresponding to coefficient name
226         @param name
227         """
228         return escript.Data(shape = getShapeOfCoefficient(name), \
229                             what = getFunctionSpaceForCoefficient(name))
230    
231       def __del__(self):
232         pass
233    
234     def getCoefficient(self,name):     def getCoefficient(self,name):
235       """       """
236       @brief return the value of the coefficient name       @brief return the value of the parameter name
237    
238       @param name       @param name
239       """       """
240       return self.__coefficient[name]       return self.COEFFICIENTS[name].getValue()
241    
242     def setValue(self,**coefficients):     def getCoefficientOfPDE(self,name):
243         """
244         @brief return the value of the coefficient name of the general PDE. This method is called by the assembling routine
245                it can be overwritten to map coefficients of a particualr PDE to the general PDE.
246         @param name
247         """
248         return self.getCoefficient(name)
249    
250       def hasCoefficient(self,name):
251        """        """
252        @brief sets new values to coefficients        @brief return true if name is the name of a coefficient
253    
254        @param coefficients        @param name
255        """        """
256        self._setValue(**coefficients)        return self.COEFFICIENTS.has_key(name)
         
257    
258     def _setValue(self,**coefficients):     def getFunctionSpaceForEquation(self):
259         """
260         @brief return true if the test functions should use reduced order
261         """
262         return self.__row_function_space
263    
264       def getFunctionSpaceForSolution(self):
265         """
266         @brief return true if the interpolation of the solution should use reduced order
267         """
268         return self.__column_function_space
269    
270       def setValue(self,**coefficients):
271        """        """
272        @brief sets new values to coefficients        @brief sets new values to coefficients
273    
274        @param coefficients        @param coefficients
275        """        """
276          self._setValue(**coefficients)
277                
       # get the dictionary of the coefficinets been altered:  
       alteredCoefficients={}  
       for i,d in coefficients.iteritems():  
          if self.hasCoefficient(i):  
             if d == None:  
                 alteredCoefficients[i]=escript.Data()  
             elif isinstance(d,escript.Data):  
                 if d.isEmpty():  
                   alteredCoefficients[i]=escript.Data()  
                 else:  
                   alteredCoefficients[i]=escript.Data(d,self.getFunctionSpaceOfCoefficient(i))  
             else:  
                 alteredCoefficients[i]=escript.Data(d,self.getFunctionSpaceOfCoefficient(i))  
          else:  
             raise ValueError,"Attempt to set undefined coefficient %s"%i  
       # if numEquations and numSolutions is undefined we try identify their values based on the coefficients:  
       if self.__numEquations<1 or self.__numSolutions<1:  
             numEquations,numSolutions=identifyNumEquationsAndSolutions(self.getDomain().getDim(),alteredCoefficients)  
             if self.__numEquations<1 and numEquations>0: self.__numEquations=numEquations  
             if self.__numSolutions<1 and numSolutions>0: self.__numSolutions=numSolutions  
             if self.debug() and self.__numEquations>0: print "PDE Debug: identified number of equations is ",self.__numEquations  
             if self.debug() and self.__numSolutions>0: print "PDE Debug: identified number of solutions is ",self.__numSolutions  
   
       # now we check the shape of the coefficient if numEquations and numSolutions are set:  
       if  self.__numEquations>0 and  self.__numSolutions>0:  
          for i in self.__coefficient.iterkeys():  
              if alteredCoefficients.has_key(i) and not alteredCoefficients[i].isEmpty():  
                  if not self.getShapeOfCoefficient(i)==alteredCoefficients[i].getShape():  
                     raise ValueError,"Expected shape for coefficient %s is %s but actual shape is %s."%(i,self.getShapeOfCoefficient(i),alteredCoefficients[i].getShape())  
              else:  
                  if not self.__coefficient[i].isEmpty():  
                     if not self.getShapeOfCoefficient(i)==self.__coefficient[i].getShape():  
                        raise ValueError,"Expected shape for coefficient %s is %s but actual shape is %s."%(i,self.getShapeOfCoefficient(i),self.__coefficient[i].getShape())  
       # overwrite new values:  
       for i,d in alteredCoefficients.iteritems():  
          if self.debug(): print "PDE Debug: Coefficient %s has been altered."%i  
          self.__coefficient[i]=d  
          self.alteredCoefficient(i)  
   
       # reset the HomogeneousConstraintFlag:  
       self.__setHomogeneousConstraintFlag()  
       if not "q" in alteredCoefficients and not self.__homogeneous_constraint: self.__rebuildSystem()  
278    
279     def cleanCoefficients(self):     def cleanCoefficients(self):
280       """       """
281       @brief resets all coefficients to default values.       @brief resets all coefficients to default values.
282       """       """
283       self.__coefficient={}       for i in self.COEFFICIENTS.iterkeys():
284       for i in _PDECoefficientTypes.iterkeys():           self.COEFFICIENTS[i].resetValue()
285           self.__coefficient[i]=escript.Data()  
286       def createNewCoefficient(self,name):
287         """
288         @brief returns a new coefficient appropriate for coefficient name:
289         """
290         return escript.Data(0,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))
291          
292    
293     def getShapeOfCoefficient(self,name):     def getShapeOfCoefficient(self,name):
294       """       """
# Line 328  class LinearPDE: Line 297  class LinearPDE:
297       @param name       @param name
298       """       """
299       if self.hasCoefficient(name):       if self.hasCoefficient(name):
300          return _PDECoefficientTypes[name].buildShape(self.getNumEquations(),self.getNumSolutions(),self.getDomain().getDim())          return self.COEFFICIENTS[name].buildShape(self.getNumEquations(),self.getNumSolutions(),self.getDomain().getDim())
301       else:       else:
302          raise ValueError,"Solution coefficient %s requested"%name          raise ValueError,"Solution coefficient %s requested"%name
303    
304     def getFunctionSpaceOfCoefficient(self,name):     def getFunctionSpaceForCoefficient(self,name):
305       """       """
306       @brief return the atoms of the coefficient name       @brief return the atoms of the coefficient name
307    
308       @param name       @param name
309       """       """
310       if self.hasCoefficient(name):       if self.hasCoefficient(name):
311          return _PDECoefficientTypes[name].getFunctionSpace(self.getDomain())          return self.COEFFICIENTS[name].getFunctionSpace(self.getDomain())
312       else:       else:
313          raise ValueError,"Solution coefficient %s requested"%name          raise ValueError,"Solution coefficient %s requested"%name
314    
# Line 350  class LinearPDE: Line 319  class LinearPDE:
319       @param name       @param name
320       """       """
321       if self.hasCoefficient(name):       if self.hasCoefficient(name):
322          if _PDECoefficientTypes[name].isAlteringOperator(): self.__rebuildOperator()          if self.COEFFICIENTS[name].isAlteringOperator(): self.__rebuildOperator()
323          if _PDECoefficientTypes[name].isAlteringRightHandSide(): self.__rebuildRightHandSide()          if self.COEFFICIENTS[name].isAlteringRightHandSide(): self.__rebuildRightHandSide()
324       else:       else:
325          raise ValueError,"Solution coefficient %s requested"%name          raise ValueError,"unknown coefficient %s requested"%name
   
    def __setHomogeneousConstraintFlag(self):  
       """  
       @brief checks if the constraints are homogeneous and sets self.__homogeneous_constraint accordingly.  
       """  
       self.__homogeneous_constraint=True  
       q=self.getCoefficient("q")  
       r=self.getCoefficient("r")  
       if not q.isEmpty() and not r.isEmpty():  
          print (q*r).Lsup(), 1.e-13*r.Lsup()  
          if (q*r).Lsup()>=1.e-13*r.Lsup(): self.__homogeneous_constraint=False  
       if self.debug():  
            if self.__homogeneous_constraint:  
                print "PDE Debug: Constraints are homogeneous."  
            else:  
                print "PDE Debug: Constraints are inhomogeneous."  
   
   
    def hasCoefficient(self,name):  
       """  
       @brief return true if name is the name of a coefficient  
   
       @param name  
       """  
       return self.__coefficient.has_key(name)  
   
    def getFunctionSpaceForEquation(self):  
      """  
      @brief return true if the test functions should use reduced order  
      """  
      return self.__row_function_space  
   
    def getFunctionSpaceForSolution(self):  
      """  
      @brief return true if the interpolation of the solution should use reduced order  
      """  
      return self.__column_function_space  
326    
327     # ===== debug ==============================================================     # ===== debug ==============================================================
328     def setDebugOn(self):     def setDebugOn(self):
# Line 417  class LinearPDE: Line 349  class LinearPDE:
349        @brief indicates to use matrix lumping        @brief indicates to use matrix lumping
350        """        """
351        if not self.isUsingLumping():        if not self.isUsingLumping():
          raise SystemError,"Lumping is not working yet! Talk to the experts"  
352           if self.debug() : print "PDE Debug: lumping is set on"           if self.debug() : print "PDE Debug: lumping is set on"
353           self.__rebuildOperator()           self.__rebuildOperator()
354           self.__lumping=True           self.__lumping=True
# Line 604  class LinearPDE: Line 535  class LinearPDE:
535       else:       else:
536          self.setReducedOrderForEquationOff()          self.setReducedOrderForEquationOff()
537                                                                                                                                                                                                                                                                                                                        
   
538     # ==== initialization =====================================================================     # ==== initialization =====================================================================
539     def __makeNewOperator(self):     def __makeNewOperator(self):
540         """         """
# Line 638  class LinearPDE: Line 568  class LinearPDE:
568             if self.debug() : print "PDE Debug: New operator allocated"             if self.debug() : print "PDE Debug: New operator allocated"
569         else:         else:
570             self.__operator.setValue(0.)             self.__operator.setValue(0.)
571               self.__operator.resetSolver()
572             if self.debug() : print "PDE Debug: Operator reset to zero"             if self.debug() : print "PDE Debug: Operator reset to zero"
573         return self.__operator         return self.__operator
574    
# Line 654  class LinearPDE: Line 585  class LinearPDE:
585             if self.debug() : print "PDE Debug: Right hand side reset to zero"             if self.debug() : print "PDE Debug: Right hand side reset to zero"
586         return  self.__righthandside         return  self.__righthandside
587    
588       #============ some serivice functions  =====================================================
589       def getDomain(self):
590         """
591         @brief returns the domain of the PDE
592         """
593         return self.__domain
594    
595       def getDim(self):
596         """
597         @brief returns the spatial dimension of the PDE
598         """
599         return self.getDomain().getDim()
600    
601       def getNumEquations(self):
602         """
603         @brief returns the number of equations
604         """
605         if self.__numEquations>0:
606             return self.__numEquations
607         else:
608             raise ValueError,"Number of equations is undefined. Please specify argument numEquations."
609    
610       def getNumSolutions(self):
611         """
612         @brief returns the number of unknowns
613         """
614         if self.__numSolutions>0:
615            return self.__numSolutions
616         else:
617            raise ValueError,"Number of solution is undefined. Please specify argument numSolutions."
618    
619    
620       def checkSymmetry(self,verbose=True):
621          """
622          @brief returns if the Operator is symmetric. This is a very expensive operation!!! The symmetry flag is not altered.
623          """
624          verbose=verbose or self.debug()
625          out=True
626          if self.getNumSolutions()!=self.getNumEquations():
627             if verbose: print "non-symmetric PDE because of different number of equations and solutions"
628             out=False
629          else:
630             A=self.getCoefficientOfPDE("A")
631             if not A.isEmpty():
632                tol=util.Lsup(A)*self.TOL
633                if self.getNumSolutions()>1:
634                   for i in range(self.getNumEquations()):
635                      for j in range(self.getDim()):
636                         for k in range(self.getNumSolutions()):
637                            for l in range(self.getDim()):
638                                if util.Lsup(A[i,j,k,l]-A[k,l,i,j])>tol:
639                                   if verbose: print "non-symmetric PDE because A[%d,%d,%d,%d]!=A[%d,%d,%d,%d]"%(i,j,k,l,k,l,i,j)
640                                   out=False
641                else:
642                   for j in range(self.getDim()):
643                      for l in range(self.getDim()):
644                         if util.Lsup(A[j,l]-A[l,j])>tol:
645                            if verbose: print "non-symmetric PDE because A[%d,%d]!=A[%d,%d]"%(j,l,l,j)
646                            out=False
647             B=self.getCoefficientOfPDE("B")
648             C=self.getCoefficientOfPDE("C")
649             if B.isEmpty() and not C.isEmpty():
650                if verbose: print "non-symmetric PDE because B is not present but C is"
651                out=False
652             elif not B.isEmpty() and C.isEmpty():
653                if verbose: print "non-symmetric PDE because C is not present but B is"
654                out=False
655             elif not B.isEmpty() and not C.isEmpty():
656                tol=(util.Lsup(B)+util.Lsup(C))*self.TOL/2.
657                if self.getNumSolutions()>1:
658                   for i in range(self.getNumEquations()):
659                       for j in range(self.getDim()):
660                          for k in range(self.getNumSolutions()):
661                             if util.Lsup(B[i,j,k]-C[k,i,j])>tol:
662                                  if verbose: print "non-symmetric PDE because B[%d,%d,%d]!=C[%d,%d,%d]"%(i,j,k,k,i,j)
663                                  out=False
664                else:
665                   for j in range(self.getDim()):
666                      if util.Lsup(B[j]-C[j])>tol:
667                         if verbose: print "non-symmetric PDE because B[%d]!=C[%d]"%(j,j)
668                         out=False
669             if self.getNumSolutions()>1:
670               D=self.getCoefficientOfPDE("D")
671               if not D.isEmpty():
672                 tol=util.Lsup(D)*self.TOL
673                 for i in range(self.getNumEquations()):
674                    for k in range(self.getNumSolutions()):
675                      if util.Lsup(D[i,k]-D[k,i])>tol:
676                          if verbose: print "non-symmetric PDE because D[%d,%d]!=D[%d,%d]"%(i,k,k,i)
677                          out=False
678                
679          return out
680    
681       def getFlux(self,u):
682           """
683           @brief returns the flux J_ij for a given u
684    
685                J_ij=A_{ijkl}u_{k,l}+B_{ijk}u_k-X_{ij}
686    
687           @param u argument of the operator
688    
689           """
690           raise SystemError,"getFlux is not implemented yet"
691           return None
692    
693       def applyOperator(self,u):
694           """
695           @brief applies the operator of the PDE to a given solution u in weak from
696    
697           @param u argument of the operator
698    
699           """
700           return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())
701                                                                                                                                                              
702       def getResidual(self,u):
703           """
704           @brief return the residual of u in the weak from
705    
706           @param u
707           """
708           return self.applyOperator(u)-self.getRightHandSide()
709    
710       def _setValue(self,**coefficients):
711          """
712          @brief sets new values to coefficient
713    
714          @param coefficients
715          """
716          # check if the coefficients are  legal:
717          for i in coefficients.iterkeys():
718             if not self.hasCoefficient(i):
719                raise ValueError,"Attempt to set unknown coefficient %s"%i
720          # if the number of unknowns or equations is still unknown we try to estimate them:
721          if self.__numEquations<1 or self.__numSolutions<1:
722             for i,d in coefficients.iteritems():
723                if hasattr(d,"shape"):
724                    s=d.shape
725                elif hasattr(d,"getShape"):
726                    s=d.getShape()
727                else:
728                    s=numarray.array(d).shape
729                if s!=None:
730                    # get number of equations and number of unknowns:
731                    res=self.COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(s,self.getDim())
732                    if res==None:
733                        raise ValueError,"Illegal shape %s of coefficient %s"%(s,i)
734                    else:
735                        if self.__numEquations<1: self.__numEquations=res[0]
736                        if self.__numSolutions<1: self.__numSolutions=res[1]
737          if self.__numEquations<1: raise ValueError,"unidententified number of equations"
738          if self.__numSolutions<1: raise ValueError,"unidententified number of solutions"
739          # now we check the shape of the coefficient if numEquations and numSolutions are set:
740          for i,d in coefficients.iteritems():
741            if d==None:
742                 d2=escript.Data()
743            elif isinstance(d,escript.Data):
744                 if d.isEmpty():
745                    d2=d
746                 else:
747                    d2=escript.Data(d,self.getFunctionSpaceForCoefficient(i))
748            else:
749                  d2=escript.Data(d,self.getFunctionSpaceForCoefficient(i))
750            if not d2.isEmpty():
751               if not self.getShapeOfCoefficient(i)==d2.getShape():
752                   raise ValueError,"Expected shape for coefficient %s is %s but actual shape is %s."%(i,self.getShapeOfCoefficient(i),d2.getShape())
753            # overwrite new values:
754            if self.debug(): print "PDE Debug: Coefficient %s has been altered."%i
755            self.COEFFICIENTS[i].setValue(d2)
756            self.alteredCoefficient(i)
757          
758          # reset the HomogeneousConstraintFlag:
759          self.__setHomogeneousConstraintFlag()
760          if len(coefficients)>0 and not self.isUsingLumping() and not self.__homogeneous_constraint: self.__rebuildSystem()
761    
762       def __setHomogeneousConstraintFlag(self):
763          """
764          @brief checks if the constraints are homogeneous and sets self.__homogeneous_constraint accordingly.
765          """
766          self.__homogeneous_constraint=True
767          q=self.getCoefficientOfPDE("q")
768          r=self.getCoefficientOfPDE("r")
769          if not q.isEmpty() and not r.isEmpty():
770             if (q*r).Lsup()>=1.e-13*r.Lsup(): self.__homogeneous_constraint=False
771          if self.debug():
772               if self.__homogeneous_constraint:
773                   print "PDE Debug: Constraints are homogeneous."
774               else:
775                   print "PDE Debug: Constraints are inhomogeneous."
776    
777    
778     # ==== rebuild switches =====================================================================     # ==== rebuild switches =====================================================================
779     def __rebuildSolution(self,deep=False):     def __rebuildSolution(self,deep=False):
780         """         """
# Line 661  class LinearPDE: Line 782  class LinearPDE:
782         """         """
783         if self.__solution_isValid and self.debug() : print "PDE Debug: PDE has to be resolved."         if self.__solution_isValid and self.debug() : print "PDE Debug: PDE has to be resolved."
784         self.__solution_isValid=False         self.__solution_isValid=False
785         if deep: self.__solution=escript.Data(deep)         if deep: self.__solution=escript.Data()
786    
787    
788     def __rebuildOperator(self,deep=False):     def __rebuildOperator(self,deep=False):
# Line 680  class LinearPDE: Line 801  class LinearPDE:
801         if self.__righthandside_isValid and self.debug() : print "PDE Debug: Right hand side has to be rebuilt."         if self.__righthandside_isValid and self.debug() : print "PDE Debug: Right hand side has to be rebuilt."
802         self.__rebuildSolution(deep)         self.__rebuildSolution(deep)
803         self.__righthandside_isValid=False         self.__righthandside_isValid=False
        if not self.__homogeneous_constraint: self.__rebuildOperator()  
804         if deep: self.__righthandside=escript.Data()         if deep: self.__righthandside=escript.Data()
805    
806     def __rebuildSystem(self,deep=False):     def __rebuildSystem(self,deep=False):
# Line 702  class LinearPDE: Line 822  class LinearPDE:
822           self.__rebuildOperator(deep=True)           self.__rebuildOperator(deep=True)
823    
824     #============ assembling =======================================================     #============ assembling =======================================================
825     def __copyConstraint(self,u):     def __copyConstraint(self):
826        """        """
827        @brief copies the constrint condition into u        @brief copies the constrint condition into u
828        """        """
829        q=self.getCoefficient("q")        if not self.__righthandside.isEmpty():
830        r=self.getCoefficient("r")           q=self.getCoefficientOfPDE("q")
831        if not q.isEmpty():           r=self.getCoefficientOfPDE("r")
832            if r.isEmpty():           if not q.isEmpty():
833               r2=escript.Data(0,u.getShape(),u.getFunctionSpace())               if r.isEmpty():
834            else:                  r2=escript.Data(0,self.__righthandside.getShape(),self.__righthandside.getFunctionSpace())
835               r2=escript.Data(r,u.getFunctionSpace())               else:
836            u.copyWithMask(r2,escript.Data(q,u.getFunctionSpace()))                  r2=escript.Data(r,self.__righthandside.getFunctionSpace())
837                 self.__righthandside.copyWithMask(r2,escript.Data(q,self.__righthandside.getFunctionSpace()))
838    
839     def __applyConstraint(self,rhs_update=True):     def __applyConstraint(self):
840         """         """
841         @brief applies the constraints  defined by q and r to the system         @brief applies the constraints defined by q and r to the system
842         """         """
843         q=self.getCoefficient("q")         q=self.getCoefficientOfPDE("q")
844         r=self.getCoefficient("r")         r=self.getCoefficientOfPDE("r")
845         if not q.isEmpty() and not self.__operator.isEmpty():         if not q.isEmpty() and not self.__operator.isEmpty():
846            # q is the row and column mask to indicate where constraints are set:            # q is the row and column mask to indicate where constraints are set:
847            row_q=escript.Data(q,self.getFunctionSpaceForEquation())            row_q=escript.Data(q,self.getFunctionSpaceForEquation())
# Line 731  class LinearPDE: Line 852  class LinearPDE:
852            else:            else:
853               r_s=escript.Data(r,self.getFunctionSpaceForSolution())               r_s=escript.Data(r,self.getFunctionSpaceForSolution())
854            u.copyWithMask(r_s,col_q)            u.copyWithMask(r_s,col_q)
855            if not self.__righthandside.isEmpty() and rhs_update: self.__righthandside-=self.__operator*u            if self.isUsingLumping():
856            self.__operator.nullifyRowsAndCols(row_q,col_q,1.)               self.__operator.copyWithMask(escript.Data(1,q.getShape(),self.getFunctionSpaceForEquation()),row_q)
857              else:
858                 if not self.__righthandside.isEmpty(): self.__righthandside-=self.__operator*u
859                 self.__operator.nullifyRowsAndCols(row_q,col_q,1.)
860    
861     def getOperator(self):     def getSystem(self):
862         """         """
863         @brief returns the operator of the PDE         @brief return the operator and right hand side of the PDE
864         """         """
865         if not self.__operator_isValid:         if not self.__operator_isValid or not self.__righthandside_isValid:
866             # some Constraints are applying for a lumpled stifness matrix:            if self.isUsingLumping():
867             if self.isUsingLumping():                if not self.__operator_isValid:
868                if self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():                   if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():
869                         raise TypeError,"Lumped matrix requires same order for equations and unknowns"                         raise TypeError,"Lumped matrix requires same order for equations and unknowns"
870                if not self.getCoefficient("A").isEmpty():                   if not self.getCoefficientOfPDE("A").isEmpty():
871                         raise Warning,"Lumped matrix does not allow coefficient A"                            raise Warning,"Lumped matrix does not allow coefficient A"
872                if not self.getCoefficient("B").isEmpty():                   if not self.getCoefficientOfPDE("B").isEmpty():
873                         raise Warning,"Lumped matrix does not allow coefficient B"                            raise Warning,"Lumped matrix does not allow coefficient B"
874                if not self.getCoefficient("C").isEmpty():                   if not self.getCoefficientOfPDE("C").isEmpty():
875                         raise Warning,"Lumped matrix does not allow coefficient C"                            raise Warning,"Lumped matrix does not allow coefficient C"
876                if not self.getCoefficient("D").isEmpty():                   if self.debug() : print "PDE Debug: New lumped operator is built."
877                         raise Warning,"Lumped matrix does not allow coefficient D"                   mat=self.__makeNewOperator()
878                if self.debug() : print "PDE Debug: New lumped operator is built."                   self.getDomain().addPDEToSystem(mat,escript.Data(), \
879                mat=self.__makeNewOperator(self)                             self.getCoefficientOfPDE("A"), \
880             else:                             self.getCoefficientOfPDE("B"), \
881                if self.debug() : print "PDE Debug: New operator is built."                             self.getCoefficientOfPDE("C"), \
882                mat=self.__getFreshOperator()                             self.getCoefficientOfPDE("D"), \
883                               escript.Data(), \
884             self.getDomain().addPDEToSystem(mat,escript.Data(), \                             escript.Data(), \
885                          self.getCoefficient("A"), \                             self.getCoefficientOfPDE("d"), \
886                          self.getCoefficient("B"), \                             escript.Data(),\
887                          self.getCoefficient("C"), \                             self.getCoefficientOfPDE("d_contact"), \
888                          self.getCoefficient("D"), \                             escript.Data())
889                          escript.Data(), \                   self.__operator=mat*escript.Data(1,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)
890                          escript.Data(), \                   self.__applyConstraint()
891                          self.getCoefficient("d"), \                   self.__operator_isValid=True
892                          escript.Data(),\                if not self.__righthandside_isValid:
893                          self.getCoefficient("d_contact"), \                   if self.debug() : print "PDE Debug: New right hand side is built."
894                          escript.Data())                   self.getDomain().addPDEToRHS(self.__getFreshRightHandSide(), \
895             if self.isUsingLumping():                                 self.getCoefficientOfPDE("X"), \
896                self.__operator=mat*escript.Data(1,(self.getNumSolutions(),),self.getFunctionSpaceOfSolution(),True)                                 self.getCoefficientOfPDE("Y"),\
897             else:                                 self.getCoefficientOfPDE("y"),\
898                self.__applyConstraint(rhs_update=False)                                 self.getCoefficientOfPDE("y_contact"))
899             self.__operator_isValid=True                   self.__copyConstraint()
900         return self.__operator                   self.__righthandside_isValid=True
901              else:
902     def getRightHandSide(self,ignoreConstraint=False):               if not self.__operator_isValid and not self.__righthandside_isValid:
903                     if self.debug() : print "PDE Debug: New system is built."
904                     self.getDomain().addPDEToSystem(self.__getFreshOperator(),self.__getFreshRightHandSide(), \
905                                   self.getCoefficientOfPDE("A"), \
906                                   self.getCoefficientOfPDE("B"), \
907                                   self.getCoefficientOfPDE("C"), \
908                                   self.getCoefficientOfPDE("D"), \
909                                   self.getCoefficientOfPDE("X"), \
910                                   self.getCoefficientOfPDE("Y"), \
911                                   self.getCoefficientOfPDE("d"), \
912                                   self.getCoefficientOfPDE("y"), \
913                                   self.getCoefficientOfPDE("d_contact"), \
914                                   self.getCoefficientOfPDE("y_contact"))
915                     self.__applyConstraint()
916                     self.__copyConstraint()
917                     self.__operator_isValid=True
918                     self.__righthandside_isValid=True
919                 elif not self.__righthandside_isValid:
920                     if self.debug() : print "PDE Debug: New right hand side is built."
921                     self.getDomain().addPDEToRHS(self.__getFreshRightHandSide(), \
922                                   self.getCoefficientOfPDE("X"), \
923                                   self.getCoefficientOfPDE("Y"),\
924                                   self.getCoefficientOfPDE("y"),\
925                                   self.getCoefficientOfPDE("y_contact"))
926                     self.__copyConstraint()
927                     self.__righthandside_isValid=True
928                 elif not self.__operator_isValid:
929                     if self.debug() : print "PDE Debug: New operator is built."
930                     self.getDomain().addPDEToSystem(self.__getFreshOperator(),escript.Data(), \
931                                self.getCoefficientOfPDE("A"), \
932                                self.getCoefficientOfPDE("B"), \
933                                self.getCoefficientOfPDE("C"), \
934                                self.getCoefficientOfPDE("D"), \
935                                escript.Data(), \
936                                escript.Data(), \
937                                self.getCoefficientOfPDE("d"), \
938                                escript.Data(),\
939                                self.getCoefficientOfPDE("d_contact"), \
940                                escript.Data())
941                     self.__applyConstraint()
942                     self.__operator_isValid=True
943           return (self.__operator,self.__righthandside)
944       def getOperator(self):
945         """         """
946         @brief returns the right hand side of the PDE         @brief returns the operator of the PDE
   
        @param ignoreConstraint  
947         """         """
948         if not self.__righthandside_isValid:         return self.getSystem()[0]
            if self.debug() : print "PDE Debug: New right hand side is built."  
            self.getDomain().addPDEToRHS(self.__getFreshRightHandSide(), \  
                          self.getCoefficient("X"), \  
                          self.getCoefficient("Y"),\  
                          self.getCoefficient("y"),\  
                          self.getCoefficient("y_contact"))  
            self.__righthandside_isValid=True  
            if ignoreConstraint: self.__copyConstraint(self.__righthandside)  
        return self.__righthandside  
949    
950     def getSystem(self):     def getRightHandSide(self):
951         """         """
952         @brief         @brief returns the right hand side of the PDE
953         """         """
954         if not self.__operator_isValid and not self.__righthandside_isValid:         return self.getSystem()[1]
           if self.debug() : print "PDE Debug: New PDE is built."  
           if self.isUsingLumping():  
               self.getRightHandSide(ignoreConstraint=True)  
               self.getOperator()  
           else:  
               self.getDomain().addPDEToSystem(self.__getFreshOperator(),self.__getFreshRightHandSide(), \  
                             self.getCoefficient("A"), \  
                             self.getCoefficient("B"), \  
                             self.getCoefficient("C"), \  
                             self.getCoefficient("D"), \  
                             self.getCoefficient("X"), \  
                             self.getCoefficient("Y"), \  
                             self.getCoefficient("d"), \  
                             self.getCoefficient("y"), \  
                             self.getCoefficient("d_contact"), \  
                             self.getCoefficient("y_contact"))  
           self.__operator_isValid=True  
           self.__righthandside_isValid=True  
           self.__applyConstraint()  
           self.__copyConstraint(self.__righthandside)  
        elif not self.__operator_isValid:  
           self.getOperator()  
        elif not self.__righthandside_isValid:  
           self.getRightHandSide()  
        return (self.__operator,self.__righthandside)  
955    
956     def solve(self,**options):     def solve(self,**options):
957        """        """
# Line 832  class LinearPDE: Line 962  class LinearPDE:
962        mat,f=self.getSystem()        mat,f=self.getSystem()
963        if self.isUsingLumping():        if self.isUsingLumping():
964           out=f/mat           out=f/mat
          self.__copyConstraint(out)  
965        else:        else:
966           options[util.TOLERANCE_KEY]=self.getTolerance()           options[util.TOLERANCE_KEY]=self.getTolerance()
967           options[util.METHOD_KEY]=self.getSolverMethod()           options[util.METHOD_KEY]=self.getSolverMethod()
# Line 852  class LinearPDE: Line 981  class LinearPDE:
981             self.__solution=self.solve(**options)             self.__solution=self.solve(**options)
982             self.__solution_isValid=True             self.__solution_isValid=True
983         return self.__solution         return self.__solution
    #============ some serivice functions  =====================================================  
    def getDomain(self):  
      """  
      @brief returns the domain of the PDE  
      """  
      return self.__domain  
984    
    def getDim(self):  
      """  
      @brief returns the spatial dimension of the PDE  
      """  
      return self.getDomain().getDim()  
985    
    def getNumEquations(self):  
      """  
      @brief returns the number of equations  
      """  
      if self.__numEquations>0:  
          return self.__numEquations  
      else:  
          raise ValueError,"Number of equations is undefined. Please specify argument numEquations."  
986    
987     def getNumSolutions(self):  def ELMAN_RAMAGE(P): return (P-1.).wherePositive()*0.5*(1.-1./(P+1.e-15))
988       """  def SIMPLIFIED_BROOK_HUGHES(P):
989       @brief returns the number of unknowns           c=(P-3.).whereNegative()
990       """           return P/6.*c+1./2.*(1.-c)
991       if self.__numSolutions>0:  def HALF(P): return escript.Scalar(0.5,P.getFunctionSpace())
         return self.__numSolutions  
      else:  
         raise ValueError,"Number of solution is undefined. Please specify argument numSolutions."  
992    
993    
994     def checkSymmetry(self):  class AdvectivePDE(LinearPDE):
995        """     """
996        @brief returns if the Operator is symmetric. This is a very expensive operation!!! The symmetry flag is not altered.     @brief Class to handel a linear PDE domineated by advective terms:
997        """    
998        raise SystemError,"checkSymmetry is not implemented yet"     class to define a linear PDE of the form
999    
1000        return None       -(A_{ijkl}u_{k,l})_{,j} -(B_{ijk}u_k)_{,j} + C_{ikl}u_{k,l} +D_{ik}u_k = - (X_{ij})_{,j} + Y_i
1001    
1002     def getFlux(self,u):       with boundary conditons:
        """  
        @brief returns the flux J_ij for a given u  
1003    
1004              J_ij=A_{ijkl}u_{k,l}+B_{ijk}u_k-X_{ij}          n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_{ik}u_k = - n_j*X_{ij} + y_i
1005    
1006         @param u argument of the operator      and contact conditions
1007    
1008         """          n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_contact_{ik}[u_k] = - n_j*X_{ij} + y_contact_i
        raise SystemError,"getFlux is not implemented yet"  
        return None  
1009    
1010     def applyOperator(self,u):      and constraints:
        """  
        @brief applies the operator of the PDE to a given solution u in weak from  
1011    
1012         @param u argument of the operator           u_i=r_i where q_i>0
1013    
1014         """     """
1015         return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())     def __init__(self,domain,numEquations=0,numSolutions=0,xi=ELMAN_RAMAGE):
1016                                                                                                                                                                    LinearPDE.__init__(self,domain,numEquations,numSolutions)
1017     def getResidual(self,u):        self.__xi=xi
1018         """        self.__Xi=escript.Data()
1019         @brief return the residual of u in the weak from  
1020       def __calculateXi(self,peclet_factor,Z,h):
1021           Z_max=util.Lsup(Z)
1022           if Z_max>0.:
1023              return h*self.__xi(Z*peclet_factor)/(Z+Z_max*self.TOL)
1024           else:
1025              return 0.
1026    
1027       def setValue(self,**args):
1028           if "A" in args.keys()   or "B" in args.keys() or "C" in args.keys(): self.__Xi=escript.Data()
1029           self._setValue(**args)
1030              
1031       def getXi(self):
1032          if self.__Xi.isEmpty():
1033             B=self.getCoefficient("B")
1034             C=self.getCoefficient("C")
1035             A=self.getCoefficient("A")
1036             h=self.getDomain().getSize()
1037             self.__Xi=escript.Scalar(0.,self.getFunctionSpaceForCoefficient("A"))
1038             if not C.isEmpty() or not B.isEmpty():
1039                if not C.isEmpty() and not B.isEmpty():
1040                    Z2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))
1041                    if self.getNumEquations()>1:
1042                       if self.getNumSolutions()>1:
1043                          for i in range(self.getNumEquations()):
1044                             for k in range(self.getNumSolutions()):
1045                                for l in range(self.getDim()): Z2+=(C[i,k,l]-B[i,l,k])**2
1046                       else:
1047                          for i in range(self.getNumEquations()):
1048                             for l in range(self.getDim()): Z2+=(C[i,l]-B[i,l])**2
1049                    else:
1050                       if self.getNumSolutions()>1:
1051                          for k in range(self.getNumSolutions()):
1052                             for l in range(self.getDim()): Z2+=(C[k,l]-B[l,k])**2
1053                       else:
1054                          for l in range(self.getDim()): Z2+=(C[l]-B[l])**2
1055                    length_of_Z=util.sqrt(Z2)
1056                elif C.isEmpty():
1057                  length_of_Z=util.length(B)
1058                else:
1059                  length_of_Z=util.length(C)
1060    
1061                Z_max=util.Lsup(length_of_Z)
1062                if Z_max>0.:
1063                   length_of_A=util.length(A)
1064                   A_max=util.Lsup(length_of_A)
1065                   if A_max>0:
1066                        inv_A=1./(length_of_A+A_max*self.TOL)
1067                   else:
1068                        inv_A=1./self.TOL
1069                   peclet_number=length_of_Z*h/2*inv_A
1070                   xi=self.__xi(peclet_number)
1071                   self.__Xi=h*xi/(length_of_Z+Z_max*self.TOL)
1072                   print "@ preclet number = %e"%util.Lsup(peclet_number),util.Lsup(xi),util.Lsup(length_of_Z)
1073          return self.__Xi
1074          
1075    
1076       def getCoefficientOfPDE(self,name):
1077         """
1078         @brief return the value of the coefficient name of the general PDE
1079         @param name
1080         """
1081         if not self.getNumEquations() == self.getNumSolutions():
1082              raise ValueError,"AdvectivePDE expects the number of solution componets and the number of equations to be equal."
1083    
1084         if name == "A" :
1085             A=self.getCoefficient("A")
1086             B=self.getCoefficient("B")
1087             C=self.getCoefficient("C")
1088             if B.isEmpty() and C.isEmpty():
1089                Aout=A
1090             else:
1091                if A.isEmpty():
1092                   Aout=self.createNewCoefficient("A")
1093                else:
1094                   Aout=A[:]
1095                Xi=self.getXi()
1096                if self.getNumEquations()>1:
1097                    for i in range(self.getNumEquations()):
1098                       for j in range(self.getDim()):
1099                          for k in range(self.getNumSolutions()):
1100                             for l in range(self.getDim()):
1101                                if not C.isEmpty() and not B.isEmpty():
1102                                   for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*(C[p,i,j]-B[p,j,i])*(C[p,k,l]-B[p,l,k])
1103                                elif C.isEmpty():
1104                                   for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*B[p,j,i]*B[p,l,k]
1105                                else:
1106                                   for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*C[p,i,j]*C[p,k,l]
1107                else:
1108                    for j in range(self.getDim()):
1109                       for l in range(self.getDim()):
1110                          if not C.isEmpty() and not B.isEmpty():
1111                              Aout[j,l]+=Xi*(C[j]-B[j])*(C[l]-B[l])
1112                          elif C.isEmpty():
1113                              Aout[j,l]+=Xi*B[j]*B[l]
1114                          else:
1115                              Aout[j,l]+=Xi*C[j]*C[l]
1116             return Aout
1117         elif name == "B" :
1118             B=self.getCoefficient("B")
1119             C=self.getCoefficient("C")
1120             D=self.getCoefficient("D")
1121             if C.isEmpty() or D.isEmpty():
1122                Bout=B
1123             else:
1124                Xi=self.getXi()
1125                if B.isEmpty():
1126                    Bout=self.createNewCoefficient("B")
1127                else:
1128                    Bout=B[:]
1129                if self.getNumEquations()>1:
1130                   for k in range(self.getNumSolutions()):
1131                      for p in range(self.getNumEquations()):
1132                         tmp=Xi*D[p,k]
1133                         for i in range(self.getNumEquations()):
1134                            for j in range(self.getDim()):
1135                               Bout[i,j,k]+=tmp*C[p,i,j]
1136                else:
1137                   tmp=Xi*D
1138                   for j in range(self.getDim()): Bout[j]+=tmp*C[j]
1139             return Bout
1140         elif name == "C" :
1141             B=self.getCoefficient("B")
1142             C=self.getCoefficient("C")
1143             D=self.getCoefficient("D")
1144             if B.isEmpty() or D.isEmpty():
1145                Cout=C
1146             else:
1147                Xi=self.getXi()
1148                if C.isEmpty():
1149                    Cout=self.createNewCoefficient("C")
1150                else:
1151                    Cout=C[:]
1152                if self.getNumEquations()>1:
1153                   for k in range(self.getNumSolutions()):
1154                       for p in range(self.getNumEquations()):
1155                          tmp=Xi*D[p,k]
1156                          for i in range(self.getNumEquations()):
1157                            for l in range(self.getDim()):
1158                                     Cout[i,k,l]+=tmp*B[p,l,i]
1159                else:
1160                   tmp=Xi*D
1161                   for j in range(self.getDim()): Cout[j]+=tmp*B[j]
1162             return Cout
1163         elif name == "D" :
1164             return self.getCoefficient("D")
1165         elif name == "X" :
1166             X=self.getCoefficient("X")
1167             Y=self.getCoefficient("Y")
1168             B=self.getCoefficient("B")
1169             C=self.getCoefficient("C")
1170             if Y.isEmpty() or (B.isEmpty() and C.isEmpty()):
1171                Xout=X
1172             else:
1173                if X.isEmpty():
1174                    Xout=self.createNewCoefficient("X")
1175                else:
1176                    Xout=X[:]
1177                Xi=self.getXi()
1178                if self.getNumEquations()>1:
1179                     for p in range(self.getNumEquations()):
1180                        tmp=Xi*Y[p]
1181                        for i in range(self.getNumEquations()):
1182                           for j in range(self.getDim()):
1183                              if not C.isEmpty() and not B.isEmpty():
1184                                 Xout[i,j]+=tmp*(C[p,i,j]-B[p,j,i])
1185                              elif C.isEmpty():
1186                                 Xout[i,j]-=tmp*B[p,j,i]
1187                              else:
1188                                 Xout[i,j]+=tmp*C[p,i,j]
1189                else:
1190                     tmp=Xi*Y
1191                     for j in range(self.getDim()):
1192                        if not C.isEmpty() and not B.isEmpty():
1193                           Xout[j]+=tmp*(C[j]-B[j])
1194                        elif C.isEmpty():
1195                           Xout[j]-=tmp*B[j]
1196                        else:
1197                           Xout[j]+=tmp*C[j]
1198             return Xout
1199         elif name == "Y" :
1200             return self.getCoefficient("Y")
1201         elif name == "d" :
1202             return self.getCoefficient("d")
1203         elif name == "y" :
1204             return self.getCoefficient("y")
1205         elif name == "d_contact" :
1206             return self.getCoefficient("d_contact")
1207         elif name == "y_contact" :
1208             return self.getCoefficient("y_contact")
1209         elif name == "r" :
1210             return self.getCoefficient("r")
1211         elif name == "q" :
1212             return self.getCoefficient("q")
1213         else:
1214             raise SystemError,"unknown PDE coefficient %s",name
1215    
        @param u  
        """  
        return self.applyOperator(u)-self.getRightHandSide()  
1216    
1217  class Poisson(LinearPDE):  class Poisson(LinearPDE):
1218     """     """
# Line 939  class Poisson(LinearPDE): Line 1232  class Poisson(LinearPDE):
1232                                                                                                                                                                                                                                                                                                                            
1233     """     """
1234    
1235     def __init__(self,domain=None,f=escript.Data(),q=escript.Data()):     def __init__(self,domain,f=escript.Data(),q=escript.Data()):
1236         LinearPDE.__init__(self,domain=identifyDomain(domain,{"f":f, "q":q}))         LinearPDE.__init__(self,domain,1,1)
1237         self._setValue(A=numarray.identity(self.getDomain().getDim()))         self.COEFFICIENTS={
1238           "f"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1239           "q"         : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.BOTH)}
1240         self.setSymmetryOn()         self.setSymmetryOn()
1241         self.setValue(f,q)         self.setValue(f,q)
1242    
1243     def setValue(self,f=escript.Data(),q=escript.Data()):     def setValue(self,f=escript.Data(),q=escript.Data()):
1244         self._setValue(Y=f,q=q)         self._setValue(f=f,q=q)
1245    
1246                                                                                                                                                                 def getCoefficientOfPDE(self,name):
1247  # $Log$       """
1248  # Revision 1.3  2004/12/17 07:43:10  jgs       @brief return the value of the coefficient name of the general PDE
1249  # *** empty log message ***       @param name
1250  #       """
1251  # Revision 1.1.2.3  2004/12/16 00:12:34  gross       if name == "A" :
1252  # __init__ of LinearPDE does not accept any coefficients anymore           return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))
1253  #       elif name == "B" :
1254  # Revision 1.1.2.2  2004/12/14 03:55:01  jgs           return escript.Data()
1255  # *** empty log message ***       elif name == "C" :
1256  #           return escript.Data()
1257  # Revision 1.1.2.1  2004/12/12 22:53:47  gross       elif name == "D" :
1258  # linearPDE has been renamed LinearPDE           return escript.Data()
1259  #       elif name == "X" :
1260  # Revision 1.1.1.1.2.7  2004/12/07 10:13:08  gross           return escript.Data()
1261  # GMRES added       elif name == "Y" :
1262  #           return self.getCoefficient("f")
1263  # Revision 1.1.1.1.2.6  2004/12/07 03:19:50  gross       elif name == "d" :
1264  # options for GMRES and PRES20 added           return escript.Data()
1265  #       elif name == "y" :
1266  # Revision 1.1.1.1.2.5  2004/12/01 06:25:15  gross           return escript.Data()
1267  # some small changes       elif name == "d_contact" :
1268  #           return escript.Data()
1269  # Revision 1.1.1.1.2.4  2004/11/24 01:50:21  gross       elif name == "y_contact" :
1270  # Finley solves 4M unknowns now           return escript.Data()
1271  #       elif name == "r" :
1272  # Revision 1.1.1.1.2.3  2004/11/15 06:05:26  gross           return escript.Data()
1273  # poisson solver added       elif name == "q" :
1274  #           return self.getCoefficient("q")
1275  # Revision 1.1.1.1.2.2  2004/11/12 06:58:15  gross       else:
1276  # a lot of changes to get the linearPDE class running: most important change is that there is no matrix format exposed to the user anymore. the format is chosen by the Domain according to the solver and symmetry           raise SystemError,"unknown PDE coefficient %s",name
 #  
 # Revision 1.1.1.1.2.1  2004/10/28 22:59:22  gross  
 # finley's RecTest.py is running now: problem in SystemMatrixAdapater fixed  
 #  
 # Revision 1.1.1.1  2004/10/26 06:53:56  jgs  
 # initial import of project esys2  
 #  
 # Revision 1.3.2.3  2004/10/26 06:43:48  jgs  
 # committing Lutz's and Paul's changes to brach jgs  
 #  
 # Revision 1.3.4.1  2004/10/20 05:32:51  cochrane  
 # Added incomplete Doxygen comments to files, or merely put the docstrings that already exist into Doxygen form.  
 #  
 # Revision 1.3  2004/09/23 00:53:23  jgs  
 # minor fixes  
 #  
 # Revision 1.1  2004/08/28 12:58:06  gross  
 # SimpleSolve is not running yet: problem with == of functionsspace  
 #  
 #  

Legend:
Removed from v.104  
changed lines
  Added in v.114

  ViewVC Help
Powered by ViewVC 1.1.26