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

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

  ViewVC Help
Powered by ViewVC 1.1.26