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

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

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

trunk/esys2/escript/py_src/linearPDEs.py revision 142 by jgs, Mon Jul 25 05:28:20 2005 UTC trunk/escript/py_src/linearPDEs.py revision 2431 by jfenwick, Wed May 20 03:39:57 2009 UTC
# Line 1  Line 1 
 # $Id$  
1    
2  ## @file linearPDEs.py  ########################################################
3    #
4    # Copyright (c) 2003-2008 by University of Queensland
5    # Earth Systems Science Computational Center (ESSCC)
6    # http://www.uq.edu.au/esscc
7    #
8    # Primary Business: Queensland, Australia
9    # Licensed under the Open Software License version 3.0
10    # http://www.opensource.org/licenses/osl-3.0.php
11    #
12    ########################################################
13    
14    __copyright__="""Copyright (c) 2003-2008 by University of Queensland
15    Earth Systems Science Computational Center (ESSCC)
16    http://www.uq.edu.au/esscc
17    Primary Business: Queensland, Australia"""
18    __license__="""Licensed under the Open Software License version 3.0
19    http://www.opensource.org/licenses/osl-3.0.php"""
20    __url__="https://launchpad.net/escript-finley"
21    
22  """  """
23  Functions and classes for linear PDEs  The module provides an interface to define and solve linear partial
24    differential equations (PDEs) and Transport problems within L{escript}.
25    L{linearPDEs} does not provide any solver capabilities in itself but hands the
26    PDE over to the PDE solver library defined through the L{Domain<escript.Domain>}
27    of the PDE. The general interface is provided through the L{LinearPDE} class.
28    L{TransportProblem} provides an interface to initial value problems dominated
29    by its advective terms.
30    
31    @var __author__: name of author
32    @var __copyright__: copyrights
33    @var __license__: licence agreement
34    @var __url__: url entry point on documentation
35    @var __version__: version
36    @var __date__: date of the version
37  """  """
38    
39    import math
40  import escript  import escript
41  import util  import util
42  import numarray  import numarray
43    
44    __author__="Lutz Gross, l.gross@uq.edu.au"
45    
46    
47  def _CompTuple2(t1,t2):  class IllegalCoefficient(ValueError):
48       """
49       Exception that is raised if an illegal coefficient of the general or
50       particular PDE is requested.
51     """     """
52     Compare two tuples     pass
53    
54     \param t1 The first tuple  class IllegalCoefficientValue(ValueError):
    \param t2 The second tuple  
55     """     """
56       Exception that is raised if an incorrect value for a coefficient is used.
57       """
58       pass
59    
60     dif=t1[0]+t1[1]-(t2[0]+t2[1])  class IllegalCoefficientFunctionSpace(ValueError):
61     if dif<0: return 1     """
62     elif dif>0: return -1     Exception that is raised if an incorrect function space for a coefficient
63     else: return 0     is used.
64       """
65    
66  def ELMAN_RAMAGE(P):  class UndefinedPDEError(ValueError):
67      return (P-1.).wherePositive()*0.5*(1.-1./(P+1.e-15))     """
68       Exception that is raised if a PDE is not fully defined yet.
69       """
70       pass
71    
72  def SIMPLIFIED_BROOK_HUGHES(P):  class PDECoef(object):
73      c=(P-3.).whereNegative()      """
74      return P/6.*c+1./2.*(1.-c)      A class for describing a PDE coefficient.
75    
76  def HALF(P):      @cvar INTERIOR: indicator that coefficient is defined on the interior of
77      return escript.Scalar(0.5,P.getFunctionSpace())                      the PDE domain
78        @cvar BOUNDARY: indicator that coefficient is defined on the boundary of
79                        the PDE domain
80        @cvar CONTACT: indicator that coefficient is defined on the contact region
81                       within the PDE domain
82        @cvar INTERIOR_REDUCED: indicator that coefficient is defined on the
83                                interior of the PDE domain using a reduced
84                                integration order
85        @cvar BOUNDARY_REDUCED: indicator that coefficient is defined on the
86                                boundary of the PDE domain using a reduced
87                                integration order
88        @cvar CONTACT_REDUCED: indicator that coefficient is defined on the contact
89                               region within the PDE domain using a reduced
90                               integration order
91        @cvar SOLUTION: indicator that coefficient is defined through a solution of
92                        the PDE
93        @cvar REDUCED: indicator that coefficient is defined through a reduced
94                       solution of the PDE
95        @cvar BY_EQUATION: indicator that the dimension of the coefficient shape is
96                           defined by the number of PDE equations
97        @cvar BY_SOLUTION: indicator that the dimension of the coefficient shape is
98                           defined by the number of PDE solutions
99        @cvar BY_DIM: indicator that the dimension of the coefficient shape is
100                      defined by the spatial dimension
101        @cvar OPERATOR: indicator that the the coefficient alters the operator of
102                        the PDE
103        @cvar RIGHTHANDSIDE: indicator that the the coefficient alters the right
104                             hand side of the PDE
105        @cvar BOTH: indicator that the the coefficient alters the operator as well
106                    as the right hand side of the PDE
107    
 class PDECoefficient:  
     """  
     A class for PDE coefficients  
108      """      """
     # identifier for location of Data objects defining COEFFICIENTS  
109      INTERIOR=0      INTERIOR=0
110      BOUNDARY=1      BOUNDARY=1
111      CONTACT=2      CONTACT=2
112      CONTINUOUS=3      SOLUTION=3
113      # identifier in the pattern of COEFFICIENTS:      REDUCED=4
114      # the pattern is a tuple of EQUATION,SOLUTION,DIM where DIM represents the spatial dimension, EQUATION the number of equations and SOLUTION the      BY_EQUATION=5
115      # number of unknowns.      BY_SOLUTION=6
116      EQUATION=3      BY_DIM=7
117      SOLUTION=4      OPERATOR=10
118      DIM=5      RIGHTHANDSIDE=11
119      # indicator for what is altered if the coefficient is altered:      BOTH=12
120      OPERATOR=5      INTERIOR_REDUCED=13
121      RIGHTHANDSIDE=6      BOUNDARY_REDUCED=14
122      BOTH=7      CONTACT_REDUCED=15
123      def __init__(self,where,pattern,altering):  
124         """      def __init__(self, where, pattern, altering):
125         Initialise a PDE Coefficient type         """
126           Initialises a PDE coefficient type.
127    
128           @param where: describes where the coefficient lives
129           @type where: one of L{INTERIOR}, L{BOUNDARY}, L{CONTACT}, L{SOLUTION},
130                        L{REDUCED}, L{INTERIOR_REDUCED}, L{BOUNDARY_REDUCED},
131                        L{CONTACT_REDUCED}
132           @param pattern: describes the shape of the coefficient and how the shape
133                           is built for a given spatial dimension and numbers of
134                           equations and solutions in then PDE. For instance,
135                           (L{BY_EQUATION},L{BY_SOLUTION},L{BY_DIM}) describes a
136                           rank 3 coefficient which is instantiated as shape (3,2,2)
137                           in case of three equations and two solution components
138                           on a 2-dimensional domain. In the case of single equation
139                           and a single solution component the shape components
140                           marked by L{BY_EQUATION} or L{BY_SOLUTION} are dropped.
141                           In this case the example would be read as (2,).
142           @type pattern: C{tuple} of L{BY_EQUATION}, L{BY_SOLUTION}, L{BY_DIM}
143           @param altering: indicates what part of the PDE is altered if the
144                            coefficient is altered
145           @type altering: one of L{OPERATOR}, L{RIGHTHANDSIDE}, L{BOTH}
146         """         """
147           super(PDECoef, self).__init__()
148         self.what=where         self.what=where
149         self.pattern=pattern         self.pattern=pattern
150         self.altering=altering         self.altering=altering
# Line 64  class PDECoefficient: Line 152  class PDECoefficient:
152    
153      def resetValue(self):      def resetValue(self):
154         """         """
155         resets coefficient value to default         Resets the coefficient value to the default.
156         """         """
157         self.value=escript.Data()         self.value=escript.Data()
158    
159      def getFunctionSpace(self,domain):      def getFunctionSpace(self,domain,reducedEquationOrder=False,reducedSolutionOrder=False):
160         """         """
161         defines the FunctionSpace of the coefficient on the domain         Returns the L{FunctionSpace<escript.FunctionSpace>} of the coefficient.
162    
163         @param domain:         @param domain: domain on which the PDE uses the coefficient
164         """         @type domain: L{Domain<escript.Domain>}
165         if self.what==self.INTERIOR: return escript.Function(domain)         @param reducedEquationOrder: True to indicate that reduced order is used
166         elif self.what==self.BOUNDARY: return escript.FunctionOnBoundary(domain)                                      to represent the equation
167         elif self.what==self.CONTACT: return escript.FunctionOnContactZero(domain)         @type reducedEquationOrder: C{bool}
168         elif self.what==self.CONTINUOUS: return escript.ContinuousFunction(domain)         @param reducedSolutionOrder: True to indicate that reduced order is used
169                                        to represent the solution
170           @type reducedSolutionOrder: C{bool}
171           @return: L{FunctionSpace<escript.FunctionSpace>} of the coefficient
172           @rtype: L{FunctionSpace<escript.FunctionSpace>}
173           """
174           if self.what==self.INTERIOR:
175                return escript.Function(domain)
176           elif self.what==self.INTERIOR_REDUCED:
177                return escript.ReducedFunction(domain)
178           elif self.what==self.BOUNDARY:
179                return escript.FunctionOnBoundary(domain)
180           elif self.what==self.BOUNDARY_REDUCED:
181                return escript.ReducedFunctionOnBoundary(domain)
182           elif self.what==self.CONTACT:
183                return escript.FunctionOnContactZero(domain)
184           elif self.what==self.CONTACT_REDUCED:
185                return escript.ReducedFunctionOnContactZero(domain)
186           elif self.what==self.SOLUTION:
187                if reducedEquationOrder and reducedSolutionOrder:
188                    return escript.ReducedSolution(domain)
189                else:
190                    return escript.Solution(domain)
191           elif self.what==self.REDUCED:
192                return escript.ReducedSolution(domain)
193    
194      def getValue(self):      def getValue(self):
195         """         """
196         returns the value of the coefficient:         Returns the value of the coefficient.
197    
198           @return: value of the coefficient
199           @rtype: L{Data<escript.Data>}
200         """         """
201         return self.value         return self.value
202        
203      def setValue(self,newValue):      def setValue(self,domain,numEquations=1,numSolutions=1,reducedEquationOrder=False,reducedSolutionOrder=False,newValue=None):
        """  
        set the value of the coefficient to new value  
204         """         """
205           Sets the value of the coefficient to a new value.
206    
207           @param domain: domain on which the PDE uses the coefficient
208           @type domain: L{Domain<escript.Domain>}
209           @param numEquations: number of equations of the PDE
210           @type numEquations: C{int}
211           @param numSolutions: number of components of the PDE solution
212           @type numSolutions: C{int}
213           @param reducedEquationOrder: True to indicate that reduced order is used
214                                        to represent the equation
215           @type reducedEquationOrder: C{bool}
216           @param reducedSolutionOrder: True to indicate that reduced order is used
217                                        to represent the solution
218           @type reducedSolutionOrder: C{bool}
219           @param newValue: number of components of the PDE solution
220           @type newValue: any object that can be converted into a
221                           L{Data<escript.Data>} object with the appropriate shape
222                           and L{FunctionSpace<escript.FunctionSpace>}
223           @raise IllegalCoefficientValue: if the shape of the assigned value does
224                                           not match the shape of the coefficient
225           @raise IllegalCoefficientFunctionSpace: if unable to interpolate value
226                                                   to appropriate function space
227           """
228           if newValue==None:
229               newValue=escript.Data()
230           elif isinstance(newValue,escript.Data):
231               if not newValue.isEmpty():
232                  if not newValue.getFunctionSpace() == self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder):
233                    try:
234                      newValue=escript.Data(newValue,self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder))
235                    except:
236                      raise IllegalCoefficientFunctionSpace,"Unable to interpolate coefficient to function space %s"%self.getFunctionSpace(domain)
237           else:
238               newValue=escript.Data(newValue,self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder))
239           if not newValue.isEmpty():
240               if not self.getShape(domain,numEquations,numSolutions)==newValue.getShape():
241                   raise IllegalCoefficientValue,"Expected shape of coefficient is %s but actual shape is %s."%(self.getShape(domain,numEquations,numSolutions),newValue.getShape())
242         self.value=newValue         self.value=newValue
243        
244      def isAlteringOperator(self):      def isAlteringOperator(self):
245          """          """
246      return true if the operator of the PDE is changed when the coefficient is changed          Checks if the coefficient alters the operator of the PDE.
247      """  
248            @return: True if the operator of the PDE is changed when the
249                     coefficient is changed
250            @rtype: C{bool}
251            """
252          if self.altering==self.OPERATOR or self.altering==self.BOTH:          if self.altering==self.OPERATOR or self.altering==self.BOTH:
253              return not None              return not None
254          else:          else:
# Line 102  class PDECoefficient: Line 256  class PDECoefficient:
256    
257      def isAlteringRightHandSide(self):      def isAlteringRightHandSide(self):
258          """          """
259      return true if the right hand side of the PDE is changed when the coefficient is changed          Checks if the coefficient alters the right hand side of the PDE.
260      """  
261            @rtype: C{bool}
262            @return: True if the right hand side of the PDE is changed when the
263                     coefficient is changed, C{None} otherwise.
264            """
265          if self.altering==self.RIGHTHANDSIDE or self.altering==self.BOTH:          if self.altering==self.RIGHTHANDSIDE or self.altering==self.BOTH:
266              return not None              return not None
267          else:          else:
268              return None              return None
269    
270      def estimateNumEquationsAndNumSolutions(self,shape=(),dim=3):      def estimateNumEquationsAndNumSolutions(self,domain,shape=()):
271         """         """
272         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 and number of solutions if
273           the coefficient has the given shape.
274    
275         @param shape:         @param domain: domain on which the PDE uses the coefficient
276         @param dim:         @type domain: L{Domain<escript.Domain>}
277           @param shape: suggested shape of the coefficient
278           @type shape: C{tuple} of C{int} values
279           @return: the number of equations and number of solutions of the PDE if
280                    the coefficient has given shape. If no appropriate numbers
281                    could be identified, C{None} is returned
282           @rtype: C{tuple} of two C{int} values or C{None}
283         """         """
284           dim=domain.getDim()
285         if len(shape)>0:         if len(shape)>0:
286             num=max(shape)+1             num=max(shape)+1
287         else:         else:
288             num=1             num=1
289         search=[]         search=[]
290         for u in range(num):         if self.definesNumEquation() and self.definesNumSolutions():
291            for e in range(num):            for u in range(num):
292               search.append((e,u))               for e in range(num):
293         search.sort(_CompTuple2)                  search.append((e,u))
294         for item in search:            search.sort(self.__CompTuple2)
295               s=self.buildShape(item[0],item[1],dim)            for item in search:
296                 s=self.getShape(domain,item[0],item[1])
297               if len(s)==0 and len(shape)==0:               if len(s)==0 and len(shape)==0:
298                   return (1,1)                   return (1,1)
299               else:               else:
300                   if s==shape: return item                   if s==shape: return item
301           elif self.definesNumEquation():
302              for e in range(num,0,-1):
303                 s=self.getShape(domain,e,0)
304                 if len(s)==0 and len(shape)==0:
305                     return (1,None)
306                 else:
307                     if s==shape: return (e,None)
308    
309           elif self.definesNumSolutions():
310              for u in range(num,0,-1):
311                 s=self.getShape(domain,0,u)
312                 if len(s)==0 and len(shape)==0:
313                     return (None,1)
314                 else:
315                     if s==shape: return (None,u)
316         return None         return None
317    
318      def buildShape(self,e=1,u=1,dim=3):      def definesNumSolutions(self):
319          """         """
320      builds the required shape for a given number of equations e, number of unknowns u and spatial dimension dim         Checks if the coefficient allows to estimate the number of solution
321           components.
322    
323      @param e:         @return: True if the coefficient allows an estimate of the number of
324      @param u:                  solution components, False otherwise
325      @param dim:         @rtype: C{bool}
326      """         """
327          s=()         for i in self.pattern:
328          for i in self.pattern:               if i==self.BY_SOLUTION: return True
329               if i==self.EQUATION:         return False
                 if e>1: s=s+(e,)  
              elif i==self.SOLUTION:  
                 if u>1: s=s+(u,)  
              else:  
                 s=s+(dim,)  
         return s  
330    
331  class LinearPDE:      def definesNumEquation(self):
332     """         """
333     Class to handle a linear PDE         Checks if the coefficient allows to estimate the number of equations.
334      
335     class to define a linear PDE of the form         @return: True if the coefficient allows an estimate of the number of
336                    equations, False otherwise
337           @rtype: C{bool}
338           """
339           for i in self.pattern:
340                 if i==self.BY_EQUATION: return True
341           return False
342    
343     \f[      def __CompTuple2(self,t1,t2):
344       -(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        """
345     \f]        Compares two tuples of possible number of equations and number of
346          solutions.
347    
348     with boundary conditons:        @param t1: the first tuple
349          @param t2: the second tuple
350          @return: 0, 1, or -1
351          """
352    
353     \f[        dif=t1[0]+t1[1]-(t2[0]+t2[1])
354     n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_{ik}u_k = - n_j*X_{ij} + y_i        if dif<0: return 1
355     \f]        elif dif>0: return -1
356          else: return 0
357    
358     and contact conditions      def getShape(self,domain,numEquations=1,numSolutions=1):
359           """
360           Builds the required shape of the coefficient.
361    
362     \f[         @param domain: domain on which the PDE uses the coefficient
363     n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_contact_{ik}[u_k] = - n_j*X_{ij} + y_contact_i         @type domain: L{Domain<escript.Domain>}
364     \f]         @param numEquations: number of equations of the PDE
365           @type numEquations: C{int}
366           @param numSolutions: number of components of the PDE solution
367           @type numSolutions: C{int}
368           @return: shape of the coefficient
369           @rtype: C{tuple} of C{int} values
370           """
371           dim=domain.getDim()
372           s=()
373           for i in self.pattern:
374                 if i==self.BY_EQUATION:
375                    if numEquations>1: s=s+(numEquations,)
376                 elif i==self.BY_SOLUTION:
377                    if numSolutions>1: s=s+(numSolutions,)
378                 else:
379                    s=s+(dim,)
380           return s
381    
382     and constraints:  class LinearProblem(object):
383       """
384       This is the base class to define a general linear PDE-type problem for
385       for an unknown function M{u} on a given domain defined through a
386       L{Domain<escript.Domain>} object. The problem can be given as a single
387       equation or as a system of equations.
388    
389       The class assumes that some sort of assembling process is required to form
390       a problem of the form
391    
392       M{L u=f}
393    
394       where M{L} is an operator and M{f} is the right hand side. This operator
395       problem will be solved to get the unknown M{u}.
396    
397       @cvar DEFAULT: The default method used to solve the system of linear
398                      equations
399       @cvar DIRECT: The direct solver based on LDU factorization
400       @cvar CHOLEVSKY: The direct solver based on LDLt factorization (can only be
401                        applied for symmetric PDEs)
402       @cvar PCG: The preconditioned conjugate gradient method (can only be applied
403                  for symmetric PDEs)
404       @cvar CR: The conjugate residual method
405       @cvar CGS: The conjugate gradient square method
406       @cvar BICGSTAB: The stabilized BiConjugate Gradient method
407       @cvar TFQMR: Transport Free Quasi Minimal Residual method
408       @cvar MINRES: Minimum residual method
409       @cvar SSOR: The symmetric overrelaxation method
410       @cvar ILU0: The incomplete LU factorization preconditioner with no fill-in
411       @cvar ILUT: The incomplete LU factorization preconditioner with fill-in
412       @cvar JACOBI: The Jacobi preconditioner
413       @cvar GMRES: The Gram-Schmidt minimum residual method
414       @cvar PRES20: Special GMRES with restart after 20 steps and truncation after
415                     5 residuals
416       @cvar LUMPING: Matrix lumping
417       @cvar NO_REORDERING: No matrix reordering allowed
418       @cvar MINIMUM_FILL_IN: Reorder matrix to reduce fill-in during factorization
419       @cvar NESTED_DISSECTION: Reorder matrix to improve load balancing during
420                                factorization
421       @cvar PASO: PASO solver package
422       @cvar SCSL: SGI SCSL solver library
423       @cvar MKL: Intel's MKL solver library
424       @cvar UMFPACK: The UMFPACK library
425       @cvar TRILINOS: The TRILINOS parallel solver class library from Sandia Natl
426                       Labs
427       @cvar ITERATIVE: The default iterative solver
428       @cvar AMG: Algebraic Multi Grid
429       @cvar RILU: Recursive ILU
430       @cvar GS: Gauss-Seidel solver
431    
432     \f[     """
433     u_i=r_i \quad \mathrm{where} \quad q_i>0     DEFAULT= 0
434     \f]     DIRECT= 1
435       CHOLEVSKY= 2
436     """     PCG= 3
437     TOL=1.e-13     CR= 4
438     DEFAULT_METHOD=util.DEFAULT_METHOD     CGS= 5
439     DIRECT=util.DIRECT     BICGSTAB= 6
440     CHOLEVSKY=util.CHOLEVSKY     SSOR= 7
441     PCG=util.PCG     ILU0= 8
442     CR=util.CR     ILUT= 9
443     CGS=util.CGS     JACOBI= 10
444     BICGSTAB=util.BICGSTAB     GMRES= 11
445     SSOR=util.SSOR     PRES20= 12
446     GMRES=util.GMRES     LUMPING= 13
447     PRES20=util.PRES20     NO_REORDERING= 17
448     ILU0=util.ILU0     MINIMUM_FILL_IN= 18
449     JACOBI=util.JACOBI     NESTED_DISSECTION= 19
450       SCSL= 14
451     def __init__(self,domain,numEquations=0,numSolutions=0):     MKL= 15
452       """     UMFPACK= 16
453       initializes a new linear PDE.     ITERATIVE= 20
454       PASO= 21
455       @param args:     AMG= 22
456       """     RILU = 23
457       # COEFFICIENTS can be overwritten by subclasses:     TRILINOS = 24
458       self.COEFFICIENTS={     NONLINEAR_GMRES = 25
459         "A"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM,PDECoefficient.SOLUTION,PDECoefficient.DIM),PDECoefficient.OPERATOR),     TFQMR = 26
460         "B"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),     MINRES = 27
461         "C"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION,PDECoefficient.DIM),PDECoefficient.OPERATOR),     GS=28
462         "D"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),  
463         "X"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM),PDECoefficient.RIGHTHANDSIDE),     SMALL_TOLERANCE=1.e-13
464         "Y"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),     PACKAGE_KEY="package"
465         "d"         : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),     METHOD_KEY="method"
466         "y"         : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),     SYMMETRY_KEY="symmetric"
467         "d_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),     TOLERANCE_KEY="tolerance"
468         "y_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),     PRECONDITIONER_KEY="preconditioner"
469         "r"         : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),  
470         "q"         : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.SOLUTION,),PDECoefficient.BOTH)}  
471       def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):
472         """
473         Initializes a linear problem.
474    
475         @param domain: domain of the PDE
476         @type domain: L{Domain<escript.Domain>}
477         @param numEquations: number of equations. If C{None} the number of
478                              equations is extracted from the coefficients.
479         @param numSolutions: number of solution components. If C{None} the number
480                              of solution components is extracted from the
481                              coefficients.
482         @param debug: if True debug information is printed
483    
484       # initialize attributes       """
485       self.__debug=None       super(LinearProblem, self).__init__()
486    
487         self.__debug=debug
488       self.__domain=domain       self.__domain=domain
489       self.__numEquations=numEquations       self.__numEquations=numEquations
490       self.__numSolutions=numSolutions       self.__numSolutions=numSolutions
491       self.cleanCoefficients()       self.__altered_coefficients=False
492         self.__reduce_equation_order=False
493       self.__operator=escript.Operator()       self.__reduce_solution_order=False
      self.__operator_isValid=False  
      self.__righthandside=escript.Data()  
      self.__righthandside_isValid=False  
      self.__solution=escript.Data()  
      self.__solution_isValid=False  
   
      # set some default values:  
      self.__homogeneous_constraint=True  
      self.__row_function_space=escript.Solution(self.__domain)  
      self.__column_function_space=escript.Solution(self.__domain)  
494       self.__tolerance=1.e-8       self.__tolerance=1.e-8
495       self.__solver_method=util.DEFAULT_METHOD       self.__solver_method=self.DEFAULT
496       self.__matrix_type=self.__domain.getSystemMatrixTypeId(util.DEFAULT_METHOD,False)       self.__solver_package=self.DEFAULT
497         self.__preconditioner=self.DEFAULT
498         self.__is_RHS_valid=False
499         self.__is_operator_valid=False
500       self.__sym=False       self.__sym=False
501       self.__lumping=False       self.__COEFFICIENTS={}
502         # initialize things:
503     def createCoefficient(self, name):       self.resetAllCoefficients()
504         self.__system_type=None
505         self.updateSystemType()
506       # ==========================================================================
507       #    general stuff:
508       # ==========================================================================
509       def __str__(self):
510         """
511         Returns a string representation of the PDE.
512    
513         @return: a simple representation of the PDE
514         @rtype: C{str}
515         """
516         return "<LinearProblem %d>"%id(self)
517       # ==========================================================================
518       #    debug :
519       # ==========================================================================
520       def setDebugOn(self):
521       """       """
522       create a data object corresponding to coefficient name       Switches debug output on.
      @param name:  
523       """       """
524       return escript.Data(shape = getShapeOfCoefficient(name), \       self.__debug=not None
                          what = getFunctionSpaceForCoefficient(name))  
525    
526     def __del__(self):     def setDebugOff(self):
527       pass       """
528         Switches debug output off.
529         """
530         self.__debug=None
531    
532     def getCoefficient(self,name):     def setDebug(self, flag):
533       """       """
534       return the value of the parameter name       Switches debug output on if C{flag} is True otherwise it is switched off.
535    
536       @param name:       @param flag: desired debug status
537         @type flag: C{bool}
538       """       """
539       return self.COEFFICIENTS[name].getValue()       if flag:
540             self.setDebugOn()
541         else:
542             self.setDebugOff()
543    
544     def getCoefficientOfPDE(self,name):     def trace(self,text):
545       """       """
546       return the value of the coefficient name of the general PDE.       Prints the text message if debug mode is switched on.
      This method is called by the assembling routine it can be  
      overwritten to map coefficients of a particualr PDE to the general PDE.  
547    
548       @param name:       @param text: message to be printed
549         @type text: C{string}
550       """       """
551       return self.getCoefficient(name)       if self.__debug: print "%s: %s"%(str(self),text)
552    
553     def hasCoefficient(self,name):     # ==========================================================================
554        """     # some service functions:
555        return true if name is the name of a coefficient     # ==========================================================================
556       def introduceCoefficients(self,**coeff):
557           """
558           Introduces new coefficients into the problem.
559    
560        @param name:         Use::
       """  
       return self.COEFFICIENTS.has_key(name)  
561    
562     def getFunctionSpaceForEquation(self):         p.introduceCoefficients(A=PDECoef(...), B=PDECoef(...))
563    
564           to introduce the coefficients M{A} ans M{B}.
565           """
566           for name, type in coeff.items():
567               if not isinstance(type,PDECoef):
568                  raise ValueError,"coefficient %s has no type."%name
569               self.__COEFFICIENTS[name]=type
570               self.__COEFFICIENTS[name].resetValue()
571               self.trace("coefficient %s has been introduced."%name)
572    
573       def getDomain(self):
574       """       """
575       return true if the test functions should use reduced order       Returns the domain of the PDE.
576    
577         @return: the domain of the PDE
578         @rtype: L{Domain<escript.Domain>}
579       """       """
580       return self.__row_function_space       return self.__domain
581    
582     def getFunctionSpaceForSolution(self):     def getDim(self):
583       """       """
584       return true if the interpolation of the solution should use reduced order       Returns the spatial dimension of the PDE.
585    
586         @return: the spatial dimension of the PDE domain
587         @rtype: C{int}
588       """       """
589       return self.__column_function_space       return self.getDomain().getDim()
590    
591     def setValue(self,**coefficients):     def getNumEquations(self):
592        """       """
593        sets new values to coefficients       Returns the number of equations.
594    
595        @param coefficients:       @return: the number of equations
596        """       @rtype: C{int}
597        self.__setValue(**coefficients)       @raise UndefinedPDEError: if the number of equations is not specified yet
598               """
599         if self.__numEquations==None:
600             if self.__numSolutions==None:
601                raise UndefinedPDEError,"Number of equations is undefined. Please specify argument numEquations."
602             else:
603                self.__numEquations=self.__numSolutions
604         return self.__numEquations
605    
606     def cleanCoefficients(self):     def getNumSolutions(self):
      """  
      resets all coefficients to default values.  
607       """       """
608       for i in self.COEFFICIENTS.iterkeys():       Returns the number of unknowns.
          self.COEFFICIENTS[i].resetValue()  
609    
610     def createNewCoefficient(self,name):       @return: the number of unknowns
611         @rtype: C{int}
612         @raise UndefinedPDEError: if the number of unknowns is not specified yet
613         """
614         if self.__numSolutions==None:
615            if self.__numEquations==None:
616                raise UndefinedPDEError,"Number of solution is undefined. Please specify argument numSolutions."
617            else:
618                self.__numSolutions=self.__numEquations
619         return self.__numSolutions
620    
621       def reduceEquationOrder(self):
622       """       """
623       returns a new coefficient appropriate for coefficient name:       Returns the status of order reduction for the equation.
624    
625         @return: True if reduced interpolation order is used for the
626                  representation of the equation, False otherwise
627         @rtype: L{bool}
628       """       """
629       return escript.Data(0,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))       return self.__reduce_equation_order
         
630    
631     def getShapeOfCoefficient(self,name):     def reduceSolutionOrder(self):
632       """       """
633       return the shape of the coefficient name       Returns the status of order reduction for the solution.
634    
635       @param name:       @return: True if reduced interpolation order is used for the
636                  representation of the solution, False otherwise
637         @rtype: L{bool}
638       """       """
639       if self.hasCoefficient(name):       return self.__reduce_solution_order
         return self.COEFFICIENTS[name].buildShape(self.getNumEquations(),self.getNumSolutions(),self.getDomain().getDim())  
      else:  
         raise ValueError,"Solution coefficient %s requested"%name  
640    
641     def getFunctionSpaceForCoefficient(self,name):     def getFunctionSpaceForEquation(self):
642       """       """
643       return the atoms of the coefficient name       Returns the L{FunctionSpace<escript.FunctionSpace>} used to discretize
644         the equation.
645    
646       @param name:       @return: representation space of equation
647         @rtype: L{FunctionSpace<escript.FunctionSpace>}
648       """       """
649       if self.hasCoefficient(name):       if self.reduceEquationOrder():
650          return self.COEFFICIENTS[name].getFunctionSpace(self.getDomain())           return escript.ReducedSolution(self.getDomain())
651       else:       else:
652          raise ValueError,"Solution coefficient %s requested"%name           return escript.Solution(self.getDomain())
653    
654     def alteredCoefficient(self,name):     def getFunctionSpaceForSolution(self):
655       """       """
656       announce that coefficient name has been changed       Returns the L{FunctionSpace<escript.FunctionSpace>} used to represent
657         the solution.
658    
659       @param name:       @return: representation space of solution
660         @rtype: L{FunctionSpace<escript.FunctionSpace>}
661       """       """
662       if self.hasCoefficient(name):       if self.reduceSolutionOrder():
663          if self.COEFFICIENTS[name].isAlteringOperator(): self.__rebuildOperator()           return escript.ReducedSolution(self.getDomain())
         if self.COEFFICIENTS[name].isAlteringRightHandSide(): self.__rebuildRightHandSide()  
664       else:       else:
665          raise ValueError,"unknown coefficient %s requested"%name           return escript.Solution(self.getDomain())
666    
667     # ===== debug ==============================================================     # ==========================================================================
668     def setDebugOn(self):     #   solver settings:
669         """     # ==========================================================================
670         """     def setSolverMethod(self,solver=None,preconditioner=None):
671         self.__debug=not None         """
672           Sets a new solver method and/or preconditioner.
673    
674           @param solver: new solver method to use
675           @type solver: one of L{DEFAULT}, L{ITERATIVE} L{DIRECT}, L{CHOLEVSKY},
676                         L{PCG}, L{CR}, L{CGS}, L{BICGSTAB}, L{SSOR}, L{GMRES},
677                         L{TFQMR}, L{MINRES}, L{PRES20}, L{LUMPING}, L{AMG}
678           @param preconditioner: new preconditioner to use
679           @type preconditioner: one of L{DEFAULT}, L{JACOBI} L{ILU0}, L{ILUT},
680                                 L{SSOR}, L{RILU}, L{GS}
681    
682           @note: the solver method chosen may not be available by the selected
683                  package.
684           """
685           if solver==None: solver=self.__solver_method
686           if preconditioner==None: preconditioner=self.__preconditioner
687           if solver==None: solver=self.DEFAULT
688           if preconditioner==None: preconditioner=self.DEFAULT
689           if not (solver,preconditioner)==self.getSolverMethod():
690               self.__solver_method=solver
691               self.__preconditioner=preconditioner
692               self.updateSystemType()
693               self.trace("New solver is %s"%self.getSolverMethodName())
694    
695       def getSolverMethodName(self):
696           """
697           Returns the name of the solver currently used.
698    
699           @return: the name of the solver currently used
700           @rtype: C{string}
701           """
702    
703           m=self.getSolverMethod()
704           p=self.getSolverPackage()
705           method=""
706           if m[0]==self.DEFAULT: method="DEFAULT"
707           elif m[0]==self.DIRECT: method= "DIRECT"
708           elif m[0]==self.ITERATIVE: method= "ITERATIVE"
709           elif m[0]==self.CHOLEVSKY: method= "CHOLEVSKY"
710           elif m[0]==self.PCG: method= "PCG"
711           elif m[0]==self.TFQMR: method= "TFQMR"
712           elif m[0]==self.MINRES: method= "MINRES"
713           elif m[0]==self.CR: method= "CR"
714           elif m[0]==self.CGS: method= "CGS"
715           elif m[0]==self.BICGSTAB: method= "BICGSTAB"
716           elif m[0]==self.SSOR: method= "SSOR"
717           elif m[0]==self.GMRES: method= "GMRES"
718           elif m[0]==self.PRES20: method= "PRES20"
719           elif m[0]==self.LUMPING: method= "LUMPING"
720           elif m[0]==self.AMG: method= "AMG"
721           if m[1]==self.DEFAULT: method+="+DEFAULT"
722           elif m[1]==self.JACOBI: method+= "+JACOBI"
723           elif m[1]==self.ILU0: method+= "+ILU0"
724           elif m[1]==self.ILUT: method+= "+ILUT"
725           elif m[1]==self.SSOR: method+= "+SSOR"
726           elif m[1]==self.AMG: method+= "+AMG"
727           elif m[1]==self.RILU: method+= "+RILU"
728           elif m[1]==self.GS: method+= "+GS"
729           if p==self.DEFAULT: package="DEFAULT"
730           elif p==self.PASO: package= "PASO"
731           elif p==self.MKL: package= "MKL"
732           elif p==self.SCSL: package= "SCSL"
733           elif p==self.UMFPACK: package= "UMFPACK"
734           elif p==self.TRILINOS: package= "TRILINOS"
735           else : method="unknown"
736           return "%s solver of %s package"%(method,package)
737    
738     def setDebugOff(self):     def getSolverMethod(self):
        """  
739         """         """
740         self.__debug=None         Returns the solver method and preconditioner used.
741    
742     def debug(self):         @return: the solver method currently used.
743           @rtype: C{tuple} of C{int}
744         """         """
745         returns true if the PDE is in the debug mode         return self.__solver_method,self.__preconditioner
746    
747       def setSolverPackage(self,package=None):
748         """         """
749         return self.__debug         Sets a new solver package.
750    
751     #===== Lumping ===========================         @param package: the new solver package
752     def setLumpingOn(self):         @type package: one of L{DEFAULT}, L{PASO} L{SCSL}, L{MKL}, L{UMFPACK},
753        """                        L{TRILINOS}
754        indicates to use matrix lumping         """
755        """         if package==None: package=self.DEFAULT
756        if not self.isUsingLumping():         if not package==self.getSolverPackage():
757           if self.debug() : print "PDE Debug: lumping is set on"             self.__solver_package=package
758           self.__rebuildOperator()             self.updateSystemType()
759           self.__lumping=True             self.trace("New solver is %s"%self.getSolverMethodName())
760    
761     def setLumpingOff(self):     def getSolverPackage(self):
762        """         """
763        switches off matrix lumping         Returns the package of the solver.
       """  
       if self.isUsingLumping():  
          if self.debug() : print "PDE Debug: lumping is set off"  
          self.__rebuildOperator()  
          self.__lumping=False  
764    
765     def setLumping(self,flag=False):         @return: the solver package currently being used
766        """         @rtype: C{int}
767        set the matrix lumping flag to flag         """
768        """         return self.__solver_package
       if flag:  
          self.setLumpingOn()  
       else:  
          self.setLumpingOff()  
769    
770     def isUsingLumping(self):     def isUsingLumping(self):
771        """        """
772                Checks if matrix lumping is the current solver method.
       """  
       return self.__lumping  
   
    #============ method business =========================================================  
    def setSolverMethod(self,solver=util.DEFAULT_METHOD):  
        """  
        sets a new solver  
        """  
        if not solver==self.getSolverMethod():  
            self.__solver_method=solver  
            if self.debug() : print "PDE Debug: New solver is %s"%solver  
            self.__checkMatrixType()  
773    
774     def getSolverMethod(self):        @return: True if the current solver method is lumping
775         """        @rtype: C{bool}
776         returns the solver method        """
777         """        return self.getSolverMethod()[0]==self.LUMPING
        return self.__solver_method  
778    
    #============ tolerance business =========================================================  
779     def setTolerance(self,tol=1.e-8):     def setTolerance(self,tol=1.e-8):
780         """         """
781         resets the tolerance to tol.         Resets the tolerance for the solver method to C{tol}.
782    
783           @param tol: new tolerance for the solver. If C{tol} is lower than the
784                       current tolerence the system will be resolved.
785           @type tol: positive C{float}
786           @raise ValueError: if tolerance is not positive
787         """         """
788         if not tol>0:         if not tol>0:
789             raise ValueException,"Tolerance as to be positive"             raise ValueError,"Tolerance has to be positive"
790         if tol<self.getTolerance(): self.__rebuildSolution()         if tol<self.getTolerance(): self.invalidateSolution()
791         if self.debug() : print "PDE Debug: New tolerance %e",tol         self.trace("New tolerance %e"%tol)
792         self.__tolerance=tol         self.__tolerance=tol
793         return         return
794    
795     def getTolerance(self):     def getTolerance(self):
796         """         """
797         returns the tolerance set for the solution         Returns the tolerance currently set for the solution.
798    
799           @return: tolerance currently used
800           @rtype: C{float}
801         """         """
802         return self.__tolerance         return self.__tolerance
803    
804     #===== symmetry  flag ==========================     # ==========================================================================
805       #    symmetry  flag:
806       # ==========================================================================
807     def isSymmetric(self):     def isSymmetric(self):
808        """        """
809        returns true is the operator is considered to be symmetric        Checks if symmetry is indicated.
810    
811          @return: True if a symmetric PDE is indicated, False otherwise
812          @rtype: C{bool}
813        """        """
814        return self.__sym        return self.__sym
815    
816     def setSymmetryOn(self):     def setSymmetryOn(self):
817        """        """
818        sets the symmetry flag to true        Sets the symmetry flag.
819        """        """
820        if not self.isSymmetric():        if not self.isSymmetric():
821           if self.debug() : print "PDE Debug: Operator is set to be symmetric"           self.trace("PDE is set to be symmetric")
822           self.__sym=True           self.__sym=True
823           self.__checkMatrixType()           self.updateSystemType()
824    
825     def setSymmetryOff(self):     def setSymmetryOff(self):
826        """        """
827        sets the symmetry flag to false        Clears the symmetry flag.
828        """        """
829        if self.isSymmetric():        if self.isSymmetric():
830           if self.debug() : print "PDE Debug: Operator is set to be unsymmetric"           self.trace("PDE is set to be nonsymmetric")
831           self.__sym=False           self.__sym=False
832           self.__checkMatrixType()           self.updateSystemType()
833    
834     def setSymmetryTo(self,flag=False):     def setSymmetryTo(self,flag=False):
835       """        """
836       sets the symmetry flag to flag        Sets the symmetry flag to C{flag}.
837    
838       @param flag:        @param flag: If True, the symmetry flag is set otherwise reset.
839       """        @type flag: C{bool}
840       if flag:        """
841          self.setSymmetryOn()        if flag:
842       else:           self.setSymmetryOn()
843          self.setSymmetryOff()        else:
844             self.setSymmetryOff()
845    
846     #===== order reduction ==========================     # ==========================================================================
847       # function space handling for the equation as well as the solution
848       # ==========================================================================
849     def setReducedOrderOn(self):     def setReducedOrderOn(self):
850       """       """
851       switches to on reduced order       Switches reduced order on for solution and equation representation.
852    
853         @raise RuntimeError: if order reduction is altered after a coefficient has
854                              been set
855       """       """
856       self.setReducedOrderForSolutionOn()       self.setReducedOrderForSolutionOn()
857       self.setReducedOrderForEquationOn()       self.setReducedOrderForEquationOn()
858    
859     def setReducedOrderOff(self):     def setReducedOrderOff(self):
860       """       """
861       switches to full order       Switches reduced order off for solution and equation representation
862    
863         @raise RuntimeError: if order reduction is altered after a coefficient has
864                              been set
865       """       """
866       self.setReducedOrderForSolutionOff()       self.setReducedOrderForSolutionOff()
867       self.setReducedOrderForEquationOff()       self.setReducedOrderForEquationOff()
868    
869     def setReducedOrderTo(self,flag=False):     def setReducedOrderTo(self,flag=False):
870       """       """
871       sets order according to flag       Sets order reduction state for both solution and equation representation
872         according to flag.
873    
874       @param flag:       @param flag: if True, the order reduction is switched on for both solution
875                      and equation representation, otherwise or if flag is not
876                      present order reduction is switched off
877         @type flag: C{bool}
878         @raise RuntimeError: if order reduction is altered after a coefficient has
879                              been set
880       """       """
881       self.setReducedOrderForSolutionTo(flag)       self.setReducedOrderForSolutionTo(flag)
882       self.setReducedOrderForEquationTo(flag)       self.setReducedOrderForEquationTo(flag)
                                                                                                                                                             
883    
884     #===== order reduction solution ==========================  
885     def setReducedOrderForSolutionOn(self):     def setReducedOrderForSolutionOn(self):
886       """       """
887       switches to reduced order to interpolate solution       Switches reduced order on for solution representation.
888    
889         @raise RuntimeError: if order reduction is altered after a coefficient has
890                              been set
891       """       """
892       new_fs=escript.ReducedSolution(self.getDomain())       if not self.__reduce_solution_order:
893       if self.getFunctionSpaceForSolution()!=new_fs:           if self.__altered_coefficients:
894           if self.debug() : print "PDE Debug: Reduced order is used to interpolate solution."                raise RuntimeError,"order cannot be altered after coefficients have been defined."
895           self.__column_function_space=new_fs           self.trace("Reduced order is used for solution representation.")
896           self.__rebuildSystem(deep=True)           self.__reduce_solution_order=True
897             self.initializeSystem()
898    
899     def setReducedOrderForSolutionOff(self):     def setReducedOrderForSolutionOff(self):
900       """       """
901       switches to full order to interpolate solution       Switches reduced order off for solution representation
902    
903         @raise RuntimeError: if order reduction is altered after a coefficient has
904                              been set.
905       """       """
906       new_fs=escript.Solution(self.getDomain())       if self.__reduce_solution_order:
907       if self.getFunctionSpaceForSolution()!=new_fs:           if self.__altered_coefficients:
908           if self.debug() : print "PDE Debug: Full order is used to interpolate solution."                raise RuntimeError,"order cannot be altered after coefficients have been defined."
909           self.__column_function_space=new_fs           self.trace("Full order is used to interpolate solution.")
910           self.__rebuildSystem(deep=True)           self.__reduce_solution_order=False
911             self.initializeSystem()
912    
913     def setReducedOrderForSolutionTo(self,flag=False):     def setReducedOrderForSolutionTo(self,flag=False):
914       """       """
915       sets order for test functions according to flag       Sets order reduction state for solution representation according to flag.
916    
917       @param flag:       @param flag: if flag is True, the order reduction is switched on for
918                      solution representation, otherwise or if flag is not present
919                      order reduction is switched off
920         @type flag: C{bool}
921         @raise RuntimeError: if order reduction is altered after a coefficient has
922                              been set
923       """       """
924       if flag:       if flag:
925          self.setReducedOrderForSolutionOn()          self.setReducedOrderForSolutionOn()
926       else:       else:
927          self.setReducedOrderForSolutionOff()          self.setReducedOrderForSolutionOff()
928                                                                                                                                                              
    #===== order reduction equation ==========================  
929     def setReducedOrderForEquationOn(self):     def setReducedOrderForEquationOn(self):
930       """       """
931       switches to reduced order for test functions       Switches reduced order on for equation representation.
932    
933         @raise RuntimeError: if order reduction is altered after a coefficient has
934                              been set
935       """       """
936       new_fs=escript.ReducedSolution(self.getDomain())       if not self.__reduce_equation_order:
937       if self.getFunctionSpaceForEquation()!=new_fs:           if self.__altered_coefficients:
938           if self.debug() : print "PDE Debug: Reduced order is used for test functions."                raise RuntimeError,"order cannot be altered after coefficients have been defined."
939           self.__row_function_space=new_fs           self.trace("Reduced order is used for test functions.")
940           self.__rebuildSystem(deep=True)           self.__reduce_equation_order=True
941             self.initializeSystem()
942    
943     def setReducedOrderForEquationOff(self):     def setReducedOrderForEquationOff(self):
944       """       """
945       switches to full order for test functions       Switches reduced order off for equation representation.
946    
947         @raise RuntimeError: if order reduction is altered after a coefficient has
948                              been set
949       """       """
950       new_fs=escript.Solution(self.getDomain())       if self.__reduce_equation_order:
951       if self.getFunctionSpaceForEquation()!=new_fs:           if self.__altered_coefficients:
952           if self.debug() : print "PDE Debug: Full order is used for test functions."                raise RuntimeError,"order cannot be altered after coefficients have been defined."
953           self.__row_function_space=new_fs           self.trace("Full order is used for test functions.")
954           self.__rebuildSystem(deep=True)           self.__reduce_equation_order=False
955             self.initializeSystem()
956    
957     def setReducedOrderForEquationTo(self,flag=False):     def setReducedOrderForEquationTo(self,flag=False):
958       """       """
959       sets order for test functions according to flag       Sets order reduction state for equation representation according to flag.
960    
961       @param flag:       @param flag: if flag is True, the order reduction is switched on for
962                      equation representation, otherwise or if flag is not present
963                      order reduction is switched off
964         @type flag: C{bool}
965         @raise RuntimeError: if order reduction is altered after a coefficient has
966                              been set
967       """       """
968       if flag:       if flag:
969          self.setReducedOrderForEquationOn()          self.setReducedOrderForEquationOn()
970       else:       else:
971          self.setReducedOrderForEquationOff()          self.setReducedOrderForEquationOff()
972                                                                                                                                                              
973     # ==== initialization =====================================================================     def updateSystemType(self):
974     def __getNewOperator(self):       """
975         Reassesses the system type and, if a new matrix is needed, resets the
976         system.
977         """
978         new_system_type=self.getRequiredSystemType()
979         if not new_system_type==self.__system_type:
980             self.trace("Matrix type is now %d."%new_system_type)
981             self.__system_type=new_system_type
982             self.initializeSystem()
983    
984       def getSystemType(self):
985          """
986          Returns the current system type.
987          """
988          return self.__system_type
989    
990       def checkSymmetricTensor(self,name,verbose=True):
991          """
992          Tests a coefficient for symmetry.
993    
994          @param name: name of the coefficient
995          @type name: C{str}
996          @param verbose: if set to True or not present a report on coefficients
997                          which break the symmetry is printed.
998          @type verbose: C{bool}
999          @return: True if coefficient C{name} is symmetric
1000          @rtype: C{bool}
1001          """
1002          A=self.getCoefficient(name)
1003          verbose=verbose or self.__debug
1004          out=True
1005          if not A.isEmpty():
1006             tol=util.Lsup(A)*self.SMALL_TOLERANCE
1007             s=A.getShape()
1008             if A.getRank() == 4:
1009                if s[0]==s[2] and s[1] == s[3]:
1010                   for i in range(s[0]):
1011                      for j in range(s[1]):
1012                         for k in range(s[2]):
1013                            for l in range(s[3]):
1014                                if util.Lsup(A[i,j,k,l]-A[k,l,i,j])>tol:
1015                                   if verbose: print "non-symmetric problem as %s[%d,%d,%d,%d]!=%s[%d,%d,%d,%d]"%(name,i,j,k,l,name,k,l,i,j)
1016                                   out=False
1017                else:
1018                     if verbose: print "non-symmetric problem because of inappropriate shape %s of coefficient %s."%(s,name)
1019                     out=False
1020             elif A.getRank() == 2:
1021                if s[0]==s[1]:
1022                   for j in range(s[0]):
1023                      for l in range(s[1]):
1024                         if util.Lsup(A[j,l]-A[l,j])>tol:
1025                            if verbose: print "non-symmetric problem because %s[%d,%d]!=%s[%d,%d]"%(name,j,l,name,l,j)
1026                            out=False
1027                else:
1028                     if verbose: print "non-symmetric problem because of inappropriate shape %s of coefficient %s."%(s,name)
1029                     out=False
1030             elif A.getRank() == 0:
1031                pass
1032             else:
1033                 raise ValueError,"Cannot check rank %s of %s."%(A.getRank(),name)
1034          return out
1035    
1036       def checkReciprocalSymmetry(self,name0,name1,verbose=True):
1037          """
1038          Tests two coefficients for reciprocal symmetry.
1039    
1040          @param name0: name of the first coefficient
1041          @type name0: C{str}
1042          @param name1: name of the second coefficient
1043          @type name1: C{str}
1044          @param verbose: if set to True or not present a report on coefficients
1045                          which break the symmetry is printed
1046          @type verbose: C{bool}
1047          @return: True if coefficients C{name0} and C{name1} are reciprocally
1048                   symmetric.
1049          @rtype: C{bool}
1050          """
1051          B=self.getCoefficient(name0)
1052          C=self.getCoefficient(name1)
1053          verbose=verbose or self.__debug
1054          out=True
1055          if B.isEmpty() and not C.isEmpty():
1056             if verbose: print "non-symmetric problem because %s is not present but %s is"%(name0,name1)
1057             out=False
1058          elif not B.isEmpty() and C.isEmpty():
1059             if verbose: print "non-symmetric problem because %s is not present but %s is"%(name0,name1)
1060             out=False
1061          elif not B.isEmpty() and not C.isEmpty():
1062             sB=B.getShape()
1063             sC=C.getShape()
1064             tol=(util.Lsup(B)+util.Lsup(C))*self.SMALL_TOLERANCE/2.
1065             if len(sB) != len(sC):
1066                 if verbose: print "non-symmetric problem because ranks of %s (=%s) and %s (=%s) are different."%(name0,len(sB),name1,len(sC))
1067                 out=False
1068             else:
1069                 if len(sB)==0:
1070                   if util.Lsup(B-C)>tol:
1071                      if verbose: print "non-symmetric problem because %s!=%s"%(name0,name1)
1072                      out=False
1073                 elif len(sB)==1:
1074                   if sB[0]==sC[0]:
1075                      for j in range(sB[0]):
1076                         if util.Lsup(B[j]-C[j])>tol:
1077                            if verbose: print "non-symmetric PDE because %s[%d]!=%s[%d]"%(name0,j,name1,j)
1078                            out=False
1079                   else:
1080                     if verbose: print "non-symmetric problem because of inappropriate shapes %s and %s of coefficients %s and %s, respectively."%(sB,sC,name0,name1)
1081                 elif len(sB)==3:
1082                   if sB[0]==sC[1] and sB[1]==sC[2] and sB[2]==sC[0]:
1083                       for i in range(sB[0]):
1084                          for j in range(sB[1]):
1085                             for k in range(sB[2]):
1086                                if util.Lsup(B[i,j,k]-C[k,i,j])>tol:
1087                                     if verbose: print "non-symmetric problem because %s[%d,%d,%d]!=%s[%d,%d,%d]"%(name0,i,j,k,name1,k,i,j)
1088                                     out=False
1089                   else:
1090                     if verbose: print "non-symmetric problem because of inappropriate shapes %s and %s of coefficients %s and %s, respectively."%(sB,sC,name0,name1)
1091                 else:
1092                     raise ValueError,"Cannot check rank %s of %s and %s."%(len(sB),name0,name1)
1093          return out
1094    
1095       def getCoefficient(self,name):
1096         """
1097         Returns the value of the coefficient C{name}.
1098    
1099         @param name: name of the coefficient requested
1100         @type name: C{string}
1101         @return: the value of the coefficient
1102         @rtype: L{Data<escript.Data>}
1103         @raise IllegalCoefficient: if C{name} is not a coefficient of the PDE
1104         """
1105         if self.hasCoefficient(name):
1106             return self.__COEFFICIENTS[name].getValue()
1107         else:
1108            raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1109    
1110       def hasCoefficient(self,name):
1111         """
1112         Returns True if C{name} is the name of a coefficient.
1113    
1114         @param name: name of the coefficient enquired
1115         @type name: C{string}
1116         @return: True if C{name} is the name of a coefficient of the general PDE,
1117                  False otherwise
1118         @rtype: C{bool}
1119         """
1120         return self.__COEFFICIENTS.has_key(name)
1121    
1122       def createCoefficient(self, name):
1123         """
1124         Creates a L{Data<escript.Data>} object corresponding to coefficient
1125         C{name}.
1126    
1127         @return: the coefficient C{name} initialized to 0
1128         @rtype: L{Data<escript.Data>}
1129         @raise IllegalCoefficient: if C{name} is not a coefficient of the PDE
1130         """
1131         if self.hasCoefficient(name):
1132            return escript.Data(0.,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))
1133         else:
1134            raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1135    
1136       def getFunctionSpaceForCoefficient(self,name):
1137         """
1138         Returns the L{FunctionSpace<escript.FunctionSpace>} to be used for
1139         coefficient C{name}.
1140    
1141         @param name: name of the coefficient enquired
1142         @type name: C{string}
1143         @return: the function space to be used for coefficient C{name}
1144         @rtype: L{FunctionSpace<escript.FunctionSpace>}
1145         @raise IllegalCoefficient: if C{name} is not a coefficient of the PDE
1146         """
1147         if self.hasCoefficient(name):
1148            return self.__COEFFICIENTS[name].getFunctionSpace(self.getDomain())
1149         else:
1150            raise ValueError,"unknown coefficient %s requested"%name
1151    
1152       def getShapeOfCoefficient(self,name):
1153         """
1154         Returns the shape of the coefficient C{name}.
1155    
1156         @param name: name of the coefficient enquired
1157         @type name: C{string}
1158         @return: the shape of the coefficient C{name}
1159         @rtype: C{tuple} of C{int}
1160         @raise IllegalCoefficient: if C{name} is not a coefficient of the PDE
1161         """
1162         if self.hasCoefficient(name):
1163            return self.__COEFFICIENTS[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())
1164         else:
1165            raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1166    
1167       def resetAllCoefficients(self):
1168         """
1169         Resets all coefficients to their default values.
1170         """
1171         for i in self.__COEFFICIENTS.iterkeys():
1172             self.__COEFFICIENTS[i].resetValue()
1173    
1174       def alteredCoefficient(self,name):
1175         """
1176         Announces that coefficient C{name} has been changed.
1177    
1178         @param name: name of the coefficient affected
1179         @type name: C{string}
1180         @raise IllegalCoefficient: if C{name} is not a coefficient of the PDE
1181         @note: if C{name} is q or r, the method will not trigger a rebuild of the
1182                system as constraints are applied to the solved system.
1183         """
1184         if self.hasCoefficient(name):
1185            self.trace("Coefficient %s has been altered."%name)
1186            if not ((name=="q" or name=="r") and self.isUsingLumping()):
1187               if self.__COEFFICIENTS[name].isAlteringOperator(): self.invalidateOperator()
1188               if self.__COEFFICIENTS[name].isAlteringRightHandSide(): self.invalidateRightHandSide()
1189         else:
1190            raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1191    
1192       def validSolution(self):
1193         """         """
1194           Marks the solution as valid.
1195         """         """
1196         return self.getDomain().newOperator( \         self.__is_solution_valid=True
                            self.getNumEquations(), \  
                            self.getFunctionSpaceForEquation(), \  
                            self.getNumSolutions(), \  
                            self.getFunctionSpaceForSolution(), \  
                            self.__matrix_type)  
1197    
1198     def __makeFreshRightHandSide(self):     def invalidateSolution(self):
1199         """         """
1200           Indicates the PDE has to be resolved if the solution is requested.
1201         """         """
1202         if self.debug() : print "PDE Debug: New right hand side allocated"         self.trace("System will be resolved.")
1203         if self.getNumEquations()>1:         self.__is_solution_valid=False
            self.__righthandside=escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)  
        else:  
            self.__righthandside=escript.Data(0.,(),self.getFunctionSpaceForEquation(),True)  
        return self.__righthandside  
1204    
1205     def __getNewSolution(self):     def isSolutionValid(self):
1206         """         """
1207           Returns True if the solution is still valid.
1208         """         """
1209         if self.debug() : print "PDE Debug: New right hand side allocated"         return self.__is_solution_valid
        if self.getNumSolutions()>1:  
            return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)  
        else:  
            return escript.Data(0.,(),self.getFunctionSpaceForSolution(),True)  
1210    
1211     def __makeFreshOperator(self):     def validOperator(self):
1212         """         """
1213           Marks the operator as valid.
1214         """         """
1215         if self.__operator.isEmpty():         self.__is_operator_valid=True
            self.__operator=self.__getNewOperator()  
            if self.debug() : print "PDE Debug: New operator allocated"  
        else:  
            self.__operator.setValue(0.)  
            self.__operator.resetSolver()  
            if self.debug() : print "PDE Debug: Operator reset to zero"  
        return self.__operator  
1216    
1217     #============ some serivice functions  =====================================================     def invalidateOperator(self):
1218     def getDomain(self):         """
1219       """         Indicates the operator has to be rebuilt next time it is used.
1220       returns the domain of the PDE         """
1221       """         self.trace("Operator will be rebuilt.")
1222       return self.__domain         self.invalidateSolution()
1223           self.__is_operator_valid=False
1224    
1225     def getDim(self):     def isOperatorValid(self):
1226       """         """
1227       returns the spatial dimension of the PDE         Returns True if the operator is still valid.
1228       """         """
1229       return self.getDomain().getDim()         return self.__is_operator_valid
1230    
1231     def getNumEquations(self):     def validRightHandSide(self):
1232           """
1233           Marks the right hand side as valid.
1234           """
1235           self.__is_RHS_valid=True
1236    
1237       def invalidateRightHandSide(self):
1238           """
1239           Indicates the right hand side has to be rebuilt next time it is used.
1240           """
1241           if self.isRightHandSideValid(): self.trace("Right hand side has to be rebuilt.")
1242           self.invalidateSolution()
1243           self.__is_RHS_valid=False
1244    
1245       def isRightHandSideValid(self):
1246           """
1247           Returns True if the operator is still valid.
1248           """
1249           return self.__is_RHS_valid
1250    
1251       def invalidateSystem(self):
1252           """
1253           Announces that everything has to be rebuilt.
1254           """
1255           if self.isRightHandSideValid(): self.trace("System has to be rebuilt.")
1256           self.invalidateSolution()
1257           self.invalidateOperator()
1258           self.invalidateRightHandSide()
1259    
1260       def isSystemValid(self):
1261           """
1262           Returns True if the system (including solution) is still vaild.
1263           """
1264           return self.isSolutionValid() and self.isOperatorValid() and self.isRightHandSideValid()
1265    
1266       def initializeSystem(self):
1267           """
1268           Resets the system clearing the operator, right hand side and solution.
1269           """
1270           self.trace("New System has been created.")
1271           self.__operator=escript.Operator()
1272           self.__righthandside=escript.Data()
1273           self.__solution=escript.Data()
1274           self.invalidateSystem()
1275    
1276       def getOperator(self):
1277       """       """
1278       returns the number of equations       Returns the operator of the linear problem.
1279    
1280         @return: the operator of the problem
1281       """       """
1282       if self.__numEquations>0:       return self.getSystem()[0]
          return self.__numEquations  
      else:  
          raise ValueError,"Number of equations is undefined. Please specify argument numEquations."  
1283    
1284     def getNumSolutions(self):     def getRightHandSide(self):
1285       """       """
1286       returns the number of unknowns       Returns the right hand side of the linear problem.
1287    
1288         @return: the right hand side of the problem
1289         @rtype: L{Data<escript.Data>}
1290       """       """
1291       if self.__numSolutions>0:       return self.getSystem()[1]
         return self.__numSolutions  
      else:  
         raise ValueError,"Number of solution is undefined. Please specify argument numSolutions."  
1292    
1293       def createRightHandSide(self):
1294           """
1295           Returns an instance of a new right hand side.
1296           """
1297           self.trace("New right hand side is allocated.")
1298           if self.getNumEquations()>1:
1299               return escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)
1300           else:
1301               return escript.Data(0.,(),self.getFunctionSpaceForEquation(),True)
1302    
1303     def checkSymmetry(self,verbose=True):     def createSolution(self):
1304        """         """
1305        returns if the Operator is symmetric. This is a very expensive operation!!! The symmetry flag is not altered.         Returns an instance of a new solution.
1306        """         """
1307        verbose=verbose or self.debug()         self.trace("New solution is allocated.")
1308        out=True         if self.getNumSolutions()>1:
1309        if self.getNumSolutions()!=self.getNumEquations():             return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)
1310           if verbose: print "non-symmetric PDE because of different number of equations and solutions"         else:
1311           out=False             return escript.Data(0.,(),self.getFunctionSpaceForSolution(),True)
       else:  
          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  
1312    
1313     def getFlux(self,u):     def resetSolution(self):
1314         """         """
1315         returns the flux J_ij for a given u         Sets the solution to zero.
1316           """
1317           if self.__solution.isEmpty():
1318               self.__solution=self.createSolution()
1319           else:
1320               self.__solution.setToZero()
1321               self.trace("Solution is reset to zero.")
1322    
1323         \f[     def setSolution(self,u):
1324         J_ij=A_{ijkl}u_{k,l}+B_{ijk}u_k-X_{ij}         """
1325         \f]         Sets the solution assuming that makes the system valid.
1326           """
1327           self.__solution=u
1328           self.validSolution()
1329    
1330         @param u: argument of the operator     def getCurrentSolution(self):
1331         """         """
1332         raise SystemError,"getFlux is not implemented yet"         Returns the solution in its current state.
1333         return None         """
1334           if self.__solution.isEmpty(): self.__solution=self.createSolution()
1335           return self.__solution
1336    
1337       def resetRightHandSide(self):
1338           """
1339           Sets the right hand side to zero.
1340           """
1341           if self.__righthandside.isEmpty():
1342               self.__righthandside=self.createRightHandSide()
1343           else:
1344               self.__righthandside.setToZero()
1345               self.trace("Right hand side is reset to zero.")
1346    
1347     def applyOperator(self,u):     def getCurrentRightHandSide(self):
1348         """         """
1349         applies the operator of the PDE to a given solution u in weak from         Returns the right hand side in its current state.
1350           """
1351           if self.__righthandside.isEmpty(): self.__righthandside=self.createRightHandSide()
1352           return self.__righthandside
1353    
1354         @param u: argument of the operator     def resetOperator(self):
1355         """         """
1356         return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())         Makes sure that the operator is instantiated and returns it initialized
1357                                                                                                                                                                     with zeros.
    def getResidual(self,u):  
1358         """         """
1359         return the residual of u in the weak from         if self.__operator.isEmpty():
1360               if self.isUsingLumping():
1361                   self.__operator=self.createSolution()
1362               else:
1363                   self.__operator=self.createOperator()
1364           else:
1365               if self.isUsingLumping():
1366                   self.__operator.setToZero()
1367               else:
1368                   self.__operator.resetValues()
1369               self.trace("Operator reset to zero")
1370    
1371         @param u:     def getCurrentOperator(self):
1372           """
1373           Returns the operator in its current state.
1374         """         """
1375         return self.applyOperator(u)-self.getRightHandSide()         if self.__operator.isEmpty(): self.resetOperator()
1376           return self.__operator
1377    
1378     def __setValue(self,**coefficients):     def setValue(self,**coefficients):
1379        """        """
1380        sets new values to coefficient        Sets new values to coefficients.
1381    
1382        @param coefficients:        @raise IllegalCoefficient: if an unknown coefficient keyword is used
1383        """        """
1384        # check if the coefficients are  legal:        # check if the coefficients are  legal:
1385        for i in coefficients.iterkeys():        for i in coefficients.iterkeys():
1386           if not self.hasCoefficient(i):           if not self.hasCoefficient(i):
1387              raise ValueError,"Attempt to set unknown coefficient %s"%i              raise IllegalCoefficient,"Attempt to set unknown coefficient %s"%i
1388        # 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:
1389        if self.__numEquations<1 or self.__numSolutions<1:        if self.__numEquations==None or self.__numSolutions==None:
1390           for i,d in coefficients.iteritems():           for i,d in coefficients.iteritems():
1391              if hasattr(d,"shape"):              if hasattr(d,"shape"):
1392                  s=d.shape                  s=d.shape
# Line 739  class LinearPDE: Line 1396  class LinearPDE:
1396                  s=numarray.array(d).shape                  s=numarray.array(d).shape
1397              if s!=None:              if s!=None:
1398                  # get number of equations and number of unknowns:                  # get number of equations and number of unknowns:
1399                  res=self.COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(s,self.getDim())                  res=self.__COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(self.getDomain(),s)
1400                  if res==None:                  if res==None:
1401                      raise ValueError,"Illegal shape %s of coefficient %s"%(s,i)                      raise IllegalCoefficientValue,"Illegal shape %s of coefficient %s"%(s,i)
1402                  else:                  else:
1403                      if self.__numEquations<1: self.__numEquations=res[0]                      if self.__numEquations==None: self.__numEquations=res[0]
1404                      if self.__numSolutions<1: self.__numSolutions=res[1]                      if self.__numSolutions==None: self.__numSolutions=res[1]
1405        if self.__numEquations<1: raise ValueError,"unidententified number of equations"        if self.__numEquations==None: raise UndefinedPDEError,"unidentified number of equations"
1406        if self.__numSolutions<1: raise ValueError,"unidententified number of solutions"        if self.__numSolutions==None: raise UndefinedPDEError,"unidentified number of solutions"
1407        # 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:
1408        for i,d in coefficients.iteritems():        for i,d in coefficients.iteritems():
1409          if d==None:          try:
1410               d2=escript.Data()             self.__COEFFICIENTS[i].setValue(self.getDomain(),
1411          elif isinstance(d,escript.Data):                       self.getNumEquations(),self.getNumSolutions(),
1412               if d.isEmpty():                       self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
1413                  d2=d             self.alteredCoefficient(i)
1414               else:          except IllegalCoefficientFunctionSpace,m:
1415                  d2=escript.Data(d,self.getFunctionSpaceForCoefficient(i))              # if the function space is wrong then we try the reduced version:
1416          else:              i_red=i+"_reduced"
1417                d2=escript.Data(d,self.getFunctionSpaceForCoefficient(i))              if (not i_red in coefficients.keys()) and i_red in self.__COEFFICIENTS.keys():
1418          if not d2.isEmpty():                  try:
1419             if not self.getShapeOfCoefficient(i)==d2.getShape():                      self.__COEFFICIENTS[i_red].setValue(self.getDomain(),
1420                 raise ValueError,"Expected shape for coefficient %s is %s but actual shape is %s."%(i,self.getShapeOfCoefficient(i),d2.getShape())                                                        self.getNumEquations(),self.getNumSolutions(),
1421          # overwrite new values:                                                        self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
1422          if self.debug(): print "PDE Debug: Coefficient %s has been altered."%i                      self.alteredCoefficient(i_red)
1423          self.COEFFICIENTS[i].setValue(d2)                  except IllegalCoefficientValue,m:
1424          self.alteredCoefficient(i)                      raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
1425                          except IllegalCoefficientFunctionSpace,m:
1426        # reset the HomogeneousConstraintFlag:                      raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m))
1427        self.__setHomogeneousConstraintFlag()              else:
1428        if len(coefficients)>0 and not self.isUsingLumping() and not self.__homogeneous_constraint: self.__rebuildSystem()                  raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m))
1429            except IllegalCoefficientValue,m:
1430     def __setHomogeneousConstraintFlag(self):             raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
1431        """        self.__altered_coefficients=True
       checks if the constraints are homogeneous and sets self.__homogeneous_constraint accordingly.  
       """  
       self.__homogeneous_constraint=True  
       q=self.getCoefficientOfPDE("q")  
       r=self.getCoefficientOfPDE("r")  
       if not q.isEmpty() and not r.isEmpty():  
          if (q*r).Lsup()>=1.e-13*r.Lsup(): self.__homogeneous_constraint=False  
       if self.debug():  
            if self.__homogeneous_constraint:  
                print "PDE Debug: Constraints are homogeneous."  
            else:  
                print "PDE Debug: Constraints are inhomogeneous."  
   
1432    
1433     # ==== rebuild switches =====================================================================     # ==========================================================================
1434     def __rebuildSolution(self,deep=False):     # methods that are typically overwritten when implementing a particular
1435         """     # linear problem
1436         indicates the PDE has to be reolved if the solution is requested     # ==========================================================================
1437         """     def getRequiredSystemType(self):
1438         if self.__solution_isValid and self.debug() : print "PDE Debug: PDE has to be resolved."        """
1439         self.__solution_isValid=False        Returns the system type which needs to be used by the current set up.
        if deep: self.__solution=escript.Data()  
1440    
1441          @note: Typically this method is overwritten when implementing a
1442                 particular linear problem.
1443          """
1444          return 0
1445    
1446     def __rebuildOperator(self,deep=False):     def createOperator(self):
1447         """         """
1448         indicates the operator has to be rebuilt next time it is used         Returns an instance of a new operator.
1449    
1450           @note: This method is overwritten when implementing a particular
1451                  linear problem.
1452         """         """
1453         if self.__operator_isValid and self.debug() : print "PDE Debug: Operator has to be rebuilt."         return escript.Operator()
        self.__rebuildSolution(deep)  
        self.__operator_isValid=False  
        if deep: self.__operator=escript.Operator()  
1454    
1455     def __rebuildRightHandSide(self,deep=False):     def checkSymmetry(self,verbose=True):
1456          """
1457          Tests the PDE for symmetry.
1458    
1459          @param verbose: if set to True or not present a report on coefficients
1460                          which break the symmetry is printed
1461          @type verbose: C{bool}
1462          @return: True if the problem is symmetric
1463          @rtype: C{bool}
1464          @note: Typically this method is overwritten when implementing a
1465                 particular linear problem.
1466          """
1467          out=True
1468          return out
1469    
1470       def getSolution(self,**options):
1471         """         """
1472         indicates the right hand side has to be rebuild next time it is used         Returns the solution of the problem.
1473    
1474           @return: the solution
1475           @rtype: L{Data<escript.Data>}
1476    
1477           @note: This method is overwritten when implementing a particular
1478                  linear problem.
1479         """         """
1480         if self.__righthandside_isValid and self.debug() : print "PDE Debug: Right hand side has to be rebuilt."         return self.getCurrentSolution()
        self.__rebuildSolution(deep)  
        self.__righthandside_isValid=False  
        if deep: self.__righthandside=escript.Data()  
1481    
1482     def __rebuildSystem(self,deep=False):     def getSystem(self):
1483         """         """
1484         annonced that all coefficient name has been changed         Returns the operator and right hand side of the PDE.
1485    
1486           @return: the discrete version of the PDE
1487           @rtype: C{tuple} of L{Operator,<escript.Operator>} and L{Data<escript.Data>}.
1488    
1489           @note: This method is overwritten when implementing a particular
1490                  linear problem.
1491         """         """
1492         self.__rebuildSolution(deep)         return (self.getCurrentOperator(), self.getCurrentRightHandSide())
1493         self.__rebuildOperator(deep)  #=====================
1494         self.__rebuildRightHandSide(deep)  
1495      class LinearPDE(LinearProblem):
1496     def __checkMatrixType(self):     """
1497       This class is used to define a general linear, steady, second order PDE
1498       for an unknown function M{u} on a given domain defined through a
1499       L{Domain<escript.Domain>} object.
1500    
1501       For a single PDE having a solution with a single component the linear PDE
1502       is defined in the following form:
1503    
1504       M{-(grad(A[j,l]+A_reduced[j,l])*grad(u)[l]+(B[j]+B_reduced[j])u)[j]+(C[l]+C_reduced[l])*grad(u)[l]+(D+D_reduced)=-grad(X+X_reduced)[j,j]+(Y+Y_reduced)}
1505    
1506       where M{grad(F)} denotes the spatial derivative of M{F}. Einstein's
1507       summation convention, ie. summation over indexes appearing twice in a term
1508       of a sum performed, is used.
1509       The coefficients M{A}, M{B}, M{C}, M{D}, M{X} and M{Y} have to be specified
1510       through L{Data<escript.Data>} objects in L{Function<escript.Function>} and
1511       the coefficients M{A_reduced}, M{B_reduced}, M{C_reduced}, M{D_reduced},
1512       M{X_reduced} and M{Y_reduced} have to be specified through
1513       L{Data<escript.Data>} objects in L{ReducedFunction<escript.ReducedFunction>}.
1514       It is also allowed to use objects that can be converted into such
1515       L{Data<escript.Data>} objects. M{A} and M{A_reduced} are rank two, M{B},
1516       M{C}, M{X}, M{B_reduced}, M{C_reduced} and M{X_reduced} are rank one and
1517       M{D}, M{D_reduced}, M{Y} and M{Y_reduced} are scalar.
1518    
1519       The following natural boundary conditions are considered:
1520    
1521       M{n[j]*((A[i,j]+A_reduced[i,j])*grad(u)[l]+(B+B_reduced)[j]*u)+(d+d_reduced)*u=n[j]*(X[j]+X_reduced[j])+y}
1522    
1523       where M{n} is the outer normal field. Notice that the coefficients M{A},
1524       M{A_reduced}, M{B}, M{B_reduced}, M{X} and M{X_reduced} are defined in the
1525       PDE. The coefficients M{d} and M{y} are each a scalar in
1526       L{FunctionOnBoundary<escript.FunctionOnBoundary>} and the coefficients
1527       M{d_reduced} and M{y_reduced} are each a scalar in
1528       L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
1529    
1530       Constraints for the solution prescribe the value of the solution at certain
1531       locations in the domain. They have the form
1532    
1533       M{u=r} where M{q>0}
1534    
1535       M{r} and M{q} are each scalar where M{q} is the characteristic function
1536       defining where the constraint is applied. The constraints override any
1537       other condition set by the PDE or the boundary condition.
1538    
1539       The PDE is symmetrical if
1540    
1541       M{A[i,j]=A[j,i]}  and M{B[j]=C[j]} and M{A_reduced[i,j]=A_reduced[j,i]}
1542       and M{B_reduced[j]=C_reduced[j]}
1543    
1544       For a system of PDEs and a solution with several components the PDE has the
1545       form
1546    
1547       M{-grad((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k])[j]+(C[i,k,l]+C_reduced[i,k,l])*grad(u[k])[l]+(D[i,k]+D_reduced[i,k]*u[k] =-grad(X[i,j]+X_reduced[i,j])[j]+Y[i]+Y_reduced[i] }
1548    
1549       M{A} and M{A_reduced} are of rank four, M{B}, M{B_reduced}, M{C} and
1550       M{C_reduced} are each of rank three, M{D}, M{D_reduced}, M{X_reduced} and
1551       M{X} are each of rank two and M{Y} and M{Y_reduced} are of rank one.
1552       The natural boundary conditions take the form:
1553    
1554       M{n[j]*((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k])+(d[i,k]+d_reduced[i,k])*u[k]=n[j]*(X[i,j]+X_reduced[i,j])+y[i]+y_reduced[i]}
1555    
1556       The coefficient M{d} is of rank two and M{y} is of rank one both in
1557       L{FunctionOnBoundary<escript.FunctionOnBoundary>}. The coefficients
1558       M{d_reduced} is of rank two and M{y_reduced} is of rank one both in
1559       L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
1560    
1561       Constraints take the form
1562    
1563       M{u[i]=r[i]}  where  M{q[i]>0}
1564    
1565       M{r} and M{q} are each rank one. Notice that at some locations not
1566       necessarily all components must have a constraint.
1567    
1568       The system of PDEs is symmetrical if
1569    
1570          - M{A[i,j,k,l]=A[k,l,i,j]}
1571          - M{A_reduced[i,j,k,l]=A_reduced[k,l,i,j]}
1572          - M{B[i,j,k]=C[k,i,j]}
1573          - M{B_reduced[i,j,k]=C_reduced[k,i,j]}
1574          - M{D[i,k]=D[i,k]}
1575          - M{D_reduced[i,k]=D_reduced[i,k]}
1576          - M{d[i,k]=d[k,i]}
1577          - M{d_reduced[i,k]=d_reduced[k,i]}
1578    
1579       L{LinearPDE} also supports solution discontinuities over a contact region
1580       in the domain. To specify the conditions across the discontinuity we are
1581       using the generalised flux M{J} which, in the case of a system of PDEs
1582       and several components of the solution, is defined as
1583    
1584       M{J[i,j]=(A[i,j,k,l]+A_reduced[[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k]-X[i,j]-X_reduced[i,j]}
1585    
1586       For the case of single solution component and single PDE M{J} is defined as
1587    
1588       M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[j]+(B[i]+B_reduced[i])*u-X[i]-X_reduced[i]}
1589    
1590       In the context of discontinuities M{n} denotes the normal on the
1591       discontinuity pointing from side 0 towards side 1 calculated from
1592       L{getNormal<escript.FunctionSpace.getNormal>} of L{FunctionOnContactZero<escript.FunctionOnContactZero>}.
1593       For a system of PDEs the contact condition takes the form
1594    
1595       M{n[j]*J0[i,j]=n[j]*J1[i,j]=(y_contact[i]+y_contact_reduced[i])- (d_contact[i,k]+d_contact_reduced[i,k])*jump(u)[k]}
1596    
1597       where M{J0} and M{J1} are the fluxes on side 0 and side 1 of the
1598       discontinuity, respectively. M{jump(u)}, which is the difference of the
1599       solution at side 1 and at side 0, denotes the jump of M{u} across
1600       discontinuity along the normal calculated by L{jump<util.jump>}.
1601       The coefficient M{d_contact} is of rank two and M{y_contact} is of rank one
1602       both in L{FunctionOnContactZero<escript.FunctionOnContactZero>} or
1603       L{FunctionOnContactOne<escript.FunctionOnContactOne>}.
1604       The coefficient M{d_contact_reduced} is of rank two and M{y_contact_reduced}
1605       is of rank one both in L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>}
1606       or L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.
1607       In case of a single PDE and a single component solution the contact
1608       condition takes the form
1609    
1610       M{n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)}
1611    
1612       In this case the coefficient M{d_contact} and M{y_contact} are each scalar
1613       both in L{FunctionOnContactZero<escript.FunctionOnContactZero>} or
1614       L{FunctionOnContactOne<escript.FunctionOnContactOne>} and the coefficient
1615       M{d_contact_reduced} and M{y_contact_reduced} are each scalar both in
1616       L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or
1617       L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.
1618    
1619       Typical usage::
1620    
1621           p = LinearPDE(dom)
1622           p.setValue(A=kronecker(dom), D=1, Y=0.5)
1623           u = p.getSolution()
1624    
1625       """
1626    
1627       def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):
1628       """       """
1629       reassess the matrix type and, if needed, initiates an operator rebuild       Initializes a new linear PDE.
1630    
1631         @param domain: domain of the PDE
1632         @type domain: L{Domain<escript.Domain>}
1633         @param numEquations: number of equations. If C{None} the number of
1634                              equations is extracted from the PDE coefficients.
1635         @param numSolutions: number of solution components. If C{None} the number
1636                              of solution components is extracted from the PDE
1637                              coefficients.
1638         @param debug: if True debug information is printed
1639    
1640         """
1641         super(LinearPDE, self).__init__(domain,numEquations,numSolutions,debug)
1642         #
1643         #   the coefficients of the PDE:
1644         #
1645         self.introduceCoefficients(
1646           A=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
1647           B=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1648           C=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
1649           D=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1650           X=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
1651           Y=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1652           d=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1653           y=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1654           d_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1655           y_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1656           A_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
1657           B_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1658           C_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
1659           D_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1660           X_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
1661           Y_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1662           d_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1663           y_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1664           d_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1665           y_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1666           r=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.RIGHTHANDSIDE),
1667           q=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.BOTH) )
1668    
1669       def __str__(self):
1670         """
1671         Returns the string representation of the PDE.
1672    
1673         @return: a simple representation of the PDE
1674         @rtype: C{str}
1675       """       """
1676       new_matrix_type=self.getDomain().getSystemMatrixTypeId(self.getSolverMethod(),self.isSymmetric())       return "<LinearPDE %d>"%id(self)
      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)  
1677    
1678     #============ assembling =======================================================     def getRequiredSystemType(self):
    def __copyConstraint(self):  
1679        """        """
1680        copies the constrint condition into u        Returns the system type which needs to be used by the current set up.
1681        """        """
1682        if not self.__righthandside.isEmpty():        return self.getDomain().getSystemMatrixTypeId(self.getSolverMethod()[0],self.getSolverMethod()[1],self.getSolverPackage(),self.isSymmetric())
          q=self.getCoefficientOfPDE("q")  
          r=self.getCoefficientOfPDE("r")  
          if not q.isEmpty():  
              if r.isEmpty():  
                 r2=escript.Data(0,self.__righthandside.getShape(),self.__righthandside.getFunctionSpace())  
              else:  
                 r2=escript.Data(r,self.__righthandside.getFunctionSpace())  
              self.__righthandside.copyWithMask(r2,escript.Data(q,self.__righthandside.getFunctionSpace()))  
1683    
1684     def __applyConstraint(self):     def checkSymmetry(self,verbose=True):
1685          """
1686          Tests the PDE for symmetry.
1687    
1688          @param verbose: if set to True or not present a report on coefficients
1689                          which break the symmetry is printed.
1690          @type verbose: C{bool}
1691          @return: True if the PDE is symmetric
1692          @rtype: L{bool}
1693          @note: This is a very expensive operation. It should be used for
1694                 degugging only! The symmetry flag is not altered.
1695          """
1696          out=True
1697          out=out and self.checkSymmetricTensor("A", verbose)
1698          out=out and self.checkSymmetricTensor("A_reduced", verbose)
1699          out=out and self.checkReciprocalSymmetry("B","C", verbose)
1700          out=out and self.checkReciprocalSymmetry("B_reduced","C_reduced", verbose)
1701          out=out and self.checkSymmetricTensor("D", verbose)
1702          out=out and self.checkSymmetricTensor("D_reduced", verbose)
1703          out=out and self.checkSymmetricTensor("d", verbose)
1704          out=out and self.checkSymmetricTensor("d_reduced", verbose)
1705          out=out and self.checkSymmetricTensor("d_contact", verbose)
1706          out=out and self.checkSymmetricTensor("d_contact_reduced", verbose)
1707          return out
1708    
1709       def createOperator(self):
1710         """         """
1711         applies the constraints defined by q and r to the system         Returns an instance of a new operator.
1712         """         """
1713         q=self.getCoefficientOfPDE("q")         self.trace("New operator is allocated.")
1714         r=self.getCoefficientOfPDE("r")         return self.getDomain().newOperator( \
1715         if not q.isEmpty() and not self.__operator.isEmpty():                             self.getNumEquations(), \
1716            # q is the row and column mask to indicate where constraints are set:                             self.getFunctionSpaceForEquation(), \
1717            row_q=escript.Data(q,self.getFunctionSpaceForEquation())                             self.getNumSolutions(), \
1718            col_q=escript.Data(q,self.getFunctionSpaceForSolution())                             self.getFunctionSpaceForSolution(), \
1719            u=self.__getNewSolution()                             self.getSystemType())
1720            if r.isEmpty():  
1721               r_s=self.__getNewSolution()     def getSolution(self,**options):
1722            else:         """
1723               r_s=escript.Data(r,self.getFunctionSpaceForSolution())         Returns the solution of the PDE.
1724            u.copyWithMask(r_s,col_q)  
1725           @return: the solution
1726           @rtype: L{Data<escript.Data>}
1727           @param options: solver options
1728           @keyword verbose: True to get some information during PDE solution
1729           @type verbose: C{bool}
1730           @keyword reordering: reordering scheme to be used during elimination.
1731                                Allowed values are L{NO_REORDERING},
1732                                L{MINIMUM_FILL_IN} and L{NESTED_DISSECTION}
1733           @keyword iter_max: maximum number of iteration steps allowed
1734           @keyword drop_tolerance: threshold for dropping in L{ILUT}
1735           @keyword drop_storage: maximum of allowed memory in L{ILUT}
1736           @keyword truncation: maximum number of residuals in L{GMRES}
1737           @keyword restart: restart cycle length in L{GMRES}
1738           """
1739           if not self.isSolutionValid():
1740              mat,f=self.getSystem()
1741            if self.isUsingLumping():            if self.isUsingLumping():
1742               self.__operator.copyWithMask(escript.Data(1,q.getShape(),self.getFunctionSpaceForEquation()),row_q)               self.setSolution(f*1/mat)
1743            else:            else:
1744               if not self.__righthandside.isEmpty(): self.__righthandside-=self.__operator*u               options[self.TOLERANCE_KEY]=self.getTolerance()
1745               self.__operator.nullifyRowsAndCols(row_q,col_q,1.)               options[self.METHOD_KEY]=self.getSolverMethod()[0]
1746                 options[self.PRECONDITIONER_KEY]=self.getSolverMethod()[1]
1747                 options[self.PACKAGE_KEY]=self.getSolverPackage()
1748                 options[self.SYMMETRY_KEY]=self.isSymmetric()
1749                 self.trace("PDE is resolved.")
1750                 self.trace("solver options: %s"%str(options))
1751                 self.setSolution(mat.solve(f,options))
1752           return self.getCurrentSolution()
1753    
1754     def getSystem(self):     def getSystem(self):
1755         """         """
1756         return the operator and right hand side of the PDE         Returns the operator and right hand side of the PDE.
1757    
1758           @return: the discrete version of the PDE
1759           @rtype: C{tuple} of L{Operator,<escript.Operator>} and
1760                   L{Data<escript.Data>}
1761         """         """
1762         if not self.__operator_isValid or not self.__righthandside_isValid:         if not self.isOperatorValid() or not self.isRightHandSideValid():
1763            if self.isUsingLumping():            if self.isUsingLumping():
1764                if not self.__operator_isValid:                if not self.isOperatorValid():
1765                   if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():                   if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():
1766                         raise TypeError,"Lumped matrix requires same order for equations and unknowns"                        raise TypeError,"Lumped matrix requires same order for equations and unknowns"
1767                   if not self.getCoefficientOfPDE("A").isEmpty():                   if not self.getCoefficient("A").isEmpty():
1768                            raise Warning,"Lumped matrix does not allow coefficient A"                        raise ValueError,"coefficient A in lumped matrix may not be present."
1769                   if not self.getCoefficientOfPDE("B").isEmpty():                   if not self.getCoefficient("B").isEmpty():
1770                            raise Warning,"Lumped matrix does not allow coefficient B"                        raise ValueError,"coefficient B in lumped matrix may not be present."
1771                   if not self.getCoefficientOfPDE("C").isEmpty():                   if not self.getCoefficient("C").isEmpty():
1772                            raise Warning,"Lumped matrix does not allow coefficient C"                        raise ValueError,"coefficient C in lumped matrix may not be present."
1773                   if self.debug() : print "PDE Debug: New lumped operator is built."                   if not self.getCoefficient("d_contact").isEmpty():
1774                   mat=self.__getNewOperator()                        raise ValueError,"coefficient d_contact in lumped matrix may not be present."
1775                   self.getDomain().addPDEToSystem(mat,escript.Data(), \                   if not self.getCoefficient("A_reduced").isEmpty():
1776                             self.getCoefficientOfPDE("A"), \                        raise ValueError,"coefficient A_reduced in lumped matrix may not be present."
1777                             self.getCoefficientOfPDE("B"), \                   if not self.getCoefficient("B_reduced").isEmpty():
1778                             self.getCoefficientOfPDE("C"), \                        raise ValueError,"coefficient B_reduced in lumped matrix may not be present."
1779                             self.getCoefficientOfPDE("D"), \                   if not self.getCoefficient("C_reduced").isEmpty():
1780                             escript.Data(), \                        raise ValueError,"coefficient C_reduced in lumped matrix may not be present."
1781                             escript.Data(), \                   if not self.getCoefficient("d_contact_reduced").isEmpty():
1782                             self.getCoefficientOfPDE("d"), \                        raise ValueError,"coefficient d_contact_reduced in lumped matrix may not be present."
1783                             escript.Data(),\                   D=self.getCoefficient("D")
1784                             self.getCoefficientOfPDE("d_contact"), \                   d=self.getCoefficient("d")
1785                             escript.Data())                   D_reduced=self.getCoefficient("D_reduced")
1786                   self.__operator=mat*escript.Data(1,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)                   d_reduced=self.getCoefficient("d_reduced")
1787                   self.__applyConstraint()                   if not D.isEmpty():
1788                   self.__operator_isValid=True                       if self.getNumSolutions()>1:
1789                if not self.__righthandside_isValid:                          D_times_e=util.matrix_mult(D,numarray.ones((self.getNumSolutions(),)))
1790                   if self.debug() : print "PDE Debug: New right hand side is built."                       else:
1791                   self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \                          D_times_e=D
1792                                 self.getCoefficientOfPDE("X"), \                   else:
1793                                 self.getCoefficientOfPDE("Y"),\                      D_times_e=escript.Data()
1794                                 self.getCoefficientOfPDE("y"),\                   if not d.isEmpty():
1795                                 self.getCoefficientOfPDE("y_contact"))                       if self.getNumSolutions()>1:
1796                   self.__copyConstraint()                          d_times_e=util.matrix_mult(d,numarray.ones((self.getNumSolutions(),)))
1797                   self.__righthandside_isValid=True                       else:
1798                            d_times_e=d
1799                     else:
1800                        d_times_e=escript.Data()
1801    
1802                     if not D_reduced.isEmpty():
1803                         if self.getNumSolutions()>1:
1804                            D_reduced_times_e=util.matrix_mult(D_reduced,numarray.ones((self.getNumSolutions(),)))
1805                         else:
1806                            D_reduced_times_e=D_reduced
1807                     else:
1808                        D_reduced_times_e=escript.Data()
1809                     if not d_reduced.isEmpty():
1810                         if self.getNumSolutions()>1:
1811                            d_reduced_times_e=util.matrix_mult(d_reduced,numarray.ones((self.getNumSolutions(),)))
1812                         else:
1813                            d_reduced_times_e=d_reduced
1814                     else:
1815                        d_reduced_times_e=escript.Data()
1816    
1817                     self.resetOperator()
1818                     operator=self.getCurrentOperator()
1819                     if False and hasattr(self.getDomain(), "addPDEToLumpedSystem") :
1820                        self.getDomain().addPDEToLumpedSystem(operator, D_times_e, d_times_e)
1821                        self.getDomain().addPDEToLumpedSystem(operator, D_reduced_times_e, d_reduced_times_e)
1822                     else:
1823                        self.getDomain().addPDEToRHS(operator, \
1824                                                     escript.Data(), \
1825                                                     D_times_e, \
1826                                                     d_times_e,\
1827                                                     escript.Data())
1828                        self.getDomain().addPDEToRHS(operator, \
1829                                                     escript.Data(), \
1830                                                     D_reduced_times_e, \
1831                                                     d_reduced_times_e,\
1832                                                     escript.Data())
1833                     self.trace("New lumped operator has been built.")
1834                  if not self.isRightHandSideValid():
1835                     self.resetRightHandSide()
1836                     righthandside=self.getCurrentRightHandSide()
1837                     self.getDomain().addPDEToRHS(righthandside, \
1838                                   self.getCoefficient("X"), \
1839                                   self.getCoefficient("Y"),\
1840                                   self.getCoefficient("y"),\
1841                                   self.getCoefficient("y_contact"))
1842                     self.getDomain().addPDEToRHS(righthandside, \
1843                                   self.getCoefficient("X_reduced"), \
1844                                   self.getCoefficient("Y_reduced"),\
1845                                   self.getCoefficient("y_reduced"),\
1846                                   self.getCoefficient("y_contact_reduced"))
1847                     self.trace("New right hand side has been built.")
1848                     self.validRightHandSide()
1849                  self.insertConstraint(rhs_only=False)
1850                  self.validOperator()
1851            else:            else:
1852               if not self.__operator_isValid and not self.__righthandside_isValid:               if not self.isOperatorValid() and not self.isRightHandSideValid():
1853                   if self.debug() : print "PDE Debug: New system is built."                   self.resetRightHandSide()
1854                   self.getDomain().addPDEToSystem(self.__makeFreshOperator(),self.__makeFreshRightHandSide(), \                   righthandside=self.getCurrentRightHandSide()
1855                                 self.getCoefficientOfPDE("A"), \                   self.resetOperator()
1856                                 self.getCoefficientOfPDE("B"), \                   operator=self.getCurrentOperator()
1857                                 self.getCoefficientOfPDE("C"), \                   self.getDomain().addPDEToSystem(operator,righthandside, \
1858                                 self.getCoefficientOfPDE("D"), \                                 self.getCoefficient("A"), \
1859                                 self.getCoefficientOfPDE("X"), \                                 self.getCoefficient("B"), \
1860                                 self.getCoefficientOfPDE("Y"), \                                 self.getCoefficient("C"), \
1861                                 self.getCoefficientOfPDE("d"), \                                 self.getCoefficient("D"), \
1862                                 self.getCoefficientOfPDE("y"), \                                 self.getCoefficient("X"), \
1863                                 self.getCoefficientOfPDE("d_contact"), \                                 self.getCoefficient("Y"), \
1864                                 self.getCoefficientOfPDE("y_contact"))                                 self.getCoefficient("d"), \
1865                   self.__applyConstraint()                                 self.getCoefficient("y"), \
1866                   self.__copyConstraint()                                 self.getCoefficient("d_contact"), \
1867                   self.__operator_isValid=True                                 self.getCoefficient("y_contact"))
1868                   self.__righthandside_isValid=True                   self.getDomain().addPDEToSystem(operator,righthandside, \
1869               elif not self.__righthandside_isValid:                                 self.getCoefficient("A_reduced"), \
1870                   if self.debug() : print "PDE Debug: New right hand side is built."                                 self.getCoefficient("B_reduced"), \
1871                   self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \                                 self.getCoefficient("C_reduced"), \
1872                                 self.getCoefficientOfPDE("X"), \                                 self.getCoefficient("D_reduced"), \
1873                                 self.getCoefficientOfPDE("Y"),\                                 self.getCoefficient("X_reduced"), \
1874                                 self.getCoefficientOfPDE("y"),\                                 self.getCoefficient("Y_reduced"), \
1875                                 self.getCoefficientOfPDE("y_contact"))                                 self.getCoefficient("d_reduced"), \
1876                   self.__copyConstraint()                                 self.getCoefficient("y_reduced"), \
1877                   self.__righthandside_isValid=True                                 self.getCoefficient("d_contact_reduced"), \
1878               elif not self.__operator_isValid:                                 self.getCoefficient("y_contact_reduced"))
1879                   if self.debug() : print "PDE Debug: New operator is built."                   self.insertConstraint(rhs_only=False)
1880                   self.getDomain().addPDEToSystem(self.__makeFreshOperator(),escript.Data(), \                   self.trace("New system has been built.")
1881                              self.getCoefficientOfPDE("A"), \                   self.validOperator()
1882                              self.getCoefficientOfPDE("B"), \                   self.validRightHandSide()
1883                              self.getCoefficientOfPDE("C"), \               elif not self.isRightHandSideValid():
1884                              self.getCoefficientOfPDE("D"), \                   self.resetRightHandSide()
1885                     righthandside=self.getCurrentRightHandSide()
1886                     self.getDomain().addPDEToRHS(righthandside,
1887                                   self.getCoefficient("X"), \
1888                                   self.getCoefficient("Y"),\
1889                                   self.getCoefficient("y"),\
1890                                   self.getCoefficient("y_contact"))
1891                     self.getDomain().addPDEToRHS(righthandside,
1892                                   self.getCoefficient("X_reduced"), \
1893                                   self.getCoefficient("Y_reduced"),\
1894                                   self.getCoefficient("y_reduced"),\
1895                                   self.getCoefficient("y_contact_reduced"))
1896                     self.insertConstraint(rhs_only=True)
1897                     self.trace("New right hand side has been built.")
1898                     self.validRightHandSide()
1899                 elif not self.isOperatorValid():
1900                     self.resetOperator()
1901                     operator=self.getCurrentOperator()
1902                     self.getDomain().addPDEToSystem(operator,escript.Data(), \
1903                                self.getCoefficient("A"), \
1904                                self.getCoefficient("B"), \
1905                                self.getCoefficient("C"), \
1906                                self.getCoefficient("D"), \
1907                              escript.Data(), \                              escript.Data(), \
1908                              escript.Data(), \                              escript.Data(), \
1909                              self.getCoefficientOfPDE("d"), \                              self.getCoefficient("d"), \
1910                              escript.Data(),\                              escript.Data(),\
1911                              self.getCoefficientOfPDE("d_contact"), \                              self.getCoefficient("d_contact"), \
1912                              escript.Data())                              escript.Data())
1913                   self.__applyConstraint()                   self.getDomain().addPDEToSystem(operator,escript.Data(), \
1914                   self.__operator_isValid=True                              self.getCoefficient("A_reduced"), \
1915         return (self.__operator,self.__righthandside)                              self.getCoefficient("B_reduced"), \
1916     def getOperator(self):                              self.getCoefficient("C_reduced"), \
1917         """                              self.getCoefficient("D_reduced"), \
1918         returns the operator of the PDE                              escript.Data(), \
1919         """                              escript.Data(), \
1920         return self.getSystem()[0]                              self.getCoefficient("d_reduced"), \
1921                                escript.Data(),\
1922     def getRightHandSide(self):                              self.getCoefficient("d_contact_reduced"), \
1923         """                              escript.Data())
1924         returns the right hand side of the PDE                   self.insertConstraint(rhs_only=False)
1925         """                   self.trace("New operator has been built.")
1926         return self.getSystem()[1]                   self.validOperator()
1927           return (self.getCurrentOperator(), self.getCurrentRightHandSide())
1928    
1929       def insertConstraint(self, rhs_only=False):
1930          """
1931          Applies the constraints defined by q and r to the PDE.
1932    
1933          @param rhs_only: if True only the right hand side is altered by the
1934                           constraint
1935          @type rhs_only: C{bool}
1936          """
1937          q=self.getCoefficient("q")
1938          r=self.getCoefficient("r")
1939          righthandside=self.getCurrentRightHandSide()
1940          operator=self.getCurrentOperator()
1941    
1942          if not q.isEmpty():
1943             if r.isEmpty():
1944                r_s=self.createSolution()
1945             else:
1946                r_s=r
1947             if not rhs_only and not operator.isEmpty():
1948                 if self.isUsingLumping():
1949                     operator.copyWithMask(escript.Data(1.,q.getShape(),q.getFunctionSpace()),q)
1950                 else:
1951                     row_q=escript.Data(q,self.getFunctionSpaceForEquation())
1952                     col_q=escript.Data(q,self.getFunctionSpaceForSolution())
1953                     u=self.createSolution()
1954                     u.copyWithMask(r_s,col_q)
1955                     righthandside-=operator*u
1956                     operator.nullifyRowsAndCols(row_q,col_q,1.)
1957             righthandside.copyWithMask(r_s,q)
1958    
1959     def solve(self,**options):     def setValue(self,**coefficients):
1960        """        """
1961        solve the PDE        Sets new values to coefficients.
1962    
1963        @param options:        @param coefficients: new values assigned to coefficients
1964        """        @keyword A: value for coefficient C{A}
1965        mat,f=self.getSystem()        @type A: any type that can be cast to a L{Data<escript.Data>} object on
1966        if self.isUsingLumping():                 L{Function<escript.Function>}
1967           out=f/mat        @keyword A_reduced: value for coefficient C{A_reduced}
1968        else:        @type A_reduced: any type that can be cast to a L{Data<escript.Data>}
1969           options[util.TOLERANCE_KEY]=self.getTolerance()                         object on L{ReducedFunction<escript.ReducedFunction>}
1970           options[util.METHOD_KEY]=self.getSolverMethod()        @keyword B: value for coefficient C{B}
1971           options[util.SYMMETRY_KEY]=self.isSymmetric()        @type B: any type that can be cast to a L{Data<escript.Data>} object on
1972           if self.debug() : print "PDE Debug: solver options: ",options                 L{Function<escript.Function>}
1973           out=mat.solve(f,options)        @keyword B_reduced: value for coefficient C{B_reduced}
1974        return out        @type B_reduced: any type that can be cast to a L{Data<escript.Data>}
1975                           object on L{ReducedFunction<escript.ReducedFunction>}
1976          @keyword C: value for coefficient C{C}
1977          @type C: any type that can be cast to a L{Data<escript.Data>} object on
1978                   L{Function<escript.Function>}
1979          @keyword C_reduced: value for coefficient C{C_reduced}
1980          @type C_reduced: any type that can be cast to a L{Data<escript.Data>}
1981                           object on L{ReducedFunction<escript.ReducedFunction>}
1982          @keyword D: value for coefficient C{D}
1983          @type D: any type that can be cast to a L{Data<escript.Data>} object on
1984                   L{Function<escript.Function>}
1985          @keyword D_reduced: value for coefficient C{D_reduced}
1986          @type D_reduced: any type that can be cast to a L{Data<escript.Data>}
1987                           object on L{ReducedFunction<escript.ReducedFunction>}
1988          @keyword X: value for coefficient C{X}
1989          @type X: any type that can be cast to a L{Data<escript.Data>} object on
1990                   L{Function<escript.Function>}
1991          @keyword X_reduced: value for coefficient C{X_reduced}
1992          @type X_reduced: any type that can be cast to a L{Data<escript.Data>}
1993                           object on L{ReducedFunction<escript.ReducedFunction>}
1994          @keyword Y: value for coefficient C{Y}
1995          @type Y: any type that can be cast to a L{Data<escript.Data>} object on
1996                   L{Function<escript.Function>}
1997          @keyword Y_reduced: value for coefficient C{Y_reduced}
1998          @type Y_reduced: any type that can be cast to a L{Data<escript.Data>}
1999                           object on L{ReducedFunction<escript.Function>}
2000          @keyword d: value for coefficient C{d}
2001          @type d: any type that can be cast to a L{Data<escript.Data>} object on
2002                   L{FunctionOnBoundary<escript.FunctionOnBoundary>}
2003          @keyword d_reduced: value for coefficient C{d_reduced}
2004          @type d_reduced: any type that can be cast to a L{Data<escript.Data>}
2005                           object on L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}
2006          @keyword y: value for coefficient C{y}
2007          @type y: any type that can be cast to a L{Data<escript.Data>} object on
2008                   L{FunctionOnBoundary<escript.FunctionOnBoundary>}
2009          @keyword d_contact: value for coefficient C{d_contact}
2010          @type d_contact: any type that can be cast to a L{Data<escript.Data>}
2011                           object on L{FunctionOnContactOne<escript.FunctionOnContactOne>}
2012                           or L{FunctionOnContactZero<escript.FunctionOnContactZero>}
2013          @keyword d_contact_reduced: value for coefficient C{d_contact_reduced}
2014          @type d_contact_reduced: any type that can be cast to a L{Data<escript.Data>}
2015                                   object on L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}
2016                                   or L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>}
2017          @keyword y_contact: value for coefficient C{y_contact}
2018          @type y_contact: any type that can be cast to a L{Data<escript.Data>}
2019                           object on L{FunctionOnContactOne<escript.FunctionOnContactOne>}
2020                           or L{FunctionOnContactZero<escript.FunctionOnContactZero>}
2021          @keyword y_contact_reduced: value for coefficient C{y_contact_reduced}
2022          @type y_contact_reduced: any type that can be cast to a L{Data<escript.Data>}
2023                                   object on L{ReducedFunctionOnContactOne<escript.FunctionOnContactOne>}
2024                                   or L{ReducedFunctionOnContactZero<escript.FunctionOnContactZero>}
2025          @keyword r: values prescribed to the solution at the locations of
2026                      constraints
2027          @type r: any type that can be cast to a L{Data<escript.Data>} object on
2028                   L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2029                   depending on whether reduced order is used for the solution
2030          @keyword q: mask for location of constraints
2031          @type q: any type that can be cast to a L{Data<escript.Data>} object on
2032                   L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2033                   depending on whether reduced order is used for the
2034                   representation of the equation
2035          @raise IllegalCoefficient: if an unknown coefficient keyword is used
2036          """
2037          super(LinearPDE,self).setValue(**coefficients)
2038          # check if the systrem is inhomogeneous:
2039          if len(coefficients)>0 and not self.isUsingLumping():
2040             q=self.getCoefficient("q")
2041             r=self.getCoefficient("r")
2042             if not q.isEmpty() and not r.isEmpty():
2043                 if util.Lsup(q*r)>0.:
2044                   self.trace("Inhomogeneous constraint detected.")
2045                   self.invalidateSystem()
2046    
2047    
2048       def getResidual(self,u=None):
2049         """
2050         Returns the residual of u or the current solution if u is not present.
2051    
2052         @param u: argument in the residual calculation. It must be representable
2053                   in L{self.getFunctionSpaceForSolution()}. If u is not present
2054                   or equals C{None} the current solution is used.
2055         @type u: L{Data<escript.Data>} or None
2056         @return: residual of u
2057         @rtype: L{Data<escript.Data>}
2058         """
2059         if u==None:
2060            return self.getOperator()*self.getSolution()-self.getRightHandSide()
2061         else:
2062            return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())-self.getRightHandSide()
2063    
2064     def getSolution(self,**options):     def getFlux(self,u=None):
2065         """       """
2066         returns the solution of the PDE       Returns the flux M{J} for a given M{u}.
2067    
2068         @param options:       M{J[i,j]=(A[i,j,k,l]+A_reduced[A[i,j,k,l]]*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])u[k]-X[i,j]-X_reduced[i,j]}
        """  
        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  
2069    
2070         or
2071    
2072         M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[l]+(B[j]+B_reduced[j])u-X[j]-X_reduced[j]}
2073    
2074         @param u: argument in the flux. If u is not present or equals L{None} the
2075                   current solution is used.
2076         @type u: L{Data<escript.Data>} or None
2077         @return: flux
2078         @rtype: L{Data<escript.Data>}
2079         """
2080         if u==None: u=self.getSolution()
2081         return util.tensormult(self.getCoefficient("A"),util.grad(u,Funtion(self.getDomain))) \
2082               +util.matrixmult(self.getCoefficient("B"),u) \
2083               -util.self.getCoefficient("X") \
2084               +util.tensormult(self.getCoefficient("A_reduced"),util.grad(u,ReducedFuntion(self.getDomain))) \
2085               +util.matrixmult(self.getCoefficient("B_reduced"),u) \
2086               -util.self.getCoefficient("X_reduced")
2087    
2088  def ELMAN_RAMAGE(P):  
2089       """   """  class Poisson(LinearPDE):
2090       return (P-1.).wherePositive()*0.5*(1.-1./(P+1.e-15))     """
2091  def SIMPLIFIED_BROOK_HUGHES(P):     Class to define a Poisson equation problem. This is generally a
2092       """   """     L{LinearPDE} of the form
2093       c=(P-3.).whereNegative()  
2094       return P/6.*c+1./2.*(1.-c)     M{-grad(grad(u)[j])[j] = f}
2095    
2096  def HALF(P):     with natural boundary conditions
2097      """ """  
2098      return escript.Scalar(0.5,P.getFunctionSpace())     M{n[j]*grad(u)[j] = 0 }
   
 class AdvectivePDE(LinearPDE):  
    """  
    Class to handle a linear PDE dominated by advective terms:  
     
    class to define a linear PDE of the form  
   
    \f[  
    -(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  
    \f]  
   
    with boundary conditons:  
   
    \f[  
    n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_{ik}u_k = - n_j*X_{ij} + y_i  
    \f]  
   
    and contact conditions  
   
    \f[  
    n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d^{contact}_{ik}[u_k] = - n_j*X_{ij} + y^{contact}_{i}  
    \f]  
2099    
2100     and constraints:     and constraints:
2101    
2102     \f[     M{u=0} where M{q>0}
    u_i=r_i \quad \mathrm{where} \quad q_i>0  
    \f]  
    """  
    def __init__(self,domain,numEquations=0,numSolutions=0,xi=ELMAN_RAMAGE):  
       LinearPDE.__init__(self,domain,numEquations,numSolutions)  
       self.__xi=xi  
       self.__Xi=escript.Data()  
   
    def __calculateXi(self,peclet_factor,Z,h):  
        Z_max=util.Lsup(Z)  
        if Z_max>0.:  
           return h*self.__xi(Z*peclet_factor)/(Z+Z_max*self.TOL)  
        else:  
           return 0.  
2103    
2104     def setValue(self,**args):     """
        if "A" in args.keys()   or "B" in args.keys() or "C" in args.keys(): self.__Xi=escript.Data()  
        self._LinearPDE__setValue(**args)  
             
    def getXi(self):  
       if self.__Xi.isEmpty():  
          B=self.getCoefficient("B")  
          C=self.getCoefficient("C")  
          A=self.getCoefficient("A")  
          h=self.getDomain().getSize()  
          self.__Xi=escript.Scalar(0.,self.getFunctionSpaceForCoefficient("A"))  
          if not C.isEmpty() or not B.isEmpty():  
             if not C.isEmpty() and not B.isEmpty():  
                 Z2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))  
                 if self.getNumEquations()>1:  
                    if self.getNumSolutions()>1:  
                       for i in range(self.getNumEquations()):  
                          for k in range(self.getNumSolutions()):  
                             for l in range(self.getDim()): Z2+=(C[i,k,l]-B[i,l,k])**2  
                    else:  
                       for i in range(self.getNumEquations()):  
                          for l in range(self.getDim()): Z2+=(C[i,l]-B[i,l])**2  
                 else:  
                    if self.getNumSolutions()>1:  
                       for k in range(self.getNumSolutions()):  
                          for l in range(self.getDim()): Z2+=(C[k,l]-B[l,k])**2  
                    else:  
                       for l in range(self.getDim()): Z2+=(C[l]-B[l])**2  
                 length_of_Z=util.sqrt(Z2)  
             elif C.isEmpty():  
               length_of_Z=util.length(B)  
             else:  
               length_of_Z=util.length(C)  
2105    
2106              Z_max=util.Lsup(length_of_Z)     def __init__(self,domain,debug=False):
2107              if Z_max>0.:       """
2108                 length_of_A=util.length(A)       Initializes a new Poisson equation.
                A_max=util.Lsup(length_of_A)  
                if A_max>0:  
                     inv_A=1./(length_of_A+A_max*self.TOL)  
                else:  
                     inv_A=1./self.TOL  
                peclet_number=length_of_Z*h/2*inv_A  
                xi=self.__xi(peclet_number)  
                self.__Xi=h*xi/(length_of_Z+Z_max*self.TOL)  
                print "@ preclet number = %e"%util.Lsup(peclet_number),util.Lsup(xi),util.Lsup(length_of_Z)  
       return self.__Xi  
         
   
    def getCoefficientOfPDE(self,name):  
      """  
      return the value of the coefficient name of the general PDE  
   
      @param name:  
      """  
      if not self.getNumEquations() == self.getNumSolutions():  
           raise ValueError,"AdvectivePDE expects the number of solution componets and the number of equations to be equal."  
   
      if name == "A" :  
          A=self.getCoefficient("A")  
          B=self.getCoefficient("B")  
          C=self.getCoefficient("C")  
          if B.isEmpty() and C.isEmpty():  
             Aout=A  
          else:  
             if A.isEmpty():  
                Aout=self.createNewCoefficient("A")  
             else:  
                Aout=A[:]  
             Xi=self.getXi()  
             if self.getNumEquations()>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 not C.isEmpty() and not B.isEmpty():  
                                for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*(C[p,i,j]-B[p,j,i])*(C[p,k,l]-B[p,l,k])  
                             elif C.isEmpty():  
                                for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*B[p,j,i]*B[p,l,k]  
                             else:  
                                for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*C[p,i,j]*C[p,k,l]  
             else:  
                 for j in range(self.getDim()):  
                    for l in range(self.getDim()):  
                       if not C.isEmpty() and not B.isEmpty():  
                           Aout[j,l]+=Xi*(C[j]-B[j])*(C[l]-B[l])  
                       elif C.isEmpty():  
                           Aout[j,l]+=Xi*B[j]*B[l]  
                       else:  
                           Aout[j,l]+=Xi*C[j]*C[l]  
          return Aout  
      elif name == "B" :  
          B=self.getCoefficient("B")  
          C=self.getCoefficient("C")  
          D=self.getCoefficient("D")  
          if C.isEmpty() or D.isEmpty():  
             Bout=B  
          else:  
             Xi=self.getXi()  
             if B.isEmpty():  
                 Bout=self.createNewCoefficient("B")  
             else:  
                 Bout=B[:]  
             if self.getNumEquations()>1:  
                for k in range(self.getNumSolutions()):  
                   for p in range(self.getNumEquations()):  
                      tmp=Xi*D[p,k]  
                      for i in range(self.getNumEquations()):  
                         for j in range(self.getDim()):  
                            Bout[i,j,k]+=tmp*C[p,i,j]  
             else:  
                tmp=Xi*D  
                for j in range(self.getDim()): Bout[j]+=tmp*C[j]  
          return Bout  
      elif name == "C" :  
          B=self.getCoefficient("B")  
          C=self.getCoefficient("C")  
          D=self.getCoefficient("D")  
          if B.isEmpty() or D.isEmpty():  
             Cout=C  
          else:  
             Xi=self.getXi()  
             if C.isEmpty():  
                 Cout=self.createNewCoefficient("C")  
             else:  
                 Cout=C[:]  
             if self.getNumEquations()>1:  
                for k in range(self.getNumSolutions()):  
                    for p in range(self.getNumEquations()):  
                       tmp=Xi*D[p,k]  
                       for i in range(self.getNumEquations()):  
                         for l in range(self.getDim()):  
                                  Cout[i,k,l]+=tmp*B[p,l,i]  
             else:  
                tmp=Xi*D  
                for j in range(self.getDim()): Cout[j]+=tmp*B[j]  
          return Cout  
      elif name == "D" :  
          return self.getCoefficient("D")  
      elif name == "X" :  
          X=self.getCoefficient("X")  
          Y=self.getCoefficient("Y")  
          B=self.getCoefficient("B")  
          C=self.getCoefficient("C")  
          if Y.isEmpty() or (B.isEmpty() and C.isEmpty()):  
             Xout=X  
          else:  
             if X.isEmpty():  
                 Xout=self.createNewCoefficient("X")  
             else:  
                 Xout=X[:]  
             Xi=self.getXi()  
             if self.getNumEquations()>1:  
                  for p in range(self.getNumEquations()):  
                     tmp=Xi*Y[p]  
                     for i in range(self.getNumEquations()):  
                        for j in range(self.getDim()):  
                           if not C.isEmpty() and not B.isEmpty():  
                              Xout[i,j]+=tmp*(C[p,i,j]-B[p,j,i])  
                           elif C.isEmpty():  
                              Xout[i,j]-=tmp*B[p,j,i]  
                           else:  
                              Xout[i,j]+=tmp*C[p,i,j]  
             else:  
                  tmp=Xi*Y  
                  for j in range(self.getDim()):  
                     if not C.isEmpty() and not B.isEmpty():  
                        Xout[j]+=tmp*(C[j]-B[j])  
                     elif C.isEmpty():  
                        Xout[j]-=tmp*B[j]  
                     else:  
                        Xout[j]+=tmp*C[j]  
          return Xout  
      elif name == "Y" :  
          return self.getCoefficient("Y")  
      elif name == "d" :  
          return self.getCoefficient("d")  
      elif name == "y" :  
          return self.getCoefficient("y")  
      elif name == "d_contact" :  
          return self.getCoefficient("d_contact")  
      elif name == "y_contact" :  
          return self.getCoefficient("y_contact")  
      elif name == "r" :  
          return self.getCoefficient("r")  
      elif name == "q" :  
          return self.getCoefficient("q")  
      else:  
          raise SystemError,"unknown PDE coefficient %s",name  
2109    
2110         @param domain: domain of the PDE
2111         @type domain: L{Domain<escript.Domain>}
2112         @param debug: if True debug information is printed
2113    
2114  class Poisson(LinearPDE):       """
2115         super(Poisson, self).__init__(domain,1,1,debug)
2116         self.introduceCoefficients(
2117                            f=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2118                            f_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE))
2119         self.setSymmetryOn()
2120    
2121       def setValue(self,**coefficients):
2122         """
2123         Sets new values to coefficients.
2124    
2125         @param coefficients: new values assigned to coefficients
2126         @keyword f: value for right hand side M{f}
2127         @type f: any type that can be cast to a L{Scalar<escript.Scalar>} object
2128                  on L{Function<escript.Function>}
2129         @keyword q: mask for location of constraints
2130         @type q: any type that can be cast to a rank zero L{Data<escript.Data>}
2131                  object on L{Solution<escript.Solution>} or
2132                  L{ReducedSolution<escript.ReducedSolution>} depending on whether
2133                  reduced order is used for the representation of the equation
2134         @raise IllegalCoefficient: if an unknown coefficient keyword is used
2135         """
2136         super(Poisson, self).setValue(**coefficients)
2137    
2138    
2139       def getCoefficient(self,name):
2140         """
2141         Returns the value of the coefficient C{name} of the general PDE.
2142    
2143         @param name: name of the coefficient requested
2144         @type name: C{string}
2145         @return: the value of the coefficient C{name}
2146         @rtype: L{Data<escript.Data>}
2147         @raise IllegalCoefficient: invalid coefficient name
2148         @note: This method is called by the assembling routine to map the Poisson
2149                equation onto the general PDE.
2150         """
2151         if name == "A" :
2152             return escript.Data(util.kronecker(self.getDim()),escript.Function(self.getDomain()))
2153         elif name == "Y" :
2154             return self.getCoefficient("f")
2155         elif name == "Y_reduced" :
2156             return self.getCoefficient("f_reduced")
2157         else:
2158             return super(Poisson, self).getCoefficient(name)
2159    
2160    class Helmholtz(LinearPDE):
2161     """     """
2162     Class to define a Poisson equation problem:     Class to define a Helmholtz equation problem. This is generally a
2163       L{LinearPDE} of the form
2164    
2165     class to define a linear PDE of the form     M{S{omega}*u - grad(k*grad(u)[j])[j] = f}
2166     \f[  
2167     -u_{,jj} = f     with natural boundary conditions
2168     \f]  
2169       M{k*n[j]*grad(u)[j] = g- S{alpha}u }
    with boundary conditons:  
   
    \f[  
    n_j*u_{,j} = 0  
    \f]  
2170    
2171     and constraints:     and constraints:
2172    
2173     \f[     M{u=r} where M{q>0}
2174     u=0 \quad \mathrm{where} \quad q>0  
2175     \f]     """
2176     """  
2177       def __init__(self,domain,debug=False):
2178     def __init__(self,domain,f=escript.Data(),q=escript.Data()):       """
2179         LinearPDE.__init__(self,domain,1,1)       Initializes a new Helmholtz equation.
2180         self.COEFFICIENTS={  
2181         "f"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),       @param domain: domain of the PDE
2182         "q"         : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.BOTH)}       @type domain: L{Domain<escript.Domain>}
2183         self.setSymmetryOn()       @param debug: if True debug information is printed
2184         self.setValue(f,q)  
2185         """
2186     def setValue(self,f=escript.Data(),q=escript.Data()):       super(Helmholtz, self).__init__(domain,1,1,debug)
2187         """set value of PDE parameters f and q"""       self.introduceCoefficients(
2188         self._LinearPDE__setValue(f=f,q=q)                          omega=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.OPERATOR),
2189                            k=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.OPERATOR),
2190     def getCoefficientOfPDE(self,name):                          f=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2191       """                          f_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2192       return the value of the coefficient name of the general PDE                          alpha=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.OPERATOR),
2193                            g=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2194       @param name:                          g_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE))
2195       """       self.setSymmetryOn()
2196       if name == "A" :  
2197           return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))     def setValue(self,**coefficients):
2198       elif name == "B" :       """
2199           return escript.Data()       Sets new values to coefficients.
2200       elif name == "C" :  
2201           return escript.Data()       @param coefficients: new values assigned to coefficients
2202       elif name == "D" :       @keyword omega: value for coefficient M{S{omega}}
2203           return escript.Data()       @type omega: any type that can be cast to a L{Scalar<escript.Scalar>}
2204       elif name == "X" :                    object on L{Function<escript.Function>}
2205           return escript.Data()       @keyword k: value for coefficient M{k}
2206       elif name == "Y" :       @type k: any type that can be cast to a L{Scalar<escript.Scalar>} object
2207                  on L{Function<escript.Function>}
2208         @keyword f: value for right hand side M{f}
2209         @type f: any type that can be cast to a L{Scalar<escript.Scalar>} object
2210                  on L{Function<escript.Function>}
2211         @keyword alpha: value for right hand side M{S{alpha}}
2212         @type alpha: any type that can be cast to a L{Scalar<escript.Scalar>}
2213                      object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}
2214         @keyword g: value for right hand side M{g}
2215         @type g: any type that can be cast to a L{Scalar<escript.Scalar>} object
2216                  on L{FunctionOnBoundary<escript.FunctionOnBoundary>}
2217         @keyword r: prescribed values M{r} for the solution in constraints
2218         @type r: any type that can be cast to a L{Scalar<escript.Scalar>} object
2219                  on L{Solution<escript.Solution>} or
2220                  L{ReducedSolution<escript.ReducedSolution>} depending on whether
2221                  reduced order is used for the representation of the equation
2222         @keyword q: mask for the location of constraints
2223         @type q: any type that can be cast to a L{Scalar<escript.Scalar>} object
2224                  on L{Solution<escript.Solution>} or
2225                  L{ReducedSolution<escript.ReducedSolution>} depending on whether
2226                  reduced order is used for the representation of the equation
2227         @raise IllegalCoefficient: if an unknown coefficient keyword is used
2228         """
2229         super(Helmholtz, self).setValue(**coefficients)
2230    
2231       def getCoefficient(self,name):
2232         """
2233         Returns the value of the coefficient C{name} of the general PDE.
2234    
2235         @param name: name of the coefficient requested
2236         @type name: C{string}
2237         @return: the value of the coefficient C{name}
2238         @rtype: L{Data<escript.Data>}
2239         @raise IllegalCoefficient: invalid name
2240         """
2241         if name == "A" :
2242             if self.getCoefficient("k").isEmpty():
2243                  return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))
2244             else:
2245                  return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))*self.getCoefficient("k")
2246         elif name == "D" :
2247             return self.getCoefficient("omega")
2248         elif name == "Y" :
2249           return self.getCoefficient("f")           return self.getCoefficient("f")
2250       elif name == "d" :       elif name == "d" :
2251           return escript.Data()           return self.getCoefficient("alpha")
2252       elif name == "y" :       elif name == "y" :
2253           return escript.Data()           return self.getCoefficient("g")
2254       elif name == "d_contact" :       elif name == "Y_reduced" :
2255           return escript.Data()           return self.getCoefficient("f_reduced")
2256       elif name == "y_contact" :       elif name == "y_reduced" :
2257           return escript.Data()          return self.getCoefficient("g_reduced")
      elif name == "r" :  
          return escript.Data()  
      elif name == "q" :  
          return self.getCoefficient("q")  
2258       else:       else:
2259           raise SystemError,"unknown PDE coefficient %s",name          return super(Helmholtz, self).getCoefficient(name)
2260    
2261  class LameEquation(LinearPDE):  class LameEquation(LinearPDE):
2262     """     """
2263     Class to define a Lame equation problem:     Class to define a Lame equation problem. This problem is defined as:
2264    
2265       M{-grad(S{mu}*(grad(u[i])[j]+grad(u[j])[i]))[j] - grad(S{lambda}*grad(u[k])[k])[j] = F_i -grad(S{sigma}[ij])[j] }
2266    
2267     class to define a linear PDE of the form     with natural boundary conditions:
2268     \f[  
2269     -(\lambda (u_{i,j}+u_{j,i}))_{,j} - \mu u_{j,ji}} = F_i -\sigma_{ij,j}     M{n[j]*(S{mu}*(grad(u[i])[j]+grad(u[j])[i]) + S{lambda}*grad(u[k])[k]) = f_i +n[j]*S{sigma}[ij] }
    \f]  
   
    with boundary conditons:  
   
    \f[  
    n_j(\lambda(u_{i,j}+u_{j,i})-sigma_{ij}) + n_i\mu u_{j,j} = f_i  
    \f]  
2270    
2271     and constraints:     and constraints:
2272    
2273     \f[     M{u[i]=r[i]} where M{q[i]>0}
2274     u_i=r_i \quad \mathrm{where} \quad q_i>0  
2275     \f]     """
2276     """  
2277       def __init__(self,domain,debug=False):
2278     def __init__(self,domain,f=escript.Data(),q=escript.Data()):        """
2279         LinearPDE.__init__(self,domain,domain.getDim(),domain.getDim())        Initializes a new Lame equation.
2280         self.COEFFICIENTS={  
2281         "lame_lambda"  : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),        @param domain: domain of the PDE
2282         "lame_mu"      : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),        @type domain: L{Domain<escript.Domain>}
2283         "F"            : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM),PDECoefficient.RIGHTHANDSIDE),        @param debug: if True debug information is printed
2284         "sigma"        : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.EQUATION),PDECoefficient.RIGHTHANDSIDE),  
2285         "f"            : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),        """
2286         "r"            : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.BOTH),        super(LameEquation, self).__init__(domain,\
2287         "q"            : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.BOTH)}                                           domain.getDim(),domain.getDim(),debug)
2288         self.setSymmetryOn()        self.introduceCoefficients(lame_lambda=PDECoef(PDECoef.INTERIOR,(),PDECoef.OPERATOR),
2289                                     lame_mu=PDECoef(PDECoef.INTERIOR,(),PDECoef.OPERATOR),
2290     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()):                                   F=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2291         """set value of PDE parameters"""                                   sigma=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
2292         self._LinearPDE__setValue(lame_lambda=lame_lambda, \                                   f=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE))
2293                                   lame_mu=lame_mu, \        self.setSymmetryOn()
2294                                   F=F, \  
2295                                   sigma=sigma, \     def setValues(self,**coefficients):
2296                                   f=f, \       """
2297                                   r=r, \       Sets new values to coefficients.
2298                                   q=q)  
2299     def getCoefficientOfPDE(self,name):       @param coefficients: new values assigned to coefficients
2300       """       @keyword lame_mu: value for coefficient M{S{mu}}
2301       return the value of the coefficient name of the general PDE       @type lame_mu: any type that can be cast to a L{Scalar<escript.Scalar>}
2302                        object on L{Function<escript.Function>}
2303       @param name:       @keyword lame_lambda: value for coefficient M{S{lambda}}
2304       """       @type lame_lambda: any type that can be cast to a L{Scalar<escript.Scalar>}
2305       if name == "A" :                          object on L{Function<escript.Function>}
2306           A =self.createNewCoefficient("A")       @keyword F: value for internal force M{F}
2307           for i in range(self.getDim()):       @type F: any type that can be cast to a L{Vector<escript.Vector>} object
2308             for j in range(self.getDim()):                on L{Function<escript.Function>}
2309               out[i,i,j,j] += self.getCoefficient("lame_mu")       @keyword sigma: value for initial stress M{S{sigma}}
2310               out[i,j,j,i] += self.getCoefficient("lame_lambda")       @type sigma: any type that can be cast to a L{Tensor<escript.Tensor>}
2311               out[i,j,i,j] += self.getCoefficient("lame_lambda")                    object on L{Function<escript.Function>}
2312         @keyword f: value for external force M{f}
2313         @type f: any type that can be cast to a L{Vector<escript.Vector>} object
2314                  on L{FunctionOnBoundary<escript.FunctionOnBoundary>}
2315         @keyword r: prescribed values M{r} for the solution in constraints
2316         @type r: any type that can be cast to a L{Vector<escript.Vector>} object
2317                  on L{Solution<escript.Solution>} or
2318                  L{ReducedSolution<escript.ReducedSolution>} depending on whether
2319                  reduced order is used for the representation of the equation
2320         @keyword q: mask for the location of constraints
2321         @type q: any type that can be cast to a L{Vector<escript.Vector>} object
2322                  on L{Solution<escript.Solution>} or
2323                  L{ReducedSolution<escript.ReducedSolution>} depending on whether
2324                  reduced order is used for the representation of the equation
2325         @raise IllegalCoefficient: if an unknown coefficient keyword is used
2326         """
2327         super(LameEquation, self).setValues(**coefficients)
2328    
2329       def getCoefficient(self,name):
2330         """
2331         Returns the value of the coefficient C{name} of the general PDE.
2332    
2333         @param name: name of the coefficient requested
2334         @type name: C{string}
2335         @return: the value of the coefficient C{name}
2336         @rtype: L{Data<escript.Data>}
2337         @raise IllegalCoefficient: invalid coefficient name
2338         """
2339         out =self.createCoefficient("A")
2340         if name == "A" :
2341             if self.getCoefficient("lame_lambda").isEmpty():
2342                if self.getCoefficient("lame_mu").isEmpty():
2343                    pass
2344                else:
2345                    for i in range(self.getDim()):
2346                      for j in range(self.getDim()):
2347                        out[i,j,j,i] += self.getCoefficient("lame_mu")
2348                        out[i,j,i,j] += self.getCoefficient("lame_mu")
2349             else:
2350                if self.getCoefficient("lame_mu").isEmpty():
2351                    for i in range(self.getDim()):
2352                      for j in range(self.getDim()):
2353                        out[i,i,j,j] += self.getCoefficient("lame_lambda")
2354                else:
2355                    for i in range(self.getDim()):
2356                      for j in range(self.getDim()):
2357                        out[i,i,j,j] += self.getCoefficient("lame_lambda")
2358                        out[i,j,j,i] += self.getCoefficient("lame_mu")
2359                        out[i,j,i,j] += self.getCoefficient("lame_mu")
2360           return out           return out
2361       elif name == "B" :       elif name == "X" :
          return escript.Data()  
      elif name == "C" :  
          return escript.Data()  
      elif name == "D" :  
          return escript.Data()  
      elif name == "X" :  
2362           return self.getCoefficient("sigma")           return self.getCoefficient("sigma")
2363       elif name == "Y" :       elif name == "Y" :
2364           return self.getCoefficient("F")           return self.getCoefficient("F")
2365       elif name == "d" :       elif name == "y" :
          return escript.Data()  
      elif name == "y" :  
2366           return self.getCoefficient("f")           return self.getCoefficient("f")
      elif name == "d_contact" :  
          return escript.Data()  
      elif name == "y_contact" :  
          return escript.Data()  
      elif name == "r" :  
          return self.getCoefficient("r")  
      elif name == "q" :  
          return self.getCoefficient("q")  
2367       else:       else:
2368           raise SystemError,"unknown PDE coefficient %s",name          return super(LameEquation, self).getCoefficient(name)
2369    
2370  # $Log$  def LinearSinglePDE(domain,debug=False):
2371  # Revision 1.9  2005/07/25 05:28:13  jgs     """
2372  # Merge of development branch back to main trunk on 2005-07-25     Defines a single linear PDE.
2373  #  
2374  # Revision 1.8  2005/06/09 05:37:59  jgs     @param domain: domain of the PDE
2375  # Merge of development branch back to main trunk on 2005-06-09     @type domain: L{Domain<escript.Domain>}
2376  #     @param debug: if True debug information is printed
2377  # Revision 1.7  2005/05/06 04:26:10  jgs     @rtype: L{LinearPDE}
2378  # Merge of development branch back to main trunk on 2005-05-06     """
2379  #     return LinearPDE(domain,numEquations=1,numSolutions=1,debug=debug)
2380  # Revision 1.1.2.24  2005/07/22 06:37:11  gross  
2381  # some extensions to modellib and linearPDEs  def LinearPDESystem(domain,debug=False):
2382  #     """
2383  # Revision 1.1.2.23  2005/05/13 00:55:20  cochrane     Defines a system of linear PDEs.
2384  # Fixed up some docstrings.  Moved module-level functions to top of file so  
2385  # that epydoc and doxygen can pick them up properly.     @param domain: domain of the PDEs
2386  #     @type domain: L{Domain<escript.Domain>}
2387  # Revision 1.1.2.22  2005/05/12 11:41:30  gross     @param debug: if True debug information is printed
2388  # some basic Models have been added     @rtype: L{LinearPDE}
2389  #     """
2390  # Revision 1.1.2.21  2005/05/12 07:16:12  cochrane     return LinearPDE(domain,numEquations=domain.getDim(),numSolutions=domain.getDim(),debug=debug)
2391  # Moved ELMAN_RAMAGE, SIMPLIFIED_BROOK_HUGHES, and HALF functions to bottom of  
2392  # file so that the AdvectivePDE class is picked up by doxygen.  Some  
2393  # reformatting of docstrings.  Addition of code to make equations come out  class TransportPDE(LinearProblem):
2394  # as proper LaTeX.     """
2395  #     This class is used to define a transport problem given by a general linear,
2396  # Revision 1.1.2.20  2005/04/15 07:09:08  gross     time dependent, second order PDE for an unknown, non-negative,
2397  # some problems with functionspace and linearPDEs fixed.     time-dependent function M{u} on a given domain defined through a
2398  #     L{Domain<escript.Domain>} object.
2399  # Revision 1.1.2.19  2005/03/04 05:27:07  gross  
2400  # bug in SystemPattern fixed.     For a single equation with a solution with a single component the transport
2401  #     problem is defined in the following form:
2402  # Revision 1.1.2.18  2005/02/08 06:16:45  gross  
2403  # Bugs in AdvectivePDE fixed, AdvectiveTest is stable but more testing is needed     M{(M+M_reduced)*u_t=-(grad(A[j,l]+A_reduced[j,l])*grad(u)[l]+(B[j]+B_reduced[j])u)[j]+(C[l]+C_reduced[l])*grad(u)[l]+(D+D_reduced)-grad(X+X_reduced)[j,j]+(Y+Y_reduced)}
2404  #  
2405  # Revision 1.1.2.17  2005/02/08 05:56:19  gross     where M{u_t} denotes the time derivative of M{u} and M{grad(F)} denotes the
2406  # Reference Number handling added     spatial derivative of M{F}. Einstein's summation convention,  ie. summation
2407  #     over indexes appearing twice in a term of a sum performed, is used.
2408  # Revision 1.1.2.16  2005/02/07 04:41:28  gross     The coefficients M{M}, M{A}, M{B}, M{C}, M{D}, M{X} and M{Y} have to be
2409  # some function exposed to python to make mesh merging running     specified through L{Data<escript.Data>} objects in L{Function<escript.Function>}
2410  #     and the coefficients M{M_reduced}, M{A_reduced}, M{B_reduced}, M{C_reduced},
2411  # Revision 1.1.2.15  2005/02/03 00:14:44  gross     M{D_reduced}, M{X_reduced} and M{Y_reduced} have to be specified through
2412  # timeseries add and ESySParameter.py renames esysXML.py for consistence     L{Data<escript.Data>} objects in L{ReducedFunction<escript.ReducedFunction>}.
2413  #     It is also allowed to use objects that can be converted into such
2414  # Revision 1.1.2.14  2005/02/01 06:44:10  gross     L{Data<escript.Data>} objects. M{M} and M{M_reduced} are scalar, M{A} and
2415  # new implementation of AdvectivePDE which now also updates right hand side. systems of PDEs are still not working     M{A_reduced} are rank two, M{B}, M{C}, M{X}, M{B_reduced}, M{C_reduced} and
2416  #     M{X_reduced} are rank one and M{D}, M{D_reduced}, M{Y} and M{Y_reduced} are
2417  # Revision 1.1.2.13  2005/01/25 00:47:07  gross     scalar.
2418  # updates in the documentation  
2419  #     The following natural boundary conditions are considered:
2420  # Revision 1.1.2.12  2005/01/12 01:28:04  matt  
2421  # Added createCoefficient method for linearPDEs.     M{n[j]*((A[i,j]+A_reduced[i,j])*grad(u)[l]+(B+B_reduced)[j]*u+X[j]+X_reduced[j])+(d+d_reduced)*u+y+y_reduced=(m+m_reduced)*u_t}
2422  #  
2423  # Revision 1.1.2.11  2005/01/11 01:55:34  gross     where M{n} is the outer normal field. Notice that the coefficients M{A},
2424  # a problem in linearPDE class fixed     M{A_reduced}, M{B}, M{B_reduced}, M{X} and M{X_reduced} are defined in the
2425  #     transport problem. The coefficients M{m}, M{d} and M{y} are each a scalar in
2426  # Revision 1.1.2.10  2005/01/07 01:13:29  gross     L{FunctionOnBoundary<escript.FunctionOnBoundary>} and the coefficients
2427  # some bugs in linearPDE fixed     M{m_reduced}, M{d_reduced} and M{y_reduced} are each a scalar in
2428  #     L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
2429  # Revision 1.1.2.9  2005/01/06 06:24:58  gross  
2430  # some bugs in slicing fixed     Constraints for the solution prescribing the value of the solution at
2431  #     certain locations in the domain have the form
2432  # Revision 1.1.2.8  2005/01/05 04:21:40  gross  
2433  # FunctionSpace checking/matchig in slicing added     M{u_t=r} where M{q>0}
2434  #  
2435  # Revision 1.1.2.7  2004/12/29 10:03:41  gross     M{r} and M{q} are each scalar where M{q} is the characteristic function
2436  # bug in setValue fixed     defining where the constraint is applied. The constraints override any other
2437  #     condition set by the transport problem or the boundary condition.
2438  # Revision 1.1.2.6  2004/12/29 05:29:59  gross  
2439  # AdvectivePDE successfully tested for Peclet number 1000000. there is still a problem with setValue and Data()     The transport problem is symmetrical if
2440  #  
2441  # Revision 1.1.2.5  2004/12/29 00:18:41  gross     M{A[i,j]=A[j,i]} and M{B[j]=C[j]} and M{A_reduced[i,j]=A_reduced[j,i]} and
2442  # AdvectivePDE added     M{B_reduced[j]=C_reduced[j]}
2443  #  
2444  # Revision 1.1.2.4  2004/12/24 06:05:41  gross     For a system and a solution with several components the transport problem
2445  # some changes in linearPDEs to add AdevectivePDE     has the form
2446  #  
2447  # Revision 1.1.2.3  2004/12/16 00:12:34  gross     M{(M[i,k]+M_reduced[i,k])*u[k]_t=-grad((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k])[j]+(C[i,k,l]+C_reduced[i,k,l])*grad(u[k])[l]+(D[i,k]+D_reduced[i,k]*u[k]-grad(X[i,j]+X_reduced[i,j])[j]+Y[i]+Y_reduced[i] }
2448  # __init__ of LinearPDE does not accept any coefficient anymore  
2449  #     M{A} and M{A_reduced} are of rank four, M{B}, M{B_reduced}, M{C} and
2450  # Revision 1.1.2.2  2004/12/14 03:55:01  jgs     M{C_reduced} are each of rank three, M{M}, M{M_reduced}, M{D}, M{D_reduced},
2451  # *** empty log message ***     M{X_reduced} and M{X} are each of rank two and M{Y} and M{Y_reduced} are of
2452  #     rank one. The natural boundary conditions take the form:
2453  # Revision 1.1.2.1  2004/12/12 22:53:47  gross  
2454  # linearPDE has been renamed LinearPDE     M{n[j]*((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k]+X[i,j]+X_reduced[i,j])+(d[i,k]+d_reduced[i,k])*u[k]+y[i]+y_reduced[i]= (m[i,k]+m_reduced[i,k])*u[k]_t}
2455  #  
2456  # Revision 1.1.1.1.2.7  2004/12/07 10:13:08  gross     The coefficient M{d} and M{m} are of rank two and M{y} is of rank one with
2457  # GMRES added     all in L{FunctionOnBoundary<escript.FunctionOnBoundary>}. The coefficients
2458  #     M{d_reduced} and M{m_reduced} are of rank two and M{y_reduced} is of rank
2459  # Revision 1.1.1.1.2.6  2004/12/07 03:19:50  gross     one all in L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
2460  # options for GMRES and PRES20 added  
2461  #     Constraints take the form
2462  # Revision 1.1.1.1.2.5  2004/12/01 06:25:15  gross  
2463  # some small changes     M{u[i]_t=r[i]} where M{q[i]>0}
2464  #  
2465  # Revision 1.1.1.1.2.4  2004/11/24 01:50:21  gross     M{r} and M{q} are each rank one. Notice that at some locations not
2466  # Finley solves 4M unknowns now     necessarily all components must have a constraint.
2467  #  
2468  # Revision 1.1.1.1.2.3  2004/11/15 06:05:26  gross     The transport problem is symmetrical if
2469  # poisson solver added  
2470  #        - M{M[i,k]=M[i,k]}
2471  # Revision 1.1.1.1.2.2  2004/11/12 06:58:15  gross        - M{M_reduced[i,k]=M_reduced[i,k]}
2472  # a lot of changes to get the linearPDE class running: most important change is that there is no matrix format exposed to the user anymore. the format is chosen by the Domain according to the solver and symmetry        - M{A[i,j,k,l]=A[k,l,i,j]}
2473  #        - M{A_reduced[i,j,k,l]=A_reduced[k,l,i,j]}
2474  # Revision 1.1.1.1.2.1  2004/10/28 22:59:22  gross        - M{B[i,j,k]=C[k,i,j]}
2475  # finley's RecTest.py is running now: problem in SystemMatrixAdapater fixed        - M{B_reduced[i,j,k]=C_reduced[k,i,j]}
2476  #        - M{D[i,k]=D[i,k]}
2477  # Revision 1.1.1.1  2004/10/26 06:53:56  jgs        - M{D_reduced[i,k]=D_reduced[i,k]}
2478  # initial import of project esys2        - M{m[i,k]=m[k,i]}
2479  #        - M{m_reduced[i,k]=m_reduced[k,i]}
2480  # Revision 1.3.2.3  2004/10/26 06:43:48  jgs        - M{d[i,k]=d[k,i]}
2481  # committing Lutz's and Paul's changes to brach jgs        - M{d_reduced[i,k]=d_reduced[k,i]}
2482  #  
2483  # Revision 1.3.4.1  2004/10/20 05:32:51  cochrane     L{TransportPDE} also supports solution discontinuities over a contact region
2484  # Added incomplete Doxygen comments to files, or merely put the docstrings that already exist into Doxygen form.     in the domain. To specify the conditions across the discontinuity we are
2485  #     using the generalised flux M{J} which, in the case of a system of PDEs and
2486  # Revision 1.3  2004/09/23 00:53:23  jgs     several components of the solution, is defined as
2487  # minor fixes  
2488  #     M{J[i,j]=(A[i,j,k,l]+A_reduced[[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k]+X[i,j]+X_reduced[i,j]}
2489  # Revision 1.1  2004/08/28 12:58:06  gross  
2490  # SimpleSolve is not running yet: problem with == of functionsspace     For the case of single solution component and single PDE M{J} is defined as
2491  #  
2492  #     M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[j]+(B[i]+B_reduced[i])*u+X[i]+X_reduced[i]}
2493    
2494       In the context of discontinuities M{n} denotes the normal on the
2495       discontinuity pointing from side 0 towards side 1 calculated from
2496       L{getNormal<escript.FunctionSpace.getNormal>} of L{FunctionOnContactZero<escript.FunctionOnContactZero>}.
2497       For a system of transport problems the contact condition takes the form
2498    
2499       M{n[j]*J0[i,j]=n[j]*J1[i,j]=(y_contact[i]+y_contact_reduced[i])- (d_contact[i,k]+d_contact_reduced[i,k])*jump(u)[k]}
2500    
2501       where M{J0} and M{J1} are the fluxes on side 0 and side 1 of the
2502       discontinuity, respectively. M{jump(u)}, which is the difference of the
2503       solution at side 1 and at side 0, denotes the jump of M{u} across
2504       discontinuity along the normal calculated by L{jump<util.jump>}.
2505       The coefficient M{d_contact} is of rank two and M{y_contact} is of rank one
2506       both in L{FunctionOnContactZero<escript.FunctionOnContactZero>} or L{FunctionOnContactOne<escript.FunctionOnContactOne>}.
2507       The coefficient M{d_contact_reduced} is of rank two and M{y_contact_reduced}
2508       is of rank one both in L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.
2509       In case of a single PDE and a single component solution the contact
2510       condition takes the form
2511    
2512       M{n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)}
2513    
2514       In this case the coefficient M{d_contact} and M{y_contact} are each scalar
2515       both in L{FunctionOnContactZero<escript.FunctionOnContactZero>} or
2516       L{FunctionOnContactOne<escript.FunctionOnContactOne>} and the coefficient
2517       M{d_contact_reduced} and M{y_contact_reduced} are each scalar both in
2518       L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or
2519       L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.
2520    
2521       Typical usage::
2522    
2523           p = TransportPDE(dom)
2524           p.setValue(M=1., C=[-1.,0.])
2525           p.setInitialSolution(u=exp(-length(dom.getX()-[0.1,0.1])**2)
2526           t = 0
2527           dt = 0.1
2528           while (t < 1.):
2529               u = p.solve(dt)
2530    
2531       """
2532       def __init__(self,domain,numEquations=None,numSolutions=None, useBackwardEuler=False, debug=False):
2533         """
2534         Initializes a transport problem.
2535    
2536         @param domain: domain of the PDE
2537         @type domain: L{Domain<escript.Domain>}
2538         @param numEquations: number of equations. If C{None} the number of
2539                              equations is extracted from the coefficients.
2540         @param numSolutions: number of solution components. If C{None} the number
2541                              of solution components is extracted from the
2542                              coefficients.
2543         @param debug: if True debug information is printed
2544         @param useBackwardEuler: if set the backward Euler scheme is used. Otherwise the Crank-Nicholson scheme is applied. Note that backward Euler scheme will return a safe time step size which is practically infinity as the scheme is unconditional unstable. The Crank-Nicholson scheme provides a higher accuracy but requires to limit the time step size to be stable.
2545         @type useBackwardEuler: C{bool}
2546         """
2547         if useBackwardEuler:
2548             self.__useBackwardEuler=True
2549         else:
2550             self.__useBackwardEuler=False
2551         super(TransportPDE, self).__init__(domain,numEquations,numSolutions,debug)
2552    
2553         self.setConstraintWeightingFactor()
2554         #
2555         #   the coefficients of the transport problem
2556         #
2557         self.introduceCoefficients(
2558           M=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2559           A=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2560           B=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2561           C=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2562           D=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2563           X=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
2564           Y=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2565           m=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2566           d=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2567           y=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2568           d_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2569           y_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2570           M_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2571           A_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2572           B_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2573           C_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2574           D_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2575           X_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
2576           Y_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2577           m_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2578           d_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2579           y_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2580           d_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2581           y_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2582           r=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.RIGHTHANDSIDE),
2583           q=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.BOTH) )
2584    
2585       def __str__(self):
2586         """
2587         Returns the string representation of the transport problem.
2588    
2589         @return: a simple representation of the transport problem
2590         @rtype: C{str}
2591         """
2592         return "<TransportPDE %d>"%id(self)
2593    
2594       def useBackwardEuler(self):
2595          """
2596          Returns true if backward Euler is used. Otherwise false is returned.
2597          @rtype: bool
2598          """
2599          return self.__useBackwardEuler
2600    
2601    
2602       def checkSymmetry(self,verbose=True):
2603          """
2604          Tests the transport problem for symmetry.
2605    
2606          @param verbose: if set to True or not present a report on coefficients
2607                          which break the symmetry is printed.
2608          @type verbose: C{bool}
2609          @return:  True if the PDE is symmetric
2610          @rtype: C{bool}
2611          @note: This is a very expensive operation. It should be used for
2612                 degugging only! The symmetry flag is not altered.
2613          """
2614          out=True
2615          out=out and self.checkSymmetricTensor("M", verbose)
2616          out=out and self.checkSymmetricTensor("M_reduced", verbose)
2617          out=out and self.checkSymmetricTensor("A", verbose)
2618          out=out and self.checkSymmetricTensor("A_reduced", verbose)
2619          out=out and self.checkReciprocalSymmetry("B","C", verbose)
2620          out=out and self.checkReciprocalSymmetry("B_reduced","C_reduced", verbose)
2621          out=out and self.checkSymmetricTensor("D", verbose)
2622          out=out and self.checkSymmetricTensor("D_reduced", verbose)
2623          out=out and self.checkSymmetricTensor("m", verbose)
2624          out=out and self.checkSymmetricTensor("m_reduced", verbose)
2625          out=out and self.checkSymmetricTensor("d", verbose)
2626          out=out and self.checkSymmetricTensor("d_reduced", verbose)
2627          out=out and self.checkSymmetricTensor("d_contact", verbose)
2628          out=out and self.checkSymmetricTensor("d_contact_reduced", verbose)
2629          return out
2630    
2631       def setValue(self,**coefficients):
2632          """
2633          Sets new values to coefficients.
2634    
2635          @param coefficients: new values assigned to coefficients
2636          @keyword M: value for coefficient C{M}
2637          @type M: any type that can be cast to a L{Data<escript.Data>} object on
2638                   L{Function<escript.Function>}
2639          @keyword M_reduced: value for coefficient C{M_reduced}
2640          @type M_reduced: any type that can be cast to a L{Data<escript.Data>}
2641                           object on L{Function<escript.ReducedFunction>}
2642          @keyword A: value for coefficient C{A}
2643          @type A: any type that can be cast to a L{Data<escript.Data>} object on
2644                   L{Function<escript.Function>}
2645          @keyword A_reduced: value for coefficient C{A_reduced}
2646          @type A_reduced: any type that can be cast to a L{Data<escript.Data>}
2647                           object on L{ReducedFunction<escript.ReducedFunction>}
2648          @keyword B: value for coefficient C{B}
2649          @type B: any type that can be cast to a L{Data<escript.Data>} object on
2650                   L{Function<escript.Function>}
2651          @keyword B_reduced: value for coefficient C{B_reduced}
2652          @type B_reduced: any type that can be cast to a L{Data<escript.Data>}
2653                           object on L{ReducedFunction<escript.ReducedFunction>}
2654          @keyword C: value for coefficient C{C}
2655          @type C: any type that can be cast to a L{Data<escript.Data>} object on
2656                   L{Function<escript.Function>}
2657          @keyword C_reduced: value for coefficient C{C_reduced}
2658          @type C_reduced: any type that can be cast to a L{Data<escript.Data>}
2659                           object on L{ReducedFunction<escript.ReducedFunction>}
2660          @keyword D: value for coefficient C{D}
2661          @type D: any type that can be cast to a L{Data<escript.Data>} object on
2662                   L{Function<escript.Function>}
2663          @keyword D_reduced: value for coefficient C{D_reduced}
2664          @type D_reduced: any type that can be cast to a L{Data<escript.Data>}
2665                           object on L{ReducedFunction<escript.ReducedFunction>}
2666          @keyword X: value for coefficient C{X}
2667          @type X: any type that can be cast to a L{Data<escript.Data>} object on
2668                   L{Function<escript.Function>}
2669          @keyword X_reduced: value for coefficient C{X_reduced}
2670          @type X_reduced: any type that can be cast to a L{Data<escript.Data>}
2671                           object on L{ReducedFunction<escript.ReducedFunction>}
2672          @keyword Y: value for coefficient C{Y}
2673          @type Y: any type that can be cast to a L{Data<escript.Data>} object on
2674                   L{Function<escript.Function>}
2675          @keyword Y_reduced: value for coefficient C{Y_reduced}
2676          @type Y_reduced: any type that can be cast to a L{Data<escript.Data>}
2677                           object on L{ReducedFunction<escript.Function>}
2678          @keyword m: value for coefficient C{m}
2679          @type m: any type that can be cast to a L{Data<escript.Data>} object on
2680                   L{FunctionOnBoundary<escript.FunctionOnBoundary>}
2681          @keyword m_reduced: value for coefficient C{m_reduced}
2682          @type m_reduced: any type that can be cast to a L{Data<escript.Data>}
2683                           object on L{FunctionOnBoundary<escript.ReducedFunctionOnBoundary>}
2684          @keyword d: value for coefficient C{d}
2685          @type d: any type that can be cast to a L{Data<escript.Data>} object on
2686                   L{FunctionOnBoundary<escript.FunctionOnBoundary>}
2687          @keyword d_reduced: value for coefficient C{d_reduced}
2688          @type d_reduced: any type that can be cast to a L{Data<escript.Data>}
2689                           object on L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}
2690          @keyword y: value for coefficient C{y}
2691          @type y: any type that can be cast to a L{Data<escript.Data>} object on
2692                   L{FunctionOnBoundary<escript.FunctionOnBoundary>}
2693          @keyword d_contact: value for coefficient C{d_contact}
2694          @type d_contact: any type that can be cast to a L{Data<escript.Data>}
2695                           object on L{FunctionOnContactOne<escript.FunctionOnContactOne>} or L{FunctionOnContactZero<escript.FunctionOnContactZero>}
2696          @keyword d_contact_reduced: value for coefficient C{d_contact_reduced}
2697          @type d_contact_reduced: any type that can be cast to a L{Data<escript.Data>} object on L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>} or L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>}
2698          @keyword y_contact: value for coefficient C{y_contact}
2699          @type y_contact: any type that can be cast to a L{Data<escript.Data>}
2700                           object on L{FunctionOnContactOne<escript.FunctionOnContactOne>} or L{FunctionOnContactZero<escript.FunctionOnContactZero>}
2701          @keyword y_contact_reduced: value for coefficient C{y_contact_reduced}
2702          @type y_contact_reduced: any type that can be cast to a L{Data<escript.Data>} object on L{ReducedFunctionOnContactOne<escript.FunctionOnContactOne>} or L{ReducedFunctionOnContactZero<escript.FunctionOnContactZero>}
2703          @keyword r: values prescribed to the solution at the locations of constraints
2704          @type r: any type that can be cast to a L{Data<escript.Data>} object on
2705                   L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2706                   depending on whether reduced order is used for the solution
2707          @keyword q: mask for the location of constraints
2708          @type q: any type that can be cast to a L{Data<escript.Data>} object on
2709                   L{Solution<escript.Solution>} or
2710                   L{ReducedSolution<escript.ReducedSolution>} depending on whether
2711                   reduced order is used for the representation of the equation
2712          @raise IllegalCoefficient: if an unknown coefficient keyword is used
2713          """
2714          super(TransportPDE,self).setValue(**coefficients)
2715    
2716       def createOperator(self):
2717           """
2718           Returns an instance of a new transport operator.
2719           """
2720           if self.useBackwardEuler():
2721             theta=1.
2722           else:
2723             theta=0.5
2724    
2725           self.trace("New Transport problem is allocated.")
2726           return self.getDomain().newTransportProblem( \
2727                                   theta,
2728                                   self.getNumEquations(), \
2729                                   self.getFunctionSpaceForSolution(), \
2730                                   self.getSystemType())
2731    
2732       def setInitialSolution(self,u):
2733           """
2734           Sets the initial solution.
2735    
2736           @param u: new initial solution
2737