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

revision 122 by jgs, Thu Jun 9 05:38:05 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         # 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         switches off debugging
294         """
295         self.__debug=None
296
297       def trace(self,text):
298       """       """
299       return escript.Data(shape = getShapeOfCoefficient(name), \       print the text message if debugging is swiched on.
what = getFunctionSpaceForCoefficient(name))
300
301     def __del__(self):       @param name: name of the coefficient enquired.
302       pass       @type name: C{string}
303         """
304         if self.__debug: print "%s: %s"%(str(self),text)
305
306     def getCoefficient(self,name):     # =============================================================================
307       # some service functions:
308       # =============================================================================
309       def getDomain(self):
310       """       """
311       return the value of the parameter name       returns the domain of the PDE
312
313         @return : the domain of the PDE
314         @rtype : C{Domain}
315
@param name:
316       """       """
317       return self.COEFFICIENTS[name].getValue()       return self.__domain
318
319     def getCoefficientOfPDE(self,name):     def getDim(self):
320       """       """
321       return the value of the coefficient name of the general PDE.       returns the spatial dimension of the PDE
This method is called by the assembling routine it can be
overwritten to map coefficients of a particualr PDE to the general PDE.
322
323       @param name:       @return : the spatial dimension of the PDE domain
324         @rtype : C{int}
325       """       """
326       return self.getCoefficient(name)       return self.getDomain().getDim()
327
328     def hasCoefficient(self,name):     def getNumEquations(self):
329        """       """
330        return true if name is the name of a coefficient       returns the number of equations
331
332        @param name:       @return : the number of equations
333        """       @rtype : C{int}
334        return self.COEFFICIENTS.has_key(name)       @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       def getNumSolutions(self):
342         """
343         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
374
375        @param coefficients:     def getOperator(self):
"""
self.__setValue(**coefficients)

def cleanCoefficients(self):
"""
resets all coefficients to default values.
376       """       """
377       for i in self.COEFFICIENTS.iterkeys():       provides access to the operator of the PDE
self.COEFFICIENTS[i].resetValue()
378
379     def createNewCoefficient(self,name):       @return : the operator of the PDE
380       """       @rtype : L{Operator}
returns a new coefficient appropriate for coefficient name:
381       """       """
382       return escript.Data(0,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))       m=self.getSystem()[0]
383               if self.isUsingLumping():
384             return self.copyConstraint(1./m)
385         else:
386             return m
387
388     def getShapeOfCoefficient(self,name):     def getRightHandSide(self):
389       """       """
390       return the shape 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.hasCoefficient(name):       r=self.getSystem()[1]
396          return self.COEFFICIENTS[name].buildShape(self.getNumEquations(),self.getNumSolutions(),self.getDomain().getDim())       if self.isUsingLumping():
397             return self.copyConstraint(r)
398       else:       else:
399          raise ValueError,"Solution coefficient %s requested"%name           return r
400
401     def getFunctionSpaceForCoefficient(self,name):     def applyOperator(self,u=None):
402       """       """
403       return the atoms 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].getFunctionSpace(self.getDomain())            return self.getOperator()*self.getSolution()
413       else:       else:
414          raise ValueError,"Solution coefficient %s requested"%name          self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())
415
416     def alteredCoefficient(self,name):     def getResidual(self,u=None):
417       """       """
418       announce that coefficient name has been changed       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()
if self.COEFFICIENTS[name].isAlteringOperator(): self.__rebuildOperator()
if self.COEFFICIENTS[name].isAlteringRightHandSide(): self.__rebuildRightHandSide()
else:
raise ValueError,"unknown coefficient %s requested"%name
427
428     # ===== debug ==============================================================     def checkSymmetry(self,verbose=True):
429     def setDebugOn(self):        """
430         """        test the PDE for symmetry.
"""
self.__debug=not None
431
def setDebugOff(self):
"""
"""
self.__debug=None
432
433     def debug(self):       @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.
434         """       @type verbose: C{bool}
435         returns true if the PDE is in the debug mode       @return:  True if the PDE is symmetric.
436         """       @rtype : C{escript.Data}
return self.__debug
437
438     #===== Lumping ===========================        @note: This is a very expensive operation. It should be used for degugging only! The symmetry flag is not altered.
def setLumpingOn(self):
"""
indicates to use matrix lumping
439        """        """
440        if not self.isUsingLumping():        verbose=verbose or self.debug()
441           if self.debug() : print "PDE Debug: lumping is set on"        out=True
442           self.__rebuildOperator()        if self.getNumSolutions()!=self.getNumEquations():
443           self.__lumping=True           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 setLumpingOff(self):        return out
"""
switches off matrix lumping
"""
if self.isUsingLumping():
if self.debug() : print "PDE Debug: lumping is set off"
self.__rebuildOperator()
self.__lumping=False
496
497     def setLumping(self,flag=False):     def getSolution(self,**options):
498        """         """
499        set the matrix lumping flag to flag         returns the solution of the PDE. If the solution is not valid the PDE is solved.
"""
if flag:
self.setLumpingOn()
else:
self.setLumpingOff()
500
501     def isUsingLumping(self):         @return: the solution
502        """         @rtype: L{escript.Data}
503                 @param options: solver options
504        """         @keyword verbose:
505        return self.__lumping         @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           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       def getFlux(self,u=None):
528         """
529         returns the flux J_ij for a given u
530
531     #============ method business =========================================================         \f[
532     def setSolverMethod(self,solver=util.DEFAULT_METHOD):         J_ij=A_{ijkl}u_{k,l}+B_{ijk}u_k-X_{ij}
533           \f]
534
535         @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         @return : flux
538         @rtype : L{escript.Data}
539
540         """
541         if u==None: u=self.getSolution()
543
544       # =============================================================================
545       #   solver settings:
546       # =============================================================================
547       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 565  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())
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       # =============================================================================
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 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          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())
1113          return u
1114
1115     def __setValue(self,**coefficients):     def setValue(self,**coefficients):
1116        """        """
1117        sets new values to coefficient        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 739  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()
1185
1186     def __setHomogeneousConstraintFlag(self):        # check if the systrem is inhomogeneous:
1187        """        if len(coefficients)>0 and not self.isUsingLumping():
1188        checks if the constraints are homogeneous and sets self.__homogeneous_constraint accordingly.           q=self.getCoefficientOfGeneralPDE("q")
1189        """           r=self.getCoefficientOfGeneralPDE("r")
1190        self.__homogeneous_constraint=True           homogeneous_constraint=True
1191        q=self.getCoefficientOfPDE("q")           if not q.isEmpty() and not r.isEmpty():
1192        r=self.getCoefficientOfPDE("r")               if util.Lsup(q*r)>=1.e-13*util.Lsup(r):
1193        if not q.isEmpty() and not r.isEmpty():                 self.trace("Inhomogeneous constraint detected.")
1194           if (q*r).Lsup()>=1.e-13*r.Lsup(): self.__homogeneous_constraint=False                 self.__invalidateSystem()
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()
1195
1196
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())

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())
if self.isUsingLumping():
else:
if not self.__righthandside.isEmpty(): self.__righthandside-=self.__operator*u
self.__operator.nullifyRowsAndCols(row_q,col_q,1.)

1197     def getSystem(self):     def getSystem(self):
1198         """         """
1199         return the operator and right hand side of the PDE         return the operator and right hand side of the PDE
# Line 876  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()
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."
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."
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."
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."
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]

def getRightHandSide(self):
"""
returns the right hand side of the PDE
"""
return self.getSystem()[1]

def solve(self,**options):
"""
solve the PDE

@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

1275
1276
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)
1277
def HALF(P):
""" """
return escript.Scalar(0.5,P.getFunctionSpace())
1278
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[
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:
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):
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")
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
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[:]
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")
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()):
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")
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[:]
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")
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()):
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
1524
1525  class Poisson(LinearPDE):  class Poisson(LinearPDE):
1526     """     """
1527     Class to define a Poisson equstion problem:     Class to define a Poisson equation problem:
1528
1529     class to define a linear PDE of the form     class to define a linear PDE of the form
1530     \f[     \f[
# Line 1260  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
1554     def setValue(self,f=escript.Data(),q=escript.Data()):     def setValue(self,f=escript.Data(),q=escript.Data()):
1555           """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")
1588         else:
1589             raise SystemError,"unknown PDE coefficient %s",name
1590
1591    class LameEquation(LinearPDE):
1592       """
1593       Class to define a Lame equation problem:
1594
1595       class to define a linear PDE of the form
1596       \f[
1597       -(\mu (u_{i,j}+u_{j,i}))_{,j} - \lambda u_{j,ji}} = F_i -\sigma_{ij,j}
1598       \f]
1599
1600       with boundary conditons:
1601
1602       \f[
1603       n_j(\mu (u_{i,j}+u_{j,i})-sigma_{ij}) + n_i\lambda u_{j,j} = f_i
1604       \f]
1605
1606       and constraints:
1607
1608       \f[
1610       \f]
1611       """
1612
1613       def __init__(self,domain,f=escript.Data(),q=escript.Data(),debug=False):
1614           LinearPDE.__init__(self,domain,domain.getDim(),domain.getDim(),debug)
1615           self.COEFFICIENTS={ "lame_lambda"  : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),
1616                              "lame_mu"      : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),
1617                              "F"            : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1618                              "sigma"        : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM),PDECoefficient.RIGHTHANDSIDE),
1619                              "f"            : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1620                              "r"            : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.BOTH),
1621                              "q"            : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.BOTH)}
1622           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()):
1625           """set value of PDE parameters"""
1626           self._LinearPDE__setValue(lame_lambda=lame_lambda, \
1627                                     lame_mu=lame_mu, \
1628                                     F=F, \
1629                                     sigma=sigma, \
1630                                     f=f, \
1631                                     r=r, \
1632                                     q=q)
1633       def getCoefficientOfGeneralPDE(self,name):
1634         """
1635         return the value of the coefficient name of the general PDE
1636
1637         @param name:
1638         """
1639         if name == "A" :
1640             out =self.createCoefficientOfGeneralPDE("A")
1641             for i in range(self.getDim()):
1642               for j in range(self.getDim()):
1643                 out[i,i,j,j] += self.getCoefficient("lame_lambda")
1644                 out[i,j,j,i] += self.getCoefficient("lame_mu")
1645                 out[i,j,i,j] += self.getCoefficient("lame_mu")
1646             return out
1647         elif name == "B" :
1648             return escript.Data()
1649         elif name == "C" :
1650             return escript.Data()
1651         elif name == "D" :
1652             return escript.Data()
1653         elif name == "X" :
1654             return self.getCoefficient("sigma")
1655         elif name == "Y" :
1656             return self.getCoefficient("F")
1657         elif name == "d" :
1658             return escript.Data()
1659         elif name == "y" :
1660             return self.getCoefficient("f")
1661         elif name == "d_contact" :
1662             return escript.Data()
1663         elif name == "y_contact" :
1664             return escript.Data()
1665         elif name == "r" :
1666             return self.getCoefficient("r")
1667         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.8  2005/06/09 05:37:59  jgs  # Revision 1.11  2005/08/23 01:24:28  jgs
1674  # Merge of development branch back to main trunk on 2005-06-09  # 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
1677    # 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
1689    # new functions in util and a new pde type in linearPDEs
1690    #
1691    # Revision 1.1.2.25  2005/07/28 04:21:09  gross
1692    # Lame equation: (linear elastic, isotropic) added
1693  #  #
1694  # Revision 1.7  2005/05/06 04:26:10  jgs  # Revision 1.1.2.24  2005/07/22 06:37:11  gross
1695  # Merge of development branch back to main trunk on 2005-05-06  # some extensions to modellib and linearPDEs
1696  #  #
1697  # Revision 1.1.2.23  2005/05/13 00:55:20  cochrane  # Revision 1.1.2.23  2005/05/13 00:55:20  cochrane
1698  # Fixed up some docstrings.  Moved module-level functions to top of file so  # Fixed up some docstrings.  Moved module-level functions to top of file so

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