/[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 147 by jgs, Fri Aug 12 01:45:47 2005 UTC revision 148 by jgs, Tue Aug 23 01:24:31 2005 UTC
# Line 11  import util Line 11  import util
11  import numarray  import numarray
12    
13    
14  def _CompTuple2(t1,t2):  class IllegalCoefficient(ValueError):
15     """     """
16     Compare two tuples     raised if an illegal coefficient of the general ar particular PDE is requested.
   
    \param t1 The first tuple  
    \param t2 The second tuple  
17     """     """
18    
19     dif=t1[0]+t1[1]-(t2[0]+t2[1])  class IllegalCoefficientValue(ValueError):
20     if dif<0: return 1     """
21     elif dif>0: return -1     raised if an incorrect value for a coefficient is used.
22     else: return 0     """
   
 def ELMAN_RAMAGE(P):  
     return (P-1.).wherePositive()*0.5*(1.-1./(P+1.e-15))  
   
 def SIMPLIFIED_BROOK_HUGHES(P):  
     c=(P-3.).whereNegative()  
     return P/6.*c+1./2.*(1.-c)  
23    
24  def HALF(P):  class UndefinedPDEError(ValueError):
25      return escript.Scalar(0.5,P.getFunctionSpace())     """
26       raised if a PDE is not fully defined yet.
27       """
28    
29    def _CompTuple2(t1,t2):
30          """
31          Compare two tuples
32      
33          @param t1 The first tuple
34          @param t2 The second tuple
35      
36          """
37      
38          dif=t1[0]+t1[1]-(t2[0]+t2[1])
39          if dif<0: return 1
40          elif dif>0: return -1
41          else: return 0
42      
43  class PDECoefficient:  class PDECoefficient:
44      """      """
45      A class for PDE coefficients      A class for PDE coefficients
# Line 84  class PDECoefficient: Line 90  class PDECoefficient:
90         returns the value of the coefficient:         returns the value of the coefficient:
91         """         """
92         return self.value         return self.value
93        
94      def setValue(self,newValue):      def setValue(self,domain,numEquations=1,numSolutions=1,newValue=None):
95         """         """
96         set the value of the coefficient to new value         set the value of the coefficient to new value
97         """         """
98           if newValue==None:
99               newValue=escript.Data()
100           elif isinstance(newValue,escript.Data):
101               if not newValue.isEmpty():
102                  newValue=escript.Data(newValue,self.getFunctionSpace(domain))
103           else:
104               newValue=escript.Data(newValue,self.getFunctionSpace(domain))
105           if not newValue.isEmpty():
106               if not self.getShape(domain,numEquations,numSolutions)==newValue.getShape():
107                   raise IllegalCoefficientValue,"Expected shape for coefficient %s is %s but actual shape is %s."%(self.getShape(domain,numEquations,numSolutions),newValue.getShape())
108         self.value=newValue         self.value=newValue
109        
110      def isAlteringOperator(self):      def isAlteringOperator(self):
111          """          """
112      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
# Line 109  class PDECoefficient: Line 125  class PDECoefficient:
125          else:          else:
126              return None              return None
127    
128      def estimateNumEquationsAndNumSolutions(self,shape=(),dim=3):      def estimateNumEquationsAndNumSolutions(self,domain,shape=()):
129         """         """
130         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
131    
132         @param shape:         @param shape:
133         @param dim:         @param dim:
134         """         """
135           dim=domain.getDim()
136         if len(shape)>0:         if len(shape)>0:
137             num=max(shape)+1             num=max(shape)+1
138         else:         else:
# Line 126  class PDECoefficient: Line 143  class PDECoefficient:
143               search.append((e,u))               search.append((e,u))
144         search.sort(_CompTuple2)         search.sort(_CompTuple2)
145         for item in search:         for item in search:
146               s=self.buildShape(item[0],item[1],dim)               s=self.getShape(domain,item[0],item[1])
147               if len(s)==0 and len(shape)==0:               if len(s)==0 and len(shape)==0:
148                   return (1,1)                   return (1,1)
149               else:               else:
150                   if s==shape: return item                   if s==shape: return item
151         return None         return None
152    
153      def buildShape(self,e=1,u=1,dim=3):      def getShape(self,domain,numEquations=1,numSolutions=1):
154          """          """
155      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
156    
# Line 141  class PDECoefficient: Line 158  class PDECoefficient:
158      @param u:      @param u:
159      @param dim:      @param dim:
160      """      """
161            dim=domain.getDim()
162          s=()          s=()
163          for i in self.pattern:          for i in self.pattern:
164               if i==self.EQUATION:               if i==self.EQUATION:
165                  if e>1: s=s+(e,)                  if numEquations>1: s=s+(numEquations,)
166               elif i==self.SOLUTION:               elif i==self.SOLUTION:
167                  if u>1: s=s+(u,)                  if numSolutions>1: s=s+(numSolutions,)
168               else:               else:
169                  s=s+(dim,)                  s=s+(dim,)
170          return s          return s
171    
172  class LinearPDE:  class LinearPDE:
173     """     """
174     Class to handle a linear PDE     Class to define a linear PDE of the form
     
    class to define a linear PDE of the form  
175    
176     \f[     \f[
177       -(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
# Line 181  class LinearPDE: Line 197  class LinearPDE:
197    
198     """     """
199     TOL=1.e-13     TOL=1.e-13
200     DEFAULT_METHOD=util.DEFAULT_METHOD     # solver methods
201     DIRECT=util.DIRECT     UNKNOWN=-1
202     CHOLEVSKY=util.CHOLEVSKY     DEFAULT_METHOD=0
203     PCG=util.PCG     DIRECT=1
204     CR=util.CR     CHOLEVSKY=2
205     CGS=util.CGS     PCG=3
206     BICGSTAB=util.BICGSTAB     CR=4
207     SSOR=util.SSOR     CGS=5
208     GMRES=util.GMRES     BICGSTAB=6
209     PRES20=util.PRES20     SSOR=7
210     ILU0=util.ILU0     ILU0=8
211     JACOBI=util.JACOBI     ILUT=9
212       JACOBI=10
213       GMRES=11
214       PRES20=12
215       LUMPING=13
216       # matrix reordering:
217       NO_REORDERING=0
218       MINIMUM_FILL_IN=1
219       NESTED_DISSECTION=2
220       # important keys in the dictonary used to hand over options to the solver library:
221       METHOD_KEY="method"
222       SYMMETRY_KEY="symmetric"
223       TOLERANCE_KEY="tolerance"
224    
225    
226       def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):
227         """
228         initializes a new linear PDE
229    
230         @param domain: domain of the PDE
231         @type domain: L{Domain}
232         @param numEquations: number of equations. If numEquations==None the number of equations
233                              is exracted from the PDE coefficients.
234         @param numSolutions: number of solution components. If  numSolutions==None the number of solution components
235                              is exracted from the PDE coefficients.
236         @param debug: if True debug informations are printed.
237    
    def __init__(self,domain,numEquations=0,numSolutions=0):  
      """  
      initializes a new linear PDE.  
238    
      @param args:  
239       """       """
240       # COEFFICIENTS can be overwritten by subclasses:       #
241       self.__COEFFICIENTS={       #   the coefficients of the general PDE:
242         #
243         self.__COEFFICIENTS_OF_GENEARL_PDE={
244         "A"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM,PDECoefficient.SOLUTION,PDECoefficient.DIM),PDECoefficient.OPERATOR),         "A"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM,PDECoefficient.SOLUTION,PDECoefficient.DIM),PDECoefficient.OPERATOR),
245         "B"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),         "B"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),
246         "C"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION,PDECoefficient.DIM),PDECoefficient.OPERATOR),         "C"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION,PDECoefficient.DIM),PDECoefficient.OPERATOR),
# Line 215  class LinearPDE: Line 254  class LinearPDE:
254         "r"         : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),         "r"         : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
255         "q"         : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.SOLUTION,),PDECoefficient.BOTH)}         "q"         : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.SOLUTION,),PDECoefficient.BOTH)}
256    
257       self.COEFFICIENTS=self.__COEFFICIENTS       # COEFFICIENTS can be overwritten by subclasses:
258         self.COEFFICIENTS=self.__COEFFICIENTS_OF_GENEARL_PDE
259       # initialize attributes       # initialize attributes
260       self.__debug=None       self.__debug=debug
261       self.__domain=domain       self.__domain=domain
262       self.__numEquations=numEquations       self.__numEquations=numEquations
263       self.__numSolutions=numSolutions       self.__numSolutions=numSolutions
264       self.cleanCoefficients()       self.__resetSystem()
   
      self.__operator=escript.Operator()  
      self.__operator_isValid=False  
      self.__righthandside=escript.Data()  
      self.__righthandside_isValid=False  
      self.__solution=escript.Data()  
      self.__solution_isValid=False  
265    
266       # set some default values:       # set some default values:
267       self.__homogeneous_constraint=True       self.__homogeneous_constraint=True
268       self.__row_function_space=escript.Solution(self.__domain)       self.__row_function_space=escript.Solution(self.__domain)
269       self.__column_function_space=escript.Solution(self.__domain)       self.__column_function_space=escript.Solution(self.__domain)
270       self.__tolerance=1.e-8       self.__tolerance=1.e-8
271       self.__solver_method=util.DEFAULT_METHOD       self.__solver_method=self.DEFAULT_METHOD
272       self.__matrix_type=self.__domain.getSystemMatrixTypeId(util.DEFAULT_METHOD,False)       self.__matrix_type=self.__domain.getSystemMatrixTypeId(self.DEFAULT_METHOD,False)
273       self.__sym=False       self.__sym=False
      self.__lumping=False  
274    
275     def createCoefficient(self, name):       self.resetCoefficients()
276         self.trace("PDE Coeffients are %s"%str(self.COEFFICIENTS.keys()))
277       # =============================================================================
278       #    general stuff:
279       # =============================================================================
280       def __str__(self):
281           return "<LinearPDE %d>"%id(self)
282       # =============================================================================
283       #    debug :
284       # =============================================================================
285       def setDebugOn(self):
286         """
287         switches on debugging
288       """       """
289       create a data object corresponding to coefficient name       self.__debug=not None
290       @param name:  
291       def setDebugOff(self):
292       """       """
293       return escript.Data(shape = getShapeOfCoefficient(name), \       switches off debugging
294                           what = getFunctionSpaceForCoefficient(name))       """
295         self.__debug=None
296    
297     def __del__(self):     def trace(self,text):
298       pass       """
299         print the text message if debugging is swiched on.
300    
301     def getCoefficient(self,name):       @param name: name of the coefficient enquired.
302         @type name: C{string}
303       """       """
304       return the value of the parameter name       if self.__debug: print "%s: %s"%(str(self),text)
305    
306       @param name:     # =============================================================================
307       # some service functions:
308       # =============================================================================
309       def getDomain(self):
310       """       """
311       return self.COEFFICIENTS[name].getValue()       returns the domain of the PDE
312        
313         @return : the domain of the PDE
314         @rtype : C{Domain}
315    
    def getCoefficientOfPDE(self,name):  
316       """       """
317       return the value of the coefficient name of the general PDE.       return self.__domain
      This method is called by the assembling routine it can be  
      overwritten to map coefficients of a particualr PDE to the general PDE.  
318    
319       @param name:     def getDim(self):
320       """       """
321       return self.getCoefficient(name)       returns the spatial dimension of the PDE
322    
323     def hasCoefficient(self,name):       @return : the spatial dimension of the PDE domain
324        """       @rtype : C{int}
325        return true if name is the name of a coefficient       """
326         return self.getDomain().getDim()
327    
328        @param name:     def getNumEquations(self):
329        """       """
330        return self.COEFFICIENTS.has_key(name)       returns the number of equations
331    
332     def hasPDECoefficient(self,name):       @return : the number of equations
333        """       @rtype : C{int}
334        return true if name is the name of a coefficient       @raise UndefinedPDEError: if the number of equations is not be specified yet.
335         """
336         if self.__numEquations==None:
337             raise UndefinedPDEError,"Number of equations is undefined. Please specify argument numEquations."
338         else:
339             return self.__numEquations
340    
341        @param name:     def getNumSolutions(self):
342        """       """
343        return self.__COEFFICIENTS.has_key(name)       returns the number of unknowns
344    
345         @return : the number of unknowns
346         @rtype : C{int}
347         @raise UndefinedPDEError: if the number of unknowns is not be specified yet.
348         """
349         if self.__numSolutions==None:
350            raise UndefinedPDEError,"Number of solution is undefined. Please specify argument numSolutions."
351         else:
352            return self.__numSolutions
353    
354     def getFunctionSpaceForEquation(self):     def getFunctionSpaceForEquation(self):
355       """       """
356       return true if the test functions should use reduced order       returns the L{escript.FunctionSpace} used to discretize the equation
357        
358         @return : representation space of equation
359         @rtype : L{escript.FunctionSpace}
360    
361       """       """
362       return self.__row_function_space       return self.__row_function_space
363    
364     def getFunctionSpaceForSolution(self):     def getFunctionSpaceForSolution(self):
365       """       """
366       return true if the interpolation of the solution should use reduced order       returns the L{escript.FunctionSpace} used to represent the solution
367        
368         @return : representation space of solution
369         @rtype : L{escript.FunctionSpace}
370    
371       """       """
372       return self.__column_function_space       return self.__column_function_space
373    
    def setValue(self,**coefficients):  
       """  
       sets new values to coefficients  
   
       @param coefficients:  
       """  
       self.__setValue(**coefficients)  
         
   
    def cleanCoefficients(self):  
      """  
      resets all coefficients to default values.  
      """  
      for i in self.COEFFICIENTS.iterkeys():  
          self.COEFFICIENTS[i].resetValue()  
374    
375     def createNewCoefficient(self,name):     def getOperator(self):
      """  
      returns a new coefficient appropriate for coefficient name:  
      """  
      return escript.Data(0,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))  
         
    def createNewCoefficientOfPDE(self,name):  
      """  
      returns a new coefficient appropriate for coefficient name:  
      """  
      return escript.Data(0,self.getShapeOfCoefficientOfPDE(name),self.getFunctionSpaceForCoefficientOfPDE(name))  
         
    def getShapeOfCoefficientOfPDE(self,name):  
376       """       """
377       return the shape of the coefficient name       provides access to the operator of the PDE
378    
379       @param name:       @return : the operator of the PDE
380         @rtype : L{Operator}
381       """       """
382       if self.hasPDECoefficient(name):       m=self.getSystem()[0]
383          return self.__COEFFICIENTS[name].buildShape(self.getNumEquations(),self.getNumSolutions(),self.getDomain().getDim())       if self.isUsingLumping():
384             return self.copyConstraint(1./m)
385       else:       else:
386          raise ValueError,"Unknown coefficient %s requested"%name           return m
387    
388     def getFunctionSpaceForCoefficientOfPDE(self,name):     def getRightHandSide(self):
389       """       """
390       return the atoms of the coefficient name       provides access to the right hand side of the PDE
391    
392       @param name:       @return : the right hand side of the PDE
393         @rtype : L{escript.Data}
394       """       """
395       if self.hasPDECoefficient(name):       r=self.getSystem()[1]
396          return self.__COEFFICIENTS[name].getFunctionSpace(self.getDomain())       if self.isUsingLumping():
397             return self.copyConstraint(r)
398       else:       else:
399          raise ValueError,"unknown coefficient %s requested"%name           return r
400    
401     def getShapeOfCoefficient(self,name):     def applyOperator(self,u=None):
402       """       """
403       return the shape of the coefficient name       applies the operator of the PDE to a given u or the solution of PDE if u is not present.
404    
405       @param name:       @param u: argument of the operator. It must be representable in C{elf.getFunctionSpaceForSolution()}. If u is not present or equals L{None}
406                   the current solution is used.
407         @type u: L{escript.Data} or None
408         @return : image of u
409         @rtype : L{escript.Data}
410       """       """
411       if self.hasCoefficient(name):       if u==None:
412          return self.COEFFICIENTS[name].buildShape(self.getNumEquations(),self.getNumSolutions(),self.getDomain().getDim())            return self.getOperator()*self.getSolution()
413       else:       else:
414          raise ValueError,"Unknown coefficient %s requested"%name          self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())
415    
416     def getFunctionSpaceForCoefficient(self,name):     def getResidual(self,u=None):
417       """       """
418       return the atoms of the coefficient name       return the residual of u or the current solution if u is not present.
419    
420       @param name:       @param u: argument in the residual calculation. It must be representable in C{elf.getFunctionSpaceForSolution()}. If u is not present or equals L{None}
421                   the current solution is used.
422         @type u: L{escript.Data} or None
423         @return : residual of u
424         @rtype : L{escript.Data}
425       """       """
426       if self.hasCoefficient(name):       return self.applyOperator(u)-self.getRightHandSide()
         return self.COEFFICIENTS[name].getFunctionSpace(self.getDomain())  
      else:  
         raise ValueError,"unknown coefficient %s requested"%name  
427    
428     def alteredCoefficient(self,name):     def checkSymmetry(self,verbose=True):
429       """        """
430       announce that coefficient name has been changed        test the PDE for symmetry.
431    
      @param name:  
      """  
      if self.hasCoefficient(name):  
         if self.COEFFICIENTS[name].isAlteringOperator(): self.__rebuildOperator()  
         if self.COEFFICIENTS[name].isAlteringRightHandSide(): self.__rebuildRightHandSide()  
      else:  
         raise ValueError,"unknown coefficient %s requested"%name  
432    
433     # ===== debug ==============================================================       @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.
434     def setDebugOn(self):       @type verbose: C{bool}
435         """       @return:  True if the PDE is symmetric.
436         """       @rtype : C{escript.Data}
        self.__debug=not None  
437    
438     def setDebugOff(self):        @note: This is a very expensive operation. It should be used for degugging only! The symmetry flag is not altered.
439         """        """
440         """        verbose=verbose or self.debug()
441         self.__debug=None        out=True
442          if self.getNumSolutions()!=self.getNumEquations():
443             if verbose: print "non-symmetric PDE because of different number of equations and solutions"
444             out=False
445          else:
446             A=self.getCoefficientOfGeneralPDE("A")
447             if not A.isEmpty():
448                tol=util.Lsup(A)*self.TOL
449                if self.getNumSolutions()>1:
450                   for i in range(self.getNumEquations()):
451                      for j in range(self.getDim()):
452                         for k in range(self.getNumSolutions()):
453                            for l in range(self.getDim()):
454                                if util.Lsup(A[i,j,k,l]-A[k,l,i,j])>tol:
455                                   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)
456                                   out=False
457                else:
458                   for j in range(self.getDim()):
459                      for l in range(self.getDim()):
460                         if util.Lsup(A[j,l]-A[l,j])>tol:
461                            if verbose: print "non-symmetric PDE because A[%d,%d]!=A[%d,%d]"%(j,l,l,j)
462                            out=False
463             B=self.getCoefficientOfGeneralPDE("B")
464             C=self.getCoefficientOfGeneralPDE("C")
465             if B.isEmpty() and not C.isEmpty():
466                if verbose: print "non-symmetric PDE because B is not present but C is"
467                out=False
468             elif not B.isEmpty() and C.isEmpty():
469                if verbose: print "non-symmetric PDE because C is not present but B is"
470                out=False
471             elif not B.isEmpty() and not C.isEmpty():
472                tol=(util.Lsup(B)+util.Lsup(C))*self.TOL/2.
473                if self.getNumSolutions()>1:
474                   for i in range(self.getNumEquations()):
475                       for j in range(self.getDim()):
476                          for k in range(self.getNumSolutions()):
477                             if util.Lsup(B[i,j,k]-C[k,i,j])>tol:
478                                  if verbose: print "non-symmetric PDE because B[%d,%d,%d]!=C[%d,%d,%d]"%(i,j,k,k,i,j)
479                                  out=False
480                else:
481                   for j in range(self.getDim()):
482                      if util.Lsup(B[j]-C[j])>tol:
483                         if verbose: print "non-symmetric PDE because B[%d]!=C[%d]"%(j,j)
484                         out=False
485             if self.getNumSolutions()>1:
486               D=self.getCoefficientOfGeneralPDE("D")
487               if not D.isEmpty():
488                 tol=util.Lsup(D)*self.TOL
489                 for i in range(self.getNumEquations()):
490                    for k in range(self.getNumSolutions()):
491                      if util.Lsup(D[i,k]-D[k,i])>tol:
492                          if verbose: print "non-symmetric PDE because D[%d,%d]!=D[%d,%d]"%(i,k,k,i)
493                          out=False
494    
495     def debug(self):        return out
496    
497       def getSolution(self,**options):
498         """         """
499         returns true if the PDE is in the debug mode         returns the solution of the PDE. If the solution is not valid the PDE is solved.
500    
501           @return: the solution
502           @rtype: L{escript.Data}
503           @param options: solver options
504           @keyword verbose:
505           @keyword reordering: reordering scheme to be used during elimination
506           @keyword preconditioner: preconditioner method to be used
507           @keyword iter_max: maximum number of iteration steps allowed.
508           @keyword drop_tolerance:
509           @keyword drop_storage:
510           @keyword truncation:
511           @keyword restart:
512         """         """
513         return self.__debug         if not self.__solution_isValid:
514              mat,f=self.getSystem()
515              if self.isUsingLumping():
516                 self.__solution=self.copyConstraint(f*mat)
517              else:
518                 options[self.TOLERANCE_KEY]=self.getTolerance()
519                 options[self.METHOD_KEY]=self.getSolverMethod()
520                 options[self.SYMMETRY_KEY]=self.isSymmetric()
521                 self.trace("PDE is resolved.")
522                 self.trace("solver options: %s"%str(options))
523                 self.__solution=mat.solve(f,options)
524              self.__solution_isValid=True
525           return self.__solution
526    
527     #===== Lumping ===========================     def getFlux(self,u=None):
528     def setLumpingOn(self):       """
529        """       returns the flux J_ij for a given u
       indicates to use matrix lumping  
       """  
       if not self.isUsingLumping():  
          if self.debug() : print "PDE Debug: lumping is set on"  
          self.__rebuildOperator()  
          self.__lumping=True  
530    
531     def setLumpingOff(self):         \f[
532        """         J_ij=A_{ijkl}u_{k,l}+B_{ijk}u_k-X_{ij}
533        switches off matrix lumping         \f]
       """  
       if self.isUsingLumping():  
          if self.debug() : print "PDE Debug: lumping is set off"  
          self.__rebuildOperator()  
          self.__lumping=False  
534    
535     def setLumping(self,flag=False):       @param u: argument in the flux. If u is not present or equals L{None} the current solution is used.
536        """       @type u: L{escript.Data} or None
537        set the matrix lumping flag to flag       @return : flux
538        """       @rtype : L{escript.Data}
       if flag:  
          self.setLumpingOn()  
       else:  
          self.setLumpingOff()  
539    
540     def isUsingLumping(self):       """
541        """       if u==None: u=self.getSolution()
542               return util.tensormult(self.getCoefficientOfGeneralPDE("A"),util.grad(u))+util.matrixmult(self.getCoefficientOfGeneralPDE("B"),u)-util.self.getCoefficientOfGeneralPDE("X")
543        """  
544        return self.__lumping     # =============================================================================
545       #   solver settings:
546     #============ method business =========================================================     # =============================================================================
547     def setSolverMethod(self,solver=util.DEFAULT_METHOD):     def setSolverMethod(self,solver=None):
548         """         """
549         sets a new solver         sets a new solver
550    
551           @param solver: sets a new solver method.
552           @type solver: C{int}
553    
554         """         """
555           if solver==None: solve=self.DEFAULT_METHOD
556         if not solver==self.getSolverMethod():         if not solver==self.getSolverMethod():
557             self.__solver_method=solver             self.__solver_method=solver
            if self.debug() : print "PDE Debug: New solver is %s"%solver  
558             self.__checkMatrixType()             self.__checkMatrixType()
559               self.trace("New solver is %s"%self.getSolverMethodName())
560    
561       def getSolverMethodName(self):
562           """
563           returns the name of the solver currently used
564    
565           @return : the name of the solver currently used.
566           @rtype: C{string}
567           """
568    
569           m=self.getSolverMethod()
570           if m==self.DEFAULT_METHOD: return "DEFAULT_METHOD"
571           elif m==self.DIRECT: return "DIRECT"
572           elif m==self.CHOLEVSKY: return "CHOLEVSKY"
573           elif m==self.PCG: return "PCG"
574           elif m==self.CR: return "CR"
575           elif m==self.CGS: return "CGS"
576           elif m==self.BICGSTAB: return "BICGSTAB"
577           elif m==self.SSOR: return "SSOR"
578           elif m==self.GMRES: return "GMRES"
579           elif m==self.PRES20: return "PRES20"
580           elif m==self.LUMPING: return "LUMPING"
581           return None
582          
583    
584     def getSolverMethod(self):     def getSolverMethod(self):
585         """         """
586         returns the solver method         returns the solver method
587      
588           @return : the solver method currently be used.
589           @rtype : C{int}
590         """         """
591         return self.__solver_method         return self.__solver_method
592    
593     #============ tolerance business =========================================================     def isUsingLumping(self):
594          """
595          checks if matrix lumping is used a solver method
596    
597          @return : True is lumping is currently used a solver method.
598          @rtype: C{bool}
599          """
600          return self.getSolverMethod()==self.LUMPING
601    
602     def setTolerance(self,tol=1.e-8):     def setTolerance(self,tol=1.e-8):
603         """         """
604         resets the tolerance to tol.         resets the tolerance for the solver method to tol where for an appropriate norm |.|
605    
606                   |self.getResidual()|<tol*|self.getRightHandSide()|
607    
608           defines the stopping criterion.
609    
610           @param tol: new tolerance for the solver. If the tol is lower then the current tolerence
611                       the system will be resolved.
612           @type solver: C{float}
613           @raise ValueException: if tolerance is not positive.
614         """         """
615         if not tol>0:         if not tol>0:
616             raise ValueException,"Tolerance as to be positive"             raise ValueException,"Tolerance as to be positive"
617         if tol<self.getTolerance(): self.__rebuildSolution()         if tol<self.getTolerance(): self.__invalidateSolution()
618         if self.debug() : print "PDE Debug: New tolerance %e",tol         self.trace("New tolerance %e"%tol)
619         self.__tolerance=tol         self.__tolerance=tol
620         return         return
621    
622     def getTolerance(self):     def getTolerance(self):
623         """         """
624         returns the tolerance set for the solution         returns the tolerance set for the solution
625    
626           @return: tolerance currently used.
627           @rtype: C{float}
628         """         """
629         return self.__tolerance         return self.__tolerance
630    
631     #===== symmetry  flag ==========================     # =============================================================================
632       #    symmetry  flag:
633       # =============================================================================
634     def isSymmetric(self):     def isSymmetric(self):
635        """        """
636        returns true is the operator is considered to be symmetric        checks if symmetry is indicated.
637        
638          @return : True is a symmetric PDE is indicated, otherwise False is returned
639          @rtype : C{bool}
640        """        """
641        return self.__sym        return self.__sym
642    
643     def setSymmetryOn(self):     def setSymmetryOn(self):
644        """        """
645        sets the symmetry flag to true        sets the symmetry flag.
646        """        """
647        if not self.isSymmetric():        if not self.isSymmetric():
648           if self.debug() : print "PDE Debug: Operator is set to be symmetric"           self.trace("PDE is set to be symmetric")
649           self.__sym=True           self.__sym=True
650           self.__checkMatrixType()           self.__checkMatrixType()
651    
652     def setSymmetryOff(self):     def setSymmetryOff(self):
653        """        """
654        sets the symmetry flag to false        removes the symmetry flag.
655        """        """
656        if self.isSymmetric():        if self.isSymmetric():
657           if self.debug() : print "PDE Debug: Operator is set to be unsymmetric"           self.trace("PDE is set to be unsymmetric")
658           self.__sym=False           self.__sym=False
659           self.__checkMatrixType()           self.__checkMatrixType()
660    
661     def setSymmetryTo(self,flag=False):     def setSymmetryTo(self,flag=False):
662       """        """
663       sets the symmetry flag to flag        sets the symmetry flag to flag
664    
665       @param flag:        @param flag: If flag, the symmetry flag is set otherwise the symmetry flag is released.
666       """        @type flag: C{bool}
667       if flag:        """
668          self.setSymmetryOn()        if flag:
669       else:           self.setSymmetryOn()
670          self.setSymmetryOff()        else:
671             self.setSymmetryOff()
672    
673     #===== order reduction ==========================    
674       # =============================================================================
675       # function space handling for the equation as well as the solution
676       # =============================================================================
677     def setReducedOrderOn(self):     def setReducedOrderOn(self):
678       """       """
679       switches to on reduced order       switches on reduced order for solution and equation representation
680       """       """
681       self.setReducedOrderForSolutionOn()       self.setReducedOrderForSolutionOn()
682       self.setReducedOrderForEquationOn()       self.setReducedOrderForEquationOn()
683    
684     def setReducedOrderOff(self):     def setReducedOrderOff(self):
685       """       """
686       switches to full order       switches off reduced order for solution and equation representation
687       """       """
688       self.setReducedOrderForSolutionOff()       self.setReducedOrderForSolutionOff()
689       self.setReducedOrderForEquationOff()       self.setReducedOrderForEquationOff()
690    
691     def setReducedOrderTo(self,flag=False):     def setReducedOrderTo(self,flag=False):
692       """       """
693       sets order according to flag       sets order reduction for both solution and equation representation according to flag.
694    
695       @param flag:       @param flag: if flag is True, the order reduction is switched on for both  solution and equation representation, otherwise or
696                      if flag is not present order reduction is switched off
697         @type flag: C{bool}
698       """       """
699       self.setReducedOrderForSolutionTo(flag)       self.setReducedOrderForSolutionTo(flag)
700       self.setReducedOrderForEquationTo(flag)       self.setReducedOrderForEquationTo(flag)
                                                                                                                                                             
701    
702     #===== order reduction solution ==========================  
703     def setReducedOrderForSolutionOn(self):     def setReducedOrderForSolutionOn(self):
704       """       """
705       switches to reduced order to interpolate solution       switches on reduced order for solution representation
706       """       """
707       new_fs=escript.ReducedSolution(self.getDomain())       new_fs=escript.ReducedSolution(self.getDomain())
708       if self.getFunctionSpaceForSolution()!=new_fs:       if self.getFunctionSpaceForSolution()!=new_fs:
709           if self.debug() : print "PDE Debug: Reduced order is used to interpolate solution."           self.trace("Reduced order is used to solution representation.")
710           self.__column_function_space=new_fs           self.__column_function_space=new_fs
711           self.__rebuildSystem(deep=True)           self.__resetSystem()
712    
713     def setReducedOrderForSolutionOff(self):     def setReducedOrderForSolutionOff(self):
714       """       """
715       switches to full order to interpolate solution       switches off reduced order for solution representation
716       """       """
717       new_fs=escript.Solution(self.getDomain())       new_fs=escript.Solution(self.getDomain())
718       if self.getFunctionSpaceForSolution()!=new_fs:       if self.getFunctionSpaceForSolution()!=new_fs:
719           if self.debug() : print "PDE Debug: Full order is used to interpolate solution."           self.trace("Full order is used to interpolate solution.")
720           self.__column_function_space=new_fs           self.__column_function_space=new_fs
721           self.__rebuildSystem(deep=True)           self.__resetSystem()
722    
723     def setReducedOrderForSolutionTo(self,flag=False):     def setReducedOrderForSolutionTo(self,flag=False):
724       """       """
725       sets order for test functions according to flag       sets order for test functions according to flag
726    
727       @param flag:       @param flag: if flag is True, the order reduction is switched on for solution representation, otherwise or
728                      if flag is not present order reduction is switched off
729         @type flag: C{bool}
730       """       """
731       if flag:       if flag:
732          self.setReducedOrderForSolutionOn()          self.setReducedOrderForSolutionOn()
733       else:       else:
734          self.setReducedOrderForSolutionOff()          self.setReducedOrderForSolutionOff()
735                                                                                                                                                              
    #===== order reduction equation ==========================  
736     def setReducedOrderForEquationOn(self):     def setReducedOrderForEquationOn(self):
737       """       """
738       switches to reduced order for test functions       switches on reduced order for equation representation
739       """       """
740       new_fs=escript.ReducedSolution(self.getDomain())       new_fs=escript.ReducedSolution(self.getDomain())
741       if self.getFunctionSpaceForEquation()!=new_fs:       if self.getFunctionSpaceForEquation()!=new_fs:
742           if self.debug() : print "PDE Debug: Reduced order is used for test functions."           self.trace("Reduced order is used for test functions.")
743           self.__row_function_space=new_fs           self.__row_function_space=new_fs
744           self.__rebuildSystem(deep=True)           self.__resetSystem()
745    
746     def setReducedOrderForEquationOff(self):     def setReducedOrderForEquationOff(self):
747       """       """
748       switches to full order for test functions       switches off reduced order for equation representation
749       """       """
750       new_fs=escript.Solution(self.getDomain())       new_fs=escript.Solution(self.getDomain())
751       if self.getFunctionSpaceForEquation()!=new_fs:       if self.getFunctionSpaceForEquation()!=new_fs:
752           if self.debug() : print "PDE Debug: Full order is used for test functions."           self.trace("Full order is used for test functions.")
753           self.__row_function_space=new_fs           self.__row_function_space=new_fs
754           self.__rebuildSystem(deep=True)           self.__resetSystem()
755    
756     def setReducedOrderForEquationTo(self,flag=False):     def setReducedOrderForEquationTo(self,flag=False):
757       """       """
758       sets order for test functions according to flag       sets order for test functions according to flag
759    
760       @param flag:       @param flag: if flag is True, the order reduction is switched on for equation representation, otherwise or
761                      if flag is not present order reduction is switched off
762         @type flag: C{bool}
763       """       """
764       if flag:       if flag:
765          self.setReducedOrderForEquationOn()          self.setReducedOrderForEquationOn()
766       else:       else:
767          self.setReducedOrderForEquationOff()          self.setReducedOrderForEquationOff()
768                                                                                                                                                              
769     # ==== initialization =====================================================================     # =============================================================================
770       # private method:
771       # =============================================================================
772       def __checkMatrixType(self):
773         """
774         reassess the matrix type and, if a new matrix is needed, resets the system.
775         """
776         new_matrix_type=self.getDomain().getSystemMatrixTypeId(self.getSolverMethod(),self.isSymmetric())
777         if not new_matrix_type==self.__matrix_type:
778             self.trace("Matrix type is now %d."%new_matrix_type)
779             self.__matrix_type=new_matrix_type
780             self.__resetSystem()
781       #
782       #   rebuild switches :
783       #
784       def __invalidateSolution(self):
785           """
786           indicates the PDE has to be resolved if the solution is requested
787           """
788           if self.__solution_isValid: self.trace("PDE has to be resolved.")
789           self.__solution_isValid=False
790    
791       def __invalidateOperator(self):
792           """
793           indicates the operator has to be rebuilt next time it is used
794           """
795           if self.__operator_isValid: self.trace("Operator has to be rebuilt.")
796           self.__invalidateSolution()
797           self.__operator_isValid=False
798    
799       def __invalidateRightHandSide(self):
800           """
801           indicates the right hand side has to be rebuild next time it is used
802           """
803           if self.__righthandside_isValid: self.trace("Right hand side has to be rebuilt.")
804           self.__invalidateSolution()
805           self.__righthandside_isValid=False
806    
807       def __invalidateSystem(self):
808           """
809           annonced that everthing has to be rebuild:
810           """
811           if self.__righthandside_isValid: self.trace("System has to be rebuilt.")
812           self.__invalidateSolution()
813           self.__invalidateOperator()
814           self.__invalidateRightHandSide()
815    
816       def __resetSystem(self):
817           """
818           annonced that everthing has to be rebuild:
819           """
820           self.trace("New System is built from scratch.")
821           self.__operator=escript.Operator()
822           self.__operator_isValid=False
823           self.__righthandside=escript.Data()
824           self.__righthandside_isValid=False
825           self.__solution=escript.Data()
826           self.__solution_isValid=False
827       #
828       #    system initialization:
829       #
830     def __getNewOperator(self):     def __getNewOperator(self):
831         """         """
832           returns an instance of a new operator
833         """         """
834           self.trace("New operator is allocated.")
835         return self.getDomain().newOperator( \         return self.getDomain().newOperator( \
836                             self.getNumEquations(), \                             self.getNumEquations(), \
837                             self.getFunctionSpaceForEquation(), \                             self.getFunctionSpaceForEquation(), \
# Line 601  class LinearPDE: Line 839  class LinearPDE:
839                             self.getFunctionSpaceForSolution(), \                             self.getFunctionSpaceForSolution(), \
840                             self.__matrix_type)                             self.__matrix_type)
841    
842     def __makeFreshRightHandSide(self):     def __getNewRightHandSide(self):
843         """         """
844           returns an instance of a new right hand side
845         """         """
846         if self.debug() : print "PDE Debug: New right hand side allocated"         self.trace("New right hand side is allocated.")
847         if self.getNumEquations()>1:         if self.getNumEquations()>1:
848             self.__righthandside=escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)             return escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)
849         else:         else:
850             self.__righthandside=escript.Data(0.,(),self.getFunctionSpaceForEquation(),True)             return escript.Data(0.,(),self.getFunctionSpaceForEquation(),True)
        return self.__righthandside  
851    
852     def __getNewSolution(self):     def __getNewSolution(self):
853         """         """
854           returns an instance of a new solution
855         """         """
856         if self.debug() : print "PDE Debug: New right hand side allocated"         self.trace("New solution is allocated.")
857         if self.getNumSolutions()>1:         if self.getNumSolutions()>1:
858             return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)             return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)
859         else:         else:
860             return escript.Data(0.,(),self.getFunctionSpaceForSolution(),True)             return escript.Data(0.,(),self.getFunctionSpaceForSolution(),True)
861    
862       def __makeFreshSolution(self):
863           """
864           makes sure that the solution is instantiated and returns it initialized by zeros
865           """
866           if self.__solution.isEmpty():
867               self.__solution=self.__getNewSolution()
868           else:
869               self.__solution*=0
870               self.trace("Solution is reset to zero.")
871           return self.__solution
872    
873       def __makeFreshRightHandSide(self):
874           """
875           makes sure that the right hand side is instantiated and returns it initialized by zeros
876           """
877           if self.__righthandside.isEmpty():
878               self.__righthandside=self.__getNewRightHandSide()
879           else:
880               self.__righthandside*=0
881               self.trace("Right hand side is reset to zero.")
882           return self.__righthandside
883    
884     def __makeFreshOperator(self):     def __makeFreshOperator(self):
885         """         """
886           makes sure that the operator is instantiated and returns it initialized by zeros
887         """         """
888         if self.__operator.isEmpty():         if self.__operator.isEmpty():
889             self.__operator=self.__getNewOperator()             self.__operator=self.__getNewOperator()
            if self.debug() : print "PDE Debug: New operator allocated"  
890         else:         else:
891             self.__operator.setValue(0.)             self.__operator.setValue(0.)
892             self.__operator.resetSolver()             self.trace("Operator reset to zero")
            if self.debug() : print "PDE Debug: Operator reset to zero"  
893         return self.__operator         return self.__operator
894    
895     #============ some serivice functions  =====================================================     def __applyConstraint(self):
896     def getDomain(self):         """
897           applies the constraints defined by q and r to the system
898           """
899           if not self.isUsingLumping():
900              q=self.getCoefficientOfGeneralPDE("q")
901              r=self.getCoefficientOfGeneralPDE("r")
902              if not q.isEmpty() and not self.__operator.isEmpty():
903                 # q is the row and column mask to indicate where constraints are set:
904                 row_q=escript.Data(q,self.getFunctionSpaceForEquation())
905                 col_q=escript.Data(q,self.getFunctionSpaceForSolution())
906                 u=self.__getNewSolution()
907                 if r.isEmpty():
908                    r_s=self.__getNewSolution()
909                 else:
910                    r_s=escript.Data(r,self.getFunctionSpaceForSolution())
911                 u.copyWithMask(r_s,col_q)
912                 if not self.__righthandside.isEmpty():
913                    self.__righthandside-=self.__operator*u
914                    self.__righthandside=self.copyConstraint(self.__righthandside)
915                 self.__operator.nullifyRowsAndCols(row_q,col_q,1.)
916       # =============================================================================
917       # function giving access to coefficients of the general PDE:
918       # =============================================================================
919       def getCoefficientOfGeneralPDE(self,name):
920         """
921         return the value of the coefficient name of the general PDE.
922    
923         @note This method is called by the assembling routine it can be overwritten
924               to map coefficients of a particular PDE to the general PDE.
925    
926         @param name: name of the coefficient requested.
927         @type name: C{string}
928         @return : the value of the coefficient  name
929         @rtype : L{escript.Data}
930         @raise IllegalCoefficient: if name is not one of coefficients
931                      "A", "B", "C", "D", "X", "Y", "d", "y", "d_contact", "y_contact", "r" or "q".
932       """       """
933       returns the domain of the PDE       if self.hasCoefficientOfGeneralPDE(name):
934            return self.getCoefficient(name)
935         else:
936            raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
937    
938       def hasCoefficientOfGeneralPDE(self,name):
939       """       """
940       return self.__domain       checks if name is a the name of a coefficient of the general PDE.
941        
942         @param name: name of the coefficient enquired.
943         @type name: C{string}
944         @return : True if name is the name of a coefficient of the general PDE. Otherwise False.
945         @rtype : C{bool}
946        
947         """
948         return self.__COEFFICIENTS_OF_GENEARL_PDE.has_key(name)
949    
950     def getDim(self):     def createCoefficientOfGeneralPDE(self,name):
951       """       """
952       returns the spatial dimension of the PDE       returns a new instance of a coefficient for coefficient name of the general PDE
953    
954         @param name: name of the coefficient requested.
955         @type name: C{string}
956         @return : a coefficient name initialized to 0.
957         @rtype : L{escript.Data}
958         @raise IllegalCoefficient: if name is not one of coefficients
959                      "A", "B", "C", "D", "X", "Y", "d", "y", "d_contact", "y_contact", "r" or "q".
960       """       """
961       return self.getDomain().getDim()       if self.hasCoefficientOfGeneralPDE(name):
962            return escript.Data(0,self.getShapeOfCoefficientOfGeneralPDE(name),self.getFunctionSpaceForCoefficientOfGeneralPDE(name))
963         else:
964            raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
965    
966     def getNumEquations(self):     def getFunctionSpaceForCoefficientOfGeneralPDE(self,name):
967       """       """
968       returns the number of equations       return the L{escript.FunctionSpace} to be used for coefficient name of the general PDE
969    
970         @param name: name of the coefficient enquired.
971         @type name: C{string}
972         @return : the function space to be used for coefficient name
973         @rtype : L{escript.FunctionSpace}
974         @raise IllegalCoefficient: if name is not one of coefficients
975                      "A", "B", "C", "D", "X", "Y", "d", "y", "d_contact", "y_contact", "r" or "q".
976       """       """
977       if self.__numEquations>0:       if self.hasCoefficientOfGeneralPDE(name):
978           return self.__numEquations          return self.__COEFFICIENTS_OF_GENEARL_PDE[name].getFunctionSpace(self.getDomain())
979       else:       else:
980           raise ValueError,"Number of equations is undefined. Please specify argument numEquations."          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
981    
982     def getNumSolutions(self):     def getShapeOfCoefficientOfGeneralPDE(self,name):
983       """       """
984       returns the number of unknowns       return the shape of the coefficient name of the general PDE
985    
986         @param name: name of the coefficient enquired.
987         @type name: C{string}
988         @return : the shape of the coefficient name
989         @rtype : C{tuple} of C{int}
990         @raise IllegalCoefficient: if name is not one of coefficients
991                      "A", "B", "C", "D", "X", "Y", "d", "y", "d_contact", "y_contact", "r" or "q".
992    
993       """       """
994       if self.__numSolutions>0:       if self.hasCoefficientOfGeneralPDE(name):
995          return self.__numSolutions          return self.__COEFFICIENTS_OF_GENEARL_PDE[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())
996       else:       else:
997          raise ValueError,"Number of solution is undefined. Please specify argument numSolutions."          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
998    
999       # =============================================================================
1000       # functions giving access to coefficients of a particular PDE implementation:
1001       # =============================================================================
1002       def getCoefficient(self,name):
1003         """
1004         returns the value of the coefficient name
1005    
1006     def checkSymmetry(self,verbose=True):       @param name: name of the coefficient requested.
1007        """       @type name: C{string}
1008        returns if the Operator is symmetric. This is a very expensive operation!!! The symmetry flag is not altered.       @return : the value of the coefficient name
1009        """       @rtype : L{escript.Data}
1010        verbose=verbose or self.debug()       @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1011        out=True       """
1012        if self.getNumSolutions()!=self.getNumEquations():       if self.hasCoefficient(name):
1013           if verbose: print "non-symmetric PDE because of different number of equations and solutions"           return self.COEFFICIENTS[name].getValue()
1014           out=False       else:
1015        else:          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
          A=self.getCoefficientOfPDE("A")  
          if not A.isEmpty():  
             tol=util.Lsup(A)*self.TOL  
             if self.getNumSolutions()>1:  
                for i in range(self.getNumEquations()):  
                   for j in range(self.getDim()):  
                      for k in range(self.getNumSolutions()):  
                         for l in range(self.getDim()):  
                             if util.Lsup(A[i,j,k,l]-A[k,l,i,j])>tol:  
                                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)  
                                out=False  
             else:  
                for j in range(self.getDim()):  
                   for l in range(self.getDim()):  
                      if util.Lsup(A[j,l]-A[l,j])>tol:  
                         if verbose: print "non-symmetric PDE because A[%d,%d]!=A[%d,%d]"%(j,l,l,j)  
                         out=False  
          B=self.getCoefficientOfPDE("B")  
          C=self.getCoefficientOfPDE("C")  
          if B.isEmpty() and not C.isEmpty():  
             if verbose: print "non-symmetric PDE because B is not present but C is"  
             out=False  
          elif not B.isEmpty() and C.isEmpty():  
             if verbose: print "non-symmetric PDE because C is not present but B is"  
             out=False  
          elif not B.isEmpty() and not C.isEmpty():  
             tol=(util.Lsup(B)+util.Lsup(C))*self.TOL/2.  
             if self.getNumSolutions()>1:  
                for i in range(self.getNumEquations()):  
                    for j in range(self.getDim()):  
                       for k in range(self.getNumSolutions()):  
                          if util.Lsup(B[i,j,k]-C[k,i,j])>tol:  
                               if verbose: print "non-symmetric PDE because B[%d,%d,%d]!=C[%d,%d,%d]"%(i,j,k,k,i,j)  
                               out=False  
             else:  
                for j in range(self.getDim()):  
                   if util.Lsup(B[j]-C[j])>tol:  
                      if verbose: print "non-symmetric PDE because B[%d]!=C[%d]"%(j,j)  
                      out=False  
          if self.getNumSolutions()>1:  
            D=self.getCoefficientOfPDE("D")  
            if not D.isEmpty():  
              tol=util.Lsup(D)*self.TOL  
              for i in range(self.getNumEquations()):  
                 for k in range(self.getNumSolutions()):  
                   if util.Lsup(D[i,k]-D[k,i])>tol:  
                       if verbose: print "non-symmetric PDE because D[%d,%d]!=D[%d,%d]"%(i,k,k,i)  
                       out=False  
               
       return out  
1016    
1017     def getFlux(self,u):     def hasCoefficient(self,name):
1018         """       """
1019         returns the flux J_ij for a given u       return True if name is the name of a coefficient
1020    
1021         \f[       @param name: name of the coefficient enquired.
1022         J_ij=A_{ijkl}u_{k,l}+B_{ijk}u_k-X_{ij}       @type name: C{string}
1023         \f]       @return : True if name is the name of a coefficient of the general PDE. Otherwise False.
1024         @rtype : C{bool}
1025         """
1026         return self.COEFFICIENTS.has_key(name)
1027    
1028         @param u: argument of the operator     def createCoefficient(self, name):
1029         """       """
1030         raise SystemError,"getFlux is not implemented yet"       create a L{escript.Data} object corresponding to coefficient name
        return None  
1031    
1032     def applyOperator(self,u):       @return : a coefficient name initialized to 0.
1033         """       @rtype : L{escript.Data}
1034         applies the operator of the PDE to a given solution u in weak from       @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1035         """
1036         if self.hasCoefficient(name):
1037            return escript.Data(0.,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))
1038         else:
1039            raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1040    
1041         @param u: argument of the operator     def getFunctionSpaceForCoefficient(self,name):
1042         """       """
1043         return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())       return the L{escript.FunctionSpace} to be used for coefficient name
                                                                                                                                                             
    def getResidual(self,u):  
        """  
        return the residual of u in the weak from  
1044    
1045         @param u:       @param name: name of the coefficient enquired.
1046         """       @type name: C{string}
1047         return self.applyOperator(u)-self.getRightHandSide()       @return : the function space to be used for coefficient name
1048         @rtype : L{escript.FunctionSpace}
1049         @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1050         """
1051         if self.hasCoefficient(name):
1052            return self.COEFFICIENTS[name].getFunctionSpace(self.getDomain())
1053         else:
1054            raise ValueError,"unknown coefficient %s requested"%name
1055    
1056       def getShapeOfCoefficient(self,name):
1057         """
1058         return the shape of the coefficient name
1059    
1060         @param name: name of the coefficient enquired.
1061         @type name: C{string}
1062         @return : the shape of the coefficient name
1063         @rtype : C{tuple} of C{int}
1064         @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1065         """
1066         if self.hasCoefficient(name):
1067            return self.COEFFICIENTS[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())
1068         else:
1069            raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1070    
1071       def resetCoefficients(self):
1072         """
1073         resets all coefficients to there default values.
1074         """
1075         for i in self.COEFFICIENTS.iterkeys():
1076             self.COEFFICIENTS[i].resetValue()
1077    
1078     def __setValue(self,**coefficients):     def alteredCoefficient(self,name):
1079         """
1080         announce that coefficient name has been changed
1081    
1082         @param name: name of the coefficient enquired.
1083         @type name: C{string}
1084         @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1085         """
1086         if self.hasCoefficient(name):
1087            self.trace("Coefficient %s has been altered."%name)
1088            if self.COEFFICIENTS[name].isAlteringOperator(): self.__invalidateOperator()
1089            if self.COEFFICIENTS[name].isAlteringRightHandSide(): self.__invalidateRightHandSide()
1090         else:
1091            raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1092    
1093       def copyConstraint(self,u):
1094        """        """
1095        sets new values to coefficient        copies the constraint into u and returns u.
1096    
1097          @param u: a function of rank 0 is a single PDE is solved and of shape (numSolution,) for a system of PDEs
1098          @type u: L{escript.Data}
1099          @return : the input u modified by the constraints.
1100          @rtype : L{escript.Data}
1101          @warning: u is altered if it has the appropriate L{escript.FunctionSpace}
1102    
1103          """
1104          q=self.getCoefficientOfGeneralPDE("q")
1105          r=self.getCoefficientOfGeneralPDE("r")
1106          if not q.isEmpty():
1107             if u.isEmpty(): u=escript.Data(0.,q.getShape(),q.getFunctionSpace())
1108             if r.isEmpty():
1109                 r=escript.Data(0,u.getShape(),u.getFunctionSpace())
1110             else:
1111                 r=escript.Data(r,u.getFunctionSpace())
1112             u.copyWithMask(r,escript.Data(q,u.getFunctionSpace()))
1113          return u
1114    
1115       def setValue(self,**coefficients):
1116          """
1117          sets new values to coefficients
1118    
1119          @note This method is called by the assembling routine it can be overwritten
1120               to map coefficients of a particular PDE to the general PDE.
1121    
1122          @param name: name of the coefficient requested.
1123          @type name: C{string}
1124          @keyword A: value for coefficient A.
1125          @type A: any type that can be interpreted as L{escript.Data} object on L{escript.Function}.
1126          @keyword B: value for coefficient B
1127          @type B: any type that can be interpreted as L{escript.Data} object on L{escript.Function}.
1128          @keyword C: value for coefficient C
1129          @type C: any type that can be interpreted as L{escript.Data} object on L{escript.Function}.
1130          @keyword D: value for coefficient D
1131          @type D: any type that can be interpreted as L{escript.Data} object on L{escript.Function}.
1132          @keyword X: value for coefficient X
1133          @type X: any type that can be interpreted as L{escript.Data} object on L{escript.Function}.
1134          @keyword Y: value for coefficient Y
1135          @type Y: any type that can be interpreted as L{escript.Data} object on L{escript.Function}.
1136          @keyword d: value for coefficient d
1137          @type d: any type that can be interpreted as L{escript.Data} object on L{escript.FunctionOnBoundary}.
1138          @keyword y: value for coefficient y
1139          @type y: any type that can be interpreted as L{escript.Data} object on L{escript.FunctionOnBoundary}.
1140          @keyword d_contact: value for coefficient d_contact
1141          @type d_contact: any type that can be interpreted as L{escript.Data} object on L{escript.FunctionOnContactOne}.
1142                           or  L{escript.FunctionOnContactZero}.
1143          @keyword y_contact: value for coefficient y_contact
1144          @type y_contact: any type that can be interpreted as L{escript.Data} object on L{escript.FunctionOnContactOne}.
1145                           or  L{escript.FunctionOnContactZero}.
1146          @keyword r: values prescribed to the solution at the locations of constraints
1147          @type r: any type that can be interpreted as L{escript.Data} object on L{escript.Solution} or L{escript.ReducedSolution}
1148                   depending of reduced order is used for the solution.
1149          @keyword q: mask for location of constraints
1150          @type q: any type that can be interpreted as L{escript.Data} object on L{escript.Solution} or L{escript.ReducedSolution}
1151                   depending of reduced order is used for the representation of the equation.
1152          @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1153    
       @param coefficients:  
1154        """        """
1155        # check if the coefficients are  legal:        # check if the coefficients are  legal:
1156        for i in coefficients.iterkeys():        for i in coefficients.iterkeys():
1157           if not self.hasCoefficient(i):           if not self.hasCoefficient(i):
1158              raise ValueError,"Attempt to set unknown coefficient %s"%i              raise IllegalCoefficient,"Attempt to set unknown coefficient %s"%i
1159        # if the number of unknowns or equations is still unknown we try to estimate them:        # if the number of unknowns or equations is still unknown we try to estimate them:
1160        if self.__numEquations<1 or self.__numSolutions<1:        if self.__numEquations==None or self.__numSolutions==None:
1161           for i,d in coefficients.iteritems():           for i,d in coefficients.iteritems():
1162              if hasattr(d,"shape"):              if hasattr(d,"shape"):
1163                  s=d.shape                  s=d.shape
# Line 775  class LinearPDE: Line 1167  class LinearPDE:
1167                  s=numarray.array(d).shape                  s=numarray.array(d).shape
1168              if s!=None:              if s!=None:
1169                  # get number of equations and number of unknowns:                  # get number of equations and number of unknowns:
1170                  res=self.COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(s,self.getDim())                  res=self.COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(self.getDomain(),s)
1171                  if res==None:                  if res==None:
1172                      raise ValueError,"Illegal shape %s of coefficient %s"%(s,i)                      raise IllegalCoefficientValue,"Illegal shape %s of coefficient %s"%(s,i)
1173                  else:                  else:
1174                      if self.__numEquations<1: self.__numEquations=res[0]                      if self.__numEquations==None: self.__numEquations=res[0]
1175                      if self.__numSolutions<1: self.__numSolutions=res[1]                      if self.__numSolutions==None: self.__numSolutions=res[1]
1176        if self.__numEquations<1: raise ValueError,"unidententified number of equations"        if self.__numEquations==None: raise UndefinedPDEError,"unidententified number of equations"
1177        if self.__numSolutions<1: raise ValueError,"unidententified number of solutions"        if self.__numSolutions==None: raise UndefinedPDEError,"unidententified number of solutions"
1178        # now we check the shape of the coefficient if numEquations and numSolutions are set:        # now we check the shape of the coefficient if numEquations and numSolutions are set:
1179        for i,d in coefficients.iteritems():        for i,d in coefficients.iteritems():
1180          if d==None:          try:
1181               d2=escript.Data()             self.COEFFICIENTS[i].setValue(self.getDomain(),self.getNumEquations(),self.getNumSolutions(),d)
1182          elif isinstance(d,escript.Data):          except IllegalCoefficientValue,m:
1183               if d.isEmpty():             raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
                 d2=d  
              else:  
                 d2=escript.Data(d,self.getFunctionSpaceForCoefficient(i))  
         else:  
               d2=escript.Data(d,self.getFunctionSpaceForCoefficient(i))  
         if not d2.isEmpty():  
            if not self.getShapeOfCoefficient(i)==d2.getShape():  
                raise ValueError,"Expected shape for coefficient %s is %s but actual shape is %s."%(i,self.getShapeOfCoefficient(i),d2.getShape())  
         # overwrite new values:  
         if self.debug(): print "PDE Debug: Coefficient %s has been altered."%i  
         self.COEFFICIENTS[i].setValue(d2)  
1184          self.alteredCoefficient(i)          self.alteredCoefficient(i)
         
       # reset the HomogeneousConstraintFlag:  
       self.__setHomogeneousConstraintFlag()  
       if len(coefficients)>0 and not self.isUsingLumping() and not self.__homogeneous_constraint: self.__rebuildSystem()  
   
    def __setHomogeneousConstraintFlag(self):  
       """  
       checks if the constraints are homogeneous and sets self.__homogeneous_constraint accordingly.  
       """  
       self.__homogeneous_constraint=True  
       q=self.getCoefficientOfPDE("q")  
       r=self.getCoefficientOfPDE("r")  
       if not q.isEmpty() and not r.isEmpty():  
          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."  
   
   
    # ==== rebuild switches =====================================================================  
    def __rebuildSolution(self,deep=False):  
        """  
        indicates the PDE has to be reolved if the solution is requested  
        """  
        if self.__solution_isValid and self.debug() : print "PDE Debug: PDE has to be resolved."  
        self.__solution_isValid=False  
        if deep: self.__solution=escript.Data()  
1185    
1186          # check if the systrem is inhomogeneous:
1187          if len(coefficients)>0 and not self.isUsingLumping():
1188             q=self.getCoefficientOfGeneralPDE("q")
1189             r=self.getCoefficientOfGeneralPDE("r")
1190             homogeneous_constraint=True
1191             if not q.isEmpty() and not r.isEmpty():
1192                 if util.Lsup(q*r)>=1.e-13*util.Lsup(r):
1193                   self.trace("Inhomogeneous constraint detected.")
1194                   self.__invalidateSystem()
1195    
    def __rebuildOperator(self,deep=False):  
        """  
        indicates the operator has to be rebuilt next time it is used  
        """  
        if self.__operator_isValid and self.debug() : print "PDE Debug: Operator has to be rebuilt."  
        self.__rebuildSolution(deep)  
        self.__operator_isValid=False  
        if deep: self.__operator=escript.Operator()  
   
    def __rebuildRightHandSide(self,deep=False):  
        """  
        indicates the right hand side has to be rebuild next time it is used  
        """  
        if self.__righthandside_isValid and self.debug() : print "PDE Debug: Right hand side has to be rebuilt."  
        self.__rebuildSolution(deep)  
        self.__righthandside_isValid=False  
        if deep: self.__righthandside=escript.Data()  
   
    def __rebuildSystem(self,deep=False):  
        """  
        annonced that all coefficient name has been changed  
        """  
        self.__rebuildSolution(deep)  
        self.__rebuildOperator(deep)  
        self.__rebuildRightHandSide(deep)  
     
    def __checkMatrixType(self):  
      """  
      reassess the matrix type and, if needed, initiates an operator rebuild  
      """  
      new_matrix_type=self.getDomain().getSystemMatrixTypeId(self.getSolverMethod(),self.isSymmetric())  
      if not new_matrix_type==self.__matrix_type:  
          if self.debug() : print "PDE Debug: Matrix type is now %d."%new_matrix_type  
          self.__matrix_type=new_matrix_type  
          self.__rebuildOperator(deep=True)  
   
    #============ assembling =======================================================  
    def __copyConstraint(self):  
       """  
       copies the constrint condition into u  
       """  
       if not self.__righthandside.isEmpty():  
          q=self.getCoefficientOfPDE("q")  
          r=self.getCoefficientOfPDE("r")  
          if not q.isEmpty():  
              if r.isEmpty():  
                 r2=escript.Data(0,self.__righthandside.getShape(),self.__righthandside.getFunctionSpace())  
              else:  
                 r2=escript.Data(r,self.__righthandside.getFunctionSpace())  
              self.__righthandside.copyWithMask(r2,escript.Data(q,self.__righthandside.getFunctionSpace()))  
   
    def __applyConstraint(self):  
        """  
        applies the constraints defined by q and r to the system  
        """  
        q=self.getCoefficientOfPDE("q")  
        r=self.getCoefficientOfPDE("r")  
        if not q.isEmpty() and not self.__operator.isEmpty():  
           # q is the row and column mask to indicate where constraints are set:  
           row_q=escript.Data(q,self.getFunctionSpaceForEquation())  
           col_q=escript.Data(q,self.getFunctionSpaceForSolution())  
           u=self.__getNewSolution()  
           if r.isEmpty():  
              r_s=self.__getNewSolution()  
           else:  
              r_s=escript.Data(r,self.getFunctionSpaceForSolution())  
           u.copyWithMask(r_s,col_q)  
           if self.isUsingLumping():  
              self.__operator.copyWithMask(escript.Data(1,q.getShape(),self.getFunctionSpaceForEquation()),row_q)  
           else:  
              if not self.__righthandside.isEmpty(): self.__righthandside-=self.__operator*u  
              self.__operator.nullifyRowsAndCols(row_q,col_q,1.)  
1196    
1197     def getSystem(self):     def getSystem(self):
1198         """         """
# Line 912  class LinearPDE: Line 1201  class LinearPDE:
1201         if not self.__operator_isValid or not self.__righthandside_isValid:         if not self.__operator_isValid or not self.__righthandside_isValid:
1202            if self.isUsingLumping():            if self.isUsingLumping():
1203                if not self.__operator_isValid:                if not self.__operator_isValid:
1204                   if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():                   if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution(): raise TypeError,"Lumped matrix requires same order for equations and unknowns"
1205                         raise TypeError,"Lumped matrix requires same order for equations and unknowns"                   if not self.getCoefficientOfGeneralPDE("A").isEmpty(): raise Warning,"Lumped matrix does not allow coefficient A"
1206                   if not self.getCoefficientOfPDE("A").isEmpty():                   if not self.getCoefficientOfGeneralPDE("B").isEmpty(): raise Warning,"Lumped matrix does not allow coefficient B"
1207                            raise Warning,"Lumped matrix does not allow coefficient A"                   if not self.getCoefficientOfGeneralPDE("C").isEmpty(): raise Warning,"Lumped matrix does not allow coefficient C"
                  if not self.getCoefficientOfPDE("B").isEmpty():  
                           raise Warning,"Lumped matrix does not allow coefficient B"  
                  if not self.getCoefficientOfPDE("C").isEmpty():  
                           raise Warning,"Lumped matrix does not allow coefficient C"  
                  if self.debug() : print "PDE Debug: New lumped operator is built."  
1208                   mat=self.__getNewOperator()                   mat=self.__getNewOperator()
1209                   self.getDomain().addPDEToSystem(mat,escript.Data(), \                   self.getDomain().addPDEToSystem(mat,escript.Data(), \
1210                             self.getCoefficientOfPDE("A"), \                             self.getCoefficientOfGeneralPDE("A"), \
1211                             self.getCoefficientOfPDE("B"), \                             self.getCoefficientOfGeneralPDE("B"), \
1212                             self.getCoefficientOfPDE("C"), \                             self.getCoefficientOfGeneralPDE("C"), \
1213                             self.getCoefficientOfPDE("D"), \                             self.getCoefficientOfGeneralPDE("D"), \
1214                             escript.Data(), \                             escript.Data(), \
1215                             escript.Data(), \                             escript.Data(), \
1216                             self.getCoefficientOfPDE("d"), \                             self.getCoefficientOfGeneralPDE("d"), \
1217                             escript.Data(),\                             escript.Data(),\
1218                             self.getCoefficientOfPDE("d_contact"), \                             self.getCoefficientOfGeneralPDE("d_contact"), \
1219                             escript.Data())                             escript.Data())
1220                   self.__operator=mat*escript.Data(1,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)                   self.__operator=1./(mat*escript.Data(1,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True))
1221                   self.__applyConstraint()                   del mat
1222                     self.trace("New lumped operator has been built.")
1223                   self.__operator_isValid=True                   self.__operator_isValid=True
1224                if not self.__righthandside_isValid:                if not self.__righthandside_isValid:
                  if self.debug() : print "PDE Debug: New right hand side is built."  
1225                   self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \                   self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \
1226                                 self.getCoefficientOfPDE("X"), \                                 self.getCoefficientOfGeneralPDE("X"), \
1227                                 self.getCoefficientOfPDE("Y"),\                                 self.getCoefficientOfGeneralPDE("Y"),\
1228                                 self.getCoefficientOfPDE("y"),\                                 self.getCoefficientOfGeneralPDE("y"),\
1229                                 self.getCoefficientOfPDE("y_contact"))                                 self.getCoefficientOfGeneralPDE("y_contact"))
1230                   self.__copyConstraint()                   self.trace("New right hand side as been built.")
1231                   self.__righthandside_isValid=True                   self.__righthandside_isValid=True
1232            else:            else:
1233               if not self.__operator_isValid and not self.__righthandside_isValid:               if not self.__operator_isValid and not self.__righthandside_isValid:
                  if self.debug() : print "PDE Debug: New system is built."  
1234                   self.getDomain().addPDEToSystem(self.__makeFreshOperator(),self.__makeFreshRightHandSide(), \                   self.getDomain().addPDEToSystem(self.__makeFreshOperator(),self.__makeFreshRightHandSide(), \
1235                                 self.getCoefficientOfPDE("A"), \                                 self.getCoefficientOfGeneralPDE("A"), \
1236                                 self.getCoefficientOfPDE("B"), \                                 self.getCoefficientOfGeneralPDE("B"), \
1237                                 self.getCoefficientOfPDE("C"), \                                 self.getCoefficientOfGeneralPDE("C"), \
1238                                 self.getCoefficientOfPDE("D"), \                                 self.getCoefficientOfGeneralPDE("D"), \
1239                                 self.getCoefficientOfPDE("X"), \                                 self.getCoefficientOfGeneralPDE("X"), \
1240                                 self.getCoefficientOfPDE("Y"), \                                 self.getCoefficientOfGeneralPDE("Y"), \
1241                                 self.getCoefficientOfPDE("d"), \                                 self.getCoefficientOfGeneralPDE("d"), \
1242                                 self.getCoefficientOfPDE("y"), \                                 self.getCoefficientOfGeneralPDE("y"), \
1243                                 self.getCoefficientOfPDE("d_contact"), \                                 self.getCoefficientOfGeneralPDE("d_contact"), \
1244                                 self.getCoefficientOfPDE("y_contact"))                                 self.getCoefficientOfGeneralPDE("y_contact"))
1245                   self.__applyConstraint()                   self.__applyConstraint()
1246                   self.__copyConstraint()                   self.__righthandside=self.copyConstraint(self.__righthandside)
1247                     self.trace("New system has been built.")
1248                   self.__operator_isValid=True                   self.__operator_isValid=True
1249                   self.__righthandside_isValid=True                   self.__righthandside_isValid=True
1250               elif not self.__righthandside_isValid:               elif not self.__righthandside_isValid:
                  if self.debug() : print "PDE Debug: New right hand side is built."  
1251                   self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \                   self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \
1252                                 self.getCoefficientOfPDE("X"), \                                 self.getCoefficientOfGeneralPDE("X"), \
1253                                 self.getCoefficientOfPDE("Y"),\                                 self.getCoefficientOfGeneralPDE("Y"),\
1254                                 self.getCoefficientOfPDE("y"),\                                 self.getCoefficientOfGeneralPDE("y"),\
1255                                 self.getCoefficientOfPDE("y_contact"))                                 self.getCoefficientOfGeneralPDE("y_contact"))
1256                   self.__copyConstraint()                   self.__righthandside=self.copyConstraint(self.__righthandside)
1257                     self.trace("New right hand side has been built.")
1258                   self.__righthandside_isValid=True                   self.__righthandside_isValid=True
1259               elif not self.__operator_isValid:               elif not self.__operator_isValid:
                  if self.debug() : print "PDE Debug: New operator is built."  
1260                   self.getDomain().addPDEToSystem(self.__makeFreshOperator(),escript.Data(), \                   self.getDomain().addPDEToSystem(self.__makeFreshOperator(),escript.Data(), \
1261                              self.getCoefficientOfPDE("A"), \                              self.getCoefficientOfGeneralPDE("A"), \
1262                              self.getCoefficientOfPDE("B"), \                              self.getCoefficientOfGeneralPDE("B"), \
1263                              self.getCoefficientOfPDE("C"), \                              self.getCoefficientOfGeneralPDE("C"), \
1264                              self.getCoefficientOfPDE("D"), \                              self.getCoefficientOfGeneralPDE("D"), \
1265                              escript.Data(), \                              escript.Data(), \
1266                              escript.Data(), \                              escript.Data(), \
1267                              self.getCoefficientOfPDE("d"), \                              self.getCoefficientOfGeneralPDE("d"), \
1268                              escript.Data(),\                              escript.Data(),\
1269                              self.getCoefficientOfPDE("d_contact"), \                              self.getCoefficientOfGeneralPDE("d_contact"), \
1270                              escript.Data())                              escript.Data())
1271                   self.__applyConstraint()                   self.__applyConstraint()
1272                     self.trace("New operator has been built.")
1273                   self.__operator_isValid=True                   self.__operator_isValid=True
1274         return (self.__operator,self.__righthandside)         return (self.__operator,self.__righthandside)
    def getOperator(self):  
        """  
        returns the operator of the PDE  
        """  
        return self.getSystem()[0]  
1275    
    def getRightHandSide(self):  
        """  
        returns the right hand side of the PDE  
        """  
        return self.getSystem()[1]  
1276    
    def solve(self,**options):  
       """  
       solve the PDE  
1277    
       @param options:  
       """  
       mat,f=self.getSystem()  
       if self.isUsingLumping():  
          out=f/mat  
       else:  
          options[util.TOLERANCE_KEY]=self.getTolerance()  
          options[util.METHOD_KEY]=self.getSolverMethod()  
          options[util.SYMMETRY_KEY]=self.isSymmetric()  
          if self.debug() : print "PDE Debug: solver options: ",options  
          out=mat.solve(f,options)  
       return out  
   
    def getSolution(self,**options):  
        """  
        returns the solution of the PDE  
   
        @param options:  
        """  
        if not self.__solution_isValid:  
            if self.debug() : print "PDE Debug: PDE is resolved."  
            self.__solution=self.solve(**options)  
            self.__solution_isValid=True  
        return self.__solution  
   
   
   
 def ELMAN_RAMAGE(P):  
      """   """  
      return (P-1.).wherePositive()*0.5*(1.-1./(P+1.e-15))  
 def SIMPLIFIED_BROOK_HUGHES(P):  
      """   """  
      c=(P-3.).whereNegative()  
      return P/6.*c+1./2.*(1.-c)  
   
 def HALF(P):  
     """ """  
     return escript.Scalar(0.5,P.getFunctionSpace())  
1278    
1279  class AdvectivePDE(LinearPDE):  class AdvectivePDE(LinearPDE):
1280     """     """
1281     Class to handle a linear PDE dominated by advective terms:     Class to handle a linear PDE dominated by advective terms:
1282      
1283     class to define a linear PDE of the form     class to define a linear PDE of the form
1284    
1285     \f[     \f[
# Line 1071  class AdvectivePDE(LinearPDE): Line 1304  class AdvectivePDE(LinearPDE):
1304     u_i=r_i \quad \mathrm{where} \quad q_i>0     u_i=r_i \quad \mathrm{where} \quad q_i>0
1305     \f]     \f]
1306     """     """
1307     def __init__(self,domain,numEquations=0,numSolutions=0,xi=ELMAN_RAMAGE):     def __init__(self,domain,numEquations=0,numSolutions=0,xi=None,debug=False):
1308        LinearPDE.__init__(self,domain,numEquations,numSolutions)        LinearPDE.__init__(self,domain,numEquations,numSolutions,debug)
1309        self.__xi=xi        if xi==None:
1310             self.__xi=AdvectivePDE.ELMAN_RAMAGE
1311          else:
1312             self.__xi=xi
1313        self.__Xi=escript.Data()        self.__Xi=escript.Data()
1314    
1315     def __calculateXi(self,peclet_factor,Z,h):     def __calculateXi(self,peclet_factor,Z,h):
# Line 1085  class AdvectivePDE(LinearPDE): Line 1321  class AdvectivePDE(LinearPDE):
1321    
1322     def setValue(self,**args):     def setValue(self,**args):
1323         if "A" in args.keys()   or "B" in args.keys() or "C" in args.keys(): self.__Xi=escript.Data()         if "A" in args.keys()   or "B" in args.keys() or "C" in args.keys(): self.__Xi=escript.Data()
1324         self._LinearPDE__setValue(**args)         LinearPDE.setValue(**args)
1325              
1326       def ELMAN_RAMAGE(P):
1327         """   """
1328         return (P-1.).wherePositive()*0.5*(1.-1./(P+1.e-15))
1329       def SIMPLIFIED_BROOK_HUGHES(P):
1330         """   """
1331         c=(P-3.).whereNegative()
1332         return P/6.*c+1./2.*(1.-c)
1333    
1334       def HALF(P):
1335        """ """
1336        return escript.Scalar(0.5,P.getFunctionSpace())
1337    
1338     def getXi(self):     def getXi(self):
1339        if self.__Xi.isEmpty():        if self.__Xi.isEmpty():
1340           B=self.getCoefficient("B")           B=self.getCoefficient("B")
# Line 1130  class AdvectivePDE(LinearPDE): Line 1378  class AdvectivePDE(LinearPDE):
1378                 self.__Xi=h*xi/(length_of_Z+Z_max*self.TOL)                 self.__Xi=h*xi/(length_of_Z+Z_max*self.TOL)
1379                 print "@ preclet number = %e"%util.Lsup(peclet_number),util.Lsup(xi),util.Lsup(length_of_Z)                 print "@ preclet number = %e"%util.Lsup(peclet_number),util.Lsup(xi),util.Lsup(length_of_Z)
1380        return self.__Xi        return self.__Xi
         
1381    
1382     def getCoefficientOfPDE(self,name):  
1383       def getCoefficientOfGeneralPDE(self,name):
1384       """       """
1385       return the value of the coefficient name of the general PDE       return the value of the coefficient name of the general PDE
1386    
# Line 1141  class AdvectivePDE(LinearPDE): Line 1389  class AdvectivePDE(LinearPDE):
1389       if not self.getNumEquations() == self.getNumSolutions():       if not self.getNumEquations() == self.getNumSolutions():
1390            raise ValueError,"AdvectivePDE expects the number of solution componets and the number of equations to be equal."            raise ValueError,"AdvectivePDE expects the number of solution componets and the number of equations to be equal."
1391    
1392       if name == "A" :       if name == "A" :
1393           A=self.getCoefficient("A")           A=self.getCoefficient("A")
1394           B=self.getCoefficient("B")           B=self.getCoefficient("B")
1395           C=self.getCoefficient("C")           C=self.getCoefficient("C")
1396           if B.isEmpty() and C.isEmpty():           if B.isEmpty() and C.isEmpty():
1397              Aout=A              Aout=A
1398           else:           else:
1399              if A.isEmpty():              if A.isEmpty():
1400                 Aout=self.createNewCoefficient("A")                 Aout=self.createNewCoefficient("A")
1401              else:              else:
1402                 Aout=A[:]                 Aout=A[:]
# Line 1174  class AdvectivePDE(LinearPDE): Line 1422  class AdvectivePDE(LinearPDE):
1422                        else:                        else:
1423                            Aout[j,l]+=Xi*C[j]*C[l]                            Aout[j,l]+=Xi*C[j]*C[l]
1424           return Aout           return Aout
1425       elif name == "B" :       elif name == "B" :
1426           B=self.getCoefficient("B")           B=self.getCoefficient("B")
1427           C=self.getCoefficient("C")           C=self.getCoefficient("C")
1428           D=self.getCoefficient("D")           D=self.getCoefficient("D")
# Line 1182  class AdvectivePDE(LinearPDE): Line 1430  class AdvectivePDE(LinearPDE):
1430              Bout=B              Bout=B
1431           else:           else:
1432              Xi=self.getXi()              Xi=self.getXi()
1433              if B.isEmpty():              if B.isEmpty():
1434                  Bout=self.createNewCoefficient("B")                  Bout=self.createNewCoefficient("B")
1435              else:              else:
1436                  Bout=B[:]                  Bout=B[:]
1437              if self.getNumEquations()>1:              if self.getNumEquations()>1:
1438                 for k in range(self.getNumSolutions()):                 for k in range(self.getNumSolutions()):
1439                    for p in range(self.getNumEquations()):                    for p in range(self.getNumEquations()):
1440                       tmp=Xi*D[p,k]                       tmp=Xi*D[p,k]
1441                       for i in range(self.getNumEquations()):                       for i in range(self.getNumEquations()):
1442                          for j in range(self.getDim()):                          for j in range(self.getDim()):
# Line 1197  class AdvectivePDE(LinearPDE): Line 1445  class AdvectivePDE(LinearPDE):
1445                 tmp=Xi*D                 tmp=Xi*D
1446                 for j in range(self.getDim()): Bout[j]+=tmp*C[j]                 for j in range(self.getDim()): Bout[j]+=tmp*C[j]
1447           return Bout           return Bout
1448       elif name == "C" :       elif name == "C" :
1449           B=self.getCoefficient("B")           B=self.getCoefficient("B")
1450           C=self.getCoefficient("C")           C=self.getCoefficient("C")
1451           D=self.getCoefficient("D")           D=self.getCoefficient("D")
# Line 1205  class AdvectivePDE(LinearPDE): Line 1453  class AdvectivePDE(LinearPDE):
1453              Cout=C              Cout=C
1454           else:           else:
1455              Xi=self.getXi()              Xi=self.getXi()
1456              if C.isEmpty():              if C.isEmpty():
1457                  Cout=self.createNewCoefficient("C")                  Cout=self.createNewCoefficient("C")
1458              else:              else:
1459                  Cout=C[:]                  Cout=C[:]
# Line 1220  class AdvectivePDE(LinearPDE): Line 1468  class AdvectivePDE(LinearPDE):
1468                 tmp=Xi*D                 tmp=Xi*D
1469                 for j in range(self.getDim()): Cout[j]+=tmp*B[j]                 for j in range(self.getDim()): Cout[j]+=tmp*B[j]
1470           return Cout           return Cout
1471       elif name == "D" :       elif name == "D" :
1472           return self.getCoefficient("D")           return self.getCoefficient("D")
1473       elif name == "X" :       elif name == "X" :
1474           X=self.getCoefficient("X")           X=self.getCoefficient("X")
1475           Y=self.getCoefficient("Y")           Y=self.getCoefficient("Y")
1476           B=self.getCoefficient("B")           B=self.getCoefficient("B")
# Line 1236  class AdvectivePDE(LinearPDE): Line 1484  class AdvectivePDE(LinearPDE):
1484                  Xout=X[:]                  Xout=X[:]
1485              Xi=self.getXi()              Xi=self.getXi()
1486              if self.getNumEquations()>1:              if self.getNumEquations()>1:
1487                   for p in range(self.getNumEquations()):                   for p in range(self.getNumEquations()):
1488                      tmp=Xi*Y[p]                      tmp=Xi*Y[p]
1489                      for i in range(self.getNumEquations()):                      for i in range(self.getNumEquations()):
1490                         for j in range(self.getDim()):                         for j in range(self.getDim()):
# Line 1256  class AdvectivePDE(LinearPDE): Line 1504  class AdvectivePDE(LinearPDE):
1504                      else:                      else:
1505                         Xout[j]+=tmp*C[j]                         Xout[j]+=tmp*C[j]
1506           return Xout           return Xout
1507       elif name == "Y" :       elif name == "Y" :
1508           return self.getCoefficient("Y")           return self.getCoefficient("Y")
1509       elif name == "d" :       elif name == "d" :
1510           return self.getCoefficient("d")           return self.getCoefficient("d")
1511       elif name == "y" :       elif name == "y" :
1512           return self.getCoefficient("y")           return self.getCoefficient("y")
1513       elif name == "d_contact" :       elif name == "d_contact" :
1514           return self.getCoefficient("d_contact")           return self.getCoefficient("d_contact")
1515       elif name == "y_contact" :       elif name == "y_contact" :
1516           return self.getCoefficient("y_contact")           return self.getCoefficient("y_contact")
1517       elif name == "r" :       elif name == "r" :
1518           return self.getCoefficient("r")           return self.getCoefficient("r")
1519       elif name == "q" :       elif name == "q" :
1520           return self.getCoefficient("q")           return self.getCoefficient("q")
1521       else:       else:
1522           raise SystemError,"unknown PDE coefficient %s",name           raise SystemError,"unknown PDE coefficient %s",name
# Line 1296  class Poisson(LinearPDE): Line 1544  class Poisson(LinearPDE):
1544     \f]     \f]
1545     """     """
1546    
1547     def __init__(self,domain,f=escript.Data(),q=escript.Data()):     def __init__(self,domain,f=escript.Data(),q=escript.Data(),debug=False):
1548         LinearPDE.__init__(self,domain,1,1)         LinearPDE.__init__(self,domain,1,1,debug)
1549         self.COEFFICIENTS={         self.COEFFICIENTS={"f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1550         "f"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),                            "q": PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.BOTH)}
        "q"         : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.BOTH)}  
1551         self.setSymmetryOn()         self.setSymmetryOn()
1552         self.setValue(f,q)         self.setValue(f,q)
1553    
# Line 1308  class Poisson(LinearPDE): Line 1555  class Poisson(LinearPDE):
1555         """set value of PDE parameters f and q"""         """set value of PDE parameters f and q"""
1556         self._LinearPDE__setValue(f=f,q=q)         self._LinearPDE__setValue(f=f,q=q)
1557    
1558     def getCoefficientOfPDE(self,name):     def getCoefficientOfGeneralPDE(self,name):
1559       """       """
1560       return the value of the coefficient name of the general PDE       return the value of the coefficient name of the general PDE
1561    
1562       @param name:       @param name:
1563       """       """
1564       if name == "A" :       if name == "A" :
1565           return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))           return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))
1566       elif name == "B" :       elif name == "B" :
1567           return escript.Data()           return escript.Data()
1568       elif name == "C" :       elif name == "C" :
1569           return escript.Data()           return escript.Data()
1570       elif name == "D" :       elif name == "D" :
1571           return escript.Data()           return escript.Data()
1572       elif name == "X" :       elif name == "X" :
1573           return escript.Data()           return escript.Data()
1574       elif name == "Y" :       elif name == "Y" :
1575           return self.getCoefficient("f")           return self.getCoefficient("f")
1576       elif name == "d" :       elif name == "d" :
1577           return escript.Data()           return escript.Data()
1578       elif name == "y" :       elif name == "y" :
1579           return escript.Data()           return escript.Data()
1580       elif name == "d_contact" :       elif name == "d_contact" :
1581           return escript.Data()           return escript.Data()
1582       elif name == "y_contact" :       elif name == "y_contact" :
1583           return escript.Data()           return escript.Data()
1584       elif name == "r" :       elif name == "r" :
1585           return escript.Data()           return escript.Data()
1586       elif name == "q" :       elif name == "q" :
1587           return self.getCoefficient("q")           return self.getCoefficient("q")
1588       else:       else:
1589           raise SystemError,"unknown PDE coefficient %s",name           raise SystemError,"unknown PDE coefficient %s",name
# Line 1363  class LameEquation(LinearPDE): Line 1610  class LameEquation(LinearPDE):
1610     \f]     \f]
1611     """     """
1612    
1613     def __init__(self,domain,f=escript.Data(),q=escript.Data()):     def __init__(self,domain,f=escript.Data(),q=escript.Data(),debug=False):
1614         LinearPDE.__init__(self,domain,domain.getDim(),domain.getDim())         LinearPDE.__init__(self,domain,domain.getDim(),domain.getDim(),debug)
1615         self.COEFFICIENTS={         self.COEFFICIENTS={ "lame_lambda"  : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),
1616         "lame_lambda"  : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),                            "lame_mu"      : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),
1617         "lame_mu"      : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),                            "F"            : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1618         "F"            : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),                            "sigma"        : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM),PDECoefficient.RIGHTHANDSIDE),
1619         "sigma"        : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM),PDECoefficient.RIGHTHANDSIDE),                            "f"            : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1620         "f"            : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),                            "r"            : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.BOTH),
1621         "r"            : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.BOTH),                            "q"            : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.BOTH)}
        "q"            : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.BOTH)}  
1622         self.setSymmetryOn()         self.setSymmetryOn()
1623    
1624     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()):     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()):
# Line 1383  class LameEquation(LinearPDE): Line 1629  class LameEquation(LinearPDE):
1629                                   sigma=sigma, \                                   sigma=sigma, \
1630                                   f=f, \                                   f=f, \
1631                                   r=r, \                                   r=r, \
1632                                   q=q)                                   q=q)
1633     def getCoefficientOfPDE(self,name):     def getCoefficientOfGeneralPDE(self,name):
1634       """       """
1635       return the value of the coefficient name of the general PDE       return the value of the coefficient name of the general PDE
1636    
1637       @param name:       @param name:
1638       """       """
1639       if name == "A" :       if name == "A" :
1640           out =self.createNewCoefficientOfPDE("A")           out =self.createCoefficientOfGeneralPDE("A")
1641           for i in range(self.getDim()):           for i in range(self.getDim()):
1642             for j in range(self.getDim()):             for j in range(self.getDim()):
1643               out[i,i,j,j] += self.getCoefficient("lame_lambda")               out[i,i,j,j] += self.getCoefficient("lame_lambda")
1644               out[i,j,j,i] += self.getCoefficient("lame_mu")               out[i,j,j,i] += self.getCoefficient("lame_mu")
1645               out[i,j,i,j] += self.getCoefficient("lame_mu")               out[i,j,i,j] += self.getCoefficient("lame_mu")
1646           return out           return out
1647       elif name == "B" :       elif name == "B" :
1648           return escript.Data()           return escript.Data()
1649       elif name == "C" :       elif name == "C" :
1650           return escript.Data()           return escript.Data()
1651       elif name == "D" :       elif name == "D" :
1652           return escript.Data()           return escript.Data()
1653       elif name == "X" :       elif name == "X" :
1654           return self.getCoefficient("sigma")           return self.getCoefficient("sigma")
1655       elif name == "Y" :       elif name == "Y" :
1656           return self.getCoefficient("F")           return self.getCoefficient("F")
1657       elif name == "d" :       elif name == "d" :
1658           return escript.Data()           return escript.Data()
1659       elif name == "y" :       elif name == "y" :
1660           return self.getCoefficient("f")           return self.getCoefficient("f")
1661       elif name == "d_contact" :       elif name == "d_contact" :
1662           return escript.Data()           return escript.Data()
1663       elif name == "y_contact" :       elif name == "y_contact" :
1664           return escript.Data()           return escript.Data()
1665       elif name == "r" :       elif name == "r" :
1666           return self.getCoefficient("r")           return self.getCoefficient("r")
1667       elif name == "q" :       elif name == "q" :
1668           return self.getCoefficient("q")           return self.getCoefficient("q")
1669       else:       else:
1670           raise SystemError,"unknown PDE coefficient %s",name           raise SystemError,"unknown PDE coefficient %s",name
1671    
1672  # $Log$  # $Log$
1673    # Revision 1.11  2005/08/23 01:24:28  jgs
1674    # Merge of development branch dev-02 back to main trunk on 2005-08-23
1675    #
1676  # Revision 1.10  2005/08/12 01:45:36  jgs  # Revision 1.10  2005/08/12 01:45:36  jgs
1677  # erge of development branch dev-02 back to main trunk on 2005-08-12  # erge of development branch dev-02 back to main trunk on 2005-08-12
1678  #  #
1679    # Revision 1.9.2.4  2005/08/22 07:11:09  gross
1680    # some problems with LinearPDEs fixed.
1681    #
1682    # Revision 1.9.2.3  2005/08/18 04:48:48  gross
1683    # the methods SetLumping*() are removed. Lumping is set trough setSolverMethod(LinearPDE.LUMPING)
1684    #
1685    # Revision 1.9.2.2  2005/08/18 04:39:32  gross
1686    # the constants have been removed from util.py as they not needed anymore. PDE related constants are accessed through LinearPDE attributes now
1687    #
1688  # Revision 1.9.2.1  2005/07/29 07:10:27  gross  # Revision 1.9.2.1  2005/07/29 07:10:27  gross
1689  # new functions in util and a new pde type in linearPDEs  # new functions in util and a new pde type in linearPDEs
1690  #  #

Legend:
Removed from v.147  
changed lines
  Added in v.148

  ViewVC Help
Powered by ViewVC 1.1.26