/[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

revision 614 by elspeth, Wed Mar 22 01:37:07 2006 UTC revision 2323 by gross, Thu Mar 19 04:23:32 2009 UTC
# Line 1  Line 1 
1  # $Id$  
2    ########################################################
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  The module provides an interface to define and solve linear partial  The module provides an interface to define and solve linear partial
24  differential equations (PDEs) within L{escript}. L{linearPDEs} does not provide any  differential equations (PDEs) and Transport problems within L{escript}.
25  solver capabilities in itself but hands the PDE over to  L{linearPDEs} does not provide any solver capabilities in itself but hands the
26  the PDE solver library defined through the L{Domain<escript.Domain>} of the PDE.  PDE over to the PDE solver library defined through the L{Domain<escript.Domain>}
27  The general interface is provided through the L{LinearPDE} class. The  of the PDE. The general interface is provided through the L{LinearPDE} class.
28  L{AdvectivePDE} which is derived from the L{LinearPDE} class  L{TransportProblem} provides an interface to initial value problems dominated
29  provides an interface to PDE dominated by its advective terms. The L{Poisson},  by its advective terms.
 L{Helmholtz}, L{LameEquation}, L{AdvectionDiffusion}  
 classs which are also derived form the L{LinearPDE} class should be used  
 to define of solve these sepecial PDEs.  
30    
31  @var __author__: name of author  @var __author__: name of author
32    @var __copyright__: copyrights
33  @var __license__: licence agreement  @var __license__: licence agreement
34  @var __url__: url entry point on documentation  @var __url__: url entry point on documentation
35  @var __version__: version  @var __version__: version
36  @var __date__: date of the version  @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  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
 __copyright__="""  Copyright (c) 2006 by ACcESS MNRF  
                     http://www.access.edu.au  
                 Primary Business: Queensland, Australia"""  
 __license__="""Licensed under the Open Software License version 3.0  
              http://www.opensource.org/licenses/osl-3.0.php"""  
 __url__="http://www.iservo.edu.au/esys/escript"  
 __version__="$Revision$"  
 __date__="$Date$"  
45    
46    
47  class IllegalCoefficient(ValueError):  class IllegalCoefficient(ValueError):
48     """     """
49     raised if an illegal coefficient of the general ar particular PDE is requested.     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):  class IllegalCoefficientValue(ValueError):
55     """     """
56     raised if an incorrect value for a coefficient is used.     Exception that is raised if an incorrect value for a coefficient is used.
57       """
58       pass
59    
60    class IllegalCoefficientFunctionSpace(ValueError):
61       """
62       Exception that is raised if an incorrect function space for a coefficient
63       is used.
64     """     """
65    
66  class UndefinedPDEError(ValueError):  class UndefinedPDEError(ValueError):
67     """     """
68     raised if a PDE is not fully defined yet.     Exception that is raised if a PDE is not fully defined yet.
69     """     """
70       pass
71    
72  class PDECoefficient(object):  class PDECoef(object):
73      """      """
74      A class for describing a PDE coefficient      A class for describing a PDE coefficient.
75    
76      @cvar INTERIOR: indicator that coefficient is defined on the interior of the PDE domain      @cvar INTERIOR: indicator that coefficient is defined on the interior of
77      @cvar BOUNDARY: indicator that coefficient is defined on the boundary of the PDE domain                      the PDE domain
78      @cvar CONTACT: indicator that coefficient is defined on the contact region within the PDE domain      @cvar BOUNDARY: indicator that coefficient is defined on the boundary of
79      @cvar SOLUTION: indicator that coefficient is defined trough a solution of the PDE                      the PDE domain
80      @cvar REDUCED: indicator that coefficient is defined trough a reduced solution of the PDE      @cvar CONTACT: indicator that coefficient is defined on the contact region
81      @cvar BY_EQUATION: indicator that the dimension of the coefficient shape is defined by the number PDE equations                     within the PDE domain
82      @cvar BY_SOLUTION: indicator that the dimension of the coefficient shape is defined by the number PDE solutions      @cvar INTERIOR_REDUCED: indicator that coefficient is defined on the
83      @cvar BY_DIM: indicator that the dimension of the coefficient shape is defined by the spatial dimension                              interior of the PDE domain using a reduced
84      @cvar OPERATOR: indicator that the the coefficient alters the operator of the PDE                              integration order
85      @cvar RIGHTHANDSIDE: indicator that the the coefficient alters the right hand side of the PDE      @cvar BOUNDARY_REDUCED: indicator that coefficient is defined on the
86      @cvar BOTH: indicator that the the coefficient alters the operator as well as the right hand side of the PDE                              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      """      """
109      INTERIOR=0      INTERIOR=0
# Line 76  class PDECoefficient(object): Line 117  class PDECoefficient(object):
117      OPERATOR=10      OPERATOR=10
118      RIGHTHANDSIDE=11      RIGHTHANDSIDE=11
119      BOTH=12      BOTH=12
120        INTERIOR_REDUCED=13
121        BOUNDARY_REDUCED=14
122        CONTACT_REDUCED=15
123    
124      def __init__(self,where,pattern,altering):      def __init__(self, where, pattern, altering):
125         """         """
126         Initialise a PDE Coefficient type         Initialises a PDE coefficient type.
127    
128         @param where: describes where the coefficient lives         @param where: describes where the coefficient lives
129         @type where: one of L{INTERIOR}, L{BOUNDARY}, L{CONTACT}, L{SOLUTION}, L{REDUCED}         @type where: one of L{INTERIOR}, L{BOUNDARY}, L{CONTACT}, L{SOLUTION},
130         @param pattern: describes the shape of the coefficient and how the shape is build for a given                      L{REDUCED}, L{INTERIOR_REDUCED}, L{BOUNDARY_REDUCED},
131                spatial dimension and numbers of equation and solution in then PDE. For instance,                      L{CONTACT_REDUCED}
132                (L{BY_EQUATION},L{BY_SOLUTION},L{BY_DIM}) descrbes a rank 3 coefficient which         @param pattern: describes the shape of the coefficient and how the shape
133                is instanciated as shape (3,2,2) in case of a three equations and two solution components                         is built for a given spatial dimension and numbers of
134                on a 2-dimensional domain. In the case of single equation and a single solution component                         equations and solutions in then PDE. For instance,
135                the shape compoments marked by L{BY_EQUATION} or L{BY_SOLUTION} are dropped. In this case                         (L{BY_EQUATION},L{BY_SOLUTION},L{BY_DIM}) describes a
136                the example would be read as (2,).                         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}         @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 coefficiennt is altered         @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}         @type altering: one of L{OPERATOR}, L{RIGHTHANDSIDE}, L{BOTH}
   
146         """         """
147         super(PDECoefficient, self).__init__()         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
# Line 103  class PDECoefficient(object): Line 152  class PDECoefficient(object):
152    
153      def resetValue(self):      def resetValue(self):
154         """         """
155         resets coefficient value to default         Resets the coefficient value to the default.
156         """         """
157         self.value=escript.Data()         self.value=escript.Data()
158    
159      def getFunctionSpace(self,domain,reducedEquationOrder=False,reducedSolutionOrder=False):      def getFunctionSpace(self,domain,reducedEquationOrder=False,reducedSolutionOrder=False):
160         """         """
161         defines the L{FunctionSpace<escript.FunctionSpace>} of the coefficient         Returns the L{FunctionSpace<escript.FunctionSpace>} of the coefficient.
162    
163         @param domain: domain on which the PDE uses the coefficient         @param domain: domain on which the PDE uses the coefficient
164         @type domain: L{Domain<escript.Domain>}         @type domain: L{Domain<escript.Domain>}
165         @param reducedEquationOrder: True to indicate that reduced order is used to represent the equation         @param reducedEquationOrder: True to indicate that reduced order is used
166         @type domain: C{bool}                                      to represent the equation
167         @param reducedSolutionOrder: True to indicate that reduced order is used to represent the solution         @type reducedEquationOrder: C{bool}
168         @type domain: C{bool}         @param reducedSolutionOrder: True to indicate that reduced order is used
169         @return:  L{FunctionSpace<escript.FunctionSpace>} of the coefficient                                      to represent the solution
170         @rtype:  L{FunctionSpace<escript.FunctionSpace>}         @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:         if self.what==self.INTERIOR:
175              return escript.Function(domain)              return escript.Function(domain)
176           elif self.what==self.INTERIOR_REDUCED:
177                return escript.ReducedFunction(domain)
178         elif self.what==self.BOUNDARY:         elif self.what==self.BOUNDARY:
179              return escript.FunctionOnBoundary(domain)              return escript.FunctionOnBoundary(domain)
180           elif self.what==self.BOUNDARY_REDUCED:
181                return escript.ReducedFunctionOnBoundary(domain)
182         elif self.what==self.CONTACT:         elif self.what==self.CONTACT:
183              return escript.FunctionOnContactZero(domain)              return escript.FunctionOnContactZero(domain)
184           elif self.what==self.CONTACT_REDUCED:
185                return escript.ReducedFunctionOnContactZero(domain)
186         elif self.what==self.SOLUTION:         elif self.what==self.SOLUTION:
187              if reducedEquationOrder and reducedSolutionOrder:              if reducedEquationOrder and reducedSolutionOrder:
188                  return escript.ReducedSolution(domain)                  return escript.ReducedSolution(domain)
# Line 136  class PDECoefficient(object): Line 193  class PDECoefficient(object):
193    
194      def getValue(self):      def getValue(self):
195         """         """
196         returns the value of the coefficient         Returns the value of the coefficient.
197    
198         @return:  value of the coefficient         @return: value of the coefficient
199         @rtype:  L{Data<escript.Data>}         @rtype: L{Data<escript.Data>}
200         """         """
201         return self.value         return self.value
202    
203      def setValue(self,domain,numEquations=1,numSolutions=1,reducedEquationOrder=False,reducedSolutionOrder=False,newValue=None):      def setValue(self,domain,numEquations=1,numSolutions=1,reducedEquationOrder=False,reducedSolutionOrder=False,newValue=None):
204         """         """
205         set the value of the coefficient to a new value         Sets the value of the coefficient to a new value.
206    
207         @param domain: domain on which the PDE uses the coefficient         @param domain: domain on which the PDE uses the coefficient
208         @type domain: L{Domain<escript.Domain>}         @type domain: L{Domain<escript.Domain>}
# Line 153  class PDECoefficient(object): Line 210  class PDECoefficient(object):
210         @type numEquations: C{int}         @type numEquations: C{int}
211         @param numSolutions: number of components of the PDE solution         @param numSolutions: number of components of the PDE solution
212         @type numSolutions: C{int}         @type numSolutions: C{int}
213         @param reducedEquationOrder: True to indicate that reduced order is used to represent the equation         @param reducedEquationOrder: True to indicate that reduced order is used
214         @type domain: C{bool}                                      to represent the equation
215         @param reducedSolutionOrder: True to indicate that reduced order is used to represent the solution         @type reducedEquationOrder: C{bool}
216         @type domain: C{bool}         @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         @param newValue: number of components of the PDE solution
220         @type newValue: any object that can be converted into a L{Data<escript.Data>} object with the appropriate shape and L{FunctionSpace<escript.FunctionSpace>}         @type newValue: any object that can be converted into a
221         @raise IllegalCoefficientValue: if the shape of the assigned value does not match the shape of the coefficient                         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:         if newValue==None:
229             newValue=escript.Data()             newValue=escript.Data()
230         elif isinstance(newValue,escript.Data):         elif isinstance(newValue,escript.Data):
231             if not newValue.isEmpty():             if not newValue.isEmpty():
232                try:                if not newValue.getFunctionSpace() == self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder):
233                   newValue=escript.Data(newValue,self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder))                  try:
234                except:                    newValue=escript.Data(newValue,self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder))
235                   raise IllegalCoefficientValue,"Unable to interpolate coefficient to function space %s"%self.getFunctionSpace(domain)                  except:
236                      raise IllegalCoefficientFunctionSpace,"Unable to interpolate coefficient to function space %s"%self.getFunctionSpace(domain)
237         else:         else:
238             newValue=escript.Data(newValue,self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder))             newValue=escript.Data(newValue,self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder))
239         if not newValue.isEmpty():         if not newValue.isEmpty():
# Line 178  class PDECoefficient(object): Line 243  class PDECoefficient(object):
243    
244      def isAlteringOperator(self):      def isAlteringOperator(self):
245          """          """
246          checks if the coefficient alters the operator of the PDE          Checks if the coefficient alters the operator of the PDE.
247    
248          @return:  True if the operator of the PDE is changed when the coefficient is changed          @return: True if the operator of the PDE is changed when the
249          @rtype:  C{bool}                   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 190  class PDECoefficient(object): Line 256  class PDECoefficient(object):
256    
257      def isAlteringRightHandSide(self):      def isAlteringRightHandSide(self):
258          """          """
259          checks if the coefficeint alters the right hand side of the PDE          Checks if the coefficient alters the right hand side of the PDE.
260    
261      @rtype:  C{bool}          @rtype: C{bool}
262          @return:  True if the right hand side of the PDE is changed when the coefficient is changed          @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:
# Line 202  class PDECoefficient(object): Line 269  class PDECoefficient(object):
269    
270      def estimateNumEquationsAndNumSolutions(self,domain,shape=()):      def estimateNumEquationsAndNumSolutions(self,domain,shape=()):
271         """         """
272         tries to estimate the number of equations and number of solutions if the coefficient has the given shape         Tries to estimate the number of equations and number of solutions if
273           the coefficient has the given shape.
274    
275         @param domain: domain on which the PDE uses the coefficient         @param domain: domain on which the PDE uses the coefficient
276         @type domain: L{Domain<escript.Domain>}         @type domain: L{Domain<escript.Domain>}
277         @param shape: suggested shape of the coefficient         @param shape: suggested shape of the coefficient
278         @type shape: C{tuple} of C{int} values         @type shape: C{tuple} of C{int} values
279         @return: the number of equations and number of solutions of the PDE is the coefficient has shape s.         @return: the number of equations and number of solutions of the PDE if
280                   If no appropriate numbers could be identified, C{None} is returned                  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}         @rtype: C{tuple} of two C{int} values or C{None}
283         """         """
284         dim=domain.getDim()         dim=domain.getDim()
# Line 245  class PDECoefficient(object): Line 314  class PDECoefficient(object):
314               else:               else:
315                   if s==shape: return (None,u)                   if s==shape: return (None,u)
316         return None         return None
317    
318      def definesNumSolutions(self):      def definesNumSolutions(self):
319         """         """
320         checks if the coefficient allows to estimate the number of solution components         Checks if the coefficient allows to estimate the number of solution
321           components.
322    
323         @return: True if the coefficient allows an estimate of the number of solution components         @return: True if the coefficient allows an estimate of the number of
324                    solution components, False otherwise
325         @rtype: C{bool}         @rtype: C{bool}
326         """         """
327         for i in self.pattern:         for i in self.pattern:
# Line 258  class PDECoefficient(object): Line 330  class PDECoefficient(object):
330    
331      def definesNumEquation(self):      def definesNumEquation(self):
332         """         """
333         checks if the coefficient allows to estimate the number of equations         Checks if the coefficient allows to estimate the number of equations.
334    
335         @return: True if the coefficient allows an estimate of the number of equations         @return: True if the coefficient allows an estimate of the number of
336                    equations, False otherwise
337         @rtype: C{bool}         @rtype: C{bool}
338         """         """
339         for i in self.pattern:         for i in self.pattern:
# Line 269  class PDECoefficient(object): Line 342  class PDECoefficient(object):
342    
343      def __CompTuple2(self,t1,t2):      def __CompTuple2(self,t1,t2):
344        """        """
345        Compare two tuples of possible number of equations and number of solutions        Compares two tuples of possible number of equations and number of
346          solutions.
       @param t1: The first tuple  
       @param t2: The second tuple  
347    
348          @param t1: the first tuple
349          @param t2: the second tuple
350          @return: 0, 1, or -1
351        """        """
352    
353        dif=t1[0]+t1[1]-(t2[0]+t2[1])        dif=t1[0]+t1[1]-(t2[0]+t2[1])
# Line 283  class PDECoefficient(object): Line 357  class PDECoefficient(object):
357    
358      def getShape(self,domain,numEquations=1,numSolutions=1):      def getShape(self,domain,numEquations=1,numSolutions=1):
359         """         """
360         builds the required shape of the coefficient         Builds the required shape of the coefficient.
361    
362         @param domain: domain on which the PDE uses the coefficient         @param domain: domain on which the PDE uses the coefficient
363         @type domain: L{Domain<escript.Domain>}         @type domain: L{Domain<escript.Domain>}
# Line 305  class PDECoefficient(object): Line 379  class PDECoefficient(object):
379                  s=s+(dim,)                  s=s+(dim,)
380         return s         return s
381    
382  class LinearPDE(object):  class LinearProblem(object):
383     """     """
384     This class is used to define a general linear, steady, second order PDE     This is the base class to define a general linear PDE-type problem for
385     for an unknown function M{u} on a given domain defined through a L{Domain<escript.Domain>} object.     for an unknown function M{u} on a given domain defined through a
386       L{Domain<escript.Domain>} object. The problem can be given as a single
387     For a single PDE with a solution with a single component the linear PDE is defined in the following form:     equation or as a system of equations.
   
    M{-grad(A[j,l]*grad(u)[l]+B[j]u)[j]+C[l]*grad(u)[l]+D*u =-grad(X)[j,j]+Y}  
   
    where M{grad(F)} denotes the spatial derivative of M{F}. Einstein's summation convention,  
    ie. summation over indexes appearing twice in a term of a sum is performed, is used.  
    The coefficients M{A}, M{B}, M{C}, M{D}, M{X} and M{Y} have to be specified through L{Data<escript.Data>} objects in the  
    L{Function<escript.Function>} on the PDE or objects that can be converted into such L{Data<escript.Data>} objects.  
    M{A} is a rank two, M{B}, M{C} and M{X} are rank one and M{D} and M{Y} are scalar.  
   
    The following natural boundary conditions are considered:  
   
    M{n[j]*(A[i,j]*grad(u)[l]+B[j]*u)+d*u=n[j]*X[j]+y}  
   
    where M{n} is the outer normal field calculated by L{getNormal<escript.FunctionSpace.getNormal>} of L{FunctionOnBoundary<escript.FunctionOnBoundary>}.  
    Notice that the coefficients M{A}, M{B} and M{X} are defined in the PDE. The coefficients M{d} and M{y} are  
    each a scalar in the L{FunctionOnBoundary<escript.FunctionOnBoundary>}.  
   
   
    Constraints for the solution prescribing the value of the solution at certain locations in the domain. They have the form  
   
    M{u=r}  where M{q>0}  
   
    M{r} and M{q} are each scalar where M{q} is the characteristic function defining where the constraint is applied.  
    The constraints override any other condition set by the PDE or the boundary condition.  
   
    The PDE is symmetrical if  
   
    M{A[i,j]=A[j,i]}  and M{B[j]=C[j]}  
   
    For a system of PDEs and a solution with several components the PDE has the form  
   
    M{-grad(A[i,j,k,l]*grad(u[k])[l]+B[i,j,k]*u[k])[j]+C[i,k,l]*grad(u[k])[l]+D[i,k]*u[k] =-grad(X[i,j])[j]+Y[i] }  
   
    M{A} is a ramk four, M{B} and M{C} are each a rank three, M{D} and M{X} are each a rank two and M{Y} is a rank one.  
    The natural boundary conditions take the form:  
   
    M{n[j]*(A[i,j,k,l]*grad(u[k])[l]+B[i,j,k]*u[k])+d[i,k]*u[k]=n[j]*X[i,j]+y[i]}  
   
   
    The coefficient M{d} is a rank two and M{y} is a  rank one both in the L{FunctionOnBoundary<escript.FunctionOnBoundary>}. Constraints take the form  
   
   
    M{u[i]=r[i]}  where  M{q[i]>0}  
   
    M{r} and M{q} are each rank one. Notice that at some locations not necessarily all components must have a constraint.  
   
    The system of PDEs is symmetrical if  
   
         - M{A[i,j,k,l]=A[k,l,i,j]}  
         - M{B[i,j,k]=C[k,i,j]}  
         - M{D[i,k]=D[i,k]}  
         - M{d[i,k]=d[k,i]}  
   
    L{LinearPDE} also supports solution discontinuities over a contact region in the domain. To specify the conditions across the  
    discontinuity we are using the generalised flux M{J} which is in the case of a systems of PDEs and several components of the solution  
    defined as  
   
    M{J[i,j]=A[i,j,k,l]*grad(u[k])[l]+B[i,j,k]*u[k]-X[i,j]}  
   
    For the case of single solution component and single PDE M{J} is defined  
   
    M{J_{j}=A[i,j]*grad(u)[j]+B[i]*u-X[i]}  
388    
389     In the context of discontinuities M{n} denotes the normal on the discontinuity pointing from side 0 towards side 1     The class assumes that some sort of assembling process is required to form
390     calculated from L{getNormal<escript.FunctionSpace.getNormal>} of L{FunctionOnContactZero<escript.FunctionOnContactZero>}. For a system of PDEs     a problem of the form
    the contact condition takes the form  
391    
392     M{n[j]*J0[i,j]=n[j]*J1[i,j]=y_contact[i]- d_contact[i,k]*jump(u)[k]}     M{L u=f}
393    
394     where M{J0} and M{J1} are the fluxes on side 0 and side 1 of the discontinuity, respectively. M{jump(u)}, which is the difference     where M{L} is an operator and M{f} is the right hand side. This operator
395     of the solution at side 1 and at side 0, denotes the jump of M{u} across discontinuity along the normal calcualted by     problem will be solved to get the unknown M{u}.
    L{jump<util.jump>}.  
    The coefficient M{d_contact} is a rank two and M{y_contact} is a rank one both in the L{FunctionOnContactZero<escript.FunctionOnContactZero>} or L{FunctionOnContactOne<escript.FunctionOnContactOne>}.  
    In case of a single PDE and a single component solution the contact condition takes the form  
396    
397     M{n[j]*J0_{j}=n[j]*J1_{j}=y_contact-d_contact*jump(u)}     @cvar DEFAULT: The default method used to solve the system of linear
398                      equations
    In this case the the coefficient M{d_contact} and M{y_contact} are eaach scalar  
    both in the L{FunctionOnContactZero<escript.FunctionOnContactZero>} or L{FunctionOnContactOne<escript.FunctionOnContactOne>}.  
   
    @cvar DEFAULT: The default method used to solve the system of linear equations  
399     @cvar DIRECT: The direct solver based on LDU factorization     @cvar DIRECT: The direct solver based on LDU factorization
400     @cvar CHOLEVSKY: The direct solver based on LDLt factorization (can only be applied for symmetric PDEs)     @cvar CHOLEVSKY: The direct solver based on LDLt factorization (can only be
401     @cvar PCG: The preconditioned conjugate gradient method (can only be applied for symmetric PDEs)                      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     @cvar CR: The conjugate residual method
405     @cvar CGS: The conjugate gardient square method     @cvar CGS: The conjugate gradient square method
406     @cvar BICGSTAB: The stabilized BiConjugate Gradient method.     @cvar BICGSTAB: The stabilized BiConjugate Gradient method
407     @cvar SSOR: The symmetric overrealaxtion method     @cvar TFQMR: Transport Free Quasi Minimal Residual method
408     @cvar ILU0: The incomplete LU factorization preconditioner  with no fill in     @cvar MINRES: Minimum residual method
409     @cvar ILUT: The incomplete LU factorization preconditioner with will in     @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     @cvar JACOBI: The Jacobi preconditioner
413     @cvar GMRES: The Gram-Schmidt minimum residual method     @cvar GMRES: The Gram-Schmidt minimum residual method
414     @cvar PRES20: Special GMRES with restart after 20 steps and truncation after 5 residuals     @cvar PRES20: Special GMRES with restart after 20 steps and truncation after
415     @cvar LUMPING: Matrix lumping.                   5 residuals
416       @cvar LUMPING: Matrix lumping
417     @cvar NO_REORDERING: No matrix reordering allowed     @cvar NO_REORDERING: No matrix reordering allowed
418     @cvar MINIMUM_FILL_IN: Reorder matrix to reduce fill-in during factorization     @cvar MINIMUM_FILL_IN: Reorder matrix to reduce fill-in during factorization
419     @cvar NESTED_DISSECTION: Reorder matrix to improve load balancing during factorization     @cvar NESTED_DISSECTION: Reorder matrix to improve load balancing during
420                                factorization
421     @cvar PASO: PASO solver package     @cvar PASO: PASO solver package
422     @cvar SCSL: SGI SCSL solver library     @cvar SCSL: SGI SCSL solver library
423     @cvar MKL: Intel's MKL solver library     @cvar MKL: Intel's MKL solver library
424     @cvar UMFPACK: the UMFPACK library     @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     @cvar ITERATIVE: The default iterative solver
428     @cvar AMG: algebraic multi grid     @cvar AMG: Algebraic Multi Grid
429     @cvar RILU: recursive ILU     @cvar RILU: Recursive ILU
430       @cvar GS: Gauss-Seidel solver
431    
432     """     """
433     DEFAULT= 0     DEFAULT= 0
# Line 441  class LinearPDE(object): Line 454  class LinearPDE(object):
454     PASO= 21     PASO= 21
455     AMG= 22     AMG= 22
456     RILU = 23     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     SMALL_TOLERANCE=1.e-13
464     __PACKAGE_KEY="package"     PACKAGE_KEY="package"
465     __METHOD_KEY="method"     METHOD_KEY="method"
466     __SYMMETRY_KEY="symmetric"     SYMMETRY_KEY="symmetric"
467     __TOLERANCE_KEY="tolerance"     TOLERANCE_KEY="tolerance"
468     __PRECONDITIONER_KEY="preconditioner"     PRECONDITIONER_KEY="preconditioner"
469    
470    
471     def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):     def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):
472       """       """
473       initializes a new linear PDE       Initializes a linear problem.
474    
475       @param domain: domain of the PDE       @param domain: domain of the PDE
476       @type domain: L{Domain<escript.Domain>}       @type domain: L{Domain<escript.Domain>}
477       @param numEquations: number of equations. If numEquations==None the number of equations       @param numEquations: number of equations. If C{None} the number of
478                            is exracted from the PDE coefficients.                            equations is extracted from the coefficients.
479       @param numSolutions: number of solution components. If  numSolutions==None the number of solution components       @param numSolutions: number of solution components. If C{None} the number
480                            is exracted from the PDE coefficients.                            of solution components is extracted from the
481       @param debug: if True debug informations are printed.                            coefficients.
482         @param debug: if True debug information is printed
483    
484       """       """
485       super(LinearPDE, self).__init__()       super(LinearProblem, self).__init__()
      #  
      #   the coefficients of the general PDE:  
      #  
      self.__COEFFICIENTS_OF_GENEARL_PDE={  
        "A"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM,PDECoefficient.BY_SOLUTION,PDECoefficient.BY_DIM),PDECoefficient.OPERATOR),  
        "B"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),  
        "C"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION,PDECoefficient.BY_DIM),PDECoefficient.OPERATOR),  
        "D"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),  
        "X"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM),PDECoefficient.RIGHTHANDSIDE),  
        "Y"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),  
        "d"         : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),  
        "y"         : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),  
        "d_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),  
        "y_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),  
        "r"         : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_SOLUTION,),PDECoefficient.RIGHTHANDSIDE),  
        "q"         : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_SOLUTION,),PDECoefficient.BOTH)}  
486    
      # COEFFICIENTS can be overwritten by subclasses:  
      self.COEFFICIENTS=self.__COEFFICIENTS_OF_GENEARL_PDE  
      self.__altered_coefficients=False  
      # initialize attributes  
487       self.__debug=debug       self.__debug=debug
488       self.__domain=domain       self.__domain=domain
489       self.__numEquations=numEquations       self.__numEquations=numEquations
490       self.__numSolutions=numSolutions       self.__numSolutions=numSolutions
491       self.__resetSystem()       self.__altered_coefficients=False
   
      # set some default values:  
492       self.__reduce_equation_order=False       self.__reduce_equation_order=False
493       self.__reduce_solution_order=False       self.__reduce_solution_order=False
494       self.__tolerance=1.e-8       self.__tolerance=1.e-8
495       self.__solver_method=self.DEFAULT       self.__solver_method=self.DEFAULT
496       self.__solver_package=self.DEFAULT       self.__solver_package=self.DEFAULT
497       self.__preconditioner=self.DEFAULT       self.__preconditioner=self.DEFAULT
498       self.__matrix_type=self.__domain.getSystemMatrixTypeId(self.DEFAULT,self.DEFAULT,False)       self.__is_RHS_valid=False
499         self.__is_operator_valid=False
500       self.__sym=False       self.__sym=False
501         self.__COEFFICIENTS={}
502       self.resetCoefficients()       # initialize things:
503       self.trace("PDE Coeffients are %s"%str(self.COEFFICIENTS.keys()))       self.resetAllCoefficients()
504     # =============================================================================       self.__system_type=None
505         self.updateSystemType()
506       # ==========================================================================
507     #    general stuff:     #    general stuff:
508     # =============================================================================     # ==========================================================================
509     def __str__(self):     def __str__(self):
510       """       """
511       returns string representation of the PDE       Returns a string representation of the PDE.
512    
513       @return: a simple representation of the PDE       @return: a simple representation of the PDE
514       @rtype: C{str}       @rtype: C{str}
515       """       """
516       return "<LinearPDE %d>"%id(self)       return "<LinearProblem %d>"%id(self)
517     # =============================================================================     # ==========================================================================
518     #    debug :     #    debug :
519     # =============================================================================     # ==========================================================================
520     def setDebugOn(self):     def setDebugOn(self):
521       """       """
522       switches on debugging       Switches debug output on.
523       """       """
524       self.__debug=not None       self.__debug=not None
525    
526     def setDebugOff(self):     def setDebugOff(self):
527       """       """
528       switches off debugging       Switches debug output off.
529       """       """
530       self.__debug=None       self.__debug=None
531    
532       def setDebug(self, flag):
533         """
534         Switches debug output on if C{flag} is True otherwise it is switched off.
535    
536         @param flag: desired debug status
537         @type flag: C{bool}
538         """
539         if flag:
540             self.setDebugOn()
541         else:
542             self.setDebugOff()
543    
544     def trace(self,text):     def trace(self,text):
545       """       """
546       print the text message if debugging is swiched on.       Prints the text message if debug mode is switched on.
547       @param text: message  
548         @param text: message to be printed
549       @type text: C{string}       @type text: C{string}
550       """       """
551       if self.__debug: print "%s: %s"%(str(self),text)       if self.__debug: print "%s: %s"%(str(self),text)
552    
553     # =============================================================================     # ==========================================================================
554     # some service functions:     # some service functions:
555     # =============================================================================     # ==========================================================================
556       def introduceCoefficients(self,**coeff):
557           """
558           Introduces new coefficients into the problem.
559    
560           Use::
561    
562           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):     def getDomain(self):
574       """       """
575       returns the domain of the PDE       Returns the domain of the PDE.
576    
577       @return: the domain of the PDE       @return: the domain of the PDE
578       @rtype: L{Domain<escript.Domain>}       @rtype: L{Domain<escript.Domain>}
# Line 551  class LinearPDE(object): Line 581  class LinearPDE(object):
581    
582     def getDim(self):     def getDim(self):
583       """       """
584       returns the spatial dimension of the PDE       Returns the spatial dimension of the PDE.
585    
586       @return: the spatial dimension of the PDE domain       @return: the spatial dimension of the PDE domain
587       @rtype: C{int}       @rtype: C{int}
# Line 560  class LinearPDE(object): Line 590  class LinearPDE(object):
590    
591     def getNumEquations(self):     def getNumEquations(self):
592       """       """
593       returns the number of equations       Returns the number of equations.
594    
595       @return: the number of equations       @return: the number of equations
596       @rtype: C{int}       @rtype: C{int}
597       @raise UndefinedPDEError: if the number of equations is not be specified yet.       @raise UndefinedPDEError: if the number of equations is not specified yet
598       """       """
599       if self.__numEquations==None:       if self.__numEquations==None:
600           raise UndefinedPDEError,"Number of equations is undefined. Please specify argument numEquations."           if self.__numSolutions==None:
601       else:              raise UndefinedPDEError,"Number of equations is undefined. Please specify argument numEquations."
602           return self.__numEquations           else:
603                self.__numEquations=self.__numSolutions
604         return self.__numEquations
605    
606     def getNumSolutions(self):     def getNumSolutions(self):
607       """       """
608       returns the number of unknowns       Returns the number of unknowns.
609    
610       @return: the number of unknowns       @return: the number of unknowns
611       @rtype: C{int}       @rtype: C{int}
612       @raise UndefinedPDEError: if the number of unknowns is not be specified yet.       @raise UndefinedPDEError: if the number of unknowns is not specified yet
613       """       """
614       if self.__numSolutions==None:       if self.__numSolutions==None:
615          raise UndefinedPDEError,"Number of solution is undefined. Please specify argument numSolutions."          if self.__numEquations==None:
616       else:              raise UndefinedPDEError,"Number of solution is undefined. Please specify argument numSolutions."
617          return self.__numSolutions          else:
618                self.__numSolutions=self.__numEquations
619         return self.__numSolutions
620    
621     def reduceEquationOrder(self):     def reduceEquationOrder(self):
622       """       """
623       return status for order reduction for equation       Returns the status of order reduction for the equation.
624    
625       @return: return True is reduced interpolation order is used for the represenation of the equation       @return: True if reduced interpolation order is used for the
626                  representation of the equation, False otherwise
627       @rtype: L{bool}       @rtype: L{bool}
628       """       """
629       return self.__reduce_equation_order       return self.__reduce_equation_order
630    
631     def reduceSolutionOrder(self):     def reduceSolutionOrder(self):
632       """       """
633       return status for order reduction for the solution       Returns the status of order reduction for the solution.
634    
635       @return: return True is reduced interpolation order is used for the represenation of the solution       @return: True if reduced interpolation order is used for the
636                  representation of the solution, False otherwise
637       @rtype: L{bool}       @rtype: L{bool}
638       """       """
639       return self.__reduce_solution_order       return self.__reduce_solution_order
640    
641     def getFunctionSpaceForEquation(self):     def getFunctionSpaceForEquation(self):
642       """       """
643       returns the L{FunctionSpace<escript.FunctionSpace>} used to discretize the equation       Returns the L{FunctionSpace<escript.FunctionSpace>} used to discretize
644         the equation.
645    
646       @return: representation space of equation       @return: representation space of equation
647       @rtype: L{FunctionSpace<escript.FunctionSpace>}       @rtype: L{FunctionSpace<escript.FunctionSpace>}
# Line 616  class LinearPDE(object): Line 653  class LinearPDE(object):
653    
654     def getFunctionSpaceForSolution(self):     def getFunctionSpaceForSolution(self):
655       """       """
656       returns the L{FunctionSpace<escript.FunctionSpace>} used to represent the solution       Returns the L{FunctionSpace<escript.FunctionSpace>} used to represent
657         the solution.
658    
659       @return: representation space of solution       @return: representation space of solution
660       @rtype: L{FunctionSpace<escript.FunctionSpace>}       @rtype: L{FunctionSpace<escript.FunctionSpace>}
# Line 626  class LinearPDE(object): Line 664  class LinearPDE(object):
664       else:       else:
665           return escript.Solution(self.getDomain())           return escript.Solution(self.getDomain())
666    
667       # ==========================================================================
    def getOperator(self):  
      """  
      provides access to the operator of the PDE  
   
      @return: the operator of the PDE  
      @rtype: L{Operator<escript.Operator>}  
      """  
      m=self.getSystem()[0]  
      if self.isUsingLumping():  
          return self.copyConstraint(1./m)  
      else:  
          return m  
   
    def getRightHandSide(self):  
      """  
      provides access to the right hand side of the PDE  
      @return: the right hand side of the PDE  
      @rtype: L{Data<escript.Data>}  
      """  
      r=self.getSystem()[1]  
      if self.isUsingLumping():  
          return self.copyConstraint(r)  
      else:  
          return r  
   
    def applyOperator(self,u=None):  
      """  
      applies the operator of the PDE to a given u or the solution of PDE if u is not present.  
   
      @param u: argument of the operator. It must be representable in C{elf.getFunctionSpaceForSolution()}. If u is not present or equals L{None}  
                the current solution is used.  
      @type u: L{Data<escript.Data>} or None  
      @return: image of u  
      @rtype: L{Data<escript.Data>}  
      """  
      if u==None:  
           return self.getOperator()*self.getSolution()  
      else:  
         self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())  
   
    def getResidual(self,u=None):  
      """  
      return the residual of u or the current solution if u is not present.  
   
      @param u: argument in the residual calculation. It must be representable in C{elf.getFunctionSpaceForSolution()}. If u is not present or equals L{None}  
                the current solution is used.  
      @type u: L{Data<escript.Data>} or None  
      @return: residual of u  
      @rtype: L{Data<escript.Data>}  
      """  
      return self.applyOperator(u)-self.getRightHandSide()  
   
    def checkSymmetry(self,verbose=True):  
       """  
       test the PDE for symmetry.  
   
       @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.  
       @type verbose: C{bool}  
       @return:  True if the PDE is symmetric.  
       @rtype: L{Data<escript.Data>}  
       @note: This is a very expensive operation. It should be used for degugging only! The symmetry flag is not altered.  
       """  
       verbose=verbose or self.__debug  
       out=True  
       if self.getNumSolutions()!=self.getNumEquations():  
          if verbose: print "non-symmetric PDE because of different number of equations and solutions"  
          out=False  
       else:  
          A=self.getCoefficientOfGeneralPDE("A")  
          if not A.isEmpty():  
             tol=util.Lsup(A)*self.SMALL_TOLERANCE  
             if self.getNumSolutions()>1:  
                for i in range(self.getNumEquations()):  
                   for j in range(self.getDim()):  
                      for k in range(self.getNumSolutions()):  
                         for l in range(self.getDim()):  
                             if util.Lsup(A[i,j,k,l]-A[k,l,i,j])>tol:  
                                if verbose: print "non-symmetric PDE because A[%d,%d,%d,%d]!=A[%d,%d,%d,%d]"%(i,j,k,l,k,l,i,j)  
                                out=False  
             else:  
                for j in range(self.getDim()):  
                   for l in range(self.getDim()):  
                      if util.Lsup(A[j,l]-A[l,j])>tol:  
                         if verbose: print "non-symmetric PDE because A[%d,%d]!=A[%d,%d]"%(j,l,l,j)  
                         out=False  
          B=self.getCoefficientOfGeneralPDE("B")  
          C=self.getCoefficientOfGeneralPDE("C")  
          if B.isEmpty() and not C.isEmpty():  
             if verbose: print "non-symmetric PDE because B is not present but C is"  
             out=False  
          elif not B.isEmpty() and C.isEmpty():  
             if verbose: print "non-symmetric PDE because C is not present but B is"  
             out=False  
          elif not B.isEmpty() and not C.isEmpty():  
             tol=(util.Lsup(B)+util.Lsup(C))*self.SMALL_TOLERANCE/2.  
             if self.getNumSolutions()>1:  
                for i in range(self.getNumEquations()):  
                    for j in range(self.getDim()):  
                       for k in range(self.getNumSolutions()):  
                          if util.Lsup(B[i,j,k]-C[k,i,j])>tol:  
                               if verbose: print "non-symmetric PDE because B[%d,%d,%d]!=C[%d,%d,%d]"%(i,j,k,k,i,j)  
                               out=False  
             else:  
                for j in range(self.getDim()):  
                   if util.Lsup(B[j]-C[j])>tol:  
                      if verbose: print "non-symmetric PDE because B[%d]!=C[%d]"%(j,j)  
                      out=False  
          if self.getNumSolutions()>1:  
            D=self.getCoefficientOfGeneralPDE("D")  
            if not D.isEmpty():  
              tol=util.Lsup(D)*self.SMALL_TOLERANCE  
              for i in range(self.getNumEquations()):  
                 for k in range(self.getNumSolutions()):  
                   if util.Lsup(D[i,k]-D[k,i])>tol:  
                       if verbose: print "non-symmetric PDE because D[%d,%d]!=D[%d,%d]"%(i,k,k,i)  
                       out=False  
            d=self.getCoefficientOfGeneralPDE("d")  
            if not d.isEmpty():  
              tol=util.Lsup(d)*self.SMALL_TOLERANCE  
              for i in range(self.getNumEquations()):  
                 for k in range(self.getNumSolutions()):  
                   if util.Lsup(d[i,k]-d[k,i])>tol:  
                       if verbose: print "non-symmetric PDE because d[%d,%d]!=d[%d,%d]"%(i,k,k,i)  
                       out=False  
            d_contact=self.getCoefficientOfGeneralPDE("d_contact")  
            if not d_contact.isEmpty():  
              tol=util.Lsup(d_contact)*self.SMALL_TOLERANCE  
              for i in range(self.getNumEquations()):  
                 for k in range(self.getNumSolutions()):  
                   if util.Lsup(d_contact[i,k]-d_contact[k,i])>tol:  
                       if verbose: print "non-symmetric PDE because d_contact[%d,%d]!=d_contact[%d,%d]"%(i,k,k,i)  
                       out=False  
       return out  
   
    def getSolution(self,**options):  
        """  
        returns the solution of the PDE. If the solution is not valid the PDE is solved.  
   
        @return: the solution  
        @rtype: L{Data<escript.Data>}  
        @param options: solver options  
        @keyword verbose: True to get some information during PDE solution  
        @type verbose: C{bool}  
        @keyword reordering: reordering scheme to be used during elimination. Allowed values are  
                             L{NO_REORDERING}, L{MINIMUM_FILL_IN}, L{NESTED_DISSECTION}  
        @keyword iter_max: maximum number of iteration steps allowed.  
        @keyword drop_tolerance: threshold for drupping in L{ILUT}  
        @keyword drop_storage: maximum of allowed memory in L{ILUT}  
        @keyword truncation: maximum number of residuals in L{GMRES}  
        @keyword restart: restart cycle length in L{GMRES}  
        """  
        if not self.__solution_isValid:  
           mat,f=self.getSystem()  
           if self.isUsingLumping():  
              self.__solution=self.copyConstraint(f*mat)  
           else:  
              options[self.__TOLERANCE_KEY]=self.getTolerance()  
              options[self.__METHOD_KEY]=self.getSolverMethod()[0]  
              options[self.__PRECONDITIONER_KEY]=self.getSolverMethod()[1]  
              options[self.__PACKAGE_KEY]=self.getSolverPackage()  
              options[self.__SYMMETRY_KEY]=self.isSymmetric()  
              self.trace("PDE is resolved.")  
              self.trace("solver options: %s"%str(options))  
              self.__solution=mat.solve(f,options)  
           self.__solution_isValid=True  
        return self.__solution  
   
    def getFlux(self,u=None):  
      """  
      returns the flux M{J} for a given M{u}  
   
      M{J[i,j]=A[i,j,k,l]*grad(u[k])[l]+B[i,j,k]u[k]-X[i,j]}  
   
      or  
   
      M{J[j]=A[i,j]*grad(u)[l]+B[j]u-X[j]}  
   
      @param u: argument in the flux. If u is not present or equals L{None} the current solution is used.  
      @type u: L{Data<escript.Data>} or None  
      @return: flux  
      @rtype: L{Data<escript.Data>}  
      """  
      if u==None: u=self.getSolution()  
      return util.tensormult(self.getCoefficientOfGeneralPDE("A"),util.grad(u))+util.matrixmult(self.getCoefficientOfGeneralPDE("B"),u)-util.self.getCoefficientOfGeneralPDE("X")  
    # =============================================================================  
668     #   solver settings:     #   solver settings:
669     # =============================================================================     # ==========================================================================
670     def setSolverMethod(self,solver=None,preconditioner=None):     def setSolverMethod(self,solver=None,preconditioner=None):
671         """         """
672         sets a new solver         Sets a new solver method and/or preconditioner.
673    
674         @param solver: sets a new solver method.         @param solver: new solver method to use
675         @type solver: one of L{DEFAULT}, L{ITERATIVE} L{DIRECT}, L{CHOLEVSKY}, L{PCG}, L{CR}, L{CGS}, L{BICGSTAB}, L{SSOR}, L{GMRES}, L{PRES20}, L{LUMPING}, L{AMG}         @type solver: one of L{DEFAULT}, L{ITERATIVE} L{DIRECT}, L{CHOLEVSKY},
676         @param preconditioner: sets a new solver method.                       L{PCG}, L{CR}, L{CGS}, L{BICGSTAB}, L{SSOR}, L{GMRES},
677         @type preconditioner: one of L{DEFAULT}, L{JACOBI} L{ILU0}, L{ILUT},L{SSOR}, L{RILU}                       L{TFQMR}, L{MINRES}, L{PRES20}, L{LUMPING}, L{AMG}
678         """         @param preconditioner: new preconditioner to use
679         if solver==None: solve=self.DEFAULT         @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         if preconditioner==None: preconditioner=self.DEFAULT
689         if not (solver,preconditioner)==self.getSolverMethod():         if not (solver,preconditioner)==self.getSolverMethod():
690             self.__solver_method=solver             self.__solver_method=solver
691             self.__preconditioner=preconditioner             self.__preconditioner=preconditioner
692             self.__checkMatrixType()             self.updateSystemType()
693             self.trace("New solver is %s"%self.getSolverMethodName())             self.trace("New solver is %s"%self.getSolverMethodName())
694    
695     def getSolverMethodName(self):     def getSolverMethodName(self):
696         """         """
697         returns the name of the solver currently used         Returns the name of the solver currently used.
698    
699         @return: the name of the solver currently used.         @return: the name of the solver currently used
700         @rtype: C{string}         @rtype: C{string}
701         """         """
702    
# Line 847  class LinearPDE(object): Line 708  class LinearPDE(object):
708         elif m[0]==self.ITERATIVE: method= "ITERATIVE"         elif m[0]==self.ITERATIVE: method= "ITERATIVE"
709         elif m[0]==self.CHOLEVSKY: method= "CHOLEVSKY"         elif m[0]==self.CHOLEVSKY: method= "CHOLEVSKY"
710         elif m[0]==self.PCG: method= "PCG"         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"         elif m[0]==self.CR: method= "CR"
714         elif m[0]==self.CGS: method= "CGS"         elif m[0]==self.CGS: method= "CGS"
715         elif m[0]==self.BICGSTAB: method= "BICGSTAB"         elif m[0]==self.BICGSTAB: method= "BICGSTAB"
# Line 862  class LinearPDE(object): Line 725  class LinearPDE(object):
725         elif m[1]==self.SSOR: method+= "+SSOR"         elif m[1]==self.SSOR: method+= "+SSOR"
726         elif m[1]==self.AMG: method+= "+AMG"         elif m[1]==self.AMG: method+= "+AMG"
727         elif m[1]==self.RILU: method+= "+RILU"         elif m[1]==self.RILU: method+= "+RILU"
728           elif m[1]==self.GS: method+= "+GS"
729         if p==self.DEFAULT: package="DEFAULT"         if p==self.DEFAULT: package="DEFAULT"
730         elif p==self.PASO: package= "PASO"         elif p==self.PASO: package= "PASO"
731         elif p==self.MKL: package= "MKL"         elif p==self.MKL: package= "MKL"
732         elif p==self.SCSL: package= "SCSL"         elif p==self.SCSL: package= "SCSL"
733         elif p==self.UMFPACK: package= "UMFPACK"         elif p==self.UMFPACK: package= "UMFPACK"
734           elif p==self.TRILINOS: package= "TRILINOS"
735         else : method="unknown"         else : method="unknown"
736         return "%s solver of %s package"%(method,package)         return "%s solver of %s package"%(method,package)
737    
   
738     def getSolverMethod(self):     def getSolverMethod(self):
739         """         """
740         returns the solver method         Returns the solver method and preconditioner used.
741    
742         @return: the solver method currently be used.         @return: the solver method currently used.
743         @rtype: C{int}         @rtype: C{tuple} of C{int}
744         """         """
745         return self.__solver_method,self.__preconditioner         return self.__solver_method,self.__preconditioner
746    
747     def setSolverPackage(self,package=None):     def setSolverPackage(self,package=None):
748         """         """
749         sets a new solver package         Sets a new solver package.
750    
751         @param solver: sets a new solver method.         @param package: the new solver package
752         @type solver: one of L{DEFAULT}, L{PASO} L{SCSL}, L{MKL}, L{UMLPACK}         @type package: one of L{DEFAULT}, L{PASO} L{SCSL}, L{MKL}, L{UMFPACK},
753                          L{TRILINOS}
754         """         """
755         if package==None: package=self.DEFAULT         if package==None: package=self.DEFAULT
756         if not package==self.getSolverPackage():         if not package==self.getSolverPackage():
757             self.__solver_method=solver             self.__solver_package=package
758             self.__checkMatrixType()             self.updateSystemType()
759             self.trace("New solver is %s"%self.getSolverMethodName())             self.trace("New solver is %s"%self.getSolverMethodName())
760    
761     def getSolverPackage(self):     def getSolverPackage(self):
762         """         """
763         returns the package of the solver         Returns the package of the solver.
764    
765         @return: the solver package currently being used.         @return: the solver package currently being used
766         @rtype: C{int}         @rtype: C{int}
767         """         """
768         return self.__solver_package         return self.__solver_package
769    
770     def isUsingLumping(self):     def isUsingLumping(self):
771        """        """
772        checks if matrix lumping is used a solver method        Checks if matrix lumping is the current solver method.
773    
774        @return: True is lumping is currently used a solver method.        @return: True if the current solver method is lumping
775        @rtype: C{bool}        @rtype: C{bool}
776        """        """
777        return self.getSolverMethod()[0]==self.LUMPING        return self.getSolverMethod()[0]==self.LUMPING
778    
779     def setTolerance(self,tol=1.e-8):     def setTolerance(self,tol=1.e-8):
780         """         """
781         resets the tolerance for the solver method to tol where for an appropriate norm M{|.|}         Resets the tolerance for the solver method to C{tol}.
782    
783         M{|L{getResidual}()|<tol*|L{getRightHandSide}()|}         @param tol: new tolerance for the solver. If C{tol} is lower than the
784                       current tolerence the system will be resolved.
        defines the stopping criterion.  
   
        @param tol: new tolerance for the solver. If the tol is lower then the current tolerence  
                    the system will be resolved.  
785         @type tol: positive C{float}         @type tol: positive C{float}
786         @raise ValueException: if tolerance is not positive.         @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.__invalidateSolution()         if tol<self.getTolerance(): self.invalidateSolution()
791         self.trace("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         returns the tolerance set for the solution         Returns the tolerance currently set for the solution.
798    
799         @return: tolerance currently used.         @return: tolerance currently used
800         @rtype: C{float}         @rtype: C{float}
801         """         """
802         return self.__tolerance         return self.__tolerance
803    
804     # =============================================================================     # ==========================================================================
805     #    symmetry  flag:     #    symmetry  flag:
806     # =============================================================================     # ==========================================================================
807     def isSymmetric(self):     def isSymmetric(self):
808        """        """
809        checks if symmetry is indicated.        Checks if symmetry is indicated.
810    
811        @return: True is a symmetric PDE is indicated, otherwise False is returned        @return: True if a symmetric PDE is indicated, False otherwise
812        @rtype: C{bool}        @rtype: C{bool}
813        """        """
814        return self.__sym        return self.__sym
815    
816     def setSymmetryOn(self):     def setSymmetryOn(self):
817        """        """
818        sets the symmetry flag.        Sets the symmetry flag.
819        """        """
820        if not self.isSymmetric():        if not self.isSymmetric():
821           self.trace("PDE 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        removes the symmetry flag.        Clears the symmetry flag.
828        """        """
829        if self.isSymmetric():        if self.isSymmetric():
830           self.trace("PDE 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        sets the symmetry flag to flag        Sets the symmetry flag to C{flag}.
837    
838        @param flag: If flag, the symmetry flag is set otherwise the symmetry flag is released.        @param flag: If True, the symmetry flag is set otherwise reset.
839        @type flag: C{bool}        @type flag: C{bool}
840        """        """
841        if flag:        if flag:
# Line 982  class LinearPDE(object): Line 843  class LinearPDE(object):
843        else:        else:
844           self.setSymmetryOff()           self.setSymmetryOff()
845    
846     # =============================================================================     # ==========================================================================
847     # function space handling for the equation as well as the solution     # function space handling for the equation as well as the solution
848     # =============================================================================     # ==========================================================================
849     def setReducedOrderOn(self):     def setReducedOrderOn(self):
850       """       """
851       switches on reduced order for solution and equation representation       Switches reduced order on for solution and equation representation.
852    
853       @raise RuntimeError: if order reduction is altered after a coefficient has been set.       @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       switches off reduced order for solution and equation representation       Switches reduced order off for solution and equation representation
862    
863       @raise RuntimeError: if order reduction is altered after a coefficient has been set.       @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       sets order reduction for both solution and equation representation according to flag.       Sets order reduction state for both solution and equation representation
872       @param flag: if flag is True, the order reduction is switched on for both  solution and equation representation, otherwise or       according to flag.
873                    if flag is not present order reduction is switched off  
874         @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}       @type flag: C{bool}
878       @raise RuntimeError: if order reduction is altered after a coefficient has been set.       @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)
# Line 1017  class LinearPDE(object): Line 884  class LinearPDE(object):
884    
885     def setReducedOrderForSolutionOn(self):     def setReducedOrderForSolutionOn(self):
886       """       """
887       switches on reduced order for solution representation       Switches reduced order on for solution representation.
888    
889       @raise RuntimeError: if order reduction is altered after a coefficient has been set.       @raise RuntimeError: if order reduction is altered after a coefficient has
890                              been set
891       """       """
892       if not self.__reduce_solution_order:       if not self.__reduce_solution_order:
893           if self.__altered_coefficients:           if self.__altered_coefficients:
894                raise RuntimeError,"order cannot be altered after coefficients have been defined."                raise RuntimeError,"order cannot be altered after coefficients have been defined."
895           self.trace("Reduced order is used to solution representation.")           self.trace("Reduced order is used for solution representation.")
896           self.__reduce_solution_order=True           self.__reduce_solution_order=True
897           self.__resetSystem()           self.initializeSystem()
898    
899     def setReducedOrderForSolutionOff(self):     def setReducedOrderForSolutionOff(self):
900       """       """
901       switches off reduced order for solution representation       Switches reduced order off for solution representation
902    
903       @raise RuntimeError: if order reduction is altered after a coefficient has been set.       @raise RuntimeError: if order reduction is altered after a coefficient has
904                              been set.
905       """       """
906       if self.__reduce_solution_order:       if self.__reduce_solution_order:
907           if self.__altered_coefficients:           if self.__altered_coefficients:
908                raise RuntimeError,"order cannot be altered after coefficients have been defined."                raise RuntimeError,"order cannot be altered after coefficients have been defined."
909           self.trace("Full order is used to interpolate solution.")           self.trace("Full order is used to interpolate solution.")
910           self.__reduce_solution_order=False           self.__reduce_solution_order=False
911           self.__resetSystem()           self.initializeSystem()
912    
913     def setReducedOrderForSolutionTo(self,flag=False):     def setReducedOrderForSolutionTo(self,flag=False):
914       """       """
915       sets order for test functions according to flag       Sets order reduction state for solution representation according to flag.
916    
917       @param flag: if flag is True, the order reduction is switched on for solution representation, otherwise or       @param flag: if flag is True, the order reduction is switched on for
918                    if flag is not present order reduction is switched off                    solution representation, otherwise or if flag is not present
919                      order reduction is switched off
920       @type flag: C{bool}       @type flag: C{bool}
921       @raise RuntimeError: if order reduction is altered after a coefficient has been set.       @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()
# Line 1057  class LinearPDE(object): Line 928  class LinearPDE(object):
928    
929     def setReducedOrderForEquationOn(self):     def setReducedOrderForEquationOn(self):
930       """       """
931       switches on reduced order for equation representation       Switches reduced order on for equation representation.
932    
933       @raise RuntimeError: if order reduction is altered after a coefficient has been set.       @raise RuntimeError: if order reduction is altered after a coefficient has
934                              been set
935       """       """
936       if not self.__reduce_equation_order:       if not self.__reduce_equation_order:
937           if self.__altered_coefficients:           if self.__altered_coefficients:
938                raise RuntimeError,"order cannot be altered after coefficients have been defined."                raise RuntimeError,"order cannot be altered after coefficients have been defined."
939           self.trace("Reduced order is used for test functions.")           self.trace("Reduced order is used for test functions.")
940           self.__reduce_equation_order=True           self.__reduce_equation_order=True
941           self.__resetSystem()           self.initializeSystem()
942    
943     def setReducedOrderForEquationOff(self):     def setReducedOrderForEquationOff(self):
944       """       """
945       switches off reduced order for equation representation       Switches reduced order off for equation representation.
946    
947       @raise RuntimeError: if order reduction is altered after a coefficient has been set.       @raise RuntimeError: if order reduction is altered after a coefficient has
948                              been set
949       """       """
950       if self.__reduce_equation_order:       if self.__reduce_equation_order:
951           if self.__altered_coefficients:           if self.__altered_coefficients:
952                raise RuntimeError,"order cannot be altered after coefficients have been defined."                raise RuntimeError,"order cannot be altered after coefficients have been defined."
953           self.trace("Full order is used for test functions.")           self.trace("Full order is used for test functions.")
954           self.__reduce_equation_order=False           self.__reduce_equation_order=False
955           self.__resetSystem()           self.initializeSystem()
956    
957     def setReducedOrderForEquationTo(self,flag=False):     def setReducedOrderForEquationTo(self,flag=False):
958       """       """
959       sets order for test functions according to flag       Sets order reduction state for equation representation according to flag.
960    
961       @param flag: if flag is True, the order reduction is switched on for equation representation, otherwise or       @param flag: if flag is True, the order reduction is switched on for
962                    if flag is not present order reduction is switched off                    equation representation, otherwise or if flag is not present
963                      order reduction is switched off
964       @type flag: C{bool}       @type flag: C{bool}
965       @raise RuntimeError: if order reduction is altered after a coefficient has been set.       @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     # =============================================================================     def updateSystemType(self):
    # private method:  
    # =============================================================================  
    def __checkMatrixType(self):  
974       """       """
975       reassess the matrix type and, if a new matrix is needed, resets the system.       Reassesses the system type and, if a new matrix is needed, resets the
976         system.
977       """       """
978       new_matrix_type=self.getDomain().getSystemMatrixTypeId(self.getSolverMethod()[0],self.getSolverPackage(),self.isSymmetric())       new_system_type=self.getRequiredSystemType()
979       if not new_matrix_type==self.__matrix_type:       if not new_system_type==self.__system_type:
980           self.trace("Matrix type is now %d."%new_matrix_type)           self.trace("Matrix type is now %d."%new_system_type)
981           self.__matrix_type=new_matrix_type           self.__system_type=new_system_type
982           self.__resetSystem()           self.initializeSystem()
    #  
    #   rebuild switches :  
    #  
    def __invalidateSolution(self):  
        """  
        indicates the PDE has to be resolved if the solution is requested  
        """  
        if self.__solution_isValid: self.trace("PDE has to be resolved.")  
        self.__solution_isValid=False  
   
    def __invalidateOperator(self):  
        """  
        indicates the operator has to be rebuilt next time it is used  
        """  
        if self.__operator_is_Valid: self.trace("Operator has to be rebuilt.")  
        self.__invalidateSolution()  
        self.__operator_is_Valid=False  
983    
984     def __invalidateRightHandSide(self):     def getSystemType(self):
985         """        """
986         indicates the right hand side has to be rebuild next time it is used        Returns the current system type.
987         """        """
988         if self.__righthandside_isValid: self.trace("Right hand side has to be rebuilt.")        return self.__system_type
        self.__invalidateSolution()  
        self.__righthandside_isValid=False  
   
    def __invalidateSystem(self):  
        """  
        annonced that everthing has to be rebuild:  
        """  
        if self.__righthandside_isValid: self.trace("System has to be rebuilt.")  
        self.__invalidateSolution()  
        self.__invalidateOperator()  
        self.__invalidateRightHandSide()  
   
    def __resetSystem(self):  
        """  
        annonced that everthing has to be rebuild:  
        """  
        self.trace("New System is built from scratch.")  
        self.__operator=escript.Operator()  
        self.__operator_is_Valid=False  
        self.__righthandside=escript.Data()  
        self.__righthandside_isValid=False  
        self.__solution=escript.Data()  
        self.__solution_isValid=False  
    #  
    #    system initialization:  
    #  
    def __getNewOperator(self):  
        """  
        returns an instance of a new operator  
        """  
        self.trace("New operator is allocated.")  
        return self.getDomain().newOperator( \  
                            self.getNumEquations(), \  
                            self.getFunctionSpaceForEquation(), \  
                            self.getNumSolutions(), \  
                            self.getFunctionSpaceForSolution(), \  
                            self.__matrix_type)  
   
    def __getNewRightHandSide(self):  
        """  
        returns an instance of a new right hand side  
        """  
        self.trace("New right hand side is allocated.")  
        if self.getNumEquations()>1:  
            return escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)  
        else:  
            return escript.Data(0.,(),self.getFunctionSpaceForEquation(),True)  
   
    def __getNewSolution(self):  
        """  
        returns an instance of a new solution  
        """  
        self.trace("New solution is allocated.")  
        if self.getNumSolutions()>1:  
            return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)  
        else:  
            return escript.Data(0.,(),self.getFunctionSpaceForSolution(),True)  
989    
990     def __makeFreshSolution(self):     def checkSymmetricTensor(self,name,verbose=True):
991         """        """
992         makes sure that the solution is instantiated and returns it initialized by zeros        Tests a coefficient for symmetry.
        """  
        if self.__solution.isEmpty():  
            self.__solution=self.__getNewSolution()  
        else:  
            self.__solution*=0  
            self.trace("Solution is reset to zero.")  
        return self.__solution  
993    
994     def __makeFreshRightHandSide(self):        @param name: name of the coefficient
995         """        @type name: C{str}
996         makes sure that the right hand side is instantiated and returns it initialized by zeros        @param verbose: if set to True or not present a report on coefficients
997         """                        which break the symmetry is printed.
998         if self.__righthandside.isEmpty():        @type verbose: C{bool}
999             self.__righthandside=self.__getNewRightHandSide()        @return: True if coefficient C{name} is symmetric
1000         else:        @rtype: C{bool}
1001             self.__righthandside*=0        """
1002             self.trace("Right hand side is reset to zero.")        A=self.getCoefficient(name)
1003         return self.__righthandside        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 __makeFreshOperator(self):     def checkReciprocalSymmetry(self,name0,name1,verbose=True):
1037         """        """
1038         makes sure that the operator is instantiated and returns it initialized by zeros        Tests two coefficients for reciprocal symmetry.
        """  
        if self.__operator.isEmpty():  
            self.__operator=self.__getNewOperator()  
        else:  
            self.__operator.resetValues()  
            self.trace("Operator reset to zero")  
        return self.__operator  
1039    
1040     def __applyConstraint(self):        @param name0: name of the first coefficient
1041         """        @type name0: C{str}
1042         applies the constraints defined by q and r to the system        @param name1: name of the second coefficient
1043         """        @type name1: C{str}
1044         if not self.isUsingLumping():        @param verbose: if set to True or not present a report on coefficients
1045            q=self.getCoefficientOfGeneralPDE("q")                        which break the symmetry is printed
1046            r=self.getCoefficientOfGeneralPDE("r")        @type verbose: C{bool}
1047            if not q.isEmpty() and not self.__operator.isEmpty():        @return: True if coefficients C{name0} and C{name1} are reciprocally
1048               # q is the row and column mask to indicate where constraints are set:                 symmetric.
1049               row_q=escript.Data(q,self.getFunctionSpaceForEquation())        @rtype: C{bool}
1050               col_q=escript.Data(q,self.getFunctionSpaceForSolution())        """
1051               u=self.__getNewSolution()        B=self.getCoefficient(name0)
1052               if r.isEmpty():        C=self.getCoefficient(name1)
1053                  r_s=self.__getNewSolution()        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:               else:
1092                  r_s=escript.Data(r,self.getFunctionSpaceForSolution())                   raise ValueError,"Cannot check rank %s of %s and %s."%(len(sB),name0,name1)
1093               u.copyWithMask(r_s,col_q)        return out
              if not self.__righthandside.isEmpty():  
                 self.__righthandside-=self.__operator*u  
                 self.__righthandside=self.copyConstraint(self.__righthandside)  
              self.__operator.nullifyRowsAndCols(row_q,col_q,1.)  
    # =============================================================================  
    # function giving access to coefficients of the general PDE:  
    # =============================================================================  
    def getCoefficientOfGeneralPDE(self,name):  
      """  
      return the value of the coefficient name of the general PDE.  
   
      @note: This method is called by the assembling routine it can be overwritten  
            to map coefficients of a particular PDE to the general PDE.  
      @param name: name of the coefficient requested.  
      @type name: C{string}  
      @return: the value of the coefficient  name  
      @rtype: L{Data<escript.Data>}  
      @raise IllegalCoefficient: if name is not one of coefficients  
                   M{A}, M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact}, M{r} or M{q}.  
      """  
      if self.hasCoefficientOfGeneralPDE(name):  
         return self.getCoefficient(name)  
      else:  
         raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name  
   
    def hasCoefficientOfGeneralPDE(self,name):  
      """  
      checks if name is a the name of a coefficient of the general PDE.  
   
      @param name: name of the coefficient enquired.  
      @type name: C{string}  
      @return: True if name is the name of a coefficient of the general PDE. Otherwise False.  
      @rtype: C{bool}  
   
      """  
      return self.__COEFFICIENTS_OF_GENEARL_PDE.has_key(name)  
   
    def createCoefficientOfGeneralPDE(self,name):  
      """  
      returns a new instance of a coefficient for coefficient name of the general PDE  
   
      @param name: name of the coefficient requested.  
      @type name: C{string}  
      @return: a coefficient name initialized to 0.  
      @rtype: L{Data<escript.Data>}  
      @raise IllegalCoefficient: if name is not one of coefficients  
                   M{A}, M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact}, M{r} or M{q}.  
      """  
      if self.hasCoefficientOfGeneralPDE(name):  
         return escript.Data(0,self.getShapeOfCoefficientOfGeneralPDE(name),self.getFunctionSpaceForCoefficientOfGeneralPDE(name))  
      else:  
         raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name  
   
    def getFunctionSpaceForCoefficientOfGeneralPDE(self,name):  
      """  
      return the L{FunctionSpace<escript.FunctionSpace>} to be used for coefficient name of the general PDE  
   
      @param name: name of the coefficient enquired.  
      @type name: C{string}  
      @return: the function space to be used for coefficient name  
      @rtype: L{FunctionSpace<escript.FunctionSpace>}  
      @raise IllegalCoefficient: if name is not one of coefficients  
                   M{A}, M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact}, M{r} or M{q}.  
      """  
      if self.hasCoefficientOfGeneralPDE(name):  
         return self.__COEFFICIENTS_OF_GENEARL_PDE[name].getFunctionSpace(self.getDomain())  
      else:  
         raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name  
   
    def getShapeOfCoefficientOfGeneralPDE(self,name):  
      """  
      return the shape of the coefficient name of the general PDE  
   
      @param name: name of the coefficient enquired.  
      @type name: C{string}  
      @return: the shape of the coefficient name  
      @rtype: C{tuple} of C{int}  
      @raise IllegalCoefficient: if name is not one of coefficients  
                   M{A}, M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact}, M{r} or M{q}.  
      """  
      if self.hasCoefficientOfGeneralPDE(name):  
         return self.__COEFFICIENTS_OF_GENEARL_PDE[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())  
      else:  
         raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name  
1094    
    # =============================================================================  
    # functions giving access to coefficients of a particular PDE implementation:  
    # =============================================================================  
1095     def getCoefficient(self,name):     def getCoefficient(self,name):
1096       """       """
1097       returns the value of the coefficient name       Returns the value of the coefficient C{name}.
1098    
1099       @param name: name of the coefficient requested.       @param name: name of the coefficient requested
1100       @type name: C{string}       @type name: C{string}
1101       @return: the value of the coefficient name       @return: the value of the coefficient
1102       @rtype: L{Data<escript.Data>}       @rtype: L{Data<escript.Data>}
1103       @raise IllegalCoefficient: if name is not a coefficient of the PDE.       @raise IllegalCoefficient: if C{name} is not a coefficient of the PDE
1104       """       """
1105       if self.hasCoefficient(name):       if self.hasCoefficient(name):
1106           return self.COEFFICIENTS[name].getValue()           return self.__COEFFICIENTS[name].getValue()
1107       else:       else:
1108          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1109    
1110     def hasCoefficient(self,name):     def hasCoefficient(self,name):
1111       """       """
1112       return True if name is the name of a coefficient       Returns True if C{name} is the name of a coefficient.
1113    
1114       @param name: name of the coefficient enquired.       @param name: name of the coefficient enquired
1115       @type name: C{string}       @type name: C{string}
1116       @return: True if name is the name of a coefficient of the general PDE. Otherwise False.       @return: True if C{name} is the name of a coefficient of the general PDE,
1117                  False otherwise
1118       @rtype: C{bool}       @rtype: C{bool}
1119       """       """
1120       return self.COEFFICIENTS.has_key(name)       return self.__COEFFICIENTS.has_key(name)
1121    
1122     def createCoefficient(self, name):     def createCoefficient(self, name):
1123       """       """
1124       create a L{Data<escript.Data>} object corresponding to coefficient name       Creates a L{Data<escript.Data>} object corresponding to coefficient
1125         C{name}.
1126    
1127       @return: a coefficient name initialized to 0.       @return: the coefficient C{name} initialized to 0
1128       @rtype: L{Data<escript.Data>}       @rtype: L{Data<escript.Data>}
1129       @raise IllegalCoefficient: if name is not a coefficient of the PDE.       @raise IllegalCoefficient: if C{name} is not a coefficient of the PDE
1130       """       """
1131       if self.hasCoefficient(name):       if self.hasCoefficient(name):
1132          return escript.Data(0.,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))          return escript.Data(0.,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))
# Line 1367  class LinearPDE(object): Line 1135  class LinearPDE(object):
1135    
1136     def getFunctionSpaceForCoefficient(self,name):     def getFunctionSpaceForCoefficient(self,name):
1137       """       """
1138       return the L{FunctionSpace<escript.FunctionSpace>} to be used for coefficient name       Returns the L{FunctionSpace<escript.FunctionSpace>} to be used for
1139         coefficient C{name}.
1140    
1141       @param name: name of the coefficient enquired.       @param name: name of the coefficient enquired
1142       @type name: C{string}       @type name: C{string}
1143       @return: the function space to be used for coefficient name       @return: the function space to be used for coefficient C{name}
1144       @rtype: L{FunctionSpace<escript.FunctionSpace>}       @rtype: L{FunctionSpace<escript.FunctionSpace>}
1145       @raise IllegalCoefficient: if name is not a coefficient of the PDE.       @raise IllegalCoefficient: if C{name} is not a coefficient of the PDE
1146       """       """
1147       if self.hasCoefficient(name):       if self.hasCoefficient(name):
1148          return self.COEFFICIENTS[name].getFunctionSpace(self.getDomain())          return self.__COEFFICIENTS[name].getFunctionSpace(self.getDomain())
1149       else:       else:
1150          raise ValueError,"unknown coefficient %s requested"%name          raise ValueError,"unknown coefficient %s requested"%name
1151    
1152     def getShapeOfCoefficient(self,name):     def getShapeOfCoefficient(self,name):
1153       """       """
1154       return the shape of the coefficient name       Returns the shape of the coefficient C{name}.
1155    
1156       @param name: name of the coefficient enquired.       @param name: name of the coefficient enquired
1157       @type name: C{string}       @type name: C{string}
1158       @return: the shape of the coefficient name       @return: the shape of the coefficient C{name}
1159       @rtype: C{tuple} of C{int}       @rtype: C{tuple} of C{int}
1160       @raise IllegalCoefficient: if name is not a coefficient of the PDE.       @raise IllegalCoefficient: if C{name} is not a coefficient of the PDE
1161       """       """
1162       if self.hasCoefficient(name):       if self.hasCoefficient(name):
1163          return self.COEFFICIENTS[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())          return self.__COEFFICIENTS[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())
1164       else:       else:
1165          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1166    
1167     def resetCoefficients(self):     def resetAllCoefficients(self):
1168       """       """
1169       resets all coefficients to there default values.       Resets all coefficients to their default values.
1170       """       """
1171       for i in self.COEFFICIENTS.iterkeys():       for i in self.__COEFFICIENTS.iterkeys():
1172           self.COEFFICIENTS[i].resetValue()           self.__COEFFICIENTS[i].resetValue()
1173    
1174     def alteredCoefficient(self,name):     def alteredCoefficient(self,name):
1175       """       """
1176       announce that coefficient name has been changed       Announces that coefficient C{name} has been changed.
1177    
1178       @param name: name of the coefficient enquired.       @param name: name of the coefficient affected
1179       @type name: C{string}       @type name: C{string}
1180       @raise IllegalCoefficient: if name is not a coefficient of the PDE.       @raise IllegalCoefficient: if C{name} is not a coefficient of the PDE
1181       @note: if name is q or r, the method will not trigger a rebuilt of the system as constraints are applied to the solved system.       @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):       if self.hasCoefficient(name):
1185          self.trace("Coefficient %s has been altered."%name)          self.trace("Coefficient %s has been altered."%name)
1186          if not ((name=="q" or name=="r") and self.isUsingLumping()):          if not ((name=="q" or name=="r") and self.isUsingLumping()):
1187             if self.COEFFICIENTS[name].isAlteringOperator(): self.__invalidateOperator()             if self.__COEFFICIENTS[name].isAlteringOperator(): self.invalidateOperator()
1188             if self.COEFFICIENTS[name].isAlteringRightHandSide(): self.__invalidateRightHandSide()             if self.__COEFFICIENTS[name].isAlteringRightHandSide(): self.invalidateRightHandSide()
1189       else:       else:
1190          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1191    
1192     def copyConstraint(self,u):     def validSolution(self):
1193        """         """
1194        copies the constraint into u and returns u.         Marks the solution as valid.
1195           """
1196           self.__is_solution_valid=True
1197    
1198        @param u: a function of rank 0 is a single PDE is solved and of shape (numSolution,) for a system of PDEs     def invalidateSolution(self):
1199        @type u: L{Data<escript.Data>}         """
1200        @return: the input u modified by the constraints.         Indicates the PDE has to be resolved if the solution is requested.
1201        @rtype: L{Data<escript.Data>}         """
1202        @warning: u is altered if it has the appropriate L{FunctionSpace<escript.FunctionSpace>}         self.trace("System will be resolved.")
1203        """         self.__is_solution_valid=False
1204        q=self.getCoefficientOfGeneralPDE("q")  
1205        r=self.getCoefficientOfGeneralPDE("r")     def isSolutionValid(self):
1206        if not q.isEmpty():         """
1207           if u.isEmpty(): u=escript.Data(0.,q.getShape(),q.getFunctionSpace())         Returns True if the solution is still valid.
1208           if r.isEmpty():         """
1209               r=escript.Data(0,u.getShape(),u.getFunctionSpace())         return self.__is_solution_valid
1210           else:  
1211               r=escript.Data(r,u.getFunctionSpace())     def validOperator(self):
1212           u.copyWithMask(r,escript.Data(q,u.getFunctionSpace()))         """
1213        return u         Marks the operator as valid.
1214           """
1215           self.__is_operator_valid=True
1216    
1217       def invalidateOperator(self):
1218           """
1219           Indicates the operator has to be rebuilt next time it is used.
1220           """
1221           self.trace("Operator will be rebuilt.")
1222           self.invalidateSolution()
1223           self.__is_operator_valid=False
1224    
1225       def isOperatorValid(self):
1226           """
1227           Returns True if the operator is still valid.
1228           """
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       def invalidateRightHandSide(self):
1238           """
1239           Indicates the right hand side has to be rebuilt next time it is used.
1240           """
1241           if self.isRightHandSideValid(): self.trace("Right hand side has to be rebuilt.")
1242           self.invalidateSolution()
1243           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 invalidateSystem(self):
1252           """
1253           Announces that everything has to be rebuilt.
1254           """
1255           if self.isRightHandSideValid(): self.trace("System has to be rebuilt.")
1256           self.invalidateSolution()
1257           self.invalidateOperator()
1258           self.invalidateRightHandSide()
1259    
1260       def isSystemValid(self):
1261           """
1262           Returns True if the system (including solution) is still vaild.
1263           """
1264           return self.isSolutionValid() and self.isOperatorValid() and self.isRightHandSideValid()
1265    
1266       def initializeSystem(self):
1267           """
1268           Resets the system clearing the operator, right hand side and solution.
1269           """
1270           self.trace("New System has been created.")
1271           self.__operator=escript.Operator()
1272           self.__righthandside=escript.Data()
1273           self.__solution=escript.Data()
1274           self.invalidateSystem()
1275    
1276       def getOperator(self):
1277         """
1278         Returns the operator of the linear problem.
1279    
1280         @return: the operator of the problem
1281         """
1282         return self.getSystem()[0]
1283    
1284       def getRightHandSide(self):
1285         """
1286         Returns the right hand side of the linear problem.
1287    
1288         @return: the right hand side of the problem
1289         @rtype: L{Data<escript.Data>}
1290         """
1291         return self.getSystem()[1]
1292    
1293       def createRightHandSide(self):
1294           """
1295           Returns an instance of a new right hand side.
1296           """
1297           self.trace("New right hand side is allocated.")
1298           if self.getNumEquations()>1:
1299               return escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)
1300           else:
1301               return escript.Data(0.,(),self.getFunctionSpaceForEquation(),True)
1302    
1303       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           Sets the solution to zero.
1316           """
1317           if self.__solution.isEmpty():
1318               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():
1361                   self.__operator=self.createSolution()
1362               else:
1363                   self.__operator=self.createOperator()
1364           else:
1365               if self.isUsingLumping():
1366                   self.__operator.setToZero()
1367               else:
1368                   self.__operator.resetValues()
1369               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
1377    
1378     def setValue(self,**coefficients):     def setValue(self,**coefficients):
1379        """        """
1380        sets new values to coefficients        Sets new values to coefficients.
1381    
1382        @param coefficients: new values assigned to coefficients        @raise IllegalCoefficient: if an unknown coefficient keyword is used
       @keyword A: value for coefficient A.  
       @type A: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.  
       @keyword B: value for coefficient B  
       @type B: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.  
       @keyword C: value for coefficient C  
       @type C: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.  
       @keyword D: value for coefficient D  
       @type D: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.  
       @keyword X: value for coefficient X  
       @type X: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.  
       @keyword Y: value for coefficient Y  
       @type Y: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.  
       @keyword d: value for coefficient d  
       @type d: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.  
       @keyword y: value for coefficient y  
       @type y: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.  
       @keyword d_contact: value for coefficient d_contact  
       @type d_contact: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnContactOne<escript.FunctionOnContactOne>}.  
                        or  L{FunctionOnContactZero<escript.FunctionOnContactZero>}.  
       @keyword y_contact: value for coefficient y_contact  
       @type y_contact: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnContactOne<escript.FunctionOnContactOne>}.  
                        or  L{FunctionOnContactZero<escript.FunctionOnContactZero>}.  
       @keyword r: values prescribed to the solution at the locations of constraints  
       @type r: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}  
                depending of reduced order is used for the solution.  
       @keyword q: mask for location of constraints  
       @type q: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}  
                depending of reduced order is used for the representation of the equation.  
       @raise IllegalCoefficient: if an unknown coefficient keyword is used.  
1383        """        """
1384        # check if the coefficients are  legal:        # check if the coefficients are  legal:
1385        for i in coefficients.iterkeys():        for i in coefficients.iterkeys():
# Line 1489  class LinearPDE(object): Line 1396  class LinearPDE(object):
1396                  s=numarray.array(d).shape                  s=numarray.array(d).shape
1397              if s!=None:              if s!=None:
1398                  # get number of equations and number of unknowns:                  # get number of equations and number of unknowns:
1399                  res=self.COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(self.getDomain(),s)                  res=self.__COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(self.getDomain(),s)
1400                  if res==None:                  if res==None:
1401                      raise IllegalCoefficientValue,"Illegal shape %s of coefficient %s"%(s,i)                      raise IllegalCoefficientValue,"Illegal shape %s of coefficient %s"%(s,i)
1402                  else:                  else:
1403                      if self.__numEquations==None: self.__numEquations=res[0]                      if self.__numEquations==None: self.__numEquations=res[0]
1404                      if self.__numSolutions==None: self.__numSolutions=res[1]                      if self.__numSolutions==None: self.__numSolutions=res[1]
1405        if self.__numEquations==None: raise UndefinedPDEError,"unidententified number of equations"        if self.__numEquations==None: raise UndefinedPDEError,"unidentified number of equations"
1406        if self.__numSolutions==None: raise UndefinedPDEError,"unidententified number of solutions"        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:        # now we check the shape of the coefficient if numEquations and numSolutions are set:
1408        for i,d in coefficients.iteritems():        for i,d in coefficients.iteritems():
1409          try:          try:
1410             self.COEFFICIENTS[i].setValue(self.getDomain(),self.getNumEquations(),self.getNumSolutions(),self.reduceEquationOrder(),self.reduceSolutionOrder(),d)             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:          except IllegalCoefficientValue,m:
1430             raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))             raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
         self.alteredCoefficient(i)  
   
1431        self.__altered_coefficients=True        self.__altered_coefficients=True
1432        # check if the systrem is inhomogeneous:  
1433        if len(coefficients)>0 and not self.isUsingLumping():     # ==========================================================================
1434           q=self.getCoefficientOfGeneralPDE("q")     # methods that are typically overwritten when implementing a particular
1435           r=self.getCoefficientOfGeneralPDE("r")     # linear problem
1436           homogeneous_constraint=True     # ==========================================================================
1437           if not q.isEmpty() and not r.isEmpty():     def getRequiredSystemType(self):
1438               if util.Lsup(q*r)>=1.e-13*util.Lsup(r):        """
1439                 self.trace("Inhomogeneous constraint detected.")        Returns the system type which needs to be used by the current set up.
1440                 self.__invalidateSystem()  
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           Returns an instance of a new operator.
1449    
1450           @note: This method is overwritten when implementing a particular
1451                  linear problem.
1452           """
1453           return escript.Operator()
1454    
1455       def checkSymmetry(self,verbose=True):
1456          """
1457          Tests the PDE for symmetry.
1458    
1459          @param verbose: if set to True or not present a report on coefficients
1460                          which break the symmetry is printed
1461          @type verbose: C{bool}
1462          @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):     def getSystem(self):
1483         """         """
1484         return the operator and right hand side of the PDE         Returns the operator and right hand side of the PDE.
1485    
1486         @return: the discrete version of the PDE         @return: the discrete version of the PDE
1487         @rtype: C{tuple} of L{Operator,<escript.Operator>} and L{Data<escript.Data>}.         @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         if not self.__operator_is_Valid or not self.__righthandside_isValid:         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():            if self.isUsingLumping():
1742                if not self.__operator_is_Valid:               self.setSolution(f*1/mat)
1743                   if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():            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):
1755           """
1756           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.isOperatorValid() or not self.isRightHandSideValid():
1763              if self.isUsingLumping():
1764                  if not self.isOperatorValid():
1765                     if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():
1766                        raise TypeError,"Lumped matrix requires same order for equations and unknowns"                        raise TypeError,"Lumped matrix requires same order for equations and unknowns"
1767                   if not self.getCoefficientOfGeneralPDE("A").isEmpty():                   if not self.getCoefficient("A").isEmpty():
                       raise ValueError,"coefficient A in lumped matrix may not be present."  
                  if not self.getCoefficientOfGeneralPDE("B").isEmpty():  
1768                        raise ValueError,"coefficient A in lumped matrix may not be present."                        raise ValueError,"coefficient A in lumped matrix may not be present."
1769                   if not self.getCoefficientOfGeneralPDE("C").isEmpty():                   if not self.getCoefficient("B").isEmpty():
1770                        raise ValueError,"coefficient A in lumped matrix may not be present."                        raise ValueError,"coefficient B in lumped matrix may not be present."
1771                   D=self.getCoefficientOfGeneralPDE("D")                   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():                   if not D.isEmpty():
1788                       if self.getNumSolutions()>1:                       if self.getNumSolutions()>1:
1789                          D_times_e=util.matrixmult(D,numarray.ones((self.getNumSolutions(),)))                          D_times_e=util.matrix_mult(D,numarray.ones((self.getNumSolutions(),)))
1790                       else:                       else:
1791                          D_times_e=D                          D_times_e=D
1792                   else:                   else:
1793                      D_times_e=escript.Data()                      D_times_e=escript.Data()
                  d=self.getCoefficientOfGeneralPDE("d")  
1794                   if not d.isEmpty():                   if not d.isEmpty():
1795                       if self.getNumSolutions()>1:                       if self.getNumSolutions()>1:
1796                          d_times_e=util.matrixmult(d,numarray.ones((self.getNumSolutions(),)))                          d_times_e=util.matrix_mult(d,numarray.ones((self.getNumSolutions(),)))
1797                       else:                       else:
1798                          d_times_e=d                          d_times_e=d
1799                   else:                   else:
1800                      d_times_e=escript.Data()                      d_times_e=escript.Data()
1801                   d_contact=self.getCoefficientOfGeneralPDE("d_contact")  
1802                   if not d_contact.isEmpty():                   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:                       if self.getNumSolutions()>1:
1811                          d_contact_times_e=util.matrixmult(d_contact,numarray.ones((self.getNumSolutions(),)))                          d_reduced_times_e=util.matrix_mult(d_reduced,numarray.ones((self.getNumSolutions(),)))
1812                       else:                       else:
1813                          d_contact_times_e=d_contact                          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:                   else:
1823                      d_contact_times_e=escript.Data()                      self.getDomain().addPDEToRHS(operator, \
1824                                                         escript.Data(), \
1825                   self.__operator=self.__getNewRightHandSide()                                                   D_times_e, \
1826                   self.getDomain().addPDEToRHS(self.__operator, \                                                   d_times_e,\
1827                                                escript.Data(), \                                                   escript.Data())
1828                                                D_times_e, \                      self.getDomain().addPDEToRHS(operator, \
1829                                                d_times_e,\                                                   escript.Data(), \
1830                                                d_contact_times_e)                                                   D_reduced_times_e, \
1831                   self.__operator=1./self.__operator                                                   d_reduced_times_e,\
1832                                                     escript.Data())
1833                   self.trace("New lumped operator has been built.")                   self.trace("New lumped operator has been built.")
1834                   self.__operator_is_Valid=True                if not self.isRightHandSideValid():
1835                if not self.__righthandside_isValid:                   self.resetRightHandSide()
1836                   self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \                   righthandside=self.getCurrentRightHandSide()
1837                                 self.getCoefficientOfGeneralPDE("X"), \                   self.getDomain().addPDEToRHS(righthandside, \
1838                                 self.getCoefficientOfGeneralPDE("Y"),\                                 self.getCoefficient("X"), \
1839                                 self.getCoefficientOfGeneralPDE("y"),\                                 self.getCoefficient("Y"),\
1840                                 self.getCoefficientOfGeneralPDE("y_contact"))                                 self.getCoefficient("y"),\
1841                   self.trace("New right hand side as been built.")                                 self.getCoefficient("y_contact"))
1842                   self.__righthandside_isValid=True                   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               if not self.__operator_is_Valid and not self.__righthandside_isValid:               if not self.isOperatorValid() and not self.isRightHandSideValid():
1853                   self.getDomain().addPDEToSystem(self.__makeFreshOperator(),self.__makeFreshRightHandSide(), \                   self.resetRightHandSide()
1854                                 self.getCoefficientOfGeneralPDE("A"), \                   righthandside=self.getCurrentRightHandSide()
1855                                 self.getCoefficientOfGeneralPDE("B"), \                   self.resetOperator()
1856                                 self.getCoefficientOfGeneralPDE("C"), \                   operator=self.getCurrentOperator()
1857                                 self.getCoefficientOfGeneralPDE("D"), \                   self.getDomain().addPDEToSystem(operator,righthandside, \
1858                                 self.getCoefficientOfGeneralPDE("X"), \                                 self.getCoefficient("A"), \
1859                                 self.getCoefficientOfGeneralPDE("Y"), \                                 self.getCoefficient("B"), \
1860                                 self.getCoefficientOfGeneralPDE("d"), \                                 self.getCoefficient("C"), \
1861                                 self.getCoefficientOfGeneralPDE("y"), \                                 self.getCoefficient("D"), \
1862                                 self.getCoefficientOfGeneralPDE("d_contact"), \                                 self.getCoefficient("X"), \
1863                                 self.getCoefficientOfGeneralPDE("y_contact"))                                 self.getCoefficient("Y"), \
1864                   self.__applyConstraint()                                 self.getCoefficient("d"), \
1865                   self.__righthandside=self.copyConstraint(self.__righthandside)                                 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.")                   self.trace("New system has been built.")
1881                   self.__operator_is_Valid=True                   self.validOperator()
1882                   self.__righthandside_isValid=True                   self.validRightHandSide()
1883               elif not self.__righthandside_isValid:               elif not self.isRightHandSideValid():
1884                   self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \                   self.resetRightHandSide()
1885                                 self.getCoefficientOfGeneralPDE("X"), \                   righthandside=self.getCurrentRightHandSide()
1886                                 self.getCoefficientOfGeneralPDE("Y"),\                   self.getDomain().addPDEToRHS(righthandside,
1887                                 self.getCoefficientOfGeneralPDE("y"),\                                 self.getCoefficient("X"), \
1888                                 self.getCoefficientOfGeneralPDE("y_contact"))                                 self.getCoefficient("Y"),\
1889                   self.__righthandside=self.copyConstraint(self.__righthandside)                                 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.")                   self.trace("New right hand side has been built.")
1898                   self.__righthandside_isValid=True                   self.validRightHandSide()
1899               elif not self.__operator_is_Valid:               elif not self.isOperatorValid():
1900                   self.getDomain().addPDEToSystem(self.__makeFreshOperator(),escript.Data(), \                   self.resetOperator()
1901                              self.getCoefficientOfGeneralPDE("A"), \                   operator=self.getCurrentOperator()
1902                              self.getCoefficientOfGeneralPDE("B"), \                   self.getDomain().addPDEToSystem(operator,escript.Data(), \
1903                              self.getCoefficientOfGeneralPDE("C"), \                              self.getCoefficient("A"), \
1904                              self.getCoefficientOfGeneralPDE("D"), \                              self.getCoefficient("B"), \
1905                                self.getCoefficient("C"), \
1906                                self.getCoefficient("D"), \
1907                              escript.Data(), \                              escript.Data(), \
1908                              escript.Data(), \                              escript.Data(), \
1909                              self.getCoefficientOfGeneralPDE("d"), \                              self.getCoefficient("d"), \
1910                              escript.Data(),\                              escript.Data(),\
1911                              self.getCoefficientOfGeneralPDE("d_contact"), \                              self.getCoefficient("d_contact"), \
1912                              escript.Data())                              escript.Data())
1913                   self.__applyConstraint()                   self.getDomain().addPDEToSystem(operator,escript.Data(), \
1914                                self.getCoefficient("A_reduced"), \
1915                                self.getCoefficient("B_reduced"), \
1916                                self.getCoefficient("C_reduced"), \
1917                                self.getCoefficient("D_reduced"), \
1918                                escript.Data(), \
1919                                escript.Data(), \
1920                                self.getCoefficient("d_reduced"), \
1921                                escript.Data(),\
1922                                self.getCoefficient("d_contact_reduced"), \
1923                                escript.Data())
1924                     self.insertConstraint(rhs_only=False)
1925                   self.trace("New operator has been built.")                   self.trace("New operator has been built.")
1926                   self.__operator_is_Valid=True                   self.validOperator()
1927         return (self.__operator,self.__righthandside)         return (self.getCurrentOperator(), self.getCurrentRightHandSide())
1928    
1929       def insertConstraint(self, rhs_only=False):
1930          """
1931          Applies the constraints defined by q and r to the PDE.
1932    
1933          @param rhs_only: if True only the right hand side is altered by the
1934                           constraint
1935          @type rhs_only: C{bool}
1936          """
1937          q=self.getCoefficient("q")
1938          r=self.getCoefficient("r")
1939          righthandside=self.getCurrentRightHandSide()
1940          operator=self.getCurrentOperator()
1941    
1942          if not q.isEmpty():
1943             if r.isEmpty():
1944                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       def setValue(self,**coefficients):
1960          """
1961          Sets new values to coefficients.
1962    
1963          @param coefficients: new values assigned to coefficients
1964          @keyword A: value for coefficient C{A}
1965          @type A: any type that can be cast to a L{Data<escript.Data>} object on
1966                   L{Function<escript.Function>}
1967          @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         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         Returns the flux M{J} for a given M{u}.
2067    
2068         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):  class Poisson(LinearPDE):
2090     """     """
2091     Class to define a Poisson equation problem, which is genear L{LinearPDE} of the form     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}     M{-grad(grad(u)[j])[j] = f}
2095    
2096     with natural boundary conditons     with natural boundary conditions
2097    
2098     M{n[j]*grad(u)[j] = 0 }     M{n[j]*grad(u)[j] = 0 }
2099    
# Line 1639  class Poisson(LinearPDE): Line 2105  class Poisson(LinearPDE):
2105    
2106     def __init__(self,domain,debug=False):     def __init__(self,domain,debug=False):
2107       """       """
2108       initializes a new Poisson equation       Initializes a new Poisson equation.
2109    
2110       @param domain: domain of the PDE       @param domain: domain of the PDE
2111       @type domain: L{Domain<escript.Domain>}       @type domain: L{Domain<escript.Domain>}
2112       @param debug: if True debug informations are printed.       @param debug: if True debug information is printed
2113    
2114       """       """
2115       super(Poisson, self).__init__(domain,1,1,debug)       super(Poisson, self).__init__(domain,1,1,debug)
2116       self.COEFFICIENTS={"f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),       self.introduceCoefficients(
2117                            "q": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}                          f=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2118                            f_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE))
2119       self.setSymmetryOn()       self.setSymmetryOn()
2120    
2121     def setValue(self,**coefficients):     def setValue(self,**coefficients):
2122       """       """
2123       sets new values to coefficients       Sets new values to coefficients.
2124    
2125       @param coefficients: new values assigned to coefficients       @param coefficients: new values assigned to coefficients
2126       @keyword f: value for right hand side M{f}       @keyword f: value for right hand side M{f}
2127       @type f: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.       @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       @keyword q: mask for location of constraints
2130       @type q: any type that can be casted to rank zeo L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}       @type q: any type that can be cast to a rank zero L{Data<escript.Data>}
2131                 depending of reduced order is used for the representation of the equation.                object on L{Solution<escript.Solution>} or
2132       @raise IllegalCoefficient: if an unknown coefficient keyword is used.                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       super(Poisson, self).setValue(**coefficients)       super(Poisson, self).setValue(**coefficients)
2137    
2138     def getCoefficientOfGeneralPDE(self,name):  
2139       def getCoefficient(self,name):
2140       """       """
2141       return the value of the coefficient name of the general PDE       Returns the value of the coefficient C{name} of the general PDE.
2142       @param name: name of the coefficient requested.  
2143         @param name: name of the coefficient requested
2144       @type name: C{string}       @type name: C{string}
2145       @return: the value of the coefficient  name       @return: the value of the coefficient C{name}
2146       @rtype: L{Data<escript.Data>}       @rtype: L{Data<escript.Data>}
2147       @raise IllegalCoefficient: if name is not one of coefficients       @raise IllegalCoefficient: invalid coefficient name
2148                    M{A}, M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact}, M{r} or M{q}.       @note: This method is called by the assembling routine to map the Poisson
2149       @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.              equation onto the general PDE.
2150       """       """
2151       if name == "A" :       if name == "A" :
2152           return escript.Data(util.kronecker(self.getDim()),escript.Function(self.getDomain()))           return escript.Data(util.kronecker(self.getDim()),escript.Function(self.getDomain()))
      elif name == "B" :  
          return escript.Data()  
      elif name == "C" :  
          return escript.Data()  
      elif name == "D" :  
          return escript.Data()  
      elif name == "X" :  
          return escript.Data()  
2153       elif name == "Y" :       elif name == "Y" :
2154           return self.getCoefficient("f")           return self.getCoefficient("f")
2155       elif name == "d" :       elif name == "Y_reduced" :
2156           return escript.Data()           return self.getCoefficient("f_reduced")
      elif name == "y" :  
          return escript.Data()  
      elif name == "d_contact" :  
          return escript.Data()  
      elif name == "y_contact" :  
          return escript.Data()  
      elif name == "r" :  
          return escript.Data()  
      elif name == "q" :  
          return self.getCoefficient("q")  
2157       else:       else:
2158          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name           return super(Poisson, self).getCoefficient(name)
2159    
2160  class Helmholtz(LinearPDE):  class Helmholtz(LinearPDE):
2161     """     """
2162     Class to define a Helmhotz equation problem, which is genear L{LinearPDE} of the form     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}     M{S{omega}*u - grad(k*grad(u)[j])[j] = f}
2166    
2167     with natural boundary conditons     with natural boundary conditions
2168    
2169     M{k*n[j]*grad(u)[j] = g- S{alpha}u }     M{k*n[j]*grad(u)[j] = g- S{alpha}u }
2170    
# Line 1721  class Helmholtz(LinearPDE): Line 2176  class Helmholtz(LinearPDE):
2176    
2177     def __init__(self,domain,debug=False):     def __init__(self,domain,debug=False):
2178       """       """
2179       initializes a new Poisson equation       Initializes a new Helmholtz equation.
2180    
2181       @param domain: domain of the PDE       @param domain: domain of the PDE
2182       @type domain: L{Domain<escript.Domain>}       @type domain: L{Domain<escript.Domain>}
2183       @param debug: if True debug informations are printed.       @param debug: if True debug information is printed
2184    
2185       """       """
2186       super(Helmholtz, self).__init__(domain,1,1,debug)       super(Helmholtz, self).__init__(domain,1,1,debug)
2187       self.COEFFICIENTS={"omega": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),       self.introduceCoefficients(
2188                          "k": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),                          omega=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.OPERATOR),
2189                          "f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),                          k=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.OPERATOR),
2190                          "alpha": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),                          f=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2191                          "g": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),                          f_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2192                          "r": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH),                          alpha=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.OPERATOR),
2193                          "q": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}                          g=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2194                            g_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE))
2195       self.setSymmetryOn()       self.setSymmetryOn()
2196    
2197     def setValue(self,**coefficients):     def setValue(self,**coefficients):
2198       """       """
2199       sets new values to coefficients       Sets new values to coefficients.
2200    
2201       @param coefficients: new values assigned to coefficients       @param coefficients: new values assigned to coefficients
2202       @keyword omega: value for coefficient M{S{omega}}       @keyword omega: value for coefficient M{S{omega}}
2203       @type omega: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.       @type omega: any type that can be cast to a L{Scalar<escript.Scalar>}
2204       @keyword k: value for coefficeint M{k}                    object on L{Function<escript.Function>}
2205       @type k: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.       @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}       @keyword f: value for right hand side M{f}
2209       @type f: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.       @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}}       @keyword alpha: value for right hand side M{S{alpha}}
2212       @type alpha: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.       @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}       @keyword g: value for right hand side M{g}
2215       @type g: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.       @type g: any type that can be cast to a L{Scalar<escript.Scalar>} object
2216       @keyword r: prescribed values M{r} for the solution in constraints.                on L{FunctionOnBoundary<escript.FunctionOnBoundary>}
2217       @type r: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}       @keyword r: prescribed values M{r} for the solution in constraints
2218                 depending of reduced order is used for the representation of the equation.       @type r: any type that can be cast to a L{Scalar<escript.Scalar>} object
2219       @keyword q: mask for location of constraints                on L{Solution<escript.Solution>} or
2220       @type q: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}                L{ReducedSolution<escript.ReducedSolution>} depending on whether
2221                 depending of reduced order is used for the representation of the equation.                reduced order is used for the representation of the equation
2222       @raise IllegalCoefficient: if an unknown coefficient keyword is used.       @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)       super(Helmholtz, self).setValue(**coefficients)
2230    
2231     def getCoefficientOfGeneralPDE(self,name):     def getCoefficient(self,name):
2232       """       """
2233       return the value of the coefficient name of the general PDE       Returns the value of the coefficient C{name} of the general PDE.
2234    
2235       @param name: name of the coefficient requested.       @param name: name of the coefficient requested
2236       @type name: C{string}       @type name: C{string}
2237       @return: the value of the coefficient  name       @return: the value of the coefficient C{name}
2238       @rtype: L{Data<escript.Data>}       @rtype: L{Data<escript.Data>}
2239       @raise IllegalCoefficient: if name is not one of coefficients       @raise IllegalCoefficient: invalid name
                   "A", M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact}, M{r} or M{q}.  
      @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.  
2240       """       """
2241       if name == "A" :       if name == "A" :
2242           return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))*self.getCoefficient("k")           if self.getCoefficient("k").isEmpty():
2243       elif name == "B" :                return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))
2244           return escript.Data()           else:
2245       elif name == "C" :                return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))*self.getCoefficient("k")
          return escript.Data()  
2246       elif name == "D" :       elif name == "D" :
2247           return self.getCoefficient("omega")           return self.getCoefficient("omega")
      elif name == "X" :  
          return escript.Data()  
2248       elif name == "Y" :       elif name == "Y" :
2249           return self.getCoefficient("f")           return self.getCoefficient("f")
2250       elif name == "d" :       elif name == "d" :
2251           return self.getCoefficient("alpha")           return self.getCoefficient("alpha")
2252       elif name == "y" :       elif name == "y" :
2253           return self.getCoefficient("g")           return self.getCoefficient("g")
2254       elif name == "d_contact" :       elif name == "Y_reduced" :
2255           return escript.Data()           return self.getCoefficient("f_reduced")
2256       elif name == "y_contact" :       elif name == "y_reduced" :
2257           return escript.Data()          return self.getCoefficient("g_reduced")
      elif name == "r" :  
          return self.getCoefficient("r")  
      elif name == "q" :  
          return self.getCoefficient("q")  
2258       else:       else:
2259          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name          return super(Helmholtz, self).getCoefficient(name)
2260    
2261  class LameEquation(LinearPDE):  class LameEquation(LinearPDE):
2262     """     """
2263     Class to define a Lame equation problem:     Class to define a Lame equation problem. This problem is defined as:
2264    
2265     M{-grad(S{mu}*(grad(u[i])[j]+grad(u[j])[i]))[j] - grad(S{lambda}*grad(u[j])[i])[j] = F_i -grad(S{sigma}[i,j])[j] }     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] }
2266    
2267     with natural boundary conditons:     with natural boundary conditions:
2268    
2269     M{n[j]*(S{mu}*(grad(u[i])[j]+grad(u[j])[i]) - S{lambda}*grad(u[j])[i]) = f_i -n[j]*S{sigma}[i,j] }     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] }
2270    
2271     and constraints:     and constraints:
2272    
# Line 1819  class LameEquation(LinearPDE): Line 2275  class LameEquation(LinearPDE):
2275     """     """
2276    
2277     def __init__(self,domain,debug=False):     def __init__(self,domain,debug=False):
2278          """
2279          Initializes a new Lame equation.
2280    
2281          @param domain: domain of the PDE
2282          @type domain: L{Domain<escript.Domain>}
2283          @param debug: if True debug information is printed
2284    
2285          """
2286        super(LameEquation, self).__init__(domain,\        super(LameEquation, self).__init__(domain,\
2287                                           domain.getDim(),domain.getDim(),debug)                                           domain.getDim(),domain.getDim(),debug)
2288        self.COEFFICIENTS={ "lame_lambda"  : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),        self.introduceCoefficients(lame_lambda=PDECoef(PDECoef.INTERIOR,(),PDECoef.OPERATOR),
2289                            "lame_mu"      : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),                                   lame_mu=PDECoef(PDECoef.INTERIOR,(),PDECoef.OPERATOR),
2290                            "F"            : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),                                   F=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2291                            "sigma"        : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM),PDECoefficient.RIGHTHANDSIDE),                                   sigma=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
2292                            "f"            : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),                                   f=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE))
                           "r"            : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH),  
                           "q"            : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}  
2293        self.setSymmetryOn()        self.setSymmetryOn()
2294    
2295     def setValues(self,**coefficients):     def setValues(self,**coefficients):
2296       """       """
2297       sets new values to coefficients       Sets new values to coefficients.
2298    
2299       @param coefficients: new values assigned to coefficients       @param coefficients: new values assigned to coefficients
2300       @keyword lame_mu: value for coefficient M{S{mu}}       @keyword lame_mu: value for coefficient M{S{mu}}
2301       @type lame_mu: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.       @type lame_mu: any type that can be cast to a L{Scalar<escript.Scalar>}
2302                        object on L{Function<escript.Function>}
2303       @keyword lame_lambda: value for coefficient M{S{lambda}}       @keyword lame_lambda: value for coefficient M{S{lambda}}
2304       @type lame_lambda: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.       @type lame_lambda: any type that can be cast to a L{Scalar<escript.Scalar>}
2305                            object on L{Function<escript.Function>}
2306       @keyword F: value for internal force M{F}       @keyword F: value for internal force M{F}
2307       @type F: any type that can be casted to L{Vector<escript.Vector>} object on L{Function<escript.Function>}       @type F: any type that can be cast to a L{Vector<escript.Vector>} object
2308                  on L{Function<escript.Function>}
2309       @keyword sigma: value for initial stress M{S{sigma}}       @keyword sigma: value for initial stress M{S{sigma}}
2310       @type sigma: any type that can be casted to L{Tensor<escript.Tensor>} object on L{Function<escript.Function>}       @type sigma: any type that can be cast to a L{Tensor<escript.Tensor>}
2311       @keyword f: value for extrenal force M{f}                    object on L{Function<escript.Function>}
2312       @type f: any type that can be casted to L{Vector<escript.Vector>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}       @keyword f: value for external force M{f}
2313       @keyword r: prescribed values M{r} for the solution in constraints.       @type f: any type that can be cast to a L{Vector<escript.Vector>} object
2314       @type r: any type that can be casted to L{Vector<escript.Vector>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}                on L{FunctionOnBoundary<escript.FunctionOnBoundary>}
2315                 depending of reduced order is used for the representation of the equation.       @keyword r: prescribed values M{r} for the solution in constraints
2316       @keyword q: mask for location of constraints       @type r: any type that can be cast to a L{Vector<escript.Vector>} object
2317       @type q: any type that can be casted to L{Vector<escript.Vector>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}                on L{Solution<escript.Solution>} or
2318                 depending of reduced order is used for the representation of the equation.                L{ReducedSolution<escript.ReducedSolution>} depending on whether
2319       @raise IllegalCoefficient: if an unknown coefficient keyword is used.                reduced order is used for the representation of the equation
2320         @keyword q: mask for the location of constraints
2321         @type q: any type that can be cast to a L{Vector<escript.Vector>} object
2322                  on L{Solution<escript.Solution>} or
2323                  L{ReducedSolution<escript.ReducedSolution>} depending on whether
2324                  reduced order is used for the representation of the equation
2325         @raise IllegalCoefficient: if an unknown coefficient keyword is used
2326       """       """
2327       super(LameEquation, self).setValues(**coefficients)       super(LameEquation, self).setValues(**coefficients)
2328    
2329     def getCoefficientOfGeneralPDE(self,name):     def getCoefficient(self,name):
2330       """       """
2331       return the value of the coefficient name of the general PDE       Returns the value of the coefficient C{name} of the general PDE.
2332    
2333       @param name: name of the coefficient requested.       @param name: name of the coefficient requested
2334       @type name: C{string}       @type name: C{string}
2335       @return: the value of the coefficient  name       @return: the value of the coefficient C{name}
2336       @rtype: L{Data<escript.Data>}       @rtype: L{Data<escript.Data>}
2337       @raise IllegalCoefficient: if name is not one of coefficients       @raise IllegalCoefficient: invalid coefficient name
                   "A", M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact}, M{r} or M{q}.  
      @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.  
2338       """       """
2339       if name == "A" :       if name == "A" :
2340           out =self.createCoefficientOfGeneralPDE("A")           out =self.createCoefficient("A")
2341           for i in range(self.getDim()):           for i in range(self.getDim()):
2342             for j in range(self.getDim()):             for j in range(self.getDim()):
2343               out[i,i,j,j] += self.getCoefficient("lame_lambda")               out[i,i,j,j] += self.getCoefficient("lame_lambda")
2344               out[i,j,j,i] += self.getCoefficient("lame_mu")               out[i,j,j,i] += self.getCoefficient("lame_mu")
2345               out[i,j,i,j] += self.getCoefficient("lame_mu")               out[i,j,i,j] += self.getCoefficient("lame_mu")
2346           return out           return out
      elif name == "B" :  
          return escript.Data()  
      elif name == "C" :  
          return escript.Data()  
      elif name == "D" :  
          return escript.Data()  
2347       elif name == "X" :       elif name == "X" :
2348           return self.getCoefficient("sigma")           return self.getCoefficient("sigma")
2349       elif name == "Y" :       elif name == "Y" :
2350           return self.getCoefficient("F")           return self.getCoefficient("F")
      elif name == "d" :  
          return escript.Data()  
2351       elif name == "y" :       elif name == "y" :
2352           return self.getCoefficient("f")           return self.getCoefficient("f")
      elif name == "d_contact" :  
          return escript.Data()  
      elif name == "y_contact" :  
          return escript.Data()  
      elif name == "r" :  
          return self.getCoefficient("r")  
      elif name == "q" :  
          return self.getCoefficient("q")  
2353       else:       else:
2354          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name          return super(LameEquation, self).getCoefficient(name)
2355    
2356  class AdvectivePDE(LinearPDE):  def LinearSinglePDE(domain,debug=False):
2357     """     """
2358     In cases of PDEs dominated by the advection terms M{B} and M{C} against the adevctive terms M{A}     Defines a single linear PDE.
    up-winding has been used.  The L{AdvectivePDE} class applies SUPG upwinding to the advective terms.  
2359    
2360     In the following we set     @param domain: domain of the PDE
2361       @type domain: L{Domain<escript.Domain>}
2362       @param debug: if True debug information is printed
2363       @rtype: L{LinearPDE}
2364       """
2365       return LinearPDE(domain,numEquations=1,numSolutions=1,debug=debug)
2366    
2367     M{Z[j]=C[j]-B[j]}  def LinearPDESystem(domain,debug=False):
2368       """
2369       Defines a system of linear PDEs.
2370    
2371     or     @param domain: domain of the PDEs
2372       @type domain: L{Domain<escript.Domain>}
2373       @param debug: if True debug information is printed
2374       @rtype: L{LinearPDE}
2375       """
2376       return LinearPDE(domain,numEquations=domain.getDim(),numSolutions=domain.getDim(),debug=debug)
2377    
    M{Z[i,k,l]=C[i,k,l]-B[i,l,k]}  
2378    
2379     To measure the dominance of the advective terms over the diffusive term M{A} the  class TransportPDE(LinearProblem):
2380     X{Pelclet number} M{P} is used. It is defined as     """
2381       This class is used to define a transport problem given by a general linear,
2382       time dependent, second order PDE for an unknown, non-negative,
2383       time-dependent function M{u} on a given domain defined through a
2384       L{Domain<escript.Domain>} object.
2385    
2386       For a single equation with a solution with a single component the transport
2387       problem is defined in the following form:
2388    
2389       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)}
2390    
2391       where M{u_t} denotes the time derivative of M{u} and M{grad(F)} denotes the
2392       spatial derivative of M{F}. Einstein's summation convention,  ie. summation
2393       over indexes appearing twice in a term of a sum performed, is used.
2394       The coefficients M{M}, M{A}, M{B}, M{C}, M{D}, M{X} and M{Y} have to be
2395       specified through L{Data<escript.Data>} objects in L{Function<escript.Function>}
2396       and the coefficients M{M_reduced}, M{A_reduced}, M{B_reduced}, M{C_reduced},
2397       M{D_reduced}, M{X_reduced} and M{Y_reduced} have to be specified through
2398       L{Data<escript.Data>} objects in L{ReducedFunction<escript.ReducedFunction>}.
2399       It is also allowed to use objects that can be converted into such
2400       L{Data<escript.Data>} objects. M{M} and M{M_reduced} are scalar, M{A} and
2401       M{A_reduced} are rank two, M{B}, M{C}, M{X}, M{B_reduced}, M{C_reduced} and
2402       M{X_reduced} are rank one and M{D}, M{D_reduced}, M{Y} and M{Y_reduced} are
2403       scalar.
2404    
2405     M{P=h|Z|/(2|A|)}     The following natural boundary conditions are considered:
2406    
2407     where M{|.|} denotes the L{length<util.length>} of the arument and M{h} is the local cell size     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}
    from L{getSize<escript.Domain.getSize>}. Where M{|A|==0} M{P} is M{S{infinity}}.  
2408    
2409     From the X{Pelclet number} the stabilization parameters M{S{Xi}} and M{S{Xi}} are calculated:     where M{n} is the outer normal field. Notice that the coefficients M{A},
2410       M{A_reduced}, M{B}, M{B_reduced}, M{X} and M{X_reduced} are defined in the
2411       transport problem. The coefficients M{m}, M{d} and M{y} are each a scalar in
2412       L{FunctionOnBoundary<escript.FunctionOnBoundary>} and the coefficients
2413       M{m_reduced}, M{d_reduced} and M{y_reduced} are each a scalar in
2414       L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
2415    
2416       Constraints for the solution prescribing the value of the solution at
2417       certain locations in the domain have the form
2418    
2419       M{u_t=r} where M{q>0}
2420    
2421       M{r} and M{q} are each scalar where M{q} is the characteristic function
2422       defining where the constraint is applied. The constraints override any other
2423       condition set by the transport problem or the boundary condition.
2424    
2425       The transport problem is symmetrical if
2426    
2427       M{A[i,j]=A[j,i]} and M{B[j]=C[j]} and M{A_reduced[i,j]=A_reduced[j,i]} and
2428       M{B_reduced[j]=C_reduced[j]}
2429    
2430       For a system and a solution with several components the transport problem
2431       has the form
2432    
2433       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] }
2434    
2435       M{A} and M{A_reduced} are of rank four, M{B}, M{B_reduced}, M{C} and
2436       M{C_reduced} are each of rank three, M{M}, M{M_reduced}, M{D}, M{D_reduced},
2437       M{X_reduced} and M{X} are each of rank two and M{Y} and M{Y_reduced} are of
2438       rank one. The natural boundary conditions take the form:
2439    
2440       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}
2441    
2442       The coefficient M{d} and M{m} are of rank two and M{y} is of rank one with
2443       all in L{FunctionOnBoundary<escript.FunctionOnBoundary>}. The coefficients
2444       M{d_reduced} and M{m_reduced} are of rank two and M{y_reduced} is of rank
2445       one all in L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
2446    
2447       Constraints take the form
2448    
2449       M{u[i]_t=r[i]} where M{q[i]>0}
2450    
2451       M{r} and M{q} are each rank one. Notice that at some locations not
2452       necessarily all components must have a constraint.
2453    
2454       The transport problem is symmetrical if
2455    
2456          - M{M[i,k]=M[i,k]}
2457          - M{M_reduced[i,k]=M_reduced[i,k]}
2458          - M{A[i,j,k,l]=A[k,l,i,j]}
2459          - M{A_reduced[i,j,k,l]=A_reduced[k,l,i,j]}
2460          - M{B[i,j,k]=C[k,i,j]}
2461          - M{B_reduced[i,j,k]=C_reduced[k,i,j]}
2462          - M{D[i,k]=D[i,k]}
2463          - M{D_reduced[i,k]=D_reduced[i,k]}
2464          - M{m[i,k]=m[k,i]}
2465          - M{m_reduced[i,k]=m_reduced[k,i]}
2466          - M{d[i,k]=d[k,i]}
2467          - M{d_reduced[i,k]=d_reduced[k,i]}
2468    
2469       L{TransportPDE} also supports solution discontinuities over a contact region
2470       in the domain. To specify the conditions across the discontinuity we are
2471       using the generalised flux M{J} which, in the case of a system of PDEs and
2472       several components of the solution, is defined as
2473    
2474       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]}
2475    
2476       For the case of single solution component and single PDE M{J} is defined as
2477    
2478       M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[j]+(B[i]+B_reduced[i])*u+X[i]+X_reduced[i]}
2479    
2480       In the context of discontinuities M{n} denotes the normal on the
2481       discontinuity pointing from side 0 towards side 1 calculated from
2482       L{getNormal<escript.FunctionSpace.getNormal>} of L{FunctionOnContactZero<escript.FunctionOnContactZero>}.
2483       For a system of transport problems the contact condition takes the form
2484    
2485       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]}
2486    
2487       where M{J0} and M{J1} are the fluxes on side 0 and side 1 of the
2488       discontinuity, respectively. M{jump(u)}, which is the difference of the
2489       solution at side 1 and at side 0, denotes the jump of M{u} across
2490       discontinuity along the normal calculated by L{jump<util.jump>}.
2491       The coefficient M{d_contact} is of rank two and M{y_contact} is of rank one
2492       both in L{FunctionOnContactZero<escript.FunctionOnContactZero>} or L{FunctionOnContactOne<escript.FunctionOnContactOne>}.
2493       The coefficient M{d_contact_reduced} is of rank two and M{y_contact_reduced}
2494       is of rank one both in L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.
2495       In case of a single PDE and a single component solution the contact
2496       condition takes the form
2497    
2498       M{n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)}
2499    
2500       In this case the coefficient M{d_contact} and M{y_contact} are each scalar
2501       both in L{FunctionOnContactZero<escript.FunctionOnContactZero>} or
2502       L{FunctionOnContactOne<escript.FunctionOnContactOne>} and the coefficient
2503       M{d_contact_reduced} and M{y_contact_reduced} are each scalar both in
2504       L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or
2505       L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.
2506    
2507       Typical usage::
2508    
2509           p = TransportPDE(dom)
2510           p.setValue(M=1., C=[-1.,0.])
2511           p.setInitialSolution(u=exp(-length(dom.getX()-[0.1,0.1])**2)
2512           t = 0
2513           dt = 0.1
2514           while (t < 1.):
2515               u = p.solve(dt)
2516    
2517     M{S{Xi}=S{xi}(P) h/|Z|}     """
2518       def __init__(self,domain,numEquations=None,numSolutions=None, theta=0.5,debug=False):
2519         """
2520         Initializes a transport problem.
2521    
2522     where M{S{xi}} is a suitable function of the Peclet number.       @param domain: domain of the PDE
2523         @type domain: L{Domain<escript.Domain>}
2524         @param numEquations: number of equations. If C{None} the number of
2525                              equations is extracted from the coefficients.
2526         @param numSolutions: number of solution components. If C{None} the number
2527                              of solution components is extracted from the
2528                              coefficients.
2529         @param debug: if True debug information is printed
2530         @param theta: time stepping control:
2531            - 1.0: backward Euler,
2532            - 0.0: forward Euler,
2533            - 0.5: Crank-Nicholson scheme, U{see here<http://en.wikipedia.org/wiki/Crank-Nicolson_method>}
2534    
2535         """
2536         if theta<0 or theta>1:
2537                  raise ValueError,"theta (=%s) needs to be between 0 and 1 (inclusive)."%theta
2538         super(TransportPDE, self).__init__(domain,numEquations,numSolutions,debug)
2539    
2540     In the case of a single PDE the coefficient are up-dated in the following way:       self.__theta=theta
2541           - M{A[i,j] S{<-} A[i,j] + S{Xi} * Z[j] * Z[l]}       self.setConstraintWeightingFactor()
2542           - M{B[j] S{<-} B[j] + S{Xi} * C[j] * D}       #
2543           - M{C[j] S{<-} C[j] + S{Xi} * B[j] * D}       #   the coefficients of the transport problem
2544           - M{X[j] S{<-} X[j] + S{Xi} * Z[j] * Y}       #
2545         self.introduceCoefficients(
2546           M=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2547           A=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2548           B=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2549           C=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2550           D=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2551           X=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
2552           Y=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2553           m=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2554           d=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2555           y=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2556           d_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2557           y_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2558           M_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2559           A_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2560           B_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2561           C_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2562           D_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2563           X_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
2564           Y_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2565           m_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2566           d_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2567           y_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2568           d_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2569           y_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2570           r=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.RIGHTHANDSIDE),
2571           q=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.BOTH) )
2572    
2573     Similar for the case of a systems of PDEs:     def __str__(self):
2574           - M{A[i,j,k,l] S{<-} A[i,j,k,l]+ S{delta}[p,m] * S{Xi} * Z[p,i,j] * Z[m,k,l]}       """
2575           - M{B[i,j,k] S{<-} B[i,j,k] +  S{delta}[p,m] * S{Xi} * D[p,k] * C[m,i,j]}       Returns the string representation of the transport problem.
          - M{C[i,k,l] S{<-} C[i,k,l] +  S{delta}[p,m] * S{Xi} * D[p,k] * B[m,l,i]}  
          - M{X[i,j] S{<-} X[i,j] + S{delta}[p,m] * S{Xi}  * Y[p] * Z[m,i,j]}  
2576    
2577     where M{S{delta}} is L{kronecker}.       @return: a simple representation of the transport problem
2578     Using upwinding in this form, introduces an additonal error which is proprtional to the cell size M{h}       @rtype: C{str}
2579     but with the intension to stabilize the solution.       """
2580         return "<TransportPDE %d>"%id(self)
2581    
2582     """     def getTheta(self):
    def __init__(self,domain,numEquations=None,numSolutions=None,xi=None,debug=False):  
2583        """        """
2584        creates a linear, steady, second order PDE on a L{Domain<escript.Domain>}        Returns the time stepping control parameter.
2585          @rtype: float
2586          """
2587          return self.__theta
2588    
2589        @param domain: domain of the PDE  
2590        @type domain: L{Domain<escript.Domain>}     def checkSymmetry(self,verbose=True):
2591        @param numEquations: number of equations. If numEquations==None the number of equations        """
2592                             is exracted from the PDE coefficients.        Tests the transport problem for symmetry.
2593        @param numSolutions: number of solution components. If  numSolutions==None the number of solution components  
2594                             is exracted from the PDE coefficients.        @param verbose: if set to True or not present a report on coefficients
2595        @param xi: defines a function which returns for any given Preclet number as L{Scalar<escript.Scalar>} object the                        which break the symmetry is printed.
2596                   M{S{xi}}-value used to define the stabilization parameters. If equal to None, L{ELMAN_RAMAGE} is used.        @type verbose: C{bool}
2597        @type xi: callable object which returns a L{Scalar<escript.Scalar>} object.        @return:  True if the PDE is symmetric
2598        @param debug: if True debug informations are printed.        @rtype: C{bool}
2599        """        @note: This is a very expensive operation. It should be used for
2600        super(AdvectivePDE, self).__init__(domain,\               degugging only! The symmetry flag is not altered.
2601                                           numEquations,numSolutions,debug)        """
2602        if xi==None:        out=True
2603           self.__xi=AdvectivePDE.ELMAN_RAMAGE        out=out and self.checkSymmetricTensor("M", verbose)
2604        else:        out=out and self.checkSymmetricTensor("M_reduced", verbose)
2605           self.__xi=xi        out=out and self.checkSymmetricTensor("A", verbose)
2606        self.__Xi=escript.Data()        out=out and self.checkSymmetricTensor("A_reduced", verbose)
2607          out=out and self.checkReciprocalSymmetry("B","C", verbose)
2608          out=out and self.checkReciprocalSymmetry("B_reduced","C_reduced", verbose)
2609          out=out and self.checkSymmetricTensor("D", verbose)
2610          out=out and self.checkSymmetricTensor("D_reduced", verbose)
2611          out=out and self.checkSymmetricTensor("m", verbose)
2612          out=out and self.checkSymmetricTensor("m_reduced", verbose)
2613          out=out and self.checkSymmetricTensor("d", verbose)
2614          out=out and self.checkSymmetricTensor("d_reduced", verbose)
2615          out=out and self.checkSymmetricTensor("d_contact", verbose)
2616          out=out and self.checkSymmetricTensor("d_contact_reduced", verbose)
2617          return out
2618    
2619     def setValue(self,**coefficients):     def setValue(self,**coefficients):
2620        """        """
2621        sets new values to coefficients        Sets new values to coefficients.
2622    
2623        @param coefficients: new values assigned to coefficients        @param coefficients: new values assigned to coefficients
2624        @keyword A: value for coefficient A.        @keyword M: value for coefficient C{M}
2625        @type A: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.        @type M: any type that can be cast to a L{Data<escript.Data>} object on
2626        @keyword B: value for coefficient B                 L{Function<escript.Function>}
2627        @type B: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.        @keyword M_reduced: value for coefficient C{M_reduced}
2628        @keyword C: value for coefficient C        @type M_reduced: any type that can be cast to a L{Data<escript.Data>}
2629        @type C: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.                         object on L{Function<escript.ReducedFunction>}
2630        @keyword D: value for coefficient D        @keyword A: value for coefficient C{A}
2631        @type D: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.        @type A: any type that can be cast to a L{Data<escript.Data>} object on
2632        @keyword X: value for coefficient X                 L{Function<escript.Function>}
2633        @type X: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.        @keyword A_reduced: value for coefficient C{A_reduced}
2634        @keyword Y: value for coefficient Y        @type A_reduced: any type that can be cast to a L{Data<escript.Data>}
2635        @type Y: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.                         object on L{ReducedFunction<escript.ReducedFunction>}
2636        @keyword d: value for coefficient d        @keyword B: value for coefficient C{B}
2637        @type d: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.        @type B: any type that can be cast to a L{Data<escript.Data>} object on
2638        @keyword y: value for coefficient y                 L{Function<escript.Function>}
2639        @type y: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.        @keyword B_reduced: value for coefficient C{B_reduced}
2640        @keyword d_contact: value for coefficient d_contact        @type B_reduced: any type that can be cast to a L{Data<escript.Data>}
2641        @type d_contact: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnContactOne<escript.FunctionOnContactOne>}.                         object on L{ReducedFunction<escript.ReducedFunction>}
2642                         or  L{FunctionOnContactZero<escript.FunctionOnContactZero>}.        @keyword C: value for coefficient C{C}
2643        @keyword y_contact: value for coefficient y_contact        @type C: any type that can be cast to a L{Data<escript.Data>} object on
2644        @type y_contact: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnContactOne<escript.FunctionOnContactOne>}.                 L{Function<escript.Function>}
2645                         or  L{FunctionOnContactZero<escript.FunctionOnContactZero>}.        @keyword C_reduced: value for coefficient C{C_reduced}
2646          @type C_reduced: any type that can be cast to a L{Data<escript.Data>}
2647                           object on L{ReducedFunction<escript.ReducedFunction>}
2648          @keyword D: value for coefficient C{D}
2649          @type D: any type that can be cast to a L{Data<escript.Data>} object on
2650                   L{Function<escript.Function>}
2651          @keyword D_reduced: value for coefficient C{D_reduced}
2652          @type D_reduced: any type that can be cast to a L{Data<escript.Data>}
2653                           object on L{ReducedFunction<escript.ReducedFunction>}
2654          @keyword X: value for coefficient C{X}
2655          @type X: any type that can be cast to a L{Data<escript.Data>} object on
2656                   L{Function<escript.Function>}
2657          @keyword X_reduced: value for coefficient C{X_reduced}
2658          @type X_reduced: any type that can be cast to a L{Data<escript.Data>}
2659                           object on L{ReducedFunction<escript.ReducedFunction>}
2660          @keyword Y: value for coefficient C{Y}
2661          @type Y: any type that can be cast to a L{Data<escript.Data>} object on
2662                   L{Function<escript.Function>}
2663          @keyword Y_reduced: value for coefficient C{Y_reduced}
2664          @type Y_reduced: any type that can be cast to a L{Data<escript.Data>}
2665                           object on L{ReducedFunction<escript.Function>}
2666          @keyword m: value for coefficient C{m}
2667          @type m: any type that can be cast to a L{Data<escript.Data>} object on
2668                   L{FunctionOnBoundary<escript.FunctionOnBoundary>}
2669          @keyword m_reduced: value for coefficient C{m_reduced}
2670          @type m_reduced: any type that can be cast to a L{Data<escript.Data>}
2671                           object on L{FunctionOnBoundary<escript.ReducedFunctionOnBoundary>}
2672          @keyword d: value for coefficient C{d}
2673          @type d: any type that can be cast to a L{Data<escript.Data>} object on
2674                   L{FunctionOnBoundary<escript.FunctionOnBoundary>}
2675          @keyword d_reduced: value for coefficient C{d_reduced}
2676          @type d_reduced: any type that can be cast to a L{Data<escript.Data>}
2677                           object on L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}
2678          @keyword y: value for coefficient C{y}
2679          @type y: any type that can be cast to a L{Data<escript.Data>} object on
2680                   L{FunctionOnBoundary<escript.FunctionOnBoundary>}
2681          @keyword d_contact: value for coefficient C{d_contact}
2682          @type d_contact: any type that can be cast to a L{Data<escript.Data>}
2683                           object on L{FunctionOnContactOne<escript.FunctionOnContactOne>} or L{FunctionOnContactZero<escript.FunctionOnContactZero>}
2684          @keyword d_contact_reduced: value for coefficient C{d_contact_reduced}
2685          @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>}
2686          @keyword y_contact: value for coefficient C{y_contact}
2687          @type y_contact: any type that can be cast to a L{Data<escript.Data>}
2688                           object on L{FunctionOnContactOne<escript.FunctionOnContactOne>} or L{FunctionOnContactZero<escript.FunctionOnContactZero>}
2689          @keyword y_contact_reduced: value for coefficient C{y_contact_reduced}
2690          @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>}
2691        @keyword r: values prescribed to the solution at the locations of constraints        @keyword r: values prescribed to the solution at the locations of constraints
2692        @type r: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}        @type r: any type that can be cast to a L{Data<escript.Data>} object on
2693                 depending of reduced order is used for the solution.                 L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2694        @keyword q: mask for location of constraints                 depending on whether reduced order is used for the solution
2695        @type q: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}        @keyword q: mask for the location of constraints
2696                 depending of reduced order is used for the representation of the equation.        @type q: any type that can be cast to a L{Data<escript.Data>} object on
2697        @raise IllegalCoefficient: if an unknown coefficient keyword is used.                 L{Solution<escript.Solution>} or
2698                   L{ReducedSolution<escript.ReducedSolution>} depending on whether
2699        """                 reduced order is used for the representation of the equation
2700        if "A" in coefficients.keys()   or "B" in coefficients.keys() or "C" in coefficients.keys(): self.__Xi=escript.Data()        @raise IllegalCoefficient: if an unknown coefficient keyword is used
2701        super(AdvectivePDE, self).setValue(**coefficients)        """
2702          super(TransportPDE,self).setValue(**coefficients)
2703     def ELMAN_RAMAGE(self,P):  
2704       """     def createOperator(self):
2705       Predefined function to set a values for M{S{xi}} from a Preclet number M{P}.         """
2706       This function uses the method suggested by H.C. Elman and A. Ramage, I{SIAM J. Numer. Anal.}, B{40} (2002)         Returns an instance of a new transport operator.
2707            - M{S{xi}(P)=0} for M{P<1}         """
2708            - M{S{xi}(P)=(1-1/P)/2} otherwise         self.trace("New Transport problem is allocated.")
2709           return self.getDomain().newTransportProblem( \
2710       @param P: Preclet number                                 self.getTheta(),
2711       @type P: L{Scalar<escript.Scalar>}                                 self.getNumEquations(), \
2712       @return: up-wind weightimg factor                                 self.getFunctionSpaceForSolution(), \
2713     &nb