/[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 1859 by gross, Wed Oct 8 03:03:37 2008 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}. L{linearPDEs} does not provide any
25  solver capabilities in itself but hands the PDE over to  solver capabilities in itself but hands the PDE over to
26  the PDE solver library defined through the L{Domain<escript.Domain>} of the PDE.  the PDE solver library defined through the L{Domain<escript.Domain>} of the PDE.
27  The general interface is provided through the L{LinearPDE} class. The  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 by its advective terms.
 provides an interface to PDE dominated by its advective terms. The L{Poisson},  
 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.  
29    
30  @var __author__: name of author  @var __author__: name of author
31    @var __copyright__: copyrights
32  @var __license__: licence agreement  @var __license__: licence agreement
33  @var __url__: url entry point on documentation  @var __url__: url entry point on documentation
34  @var __version__: version  @var __version__: version
35  @var __date__: date of the version  @var __date__: date of the version
36  """  """
37    
38    import math
39  import escript  import escript
40  import util  import util
41  import numarray  import numarray
42    
43  __author__="Lutz Gross, l.gross@uq.edu.au"  __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$"  
44    
45    
46  class IllegalCoefficient(ValueError):  class IllegalCoefficient(ValueError):
47     """     """
48     raised if an illegal coefficient of the general ar particular PDE is requested.     raised if an illegal coefficient of the general ar particular PDE is requested.
49     """     """
50       pass
51    
52  class IllegalCoefficientValue(ValueError):  class IllegalCoefficientValue(ValueError):
53     """     """
54     raised if an incorrect value for a coefficient is used.     raised if an incorrect value for a coefficient is used.
55     """     """
56       pass
57    
58    class IllegalCoefficientFunctionSpace(ValueError):
59       """
60       raised if an incorrect function space for a coefficient is used.
61       """
62    
63  class UndefinedPDEError(ValueError):  class UndefinedPDEError(ValueError):
64     """     """
65     raised if a PDE is not fully defined yet.     raised if a PDE is not fully defined yet.
66     """     """
67       pass
68    
69  class PDECoefficient(object):  class PDECoef(object):
70      """      """
71      A class for describing a PDE coefficient      A class for describing a PDE coefficient
72    
73      @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 the PDE domain
74      @cvar BOUNDARY: indicator that coefficient is defined on the boundary of the PDE domain      @cvar BOUNDARY: indicator that coefficient is defined on the boundary of the PDE domain
75      @cvar CONTACT: indicator that coefficient is defined on the contact region within the PDE domain      @cvar CONTACT: indicator that coefficient is defined on the contact region within the PDE domain
76        @cvar INTERIOR_REDUCED: indicator that coefficient is defined on the interior of the PDE domain using a reduced integration order
77        @cvar BOUNDARY_REDUCED: indicator that coefficient is defined on the boundary of the PDE domain using a reduced integration order
78        @cvar CONTACT_REDUCED: indicator that coefficient is defined on the contact region within the PDE domain using a reduced integration order
79      @cvar SOLUTION: indicator that coefficient is defined trough a solution of the PDE      @cvar SOLUTION: indicator that coefficient is defined trough a solution of the PDE
80      @cvar REDUCED: indicator that coefficient is defined trough a reduced solution of the PDE      @cvar REDUCED: indicator that coefficient is defined trough a reduced solution of the PDE
81      @cvar BY_EQUATION: indicator that the dimension of the coefficient shape is defined by the number PDE equations      @cvar BY_EQUATION: indicator that the dimension of the coefficient shape is defined by the number PDE equations
# Line 76  class PDECoefficient(object): Line 97  class PDECoefficient(object):
97      OPERATOR=10      OPERATOR=10
98      RIGHTHANDSIDE=11      RIGHTHANDSIDE=11
99      BOTH=12      BOTH=12
100        INTERIOR_REDUCED=13
101        BOUNDARY_REDUCED=14
102        CONTACT_REDUCED=15
103    
104      def __init__(self,where,pattern,altering):      def __init__(self, where, pattern, altering):
105         """         """
106         Initialise a PDE Coefficient type         Initialise a PDE Coefficient type
107    
108         @param where: describes where the coefficient lives         @param where: describes where the coefficient lives
109         @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}, L{REDUCED},
110                               L{INTERIOR_REDUCED}, L{BOUNDARY_REDUCED}, L{CONTACT_REDUCED}.
111         @param pattern: describes the shape of the coefficient and how the shape is build for a given         @param pattern: describes the shape of the coefficient and how the shape is build for a given
112                spatial dimension and numbers of equation and solution in then PDE. For instance,                spatial dimension and numbers of equation and solution in then PDE. For instance,
113                (L{BY_EQUATION},L{BY_SOLUTION},L{BY_DIM}) descrbes a rank 3 coefficient which                (L{BY_EQUATION},L{BY_SOLUTION},L{BY_DIM}) descrbes a rank 3 coefficient which
# Line 93  class PDECoefficient(object): Line 118  class PDECoefficient(object):
118         @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}
119         @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 coefficiennt is altered
120         @type altering: one of L{OPERATOR}, L{RIGHTHANDSIDE}, L{BOTH}         @type altering: one of L{OPERATOR}, L{RIGHTHANDSIDE}, L{BOTH}
121           @param reduced: indicates if reduced
122           @type reduced: C{bool}
123         """         """
124         super(PDECoefficient, self).__init__()         super(PDECoef, self).__init__()
125         self.what=where         self.what=where
126         self.pattern=pattern         self.pattern=pattern
127         self.altering=altering         self.altering=altering
# Line 114  class PDECoefficient(object): Line 140  class PDECoefficient(object):
140         @param domain: domain on which the PDE uses the coefficient         @param domain: domain on which the PDE uses the coefficient
141         @type domain: L{Domain<escript.Domain>}         @type domain: L{Domain<escript.Domain>}
142         @param reducedEquationOrder: True to indicate that reduced order is used to represent the equation         @param reducedEquationOrder: True to indicate that reduced order is used to represent the equation
143         @type domain: C{bool}         @type reducedEquationOrder: C{bool}
144         @param reducedSolutionOrder: True to indicate that reduced order is used to represent the solution         @param reducedSolutionOrder: True to indicate that reduced order is used to represent the solution
145         @type domain: C{bool}         @type reducedSolutionOrder: C{bool}
146         @return:  L{FunctionSpace<escript.FunctionSpace>} of the coefficient         @return:  L{FunctionSpace<escript.FunctionSpace>} of the coefficient
147         @rtype:  L{FunctionSpace<escript.FunctionSpace>}         @rtype:  L{FunctionSpace<escript.FunctionSpace>}
148         """         """
149         if self.what==self.INTERIOR:         if self.what==self.INTERIOR:
150              return escript.Function(domain)              return escript.Function(domain)
151           elif self.what==self.INTERIOR_REDUCED:
152                return escript.ReducedFunction(domain)
153         elif self.what==self.BOUNDARY:         elif self.what==self.BOUNDARY:
154              return escript.FunctionOnBoundary(domain)              return escript.FunctionOnBoundary(domain)
155           elif self.what==self.BOUNDARY_REDUCED:
156                return escript.ReducedFunctionOnBoundary(domain)
157         elif self.what==self.CONTACT:         elif self.what==self.CONTACT:
158              return escript.FunctionOnContactZero(domain)              return escript.FunctionOnContactZero(domain)
159           elif self.what==self.CONTACT_REDUCED:
160                return escript.ReducedFunctionOnContactZero(domain)
161         elif self.what==self.SOLUTION:         elif self.what==self.SOLUTION:
162              if reducedEquationOrder and reducedSolutionOrder:              if reducedEquationOrder and reducedSolutionOrder:
163                  return escript.ReducedSolution(domain)                  return escript.ReducedSolution(domain)
# Line 154  class PDECoefficient(object): Line 186  class PDECoefficient(object):
186         @param numSolutions: number of components of the PDE solution         @param numSolutions: number of components of the PDE solution
187         @type numSolutions: C{int}         @type numSolutions: C{int}
188         @param reducedEquationOrder: True to indicate that reduced order is used to represent the equation         @param reducedEquationOrder: True to indicate that reduced order is used to represent the equation
189         @type domain: C{bool}         @type reducedEquationOrder: C{bool}
190         @param reducedSolutionOrder: True to indicate that reduced order is used to represent the solution         @param reducedSolutionOrder: True to indicate that reduced order is used to represent the solution
191         @type domain: C{bool}         @type reducedSolutionOrder: C{bool}
192         @param newValue: number of components of the PDE solution         @param newValue: number of components of the PDE solution
193         @type newValue: any object that can be converted into a L{Data<escript.Data>} object with the appropriate shape and L{FunctionSpace<escript.FunctionSpace>}         @type newValue: any object that can be converted into a L{Data<escript.Data>} object with the appropriate shape and L{FunctionSpace<escript.FunctionSpace>}
194         @raise IllegalCoefficientValue: if the shape of the assigned value does not match the shape of the coefficient         @raise IllegalCoefficientValue: if the shape of the assigned value does not match the shape of the coefficient
195           @raise IllegalCoefficientFunctionSpace: if unable to interploate value to appropriate function space
196         """         """
197         if newValue==None:         if newValue==None:
198             newValue=escript.Data()             newValue=escript.Data()
199         elif isinstance(newValue,escript.Data):         elif isinstance(newValue,escript.Data):
200             if not newValue.isEmpty():             if not newValue.isEmpty():
201                try:                if not newValue.getFunctionSpace() == self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder):
202                   newValue=escript.Data(newValue,self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder))                  try:
203                except:                    newValue=escript.Data(newValue,self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder))
204                   raise IllegalCoefficientValue,"Unable to interpolate coefficient to function space %s"%self.getFunctionSpace(domain)                  except:
205                      raise IllegalCoefficientFunctionSpace,"Unable to interpolate coefficient to function space %s"%self.getFunctionSpace(domain)
206         else:         else:
207             newValue=escript.Data(newValue,self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder))             newValue=escript.Data(newValue,self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder))
208         if not newValue.isEmpty():         if not newValue.isEmpty():
# Line 305  class PDECoefficient(object): Line 339  class PDECoefficient(object):
339                  s=s+(dim,)                  s=s+(dim,)
340         return s         return s
341    
342  class LinearPDE(object):  class LinearProblem(object):
343     """     """
344     This class is used to define a general linear, steady, second order PDE     This is the base class is to define a general linear PDE-type problem for an u
345     for an unknown function M{u} on a given domain defined through a L{Domain<escript.Domain>} object.     for an unknown function M{u} on a given domain defined through a L{Domain<escript.Domain>} object.
346       The problem can be given as a single equations or as a system of equations.
347     For a single PDE with a solution with a single component the linear PDE is defined in the following form:  
348       The class assumes that some sort of assembling process is required to form a problem of the form
349     M{-grad(A[j,l]*grad(u)[l]+B[j]u)[j]+C[l]*grad(u)[l]+D*u =-grad(X)[j,j]+Y}  
350       M{L u=f}
351     where M{grad(F)} denotes the spatial derivative of M{F}. Einstein's summation convention,    
352     ie. summation over indexes appearing twice in a term of a sum is performed, is used.     where M{L} is an operator and M{f} is the right hand side. This operator problem will be solved to get the
353     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     unknown M{u}.
    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]}  
   
    In the context of discontinuities M{n} denotes the normal on the discontinuity pointing from side 0 towards side 1  
    calculated from L{getNormal<escript.FunctionSpace.getNormal>} of L{FunctionOnContactZero<escript.FunctionOnContactZero>}. For a system of PDEs  
    the contact condition takes the form  
   
    M{n[j]*J0[i,j]=n[j]*J1[i,j]=y_contact[i]- d_contact[i,k]*jump(u)[k]}  
   
    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  
    of the solution at side 1 and at side 0, denotes the jump of M{u} across discontinuity along the normal calcualted by  
    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  
   
    M{n[j]*J0_{j}=n[j]*J1_{j}=y_contact-d_contact*jump(u)}  
   
    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>}.  
354    
355     @cvar DEFAULT: The default method used to solve the system of linear equations     @cvar DEFAULT: The default method used to solve the system of linear equations
356     @cvar DIRECT: The direct solver based on LDU factorization     @cvar DIRECT: The direct solver based on LDU factorization
# Line 398  class LinearPDE(object): Line 359  class LinearPDE(object):
359     @cvar CR: The conjugate residual method     @cvar CR: The conjugate residual method
360     @cvar CGS: The conjugate gardient square method     @cvar CGS: The conjugate gardient square method
361     @cvar BICGSTAB: The stabilized BiConjugate Gradient method.     @cvar BICGSTAB: The stabilized BiConjugate Gradient method.
362       @cvar TFQMR: Transport Free Quasi Minimal Residual method.
363       @cvar MINRES: Minimum residual method.
364     @cvar SSOR: The symmetric overrealaxtion method     @cvar SSOR: The symmetric overrealaxtion method
365     @cvar ILU0: The incomplete LU factorization preconditioner  with no fill in     @cvar ILU0: The incomplete LU factorization preconditioner  with no fill in
366     @cvar ILUT: The incomplete LU factorization preconditioner with will in     @cvar ILUT: The incomplete LU factorization preconditioner with will in
# Line 412  class LinearPDE(object): Line 375  class LinearPDE(object):
375     @cvar SCSL: SGI SCSL solver library     @cvar SCSL: SGI SCSL solver library
376     @cvar MKL: Intel's MKL solver library     @cvar MKL: Intel's MKL solver library
377     @cvar UMFPACK: the UMFPACK library     @cvar UMFPACK: the UMFPACK library
378       @cvar TRILINOS: the TRILINOS parallel solver class library from Sandia Natl Labs
379     @cvar ITERATIVE: The default iterative solver     @cvar ITERATIVE: The default iterative solver
380     @cvar AMG: algebraic multi grid     @cvar AMG: algebraic multi grid
381     @cvar RILU: recursive ILU     @cvar RILU: recursive ILU
382       @cvar GS: Gauss-Seidel solver
383    
384     """     """
385     DEFAULT= 0     DEFAULT= 0
# Line 441  class LinearPDE(object): Line 406  class LinearPDE(object):
406     PASO= 21     PASO= 21
407     AMG= 22     AMG= 22
408     RILU = 23     RILU = 23
409       TRILINOS = 24
410       NONLINEAR_GMRES = 25
411       TFQMR = 26
412       MINRES = 27
413       GS=28
414    
415     SMALL_TOLERANCE=1.e-13     SMALL_TOLERANCE=1.e-13
416     __PACKAGE_KEY="package"     PACKAGE_KEY="package"
417     __METHOD_KEY="method"     METHOD_KEY="method"
418     __SYMMETRY_KEY="symmetric"     SYMMETRY_KEY="symmetric"
419     __TOLERANCE_KEY="tolerance"     TOLERANCE_KEY="tolerance"
420     __PRECONDITIONER_KEY="preconditioner"     PRECONDITIONER_KEY="preconditioner"
421    
422    
423     def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):     def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):
424       """       """
425       initializes a new linear PDE       initializes a linear problem
426    
427       @param domain: domain of the PDE       @param domain: domain of the PDE
428       @type domain: L{Domain<escript.Domain>}       @type domain: L{Domain<escript.Domain>}
429       @param numEquations: number of equations. If numEquations==None the number of equations       @param numEquations: number of equations. If numEquations==None the number of equations
430                            is exracted from the PDE coefficients.                            is extracted from the coefficients.
431       @param numSolutions: number of solution components. If  numSolutions==None the number of solution components       @param numSolutions: number of solution components. If  numSolutions==None the number of solution components
432                            is exracted from the PDE coefficients.                            is extracted from the coefficients.
433       @param debug: if True debug informations are printed.       @param debug: if True debug informations are printed.
434    
435       """       """
436       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)}  
437    
      # COEFFICIENTS can be overwritten by subclasses:  
      self.COEFFICIENTS=self.__COEFFICIENTS_OF_GENEARL_PDE  
      self.__altered_coefficients=False  
      # initialize attributes  
438       self.__debug=debug       self.__debug=debug
439       self.__domain=domain       self.__domain=domain
440       self.__numEquations=numEquations       self.__numEquations=numEquations
441       self.__numSolutions=numSolutions       self.__numSolutions=numSolutions
442       self.__resetSystem()       self.__altered_coefficients=False
   
      # set some default values:  
443       self.__reduce_equation_order=False       self.__reduce_equation_order=False
444       self.__reduce_solution_order=False       self.__reduce_solution_order=False
445       self.__tolerance=1.e-8       self.__tolerance=1.e-8
446       self.__solver_method=self.DEFAULT       self.__solver_method=self.DEFAULT
447       self.__solver_package=self.DEFAULT       self.__solver_package=self.DEFAULT
448       self.__preconditioner=self.DEFAULT       self.__preconditioner=self.DEFAULT
449       self.__matrix_type=self.__domain.getSystemMatrixTypeId(self.DEFAULT,self.DEFAULT,False)       self.__is_RHS_valid=False
450         self.__is_operator_valid=False
451       self.__sym=False       self.__sym=False
452         self.__COEFFICIENTS={}
453       self.resetCoefficients()       # initialize things:
454       self.trace("PDE Coeffients are %s"%str(self.COEFFICIENTS.keys()))       self.resetAllCoefficients()
455         self.__system_type=None
456         self.updateSystemType()
457     # =============================================================================     # =============================================================================
458     #    general stuff:     #    general stuff:
459     # =============================================================================     # =============================================================================
# Line 513  class LinearPDE(object): Line 464  class LinearPDE(object):
464       @return: a simple representation of the PDE       @return: a simple representation of the PDE
465       @rtype: C{str}       @rtype: C{str}
466       """       """
467       return "<LinearPDE %d>"%id(self)       return "<LinearProblem %d>"%id(self)
468     # =============================================================================     # =============================================================================
469     #    debug :     #    debug :
470     # =============================================================================     # =============================================================================
# Line 540  class LinearPDE(object): Line 491  class LinearPDE(object):
491     # =============================================================================     # =============================================================================
492     # some service functions:     # some service functions:
493     # =============================================================================     # =============================================================================
494       def introduceCoefficients(self,**coeff):
495           """
496           introduces a new coefficient into the problem.
497    
498           use
499    
500             self.introduceCoefficients(A=PDECoef(...),B=PDECoef(...))
501    
502           to introduce the coefficients M{A} ans M{B}.
503           """
504           for name, type in coeff.items():
505               if not isinstance(type,PDECoef):
506                  raise ValueError,"coefficient %s has no type."%name
507               self.__COEFFICIENTS[name]=type
508               self.__COEFFICIENTS[name].resetValue()
509               self.trace("coefficient %s has been introduced."%name)
510    
511     def getDomain(self):     def getDomain(self):
512       """       """
513       returns the domain of the PDE       returns the domain of the PDE
# Line 567  class LinearPDE(object): Line 535  class LinearPDE(object):
535       @raise UndefinedPDEError: if the number of equations is not be specified yet.       @raise UndefinedPDEError: if the number of equations is not be specified yet.
536       """       """
537       if self.__numEquations==None:       if self.__numEquations==None:
538           raise UndefinedPDEError,"Number of equations is undefined. Please specify argument numEquations."           if self.__numSolutions==None:
539       else:              raise UndefinedPDEError,"Number of equations is undefined. Please specify argument numEquations."
540           return self.__numEquations           else:
541                self.__numEquations=self.__numSolutions
542         return self.__numEquations
543    
544     def getNumSolutions(self):     def getNumSolutions(self):
545       """       """
# Line 580  class LinearPDE(object): Line 550  class LinearPDE(object):
550       @raise UndefinedPDEError: if the number of unknowns is not be specified yet.       @raise UndefinedPDEError: if the number of unknowns is not be specified yet.
551       """       """
552       if self.__numSolutions==None:       if self.__numSolutions==None:
553          raise UndefinedPDEError,"Number of solution is undefined. Please specify argument numSolutions."          if self.__numEquations==None:
554       else:              raise UndefinedPDEError,"Number of solution is undefined. Please specify argument numSolutions."
555          return self.__numSolutions          else:
556                self.__numSolutions=self.__numEquations
557         return self.__numSolutions
558    
559     def reduceEquationOrder(self):     def reduceEquationOrder(self):
560       """       """
# Line 626  class LinearPDE(object): Line 598  class LinearPDE(object):
598       else:       else:
599           return escript.Solution(self.getDomain())           return escript.Solution(self.getDomain())
600    
   
    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")  
601     # =============================================================================     # =============================================================================
602     #   solver settings:     #   solver settings:
603     # =============================================================================     # =============================================================================
# Line 819  class LinearPDE(object): Line 606  class LinearPDE(object):
606         sets a new solver         sets a new solver
607    
608         @param solver: sets a new solver method.         @param solver: sets a new solver method.
609         @type solver: one of L{DEFAULT}, L{ITERATIVE} L{DIRECT}, L{CHOLEVSKY}, L{PCG}, L{CR}, L{CGS}, L{BICGSTAB}, L{SSOR}, L{GMRES}, L{PRES20}, L{LUMPING}, L{AMG}         @type solver: one of L{DEFAULT}, L{ITERATIVE} L{DIRECT}, L{CHOLEVSKY}, L{PCG}, L{CR}, L{CGS}, L{BICGSTAB}, L{SSOR}, L{GMRES}, L{TFQMR}, L{MINRES}, L{PRES20}, L{LUMPING}, L{AMG}
610         @param preconditioner: sets a new solver method.         @param preconditioner: sets a new solver method.
611         @type preconditioner: one of L{DEFAULT}, L{JACOBI} L{ILU0}, L{ILUT},L{SSOR}, L{RILU}         @type preconditioner: one of L{DEFAULT}, L{JACOBI} L{ILU0}, L{ILUT},L{SSOR}, L{RILU},  L{GS}
612    
613           @note: the solver method choosen may not be available by the selected package.
614         """         """
615         if solver==None: solve=self.DEFAULT         if solver==None: solver=self.__solver_method
616           if preconditioner==None: preconditioner=self.__preconditioner
617           if solver==None: solver=self.DEFAULT
618         if preconditioner==None: preconditioner=self.DEFAULT         if preconditioner==None: preconditioner=self.DEFAULT
619         if not (solver,preconditioner)==self.getSolverMethod():         if not (solver,preconditioner)==self.getSolverMethod():
620             self.__solver_method=solver             self.__solver_method=solver
621             self.__preconditioner=preconditioner             self.__preconditioner=preconditioner
622             self.__checkMatrixType()             self.updateSystemType()
623             self.trace("New solver is %s"%self.getSolverMethodName())             self.trace("New solver is %s"%self.getSolverMethodName())
624    
625     def getSolverMethodName(self):     def getSolverMethodName(self):
626         """         """
627         returns the name of the solver currently used         returns the name of the solver currently used
628    
# Line 847  class LinearPDE(object): Line 638  class LinearPDE(object):
638         elif m[0]==self.ITERATIVE: method= "ITERATIVE"         elif m[0]==self.ITERATIVE: method= "ITERATIVE"
639         elif m[0]==self.CHOLEVSKY: method= "CHOLEVSKY"         elif m[0]==self.CHOLEVSKY: method= "CHOLEVSKY"
640         elif m[0]==self.PCG: method= "PCG"         elif m[0]==self.PCG: method= "PCG"
641           elif m[0]==self.TFQMR: method= "TFQMR"
642           elif m[0]==self.MINRES: method= "MINRES"
643         elif m[0]==self.CR: method= "CR"         elif m[0]==self.CR: method= "CR"
644         elif m[0]==self.CGS: method= "CGS"         elif m[0]==self.CGS: method= "CGS"
645         elif m[0]==self.BICGSTAB: method= "BICGSTAB"         elif m[0]==self.BICGSTAB: method= "BICGSTAB"
# Line 862  class LinearPDE(object): Line 655  class LinearPDE(object):
655         elif m[1]==self.SSOR: method+= "+SSOR"         elif m[1]==self.SSOR: method+= "+SSOR"
656         elif m[1]==self.AMG: method+= "+AMG"         elif m[1]==self.AMG: method+= "+AMG"
657         elif m[1]==self.RILU: method+= "+RILU"         elif m[1]==self.RILU: method+= "+RILU"
658           elif m[1]==self.GS: method+= "+GS"
659         if p==self.DEFAULT: package="DEFAULT"         if p==self.DEFAULT: package="DEFAULT"
660         elif p==self.PASO: package= "PASO"         elif p==self.PASO: package= "PASO"
661         elif p==self.MKL: package= "MKL"         elif p==self.MKL: package= "MKL"
662         elif p==self.SCSL: package= "SCSL"         elif p==self.SCSL: package= "SCSL"
663         elif p==self.UMFPACK: package= "UMFPACK"         elif p==self.UMFPACK: package= "UMFPACK"
664           elif p==self.TRILINOS: package= "TRILINOS"
665         else : method="unknown"         else : method="unknown"
666         return "%s solver of %s package"%(method,package)         return "%s solver of %s package"%(method,package)
667    
668       def getSolverMethod(self):
    def getSolverMethod(self):  
669         """         """
670         returns the solver method         returns the solver method and preconditioner used
671    
672         @return: the solver method currently be used.         @return: the solver method currently be used.
673         @rtype: C{int}         @rtype: C{tuple} of C{int}
674         """         """
675         return self.__solver_method,self.__preconditioner         return self.__solver_method,self.__preconditioner
676    
677     def setSolverPackage(self,package=None):     def setSolverPackage(self,package=None):
678         """         """
679         sets a new solver package         sets a new solver package
680    
681         @param solver: sets a new solver method.         @param package: sets a new solver method.
682         @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}, L{TRILINOS}
683         """         """
684         if package==None: package=self.DEFAULT         if package==None: package=self.DEFAULT
685         if not package==self.getSolverPackage():         if not package==self.getSolverPackage():
686             self.__solver_method=solver             self.__solver_package=package
687             self.__checkMatrixType()             self.updateSystemType()
688             self.trace("New solver is %s"%self.getSolverMethodName())             self.trace("New solver is %s"%self.getSolverMethodName())
689    
690     def getSolverPackage(self):     def getSolverPackage(self):
691         """         """
692         returns the package of the solver         returns the package of the solver
693    
# Line 904  class LinearPDE(object): Line 698  class LinearPDE(object):
698    
699     def isUsingLumping(self):     def isUsingLumping(self):
700        """        """
701        checks if matrix lumping is used a solver method        checks if matrix lumping is used as a solver method
702    
703        @return: True is lumping is currently used a solver method.        @return: True is lumping is currently used a solver method.
704        @rtype: C{bool}        @rtype: C{bool}
705        """        """
706        return self.getSolverMethod()[0]==self.LUMPING        return self.getSolverMethod()[0]==self.LUMPING
707    
708     def setTolerance(self,tol=1.e-8):     def setTolerance(self,tol=1.e-8):
709         """         """
710         resets the tolerance for the solver method to tol where for an appropriate norm M{|.|}         resets the tolerance for the solver method to tol.
   
        M{|L{getResidual}()|<tol*|L{getRightHandSide}()|}  
   
        defines the stopping criterion.  
711    
712         @param tol: new tolerance for the solver. If the tol is lower then the current tolerence         @param tol: new tolerance for the solver. If the tol is lower then the current tolerence
713                     the system will be resolved.                     the system will be resolved.
714         @type tol: positive C{float}         @type tol: positive C{float}
715         @raise ValueException: if tolerance is not positive.         @raise ValueError: if tolerance is not positive.
716         """         """
717         if not tol>0:         if not tol>0:
718             raise ValueException,"Tolerance as to be positive"             raise ValueError,"Tolerance as to be positive"
719         if tol<self.getTolerance(): self.__invalidateSolution()         if tol<self.getTolerance(): self.invalidateSolution()
720         self.trace("New tolerance %e"%tol)         self.trace("New tolerance %e"%tol)
721         self.__tolerance=tol         self.__tolerance=tol
722         return         return
723    
724     def getTolerance(self):     def getTolerance(self):
725         """         """
726         returns the tolerance set for the solution         returns the tolerance set for the solution
727    
# Line 943  class LinearPDE(object): Line 733  class LinearPDE(object):
733     # =============================================================================     # =============================================================================
734     #    symmetry  flag:     #    symmetry  flag:
735     # =============================================================================     # =============================================================================
736     def isSymmetric(self):     def isSymmetric(self):
737        """        """
738        checks if symmetry is indicated.        checks if symmetry is indicated.
739    
# Line 952  class LinearPDE(object): Line 742  class LinearPDE(object):
742        """        """
743        return self.__sym        return self.__sym
744    
745     def setSymmetryOn(self):     def setSymmetryOn(self):
746        """        """
747        sets the symmetry flag.        sets the symmetry flag.
748        """        """
749        if not self.isSymmetric():        if not self.isSymmetric():
750           self.trace("PDE is set to be symmetric")           self.trace("PDE is set to be symmetric")
751           self.__sym=True           self.__sym=True
752           self.__checkMatrixType()           self.updateSystemType()
753    
754     def setSymmetryOff(self):     def setSymmetryOff(self):
755        """        """
756        removes the symmetry flag.        removes the symmetry flag.
757        """        """
758        if self.isSymmetric():        if self.isSymmetric():
759           self.trace("PDE is set to be unsymmetric")           self.trace("PDE is set to be unsymmetric")
760           self.__sym=False           self.__sym=False
761           self.__checkMatrixType()           self.updateSystemType()
762    
763     def setSymmetryTo(self,flag=False):     def setSymmetryTo(self,flag=False):
764        """        """
765        sets the symmetry flag to flag        sets the symmetry flag to flag
766    
# Line 985  class LinearPDE(object): Line 775  class LinearPDE(object):
775     # =============================================================================     # =============================================================================
776     # function space handling for the equation as well as the solution     # function space handling for the equation as well as the solution
777     # =============================================================================     # =============================================================================
778     def setReducedOrderOn(self):     def setReducedOrderOn(self):
779       """       """
780       switches on reduced order for solution and equation representation       switches on reduced order for solution and equation representation
781    
# Line 994  class LinearPDE(object): Line 784  class LinearPDE(object):
784       self.setReducedOrderForSolutionOn()       self.setReducedOrderForSolutionOn()
785       self.setReducedOrderForEquationOn()       self.setReducedOrderForEquationOn()
786    
787     def setReducedOrderOff(self):     def setReducedOrderOff(self):
788       """       """
789       switches off reduced order for solution and equation representation       switches off reduced order for solution and equation representation
790    
# Line 1003  class LinearPDE(object): Line 793  class LinearPDE(object):
793       self.setReducedOrderForSolutionOff()       self.setReducedOrderForSolutionOff()
794       self.setReducedOrderForEquationOff()       self.setReducedOrderForEquationOff()
795    
796     def setReducedOrderTo(self,flag=False):     def setReducedOrderTo(self,flag=False):
797       """       """
798       sets order reduction for both solution and equation representation according to flag.       sets order reduction for both solution and equation representation according to flag.
799       @param flag: if flag is True, the order reduction is switched on for both  solution and equation representation, otherwise or       @param flag: if flag is True, the order reduction is switched on for both  solution and equation representation, otherwise or
# Line 1015  class LinearPDE(object): Line 805  class LinearPDE(object):
805       self.setReducedOrderForEquationTo(flag)       self.setReducedOrderForEquationTo(flag)
806    
807    
808     def setReducedOrderForSolutionOn(self):     def setReducedOrderForSolutionOn(self):
809       """       """
810       switches on reduced order for solution representation       switches on reduced order for solution representation
811    
# Line 1026  class LinearPDE(object): Line 816  class LinearPDE(object):
816                raise RuntimeError,"order cannot be altered after coefficients have been defined."                raise RuntimeError,"order cannot be altered after coefficients have been defined."
817           self.trace("Reduced order is used to solution representation.")           self.trace("Reduced order is used to solution representation.")
818           self.__reduce_solution_order=True           self.__reduce_solution_order=True
819           self.__resetSystem()           self.initializeSystem()
820    
821     def setReducedOrderForSolutionOff(self):     def setReducedOrderForSolutionOff(self):
822       """       """
823       switches off reduced order for solution representation       switches off reduced order for solution representation
824    
# Line 1039  class LinearPDE(object): Line 829  class LinearPDE(object):
829                raise RuntimeError,"order cannot be altered after coefficients have been defined."                raise RuntimeError,"order cannot be altered after coefficients have been defined."
830           self.trace("Full order is used to interpolate solution.")           self.trace("Full order is used to interpolate solution.")
831           self.__reduce_solution_order=False           self.__reduce_solution_order=False
832           self.__resetSystem()           self.initializeSystem()
833    
834     def setReducedOrderForSolutionTo(self,flag=False):     def setReducedOrderForSolutionTo(self,flag=False):
835       """       """
836       sets order for test functions according to flag       sets order for test functions according to flag
837    
# Line 1055  class LinearPDE(object): Line 845  class LinearPDE(object):
845       else:       else:
846          self.setReducedOrderForSolutionOff()          self.setReducedOrderForSolutionOff()
847    
848     def setReducedOrderForEquationOn(self):     def setReducedOrderForEquationOn(self):
849       """       """
850       switches on reduced order for equation representation       switches on reduced order for equation representation
851    
# Line 1066  class LinearPDE(object): Line 856  class LinearPDE(object):
856                raise RuntimeError,"order cannot be altered after coefficients have been defined."                raise RuntimeError,"order cannot be altered after coefficients have been defined."
857           self.trace("Reduced order is used for test functions.")           self.trace("Reduced order is used for test functions.")
858           self.__reduce_equation_order=True           self.__reduce_equation_order=True
859           self.__resetSystem()           self.initializeSystem()
860    
861     def setReducedOrderForEquationOff(self):     def setReducedOrderForEquationOff(self):
862       """       """
863       switches off reduced order for equation representation       switches off reduced order for equation representation
864    
# Line 1079  class LinearPDE(object): Line 869  class LinearPDE(object):
869                raise RuntimeError,"order cannot be altered after coefficients have been defined."                raise RuntimeError,"order cannot be altered after coefficients have been defined."
870           self.trace("Full order is used for test functions.")           self.trace("Full order is used for test functions.")
871           self.__reduce_equation_order=False           self.__reduce_equation_order=False
872           self.__resetSystem()           self.initializeSystem()
873    
874     def setReducedOrderForEquationTo(self,flag=False):     def setReducedOrderForEquationTo(self,flag=False):
875       """       """
876       sets order for test functions according to flag       sets order for test functions according to flag
877    
# Line 1095  class LinearPDE(object): Line 885  class LinearPDE(object):
885       else:       else:
886          self.setReducedOrderForEquationOff()          self.setReducedOrderForEquationOff()
887    
888     # =============================================================================     def updateSystemType(self):
    # private method:  
    # =============================================================================  
    def __checkMatrixType(self):  
889       """       """
890       reassess the matrix type and, if a new matrix is needed, resets the system.       reasses the system type and, if a new matrix is needed, resets the system.
891       """       """
892       new_matrix_type=self.getDomain().getSystemMatrixTypeId(self.getSolverMethod()[0],self.getSolverPackage(),self.isSymmetric())       new_system_type=self.getRequiredSystemType()
893       if not new_matrix_type==self.__matrix_type:       if not new_system_type==self.__system_type:
894           self.trace("Matrix type is now %d."%new_matrix_type)           self.trace("Matrix type is now %d."%new_system_type)
895           self.__matrix_type=new_matrix_type           self.__system_type=new_system_type
896           self.__resetSystem()           self.initializeSystem()
897     #     def getSystemType(self):
898     #   rebuild switches :        """
899     #        return the current system type
900     def __invalidateSolution(self):        """
901         """        return self.__system_type
        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  
   
    def __invalidateRightHandSide(self):  
        """  
        indicates the right hand side has to be rebuild next time it is used  
        """  
        if self.__righthandside_isValid: self.trace("Right hand side has to be rebuilt.")  
        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)  
902    
903     def __makeFreshSolution(self):     def checkSymmetricTensor(self,name,verbose=True):
904         """        """
905         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  
906    
907     def __makeFreshRightHandSide(self):        @param name: name of the coefficient
908         """        @type name: C{str}
909         makes sure that the right hand side is instantiated and returns it initialized by zeros        @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.
910         """        @type verbose: C{bool}
911         if self.__righthandside.isEmpty():        @return:  True if coefficient C{name} is symmetric.
912             self.__righthandside=self.__getNewRightHandSide()        @rtype: C{bool}
913         else:        """
914             self.__righthandside*=0        A=self.getCoefficient(name)
915             self.trace("Right hand side is reset to zero.")        verbose=verbose or self.__debug
916         return self.__righthandside        out=True
917          if not A.isEmpty():
918             tol=util.Lsup(A)*self.SMALL_TOLERANCE
919             s=A.getShape()
920             if A.getRank() == 4:
921                if s[0]==s[2] and s[1] == s[3]:
922                   for i in range(s[0]):
923                      for j in range(s[1]):
924                         for k in range(s[2]):
925                            for l in range(s[3]):
926                                if util.Lsup(A[i,j,k,l]-A[k,l,i,j])>tol:
927                                   if verbose: print "non-symmetric problem as %s[%d,%d,%d,%d]!=%s[%d,%d,%d,%d]"%(name,i,j,k,l,name,k,l,i,j)
928                                   out=False
929                else:
930                     if verbose: print "non-symmetric problem because of inappropriate shape %s of coefficient %s."%(s,name)
931                     out=False
932             elif A.getRank() == 2:
933                if s[0]==s[1]:
934                   for j in range(s[0]):
935                      for l in range(s[1]):
936                         if util.Lsup(A[j,l]-A[l,j])>tol:
937                            if verbose: print "non-symmetric problem because %s[%d,%d]!=%s[%d,%d]"%(name,j,l,name,l,j)
938                            out=False
939                else:
940                     if verbose: print "non-symmetric problem because of inappropriate shape %s of coefficient %s."%(s,name)
941                     out=False
942             elif A.getRank() == 0:
943                pass
944             else:
945                 raise ValueError,"Cannot check rank %s of %s."%(A.getRank(),name)
946          return out
947    
948     def __makeFreshOperator(self):     def checkReciprocalSymmetry(self,name0,name1,verbose=True):
949         """        """
950         makes sure that the operator is instantiated and returns it initialized by zeros        test two coefficients for resciprocal symmetry
        """  
        if self.__operator.isEmpty():  
            self.__operator=self.__getNewOperator()  
        else:  
            self.__operator.resetValues()  
            self.trace("Operator reset to zero")  
        return self.__operator  
951    
952     def __applyConstraint(self):        @param name0: name of the first coefficient
953         """        @type name: C{str}
954         applies the constraints defined by q and r to the system        @param name1: name of the second coefficient
955         """        @type name: C{str}
956         if not self.isUsingLumping():        @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.
957            q=self.getCoefficientOfGeneralPDE("q")        @type verbose: C{bool}
958            r=self.getCoefficientOfGeneralPDE("r")        @return:  True if coefficients C{name0} and C{name1}  are resciprocally symmetric.
959            if not q.isEmpty() and not self.__operator.isEmpty():        @rtype: C{bool}
960               # q is the row and column mask to indicate where constraints are set:        """
961               row_q=escript.Data(q,self.getFunctionSpaceForEquation())        B=self.getCoefficient(name0)
962               col_q=escript.Data(q,self.getFunctionSpaceForSolution())        C=self.getCoefficient(name1)
963               u=self.__getNewSolution()        verbose=verbose or self.__debug
964               if r.isEmpty():        out=True
965                  r_s=self.__getNewSolution()        if B.isEmpty() and not C.isEmpty():
966             if verbose: print "non-symmetric problem because %s is not present but %s is"%(name0,name1)
967             out=False
968          elif not B.isEmpty() and C.isEmpty():
969             if verbose: print "non-symmetric problem because %s is not present but %s is"%(name0,name1)
970             out=False
971          elif not B.isEmpty() and not C.isEmpty():
972             sB=B.getShape()
973             sC=C.getShape()
974             tol=(util.Lsup(B)+util.Lsup(C))*self.SMALL_TOLERANCE/2.
975             if len(sB) != len(sC):
976                 if verbose: print "non-symmetric problem because ranks of %s (=%s) and %s (=%s) are different."%(name0,len(sB),name1,len(sC))
977                 out=False
978             else:
979                 if len(sB)==0:
980                   if util.Lsup(B-C)>tol:
981                      if verbose: print "non-symmetric problem because %s!=%s"%(name0,name1)
982                      out=False
983                 elif len(sB)==1:
984                   if sB[0]==sC[0]:
985                      for j in range(sB[0]):
986                         if util.Lsup(B[j]-C[j])>tol:
987                            if verbose: print "non-symmetric PDE because %s[%d]!=%s[%d]"%(name0,j,name1,j)
988                            out=False
989                   else:
990                     if verbose: print "non-symmetric problem because of inappropriate shapes %s and %s of coefficients %s and %s, respectively."%(sB,sC,name0,name1)
991                 elif len(sB)==3:
992                   if sB[0]==sC[1] and sB[1]==sC[2] and sB[2]==sC[0]:
993                       for i in range(sB[0]):
994                          for j in range(sB[1]):
995                             for k in range(sB[2]):
996                                if util.Lsup(B[i,j,k]-C[k,i,j])>tol:
997                                     if verbose: print "non-symmetric problem because %s[%d,%d,%d]!=%s[%d,%d,%d]"%(name2,i,j,k,name1,k,i,j)
998                                     out=False
999                   else:
1000                     if verbose: print "non-symmetric problem because of inappropriate shapes %s and %s of coefficients %s and %s, respectively."%(sB,sC,name0,name1)
1001               else:               else:
1002                  r_s=escript.Data(r,self.getFunctionSpaceForSolution())                   raise ValueError,"Cannot check rank %s of %s and %s."%(len(sB),name0,name1)
1003               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  
1004    
1005     # =============================================================================     def getCoefficient(self,name):
    # functions giving access to coefficients of a particular PDE implementation:  
    # =============================================================================  
    def getCoefficient(self,name):  
1006       """       """
1007       returns the value of the coefficient name       returns the value of the coefficient name
1008    
# Line 1337  class LinearPDE(object): Line 1013  class LinearPDE(object):
1013       @raise IllegalCoefficient: if name is not a coefficient of the PDE.       @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1014       """       """
1015       if self.hasCoefficient(name):       if self.hasCoefficient(name):
1016           return self.COEFFICIENTS[name].getValue()           return self.__COEFFICIENTS[name].getValue()
1017       else:       else:
1018          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1019    
1020     def hasCoefficient(self,name):     def hasCoefficient(self,name):
1021       """       """
1022       return True if name is the name of a coefficient       return True if name is the name of a coefficient
1023    
# Line 1350  class LinearPDE(object): Line 1026  class LinearPDE(object):
1026       @return: True if name is the name of a coefficient of the general PDE. Otherwise False.       @return: True if name is the name of a coefficient of the general PDE. Otherwise False.
1027       @rtype: C{bool}       @rtype: C{bool}
1028       """       """
1029       return self.COEFFICIENTS.has_key(name)       return self.__COEFFICIENTS.has_key(name)
1030    
1031     def createCoefficient(self, name):     def createCoefficient(self, name):
1032       """       """
1033       create a L{Data<escript.Data>} object corresponding to coefficient name       create a L{Data<escript.Data>} object corresponding to coefficient name
1034    
# Line 1365  class LinearPDE(object): Line 1041  class LinearPDE(object):
1041       else:       else:
1042          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1043    
1044     def getFunctionSpaceForCoefficient(self,name):     def getFunctionSpaceForCoefficient(self,name):
1045       """       """
1046       return the L{FunctionSpace<escript.FunctionSpace>} to be used for coefficient name       return the L{FunctionSpace<escript.FunctionSpace>} to be used for coefficient name
1047    
# Line 1376  class LinearPDE(object): Line 1052  class LinearPDE(object):
1052       @raise IllegalCoefficient: if name is not a coefficient of the PDE.       @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1053       """       """
1054       if self.hasCoefficient(name):       if self.hasCoefficient(name):
1055          return self.COEFFICIENTS[name].getFunctionSpace(self.getDomain())          return self.__COEFFICIENTS[name].getFunctionSpace(self.getDomain())
1056       else:       else:
1057          raise ValueError,"unknown coefficient %s requested"%name          raise ValueError,"unknown coefficient %s requested"%name
1058     def getShapeOfCoefficient(self,name):  
1059       def getShapeOfCoefficient(self,name):
1060       """       """
1061       return the shape of the coefficient name       return the shape of the coefficient name
1062    
# Line 1390  class LinearPDE(object): Line 1067  class LinearPDE(object):
1067       @raise IllegalCoefficient: if name is not a coefficient of the PDE.       @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1068       """       """
1069       if self.hasCoefficient(name):       if self.hasCoefficient(name):
1070          return self.COEFFICIENTS[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())          return self.__COEFFICIENTS[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())
1071       else:       else:
1072          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1073    
1074     def resetCoefficients(self):     def resetAllCoefficients(self):
1075       """       """
1076       resets all coefficients to there default values.       resets all coefficients to there default values.
1077       """       """
1078       for i in self.COEFFICIENTS.iterkeys():       for i in self.__COEFFICIENTS.iterkeys():
1079           self.COEFFICIENTS[i].resetValue()           self.__COEFFICIENTS[i].resetValue()
1080    
1081     def alteredCoefficient(self,name):     def alteredCoefficient(self,name):
1082       """       """
1083       announce that coefficient name has been changed       announce that coefficient name has been changed
1084    
# Line 1413  class LinearPDE(object): Line 1090  class LinearPDE(object):
1090       if self.hasCoefficient(name):       if self.hasCoefficient(name):
1091          self.trace("Coefficient %s has been altered."%name)          self.trace("Coefficient %s has been altered."%name)
1092          if not ((name=="q" or name=="r") and self.isUsingLumping()):          if not ((name=="q" or name=="r") and self.isUsingLumping()):
1093             if self.COEFFICIENTS[name].isAlteringOperator(): self.__invalidateOperator()             if self.__COEFFICIENTS[name].isAlteringOperator(): self.invalidateOperator()
1094             if self.COEFFICIENTS[name].isAlteringRightHandSide(): self.__invalidateRightHandSide()             if self.__COEFFICIENTS[name].isAlteringRightHandSide(): self.invalidateRightHandSide()
1095       else:       else:
1096          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1097    
1098     def copyConstraint(self,u):     def validSolution(self):
1099        """         """
1100        copies the constraint into u and returns u.         markes the solution as valid
1101           """
1102           self.__is_solution_valid=True
1103    
1104        @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):
1105        @type u: L{Data<escript.Data>}         """
1106        @return: the input u modified by the constraints.         indicates the PDE has to be resolved if the solution is requested
1107        @rtype: L{Data<escript.Data>}         """
1108        @warning: u is altered if it has the appropriate L{FunctionSpace<escript.FunctionSpace>}         self.trace("System will be resolved.")
1109        """         self.__is_solution_valid=False
       q=self.getCoefficientOfGeneralPDE("q")  
       r=self.getCoefficientOfGeneralPDE("r")  
       if not q.isEmpty():  
          if u.isEmpty(): u=escript.Data(0.,q.getShape(),q.getFunctionSpace())  
          if r.isEmpty():  
              r=escript.Data(0,u.getShape(),u.getFunctionSpace())  
          else:  
              r=escript.Data(r,u.getFunctionSpace())  
          u.copyWithMask(r,escript.Data(q,u.getFunctionSpace()))  
       return u  
1110    
1111     def setValue(self,**coefficients):     def isSolutionValid(self):
1112           """
1113           returns True is the solution is still valid.
1114           """
1115           return self.__is_solution_valid
1116    
1117       def validOperator(self):
1118           """
1119           markes the  operator as valid
1120           """
1121           self.__is_operator_valid=True
1122    
1123       def invalidateOperator(self):
1124           """
1125           indicates the operator has to be rebuilt next time it is used
1126           """
1127           self.trace("Operator will be rebuilt.")
1128           self.invalidateSolution()
1129           self.__is_operator_valid=False
1130    
1131       def isOperatorValid(self):
1132           """
1133           returns True is the operator is still valid.
1134           """
1135           return self.__is_operator_valid
1136    
1137       def validRightHandSide(self):
1138           """
1139           markes the right hand side as valid
1140           """
1141           self.__is_RHS_valid=True
1142    
1143       def invalidateRightHandSide(self):
1144           """
1145           indicates the right hand side has to be rebuild next time it is used
1146           """
1147           if self.isRightHandSideValid(): self.trace("Right hand side has to be rebuilt.")
1148           self.invalidateSolution()
1149           self.__is_RHS_valid=False
1150    
1151       def isRightHandSideValid(self):
1152           """
1153           returns True is the operator is still valid.
1154           """
1155           return self.__is_RHS_valid
1156    
1157       def invalidateSystem(self):
1158           """
1159           annonced that everthing has to be rebuild:
1160           """
1161           if self.isRightHandSideValid(): self.trace("System has to be rebuilt.")
1162           self.invalidateSolution()
1163           self.invalidateOperator()
1164           self.invalidateRightHandSide()
1165    
1166       def isSystemValid(self):
1167           """
1168           returns true if the system (including solution) is still vaild:
1169           """
1170           return self.isSolutionValid() and self.isOperatorValid() and self.isRightHandSideValid()
1171    
1172       def initializeSystem(self):
1173           """
1174           resets the system
1175           """
1176           self.trace("New System has been created.")
1177           self.__operator=escript.Operator()
1178           self.__righthandside=escript.Data()
1179           self.__solution=escript.Data()
1180           self.invalidateSystem()
1181    
1182       def getOperator(self):
1183         """
1184         returns the operator of the linear problem
1185    
1186         @return: the operator of the problem
1187         """
1188         return self.getSystem()[0]
1189    
1190       def getRightHandSide(self):
1191         """
1192         returns the right hand side of the linear problem
1193    
1194         @return: the right hand side of the problem
1195         @rtype: L{Data<escript.Data>}
1196         """
1197         return self.getSystem()[1]
1198    
1199       def createRightHandSide(self):
1200           """
1201           returns an instance of a new right hand side
1202           """
1203           self.trace("New right hand side is allocated.")
1204           if self.getNumEquations()>1:
1205               return escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)
1206           else:
1207               return escript.Data(0.,(),self.getFunctionSpaceForEquation(),True)
1208    
1209       def createSolution(self):
1210           """
1211           returns an instance of a new solution
1212           """
1213           self.trace("New solution is allocated.")
1214           if self.getNumSolutions()>1:
1215               return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)
1216           else:
1217               return escript.Data(0.,(),self.getFunctionSpaceForSolution(),True)
1218    
1219       def resetSolution(self):
1220           """
1221           sets the solution to zero
1222           """
1223           if self.__solution.isEmpty():
1224               self.__solution=self.createSolution()
1225           else:
1226               self.__solution.setToZero()
1227               self.trace("Solution is reset to zero.")
1228       def setSolution(self,u):
1229           """
1230           sets the solution assuming that makes the sytem valid.
1231           """
1232           self.__solution=u
1233           self.validSolution()
1234    
1235       def getCurrentSolution(self):
1236           """
1237           returns the solution in its current state
1238           """
1239           if self.__solution.isEmpty(): self.__solution=self.createSolution()
1240           return self.__solution
1241    
1242       def resetRightHandSide(self):
1243           """
1244           sets the right hand side to zero
1245           """
1246           if self.__righthandside.isEmpty():
1247               self.__righthandside=self.createRightHandSide()
1248           else:
1249               self.__righthandside.setToZero()
1250               self.trace("Right hand side is reset to zero.")
1251    
1252       def getCurrentRightHandSide(self):
1253           """
1254           returns the right hand side  in its current state
1255           """
1256           if self.__righthandside.isEmpty(): self.__righthandside=self.createRightHandSide()
1257           return self.__righthandside
1258            
1259    
1260       def resetOperator(self):
1261           """
1262           makes sure that the operator is instantiated and returns it initialized by zeros
1263           """
1264           if self.__operator.isEmpty():
1265               if self.isUsingLumping():
1266                   self.__operator=self.createSolution()
1267               else:
1268                   self.__operator=self.createOperator()
1269           else:
1270               if self.isUsingLumping():
1271                   self.__operator.setToZero()
1272               else:
1273                   self.__operator.resetValues()
1274               self.trace("Operator reset to zero")
1275    
1276       def getCurrentOperator(self):
1277           """
1278           returns the operator in its current state
1279           """
1280           if self.__operator.isEmpty(): self.resetOperator()
1281           return self.__operator
1282    
1283       def setValue(self,**coefficients):
1284        """        """
1285        sets new values to coefficients        sets new values to coefficients
1286    
       @param coefficients: new values assigned to coefficients  
       @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.  
1287        @raise IllegalCoefficient: if an unknown coefficient keyword is used.        @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1288        """        """
1289        # check if the coefficients are  legal:        # check if the coefficients are  legal:
# Line 1489  class LinearPDE(object): Line 1301  class LinearPDE(object):
1301                  s=numarray.array(d).shape                  s=numarray.array(d).shape
1302              if s!=None:              if s!=None:
1303                  # get number of equations and number of unknowns:                  # get number of equations and number of unknowns:
1304                  res=self.COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(self.getDomain(),s)                  res=self.__COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(self.getDomain(),s)
1305                  if res==None:                  if res==None:
1306                      raise IllegalCoefficientValue,"Illegal shape %s of coefficient %s"%(s,i)                      raise IllegalCoefficientValue,"Illegal shape %s of coefficient %s"%(s,i)
1307                  else:                  else:
# Line 1500  class LinearPDE(object): Line 1312  class LinearPDE(object):
1312        # now we check the shape of the coefficient if numEquations and numSolutions are set:        # now we check the shape of the coefficient if numEquations and numSolutions are set:
1313        for i,d in coefficients.iteritems():        for i,d in coefficients.iteritems():
1314          try:          try:
1315             self.COEFFICIENTS[i].setValue(self.getDomain(),self.getNumEquations(),self.getNumSolutions(),self.reduceEquationOrder(),self.reduceSolutionOrder(),d)             self.__COEFFICIENTS[i].setValue(self.getDomain(),
1316                                             self.getNumEquations(),self.getNumSolutions(),
1317                                             self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
1318               self.alteredCoefficient(i)
1319            except IllegalCoefficientFunctionSpace,m:
1320                # if the function space is wrong then we try the reduced version:
1321                i_red=i+"_reduced"
1322                if (not i_red in coefficients.keys()) and i_red in self.__COEFFICIENTS.keys():
1323                    try:
1324                        self.__COEFFICIENTS[i_red].setValue(self.getDomain(),
1325                                                          self.getNumEquations(),self.getNumSolutions(),
1326                                                          self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
1327                        self.alteredCoefficient(i_red)
1328                    except IllegalCoefficientValue,m:
1329                        raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
1330                    except IllegalCoefficientFunctionSpace,m:
1331                        raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m))
1332                else:
1333                    raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m))
1334          except IllegalCoefficientValue,m:          except IllegalCoefficientValue,m:
1335             raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))             raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
         self.alteredCoefficient(i)  
   
1336        self.__altered_coefficients=True        self.__altered_coefficients=True
       # check if the systrem is inhomogeneous:  
       if len(coefficients)>0 and not self.isUsingLumping():  
          q=self.getCoefficientOfGeneralPDE("q")  
          r=self.getCoefficientOfGeneralPDE("r")  
          homogeneous_constraint=True  
          if not q.isEmpty() and not r.isEmpty():  
              if util.Lsup(q*r)>=1.e-13*util.Lsup(r):  
                self.trace("Inhomogeneous constraint detected.")  
                self.__invalidateSystem()  
1337    
1338     def getSystem(self):     # ====================================================================================
1339       # methods that are typically overwritten when implementing a particular linear problem
1340       # ====================================================================================
1341       def getRequiredSystemType(self):
1342          """
1343          returns the system type which needs to be used by the current set up.
1344    
1345          @note: Typically this method is overwritten when implementing a particular linear problem.
1346          """
1347          return 0
1348    
1349       def createOperator(self):
1350           """
1351           returns an instance of a new operator
1352    
1353           @note: This method is overwritten when implementing a particular linear problem.
1354           """
1355           return escript.Operator()
1356    
1357       def checkSymmetry(self,verbose=True):
1358          """
1359          test the PDE for symmetry.
1360    
1361          @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.
1362          @type verbose: C{bool}
1363          @return:  True if the problem is symmetric.
1364          @rtype: C{bool}
1365          @note: Typically this method is overwritten when implementing a particular linear problem.
1366          """
1367          out=True
1368          return out
1369    
1370       def getSolution(self,**options):
1371           """
1372           returns the solution of the problem
1373           @return: the solution
1374           @rtype: L{Data<escript.Data>}
1375    
1376           @note: This method is overwritten when implementing a particular linear problem.
1377           """
1378           return self.getCurrentSolution()
1379    
1380       def getSystem(self):
1381           """
1382           return the operator and right hand side of the PDE
1383    
1384           @return: the discrete version of the PDE
1385           @rtype: C{tuple} of L{Operator,<escript.Operator>} and L{Data<escript.Data>}.
1386    
1387           @note: This method is overwritten when implementing a particular linear problem.
1388           """
1389           return (self.getCurrentOperator(), self.getCurrentRightHandSide())
1390    #=====================
1391    
1392    class LinearPDE(LinearProblem):
1393       """
1394       This class is used to define a general linear, steady, second order PDE
1395       for an unknown function M{u} on a given domain defined through a L{Domain<escript.Domain>} object.
1396    
1397       For a single PDE with a solution with a single component the linear PDE is defined in the following form:
1398    
1399       M{-(grad(A[j,l]+A_reduced[j,l])*grad(u)[l]+(B[j]+B_reduced[j])u)[j]+(C[l]+C_reduced[l])*grad(u)[l]+(D+D_reduced)=-grad(X+X_reduced)[j,j]+(Y+Y_reduced)}
1400    
1401    
1402       where M{grad(F)} denotes the spatial derivative of M{F}. Einstein's summation convention,
1403       ie. summation over indexes appearing twice in a term of a sum is performed, is used.
1404       The coefficients M{A}, M{B}, M{C}, M{D}, M{X} and M{Y} have to be specified through L{Data<escript.Data>} objects in the
1405       L{Function<escript.Function>} and the coefficients M{A_reduced}, M{B_reduced}, M{C_reduced}, M{D_reduced}, M{X_reduced} and M{Y_reduced} have to be specified through L{Data<escript.Data>} objects in the
1406       L{ReducedFunction<escript.ReducedFunction>}. It is also allowd to use objects that can be converted into
1407       such L{Data<escript.Data>} objects. M{A} and M{A_reduced} are rank two, M{B}, M{C}, M{X},
1408       M{B_reduced}, M{C_reduced} and M{X_reduced} are rank one and M{D}, M{D_reduced} M{Y} and M{Y_reduced} are scalar.
1409    
1410       The following natural boundary conditions are considered:
1411    
1412       M{n[j]*((A[i,j]+A_reduced[i,j])*grad(u)[l]+(B+B_reduced)[j]*u)+(d+d_reduced)*u=n[j]*(X[j]+X_reduced[j])+y}
1413    
1414       where M{n} is the outer normal field. Notice that the coefficients M{A}, M{A_reduced}, M{B}, M{B_reduced}, M{X} and M{X_reduced} are defined in the PDE. The coefficients M{d} and M{y} and are each a scalar in the L{FunctionOnBoundary<escript.FunctionOnBoundary>} and the coefficients M{d_reduced} and M{y_reduced} and are each a scalar in the L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
1415    
1416    
1417       Constraints for the solution prescribing the value of the solution at certain locations in the domain. They have the form
1418    
1419       M{u=r}  where M{q>0}
1420    
1421       M{r} and M{q} are each scalar where M{q} is the characteristic function defining where the constraint is applied.
1422       The constraints override any other condition set by the PDE or the boundary condition.
1423    
1424       The PDE is symmetrical if
1425    
1426       M{A[i,j]=A[j,i]}  and M{B[j]=C[j]} and M{A_reduced[i,j]=A_reduced[j,i]}  and M{B_reduced[j]=C_reduced[j]}
1427    
1428       For a system of PDEs and a solution with several components the PDE has the form
1429    
1430       M{-grad((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k])[j]+(C[i,k,l]+C_reduced[i,k,l])*grad(u[k])[l]+(D[i,k]+D_reduced[i,k]*u[k] =-grad(X[i,j]+X_reduced[i,j])[j]+Y[i]+Y_reduced[i] }
1431    
1432       M{A} and M{A_reduced} are of rank four, M{B}, M{B_reduced}, M{C} and M{C_reduced} are each of rank three, M{D}, M{D_reduced}, M{X_reduced} and M{X} are each a rank two and M{Y} and M{Y_reduced} are of rank one.
1433       The natural boundary conditions take the form:
1434    
1435       M{n[j]*((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k])+(d[i,k]+d_reduced[i,k])*u[k]=n[j]*(X[i,j]+X_reduced[i,j])+y[i]+y_reduced[i]}
1436    
1437    
1438       The coefficient M{d} is a rank two and M{y} is a rank one both in the L{FunctionOnBoundary<escript.FunctionOnBoundary>}. Constraints take the form and the coefficients M{d_reduced} is a rank two and M{y_reduced} is a rank one both in the L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
1439    
1440       Constraints take the form
1441    
1442       M{u[i]=r[i]}  where  M{q[i]>0}
1443    
1444       M{r} and M{q} are each rank one. Notice that at some locations not necessarily all components must have a constraint.
1445    
1446       The system of PDEs is symmetrical if
1447    
1448            - M{A[i,j,k,l]=A[k,l,i,j]}
1449            - M{A_reduced[i,j,k,l]=A_reduced[k,l,i,j]}
1450            - M{B[i,j,k]=C[k,i,j]}
1451            - M{B_reduced[i,j,k]=C_reduced[k,i,j]}
1452            - M{D[i,k]=D[i,k]}
1453            - M{D_reduced[i,k]=D_reduced[i,k]}
1454            - M{d[i,k]=d[k,i]}
1455            - M{d_reduced[i,k]=d_reduced[k,i]}
1456    
1457       L{LinearPDE} also supports solution discontinuities over a contact region in the domain. To specify the conditions across the
1458       discontinuity we are using the generalised flux M{J} which is in the case of a systems of PDEs and several components of the solution
1459       defined as
1460    
1461       M{J[i,j]=(A[i,j,k,l]+A_reduced[[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k]-X[i,j]-X_reduced[i,j]}
1462    
1463       For the case of single solution component and single PDE M{J} is defined
1464    
1465       M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[j]+(B[i]+B_reduced[i])*u-X[i]-X_reduced[i]}
1466    
1467       In the context of discontinuities M{n} denotes the normal on the discontinuity pointing from side 0 towards side 1
1468       calculated from L{getNormal<escript.FunctionSpace.getNormal>} of L{FunctionOnContactZero<escript.FunctionOnContactZero>}. For a system of PDEs
1469       the contact condition takes the form
1470    
1471       M{n[j]*J0[i,j]=n[j]*J1[i,j]=(y_contact[i]+y_contact_reduced[i])- (d_contact[i,k]+d_contact_reduced[i,k])*jump(u)[k]}
1472    
1473       where M{J0} and M{J1} are the fluxes on side 0 and side 1 of the discontinuity, respectively. M{jump(u)}, which is the difference
1474       of the solution at side 1 and at side 0, denotes the jump of M{u} across discontinuity along the normal calculated by
1475       L{jump<util.jump>}.
1476       The coefficient M{d_contact} is a rank two and M{y_contact} is a rank one both in the L{FunctionOnContactZero<escript.FunctionOnContactZero>} or L{FunctionOnContactOne<escript.FunctionOnContactOne>}.
1477       The coefficient M{d_contact_reduced} is a rank two and M{y_contact_reduced} is a rank one both in the L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.
1478       In case of a single PDE and a single component solution the contact condition takes the form
1479    
1480       M{n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)}
1481    
1482       In this case the coefficient M{d_contact} and M{y_contact} are each scalar both in the L{FunctionOnContactZero<escript.FunctionOnContactZero>} or L{FunctionOnContactOne<escript.FunctionOnContactOne>} and the coefficient M{d_contact_reduced} and M{y_contact_reduced} are each scalar both in the L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}
1483    
1484       typical usage:
1485    
1486          p=LinearPDE(dom)
1487          p.setValue(A=kronecker(dom),D=1,Y=0.5)
1488          u=p.getSolution()
1489    
1490       """
1491    
1492    
1493       def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):
1494         """
1495         initializes a new linear PDE
1496    
1497         @param domain: domain of the PDE
1498         @type domain: L{Domain<escript.Domain>}
1499         @param numEquations: number of equations. If numEquations==None the number of equations
1500                              is extracted from the PDE coefficients.
1501         @param numSolutions: number of solution components. If  numSolutions==None the number of solution components
1502                              is extracted from the PDE coefficients.
1503         @param debug: if True debug informations are printed.
1504    
1505         """
1506         super(LinearPDE, self).__init__(domain,numEquations,numSolutions,debug)
1507         #
1508         #   the coefficients of the PDE:
1509         #
1510         self.introduceCoefficients(
1511           A=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
1512           B=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1513           C=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
1514           D=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1515           X=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
1516           Y=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1517           d=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1518           y=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1519           d_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1520           y_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1521           A_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
1522           B_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1523           C_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
1524           D_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1525           X_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
1526           Y_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1527           d_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1528           y_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1529           d_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1530           y_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1531           r=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.RIGHTHANDSIDE),
1532           q=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.BOTH) )
1533    
1534       def __str__(self):
1535         """
1536         returns string representation of the PDE
1537    
1538         @return: a simple representation of the PDE
1539         @rtype: C{str}
1540         """
1541         return "<LinearPDE %d>"%id(self)
1542    
1543       def getRequiredSystemType(self):
1544          """
1545          returns the system type which needs to be used by the current set up.
1546          """
1547          return self.getDomain().getSystemMatrixTypeId(self.getSolverMethod()[0],self.getSolverMethod()[1],self.getSolverPackage(),self.isSymmetric())
1548    
1549       def checkSymmetry(self,verbose=True):
1550          """
1551          test the PDE for symmetry.
1552    
1553          @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.
1554          @type verbose: C{bool}
1555          @return:  True if the PDE is symmetric.
1556          @rtype: L{bool}
1557          @note: This is a very expensive operation. It should be used for degugging only! The symmetry flag is not altered.
1558          """
1559          out=True
1560          out=out and self.checkSymmetricTensor("A", verbose)
1561          out=out and self.checkSymmetricTensor("A_reduced", verbose)
1562          out=out and self.checkReciprocalSymmetry("B","C", verbose)
1563          out=out and self.checkReciprocalSymmetry("B_reduced","C_reduced", verbose)
1564          out=out and self.checkSymmetricTensor("D", verbose)
1565          out=out and self.checkSymmetricTensor("D_reduced", verbose)
1566          out=out and self.checkSymmetricTensor("d", verbose)
1567          out=out and self.checkSymmetricTensor("d_reduced", verbose)
1568          out=out and self.checkSymmetricTensor("d_contact", verbose)
1569          out=out and self.checkSymmetricTensor("d_contact_reduced", verbose)
1570          return out
1571    
1572       def createOperator(self):
1573           """
1574           returns an instance of a new operator
1575           """
1576           self.trace("New operator is allocated.")
1577           return self.getDomain().newOperator( \
1578                               self.getNumEquations(), \
1579                               self.getFunctionSpaceForEquation(), \
1580                               self.getNumSolutions(), \
1581                               self.getFunctionSpaceForSolution(), \
1582                               self.getSystemType())
1583    
1584       def getSolution(self,**options):
1585           """
1586           returns the solution of the PDE
1587    
1588           @return: the solution
1589           @rtype: L{Data<escript.Data>}
1590           @param options: solver options
1591           @keyword verbose: True to get some information during PDE solution
1592           @type verbose: C{bool}
1593           @keyword reordering: reordering scheme to be used during elimination. Allowed values are
1594                                L{NO_REORDERING}, L{MINIMUM_FILL_IN}, L{NESTED_DISSECTION}
1595           @keyword iter_max: maximum number of iteration steps allowed.
1596           @keyword drop_tolerance: threshold for drupping in L{ILUT}
1597           @keyword drop_storage: maximum of allowed memory in L{ILUT}
1598           @keyword truncation: maximum number of residuals in L{GMRES}
1599           @keyword restart: restart cycle length in L{GMRES}
1600           """
1601           if not self.isSolutionValid():
1602              mat,f=self.getSystem()
1603              if self.isUsingLumping():
1604                 self.setSolution(f*1/mat)
1605              else:
1606                 options[self.TOLERANCE_KEY]=self.getTolerance()
1607                 options[self.METHOD_KEY]=self.getSolverMethod()[0]
1608                 options[self.PRECONDITIONER_KEY]=self.getSolverMethod()[1]
1609                 options[self.PACKAGE_KEY]=self.getSolverPackage()
1610                 options[self.SYMMETRY_KEY]=self.isSymmetric()
1611                 self.trace("PDE is resolved.")
1612                 self.trace("solver options: %s"%str(options))
1613                 self.setSolution(mat.solve(f,options))
1614           return self.getCurrentSolution()
1615    
1616       def getSystem(self):
1617         """         """
1618         return the operator and right hand side of the PDE         return the operator and right hand side of the PDE
1619    
1620         @return: the discrete version of the PDE         @return: the discrete version of the PDE
1621         @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>}.
1622         """         """
1623         if not self.__operator_is_Valid or not self.__righthandside_isValid:         if not self.isOperatorValid() or not self.isRightHandSideValid():
1624            if self.isUsingLumping():            if self.isUsingLumping():
1625                if not self.__operator_is_Valid:                if not self.isOperatorValid():
1626                   if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():                   if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():
1627                        raise TypeError,"Lumped matrix requires same order for equations and unknowns"                        raise TypeError,"Lumped matrix requires same order for equations and unknowns"
1628                   if not self.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():  
1629                        raise ValueError,"coefficient A in lumped matrix may not be present."                        raise ValueError,"coefficient A in lumped matrix may not be present."
1630                   if not self.getCoefficientOfGeneralPDE("C").isEmpty():                   if not self.getCoefficient("B").isEmpty():
1631                        raise ValueError,"coefficient A in lumped matrix may not be present."                        raise ValueError,"coefficient B in lumped matrix may not be present."
1632                   D=self.getCoefficientOfGeneralPDE("D")                   if not self.getCoefficient("C").isEmpty():
1633                          raise ValueError,"coefficient C in lumped matrix may not be present."
1634                     if not self.getCoefficient("d_contact").isEmpty():
1635                          raise ValueError,"coefficient d_contact in lumped matrix may not be present."
1636                     if not self.getCoefficient("A_reduced").isEmpty():
1637                          raise ValueError,"coefficient A_reduced in lumped matrix may not be present."
1638                     if not self.getCoefficient("B_reduced").isEmpty():
1639                          raise ValueError,"coefficient B_reduced in lumped matrix may not be present."
1640                     if not self.getCoefficient("C_reduced").isEmpty():
1641                          raise ValueError,"coefficient C_reduced in lumped matrix may not be present."
1642                     if not self.getCoefficient("d_contact_reduced").isEmpty():
1643                          raise ValueError,"coefficient d_contact_reduced in lumped matrix may not be present."
1644                     D=self.getCoefficient("D")
1645                     d=self.getCoefficient("d")
1646                     D_reduced=self.getCoefficient("D_reduced")
1647                     d_reduced=self.getCoefficient("d_reduced")
1648                   if not D.isEmpty():                   if not D.isEmpty():
1649                       if self.getNumSolutions()>1:                       if self.getNumSolutions()>1:
1650                          D_times_e=util.matrixmult(D,numarray.ones((self.getNumSolutions(),)))                          D_times_e=util.matrix_mult(D,numarray.ones((self.getNumSolutions(),)))
1651                       else:                       else:
1652                          D_times_e=D                          D_times_e=D
1653                   else:                   else:
1654                      D_times_e=escript.Data()                      D_times_e=escript.Data()
                  d=self.getCoefficientOfGeneralPDE("d")  
1655                   if not d.isEmpty():                   if not d.isEmpty():
1656                       if self.getNumSolutions()>1:                       if self.getNumSolutions()>1:
1657                          d_times_e=util.matrixmult(d,numarray.ones((self.getNumSolutions(),)))                          d_times_e=util.matrix_mult(d,numarray.ones((self.getNumSolutions(),)))
1658                       else:                       else:
1659                          d_times_e=d                          d_times_e=d
1660                   else:                   else:
1661                      d_times_e=escript.Data()                      d_times_e=escript.Data()
1662                   d_contact=self.getCoefficientOfGeneralPDE("d_contact")        
1663                   if not d_contact.isEmpty():                   if not D_reduced.isEmpty():
1664                       if self.getNumSolutions()>1:                       if self.getNumSolutions()>1:
1665                          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(),)))
1666                       else:                       else:
1667                          d_contact_times_e=d_contact                          D_reduced_times_e=D_reduced
1668                   else:                   else:
1669                      d_contact_times_e=escript.Data()                      D_reduced_times_e=escript.Data()
1670                         if not d_reduced.isEmpty():
1671                   self.__operator=self.__getNewRightHandSide()                       if self.getNumSolutions()>1:
1672                   self.getDomain().addPDEToRHS(self.__operator, \                          d_reduced_times_e=util.matrix_mult(d_reduced,numarray.ones((self.getNumSolutions(),)))
1673                                                escript.Data(), \                       else:
1674                                                D_times_e, \                          d_reduced_times_e=d_reduced
1675                                                d_times_e,\                   else:
1676                                                d_contact_times_e)                      d_reduced_times_e=escript.Data()
1677                   self.__operator=1./self.__operator  
1678                     self.resetOperator()
1679                     operator=self.getCurrentOperator()
1680                     if False and hasattr(self.getDomain(), "addPDEToLumpedSystem") :
1681                        self.getDomain().addPDEToLumpedSystem(operator, D_times_e, d_times_e)
1682                        self.getDomain().addPDEToLumpedSystem(operator, D_reduced_times_e, d_reduced_times_e)
1683                     else:
1684                        self.getDomain().addPDEToRHS(operator, \
1685                                                     escript.Data(), \
1686                                                     D_times_e, \
1687                                                     d_times_e,\
1688                                                     escript.Data())
1689                        self.getDomain().addPDEToRHS(operator, \
1690                                                     escript.Data(), \
1691                                                     D_reduced_times_e, \
1692                                                     d_reduced_times_e,\
1693                                                     escript.Data())
1694                   self.trace("New lumped operator has been built.")                   self.trace("New lumped operator has been built.")
1695                   self.__operator_is_Valid=True                if not self.isRightHandSideValid():
1696                if not self.__righthandside_isValid:                   self.resetRightHandSide()
1697                   self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \                   righthandside=self.getCurrentRightHandSide()
1698                                 self.getCoefficientOfGeneralPDE("X"), \                   self.getDomain().addPDEToRHS(righthandside, \
1699                                 self.getCoefficientOfGeneralPDE("Y"),\                                 self.getCoefficient("X"), \
1700                                 self.getCoefficientOfGeneralPDE("y"),\                                 self.getCoefficient("Y"),\
1701                                 self.getCoefficientOfGeneralPDE("y_contact"))                                 self.getCoefficient("y"),\
1702                                   self.getCoefficient("y_contact"))
1703                     self.getDomain().addPDEToRHS(righthandside, \
1704                                   self.getCoefficient("X_reduced"), \
1705                                   self.getCoefficient("Y_reduced"),\
1706                                   self.getCoefficient("y_reduced"),\
1707                                   self.getCoefficient("y_contact_reduced"))
1708                   self.trace("New right hand side as been built.")                   self.trace("New right hand side as been built.")
1709                   self.__righthandside_isValid=True                   self.validRightHandSide()
1710                  self.insertConstraint(rhs_only=False)
1711                  self.validOperator()
1712            else:            else:
1713               if not self.__operator_is_Valid and not self.__righthandside_isValid:               if not self.isOperatorValid() and not self.isRightHandSideValid():
1714                   self.getDomain().addPDEToSystem(self.__makeFreshOperator(),self.__makeFreshRightHandSide(), \                   self.resetRightHandSide()
1715                                 self.getCoefficientOfGeneralPDE("A"), \                   righthandside=self.getCurrentRightHandSide()
1716                                 self.getCoefficientOfGeneralPDE("B"), \                   self.resetOperator()
1717                                 self.getCoefficientOfGeneralPDE("C"), \                   operator=self.getCurrentOperator()
1718                                 self.getCoefficientOfGeneralPDE("D"), \                   self.getDomain().addPDEToSystem(operator,righthandside, \
1719                                 self.getCoefficientOfGeneralPDE("X"), \                                 self.getCoefficient("A"), \
1720                                 self.getCoefficientOfGeneralPDE("Y"), \                                 self.getCoefficient("B"), \
1721                                 self.getCoefficientOfGeneralPDE("d"), \                                 self.getCoefficient("C"), \
1722                                 self.getCoefficientOfGeneralPDE("y"), \                                 self.getCoefficient("D"), \
1723                                 self.getCoefficientOfGeneralPDE("d_contact"), \                                 self.getCoefficient("X"), \
1724                                 self.getCoefficientOfGeneralPDE("y_contact"))                                 self.getCoefficient("Y"), \
1725                   self.__applyConstraint()                                 self.getCoefficient("d"), \
1726                   self.__righthandside=self.copyConstraint(self.__righthandside)                                 self.getCoefficient("y"), \
1727                                   self.getCoefficient("d_contact"), \
1728                                   self.getCoefficient("y_contact"))
1729                     self.getDomain().addPDEToSystem(operator,righthandside, \
1730                                   self.getCoefficient("A_reduced"), \
1731                                   self.getCoefficient("B_reduced"), \
1732                                   self.getCoefficient("C_reduced"), \
1733                                   self.getCoefficient("D_reduced"), \
1734                                   self.getCoefficient("X_reduced"), \
1735                                   self.getCoefficient("Y_reduced"), \
1736                                   self.getCoefficient("d_reduced"), \
1737                                   self.getCoefficient("y_reduced"), \
1738                                   self.getCoefficient("d_contact_reduced"), \
1739                                   self.getCoefficient("y_contact_reduced"))
1740                     self.insertConstraint(rhs_only=False)
1741                   self.trace("New system has been built.")                   self.trace("New system has been built.")
1742                   self.__operator_is_Valid=True                   self.validOperator()
1743                   self.__righthandside_isValid=True                   self.validRightHandSide()
1744               elif not self.__righthandside_isValid:               elif not self.isRightHandSideValid():
1745                   self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \                   self.resetRightHandSide()
1746                                 self.getCoefficientOfGeneralPDE("X"), \                   righthandside=self.getCurrentRightHandSide()
1747                                 self.getCoefficientOfGeneralPDE("Y"),\                   self.getDomain().addPDEToRHS(righthandside,
1748                                 self.getCoefficientOfGeneralPDE("y"),\                                 self.getCoefficient("X"), \
1749                                 self.getCoefficientOfGeneralPDE("y_contact"))                                 self.getCoefficient("Y"),\
1750                   self.__righthandside=self.copyConstraint(self.__righthandside)                                 self.getCoefficient("y"),\
1751                                   self.getCoefficient("y_contact"))
1752                     self.getDomain().addPDEToRHS(righthandside,
1753                                   self.getCoefficient("X_reduced"), \
1754                                   self.getCoefficient("Y_reduced"),\
1755                                   self.getCoefficient("y_reduced"),\
1756                                   self.getCoefficient("y_contact_reduced"))
1757                     self.insertConstraint(rhs_only=True)
1758                   self.trace("New right hand side has been built.")                   self.trace("New right hand side has been built.")
1759                   self.__righthandside_isValid=True                   self.validRightHandSide()
1760               elif not self.__operator_is_Valid:               elif not self.isOperatorValid():
1761                   self.getDomain().addPDEToSystem(self.__makeFreshOperator(),escript.Data(), \                   self.resetOperator()
1762                              self.getCoefficientOfGeneralPDE("A"), \                   operator=self.getCurrentOperator()
1763                              self.getCoefficientOfGeneralPDE("B"), \                   self.getDomain().addPDEToSystem(operator,escript.Data(), \
1764                              self.getCoefficientOfGeneralPDE("C"), \                              self.getCoefficient("A"), \
1765                              self.getCoefficientOfGeneralPDE("D"), \                              self.getCoefficient("B"), \
1766                                self.getCoefficient("C"), \
1767                                self.getCoefficient("D"), \
1768                                escript.Data(), \
1769                                escript.Data(), \
1770                                self.getCoefficient("d"), \
1771                                escript.Data(),\
1772                                self.getCoefficient("d_contact"), \
1773                                escript.Data())
1774                     self.getDomain().addPDEToSystem(operator,escript.Data(), \
1775                                self.getCoefficient("A_reduced"), \
1776                                self.getCoefficient("B_reduced"), \
1777                                self.getCoefficient("C_reduced"), \
1778                                self.getCoefficient("D_reduced"), \
1779                              escript.Data(), \                              escript.Data(), \
1780                              escript.Data(), \                              escript.Data(), \
1781                              self.getCoefficientOfGeneralPDE("d"), \                              self.getCoefficient("d_reduced"), \
1782                              escript.Data(),\                              escript.Data(),\
1783                              self.getCoefficientOfGeneralPDE("d_contact"), \                              self.getCoefficient("d_contact_reduced"), \
1784                              escript.Data())                              escript.Data())
1785                   self.__applyConstraint()                   self.insertConstraint(rhs_only=False)
1786                   self.trace("New operator has been built.")                   self.trace("New operator has been built.")
1787                   self.__operator_is_Valid=True                   self.validOperator()
1788         return (self.__operator,self.__righthandside)         return (self.getCurrentOperator(), self.getCurrentRightHandSide())
1789    
1790       def insertConstraint(self, rhs_only=False):
1791          """
1792          applies the constraints defined by q and r to the PDE
1793          
1794          @param rhs_only: if True only the right hand side is altered by the constraint.
1795          @type rhs_only: C{bool}
1796          """
1797          q=self.getCoefficient("q")
1798          r=self.getCoefficient("r")
1799          righthandside=self.getCurrentRightHandSide()
1800          operator=self.getCurrentOperator()
1801    
1802          if not q.isEmpty():
1803             if r.isEmpty():
1804                r_s=self.createSolution()
1805             else:
1806                r_s=r
1807             if not rhs_only and not operator.isEmpty():
1808                 if self.isUsingLumping():
1809                     operator.copyWithMask(escript.Data(1.,q.getShape(),q.getFunctionSpace()),q)
1810                 else:
1811                     row_q=escript.Data(q,self.getFunctionSpaceForEquation())
1812                     col_q=escript.Data(q,self.getFunctionSpaceForSolution())
1813                     u=self.createSolution()
1814                     u.copyWithMask(r_s,col_q)
1815                     righthandside-=operator*u
1816                     operator.nullifyRowsAndCols(row_q,col_q,1.)
1817             righthandside.copyWithMask(r_s,q)
1818    
1819       def setValue(self,**coefficients):
1820          """
1821          sets new values to coefficients
1822    
1823          @param coefficients: new values assigned to coefficients
1824          @keyword A: value for coefficient A.
1825          @type A: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1826          @keyword A_reduced: value for coefficient A_reduced.
1827          @type A_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1828          @keyword B: value for coefficient B
1829          @type B: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1830          @keyword B_reduced: value for coefficient B_reduced
1831          @type B_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1832          @keyword C: value for coefficient C
1833          @type C: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1834          @keyword C_reduced: value for coefficient C_reduced
1835          @type C_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1836          @keyword D: value for coefficient D
1837          @type D: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1838          @keyword D_reduced: value for coefficient D_reduced
1839          @type D_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1840          @keyword X: value for coefficient X
1841          @type X: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1842          @keyword X_reduced: value for coefficient X_reduced
1843          @type X_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1844          @keyword Y: value for coefficient Y
1845          @type Y: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1846          @keyword Y_reduced: value for coefficient Y_reduced
1847          @type Y_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.Function>}.
1848          @keyword d: value for coefficient d
1849          @type d: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1850          @keyword d_reduced: value for coefficient d_reduced
1851          @type d_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
1852          @keyword y: value for coefficient y
1853          @type y: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1854          @keyword d_contact: value for coefficient d_contact
1855          @type d_contact: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnContactOne<escript.FunctionOnContactOne>} or  L{FunctionOnContactZero<escript.FunctionOnContactZero>}.
1856          @keyword d_contact_reduced: value for coefficient d_contact_reduced
1857          @type d_contact_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>} or  L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>}.
1858          @keyword y_contact: value for coefficient y_contact
1859          @type y_contact: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnContactOne<escript.FunctionOnContactOne>} or  L{FunctionOnContactZero<escript.FunctionOnContactZero>}.
1860          @keyword y_contact_reduced: value for coefficient y_contact_reduced
1861          @type y_contact_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunctionOnContactOne<escript.FunctionOnContactOne>} or L{ReducedFunctionOnContactZero<escript.FunctionOnContactZero>}.
1862          @keyword r: values prescribed to the solution at the locations of constraints
1863          @type r: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1864                   depending of reduced order is used for the solution.
1865          @keyword q: mask for location of constraints
1866          @type q: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1867                   depending of reduced order is used for the representation of the equation.
1868          @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1869          """
1870          super(LinearPDE,self).setValue(**coefficients)
1871          # check if the systrem is inhomogeneous:
1872          if len(coefficients)>0 and not self.isUsingLumping():
1873             q=self.getCoefficient("q")
1874             r=self.getCoefficient("r")
1875             if not q.isEmpty() and not r.isEmpty():
1876                 if util.Lsup(q*r)>0.:
1877                   self.trace("Inhomogeneous constraint detected.")
1878                   self.invalidateSystem()
1879    
1880    
1881       def getResidual(self,u=None):
1882         """
1883         return the residual of u or the current solution if u is not present.
1884    
1885         @param u: argument in the residual calculation. It must be representable in C{elf.getFunctionSpaceForSolution()}. If u is not present or equals L{None}
1886                   the current solution is used.
1887         @type u: L{Data<escript.Data>} or None
1888         @return: residual of u
1889         @rtype: L{Data<escript.Data>}
1890         """
1891         if u==None:
1892            return self.getOperator()*self.getSolution()-self.getRightHandSide()
1893         else:
1894            return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())-self.getRightHandSide()
1895    
1896       def getFlux(self,u=None):
1897         """
1898         returns the flux M{J} for a given M{u}
1899    
1900         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]}
1901    
1902         or
1903    
1904         M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[l]+(B[j]+B_reduced[j])u-X[j]-X_reduced[j]}
1905    
1906         @param u: argument in the flux. If u is not present or equals L{None} the current solution is used.
1907         @type u: L{Data<escript.Data>} or None
1908         @return: flux
1909         @rtype: L{Data<escript.Data>}
1910         """
1911         if u==None: u=self.getSolution()
1912         return util.tensormult(self.getCoefficient("A"),util.grad(u,Funtion(self.getDomain))) \
1913               +util.matrixmult(self.getCoefficient("B"),u) \
1914               -util.self.getCoefficient("X") \
1915               +util.tensormult(self.getCoefficient("A_reduced"),util.grad(u,ReducedFuntion(self.getDomain))) \
1916               +util.matrixmult(self.getCoefficient("B_reduced"),u) \
1917               -util.self.getCoefficient("X_reduced")
1918    
1919    
1920  class Poisson(LinearPDE):  class Poisson(LinearPDE):
# Line 1647  class Poisson(LinearPDE): Line 1943  class Poisson(LinearPDE):
1943    
1944       """       """
1945       super(Poisson, self).__init__(domain,1,1,debug)       super(Poisson, self).__init__(domain,1,1,debug)
1946       self.COEFFICIENTS={"f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),       self.introduceCoefficients(
1947                            "q": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}                          f=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1948                            f_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE))
1949       self.setSymmetryOn()       self.setSymmetryOn()
1950    
1951     def setValue(self,**coefficients):     def setValue(self,**coefficients):
# Line 1665  class Poisson(LinearPDE): Line 1962  class Poisson(LinearPDE):
1962       """       """
1963       super(Poisson, self).setValue(**coefficients)       super(Poisson, self).setValue(**coefficients)
1964    
1965     def getCoefficientOfGeneralPDE(self,name):  
1966       def getCoefficient(self,name):
1967       """       """
1968       return the value of the coefficient name of the general PDE       return the value of the coefficient name of the general PDE
1969       @param name: name of the coefficient requested.       @param name: name of the coefficient requested.
1970       @type name: C{string}       @type name: C{string}
1971       @return: the value of the coefficient  name       @return: the value of the coefficient  name
1972       @rtype: L{Data<escript.Data>}       @rtype: L{Data<escript.Data>}
1973       @raise IllegalCoefficient: if name is not one of coefficients       @raise IllegalCoefficient: invalid coefficient name
                   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}.  
1974       @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.       @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
1975       """       """
1976       if name == "A" :       if name == "A" :
1977           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()  
1978       elif name == "Y" :       elif name == "Y" :
1979           return self.getCoefficient("f")           return self.getCoefficient("f")
1980       elif name == "d" :       elif name == "Y_reduced" :
1981           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")  
1982       else:       else:
1983          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name          return super(Poisson, self).getCoefficient(name)
1984    
1985  class Helmholtz(LinearPDE):  class Helmholtz(LinearPDE):
1986     """     """
# Line 1729  class Helmholtz(LinearPDE): Line 2008  class Helmholtz(LinearPDE):
2008    
2009       """       """
2010       super(Helmholtz, self).__init__(domain,1,1,debug)       super(Helmholtz, self).__init__(domain,1,1,debug)
2011       self.COEFFICIENTS={"omega": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),       self.introduceCoefficients(
2012                          "k": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),                          omega=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.OPERATOR),
2013                          "f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),                          k=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.OPERATOR),
2014                          "alpha": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),                          f=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2015                          "g": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),                          f_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2016                          "r": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH),                          alpha=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.OPERATOR),
2017                          "q": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}                          g=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2018                            g_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE))
2019       self.setSymmetryOn()       self.setSymmetryOn()
2020    
2021     def setValue(self,**coefficients):     def setValue(self,**coefficients):
2022       """       """
2023       sets new values to coefficients       sets new values to coefficients
2024    
# Line 1763  class Helmholtz(LinearPDE): Line 2043  class Helmholtz(LinearPDE):
2043       """       """
2044       super(Helmholtz, self).setValue(**coefficients)       super(Helmholtz, self).setValue(**coefficients)
2045    
2046     def getCoefficientOfGeneralPDE(self,name):     def getCoefficient(self,name):
2047       """       """
2048       return the value of the coefficient name of the general PDE       return the value of the coefficient name of the general PDE
2049    
# Line 1771  class Helmholtz(LinearPDE): Line 2051  class Helmholtz(LinearPDE):
2051       @type name: C{string}       @type name: C{string}
2052       @return: the value of the coefficient  name       @return: the value of the coefficient  name
2053       @rtype: L{Data<escript.Data>}       @rtype: L{Data<escript.Data>}
2054       @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.  
2055       """       """
2056       if name == "A" :       if name == "A" :
2057           return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))*self.getCoefficient("k")           return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))*self.getCoefficient("k")
      elif name == "B" :  
          return escript.Data()  
      elif name == "C" :  
          return escript.Data()  
2058       elif name == "D" :       elif name == "D" :
2059           return self.getCoefficient("omega")           return self.getCoefficient("omega")
      elif name == "X" :  
          return escript.Data()  
2060       elif name == "Y" :       elif name == "Y" :
2061           return self.getCoefficient("f")           return self.getCoefficient("f")
2062       elif name == "d" :       elif name == "d" :
2063           return self.getCoefficient("alpha")           return self.getCoefficient("alpha")
2064       elif name == "y" :       elif name == "y" :
2065           return self.getCoefficient("g")           return self.getCoefficient("g")
2066       elif name == "d_contact" :       elif name == "Y_reduced" :
2067           return escript.Data()           return self.getCoefficient("f_reduced")
2068       elif name == "y_contact" :       elif name == "y_reduced" :
2069           return escript.Data()          return self.getCoefficient("g_reduced")
      elif name == "r" :  
          return self.getCoefficient("r")  
      elif name == "q" :  
          return self.getCoefficient("q")  
2070       else:       else:
2071          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name          return super(Helmholtz, self).getCoefficient(name)
2072    
2073  class LameEquation(LinearPDE):  class LameEquation(LinearPDE):
2074     """     """
2075     Class to define a Lame equation problem:     Class to define a Lame equation problem:
2076    
2077     M{-grad(S{mu}*(grad(u[i])[j]+grad(u[j])[i]))[j] - grad(S{lambda}*grad(u[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] }
2078    
2079     with natural boundary conditons:     with natural boundary conditons:
2080    
2081     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] }
2082    
2083     and constraints:     and constraints:
2084    
# Line 1821  class LameEquation(LinearPDE): Line 2089  class LameEquation(LinearPDE):
2089     def __init__(self,domain,debug=False):     def __init__(self,domain,debug=False):
2090        super(LameEquation, self).__init__(domain,\        super(LameEquation, self).__init__(domain,\
2091                                           domain.getDim(),domain.getDim(),debug)                                           domain.getDim(),domain.getDim(),debug)
2092        self.COEFFICIENTS={ "lame_lambda"  : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),        self.introduceCoefficients(lame_lambda=PDECoef(PDECoef.INTERIOR,(),PDECoef.OPERATOR),
2093                            "lame_mu"      : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),                                   lame_mu=PDECoef(PDECoef.INTERIOR,(),PDECoef.OPERATOR),
2094                            "F"            : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),                                   F=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2095                            "sigma"        : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM),PDECoefficient.RIGHTHANDSIDE),                                   sigma=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
2096                            "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)}  
2097        self.setSymmetryOn()        self.setSymmetryOn()
2098    
2099     def setValues(self,**coefficients):     def setValues(self,**coefficients):
2100       """       """
2101       sets new values to coefficients       sets new values to coefficients
2102    
# Line 1855  class LameEquation(LinearPDE): Line 2121  class LameEquation(LinearPDE):
2121       """       """
2122       super(LameEquation, self).setValues(**coefficients)       super(LameEquation, self).setValues(**coefficients)
2123    
2124     def getCoefficientOfGeneralPDE(self,name):     def getCoefficient(self,name):
2125       """       """
2126       return the value of the coefficient name of the general PDE       return the value of the coefficient name of the general PDE
2127    
# Line 1863  class LameEquation(LinearPDE): Line 2129  class LameEquation(LinearPDE):
2129       @type name: C{string}       @type name: C{string}
2130       @return: the value of the coefficient  name       @return: the value of the coefficient  name
2131       @rtype: L{Data<escript.Data>}       @rtype: L{Data<escript.Data>}
2132       @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.  
2133       """       """
2134       if name == "A" :       if name == "A" :
2135           out =self.createCoefficientOfGeneralPDE("A")           out =self.createCoefficient("A")
2136           for i in range(self.getDim()):           for i in range(self.getDim()):
2137             for j in range(self.getDim()):             for j in range(self.getDim()):
2138               out[i,i,j,j] += self.getCoefficient("lame_lambda")               out[i,i,j,j] += self.getCoefficient("lame_lambda")
2139               out[i,j,j,i] += self.getCoefficient("lame_mu")               out[i,j,j,i] += self.getCoefficient("lame_mu")
2140               out[i,j,i,j] += self.getCoefficient("lame_mu")               out[i,j,i,j] += self.getCoefficient("lame_mu")
2141           return out           return out
      elif name == "B" :  
          return escript.Data()  
      elif name == "C" :  
          return escript.Data()  
      elif name == "D" :  
          return escript.Data()  
2142       elif name == "X" :       elif name == "X" :
2143           return self.getCoefficient("sigma")           return self.getCoefficient("sigma")
2144       elif name == "Y" :       elif name == "Y" :
2145           return self.getCoefficient("F")           return self.getCoefficient("F")
      elif name == "d" :  
          return escript.Data()  
2146       elif name == "y" :       elif name == "y" :
2147           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")  
2148       else:       else:
2149          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name          return super(LameEquation, self).getCoefficient(name)
2150    
2151    def LinearSinglePDE(domain,debug=False):
2152       """
2153       defines a single linear PDEs
2154    
2155       @param domain: domain of the PDE
2156       @type domain: L{Domain<escript.Domain>}
2157       @param debug: if True debug informations are printed.
2158       @rtype: L{LinearPDE}
2159       """
2160       return LinearPDE(domain,numEquations=1,numSolutions=1,debug=debug)
2161    
2162    def LinearPDESystem(domain,debug=False):
2163       """
2164       defines a system of linear PDEs
2165    
2166  class AdvectivePDE(LinearPDE):     @param domain: domain of the PDE
2167       @type domain: L{Domain<escript.Domain>}
2168       @param debug: if True debug informations are printed.
2169       @rtype: L{LinearPDE}
2170     """     """
2171     In cases of PDEs dominated by the advection terms M{B} and M{C} against the adevctive terms M{A}     return LinearPDE(domain,numEquations=domain.getDim(),numSolutions=domain.getDim(),debug=debug)
2172     up-winding has been used.  The L{AdvectivePDE} class applies SUPG upwinding to the advective terms.  
2173    class TransportPDE(LinearProblem):
2174       """
2175       This class is used to define a transport problem given by a general linear, time dependend, second order PDE
2176       for an unknown, non-negative, time-dependend function M{u} on a given domain defined through a L{Domain<escript.Domain>} object.
2177    
2178       For a single equation with a solution with a single component the transport problem is defined in the following form:
2179    
2180       M{(M+M_reduced)*u_t=-(grad(A[j,l]+A_reduced[j,l])*grad(u)[l]+(B[j]+B_reduced[j])u)[j]+(C[l]+C_reduced[l])*grad(u)[l]+(D+D_reduced)-grad(X+X_reduced)[j,j]+(Y+Y_reduced)}
2181    
2182    
2183       where M{u_t} denotes the time derivative of M{u} and M{grad(F)} denotes the spatial derivative of M{F}. Einstein's summation convention,
2184       ie. summation over indexes appearing twice in a term of a sum is performed, is used.
2185       The coefficients M{M}, M{A}, M{B}, M{C}, M{D}, M{X} and M{Y} have to be specified through L{Data<escript.Data>} objects in the
2186       L{Function<escript.Function>} and the coefficients M{M_reduced}, M{A_reduced}, M{B_reduced}, M{C_reduced}, M{D_reduced}, M{X_reduced} and M{Y_reduced} have to be specified through L{Data<escript.Data>} objects in the
2187       L{ReducedFunction<escript.ReducedFunction>}. It is also allowed to use objects that can be converted into
2188       such L{Data<escript.Data>} objects. M{M} and M{M_reduced} are scalar, M{A} and M{A_reduced} are rank two,
2189       M{B}, M{C}, M{X}, M{B_reduced}, M{C_reduced} and M{X_reduced} are rank one and M{D}, M{D_reduced}, M{Y} and M{Y_reduced} are scalar.
2190    
2191       The following natural boundary conditions are considered:
2192    
2193       M{n[j]*((A[i,j]+A_reduced[i,j])*grad(u)[l]+(B+B_reduced)[j]*u+X[j]+X_reduced[j])+(d+d_reduced)*u+y+y_reduced=(m+m_reduced)*u_t}
2194    
2195       where M{n} is the outer normal field. Notice that the coefficients M{A}, M{A_reduced}, M{B}, M{B_reduced}, M{X} and M{X_reduced} are defined in the transport problem. The coefficients M{m}, M{d} and M{y} and are each a scalar in the L{FunctionOnBoundary<escript.FunctionOnBoundary>} and the coefficients M{m_reduced}, M{d_reduced} and M{y_reduced} and are each a scalar in the L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
2196    
2197       Constraints for the solution prescribing the value of the solution at certain locations in the domain. They have the form
2198    
2199       M{u_t=r}  where M{q>0}
2200    
2201     In the following we set     M{r} and M{q} are each scalar where M{q} is the characteristic function defining where the constraint is applied.
2202       The constraints override any other condition set by the transport problem or the boundary condition.
2203    
2204       The transport problem is symmetrical if
2205    
2206       M{A[i,j]=A[j,i]}  and M{B[j]=C[j]} and M{A_reduced[i,j]=A_reduced[j,i]}  and M{B_reduced[j]=C_reduced[j]}
2207    
2208       For a system and a solution with several components the transport problem has the form
2209    
2210       M{(M[i,k]+M_reduced[i,k])*u[k]_t=-grad((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k])[j]+(C[i,k,l]+C_reduced[i,k,l])*grad(u[k])[l]+(D[i,k]+D_reduced[i,k]*u[k]-grad(X[i,j]+X_reduced[i,j])[j]+Y[i]+Y_reduced[i] }
2211    
2212       M{A} and M{A_reduced} are of rank four, M{B}, M{B_reduced}, M{C} and M{C_reduced} are each of rank three, M{M}, M{M_reduced}, M{D}, M{D_reduced}, M{X_reduced} and M{X} are each a rank two and M{Y} and M{Y_reduced} are of rank one. The natural boundary conditions take the form:
2213    
2214       M{n[j]*((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k]+X[i,j]+X_reduced[i,j])+(d[i,k]+d_reduced[i,k])*u[k]+y[i]+y_reduced[i]= (m[i,k]+m_reduced[i,k])*u[k]_t}
2215    
2216       The coefficient M{d} and M{m} are of rank two and M{y} are of rank one with all in the L{FunctionOnBoundary<escript.FunctionOnBoundary>}. Constraints take the form and the coefficients M{d_reduced} and M{m_reduced} are of rank two and M{y_reduced} is of rank one with all in the L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
2217    
2218     M{Z[j]=C[j]-B[j]}     Constraints take the form
2219    
2220       M{u[i]_t=r[i]}  where  M{q[i]>0}
2221    
2222       M{r} and M{q} are each rank one. Notice that at some locations not necessarily all components must have a constraint.
2223    
2224       The transport problem is symmetrical if
2225    
2226            - M{M[i,k]=M[i,k]}
2227            - M{M_reduced[i,k]=M_reduced[i,k]}
2228            - M{A[i,j,k,l]=A[k,l,i,j]}
2229            - M{A_reduced[i,j,k,l]=A_reduced[k,l,i,j]}
2230            - M{B[i,j,k]=C[k,i,j]}
2231            - M{B_reduced[i,j,k]=C_reduced[k,i,j]}
2232            - M{D[i,k]=D[i,k]}
2233            - M{D_reduced[i,k]=D_reduced[i,k]}
2234            - M{m[i,k]=m[k,i]}
2235            - M{m_reduced[i,k]=m_reduced[k,i]}
2236            - M{d[i,k]=d[k,i]}
2237            - M{d_reduced[i,k]=d_reduced[k,i]}
2238    
2239     or     L{TransportPDE} also supports solution discontinuities over a contact region in the domain. To specify the conditions across the
2240       discontinuity we are using the generalised flux M{J} which is in the case of a systems of PDEs and several components of the solution
2241       defined as
2242    
2243     M{Z[i,k,l]=C[i,k,l]-B[i,l,k]}     M{J[i,j]=(A[i,j,k,l]+A_reduced[[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k]+X[i,j]+X_reduced[i,j]}
2244    
2245     To measure the dominance of the advective terms over the diffusive term M{A} the     For the case of single solution component and single PDE M{J} is defined
    X{Pelclet number} M{P} is used. It is defined as  
2246    
2247     M{P=h|Z|/(2|A|)}     M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[j]+(B[i]+B_reduced[i])*u+X[i]+X_reduced[i]}
2248    
2249     where M{|.|} denotes the L{length<util.length>} of the arument and M{h} is the local cell size     In the context of discontinuities M{n} denotes the normal on the discontinuity pointing from side 0 towards side 1
2250     from L{getSize<escript.Domain.getSize>}. Where M{|A|==0} M{P} is M{S{infinity}}.     calculated from L{getNormal<escript.FunctionSpace.getNormal>} of L{FunctionOnContactZero<escript.FunctionOnContactZero>}. For a system of transport problems the contact condition takes the form
2251    
2252     From the X{Pelclet number} the stabilization parameters M{S{Xi}} and M{S{Xi}} are calculated:     M{n[j]*J0[i,j]=n[j]*J1[i,j]=(y_contact[i]+y_contact_reduced[i])- (d_contact[i,k]+d_contact_reduced[i,k])*jump(u)[k]}
2253    
2254     M{S{Xi}=S{xi}(P) h/|Z|}     where M{J0} and M{J1} are the fluxes on side 0 and side 1 of the discontinuity, respectively. M{jump(u)}, which is the difference
2255       of the solution at side 1 and at side 0, denotes the jump of M{u} across discontinuity along the normal calculated by
2256       L{jump<util.jump>}.
2257       The coefficient M{d_contact} is a rank two and M{y_contact} is a rank one both in the L{FunctionOnContactZero<escript.FunctionOnContactZero>} or L{FunctionOnContactOne<escript.FunctionOnContactOne>}.
2258       The coefficient M{d_contact_reduced} is a rank two and M{y_contact_reduced} is a rank one both in the L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.
2259       In case of a single PDE and a single component solution the contact condition takes the form
2260    
2261     where M{S{xi}} is a suitable function of the Peclet number.     M{n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)}
2262    
2263     In the case of a single PDE the coefficient are up-dated in the following way:     In this case the coefficient M{d_contact} and M{y_contact} are each scalar both in the L{FunctionOnContactZero<escript.FunctionOnContactZero>} or L{FunctionOnContactOne<escript.FunctionOnContactOne>} and the coefficient M{d_contact_reduced} and M{y_contact_reduced} are each scalar both in the L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}
          - M{A[i,j] S{<-} A[i,j] + S{Xi} * Z[j] * Z[l]}  
          - M{B[j] S{<-} B[j] + S{Xi} * C[j] * D}  
          - M{C[j] S{<-} C[j] + S{Xi} * B[j] * D}  
          - M{X[j] S{<-} X[j] + S{Xi} * Z[j] * Y}  
2264    
2265     Similar for the case of a systems of PDEs:     typical usage:
          - 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]}  
          - M{B[i,j,k] S{<-} B[i,j,k] +  S{delta}[p,m] * S{Xi} * D[p,k] * C[m,i,j]}  
          - 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]}  
2266    
2267     where M{S{delta}} is L{kronecker}.        p=TransportPDE(dom)
2268     Using upwinding in this form, introduces an additonal error which is proprtional to the cell size M{h}        p.setValue(M=1.,C=[-1.,0.])
2269     but with the intension to stabilize the solution.        p.setInitialSolution(u=exp(-length(dom.getX()-[0.1,0.1])**2)
2270          t=0
2271          dt=0.1
2272          while (t<1.):
2273              u=p.solve(dt)
2274    
2275     """     """
2276     def __init__(self,domain,numEquations=None,numSolutions=None,xi=None,debug=False):     def __init__(self,domain,numEquations=None,numSolutions=None, theta=0.5,debug=False):
2277         """
2278         initializes a linear problem
2279    
2280         @param domain: domain of the PDE
2281         @type domain: L{Domain<escript.Domain>}
2282         @param numEquations: number of equations. If numEquations==None the number of equations
2283                              is extracted from the coefficients.
2284         @param numSolutions: number of solution components. If  numSolutions==None the number of solution components
2285                              is extracted from the coefficients.
2286         @param debug: if True debug informations are printed.
2287         @param theta: time stepping control: =1 backward Euler, =0 forward Euler, =0.5 Crank-Nicholson scheme, U{see here<http://en.wikipedia.org/wiki/Crank-Nicolson_method>}.
2288    
2289         """
2290         if theta<0 or theta>1:
2291                  raise ValueError,"theta (=%s) needs to be between 0 and 1 (inclusive)."%theta
2292         super(TransportPDE, self).__init__(domain,numEquations,numSolutions,debug)
2293    
2294         self.__theta=theta
2295         #
2296         #   the coefficients of the transport problem
2297         #
2298         self.introduceCoefficients(
2299           M=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2300           A=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2301           B=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2302           C=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2303           D=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2304           X=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
2305           Y=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2306           m=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2307           d=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2308           y=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2309           d_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2310           y_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2311           M_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2312           A_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2313           B_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2314           C_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2315           D_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2316           X_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
2317           Y_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2318           m_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2319           d_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2320           y_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2321           d_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2322           y_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2323           r=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.RIGHTHANDSIDE),
2324           q=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.BOTH) )
2325    
2326       def __str__(self):
2327         """
2328         returns string representation of the transport problem
2329    
2330         @return: a simple representation of the transport problem
2331         @rtype: C{str}
2332         """
2333         return "<TransportPDE %d>"%id(self)
2334    
2335       def getTheta(self):
2336          """
2337          returns the time stepping control parameter
2338          @rtype: float
2339        """        """
2340        creates a linear, steady, second order PDE on a L{Domain<escript.Domain>}        return self.__theta
2341    
2342        @param domain: domain of the PDE  
2343        @type domain: L{Domain<escript.Domain>}     def checkSymmetry(self,verbose=True):
       @param numEquations: number of equations. If numEquations==None the number of equations  
                            is exracted from the PDE coefficients.  
       @param numSolutions: number of solution components. If  numSolutions==None the number of solution components  
                            is exracted from the PDE coefficients.  
       @param xi: defines a function which returns for any given Preclet number as L{Scalar<escript.Scalar>} object the  
                  M{S{xi}}-value used to define the stabilization parameters. If equal to None, L{ELMAN_RAMAGE} is used.  
       @type xi: callable object which returns a L{Scalar<escript.Scalar>} object.  
       @param debug: if True debug informations are printed.  
2344        """        """
2345        super(AdvectivePDE, self).__init__(domain,\        test the transport problem for symmetry.
                                          numEquations,numSolutions,debug)  
       if xi==None:  
          self.__xi=AdvectivePDE.ELMAN_RAMAGE  
       else:  
          self.__xi=xi  
       self.__Xi=escript.Data()  
2346    
2347     def setValue(self,**coefficients):        @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.
2348          @type verbose: C{bool}
2349          @return:  True if the PDE is symmetric.
2350          @rtype: C{bool}
2351          @note: This is a very expensive operation. It should be used for degugging only! The symmetry flag is not altered.
2352          """
2353          out=True
2354          out=out and self.checkSymmetricTensor("M", verbose)
2355          out=out and self.checkSymmetricTensor("M_reduced", verbose)
2356          out=out and self.checkSymmetricTensor("A", verbose)
2357          out=out and self.checkSymmetricTensor("A_reduced", verbose)
2358          out=out and self.checkReciprocalSymmetry("B","C", verbose)
2359          out=out and self.checkReciprocalSymmetry("B_reduced","C_reduced", verbose)
2360          out=out and self.checkSymmetricTensor("D", verbose)
2361          out=out and self.checkSymmetricTensor("D_reduced", verbose)
2362          out=out and self.checkSymmetricTensor("m", verbose)
2363          out=out and self.checkSymmetricTensor("m_reduced", verbose)
2364          out=out and self.checkSymmetricTensor("d", verbose)
2365          out=out and self.checkSymmetricTensor("d_reduced", verbose)
2366          out=out and self.checkSymmetricTensor("d_contact", verbose)
2367          out=out and self.checkSymmetricTensor("d_contact_reduced", verbose)
2368          return out
2369    
2370       def setValue(self,**coefficients):
2371        """        """
2372        sets new values to coefficients        sets new values to coefficients
2373    
2374        @param coefficients: new values assigned to coefficients        @param coefficients: new values assigned to coefficients
2375          @keyword M: value for coefficient M.
2376          @type M: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
2377          @keyword M_reduced: value for coefficient M_reduced.
2378          @type M_reduced: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.ReducedFunction>}.
2379        @keyword A: value for coefficient A.        @keyword A: value for coefficient A.
2380        @type A: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.        @type A: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
2381          @keyword A_reduced: value for coefficient A_reduced.
2382          @type A_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
2383        @keyword B: value for coefficient B        @keyword B: value for coefficient B
2384        @type B: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.        @type B: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
2385          @keyword B_reduced: value for coefficient B_reduced
2386          @type B_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
2387        @keyword C: value for coefficient C        @keyword C: value for coefficient C
2388        @type C: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.        @type C: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
2389          @keyword C_reduced: value for coefficient C_reduced
2390          @type C_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
2391        @keyword D: value for coefficient D        @keyword D: value for coefficient D
2392        @type D: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.        @type D: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
2393          @keyword D_reduced: value for coefficient D_reduced
2394          @type D_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
2395        @keyword X: value for coefficient X        @keyword X: value for coefficient X
2396        @type X: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.        @type X: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
2397          @keyword X_reduced: value for coefficient X_reduced
2398          @type X_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
2399        @keyword Y: value for coefficient Y        @keyword Y: value for coefficient Y
2400        @type Y: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.        @type Y: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
2401          @keyword Y_reduced: value for coefficient Y_reduced
2402          @type Y_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.Function>}.
2403          @keyword m: value for coefficient m
2404          @type m: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
2405          @keyword m_reduced: value for coefficient m_reduced
2406          @type m_reduced: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
2407        @keyword d: value for coefficient d        @keyword d: value for coefficient d
2408        @type d: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.        @type d: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
2409          @keyword d_reduced: value for coefficient d_reduced
2410          @type d_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
2411        @keyword y: value for coefficient y        @keyword y: value for coefficient y
2412        @type y: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.        @type y: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
2413        @keyword d_contact: value for coefficient d_contact        @keyword d_contact: value for coefficient d_contact
2414        @type d_contact: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnContactOne<escript.FunctionOnContactOne>}.        @type d_contact: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnContactOne<escript.FunctionOnContactOne>} or  L{FunctionOnContactZero<escript.FunctionOnContactZero>}.
2415                         or  L{FunctionOnContactZero<escript.FunctionOnContactZero>}.        @keyword d_contact_reduced: value for coefficient d_contact_reduced
2416          @type d_contact_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>} or  L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>}.
2417        @keyword y_contact: value for coefficient y_contact        @keyword y_contact: value for coefficient y_contact
2418        @type y_contact: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnContactOne<escript.FunctionOnContactOne>}.        @type y_contact: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnContactOne<escript.FunctionOnContactOne>} or  L{FunctionOnContactZero<escript.FunctionOnContactZero>}.
2419                         or  L{FunctionOnContactZero<escript.FunctionOnContactZero>}.        @keyword y_contact_reduced: value for coefficient y_contact_reduced
2420          @type y_contact_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunctionOnContactOne<escript.FunctionOnContactOne>} or L{ReducedFunctionOnContactZero<escript.FunctionOnContactZero>}.
2421        @keyword r: values prescribed to the solution at the locations of constraints        @keyword r: values prescribed to the solution at the locations of constraints
2422        @type r: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}        @type r: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2423                 depending of reduced order is used for the solution.                 depending of reduced order is used for the solution.
# Line 2001  class AdvectivePDE(LinearPDE): Line 2425  class AdvectivePDE(LinearPDE):
2425        @type q: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}        @type q: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2426                 depending of reduced order is used for the representation of the equation.                 depending of reduced order is used for the representation of the equation.
2427        @raise IllegalCoefficient: if an unknown coefficient keyword is used.        @raise IllegalCoefficient: if an unknown coefficient keyword is used.
2428          """
2429          super(TransportPDE,self).setValue(**coefficients)
2430    
2431    
2432       def createOperator(self):
2433           """
2434           returns an instance of a new transport operator
2435           """
2436           self.trace("New Transport problem is allocated.")
2437           return self.getDomain().newTransportProblem( \
2438                                   self.getTheta(),
2439                                   self.getNumEquations(), \
2440                                   self.getFunctionSpaceForSolution(), \
2441                                   self.getSystemType())
2442    
2443       def setInitialSolution(self,u):
2444           """
2445           sets the initial solution
2446        
2447           @param u: new initial solution
2448           @type u: any object that can ne interpolated to a L{Data<escript.Data>} object on on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}.
2449           @note: C{u} must be non negative.
2450           """
2451           u2=util.interpolate(u,self.getFunctionSpaceForSolution())
2452           if self.getNumSolutions() == 1:
2453              if u2.getShape()!=():
2454                  raise ValueError,"Illegal shape %s of initial solution."%(u2.getShape(),)
2455           else:
2456              if u2.getShape()!=(self.getNumSolutions(),):
2457                  raise ValueError,"Illegal shape %s of initial solution."%(u2.getShape(),)
2458           self.getOperator().setInitialValue(u2)
2459       def getRequiredSystemType(self):
2460        """        """
2461        if "A" in coefficients.keys()   or "B" in coefficients.keys() or "C" in coefficients.keys(): self.__Xi=escript.Data()        returns the system type which needs to be used by the current set up.
       super(AdvectivePDE, self).setValue(**coefficients)  
2462    
2463     def ELMAN_RAMAGE(self,P):        @return: a code to indicate the type of transport problem scheme used.
2464       """        @rtype: C{float}
2465       Predefined function to set a values for M{S{xi}} from a Preclet number M{P}.        """
2466       This function uses the method suggested by H.C. Elman and A. Ramage, I{SIAM J. Numer. Anal.}, B{40} (2002)        return self.getDomain().getTransportTypeId(self.getSolverMethod()[0],self.getSolverMethod()[1],self.getSolverPackage(),self.isSymmetric())
           - M{S{xi}(P)=0} for M{P<1}  
           - M{S{xi}(P)=(1-1/P)/2} otherwise  
   
      @param P: Preclet number  
      @type P: L{Scalar<escript.Scalar>}  
      @return: up-wind weightimg factor  
      @rtype: L{Scalar<escript.Scalar>}  
      """  
      return util.wherePositive(P-1.)*0.5*(1.-1./(P+1.e-15))  
   
    def SIMPLIFIED_BROOK_HUGHES(self,P):  
      """  
      Predefined function to set a values for M{S{xi}} from a Preclet number M{P}.  
      The original methods is  
   
      M{S{xi}(P)=coth(P)-1/P}  
   
      As the evaluation of M{coth} is expensive we are using the approximation:  
   
          - M{S{xi}(P)=P/3} where M{P<3}  
          - M{S{xi}(P)=1/2} otherwise  
   
      @param P: Preclet number  
      @type P: L{Scalar<escript.Scalar>}  
      @return: up-wind weightimg factor  
      @rtype: L{Scalar<escript.Scalar>}  
      """  
      c=util.whereNegative(P-3.)  
      return P/6.*c+1./2.*(1.-c)  
   
    def HALF(self,P):  
      """  
      Predefined function to set value M{1/2} for M{S{xi}}  
   
      @param P: Preclet number  
      @type P: L{Scalar<escript.Scalar>}  
      @return: up-wind weightimg factor  
      @rtype: L{Scalar<escript.Scalar>}  
      """  
      return escript.Scalar(0.5,P.getFunctionSpace())  
   
    def __getXi(self):  
       if self.__Xi.isEmpty():  
          B=self.getCoefficient("B")  
          C=self.getCoefficient("C")  
          A=self.getCoefficient("A")  
          h=self.getDomain().getSize()  
          self.__Xi=escript.Scalar(0.,self.getFunctionSpaceForCoefficient("A"))  
          if not C.isEmpty() or not B.isEmpty():  
             if not C.isEmpty() and not B.isEmpty():  
                 if self.getNumEquations()>1:  
                    if self.getNumSolutions()>1:  
                       flux2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))  
                       for i in range(self.getNumEquations()):  
                          for k in range(self.getNumSolutions()):  
                             for l in range(self.getDim()): flux2+=(C[i,k,l]-B[i,l,k])**2  
                       length_of_flux=util.sqrt(flux2)  
                       # flux=C-util.reorderComponents(B,[0,2,1])  
                    else:  
                       flux2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))  
                       for i in range(self.getNumEquations()):  
                          for l in range(self.getDim()): flux2+=(C[i,l]-B[i,l])**2  
                       length_of_flux=util.sqrt(flux2)  
                       # flux=C-B  
                 else:  
                    if self.getNumSolutions()>1:  
                       flux2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))  
                       for k in range(self.getNumSolutions()):  
                          for l in range(self.getDim()): flux2+=(C[k,l]-B[l,k])**2  
                       # flux=C-util.reorderComponents(B,[1,0])  
                       length_of_flux=util.sqrt(flux2)  
                    else:  
                       length_of_flux=util.length(C-B)  
             elif C.isEmpty():  
               length_of_flux=util.length(B)  
             else:  
               length_of_flux=util.length(C)  
             flux_max=util.Lsup(length_of_flux)  
             if flux_max>0.:  
               if A.isEmpty():  
                   inv_A=1./self.SMALL_TOLERANCE  
                   peclet_number=escript.Scalar(inv_A,length_of_flux.getFunctionSpace())  
                   xi=self.__xi(self,peclet_number)  
               else:  
                   # length_of_A=util.inner(flux,util.tensormutiply(A,flux))  
                   length_of_A=util.length(A)  
                   A_max=util.Lsup(length_of_A)  
                   if A_max>0:  
                        inv_A=1./(length_of_A+A_max*self.SMALL_TOLERANCE)  
                   else:  
                        inv_A=1./self.SMALL_TOLERANCE  
                   peclet_number=length_of_flux*h/2*inv_A  
                   xi=self.__xi(self,peclet_number)  
               self.__Xi=h*xi/(length_of_flux+flux_max*self.SMALL_TOLERANCE)  
               self.trace("preclet number = %e"%util.Lsup(peclet_number))  
             else:  
               self.__Xi=escript.Scalar(0.,length_of_flux.getFunctionSpace())  
       return self.__Xi  
2467    
2468       def getUnlimitedTimeStepSize(self):
2469          """
2470          returns the value returned by the C{getSafeTimeStepSize} method to indicate no limit on the save time step size.
2471    
2472     def getCoefficientOfGeneralPDE(self,name):         @return: the value used to indicate that no limit is set to the time step size
2473       """         @rtype: C{float}
2474       return the value of the coefficient name of the general PDE         @note: Typically the biggest positive float is returned.
2475          """
2476          return self.getOperator().getUnlimitedTimeStepSize()
2477    
2478       @param name: name of the coefficient requested.     def getSafeTimeStepSize(self):
2479       @type name: C{string}         """
2480       @return: the value of the coefficient name         returns a safe time step size to do the next time step.
      @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}.  
      @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.  
      """  
      if not self.getNumEquations() == self.getNumSolutions():  
           raise ValueError,"AdvectivePDE expects the number of solution componets and the number of equations to be equal."  
2481    
2482       if name == "A" :         @return: safe time step size
2483           A=self.getCoefficient("A")         @rtype: C{float}
2484           B=self.getCoefficient("B")         @note: If not C{getSafeTimeStepSize()} < C{getUnlimitedTimeStepSize()} any time step size can be used.
2485           C=self.getCoefficient("C")         """
2486           if B.isEmpty() and C.isEmpty():         return self.getOperator().getSafeTimeStepSize()
2487              Aout=A     #====================================================================
2488           else:     def getSolution(self,dt,**options):
2489              if A.isEmpty():         """
2490                 Aout=self.createCoefficientOfGeneralPDE("A")         returns the solution of the problem
2491              else:         @return: the solution
2492                 Aout=A[:]         @rtype: L{Data<escript.Data>}
             Xi=self.__getXi()  
             if self.getNumEquations()>1:  
                 for i in range(self.getNumEquations()):  
                    for j in range(self.getDim()):  
                       for k in range(self.getNumSolutions()):  
                          for l in range(self.getDim()):  
                             if not C.isEmpty() and not B.isEmpty():  
                                # tmp=C-util.reorderComponents(B,[0,2,1])  
                                # Aout=Aout+Xi*util.generalTensorProduct(util.reorder(tmp,[1,2,0]),tmp,offset=1)  
                                for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*(C[p,i,j]-B[p,j,i])*(C[p,k,l]-B[p,l,k])  
                             elif C.isEmpty():  
                                for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*B[p,j,i]*B[p,l,k]  
                                # Aout=Aout+Xi*util.generalTensorProduct(util.reorder(B,[2,1,0]),util.reorder(B,[0,2,1]),offset=1)  
                             else:  
                                for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*C[p,i,j]*C[p,k,l]  
                                # Aout=Aout+Xi*util.generalTensorProduct(util.reorder(C,[1,2,0]),C,offset=1)  
             else:  
                if not C.isEmpty() and not B.isEmpty():  
                    delta=(C-B)  
                    Aout+=util.outer(Xi*delta,delta)  
                elif not B.isEmpty():  
                    Aout+=util.outer(Xi*B,B)  
                elif not C.isEmpty():  
                    Aout+=util.outer(Xi*C,C)  
          return Aout  
      elif name == "B" :  
          # return self.getCoefficient("B")  
          B=self.getCoefficient("B")  
          C=self.getCoefficient("C")  
          D=self.getCoefficient("D")  
          if C.isEmpty() or D.isEmpty():  
             Bout=B  
          else:  
             Xi=self.__getXi()  
             if B.isEmpty():  
                 Bout=self.createCoefficientOfGeneralPDE("B")  
             else:  
                 Bout=B[:]  
             if self.getNumEquations()>1:  
                for k in range(self.getNumSolutions()):  
                   for p in range(self.getNumEquations()):  
                      tmp=Xi*D[p,k]  
                      for i in range(self.getNumEquations()):  
                         for j in range(self.getDim()):  
                            Bout[i,j,k]+=tmp*C[p,i,j]  
                            # Bout=Bout+Xi*util.generalTensorProduct(util.reorder(C,[1,2,0]),D,offset=1)  
             else:  
                Bout+=(Xi*D)*C  
          return Bout  
      elif name == "C" :  
          # return self.getCoefficient("C")  
          B=self.getCoefficient("B")  
          C=self.getCoefficient("C")  
          D=self.getCoefficient("D")  
          if B.isEmpty() or D.isEmpty():  
             Cout=C  
          else:  
             Xi=self.__getXi()  
             if C.isEmpty():  
                 Cout=self.createCoefficientOfGeneralPDE("C")  
             else:  
                 Cout=C[:]  
             if self.getNumEquations()>1:  
                for k in range(self.getNumSolutions()):  
                    for p in range(self.getNumEquations()):  
                       tmp=Xi*D[p,k]  
                       for i in range(self.getNumEquations()):  
                         for l in range(self.getDim()):  
                                  Cout[i,k,l]+=tmp*B[p,l,i]  
                                  # Cout=Cout+Xi*B[p,l,i]*D[p,k]  
             else:  
                Cout+=(Xi*D)*B  
          return Cout  
      elif name == "D" :  
          return self.getCoefficient("D")  
      elif name == "X" :  
          # return self.getCoefficient("X")  
          X=self.getCoefficient("X")  
          Y=self.getCoefficient("Y")  
          B=self.getCoefficient("B")  
          C=self.getCoefficient("C")  
          if Y.isEmpty() or (B.isEmpty() and C.isEmpty()):  
             Xout=X  
          else:  
             if X.isEmpty():  
                 Xout=self.createCoefficientOfGeneralPDE("X")  
             else:  
                 Xout=X[:]  
             Xi=self.__getXi()  
             if self.getNumEquations()>1:  
                  for p in range(self.getNumEquations()):  
                     tmp=Xi*Y[p]  
                     for i in range(self.getNumEquations()):  
                        for j in range(self.getDim()):  
                           if not C.isEmpty() and not B.isEmpty():  
                              Xout[i,j]+=tmp*(C[p,i,j]-B[p,j,i])  
                              # Xout=X_out+Xi*util.inner(Y,C-util.reorderComponents(B,[0,2,1]),offset=1)  
                           elif C.isEmpty():  
                              Xout[i,j]-=tmp*B[p,j,i]  
                              # Xout=X_out-Xi*util.inner(Y,util.reorderComponents(B,[0,2,1]),offset=1)  
                           else:  
                              Xout[i,j]+=tmp*C[p,i,j]  
                              # Xout=X_out+Xi*util.inner(Y,C,offset=1)  
             else:  
               if not C.isEmpty() and not B.isEmpty():  
                 Xout+=(Xi*Y)*(C-B)  
               elif C.isEmpty():  
                 Xout-=(Xi*Y)*B  
               else:  
                 Xout+=(Xi*Y)*C  
          return Xout  
      elif name == "Y" :  
          return self.getCoefficient("Y")  
      elif name == "d" :  
          return self.getCoefficient("d")  
      elif name == "y" :  
          return self.getCoefficient("y")  
      elif name == "d_contact" :  
          return self.getCoefficient("d_contact")  
      elif name == "y_contact" :  
          return self.getCoefficient("y_contact")  
      elif name == "r" :  
          return self.getCoefficient("r")  
      elif name == "q" :  
          return self.getCoefficient("q")  
      else:  
         raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name  
2493    
2494  # $Log$         @note: This method is overwritten when implementing a particular linear problem.
2495  # Revision 1.14  2005/09/22 01:54:57  jgs         """
2496  # Merge of development branch dev-02 back to main trunk on 2005-09-22         if dt<=0:
2497  #             raise ValueError,"step size needs to be positive."
2498  # Revision 1.13  2005/09/15 03:44:19  jgs         options[self.TOLERANCE_KEY]=self.getTolerance()
2499  # Merge of development branch dev-02 back to main trunk on 2005-09-15         self.setSolution(self.getOperator().solve(self.getRightHandSide(),dt,options))
2500  #         self.validSolution()
2501  # Revision 1.12  2005/09/01 03:31:28  jgs         return self.getCurrentSolution()
2502  # Merge of development branch dev-02 back to main trunk on 2005-09-01  
2503  #     def getSystem(self):
2504  # Revision 1.11  2005/08/23 01:24:28  jgs         """
2505  # Merge of development branch dev-02 back to main trunk on 2005-08-23         return the operator and right hand side of the PDE
2506  #  
2507  # Revision 1.10  2005/08/12 01:45:36  jgs         @return: the discrete version of the PDE
2508  # erge of development branch dev-02 back to main trunk on 2005-08-12         @rtype: C{tuple} of L{Operator,<escript.Operator>} and L{Data<escript.Data>}.
2509  #  
2510  # Revision 1.9.2.17  2005/09/21 07:03:33  matt         @note: This method is overwritten when implementing a particular linear problem.
2511  # PDECoefficient and LinearPDE are now new style classes (introduced in Python         """
2512  # 2.2). Classes Poisson, Helmholtz, LameEquation and AdvectivePDE have been         if not self.isOperatorValid() or not self.isRightHandSideValid():
2513  # modified to instead use portable/cooperative "super" calls to extend base            self.resetRightHandSide()
2514  # class methods.            righthandside=self.getCurrentRightHandSide()
2515  #            self.resetOperator()
2516  # Revision 1.9.2.16  2005/09/16 01:54:37  matt            operator=self.getCurrentOperator()
2517  # Removed redundant if-loop.            self.getDomain().addPDEToTransportProblem(
2518  #                              operator,
2519  # Revision 1.9.2.15  2005/09/14 08:09:18  matt                              righthandside,
2520  # Added "REDUCED" solution PDECoefficient descriptors for LinearPDEs.                              self.getCoefficient("M"),
2521  #                              self.getCoefficient("A"),
2522  # Revision 1.9.2.14  2005/09/07 06:26:16  gross                              self.getCoefficient("B"),
2523  # the solver from finley are put into the standalone package paso now                              self.getCoefficient("C"),
2524  #                              self.getCoefficient("D"),
2525  # Revision 1.9.2.13  2005/08/31 08:45:03  gross                              self.getCoefficient("X"),
2526  # in the case of lumping no new system is allocated if the constraint is changed.                              self.getCoefficient("Y"),
2527  #                              self.getCoefficient("d"),
2528  # Revision 1.9.2.12  2005/08/31 07:10:23  gross                              self.getCoefficient("y"),
2529  # test for Lumping added                              self.getCoefficient("d_contact"),
2530  #                              self.getCoefficient("y_contact"))
2531  # Revision 1.9.2.11  2005/08/30 01:53:45  gross            self.getDomain().addPDEToTransportProblem(
2532  # bug in format fixed.                              operator,
2533  #                              righthandside,
2534  # Revision 1.9.2.10  2005/08/26 07:14:17  gross                              self.getCoefficient("M_reduced"),
2535  # a few more bugs in linearPDE fixed. remaining problem are finley problems                              self.getCoefficient("A_reduced"),
2536  #                              self.getCoefficient("B_reduced"),
2537  # Revision 1.9.2.9  2005/08/26 06:30:45  gross                              self.getCoefficient("C_reduced"),
2538  # fix for reported bug  0000004. test_linearPDE passes a few more tests                              self.getCoefficient("D_reduced"),
2539  #                              self.getCoefficient("X_reduced"),
2540  # Revision 1.9.2.8  2005/08/26 04:30:13  gross                              self.getCoefficient("Y_reduced"),
2541  # gneric unit testing for linearPDE                              self.getCoefficient("d_reduced"),
2542  #                              self.getCoefficient("y_reduced"),
2543  # Revision 1.9.2.7  2005/08/25 07:06:50  gross                              self.getCoefficient("d_contact_reduced"),
2544  # linearPDE documentation is parsed now by epydoc. there is still a problem with links into escriptcpp.so                              self.getCoefficient("y_contact_reduced"))
2545  #            operator.insertConstraint(righthandside,self.getCoefficient("q"),self.getCoefficient("r"))
2546  # Revision 1.9.2.6  2005/08/24 05:01:24  gross            self.trace("New system has been built.")
2547  # problem with resetting the matrix in case of resetting its values to 0 fixed.            self.validOperator()
2548  #            self.validRightHandSide()
2549  # Revision 1.9.2.5  2005/08/24 02:03:28  gross         return (self.getCurrentOperator(), self.getCurrentRightHandSide())
 # epydoc mark up partially fixed  
 #  
 # Revision 1.9.2.4  2005/08/22 07:11:09  gross  
 # some problems with LinearPDEs fixed.  
 #  
 # Revision 1.9.2.3  2005/08/18 04:48:48  gross  
 # the methods SetLumping*() are removed. Lumping is set trough setSolverMethod(LinearPDE.LUMPING)  
 #  
 # Revision 1.9.2.2  2005/08/18 04:39:32  gross  
 # the constants have been removed from util.py as they not needed anymore. PDE related constants are accessed through LinearPDE attributes now  
 #  
 # Revision 1.9.2.1  2005/07/29 07:10:27  gross  
 # new functions in util and a new pde type in linearPDEs  
 #  
 # Revision 1.1.2.25  2005/07/28 04:21:09  gross  
 # Lame equation: (linear elastic, isotropic) added  
 #  
 # Revision 1.1.2.24  2005/07/22 06:37:11  gross  
 # some extensions to modellib and linearPDEs  
 #  
 # Revision 1.1.2.23  2005/05/13 00:55:20  cochrane  
 # Fixed up some docstrings.  Moved module-level functions to top of file so  
 # that epydoc and doxygen can pick them up properly.  
 #  
 # Revision 1.1.2.22  2005/05/12 11:41:30  gross  
 # some basic Models have been added  
 #  
 # Revision 1.1.2.21  2005/05/12 07:16:12  cochrane  
 # Moved ELMAN_RAMAGE, SIMPLIFIED_BROOK_HUGHES, and HALF functions to bottom of  
 # file so that the AdvectivePDE class is picked up by doxygen.  Some  
 # reformatting of docstrings.  Addition of code to make equations come out  
 # as proper LaTeX.  
 #  
 # Revision 1.1.2.20  2005/04/15 07:09:08  gross  
 # some problems with functionspace and linearPDEs fixed.  
 #  
 # Revision 1.1.2.19  2005/03/04 05:27:07  gross  
 # bug in SystemPattern fixed.  
 #  
 # Revision 1.1.2.18  2005/02/08 06:16:45  gross  
 # Bugs in AdvectivePDE fixed, AdvectiveTest is stable but more testing is needed  
 #  
 # Revision 1.1.2.17  2005/02/08 05:56:19  gross  
 # Reference Number handling added  
 #  
 # Revision 1.1.2.16  2005/02/07 04:41:28  gross  
 # some function exposed to python to make mesh merging running  
 #  
 # Revision 1.1.2.15  2005/02/03 00:14:44  gross  
 # timeseries add and ESySParameter.py renames esysXML.py for consistence  
 #  
 # Revision 1.1.2.14  2005/02/01 06:44:10  gross  
 # new implementation of AdvectivePDE which now also updates right hand side. systems of PDEs are still not working  
 #  
 # Revision 1.1.2.13  2005/01/25 00:47:07  gross  
 # updates in the documentation  
 #  
 # Revision 1.1.2.12  2005/01/12 01:28:04  matt  
 # Added createCoefficient method for linearPDEs.  
 #  
 # Revision 1.1.2.11  2005/01/11 01:55:34  gross  
 # a problem in linearPDE class fixed  
 #  
 # Revision 1.1.2.10  2005/01/07 01:13:29  gross  
 # some bugs in linearPDE fixed  
 #  
 # Revision 1.1.2.9  2005/01/06 06:24:58  gross  
 # some bugs in slicing fixed  
 #  
 # Revision 1.1.2.8  2005/01/05 04:21:40  gross  
 # FunctionSpace checking/matchig in slicing added  
 #  
 # Revision 1.1.2.7  2004/12/29 10:03:41  gross  
 # bug in setValue fixed  
 #  
 # Revision 1.1.2.6  2004/12/29 05:29:59  gross  
 # AdvectivePDE successfully tested for Peclet number 1000000. there is still a problem with setValue and Data()  
 #  
 # Revision 1.1.2.5  2004/12/29 00:18:41  gross  
 # AdvectivePDE added  
 #  
 # Revision 1.1.2.4  2004/12/24 06:05:41  gross  
 # some changes in linearPDEs to add AdevectivePDE  
 #  
 # Revision 1.1.2.3  2004/12/16 00:12:34  gross  
 # __init__ of LinearPDE does not accept any coefficient anymore  
 #  
 # Revision 1.1.2.2  2004/12/14 03:55:01  jgs  
 # *** empty log message ***  
 #  
 # Revision 1.1.2.1  2004/12/12 22:53:47  gross  
 # linearPDE has been renamed LinearPDE  
 #  
 # Revision 1.1.1.1.2.7  2004/12/07 10:13:08  gross  
 # GMRES added  
 #  
 # Revision 1.1.1.1.2.6  2004/12/07 03:19:50  gross  
 # options for GMRES and PRES20 added  
 #  
 # Revision 1.1.1.1.2.5  2004/12/01 06:25:15  gross  
 # some small changes  
 #  
 # Revision 1.1.1.1.2.4  2004/11/24 01:50:21  gross  
 # Finley solves 4M unknowns now  
 #  
 # Revision 1.1.1.1.2.3  2004/11/15 06:05:26  gross  
 # poisson solver added  
 #  
 # Revision 1.1.1.1.2.2  2004/11/12 06:58:15  gross  
 # a lot of changes to get the linearPDE class running: most important change is that there is no matrix format exposed to the user anymore. the format is chosen by the Domain according to the solver and symmetry  
 #  
 # Revision 1.1.1.1.2.1  2004/10/28 22:59:22  gross  
 # finley's RecTest.py is running now: problem in SystemMatrixAdapater fixed  
 #  
 # Revision 1.1.1.1  2004/10/26 06:53:56  jgs  
 # initial import of project esys2  
 #  
 # Revision 1.3.2.3  2004/10/26 06:43:48  jgs  
 # committing Lutz's and Paul's changes to brach jgs  
 #  
 # Revision 1.3.4.1  2004/10/20 05:32:51  cochrane  
 # Added incomplete Doxygen comments to files, or merely put the docstrings that already exist into Doxygen form.  
 #  
 # Revision 1.3  2004/09/23 00:53:23  jgs  
 # minor fixes  
 #  
 # Revision 1.1  2004/08/28 12:58:06  gross  
 # SimpleSolve is not running yet: problem with == of functionsspace  
 #  

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

  ViewVC Help
Powered by ViewVC 1.1.26