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

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

  ViewVC Help
Powered by ViewVC 1.1.26