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

Legend:
Removed from v.1400  
changed lines
  Added in v.2061

  ViewVC Help
Powered by ViewVC 1.1.26