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

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

  ViewVC Help
Powered by ViewVC 1.1.26