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

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

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

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

Legend:
Removed from v.142  
changed lines
  Added in v.2100

  ViewVC Help
Powered by ViewVC 1.1.26