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

Legend:
Removed from v.1819  
changed lines
  Added in v.1841

  ViewVC Help
Powered by ViewVC 1.1.26