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

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

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

revision 104 by jgs, Fri Dec 17 07:43:12 2004 UTC revision 142 by jgs, Mon Jul 25 05:28:20 2005 UTC
# Line 3  Line 3 
3  ## @file linearPDEs.py  ## @file linearPDEs.py
4    
5  """  """
6  @brief Functions and classes for linear PDEs  Functions and classes for linear PDEs
7  """  """
8    
9  import escript  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     """     """
16     @brief     Compare two tuples
17    
18     @param t1     \param t1 The first tuple
19     @param t2     \param t2 The second tuple
20     """     """
21    
22     dif=t1[0]+t1[1]-(t2[0]+t2[1])     dif=t1[0]+t1[1]-(t2[0]+t2[1])
23     if dif<0: return 1     if dif<0: return 1
24     elif dif>0: return -1     elif dif>0: return -1
25     else: return 0     else: return 0
26    
27  class PDECoefficientType:  def ELMAN_RAMAGE(P):
28        return (P-1.).wherePositive()*0.5*(1.-1./(P+1.e-15))
29    
30    def SIMPLIFIED_BROOK_HUGHES(P):
31        c=(P-3.).whereNegative()
32        return P/6.*c+1./2.*(1.-c)
33    
34    def HALF(P):
35        return escript.Scalar(0.5,P.getFunctionSpace())
36    
37    class PDECoefficient:
38      """      """
39      @brief      A class for PDE coefficients
40      """      """
41      # identifier for location of Data objects defining coefficients      # identifier for location of Data objects defining COEFFICIENTS
42      INTERIOR=0      INTERIOR=0
43      BOUNDARY=1      BOUNDARY=1
44      CONTACT=2      CONTACT=2
45      CONTINUOUS=3      CONTINUOUS=3
46      # identifier in the pattern of coefficients:      # identifier in the pattern of COEFFICIENTS:
47      # 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
48      # number of unknowns.      # number of unknowns.
49      EQUATION=3      EQUATION=3
# Line 91  class PDECoefficientType: Line 55  class PDECoefficientType:
55      BOTH=7      BOTH=7
56      def __init__(self,where,pattern,altering):      def __init__(self,where,pattern,altering):
57         """         """
58         @brief Initialise a PDE Coefficient type         Initialise a PDE Coefficient type
59         """         """
60         self.what=where         self.what=where
61         self.pattern=pattern         self.pattern=pattern
62         self.altering=altering         self.altering=altering
63           self.resetValue()
64    
65        def resetValue(self):
66           """
67           resets coefficient value to default
68           """
69           self.value=escript.Data()
70    
71      def getFunctionSpace(self,domain):      def getFunctionSpace(self,domain):
72         """         """
73         @brief defines the FunctionSpace of the coefficient on the domain         defines the FunctionSpace of the coefficient on the domain
74    
75         @param domain         @param domain:
76         """         """
77         if self.what==self.INTERIOR: return escript.Function(domain)         if self.what==self.INTERIOR: return escript.Function(domain)
78         elif self.what==self.BOUNDARY: return escript.FunctionOnBoundary(domain)         elif self.what==self.BOUNDARY: return escript.FunctionOnBoundary(domain)
79         elif self.what==self.CONTACT: return escript.FunctionOnContactZero(domain)         elif self.what==self.CONTACT: return escript.FunctionOnContactZero(domain)
80         elif self.what==self.CONTINUOUS: return escript.ContinuousFunction(domain)         elif self.what==self.CONTINUOUS: return escript.ContinuousFunction(domain)
81    
82        def getValue(self):
83           """
84           returns the value of the coefficient:
85           """
86           return self.value
87        
88        def setValue(self,newValue):
89           """
90           set the value of the coefficient to new value
91           """
92           self.value=newValue
93        
94      def isAlteringOperator(self):      def isAlteringOperator(self):
95          """          """
96      @brief return true if the operator of the PDE is changed when the coefficient is changed      return true if the operator of the PDE is changed when the coefficient is changed
97      """      """
98          if self.altering==self.OPERATOR or self.altering==self.BOTH:          if self.altering==self.OPERATOR or self.altering==self.BOTH:
99              return not None              return not None
# Line 119  class PDECoefficientType: Line 102  class PDECoefficientType:
102    
103      def isAlteringRightHandSide(self):      def isAlteringRightHandSide(self):
104          """          """
105      @brief return true if the right hand side of the PDE is changed when the coefficient is changed      return true if the right hand side of the PDE is changed when the coefficient is changed
106      """      """
107          if self.altering==self.RIGHTHANDSIDE or self.altering==self.BOTH:          if self.altering==self.RIGHTHANDSIDE or self.altering==self.BOTH:
108              return not None              return not None
# Line 128  class PDECoefficientType: Line 111  class PDECoefficientType:
111    
112      def estimateNumEquationsAndNumSolutions(self,shape=(),dim=3):      def estimateNumEquationsAndNumSolutions(self,shape=(),dim=3):
113         """         """
114         @brief tries to estimate the number of equations in a given tensor shape for a given spatial dimension dim         tries to estimate the number of equations in a given tensor shape for a given spatial dimension dim
115    
116         @param shape         @param shape:
117         @param dim         @param dim:
118         """         """
119         if len(shape)>0:         if len(shape)>0:
120             num=max(shape)+1             num=max(shape)+1
# Line 152  class PDECoefficientType: Line 135  class PDECoefficientType:
135    
136      def buildShape(self,e=1,u=1,dim=3):      def buildShape(self,e=1,u=1,dim=3):
137          """          """
138      @brief builds the required shape for a given number of equations e, number of unknowns u and spatial dimension dim      builds the required shape for a given number of equations e, number of unknowns u and spatial dimension dim
139    
140      @param e      @param e:
141      @param u      @param u:
142      @param dim      @param dim:
143      """      """
144          s=()          s=()
145          for i in self.pattern:          for i in self.pattern:
# Line 168  class PDECoefficientType: Line 151  class PDECoefficientType:
151                  s=s+(dim,)                  s=s+(dim,)
152          return s          return s
153    
 _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),  
 }  
   
154  class LinearPDE:  class LinearPDE:
155     """     """
156     @brief Class to define a linear PDE     Class to handle a linear PDE
157        
158     class to define a linear PDE of the form     class to define a linear PDE of the form
159    
160       \f[
161       -(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       -(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
162       \f]
163    
164       with boundary conditons:     with boundary conditons:
165    
166          n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_{ik}u_k = - n_j*X_{ij} + y_i     \f[
167       n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_{ik}u_k = - n_j*X_{ij} + y_i
168       \f]
169    
170      and contact conditions     and contact conditions
171    
172          n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_contact_{ik}[u_k] = - n_j*X_{ij} + y_contact_i     \f[
173       n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_contact_{ik}[u_k] = - n_j*X_{ij} + y_contact_i
174       \f]
175    
176      and constraints:     and constraints:
177    
178           u_i=r_i where q_i>0     \f[
179       u_i=r_i \quad \mathrm{where} \quad q_i>0
180       \f]
181    
182     """     """
183       TOL=1.e-13
184     DEFAULT_METHOD=util.DEFAULT_METHOD     DEFAULT_METHOD=util.DEFAULT_METHOD
185     DIRECT=util.DIRECT     DIRECT=util.DIRECT
186     CHOLEVSKY=util.CHOLEVSKY     CHOLEVSKY=util.CHOLEVSKY
# Line 214  class LinearPDE: Line 191  class LinearPDE:
191     SSOR=util.SSOR     SSOR=util.SSOR
192     GMRES=util.GMRES     GMRES=util.GMRES
193     PRES20=util.PRES20     PRES20=util.PRES20
194       ILU0=util.ILU0
195       JACOBI=util.JACOBI
196    
197     def __init__(self,domain,numEquations=0,numSolutions=0):     def __init__(self,domain,numEquations=0,numSolutions=0):
198       """       """
199       @brief initializes a new linear PDE.       initializes a new linear PDE.
200    
201       @param args       @param args:
202       """       """
203         # COEFFICIENTS can be overwritten by subclasses:
204         self.COEFFICIENTS={
205           "A"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM,PDECoefficient.SOLUTION,PDECoefficient.DIM),PDECoefficient.OPERATOR),
206           "B"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),
207           "C"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION,PDECoefficient.DIM),PDECoefficient.OPERATOR),
208           "D"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),
209           "X"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM),PDECoefficient.RIGHTHANDSIDE),
210           "Y"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
211           "d"         : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),
212           "y"         : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
213           "d_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),
214           "y_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
215           "r"         : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
216           "q"         : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.SOLUTION,),PDECoefficient.BOTH)}
217    
218       # initialize attributes       # initialize attributes
219       self.__debug=None       self.__debug=None
# Line 246  class LinearPDE: Line 239  class LinearPDE:
239       self.__sym=False       self.__sym=False
240       self.__lumping=False       self.__lumping=False
241    
242       def createCoefficient(self, name):
243         """
244         create a data object corresponding to coefficient name
245         @param name:
246         """
247         return escript.Data(shape = getShapeOfCoefficient(name), \
248                             what = getFunctionSpaceForCoefficient(name))
249    
250       def __del__(self):
251         pass
252    
253     def getCoefficient(self,name):     def getCoefficient(self,name):
254       """       """
255       @brief return the value of the coefficient name       return the value of the parameter name
256    
257       @param name       @param name:
258       """       """
259       return self.__coefficient[name]       return self.COEFFICIENTS[name].getValue()
260    
261     def setValue(self,**coefficients):     def getCoefficientOfPDE(self,name):
262         """
263         return the value of the coefficient name of the general PDE.
264         This method is called by the assembling routine it can be
265         overwritten to map coefficients of a particualr PDE to the general PDE.
266    
267         @param name:
268         """
269         return self.getCoefficient(name)
270    
271       def hasCoefficient(self,name):
272        """        """
273        @brief sets new values to coefficients        return true if name is the name of a coefficient
274    
275        @param coefficients        @param name:
276        """        """
277        self._setValue(**coefficients)        return self.COEFFICIENTS.has_key(name)
         
278    
279     def _setValue(self,**coefficients):     def getFunctionSpaceForEquation(self):
280         """
281         return true if the test functions should use reduced order
282         """
283         return self.__row_function_space
284    
285       def getFunctionSpaceForSolution(self):
286         """
287         return true if the interpolation of the solution should use reduced order
288         """
289         return self.__column_function_space
290    
291       def setValue(self,**coefficients):
292        """        """
293        @brief sets new values to coefficients        sets new values to coefficients
294    
295        @param coefficients        @param coefficients:
296        """        """
297          self.__setValue(**coefficients)
298                
       # 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()  
299    
300     def cleanCoefficients(self):     def cleanCoefficients(self):
301       """       """
302       @brief resets all coefficients to default values.       resets all coefficients to default values.
303         """
304         for i in self.COEFFICIENTS.iterkeys():
305             self.COEFFICIENTS[i].resetValue()
306    
307       def createNewCoefficient(self,name):
308       """       """
309       self.__coefficient={}       returns a new coefficient appropriate for coefficient name:
310       for i in _PDECoefficientTypes.iterkeys():       """
311           self.__coefficient[i]=escript.Data()       return escript.Data(0,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))
312          
313    
314     def getShapeOfCoefficient(self,name):     def getShapeOfCoefficient(self,name):
315       """       """
316       @brief return the shape of the coefficient name       return the shape of the coefficient name
317    
318       @param name       @param name:
319       """       """
320       if self.hasCoefficient(name):       if self.hasCoefficient(name):
321          return _PDECoefficientTypes[name].buildShape(self.getNumEquations(),self.getNumSolutions(),self.getDomain().getDim())          return self.COEFFICIENTS[name].buildShape(self.getNumEquations(),self.getNumSolutions(),self.getDomain().getDim())
322       else:       else:
323          raise ValueError,"Solution coefficient %s requested"%name          raise ValueError,"Solution coefficient %s requested"%name
324    
325     def getFunctionSpaceOfCoefficient(self,name):     def getFunctionSpaceForCoefficient(self,name):
326       """       """
327       @brief return the atoms of the coefficient name       return the atoms of the coefficient name
328    
329       @param name       @param name:
330       """       """
331       if self.hasCoefficient(name):       if self.hasCoefficient(name):
332          return _PDECoefficientTypes[name].getFunctionSpace(self.getDomain())          return self.COEFFICIENTS[name].getFunctionSpace(self.getDomain())
333       else:       else:
334          raise ValueError,"Solution coefficient %s requested"%name          raise ValueError,"Solution coefficient %s requested"%name
335    
336     def alteredCoefficient(self,name):     def alteredCoefficient(self,name):
337       """       """
338       @brief annonced that coefficient name has been changed       announce that coefficient name has been changed
339    
340       @param name       @param name:
341       """       """
342       if self.hasCoefficient(name):       if self.hasCoefficient(name):
343          if _PDECoefficientTypes[name].isAlteringOperator(): self.__rebuildOperator()          if self.COEFFICIENTS[name].isAlteringOperator(): self.__rebuildOperator()
344          if _PDECoefficientTypes[name].isAlteringRightHandSide(): self.__rebuildRightHandSide()          if self.COEFFICIENTS[name].isAlteringRightHandSide(): self.__rebuildRightHandSide()
345       else:       else:
346          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  
347    
348     # ===== debug ==============================================================     # ===== debug ==============================================================
349     def setDebugOn(self):     def setDebugOn(self):
350         """         """
        @brief  
351         """         """
352         self.__debug=not None         self.__debug=not None
353    
354     def setDebugOff(self):     def setDebugOff(self):
355         """         """
        @brief  
356         """         """
357         self.__debug=None         self.__debug=None
358    
359     def debug(self):     def debug(self):
360         """         """
361         @brief returns true if the PDE is in the debug mode         returns true if the PDE is in the debug mode
362         """         """
363         return self.__debug         return self.__debug
364    
365     #===== Lumping ===========================     #===== Lumping ===========================
366     def setLumpingOn(self):     def setLumpingOn(self):
367        """        """
368        @brief indicates to use matrix lumping        indicates to use matrix lumping
369        """        """
370        if not self.isUsingLumping():        if not self.isUsingLumping():
          raise SystemError,"Lumping is not working yet! Talk to the experts"  
371           if self.debug() : print "PDE Debug: lumping is set on"           if self.debug() : print "PDE Debug: lumping is set on"
372           self.__rebuildOperator()           self.__rebuildOperator()
373           self.__lumping=True           self.__lumping=True
374    
375     def setLumpingOff(self):     def setLumpingOff(self):
376        """        """
377        @brief switches off matrix lumping        switches off matrix lumping
378        """        """
379        if self.isUsingLumping():        if self.isUsingLumping():
380           if self.debug() : print "PDE Debug: lumping is set off"           if self.debug() : print "PDE Debug: lumping is set off"
# Line 433  class LinearPDE: Line 383  class LinearPDE:
383    
384     def setLumping(self,flag=False):     def setLumping(self,flag=False):
385        """        """
386        @brief set the matrix lumping flag to flag        set the matrix lumping flag to flag
387        """        """
388        if flag:        if flag:
389           self.setLumpingOn()           self.setLumpingOn()
# Line 442  class LinearPDE: Line 392  class LinearPDE:
392    
393     def isUsingLumping(self):     def isUsingLumping(self):
394        """        """
395        @brief        
396        """        """
397        return self.__lumping        return self.__lumping
398    
399     #============ method business =========================================================     #============ method business =========================================================
400     def setSolverMethod(self,solver=util.DEFAULT_METHOD):     def setSolverMethod(self,solver=util.DEFAULT_METHOD):
401         """         """
402         @brief sets a new solver         sets a new solver
403         """         """
404         if not solver==self.getSolverMethod():         if not solver==self.getSolverMethod():
405             self.__solver_method=solver             self.__solver_method=solver
# Line 458  class LinearPDE: Line 408  class LinearPDE:
408    
409     def getSolverMethod(self):     def getSolverMethod(self):
410         """         """
411         @brief returns the solver method         returns the solver method
412         """         """
413         return self.__solver_method         return self.__solver_method
414    
415     #============ tolerance business =========================================================     #============ tolerance business =========================================================
416     def setTolerance(self,tol=1.e-8):     def setTolerance(self,tol=1.e-8):
417         """         """
418         @brief resets the tolerance to tol.         resets the tolerance to tol.
419         """         """
420         if not tol>0:         if not tol>0:
421             raise ValueException,"Tolerance as to be positive"             raise ValueException,"Tolerance as to be positive"
# Line 475  class LinearPDE: Line 425  class LinearPDE:
425         return         return
426     def getTolerance(self):     def getTolerance(self):
427         """         """
428         @brief returns the tolerance set for the solution         returns the tolerance set for the solution
429         """         """
430         return self.__tolerance         return self.__tolerance
431    
432     #===== symmetry  flag ==========================     #===== symmetry  flag ==========================
433     def isSymmetric(self):     def isSymmetric(self):
434        """        """
435        @brief returns true is the operator is considered to be symmetric        returns true is the operator is considered to be symmetric
436        """        """
437        return self.__sym        return self.__sym
438    
439     def setSymmetryOn(self):     def setSymmetryOn(self):
440        """        """
441        @brief sets the symmetry flag to true        sets the symmetry flag to true
442        """        """
443        if not self.isSymmetric():        if not self.isSymmetric():
444           if self.debug() : print "PDE Debug: Operator is set to be symmetric"           if self.debug() : print "PDE Debug: Operator is set to be symmetric"
# Line 497  class LinearPDE: Line 447  class LinearPDE:
447    
448     def setSymmetryOff(self):     def setSymmetryOff(self):
449        """        """
450        @brief sets the symmetry flag to false        sets the symmetry flag to false
451        """        """
452        if self.isSymmetric():        if self.isSymmetric():
453           if self.debug() : print "PDE Debug: Operator is set to be unsymmetric"           if self.debug() : print "PDE Debug: Operator is set to be unsymmetric"
# Line 506  class LinearPDE: Line 456  class LinearPDE:
456    
457     def setSymmetryTo(self,flag=False):     def setSymmetryTo(self,flag=False):
458       """       """
459       @brief sets the symmetry flag to flag       sets the symmetry flag to flag
460    
461       @param flag       @param flag:
462       """       """
463       if flag:       if flag:
464          self.setSymmetryOn()          self.setSymmetryOn()
# Line 518  class LinearPDE: Line 468  class LinearPDE:
468     #===== order reduction ==========================     #===== order reduction ==========================
469     def setReducedOrderOn(self):     def setReducedOrderOn(self):
470       """       """
471       @brief switches to on reduced order       switches to on reduced order
472       """       """
473       self.setReducedOrderForSolutionOn()       self.setReducedOrderForSolutionOn()
474       self.setReducedOrderForEquationOn()       self.setReducedOrderForEquationOn()
475    
476     def setReducedOrderOff(self):     def setReducedOrderOff(self):
477       """       """
478       @brief switches to full order       switches to full order
479       """       """
480       self.setReducedOrderForSolutionOff()       self.setReducedOrderForSolutionOff()
481       self.setReducedOrderForEquationOff()       self.setReducedOrderForEquationOff()
482    
483     def setReducedOrderTo(self,flag=False):     def setReducedOrderTo(self,flag=False):
484       """       """
485       @brief sets order according to flag       sets order according to flag
486    
487       @param flag       @param flag:
488       """       """
489       self.setReducedOrderForSolutionTo(flag)       self.setReducedOrderForSolutionTo(flag)
490       self.setReducedOrderForEquationTo(flag)       self.setReducedOrderForEquationTo(flag)
# Line 543  class LinearPDE: Line 493  class LinearPDE:
493     #===== order reduction solution ==========================     #===== order reduction solution ==========================
494     def setReducedOrderForSolutionOn(self):     def setReducedOrderForSolutionOn(self):
495       """       """
496       @brief switches to reduced order to interpolate solution       switches to reduced order to interpolate solution
497       """       """
498       new_fs=escript.ReducedSolution(self.getDomain())       new_fs=escript.ReducedSolution(self.getDomain())
499       if self.getFunctionSpaceForSolution()!=new_fs:       if self.getFunctionSpaceForSolution()!=new_fs:
# Line 553  class LinearPDE: Line 503  class LinearPDE:
503    
504     def setReducedOrderForSolutionOff(self):     def setReducedOrderForSolutionOff(self):
505       """       """
506       @brief switches to full order to interpolate solution       switches to full order to interpolate solution
507       """       """
508       new_fs=escript.Solution(self.getDomain())       new_fs=escript.Solution(self.getDomain())
509       if self.getFunctionSpaceForSolution()!=new_fs:       if self.getFunctionSpaceForSolution()!=new_fs:
# Line 563  class LinearPDE: Line 513  class LinearPDE:
513    
514     def setReducedOrderForSolutionTo(self,flag=False):     def setReducedOrderForSolutionTo(self,flag=False):
515       """       """
516       @brief sets order for test functions according to flag       sets order for test functions according to flag
517    
518       @param flag       @param flag:
519       """       """
520       if flag:       if flag:
521          self.setReducedOrderForSolutionOn()          self.setReducedOrderForSolutionOn()
# Line 575  class LinearPDE: Line 525  class LinearPDE:
525     #===== order reduction equation ==========================     #===== order reduction equation ==========================
526     def setReducedOrderForEquationOn(self):     def setReducedOrderForEquationOn(self):
527       """       """
528       @brief switches to reduced order for test functions       switches to reduced order for test functions
529       """       """
530       new_fs=escript.ReducedSolution(self.getDomain())       new_fs=escript.ReducedSolution(self.getDomain())
531       if self.getFunctionSpaceForEquation()!=new_fs:       if self.getFunctionSpaceForEquation()!=new_fs:
# Line 585  class LinearPDE: Line 535  class LinearPDE:
535    
536     def setReducedOrderForEquationOff(self):     def setReducedOrderForEquationOff(self):
537       """       """
538       @brief switches to full order for test functions       switches to full order for test functions
539       """       """
540       new_fs=escript.Solution(self.getDomain())       new_fs=escript.Solution(self.getDomain())
541       if self.getFunctionSpaceForEquation()!=new_fs:       if self.getFunctionSpaceForEquation()!=new_fs:
# Line 595  class LinearPDE: Line 545  class LinearPDE:
545    
546     def setReducedOrderForEquationTo(self,flag=False):     def setReducedOrderForEquationTo(self,flag=False):
547       """       """
548       @brief sets order for test functions according to flag       sets order for test functions according to flag
549    
550       @param flag       @param flag:
551       """       """
552       if flag:       if flag:
553          self.setReducedOrderForEquationOn()          self.setReducedOrderForEquationOn()
554       else:       else:
555          self.setReducedOrderForEquationOff()          self.setReducedOrderForEquationOff()
556                                                                                                                                                                                                                                                                                                                        
   
557     # ==== initialization =====================================================================     # ==== initialization =====================================================================
558     def __makeNewOperator(self):     def __getNewOperator(self):
559         """         """
        @brief  
560         """         """
561         return self.getDomain().newOperator( \         return self.getDomain().newOperator( \
562                             self.getNumEquations(), \                             self.getNumEquations(), \
# Line 617  class LinearPDE: Line 565  class LinearPDE:
565                             self.getFunctionSpaceForSolution(), \                             self.getFunctionSpaceForSolution(), \
566                             self.__matrix_type)                             self.__matrix_type)
567    
568     def __makeNewRightHandSide(self):     def __makeFreshRightHandSide(self):
569         """         """
        @brief  
570         """         """
571         return escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)         if self.debug() : print "PDE Debug: New right hand side allocated"
572           if self.getNumEquations()>1:
573               self.__righthandside=escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)
574           else:
575               self.__righthandside=escript.Data(0.,(),self.getFunctionSpaceForEquation(),True)
576           return self.__righthandside
577    
578     def __makeNewSolution(self):     def __getNewSolution(self):
579         """         """
        @brief  
580         """         """
581         return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)         if self.debug() : print "PDE Debug: New right hand side allocated"
582           if self.getNumSolutions()>1:
583               return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)
584           else:
585               return escript.Data(0.,(),self.getFunctionSpaceForSolution(),True)
586    
587     def __getFreshOperator(self):     def __makeFreshOperator(self):
588         """         """
        @brief  
589         """         """
590         if self.__operator.isEmpty():         if self.__operator.isEmpty():
591             self.__operator=self.__makeNewOperator()             self.__operator=self.__getNewOperator()
592             if self.debug() : print "PDE Debug: New operator allocated"             if self.debug() : print "PDE Debug: New operator allocated"
593         else:         else:
594             self.__operator.setValue(0.)             self.__operator.setValue(0.)
595               self.__operator.resetSolver()
596             if self.debug() : print "PDE Debug: Operator reset to zero"             if self.debug() : print "PDE Debug: Operator reset to zero"
597         return self.__operator         return self.__operator
598    
599     def __getFreshRightHandSide(self):     #============ some serivice functions  =====================================================
600       def getDomain(self):
601         """
602         returns the domain of the PDE
603         """
604         return self.__domain
605    
606       def getDim(self):
607         """
608         returns the spatial dimension of the PDE
609         """
610         return self.getDomain().getDim()
611    
612       def getNumEquations(self):
613         """
614         returns the number of equations
615         """
616         if self.__numEquations>0:
617             return self.__numEquations
618         else:
619             raise ValueError,"Number of equations is undefined. Please specify argument numEquations."
620    
621       def getNumSolutions(self):
622         """
623         returns the number of unknowns
624         """
625         if self.__numSolutions>0:
626            return self.__numSolutions
627         else:
628            raise ValueError,"Number of solution is undefined. Please specify argument numSolutions."
629    
630    
631       def checkSymmetry(self,verbose=True):
632          """
633          returns if the Operator is symmetric. This is a very expensive operation!!! The symmetry flag is not altered.
634          """
635          verbose=verbose or self.debug()
636          out=True
637          if self.getNumSolutions()!=self.getNumEquations():
638             if verbose: print "non-symmetric PDE because of different number of equations and solutions"
639             out=False
640          else:
641             A=self.getCoefficientOfPDE("A")
642             if not A.isEmpty():
643                tol=util.Lsup(A)*self.TOL
644                if self.getNumSolutions()>1:
645                   for i in range(self.getNumEquations()):
646                      for j in range(self.getDim()):
647                         for k in range(self.getNumSolutions()):
648                            for l in range(self.getDim()):
649                                if util.Lsup(A[i,j,k,l]-A[k,l,i,j])>tol:
650                                   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)
651                                   out=False
652                else:
653                   for j in range(self.getDim()):
654                      for l in range(self.getDim()):
655                         if util.Lsup(A[j,l]-A[l,j])>tol:
656                            if verbose: print "non-symmetric PDE because A[%d,%d]!=A[%d,%d]"%(j,l,l,j)
657                            out=False
658             B=self.getCoefficientOfPDE("B")
659             C=self.getCoefficientOfPDE("C")
660             if B.isEmpty() and not C.isEmpty():
661                if verbose: print "non-symmetric PDE because B is not present but C is"
662                out=False
663             elif not B.isEmpty() and C.isEmpty():
664                if verbose: print "non-symmetric PDE because C is not present but B is"
665                out=False
666             elif not B.isEmpty() and not C.isEmpty():
667                tol=(util.Lsup(B)+util.Lsup(C))*self.TOL/2.
668                if self.getNumSolutions()>1:
669                   for i in range(self.getNumEquations()):
670                       for j in range(self.getDim()):
671                          for k in range(self.getNumSolutions()):
672                             if util.Lsup(B[i,j,k]-C[k,i,j])>tol:
673                                  if verbose: print "non-symmetric PDE because B[%d,%d,%d]!=C[%d,%d,%d]"%(i,j,k,k,i,j)
674                                  out=False
675                else:
676                   for j in range(self.getDim()):
677                      if util.Lsup(B[j]-C[j])>tol:
678                         if verbose: print "non-symmetric PDE because B[%d]!=C[%d]"%(j,j)
679                         out=False
680             if self.getNumSolutions()>1:
681               D=self.getCoefficientOfPDE("D")
682               if not D.isEmpty():
683                 tol=util.Lsup(D)*self.TOL
684                 for i in range(self.getNumEquations()):
685                    for k in range(self.getNumSolutions()):
686                      if util.Lsup(D[i,k]-D[k,i])>tol:
687                          if verbose: print "non-symmetric PDE because D[%d,%d]!=D[%d,%d]"%(i,k,k,i)
688                          out=False
689                
690          return out
691    
692       def getFlux(self,u):
693         """         """
694         @brief         returns the flux J_ij for a given u
695    
696           \f[
697           J_ij=A_{ijkl}u_{k,l}+B_{ijk}u_k-X_{ij}
698           \f]
699    
700           @param u: argument of the operator
701         """         """
702         if self.__righthandside.isEmpty():         raise SystemError,"getFlux is not implemented yet"
703             self.__righthandside=self.__makeNewRightHandSide()         return None
704             if self.debug() : print "PDE Debug: New right hand side allocated"  
705         else:     def applyOperator(self,u):
706             print "fix self.__righthandside*=0"         """
707             self.__righthandside*=0.         applies the operator of the PDE to a given solution u in weak from
708             if self.debug() : print "PDE Debug: Right hand side reset to zero"  
709         return  self.__righthandside         @param u: argument of the operator
710           """
711           return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())
712                                                                                                                                                              
713       def getResidual(self,u):
714           """
715           return the residual of u in the weak from
716    
717           @param u:
718           """
719           return self.applyOperator(u)-self.getRightHandSide()
720    
721       def __setValue(self,**coefficients):
722          """
723          sets new values to coefficient
724    
725          @param coefficients:
726          """
727          # check if the coefficients are  legal:
728          for i in coefficients.iterkeys():
729             if not self.hasCoefficient(i):
730                raise ValueError,"Attempt to set unknown coefficient %s"%i
731          # if the number of unknowns or equations is still unknown we try to estimate them:
732          if self.__numEquations<1 or self.__numSolutions<1:
733             for i,d in coefficients.iteritems():
734                if hasattr(d,"shape"):
735                    s=d.shape
736                elif hasattr(d,"getShape"):
737                    s=d.getShape()
738                else:
739                    s=numarray.array(d).shape
740                if s!=None:
741                    # get number of equations and number of unknowns:
742                    res=self.COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(s,self.getDim())
743                    if res==None:
744                        raise ValueError,"Illegal shape %s of coefficient %s"%(s,i)
745                    else:
746                        if self.__numEquations<1: self.__numEquations=res[0]
747                        if self.__numSolutions<1: self.__numSolutions=res[1]
748          if self.__numEquations<1: raise ValueError,"unidententified number of equations"
749          if self.__numSolutions<1: raise ValueError,"unidententified number of solutions"
750          # now we check the shape of the coefficient if numEquations and numSolutions are set:
751          for i,d in coefficients.iteritems():
752            if d==None:
753                 d2=escript.Data()
754            elif isinstance(d,escript.Data):
755                 if d.isEmpty():
756                    d2=d
757                 else:
758                    d2=escript.Data(d,self.getFunctionSpaceForCoefficient(i))
759            else:
760                  d2=escript.Data(d,self.getFunctionSpaceForCoefficient(i))
761            if not d2.isEmpty():
762               if not self.getShapeOfCoefficient(i)==d2.getShape():
763                   raise ValueError,"Expected shape for coefficient %s is %s but actual shape is %s."%(i,self.getShapeOfCoefficient(i),d2.getShape())
764            # overwrite new values:
765            if self.debug(): print "PDE Debug: Coefficient %s has been altered."%i
766            self.COEFFICIENTS[i].setValue(d2)
767            self.alteredCoefficient(i)
768          
769          # reset the HomogeneousConstraintFlag:
770          self.__setHomogeneousConstraintFlag()
771          if len(coefficients)>0 and not self.isUsingLumping() and not self.__homogeneous_constraint: self.__rebuildSystem()
772    
773       def __setHomogeneousConstraintFlag(self):
774          """
775          checks if the constraints are homogeneous and sets self.__homogeneous_constraint accordingly.
776          """
777          self.__homogeneous_constraint=True
778          q=self.getCoefficientOfPDE("q")
779          r=self.getCoefficientOfPDE("r")
780          if not q.isEmpty() and not r.isEmpty():
781             if (q*r).Lsup()>=1.e-13*r.Lsup(): self.__homogeneous_constraint=False
782          if self.debug():
783               if self.__homogeneous_constraint:
784                   print "PDE Debug: Constraints are homogeneous."
785               else:
786                   print "PDE Debug: Constraints are inhomogeneous."
787    
788    
789     # ==== rebuild switches =====================================================================     # ==== rebuild switches =====================================================================
790     def __rebuildSolution(self,deep=False):     def __rebuildSolution(self,deep=False):
791         """         """
792         @brief indicates the PDE has to be reolved if the solution is requested         indicates the PDE has to be reolved if the solution is requested
793         """         """
794         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."
795         self.__solution_isValid=False         self.__solution_isValid=False
796         if deep: self.__solution=escript.Data(deep)         if deep: self.__solution=escript.Data()
797    
798    
799     def __rebuildOperator(self,deep=False):     def __rebuildOperator(self,deep=False):
800         """         """
801         @brief indicates the operator has to be rebuilt next time it is used         indicates the operator has to be rebuilt next time it is used
802         """         """
803         if self.__operator_isValid and self.debug() : print "PDE Debug: Operator has to be rebuilt."         if self.__operator_isValid and self.debug() : print "PDE Debug: Operator has to be rebuilt."
804         self.__rebuildSolution(deep)         self.__rebuildSolution(deep)
# Line 675  class LinearPDE: Line 807  class LinearPDE:
807    
808     def __rebuildRightHandSide(self,deep=False):     def __rebuildRightHandSide(self,deep=False):
809         """         """
810         @brief indicates the right hand side has to be rebuild next time it is used         indicates the right hand side has to be rebuild next time it is used
811         """         """
812         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."
813         self.__rebuildSolution(deep)         self.__rebuildSolution(deep)
814         self.__righthandside_isValid=False         self.__righthandside_isValid=False
        if not self.__homogeneous_constraint: self.__rebuildOperator()  
815         if deep: self.__righthandside=escript.Data()         if deep: self.__righthandside=escript.Data()
816    
817     def __rebuildSystem(self,deep=False):     def __rebuildSystem(self,deep=False):
818         """         """
819         @brief annonced that all coefficient name has been changed         annonced that all coefficient name has been changed
820         """         """
821         self.__rebuildSolution(deep)         self.__rebuildSolution(deep)
822         self.__rebuildOperator(deep)         self.__rebuildOperator(deep)
# Line 693  class LinearPDE: Line 824  class LinearPDE:
824        
825     def __checkMatrixType(self):     def __checkMatrixType(self):
826       """       """
827       @brief reassess the matrix type and, if needed, initiates an operator rebuild       reassess the matrix type and, if needed, initiates an operator rebuild
828       """       """
829       new_matrix_type=self.getDomain().getSystemMatrixTypeId(self.getSolverMethod(),self.isSymmetric())       new_matrix_type=self.getDomain().getSystemMatrixTypeId(self.getSolverMethod(),self.isSymmetric())
830       if not new_matrix_type==self.__matrix_type:       if not new_matrix_type==self.__matrix_type:
# Line 702  class LinearPDE: Line 833  class LinearPDE:
833           self.__rebuildOperator(deep=True)           self.__rebuildOperator(deep=True)
834    
835     #============ assembling =======================================================     #============ assembling =======================================================
836     def __copyConstraint(self,u):     def __copyConstraint(self):
837        """        """
838        @brief copies the constrint condition into u        copies the constrint condition into u
839        """        """
840        q=self.getCoefficient("q")        if not self.__righthandside.isEmpty():
841        r=self.getCoefficient("r")           q=self.getCoefficientOfPDE("q")
842        if not q.isEmpty():           r=self.getCoefficientOfPDE("r")
843            if r.isEmpty():           if not q.isEmpty():
844               r2=escript.Data(0,u.getShape(),u.getFunctionSpace())               if r.isEmpty():
845            else:                  r2=escript.Data(0,self.__righthandside.getShape(),self.__righthandside.getFunctionSpace())
846               r2=escript.Data(r,u.getFunctionSpace())               else:
847            u.copyWithMask(r2,escript.Data(q,u.getFunctionSpace()))                  r2=escript.Data(r,self.__righthandside.getFunctionSpace())
848                 self.__righthandside.copyWithMask(r2,escript.Data(q,self.__righthandside.getFunctionSpace()))
849    
850     def __applyConstraint(self,rhs_update=True):     def __applyConstraint(self):
851         """         """
852         @brief applies the constraints  defined by q and r to the system         applies the constraints defined by q and r to the system
853         """         """
854         q=self.getCoefficient("q")         q=self.getCoefficientOfPDE("q")
855         r=self.getCoefficient("r")         r=self.getCoefficientOfPDE("r")
856         if not q.isEmpty() and not self.__operator.isEmpty():         if not q.isEmpty() and not self.__operator.isEmpty():
857            # 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:
858            row_q=escript.Data(q,self.getFunctionSpaceForEquation())            row_q=escript.Data(q,self.getFunctionSpaceForEquation())
859            col_q=escript.Data(q,self.getFunctionSpaceForSolution())            col_q=escript.Data(q,self.getFunctionSpaceForSolution())
860            u=self.__makeNewSolution()            u=self.__getNewSolution()
861            if r.isEmpty():            if r.isEmpty():
862               r_s=self.__makeNewSolution()               r_s=self.__getNewSolution()
863            else:            else:
864               r_s=escript.Data(r,self.getFunctionSpaceForSolution())               r_s=escript.Data(r,self.getFunctionSpaceForSolution())
865            u.copyWithMask(r_s,col_q)            u.copyWithMask(r_s,col_q)
866            if not self.__righthandside.isEmpty() and rhs_update: self.__righthandside-=self.__operator*u            if self.isUsingLumping():
867            self.__operator.nullifyRowsAndCols(row_q,col_q,1.)               self.__operator.copyWithMask(escript.Data(1,q.getShape(),self.getFunctionSpaceForEquation()),row_q)
868              else:
869                 if not self.__righthandside.isEmpty(): self.__righthandside-=self.__operator*u
870                 self.__operator.nullifyRowsAndCols(row_q,col_q,1.)
871    
872     def getOperator(self):     def getSystem(self):
873         """         """
874         @brief returns the operator of the PDE         return the operator and right hand side of the PDE
875         """         """
876         if not self.__operator_isValid:         if not self.__operator_isValid or not self.__righthandside_isValid:
877             # some Constraints are applying for a lumpled stifness matrix:            if self.isUsingLumping():
878             if self.isUsingLumping():                if not self.__operator_isValid:
879                if self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():                   if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():
880                         raise TypeError,"Lumped matrix requires same order for equations and unknowns"                         raise TypeError,"Lumped matrix requires same order for equations and unknowns"
881                if not self.getCoefficient("A").isEmpty():                   if not self.getCoefficientOfPDE("A").isEmpty():
882                         raise Warning,"Lumped matrix does not allow coefficient A"                            raise Warning,"Lumped matrix does not allow coefficient A"
883                if not self.getCoefficient("B").isEmpty():                   if not self.getCoefficientOfPDE("B").isEmpty():
884                         raise Warning,"Lumped matrix does not allow coefficient B"                            raise Warning,"Lumped matrix does not allow coefficient B"
885                if not self.getCoefficient("C").isEmpty():                   if not self.getCoefficientOfPDE("C").isEmpty():
886                         raise Warning,"Lumped matrix does not allow coefficient C"                            raise Warning,"Lumped matrix does not allow coefficient C"
887                if not self.getCoefficient("D").isEmpty():                   if self.debug() : print "PDE Debug: New lumped operator is built."
888                         raise Warning,"Lumped matrix does not allow coefficient D"                   mat=self.__getNewOperator()
889                if self.debug() : print "PDE Debug: New lumped operator is built."                   self.getDomain().addPDEToSystem(mat,escript.Data(), \
890                mat=self.__makeNewOperator(self)                             self.getCoefficientOfPDE("A"), \
891             else:                             self.getCoefficientOfPDE("B"), \
892                if self.debug() : print "PDE Debug: New operator is built."                             self.getCoefficientOfPDE("C"), \
893                mat=self.__getFreshOperator()                             self.getCoefficientOfPDE("D"), \
894                               escript.Data(), \
895             self.getDomain().addPDEToSystem(mat,escript.Data(), \                             escript.Data(), \
896                          self.getCoefficient("A"), \                             self.getCoefficientOfPDE("d"), \
897                          self.getCoefficient("B"), \                             escript.Data(),\
898                          self.getCoefficient("C"), \                             self.getCoefficientOfPDE("d_contact"), \
899                          self.getCoefficient("D"), \                             escript.Data())
900                          escript.Data(), \                   self.__operator=mat*escript.Data(1,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)
901                          escript.Data(), \                   self.__applyConstraint()
902                          self.getCoefficient("d"), \                   self.__operator_isValid=True
903                          escript.Data(),\                if not self.__righthandside_isValid:
904                          self.getCoefficient("d_contact"), \                   if self.debug() : print "PDE Debug: New right hand side is built."
905                          escript.Data())                   self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \
906             if self.isUsingLumping():                                 self.getCoefficientOfPDE("X"), \
907                self.__operator=mat*escript.Data(1,(self.getNumSolutions(),),self.getFunctionSpaceOfSolution(),True)                                 self.getCoefficientOfPDE("Y"),\
908             else:                                 self.getCoefficientOfPDE("y"),\
909                self.__applyConstraint(rhs_update=False)                                 self.getCoefficientOfPDE("y_contact"))
910             self.__operator_isValid=True                   self.__copyConstraint()
911         return self.__operator                   self.__righthandside_isValid=True
912              else:
913     def getRightHandSide(self,ignoreConstraint=False):               if not self.__operator_isValid and not self.__righthandside_isValid:
914                     if self.debug() : print "PDE Debug: New system is built."
915                     self.getDomain().addPDEToSystem(self.__makeFreshOperator(),self.__makeFreshRightHandSide(), \
916                                   self.getCoefficientOfPDE("A"), \
917                                   self.getCoefficientOfPDE("B"), \
918                                   self.getCoefficientOfPDE("C"), \
919                                   self.getCoefficientOfPDE("D"), \
920                                   self.getCoefficientOfPDE("X"), \
921                                   self.getCoefficientOfPDE("Y"), \
922                                   self.getCoefficientOfPDE("d"), \
923                                   self.getCoefficientOfPDE("y"), \
924                                   self.getCoefficientOfPDE("d_contact"), \
925                                   self.getCoefficientOfPDE("y_contact"))
926                     self.__applyConstraint()
927                     self.__copyConstraint()
928                     self.__operator_isValid=True
929                     self.__righthandside_isValid=True
930                 elif not self.__righthandside_isValid:
931                     if self.debug() : print "PDE Debug: New right hand side is built."
932                     self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \
933                                   self.getCoefficientOfPDE("X"), \
934                                   self.getCoefficientOfPDE("Y"),\
935                                   self.getCoefficientOfPDE("y"),\
936                                   self.getCoefficientOfPDE("y_contact"))
937                     self.__copyConstraint()
938                     self.__righthandside_isValid=True
939                 elif not self.__operator_isValid:
940                     if self.debug() : print "PDE Debug: New operator is built."
941                     self.getDomain().addPDEToSystem(self.__makeFreshOperator(),escript.Data(), \
942                                self.getCoefficientOfPDE("A"), \
943                                self.getCoefficientOfPDE("B"), \
944                                self.getCoefficientOfPDE("C"), \
945                                self.getCoefficientOfPDE("D"), \
946                                escript.Data(), \
947                                escript.Data(), \
948                                self.getCoefficientOfPDE("d"), \
949                                escript.Data(),\
950                                self.getCoefficientOfPDE("d_contact"), \
951                                escript.Data())
952                     self.__applyConstraint()
953                     self.__operator_isValid=True
954           return (self.__operator,self.__righthandside)
955       def getOperator(self):
956         """         """
957         @brief returns the right hand side of the PDE         returns the operator of the PDE
   
        @param ignoreConstraint  
958         """         """
959         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  
960    
961     def getSystem(self):     def getRightHandSide(self):
962         """         """
963         @brief         returns the right hand side of the PDE
964         """         """
965         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)  
966    
967     def solve(self,**options):     def solve(self,**options):
968        """        """
969        @brief solve the PDE        solve the PDE
970    
971        @param options        @param options:
972        """        """
973        mat,f=self.getSystem()        mat,f=self.getSystem()
974        if self.isUsingLumping():        if self.isUsingLumping():
975           out=f/mat           out=f/mat
          self.__copyConstraint(out)  
976        else:        else:
977           options[util.TOLERANCE_KEY]=self.getTolerance()           options[util.TOLERANCE_KEY]=self.getTolerance()
978           options[util.METHOD_KEY]=self.getSolverMethod()           options[util.METHOD_KEY]=self.getSolverMethod()
# Line 843  class LinearPDE: Line 983  class LinearPDE:
983    
984     def getSolution(self,**options):     def getSolution(self,**options):
985         """         """
986         @brief returns the solution of the PDE         returns the solution of the PDE
987    
988         @param options         @param options:
989         """         """
990         if not self.__solution_isValid:         if not self.__solution_isValid:
991             if self.debug() : print "PDE Debug: PDE is resolved."             if self.debug() : print "PDE Debug: PDE is resolved."
992             self.__solution=self.solve(**options)             self.__solution=self.solve(**options)
993             self.__solution_isValid=True             self.__solution_isValid=True
994         return self.__solution         return self.__solution
    #============ some serivice functions  =====================================================  
    def getDomain(self):  
      """  
      @brief returns the domain of the PDE  
      """  
      return self.__domain  
995    
    def getDim(self):  
      """  
      @brief returns the spatial dimension of the PDE  
      """  
      return self.getDomain().getDim()  
996    
    def getNumEquations(self):  
      """  
      @brief returns the number of equations  
      """  
      if self.__numEquations>0:  
          return self.__numEquations  
      else:  
          raise ValueError,"Number of equations is undefined. Please specify argument numEquations."  
997    
998     def getNumSolutions(self):  def ELMAN_RAMAGE(P):
999       """       """   """
1000       @brief returns the number of unknowns       return (P-1.).wherePositive()*0.5*(1.-1./(P+1.e-15))
1001       """  def SIMPLIFIED_BROOK_HUGHES(P):
1002       if self.__numSolutions>0:       """   """
1003          return self.__numSolutions       c=(P-3.).whereNegative()
1004       else:       return P/6.*c+1./2.*(1.-c)
1005          raise ValueError,"Number of solution is undefined. Please specify argument numSolutions."  
1006    def HALF(P):
1007        """ """
1008        return escript.Scalar(0.5,P.getFunctionSpace())
1009    
1010    class AdvectivePDE(LinearPDE):
1011       """
1012       Class to handle a linear PDE dominated by advective terms:
1013      
1014       class to define a linear PDE of the form
1015    
1016     def checkSymmetry(self):     \f[
1017        """     -(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
1018        @brief returns if the Operator is symmetric. This is a very expensive operation!!! The symmetry flag is not altered.     \f]
       """  
       raise SystemError,"checkSymmetry is not implemented yet"  
1019    
1020        return None     with boundary conditons:
1021    
1022     def getFlux(self,u):     \f[
1023         """     n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_{ik}u_k = - n_j*X_{ij} + y_i
1024         @brief returns the flux J_ij for a given u     \f]
1025    
1026              J_ij=A_{ijkl}u_{k,l}+B_{ijk}u_k-X_{ij}     and contact conditions
1027    
1028         @param u argument of the operator     \f[
1029       n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d^{contact}_{ik}[u_k] = - n_j*X_{ij} + y^{contact}_{i}
1030       \f]
1031    
1032         """     and constraints:
        raise SystemError,"getFlux is not implemented yet"  
        return None  
1033    
1034     def applyOperator(self,u):     \f[
1035         """     u_i=r_i \quad \mathrm{where} \quad q_i>0
1036         @brief applies the operator of the PDE to a given solution u in weak from     \f]
1037       """
1038       def __init__(self,domain,numEquations=0,numSolutions=0,xi=ELMAN_RAMAGE):
1039          LinearPDE.__init__(self,domain,numEquations,numSolutions)
1040          self.__xi=xi
1041          self.__Xi=escript.Data()
1042    
1043       def __calculateXi(self,peclet_factor,Z,h):
1044           Z_max=util.Lsup(Z)
1045           if Z_max>0.:
1046              return h*self.__xi(Z*peclet_factor)/(Z+Z_max*self.TOL)
1047           else:
1048              return 0.
1049    
1050       def setValue(self,**args):
1051           if "A" in args.keys()   or "B" in args.keys() or "C" in args.keys(): self.__Xi=escript.Data()
1052           self._LinearPDE__setValue(**args)
1053              
1054       def getXi(self):
1055          if self.__Xi.isEmpty():
1056             B=self.getCoefficient("B")
1057             C=self.getCoefficient("C")
1058             A=self.getCoefficient("A")
1059             h=self.getDomain().getSize()
1060             self.__Xi=escript.Scalar(0.,self.getFunctionSpaceForCoefficient("A"))
1061             if not C.isEmpty() or not B.isEmpty():
1062                if not C.isEmpty() and not B.isEmpty():
1063                    Z2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))
1064                    if self.getNumEquations()>1:
1065                       if self.getNumSolutions()>1:
1066                          for i in range(self.getNumEquations()):
1067                             for k in range(self.getNumSolutions()):
1068                                for l in range(self.getDim()): Z2+=(C[i,k,l]-B[i,l,k])**2
1069                       else:
1070                          for i in range(self.getNumEquations()):
1071                             for l in range(self.getDim()): Z2+=(C[i,l]-B[i,l])**2
1072                    else:
1073                       if self.getNumSolutions()>1:
1074                          for k in range(self.getNumSolutions()):
1075                             for l in range(self.getDim()): Z2+=(C[k,l]-B[l,k])**2
1076                       else:
1077                          for l in range(self.getDim()): Z2+=(C[l]-B[l])**2
1078                    length_of_Z=util.sqrt(Z2)
1079                elif C.isEmpty():
1080                  length_of_Z=util.length(B)
1081                else:
1082                  length_of_Z=util.length(C)
1083    
1084         @param u argument of the operator              Z_max=util.Lsup(length_of_Z)
1085                if Z_max>0.:
1086                   length_of_A=util.length(A)
1087                   A_max=util.Lsup(length_of_A)
1088                   if A_max>0:
1089                        inv_A=1./(length_of_A+A_max*self.TOL)
1090                   else:
1091                        inv_A=1./self.TOL
1092                   peclet_number=length_of_Z*h/2*inv_A
1093                   xi=self.__xi(peclet_number)
1094                   self.__Xi=h*xi/(length_of_Z+Z_max*self.TOL)
1095                   print "@ preclet number = %e"%util.Lsup(peclet_number),util.Lsup(xi),util.Lsup(length_of_Z)
1096          return self.__Xi
1097          
1098    
1099         """     def getCoefficientOfPDE(self,name):
1100         return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())       """
1101                                                                                                                                                                   return the value of the coefficient name of the general PDE
1102     def getResidual(self,u):  
1103         """       @param name:
1104         @brief return the residual of u in the weak from       """
1105         if not self.getNumEquations() == self.getNumSolutions():
1106              raise ValueError,"AdvectivePDE expects the number of solution componets and the number of equations to be equal."
1107    
1108         if name == "A" :
1109             A=self.getCoefficient("A")
1110             B=self.getCoefficient("B")
1111             C=self.getCoefficient("C")
1112             if B.isEmpty() and C.isEmpty():
1113                Aout=A
1114             else:
1115                if A.isEmpty():
1116                   Aout=self.createNewCoefficient("A")
1117                else:
1118                   Aout=A[:]
1119                Xi=self.getXi()
1120                if self.getNumEquations()>1:
1121                    for i in range(self.getNumEquations()):
1122                       for j in range(self.getDim()):
1123                          for k in range(self.getNumSolutions()):
1124                             for l in range(self.getDim()):
1125                                if not C.isEmpty() and not B.isEmpty():
1126                                   for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*(C[p,i,j]-B[p,j,i])*(C[p,k,l]-B[p,l,k])
1127                                elif C.isEmpty():
1128                                   for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*B[p,j,i]*B[p,l,k]
1129                                else:
1130                                   for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*C[p,i,j]*C[p,k,l]
1131                else:
1132                    for j in range(self.getDim()):
1133                       for l in range(self.getDim()):
1134                          if not C.isEmpty() and not B.isEmpty():
1135                              Aout[j,l]+=Xi*(C[j]-B[j])*(C[l]-B[l])
1136                          elif C.isEmpty():
1137                              Aout[j,l]+=Xi*B[j]*B[l]
1138                          else:
1139                              Aout[j,l]+=Xi*C[j]*C[l]
1140             return Aout
1141         elif name == "B" :
1142             B=self.getCoefficient("B")
1143             C=self.getCoefficient("C")
1144             D=self.getCoefficient("D")
1145             if C.isEmpty() or D.isEmpty():
1146                Bout=B
1147             else:
1148                Xi=self.getXi()
1149                if B.isEmpty():
1150                    Bout=self.createNewCoefficient("B")
1151                else:
1152                    Bout=B[:]
1153                if self.getNumEquations()>1:
1154                   for k in range(self.getNumSolutions()):
1155                      for p in range(self.getNumEquations()):
1156                         tmp=Xi*D[p,k]
1157                         for i in range(self.getNumEquations()):
1158                            for j in range(self.getDim()):
1159                               Bout[i,j,k]+=tmp*C[p,i,j]
1160                else:
1161                   tmp=Xi*D
1162                   for j in range(self.getDim()): Bout[j]+=tmp*C[j]
1163             return Bout
1164         elif name == "C" :
1165             B=self.getCoefficient("B")
1166             C=self.getCoefficient("C")
1167             D=self.getCoefficient("D")
1168             if B.isEmpty() or D.isEmpty():
1169                Cout=C
1170             else:
1171                Xi=self.getXi()
1172                if C.isEmpty():
1173                    Cout=self.createNewCoefficient("C")
1174                else:
1175                    Cout=C[:]
1176                if self.getNumEquations()>1:
1177                   for k in range(self.getNumSolutions()):
1178                       for p in range(self.getNumEquations()):
1179                          tmp=Xi*D[p,k]
1180                          for i in range(self.getNumEquations()):
1181                            for l in range(self.getDim()):
1182                                     Cout[i,k,l]+=tmp*B[p,l,i]
1183                else:
1184                   tmp=Xi*D
1185                   for j in range(self.getDim()): Cout[j]+=tmp*B[j]
1186             return Cout
1187         elif name == "D" :
1188             return self.getCoefficient("D")
1189         elif name == "X" :
1190             X=self.getCoefficient("X")
1191             Y=self.getCoefficient("Y")
1192             B=self.getCoefficient("B")
1193             C=self.getCoefficient("C")
1194             if Y.isEmpty() or (B.isEmpty() and C.isEmpty()):
1195                Xout=X
1196             else:
1197                if X.isEmpty():
1198                    Xout=self.createNewCoefficient("X")
1199                else:
1200                    Xout=X[:]
1201                Xi=self.getXi()
1202                if self.getNumEquations()>1:
1203                     for p in range(self.getNumEquations()):
1204                        tmp=Xi*Y[p]
1205                        for i in range(self.getNumEquations()):
1206                           for j in range(self.getDim()):
1207                              if not C.isEmpty() and not B.isEmpty():
1208                                 Xout[i,j]+=tmp*(C[p,i,j]-B[p,j,i])
1209                              elif C.isEmpty():
1210                                 Xout[i,j]-=tmp*B[p,j,i]
1211                              else:
1212                                 Xout[i,j]+=tmp*C[p,i,j]
1213                else:
1214                     tmp=Xi*Y
1215                     for j in range(self.getDim()):
1216                        if not C.isEmpty() and not B.isEmpty():
1217                           Xout[j]+=tmp*(C[j]-B[j])
1218                        elif C.isEmpty():
1219                           Xout[j]-=tmp*B[j]
1220                        else:
1221                           Xout[j]+=tmp*C[j]
1222             return Xout
1223         elif name == "Y" :
1224             return self.getCoefficient("Y")
1225         elif name == "d" :
1226             return self.getCoefficient("d")
1227         elif name == "y" :
1228             return self.getCoefficient("y")
1229         elif name == "d_contact" :
1230             return self.getCoefficient("d_contact")
1231         elif name == "y_contact" :
1232             return self.getCoefficient("y_contact")
1233         elif name == "r" :
1234             return self.getCoefficient("r")
1235         elif name == "q" :
1236             return self.getCoefficient("q")
1237         else:
1238             raise SystemError,"unknown PDE coefficient %s",name
1239    
        @param u  
        """  
        return self.applyOperator(u)-self.getRightHandSide()  
1240    
1241  class Poisson(LinearPDE):  class Poisson(LinearPDE):
1242     """     """
1243     @brief Class to define a Poisson equstion problem:     Class to define a Poisson equation problem:
1244                                                                                                                                                                
1245     class to define a linear PDE of the form     class to define a linear PDE of the form
1246                                                                                                                                                                   \f[
1247          -u_{,jj} = f     -u_{,jj} = f
1248                                                                                                                                                                   \f]
1249       with boundary conditons:  
1250                                                                                                                                                                   with boundary conditons:
1251          n_j*u_{,j} = 0  
1252                                                                                                                                                                   \f[
1253      and constraints:     n_j*u_{,j} = 0
1254                                                                                                                                                                   \f]
1255           u=0 where q>0  
1256                                                                                                                                                                   and constraints:
1257    
1258       \f[
1259       u=0 \quad \mathrm{where} \quad q>0
1260       \f]
1261     """     """
1262    
1263     def __init__(self,domain=None,f=escript.Data(),q=escript.Data()):     def __init__(self,domain,f=escript.Data(),q=escript.Data()):
1264         LinearPDE.__init__(self,domain=identifyDomain(domain,{"f":f, "q":q}))         LinearPDE.__init__(self,domain,1,1)
1265         self._setValue(A=numarray.identity(self.getDomain().getDim()))         self.COEFFICIENTS={
1266           "f"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1267           "q"         : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.BOTH)}
1268         self.setSymmetryOn()         self.setSymmetryOn()
1269         self.setValue(f,q)         self.setValue(f,q)
1270    
1271     def setValue(self,f=escript.Data(),q=escript.Data()):     def setValue(self,f=escript.Data(),q=escript.Data()):
1272         self._setValue(Y=f,q=q)         """set value of PDE parameters f and q"""
1273           self._LinearPDE__setValue(f=f,q=q)
1274                                                                                                                                                              
1275       def getCoefficientOfPDE(self,name):
1276         """
1277         return the value of the coefficient name of the general PDE
1278    
1279         @param name:
1280         """
1281         if name == "A" :
1282             return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))
1283         elif name == "B" :
1284             return escript.Data()
1285         elif name == "C" :
1286             return escript.Data()
1287         elif name == "D" :
1288             return escript.Data()
1289         elif name == "X" :
1290             return escript.Data()
1291         elif name == "Y" :
1292             return self.getCoefficient("f")
1293         elif name == "d" :
1294             return escript.Data()
1295         elif name == "y" :
1296             return escript.Data()
1297         elif name == "d_contact" :
1298             return escript.Data()
1299         elif name == "y_contact" :
1300             return escript.Data()
1301         elif name == "r" :
1302             return escript.Data()
1303         elif name == "q" :
1304             return self.getCoefficient("q")
1305         else:
1306             raise SystemError,"unknown PDE coefficient %s",name
1307    
1308    class LameEquation(LinearPDE):
1309       """
1310       Class to define a Lame equation problem:
1311    
1312       class to define a linear PDE of the form
1313       \f[
1314       -(\lambda (u_{i,j}+u_{j,i}))_{,j} - \mu u_{j,ji}} = F_i -\sigma_{ij,j}
1315       \f]
1316    
1317       with boundary conditons:
1318    
1319       \f[
1320       n_j(\lambda(u_{i,j}+u_{j,i})-sigma_{ij}) + n_i\mu u_{j,j} = f_i
1321       \f]
1322    
1323       and constraints:
1324    
1325       \f[
1326       u_i=r_i \quad \mathrm{where} \quad q_i>0
1327       \f]
1328       """
1329    
1330       def __init__(self,domain,f=escript.Data(),q=escript.Data()):
1331           LinearPDE.__init__(self,domain,domain.getDim(),domain.getDim())
1332           self.COEFFICIENTS={
1333           "lame_lambda"  : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),
1334           "lame_mu"      : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),
1335           "F"            : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM),PDECoefficient.RIGHTHANDSIDE),
1336           "sigma"        : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.EQUATION),PDECoefficient.RIGHTHANDSIDE),
1337           "f"            : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1338           "r"            : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.BOTH),
1339           "q"            : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.BOTH)}
1340           self.setSymmetryOn()
1341    
1342       def setValue(self,lame_lambda=escript.Data(),lame_mu=escript.Data(),F=escript.Data(),sigma=escript.Data(),f=escript.Data(),r=escript.Data(),q=escript.Data()):
1343           """set value of PDE parameters"""
1344           self._LinearPDE__setValue(lame_lambda=lame_lambda, \
1345                                     lame_mu=lame_mu, \
1346                                     F=F, \
1347                                     sigma=sigma, \
1348                                     f=f, \
1349                                     r=r, \
1350                                     q=q)
1351       def getCoefficientOfPDE(self,name):
1352         """
1353         return the value of the coefficient name of the general PDE
1354    
1355         @param name:
1356         """
1357         if name == "A" :
1358             A =self.createNewCoefficient("A")
1359             for i in range(self.getDim()):
1360               for j in range(self.getDim()):
1361                 out[i,i,j,j] += self.getCoefficient("lame_mu")
1362                 out[i,j,j,i] += self.getCoefficient("lame_lambda")
1363                 out[i,j,i,j] += self.getCoefficient("lame_lambda")
1364             return out
1365         elif name == "B" :
1366             return escript.Data()
1367         elif name == "C" :
1368             return escript.Data()
1369         elif name == "D" :
1370             return escript.Data()
1371         elif name == "X" :
1372             return self.getCoefficient("sigma")
1373         elif name == "Y" :
1374             return self.getCoefficient("F")
1375         elif name == "d" :
1376             return escript.Data()
1377         elif name == "y" :
1378             return self.getCoefficient("f")
1379         elif name == "d_contact" :
1380             return escript.Data()
1381         elif name == "y_contact" :
1382             return escript.Data()
1383         elif name == "r" :
1384             return self.getCoefficient("r")
1385         elif name == "q" :
1386             return self.getCoefficient("q")
1387         else:
1388             raise SystemError,"unknown PDE coefficient %s",name
1389    
1390  # $Log$  # $Log$
1391  # Revision 1.3  2004/12/17 07:43:10  jgs  # Revision 1.9  2005/07/25 05:28:13  jgs
1392  # *** empty log message ***  # Merge of development branch back to main trunk on 2005-07-25
1393    #
1394    # Revision 1.8  2005/06/09 05:37:59  jgs
1395    # Merge of development branch back to main trunk on 2005-06-09
1396    #
1397    # Revision 1.7  2005/05/06 04:26:10  jgs
1398    # Merge of development branch back to main trunk on 2005-05-06
1399    #
1400    # Revision 1.1.2.24  2005/07/22 06:37:11  gross
1401    # some extensions to modellib and linearPDEs
1402    #
1403    # Revision 1.1.2.23  2005/05/13 00:55:20  cochrane
1404    # Fixed up some docstrings.  Moved module-level functions to top of file so
1405    # that epydoc and doxygen can pick them up properly.
1406    #
1407    # Revision 1.1.2.22  2005/05/12 11:41:30  gross
1408    # some basic Models have been added
1409    #
1410    # Revision 1.1.2.21  2005/05/12 07:16:12  cochrane
1411    # Moved ELMAN_RAMAGE, SIMPLIFIED_BROOK_HUGHES, and HALF functions to bottom of
1412    # file so that the AdvectivePDE class is picked up by doxygen.  Some
1413    # reformatting of docstrings.  Addition of code to make equations come out
1414    # as proper LaTeX.
1415    #
1416    # Revision 1.1.2.20  2005/04/15 07:09:08  gross
1417    # some problems with functionspace and linearPDEs fixed.
1418    #
1419    # Revision 1.1.2.19  2005/03/04 05:27:07  gross
1420    # bug in SystemPattern fixed.
1421    #
1422    # Revision 1.1.2.18  2005/02/08 06:16:45  gross
1423    # Bugs in AdvectivePDE fixed, AdvectiveTest is stable but more testing is needed
1424    #
1425    # Revision 1.1.2.17  2005/02/08 05:56:19  gross
1426    # Reference Number handling added
1427    #
1428    # Revision 1.1.2.16  2005/02/07 04:41:28  gross
1429    # some function exposed to python to make mesh merging running
1430    #
1431    # Revision 1.1.2.15  2005/02/03 00:14:44  gross
1432    # timeseries add and ESySParameter.py renames esysXML.py for consistence
1433    #
1434    # Revision 1.1.2.14  2005/02/01 06:44:10  gross
1435    # new implementation of AdvectivePDE which now also updates right hand side. systems of PDEs are still not working
1436    #
1437    # Revision 1.1.2.13  2005/01/25 00:47:07  gross
1438    # updates in the documentation
1439    #
1440    # Revision 1.1.2.12  2005/01/12 01:28:04  matt
1441    # Added createCoefficient method for linearPDEs.
1442    #
1443    # Revision 1.1.2.11  2005/01/11 01:55:34  gross
1444    # a problem in linearPDE class fixed
1445    #
1446    # Revision 1.1.2.10  2005/01/07 01:13:29  gross
1447    # some bugs in linearPDE fixed
1448    #
1449    # Revision 1.1.2.9  2005/01/06 06:24:58  gross
1450    # some bugs in slicing fixed
1451    #
1452    # Revision 1.1.2.8  2005/01/05 04:21:40  gross
1453    # FunctionSpace checking/matchig in slicing added
1454    #
1455    # Revision 1.1.2.7  2004/12/29 10:03:41  gross
1456    # bug in setValue fixed
1457    #
1458    # Revision 1.1.2.6  2004/12/29 05:29:59  gross
1459    # AdvectivePDE successfully tested for Peclet number 1000000. there is still a problem with setValue and Data()
1460    #
1461    # Revision 1.1.2.5  2004/12/29 00:18:41  gross
1462    # AdvectivePDE added
1463    #
1464    # Revision 1.1.2.4  2004/12/24 06:05:41  gross
1465    # some changes in linearPDEs to add AdevectivePDE
1466  #  #
1467  # Revision 1.1.2.3  2004/12/16 00:12:34  gross  # Revision 1.1.2.3  2004/12/16 00:12:34  gross
1468  # __init__ of LinearPDE does not accept any coefficients anymore  # __init__ of LinearPDE does not accept any coefficient anymore
1469  #  #
1470  # Revision 1.1.2.2  2004/12/14 03:55:01  jgs  # Revision 1.1.2.2  2004/12/14 03:55:01  jgs
1471  # *** empty log message ***  # *** empty log message ***

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

  ViewVC Help
Powered by ViewVC 1.1.26