/[escript]/trunk-mpi-branch/escript/py_src/linearPDEs.py
ViewVC logotype

Diff of /trunk-mpi-branch/escript/py_src/linearPDEs.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/esys2/escript/py_src/linearPDEs.py revision 104 by jgs, Fri Dec 17 07:43:12 2004 UTC trunk-mpi-branch/escript/py_src/linearPDEs.py revision 1302 by ksteube, Mon Sep 17 06:57:51 2007 UTC
# Line 1  Line 1 
1  # $Id$  # $Id$
   
 ## @file linearPDEs.py  
   
2  """  """
3  @brief Functions and classes for linear PDEs  The module provides an interface to define and solve linear partial
4    differential equations (PDEs) within L{escript}. L{linearPDEs} does not provide any
5    solver capabilities in itself but hands the PDE over to
6    the PDE solver library defined through the L{Domain<escript.Domain>} of the PDE.
7    The general interface is provided through the L{LinearPDE} class. The
8    L{AdvectivePDE} which is derived from the L{LinearPDE} class
9    provides an interface to PDE dominated by its advective terms. The L{Poisson},
10    L{Helmholtz}, L{LameEquation}, L{AdvectivePDE}
11    classs which are also derived form the L{LinearPDE} class should be used
12    to define of solve these sepecial PDEs.
13    
14    @var __author__: name of author
15    @var __copyright__: copyrights
16    @var __license__: licence agreement
17    @var __url__: url entry point on documentation
18    @var __version__: version
19    @var __date__: date of the version
20  """  """
21    
22  import escript  import escript
23  import util  import util
24  import numarray  import numarray
25    
26  def identifyDomain(domain=None,data={}):  __author__="Lutz Gross, l.gross@uq.edu.au"
27       """  __copyright__="""  Copyright (c) 2006 by ACcESS MNRF
28       @brief Return the Domain which is equal to the input domain (if not None)                      http://www.access.edu.au
29       and is the domain of all Data objects in the dictionary data.                  Primary Business: Queensland, Australia"""
30       An exception is raised if this is not possible  __license__="""Licensed under the Open Software License version 3.0
31                 http://www.opensource.org/licenses/osl-3.0.php"""
32       @param domain  __url__="http://www.iservo.edu.au/esys"
33       @param data  __version__="$Revision$"
34       """  __date__="$Date$"
35       # get the domain used by any Data object in the list data:  
36       data_domain=None  
37       for d in data.itervalues():  class IllegalCoefficient(ValueError):
38            if isinstance(d,escript.Data):     """
39               if not d.isEmpty(): data_domain=d.getDomain()     raised if an illegal coefficient of the general ar particular PDE is requested.
40       # check if domain and data_domain are identical?     """
41       if domain == None:     pass
42           if data_domain == None:  
43                raise ValueError,"Undefined PDE domain. Specify a domain or use a Data class object as coefficient"  class IllegalCoefficientValue(ValueError):
44       else:     """
45           if data_domain == None:     raised if an incorrect value for a coefficient is used.
46                data_domain=domain     """
47           else:     pass
48             if not data_domain == domain:  
49                   raise ValueError,"Domain of coefficients doesnot match specified domain"  class IllegalCoefficientFunctionSpace(ValueError):
50       # now we check if all Data class object coefficients are defined on data_domain:     """
51       for i,d in data.iteritems():     raised if an incorrect function space for a coefficient is used.
52           if isinstance(d,escript.Data):     """
             if not d.isEmpty():  
                if not data_domain==d.getDomain():  
                  raise ValueError,"Illegal domain for coefficient %s."%i  
      # done:  
      return data_domain  
   
 def identifyNumEquationsAndSolutions(dim,coef={}):  
      # get number of equations and number of unknowns:  
      numEquations=0  
      numSolutions=0  
      for i in coef.iterkeys():  
         if not coef[i].isEmpty():  
            res=_PDECoefficientTypes[i].estimateNumEquationsAndNumSolutions(coef[i].getShape(),dim)  
            if res==None:  
                raise ValueError,"Illegal shape %s of coefficient %s"%(coef[i].getShape().__str__(),i)  
            else:  
                numEquations=max(numEquations,res[0])  
                numSolutions=max(numSolutions,res[1])  
      return numEquations,numSolutions  
   
   
 def _CompTuple2(t1,t2):  
    """  
    @brief  
   
    @param t1  
    @param t2  
    """  
    dif=t1[0]+t1[1]-(t2[0]+t2[1])  
    if dif<0: return 1  
    elif dif>0: return -1  
    else: return 0  
53    
54  class PDECoefficientType:  class UndefinedPDEError(ValueError):
55       """
56       raised if a PDE is not fully defined yet.
57       """
58       pass
59    
60    class PDECoefficient(object):
61      """      """
62      @brief      A class for describing a PDE coefficient
63    
64        @cvar INTERIOR: indicator that coefficient is defined on the interior of the PDE domain
65        @cvar BOUNDARY: indicator that coefficient is defined on the boundary of the PDE domain
66        @cvar CONTACT: indicator that coefficient is defined on the contact region within the PDE domain
67        @cvar INTERIOR_REDUCED: indicator that coefficient is defined on the interior of the PDE domain using a reduced integration order
68        @cvar BOUNDARY_REDUCED: indicator that coefficient is defined on the boundary of the PDE domain using a reduced integration order
69        @cvar CONTACT_REDUCED: indicator that coefficient is defined on the contact region within the PDE domain using a reduced integration order
70        @cvar SOLUTION: indicator that coefficient is defined trough a solution of the PDE
71        @cvar REDUCED: indicator that coefficient is defined trough a reduced solution of the PDE
72        @cvar BY_EQUATION: indicator that the dimension of the coefficient shape is defined by the number PDE equations
73        @cvar BY_SOLUTION: indicator that the dimension of the coefficient shape is defined by the number PDE solutions
74        @cvar BY_DIM: indicator that the dimension of the coefficient shape is defined by the spatial dimension
75        @cvar OPERATOR: indicator that the the coefficient alters the operator of the PDE
76        @cvar RIGHTHANDSIDE: indicator that the the coefficient alters the right hand side of the PDE
77        @cvar BOTH: indicator that the the coefficient alters the operator as well as the right hand side of the PDE
78    
79      """      """
     # identifier for location of Data objects defining coefficients  
80      INTERIOR=0      INTERIOR=0
81      BOUNDARY=1      BOUNDARY=1
82      CONTACT=2      CONTACT=2
83      CONTINUOUS=3      SOLUTION=3
84      # identifier in the pattern of coefficients:      REDUCED=4
85      # the pattern is a tuple of EQUATION,SOLUTION,DIM where DIM represents the spatial dimension, EQUATION the number of equations and SOLUTION the      BY_EQUATION=5
86      # number of unknowns.      BY_SOLUTION=6
87      EQUATION=3      BY_DIM=7
88      SOLUTION=4      OPERATOR=10
89      DIM=5      RIGHTHANDSIDE=11
90      # indicator for what is altered if the coefficient is altered:      BOTH=12
91      OPERATOR=5      INTERIOR_REDUCED=13
92      RIGHTHANDSIDE=6      BOUNDARY_REDUCED=14
93      BOTH=7      CONTACT_REDUCED=15
94      def __init__(self,where,pattern,altering):  
95         """      def __init__(self, where, pattern, altering):
96         @brief Initialise a PDE Coefficient type         """
97           Initialise a PDE Coefficient type
98    
99           @param where: describes where the coefficient lives
100           @type where: one of L{INTERIOR}, L{BOUNDARY}, L{CONTACT}, L{SOLUTION}, L{REDUCED},
101                               L{INTERIOR_REDUCED}, L{BOUNDARY_REDUCED}, L{CONTACT_REDUCED}.
102           @param pattern: describes the shape of the coefficient and how the shape is build for a given
103                  spatial dimension and numbers of equation and solution in then PDE. For instance,
104                  (L{BY_EQUATION},L{BY_SOLUTION},L{BY_DIM}) descrbes a rank 3 coefficient which
105                  is instanciated as shape (3,2,2) in case of a three equations and two solution components
106                  on a 2-dimensional domain. In the case of single equation and a single solution component
107                  the shape compoments marked by L{BY_EQUATION} or L{BY_SOLUTION} are dropped. In this case
108                  the example would be read as (2,).
109           @type pattern: C{tuple} of L{BY_EQUATION}, L{BY_SOLUTION}, L{BY_DIM}
110           @param altering: indicates what part of the PDE is altered if the coefficiennt is altered
111           @type altering: one of L{OPERATOR}, L{RIGHTHANDSIDE}, L{BOTH}
112           @param reduced: indicates if reduced
113           @type reduced: C{bool}
114         """         """
115           super(PDECoefficient, self).__init__()
116         self.what=where         self.what=where
117         self.pattern=pattern         self.pattern=pattern
118         self.altering=altering         self.altering=altering
119           self.resetValue()
120    
121      def getFunctionSpace(self,domain):      def resetValue(self):
122         """         """
123         @brief defines the FunctionSpace of the coefficient on the domain         resets coefficient value to default
   
        @param domain  
124         """         """
125         if self.what==self.INTERIOR: return escript.Function(domain)         self.value=escript.Data()
126         elif self.what==self.BOUNDARY: return escript.FunctionOnBoundary(domain)  
127         elif self.what==self.CONTACT: return escript.FunctionOnContactZero(domain)      def getFunctionSpace(self,domain,reducedEquationOrder=False,reducedSolutionOrder=False):
128         elif self.what==self.CONTINUOUS: return escript.ContinuousFunction(domain)         """
129           defines the L{FunctionSpace<escript.FunctionSpace>} of the coefficient
130    
131           @param domain: domain on which the PDE uses the coefficient
132           @type domain: L{Domain<escript.Domain>}
133           @param reducedEquationOrder: True to indicate that reduced order is used to represent the equation
134           @type reducedEquationOrder: C{bool}
135           @param reducedSolutionOrder: True to indicate that reduced order is used to represent the solution
136           @type reducedSolutionOrder: C{bool}
137           @return:  L{FunctionSpace<escript.FunctionSpace>} of the coefficient
138           @rtype:  L{FunctionSpace<escript.FunctionSpace>}
139           """
140           if self.what==self.INTERIOR:
141                return escript.Function(domain)
142           elif self.what==self.INTERIOR_REDUCED:
143                return escript.ReducedFunction(domain)
144           elif self.what==self.BOUNDARY:
145                return escript.FunctionOnBoundary(domain)
146           elif self.what==self.BOUNDARY_REDUCED:
147                return escript.ReducedFunctionOnBoundary(domain)
148           elif self.what==self.CONTACT:
149                return escript.FunctionOnContactZero(domain)
150           elif self.what==self.CONTACT_REDUCED:
151                return escript.ReducedFunctionOnContactZero(domain)
152           elif self.what==self.SOLUTION:
153                if reducedEquationOrder and reducedSolutionOrder:
154                    return escript.ReducedSolution(domain)
155                else:
156                    return escript.Solution(domain)
157           elif self.what==self.REDUCED:
158                return escript.ReducedSolution(domain)
159    
160        def getValue(self):
161           """
162           returns the value of the coefficient
163    
164           @return:  value of the coefficient
165           @rtype:  L{Data<escript.Data>}
166           """
167           return self.value
168    
169        def setValue(self,domain,numEquations=1,numSolutions=1,reducedEquationOrder=False,reducedSolutionOrder=False,newValue=None):
170           """
171           set the value of the coefficient to a new value
172    
173           @param domain: domain on which the PDE uses the coefficient
174           @type domain: L{Domain<escript.Domain>}
175           @param numEquations: number of equations of the PDE
176           @type numEquations: C{int}
177           @param numSolutions: number of components of the PDE solution
178           @type numSolutions: C{int}
179           @param reducedEquationOrder: True to indicate that reduced order is used to represent the equation
180           @type reducedEquationOrder: C{bool}
181           @param reducedSolutionOrder: True to indicate that reduced order is used to represent the solution
182           @type reducedSolutionOrder: C{bool}
183           @param newValue: number of components of the PDE solution
184           @type newValue: any object that can be converted into a L{Data<escript.Data>} object with the appropriate shape and L{FunctionSpace<escript.FunctionSpace>}
185           @raise IllegalCoefficientValue: if the shape of the assigned value does not match the shape of the coefficient
186           @raise IllegalCoefficientFunctionSpace: if unable to interploate value to appropriate function space
187           """
188           if newValue==None:
189               newValue=escript.Data()
190           elif isinstance(newValue,escript.Data):
191               if not newValue.isEmpty():
192                  if not newValue.getFunctionSpace() == self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder):
193                    try:
194                      newValue=escript.Data(newValue,self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder))
195                    except:
196                      raise IllegalCoefficientFunctionSpace,"Unable to interpolate coefficient to function space %s"%self.getFunctionSpace(domain)
197           else:
198               newValue=escript.Data(newValue,self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder))
199           if not newValue.isEmpty():
200               if not self.getShape(domain,numEquations,numSolutions)==newValue.getShape():
201                   raise IllegalCoefficientValue,"Expected shape of coefficient is %s but actual shape is %s."%(self.getShape(domain,numEquations,numSolutions),newValue.getShape())
202           self.value=newValue
203    
204      def isAlteringOperator(self):      def isAlteringOperator(self):
205          """          """
206      @brief return true if the operator of the PDE is changed when the coefficient is changed          checks if the coefficient alters the operator of the PDE
207    
208            @return:  True if the operator of the PDE is changed when the coefficient is changed
209            @rtype:  C{bool}
210      """      """
211          if self.altering==self.OPERATOR or self.altering==self.BOTH:          if self.altering==self.OPERATOR or self.altering==self.BOTH:
212              return not None              return not None
# Line 119  class PDECoefficientType: Line 215  class PDECoefficientType:
215    
216      def isAlteringRightHandSide(self):      def isAlteringRightHandSide(self):
217          """          """
218      @brief return true if the right hand side of the PDE is changed when the coefficient is changed          checks if the coefficeint alters the right hand side of the PDE
219    
220        @rtype:  C{bool}
221            @return:  True if the right hand side of the PDE is changed when the coefficient is changed
222      """      """
223          if self.altering==self.RIGHTHANDSIDE or self.altering==self.BOTH:          if self.altering==self.RIGHTHANDSIDE or self.altering==self.BOTH:
224              return not None              return not None
225          else:          else:
226              return None              return None
227    
228      def estimateNumEquationsAndNumSolutions(self,shape=(),dim=3):      def estimateNumEquationsAndNumSolutions(self,domain,shape=()):
229         """         """
230         @brief tries to estimate the number of equations in a given tensor shape for a given spatial dimension dim         tries to estimate the number of equations and number of solutions if the coefficient has the given shape
231    
232         @param shape         @param domain: domain on which the PDE uses the coefficient
233         @param dim         @type domain: L{Domain<escript.Domain>}
234           @param shape: suggested shape of the coefficient
235           @type shape: C{tuple} of C{int} values
236           @return: the number of equations and number of solutions of the PDE is the coefficient has shape s.
237                     If no appropriate numbers could be identified, C{None} is returned
238           @rtype: C{tuple} of two C{int} values or C{None}
239         """         """
240           dim=domain.getDim()
241         if len(shape)>0:         if len(shape)>0:
242             num=max(shape)+1             num=max(shape)+1
243         else:         else:
244             num=1             num=1
245         search=[]         search=[]
246         for u in range(num):         if self.definesNumEquation() and self.definesNumSolutions():
247            for e in range(num):            for u in range(num):
248               search.append((e,u))               for e in range(num):
249         search.sort(_CompTuple2)                  search.append((e,u))
250         for item in search:            search.sort(self.__CompTuple2)
251               s=self.buildShape(item[0],item[1],dim)            for item in search:
252                 s=self.getShape(domain,item[0],item[1])
253               if len(s)==0 and len(shape)==0:               if len(s)==0 and len(shape)==0:
254                   return (1,1)                   return (1,1)
255               else:               else:
256                   if s==shape: return item                   if s==shape: return item
257           elif self.definesNumEquation():
258              for e in range(num,0,-1):
259                 s=self.getShape(domain,e,0)
260                 if len(s)==0 and len(shape)==0:
261                     return (1,None)
262                 else:
263                     if s==shape: return (e,None)
264    
265           elif self.definesNumSolutions():
266              for u in range(num,0,-1):
267                 s=self.getShape(domain,0,u)
268                 if len(s)==0 and len(shape)==0:
269                     return (None,1)
270                 else:
271                     if s==shape: return (None,u)
272         return None         return None
273        def definesNumSolutions(self):
274           """
275           checks if the coefficient allows to estimate the number of solution components
276    
277      def buildShape(self,e=1,u=1,dim=3):         @return: True if the coefficient allows an estimate of the number of solution components
278          """         @rtype: C{bool}
279      @brief builds the required shape for a given number of equations e, number of unknowns u and spatial dimension dim         """
280           for i in self.pattern:
281                 if i==self.BY_SOLUTION: return True
282           return False
283    
284      @param e      def definesNumEquation(self):
285      @param u         """
286      @param dim         checks if the coefficient allows to estimate the number of equations
287      """  
288          s=()         @return: True if the coefficient allows an estimate of the number of equations
289          for i in self.pattern:         @rtype: C{bool}
290               if i==self.EQUATION:         """
291                  if e>1: s=s+(e,)         for i in self.pattern:
292               elif i==self.SOLUTION:               if i==self.BY_EQUATION: return True
293                  if u>1: s=s+(u,)         return False
294    
295        def __CompTuple2(self,t1,t2):
296          """
297          Compare two tuples of possible number of equations and number of solutions
298    
299          @param t1: The first tuple
300          @param t2: The second tuple
301    
302          """
303    
304          dif=t1[0]+t1[1]-(t2[0]+t2[1])
305          if dif<0: return 1
306          elif dif>0: return -1
307          else: return 0
308    
309        def getShape(self,domain,numEquations=1,numSolutions=1):
310           """
311           builds the required shape of the coefficient
312    
313           @param domain: domain on which the PDE uses the coefficient
314           @type domain: L{Domain<escript.Domain>}
315           @param numEquations: number of equations of the PDE
316           @type numEquations: C{int}
317           @param numSolutions: number of components of the PDE solution
318           @type numSolutions: C{int}
319           @return: shape of the coefficient
320           @rtype: C{tuple} of C{int} values
321           """
322           dim=domain.getDim()
323           s=()
324           for i in self.pattern:
325                 if i==self.BY_EQUATION:
326                    if numEquations>1: s=s+(numEquations,)
327                 elif i==self.BY_SOLUTION:
328                    if numSolutions>1: s=s+(numSolutions,)
329               else:               else:
330                  s=s+(dim,)                  s=s+(dim,)
331          return s         return s
   
 _PDECoefficientTypes={  
 "A"         : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,PDECoefficientType.DIM,PDECoefficientType.SOLUTION,PDECoefficientType.DIM),PDECoefficientType.OPERATOR),  
 "B"         : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,PDECoefficientType.DIM,PDECoefficientType.SOLUTION),PDECoefficientType.OPERATOR),  
 "C"         : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,PDECoefficientType.SOLUTION,PDECoefficientType.DIM),PDECoefficientType.OPERATOR),  
 "D"         : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,PDECoefficientType.SOLUTION),PDECoefficientType.OPERATOR),  
 "X"         : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,PDECoefficientType.DIM),PDECoefficientType.RIGHTHANDSIDE),  
 "Y"         : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,),PDECoefficientType.RIGHTHANDSIDE),  
 "d"         : PDECoefficientType(PDECoefficientType.BOUNDARY,(PDECoefficientType.EQUATION,PDECoefficientType.SOLUTION),PDECoefficientType.OPERATOR),  
 "y"         : PDECoefficientType(PDECoefficientType.BOUNDARY,(PDECoefficientType.EQUATION,),PDECoefficientType.RIGHTHANDSIDE),  
 "d_contact" : PDECoefficientType(PDECoefficientType.CONTACT,(PDECoefficientType.EQUATION,PDECoefficientType.SOLUTION),PDECoefficientType.OPERATOR),  
 "y_contact" : PDECoefficientType(PDECoefficientType.CONTACT,(PDECoefficientType.EQUATION,),PDECoefficientType.RIGHTHANDSIDE),  
 "r"         : PDECoefficientType(PDECoefficientType.CONTINUOUS,(PDECoefficientType.EQUATION,),PDECoefficientType.RIGHTHANDSIDE),  
 "q"         : PDECoefficientType(PDECoefficientType.CONTINUOUS,(PDECoefficientType.SOLUTION,),PDECoefficientType.BOTH),  
 }  
332    
333  class LinearPDE:  class LinearPDE(object):
334     """     """
335     @brief Class to define a linear PDE     This class is used to define a general linear, steady, second order PDE
336         for an unknown function M{u} on a given domain defined through a L{Domain<escript.Domain>} object.
    class to define a linear PDE of the form  
337    
338       -(A_{ijkl}u_{k,l})_{,j} -(B_{ijk}u_k)_{,j} + C_{ikl}u_{k,l} +D_{ik}u_k = - (X_{ij})_{,j} + Y_i     For a single PDE with a solution with a single component the linear PDE is defined in the following form:
339    
340       with boundary conditons:     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)}
341    
         n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_{ik}u_k = - n_j*X_{ij} + y_i  
342    
343      and contact conditions     where M{grad(F)} denotes the spatial derivative of M{F}. Einstein's summation convention,
344       ie. summation over indexes appearing twice in a term of a sum is performed, is used.
345       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
346       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
347       L{ReducedFunction<escript.ReducedFunction>}. It is also allowd to use objects that can be converted into
348       such L{Data<escript.Data>} objects. M{A} and M{A_reduced} are rank two, M{B_reduced}, M{C_reduced}, M{X_reduced}
349       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.
350    
351          n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_contact_{ik}[u_k] = - n_j*X_{ij} + y_contact_i     The following natural boundary conditions are considered:
352    
353      and constraints:     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}
354    
355           u_i=r_i where q_i>0     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>}.
356    
    """  
    DEFAULT_METHOD=util.DEFAULT_METHOD  
    DIRECT=util.DIRECT  
    CHOLEVSKY=util.CHOLEVSKY  
    PCG=util.PCG  
    CR=util.CR  
    CGS=util.CGS  
    BICGSTAB=util.BICGSTAB  
    SSOR=util.SSOR  
    GMRES=util.GMRES  
    PRES20=util.PRES20  
357    
358     def __init__(self,domain,numEquations=0,numSolutions=0):     Constraints for the solution prescribing the value of the solution at certain locations in the domain. They have the form
      """  
      @brief initializes a new linear PDE.  
359    
360       @param args     M{u=r}  where M{q>0}
361       """  
362       M{r} and M{q} are each scalar where M{q} is the characteristic function defining where the constraint is applied.
363       The constraints override any other condition set by the PDE or the boundary condition.
364    
365       The PDE is symmetrical if
366    
367       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]}
368    
369       For a system of PDEs and a solution with several components the PDE has the form
370    
371       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] }
372    
373       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.
374       The natural boundary conditions take the form:
375    
376       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]}
377    
378    
379       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>}.
380    
381       Constraints take the form
382    
383       M{u[i]=r[i]}  where  M{q[i]>0}
384    
385       M{r} and M{q} are each rank one. Notice that at some locations not necessarily all components must have a constraint.
386    
387       The system of PDEs is symmetrical if
388    
389            - M{A[i,j,k,l]=A[k,l,i,j]}
390            - M{A_reduced[i,j,k,l]=A_reduced[k,l,i,j]}
391            - M{B[i,j,k]=C[k,i,j]}
392            - M{B_reduced[i,j,k]=C_reduced[k,i,j]}
393            - M{D[i,k]=D[i,k]}
394            - M{D_reduced[i,k]=D_reduced[i,k]}
395            - M{d[i,k]=d[k,i]}
396            - M{d_reduced[i,k]=d_reduced[k,i]}
397    
398       L{LinearPDE} also supports solution discontinuities over a contact region in the domain. To specify the conditions across the
399       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
400       defined as
401    
402       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]}
403    
404       For the case of single solution component and single PDE M{J} is defined
405    
406       M{J_{j}=(A[i,j]+A_reduced[i,j])*grad(u)[j]+(B[i]+B_reduced[i])*u-X[i]-X_reduced[i]}
407    
408       In the context of discontinuities M{n} denotes the normal on the discontinuity pointing from side 0 towards side 1
409       calculated from L{getNormal<escript.FunctionSpace.getNormal>} of L{FunctionOnContactZero<escript.FunctionOnContactZero>}. For a system of PDEs
410       the contact condition takes the form
411    
412       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]}
413    
414       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
415       of the solution at side 1 and at side 0, denotes the jump of M{u} across discontinuity along the normal calcualted by
416       L{jump<util.jump>}.
417       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>}.
418       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>}.
419       In case of a single PDE and a single component solution the contact condition takes the form
420    
421       M{n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)}
422    
423       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>}
424    
425       @cvar DEFAULT: The default method used to solve the system of linear equations
426       @cvar DIRECT: The direct solver based on LDU factorization
427       @cvar CHOLEVSKY: The direct solver based on LDLt factorization (can only be applied for symmetric PDEs)
428       @cvar PCG: The preconditioned conjugate gradient method (can only be applied for symmetric PDEs)
429       @cvar CR: The conjugate residual method
430       @cvar CGS: The conjugate gardient square method
431       @cvar BICGSTAB: The stabilized BiConjugate Gradient method.
432       @cvar SSOR: The symmetric overrealaxtion method
433       @cvar ILU0: The incomplete LU factorization preconditioner  with no fill in
434       @cvar ILUT: The incomplete LU factorization preconditioner with will in
435       @cvar JACOBI: The Jacobi preconditioner
436       @cvar GMRES: The Gram-Schmidt minimum residual method
437       @cvar PRES20: Special GMRES with restart after 20 steps and truncation after 5 residuals
438       @cvar LUMPING: Matrix lumping.
439       @cvar NO_REORDERING: No matrix reordering allowed
440       @cvar MINIMUM_FILL_IN: Reorder matrix to reduce fill-in during factorization
441       @cvar NESTED_DISSECTION: Reorder matrix to improve load balancing during factorization
442       @cvar PASO: PASO solver package
443       @cvar SCSL: SGI SCSL solver library
444       @cvar MKL: Intel's MKL solver library
445       @cvar UMFPACK: the UMFPACK library
446       @cvar TRILINOS: the TRILINOS parallel solver class library from Sandia Natl Labs
447       @cvar ITERATIVE: The default iterative solver
448       @cvar AMG: algebraic multi grid
449       @cvar RILU: recursive ILU
450    
451       """
452       DEFAULT= 0
453       DIRECT= 1
454       CHOLEVSKY= 2
455       PCG= 3
456       CR= 4
457       CGS= 5
458       BICGSTAB= 6
459       SSOR= 7
460       ILU0= 8
461       ILUT= 9
462       JACOBI= 10
463       GMRES= 11
464       PRES20= 12
465       LUMPING= 13
466       NO_REORDERING= 17
467       MINIMUM_FILL_IN= 18
468       NESTED_DISSECTION= 19
469       SCSL= 14
470       MKL= 15
471       UMFPACK= 16
472       ITERATIVE= 20
473       PASO= 21
474       AMG= 22
475       RILU = 23
476       TRILINOS = 24
477    
478       SMALL_TOLERANCE=1.e-13
479       __PACKAGE_KEY="package"
480       __METHOD_KEY="method"
481       __SYMMETRY_KEY="symmetric"
482       __TOLERANCE_KEY="tolerance"
483       __PRECONDITIONER_KEY="preconditioner"
484    
485    
486       def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):
487         """
488         initializes a new linear PDE
489    
490         @param domain: domain of the PDE
491         @type domain: L{Domain<escript.Domain>}
492         @param numEquations: number of equations. If numEquations==None the number of equations
493                              is exracted from the PDE coefficients.
494         @param numSolutions: number of solution components. If  numSolutions==None the number of solution components
495                              is exracted from the PDE coefficients.
496         @param debug: if True debug informations are printed.
497    
498         """
499         super(LinearPDE, self).__init__()
500         #
501         #   the coefficients of the general PDE:
502         #
503         self.__COEFFICIENTS_OF_GENEARL_PDE={
504           "A"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM,PDECoefficient.BY_SOLUTION,PDECoefficient.BY_DIM),PDECoefficient.OPERATOR),
505           "B"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
506           "C"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION,PDECoefficient.BY_DIM),PDECoefficient.OPERATOR),
507           "D"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
508           "X"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM),PDECoefficient.RIGHTHANDSIDE),
509           "Y"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
510           "d"         : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
511           "y"         : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
512           "d_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
513           "y_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
514           "A_reduced"         : PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM,PDECoefficient.BY_SOLUTION,PDECoefficient.BY_DIM),PDECoefficient.OPERATOR),
515           "B_reduced"         : PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
516           "C_reduced"         : PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION,PDECoefficient.BY_DIM),PDECoefficient.OPERATOR),
517           "D_reduced"         : PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
518           "X_reduced"         : PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM),PDECoefficient.RIGHTHANDSIDE),
519           "Y_reduced"         : PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
520           "d_reduced"         : PDECoefficient(PDECoefficient.BOUNDARY_REDUCED,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
521           "y_reduced"         : PDECoefficient(PDECoefficient.BOUNDARY_REDUCED,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
522           "d_contact_reduced" : PDECoefficient(PDECoefficient.CONTACT_REDUCED,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
523           "y_contact_reduced" : PDECoefficient(PDECoefficient.CONTACT_REDUCED,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
524           "r"         : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_SOLUTION,),PDECoefficient.RIGHTHANDSIDE),
525           "q"         : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_SOLUTION,),PDECoefficient.BOTH)}
526    
527         # COEFFICIENTS can be overwritten by subclasses:
528         self.COEFFICIENTS=self.__COEFFICIENTS_OF_GENEARL_PDE
529         self.__altered_coefficients=False
530       # initialize attributes       # initialize attributes
531       self.__debug=None       self.__debug=debug
532       self.__domain=domain       self.__domain=domain
533       self.__numEquations=numEquations       self.__numEquations=numEquations
534       self.__numSolutions=numSolutions       self.__numSolutions=numSolutions
535       self.cleanCoefficients()       self.__resetSystem()
   
      self.__operator=escript.Operator()  
      self.__operator_isValid=False  
      self.__righthandside=escript.Data()  
      self.__righthandside_isValid=False  
      self.__solution=escript.Data()  
      self.__solution_isValid=False  
536    
537       # set some default values:       # set some default values:
538       self.__homogeneous_constraint=True       self.__reduce_equation_order=False
539       self.__row_function_space=escript.Solution(self.__domain)       self.__reduce_solution_order=False
      self.__column_function_space=escript.Solution(self.__domain)  
540       self.__tolerance=1.e-8       self.__tolerance=1.e-8
541       self.__solver_method=util.DEFAULT_METHOD       self.__solver_method=self.DEFAULT
542       self.__matrix_type=self.__domain.getSystemMatrixTypeId(util.DEFAULT_METHOD,False)       self.__solver_package=self.DEFAULT
543         self.__preconditioner=self.DEFAULT
544         self.__matrix_type=self.__domain.getSystemMatrixTypeId(self.DEFAULT,self.DEFAULT,False)
545       self.__sym=False       self.__sym=False
      self.__lumping=False  
546    
547     def getCoefficient(self,name):       self.resetCoefficients()
548         self.trace("PDE Coeffients are %s"%str(self.COEFFICIENTS.keys()))
549       # =============================================================================
550       #    general stuff:
551       # =============================================================================
552       def __str__(self):
553         """
554         returns string representation of the PDE
555    
556         @return: a simple representation of the PDE
557         @rtype: C{str}
558         """
559         return "<LinearPDE %d>"%id(self)
560       # =============================================================================
561       #    debug :
562       # =============================================================================
563       def setDebugOn(self):
564         """
565         switches on debugging
566       """       """
567       @brief return the value of the coefficient name       self.__debug=not None
568    
569       @param name     def setDebugOff(self):
570       """       """
571       return self.__coefficient[name]       switches off debugging
572         """
573         self.__debug=None
574    
575     def setValue(self,**coefficients):     def trace(self,text):
576        """       """
577        @brief sets new values to coefficients       print the text message if debugging is swiched on.
578         @param text: message
579         @type text: C{string}
580         """
581         if self.__debug: print "%s: %s"%(str(self),text)
582    
583        @param coefficients     # =============================================================================
584        """     # some service functions:
585        self._setValue(**coefficients)     # =============================================================================
586             def getDomain(self):
587         """
588         returns the domain of the PDE
589    
590     def _setValue(self,**coefficients):       @return: the domain of the PDE
591        """       @rtype: L{Domain<escript.Domain>}
592        @brief sets new values to coefficients       """
593         return self.__domain
594    
595        @param coefficients     def getDim(self):
596        """       """
597               returns the spatial dimension of the PDE
       # get the dictionary of the coefficinets been altered:  
       alteredCoefficients={}  
       for i,d in coefficients.iteritems():  
          if self.hasCoefficient(i):  
             if d == None:  
                 alteredCoefficients[i]=escript.Data()  
             elif isinstance(d,escript.Data):  
                 if d.isEmpty():  
                   alteredCoefficients[i]=escript.Data()  
                 else:  
                   alteredCoefficients[i]=escript.Data(d,self.getFunctionSpaceOfCoefficient(i))  
             else:  
                 alteredCoefficients[i]=escript.Data(d,self.getFunctionSpaceOfCoefficient(i))  
          else:  
             raise ValueError,"Attempt to set undefined coefficient %s"%i  
       # if numEquations and numSolutions is undefined we try identify their values based on the coefficients:  
       if self.__numEquations<1 or self.__numSolutions<1:  
             numEquations,numSolutions=identifyNumEquationsAndSolutions(self.getDomain().getDim(),alteredCoefficients)  
             if self.__numEquations<1 and numEquations>0: self.__numEquations=numEquations  
             if self.__numSolutions<1 and numSolutions>0: self.__numSolutions=numSolutions  
             if self.debug() and self.__numEquations>0: print "PDE Debug: identified number of equations is ",self.__numEquations  
             if self.debug() and self.__numSolutions>0: print "PDE Debug: identified number of solutions is ",self.__numSolutions  
598    
599        # now we check the shape of the coefficient if numEquations and numSolutions are set:       @return: the spatial dimension of the PDE domain
600        if  self.__numEquations>0 and  self.__numSolutions>0:       @rtype: C{int}
601           for i in self.__coefficient.iterkeys():       """
602               if alteredCoefficients.has_key(i) and not alteredCoefficients[i].isEmpty():       return self.getDomain().getDim()
                  if not self.getShapeOfCoefficient(i)==alteredCoefficients[i].getShape():  
                     raise ValueError,"Expected shape for coefficient %s is %s but actual shape is %s."%(i,self.getShapeOfCoefficient(i),alteredCoefficients[i].getShape())  
              else:  
                  if not self.__coefficient[i].isEmpty():  
                     if not self.getShapeOfCoefficient(i)==self.__coefficient[i].getShape():  
                        raise ValueError,"Expected shape for coefficient %s is %s but actual shape is %s."%(i,self.getShapeOfCoefficient(i),self.__coefficient[i].getShape())  
       # overwrite new values:  
       for i,d in alteredCoefficients.iteritems():  
          if self.debug(): print "PDE Debug: Coefficient %s has been altered."%i  
          self.__coefficient[i]=d  
          self.alteredCoefficient(i)  
   
       # reset the HomogeneousConstraintFlag:  
       self.__setHomogeneousConstraintFlag()  
       if not "q" in alteredCoefficients and not self.__homogeneous_constraint: self.__rebuildSystem()  
   
    def cleanCoefficients(self):  
      """  
      @brief resets all coefficients to default values.  
      """  
      self.__coefficient={}  
      for i in _PDECoefficientTypes.iterkeys():  
          self.__coefficient[i]=escript.Data()  
603    
604     def getShapeOfCoefficient(self,name):     def getNumEquations(self):
605       """       """
606       @brief return the shape of the coefficient name       returns the number of equations
607    
608       @param name       @return: the number of equations
609         @rtype: C{int}
610         @raise UndefinedPDEError: if the number of equations is not be specified yet.
611       """       """
612       if self.hasCoefficient(name):       if self.__numEquations==None:
613          return _PDECoefficientTypes[name].buildShape(self.getNumEquations(),self.getNumSolutions(),self.getDomain().getDim())           raise UndefinedPDEError,"Number of equations is undefined. Please specify argument numEquations."
614       else:       else:
615          raise ValueError,"Solution coefficient %s requested"%name           return self.__numEquations
616    
617     def getFunctionSpaceOfCoefficient(self,name):     def getNumSolutions(self):
618       """       """
619       @brief return the atoms of the coefficient name       returns the number of unknowns
620    
621       @param name       @return: the number of unknowns
622         @rtype: C{int}
623         @raise UndefinedPDEError: if the number of unknowns is not be specified yet.
624       """       """
625       if self.hasCoefficient(name):       if self.__numSolutions==None:
626          return _PDECoefficientTypes[name].getFunctionSpace(self.getDomain())          raise UndefinedPDEError,"Number of solution is undefined. Please specify argument numSolutions."
627       else:       else:
628          raise ValueError,"Solution coefficient %s requested"%name          return self.__numSolutions
629    
630     def alteredCoefficient(self,name):     def reduceEquationOrder(self):
631       """       """
632       @brief annonced that coefficient name has been changed       return status for order reduction for equation
633    
634       @param name       @return: return True is reduced interpolation order is used for the represenation of the equation
635         @rtype: L{bool}
636       """       """
637       if self.hasCoefficient(name):       return self.__reduce_equation_order
         if _PDECoefficientTypes[name].isAlteringOperator(): self.__rebuildOperator()  
         if _PDECoefficientTypes[name].isAlteringRightHandSide(): self.__rebuildRightHandSide()  
      else:  
         raise ValueError,"Solution coefficient %s requested"%name  
   
    def __setHomogeneousConstraintFlag(self):  
       """  
       @brief checks if the constraints are homogeneous and sets self.__homogeneous_constraint accordingly.  
       """  
       self.__homogeneous_constraint=True  
       q=self.getCoefficient("q")  
       r=self.getCoefficient("r")  
       if not q.isEmpty() and not r.isEmpty():  
          print (q*r).Lsup(), 1.e-13*r.Lsup()  
          if (q*r).Lsup()>=1.e-13*r.Lsup(): self.__homogeneous_constraint=False  
       if self.debug():  
            if self.__homogeneous_constraint:  
                print "PDE Debug: Constraints are homogeneous."  
            else:  
                print "PDE Debug: Constraints are inhomogeneous."  
   
638    
639     def hasCoefficient(self,name):     def reduceSolutionOrder(self):
640        """       """
641        @brief return true if name is the name of a coefficient       return status for order reduction for the solution
642    
643        @param name       @return: return True is reduced interpolation order is used for the represenation of the solution
644        """       @rtype: L{bool}
645        return self.__coefficient.has_key(name)       """
646         return self.__reduce_solution_order
647    
648     def getFunctionSpaceForEquation(self):     def getFunctionSpaceForEquation(self):
649       """       """
650       @brief return true if the test functions should use reduced order       returns the L{FunctionSpace<escript.FunctionSpace>} used to discretize the equation
651    
652         @return: representation space of equation
653         @rtype: L{FunctionSpace<escript.FunctionSpace>}
654       """       """
655       return self.__row_function_space       if self.reduceEquationOrder():
656             return escript.ReducedSolution(self.getDomain())
657         else:
658             return escript.Solution(self.getDomain())
659    
660     def getFunctionSpaceForSolution(self):     def getFunctionSpaceForSolution(self):
661       """       """
662       @brief return true if the interpolation of the solution should use reduced order       returns the L{FunctionSpace<escript.FunctionSpace>} used to represent the solution
663    
664         @return: representation space of solution
665         @rtype: L{FunctionSpace<escript.FunctionSpace>}
666       """       """
667       return self.__column_function_space       if self.reduceSolutionOrder():
668             return escript.ReducedSolution(self.getDomain())
669         else:
670             return escript.Solution(self.getDomain())
671    
    # ===== debug ==============================================================  
    def setDebugOn(self):  
        """  
        @brief  
        """  
        self.__debug=not None  
672    
673     def setDebugOff(self):     def getOperator(self):
674         """       """
675         @brief       provides access to the operator of the PDE
        """  
        self.__debug=None  
676    
677     def debug(self):       @return: the operator of the PDE
678         """       @rtype: L{Operator<escript.Operator>}
679         @brief returns true if the PDE is in the debug mode       """
680         """       m=self.getSystem()[0]
681         return self.__debug       if self.isUsingLumping():
682             return self.copyConstraint(1./m)
683         else:
684             return m
685    
686     #===== Lumping ===========================     def getRightHandSide(self):
687     def setLumpingOn(self):       """
688        """       provides access to the right hand side of the PDE
689        @brief indicates to use matrix lumping       @return: the right hand side of the PDE
690        """       @rtype: L{Data<escript.Data>}
691        if not self.isUsingLumping():       """
692           raise SystemError,"Lumping is not working yet! Talk to the experts"       r=self.getSystem()[1]
693           if self.debug() : print "PDE Debug: lumping is set on"       if self.isUsingLumping():
694           self.__rebuildOperator()           return self.copyConstraint(r)
695           self.__lumping=True       else:
696             return r
697    
698     def setLumpingOff(self):     def applyOperator(self,u=None):
699        """       """
700        @brief switches off matrix lumping       applies the operator of the PDE to a given u or the solution of PDE if u is not present.
       """  
       if self.isUsingLumping():  
          if self.debug() : print "PDE Debug: lumping is set off"  
          self.__rebuildOperator()  
          self.__lumping=False  
701    
702     def setLumping(self,flag=False):       @param u: argument of the operator. It must be representable in C{elf.getFunctionSpaceForSolution()}. If u is not present or equals L{None}
703        """                 the current solution is used.
704        @brief set the matrix lumping flag to flag       @type u: L{Data<escript.Data>} or None
705        """       @return: image of u
706        if flag:       @rtype: L{Data<escript.Data>}
707           self.setLumpingOn()       """
708        else:       if u==None:
709           self.setLumpingOff()          return self.getOperator()*self.getSolution()
710         else:
711            return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())
712    
713     def isUsingLumping(self):     def getResidual(self,u=None):
714         """
715         return the residual of u or the current solution if u is not present.
716    
717         @param u: argument in the residual calculation. It must be representable in C{elf.getFunctionSpaceForSolution()}. If u is not present or equals L{None}
718                   the current solution is used.
719         @type u: L{Data<escript.Data>} or None
720         @return: residual of u
721         @rtype: L{Data<escript.Data>}
722         """
723         return self.applyOperator(u)-self.getRightHandSide()
724    
725       def checkSymmetry(self,verbose=True):
726        """        """
727        @brief        test the PDE for symmetry.
728    
729          @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.
730          @type verbose: C{bool}
731          @return:  True if the PDE is symmetric.
732          @rtype: L{Data<escript.Data>}
733          @note: This is a very expensive operation. It should be used for degugging only! The symmetry flag is not altered.
734        """        """
735        return self.__lumping        verbose=verbose or self.__debug
736          out=True
737          if self.getNumSolutions()!=self.getNumEquations():
738             if verbose: print "non-symmetric PDE because of different number of equations and solutions"
739             out=False
740          else:
741             A=self.getCoefficientOfGeneralPDE("A")
742             if not A.isEmpty():
743                tol=util.Lsup(A)*self.SMALL_TOLERANCE
744                if self.getNumSolutions()>1:
745                   for i in range(self.getNumEquations()):
746                      for j in range(self.getDim()):
747                         for k in range(self.getNumSolutions()):
748                            for l in range(self.getDim()):
749                                if util.Lsup(A[i,j,k,l]-A[k,l,i,j])>tol:
750                                   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)
751                                   out=False
752                else:
753                   for j in range(self.getDim()):
754                      for l in range(self.getDim()):
755                         if util.Lsup(A[j,l]-A[l,j])>tol:
756                            if verbose: print "non-symmetric PDE because A[%d,%d]!=A[%d,%d]"%(j,l,l,j)
757                            out=False
758             B=self.getCoefficientOfGeneralPDE("B")
759             C=self.getCoefficientOfGeneralPDE("C")
760             if B.isEmpty() and not C.isEmpty():
761                if verbose: print "non-symmetric PDE because B is not present but C is"
762                out=False
763             elif not B.isEmpty() and C.isEmpty():
764                if verbose: print "non-symmetric PDE because C is not present but B is"
765                out=False
766             elif not B.isEmpty() and not C.isEmpty():
767                tol=(util.Lsup(B)+util.Lsup(C))*self.SMALL_TOLERANCE/2.
768                if self.getNumSolutions()>1:
769                   for i in range(self.getNumEquations()):
770                       for j in range(self.getDim()):
771                          for k in range(self.getNumSolutions()):
772                             if util.Lsup(B[i,j,k]-C[k,i,j])>tol:
773                                  if verbose: print "non-symmetric PDE because B[%d,%d,%d]!=C[%d,%d,%d]"%(i,j,k,k,i,j)
774                                  out=False
775                else:
776                   for j in range(self.getDim()):
777                      if util.Lsup(B[j]-C[j])>tol:
778                         if verbose: print "non-symmetric PDE because B[%d]!=C[%d]"%(j,j)
779                         out=False
780             if self.getNumSolutions()>1:
781               D=self.getCoefficientOfGeneralPDE("D")
782               if not D.isEmpty():
783                 tol=util.Lsup(D)*self.SMALL_TOLERANCE
784                 for i in range(self.getNumEquations()):
785                    for k in range(self.getNumSolutions()):
786                      if util.Lsup(D[i,k]-D[k,i])>tol:
787                          if verbose: print "non-symmetric PDE because D[%d,%d]!=D[%d,%d]"%(i,k,k,i)
788                          out=False
789               d=self.getCoefficientOfGeneralPDE("d")
790               if not d.isEmpty():
791                 tol=util.Lsup(d)*self.SMALL_TOLERANCE
792                 for i in range(self.getNumEquations()):
793                    for k in range(self.getNumSolutions()):
794                      if util.Lsup(d[i,k]-d[k,i])>tol:
795                          if verbose: print "non-symmetric PDE because d[%d,%d]!=d[%d,%d]"%(i,k,k,i)
796                          out=False
797               d_contact=self.getCoefficientOfGeneralPDE("d_contact")
798               if not d_contact.isEmpty():
799                 tol=util.Lsup(d_contact)*self.SMALL_TOLERANCE
800                 for i in range(self.getNumEquations()):
801                    for k in range(self.getNumSolutions()):
802                      if util.Lsup(d_contact[i,k]-d_contact[k,i])>tol:
803                          if verbose: print "non-symmetric PDE because d_contact[%d,%d]!=d_contact[%d,%d]"%(i,k,k,i)
804                          out=False
805             # and now the reduced coefficients
806             A_reduced=self.getCoefficientOfGeneralPDE("A_reduced")
807             if not A_reduced.isEmpty():
808                tol=util.Lsup(A_reduced)*self.SMALL_TOLERANCE
809                if self.getNumSolutions()>1:
810                   for i in range(self.getNumEquations()):
811                      for j in range(self.getDim()):
812                         for k in range(self.getNumSolutions()):
813                            for l in range(self.getDim()):
814                                if util.Lsup(A_reduced[i,j,k,l]-A_reduced[k,l,i,j])>tol:
815                                   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)
816                                   out=False
817                else:
818                   for j in range(self.getDim()):
819                      for l in range(self.getDim()):
820                         if util.Lsup(A_reduced[j,l]-A_reduced[l,j])>tol:
821                            if verbose: print "non-symmetric PDE because A_reduced[%d,%d]!=A_reduced[%d,%d]"%(j,l,l,j)
822                            out=False
823             B_reduced=self.getCoefficientOfGeneralPDE("B_reduced")
824             C_reduced=self.getCoefficientOfGeneralPDE("C_reduced")
825             if B_reduced.isEmpty() and not C_reduced.isEmpty():
826                if verbose: print "non-symmetric PDE because B_reduced is not present but C_reduced is"
827                out=False
828             elif not B_reduced.isEmpty() and C_reduced.isEmpty():
829                if verbose: print "non-symmetric PDE because C_reduced is not present but B_reduced is"
830                out=False
831             elif not B_reduced.isEmpty() and not C_reduced.isEmpty():
832                tol=(util.Lsup(B_reduced)+util.Lsup(C_reduced))*self.SMALL_TOLERANCE/2.
833                if self.getNumSolutions()>1:
834                   for i in range(self.getNumEquations()):
835                       for j in range(self.getDim()):
836                          for k in range(self.getNumSolutions()):
837                             if util.Lsup(B_reduced[i,j,k]-C_reduced[k,i,j])>tol:
838                                  if verbose: print "non-symmetric PDE because B_reduced[%d,%d,%d]!=C_reduced[%d,%d,%d]"%(i,j,k,k,i,j)
839                                  out=False
840                else:
841                   for j in range(self.getDim()):
842                      if util.Lsup(B_reduced[j]-C_reduced[j])>tol:
843                         if verbose: print "non-symmetric PDE because B_reduced[%d]!=C_reduced[%d]"%(j,j)
844                         out=False
845             if self.getNumSolutions()>1:
846               D_reduced=self.getCoefficientOfGeneralPDE("D_reduced")
847               if not D_reduced.isEmpty():
848                 tol=util.Lsup(D_reduced)*self.SMALL_TOLERANCE
849                 for i in range(self.getNumEquations()):
850                    for k in range(self.getNumSolutions()):
851                      if util.Lsup(D_reduced[i,k]-D_reduced[k,i])>tol:
852                          if verbose: print "non-symmetric PDE because D_reduced[%d,%d]!=D_reduced[%d,%d]"%(i,k,k,i)
853                          out=False
854               d_reduced=self.getCoefficientOfGeneralPDE("d_reduced")
855               if not d_reduced.isEmpty():
856                 tol=util.Lsup(d_reduced)*self.SMALL_TOLERANCE
857                 for i in range(self.getNumEquations()):
858                    for k in range(self.getNumSolutions()):
859                      if util.Lsup(d_reduced[i,k]-d_reduced[k,i])>tol:
860                          if verbose: print "non-symmetric PDE because d_reduced[%d,%d]!=d_reduced[%d,%d]"%(i,k,k,i)
861                          out=False
862               d_contact_reduced=self.getCoefficientOfGeneralPDE("d_contact_reduced")
863               if not d_contact_reduced.isEmpty():
864                 tol=util.Lsup(d_contact_reduced)*self.SMALL_TOLERANCE
865                 for i in range(self.getNumEquations()):
866                    for k in range(self.getNumSolutions()):
867                      if util.Lsup(d_contact_reduced[i,k]-d_contact_reduced[k,i])>tol:
868                          if verbose: print "non-symmetric PDE because d_contact_reduced[%d,%d]!=d_contact_reduced[%d,%d]"%(i,k,k,i)
869                          out=False
870          return out
871    
872     #============ method business =========================================================     def getSolution(self,**options):
    def setSolverMethod(self,solver=util.DEFAULT_METHOD):  
873         """         """
874         @brief sets a new solver         returns the solution of the PDE. If the solution is not valid the PDE is solved.
875    
876           @return: the solution
877           @rtype: L{Data<escript.Data>}
878           @param options: solver options
879           @keyword verbose: True to get some information during PDE solution
880           @type verbose: C{bool}
881           @keyword reordering: reordering scheme to be used during elimination. Allowed values are
882                                L{NO_REORDERING}, L{MINIMUM_FILL_IN}, L{NESTED_DISSECTION}
883           @keyword iter_max: maximum number of iteration steps allowed.
884           @keyword drop_tolerance: threshold for drupping in L{ILUT}
885           @keyword drop_storage: maximum of allowed memory in L{ILUT}
886           @keyword truncation: maximum number of residuals in L{GMRES}
887           @keyword restart: restart cycle length in L{GMRES}
888         """         """
889         if not solver==self.getSolverMethod():         if not self.__solution_isValid:
890              mat,f=self.getSystem()
891              if self.isUsingLumping():
892                 self.__solution=self.copyConstraint(f*mat)
893              else:
894                 options[self.__TOLERANCE_KEY]=self.getTolerance()
895                 options[self.__METHOD_KEY]=self.getSolverMethod()[0]
896                 options[self.__PRECONDITIONER_KEY]=self.getSolverMethod()[1]
897                 options[self.__PACKAGE_KEY]=self.getSolverPackage()
898                 options[self.__SYMMETRY_KEY]=self.isSymmetric()
899                 self.trace("PDE is resolved.")
900                 self.trace("solver options: %s"%str(options))
901                 self.__solution=mat.solve(f,options)
902              self.__solution_isValid=True
903           return self.__solution
904    
905       def getFlux(self,u=None):
906         """
907         returns the flux M{J} for a given M{u}
908    
909         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]}
910    
911         or
912    
913         M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[l]+(B[j]+B_reduced[j])u-X[j]-X_reduced[j]}
914    
915         @param u: argument in the flux. If u is not present or equals L{None} the current solution is used.
916         @type u: L{Data<escript.Data>} or None
917         @return: flux
918         @rtype: L{Data<escript.Data>}
919         """
920         if u==None: u=self.getSolution()
921         return util.tensormult(self.getCoefficientOfGeneralPDE("A"),util.grad(u,Funtion(self.getDomain))) \
922               +util.matrixmult(self.getCoefficientOfGeneralPDE("B"),u) \
923               -util.self.getCoefficientOfGeneralPDE("X") \
924               +util.tensormult(self.getCoefficientOfGeneralPDE("A_reduced"),util.grad(u,ReducedFuntion(self.getDomain))) \
925               +util.matrixmult(self.getCoefficientOfGeneralPDE("B_reduced"),u) \
926               -util.self.getCoefficientOfGeneralPDE("X_reduced")
927       # =============================================================================
928       #   solver settings:
929       # =============================================================================
930       def setSolverMethod(self,solver=None,preconditioner=None):
931           """
932           sets a new solver
933    
934           @param solver: sets a new solver method.
935           @type solver: one of L{DEFAULT}, L{ITERATIVE} L{DIRECT}, L{CHOLEVSKY}, L{PCG}, L{CR}, L{CGS}, L{BICGSTAB}, L{SSOR}, L{GMRES}, L{PRES20}, L{LUMPING}, L{AMG}
936           @param preconditioner: sets a new solver method.
937           @type preconditioner: one of L{DEFAULT}, L{JACOBI} L{ILU0}, L{ILUT},L{SSOR}, L{RILU}
938           """
939           if solver==None: solver=self.__solver_method
940           if preconditioner==None: preconditioner=self.__preconditioner
941           if solver==None: solver=self.DEFAULT
942           if preconditioner==None: preconditioner=self.DEFAULT
943           if not (solver,preconditioner)==self.getSolverMethod():
944             self.__solver_method=solver             self.__solver_method=solver
945             if self.debug() : print "PDE Debug: New solver is %s"%solver             self.__preconditioner=preconditioner
946             self.__checkMatrixType()             self.__checkMatrixType()
947               self.trace("New solver is %s"%self.getSolverMethodName())
948    
949       def getSolverMethodName(self):
950           """
951           returns the name of the solver currently used
952    
953           @return: the name of the solver currently used.
954           @rtype: C{string}
955           """
956    
957           m=self.getSolverMethod()
958           p=self.getSolverPackage()
959           method=""
960           if m[0]==self.DEFAULT: method="DEFAULT"
961           elif m[0]==self.DIRECT: method= "DIRECT"
962           elif m[0]==self.ITERATIVE: method= "ITERATIVE"
963           elif m[0]==self.CHOLEVSKY: method= "CHOLEVSKY"
964           elif m[0]==self.PCG: method= "PCG"
965           elif m[0]==self.CR: method= "CR"
966           elif m[0]==self.CGS: method= "CGS"
967           elif m[0]==self.BICGSTAB: method= "BICGSTAB"
968           elif m[0]==self.SSOR: method= "SSOR"
969           elif m[0]==self.GMRES: method= "GMRES"
970           elif m[0]==self.PRES20: method= "PRES20"
971           elif m[0]==self.LUMPING: method= "LUMPING"
972           elif m[0]==self.AMG: method= "AMG"
973           if m[1]==self.DEFAULT: method+="+DEFAULT"
974           elif m[1]==self.JACOBI: method+= "+JACOBI"
975           elif m[1]==self.ILU0: method+= "+ILU0"
976           elif m[1]==self.ILUT: method+= "+ILUT"
977           elif m[1]==self.SSOR: method+= "+SSOR"
978           elif m[1]==self.AMG: method+= "+AMG"
979           elif m[1]==self.RILU: method+= "+RILU"
980           if p==self.DEFAULT: package="DEFAULT"
981           elif p==self.PASO: package= "PASO"
982           elif p==self.MKL: package= "MKL"
983           elif p==self.SCSL: package= "SCSL"
984           elif p==self.UMFPACK: package= "UMFPACK"
985           elif p==self.TRILINOS: package= "TRILINOS"
986           else : method="unknown"
987           return "%s solver of %s package"%(method,package)
988    
989    
990     def getSolverMethod(self):     def getSolverMethod(self):
991         """         """
992         @brief returns the solver method         returns the solver method
993    
994           @return: the solver method currently be used.
995           @rtype: C{int}
996           """
997           return self.__solver_method,self.__preconditioner
998    
999       def setSolverPackage(self,package=None):
1000         """         """
1001         return self.__solver_method         sets a new solver package
1002    
1003           @param package: sets a new solver method.
1004           @type package: one of L{DEFAULT}, L{PASO} L{SCSL}, L{MKL}, L{UMFPACK}, L{TRILINOS}
1005           """
1006           if package==None: package=self.DEFAULT
1007           if not package==self.getSolverPackage():
1008               self.__solver_package=package
1009               self.__checkMatrixType()
1010               self.trace("New solver is %s"%self.getSolverMethodName())
1011    
1012       def getSolverPackage(self):
1013           """
1014           returns the package of the solver
1015    
1016           @return: the solver package currently being used.
1017           @rtype: C{int}
1018           """
1019           return self.__solver_package
1020    
1021       def isUsingLumping(self):
1022          """
1023          checks if matrix lumping is used a solver method
1024    
1025          @return: True is lumping is currently used a solver method.
1026          @rtype: C{bool}
1027          """
1028          return self.getSolverMethod()[0]==self.LUMPING
1029    
    #============ tolerance business =========================================================  
1030     def setTolerance(self,tol=1.e-8):     def setTolerance(self,tol=1.e-8):
1031         """         """
1032         @brief resets the tolerance to tol.         resets the tolerance for the solver method to tol where for an appropriate norm M{|.|}
1033    
1034           M{|L{getResidual}()|<tol*|L{getRightHandSide}()|}
1035    
1036           defines the stopping criterion.
1037    
1038           @param tol: new tolerance for the solver. If the tol is lower then the current tolerence
1039                       the system will be resolved.
1040           @type tol: positive C{float}
1041           @raise ValueError: if tolerance is not positive.
1042         """         """
1043         if not tol>0:         if not tol>0:
1044             raise ValueException,"Tolerance as to be positive"             raise ValueError,"Tolerance as to be positive"
1045         if tol<self.getTolerance(): self.__rebuildSolution()         if tol<self.getTolerance(): self.__invalidateSolution()
1046         if self.debug() : print "PDE Debug: New tolerance %e",tol         self.trace("New tolerance %e"%tol)
1047         self.__tolerance=tol         self.__tolerance=tol
1048         return         return
1049    
1050     def getTolerance(self):     def getTolerance(self):
1051         """         """
1052         @brief returns the tolerance set for the solution         returns the tolerance set for the solution
1053    
1054           @return: tolerance currently used.
1055           @rtype: C{float}
1056         """         """
1057         return self.__tolerance         return self.__tolerance
1058    
1059     #===== symmetry  flag ==========================     # =============================================================================
1060       #    symmetry  flag:
1061       # =============================================================================
1062     def isSymmetric(self):     def isSymmetric(self):
1063        """        """
1064        @brief returns true is the operator is considered to be symmetric        checks if symmetry is indicated.
1065    
1066          @return: True is a symmetric PDE is indicated, otherwise False is returned
1067          @rtype: C{bool}
1068        """        """
1069        return self.__sym        return self.__sym
1070    
1071     def setSymmetryOn(self):     def setSymmetryOn(self):
1072        """        """
1073        @brief sets the symmetry flag to true        sets the symmetry flag.
1074        """        """
1075        if not self.isSymmetric():        if not self.isSymmetric():
1076           if self.debug() : print "PDE Debug: Operator is set to be symmetric"           self.trace("PDE is set to be symmetric")
1077           self.__sym=True           self.__sym=True
1078           self.__checkMatrixType()           self.__checkMatrixType()
1079    
1080     def setSymmetryOff(self):     def setSymmetryOff(self):
1081        """        """
1082        @brief sets the symmetry flag to false        removes the symmetry flag.
1083        """        """
1084        if self.isSymmetric():        if self.isSymmetric():
1085           if self.debug() : print "PDE Debug: Operator is set to be unsymmetric"           self.trace("PDE is set to be unsymmetric")
1086           self.__sym=False           self.__sym=False
1087           self.__checkMatrixType()           self.__checkMatrixType()
1088    
1089     def setSymmetryTo(self,flag=False):     def setSymmetryTo(self,flag=False):
1090       """        """
1091       @brief sets the symmetry flag to flag        sets the symmetry flag to flag
1092    
1093       @param flag        @param flag: If flag, the symmetry flag is set otherwise the symmetry flag is released.
1094       """        @type flag: C{bool}
1095       if flag:        """
1096          self.setSymmetryOn()        if flag:
1097       else:           self.setSymmetryOn()
1098          self.setSymmetryOff()        else:
1099             self.setSymmetryOff()
1100    
1101     #===== order reduction ==========================     # =============================================================================
1102       # function space handling for the equation as well as the solution
1103       # =============================================================================
1104     def setReducedOrderOn(self):     def setReducedOrderOn(self):
1105       """       """
1106       @brief switches to on reduced order       switches on reduced order for solution and equation representation
1107    
1108         @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1109       """       """
1110       self.setReducedOrderForSolutionOn()       self.setReducedOrderForSolutionOn()
1111       self.setReducedOrderForEquationOn()       self.setReducedOrderForEquationOn()
1112    
1113     def setReducedOrderOff(self):     def setReducedOrderOff(self):
1114       """       """
1115       @brief switches to full order       switches off reduced order for solution and equation representation
1116    
1117         @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1118       """       """
1119       self.setReducedOrderForSolutionOff()       self.setReducedOrderForSolutionOff()
1120       self.setReducedOrderForEquationOff()       self.setReducedOrderForEquationOff()
1121    
1122     def setReducedOrderTo(self,flag=False):     def setReducedOrderTo(self,flag=False):
1123       """       """
1124       @brief sets order according to flag       sets order reduction for both solution and equation representation according to flag.
1125         @param flag: if flag is True, the order reduction is switched on for both  solution and equation representation, otherwise or
1126       @param flag                    if flag is not present order reduction is switched off
1127         @type flag: C{bool}
1128         @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1129       """       """
1130       self.setReducedOrderForSolutionTo(flag)       self.setReducedOrderForSolutionTo(flag)
1131       self.setReducedOrderForEquationTo(flag)       self.setReducedOrderForEquationTo(flag)
                                                                                                                                                             
1132    
1133     #===== order reduction solution ==========================  
1134     def setReducedOrderForSolutionOn(self):     def setReducedOrderForSolutionOn(self):
1135       """       """
1136       @brief switches to reduced order to interpolate solution       switches on reduced order for solution representation
1137    
1138         @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1139       """       """
1140       new_fs=escript.ReducedSolution(self.getDomain())       if not self.__reduce_solution_order:
1141       if self.getFunctionSpaceForSolution()!=new_fs:           if self.__altered_coefficients:
1142           if self.debug() : print "PDE Debug: Reduced order is used to interpolate solution."                raise RuntimeError,"order cannot be altered after coefficients have been defined."
1143           self.__column_function_space=new_fs           self.trace("Reduced order is used to solution representation.")
1144           self.__rebuildSystem(deep=True)           self.__reduce_solution_order=True
1145             self.__resetSystem()
1146    
1147     def setReducedOrderForSolutionOff(self):     def setReducedOrderForSolutionOff(self):
1148       """       """
1149       @brief switches to full order to interpolate solution       switches off reduced order for solution representation
1150    
1151         @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1152       """       """
1153       new_fs=escript.Solution(self.getDomain())       if self.__reduce_solution_order:
1154       if self.getFunctionSpaceForSolution()!=new_fs:           if self.__altered_coefficients:
1155           if self.debug() : print "PDE Debug: Full order is used to interpolate solution."                raise RuntimeError,"order cannot be altered after coefficients have been defined."
1156           self.__column_function_space=new_fs           self.trace("Full order is used to interpolate solution.")
1157           self.__rebuildSystem(deep=True)           self.__reduce_solution_order=False
1158             self.__resetSystem()
1159    
1160     def setReducedOrderForSolutionTo(self,flag=False):     def setReducedOrderForSolutionTo(self,flag=False):
1161       """       """
1162       @brief sets order for test functions according to flag       sets order for test functions according to flag
1163    
1164       @param flag       @param flag: if flag is True, the order reduction is switched on for solution representation, otherwise or
1165                      if flag is not present order reduction is switched off
1166         @type flag: C{bool}
1167         @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1168       """       """
1169       if flag:       if flag:
1170          self.setReducedOrderForSolutionOn()          self.setReducedOrderForSolutionOn()
1171       else:       else:
1172          self.setReducedOrderForSolutionOff()          self.setReducedOrderForSolutionOff()
1173                                                                                                                                                              
    #===== order reduction equation ==========================  
1174     def setReducedOrderForEquationOn(self):     def setReducedOrderForEquationOn(self):
1175       """       """
1176       @brief switches to reduced order for test functions       switches on reduced order for equation representation
1177    
1178         @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1179       """       """
1180       new_fs=escript.ReducedSolution(self.getDomain())       if not self.__reduce_equation_order:
1181       if self.getFunctionSpaceForEquation()!=new_fs:           if self.__altered_coefficients:
1182           if self.debug() : print "PDE Debug: Reduced order is used for test functions."                raise RuntimeError,"order cannot be altered after coefficients have been defined."
1183           self.__row_function_space=new_fs           self.trace("Reduced order is used for test functions.")
1184           self.__rebuildSystem(deep=True)           self.__reduce_equation_order=True
1185             self.__resetSystem()
1186    
1187     def setReducedOrderForEquationOff(self):     def setReducedOrderForEquationOff(self):
1188       """       """
1189       @brief switches to full order for test functions       switches off reduced order for equation representation
1190    
1191         @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1192       """       """
1193       new_fs=escript.Solution(self.getDomain())       if self.__reduce_equation_order:
1194       if self.getFunctionSpaceForEquation()!=new_fs:           if self.__altered_coefficients:
1195           if self.debug() : print "PDE Debug: Full order is used for test functions."                raise RuntimeError,"order cannot be altered after coefficients have been defined."
1196           self.__row_function_space=new_fs           self.trace("Full order is used for test functions.")
1197           self.__rebuildSystem(deep=True)           self.__reduce_equation_order=False
1198             self.__resetSystem()
1199    
1200     def setReducedOrderForEquationTo(self,flag=False):     def setReducedOrderForEquationTo(self,flag=False):
1201       """       """
1202       @brief sets order for test functions according to flag       sets order for test functions according to flag
1203    
1204       @param flag       @param flag: if flag is True, the order reduction is switched on for equation representation, otherwise or
1205                      if flag is not present order reduction is switched off
1206         @type flag: C{bool}
1207         @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1208       """       """
1209       if flag:       if flag:
1210          self.setReducedOrderForEquationOn()          self.setReducedOrderForEquationOn()
1211       else:       else:
1212          self.setReducedOrderForEquationOff()          self.setReducedOrderForEquationOff()
                                                                                                                                                             
1213    
1214     # ==== initialization =====================================================================     # =============================================================================
1215     def __makeNewOperator(self):     # private method:
1216       # =============================================================================
1217       def __checkMatrixType(self):
1218         """
1219         reassess the matrix type and, if a new matrix is needed, resets the system.
1220         """
1221         new_matrix_type=self.getDomain().getSystemMatrixTypeId(self.getSolverMethod()[0],self.getSolverPackage(),self.isSymmetric())
1222         if not new_matrix_type==self.__matrix_type:
1223             self.trace("Matrix type is now %d."%new_matrix_type)
1224             self.__matrix_type=new_matrix_type
1225             self.__resetSystem()
1226       #
1227       #   rebuild switches :
1228       #
1229       def __invalidateSolution(self):
1230         """         """
1231         @brief         indicates the PDE has to be resolved if the solution is requested
1232         """         """
1233           if self.__solution_isValid: self.trace("PDE has to be resolved.")
1234           self.__solution_isValid=False
1235    
1236       def __invalidateOperator(self):
1237           """
1238           indicates the operator has to be rebuilt next time it is used
1239           """
1240           if self.__operator_is_Valid: self.trace("Operator has to be rebuilt.")
1241           self.__invalidateSolution()
1242           self.__operator_is_Valid=False
1243    
1244       def __invalidateRightHandSide(self):
1245           """
1246           indicates the right hand side has to be rebuild next time it is used
1247           """
1248           if self.__righthandside_isValid: self.trace("Right hand side has to be rebuilt.")
1249           self.__invalidateSolution()
1250           self.__righthandside_isValid=False
1251    
1252       def __invalidateSystem(self):
1253           """
1254           annonced that everthing has to be rebuild:
1255           """
1256           if self.__righthandside_isValid: self.trace("System has to be rebuilt.")
1257           self.__invalidateSolution()
1258           self.__invalidateOperator()
1259           self.__invalidateRightHandSide()
1260    
1261       def __resetSystem(self):
1262           """
1263           annonced that everthing has to be rebuild:
1264           """
1265           self.trace("New System is built from scratch.")
1266           self.__operator=escript.Operator()
1267           self.__operator_is_Valid=False
1268           self.__righthandside=escript.Data()
1269           self.__righthandside_isValid=False
1270           self.__solution=escript.Data()
1271           self.__solution_isValid=False
1272       #
1273       #    system initialization:
1274       #
1275       def __getNewOperator(self):
1276           """
1277           returns an instance of a new operator
1278           """
1279           self.trace("New operator is allocated.")
1280         return self.getDomain().newOperator( \         return self.getDomain().newOperator( \
1281                             self.getNumEquations(), \                             self.getNumEquations(), \
1282                             self.getFunctionSpaceForEquation(), \                             self.getFunctionSpaceForEquation(), \
# Line 617  class LinearPDE: Line 1284  class LinearPDE:
1284                             self.getFunctionSpaceForSolution(), \                             self.getFunctionSpaceForSolution(), \
1285                             self.__matrix_type)                             self.__matrix_type)
1286    
1287     def __makeNewRightHandSide(self):     def __getNewRightHandSide(self):
1288         """         """
1289         @brief         returns an instance of a new right hand side
1290         """         """
1291         return escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)         self.trace("New right hand side is allocated.")
1292           if self.getNumEquations()>1:
1293               return escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)
1294           else:
1295               return escript.Data(0.,(),self.getFunctionSpaceForEquation(),True)
1296    
1297     def __makeNewSolution(self):     def __getNewSolution(self):
1298         """         """
1299         @brief         returns an instance of a new solution
1300         """         """
1301         return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)         self.trace("New solution is allocated.")
1302           if self.getNumSolutions()>1:
1303               return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)
1304           else:
1305               return escript.Data(0.,(),self.getFunctionSpaceForSolution(),True)
1306    
1307     def __getFreshOperator(self):     def __makeFreshSolution(self):
1308         """         """
1309         @brief         makes sure that the solution is instantiated and returns it initialized by zeros
1310         """         """
1311         if self.__operator.isEmpty():         if self.__solution.isEmpty():
1312             self.__operator=self.__makeNewOperator()             self.__solution=self.__getNewSolution()
            if self.debug() : print "PDE Debug: New operator allocated"  
1313         else:         else:
1314             self.__operator.setValue(0.)             self.__solution*=0
1315             if self.debug() : print "PDE Debug: Operator reset to zero"             self.trace("Solution is reset to zero.")
1316         return self.__operator         return self.__solution
1317    
1318     def __getFreshRightHandSide(self):     def __makeFreshRightHandSide(self):
1319         """         """
1320         @brief         makes sure that the right hand side is instantiated and returns it initialized by zeros
1321         """         """
1322         if self.__righthandside.isEmpty():         if self.__righthandside.isEmpty():
1323             self.__righthandside=self.__makeNewRightHandSide()             self.__righthandside=self.__getNewRightHandSide()
            if self.debug() : print "PDE Debug: New right hand side allocated"  
1324         else:         else:
1325             print "fix self.__righthandside*=0"             self.__righthandside.setToZero()
1326             self.__righthandside*=0.             self.trace("Right hand side is reset to zero.")
1327             if self.debug() : print "PDE Debug: Right hand side reset to zero"         return self.__righthandside
        return  self.__righthandside  
   
    # ==== rebuild switches =====================================================================  
    def __rebuildSolution(self,deep=False):  
        """  
        @brief indicates the PDE has to be reolved if the solution is requested  
        """  
        if self.__solution_isValid and self.debug() : print "PDE Debug: PDE has to be resolved."  
        self.__solution_isValid=False  
        if deep: self.__solution=escript.Data(deep)  
   
1328    
1329     def __rebuildOperator(self,deep=False):     def __makeFreshOperator(self):
1330         """         """
1331         @brief indicates the operator has to be rebuilt next time it is used         makes sure that the operator is instantiated and returns it initialized by zeros
1332         """         """
1333         if self.__operator_isValid and self.debug() : print "PDE Debug: Operator has to be rebuilt."         if self.__operator.isEmpty():
1334         self.__rebuildSolution(deep)             self.__operator=self.__getNewOperator()
1335         self.__operator_isValid=False         else:
1336         if deep: self.__operator=escript.Operator()             self.__operator.resetValues()
1337               self.trace("Operator reset to zero")
1338           return self.__operator
1339    
1340     def __rebuildRightHandSide(self,deep=False):     def __applyConstraint(self):
1341         """         """
1342         @brief indicates the right hand side has to be rebuild next time it is used         applies the constraints defined by q and r to the system
1343         """         """
1344         if self.__righthandside_isValid and self.debug() : print "PDE Debug: Right hand side has to be rebuilt."         if not self.isUsingLumping():
1345         self.__rebuildSolution(deep)            q=self.getCoefficientOfGeneralPDE("q")
1346         self.__righthandside_isValid=False            r=self.getCoefficientOfGeneralPDE("r")
1347         if not self.__homogeneous_constraint: self.__rebuildOperator()            if not q.isEmpty() and not self.__operator.isEmpty():
1348         if deep: self.__righthandside=escript.Data()               # q is the row and column mask to indicate where constraints are set:
1349                 row_q=escript.Data(q,self.getFunctionSpaceForEquation())
1350                 col_q=escript.Data(q,self.getFunctionSpaceForSolution())
1351                 u=self.__getNewSolution()
1352                 if r.isEmpty():
1353                    r_s=self.__getNewSolution()
1354                 else:
1355                    r_s=escript.Data(r,self.getFunctionSpaceForSolution())
1356                 u.copyWithMask(r_s,col_q)
1357                 if not self.__righthandside.isEmpty():
1358                    self.__righthandside-=self.__operator*u
1359                    self.__righthandside=self.copyConstraint(self.__righthandside)
1360                 self.__operator.nullifyRowsAndCols(row_q,col_q,1.)
1361       # =============================================================================
1362       # function giving access to coefficients of the general PDE:
1363       # =============================================================================
1364       def getCoefficientOfGeneralPDE(self,name):
1365         """
1366         return the value of the coefficient name of the general PDE.
1367    
1368         @note: This method is called by the assembling routine it can be overwritten
1369               to map coefficients of a particular PDE to the general PDE.
1370         @param name: name of the coefficient requested.
1371         @type name: C{string}
1372         @return: the value of the coefficient  name
1373         @rtype: L{Data<escript.Data>}
1374         @raise IllegalCoefficient: if name is not one of coefficients
1375                      M{A}, M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact},
1376                      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}.
1377         """
1378         if self.hasCoefficientOfGeneralPDE(name):
1379            return self.getCoefficient(name)
1380         else:
1381            raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1382    
1383     def __rebuildSystem(self,deep=False):     def hasCoefficientOfGeneralPDE(self,name):
        """  
        @brief annonced that all coefficient name has been changed  
        """  
        self.__rebuildSolution(deep)  
        self.__rebuildOperator(deep)  
        self.__rebuildRightHandSide(deep)  
     
    def __checkMatrixType(self):  
1384       """       """
1385       @brief reassess the matrix type and, if needed, initiates an operator rebuild       checks if name is a the name of a coefficient of the general PDE.
1386    
1387         @param name: name of the coefficient enquired.
1388         @type name: C{string}
1389         @return: True if name is the name of a coefficient of the general PDE. Otherwise False.
1390         @rtype: C{bool}
1391    
1392       """       """
1393       new_matrix_type=self.getDomain().getSystemMatrixTypeId(self.getSolverMethod(),self.isSymmetric())       return self.__COEFFICIENTS_OF_GENEARL_PDE.has_key(name)
      if not new_matrix_type==self.__matrix_type:  
          if self.debug() : print "PDE Debug: Matrix type is now %d."%new_matrix_type  
          self.__matrix_type=new_matrix_type  
          self.__rebuildOperator(deep=True)  
1394    
1395     #============ assembling =======================================================     def createCoefficientOfGeneralPDE(self,name):
1396     def __copyConstraint(self,u):       """
1397        """       returns a new instance of a coefficient for coefficient name of the general PDE
       @brief copies the constrint condition into u  
       """  
       q=self.getCoefficient("q")  
       r=self.getCoefficient("r")  
       if not q.isEmpty():  
           if r.isEmpty():  
              r2=escript.Data(0,u.getShape(),u.getFunctionSpace())  
           else:  
              r2=escript.Data(r,u.getFunctionSpace())  
           u.copyWithMask(r2,escript.Data(q,u.getFunctionSpace()))  
1398    
1399     def __applyConstraint(self,rhs_update=True):       @param name: name of the coefficient requested.
1400         """       @type name: C{string}
1401         @brief applies the constraints  defined by q and r to the system       @return: a coefficient name initialized to 0.
1402         """       @rtype: L{Data<escript.Data>}
1403         q=self.getCoefficient("q")       @raise IllegalCoefficient: if name is not one of coefficients
1404         r=self.getCoefficient("r")                    M{A}, M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact},
1405         if not q.isEmpty() and not self.__operator.isEmpty():                    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}.
1406            # q is the row and column mask to indicate where constraints are set:       """
1407            row_q=escript.Data(q,self.getFunctionSpaceForEquation())       if self.hasCoefficientOfGeneralPDE(name):
1408            col_q=escript.Data(q,self.getFunctionSpaceForSolution())          return escript.Data(0,self.getShapeOfCoefficientOfGeneralPDE(name),self.getFunctionSpaceForCoefficientOfGeneralPDE(name))
1409            u=self.__makeNewSolution()       else:
1410            if r.isEmpty():          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
              r_s=self.__makeNewSolution()  
           else:  
              r_s=escript.Data(r,self.getFunctionSpaceForSolution())  
           u.copyWithMask(r_s,col_q)  
           if not self.__righthandside.isEmpty() and rhs_update: self.__righthandside-=self.__operator*u  
           self.__operator.nullifyRowsAndCols(row_q,col_q,1.)  
1411    
1412     def getOperator(self):     def getFunctionSpaceForCoefficientOfGeneralPDE(self,name):
1413         """       """
1414         @brief returns the operator of the PDE       return the L{FunctionSpace<escript.FunctionSpace>} to be used for coefficient name of the general PDE
        """  
        if not self.__operator_isValid:  
            # some Constraints are applying for a lumpled stifness matrix:  
            if self.isUsingLumping():  
               if self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():  
                        raise TypeError,"Lumped matrix requires same order for equations and unknowns"  
               if not self.getCoefficient("A").isEmpty():  
                        raise Warning,"Lumped matrix does not allow coefficient A"  
               if not self.getCoefficient("B").isEmpty():  
                        raise Warning,"Lumped matrix does not allow coefficient B"  
               if not self.getCoefficient("C").isEmpty():  
                        raise Warning,"Lumped matrix does not allow coefficient C"  
               if not self.getCoefficient("D").isEmpty():  
                        raise Warning,"Lumped matrix does not allow coefficient D"  
               if self.debug() : print "PDE Debug: New lumped operator is built."  
               mat=self.__makeNewOperator(self)  
            else:  
               if self.debug() : print "PDE Debug: New operator is built."  
               mat=self.__getFreshOperator()  
   
            self.getDomain().addPDEToSystem(mat,escript.Data(), \  
                         self.getCoefficient("A"), \  
                         self.getCoefficient("B"), \  
                         self.getCoefficient("C"), \  
                         self.getCoefficient("D"), \  
                         escript.Data(), \  
                         escript.Data(), \  
                         self.getCoefficient("d"), \  
                         escript.Data(),\  
                         self.getCoefficient("d_contact"), \  
                         escript.Data())  
            if self.isUsingLumping():  
               self.__operator=mat*escript.Data(1,(self.getNumSolutions(),),self.getFunctionSpaceOfSolution(),True)  
            else:  
               self.__applyConstraint(rhs_update=False)  
            self.__operator_isValid=True  
        return self.__operator  
1415    
1416     def getRightHandSide(self,ignoreConstraint=False):       @param name: name of the coefficient enquired.
1417         """       @type name: C{string}
1418         @brief returns the right hand side of the PDE       @return: the function space to be used for coefficient name
1419         @rtype: L{FunctionSpace<escript.FunctionSpace>}
1420         @raise IllegalCoefficient: if name is not one of coefficients
1421                      M{A}, M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact},
1422                      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}.
1423         """
1424         if self.hasCoefficientOfGeneralPDE(name):
1425            return self.__COEFFICIENTS_OF_GENEARL_PDE[name].getFunctionSpace(self.getDomain())
1426         else:
1427            raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1428    
1429         @param ignoreConstraint     def getShapeOfCoefficientOfGeneralPDE(self,name):
1430         """       """
1431         if not self.__righthandside_isValid:       return the shape of the coefficient name of the general PDE
            if self.debug() : print "PDE Debug: New right hand side is built."  
            self.getDomain().addPDEToRHS(self.__getFreshRightHandSide(), \  
                          self.getCoefficient("X"), \  
                          self.getCoefficient("Y"),\  
                          self.getCoefficient("y"),\  
                          self.getCoefficient("y_contact"))  
            self.__righthandside_isValid=True  
            if ignoreConstraint: self.__copyConstraint(self.__righthandside)  
        return self.__righthandside  
1432    
1433     def getSystem(self):       @param name: name of the coefficient enquired.
1434         """       @type name: C{string}
1435         @brief       @return: the shape of the coefficient name
1436         """       @rtype: C{tuple} of C{int}
1437         if not self.__operator_isValid and not self.__righthandside_isValid:       @raise IllegalCoefficient: if name is not one of coefficients
1438            if self.debug() : print "PDE Debug: New PDE is built."                    M{A}, M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact},
1439            if self.isUsingLumping():                    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}.
1440                self.getRightHandSide(ignoreConstraint=True)       """
1441                self.getOperator()       if self.hasCoefficientOfGeneralPDE(name):
1442            else:          return self.__COEFFICIENTS_OF_GENEARL_PDE[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())
1443                self.getDomain().addPDEToSystem(self.__getFreshOperator(),self.__getFreshRightHandSide(), \       else:
1444                              self.getCoefficient("A"), \          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
                             self.getCoefficient("B"), \  
                             self.getCoefficient("C"), \  
                             self.getCoefficient("D"), \  
                             self.getCoefficient("X"), \  
                             self.getCoefficient("Y"), \  
                             self.getCoefficient("d"), \  
                             self.getCoefficient("y"), \  
                             self.getCoefficient("d_contact"), \  
                             self.getCoefficient("y_contact"))  
           self.__operator_isValid=True  
           self.__righthandside_isValid=True  
           self.__applyConstraint()  
           self.__copyConstraint(self.__righthandside)  
        elif not self.__operator_isValid:  
           self.getOperator()  
        elif not self.__righthandside_isValid:  
           self.getRightHandSide()  
        return (self.__operator,self.__righthandside)  
1445    
1446     def solve(self,**options):     # =============================================================================
1447        """     # functions giving access to coefficients of a particular PDE implementation:
1448        @brief solve the PDE     # =============================================================================
1449       def getCoefficient(self,name):
1450         """
1451         returns the value of the coefficient name
1452    
1453        @param options       @param name: name of the coefficient requested.
1454        """       @type name: C{string}
1455        mat,f=self.getSystem()       @return: the value of the coefficient name
1456        if self.isUsingLumping():       @rtype: L{Data<escript.Data>}
1457           out=f/mat       @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1458           self.__copyConstraint(out)       """
1459        else:       if self.hasCoefficient(name):
1460           options[util.TOLERANCE_KEY]=self.getTolerance()           return self.COEFFICIENTS[name].getValue()
1461           options[util.METHOD_KEY]=self.getSolverMethod()       else:
1462           options[util.SYMMETRY_KEY]=self.isSymmetric()          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
          if self.debug() : print "PDE Debug: solver options: ",options  
          out=mat.solve(f,options)  
       return out  
1463    
1464     def getSolution(self,**options):     def hasCoefficient(self,name):
1465         """       """
1466         @brief returns the solution of the PDE       return True if name is the name of a coefficient
1467    
1468         @param options       @param name: name of the coefficient enquired.
1469         """       @type name: C{string}
1470         if not self.__solution_isValid:       @return: True if name is the name of a coefficient of the general PDE. Otherwise False.
1471             if self.debug() : print "PDE Debug: PDE is resolved."       @rtype: C{bool}
            self.__solution=self.solve(**options)  
            self.__solution_isValid=True  
        return self.__solution  
    #============ some serivice functions  =====================================================  
    def getDomain(self):  
1472       """       """
1473       @brief returns the domain of the PDE       return self.COEFFICIENTS.has_key(name)
1474    
1475       def createCoefficient(self, name):
1476       """       """
1477       return self.__domain       create a L{Data<escript.Data>} object corresponding to coefficient name
1478    
1479     def getDim(self):       @return: a coefficient name initialized to 0.
1480         @rtype: L{Data<escript.Data>}
1481         @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1482       """       """
1483       @brief returns the spatial dimension of the PDE       if self.hasCoefficient(name):
1484            return escript.Data(0.,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))
1485         else:
1486            raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1487    
1488       def getFunctionSpaceForCoefficient(self,name):
1489       """       """
1490       return self.getDomain().getDim()       return the L{FunctionSpace<escript.FunctionSpace>} to be used for coefficient name
1491    
1492     def getNumEquations(self):       @param name: name of the coefficient enquired.
1493         @type name: C{string}
1494         @return: the function space to be used for coefficient name
1495         @rtype: L{FunctionSpace<escript.FunctionSpace>}
1496         @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1497       """       """
1498       @brief returns the number of equations       if self.hasCoefficient(name):
1499            return self.COEFFICIENTS[name].getFunctionSpace(self.getDomain())
1500         else:
1501            raise ValueError,"unknown coefficient %s requested"%name
1502       def getShapeOfCoefficient(self,name):
1503       """       """
1504       if self.__numEquations>0:       return the shape of the coefficient name
1505           return self.__numEquations  
1506         @param name: name of the coefficient enquired.
1507         @type name: C{string}
1508         @return: the shape of the coefficient name
1509         @rtype: C{tuple} of C{int}
1510         @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1511         """
1512         if self.hasCoefficient(name):
1513            return self.COEFFICIENTS[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())
1514       else:       else:
1515           raise ValueError,"Number of equations is undefined. Please specify argument numEquations."          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1516    
1517     def getNumSolutions(self):     def resetCoefficients(self):
1518       """       """
1519       @brief returns the number of unknowns       resets all coefficients to there default values.
1520       """       """
1521       if self.__numSolutions>0:       for i in self.COEFFICIENTS.iterkeys():
1522          return self.__numSolutions           self.COEFFICIENTS[i].resetValue()
1523    
1524       def alteredCoefficient(self,name):
1525         """
1526         announce that coefficient name has been changed
1527    
1528         @param name: name of the coefficient enquired.
1529         @type name: C{string}
1530         @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1531         @note: if name is q or r, the method will not trigger a rebuilt of the system as constraints are applied to the solved system.
1532         """
1533         if self.hasCoefficient(name):
1534            self.trace("Coefficient %s has been altered."%name)
1535            if not ((name=="q" or name=="r") and self.isUsingLumping()):
1536               if self.COEFFICIENTS[name].isAlteringOperator(): self.__invalidateOperator()
1537               if self.COEFFICIENTS[name].isAlteringRightHandSide(): self.__invalidateRightHandSide()
1538       else:       else:
1539          raise ValueError,"Number of solution is undefined. Please specify argument numSolutions."          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1540    
1541       def copyConstraint(self,u):
1542          """
1543          copies the constraint into u and returns u.
1544    
1545     def checkSymmetry(self):        @param u: a function of rank 0 is a single PDE is solved and of shape (numSolution,) for a system of PDEs
1546          @type u: L{Data<escript.Data>}
1547          @return: the input u modified by the constraints.
1548          @rtype: L{Data<escript.Data>}
1549          @warning: u is altered if it has the appropriate L{FunctionSpace<escript.FunctionSpace>}
1550        """        """
1551        @brief returns if the Operator is symmetric. This is a very expensive operation!!! The symmetry flag is not altered.        q=self.getCoefficientOfGeneralPDE("q")
1552          r=self.getCoefficientOfGeneralPDE("r")
1553          if not q.isEmpty():
1554             if u.isEmpty(): u=escript.Data(0.,q.getShape(),q.getFunctionSpace())
1555             if r.isEmpty():
1556                 r=escript.Data(0,u.getShape(),u.getFunctionSpace())
1557             else:
1558                 r=escript.Data(r,u.getFunctionSpace())
1559             u.copyWithMask(r,escript.Data(q,u.getFunctionSpace()))
1560          return u
1561    
1562       def setValue(self,**coefficients):
1563        """        """
1564        raise SystemError,"checkSymmetry is not implemented yet"        sets new values to coefficients
1565    
1566        return None        @param coefficients: new values assigned to coefficients
1567          @keyword A: value for coefficient A.
1568          @type A: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1569          @keyword A_reduced: value for coefficient A_reduced.
1570          @type A_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1571          @keyword B: value for coefficient B
1572          @type B: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1573          @keyword B_reduced: value for coefficient B_reduced
1574          @type B_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1575          @keyword C: value for coefficient C
1576          @type C: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1577          @keyword C_reduced: value for coefficient C_reduced
1578          @type C_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1579          @keyword D: value for coefficient D
1580          @type D: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1581          @keyword D_reduced: value for coefficient D_reduced
1582          @type D_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1583          @keyword X: value for coefficient X
1584          @type X: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1585          @keyword X_reduced: value for coefficient X_reduced
1586          @type X_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1587          @keyword Y: value for coefficient Y
1588          @type Y: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1589          @keyword Y_reduced: value for coefficient Y_reduced
1590          @type Y_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.Function>}.
1591          @keyword d: value for coefficient d
1592          @type d: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1593          @keyword d_reduced: value for coefficient d_reduced
1594          @type d_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
1595          @keyword y: value for coefficient y
1596          @type y: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1597          @keyword d_contact: value for coefficient d_contact
1598          @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>}.
1599          @keyword d_contact_reduced: value for coefficient d_contact_reduced
1600          @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>}.
1601          @keyword y_contact: value for coefficient y_contact
1602          @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>}.
1603          @keyword y_contact_reduced: value for coefficient y_contact_reduced
1604          @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>}.
1605          @keyword r: values prescribed to the solution at the locations of constraints
1606          @type r: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1607                   depending of reduced order is used for the solution.
1608          @keyword q: mask for location of constraints
1609          @type q: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1610                   depending of reduced order is used for the representation of the equation.
1611          @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1612          """
1613          # check if the coefficients are  legal:
1614          for i in coefficients.iterkeys():
1615             if not self.hasCoefficient(i):
1616                raise IllegalCoefficient,"Attempt to set unknown coefficient %s"%i
1617          # if the number of unknowns or equations is still unknown we try to estimate them:
1618          if self.__numEquations==None or self.__numSolutions==None:
1619             for i,d in coefficients.iteritems():
1620                if hasattr(d,"shape"):
1621                    s=d.shape
1622                elif hasattr(d,"getShape"):
1623                    s=d.getShape()
1624                else:
1625                    s=numarray.array(d).shape
1626                if s!=None:
1627                    # get number of equations and number of unknowns:
1628                    res=self.COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(self.getDomain(),s)
1629                    if res==None:
1630                        raise IllegalCoefficientValue,"Illegal shape %s of coefficient %s"%(s,i)
1631                    else:
1632                        if self.__numEquations==None: self.__numEquations=res[0]
1633                        if self.__numSolutions==None: self.__numSolutions=res[1]
1634          if self.__numEquations==None: raise UndefinedPDEError,"unidententified number of equations"
1635          if self.__numSolutions==None: raise UndefinedPDEError,"unidententified number of solutions"
1636          # now we check the shape of the coefficient if numEquations and numSolutions are set:
1637          for i,d in coefficients.iteritems():
1638            try:
1639               self.COEFFICIENTS[i].setValue(self.getDomain(),
1640                                             self.getNumEquations(),self.getNumSolutions(),
1641                                             self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
1642               self.alteredCoefficient(i)
1643            except IllegalCoefficientFunctionSpace,m:
1644                # if the function space is wrong then we try the reduced version:
1645                i_red=i+"_reduced"
1646                if (not i_red in coefficients.keys()) and i_red in self.COEFFICIENTS.keys():
1647                    try:
1648                        self.COEFFICIENTS[i_red].setValue(self.getDomain(),
1649                                                          self.getNumEquations(),self.getNumSolutions(),
1650                                                          self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
1651                        self.alteredCoefficient(i_red)
1652                    except IllegalCoefficientValue,m:
1653                        raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
1654                    except IllegalCoefficientFunctionSpace,m:
1655                        raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m))
1656                else:
1657                    raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m))
1658            except IllegalCoefficientValue,m:
1659               raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
1660          self.__altered_coefficients=True
1661          # check if the systrem is inhomogeneous:
1662          if len(coefficients)>0 and not self.isUsingLumping():
1663             q=self.getCoefficientOfGeneralPDE("q")
1664             r=self.getCoefficientOfGeneralPDE("r")
1665             homogeneous_constraint=True
1666             if not q.isEmpty() and not r.isEmpty():
1667                 if util.Lsup(q*r)>0.:
1668                   self.trace("Inhomogeneous constraint detected.")
1669                   self.__invalidateSystem()
1670    
1671     def getFlux(self,u):     def getSystem(self):
1672         """         """
1673         @brief returns the flux J_ij for a given u         return the operator and right hand side of the PDE
1674    
1675              J_ij=A_{ijkl}u_{k,l}+B_{ijk}u_k-X_{ij}         @return: the discrete version of the PDE
1676           @rtype: C{tuple} of L{Operator,<escript.Operator>} and L{Data<escript.Data>}.
1677           """
1678           if not self.__operator_is_Valid or not self.__righthandside_isValid:
1679              if self.isUsingLumping():
1680                  if not self.__operator_is_Valid:
1681                     if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():
1682                          raise TypeError,"Lumped matrix requires same order for equations and unknowns"
1683                     if not self.getCoefficientOfGeneralPDE("A").isEmpty():
1684                          raise ValueError,"coefficient A in lumped matrix may not be present."
1685                     if not self.getCoefficientOfGeneralPDE("B").isEmpty():
1686                          raise ValueError,"coefficient B in lumped matrix may not be present."
1687                     if not self.getCoefficientOfGeneralPDE("C").isEmpty():
1688                          raise ValueError,"coefficient C in lumped matrix may not be present."
1689                     if not self.getCoefficientOfGeneralPDE("d_contact").isEmpty():
1690                          raise ValueError,"coefficient d_contact in lumped matrix may not be present."
1691                     if not self.getCoefficientOfGeneralPDE("A_reduced").isEmpty():
1692                          raise ValueError,"coefficient A_reduced in lumped matrix may not be present."
1693                     if not self.getCoefficientOfGeneralPDE("B_reduced").isEmpty():
1694                          raise ValueError,"coefficient B_reduced in lumped matrix may not be present."
1695                     if not self.getCoefficientOfGeneralPDE("C_reduced").isEmpty():
1696                          raise ValueError,"coefficient C_reduced in lumped matrix may not be present."
1697                     if not self.getCoefficientOfGeneralPDE("d_contact_reduced").isEmpty():
1698                          raise ValueError,"coefficient d_contact_reduced in lumped matrix may not be present."
1699                     D=self.getCoefficientOfGeneralPDE("D")
1700                     d=self.getCoefficientOfGeneralPDE("d")
1701                     D_reduced=self.getCoefficientOfGeneralPDE("D_reduced")
1702                     d_reduced=self.getCoefficientOfGeneralPDE("d_reduced")
1703                     if not D.isEmpty():
1704                         if self.getNumSolutions()>1:
1705                            D_times_e=util.matrix_mult(D,numarray.ones((self.getNumSolutions(),)))
1706                         else:
1707                            D_times_e=D
1708                     else:
1709                        D_times_e=escript.Data()
1710                     if not d.isEmpty():
1711                         if self.getNumSolutions()>1:
1712                            d_times_e=util.matrix_mult(d,numarray.ones((self.getNumSolutions(),)))
1713                         else:
1714                            d_times_e=d
1715                     else:
1716                        d_times_e=escript.Data()
1717          
1718                     if not D_reduced.isEmpty():
1719                         if self.getNumSolutions()>1:
1720                            D_reduced_times_e=util.matrix_mult(D_reduced,numarray.ones((self.getNumSolutions(),)))
1721                         else:
1722                            D_reduced_times_e=D_reduced
1723                     else:
1724                        D_reduced_times_e=escript.Data()
1725                     if not d_reduced.isEmpty():
1726                         if self.getNumSolutions()>1:
1727                            d_reduced_times_e=util.matrix_mult(d_reduced,numarray.ones((self.getNumSolutions(),)))
1728                         else:
1729                            d_reduced_times_e=d_reduced
1730                     else:
1731                        d_reduced_times_e=escript.Data()
1732    
1733                     self.__operator=self.__getNewRightHandSide()
1734                     if hasattr(self.getDomain(), "addPDEToLumpedSystem") :
1735                        self.getDomain().addPDEToLumpedSystem(self.__operator, D_times_e, d_times_e)
1736                        self.getDomain().addPDEToLumpedSystem(self.__operator, D_reduced_times_e, d_reduced_times_e)
1737                     else:
1738                        self.getDomain().addPDEToRHS(self.__operator, \
1739                                                     escript.Data(), \
1740                                                     D_times_e, \
1741                                                     d_times_e,\
1742                                                     escript.Data())
1743                        self.getDomain().addPDEToRHS(self.__operator, \
1744                                                     escript.Data(), \
1745                                                     D_reduced_times_e, \
1746                                                     d_reduced_times_e,\
1747                                                     escript.Data())
1748                     self.__operator=1./self.__operator
1749                     self.trace("New lumped operator has been built.")
1750                     self.__operator_is_Valid=True
1751                  if not self.__righthandside_isValid:
1752                     self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \
1753                                   self.getCoefficientOfGeneralPDE("X"), \
1754                                   self.getCoefficientOfGeneralPDE("Y"),\
1755                                   self.getCoefficientOfGeneralPDE("y"),\
1756                                   self.getCoefficientOfGeneralPDE("y_contact"))
1757                     self.getDomain().addPDEToRHS(self.__righthandside, \
1758                                   self.getCoefficientOfGeneralPDE("X_reduced"), \
1759                                   self.getCoefficientOfGeneralPDE("Y_reduced"),\
1760                                   self.getCoefficientOfGeneralPDE("y_reduced"),\
1761                                   self.getCoefficientOfGeneralPDE("y_contact_reduced"))
1762                     self.trace("New right hand side as been built.")
1763                     self.__righthandside_isValid=True
1764              else:
1765                 if not self.__operator_is_Valid and not self.__righthandside_isValid:
1766                     self.getDomain().addPDEToSystem(self.__makeFreshOperator(),self.__makeFreshRightHandSide(), \
1767                                   self.getCoefficientOfGeneralPDE("A"), \
1768                                   self.getCoefficientOfGeneralPDE("B"), \
1769                                   self.getCoefficientOfGeneralPDE("C"), \
1770                                   self.getCoefficientOfGeneralPDE("D"), \
1771                                   self.getCoefficientOfGeneralPDE("X"), \
1772                                   self.getCoefficientOfGeneralPDE("Y"), \
1773                                   self.getCoefficientOfGeneralPDE("d"), \
1774                                   self.getCoefficientOfGeneralPDE("y"), \
1775                                   self.getCoefficientOfGeneralPDE("d_contact"), \
1776                                   self.getCoefficientOfGeneralPDE("y_contact"))
1777                     self.getDomain().addPDEToSystem(self.__operator,self.__righthandside, \
1778                                   self.getCoefficientOfGeneralPDE("A_reduced"), \
1779                                   self.getCoefficientOfGeneralPDE("B_reduced"), \
1780                                   self.getCoefficientOfGeneralPDE("C_reduced"), \
1781                                   self.getCoefficientOfGeneralPDE("D_reduced"), \
1782                                   self.getCoefficientOfGeneralPDE("X_reduced"), \
1783                                   self.getCoefficientOfGeneralPDE("Y_reduced"), \
1784                                   self.getCoefficientOfGeneralPDE("d_reduced"), \
1785                                   self.getCoefficientOfGeneralPDE("y_reduced"), \
1786                                   self.getCoefficientOfGeneralPDE("d_contact_reduced"), \
1787                                   self.getCoefficientOfGeneralPDE("y_contact_reduced"))
1788                     self.__applyConstraint()
1789                     self.__righthandside=self.copyConstraint(self.__righthandside)
1790                     self.trace("New system has been built.")
1791                     self.__operator_is_Valid=True
1792                     self.__righthandside_isValid=True
1793                 elif not self.__righthandside_isValid:
1794                     self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \
1795                                   self.getCoefficientOfGeneralPDE("X"), \
1796                                   self.getCoefficientOfGeneralPDE("Y"),\
1797                                   self.getCoefficientOfGeneralPDE("y"),\
1798                                   self.getCoefficientOfGeneralPDE("y_contact"))
1799                     self.getDomain().addPDEToRHS(self.__righthandside, \
1800                                   self.getCoefficientOfGeneralPDE("X_reduced"), \
1801                                   self.getCoefficientOfGeneralPDE("Y_reduced"),\
1802                                   self.getCoefficientOfGeneralPDE("y_reduced"),\
1803                                   self.getCoefficientOfGeneralPDE("y_contact_reduced"))
1804                     self.__righthandside=self.copyConstraint(self.__righthandside)
1805                     self.trace("New right hand side has been built.")
1806                     self.__righthandside_isValid=True
1807                 elif not self.__operator_is_Valid:
1808                     self.getDomain().addPDEToSystem(self.__makeFreshOperator(),escript.Data(), \
1809                                self.getCoefficientOfGeneralPDE("A"), \
1810                                self.getCoefficientOfGeneralPDE("B"), \
1811                                self.getCoefficientOfGeneralPDE("C"), \
1812                                self.getCoefficientOfGeneralPDE("D"), \
1813                                escript.Data(), \
1814                                escript.Data(), \
1815                                self.getCoefficientOfGeneralPDE("d"), \
1816                                escript.Data(),\
1817                                self.getCoefficientOfGeneralPDE("d_contact"), \
1818                                escript.Data())
1819                     self.getDomain().addPDEToSystem(self.__operator,escript.Data(), \
1820                                self.getCoefficientOfGeneralPDE("A_reduced"), \
1821                                self.getCoefficientOfGeneralPDE("B_reduced"), \
1822                                self.getCoefficientOfGeneralPDE("C_reduced"), \
1823                                self.getCoefficientOfGeneralPDE("D_reduced"), \
1824                                escript.Data(), \
1825                                escript.Data(), \
1826                                self.getCoefficientOfGeneralPDE("d_reduced"), \
1827                                escript.Data(),\
1828                                self.getCoefficientOfGeneralPDE("d_contact_reduced"), \
1829                                escript.Data())
1830                     self.__applyConstraint()
1831                     self.trace("New operator has been built.")
1832                     self.__operator_is_Valid=True
1833           return (self.__operator,self.__righthandside)
1834    
        @param u argument of the operator  
1835    
1836         """  class Poisson(LinearPDE):
1837         raise SystemError,"getFlux is not implemented yet"     """
1838         return None     Class to define a Poisson equation problem, which is genear L{LinearPDE} of the form
1839    
1840     def applyOperator(self,u):     M{-grad(grad(u)[j])[j] = f}
        """  
        @brief applies the operator of the PDE to a given solution u in weak from  
1841    
1842         @param u argument of the operator     with natural boundary conditons
1843    
1844         """     M{n[j]*grad(u)[j] = 0 }
        return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())  
                                                                                                                                                             
    def getResidual(self,u):  
        """  
        @brief return the residual of u in the weak from  
1845    
1846         @param u     and constraints:
1847         """  
1848         return self.applyOperator(u)-self.getRightHandSide()     M{u=0} where M{q>0}
1849    
1850       """
1851    
1852       def __init__(self,domain,debug=False):
1853         """
1854         initializes a new Poisson equation
1855    
1856         @param domain: domain of the PDE
1857         @type domain: L{Domain<escript.Domain>}
1858         @param debug: if True debug informations are printed.
1859    
1860         """
1861         super(Poisson, self).__init__(domain,1,1,debug)
1862         self.COEFFICIENTS={"f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1863                            "f_reduced": PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1864                            "q": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
1865         self.setSymmetryOn()
1866    
1867       def setValue(self,**coefficients):
1868         """
1869         sets new values to coefficients
1870    
1871         @param coefficients: new values assigned to coefficients
1872         @keyword f: value for right hand side M{f}
1873         @type f: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1874         @keyword q: mask for location of constraints
1875         @type q: any type that can be casted to rank zeo L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1876                   depending of reduced order is used for the representation of the equation.
1877         @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1878         """
1879         super(Poisson, self).setValue(**coefficients)
1880    
1881       def getCoefficientOfGeneralPDE(self,name):
1882         """
1883         return the value of the coefficient name of the general PDE
1884         @param name: name of the coefficient requested.
1885         @type name: C{string}
1886         @return: the value of the coefficient  name
1887         @rtype: L{Data<escript.Data>}
1888         @raise IllegalCoefficient: if name is not one of coefficients
1889                      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}.
1890         @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
1891         """
1892         if name == "A" :
1893             return escript.Data(util.kronecker(self.getDim()),escript.Function(self.getDomain()))
1894         elif name == "B" :
1895             return escript.Data()
1896         elif name == "C" :
1897             return escript.Data()
1898         elif name == "D" :
1899             return escript.Data()
1900         elif name == "X" :
1901             return escript.Data()
1902         elif name == "Y" :
1903             return self.getCoefficient("f")
1904         elif name == "d" :
1905             return escript.Data()
1906         elif name == "y" :
1907             return escript.Data()
1908         elif name == "d_contact" :
1909             return escript.Data()
1910         elif name == "y_contact" :
1911             return escript.Data()
1912         elif name == "A_reduced" :
1913             return escript.Data()
1914         elif name == "B_reduced" :
1915             return escript.Data()
1916         elif name == "C_reduced" :
1917             return escript.Data()
1918         elif name == "D_reduced" :
1919             return escript.Data()
1920         elif name == "X_reduced" :
1921             return escript.Data()
1922         elif name == "Y_reduced" :
1923             return self.getCoefficient("f_reduced")
1924         elif name == "d_reduced" :
1925             return escript.Data()
1926         elif name == "y_reduced" :
1927             return escript.Data()
1928         elif name == "d_contact_reduced" :
1929             return escript.Data()
1930         elif name == "y_contact_reduced" :
1931             return escript.Data()
1932         elif name == "r" :
1933             return escript.Data()
1934         elif name == "q" :
1935             return self.getCoefficient("q")
1936         else:
1937            raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1938    
1939    class Helmholtz(LinearPDE):
1940       """
1941       Class to define a Helmhotz equation problem, which is genear L{LinearPDE} of the form
1942    
1943       M{S{omega}*u - grad(k*grad(u)[j])[j] = f}
1944    
1945       with natural boundary conditons
1946    
1947       M{k*n[j]*grad(u)[j] = g- S{alpha}u }
1948    
1949       and constraints:
1950    
1951       M{u=r} where M{q>0}
1952    
 class Poisson(LinearPDE):  
1953     """     """
1954     @brief Class to define a Poisson equstion problem:  
1955                                                                                                                                                                   def __init__(self,domain,debug=False):
1956     class to define a linear PDE of the form       """
1957                                                                                                                                                                     initializes a new Poisson equation
1958          -u_{,jj} = f  
1959                                                                                                                                                                     @param domain: domain of the PDE
1960       with boundary conditons:       @type domain: L{Domain<escript.Domain>}
1961                                                                                                                                                                     @param debug: if True debug informations are printed.
1962          n_j*u_{,j} = 0  
1963                                                                                                                                                                     """
1964      and constraints:       super(Helmholtz, self).__init__(domain,1,1,debug)
1965                                                                                                                                                                     self.COEFFICIENTS={"omega": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
1966           u=0 where q>0                          "k": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
1967                                                                                                                                                                                        "f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1968     """                          "f_reduced": PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1969                            "alpha": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
1970     def __init__(self,domain=None,f=escript.Data(),q=escript.Data()):                          "g": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1971         LinearPDE.__init__(self,domain=identifyDomain(domain,{"f":f, "q":q}))                          "g_reduced": PDECoefficient(PDECoefficient.BOUNDARY_REDUCED,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1972         self._setValue(A=numarray.identity(self.getDomain().getDim()))                          "r": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH),
1973         self.setSymmetryOn()                          "q": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
1974         self.setValue(f,q)       self.setSymmetryOn()
1975    
1976     def setValue(self,f=escript.Data(),q=escript.Data()):     def setValue(self,**coefficients):
1977         self._setValue(Y=f,q=q)       """
1978         sets new values to coefficients
1979                                                                                                                                                              
1980  # $Log$       @param coefficients: new values assigned to coefficients
1981  # Revision 1.3  2004/12/17 07:43:10  jgs       @keyword omega: value for coefficient M{S{omega}}
1982  # *** empty log message ***       @type omega: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1983  #       @keyword k: value for coefficeint M{k}
1984  # Revision 1.1.2.3  2004/12/16 00:12:34  gross       @type k: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1985  # __init__ of LinearPDE does not accept any coefficients anymore       @keyword f: value for right hand side M{f}
1986  #       @type f: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1987  # Revision 1.1.2.2  2004/12/14 03:55:01  jgs       @keyword alpha: value for right hand side M{S{alpha}}
1988  # *** empty log message ***       @type alpha: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1989  #       @keyword g: value for right hand side M{g}
1990  # Revision 1.1.2.1  2004/12/12 22:53:47  gross       @type g: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1991  # linearPDE has been renamed LinearPDE       @keyword r: prescribed values M{r} for the solution in constraints.
1992  #       @type r: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1993  # Revision 1.1.1.1.2.7  2004/12/07 10:13:08  gross                 depending of reduced order is used for the representation of the equation.
1994  # GMRES added       @keyword q: mask for location of constraints
1995  #       @type q: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1996  # Revision 1.1.1.1.2.6  2004/12/07 03:19:50  gross                 depending of reduced order is used for the representation of the equation.
1997  # options for GMRES and PRES20 added       @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1998  #       """
1999  # Revision 1.1.1.1.2.5  2004/12/01 06:25:15  gross       super(Helmholtz, self).setValue(**coefficients)
2000  # some small changes  
2001  #     def getCoefficientOfGeneralPDE(self,name):
2002  # Revision 1.1.1.1.2.4  2004/11/24 01:50:21  gross       """
2003  # Finley solves 4M unknowns now       return the value of the coefficient name of the general PDE
2004  #  
2005  # Revision 1.1.1.1.2.3  2004/11/15 06:05:26  gross       @param name: name of the coefficient requested.
2006  # poisson solver added       @type name: C{string}
2007  #       @return: the value of the coefficient  name
2008  # Revision 1.1.1.1.2.2  2004/11/12 06:58:15  gross       @rtype: L{Data<escript.Data>}
2009  # a lot of changes to get the linearPDE class running: most important change is that there is no matrix format exposed to the user anymore. the format is chosen by the Domain according to the solver and symmetry       @raise IllegalCoefficient: if name is not one of coefficients
2010  #                    "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}.
2011  # Revision 1.1.1.1.2.1  2004/10/28 22:59:22  gross       @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
2012  # finley's RecTest.py is running now: problem in SystemMatrixAdapater fixed       """
2013  #       if name == "A" :
2014  # Revision 1.1.1.1  2004/10/26 06:53:56  jgs           return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))*self.getCoefficient("k")
2015  # initial import of project esys2       elif name == "B" :
2016  #           return escript.Data()
2017  # Revision 1.3.2.3  2004/10/26 06:43:48  jgs       elif name == "C" :
2018  # committing Lutz's and Paul's changes to brach jgs           return escript.Data()
2019  #       elif name == "D" :
2020  # Revision 1.3.4.1  2004/10/20 05:32:51  cochrane           return self.getCoefficient("omega")
2021  # Added incomplete Doxygen comments to files, or merely put the docstrings that already exist into Doxygen form.       elif name == "X" :
2022  #           return escript.Data()
2023  # Revision 1.3  2004/09/23 00:53:23  jgs       elif name == "Y" :
2024  # minor fixes           return self.getCoefficient("f")
2025  #       elif name == "d" :
2026  # Revision 1.1  2004/08/28 12:58:06  gross           return self.getCoefficient("alpha")
2027  # SimpleSolve is not running yet: problem with == of functionsspace       elif name == "y" :
2028  #           return self.getCoefficient("g")
2029  #       elif name == "d_contact" :
2030             return escript.Data()
2031         elif name == "y_contact" :
2032             return escript.Data()
2033         elif name == "A_reduced" :
2034             return escript.Data()
2035         elif name == "B_reduced" :
2036             return escript.Data()
2037         elif name == "C_reduced" :
2038             return escript.Data()
2039         elif name == "D_reduced" :
2040             return escript.Data()
2041         elif name == "X_reduced" :
2042             return escript.Data()
2043         elif name == "Y_reduced" :
2044             return self.getCoefficient("f_reduced")
2045         elif name == "d_reduced" :
2046             return escript.Data()
2047         elif name == "y_reduced" :
2048            return self.getCoefficient("g_reduced")
2049         elif name == "d_contact_reduced" :
2050             return escript.Data()
2051         elif name == "y_contact_reduced" :
2052             return escript.Data()
2053         elif name == "r" :
2054             return self.getCoefficient("r")
2055         elif name == "q" :
2056             return self.getCoefficient("q")
2057         else:
2058            raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
2059    
2060    class LameEquation(LinearPDE):
2061       """
2062       Class to define a Lame equation problem:
2063    
2064       M{-grad(S{mu}*(grad(u[i])[j]+grad(u[j])[i]))[j] - grad(S{lambda}*grad(u[k])[k])[j] = F_i -grad(S{sigma}[ij])[j] }
2065    
2066       with natural boundary conditons:
2067    
2068       M{n[j]*(S{mu}*(grad(u[i])[j]+grad(u[j])[i]) + S{lambda}*grad(u[k])[k]) = f_i +n[j]*S{sigma}[ij] }
2069    
2070       and constraints:
2071    
2072       M{u[i]=r[i]} where M{q[i]>0}
2073    
2074       """
2075    
2076       def __init__(self,domain,debug=False):
2077          super(LameEquation, self).__init__(domain,\
2078                                             domain.getDim(),domain.getDim(),debug)
2079          self.COEFFICIENTS={ "lame_lambda"  : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),
2080                              "lame_mu"      : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),
2081                              "F"            : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
2082                              "sigma"        : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM),PDECoefficient.RIGHTHANDSIDE),
2083                              "f"            : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
2084                              "r"            : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH),
2085                              "q"            : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
2086          self.setSymmetryOn()
2087    
2088       def setValues(self,**coefficients):
2089         """
2090         sets new values to coefficients
2091    
2092         @param coefficients: new values assigned to coefficients
2093         @keyword lame_mu: value for coefficient M{S{mu}}
2094         @type lame_mu: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2095         @keyword lame_lambda: value for coefficient M{S{lambda}}
2096         @type lame_lambda: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2097         @keyword F: value for internal force M{F}
2098         @type F: any type that can be casted to L{Vector<escript.Vector>} object on L{Function<escript.Function>}
2099         @keyword sigma: value for initial stress M{S{sigma}}
2100         @type sigma: any type that can be casted to L{Tensor<escript.Tensor>} object on L{Function<escript.Function>}
2101         @keyword f: value for extrenal force M{f}
2102         @type f: any type that can be casted to L{Vector<escript.Vector>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}
2103         @keyword r: prescribed values M{r} for the solution in constraints.
2104         @type r: any type that can be casted to L{Vector<escript.Vector>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2105                   depending of reduced order is used for the representation of the equation.
2106         @keyword q: mask for location of constraints
2107         @type q: any type that can be casted to L{Vector<escript.Vector>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2108                   depending of reduced order is used for the representation of the equation.
2109         @raise IllegalCoefficient: if an unknown coefficient keyword is used.
2110         """
2111         super(LameEquation, self).setValues(**coefficients)
2112    
2113       def getCoefficientOfGeneralPDE(self,name):
2114         """
2115         return the value of the coefficient name of the general PDE
2116    
2117         @param name: name of the coefficient requested.
2118         @type name: C{string}
2119         @return: the value of the coefficient  name
2120         @rtype: L{Data<escript.Data>}
2121         @raise IllegalCoefficient: if name is not one of coefficients
2122                      "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}.
2123         @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
2124         """
2125         if name == "A" :
2126             out =self.createCoefficientOfGeneralPDE("A")
2127             for i in range(self.getDim()):
2128               for j in range(self.getDim()):
2129                 out[i,i,j,j] += self.getCoefficient("lame_lambda")
2130                 out[i,j,j,i] += self.getCoefficient("lame_mu")
2131                 out[i,j,i,j] += self.getCoefficient("lame_mu")
2132             return out
2133         elif name == "B" :
2134             return escript.Data()
2135         elif name == "C" :
2136             return escript.Data()
2137         elif name == "D" :
2138             return escript.Data()
2139         elif name == "X" :
2140             return self.getCoefficient("sigma")
2141         elif name == "Y" :
2142             return self.getCoefficient("F")
2143         elif name == "d" :
2144             return escript.Data()
2145         elif name == "y" :
2146             return self.getCoefficient("f")
2147         elif name == "d_contact" :
2148             return escript.Data()
2149         elif name == "y_contact" :
2150             return escript.Data()
2151         elif name == "A_reduced" :
2152             return escript.Data()
2153         elif name == "B_reduced" :
2154             return escript.Data()
2155         elif name == "C_reduced" :
2156             return escript.Data()
2157         elif name == "D_reduced" :
2158             return escript.Data()
2159         elif name == "X_reduced" :
2160             return escript.Data()
2161         elif name == "Y_reduced" :
2162             return escript.Data()
2163         elif name == "d_reduced" :
2164             return escript.Data()
2165         elif name == "y_reduced" :
2166             return escript.Data()
2167         elif name == "d_contact_reduced" :
2168             return escript.Data()
2169         elif name == "y_contact_reduced" :
2170             return escript.Data()
2171         elif name == "r" :
2172             return self.getCoefficient("r")
2173         elif name == "q" :
2174             return self.getCoefficient("q")
2175         else:
2176            raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
2177    

Legend:
Removed from v.104  
changed lines
  Added in v.1302

  ViewVC Help
Powered by ViewVC 1.1.26