/[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 425 by gross, Tue Jan 10 04:10:39 2006 UTC revision 1841 by gross, Fri Oct 3 03:57:52 2008 UTC
# Line 1  Line 1 
 # $Id$  
1    
2    ########################################################
3  #  #
4  #      COPYRIGHT ACcESS 2004 -  All Rights Reserved  # Copyright (c) 2003-2008 by University of Queensland
5  #  # Earth Systems Science Computational Center (ESSCC)
6  #   This software is the property of ACcESS.  No part of this code  # http://www.uq.edu.au/esscc
7  #   may be copied in any form or by any means without the expressed written  #
8  #   consent of ACcESS.  Copying, use or modification of this software  # Primary Business: Queensland, Australia
9  #   by any unauthorised person is illegal unless that  # Licensed under the Open Software License version 3.0
10  #   person has a software license agreement with ACcESS.  # 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 __licence__: licence agreement  @var __copyright__: copyrights
32    @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"
 __licence__="contact: esys@access.uq.edu.au"  
 __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 82  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 99  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 120  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 160  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 311  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 404  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 418  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
381       @cvar RILU: recursive ILU
382       @cvar GS: Gauss-Seidel solver
383    
384     """     """
385     DEFAULT= 0     DEFAULT= 0
# Line 443  class LinearPDE(object): Line 404  class LinearPDE(object):
404     UMFPACK= 16     UMFPACK= 16
405     ITERATIVE= 20     ITERATIVE= 20
406     PASO= 21     PASO= 21
407       AMG= 22
408     __TOL=1.e-13     RILU = 23
409     __PACKAGE_KEY="package"     TRILINOS = 24
410     __METHOD_KEY="method"     NONLINEAR_GMRES = 25
411     __SYMMETRY_KEY="symmetric"     TFQMR = 26
412     __TOLERANCE_KEY="tolerance"     MINRES = 27
413     __PRECONDITIONER_KEY="preconditioner"     GS=28
414    
415       SMALL_TOLERANCE=1.e-13
416       PACKAGE_KEY="package"
417       METHOD_KEY="method"
418       SYMMETRY_KEY="symmetric"
419       TOLERANCE_KEY="tolerance"
420       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 515  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 542  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 628  class LinearPDE(object): Line 594  class LinearPDE(object):
594       else:       else:
595           return escript.Solution(self.getDomain())           return escript.Solution(self.getDomain())
596    
   
    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.__TOL  
             if self.getNumSolutions()>1:  
                for i in range(self.getNumEquations()):  
                   for j in range(self.getDim()):  
                      for k in range(self.getNumSolutions()):  
                         for l in range(self.getDim()):  
                             if util.Lsup(A[i,j,k,l]-A[k,l,i,j])>tol:  
                                if verbose: print "non-symmetric PDE because A[%d,%d,%d,%d]!=A[%d,%d,%d,%d]"%(i,j,k,l,k,l,i,j)  
                                out=False  
             else:  
                for j in range(self.getDim()):  
                   for l in range(self.getDim()):  
                      if util.Lsup(A[j,l]-A[l,j])>tol:  
                         if verbose: print "non-symmetric PDE because A[%d,%d]!=A[%d,%d]"%(j,l,l,j)  
                         out=False  
          B=self.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.__TOL/2.  
             if self.getNumSolutions()>1:  
                for i in range(self.getNumEquations()):  
                    for j in range(self.getDim()):  
                       for k in range(self.getNumSolutions()):  
                          if util.Lsup(B[i,j,k]-C[k,i,j])>tol:  
                               if verbose: print "non-symmetric PDE because B[%d,%d,%d]!=C[%d,%d,%d]"%(i,j,k,k,i,j)  
                               out=False  
             else:  
                for j in range(self.getDim()):  
                   if util.Lsup(B[j]-C[j])>tol:  
                      if verbose: print "non-symmetric PDE because B[%d]!=C[%d]"%(j,j)  
                      out=False  
          if self.getNumSolutions()>1:  
            D=self.getCoefficientOfGeneralPDE("D")  
            if not D.isEmpty():  
              tol=util.Lsup(D)*self.__TOL  
              for i in range(self.getNumEquations()):  
                 for k in range(self.getNumSolutions()):  
                   if util.Lsup(D[i,k]-D[k,i])>tol:  
                       if verbose: print "non-symmetric PDE because D[%d,%d]!=D[%d,%d]"%(i,k,k,i)  
                       out=False  
            d=self.getCoefficientOfGeneralPDE("d")  
            if not d.isEmpty():  
              tol=util.Lsup(d)*self.__TOL  
              for i in range(self.getNumEquations()):  
                 for k in range(self.getNumSolutions()):  
                   if util.Lsup(d[i,k]-d[k,i])>tol:  
                       if verbose: print "non-symmetric PDE because d[%d,%d]!=d[%d,%d]"%(i,k,k,i)  
                       out=False  
            d_contact=self.getCoefficientOfGeneralPDE("d_contact")  
            if not d_contact.isEmpty():  
              tol=util.Lsup(d_contact)*self.__TOL  
              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")  
597     # =============================================================================     # =============================================================================
598     #   solver settings:     #   solver settings:
599     # =============================================================================     # =============================================================================
# Line 821  class LinearPDE(object): Line 602  class LinearPDE(object):
602         sets a new solver         sets a new solver
603    
604         @param solver: sets a new solver method.         @param solver: sets a new solver method.
605         @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}.         @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}
606         @param preconditioner: sets a new solver method.         @param preconditioner: sets a new solver method.
607         @type solver: one of L{DEFAULT}, L{JACOBI} L{ILU0}, L{ILUT},L{SSOR}         @type preconditioner: one of L{DEFAULT}, L{JACOBI} L{ILU0}, L{ILUT},L{SSOR}, L{RILU},  L{GS}
608    
609           @note: the solver method choosen may not be available by the selected package.
610         """         """
611         if solver==None: solve=self.DEFAULT         if solver==None: solver=self.__solver_method
612           if preconditioner==None: preconditioner=self.__preconditioner
613           if solver==None: solver=self.DEFAULT
614         if preconditioner==None: preconditioner=self.DEFAULT         if preconditioner==None: preconditioner=self.DEFAULT
615         if not (solver,preconditioner)==self.getSolverMethod():         if not (solver,preconditioner)==self.getSolverMethod():
616             self.__solver_method=solver             self.__solver_method=solver
617             self.__preconditioner=preconditioner             self.__preconditioner=preconditioner
618             self.__checkMatrixType()             self.updateSystemType()
619             self.trace("New solver is %s"%self.getSolverMethodName())             self.trace("New solver is %s"%self.getSolverMethodName())
620    
621     def getSolverMethodName(self):     def getSolverMethodName(self):
622         """         """
623         returns the name of the solver currently used         returns the name of the solver currently used
624    
# Line 849  class LinearPDE(object): Line 634  class LinearPDE(object):
634         elif m[0]==self.ITERATIVE: method= "ITERATIVE"         elif m[0]==self.ITERATIVE: method= "ITERATIVE"
635         elif m[0]==self.CHOLEVSKY: method= "CHOLEVSKY"         elif m[0]==self.CHOLEVSKY: method= "CHOLEVSKY"
636         elif m[0]==self.PCG: method= "PCG"         elif m[0]==self.PCG: method= "PCG"
637           elif m[0]==self.TFQMR: method= "TFQMR"
638           elif m[0]==self.MINRES: method= "MINRES"
639         elif m[0]==self.CR: method= "CR"         elif m[0]==self.CR: method= "CR"
640         elif m[0]==self.CGS: method= "CGS"         elif m[0]==self.CGS: method= "CGS"
641         elif m[0]==self.BICGSTAB: method= "BICGSTAB"         elif m[0]==self.BICGSTAB: method= "BICGSTAB"
# Line 856  class LinearPDE(object): Line 643  class LinearPDE(object):
643         elif m[0]==self.GMRES: method= "GMRES"         elif m[0]==self.GMRES: method= "GMRES"
644         elif m[0]==self.PRES20: method= "PRES20"         elif m[0]==self.PRES20: method= "PRES20"
645         elif m[0]==self.LUMPING: method= "LUMPING"         elif m[0]==self.LUMPING: method= "LUMPING"
646           elif m[0]==self.AMG: method= "AMG"
647         if m[1]==self.DEFAULT: method+="+DEFAULT"         if m[1]==self.DEFAULT: method+="+DEFAULT"
648         elif m[1]==self.JACOBI: method+= "+JACOBI"         elif m[1]==self.JACOBI: method+= "+JACOBI"
649         elif m[1]==self.ILU0: method+= "+ILU0"         elif m[1]==self.ILU0: method+= "+ILU0"
650         elif m[1]==self.ILUT: method+= "+ILUT"         elif m[1]==self.ILUT: method+= "+ILUT"
651         elif m[1]==self.SSOR: method+= "+SSOR"         elif m[1]==self.SSOR: method+= "+SSOR"
652           elif m[1]==self.AMG: method+= "+AMG"
653           elif m[1]==self.RILU: method+= "+RILU"
654           elif m[1]==self.GS: method+= "+GS"
655         if p==self.DEFAULT: package="DEFAULT"         if p==self.DEFAULT: package="DEFAULT"
656         elif p==self.PASO: package= "PASO"         elif p==self.PASO: package= "PASO"
657         elif p==self.MKL: package= "MKL"         elif p==self.MKL: package= "MKL"
658         elif p==self.SCSL: package= "SCSL"         elif p==self.SCSL: package= "SCSL"
659         elif p==self.UMFPACK: package= "UMFPACK"         elif p==self.UMFPACK: package= "UMFPACK"
660           elif p==self.TRILINOS: package= "TRILINOS"
661         else : method="unknown"         else : method="unknown"
662         return "%s solver of %s package"%(method,package)         return "%s solver of %s package"%(method,package)
663    
664       def getSolverMethod(self):
    def getSolverMethod(self):  
665         """         """
666         returns the solver method         returns the solver method and preconditioner used
667    
668         @return: the solver method currently be used.         @return: the solver method currently be used.
669         @rtype: C{int}         @rtype: C{tuple} of C{int}
670         """         """
671         return self.__solver_method,self.__preconditioner         return self.__solver_method,self.__preconditioner
672    
673     def setSolverPackage(self,package=None):     def setSolverPackage(self,package=None):
674         """         """
675         sets a new solver package         sets a new solver package
676    
677         @param solver: sets a new solver method.         @param package: sets a new solver method.
678         @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}
679         """         """
680         if package==None: package=self.DEFAULT         if package==None: package=self.DEFAULT
681         if not package==self.getSolverPackage():         if not package==self.getSolverPackage():
682             self.__solver_method=solver             self.__solver_package=package
683             self.__checkMatrixType()             self.updateSystemType()
684             self.trace("New solver is %s"%self.getSolverMethodName())             self.trace("New solver is %s"%self.getSolverMethodName())
685    
686     def getSolverPackage(self):     def getSolverPackage(self):
687         """         """
688         returns the package of the solver         returns the package of the solver
689    
# Line 903  class LinearPDE(object): Line 694  class LinearPDE(object):
694    
695     def isUsingLumping(self):     def isUsingLumping(self):
696        """        """
697        checks if matrix lumping is used a solver method        checks if matrix lumping is used as a solver method
698    
699        @return: True is lumping is currently used a solver method.        @return: True is lumping is currently used a solver method.
700        @rtype: C{bool}        @rtype: C{bool}
701        """        """
702        return self.getSolverMethod()[0]==self.LUMPING        return self.getSolverMethod()[0]==self.LUMPING
703    
704     def setTolerance(self,tol=1.e-8):     def setTolerance(self,tol=1.e-8):
705         """         """
706         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.  
707    
708         @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
709                     the system will be resolved.                     the system will be resolved.
710         @type tol: positive C{float}         @type tol: positive C{float}
711         @raise ValueException: if tolerance is not positive.         @raise ValueError: if tolerance is not positive.
712         """         """
713         if not tol>0:         if not tol>0:
714             raise ValueException,"Tolerance as to be positive"             raise ValueError,"Tolerance as to be positive"
715         if tol<self.getTolerance(): self.__invalidateSolution()         if tol<self.getTolerance(): self.invalidateSolution()
716         self.trace("New tolerance %e"%tol)         self.trace("New tolerance %e"%tol)
717         self.__tolerance=tol         self.__tolerance=tol
718         return         return
719    
720     def getTolerance(self):     def getTolerance(self):
721         """         """
722         returns the tolerance set for the solution         returns the tolerance set for the solution
723    
# Line 942  class LinearPDE(object): Line 729  class LinearPDE(object):
729     # =============================================================================     # =============================================================================
730     #    symmetry  flag:     #    symmetry  flag:
731     # =============================================================================     # =============================================================================
732     def isSymmetric(self):     def isSymmetric(self):
733        """        """
734        checks if symmetry is indicated.        checks if symmetry is indicated.
735    
# Line 951  class LinearPDE(object): Line 738  class LinearPDE(object):
738        """        """
739        return self.__sym        return self.__sym
740    
741     def setSymmetryOn(self):     def setSymmetryOn(self):
742        """        """
743        sets the symmetry flag.        sets the symmetry flag.
744        """        """
745        if not self.isSymmetric():        if not self.isSymmetric():
746           self.trace("PDE is set to be symmetric")           self.trace("PDE is set to be symmetric")
747           self.__sym=True           self.__sym=True
748           self.__checkMatrixType()           self.updateSystemType()
749    
750     def setSymmetryOff(self):     def setSymmetryOff(self):
751        """        """
752        removes the symmetry flag.        removes the symmetry flag.
753        """        """
754        if self.isSymmetric():        if self.isSymmetric():
755           self.trace("PDE is set to be unsymmetric")           self.trace("PDE is set to be unsymmetric")
756           self.__sym=False           self.__sym=False
757           self.__checkMatrixType()           self.updateSystemType()
758    
759     def setSymmetryTo(self,flag=False):     def setSymmetryTo(self,flag=False):
760        """        """
761        sets the symmetry flag to flag        sets the symmetry flag to flag
762    
# Line 984  class LinearPDE(object): Line 771  class LinearPDE(object):
771     # =============================================================================     # =============================================================================
772     # function space handling for the equation as well as the solution     # function space handling for the equation as well as the solution
773     # =============================================================================     # =============================================================================
774     def setReducedOrderOn(self):     def setReducedOrderOn(self):
775       """       """
776       switches on reduced order for solution and equation representation       switches on reduced order for solution and equation representation
777    
# Line 993  class LinearPDE(object): Line 780  class LinearPDE(object):
780       self.setReducedOrderForSolutionOn()       self.setReducedOrderForSolutionOn()
781       self.setReducedOrderForEquationOn()       self.setReducedOrderForEquationOn()
782    
783     def setReducedOrderOff(self):     def setReducedOrderOff(self):
784       """       """
785       switches off reduced order for solution and equation representation       switches off reduced order for solution and equation representation
786    
# Line 1002  class LinearPDE(object): Line 789  class LinearPDE(object):
789       self.setReducedOrderForSolutionOff()       self.setReducedOrderForSolutionOff()
790       self.setReducedOrderForEquationOff()       self.setReducedOrderForEquationOff()
791    
792     def setReducedOrderTo(self,flag=False):     def setReducedOrderTo(self,flag=False):
793       """       """
794       sets order reduction for both solution and equation representation according to flag.       sets order reduction for both solution and equation representation according to flag.
795       @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 1014  class LinearPDE(object): Line 801  class LinearPDE(object):
801       self.setReducedOrderForEquationTo(flag)       self.setReducedOrderForEquationTo(flag)
802    
803    
804     def setReducedOrderForSolutionOn(self):     def setReducedOrderForSolutionOn(self):
805       """       """
806       switches on reduced order for solution representation       switches on reduced order for solution representation
807    
# Line 1025  class LinearPDE(object): Line 812  class LinearPDE(object):
812                raise RuntimeError,"order cannot be altered after coefficients have been defined."                raise RuntimeError,"order cannot be altered after coefficients have been defined."
813           self.trace("Reduced order is used to solution representation.")           self.trace("Reduced order is used to solution representation.")
814           self.__reduce_solution_order=True           self.__reduce_solution_order=True
815           self.__resetSystem()           self.initializeSystem()
816    
817     def setReducedOrderForSolutionOff(self):     def setReducedOrderForSolutionOff(self):
818       """       """
819       switches off reduced order for solution representation       switches off reduced order for solution representation
820    
# Line 1038  class LinearPDE(object): Line 825  class LinearPDE(object):
825                raise RuntimeError,"order cannot be altered after coefficients have been defined."                raise RuntimeError,"order cannot be altered after coefficients have been defined."
826           self.trace("Full order is used to interpolate solution.")           self.trace("Full order is used to interpolate solution.")
827           self.__reduce_solution_order=False           self.__reduce_solution_order=False
828           self.__resetSystem()           self.initializeSystem()
829    
830     def setReducedOrderForSolutionTo(self,flag=False):     def setReducedOrderForSolutionTo(self,flag=False):
831       """       """
832       sets order for test functions according to flag       sets order for test functions according to flag
833    
# Line 1054  class LinearPDE(object): Line 841  class LinearPDE(object):
841       else:       else:
842          self.setReducedOrderForSolutionOff()          self.setReducedOrderForSolutionOff()
843    
844     def setReducedOrderForEquationOn(self):     def setReducedOrderForEquationOn(self):
845       """       """
846       switches on reduced order for equation representation       switches on reduced order for equation representation
847    
# Line 1065  class LinearPDE(object): Line 852  class LinearPDE(object):
852                raise RuntimeError,"order cannot be altered after coefficients have been defined."                raise RuntimeError,"order cannot be altered after coefficients have been defined."
853           self.trace("Reduced order is used for test functions.")           self.trace("Reduced order is used for test functions.")
854           self.__reduce_equation_order=True           self.__reduce_equation_order=True
855           self.__resetSystem()           self.initializeSystem()
856    
857     def setReducedOrderForEquationOff(self):     def setReducedOrderForEquationOff(self):
858       """       """
859       switches off reduced order for equation representation       switches off reduced order for equation representation
860    
# Line 1078  class LinearPDE(object): Line 865  class LinearPDE(object):
865                raise RuntimeError,"order cannot be altered after coefficients have been defined."                raise RuntimeError,"order cannot be altered after coefficients have been defined."
866           self.trace("Full order is used for test functions.")           self.trace("Full order is used for test functions.")
867           self.__reduce_equation_order=False           self.__reduce_equation_order=False
868           self.__resetSystem()           self.initializeSystem()
869    
870     def setReducedOrderForEquationTo(self,flag=False):     def setReducedOrderForEquationTo(self,flag=False):
871       """       """
872       sets order for test functions according to flag       sets order for test functions according to flag
873    
# Line 1094  class LinearPDE(object): Line 881  class LinearPDE(object):
881       else:       else:
882          self.setReducedOrderForEquationOff()          self.setReducedOrderForEquationOff()
883    
884     # =============================================================================     def updateSystemType(self):
    # private method:  
    # =============================================================================  
    def __checkMatrixType(self):  
885       """       """
886       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.
887       """       """
888       new_matrix_type=self.getDomain().getSystemMatrixTypeId(self.getSolverMethod()[0],self.getSolverPackage(),self.isSymmetric())       new_system_type=self.getRequiredSystemType()
889       if not new_matrix_type==self.__matrix_type:       if not new_system_type==self.__system_type:
890           self.trace("Matrix type is now %d."%new_matrix_type)           self.trace("Matrix type is now %d."%new_system_type)
891           self.__matrix_type=new_matrix_type           self.__system_type=new_system_type
892           self.__resetSystem()           self.initializeSystem()
893     #     def getSystemType(self):
894     #   rebuild switches :        """
895     #        return the current system type
896     def __invalidateSolution(self):        """
897         """        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)  
898    
899     def __makeFreshSolution(self):     def checkSymmetricTensor(self,name,verbose=True):
900         """        """
901         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  
902    
903     def __makeFreshRightHandSide(self):        @param name: name of the coefficient
904         """        @type name: C{str}
905         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.
906         """        @type verbose: C{bool}
907         if self.__righthandside.isEmpty():        @return:  True if coefficient C{name} is symmetric.
908             self.__righthandside=self.__getNewRightHandSide()        @rtype: C{bool}
909         else:        """
910             self.__righthandside*=0        A=self.getCoefficient(name)
911             self.trace("Right hand side is reset to zero.")        verbose=verbose or self.__debug
912         return self.__righthandside        out=True
913          if not A.isEmpty():
914             tol=util.Lsup(A)*self.SMALL_TOLERANCE
915             s=A.getShape()
916             if A.getRank() == 4:
917                if s[0]==s[2] and s[1] == s[3]:
918                   for i in range(s[0]):
919                      for j in range(s[1]):
920                         for k in range(s[2]):
921                            for l in range(s[3]):
922                                if util.Lsup(A[i,j,k,l]-A[k,l,i,j])>tol:
923                                   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)
924                                   out=False
925                else:
926                     if verbose: print "non-symmetric problem because of inappropriate shape %s of coefficient %s."%(s,name)
927                     out=False
928             elif A.getRank() == 2:
929                if s[0]==s[1]:
930                   for j in range(s[0]):
931                      for l in range(s[1]):
932                         if util.Lsup(A[j,l]-A[l,j])>tol:
933                            if verbose: print "non-symmetric problem because %s[%d,%d]!=%s[%d,%d]"%(name,j,l,name,l,j)
934                            out=False
935                else:
936                     if verbose: print "non-symmetric problem because of inappropriate shape %s of coefficient %s."%(s,name)
937                     out=False
938             elif A.getRank() == 0:
939                pass
940             else:
941                 raise ValueError,"Cannot check rank %s of %s."%(A.getRank(),name)
942          return out
943    
944     def __makeFreshOperator(self):     def checkReciprocalSymmetry(self,name0,name1,verbose=True):
945         """        """
946         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  
947    
948     def __applyConstraint(self):        @param name0: name of the first coefficient
949         """        @type name: C{str}
950         applies the constraints defined by q and r to the system        @param name1: name of the second coefficient
951         """        @type name: C{str}
952         if not self.isUsingLumping():        @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.
953            q=self.getCoefficientOfGeneralPDE("q")        @type verbose: C{bool}
954            r=self.getCoefficientOfGeneralPDE("r")        @return:  True if coefficients C{name0} and C{name1}  are resciprocally symmetric.
955            if not q.isEmpty() and not self.__operator.isEmpty():        @rtype: C{bool}
956               # q is the row and column mask to indicate where constraints are set:        """
957               row_q=escript.Data(q,self.getFunctionSpaceForEquation())        B=self.getCoefficient(name0)
958               col_q=escript.Data(q,self.getFunctionSpaceForSolution())        C=self.getCoefficient(name1)
959               u=self.__getNewSolution()        verbose=verbose or self.__debug
960               if r.isEmpty():        out=True
961                  r_s=self.__getNewSolution()        if B.isEmpty() and not C.isEmpty():
962             if verbose: print "non-symmetric problem because %s is not present but %s is"%(name0,name1)
963             out=False
964          elif not B.isEmpty() and C.isEmpty():
965             if verbose: print "non-symmetric problem because %s is not present but %s is"%(name0,name1)
966             out=False
967          elif not B.isEmpty() and not C.isEmpty():
968             sB=B.getShape()
969             sC=C.getShape()
970             tol=(util.Lsup(B)+util.Lsup(C))*self.SMALL_TOLERANCE/2.
971             if len(sB) != len(sC):
972                 if verbose: print "non-symmetric problem because ranks of %s (=%s) and %s (=%s) are different."%(name0,len(sB),name1,len(sC))
973                 out=False
974             else:
975                 if len(sB)==0:
976                   if util.Lsup(B-C)>tol:
977                      if verbose: print "non-symmetric problem because %s!=%s"%(name0,name1)
978                      out=False
979                 elif len(sB)==1:
980                   if sB[0]==sC[0]:
981                      for j in range(sB[0]):
982                         if util.Lsup(B[j]-C[j])>tol:
983                            if verbose: print "non-symmetric PDE because %s[%d]!=%s[%d]"%(name0,j,name1,j)
984                            out=False
985                   else:
986                     if verbose: print "non-symmetric problem because of inappropriate shapes %s and %s of coefficients %s and %s, respectively."%(sB,sC,name0,name1)
987                 elif len(sB)==3:
988                   if sB[0]==sC[1] and sB[1]==sC[2] and sB[2]==sC[0]:
989                       for i in range(sB[0]):
990                          for j in range(sB[1]):
991                             for k in range(sB[2]):
992                                if util.Lsup(B[i,j,k]-C[k,i,j])>tol:
993                                     if verbose: print "non-symmetric problem because %s[%d,%d,%d]!=%s[%d,%d,%d]"%(name2,i,j,k,name1,k,i,j)
994                                     out=False
995                   else:
996                     if verbose: print "non-symmetric problem because of inappropriate shapes %s and %s of coefficients %s and %s, respectively."%(sB,sC,name0,name1)
997               else:               else:
998                  r_s=escript.Data(r,self.getFunctionSpaceForSolution())                   raise ValueError,"Cannot check rank %s of %s and %s."%(len(sB),name0,name1)
999               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  
1000    
1001     # =============================================================================     def getCoefficient(self,name):
    # functions giving access to coefficients of a particular PDE implementation:  
    # =============================================================================  
    def getCoefficient(self,name):  
1002       """       """
1003       returns the value of the coefficient name       returns the value of the coefficient name
1004    
# Line 1336  class LinearPDE(object): Line 1009  class LinearPDE(object):
1009       @raise IllegalCoefficient: if name is not a coefficient of the PDE.       @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1010       """       """
1011       if self.hasCoefficient(name):       if self.hasCoefficient(name):
1012           return self.COEFFICIENTS[name].getValue()           return self.__COEFFICIENTS[name].getValue()
1013       else:       else:
1014          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1015    
1016     def hasCoefficient(self,name):     def hasCoefficient(self,name):
1017       """       """
1018       return True if name is the name of a coefficient       return True if name is the name of a coefficient
1019    
# Line 1349  class LinearPDE(object): Line 1022  class LinearPDE(object):
1022       @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.
1023       @rtype: C{bool}       @rtype: C{bool}
1024       """       """
1025       return self.COEFFICIENTS.has_key(name)       return self.__COEFFICIENTS.has_key(name)
1026    
1027     def createCoefficient(self, name):     def createCoefficient(self, name):
1028       """       """
1029       create a L{Data<escript.Data>} object corresponding to coefficient name       create a L{Data<escript.Data>} object corresponding to coefficient name
1030    
# Line 1364  class LinearPDE(object): Line 1037  class LinearPDE(object):
1037       else:       else:
1038          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1039    
1040     def getFunctionSpaceForCoefficient(self,name):     def getFunctionSpaceForCoefficient(self,name):
1041       """       """
1042       return the L{FunctionSpace<escript.FunctionSpace>} to be used for coefficient name       return the L{FunctionSpace<escript.FunctionSpace>} to be used for coefficient name
1043    
# Line 1375  class LinearPDE(object): Line 1048  class LinearPDE(object):
1048       @raise IllegalCoefficient: if name is not a coefficient of the PDE.       @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1049       """       """
1050       if self.hasCoefficient(name):       if self.hasCoefficient(name):
1051          return self.COEFFICIENTS[name].getFunctionSpace(self.getDomain())          return self.__COEFFICIENTS[name].getFunctionSpace(self.getDomain())
1052       else:       else:
1053          raise ValueError,"unknown coefficient %s requested"%name          raise ValueError,"unknown coefficient %s requested"%name
1054     def getShapeOfCoefficient(self,name):  
1055       def getShapeOfCoefficient(self,name):
1056       """       """
1057       return the shape of the coefficient name       return the shape of the coefficient name
1058    
# Line 1389  class LinearPDE(object): Line 1063  class LinearPDE(object):
1063       @raise IllegalCoefficient: if name is not a coefficient of the PDE.       @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1064       """       """
1065       if self.hasCoefficient(name):       if self.hasCoefficient(name):
1066          return self.COEFFICIENTS[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())          return self.__COEFFICIENTS[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())
1067       else:       else:
1068          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1069    
1070     def resetCoefficients(self):     def resetAllCoefficients(self):
1071       """       """
1072       resets all coefficients to there default values.       resets all coefficients to there default values.
1073       """       """
1074       for i in self.COEFFICIENTS.iterkeys():       for i in self.__COEFFICIENTS.iterkeys():
1075           self.COEFFICIENTS[i].resetValue()           self.__COEFFICIENTS[i].resetValue()
1076    
1077     def alteredCoefficient(self,name):     def alteredCoefficient(self,name):
1078       """       """
1079       announce that coefficient name has been changed       announce that coefficient name has been changed
1080    
# Line 1412  class LinearPDE(object): Line 1086  class LinearPDE(object):
1086       if self.hasCoefficient(name):       if self.hasCoefficient(name):
1087          self.trace("Coefficient %s has been altered."%name)          self.trace("Coefficient %s has been altered."%name)
1088          if not ((name=="q" or name=="r") and self.isUsingLumping()):          if not ((name=="q" or name=="r") and self.isUsingLumping()):
1089             if self.COEFFICIENTS[name].isAlteringOperator(): self.__invalidateOperator()             if self.__COEFFICIENTS[name].isAlteringOperator(): self.invalidateOperator()
1090             if self.COEFFICIENTS[name].isAlteringRightHandSide(): self.__invalidateRightHandSide()             if self.__COEFFICIENTS[name].isAlteringRightHandSide(): self.invalidateRightHandSide()
1091       else:       else:
1092          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1093    
1094     def copyConstraint(self,u):     def validSolution(self):
1095        """         """
1096        copies the constraint into u and returns u.         markes the solution as valid
1097           """
1098           self.__is_solution_valid=True
1099    
1100        @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):
1101        @type u: L{Data<escript.Data>}         """
1102        @return: the input u modified by the constraints.         indicates the PDE has to be resolved if the solution is requested
1103        @rtype: L{Data<escript.Data>}         """
1104        @warning: u is altered if it has the appropriate L{FunctionSpace<escript.FunctionSpace>}         self.trace("System will be resolved.")
1105        """         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  
1106    
1107     def setValue(self,**coefficients):     def isSolutionValid(self):
1108           """
1109           returns True is the solution is still valid.
1110           """
1111           return self.__is_solution_valid
1112    
1113       def validOperator(self):
1114           """
1115           markes the  operator as valid
1116           """
1117           self.__is_operator_valid=True
1118    
1119       def invalidateOperator(self):
1120           """
1121           indicates the operator has to be rebuilt next time it is used
1122           """
1123           self.trace("Operator will be rebuilt.")
1124           self.invalidateSolution()
1125           self.__is_operator_valid=False
1126    
1127       def isOperatorValid(self):
1128           """
1129           returns True is the operator is still valid.
1130           """
1131           return self.__is_operator_valid
1132    
1133       def validRightHandSide(self):
1134           """
1135           markes the right hand side as valid
1136           """
1137           self.__is_RHS_valid=True
1138    
1139       def invalidateRightHandSide(self):
1140           """
1141           indicates the right hand side has to be rebuild next time it is used
1142           """
1143           if self.isRightHandSideValid(): self.trace("Right hand side has to be rebuilt.")
1144           self.invalidateSolution()
1145           self.__is_RHS_valid=False
1146    
1147       def isRightHandSideValid(self):
1148           """
1149           returns True is the operator is still valid.
1150           """
1151           return self.__is_RHS_valid
1152    
1153       def invalidateSystem(self):
1154           """
1155           annonced that everthing has to be rebuild:
1156           """
1157           if self.isRightHandSideValid(): self.trace("System has to be rebuilt.")
1158           self.invalidateSolution()
1159           self.invalidateOperator()
1160           self.invalidateRightHandSide()
1161    
1162       def isSystemValid(self):
1163           """
1164           returns true if the system (including solution) is still vaild:
1165           """
1166           return self.isSolutionValid() and self.isOperatorValid() and self.isRightHandSideValid()
1167    
1168       def initializeSystem(self):
1169           """
1170           resets the system
1171           """
1172           self.trace("New System has been created.")
1173           self.__operator=escript.Operator()
1174           self.__righthandside=escript.Data()
1175           self.__solution=escript.Data()
1176           self.invalidateSystem()
1177    
1178       def getOperator(self):
1179         """
1180         returns the operator of the linear problem
1181    
1182         @return: the operator of the problem
1183         """
1184         return self.getSystem()[0]
1185    
1186       def getRightHandSide(self):
1187         """
1188         returns the right hand side of the linear problem
1189    
1190         @return: the right hand side of the problem
1191         @rtype: L{Data<escript.Data>}
1192         """
1193         return self.getSystem()[1]
1194    
1195       def createRightHandSide(self):
1196           """
1197           returns an instance of a new right hand side
1198           """
1199           self.trace("New right hand side is allocated.")
1200           if self.getNumEquations()>1:
1201               return escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)
1202           else:
1203               return escript.Data(0.,(),self.getFunctionSpaceForEquation(),True)
1204    
1205       def createSolution(self):
1206           """
1207           returns an instance of a new solution
1208           """
1209           self.trace("New solution is allocated.")
1210           if self.getNumSolutions()>1:
1211               return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)
1212           else:
1213               return escript.Data(0.,(),self.getFunctionSpaceForSolution(),True)
1214    
1215       def resetSolution(self):
1216           """
1217           sets the solution to zero
1218           """
1219           if self.__solution.isEmpty():
1220               self.__solution=self.createSolution()
1221           else:
1222               self.__solution.setToZero()
1223               self.trace("Solution is reset to zero.")
1224       def setSolution(self,u):
1225           """
1226           sets the solution assuming that makes the sytem valid.
1227           """
1228           self.__solution=u
1229           self.validSolution()
1230    
1231       def getCurrentSolution(self):
1232           """
1233           returns the solution in its current state
1234           """
1235           if self.__solution.isEmpty(): self.__solution=self.createSolution()
1236           return self.__solution
1237    
1238       def resetRightHandSide(self):
1239           """
1240           sets the right hand side to zero
1241           """
1242           if self.__righthandside.isEmpty():
1243               self.__righthandside=self.createRightHandSide()
1244           else:
1245               self.__righthandside.setToZero()
1246               self.trace("Right hand side is reset to zero.")
1247    
1248       def getCurrentRightHandSide(self):
1249           """
1250           returns the right hand side  in its current state
1251           """
1252           if self.__righthandside.isEmpty(): self.__righthandside=self.createRightHandSide()
1253           return self.__righthandside
1254            
1255    
1256       def resetOperator(self):
1257           """
1258           makes sure that the operator is instantiated and returns it initialized by zeros
1259           """
1260           if self.__operator.isEmpty():
1261               if self.isUsingLumping():
1262                   self.__operator=self.createSolution()
1263               else:
1264                   self.__operator=self.createOperator()
1265           else:
1266               if self.isUsingLumping():
1267                   self.__operator.setToZero()
1268               else:
1269                   self.__operator.resetValues()
1270               self.trace("Operator reset to zero")
1271    
1272       def getCurrentOperator(self):
1273           """
1274           returns the operator in its current state
1275           """
1276           if self.__operator.isEmpty(): self.resetOperator()
1277           return self.__operator
1278    
1279       def setValue(self,**coefficients):
1280        """        """
1281        sets new values to coefficients        sets new values to coefficients
1282    
       @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.  
1283        @raise IllegalCoefficient: if an unknown coefficient keyword is used.        @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1284        """        """
1285        # check if the coefficients are  legal:        # check if the coefficients are  legal:
# Line 1488  class LinearPDE(object): Line 1297  class LinearPDE(object):
1297                  s=numarray.array(d).shape                  s=numarray.array(d).shape
1298              if s!=None:              if s!=None:
1299                  # get number of equations and number of unknowns:                  # get number of equations and number of unknowns:
1300                  res=self.COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(self.getDomain(),s)                  res=self.__COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(self.getDomain(),s)
1301                  if res==None:                  if res==None:
1302                      raise IllegalCoefficientValue,"Illegal shape %s of coefficient %s"%(s,i)                      raise IllegalCoefficientValue,"Illegal shape %s of coefficient %s"%(s,i)
1303                  else:                  else:
# Line 1499  class LinearPDE(object): Line 1308  class LinearPDE(object):
1308        # 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:
1309        for i,d in coefficients.iteritems():        for i,d in coefficients.iteritems():
1310          try:          try:
1311             self.COEFFICIENTS[i].setValue(self.getDomain(),self.getNumEquations(),self.getNumSolutions(),self.reduceEquationOrder(),self.reduceSolutionOrder(),d)             self.__COEFFICIENTS[i].setValue(self.getDomain(),
1312                                             self.getNumEquations(),self.getNumSolutions(),
1313                                             self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
1314               self.alteredCoefficient(i)
1315            except IllegalCoefficientFunctionSpace,m:
1316                # if the function space is wrong then we try the reduced version:
1317                i_red=i+"_reduced"
1318                if (not i_red in coefficients.keys()) and i_red in self.__COEFFICIENTS.keys():
1319                    try:
1320                        self.__COEFFICIENTS[i_red].setValue(self.getDomain(),
1321                                                          self.getNumEquations(),self.getNumSolutions(),
1322                                                          self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
1323                        self.alteredCoefficient(i_red)
1324                    except IllegalCoefficientValue,m:
1325                        raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
1326                    except IllegalCoefficientFunctionSpace,m:
1327                        raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m))
1328                else:
1329                    raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m))
1330          except IllegalCoefficientValue,m:          except IllegalCoefficientValue,m:
1331             raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))             raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
         self.alteredCoefficient(i)  
   
1332        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()  
1333    
1334     def getSystem(self):     # ====================================================================================
1335       # methods that are typically overwritten when implementing a particular linear problem
1336       # ====================================================================================
1337       def getRequiredSystemType(self):
1338          """
1339          returns the system type which needs to be used by the current set up.
1340    
1341          @note: Typically this method is overwritten when implementing a particular linear problem.
1342          """
1343          return 0
1344    
1345       def createOperator(self):
1346           """
1347           returns an instance of a new operator
1348    
1349           @note: This method is overwritten when implementing a particular linear problem.
1350           """
1351           return escript.Operator()
1352    
1353       def checkSymmetry(self,verbose=True):
1354          """
1355          test the PDE for symmetry.
1356    
1357          @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.
1358          @type verbose: C{bool}
1359          @return:  True if the problem is symmetric.
1360          @rtype: C{bool}
1361          @note: Typically this method is overwritten when implementing a particular linear problem.
1362          """
1363          out=True
1364          return out
1365    
1366       def getSolution(self,**options):
1367           """
1368           returns the solution of the problem
1369           @return: the solution
1370           @rtype: L{Data<escript.Data>}
1371    
1372           @note: This method is overwritten when implementing a particular linear problem.
1373           """
1374           return self.getCurrentSolution()
1375    
1376       def getSystem(self):
1377         """         """
1378         return the operator and right hand side of the PDE         return the operator and right hand side of the PDE
1379    
1380         @return: the discrete version of the PDE         @return: the discrete version of the PDE
1381         @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>}.
1382    
1383           @note: This method is overwritten when implementing a particular linear problem.
1384         """         """
1385         if not self.__operator_is_Valid or not self.__righthandside_isValid:         return (self.getCurrentOperator(), self.getCurrentRightHandSide())
1386    #=====================
1387    
1388    class LinearPDE(LinearProblem):
1389       """
1390       This class is used to define a general linear, steady, second order PDE
1391       for an unknown function M{u} on a given domain defined through a L{Domain<escript.Domain>} object.
1392    
1393       For a single PDE with a solution with a single component the linear PDE is defined in the following form:
1394    
1395       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)}
1396    
1397    
1398       where M{grad(F)} denotes the spatial derivative of M{F}. Einstein's summation convention,
1399       ie. summation over indexes appearing twice in a term of a sum is performed, is used.
1400       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
1401       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
1402       L{ReducedFunction<escript.ReducedFunction>}. It is also allowd to use objects that can be converted into
1403       such L{Data<escript.Data>} objects. M{A} and M{A_reduced} are rank two, M{B}, M{C}, M{X},
1404       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.
1405    
1406       The following natural boundary conditions are considered:
1407    
1408       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}
1409    
1410       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>}.
1411    
1412    
1413       Constraints for the solution prescribing the value of the solution at certain locations in the domain. They have the form
1414    
1415       M{u=r}  where M{q>0}
1416    
1417       M{r} and M{q} are each scalar where M{q} is the characteristic function defining where the constraint is applied.
1418       The constraints override any other condition set by the PDE or the boundary condition.
1419    
1420       The PDE is symmetrical if
1421    
1422       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]}
1423    
1424       For a system of PDEs and a solution with several components the PDE has the form
1425    
1426       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] }
1427    
1428       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.
1429       The natural boundary conditions take the form:
1430    
1431       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]}
1432    
1433    
1434       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>}.
1435    
1436       Constraints take the form
1437    
1438       M{u[i]=r[i]}  where  M{q[i]>0}
1439    
1440       M{r} and M{q} are each rank one. Notice that at some locations not necessarily all components must have a constraint.
1441    
1442       The system of PDEs is symmetrical if
1443    
1444            - M{A[i,j,k,l]=A[k,l,i,j]}
1445            - M{A_reduced[i,j,k,l]=A_reduced[k,l,i,j]}
1446            - M{B[i,j,k]=C[k,i,j]}
1447            - M{B_reduced[i,j,k]=C_reduced[k,i,j]}
1448            - M{D[i,k]=D[i,k]}
1449            - M{D_reduced[i,k]=D_reduced[i,k]}
1450            - M{d[i,k]=d[k,i]}
1451            - M{d_reduced[i,k]=d_reduced[k,i]}
1452    
1453       L{LinearPDE} also supports solution discontinuities over a contact region in the domain. To specify the conditions across the
1454       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
1455       defined as
1456    
1457       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]}
1458    
1459       For the case of single solution component and single PDE M{J} is defined
1460    
1461       M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[j]+(B[i]+B_reduced[i])*u-X[i]-X_reduced[i]}
1462    
1463       In the context of discontinuities M{n} denotes the normal on the discontinuity pointing from side 0 towards side 1
1464       calculated from L{getNormal<escript.FunctionSpace.getNormal>} of L{FunctionOnContactZero<escript.FunctionOnContactZero>}. For a system of PDEs
1465       the contact condition takes the form
1466    
1467       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]}
1468    
1469       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
1470       of the solution at side 1 and at side 0, denotes the jump of M{u} across discontinuity along the normal calculated by
1471       L{jump<util.jump>}.
1472       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>}.
1473       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>}.
1474       In case of a single PDE and a single component solution the contact condition takes the form
1475    
1476       M{n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)}
1477    
1478       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>}
1479    
1480       typical usage:
1481    
1482          p=LinearPDE(dom)
1483          p.setValue(A=kronecker(dom),D=1,Y=0.5)
1484          u=p.getSolution()
1485    
1486       """
1487    
1488    
1489       def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):
1490         """
1491         initializes a new linear PDE
1492    
1493         @param domain: domain of the PDE
1494         @type domain: L{Domain<escript.Domain>}
1495         @param numEquations: number of equations. If numEquations==None the number of equations
1496                              is extracted from the PDE coefficients.
1497         @param numSolutions: number of solution components. If  numSolutions==None the number of solution components
1498                              is extracted from the PDE coefficients.
1499         @param debug: if True debug informations are printed.
1500    
1501         """
1502         super(LinearPDE, self).__init__(domain,numEquations,numSolutions,debug)
1503         #
1504         #   the coefficients of the PDE:
1505         #
1506         self.introduceCoefficients(
1507           A=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
1508           B=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1509           C=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
1510           D=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1511           X=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
1512           Y=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1513           d=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1514           y=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1515           d_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1516           y_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1517           A_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
1518           B_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1519           C_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
1520           D_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1521           X_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
1522           Y_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1523           d_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1524           y_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1525           d_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1526           y_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1527           r=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.RIGHTHANDSIDE),
1528           q=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.BOTH) )
1529    
1530       def __str__(self):
1531         """
1532         returns string representation of the PDE
1533    
1534         @return: a simple representation of the PDE
1535         @rtype: C{str}
1536         """
1537         return "<LinearPDE %d>"%id(self)
1538    
1539       def getRequiredSystemType(self):
1540          """
1541          returns the system type which needs to be used by the current set up.
1542          """
1543          return self.getDomain().getSystemMatrixTypeId(self.getSolverMethod()[0],self.getSolverPackage(),self.isSymmetric())
1544    
1545       def checkSymmetry(self,verbose=True):
1546          """
1547          test the PDE for symmetry.
1548    
1549          @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.
1550          @type verbose: C{bool}
1551          @return:  True if the PDE is symmetric.
1552          @rtype: L{bool}
1553          @note: This is a very expensive operation. It should be used for degugging only! The symmetry flag is not altered.
1554          """
1555          out=True
1556          out=out and self.checkSymmetricTensor("A", verbose)
1557          out=out and self.checkSymmetricTensor("A_reduced", verbose)
1558          out=out and self.checkReciprocalSymmetry("B","C", verbose)
1559          out=out and self.checkReciprocalSymmetry("B_reduced","C_reduced", verbose)
1560          out=out and self.checkSymmetricTensor("D", verbose)
1561          out=out and self.checkSymmetricTensor("D_reduced", verbose)
1562          out=out and self.checkSymmetricTensor("d", verbose)
1563          out=out and self.checkSymmetricTensor("d_reduced", verbose)
1564          out=out and self.checkSymmetricTensor("d_contact", verbose)
1565          out=out and self.checkSymmetricTensor("d_contact_reduced", verbose)
1566          return out
1567    
1568       def createOperator(self):
1569           """
1570           returns an instance of a new operator
1571           """
1572           self.trace("New operator is allocated.")
1573           return self.getDomain().newOperator( \
1574                               self.getNumEquations(), \
1575                               self.getFunctionSpaceForEquation(), \
1576                               self.getNumSolutions(), \
1577                               self.getFunctionSpaceForSolution(), \
1578                               self.getSystemType())
1579    
1580       def getSolution(self,**options):
1581           """
1582           returns the solution of the PDE
1583    
1584           @return: the solution
1585           @rtype: L{Data<escript.Data>}
1586           @param options: solver options
1587           @keyword verbose: True to get some information during PDE solution
1588           @type verbose: C{bool}
1589           @keyword reordering: reordering scheme to be used during elimination. Allowed values are
1590                                L{NO_REORDERING}, L{MINIMUM_FILL_IN}, L{NESTED_DISSECTION}
1591           @keyword iter_max: maximum number of iteration steps allowed.
1592           @keyword drop_tolerance: threshold for drupping in L{ILUT}
1593           @keyword drop_storage: maximum of allowed memory in L{ILUT}
1594           @keyword truncation: maximum number of residuals in L{GMRES}
1595           @keyword restart: restart cycle length in L{GMRES}
1596           """
1597           if not self.isSolutionValid():
1598              mat,f=self.getSystem()
1599              if self.isUsingLumping():
1600                 self.setSolution(f*1/mat)
1601              else:
1602                 options[self.TOLERANCE_KEY]=self.getTolerance()
1603                 options[self.METHOD_KEY]=self.getSolverMethod()[0]
1604                 options[self.PRECONDITIONER_KEY]=self.getSolverMethod()[1]
1605                 options[self.PACKAGE_KEY]=self.getSolverPackage()
1606                 options[self.SYMMETRY_KEY]=self.isSymmetric()
1607                 self.trace("PDE is resolved.")
1608                 self.trace("solver options: %s"%str(options))
1609                 self.setSolution(mat.solve(f,options))
1610           return self.getCurrentSolution()
1611    
1612       def getSystem(self):
1613           """
1614           return the operator and right hand side of the PDE
1615    
1616           @return: the discrete version of the PDE
1617           @rtype: C{tuple} of L{Operator,<escript.Operator>} and L{Data<escript.Data>}.
1618           """
1619           if not self.isOperatorValid() or not self.isRightHandSideValid():
1620            if self.isUsingLumping():            if self.isUsingLumping():
1621                if not self.__operator_is_Valid:                if not self.isOperatorValid():
1622                   if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution(): raise TypeError,"Lumped matrix requires same order for equations and unknowns"                   if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():
1623                   if not self.getCoefficientOfGeneralPDE("A").isEmpty(): raise Warning,"Using coefficient A in lumped matrix can produce wrong results"                        raise TypeError,"Lumped matrix requires same order for equations and unknowns"
1624                   if not self.getCoefficientOfGeneralPDE("B").isEmpty(): raise Warning,"Using coefficient B in lumped matrix can produce wrong results"                   if not self.getCoefficient("A").isEmpty():
1625                   if not self.getCoefficientOfGeneralPDE("C").isEmpty(): raise Warning,"Using coefficient C in lumped matrix can produce wrong results"                        raise ValueError,"coefficient A in lumped matrix may not be present."
1626                   mat=self.__getNewOperator()                   if not self.getCoefficient("B").isEmpty():
1627                   self.getDomain().addPDEToSystem(mat,escript.Data(), \                        raise ValueError,"coefficient B in lumped matrix may not be present."
1628                             self.getCoefficientOfGeneralPDE("A"), \                   if not self.getCoefficient("C").isEmpty():
1629                             self.getCoefficientOfGeneralPDE("B"), \                        raise ValueError,"coefficient C in lumped matrix may not be present."
1630                             self.getCoefficientOfGeneralPDE("C"), \                   if not self.getCoefficient("d_contact").isEmpty():
1631                             self.getCoefficientOfGeneralPDE("D"), \                        raise ValueError,"coefficient d_contact in lumped matrix may not be present."
1632                             escript.Data(), \                   if not self.getCoefficient("A_reduced").isEmpty():
1633                             escript.Data(), \                        raise ValueError,"coefficient A_reduced in lumped matrix may not be present."
1634                             self.getCoefficientOfGeneralPDE("d"), \                   if not self.getCoefficient("B_reduced").isEmpty():
1635                             escript.Data(),\                        raise ValueError,"coefficient B_reduced in lumped matrix may not be present."
1636                             self.getCoefficientOfGeneralPDE("d_contact"), \                   if not self.getCoefficient("C_reduced").isEmpty():
1637                             escript.Data())                        raise ValueError,"coefficient C_reduced in lumped matrix may not be present."
1638                   self.__operator=1./(mat*escript.Data(1,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True))                   if not self.getCoefficient("d_contact_reduced").isEmpty():
1639                   del mat                        raise ValueError,"coefficient d_contact_reduced in lumped matrix may not be present."
1640                     D=self.getCoefficient("D")
1641                     d=self.getCoefficient("d")
1642                     D_reduced=self.getCoefficient("D_reduced")
1643                     d_reduced=self.getCoefficient("d_reduced")
1644                     if not D.isEmpty():
1645                         if self.getNumSolutions()>1:
1646                            D_times_e=util.matrix_mult(D,numarray.ones((self.getNumSolutions(),)))
1647                         else:
1648                            D_times_e=D
1649                     else:
1650                        D_times_e=escript.Data()
1651                     if not d.isEmpty():
1652                         if self.getNumSolutions()>1:
1653                            d_times_e=util.matrix_mult(d,numarray.ones((self.getNumSolutions(),)))
1654                         else:
1655                            d_times_e=d
1656                     else:
1657                        d_times_e=escript.Data()
1658          
1659                     if not D_reduced.isEmpty():
1660                         if self.getNumSolutions()>1:
1661                            D_reduced_times_e=util.matrix_mult(D_reduced,numarray.ones((self.getNumSolutions(),)))
1662                         else:
1663                            D_reduced_times_e=D_reduced
1664                     else:
1665                        D_reduced_times_e=escript.Data()
1666                     if not d_reduced.isEmpty():
1667                         if self.getNumSolutions()>1:
1668                            d_reduced_times_e=util.matrix_mult(d_reduced,numarray.ones((self.getNumSolutions(),)))
1669                         else:
1670                            d_reduced_times_e=d_reduced
1671                     else:
1672                        d_reduced_times_e=escript.Data()
1673    
1674                     self.resetOperator()
1675                     operator=self.getCurrentOperator()
1676                     if False and hasattr(self.getDomain(), "addPDEToLumpedSystem") :
1677                        self.getDomain().addPDEToLumpedSystem(operator, D_times_e, d_times_e)
1678                        self.getDomain().addPDEToLumpedSystem(operator, D_reduced_times_e, d_reduced_times_e)
1679                     else:
1680                        self.getDomain().addPDEToRHS(operator, \
1681                                                     escript.Data(), \
1682                                                     D_times_e, \
1683                                                     d_times_e,\
1684                                                     escript.Data())
1685                        self.getDomain().addPDEToRHS(operator, \
1686                                                     escript.Data(), \
1687                                                     D_reduced_times_e, \
1688                                                     d_reduced_times_e,\
1689                                                     escript.Data())
1690                   self.trace("New lumped operator has been built.")                   self.trace("New lumped operator has been built.")
1691                   self.__operator_is_Valid=True                if not self.isRightHandSideValid():
1692                if not self.__righthandside_isValid:                   self.resetRightHandSide()
1693                   self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \                   righthandside=self.getCurrentRightHandSide()
1694                                 self.getCoefficientOfGeneralPDE("X"), \                   self.getDomain().addPDEToRHS(righthandside, \
1695                                 self.getCoefficientOfGeneralPDE("Y"),\                                 self.getCoefficient("X"), \
1696                                 self.getCoefficientOfGeneralPDE("y"),\                                 self.getCoefficient("Y"),\
1697                                 self.getCoefficientOfGeneralPDE("y_contact"))                                 self.getCoefficient("y"),\
1698                                   self.getCoefficient("y_contact"))
1699                     self.getDomain().addPDEToRHS(righthandside, \
1700                                   self.getCoefficient("X_reduced"), \
1701                                   self.getCoefficient("Y_reduced"),\
1702                                   self.getCoefficient("y_reduced"),\
1703                                   self.getCoefficient("y_contact_reduced"))
1704                   self.trace("New right hand side as been built.")                   self.trace("New right hand side as been built.")
1705                   self.__righthandside_isValid=True                   self.validRightHandSide()
1706                  self.insertConstraint(rhs_only=False)
1707                  self.validOperator()
1708            else:            else:
1709               if not self.__operator_is_Valid and not self.__righthandside_isValid:               if not self.isOperatorValid() and not self.isRightHandSideValid():
1710                   self.getDomain().addPDEToSystem(self.__makeFreshOperator(),self.__makeFreshRightHandSide(), \                   self.resetRightHandSide()
1711                                 self.getCoefficientOfGeneralPDE("A"), \                   righthandside=self.getCurrentRightHandSide()
1712                                 self.getCoefficientOfGeneralPDE("B"), \                   self.resetOperator()
1713                                 self.getCoefficientOfGeneralPDE("C"), \                   operator=self.getCurrentOperator()
1714                                 self.getCoefficientOfGeneralPDE("D"), \                   self.getDomain().addPDEToSystem(operator,righthandside, \
1715                                 self.getCoefficientOfGeneralPDE("X"), \                                 self.getCoefficient("A"), \
1716                                 self.getCoefficientOfGeneralPDE("Y"), \                                 self.getCoefficient("B"), \
1717                                 self.getCoefficientOfGeneralPDE("d"), \                                 self.getCoefficient("C"), \
1718                                 self.getCoefficientOfGeneralPDE("y"), \                                 self.getCoefficient("D"), \
1719                                 self.getCoefficientOfGeneralPDE("d_contact"), \                                 self.getCoefficient("X"), \
1720                                 self.getCoefficientOfGeneralPDE("y_contact"))                                 self.getCoefficient("Y"), \
1721                   self.__applyConstraint()                                 self.getCoefficient("d"), \
1722                   self.__righthandside=self.copyConstraint(self.__righthandside)                                 self.getCoefficient("y"), \
1723                                   self.getCoefficient("d_contact"), \
1724                                   self.getCoefficient("y_contact"))
1725                     self.getDomain().addPDEToSystem(operator,righthandside, \
1726                                   self.getCoefficient("A_reduced"), \
1727                                   self.getCoefficient("B_reduced"), \
1728                                   self.getCoefficient("C_reduced"), \
1729                                   self.getCoefficient("D_reduced"), \
1730                                   self.getCoefficient("X_reduced"), \
1731                                   self.getCoefficient("Y_reduced"), \
1732                                   self.getCoefficient("d_reduced"), \
1733                                   self.getCoefficient("y_reduced"), \
1734                                   self.getCoefficient("d_contact_reduced"), \
1735                                   self.getCoefficient("y_contact_reduced"))
1736                     self.insertConstraint(rhs_only=False)
1737                   self.trace("New system has been built.")                   self.trace("New system has been built.")
1738                   self.__operator_is_Valid=True                   self.validOperator()
1739                   self.__righthandside_isValid=True                   self.validRightHandSide()
1740               elif not self.__righthandside_isValid:               elif not self.isRightHandSideValid():
1741                   self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \                   self.resetRightHandSide()
1742                                 self.getCoefficientOfGeneralPDE("X"), \                   righthandside=self.getCurrentRightHandSide()
1743                                 self.getCoefficientOfGeneralPDE("Y"),\                   self.getDomain().addPDEToRHS(righthandside,
1744                                 self.getCoefficientOfGeneralPDE("y"),\                                 self.getCoefficient("X"), \
1745                                 self.getCoefficientOfGeneralPDE("y_contact"))                                 self.getCoefficient("Y"),\
1746                   self.__righthandside=self.copyConstraint(self.__righthandside)                                 self.getCoefficient("y"),\
1747                                   self.getCoefficient("y_contact"))
1748                     self.getDomain().addPDEToRHS(righthandside,
1749                                   self.getCoefficient("X_reduced"), \
1750                                   self.getCoefficient("Y_reduced"),\
1751                                   self.getCoefficient("y_reduced"),\
1752                                   self.getCoefficient("y_contact_reduced"))
1753                     self.insertConstraint(rhs_only=True)
1754                   self.trace("New right hand side has been built.")                   self.trace("New right hand side has been built.")
1755                   self.__righthandside_isValid=True                   self.validRightHandSide()
1756               elif not self.__operator_is_Valid:               elif not self.isOperatorValid():
1757                   self.getDomain().addPDEToSystem(self.__makeFreshOperator(),escript.Data(), \                   self.resetOperator()
1758                              self.getCoefficientOfGeneralPDE("A"), \                   operator=self.getCurrentOperator()
1759                              self.getCoefficientOfGeneralPDE("B"), \                   self.getDomain().addPDEToSystem(operator,escript.Data(), \
1760                              self.getCoefficientOfGeneralPDE("C"), \                              self.getCoefficient("A"), \
1761                              self.getCoefficientOfGeneralPDE("D"), \                              self.getCoefficient("B"), \
1762                                self.getCoefficient("C"), \
1763                                self.getCoefficient("D"), \
1764                                escript.Data(), \
1765                                escript.Data(), \
1766                                self.getCoefficient("d"), \
1767                                escript.Data(),\
1768                                self.getCoefficient("d_contact"), \
1769                                escript.Data())
1770                     self.getDomain().addPDEToSystem(operator,escript.Data(), \
1771                                self.getCoefficient("A_reduced"), \
1772                                self.getCoefficient("B_reduced"), \
1773                                self.getCoefficient("C_reduced"), \
1774                                self.getCoefficient("D_reduced"), \
1775                              escript.Data(), \                              escript.Data(), \
1776                              escript.Data(), \                              escript.Data(), \
1777                              self.getCoefficientOfGeneralPDE("d"), \                              self.getCoefficient("d_reduced"), \
1778                              escript.Data(),\                              escript.Data(),\
1779                              self.getCoefficientOfGeneralPDE("d_contact"), \                              self.getCoefficient("d_contact_reduced"), \
1780                              escript.Data())                              escript.Data())
1781                   self.__applyConstraint()                   self.insertConstraint(rhs_only=False)
1782                   self.trace("New operator has been built.")                   self.trace("New operator has been built.")
1783                   self.__operator_is_Valid=True                   self.validOperator()
1784         return (self.__operator,self.__righthandside)         return (self.getCurrentOperator(), self.getCurrentRightHandSide())
1785    
1786       def insertConstraint(self, rhs_only=False):
1787          """
1788          applies the constraints defined by q and r to the PDE
1789          
1790          @param rhs_only: if True only the right hand side is altered by the constraint.
1791          @type rhs_only: C{bool}
1792          """
1793          q=self.getCoefficient("q")
1794          r=self.getCoefficient("r")
1795          righthandside=self.getCurrentRightHandSide()
1796          operator=self.getCurrentOperator()
1797    
1798          if not q.isEmpty():
1799             if r.isEmpty():
1800                r_s=self.createSolution()
1801             else:
1802                r_s=r
1803             if not rhs_only and not operator.isEmpty():
1804                 if self.isUsingLumping():
1805                     operator.copyWithMask(escript.Data(1.,q.getShape(),q.getFunctionSpace()),q)
1806                 else:
1807                     row_q=escript.Data(q,self.getFunctionSpaceForEquation())
1808                     col_q=escript.Data(q,self.getFunctionSpaceForSolution())
1809                     u=self.createSolution()
1810                     u.copyWithMask(r_s,col_q)
1811                     righthandside-=operator*u
1812                     operator.nullifyRowsAndCols(row_q,col_q,1.)
1813             righthandside.copyWithMask(r_s,q)
1814    
1815       def setValue(self,**coefficients):
1816          """
1817          sets new values to coefficients
1818    
1819          @param coefficients: new values assigned to coefficients
1820          @keyword A: value for coefficient A.
1821          @type A: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1822          @keyword A_reduced: value for coefficient A_reduced.
1823          @type A_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1824          @keyword B: value for coefficient B
1825          @type B: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1826          @keyword B_reduced: value for coefficient B_reduced
1827          @type B_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1828          @keyword C: value for coefficient C
1829          @type C: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1830          @keyword C_reduced: value for coefficient C_reduced
1831          @type C_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1832          @keyword D: value for coefficient D
1833          @type D: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1834          @keyword D_reduced: value for coefficient D_reduced
1835          @type D_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1836          @keyword X: value for coefficient X
1837          @type X: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1838          @keyword X_reduced: value for coefficient X_reduced
1839          @type X_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1840          @keyword Y: value for coefficient Y
1841          @type Y: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1842          @keyword Y_reduced: value for coefficient Y_reduced
1843          @type Y_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.Function>}.
1844          @keyword d: value for coefficient d
1845          @type d: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1846          @keyword d_reduced: value for coefficient d_reduced
1847          @type d_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
1848          @keyword y: value for coefficient y
1849          @type y: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1850          @keyword d_contact: value for coefficient d_contact
1851          @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>}.
1852          @keyword d_contact_reduced: value for coefficient d_contact_reduced
1853          @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>}.
1854          @keyword y_contact: value for coefficient y_contact
1855          @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>}.
1856          @keyword y_contact_reduced: value for coefficient y_contact_reduced
1857          @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>}.
1858          @keyword r: values prescribed to the solution at the locations of constraints
1859          @type r: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1860                   depending of reduced order is used for the solution.
1861          @keyword q: mask for location of constraints
1862          @type q: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1863                   depending of reduced order is used for the representation of the equation.
1864          @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1865          """
1866          super(LinearPDE,self).setValue(**coefficients)
1867          # check if the systrem is inhomogeneous:
1868          if len(coefficients)>0 and not self.isUsingLumping():
1869             q=self.getCoefficient("q")
1870             r=self.getCoefficient("r")
1871             if not q.isEmpty() and not r.isEmpty():
1872                 if util.Lsup(q*r)>0.:
1873                   self.trace("Inhomogeneous constraint detected.")
1874                   self.invalidateSystem()
1875    
1876    
1877       def getResidual(self,u=None):
1878         """
1879         return the residual of u or the current solution if u is not present.
1880    
1881         @param u: argument in the residual calculation. It must be representable in C{elf.getFunctionSpaceForSolution()}. If u is not present or equals L{None}
1882                   the current solution is used.
1883         @type u: L{Data<escript.Data>} or None
1884         @return: residual of u
1885         @rtype: L{Data<escript.Data>}
1886         """
1887         if u==None:
1888            return self.getOperator()*self.getSolution()-self.getRightHandSide()
1889         else:
1890            return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())-self.getRightHandSide()
1891    
1892       def getFlux(self,u=None):
1893         """
1894         returns the flux M{J} for a given M{u}
1895    
1896         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]}
1897    
1898         or
1899    
1900         M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[l]+(B[j]+B_reduced[j])u-X[j]-X_reduced[j]}
1901    
1902         @param u: argument in the flux. If u is not present or equals L{None} the current solution is used.
1903         @type u: L{Data<escript.Data>} or None
1904         @return: flux
1905         @rtype: L{Data<escript.Data>}
1906         """
1907         if u==None: u=self.getSolution()
1908         return util.tensormult(self.getCoefficient("A"),util.grad(u,Funtion(self.getDomain))) \
1909               +util.matrixmult(self.getCoefficient("B"),u) \
1910               -util.self.getCoefficient("X") \
1911               +util.tensormult(self.getCoefficient("A_reduced"),util.grad(u,ReducedFuntion(self.getDomain))) \
1912               +util.matrixmult(self.getCoefficient("B_reduced"),u) \
1913               -util.self.getCoefficient("X_reduced")
1914    
1915    
1916  class Poisson(LinearPDE):  class Poisson(LinearPDE):
# Line 1624  class Poisson(LinearPDE): Line 1939  class Poisson(LinearPDE):
1939    
1940       """       """
1941       super(Poisson, self).__init__(domain,1,1,debug)       super(Poisson, self).__init__(domain,1,1,debug)
1942       self.COEFFICIENTS={"f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),       self.introduceCoefficients(
1943                            "q": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}                          f=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1944                            f_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE))
1945       self.setSymmetryOn()       self.setSymmetryOn()
1946    
1947     def setValue(self,**coefficients):     def setValue(self,**coefficients):
# Line 1642  class Poisson(LinearPDE): Line 1958  class Poisson(LinearPDE):
1958       """       """
1959       super(Poisson, self).setValue(**coefficients)       super(Poisson, self).setValue(**coefficients)
1960    
1961     def getCoefficientOfGeneralPDE(self,name):  
1962       def getCoefficient(self,name):
1963       """       """
1964       return the value of the coefficient name of the general PDE       return the value of the coefficient name of the general PDE
1965       @param name: name of the coefficient requested.       @param name: name of the coefficient requested.
1966       @type name: C{string}       @type name: C{string}
1967       @return: the value of the coefficient  name       @return: the value of the coefficient  name
1968       @rtype: L{Data<escript.Data>}       @rtype: L{Data<escript.Data>}
1969       @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}.  
1970       @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.
1971       """       """
1972       if name == "A" :       if name == "A" :
1973           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()  
1974       elif name == "Y" :       elif name == "Y" :
1975           return self.getCoefficient("f")           return self.getCoefficient("f")
1976       elif name == "d" :       elif name == "Y_reduced" :
1977           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")  
1978       else:       else:
1979          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name          return super(Poisson, self).getCoefficient(name)
1980    
1981  class Helmholtz(LinearPDE):  class Helmholtz(LinearPDE):
1982     """     """
# Line 1706  class Helmholtz(LinearPDE): Line 2004  class Helmholtz(LinearPDE):
2004    
2005       """       """
2006       super(Helmholtz, self).__init__(domain,1,1,debug)       super(Helmholtz, self).__init__(domain,1,1,debug)
2007       self.COEFFICIENTS={"omega": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),       self.introduceCoefficients(
2008                          "k": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),                          omega=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.OPERATOR),
2009                          "f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),                          k=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.OPERATOR),
2010                          "alpha": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),                          f=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2011                          "g": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),                          f_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2012                          "r": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH),                          alpha=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.OPERATOR),
2013                          "q": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}                          g=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2014                            g_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE))
2015       self.setSymmetryOn()       self.setSymmetryOn()
2016    
2017     def setValue(self,**coefficients):     def setValue(self,**coefficients):
2018       """       """
2019       sets new values to coefficients       sets new values to coefficients
2020    
# Line 1740  class Helmholtz(LinearPDE): Line 2039  class Helmholtz(LinearPDE):
2039       """       """
2040       super(Helmholtz, self).setValue(**coefficients)       super(Helmholtz, self).setValue(**coefficients)
2041    
2042     def getCoefficientOfGeneralPDE(self,name):     def getCoefficient(self,name):
2043       """       """
2044       return the value of the coefficient name of the general PDE       return the value of the coefficient name of the general PDE
2045    
# Line 1748  class Helmholtz(LinearPDE): Line 2047  class Helmholtz(LinearPDE):
2047       @type name: C{string}       @type name: C{string}
2048       @return: the value of the coefficient  name       @return: the value of the coefficient  name
2049       @rtype: L{Data<escript.Data>}       @rtype: L{Data<escript.Data>}
2050       @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.  
2051       """       """
2052       if name == "A" :       if name == "A" :
2053           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()  
2054       elif name == "D" :       elif name == "D" :
2055           return self.getCoefficient("omega")           return self.getCoefficient("omega")
      elif name == "X" :  
          return escript.Data()  
2056       elif name == "Y" :       elif name == "Y" :
2057           return self.getCoefficient("f")           return self.getCoefficient("f")
2058       elif name == "d" :       elif name == "d" :
2059           return self.getCoefficient("alpha")           return self.getCoefficient("alpha")
2060       elif name == "y" :       elif name == "y" :
2061           return self.getCoefficient("g")           return self.getCoefficient("g")
2062       elif name == "d_contact" :       elif name == "Y_reduced" :
2063           return escript.Data()           return self.getCoefficient("f_reduced")
2064       elif name == "y_contact" :       elif name == "y_reduced" :
2065           return escript.Data()          return self.getCoefficient("g_reduced")
      elif name == "r" :  
          return self.getCoefficient("r")  
      elif name == "q" :  
          return self.getCoefficient("q")  
2066       else:       else:
2067          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name          return super(Helmholtz, self).getCoefficient(name)
2068    
2069  class LameEquation(LinearPDE):  class LameEquation(LinearPDE):
2070     """     """
2071     Class to define a Lame equation problem:     Class to define a Lame equation problem:
2072    
2073     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] }
2074    
2075     with natural boundary conditons:     with natural boundary conditons:
2076    
2077     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] }
2078    
2079     and constraints:     and constraints:
2080    
# Line 1798  class LameEquation(LinearPDE): Line 2085  class LameEquation(LinearPDE):
2085     def __init__(self,domain,debug=False):     def __init__(self,domain,debug=False):
2086        super(LameEquation, self).__init__(domain,\        super(LameEquation, self).__init__(domain,\
2087                                           domain.getDim(),domain.getDim(),debug)                                           domain.getDim(),domain.getDim(),debug)
2088        self.COEFFICIENTS={ "lame_lambda"  : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),        self.introduceCoefficients(lame_lambda=PDECoef(PDECoef.INTERIOR,(),PDECoef.OPERATOR),
2089                            "lame_mu"      : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),                                   lame_mu=PDECoef(PDECoef.INTERIOR,(),PDECoef.OPERATOR),
2090                            "F"            : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),                                   F=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2091                            "sigma"        : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM),PDECoefficient.RIGHTHANDSIDE),                                   sigma=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
2092                            "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)}  
2093        self.setSymmetryOn()        self.setSymmetryOn()
2094    
2095     def setValue(self,**coefficients):     def setValues(self,**coefficients):
2096       """       """
2097       sets new values to coefficients       sets new values to coefficients
2098    
# Line 1830  class LameEquation(LinearPDE): Line 2115  class LameEquation(LinearPDE):
2115                 depending of reduced order is used for the representation of the equation.                 depending of reduced order is used for the representation of the equation.
2116       @raise IllegalCoefficient: if an unknown coefficient keyword is used.       @raise IllegalCoefficient: if an unknown coefficient keyword is used.
2117       """       """
2118       super(LameEquation, self).setValue(**coefficients)       super(LameEquation, self).setValues(**coefficients)
2119    
2120     def getCoefficientOfGeneralPDE(self,name):     def getCoefficient(self,name):
2121       """       """
2122       return the value of the coefficient name of the general PDE       return the value of the coefficient name of the general PDE
2123    
# Line 1840  class LameEquation(LinearPDE): Line 2125  class LameEquation(LinearPDE):
2125       @type name: C{string}       @type name: C{string}
2126       @return: the value of the coefficient  name       @return: the value of the coefficient  name
2127       @rtype: L{Data<escript.Data>}       @rtype: L{Data<escript.Data>}
2128       @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.  
2129       """       """
2130       if name == "A" :       if name == "A" :
2131           out =self.createCoefficientOfGeneralPDE("A")           out =self.createCoefficient("A")
2132           for i in range(self.getDim()):           for i in range(self.getDim()):
2133             for j in range(self.getDim()):             for j in range(self.getDim()):
2134               out[i,i,j,j] += self.getCoefficient("lame_lambda")               out[i,i,j,j] += self.getCoefficient("lame_lambda")
2135               out[i,j,j,i] += self.getCoefficient("lame_mu")               out[i,j,j,i] += self.getCoefficient("lame_mu")
2136               out[i,j,i,j] += self.getCoefficient("lame_mu")               out[i,j,i,j] += self.getCoefficient("lame_mu")
2137           return out           return out
      elif name == "B" :  
          return escript.Data()  
      elif name == "C" :  
          return escript.Data()  
      elif name == "D" :  
          return escript.Data()  
2138       elif name == "X" :       elif name == "X" :
2139           return self.getCoefficient("sigma")           return self.getCoefficient("sigma")
2140       elif name == "Y" :       elif name == "Y" :
2141           return self.getCoefficient("F")           return self.getCoefficient("F")
      elif name == "d" :  
          return escript.Data()  
2142       elif name == "y" :       elif name == "y" :
2143           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")  
2144       else:       else:
2145          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name          return super(LameEquation, self).getCoefficient(name)
2146    
2147  class AdvectivePDE(LinearPDE):  def LinearSinglePDE(domain,debug=False):
2148     """     """
2149     In cases of PDEs dominated by the advection terms M{B} and M{C} against the adevctive terms M{A}     defines a single linear PDEs
    up-winding has been used.  The L{AdvectivePDE} class applies SUPG upwinding to the advective terms.  
2150    
2151     In the following we set     @param domain: domain of the PDE
2152       @type domain: L{Domain<escript.Domain>}
2153       @param debug: if True debug informations are printed.
2154       @rtype: L{LinearPDE}
2155       """
2156       return LinearPDE(domain,numEquations=1,numSolutions=1,debug=debug)
2157    
2158     M{Z[j]=C[j]-B[j]}  def LinearPDESystem(domain,debug=False):
2159       """
2160       defines a system of linear PDEs
2161    
2162     or     @param domain: domain of the PDE
2163       @type domain: L{Domain<escript.Domain>}
2164       @param debug: if True debug informations are printed.
2165       @rtype: L{LinearPDE}
2166       """
2167       return LinearPDE(domain,numEquations=domain.getDim(),numSolutions=domain.getDim(),debug=debug)
2168    
2169    class TransportPDE(LinearProblem):
2170       """
2171       This class is used to define a transport problem given by a general linear, time dependend, second order PDE
2172       for an unknown, non-negative, time-dependend function M{u} on a given domain defined through a L{Domain<escript.Domain>} object.
2173    
2174       For a single equation with a solution with a single component the transport problem is defined in the following form:
2175    
2176     M{Z[i,k,l]=C[i,k,l]-B[i,l,k]}     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)}
2177    
    To measure the dominance of the advective terms over the diffusive term M{A} the  
    X{Pelclet number} M{P} is used. It is defined as  
2178    
2179     M{P=h|Z|/(2|A|)}     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,
2180       ie. summation over indexes appearing twice in a term of a sum is performed, is used.
2181       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
2182       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
2183       L{ReducedFunction<escript.ReducedFunction>}. It is also allowed to use objects that can be converted into
2184       such L{Data<escript.Data>} objects. M{M} and M{M_reduced} are scalar, M{A} and M{A_reduced} are rank two,
2185       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.
2186    
2187       The following natural boundary conditions are considered:
2188    
2189       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}
2190    
2191       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>}.
2192    
2193       Constraints for the solution prescribing the value of the solution at certain locations in the domain. They have the form
2194    
2195       M{u_t=r}  where M{q>0}
2196    
2197       M{r} and M{q} are each scalar where M{q} is the characteristic function defining where the constraint is applied.
2198       The constraints override any other condition set by the transport problem or the boundary condition.
2199    
2200       The transport problem is symmetrical if
2201    
2202     where M{|.|} denotes the L{length<util.length>} of the arument and M{h} is the local cell size     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]}
    from L{getSize<escript.Domain.getSize>}. Where M{|A|==0} M{P} is M{S{infinity}}.  
2203    
2204     From the X{Pelclet number} the stabilization parameters M{S{Xi}} and M{S{Xi}} are calculated:     For a system and a solution with several components the transport problem has the form
2205    
2206     M{S{Xi}=S{xi}(P) h/|Z|}     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] }
2207    
2208     where M{S{xi}} is a suitable function of the Peclet number.     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:
2209    
2210     In the case of a single PDE the coefficient are up-dated in the following way:     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}
          - 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}  
2211    
2212     Similar for the case of a systems of PDEs:     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>}.
          - 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]}  
2213    
2214     where M{S{delta}} is L{kronecker}.     Constraints take the form
2215     Using upwinding in this form, introduces an additonal error which is proprtional to the cell size M{h}  
2216     but with the intension to stabilize the solution.     M{u[i]_t=r[i]}  where  M{q[i]>0}
2217    
2218       M{r} and M{q} are each rank one. Notice that at some locations not necessarily all components must have a constraint.
2219    
2220       The transport problem is symmetrical if
2221    
2222            - M{M[i,k]=M[i,k]}
2223            - M{M_reduced[i,k]=M_reduced[i,k]}
2224            - M{A[i,j,k,l]=A[k,l,i,j]}
2225            - M{A_reduced[i,j,k,l]=A_reduced[k,l,i,j]}
2226            - M{B[i,j,k]=C[k,i,j]}
2227            - M{B_reduced[i,j,k]=C_reduced[k,i,j]}
2228            - M{D[i,k]=D[i,k]}
2229            - M{D_reduced[i,k]=D_reduced[i,k]}
2230            - M{m[i,k]=m[k,i]}
2231            - M{m_reduced[i,k]=m_reduced[k,i]}
2232            - M{d[i,k]=d[k,i]}
2233            - M{d_reduced[i,k]=d_reduced[k,i]}
2234    
2235       L{TransportPDE} also supports solution discontinuities over a contact region in the domain. To specify the conditions across the
2236       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
2237       defined as
2238    
2239       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]}
2240    
2241       For the case of single solution component and single PDE M{J} is defined
2242    
2243       M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[j]+(B[i]+B_reduced[i])*u+X[i]+X_reduced[i]}
2244    
2245       In the context of discontinuities M{n} denotes the normal on the discontinuity pointing from side 0 towards side 1
2246       calculated from L{getNormal<escript.FunctionSpace.getNormal>} of L{FunctionOnContactZero<escript.FunctionOnContactZero>}. For a system of transport problems the contact condition takes the form
2247    
2248       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]}
2249    
2250       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
2251       of the solution at side 1 and at side 0, denotes the jump of M{u} across discontinuity along the normal calculated by
2252       L{jump<util.jump>}.
2253       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>}.
2254       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>}.
2255       In case of a single PDE and a single component solution the contact condition takes the form
2256    
2257       M{n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)}
2258    
2259       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>}
2260    
2261       typical usage:
2262    
2263          p=TransportPDE(dom)
2264          p.setValue(M=1.,C=[-1.,0.])
2265          p.setInitialSolution(u=exp(-length(dom.getX()-[0.1,0.1])**2)
2266          t=0
2267          dt=0.1
2268          while (t<1.):
2269              u=p.solve(dt)
2270    
2271     """     """
2272     def __init__(self,domain,numEquations=None,numSolutions=None,xi=None,debug=False):     def __init__(self,domain,num_equations=None,numSolutions=None, theta=0.5,debug=False):
2273         """
2274         initializes a linear problem
2275    
2276         @param domain: domain of the PDE
2277         @type domain: L{Domain<escript.Domain>}
2278         @param numEquations: number of equations. If numEquations==None the number of equations
2279                              is extracted from the coefficients.
2280         @param numSolutions: number of solution components. If  numSolutions==None the number of solution components
2281                              is extracted from the coefficients.
2282         @param debug: if True debug informations are printed.
2283         @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>}.
2284    
2285         """
2286         if theta<0 or theta>1:
2287                  raise ValueError,"theta (=%s) needs to be between 0 and 1 (inclusive)."
2288         super(TransportPDE, self).__init__(domain,num_equations,numSolutions,debug)
2289    
2290         self.__theta=theta
2291         #
2292         #   the coefficients of the transport problem
2293         #
2294         self.introduceCoefficients(
2295           M=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2296           A=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2297           B=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2298           C=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2299           D=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2300           X=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
2301           Y=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2302           m=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2303           d=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2304           y=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2305           d_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2306           y_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2307           M_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2308           A_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2309           B_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2310           C_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2311           D_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2312           X_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
2313           Y_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2314           m_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2315           d_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2316           y_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2317           d_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2318           y_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2319           r=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.RIGHTHANDSIDE),
2320           q=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.BOTH) )
2321    
2322       def __str__(self):
2323         """
2324         returns string representation of the transport problem
2325    
2326         @return: a simple representation of the transport problem
2327         @rtype: C{str}
2328         """
2329         return "<TransportPDE %d>"%id(self)
2330    
2331       def getTheta(self):
2332          """
2333          returns the time stepping control parameter
2334          @rtype: float
2335        """        """
2336        creates a linear, steady, second order PDE on a L{Domain<escript.Domain>}        return self.__theta
2337       def getRequiredSystemType(self):
2338          """
2339          returns the system type which needs to be used by the current set up.
2340          """
2341          return 0
2342          return self.getDomain().getTransportMatrixTypeId(self.getSolverMethod()[0],self.getSolverMethod()[1],self.getSolverPackage(),self.isSymmetric())
2343    
2344        @param domain: domain of the PDE     def checkSymmetry(self,verbose=True):
       @type domain: L{Domain<escript.Domain>}  
       @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.  
2345        """        """
2346        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()  
2347    
2348     def setValue(**coefficients):        @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.
2349          @type verbose: C{bool}
2350          @return:  True if the PDE is symmetric.
2351          @rtype: C{bool}
2352          @note: This is a very expensive operation. It should be used for degugging only! The symmetry flag is not altered.
2353          """
2354          out=True
2355          out=out and self.checkSymmetricTensor("M", verbose)
2356          out=out and self.checkSymmetricTensor("M_reduced", verbose)
2357          out=out and self.checkSymmetricTensor("A", verbose)
2358          out=out and self.checkSymmetricTensor("A_reduced", verbose)
2359          out=out and self.checkReciprocalSymmetry("B","C", verbose)
2360          out=out and self.checkReciprocalSymmetry("B_reduced","C_reduced", verbose)
2361          out=out and self.checkSymmetricTensor("D", verbose)
2362          out=out and self.checkSymmetricTensor("D_reduced", verbose)
2363          out=out and self.checkSymmetricTensor("m", verbose)
2364          out=out and self.checkSymmetricTensor("m_reduced", verbose)
2365          out=out and self.checkSymmetricTensor("d", verbose)
2366          out=out and self.checkSymmetricTensor("d_reduced", verbose)
2367          out=out and self.checkSymmetricTensor("d_contact", verbose)
2368          out=out and self.checkSymmetricTensor("d_contact_reduced", verbose)
2369          return out
2370    
2371       def setValue(self,**coefficients):
2372        """        """
2373        sets new values to coefficients        sets new values to coefficients
2374    
2375        @param coefficients: new values assigned to coefficients        @param coefficients: new values assigned to coefficients
2376          @keyword M: value for coefficient M.
2377          @type M: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
2378          @keyword M_reduced: value for coefficient M_reduced.
2379          @type M_reduced: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.ReducedFunction>}.
2380        @keyword A: value for coefficient A.        @keyword A: value for coefficient A.
2381        @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>}.
2382          @keyword A_reduced: value for coefficient A_reduced.
2383          @type A_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
2384        @keyword B: value for coefficient B        @keyword B: value for coefficient B
2385        @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>}.
2386          @keyword B_reduced: value for coefficient B_reduced
2387          @type B_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
2388        @keyword C: value for coefficient C        @keyword C: value for coefficient C
2389        @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>}.
2390          @keyword C_reduced: value for coefficient C_reduced
2391          @type C_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
2392        @keyword D: value for coefficient D        @keyword D: value for coefficient D
2393        @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>}.
2394          @keyword D_reduced: value for coefficient D_reduced
2395          @type D_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
2396        @keyword X: value for coefficient X        @keyword X: value for coefficient X
2397        @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>}.
2398          @keyword X_reduced: value for coefficient X_reduced
2399          @type X_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
2400        @keyword Y: value for coefficient Y        @keyword Y: value for coefficient Y
2401        @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>}.
2402          @keyword Y_reduced: value for coefficient Y_reduced
2403          @type Y_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.Function>}.
2404          @keyword m: value for coefficient m
2405          @type m: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
2406          @keyword m_reduced: value for coefficient m_reduced
2407          @type m_reduced: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
2408        @keyword d: value for coefficient d        @keyword d: value for coefficient d
2409        @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>}.
2410          @keyword d_reduced: value for coefficient d_reduced
2411          @type d_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
2412        @keyword y: value for coefficient y        @keyword y: value for coefficient y
2413        @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>}.
2414        @keyword d_contact: value for coefficient d_contact        @keyword d_contact: value for coefficient d_contact
2415        @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>}.
2416                         or  L{FunctionOnContactZero<escript.FunctionOnContactZero>}.        @keyword d_contact_reduced: value for coefficient d_contact_reduced
2417          @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>}.
2418        @keyword y_contact: value for coefficient y_contact        @keyword y_contact: value for coefficient y_contact
2419        @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>}.
2420                         or  L{FunctionOnContactZero<escript.FunctionOnContactZero>}.        @keyword y_contact_reduced: value for coefficient y_contact_reduced
2421          @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>}.
2422        @keyword r: values prescribed to the solution at the locations of constraints        @keyword r: values prescribed to the solution at the locations of constraints
2423        @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>}
2424                 depending of reduced order is used for the solution.                 depending of reduced order is used for the solution.
# Line 1978  class AdvectivePDE(LinearPDE): Line 2426  class AdvectivePDE(LinearPDE):
2426        @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>}
2427                 depending of reduced order is used for the representation of the equation.                 depending of reduced order is used for the representation of the equation.
2428        @raise IllegalCoefficient: if an unknown coefficient keyword is used.        @raise IllegalCoefficient: if an unknown coefficient keyword is used.
   
2429        """        """
2430        if "A" in coefficients.keys()   or "B" in coefficients.keys() or "C" in coefficients.keys(): self.__Xi=escript.Data()        super(TransportPDE,self).setValue(**coefficients)
       super(AdvectivePDE, self).setValue(**coefficients)  
   
    def ELMAN_RAMAGE(self,P):  
      """  
      Predefined function to set a values for M{S{xi}} from a Preclet number M{P}.  
      This function uses the method suggested by H.C. Elman and A. Ramage, I{SIAM J. Numer. Anal.}, B{40} (2002)  
           - 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 __calculateXi(self,peclet_factor,flux,h):  
        flux=util.Lsup(flux)  
        if flux_max>0.:  
           return h*self.__xi(flux*peclet_factor)/(flux+flux_max*self.__TOL)  
        else:  
           return 0.  
   
    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():  
                 flux2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))  
                 if self.getNumEquations()>1:  
                    if self.getNumSolutions()>1:  
                       for i in range(self.getNumEquations()):  
                          for k in range(self.getNumSolutions()):  
                             for l in range(self.getDim()): flux2+=(C[i,k,l]-B[i,l,k])**2  
                       # flux=C-util.reorderComponents(B,[0,2,1])  
                    else:  
                       for i in range(self.getNumEquations()):  
                          for l in range(self.getDim()): flux2+=(C[i,l]-B[i,l])**2  
                       # flux=C-B  
                 else:  
                    if self.getNumSolutions()>1:  
                       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])  
                    else:  
                       for l in range(self.getDim()): flux2+=(C[l]-B[l])**2  
                       #flux=C-B  
                 length_of_flux=util.sqrt(flux2)  
             elif C.isEmpty():  
               length_of_flux=util.length(B)  
               #flux=B  
             else:  
               length_of_flux=util.length(C)  
               #flux=C  
   
             #length_of_flux=util.length(flux)  
             flux_max=util.Lsup(length_of_flux)  
             if flux_max>0.:  
                # 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.__TOL)  
                else:  
                     inv_A=1./self.__TOL  
                peclet_number=length_of_flux*h/2*inv_A  
                xi=self.__xi(peclet_number)  
                self.__Xi=h*xi/(length_of_flux+flux_max*self.__TOL)  
                self.trace("preclet number = %e"%util.Lsup(peclet_number))  
       return self.__Xi  
   
   
    def getCoefficientOfGeneralPDE(self,name):  
      """  
      return the value of the coefficient name of 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}.  
      @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."  
   
      if name == "A" :  
          A=self.getCoefficient("A")  
          B=self.getCoefficient("B")  
          C=self.getCoefficient("C")  
          if B.isEmpty() and C.isEmpty():  
             Aout=A  
          else:  
             if A.isEmpty():  
                Aout=self.createNewCoefficient("A")  
             else:  
                Aout=A[:]  
             Xi=self.__getXi()  
             if self.getNumEquations()>1:  
                 for i in range(self.getNumEquations()):  
                    for j in range(self.getDim()):  
                       for k in range(self.getNumSolutions()):  
                          for l in range(self.getDim()):  
                             if not C.isEmpty() and not B.isEmpty():  
                                # 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:  
                 for j in range(self.getDim()):  
                    for l in range(self.getDim()):  
                       if not C.isEmpty() and not B.isEmpty():  
                           Aout[j,l]+=Xi*(C[j]-B[j])*(C[l]-B[l])  
                       elif C.isEmpty():  
                           Aout[j,l]+=Xi*B[j]*B[l]  
                       else:  
                           Aout[j,l]+=Xi*C[j]*C[l]  
                  # if not C.isEmpty() and not B.isEmpty():  
                  #    tmp=C-B  
                  #    Aout=Aout+Xi*util.outer(tmp,tmp)  
                  # elif C.isEmpty():  
                  #    Aout=Aout+Xi*util.outer(B,B)  
                  # else:  
                  # Aout=Aout+Xi*util.outer(C,C)  
          return Aout  
      elif name == "B" :  
          B=self.getCoefficient("B")  
          C=self.getCoefficient("C")  
          D=self.getCoefficient("D")  
          if C.isEmpty() or D.isEmpty():  
             Bout=B  
          else:  
             Xi=self.__getXi()  
             if B.isEmpty():  
                 Bout=self.createNewCoefficient("B")  
             else:  
                 Bout=B[:]  
             if self.getNumEquations()>1:  
                for k in range(self.getNumSolutions()):  
                   for p in range(self.getNumEquations()):  
                      tmp=Xi*D[p,k]  
                      for i in range(self.getNumEquations()):  
                         for j in range(self.getDim()):  
                            Bout[i,j,k]+=tmp*C[p,i,j]  
                            # Bout=Bout+Xi*util.generalTensorProduct(util.reorder(C,[1,2,0]),D,offset=1)  
             else:  
                tmp=Xi*D  
                for j in range(self.getDim()): Bout[j]+=tmp*C[j]  
                # Bout=Bout+Xi*D*C  
          return Bout  
      elif name == "C" :  
          B=self.getCoefficient("B")  
          C=self.getCoefficient("C")  
          D=self.getCoefficient("D")  
          if B.isEmpty() or D.isEmpty():  
             Cout=C  
          else:  
             Xi=self.__getXi()  
             if C.isEmpty():  
                 Cout=self.createNewCoefficient("C")  
             else:  
                 Cout=C[:]  
             if self.getNumEquations()>1:  
                for k in range(self.getNumSolutions()):  
                    for p in range(self.getNumEquations()):  
                       tmp=Xi*D[p,k]  
                       for i in range(self.getNumEquations()):  
                         for l in range(self.getDim()):  
                                  Cout[i,k,l]+=tmp*B[p,l,i]  
                                  # Cout=Cout+Xi*B[p,l,i]*D[p,k]  
             else:  
                tmp=Xi*D  
                for j in range(self.getDim()): Cout[j]+=tmp*B[j]  
                # Cout=Cout+tmp*D*B  
          return Cout  
      elif name == "D" :  
          return self.getCoefficient("D")  
      elif name == "X" :  
          X=self.getCoefficient("X")  
          Y=self.getCoefficient("Y")  
          B=self.getCoefficient("B")  
          C=self.getCoefficient("C")  
          if Y.isEmpty() or (B.isEmpty() and C.isEmpty()):  
             Xout=X  
          else:  
             if X.isEmpty():  
                 Xout=self.createNewCoefficient("X")  
             else:  
                 Xout=X[:]  
             Xi=self.__getXi()  
             if self.getNumEquations()>1:  
                  for p in range(self.getNumEquations()):  
                     tmp=Xi*Y[p]  
                     for i in range(self.getNumEquations()):  
                        for j in range(self.getDim()):  
                           if not C.isEmpty() and not B.isEmpty():  
                              Xout[i,j]+=tmp*(C[p,i,j]-B[p,j,i])  
                              # 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:  
                  tmp=Xi*Y  
                  for j in range(self.getDim()):  
                     if not C.isEmpty() and not B.isEmpty():  
                        Xout[j]+=tmp*(C[j]-B[j])  
                        # Xout=Xout+Xi*Y*(C-B)  
                     elif C.isEmpty():  
                        Xout[j]-=tmp*B[j]  
                        # Xout=Xout-Xi*Y*B  
                     else:  
                        Xout[j]+=tmp*C[j]  
                        # Xout=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  
   
 class AdvectionDiffusion(LinearPDE):  
    """  
    Class to define PDE equation of the unisotropic advection-diffusion problem, which is genear L{LinearPDE} of the form  
   
    M{S{omega}*u + inner(v,grad(u))- grad(matrixmult(k_bar,grad(u))[j])[j] = f}  
   
    with natural boundary conditons  
   
    M{n[j]*matrixmult(k,grad(u))[j] = g- S{alpha}u }  
   
    and constraints:  
   
    M{u=r} where M{q>0}  
   
    and  
   
    M{k_bar[i,j]=k[i,j]+upwind[i]*upwind[j]}  
   
    """  
   
    def __init__(self,domain,debug=False):  
      """  
      initializes a new Poisson equation  
   
      @param domain: domain of the PDE  
      @type domain: L{Domain<escript.Domain>}  
      @param debug: if True debug informations are printed.  
   
      """  
      super(AdvectionDiffusion, self).__init__(domain,1,1,debug)  
      self.COEFFICIENTS={"omega": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),  
                         "k": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_DIM,PDECoefficient.BY_DIM),PDECoefficient.OPERATOR),  
                         "f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),  
                         "v": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_DIM,),PDECoefficient.OPERATOR),  
                         "upwind": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_DIM,),PDECoefficient.OPERATOR),  
                         "alpha": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),  
                         "g": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),  
                         "r": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH),  
                         "q": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}  
2431    
2432     def setValue(self,**coefficients):     def createOperator(self):
2433       """         """
2434       sets new values to coefficients         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       @param coefficients: new values assigned to coefficients     def getSolution(self,dt,**options):
2444       @keyword omega: value for coefficient M{S{omega}}         """
2445       @type omega: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.         returns the solution of the problem
2446       @keyword k: value for coefficient M{k}         @return: the solution
2447       @type k: any type that can be casted to L{Tensor<escript.Tensor>} object on L{Function<escript.Function>}.         @rtype: L{Data<escript.Data>}
      @keyword v: value for coefficient M{v}  
      @type v: any type that can be casted to L{Vector<escript.Vector>} object on L{Function<escript.Function>}.  
      @keyword upwind: value for upwind term M{upwind}  
      @type upwind: any type that can be casted to L{Vector<escript.Vector>} object on L{Function<escript.Function>}.  
      @keyword f: value for right hand side M{f}  
      @type f: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.  
      @keyword alpha: value for right hand side M{S{alpha}}  
      @type alpha: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.  
      @keyword g: value for right hand side M{g}  
      @type g: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.  
      @keyword r: prescribed values M{r} for the solution in constraints.  
      @type r: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}  
                depending of reduced order is used for the representation of the equation.  
      @keyword q: mask for location of constraints  
      @type q: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}  
                depending of reduced order is used for the representation of the equation.  
      @raise IllegalCoefficient: if an unknown coefficient keyword is used.  
      """  
      super(AdvectionDiffusion, self).setValue(**coefficients)  
2448    
2449     def getCoefficientOfGeneralPDE(self,name):         @note: This method is overwritten when implementing a particular linear problem.
2450       """         """
2451       return the value of the coefficient name of the general PDE         if not self.isSolutionValid():
2452               options[self.TOLERANCE_KEY]=self.getTolerance()
2453               self.setSolution(self.getOperator().solve(self.getRightHandSide(),dt,options))
2454               self.validSolution()
2455           return self.getCurrentSolution()
2456       def setInitialSolution(self,u):
2457           self.getOperator().setInitialValue(util.interpolate(u,self.getFunctionSpaceForSolution()))
2458       def getSafeTimeStepSize(self):
2459           return self.getOperator().getSafeTimeStepSize()
2460    
2461       @param name: name of the coefficient requested.     def getSystem(self):
2462       @type name: C{string}         """
2463       @return: the value of the coefficient  name         return the operator and right hand side of the PDE
      @rtype: L{Data<escript.Data>}  
      @raise IllegalCoefficient: if name is not one of coefficients  
                   "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 name == "A" :  
          return self.getCoefficient("k")+util.outer(self.getCoefficient("upwind"),self.getCoefficient("upwind"))  
      elif name == "B" :  
          return escript.Data()  
      elif name == "C" :  
          return self.getCoefficient("v")  
      elif name == "D" :  
          return self.getCoefficient("omega")  
      elif name == "X" :  
          return escript.Data()  
      elif name == "Y" :  
          return self.getCoefficient("f")  
      elif name == "d" :  
          return self.getCoefficient("alpha")  
      elif name == "y" :  
          return self.getCoefficient("g")  
      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")  
      else:  
         raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name  
2464    
2465           @return: the discrete version of the PDE
2466           @rtype: C{tuple} of L{Operator,<escript.Operator>} and L{Data<escript.Data>}.
2467    
2468  # $Log$         @note: This method is overwritten when implementing a particular linear problem.
2469  # Revision 1.14  2005/09/22 01:54:57  jgs         """
2470  # Merge of development branch dev-02 back to main trunk on 2005-09-22         if not self.isOperatorValid() or not self.isRightHandSideValid():
2471  #            self.resetRightHandSide()
2472  # Revision 1.13  2005/09/15 03:44:19  jgs            righthandside=self.getCurrentRightHandSide()
2473  # Merge of development branch dev-02 back to main trunk on 2005-09-15            self.resetOperator()
2474  #            operator=self.getCurrentOperator()
2475  # Revision 1.12  2005/09/01 03:31:28  jgs            self.getDomain().addPDEToTransportProblem(
2476  # Merge of development branch dev-02 back to main trunk on 2005-09-01                              operator,
2477  #                              righthandside,
2478  # Revision 1.11  2005/08/23 01:24:28  jgs                              self.getCoefficient("M"),
2479  # Merge of development branch dev-02 back to main trunk on 2005-08-23                              self.getCoefficient("A"),
2480  #                              self.getCoefficient("B"),
2481  # Revision 1.10  2005/08/12 01:45:36  jgs                              self.getCoefficient("C"),
2482  # erge of development branch dev-02 back to main trunk on 2005-08-12                              self.getCoefficient("D"),
2483  #                              self.getCoefficient("X"),
2484  # Revision 1.9.2.17  2005/09/21 07:03:33  matt                              self.getCoefficient("Y"),
2485  # PDECoefficient and LinearPDE are now new style classes (introduced in Python                              self.getCoefficient("d"),
2486  # 2.2). Classes Poisson, Helmholtz, LameEquation and AdvectivePDE have been                              self.getCoefficient("y"),
2487  # modified to instead use portable/cooperative "super" calls to extend base                              self.getCoefficient("d_contact"),
2488  # class methods.                              self.getCoefficient("y_contact"))
2489  #            operator.insertConstraint(righthandside,self.getCoefficient("q"),self.getCoefficient("r"))
2490  # Revision 1.9.2.16  2005/09/16 01:54:37  matt            self.trace("New system has been built.")
2491  # Removed redundant if-loop.            self.validOperator()
2492  #            self.validRightHandSide()
2493  # Revision 1.9.2.15  2005/09/14 08:09:18  matt         return (self.getCurrentOperator(), self.getCurrentRightHandSide())
 # Added "REDUCED" solution PDECoefficient descriptors for LinearPDEs.  
 #  
 # Revision 1.9.2.14  2005/09/07 06:26:16  gross  
 # the solver from finley are put into the standalone package paso now  
 #  
 # Revision 1.9.2.13  2005/08/31 08:45:03  gross  
 # in the case of lumping no new system is allocated if the constraint is changed.  
 #  
 # Revision 1.9.2.12  2005/08/31 07:10:23  gross  
 # test for Lumping added  
 #  
 # Revision 1.9.2.11  2005/08/30 01:53:45  gross  
 # bug in format fixed.  
 #  
 # Revision 1.9.2.10  2005/08/26 07:14:17  gross  
 # a few more bugs in linearPDE fixed. remaining problem are finley problems  
 #  
 # Revision 1.9.2.9  2005/08/26 06:30:45  gross  
 # fix for reported bug  0000004. test_linearPDE passes a few more tests  
 #  
 # Revision 1.9.2.8  2005/08/26 04:30:13  gross  
 # gneric unit testing for linearPDE  
 #  
 # Revision 1.9.2.7  2005/08/25 07:06:50  gross  
 # linearPDE documentation is parsed now by epydoc. there is still a problem with links into escriptcpp.so  
 #  
 # Revision 1.9.2.6  2005/08/24 05:01:24  gross  
 # problem with resetting the matrix in case of resetting its values to 0 fixed.  
 #  
 # Revision 1.9.2.5  2005/08/24 02:03:28  gross  
 # 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.425  
changed lines
  Added in v.1841

  ViewVC Help
Powered by ViewVC 1.1.26