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

Legend:
Removed from v.720  
changed lines
  Added in v.2100

  ViewVC Help
Powered by ViewVC 1.1.26