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

  ViewVC Help
Powered by ViewVC 1.1.26