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

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

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

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

Legend:
Removed from v.104  
changed lines
  Added in v.2197

  ViewVC Help
Powered by ViewVC 1.1.26