/[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 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 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 = 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._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[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 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    
983     def getDim(self):  class AdvectivePDE(LinearPDE):
984       """     """
985       @brief returns the spatial dimension of the PDE     @brief Class to handel a linear PDE domineated by advective terms:
986       """    
987       return self.getDomain().getDim()     class to define a linear PDE of the form
988    
989     def getNumEquations(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 equations  
      """  
      if self.__numEquations>0:  
          return self.__numEquations  
      else:  
          raise ValueError,"Number of equations is undefined. Please specify argument numEquations."  
990    
991     def getNumSolutions(self):       with boundary conditons:
      """  
      @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."  
992    
993            n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_{ik}u_k = - n_j*X_{ij} + y_i
994    
995     def checkSymmetry(self):      and contact conditions
       """  
       @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"  
996    
997        return None          n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_contact_{ik}[u_k] = - n_j*X_{ij} + y_contact_i
998    
999     def getFlux(self,u):      and constraints:
        """  
        @brief returns the flux J_ij for a given u  
1000    
1001              J_ij=A_{ijkl}u_{k,l}+B_{ijk}u_k-X_{ij}           u_i=r_i where q_i>0
1002    
1003         @param u argument of the operator      The PDE is solved by stabilizing the advective terms using SUPG approach:
1004    
1005         """         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)
        raise SystemError,"getFlux is not implemented yet"  
        return None  
1006    
1007     def applyOperator(self,u):      where  
        """  
        @brief applies the operator of the PDE to a given solution u in weak from  
1008    
1009         @param u argument of the operator             b_{ik}=length(B_{i:k})*h/2/length(A_{i:k:})
1010               c_{ik}=length(C_{i:k})*h/2/length(A_{i:k:})
1011    
1012         """                        alpha/3        alpha<3
1013         return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())             xi(alpha)=          for                  approximating cotanh(alpha)-1/alpha
1014                                                                                                                                                                                     1             alpha>=3
1015     def getResidual(self,u):     """
1016         """     def __getXi(self,alpha):
1017         @brief return the residual of u in the weak from           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 939  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._setValue(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._setValue(Y=f,q=q)         self._setValue(f=f,q=q)
1167    
1168                                                                                                                                                                 def getCoefficientOfPDE(self,name):
1169  # $Log$       """
1170  # Revision 1.3  2004/12/17 07:43:10  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.3  2004/12/16 00:12:34  gross       if name == "A" :
1174  # __init__ of LinearPDE does not accept any coefficients anymore           return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))
1175  #       elif name == "B" :
1176  # Revision 1.1.2.2  2004/12/14 03:55:01  jgs           return escript.Data()
1177  # *** empty log message ***       elif name == "C" :
1178  #           return escript.Data()
1179  # Revision 1.1.2.1  2004/12/12 22:53:47  gross       elif name == "D" :
1180  # linearPDE has been renamed LinearPDE           return escript.Data()
1181  #       elif name == "X" :
1182  # Revision 1.1.1.1.2.7  2004/12/07 10:13:08  gross           return escript.Data()
1183  # GMRES added       elif name == "Y" :
1184  #           return self.getCoefficient("f")
1185  # Revision 1.1.1.1.2.6  2004/12/07 03:19:50  gross       elif name == "d" :
1186  # options for GMRES and PRES20 added           return escript.Data()
1187  #       elif name == "y" :
1188  # Revision 1.1.1.1.2.5  2004/12/01 06:25:15  gross           return escript.Data()
1189  # some small changes       elif name == "d_contact" :
1190  #           return escript.Data()
1191  # Revision 1.1.1.1.2.4  2004/11/24 01:50:21  gross       elif name == "y_contact" :
1192  # Finley solves 4M unknowns now           return escript.Data()
1193  #       elif name == "r" :
1194  # Revision 1.1.1.1.2.3  2004/11/15 06:05:26  gross           return escript.Data()
1195  # poisson solver added       elif name == "q" :
1196  #           return self.getCoefficient("q")
1197  # Revision 1.1.1.1.2.2  2004/11/12 06:58:15  gross       else:
1198  # 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.108
  ViewVC Help
Powered by ViewVC 1.1.26