/[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 104 by jgs, Fri Dec 17 07:43:12 2004 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  @brief 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  def identifyDomain(domain=None,data={}):  __author__="Lutz Gross, l.gross@uq.edu.au"
      """  
      @brief Return the Domain which is equal to the input domain (if not None)  
      and is the domain of all Data objects in the dictionary data.  
      An exception is raised if this is not possible  
   
      @param domain  
      @param data  
      """  
      # get the domain used by any Data object in the list data:  
      data_domain=None  
      for d in data.itervalues():  
           if isinstance(d,escript.Data):  
              if not d.isEmpty(): data_domain=d.getDomain()  
      # check if domain and data_domain are identical?  
      if domain == None:  
          if data_domain == None:  
               raise ValueError,"Undefined PDE domain. Specify a domain or use a Data class object as coefficient"  
      else:  
          if data_domain == None:  
               data_domain=domain  
          else:  
            if not data_domain == domain:  
                  raise ValueError,"Domain of coefficients doesnot match specified domain"  
      # now we check if all Data class object coefficients are defined on data_domain:  
      for i,d in data.iteritems():  
          if isinstance(d,escript.Data):  
             if not d.isEmpty():  
                if not data_domain==d.getDomain():  
                  raise ValueError,"Illegal domain for coefficient %s."%i  
      # done:  
      return data_domain  
   
 def identifyNumEquationsAndSolutions(dim,coef={}):  
      # get number of equations and number of unknowns:  
      numEquations=0  
      numSolutions=0  
      for i in coef.iterkeys():  
         if not coef[i].isEmpty():  
            res=_PDECoefficientTypes[i].estimateNumEquationsAndNumSolutions(coef[i].getShape(),dim)  
            if res==None:  
                raise ValueError,"Illegal shape %s of coefficient %s"%(coef[i].getShape().__str__(),i)  
            else:  
                numEquations=max(numEquations,res[0])  
                numSolutions=max(numSolutions,res[1])  
      return numEquations,numSolutions  
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       pass
51    
52    class IllegalCoefficientValue(ValueError):
53       """
54       raised if an incorrect value for a coefficient is used.
55       """
56       pass
57    
58    class IllegalCoefficientFunctionSpace(ValueError):
59       """
60       raised if an incorrect function space for a coefficient is used.
61     """     """
    @brief  
62    
63     @param t1  class UndefinedPDEError(ValueError):
    @param t2  
64     """     """
65     dif=t1[0]+t1[1]-(t2[0]+t2[1])     raised if a PDE is not fully defined yet.
66     if dif<0: return 1     """
67     elif dif>0: return -1     pass
    else: return 0  
68    
69  class PDECoefficientType:  class PDECoef(object):
70      """      """
71      @brief      A class for describing a PDE coefficient
72    
73        @cvar INTERIOR: indicator that coefficient is defined on the interior of the PDE domain
74        @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    
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         @brief Initialise a PDE Coefficient type         """
106           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
128           self.resetValue()
129    
130      def getFunctionSpace(self,domain):      def resetValue(self):
131         """         """
132         @brief defines the FunctionSpace of the coefficient on the domain         resets coefficient value to default
   
        @param domain  
133         """         """
134         if self.what==self.INTERIOR: return escript.Function(domain)         self.value=escript.Data()
135         elif self.what==self.BOUNDARY: return escript.FunctionOnBoundary(domain)  
136         elif self.what==self.CONTACT: return escript.FunctionOnContactZero(domain)      def getFunctionSpace(self,domain,reducedEquationOrder=False,reducedSolutionOrder=False):
137         elif self.what==self.CONTINUOUS: return escript.ContinuousFunction(domain)         """
138           defines the L{FunctionSpace<escript.FunctionSpace>} of the coefficient
139    
140           @param domain: domain on which the PDE uses the coefficient
141           @type domain: L{Domain<escript.Domain>}
142           @param reducedEquationOrder: True to indicate that reduced order is used to represent the equation
143           @type reducedEquationOrder: C{bool}
144           @param reducedSolutionOrder: True to indicate that reduced order is used to represent the solution
145           @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):
170           """
171           returns the value of the coefficient
172    
173           @return:  value of the coefficient
174           @rtype:  L{Data<escript.Data>}
175           """
176           return self.value
177    
178        def setValue(self,domain,numEquations=1,numSolutions=1,reducedEquationOrder=False,reducedSolutionOrder=False,newValue=None):
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
212    
213      def isAlteringOperator(self):      def isAlteringOperator(self):
214          """          """
215      @brief 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 119  class PDECoefficientType: Line 224  class PDECoefficientType:
224    
225      def isAlteringRightHandSide(self):      def isAlteringRightHandSide(self):
226          """          """
227      @brief 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         @brief 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         return None         elif self.definesNumEquation():
267              for e in range(num,0,-1):
268      def buildShape(self,e=1,u=1,dim=3):               s=self.getShape(domain,e,0)
269          """               if len(s)==0 and len(shape)==0:
270      @brief builds the required shape for a given number of equations e, number of unknowns u and spatial dimension dim                   return (1,None)
271                 else:
272                     if s==shape: return (e,None)
273    
274      @param e         elif self.definesNumSolutions():
275      @param u            for u in range(num,0,-1):
276      @param dim               s=self.getShape(domain,0,u)
277      """               if len(s)==0 and len(shape)==0:
278          s=()                   return (None,1)
         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,)  
279               else:               else:
280                  s=s+(dim,)                   if s==shape: return (None,u)
281          return s         return None
282        def definesNumSolutions(self):
283           """
284           checks if the coefficient allows to estimate the number of solution components
285    
286  _PDECoefficientTypes={         @return: True if the coefficient allows an estimate of the number of solution components
287  "A"         : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,PDECoefficientType.DIM,PDECoefficientType.SOLUTION,PDECoefficientType.DIM),PDECoefficientType.OPERATOR),         @rtype: C{bool}
288  "B"         : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,PDECoefficientType.DIM,PDECoefficientType.SOLUTION),PDECoefficientType.OPERATOR),         """
289  "C"         : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,PDECoefficientType.SOLUTION,PDECoefficientType.DIM),PDECoefficientType.OPERATOR),         for i in self.pattern:
290  "D"         : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,PDECoefficientType.SOLUTION),PDECoefficientType.OPERATOR),               if i==self.BY_SOLUTION: return True
291  "X"         : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,PDECoefficientType.DIM),PDECoefficientType.RIGHTHANDSIDE),         return False
 "Y"         : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,),PDECoefficientType.RIGHTHANDSIDE),  
 "d"         : PDECoefficientType(PDECoefficientType.BOUNDARY,(PDECoefficientType.EQUATION,PDECoefficientType.SOLUTION),PDECoefficientType.OPERATOR),  
 "y"         : PDECoefficientType(PDECoefficientType.BOUNDARY,(PDECoefficientType.EQUATION,),PDECoefficientType.RIGHTHANDSIDE),  
 "d_contact" : PDECoefficientType(PDECoefficientType.CONTACT,(PDECoefficientType.EQUATION,PDECoefficientType.SOLUTION),PDECoefficientType.OPERATOR),  
 "y_contact" : PDECoefficientType(PDECoefficientType.CONTACT,(PDECoefficientType.EQUATION,),PDECoefficientType.RIGHTHANDSIDE),  
 "r"         : PDECoefficientType(PDECoefficientType.CONTINUOUS,(PDECoefficientType.EQUATION,),PDECoefficientType.RIGHTHANDSIDE),  
 "q"         : PDECoefficientType(PDECoefficientType.CONTINUOUS,(PDECoefficientType.SOLUTION,),PDECoefficientType.BOTH),  
 }  
292    
293  class LinearPDE:      def definesNumEquation(self):
294     """         """
295     @brief Class to define a linear PDE         checks if the coefficient allows to estimate the number of equations
     
    class to define a linear PDE of the form  
296    
297       -(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         @return: True if the coefficient allows an estimate of the number of equations
298           @rtype: C{bool}
299           """
300           for i in self.pattern:
301                 if i==self.BY_EQUATION: return True
302           return False
303    
304       with boundary conditons:      def __CompTuple2(self,t1,t2):
305          """
306          Compare two tuples of possible number of equations and number of solutions
307    
308          n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_{ik}u_k = - n_j*X_{ij} + y_i        @param t1: The first tuple
309          @param t2: The second tuple
310    
311      and contact conditions        """
312    
313          n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_contact_{ik}[u_k] = - n_j*X_{ij} + y_contact_i        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      and constraints:      def getShape(self,domain,numEquations=1,numSolutions=1):
319           """
320           builds the required shape of the coefficient
321    
322           u_i=r_i where q_i>0         @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    class LinearProblem(object):
343     """     """
344     DEFAULT_METHOD=util.DEFAULT_METHOD     This is the base class is to define a general linear PDE-type problem for an u
345     DIRECT=util.DIRECT     for an unknown function M{u} on a given domain defined through a L{Domain<escript.Domain>} object.
346     CHOLEVSKY=util.CHOLEVSKY     The problem can be given as a single equations or as a system of equations.
347     PCG=util.PCG  
348     CR=util.CR     The class assumes that some sort of assembling process is required to form a problem of the form
349     CGS=util.CGS  
350     BICGSTAB=util.BICGSTAB     M{L u=f}
351     SSOR=util.SSOR    
352     GMRES=util.GMRES     where M{L} is an operator and M{f} is the right hand side. This operator problem will be solved to get the
353     PRES20=util.PRES20     unknown M{u}.
354    
355       @cvar DEFAULT: The default method used to solve the system of linear equations
356       @cvar DIRECT: The direct solver based on LDU factorization
357       @cvar CHOLEVSKY: The direct solver based on LDLt factorization (can only be applied for symmetric PDEs)
358       @cvar PCG: The preconditioned conjugate gradient method (can only be applied for symmetric PDEs)
359       @cvar CR: The conjugate residual method
360       @cvar CGS: The conjugate gardient square method
361       @cvar BICGSTAB: The stabilized BiConjugate Gradient method.
362       @cvar TFQMR: Transport Free Quasi Minimal Residual method.
363       @cvar MINRES: Minimum residual method.
364       @cvar SSOR: The symmetric overrealaxtion method
365       @cvar ILU0: The incomplete LU factorization preconditioner  with no fill in
366       @cvar ILUT: The incomplete LU factorization preconditioner with will in
367       @cvar JACOBI: The Jacobi preconditioner
368       @cvar GMRES: The Gram-Schmidt minimum residual method
369       @cvar PRES20: Special GMRES with restart after 20 steps and truncation after 5 residuals
370       @cvar LUMPING: Matrix lumping.
371       @cvar NO_REORDERING: No matrix reordering allowed
372       @cvar MINIMUM_FILL_IN: Reorder matrix to reduce fill-in during factorization
373       @cvar NESTED_DISSECTION: Reorder matrix to improve load balancing during factorization
374       @cvar PASO: PASO solver package
375       @cvar SCSL: SGI SCSL solver library
376       @cvar MKL: Intel's MKL solver library
377       @cvar UMFPACK: the UMFPACK library
378       @cvar TRILINOS: the TRILINOS parallel solver class library from Sandia Natl Labs
379       @cvar ITERATIVE: The default iterative solver
380       @cvar AMG: algebraic multi grid
381       @cvar RILU: recursive ILU
382       @cvar GS: Gauss-Seidel solver
383    
384     def __init__(self,domain,numEquations=0,numSolutions=0):     """
385       """     DEFAULT= 0
386       @brief initializes a new linear PDE.     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    
      @param args  
435       """       """
436         super(LinearProblem, self).__init__()
437    
438       # initialize attributes       self.__debug=debug
      self.__debug=None  
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         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         switches on debugging
474         """
475         self.__debug=not None
476    
477     def getCoefficient(self,name):     def setDebugOff(self):
478       """       """
479       @brief return the value of the coefficient name       switches off debugging
480         """
481         self.__debug=None
482    
483       @param name     def trace(self,text):
484       """       """
485       return self.__coefficient[name]       print the text message if debugging is swiched on.
486         @param text: message
487         @type text: C{string}
488         """
489         if self.__debug: print "%s: %s"%(str(self),text)
490    
491     def setValue(self,**coefficients):     # =============================================================================
492        """     # some service functions:
493        @brief sets new values to coefficients     # =============================================================================
494       def introduceCoefficients(self,**coeff):
495           """
496           introduces a new coefficient into the problem.
497    
498        @param coefficients         use
       """  
       self._setValue(**coefficients)  
         
499    
500     def _setValue(self,**coefficients):           self.introduceCoefficients(A=PDECoef(...),B=PDECoef(...))
       """  
       @brief sets new values to coefficients  
501    
502        @param coefficients         to introduce the coefficients M{A} ans M{B}.
503        """         """
504                 for name, type in coeff.items():
505        # get the dictionary of the coefficinets been altered:             if not isinstance(type,PDECoef):
506        alteredCoefficients={}                raise ValueError,"coefficient %s has no type."%name
507        for i,d in coefficients.iteritems():             self.__COEFFICIENTS[name]=type
508           if self.hasCoefficient(i):             self.__COEFFICIENTS[name].resetValue()
509              if d == None:             self.trace("coefficient %s has been introduced."%name)
                 alteredCoefficients[i]=escript.Data()  
             elif isinstance(d,escript.Data):  
                 if d.isEmpty():  
                   alteredCoefficients[i]=escript.Data()  
                 else:  
                   alteredCoefficients[i]=escript.Data(d,self.getFunctionSpaceOfCoefficient(i))  
             else:  
                 alteredCoefficients[i]=escript.Data(d,self.getFunctionSpaceOfCoefficient(i))  
          else:  
             raise ValueError,"Attempt to set undefined coefficient %s"%i  
       # if numEquations and numSolutions is undefined we try identify their values based on the coefficients:  
       if self.__numEquations<1 or self.__numSolutions<1:  
             numEquations,numSolutions=identifyNumEquationsAndSolutions(self.getDomain().getDim(),alteredCoefficients)  
             if self.__numEquations<1 and numEquations>0: self.__numEquations=numEquations  
             if self.__numSolutions<1 and numSolutions>0: self.__numSolutions=numSolutions  
             if self.debug() and self.__numEquations>0: print "PDE Debug: identified number of equations is ",self.__numEquations  
             if self.debug() and self.__numSolutions>0: print "PDE Debug: identified number of solutions is ",self.__numSolutions  
510    
511        # now we check the shape of the coefficient if numEquations and numSolutions are set:     def getDomain(self):
512        if  self.__numEquations>0 and  self.__numSolutions>0:       """
513           for i in self.__coefficient.iterkeys():       returns the domain of the PDE
              if alteredCoefficients.has_key(i) and not alteredCoefficients[i].isEmpty():  
                  if not self.getShapeOfCoefficient(i)==alteredCoefficients[i].getShape():  
                     raise ValueError,"Expected shape for coefficient %s is %s but actual shape is %s."%(i,self.getShapeOfCoefficient(i),alteredCoefficients[i].getShape())  
              else:  
                  if not self.__coefficient[i].isEmpty():  
                     if not self.getShapeOfCoefficient(i)==self.__coefficient[i].getShape():  
                        raise ValueError,"Expected shape for coefficient %s is %s but actual shape is %s."%(i,self.getShapeOfCoefficient(i),self.__coefficient[i].getShape())  
       # overwrite new values:  
       for i,d in alteredCoefficients.iteritems():  
          if self.debug(): print "PDE Debug: Coefficient %s has been altered."%i  
          self.__coefficient[i]=d  
          self.alteredCoefficient(i)  
514    
515        # reset the HomogeneousConstraintFlag:       @return: the domain of the PDE
516        self.__setHomogeneousConstraintFlag()       @rtype: L{Domain<escript.Domain>}
517        if not "q" in alteredCoefficients and not self.__homogeneous_constraint: self.__rebuildSystem()       """
518         return self.__domain
519    
520     def cleanCoefficients(self):     def getDim(self):
521       """       """
522       @brief resets all coefficients to default values.       returns the spatial dimension of the PDE
523    
524         @return: the spatial dimension of the PDE domain
525         @rtype: C{int}
526       """       """
527       self.__coefficient={}       return self.getDomain().getDim()
      for i in _PDECoefficientTypes.iterkeys():  
          self.__coefficient[i]=escript.Data()  
528    
529     def getShapeOfCoefficient(self,name):     def getNumEquations(self):
530       """       """
531       @brief return the shape of the coefficient name       returns the number of equations
532    
533       @param name       @return: the number of equations
534         @rtype: C{int}
535         @raise UndefinedPDEError: if the number of equations is not be specified yet.
536       """       """
537       if self.hasCoefficient(name):       if self.__numEquations==None:
538          return _PDECoefficientTypes[name].buildShape(self.getNumEquations(),self.getNumSolutions(),self.getDomain().getDim())           raise UndefinedPDEError,"Number of equations is undefined. Please specify argument numEquations."
539       else:       else:
540          raise ValueError,"Solution coefficient %s requested"%name           return self.__numEquations
541    
542     def getFunctionSpaceOfCoefficient(self,name):     def getNumSolutions(self):
543       """       """
544       @brief return the atoms of the coefficient name       returns the number of unknowns
545    
546       @param name       @return: the number of unknowns
547         @rtype: C{int}
548         @raise UndefinedPDEError: if the number of unknowns is not be specified yet.
549       """       """
550       if self.hasCoefficient(name):       if self.__numSolutions==None:
551          return _PDECoefficientTypes[name].getFunctionSpace(self.getDomain())          raise UndefinedPDEError,"Number of solution is undefined. Please specify argument numSolutions."
552       else:       else:
553          raise ValueError,"Solution coefficient %s requested"%name          return self.__numSolutions
554    
555     def alteredCoefficient(self,name):     def reduceEquationOrder(self):
556       """       """
557       @brief annonced that coefficient name has been changed       return status for order reduction for equation
558    
559       @param name       @return: return True is reduced interpolation order is used for the represenation of the equation
560         @rtype: L{bool}
561       """       """
562       if self.hasCoefficient(name):       return self.__reduce_equation_order
         if _PDECoefficientTypes[name].isAlteringOperator(): self.__rebuildOperator()  
         if _PDECoefficientTypes[name].isAlteringRightHandSide(): self.__rebuildRightHandSide()  
      else:  
         raise ValueError,"Solution coefficient %s requested"%name  
563    
564     def __setHomogeneousConstraintFlag(self):     def reduceSolutionOrder(self):
565        """       """
566        @brief checks if the constraints are homogeneous and sets self.__homogeneous_constraint accordingly.       return status for order reduction for the solution
       """  
       self.__homogeneous_constraint=True  
       q=self.getCoefficient("q")  
       r=self.getCoefficient("r")  
       if not q.isEmpty() and not r.isEmpty():  
          print (q*r).Lsup(), 1.e-13*r.Lsup()  
          if (q*r).Lsup()>=1.e-13*r.Lsup(): self.__homogeneous_constraint=False  
       if self.debug():  
            if self.__homogeneous_constraint:  
                print "PDE Debug: Constraints are homogeneous."  
            else:  
                print "PDE Debug: Constraints are inhomogeneous."  
   
   
    def hasCoefficient(self,name):  
       """  
       @brief return true if name is the name of a coefficient  
567    
568        @param name       @return: return True is reduced interpolation order is used for the represenation of the solution
569        """       @rtype: L{bool}
570        return self.__coefficient.has_key(name)       """
571         return self.__reduce_solution_order
572    
573     def getFunctionSpaceForEquation(self):     def getFunctionSpaceForEquation(self):
574       """       """
575       @brief return true if the test functions should use reduced order       returns the L{FunctionSpace<escript.FunctionSpace>} used to discretize the equation
576    
577         @return: representation space of equation
578         @rtype: L{FunctionSpace<escript.FunctionSpace>}
579       """       """
580       return self.__row_function_space       if self.reduceEquationOrder():
581             return escript.ReducedSolution(self.getDomain())
582         else:
583             return escript.Solution(self.getDomain())
584    
585     def getFunctionSpaceForSolution(self):     def getFunctionSpaceForSolution(self):
586       """       """
587       @brief return true if the interpolation of the solution should use reduced order       returns the L{FunctionSpace<escript.FunctionSpace>} used to represent the solution
588    
589         @return: representation space of solution
590         @rtype: L{FunctionSpace<escript.FunctionSpace>}
591       """       """
592       return self.__column_function_space       if self.reduceSolutionOrder():
593             return escript.ReducedSolution(self.getDomain())
594         else:
595             return escript.Solution(self.getDomain())
596    
597     # ===== debug ==============================================================     # =============================================================================
598     def setDebugOn(self):     #   solver settings:
599         """     # =============================================================================
600         @brief     def setSolverMethod(self,solver=None,preconditioner=None):
601         """         """
602         self.__debug=not None         sets a new solver
603    
604           @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           @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 setDebugOff(self):     def getSolverPackage(self):
687         """         """
688         @brief         returns the package of the solver
        """  
        self.__debug=None  
689    
690     def debug(self):         @return: the solver package currently being used.
691         """         @rtype: C{int}
        @brief returns true if the PDE is in the debug mode  
692         """         """
693         return self.__debug         return self.__solver_package
   
    #===== Lumping ===========================  
    def setLumpingOn(self):  
       """  
       @brief indicates to use matrix lumping  
       """  
       if not self.isUsingLumping():  
          raise SystemError,"Lumping is not working yet! Talk to the experts"  
          if self.debug() : print "PDE Debug: lumping is set on"  
          self.__rebuildOperator()  
          self.__lumping=True  
   
    def setLumpingOff(self):  
       """  
       @brief 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):  
       """  
       @brief set the matrix lumping flag to flag  
       """  
       if flag:  
          self.setLumpingOn()  
       else:  
          self.setLumpingOff()  
694    
695     def isUsingLumping(self):     def isUsingLumping(self):
696        """        """
697        @brief        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         @brief 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         @brief 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         @brief 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         @brief 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        @brief 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        @brief 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        @brief 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       @brief 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       @brief 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       @brief 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       @brief 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       @brief 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       @brief 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       @brief 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       @brief 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       @brief 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       @brief 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 __makeNewOperator(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         @brief         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 __makeNewRightHandSide(self):     def invalidateSolution(self):
1101         """         """
1102         @brief         indicates the PDE has to be resolved if the solution is requested
1103         """         """
1104         return escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)         self.trace("System will be resolved.")
1105           self.__is_solution_valid=False
1106    
1107     def __makeNewSolution(self):     def isSolutionValid(self):
1108         """         """
1109         @brief         returns True is the solution is still valid.
1110         """         """
1111         return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)         return self.__is_solution_valid
1112    
1113     def __getFreshOperator(self):     def validOperator(self):
1114         """         """
1115         @brief         markes the  operator as valid
1116         """         """
1117         if self.__operator.isEmpty():         self.__is_operator_valid=True
            self.__operator=self.__makeNewOperator()  
            if self.debug() : print "PDE Debug: New operator allocated"  
        else:  
            self.__operator.setValue(0.)  
            if self.debug() : print "PDE Debug: Operator reset to zero"  
        return self.__operator  
1118    
1119     def __getFreshRightHandSide(self):     def invalidateOperator(self):
1120         """         """
1121         @brief         indicates the operator has to be rebuilt next time it is used
1122         """         """
1123         if self.__righthandside.isEmpty():         self.trace("Operator will be rebuilt.")
1124             self.__righthandside=self.__makeNewRightHandSide()         self.invalidateSolution()
1125             if self.debug() : print "PDE Debug: New right hand side allocated"         self.__is_operator_valid=False
        else:  
            print "fix self.__righthandside*=0"  
            self.__righthandside*=0.  
            if self.debug() : print "PDE Debug: Right hand side reset to zero"  
        return  self.__righthandside  
1126    
1127     # ==== rebuild switches =====================================================================     def isOperatorValid(self):
    def __rebuildSolution(self,deep=False):  
1128         """         """
1129         @brief indicates the PDE has to be reolved if the solution is requested         returns True is the operator is still valid.
1130         """         """
1131         if self.__solution_isValid and self.debug() : print "PDE Debug: PDE has to be resolved."         return self.__is_operator_valid
        self.__solution_isValid=False  
        if deep: self.__solution=escript.Data(deep)  
1132    
1133       def validRightHandSide(self):
1134           """
1135           markes the right hand side as valid
1136           """
1137           self.__is_RHS_valid=True
1138    
1139     def __rebuildOperator(self,deep=False):     def invalidateRightHandSide(self):
1140         """         """
1141         @brief indicates the operator has to be rebuilt next time it is used         indicates the right hand side has to be rebuild next time it is used
1142         """         """
1143         if self.__operator_isValid and self.debug() : print "PDE Debug: Operator has to be rebuilt."         if self.isRightHandSideValid(): self.trace("Right hand side has to be rebuilt.")
1144         self.__rebuildSolution(deep)         self.invalidateSolution()
1145         self.__operator_isValid=False         self.__is_RHS_valid=False
        if deep: self.__operator=escript.Operator()  
1146    
1147     def __rebuildRightHandSide(self,deep=False):     def isRightHandSideValid(self):
1148         """         """
1149         @brief indicates the right hand side has to be rebuild next time it is used         returns True is the operator is still valid.
1150         """         """
1151         if self.__righthandside_isValid and self.debug() : print "PDE Debug: Right hand side has to be rebuilt."         return self.__is_RHS_valid
        self.__rebuildSolution(deep)  
        self.__righthandside_isValid=False  
        if not self.__homogeneous_constraint: self.__rebuildOperator()  
        if deep: self.__righthandside=escript.Data()  
1152    
1153     def __rebuildSystem(self,deep=False):     def invalidateSystem(self):
1154         """         """
1155         @brief annonced that all coefficient name has been changed         annonced that everthing has to be rebuild:
1156         """         """
1157         self.__rebuildSolution(deep)         if self.isRightHandSideValid(): self.trace("System has to be rebuilt.")
1158         self.__rebuildOperator(deep)         self.invalidateSolution()
1159         self.__rebuildRightHandSide(deep)         self.invalidateOperator()
1160             self.invalidateRightHandSide()
1161     def __checkMatrixType(self):  
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       @brief reassess the matrix type and, if needed, initiates an operator rebuild       returns the operator of the linear problem
1181    
1182         @return: the operator of the problem
1183       """       """
1184       new_matrix_type=self.getDomain().getSystemMatrixTypeId(self.getSolverMethod(),self.isSymmetric())       return self.getSystem()[0]
      if not new_matrix_type==self.__matrix_type:  
          if self.debug() : print "PDE Debug: Matrix type is now %d."%new_matrix_type  
          self.__matrix_type=new_matrix_type  
          self.__rebuildOperator(deep=True)  
1185    
1186     #============ assembling =======================================================     def getRightHandSide(self):
1187     def __copyConstraint(self,u):       """
1188        """       returns the right hand side of the linear problem
       @brief copies the constrint condition into u  
       """  
       q=self.getCoefficient("q")  
       r=self.getCoefficient("r")  
       if not q.isEmpty():  
           if r.isEmpty():  
              r2=escript.Data(0,u.getShape(),u.getFunctionSpace())  
           else:  
              r2=escript.Data(r,u.getFunctionSpace())  
           u.copyWithMask(r2,escript.Data(q,u.getFunctionSpace()))  
1189    
1190     def __applyConstraint(self,rhs_update=True):       @return: the right hand side of the problem
1191         @rtype: L{Data<escript.Data>}
1192         """
1193         return self.getSystem()[1]
1194    
1195       def createRightHandSide(self):
1196         """         """
1197         @brief applies the constraints  defined by q and r to the system         returns an instance of a new right hand side
1198         """         """
1199         q=self.getCoefficient("q")         self.trace("New right hand side is allocated.")
1200         r=self.getCoefficient("r")         if self.getNumEquations()>1:
1201         if not q.isEmpty() and not self.__operator.isEmpty():             return escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)
1202            # q is the row and column mask to indicate where constraints are set:         else:
1203            row_q=escript.Data(q,self.getFunctionSpaceForEquation())             return escript.Data(0.,(),self.getFunctionSpaceForEquation(),True)
1204            col_q=escript.Data(q,self.getFunctionSpaceForSolution())  
1205            u=self.__makeNewSolution()     def createSolution(self):
1206            if r.isEmpty():         """
1207               r_s=self.__makeNewSolution()         returns an instance of a new solution
1208            else:         """
1209               r_s=escript.Data(r,self.getFunctionSpaceForSolution())         self.trace("New solution is allocated.")
1210            u.copyWithMask(r_s,col_q)         if self.getNumSolutions()>1:
1211            if not self.__righthandside.isEmpty() and rhs_update: self.__righthandside-=self.__operator*u             return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)
1212            self.__operator.nullifyRowsAndCols(row_q,col_q,1.)         else:
1213               return escript.Data(0.,(),self.getFunctionSpaceForSolution(),True)
1214    
1215       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           self.__solution=u
1229           self.validSolution()
1230    
1231       def getCurrentSolution(self):
1232           """
1233           returns the solution in its current state
1234           """
1235           if self.__solution.isEmpty(): self.__solution=self.createSolution()
1236           return self.__solution
1237    
1238       def resetRightHandSide(self):
1239           """
1240           sets the right hand side to zero
1241           """
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 getOperator(self):     def getCurrentRightHandSide(self):
1249         """         """
1250         @brief returns the operator of the PDE         returns the right hand side  in its current state
1251         """         """
1252         if not self.__operator_isValid:         if self.__righthandside.isEmpty(): self.__righthandside=self.createRightHandSide()
1253             # some Constraints are applying for a lumpled stifness matrix:         return self.__righthandside
1254            
1255    
1256       def resetOperator(self):
1257           """
1258           makes sure that the operator is instantiated and returns it initialized by zeros
1259           """
1260           if self.__operator.isEmpty():
1261             if self.isUsingLumping():             if self.isUsingLumping():
1262                if self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():                 self.__operator=self.createSolution()
                        raise TypeError,"Lumped matrix requires same order for equations and unknowns"  
               if not self.getCoefficient("A").isEmpty():  
                        raise Warning,"Lumped matrix does not allow coefficient A"  
               if not self.getCoefficient("B").isEmpty():  
                        raise Warning,"Lumped matrix does not allow coefficient B"  
               if not self.getCoefficient("C").isEmpty():  
                        raise Warning,"Lumped matrix does not allow coefficient C"  
               if not self.getCoefficient("D").isEmpty():  
                        raise Warning,"Lumped matrix does not allow coefficient D"  
               if self.debug() : print "PDE Debug: New lumped operator is built."  
               mat=self.__makeNewOperator(self)  
1263             else:             else:
1264                if self.debug() : print "PDE Debug: New operator is built."                 self.__operator=self.createOperator()
1265                mat=self.__getFreshOperator()         else:
   
            self.getDomain().addPDEToSystem(mat,escript.Data(), \  
                         self.getCoefficient("A"), \  
                         self.getCoefficient("B"), \  
                         self.getCoefficient("C"), \  
                         self.getCoefficient("D"), \  
                         escript.Data(), \  
                         escript.Data(), \  
                         self.getCoefficient("d"), \  
                         escript.Data(),\  
                         self.getCoefficient("d_contact"), \  
                         escript.Data())  
1266             if self.isUsingLumping():             if self.isUsingLumping():
1267                self.__operator=mat*escript.Data(1,(self.getNumSolutions(),),self.getFunctionSpaceOfSolution(),True)                 self.__operator.setToZero()
1268             else:             else:
1269                self.__applyConstraint(rhs_update=False)                 self.__operator.resetValues()
1270             self.__operator_isValid=True             self.trace("Operator reset to zero")
1271    
1272       def getCurrentOperator(self):
1273           """
1274           returns the operator in its current state
1275           """
1276           if self.__operator.isEmpty(): self.resetOperator()
1277         return self.__operator         return self.__operator
1278    
1279     def getRightHandSide(self,ignoreConstraint=False):     def setValue(self,**coefficients):
1280          """
1281          sets new values to coefficients
1282    
1283          @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1284          """
1285          # check if the coefficients are  legal:
1286          for i in coefficients.iterkeys():
1287             if not self.hasCoefficient(i):
1288                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:
1290          if self.__numEquations==None or self.__numSolutions==None:
1291             for i,d in coefficients.iteritems():
1292                if hasattr(d,"shape"):
1293                    s=d.shape
1294                elif hasattr(d,"getShape"):
1295                    s=d.getShape()
1296                else:
1297                    s=numarray.array(d).shape
1298                if s!=None:
1299                    # get number of equations and number of unknowns:
1300                    res=self.__COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(self.getDomain(),s)
1301                    if res==None:
1302                        raise IllegalCoefficientValue,"Illegal shape %s of coefficient %s"%(s,i)
1303                    else:
1304                        if self.__numEquations==None: self.__numEquations=res[0]
1305                        if self.__numSolutions==None: self.__numSolutions=res[1]
1306          if self.__numEquations==None: raise UndefinedPDEError,"unidententified number of equations"
1307          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:
1309          for i,d in coefficients.iteritems():
1310            try:
1311               self.__COEFFICIENTS[i].setValue(self.getDomain(),
1312                                             self.getNumEquations(),self.getNumSolutions(),
1313                                             self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
1314               self.alteredCoefficient(i)
1315            except IllegalCoefficientFunctionSpace,m:
1316                # if the function space is wrong then we try the reduced version:
1317                i_red=i+"_reduced"
1318                if (not i_red in coefficients.keys()) and i_red in self.__COEFFICIENTS.keys():
1319                    try:
1320                        self.__COEFFICIENTS[i_red].setValue(self.getDomain(),
1321                                                          self.getNumEquations(),self.getNumSolutions(),
1322                                                          self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
1323                        self.alteredCoefficient(i_red)
1324                    except IllegalCoefficientValue,m:
1325                        raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
1326                    except IllegalCoefficientFunctionSpace,m:
1327                        raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m))
1328                else:
1329                    raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m))
1330            except IllegalCoefficientValue,m:
1331               raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
1332          self.__altered_coefficients=True
1333    
1334       # ====================================================================================
1335       # methods that are typically overwritten when implementing a particular linear problem
1336       # ====================================================================================
1337       def getRequiredSystemType(self):
1338          """
1339          returns the system type which needs to be used by the current set up.
1340    
1341          @note: Typically this method is overwritten when implementing a particular linear problem.
1342          """
1343          return 0
1344    
1345       def createOperator(self):
1346         """         """
1347         @brief returns the right hand side of the PDE         returns an instance of a new operator
1348    
1349         @param ignoreConstraint         @note: This method is overwritten when implementing a particular linear problem.
1350         """         """
1351         if not self.__righthandside_isValid:         return escript.Operator()
1352             if self.debug() : print "PDE Debug: New right hand side is built."  
1353             self.getDomain().addPDEToRHS(self.__getFreshRightHandSide(), \     def checkSymmetry(self,verbose=True):
1354                           self.getCoefficient("X"), \        """
1355                           self.getCoefficient("Y"),\        test the PDE for symmetry.
                          self.getCoefficient("y"),\  
                          self.getCoefficient("y_contact"))  
            self.__righthandside_isValid=True  
            if ignoreConstraint: self.__copyConstraint(self.__righthandside)  
        return self.__righthandside  
1356    
1357     def getSystem(self):        @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         @brief         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 not self.__operator_isValid and not self.__righthandside_isValid:         return self.getCurrentSolution()
1375            if self.debug() : print "PDE Debug: New PDE is built."  
1376       def getSystem(self):
1377           """
1378           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           return (self.getCurrentOperator(), self.getCurrentRightHandSide())
1386    #=====================
1387    
1388    class LinearPDE(LinearProblem):
1389       """
1390       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    
1393       For a single PDE with a solution with a single component the linear PDE is defined in the following form:
1394    
1395       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    
1397    
1398       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       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       L{ReducedFunction<escript.ReducedFunction>}. It is also allowd to use objects that can be converted into
1403       such L{Data<escript.Data>} objects. M{A} and M{A_reduced} are rank two, M{B}, M{C}, M{X},
1404       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    
1406       The following natural boundary conditions are considered:
1407    
1408       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    
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    
1412    
1413       Constraints for the solution prescribing the value of the solution at certain locations in the domain. They have the form
1414    
1415       M{u=r}  where M{q>0}
1416    
1417       M{r} and M{q} are each scalar where M{q} is the characteristic function defining where the constraint is applied.
1418       The constraints override any other condition set by the PDE or the boundary condition.
1419    
1420       The PDE is symmetrical if
1421    
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       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           returns an instance of a new operator
1571           """
1572           self.trace("New operator is allocated.")
1573           return self.getDomain().newOperator( \
1574                               self.getNumEquations(), \
1575                               self.getFunctionSpaceForEquation(), \
1576                               self.getNumSolutions(), \
1577                               self.getFunctionSpaceForSolution(), \
1578                               self.getSystemType())
1579    
1580       def getSolution(self,**options):
1581           """
1582           returns the solution of the PDE
1583    
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.getRightHandSide(ignoreConstraint=True)               self.setSolution(f*1/mat)
               self.getOperator()  
1601            else:            else:
1602                self.getDomain().addPDEToSystem(self.__getFreshOperator(),self.__getFreshRightHandSide(), \               options[self.TOLERANCE_KEY]=self.getTolerance()
1603                 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):
1613           """
1614           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.isOperatorValid() or not self.isRightHandSideValid():
1620              if self.isUsingLumping():
1621                  if not self.isOperatorValid():
1622                     if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():
1623                          raise TypeError,"Lumped matrix requires same order for equations and unknowns"
1624                     if not self.getCoefficient("A").isEmpty():
1625                          raise ValueError,"coefficient A in lumped matrix may not be present."
1626                     if not self.getCoefficient("B").isEmpty():
1627                          raise ValueError,"coefficient B in lumped matrix may not be present."
1628                     if not self.getCoefficient("C").isEmpty():
1629                          raise ValueError,"coefficient C in lumped matrix may not be present."
1630                     if not self.getCoefficient("d_contact").isEmpty():
1631                          raise ValueError,"coefficient d_contact in lumped matrix may not be present."
1632                     if not self.getCoefficient("A_reduced").isEmpty():
1633                          raise ValueError,"coefficient A_reduced in lumped matrix may not be present."
1634                     if not self.getCoefficient("B_reduced").isEmpty():
1635                          raise ValueError,"coefficient B_reduced in lumped matrix may not be present."
1636                     if not self.getCoefficient("C_reduced").isEmpty():
1637                          raise ValueError,"coefficient C_reduced in lumped matrix may not be present."
1638                     if not self.getCoefficient("d_contact_reduced").isEmpty():
1639                          raise ValueError,"coefficient d_contact_reduced in lumped matrix may not be present."
1640                     D=self.getCoefficient("D")
1641                     d=self.getCoefficient("d")
1642                     D_reduced=self.getCoefficient("D_reduced")
1643                     d_reduced=self.getCoefficient("d_reduced")
1644                     if not D.isEmpty():
1645                         if self.getNumSolutions()>1:
1646                            D_times_e=util.matrix_mult(D,numarray.ones((self.getNumSolutions(),)))
1647                         else:
1648                            D_times_e=D
1649                     else:
1650                        D_times_e=escript.Data()
1651                     if not d.isEmpty():
1652                         if self.getNumSolutions()>1:
1653                            d_times_e=util.matrix_mult(d,numarray.ones((self.getNumSolutions(),)))
1654                         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:
1709                 if not self.isOperatorValid() and not self.isRightHandSideValid():
1710                     self.resetRightHandSide()
1711                     righthandside=self.getCurrentRightHandSide()
1712                     self.resetOperator()
1713                     operator=self.getCurrentOperator()
1714                     self.getDomain().addPDEToSystem(operator,righthandside, \
1715                                   self.getCoefficient("A"), \
1716                                   self.getCoefficient("B"), \
1717                                   self.getCoefficient("C"), \
1718                                   self.getCoefficient("D"), \
1719                                   self.getCoefficient("X"), \
1720                                   self.getCoefficient("Y"), \
1721                                   self.getCoefficient("d"), \
1722                                   self.getCoefficient("y"), \
1723                                   self.getCoefficient("d_contact"), \
1724                                   self.getCoefficient("y_contact"))
1725                     self.getDomain().addPDEToSystem(operator,righthandside, \
1726                                   self.getCoefficient("A_reduced"), \
1727                                   self.getCoefficient("B_reduced"), \
1728                                   self.getCoefficient("C_reduced"), \
1729                                   self.getCoefficient("D_reduced"), \
1730                                   self.getCoefficient("X_reduced"), \
1731                                   self.getCoefficient("Y_reduced"), \
1732                                   self.getCoefficient("d_reduced"), \
1733                                   self.getCoefficient("y_reduced"), \
1734                                   self.getCoefficient("d_contact_reduced"), \
1735                                   self.getCoefficient("y_contact_reduced"))
1736                     self.insertConstraint(rhs_only=False)
1737                     self.trace("New system has been built.")
1738                     self.validOperator()
1739                     self.validRightHandSide()
1740                 elif not self.isRightHandSideValid():
1741                     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"), \                              self.getCoefficient("A"), \
1761                              self.getCoefficient("B"), \                              self.getCoefficient("B"), \
1762                              self.getCoefficient("C"), \                              self.getCoefficient("C"), \
1763                              self.getCoefficient("D"), \                              self.getCoefficient("D"), \
1764                              self.getCoefficient("X"), \                              escript.Data(), \
1765                              self.getCoefficient("Y"), \                              escript.Data(), \
1766                              self.getCoefficient("d"), \                              self.getCoefficient("d"), \
1767                              self.getCoefficient("y"), \                              escript.Data(),\
1768                              self.getCoefficient("d_contact"), \                              self.getCoefficient("d_contact"), \
1769                              self.getCoefficient("y_contact"))                              escript.Data())
1770            self.__operator_isValid=True                   self.getDomain().addPDEToSystem(operator,escript.Data(), \
1771            self.__righthandside_isValid=True                              self.getCoefficient("A_reduced"), \
1772            self.__applyConstraint()                              self.getCoefficient("B_reduced"), \
1773            self.__copyConstraint(self.__righthandside)                              self.getCoefficient("C_reduced"), \
1774         elif not self.__operator_isValid:                              self.getCoefficient("D_reduced"), \
1775            self.getOperator()                              escript.Data(), \
1776         elif not self.__righthandside_isValid:                              escript.Data(), \
1777            self.getRightHandSide()                              self.getCoefficient("d_reduced"), \
1778         return (self.__operator,self.__righthandside)                              escript.Data(),\
1779                                self.getCoefficient("d_contact_reduced"), \
1780     def solve(self,**options):                              escript.Data())
1781        """                   self.insertConstraint(rhs_only=False)
1782        @brief solve the PDE                   self.trace("New operator has been built.")
1783                     self.validOperator()
1784        @param options         return (self.getCurrentOperator(), self.getCurrentRightHandSide())
       """  
       mat,f=self.getSystem()  
       if self.isUsingLumping():  
          out=f/mat  
          self.__copyConstraint(out)  
       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  
1785    
1786     def getSolution(self,**options):     def insertConstraint(self, rhs_only=False):
1787         """        """
1788         @brief returns the solution of the PDE        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          q=self.getCoefficient("q")
1794          r=self.getCoefficient("r")
1795          righthandside=self.getCurrentRightHandSide()
1796          operator=self.getCurrentOperator()
1797    
1798         @param options        if not q.isEmpty():
1799         """           if r.isEmpty():
1800         if not self.__solution_isValid:              r_s=self.createSolution()
1801             if self.debug() : print "PDE Debug: PDE is resolved."           else:
1802             self.__solution=self.solve(**options)              r_s=r
1803             self.__solution_isValid=True           if not rhs_only and not operator.isEmpty():
1804         return self.__solution               if self.isUsingLumping():
1805     #============ some serivice functions  =====================================================                   operator.copyWithMask(escript.Data(1.,q.getShape(),q.getFunctionSpace()),q)
1806     def getDomain(self):               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       def setValue(self,**coefficients):
1816          """
1817          sets new values to coefficients
1818    
1819          @param coefficients: new values assigned to coefficients
1820          @keyword A: value for coefficient A.
1821          @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       @brief returns the domain of the PDE       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       def getFlux(self,u=None):
1893       """       """
1894       return self.__domain       returns the flux M{J} for a given M{u}
1895    
1896     def getDim(self):       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    
1898         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       @brief returns the spatial dimension of the PDE       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       M{n[j]*grad(u)[j] = 0 }
1925    
1926       and constraints:
1927    
1928       M{u=0} where M{q>0}
1929    
1930       """
1931    
1932       def __init__(self,domain,debug=False):
1933       """       """
1934       return self.getDomain().getDim()       initializes a new Poisson equation
1935    
1936         @param domain: domain of the PDE
1937         @type domain: L{Domain<escript.Domain>}
1938         @param debug: if True debug informations are printed.
1939    
    def getNumEquations(self):  
1940       """       """
1941       @brief returns the number of equations       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       def setValue(self,**coefficients):
1948       """       """
1949       if self.__numEquations>0:       sets new values to coefficients
1950           return self.__numEquations  
1951         @param coefficients: new values assigned to coefficients
1952         @keyword f: value for right hand side M{f}
1953         @type f: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1954         @keyword q: mask for location of constraints
1955         @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                   depending of reduced order is used for the representation of the equation.
1957         @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1958         """
1959         super(Poisson, self).setValue(**coefficients)
1960    
1961    
1962       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:       else:
1979           raise ValueError,"Number of equations is undefined. Please specify argument numEquations."          return super(Poisson, self).getCoefficient(name)
1980    
1981     def getNumSolutions(self):  class Helmholtz(LinearPDE):
1982       """     """
1983       @brief returns the number of unknowns     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       M{k*n[j]*grad(u)[j] = g- S{alpha}u }
1990    
1991       and constraints:
1992    
1993       M{u=r} where M{q>0}
1994    
1995       """
1996    
1997       def __init__(self,domain,debug=False):
1998       """       """
1999       if self.__numSolutions>0:       initializes a new Poisson equation
2000          return self.__numSolutions  
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 getCoefficient(self,name):
2043         """
2044         return the value of the coefficient name of the general PDE
2045    
2046         @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:       else:
2067          raise ValueError,"Number of solution is undefined. Please specify argument numSolutions."          return super(Helmholtz, self).getCoefficient(name)
2068    
2069    class LameEquation(LinearPDE):
2070       """
2071       Class to define a Lame equation problem:
2072    
2073     def checkSymmetry(self):     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] }
       """  
       @brief returns if the Operator is symmetric. This is a very expensive operation!!! The symmetry flag is not altered.  
       """  
       raise SystemError,"checkSymmetry is not implemented yet"  
2074    
2075        return None     with natural boundary conditons:
2076    
2077     def getFlux(self,u):     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] }
        """  
        @brief returns the flux J_ij for a given u  
2078    
2079              J_ij=A_{ijkl}u_{k,l}+B_{ijk}u_k-X_{ij}     and constraints:
2080    
2081         @param u argument of the operator     M{u[i]=r[i]} where M{q[i]>0}
2082    
2083         """     """
        raise SystemError,"getFlux is not implemented yet"  
        return None  
2084    
2085     def applyOperator(self,u):     def __init__(self,domain,debug=False):
2086         """        super(LameEquation, self).__init__(domain,\
2087         @brief applies the operator of the PDE to a given solution u in weak from                                           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         return the value of the coefficient name of the general PDE
2123    
2124         @param name: name of the coefficient requested.
2125         @type name: C{string}
2126         @return: the value of the coefficient  name
2127         @rtype: L{Data<escript.Data>}
2128         @raise IllegalCoefficient: invalid coefficient name
2129         """
2130         if name == "A" :
2131             out =self.createCoefficient("A")
2132             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")
2144         else:
2145            return super(LameEquation, self).getCoefficient(name)
2146    
2147         @param u argument of the operator  def LinearSinglePDE(domain,debug=False):
2148       """
2149       defines a single linear PDEs
2150    
2151         """     @param domain: domain of the PDE
2152         return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())     @type domain: L{Domain<escript.Domain>}
2153                                                                                                                                                                 @param debug: if True debug informations are printed.
2154     def getResidual(self,u):     @rtype: L{LinearPDE}
2155         """     """
2156         @brief return the residual of u in the weak from     return LinearPDE(domain,numEquations=1,numSolutions=1,debug=debug)
2157    
2158         @param u  def LinearPDESystem(domain,debug=False):
2159         """     """
2160         return self.applyOperator(u)-self.getRightHandSide()     defines a system of linear PDEs
2161    
2162  class Poisson(LinearPDE):     @param domain: domain of the PDE
2163       @type domain: L{Domain<escript.Domain>}
2164       @param debug: if True debug informations are printed.
2165       @rtype: L{LinearPDE}
2166     """     """
2167     @brief Class to define a Poisson equstion problem:     return LinearPDE(domain,numEquations=domain.getDim(),numSolutions=domain.getDim(),debug=debug)
2168                                                                                                                                                                
2169     class to define a linear PDE of the form  class TransportPDE(LinearProblem):
2170                                                                                                                                                                   """
2171          -u_{,jj} = f     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       with boundary conditons:  
2174                                                                                                                                                                   For a single equation with a solution with a single component the transport problem is defined in the following form:
2175          n_j*u_{,j} = 0  
2176                                                                                                                                                                   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      and constraints:  
2178                                                                                                                                                                
2179           u=0 where q>0     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                                                                                                                                                                   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       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     def __init__(self,domain=None,f=escript.Data(),q=escript.Data()):     L{ReducedFunction<escript.ReducedFunction>}. It is also allowed to use objects that can be converted into
2184         LinearPDE.__init__(self,domain=identifyDomain(domain,{"f":f, "q":q}))     such L{Data<escript.Data>} objects. M{M} and M{M_reduced} are scalar, M{A} and M{A_reduced} are rank two,
2185         self._setValue(A=numarray.identity(self.getDomain().getDim()))     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         self.setSymmetryOn()  
2187         self.setValue(f,q)     The following natural boundary conditions are considered:
2188    
2189     def setValue(self,f=escript.Data(),q=escript.Data()):     M{n[j]*((A[i,j]+A_reduced[i,j])*grad(u)[l]+(B+B_reduced)[j]*u+X[j]+X_reduced[j])+(d+d_reduced)*u+y+y_reduced=(m+m_reduced)*u_t}
2190         self._setValue(Y=f,q=q)  
2191       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                                                                                                                                                              
2193  # $Log$     Constraints for the solution prescribing the value of the solution at certain locations in the domain. They have the form
2194  # Revision 1.3  2004/12/17 07:43:10  jgs  
2195  # *** empty log message ***     M{u_t=r}  where M{q>0}
2196  #  
2197  # Revision 1.1.2.3  2004/12/16 00:12:34  gross     M{r} and M{q} are each scalar where M{q} is the characteristic function defining where the constraint is applied.
2198  # __init__ of LinearPDE does not accept any coefficients anymore     The constraints override any other condition set by the transport problem or the boundary condition.
2199  #  
2200  # Revision 1.1.2.2  2004/12/14 03:55:01  jgs     The transport problem is symmetrical if
2201  # *** empty log message ***  
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.1  2004/12/12 22:53:47  gross  
2204  # linearPDE has been renamed LinearPDE     For a system and a solution with several components the transport problem has the form
2205  #  
2206  # Revision 1.1.1.1.2.7  2004/12/07 10:13:08  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  # GMRES added  
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.1.1.2.6  2004/12/07 03:19:50  gross  
2210  # options for GMRES and PRES20 added     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.1.1.2.5  2004/12/01 06:25:15  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  # some small changes  
2214  #     Constraints take the form
2215  # Revision 1.1.1.1.2.4  2004/11/24 01:50:21  gross  
2216  # Finley solves 4M unknowns now     M{u[i]_t=r[i]}  where  M{q[i]>0}
2217  #  
2218  # Revision 1.1.1.1.2.3  2004/11/15 06:05:26  gross     M{r} and M{q} are each rank one. Notice that at some locations not necessarily all components must have a constraint.
2219  # poisson solver added  
2220  #     The transport problem is symmetrical if
2221  # Revision 1.1.1.1.2.2  2004/11/12 06:58:15  gross  
2222  # 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{M[i,k]=M[i,k]}
2223  #          - M{M_reduced[i,k]=M_reduced[i,k]}
2224  # Revision 1.1.1.1.2.1  2004/10/28 22:59:22  gross          - M{A[i,j,k,l]=A[k,l,i,j]}
2225  # finley's RecTest.py is running now: problem in SystemMatrixAdapater fixed          - 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  2004/10/26 06:53:56  jgs          - M{B_reduced[i,j,k]=C_reduced[k,i,j]}
2228  # initial import of project esys2          - M{D[i,k]=D[i,k]}
2229  #          - M{D_reduced[i,k]=D_reduced[i,k]}
2230  # Revision 1.3.2.3  2004/10/26 06:43:48  jgs          - M{m[i,k]=m[k,i]}
2231  # committing Lutz's and Paul's changes to brach jgs          - M{m_reduced[i,k]=m_reduced[k,i]}
2232  #          - M{d[i,k]=d[k,i]}
2233  # Revision 1.3.4.1  2004/10/20 05:32:51  cochrane          - M{d_reduced[i,k]=d_reduced[k,i]}
2234  # Added incomplete Doxygen comments to files, or merely put the docstrings that already exist into Doxygen form.  
2235  #     L{TransportPDE} also supports solution discontinuities over a contact region in the domain. To specify the conditions across the
2236  # Revision 1.3  2004/09/23 00:53:23  jgs     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  # minor fixes     defined as
2238  #  
2239  # Revision 1.1  2004/08/28 12:58:06  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  # SimpleSolve is not running yet: problem with == of functionsspace  
2241  #     For the case of single solution component and single PDE M{J} is defined
2242  #  
2243       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       In the context of discontinuities M{n} denotes the normal on the discontinuity pointing from side 0 towards side 1
2246       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       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    
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       of the solution at side 1 and at side 0, denotes the jump of M{u} across discontinuity along the normal calculated by
2252       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       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       In case of a single PDE and a single component solution the contact condition takes the form
2256    
2257       M{n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)}
2258    
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    
2261       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.104  
changed lines
  Added in v.1841

  ViewVC Help
Powered by ViewVC 1.1.26