/[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 102 by jgs, Wed Dec 15 07:08:39 2004 UTC revision 108 by jgs, Thu Jan 27 06:21:59 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 215  class LinearPDE: Line 173  class LinearPDE:
173     GMRES=util.GMRES     GMRES=util.GMRES
174     PRES20=util.PRES20     PRES20=util.PRES20
175    
176     def __init__(self,**args):     def __init__(self,domain,numEquations=0,numSolutions=0):
177       """       """
178       @brief initializes a new linear PDE.       @brief initializes a new linear PDE.
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
199       self.__domain=None       self.__domain=domain
200       self.__numEquations=0       self.__numEquations=numEquations
201       self.__numSolutions=0       self.__numSolutions=numSolutions
202       self.cleanCoefficients()       self.cleanCoefficients()
203    
204       self.__operator=escript.Operator()       self.__operator=escript.Operator()
# Line 236  class LinearPDE: Line 208  class LinearPDE:
208       self.__solution=escript.Data()       self.__solution=escript.Data()
209       self.__solution_isValid=False       self.__solution_isValid=False
210    
      # check the arguments  
      coef={}  
      for arg in args.iterkeys():  
           if arg=="domain":  
               self.__domain=args[arg]  
           elif arg=="numEquations":  
               self.__numEquations=args[arg]  
           elif arg=="numSolutions":  
               self.__numSolutions=args[arg]  
           elif _PDECoefficientTypes.has_key(arg):  
               coef[arg]=args[arg]  
           else:  
               raise ValueError,"Illegal argument %s"%arg  
   
      # get the domain of the PDE  
      self.__domain=identifyDomain(self.__domain,coef)  
   
211       # set some default values:       # set some default values:
212       self.__homogeneous_constraint=True       self.__homogeneous_constraint=True
213       self.__row_function_space=escript.Solution(self.__domain)       self.__row_function_space=escript.Solution(self.__domain)
# Line 262  class LinearPDE: Line 217  class LinearPDE:
217       self.__matrix_type=self.__domain.getSystemMatrixTypeId(util.DEFAULT_METHOD,False)       self.__matrix_type=self.__domain.getSystemMatrixTypeId(util.DEFAULT_METHOD,False)
218       self.__sym=False       self.__sym=False
219       self.__lumping=False       self.__lumping=False
220       self.__numEquations=0  
221       self.__numSolutions=0     def createCoefficient(self, name):
222       # now we can set the ceofficients:       """
223       self._setCoefficient(**coef)       @brief create a data object corresponding to coefficient name
224         @param name
225         """
226         return escript.Data(shape = getShapeOfCoefficient(name), \
227                             what = getFunctionSpaceOfCoefficient(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._setCoefficient(**coefficients)        return self.COEFFICIENTS.has_key(name)
         
255    
256     def _setCoefficient(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:  
                 if self.__numEquations>0 and  self.__numSolutions>0:  
                    alteredCoefficients[i]=escript.Data(d,self.getShapeOfCoefficient(i),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 352  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 374  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 441  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 628  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 662  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 678  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[i,j,k])>tol:
660                                  if verbose: print "non-symmetric PDE because B[%d,%d,%d]!=C[%d,%d,%d]"%(i,j,k,i,j,k)
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 685  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 704  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 726  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 755  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 856  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 876  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    
983     def getNumEquations(self):  class AdvectivePDE(LinearPDE):
984       """     """
985       @brief returns the number of equations     @brief Class to handel a linear PDE domineated by advective terms:
986       """    
987       if self.__numEquations>0:     class to define a linear PDE of the form
          return self.__numEquations  
      else:  
          raise ValueError,"Number of equations is undefined. Please specify argument numEquations."  
988    
989     def getNumSolutions(self):       -(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
      """  
      @brief returns the number of unknowns  
      """  
      if self.__numSolutions>0:  
         return self.__numSolutions  
      else:  
         raise ValueError,"Number of solution is undefined. Please specify argument numSolutions."  
990    
991         with boundary conditons:
992    
993     def checkSymmetry(self):          n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_{ik}u_k = - n_j*X_{ij} + y_i
       """  
       @brief returns if the Operator is symmetric. This is a very expensive operation!!! The symmetry flag is not altered.  
       """  
       raise SystemError,"checkSymmetry is not implemented yet"  
994    
995        return None      and contact conditions
996    
997     def getFlux(self,u):          n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_contact_{ik}[u_k] = - n_j*X_{ij} + y_contact_i
        """  
        @brief returns the flux J_ij for a given u  
998    
999              J_ij=A_{ijkl}u_{k,l}+B_{ijk}u_k-X_{ij}      and constraints:
1000    
1001         @param u argument of the operator           u_i=r_i where q_i>0
1002    
1003         """      The PDE is solved by stabilizing the advective terms using SUPG approach:
        raise SystemError,"getFlux is not implemented yet"  
        return None  
1004    
1005     def applyOperator(self,u):         A_{ijkl}<-A_{ijkl}+0.5*h*(xi(b_{ik})*B_{ijk}*B_{ilk}/length(B_{i:k})^2)+0.5*h*xi_{c_{ik}}*(C_{ikj}*C_{ikl}/length(C_{ik:})^2)
        """  
        @brief applies the operator of the PDE to a given solution u in weak from  
1006    
1007         @param u argument of the operator      where  
1008    
1009         """             b_{ik}=length(B_{i:k})*h/2/length(A_{i:k:})
1010         return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())             c_{ik}=length(C_{i:k})*h/2/length(A_{i:k:})
1011                                                                                                                                                              
1012     def getResidual(self,u):                        alpha/3        alpha<3
1013         """             xi(alpha)=          for                  approximating cotanh(alpha)-1/alpha
1014         @brief return the residual of u in the weak from                         1             alpha>=3
1015       """
1016       def __getXi(self,alpha):
1017             c=alpha-3.
1018             return c*c.whereNegative()/3.+1.
1019    
1020       def __getUpdateVector(self,V,hover2,alphaByU):
1021         v=util.length(V)
1022         v_max=util.Lsup(v)
1023         if v_max>0:
1024             V/=v+v_max*self.TOL
1025             alpha=alphaByU*v
1026             A_bar=v*hover2*self.__getXi(alpha)
1027             print "-------------"
1028             print "@ max alpha ",util.Lsup(alpha)
1029             print "-------------"
1030         else:
1031             A_bar=1.
1032         return V,A_bar
1033    
1034       def __getAlphaByU(self,A,hover2):
1035          a=util.length(A)
1036          a_max=util.Lsup(a)
1037          if a_max>0:
1038             return hover2/(a+a_max*self.TOL)
1039          else:
1040             return 1./self.TOL
1041    
1042    
1043       def getCoefficientOfPDE(self,name):
1044         """
1045         @brief return the value of the coefficient name of the general PDE
1046         @param name
1047         """
1048         if name == "A" :
1049             A=self.getCoefficient("A")
1050             B=self.getCoefficient("B")
1051             C=self.getCoefficient("C")
1052             if not B.isEmpty() or not C.isEmpty():
1053                 if A.isEmpty():
1054                     A=self.createNewCoefficient("A")
1055                 else:
1056                     A=A[:]
1057                 hover2=self.getDomain().getSize()/2.
1058                 if self.getNumEquations()>1:
1059                    if self.getNumSolutions()>1:
1060                       for i in range(self.getNumEquations()):
1061                          for k in range(self.getNumSolutions()):
1062                             alphaByU=self.__getAlphaByU(A[i,:,k,:],hover2)
1063                             if not B.isEmpty():
1064                                 b_sub,f=self.__getUpdateVector(B[i,:,k],hover2,alphaByU)
1065                                 for j in range(self.getDim()):
1066                                    for l in range(self.getDim()):
1067                                       A[i,j,k,l]+=f*b_sub[j]*b_sub[l]
1068                             if not C.isEmpty():
1069                                 c_sub,f=self.__getUpdateVector(C[i,k,:],hover2,alphaByU)
1070                                 for j in range(self.getDim()):
1071                                    for l in range(self.getDim()):
1072                                       A[i,j,k,l]+=f*c_sub[j]*c_sub[l]
1073                    else:  
1074                       for i in range(self.getNumEquations()):
1075                          alphaByU=self.__getAlphaByU(A[i,:,:],hover2)
1076                          if not B.isEmpty():
1077                              b_sub,f=self.__getUpdateVector(B[i,:],hover2,alphaByU)
1078                              for j in range(self.getDim()):
1079                                 for l in range(self.getDim()):
1080                                     A[i,j,l]+=f*b_sub[j]*b_sub[l]
1081                          if not C.isEmpty():
1082                               c_sub,f=self.__getUpdateVector(C[i,:],hover2,alphaByU)
1083                               for j in range(self.getDim()):
1084                                  for l in range(self.getDim()):
1085                                     A[i,j,l]+=f*c_sub[j]*c_sub[l]
1086                 else:
1087                    if self.getNumSolutions()>1:
1088                       for k in range(self.getNumSolutions()):
1089                          alphaByU=self.__getAlphaByU(A[:,k,:],hover2)
1090                          if not B.isEmpty():
1091                             b_sub,f=self.__getUpdateVector(B[:,k],hover2,alphaByU)
1092                             for j in range(self.getDim()):
1093                                for l in range(self.getDim()):
1094                                       A[j,k,l]+=f*b_sub[j]*b_sub[l]
1095                          if not C.isEmpty():
1096                             c_sub,f=self.__getUpdateVector(C[k,:],hover2,alphaByU)
1097                             for j in range(self.getDim()):
1098                                for l in range(self.getDim()):
1099                                   A[j,k,l]+=f*c_sub[j]*c_sub[l]
1100                    else:  
1101                       alphaByU=self.__getAlphaByU(A[:,:],hover2)
1102                       if not B.isEmpty():
1103                           b_sub,f=self.__getUpdateVector(B[:],hover2,alphaByU)
1104                           for j in range(self.getDim()):
1105                              for l in range(self.getDim()):
1106                                 A[j,l]+=f*b_sub[j]*b_sub[l]
1107                       if not C.isEmpty():
1108                          c_sub,f=self.__getUpdateVector(C[:],hover2,alphaByU)
1109                          for j in range(self.getDim()):
1110                              for l in range(self.getDim()):
1111                                 A[j,l]+=f*c_sub[j]*c_sub[l]
1112             return A
1113         elif name == "B" :
1114             return self.getCoefficient("B")
1115         elif name == "C" :
1116             return self.getCoefficient("C")
1117         elif name == "D" :
1118             return self.getCoefficient("D")
1119         elif name == "X" :
1120             return self.getCoefficient("X")
1121         elif name == "Y" :
1122             return self.getCoefficient("Y")
1123         elif name == "d" :
1124             return self.getCoefficient("d")
1125         elif name == "y" :
1126             return self.getCoefficient("y")
1127         elif name == "d_contact" :
1128             return self.getCoefficient("d_contact")
1129         elif name == "y_contact" :
1130             return self.getCoefficient("y_contact")
1131         elif name == "r" :
1132             return self.getCoefficient("r")
1133         elif name == "q" :
1134             return self.getCoefficient("q")
1135         else:
1136             raise SystemError,"unknown PDE coefficient %s",name
1137    
        @param u  
        """  
        return self.applyOperator(u)-self.getRightHandSide()  
1138    
1139  class Poisson(LinearPDE):  class Poisson(LinearPDE):
1140     """     """
# Line 957  class Poisson(LinearPDE): Line 1154  class Poisson(LinearPDE):
1154                                                                                                                                                                                                                                                                                                                            
1155     """     """
1156    
1157     def __init__(self,domain=None,f=escript.Data(),q=escript.Data()):     def __init__(self,domain,f=escript.Data(),q=escript.Data()):
1158         LinearPDE.__init__(self,domain=identifyDomain(domain,{"f":f, "q":q}))         LinearPDE.__init__(self,domain,1,1)
1159         self._setCoefficient(A=numarray.identity(self.getDomain().getDim()))         self.COEFFICIENTS={
1160           "f"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1161           "q"         : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.BOTH)}
1162         self.setSymmetryOn()         self.setSymmetryOn()
1163         self.setValue(f,q)         self.setValue(f,q)
1164    
1165     def setValue(self,f=escript.Data(),q=escript.Data()):     def setValue(self,f=escript.Data(),q=escript.Data()):
1166         self._setCoefficient(Y=f,q=q)         self._setValue(f=f,q=q)
1167    
1168                                                                                                                                                                 def getCoefficientOfPDE(self,name):
1169  # $Log$       """
1170  # Revision 1.2  2004/12/15 07:08:27  jgs       @brief return the value of the coefficient name of the general PDE
1171  # *** empty log message ***       @param name
1172  #       """
1173  # Revision 1.1.2.2  2004/12/14 03:55:01  jgs       if name == "A" :
1174  # *** empty log message ***           return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))
1175  #       elif name == "B" :
1176  # Revision 1.1.2.1  2004/12/12 22:53:47  gross           return escript.Data()
1177  # linearPDE has been renamed LinearPDE       elif name == "C" :
1178  #           return escript.Data()
1179  # Revision 1.1.1.1.2.7  2004/12/07 10:13:08  gross       elif name == "D" :
1180  # GMRES added           return escript.Data()
1181  #       elif name == "X" :
1182  # Revision 1.1.1.1.2.6  2004/12/07 03:19:50  gross           return escript.Data()
1183  # options for GMRES and PRES20 added       elif name == "Y" :
1184  #           return self.getCoefficient("f")
1185  # Revision 1.1.1.1.2.5  2004/12/01 06:25:15  gross       elif name == "d" :
1186  # some small changes           return escript.Data()
1187  #       elif name == "y" :
1188  # Revision 1.1.1.1.2.4  2004/11/24 01:50:21  gross           return escript.Data()
1189  # Finley solves 4M unknowns now       elif name == "d_contact" :
1190  #           return escript.Data()
1191  # Revision 1.1.1.1.2.3  2004/11/15 06:05:26  gross       elif name == "y_contact" :
1192  # poisson solver added           return escript.Data()
1193  #       elif name == "r" :
1194  # Revision 1.1.1.1.2.2  2004/11/12 06:58:15  gross           return escript.Data()
1195  # 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       elif name == "q" :
1196  #           return self.getCoefficient("q")
1197  # Revision 1.1.1.1.2.1  2004/10/28 22:59:22  gross       else:
1198  # finley's RecTest.py is running now: problem in SystemMatrixAdapater fixed           raise SystemError,"unknown PDE coefficient %s",name
 #  
 # 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.102  
changed lines
  Added in v.108

  ViewVC Help
Powered by ViewVC 1.1.26