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

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

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

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

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

  ViewVC Help
Powered by ViewVC 1.1.26