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

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