/[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 122 by jgs, Thu Jun 9 05:38:05 2005 UTC trunk/escript/py_src/linearPDEs.py revision 2344 by jfenwick, Mon Mar 30 02:13:58 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        def __CompTuple2(self,t1,t2):
344          """
345          Compares two tuples of possible number of equations and number of
346          solutions.
347    
348     \f[        @param t1: the first tuple
349       -(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        @param t2: the second tuple
350     \f]        @return: 0, 1, or -1
351          """
352    
353     with boundary conditons:        dif=t1[0]+t1[1]-(t2[0]+t2[1])
354          if dif<0: return 1
355          elif dif>0: return -1
356          else: return 0
357    
358     \f[      def getShape(self,domain,numEquations=1,numSolutions=1):
359     n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_{ik}u_k = - n_j*X_{ij} + y_i         """
360     \f]         Builds the required shape of the coefficient.
361    
362     and contact conditions         @param domain: domain on which the PDE uses the coefficient
363           @type domain: L{Domain<escript.Domain>}
364           @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     \f[  class LinearProblem(object):
383     n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_contact_{ik}[u_k] = - n_j*X_{ij} + y_contact_i     """
384     \f]     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     and constraints:     """
433       DEFAULT= 0
434       DIRECT= 1
435       CHOLEVSKY= 2
436       PCG= 3
437       CR= 4
438       CGS= 5
439       BICGSTAB= 6
440       SSOR= 7
441       ILU0= 8
442       ILUT= 9
443       JACOBI= 10
444       GMRES= 11
445       PRES20= 12
446       LUMPING= 13
447       NO_REORDERING= 17
448       MINIMUM_FILL_IN= 18
449       NESTED_DISSECTION= 19
450       SCSL= 14
451       MKL= 15
452       UMFPACK= 16
453       ITERATIVE= 20
454       PASO= 21
455       AMG= 22
456       RILU = 23
457       TRILINOS = 24
458       NONLINEAR_GMRES = 25
459       TFQMR = 26
460       MINRES = 27
461       GS=28
462    
463       SMALL_TOLERANCE=1.e-13
464       PACKAGE_KEY="package"
465       METHOD_KEY="method"
466       SYMMETRY_KEY="symmetric"
467       TOLERANCE_KEY="tolerance"
468       PRECONDITIONER_KEY="preconditioner"
469    
470    
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     \f[       """
485     u_i=r_i \quad \mathrm{where} \quad q_i>0       super(LinearProblem, self).__init__()
    \f]  
   
    """  
    TOL=1.e-13  
    DEFAULT_METHOD=util.DEFAULT_METHOD  
    DIRECT=util.DIRECT  
    CHOLEVSKY=util.CHOLEVSKY  
    PCG=util.PCG  
    CR=util.CR  
    CGS=util.CGS  
    BICGSTAB=util.BICGSTAB  
    SSOR=util.SSOR  
    GMRES=util.GMRES  
    PRES20=util.PRES20  
    ILU0=util.ILU0  
    JACOBI=util.JACOBI  
   
    def __init__(self,domain,numEquations=0,numSolutions=0):  
      """  
      initializes a new linear PDE.  
   
      @param args:  
      """  
      # COEFFICIENTS can be overwritten by subclasses:  
      self.COEFFICIENTS={  
        "A"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM,PDECoefficient.SOLUTION,PDECoefficient.DIM),PDECoefficient.OPERATOR),  
        "B"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),  
        "C"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION,PDECoefficient.DIM),PDECoefficient.OPERATOR),  
        "D"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),  
        "X"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM),PDECoefficient.RIGHTHANDSIDE),  
        "Y"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),  
        "d"         : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),  
        "y"         : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),  
        "d_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),  
        "y_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),  
        "r"         : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),  
        "q"         : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.SOLUTION,),PDECoefficient.BOTH)}  
486    
487       # initialize attributes       self.__debug=debug
      self.__debug=None  
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           Sets the solution to zero.
1316         """         """
1317         returns the flux J_ij for a given u         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 applyOperator(self,u):     def resetRightHandSide(self):
1338         """         """
1339         applies the operator of the PDE to a given solution u in weak from         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         @param u: argument of the operator     def getCurrentRightHandSide(self):
1348         """         """
1349         return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())         Returns the right hand side in its current state.
                                                                                                                                                             
    def getResidual(self,u):  
1350         """         """
1351         return the residual of u in the weak from         if self.__righthandside.isEmpty(): self.__righthandside=self.createRightHandSide()
1352           return self.__righthandside
1353    
1354         @param u:     def resetOperator(self):
1355         """         """
1356         return self.applyOperator(u)-self.getRightHandSide()         Makes sure that the operator is instantiated and returns it initialized
1357           with zeros.
1358           """
1359           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     def __setValue(self,**coefficients):     def getCurrentOperator(self):
1372           """
1373           Returns the operator in its current state.
1374           """
1375           if self.__operator.isEmpty(): self.resetOperator()
1376           return self.__operator
1377    
1378       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()
1454         self.__rebuildSolution(deep)  
1455         self.__operator_isValid=False     def checkSymmetry(self,verbose=True):
1456         if deep: self.__operator=escript.Operator()        """
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 __rebuildRightHandSide(self,deep=False):     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())
1683           q=self.getCoefficientOfPDE("q")  
1684           r=self.getCoefficientOfPDE("r")     def checkSymmetry(self,verbose=True):
1685           if not q.isEmpty():        """
1686               if r.isEmpty():        Tests the PDE for symmetry.
1687                  r2=escript.Data(0,self.__righthandside.getShape(),self.__righthandside.getFunctionSpace())  
1688               else:        @param verbose: if set to True or not present a report on coefficients
1689                  r2=escript.Data(r,self.__righthandside.getFunctionSpace())                        which break the symmetry is printed.
1690               self.__righthandside.copyWithMask(r2,escript.Data(q,self.__righthandside.getFunctionSpace()))        @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 __applyConstraint(self):     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]}
2069         """  
2070         if not self.__solution_isValid:       or
2071             if self.debug() : print "PDE Debug: PDE is resolved."  
2072             self.__solution=self.solve(**options)       M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[l]+(B[j]+B_reduced[j])u-X[j]-X_reduced[j]}
            self.__solution_isValid=True  
        return self.__solution  
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    
2089    class Poisson(LinearPDE):
2090       """
2091       Class to define a Poisson equation problem. This is generally a
2092       L{LinearPDE} of the form
2093    
2094       M{-grad(grad(u)[j])[j] = f}
2095    
2096  def ELMAN_RAMAGE(P):     with natural boundary conditions
2097       """   """  
2098       return (P-1.).wherePositive()*0.5*(1.-1./(P+1.e-15))     M{n[j]*grad(u)[j] = 0 }
 def SIMPLIFIED_BROOK_HUGHES(P):  
      """   """  
      c=(P-3.).whereNegative()  
      return P/6.*c+1./2.*(1.-c)  
   
 def HALF(P):  
     """ """  
     return escript.Scalar(0.5,P.getFunctionSpace())  
   
 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.
2109                 A_max=util.Lsup(length_of_A)  
2110                 if A_max>0:       @param domain: domain of the PDE
2111                      inv_A=1./(length_of_A+A_max*self.TOL)       @type domain: L{Domain<escript.Domain>}
2112                 else:       @param debug: if True debug information is printed
2113                      inv_A=1./self.TOL  
2114                 peclet_number=length_of_Z*h/2*inv_A       """
2115                 xi=self.__xi(peclet_number)       super(Poisson, self).__init__(domain,1,1,debug)
2116                 self.__Xi=h*xi/(length_of_Z+Z_max*self.TOL)       self.introduceCoefficients(
2117                 print "@ preclet number = %e"%util.Lsup(peclet_number),util.Lsup(xi),util.Lsup(length_of_Z)                          f=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2118        return self.__Xi                          f_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE))
2119               self.setSymmetryOn()
2120    
2121     def getCoefficientOfPDE(self,name):     def setValue(self,**coefficients):
2122       """       """
2123       return the value of the coefficient name of the general PDE       Sets new values to coefficients.
2124    
2125       @param name:       @param coefficients: new values assigned to coefficients
2126       """       @keyword f: value for right hand side M{f}
2127       if not self.getNumEquations() == self.getNumSolutions():       @type f: any type that can be cast to a L{Scalar<escript.Scalar>} object
2128            raise ValueError,"AdvectivePDE expects the number of solution componets and the number of equations to be equal."                on L{Function<escript.Function>}
2129         @keyword q: mask for location of constraints
2130       if name == "A" :       @type q: any type that can be cast to a rank zero L{Data<escript.Data>}
2131           A=self.getCoefficient("A")                object on L{Solution<escript.Solution>} or
2132           B=self.getCoefficient("B")                L{ReducedSolution<escript.ReducedSolution>} depending on whether
2133           C=self.getCoefficient("C")                reduced order is used for the representation of the equation
2134           if B.isEmpty() and C.isEmpty():       @raise IllegalCoefficient: if an unknown coefficient keyword is used
2135              Aout=A       """
2136           else:       super(Poisson, self).setValue(**coefficients)
2137              if A.isEmpty():  
2138                 Aout=self.createNewCoefficient("A")  
2139              else:     def getCoefficient(self,name):
2140                 Aout=A[:]       """
2141              Xi=self.getXi()       Returns the value of the coefficient C{name} of the general PDE.
2142              if self.getNumEquations()>1:  
2143                  for i in range(self.getNumEquations()):       @param name: name of the coefficient requested
2144                     for j in range(self.getDim()):       @type name: C{string}
2145                        for k in range(self.getNumSolutions()):       @return: the value of the coefficient C{name}
2146                           for l in range(self.getDim()):       @rtype: L{Data<escript.Data>}
2147                              if not C.isEmpty() and not B.isEmpty():       @raise IllegalCoefficient: invalid coefficient name
2148                                 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])       @note: This method is called by the assembling routine to map the Poisson
2149                              elif C.isEmpty():              equation onto the general PDE.
2150                                 for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*B[p,j,i]*B[p,l,k]       """
2151                              else:       if name == "A" :
2152                                 for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*C[p,i,j]*C[p,k,l]           return escript.Data(util.kronecker(self.getDim()),escript.Function(self.getDomain()))
2153              else:       elif name == "Y" :
2154                  for j in range(self.getDim()):           return self.getCoefficient("f")
2155                     for l in range(self.getDim()):       elif name == "Y_reduced" :
2156                        if not C.isEmpty() and not B.isEmpty():           return self.getCoefficient("f_reduced")
2157                            Aout[j,l]+=Xi*(C[j]-B[j])*(C[l]-B[l])       else:
2158                        elif C.isEmpty():           return super(Poisson, self).getCoefficient(name)
2159                            Aout[j,l]+=Xi*B[j]*B[l]  
2160                        else:  class Helmholtz(LinearPDE):
2161                            Aout[j,l]+=Xi*C[j]*C[l]     """
2162           return Aout     Class to define a Helmholtz equation problem. This is generally a
2163       elif name == "B" :     L{LinearPDE} of the form
2164           B=self.getCoefficient("B")  
2165           C=self.getCoefficient("C")     M{S{omega}*u - grad(k*grad(u)[j])[j] = f}
2166           D=self.getCoefficient("D")  
2167           if C.isEmpty() or D.isEmpty():     with natural boundary conditions
2168              Bout=B  
2169           else:     M{k*n[j]*grad(u)[j] = g- S{alpha}u }
2170              Xi=self.getXi()  
2171              if B.isEmpty():     and constraints:
2172                  Bout=self.createNewCoefficient("B")  
2173              else:     M{u=r} where M{q>0}
2174                  Bout=B[:]  
2175              if self.getNumEquations()>1:     """
2176                 for k in range(self.getNumSolutions()):  
2177                    for p in range(self.getNumEquations()):     def __init__(self,domain,debug=False):
2178                       tmp=Xi*D[p,k]       """
2179                       for i in range(self.getNumEquations()):       Initializes a new Helmholtz equation.
2180                          for j in range(self.getDim()):  
2181                             Bout[i,j,k]+=tmp*C[p,i,j]       @param domain: domain of the PDE
2182              else:       @type domain: L{Domain<escript.Domain>}
2183                 tmp=Xi*D       @param debug: if True debug information is printed
2184                 for j in range(self.getDim()): Bout[j]+=tmp*C[j]  
2185           return Bout       """
2186       elif name == "C" :       super(Helmholtz, self).__init__(domain,1,1,debug)
2187           B=self.getCoefficient("B")       self.introduceCoefficients(
2188           C=self.getCoefficient("C")                          omega=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.OPERATOR),
2189           D=self.getCoefficient("D")                          k=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.OPERATOR),
2190           if B.isEmpty() or D.isEmpty():                          f=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2191              Cout=C                          f_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2192           else:                          alpha=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.OPERATOR),
2193              Xi=self.getXi()                          g=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2194              if C.isEmpty():                          g_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE))
2195                  Cout=self.createNewCoefficient("C")       self.setSymmetryOn()
2196              else:  
2197                  Cout=C[:]     def setValue(self,**coefficients):
2198              if self.getNumEquations()>1:       """
2199                 for k in range(self.getNumSolutions()):       Sets new values to coefficients.
2200                     for p in range(self.getNumEquations()):  
2201                        tmp=Xi*D[p,k]       @param coefficients: new values assigned to coefficients
2202                        for i in range(self.getNumEquations()):       @keyword omega: value for coefficient M{S{omega}}
2203                          for l in range(self.getDim()):       @type omega: any type that can be cast to a L{Scalar<escript.Scalar>}
2204                                   Cout[i,k,l]+=tmp*B[p,l,i]                    object on L{Function<escript.Function>}
2205              else:       @keyword k: value for coefficient M{k}
2206                 tmp=Xi*D       @type k: any type that can be cast to a L{Scalar<escript.Scalar>} object
2207                 for j in range(self.getDim()): Cout[j]+=tmp*B[j]                on L{Function<escript.Function>}
2208           return Cout       @keyword f: value for right hand side M{f}
2209       elif name == "D" :       @type f: any type that can be cast to a L{Scalar<escript.Scalar>} object
2210           return self.getCoefficient("D")                on L{Function<escript.Function>}
2211       elif name == "X" :       @keyword alpha: value for right hand side M{S{alpha}}
2212           X=self.getCoefficient("X")       @type alpha: any type that can be cast to a L{Scalar<escript.Scalar>}
2213           Y=self.getCoefficient("Y")                    object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}
2214           B=self.getCoefficient("B")       @keyword g: value for right hand side M{g}
2215           C=self.getCoefficient("C")       @type g: any type that can be cast to a L{Scalar<escript.Scalar>} object
2216           if Y.isEmpty() or (B.isEmpty() and C.isEmpty()):                on L{FunctionOnBoundary<escript.FunctionOnBoundary>}
2217              Xout=X       @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:           else:
2245              if X.isEmpty():                return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))*self.getCoefficient("k")
2246                  Xout=self.createNewCoefficient("X")       elif name == "D" :
2247             return self.getCoefficient("omega")
2248         elif name == "Y" :
2249             return self.getCoefficient("f")
2250         elif name == "d" :
2251             return self.getCoefficient("alpha")
2252         elif name == "y" :
2253             return self.getCoefficient("g")
2254         elif name == "Y_reduced" :
2255             return self.getCoefficient("f_reduced")
2256         elif name == "y_reduced" :
2257            return self.getCoefficient("g_reduced")
2258         else:
2259            return super(Helmholtz, self).getCoefficient(name)
2260    
2261    class LameEquation(LinearPDE):
2262       """
2263       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       with natural boundary conditions:
2268    
2269       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] }
2270    
2271       and constraints:
2272    
2273       M{u[i]=r[i]} where M{q[i]>0}
2274    
2275       """
2276    
2277       def __init__(self,domain,debug=False):
2278          """
2279          Initializes a new Lame equation.
2280    
2281          @param domain: domain of the PDE
2282          @type domain: L{Domain<escript.Domain>}
2283          @param debug: if True debug information is printed
2284    
2285          """
2286          super(LameEquation, self).__init__(domain,\
2287                                             domain.getDim(),domain.getDim(),debug)
2288          self.introduceCoefficients(lame_lambda=PDECoef(PDECoef.INTERIOR,(),PDECoef.OPERATOR),
2289                                     lame_mu=PDECoef(PDECoef.INTERIOR,(),PDECoef.OPERATOR),
2290                                     F=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2291                                     sigma=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
2292                                     f=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE))
2293          self.setSymmetryOn()
2294    
2295       def setValues(self,**coefficients):
2296         """
2297         Sets new values to coefficients.
2298    
2299         @param coefficients: new values assigned to coefficients
2300         @keyword lame_mu: value for coefficient M{S{mu}}
2301         @type lame_mu: any type that can be cast to a L{Scalar<escript.Scalar>}
2302                        object on L{Function<escript.Function>}
2303         @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                            object on L{Function<escript.Function>}
2306         @keyword F: value for internal force M{F}
2307         @type F: any type that can be cast to a L{Vector<escript.Vector>} object
2308                  on L{Function<escript.Function>}
2309         @keyword sigma: value for initial stress M{S{sigma}}
2310         @type sigma: any type that can be cast to a L{Tensor<escript.Tensor>}
2311                      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:              else:
2345                  Xout=X[:]                  for i in range(self.getDim()):
2346              Xi=self.getXi()                    for j in range(self.getDim()):
2347              if self.getNumEquations()>1:                      out[i,j,j,i] += self.getCoefficient("lame_mu")
2348                   for p in range(self.getNumEquations()):                      out[i,j,i,j] += self.getCoefficient("lame_mu")
2349                      tmp=Xi*Y[p]           else:
2350                      for i in range(self.getNumEquations()):              if self.getCoefficient("lame_mu").isEmpty():
2351                         for j in range(self.getDim()):                  for i in range(self.getDim()):
2352                            if not C.isEmpty() and not B.isEmpty():                    for j in range(self.getDim()):
2353                               Xout[i,j]+=tmp*(C[p,i,j]-B[p,j,i])                      out[i,i,j,j] += self.getCoefficient("lame_lambda")
                           elif C.isEmpty():  
                              Xout[i,j]-=tmp*B[p,j,i]  
                           else:  
                              Xout[i,j]+=tmp*C[p,i,j]  
2354              else:              else:
2355                   tmp=Xi*Y                  for i in range(self.getDim()):
2356                   for j in range(self.getDim()):                    for j in range(self.getDim()):
2357                      if not C.isEmpty() and not B.isEmpty():                      out[i,i,j,j] += self.getCoefficient("lame_lambda")
2358                         Xout[j]+=tmp*(C[j]-B[j])                      out[i,j,j,i] += self.getCoefficient("lame_mu")
2359                      elif C.isEmpty():                      out[i,j,i,j] += self.getCoefficient("lame_mu")
2360                         Xout[j]-=tmp*B[j]           return out
2361                      else:       elif name == "X" :
2362                         Xout[j]+=tmp*C[j]           return self.getCoefficient("sigma")
2363           return Xout       elif name == "Y" :
2364       elif name == "Y" :           return self.getCoefficient("F")
2365           return self.getCoefficient("Y")       elif name == "y" :
2366       elif name == "d" :           return self.getCoefficient("f")
          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")  
2367       else:       else:
2368           raise SystemError,"unknown PDE coefficient %s",name          return super(LameEquation, self).getCoefficient(name)
2369    
2370    def LinearSinglePDE(domain,debug=False):
2371       """
2372       Defines a single linear PDE.
2373    
2374  class Poisson(LinearPDE):     @param domain: domain of the PDE
2375       @type domain: L{Domain<escript.Domain>}
2376       @param debug: if True debug information is printed
2377       @rtype: L{LinearPDE}
2378     """     """
2379     Class to define a Poisson equstion problem:     return LinearPDE(domain,numEquations=1,numSolutions=1,debug=debug)
2380    
2381     class to define a linear PDE of the form  def LinearPDESystem(domain,debug=False):
2382     \f[     """
2383     -u_{,jj} = f     Defines a system of linear PDEs.
    \f]  
   
    with boundary conditons:  
   
    \f[  
    n_j*u_{,j} = 0  
    \f]  
2384    
2385     and constraints:     @param domain: domain of the PDEs
2386       @type domain: L{Domain<escript.Domain>}
2387       @param debug: if True debug information is printed
2388       @rtype: L{LinearPDE}
2389       """
2390       return LinearPDE(domain,numEquations=domain.getDim(),numSolutions=domain.getDim(),debug=debug)
2391    
2392     \f[  
2393     u=0 \quad \mathrm{where} \quad q>0  class TransportPDE(LinearProblem):
2394     \f]     """
2395     """     This class is used to define a transport problem given by a general linear,
2396       time dependent, second order PDE for an unknown, non-negative,
2397     def __init__(self,domain,f=escript.Data(),q=escript.Data()):     time-dependent function M{u} on a given domain defined through a
2398         LinearPDE.__init__(self,domain,1,1)     L{Domain<escript.Domain>} object.
2399         self.COEFFICIENTS={  
2400         "f"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),     For a single equation with a solution with a single component the transport
2401         "q"         : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.BOTH)}     problem is defined in the following form:
2402         self.setSymmetryOn()  
2403         self.setValue(f,q)     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     def setValue(self,f=escript.Data(),q=escript.Data()):     where M{u_t} denotes the time derivative of M{u} and M{grad(F)} denotes the
2406         self._LinearPDE__setValue(f=f,q=q)     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     def getCoefficientOfPDE(self,name):     The coefficients M{M}, M{A}, M{B}, M{C}, M{D}, M{X} and M{Y} have to be
2409       """     specified through L{Data<escript.Data>} objects in L{Function<escript.Function>}
2410       return the value of the coefficient name of the general PDE     and the coefficients M{M_reduced}, M{A_reduced}, M{B_reduced}, M{C_reduced},
2411       M{D_reduced}, M{X_reduced} and M{Y_reduced} have to be specified through
2412       @param name:     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       if name == "A" :     L{Data<escript.Data>} objects. M{M} and M{M_reduced} are scalar, M{A} and
2415           return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))     M{A_reduced} are rank two, M{B}, M{C}, M{X}, M{B_reduced}, M{C_reduced} and
2416       elif name == "B" :     M{X_reduced} are rank one and M{D}, M{D_reduced}, M{Y} and M{Y_reduced} are
2417           return escript.Data()     scalar.
2418       elif name == "C" :  
2419           return escript.Data()     The following natural boundary conditions are considered:
2420       elif name == "D" :  
2421           return escript.Data()     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       elif name == "X" :  
2423           return escript.Data()     where M{n} is the outer normal field. Notice that the coefficients M{A},
2424       elif name == "Y" :     M{A_reduced}, M{B}, M{B_reduced}, M{X} and M{X_reduced} are defined in the
2425           return self.getCoefficient("f")     transport problem. The coefficients M{m}, M{d} and M{y} are each a scalar in
2426       elif name == "d" :     L{FunctionOnBoundary<escript.FunctionOnBoundary>} and the coefficients
2427           return escript.Data()     M{m_reduced}, M{d_reduced} and M{y_reduced} are each a scalar in
2428       elif name == "y" :     L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
2429           return escript.Data()  
2430       elif name == "d_contact" :     Constraints for the solution prescribing the value of the solution at
2431           return escript.Data()     certain locations in the domain have the form
2432       elif name == "y_contact" :  
2433           return escript.Data()     M{u_t=r} where M{q>0}
2434       elif name == "r" :  
2435           return escript.Data()     M{r} and M{q} are each scalar where M{q} is the characteristic function
2436       elif name == "q" :     defining where the constraint is applied. The constraints override any other
2437           return self.getCoefficient("q")     condition set by the transport problem or the boundary condition.
2438    
2439       The transport problem is symmetrical if
2440    
2441       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       M{B_reduced[j]=C_reduced[j]}
2443    
2444       For a system and a solution with several components the transport problem
2445       has the form
2446    
2447       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    
2449       M{A} and M{A_reduced} are of rank four, M{B}, M{B_reduced}, M{C} and
2450       M{C_reduced} are each of rank three, M{M}, M{M_reduced}, M{D}, M{D_reduced},
2451       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    
2454       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       The coefficient M{d} and M{m} are of rank two and M{y} is of rank one with
2457       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       one all in L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
2460    
2461       Constraints take the form
2462    
2463       M{u[i]_t=r[i]} where M{q[i]>0}
2464    
2465       M{r} and M{q} are each rank one. Notice that at some locations not
2466       necessarily all components must have a constraint.
2467    
2468       The transport problem is symmetrical if
2469    
2470          - M{M[i,k]=M[i,k]}
2471          - M{M_reduced[i,k]=M_reduced[i,k]}
2472          - 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          - M{B[i,j,k]=C[k,i,j]}
2475          - M{B_reduced[i,j,k]=C_reduced[k,i,j]}
2476          - M{D[i,k]=D[i,k]}
2477          - M{D_reduced[i,k]=D_reduced[i,k]}
2478          - M{m[i,k]=m[k,i]}
2479          - M{m_reduced[i,k]=m_reduced[k,i]}
2480          - M{d[i,k]=d[k,i]}
2481          - M{d_reduced[i,k]=d_reduced[k,i]}
2482    
2483       L{TransportPDE} also supports solution discontinuities over a contact region
2484       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       several components of the solution, is defined as
2487    
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    
2490       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:       else:
2550           raise SystemError,"unknown PDE coefficient %s",name           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           @type u: any object that can be interpolated to a L{Data<escript.Data>}
2738                    object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2739           @note: C{u} must be non-negative
2740           """
2741           u2=util.interpolate(u,self.getFunctionSpaceForSolution())
2742           if self.getNumSolutions() == 1:
2743              if u2.getShape()!=():
2744                  raise ValueError,"Illegal shape %s of initial solution."%(u2.getShape(),)
2745           else:
2746              if u2.getShape()!=(self.getNumSolutions(),):
2747                  raise ValueError,"Illegal shape %s of initial solution."%(u2.getShape(),)
2748           self.getOperator().setInitialValue(u2)
2749    
2750       def getRequiredSystemType(self):
2751          """
2752          Returns the system type which needs to be used by the current set up.
2753    
2754          @return: a code to indicate the type of transport problem scheme used
2755          @rtype: C{float}
2756          """
2757          return self.getDomain().getTransportTypeId(self.getSolverMethod()[0],self.getSolverMethod()[1],self.getSolverPackage(),self.isSymmetric())
2758    
2759       def getUnlimitedTimeStepSize(self):
2760          """
2761          Returns the value returned by the C{getSafeTimeStepSize} method to
2762          indicate no limit on the safe time step size.
2763    
2764           @return: the value used to indicate that no limit is set to the time
2765                    step size
2766           @rtype: C{float}
2767           @note: Typically the biggest positive float is returned
2768          """
2769          return self.getOperator().getUnlimitedTimeStepSize()
2770    
2771       def getSafeTimeStepSize(self):
2772           """
2773           Returns a safe time step size to do the next time step.
2774    
2775           @return: safe time step size
2776           @rtype: C{float}
2777           @note: If not C{getSafeTimeStepSize()} < C{getUnlimitedTimeStepSize()}
2778                  any time step size can be used.
2779           """
2780           return self.getOperator().getSafeTimeStepSize()
2781    
2782       def setConstraintWeightingFactor(self,value=1./util.sqrt(util.EPSILON)):
2783           """
2784           Sets the weighting factor used to insert the constraints into the problem
2785    
2786           @param value: value for the weighting factor
2787           @type value: large positive C{float}
2788           """
2789           if not value>0:
2790             raise ValueError,"weighting factor needs to be positive."
2791           self.__constraint_factor=value
2792           self.trace("Weighting factor for constraints is set to %e."%value)
2793    
2794       def getConstraintWeightingFactor(self):
2795           """
2796           returns the weighting factor used to insert the constraints into the problem
2797           @return: value for the weighting factor
2798           @rtype: C{float}
2799           """
2800           return self.__constraint_factor
2801       #====================================================================
2802       def getSolution(self,dt,**options):
2803           """