/[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 147 by jgs, Fri Aug 12 01:45:47 2005 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  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    __author__="Lutz Gross, l.gross@uq.edu.au"
27    __copyright__="""  Copyright (c) 2006 by ACcESS MNRF
28                        http://www.access.edu.au
29                    Primary Business: Queensland, Australia"""
30    __license__="""Licensed under the Open Software License version 3.0
31                 http://www.opensource.org/licenses/osl-3.0.php"""
32    __url__="http://www.iservo.edu.au/esys"
33    __version__="$Revision$"
34    __date__="$Date$"
35    
36    
37  def _CompTuple2(t1,t2):  class IllegalCoefficient(ValueError):
38     """     """
39     Compare two tuples     raised if an illegal coefficient of the general ar particular PDE is requested.
40       """
41       pass
42    
43     \param t1 The first tuple  class IllegalCoefficientValue(ValueError):
44     \param t2 The second tuple     """
45       raised if an incorrect value for a coefficient is used.
46     """     """
47       pass
48    
49     dif=t1[0]+t1[1]-(t2[0]+t2[1])  class IllegalCoefficientFunctionSpace(ValueError):
50     if dif<0: return 1     """
51     elif dif>0: return -1     raised if an incorrect function space for a coefficient is used.
52     else: return 0     """
   
 def ELMAN_RAMAGE(P):  
     return (P-1.).wherePositive()*0.5*(1.-1./(P+1.e-15))  
   
 def SIMPLIFIED_BROOK_HUGHES(P):  
     c=(P-3.).whereNegative()  
     return P/6.*c+1./2.*(1.-c)  
53    
54  def HALF(P):  class UndefinedPDEError(ValueError):
55      return escript.Scalar(0.5,P.getFunctionSpace())     """
56       raised if a PDE is not fully defined yet.
57       """
58       pass
59    
60  class PDECoefficient:  class PDECoefficient(object):
61      """      """
62      A class for PDE coefficients      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         """         """
97         Initialise a PDE Coefficient type         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
# Line 68  class PDECoefficient: Line 124  class PDECoefficient:
124         """         """
125         self.value=escript.Data()         self.value=escript.Data()
126    
127      def getFunctionSpace(self,domain):      def getFunctionSpace(self,domain,reducedEquationOrder=False,reducedSolutionOrder=False):
128         """         """
129         defines the FunctionSpace of the coefficient on the domain         defines the L{FunctionSpace<escript.FunctionSpace>} of the coefficient
130    
131         @param domain:         @param domain: domain on which the PDE uses the coefficient
132         """         @type domain: L{Domain<escript.Domain>}
133         if self.what==self.INTERIOR: return escript.Function(domain)         @param reducedEquationOrder: True to indicate that reduced order is used to represent the equation
134         elif self.what==self.BOUNDARY: return escript.FunctionOnBoundary(domain)         @type reducedEquationOrder: C{bool}
135         elif self.what==self.CONTACT: return escript.FunctionOnContactZero(domain)         @param reducedSolutionOrder: True to indicate that reduced order is used to represent the solution
136         elif self.what==self.CONTINUOUS: return escript.ContinuousFunction(domain)         @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):      def getValue(self):
161         """         """
162         returns the value of the coefficient:         returns the value of the coefficient
163    
164           @return:  value of the coefficient
165           @rtype:  L{Data<escript.Data>}
166         """         """
167         return self.value         return self.value
168        
169      def setValue(self,newValue):      def setValue(self,domain,numEquations=1,numSolutions=1,reducedEquationOrder=False,reducedSolutionOrder=False,newValue=None):
        """  
        set the value of the coefficient to new value  
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         self.value=newValue
203        
204      def isAlteringOperator(self):      def isAlteringOperator(self):
205          """          """
206      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 102  class PDECoefficient: Line 215  class PDECoefficient:
215    
216      def isAlteringRightHandSide(self):      def isAlteringRightHandSide(self):
217          """          """
218      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         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      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
332    
333  class LinearPDE:  class LinearPDE(object):
334     """     """
335     Class to handle 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     \f[     For a single PDE with a solution with a single component the linear PDE is defined in the following form:
      -(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  
    \f]  
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    
    \f[  
    n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_{ik}u_k = - n_j*X_{ij} + y_i  
    \f]  
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     \f[     The following natural boundary conditions are considered:
    n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_contact_{ik}[u_k] = - n_j*X_{ij} + y_contact_i  
    \f]  
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     \f[     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>}.
    u_i=r_i \quad \mathrm{where} \quad q_i>0  
    \f]  
356    
    """  
    TOL=1.e-13  
    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  
    ILU0=util.ILU0  
    JACOBI=util.JACOBI  
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
      """  
      initializes a new linear PDE.  
359    
360       @param args:     M{u=r}  where M{q>0}
361       """  
362       # COEFFICIENTS can be overwritten by subclasses:     M{r} and M{q} are each scalar where M{q} is the characteristic function defining where the constraint is applied.
363       self.__COEFFICIENTS={     The constraints override any other condition set by the PDE or the boundary condition.
364         "A"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM,PDECoefficient.SOLUTION,PDECoefficient.DIM),PDECoefficient.OPERATOR),  
365         "B"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),     The PDE is symmetrical if
366         "C"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION,PDECoefficient.DIM),PDECoefficient.OPERATOR),  
367         "D"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),     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         "X"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM),PDECoefficient.RIGHTHANDSIDE),  
369         "Y"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),     For a system of PDEs and a solution with several components the PDE has the form
370         "d"         : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),  
371         "y"         : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),     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         "d_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),  
373         "y_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),     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         "r"         : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),     The natural boundary conditions take the form:
375         "q"         : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.SOLUTION,),PDECoefficient.BOTH)}  
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       self.COEFFICIENTS=self.__COEFFICIENTS     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 createCoefficient(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       create a data object corresponding to coefficient name       switches on debugging
      @param name:  
566       """       """
567       return escript.Data(shape = getShapeOfCoefficient(name), \       self.__debug=not None
                          what = getFunctionSpaceForCoefficient(name))  
   
    def __del__(self):  
      pass  
568    
569     def getCoefficient(self,name):     def setDebugOff(self):
570       """       """
571       return the value of the parameter name       switches off debugging
   
      @param name:  
572       """       """
573       return self.COEFFICIENTS[name].getValue()       self.__debug=None
574    
575     def getCoefficientOfPDE(self,name):     def trace(self,text):
576       """       """
577       return the value of the coefficient name of the general PDE.       print the text message if debugging is swiched on.
578       This method is called by the assembling routine it can be       @param text: message
579       overwritten to map coefficients of a particualr PDE to the general PDE.       @type text: C{string}
   
      @param name:  
580       """       """
581       return self.getCoefficient(name)       if self.__debug: print "%s: %s"%(str(self),text)
   
    def hasCoefficient(self,name):  
       """  
       return true if name is the name of a coefficient  
582    
583        @param name:     # =============================================================================
584        """     # some service functions:
585        return self.COEFFICIENTS.has_key(name)     # =============================================================================
586       def getDomain(self):
    def hasPDECoefficient(self,name):  
       """  
       return true if name is the name of a coefficient  
   
       @param name:  
       """  
       return self.__COEFFICIENTS.has_key(name)  
   
    def getFunctionSpaceForEquation(self):  
587       """       """
588       return true if the test functions should use reduced order       returns the domain of the PDE
589    
590         @return: the domain of the PDE
591         @rtype: L{Domain<escript.Domain>}
592       """       """
593       return self.__row_function_space       return self.__domain
594    
595     def getFunctionSpaceForSolution(self):     def getDim(self):
596       """       """
597       return true if the interpolation of the solution should use reduced order       returns the spatial dimension of the PDE
598    
599         @return: the spatial dimension of the PDE domain
600         @rtype: C{int}
601       """       """
602       return self.__column_function_space       return self.getDomain().getDim()
603    
604     def setValue(self,**coefficients):     def getNumEquations(self):
605        """       """
606        sets new values to coefficients       returns the number of equations
607    
608        @param coefficients:       @return: the number of equations
609        """       @rtype: C{int}
610        self.__setValue(**coefficients)       @raise UndefinedPDEError: if the number of equations is not be specified yet.
611               """
612         if self.__numEquations==None:
613             raise UndefinedPDEError,"Number of equations is undefined. Please specify argument numEquations."
614         else:
615             return self.__numEquations
616    
617     def cleanCoefficients(self):     def getNumSolutions(self):
618       """       """
619       resets all coefficients to default values.       returns the number of unknowns
620    
621         @return: the number of unknowns
622         @rtype: C{int}
623         @raise UndefinedPDEError: if the number of unknowns is not be specified yet.
624       """       """
625       for i in self.COEFFICIENTS.iterkeys():       if self.__numSolutions==None:
626           self.COEFFICIENTS[i].resetValue()          raise UndefinedPDEError,"Number of solution is undefined. Please specify argument numSolutions."
627         else:
628            return self.__numSolutions
629    
630     def createNewCoefficient(self,name):     def reduceEquationOrder(self):
631       """       """
632       returns a new coefficient appropriate for coefficient name:       return status for order reduction for equation
633    
634         @return: return True is reduced interpolation order is used for the represenation of the equation
635         @rtype: L{bool}
636       """       """
637       return escript.Data(0,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))       return self.__reduce_equation_order
638          
639     def createNewCoefficientOfPDE(self,name):     def reduceSolutionOrder(self):
640       """       """
641       returns a new coefficient appropriate for coefficient name:       return status for order reduction for the solution
642    
643         @return: return True is reduced interpolation order is used for the represenation of the solution
644         @rtype: L{bool}
645       """       """
646       return escript.Data(0,self.getShapeOfCoefficientOfPDE(name),self.getFunctionSpaceForCoefficientOfPDE(name))       return self.__reduce_solution_order
647          
648     def getShapeOfCoefficientOfPDE(self,name):     def getFunctionSpaceForEquation(self):
649       """       """
650       return the shape of the coefficient name       returns the L{FunctionSpace<escript.FunctionSpace>} used to discretize the equation
651    
652       @param name:       @return: representation space of equation
653         @rtype: L{FunctionSpace<escript.FunctionSpace>}
654       """       """
655       if self.hasPDECoefficient(name):       if self.reduceEquationOrder():
656          return self.__COEFFICIENTS[name].buildShape(self.getNumEquations(),self.getNumSolutions(),self.getDomain().getDim())           return escript.ReducedSolution(self.getDomain())
657       else:       else:
658          raise ValueError,"Unknown coefficient %s requested"%name           return escript.Solution(self.getDomain())
659    
660     def getFunctionSpaceForCoefficientOfPDE(self,name):     def getFunctionSpaceForSolution(self):
661       """       """
662       return the atoms of the coefficient name       returns the L{FunctionSpace<escript.FunctionSpace>} used to represent the solution
663    
664       @param name:       @return: representation space of solution
665         @rtype: L{FunctionSpace<escript.FunctionSpace>}
666       """       """
667       if self.hasPDECoefficient(name):       if self.reduceSolutionOrder():
668          return self.__COEFFICIENTS[name].getFunctionSpace(self.getDomain())           return escript.ReducedSolution(self.getDomain())
669       else:       else:
670          raise ValueError,"unknown coefficient %s requested"%name           return escript.Solution(self.getDomain())
671    
672     def getShapeOfCoefficient(self,name):  
673       def getOperator(self):
674       """       """
675       return the shape of the coefficient name       provides access to the operator of the PDE
676    
677       @param name:       @return: the operator of the PDE
678         @rtype: L{Operator<escript.Operator>}
679       """       """
680       if self.hasCoefficient(name):       m=self.getSystem()[0]
681          return self.COEFFICIENTS[name].buildShape(self.getNumEquations(),self.getNumSolutions(),self.getDomain().getDim())       if self.isUsingLumping():
682             return self.copyConstraint(1./m)
683       else:       else:
684          raise ValueError,"Unknown coefficient %s requested"%name           return m
685    
686     def getFunctionSpaceForCoefficient(self,name):     def getRightHandSide(self):
687       """       """
688       return the atoms of the coefficient name       provides access to the right hand side of the PDE
689         @return: the right hand side of the PDE
690         @rtype: L{Data<escript.Data>}
691         """
692         r=self.getSystem()[1]
693         if self.isUsingLumping():
694             return self.copyConstraint(r)
695         else:
696             return r
697    
698       @param name:     def applyOperator(self,u=None):
699       """       """
700       if self.hasCoefficient(name):       applies the operator of the PDE to a given u or the solution of PDE if u is not present.
701          return self.COEFFICIENTS[name].getFunctionSpace(self.getDomain())  
702         @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         @type u: L{Data<escript.Data>} or None
705         @return: image of u
706         @rtype: L{Data<escript.Data>}
707         """
708         if u==None:
709            return self.getOperator()*self.getSolution()
710       else:       else:
711          raise ValueError,"unknown coefficient %s requested"%name          return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())
712    
713     def alteredCoefficient(self,name):     def getResidual(self,u=None):
714       """       """
715       announce that coefficient name has been changed       return the residual of u or the current solution if u is not present.
716    
717       @param name:       @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       if self.hasCoefficient(name):       return self.applyOperator(u)-self.getRightHandSide()
         if self.COEFFICIENTS[name].isAlteringOperator(): self.__rebuildOperator()  
         if self.COEFFICIENTS[name].isAlteringRightHandSide(): self.__rebuildRightHandSide()  
      else:  
         raise ValueError,"unknown coefficient %s requested"%name  
724    
725     # ===== debug ==============================================================     def checkSymmetry(self,verbose=True):
726     def setDebugOn(self):        """
727         """        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          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       def getSolution(self,**options):
873         """         """
874         self.__debug=not None         returns the solution of the PDE. If the solution is not valid the PDE is solved.
875    
876     def setDebugOff(self):         @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 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         self.__debug=None         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
945               self.__preconditioner=preconditioner
946               self.__checkMatrixType()
947               self.trace("New solver is %s"%self.getSolverMethodName())
948    
949     def debug(self):     def getSolverMethodName(self):
950         """         """
951         returns true if the PDE is in the debug mode         returns the name of the solver currently used
952    
953           @return: the name of the solver currently used.
954           @rtype: C{string}
955         """         """
        return self.__debug  
956    
957     #===== Lumping ===========================         m=self.getSolverMethod()
958     def setLumpingOn(self):         p=self.getSolverPackage()
959        """         method=""
960        indicates to use matrix lumping         if m[0]==self.DEFAULT: method="DEFAULT"
961        """         elif m[0]==self.DIRECT: method= "DIRECT"
962        if not self.isUsingLumping():         elif m[0]==self.ITERATIVE: method= "ITERATIVE"
963           if self.debug() : print "PDE Debug: lumping is set on"         elif m[0]==self.CHOLEVSKY: method= "CHOLEVSKY"
964           self.__rebuildOperator()         elif m[0]==self.PCG: method= "PCG"
965           self.__lumping=True         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    
    def setLumpingOff(self):  
       """  
       switches off matrix lumping  
       """  
       if self.isUsingLumping():  
          if self.debug() : print "PDE Debug: lumping is set off"  
          self.__rebuildOperator()  
          self.__lumping=False  
989    
990     def setLumping(self,flag=False):     def getSolverMethod(self):
991        """         """
992        set the matrix lumping flag to flag         returns the solver method
       """  
       if flag:  
          self.setLumpingOn()  
       else:  
          self.setLumpingOff()  
993    
994     def isUsingLumping(self):         @return: the solver method currently be used.
995        """         @rtype: C{int}
996                 """
997        """         return self.__solver_method,self.__preconditioner
       return self.__lumping  
998    
999     #============ method business =========================================================     def setSolverPackage(self,package=None):
    def setSolverMethod(self,solver=util.DEFAULT_METHOD):  
1000         """         """
1001         sets a new solver         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 not solver==self.getSolverMethod():         if package==None: package=self.DEFAULT
1007             self.__solver_method=solver         if not package==self.getSolverPackage():
1008             if self.debug() : print "PDE Debug: New solver is %s"%solver             self.__solver_package=package
1009             self.__checkMatrixType()             self.__checkMatrixType()
1010               self.trace("New solver is %s"%self.getSolverMethodName())
1011    
1012     def getSolverMethod(self):     def getSolverPackage(self):
1013         """         """
1014         returns the solver method         returns the package of the solver
1015    
1016           @return: the solver package currently being used.
1017           @rtype: C{int}
1018         """         """
1019         return self.__solver_method         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         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         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        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        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        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       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       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       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       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       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       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       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       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       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       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       # 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           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):     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 601  class LinearPDE: Line 1284  class LinearPDE:
1284                             self.getFunctionSpaceForSolution(), \                             self.getFunctionSpaceForSolution(), \
1285                             self.__matrix_type)                             self.__matrix_type)
1286    
1287     def __makeFreshRightHandSide(self):     def __getNewRightHandSide(self):
1288         """         """
1289           returns an instance of a new right hand side
1290         """         """
1291         if self.debug() : print "PDE Debug: New right hand side allocated"         self.trace("New right hand side is allocated.")
1292         if self.getNumEquations()>1:         if self.getNumEquations()>1:
1293             self.__righthandside=escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)             return escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)
1294         else:         else:
1295             self.__righthandside=escript.Data(0.,(),self.getFunctionSpaceForEquation(),True)             return escript.Data(0.,(),self.getFunctionSpaceForEquation(),True)
        return self.__righthandside  
1296    
1297     def __getNewSolution(self):     def __getNewSolution(self):
1298         """         """
1299           returns an instance of a new solution
1300         """         """
1301         if self.debug() : print "PDE Debug: New right hand side allocated"         self.trace("New solution is allocated.")
1302         if self.getNumSolutions()>1:         if self.getNumSolutions()>1:
1303             return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)             return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)
1304         else:         else:
1305             return escript.Data(0.,(),self.getFunctionSpaceForSolution(),True)             return escript.Data(0.,(),self.getFunctionSpaceForSolution(),True)
1306    
1307       def __makeFreshSolution(self):
1308           """
1309           makes sure that the solution is instantiated and returns it initialized by zeros
1310           """
1311           if self.__solution.isEmpty():
1312               self.__solution=self.__getNewSolution()
1313           else:
1314               self.__solution*=0
1315               self.trace("Solution is reset to zero.")
1316           return self.__solution
1317    
1318       def __makeFreshRightHandSide(self):
1319           """
1320           makes sure that the right hand side is instantiated and returns it initialized by zeros
1321           """
1322           if self.__righthandside.isEmpty():
1323               self.__righthandside=self.__getNewRightHandSide()
1324           else:
1325               self.__righthandside.setToZero()
1326               self.trace("Right hand side is reset to zero.")
1327           return self.__righthandside
1328    
1329     def __makeFreshOperator(self):     def __makeFreshOperator(self):
1330         """         """
1331           makes sure that the operator is instantiated and returns it initialized by zeros
1332         """         """
1333         if self.__operator.isEmpty():         if self.__operator.isEmpty():
1334             self.__operator=self.__getNewOperator()             self.__operator=self.__getNewOperator()
            if self.debug() : print "PDE Debug: New operator allocated"  
1335         else:         else:
1336             self.__operator.setValue(0.)             self.__operator.resetValues()
1337             self.__operator.resetSolver()             self.trace("Operator reset to zero")
            if self.debug() : print "PDE Debug: Operator reset to zero"  
1338         return self.__operator         return self.__operator
1339    
1340     #============ some serivice functions  =====================================================     def __applyConstraint(self):
1341     def getDomain(self):         """
1342           applies the constraints defined by q and r to the system
1343           """
1344           if not self.isUsingLumping():
1345              q=self.getCoefficientOfGeneralPDE("q")
1346              r=self.getCoefficientOfGeneralPDE("r")
1347              if not q.isEmpty() and not self.__operator.isEmpty():
1348                 # 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       returns the domain of the PDE       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 hasCoefficientOfGeneralPDE(self,name):
1384       """       """
1385       return self.__domain       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    
    def getDim(self):  
1392       """       """
1393       returns the spatial dimension of the PDE       return self.__COEFFICIENTS_OF_GENEARL_PDE.has_key(name)
1394    
1395       def createCoefficientOfGeneralPDE(self,name):
1396       """       """
1397       return self.getDomain().getDim()       returns a new instance of a coefficient for coefficient name of the general PDE
1398    
1399     def getNumEquations(self):       @param name: name of the coefficient requested.
1400         @type name: C{string}
1401         @return: a coefficient name initialized to 0.
1402         @rtype: L{Data<escript.Data>}
1403         @raise IllegalCoefficient: if name is not one of coefficients
1404                      M{A}, M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact},
1405                      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       """       """
1407       returns the number of equations       if self.hasCoefficientOfGeneralPDE(name):
1408            return escript.Data(0,self.getShapeOfCoefficientOfGeneralPDE(name),self.getFunctionSpaceForCoefficientOfGeneralPDE(name))
1409         else:
1410            raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1411    
1412       def getFunctionSpaceForCoefficientOfGeneralPDE(self,name):
1413       """       """
1414       if self.__numEquations>0:       return the L{FunctionSpace<escript.FunctionSpace>} to be used for coefficient name of the general PDE
1415           return self.__numEquations  
1416         @param name: name of the coefficient enquired.
1417         @type name: C{string}
1418         @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:       else:
1427           raise ValueError,"Number of equations is undefined. Please specify argument numEquations."          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1428    
1429     def getNumSolutions(self):     def getShapeOfCoefficientOfGeneralPDE(self,name):
1430       """       """
1431       returns the number of unknowns       return the shape of the coefficient name of the general PDE
1432    
1433         @param name: name of the coefficient enquired.
1434         @type name: C{string}
1435         @return: the shape of the coefficient name
1436         @rtype: C{tuple} of C{int}
1437         @raise IllegalCoefficient: if name is not one of coefficients
1438                      M{A}, M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact},
1439                      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       """       """
1441       if self.__numSolutions>0:       if self.hasCoefficientOfGeneralPDE(name):
1442          return self.__numSolutions          return self.__COEFFICIENTS_OF_GENEARL_PDE[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())
1443       else:       else:
1444          raise ValueError,"Number of solution is undefined. Please specify argument numSolutions."          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1445    
1446       # =============================================================================
1447       # functions giving access to coefficients of a particular PDE implementation:
1448       # =============================================================================
1449       def getCoefficient(self,name):
1450         """
1451         returns the value of the coefficient name
1452    
1453     def checkSymmetry(self,verbose=True):       @param name: name of the coefficient requested.
1454        """       @type name: C{string}
1455        returns if the Operator is symmetric. This is a very expensive operation!!! The symmetry flag is not altered.       @return: the value of the coefficient name
1456        """       @rtype: L{Data<escript.Data>}
1457        verbose=verbose or self.debug()       @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1458        out=True       """
1459        if self.getNumSolutions()!=self.getNumEquations():       if self.hasCoefficient(name):
1460           if verbose: print "non-symmetric PDE because of different number of equations and solutions"           return self.COEFFICIENTS[name].getValue()
1461           out=False       else:
1462        else:          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
          A=self.getCoefficientOfPDE("A")  
          if not A.isEmpty():  
             tol=util.Lsup(A)*self.TOL  
             if self.getNumSolutions()>1:  
                for i in range(self.getNumEquations()):  
                   for j in range(self.getDim()):  
                      for k in range(self.getNumSolutions()):  
                         for l in range(self.getDim()):  
                             if util.Lsup(A[i,j,k,l]-A[k,l,i,j])>tol:  
                                if verbose: print "non-symmetric PDE because A[%d,%d,%d,%d]!=A[%d,%d,%d,%d]"%(i,j,k,l,k,l,i,j)  
                                out=False  
             else:  
                for j in range(self.getDim()):  
                   for l in range(self.getDim()):  
                      if util.Lsup(A[j,l]-A[l,j])>tol:  
                         if verbose: print "non-symmetric PDE because A[%d,%d]!=A[%d,%d]"%(j,l,l,j)  
                         out=False  
          B=self.getCoefficientOfPDE("B")  
          C=self.getCoefficientOfPDE("C")  
          if B.isEmpty() and not C.isEmpty():  
             if verbose: print "non-symmetric PDE because B is not present but C is"  
             out=False  
          elif not B.isEmpty() and C.isEmpty():  
             if verbose: print "non-symmetric PDE because C is not present but B is"  
             out=False  
          elif not B.isEmpty() and not C.isEmpty():  
             tol=(util.Lsup(B)+util.Lsup(C))*self.TOL/2.  
             if self.getNumSolutions()>1:  
                for i in range(self.getNumEquations()):  
                    for j in range(self.getDim()):  
                       for k in range(self.getNumSolutions()):  
                          if util.Lsup(B[i,j,k]-C[k,i,j])>tol:  
                               if verbose: print "non-symmetric PDE because B[%d,%d,%d]!=C[%d,%d,%d]"%(i,j,k,k,i,j)  
                               out=False  
             else:  
                for j in range(self.getDim()):  
                   if util.Lsup(B[j]-C[j])>tol:  
                      if verbose: print "non-symmetric PDE because B[%d]!=C[%d]"%(j,j)  
                      out=False  
          if self.getNumSolutions()>1:  
            D=self.getCoefficientOfPDE("D")  
            if not D.isEmpty():  
              tol=util.Lsup(D)*self.TOL  
              for i in range(self.getNumEquations()):  
                 for k in range(self.getNumSolutions()):  
                   if util.Lsup(D[i,k]-D[k,i])>tol:  
                       if verbose: print "non-symmetric PDE because D[%d,%d]!=D[%d,%d]"%(i,k,k,i)  
                       out=False  
               
       return out  
1463    
1464     def getFlux(self,u):     def hasCoefficient(self,name):
1465         """       """
1466         returns the flux J_ij for a given u       return True if name is the name of a coefficient
1467    
1468         \f[       @param name: name of the coefficient enquired.
1469         J_ij=A_{ijkl}u_{k,l}+B_{ijk}u_k-X_{ij}       @type name: C{string}
1470         \f]       @return: True if name is the name of a coefficient of the general PDE. Otherwise False.
1471         @rtype: C{bool}
1472         """
1473         return self.COEFFICIENTS.has_key(name)
1474    
1475         @param u: argument of the operator     def createCoefficient(self, name):
1476         """       """
1477         raise SystemError,"getFlux is not implemented yet"       create a L{Data<escript.Data>} object corresponding to coefficient name
        return None  
1478    
1479     def applyOperator(self,u):       @return: a coefficient name initialized to 0.
1480         """       @rtype: L{Data<escript.Data>}
1481         applies the operator of the PDE to a given solution u in weak from       @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1482         """
1483         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         @param u: argument of the operator     def getFunctionSpaceForCoefficient(self,name):
1489         """       """
1490         return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())       return the L{FunctionSpace<escript.FunctionSpace>} to be used for coefficient name
                                                                                                                                                             
    def getResidual(self,u):  
        """  
        return the residual of u in the weak from  
1491    
1492         @param u:       @param name: name of the coefficient enquired.
1493         """       @type name: C{string}
1494         return self.applyOperator(u)-self.getRightHandSide()       @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         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         return the shape of the coefficient name
1505    
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:
1515            raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1516    
1517     def __setValue(self,**coefficients):     def resetCoefficients(self):
1518         """
1519         resets all coefficients to there default values.
1520         """
1521         for i in self.COEFFICIENTS.iterkeys():
1522             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:
1539            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          @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          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        sets new values to coefficient        sets new values to coefficients
1565    
1566        @param coefficients:        @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:        # check if the coefficients are  legal:
1614        for i in coefficients.iterkeys():        for i in coefficients.iterkeys():
1615           if not self.hasCoefficient(i):           if not self.hasCoefficient(i):
1616              raise ValueError,"Attempt to set unknown coefficient %s"%i              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:        # if the number of unknowns or equations is still unknown we try to estimate them:
1618        if self.__numEquations<1 or self.__numSolutions<1:        if self.__numEquations==None or self.__numSolutions==None:
1619           for i,d in coefficients.iteritems():           for i,d in coefficients.iteritems():
1620              if hasattr(d,"shape"):              if hasattr(d,"shape"):
1621                  s=d.shape                  s=d.shape
# Line 775  class LinearPDE: Line 1625  class LinearPDE:
1625                  s=numarray.array(d).shape                  s=numarray.array(d).shape
1626              if s!=None:              if s!=None:
1627                  # get number of equations and number of unknowns:                  # get number of equations and number of unknowns:
1628                  res=self.COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(s,self.getDim())                  res=self.COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(self.getDomain(),s)
1629                  if res==None:                  if res==None:
1630                      raise ValueError,"Illegal shape %s of coefficient %s"%(s,i)                      raise IllegalCoefficientValue,"Illegal shape %s of coefficient %s"%(s,i)
1631                  else:                  else:
1632                      if self.__numEquations<1: self.__numEquations=res[0]                      if self.__numEquations==None: self.__numEquations=res[0]
1633                      if self.__numSolutions<1: self.__numSolutions=res[1]                      if self.__numSolutions==None: self.__numSolutions=res[1]
1634        if self.__numEquations<1: raise ValueError,"unidententified number of equations"        if self.__numEquations==None: raise UndefinedPDEError,"unidententified number of equations"
1635        if self.__numSolutions<1: raise ValueError,"unidententified number of solutions"        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:        # now we check the shape of the coefficient if numEquations and numSolutions are set:
1637        for i,d in coefficients.iteritems():        for i,d in coefficients.iteritems():
1638          if d==None:          try:
1639               d2=escript.Data()             self.COEFFICIENTS[i].setValue(self.getDomain(),
1640          elif isinstance(d,escript.Data):                                           self.getNumEquations(),self.getNumSolutions(),
1641               if d.isEmpty():                                           self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
1642                  d2=d             self.alteredCoefficient(i)
1643               else:          except IllegalCoefficientFunctionSpace,m:
1644                  d2=escript.Data(d,self.getFunctionSpaceForCoefficient(i))              # if the function space is wrong then we try the reduced version:
1645          else:              i_red=i+"_reduced"
1646                d2=escript.Data(d,self.getFunctionSpaceForCoefficient(i))              if (not i_red in coefficients.keys()) and i_red in self.COEFFICIENTS.keys():
1647          if not d2.isEmpty():                  try:
1648             if not self.getShapeOfCoefficient(i)==d2.getShape():                      self.COEFFICIENTS[i_red].setValue(self.getDomain(),
1649                 raise ValueError,"Expected shape for coefficient %s is %s but actual shape is %s."%(i,self.getShapeOfCoefficient(i),d2.getShape())                                                        self.getNumEquations(),self.getNumSolutions(),
1650          # overwrite new values:                                                        self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
1651          if self.debug(): print "PDE Debug: Coefficient %s has been altered."%i                      self.alteredCoefficient(i_red)
1652          self.COEFFICIENTS[i].setValue(d2)                  except IllegalCoefficientValue,m:
1653          self.alteredCoefficient(i)                      raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
1654                          except IllegalCoefficientFunctionSpace,m:
1655        # reset the HomogeneousConstraintFlag:                      raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m))
1656        self.__setHomogeneousConstraintFlag()              else:
1657        if len(coefficients)>0 and not self.isUsingLumping() and not self.__homogeneous_constraint: self.__rebuildSystem()                  raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m))
1658            except IllegalCoefficientValue,m:
1659     def __setHomogeneousConstraintFlag(self):             raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
1660        """        self.__altered_coefficients=True
1661        checks if the constraints are homogeneous and sets self.__homogeneous_constraint accordingly.        # check if the systrem is inhomogeneous:
1662        """        if len(coefficients)>0 and not self.isUsingLumping():
1663        self.__homogeneous_constraint=True           q=self.getCoefficientOfGeneralPDE("q")
1664        q=self.getCoefficientOfPDE("q")           r=self.getCoefficientOfGeneralPDE("r")
1665        r=self.getCoefficientOfPDE("r")           homogeneous_constraint=True
1666        if not q.isEmpty() and not r.isEmpty():           if not q.isEmpty() and not r.isEmpty():
1667           if (q*r).Lsup()>=1.e-13*r.Lsup(): self.__homogeneous_constraint=False               if util.Lsup(q*r)>0.:
1668        if self.debug():                 self.trace("Inhomogeneous constraint detected.")
1669             if self.__homogeneous_constraint:                 self.__invalidateSystem()
                print "PDE Debug: Constraints are homogeneous."  
            else:  
                print "PDE Debug: Constraints are inhomogeneous."  
   
   
    # ==== rebuild switches =====================================================================  
    def __rebuildSolution(self,deep=False):  
        """  
        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()  
   
   
    def __rebuildOperator(self,deep=False):  
        """  
        indicates the operator has to be rebuilt next time it is used  
        """  
        if self.__operator_isValid and self.debug() : print "PDE Debug: Operator has to be rebuilt."  
        self.__rebuildSolution(deep)  
        self.__operator_isValid=False  
        if deep: self.__operator=escript.Operator()  
   
    def __rebuildRightHandSide(self,deep=False):  
        """  
        indicates the right hand side has to be rebuild next time it is used  
        """  
        if self.__righthandside_isValid and self.debug() : print "PDE Debug: Right hand side has to be rebuilt."  
        self.__rebuildSolution(deep)  
        self.__righthandside_isValid=False  
        if deep: self.__righthandside=escript.Data()  
   
    def __rebuildSystem(self,deep=False):  
        """  
        annonced that all coefficient name has been changed  
        """  
        self.__rebuildSolution(deep)  
        self.__rebuildOperator(deep)  
        self.__rebuildRightHandSide(deep)  
     
    def __checkMatrixType(self):  
      """  
      reassess the matrix type and, if needed, initiates an operator rebuild  
      """  
      new_matrix_type=self.getDomain().getSystemMatrixTypeId(self.getSolverMethod(),self.isSymmetric())  
      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)  
   
    #============ assembling =======================================================  
    def __copyConstraint(self):  
       """  
       copies the constrint condition into u  
       """  
       if not self.__righthandside.isEmpty():  
          q=self.getCoefficientOfPDE("q")  
          r=self.getCoefficientOfPDE("r")  
          if not q.isEmpty():  
              if r.isEmpty():  
                 r2=escript.Data(0,self.__righthandside.getShape(),self.__righthandside.getFunctionSpace())  
              else:  
                 r2=escript.Data(r,self.__righthandside.getFunctionSpace())  
              self.__righthandside.copyWithMask(r2,escript.Data(q,self.__righthandside.getFunctionSpace()))  
   
    def __applyConstraint(self):  
        """  
        applies the constraints defined by q and r to the system  
        """  
        q=self.getCoefficientOfPDE("q")  
        r=self.getCoefficientOfPDE("r")  
        if not q.isEmpty() and not self.__operator.isEmpty():  
           # q is the row and column mask to indicate where constraints are set:  
           row_q=escript.Data(q,self.getFunctionSpaceForEquation())  
           col_q=escript.Data(q,self.getFunctionSpaceForSolution())  
           u=self.__getNewSolution()  
           if r.isEmpty():  
              r_s=self.__getNewSolution()  
           else:  
              r_s=escript.Data(r,self.getFunctionSpaceForSolution())  
           u.copyWithMask(r_s,col_q)  
           if self.isUsingLumping():  
              self.__operator.copyWithMask(escript.Data(1,q.getShape(),self.getFunctionSpaceForEquation()),row_q)  
           else:  
              if not self.__righthandside.isEmpty(): self.__righthandside-=self.__operator*u  
              self.__operator.nullifyRowsAndCols(row_q,col_q,1.)  
1670    
1671     def getSystem(self):     def getSystem(self):
1672         """         """
1673         return the operator and right hand side of the PDE         return the operator and right hand side of the PDE
1674    
1675           @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_isValid or not self.__righthandside_isValid:         if not self.__operator_is_Valid or not self.__righthandside_isValid:
1679            if self.isUsingLumping():            if self.isUsingLumping():
1680                if not self.__operator_isValid:                if not self.__operator_is_Valid:
1681                   if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():                   if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():
1682                         raise TypeError,"Lumped matrix requires same order for equations and unknowns"                        raise TypeError,"Lumped matrix requires same order for equations and unknowns"
1683                   if not self.getCoefficientOfPDE("A").isEmpty():                   if not self.getCoefficientOfGeneralPDE("A").isEmpty():
1684                            raise Warning,"Lumped matrix does not allow coefficient A"                        raise ValueError,"coefficient A in lumped matrix may not be present."
1685                   if not self.getCoefficientOfPDE("B").isEmpty():                   if not self.getCoefficientOfGeneralPDE("B").isEmpty():
1686                            raise Warning,"Lumped matrix does not allow coefficient B"                        raise ValueError,"coefficient B in lumped matrix may not be present."
1687                   if not self.getCoefficientOfPDE("C").isEmpty():                   if not self.getCoefficientOfGeneralPDE("C").isEmpty():
1688                            raise Warning,"Lumped matrix does not allow coefficient C"                        raise ValueError,"coefficient C in lumped matrix may not be present."
1689                   if self.debug() : print "PDE Debug: New lumped operator is built."                   if not self.getCoefficientOfGeneralPDE("d_contact").isEmpty():
1690                   mat=self.__getNewOperator()                        raise ValueError,"coefficient d_contact in lumped matrix may not be present."
1691                   self.getDomain().addPDEToSystem(mat,escript.Data(), \                   if not self.getCoefficientOfGeneralPDE("A_reduced").isEmpty():
1692                             self.getCoefficientOfPDE("A"), \                        raise ValueError,"coefficient A_reduced in lumped matrix may not be present."
1693                             self.getCoefficientOfPDE("B"), \                   if not self.getCoefficientOfGeneralPDE("B_reduced").isEmpty():
1694                             self.getCoefficientOfPDE("C"), \                        raise ValueError,"coefficient B_reduced in lumped matrix may not be present."
1695                             self.getCoefficientOfPDE("D"), \                   if not self.getCoefficientOfGeneralPDE("C_reduced").isEmpty():
1696                             escript.Data(), \                        raise ValueError,"coefficient C_reduced in lumped matrix may not be present."
1697                             escript.Data(), \                   if not self.getCoefficientOfGeneralPDE("d_contact_reduced").isEmpty():
1698                             self.getCoefficientOfPDE("d"), \                        raise ValueError,"coefficient d_contact_reduced in lumped matrix may not be present."
1699                             escript.Data(),\                   D=self.getCoefficientOfGeneralPDE("D")
1700                             self.getCoefficientOfPDE("d_contact"), \                   d=self.getCoefficientOfGeneralPDE("d")
1701                             escript.Data())                   D_reduced=self.getCoefficientOfGeneralPDE("D_reduced")
1702                   self.__operator=mat*escript.Data(1,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)                   d_reduced=self.getCoefficientOfGeneralPDE("d_reduced")
1703                   self.__applyConstraint()                   if not D.isEmpty():
1704                   self.__operator_isValid=True                       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:                if not self.__righthandside_isValid:
                  if self.debug() : print "PDE Debug: New right hand side is built."  
1752                   self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \                   self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \
1753                                 self.getCoefficientOfPDE("X"), \                                 self.getCoefficientOfGeneralPDE("X"), \
1754                                 self.getCoefficientOfPDE("Y"),\                                 self.getCoefficientOfGeneralPDE("Y"),\
1755                                 self.getCoefficientOfPDE("y"),\                                 self.getCoefficientOfGeneralPDE("y"),\
1756                                 self.getCoefficientOfPDE("y_contact"))                                 self.getCoefficientOfGeneralPDE("y_contact"))
1757                   self.__copyConstraint()                   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                   self.__righthandside_isValid=True
1764            else:            else:
1765               if not self.__operator_isValid and not self.__righthandside_isValid:               if not self.__operator_is_Valid and not self.__righthandside_isValid:
                  if self.debug() : print "PDE Debug: New system is built."  
1766                   self.getDomain().addPDEToSystem(self.__makeFreshOperator(),self.__makeFreshRightHandSide(), \                   self.getDomain().addPDEToSystem(self.__makeFreshOperator(),self.__makeFreshRightHandSide(), \
1767                                 self.getCoefficientOfPDE("A"), \                                 self.getCoefficientOfGeneralPDE("A"), \
1768                                 self.getCoefficientOfPDE("B"), \                                 self.getCoefficientOfGeneralPDE("B"), \
1769                                 self.getCoefficientOfPDE("C"), \                                 self.getCoefficientOfGeneralPDE("C"), \
1770                                 self.getCoefficientOfPDE("D"), \                                 self.getCoefficientOfGeneralPDE("D"), \
1771                                 self.getCoefficientOfPDE("X"), \                                 self.getCoefficientOfGeneralPDE("X"), \
1772                                 self.getCoefficientOfPDE("Y"), \                                 self.getCoefficientOfGeneralPDE("Y"), \
1773                                 self.getCoefficientOfPDE("d"), \                                 self.getCoefficientOfGeneralPDE("d"), \
1774                                 self.getCoefficientOfPDE("y"), \                                 self.getCoefficientOfGeneralPDE("y"), \
1775                                 self.getCoefficientOfPDE("d_contact"), \                                 self.getCoefficientOfGeneralPDE("d_contact"), \
1776                                 self.getCoefficientOfPDE("y_contact"))                                 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()                   self.__applyConstraint()
1789                   self.__copyConstraint()                   self.__righthandside=self.copyConstraint(self.__righthandside)
1790                   self.__operator_isValid=True                   self.trace("New system has been built.")
1791                     self.__operator_is_Valid=True
1792                   self.__righthandside_isValid=True                   self.__righthandside_isValid=True
1793               elif not self.__righthandside_isValid:               elif not self.__righthandside_isValid:
                  if self.debug() : print "PDE Debug: New right hand side is built."  
1794                   self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \                   self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \
1795                                 self.getCoefficientOfPDE("X"), \                                 self.getCoefficientOfGeneralPDE("X"), \
1796                                 self.getCoefficientOfPDE("Y"),\                                 self.getCoefficientOfGeneralPDE("Y"),\
1797                                 self.getCoefficientOfPDE("y"),\                                 self.getCoefficientOfGeneralPDE("y"),\
1798                                 self.getCoefficientOfPDE("y_contact"))                                 self.getCoefficientOfGeneralPDE("y_contact"))
1799                   self.__copyConstraint()                   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                   self.__righthandside_isValid=True
1807               elif not self.__operator_isValid:               elif not self.__operator_is_Valid:
                  if self.debug() : print "PDE Debug: New operator is built."  
1808                   self.getDomain().addPDEToSystem(self.__makeFreshOperator(),escript.Data(), \                   self.getDomain().addPDEToSystem(self.__makeFreshOperator(),escript.Data(), \
1809                              self.getCoefficientOfPDE("A"), \                              self.getCoefficientOfGeneralPDE("A"), \
1810                              self.getCoefficientOfPDE("B"), \                              self.getCoefficientOfGeneralPDE("B"), \
1811                              self.getCoefficientOfPDE("C"), \                              self.getCoefficientOfGeneralPDE("C"), \
1812                              self.getCoefficientOfPDE("D"), \                              self.getCoefficientOfGeneralPDE("D"), \
1813                              escript.Data(), \                              escript.Data(), \
1814                              escript.Data(), \                              escript.Data(), \
1815                              self.getCoefficientOfPDE("d"), \                              self.getCoefficientOfGeneralPDE("d"), \
1816                              escript.Data(),\                              escript.Data(),\
1817                              self.getCoefficientOfPDE("d_contact"), \                              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())                              escript.Data())
1830                   self.__applyConstraint()                   self.__applyConstraint()
1831                   self.__operator_isValid=True                   self.trace("New operator has been built.")
1832                     self.__operator_is_Valid=True
1833         return (self.__operator,self.__righthandside)         return (self.__operator,self.__righthandside)
    def getOperator(self):  
        """  
        returns the operator of the PDE  
        """  
        return self.getSystem()[0]  
   
    def getRightHandSide(self):  
        """  
        returns the right hand side of the PDE  
        """  
        return self.getSystem()[1]  
   
    def solve(self,**options):  
       """  
       solve the PDE  
   
       @param options:  
       """  
       mat,f=self.getSystem()  
       if self.isUsingLumping():  
          out=f/mat  
       else:  
          options[util.TOLERANCE_KEY]=self.getTolerance()  
          options[util.METHOD_KEY]=self.getSolverMethod()  
          options[util.SYMMETRY_KEY]=self.isSymmetric()  
          if self.debug() : print "PDE Debug: solver options: ",options  
          out=mat.solve(f,options)  
       return out  
   
    def getSolution(self,**options):  
        """  
        returns the solution of the PDE  
   
        @param options:  
        """  
        if not self.__solution_isValid:  
            if self.debug() : print "PDE Debug: PDE is resolved."  
            self.__solution=self.solve(**options)  
            self.__solution_isValid=True  
        return self.__solution  
   
1834    
1835    
1836  def ELMAN_RAMAGE(P):  class Poisson(LinearPDE):
      """   """  
      return (P-1.).wherePositive()*0.5*(1.-1./(P+1.e-15))  
 def SIMPLIFIED_BROOK_HUGHES(P):  
      """   """  
      c=(P-3.).whereNegative()  
      return P/6.*c+1./2.*(1.-c)  
   
 def HALF(P):  
     """ """  
     return escript.Scalar(0.5,P.getFunctionSpace())  
   
 class AdvectivePDE(LinearPDE):  
1837     """     """
1838     Class to handle a linear PDE dominated by advective terms:     Class to define a Poisson equation problem, which is genear L{LinearPDE} of the form
     
    class to define a linear PDE of the form  
   
    \f[  
    -(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  
    \f]  
1839    
1840     with boundary conditons:     M{-grad(grad(u)[j])[j] = f}
1841    
1842     \f[     with natural boundary conditons
    n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_{ik}u_k = - n_j*X_{ij} + y_i  
    \f]  
1843    
1844     and contact conditions     M{n[j]*grad(u)[j] = 0 }
   
    \f[  
    n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d^{contact}_{ik}[u_k] = - n_j*X_{ij} + y^{contact}_{i}  
    \f]  
1845    
1846     and constraints:     and constraints:
1847    
1848     \f[     M{u=0} where M{q>0}
1849     u_i=r_i \quad \mathrm{where} \quad q_i>0  
    \f]  
1850     """     """
    def __init__(self,domain,numEquations=0,numSolutions=0,xi=ELMAN_RAMAGE):  
       LinearPDE.__init__(self,domain,numEquations,numSolutions)  
       self.__xi=xi  
       self.__Xi=escript.Data()  
   
    def __calculateXi(self,peclet_factor,Z,h):  
        Z_max=util.Lsup(Z)  
        if Z_max>0.:  
           return h*self.__xi(Z*peclet_factor)/(Z+Z_max*self.TOL)  
        else:  
           return 0.  
1851    
1852     def setValue(self,**args):     def __init__(self,domain,debug=False):
1853         if "A" in args.keys()   or "B" in args.keys() or "C" in args.keys(): self.__Xi=escript.Data()       """
1854         self._LinearPDE__setValue(**args)       initializes a new Poisson equation
             
    def getXi(self):  
       if self.__Xi.isEmpty():  
          B=self.getCoefficient("B")  
          C=self.getCoefficient("C")  
          A=self.getCoefficient("A")  
          h=self.getDomain().getSize()  
          self.__Xi=escript.Scalar(0.,self.getFunctionSpaceForCoefficient("A"))  
          if not C.isEmpty() or not B.isEmpty():  
             if not C.isEmpty() and not B.isEmpty():  
                 Z2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))  
                 if self.getNumEquations()>1:  
                    if self.getNumSolutions()>1:  
                       for i in range(self.getNumEquations()):  
                          for k in range(self.getNumSolutions()):  
                             for l in range(self.getDim()): Z2+=(C[i,k,l]-B[i,l,k])**2  
                    else:  
                       for i in range(self.getNumEquations()):  
                          for l in range(self.getDim()): Z2+=(C[i,l]-B[i,l])**2  
                 else:  
                    if self.getNumSolutions()>1:  
                       for k in range(self.getNumSolutions()):  
                          for l in range(self.getDim()): Z2+=(C[k,l]-B[l,k])**2  
                    else:  
                       for l in range(self.getDim()): Z2+=(C[l]-B[l])**2  
                 length_of_Z=util.sqrt(Z2)  
             elif C.isEmpty():  
               length_of_Z=util.length(B)  
             else:  
               length_of_Z=util.length(C)  
1855    
1856              Z_max=util.Lsup(length_of_Z)       @param domain: domain of the PDE
1857              if Z_max>0.:       @type domain: L{Domain<escript.Domain>}
1858                 length_of_A=util.length(A)       @param debug: if True debug informations are printed.
                A_max=util.Lsup(length_of_A)  
                if A_max>0:  
                     inv_A=1./(length_of_A+A_max*self.TOL)  
                else:  
                     inv_A=1./self.TOL  
                peclet_number=length_of_Z*h/2*inv_A  
                xi=self.__xi(peclet_number)  
                self.__Xi=h*xi/(length_of_Z+Z_max*self.TOL)  
                print "@ preclet number = %e"%util.Lsup(peclet_number),util.Lsup(xi),util.Lsup(length_of_Z)  
       return self.__Xi  
         
1859    
    def getCoefficientOfPDE(self,name):  
1860       """       """
1861       return the value of the coefficient name of the general PDE       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       @param name:     def setValue(self,**coefficients):
1868       """       """
1869       if not self.getNumEquations() == self.getNumSolutions():       sets new values to coefficients
           raise ValueError,"AdvectivePDE expects the number of solution componets and the number of equations to be equal."  
1870    
1871       if name == "A" :       @param coefficients: new values assigned to coefficients
1872           A=self.getCoefficient("A")       @keyword f: value for right hand side M{f}
1873           B=self.getCoefficient("B")       @type f: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1874           C=self.getCoefficient("C")       @keyword q: mask for location of constraints
1875           if B.isEmpty() and C.isEmpty():       @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              Aout=A                 depending of reduced order is used for the representation of the equation.
1877           else:       @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1878              if A.isEmpty():       """
1879                 Aout=self.createNewCoefficient("A")       super(Poisson, self).setValue(**coefficients)
1880              else:  
1881                 Aout=A[:]     def getCoefficientOfGeneralPDE(self,name):
1882              Xi=self.getXi()       """
1883              if self.getNumEquations()>1:       return the value of the coefficient name of the general PDE
1884                  for i in range(self.getNumEquations()):       @param name: name of the coefficient requested.
1885                     for j in range(self.getDim()):       @type name: C{string}
1886                        for k in range(self.getNumSolutions()):       @return: the value of the coefficient  name
1887                           for l in range(self.getDim()):       @rtype: L{Data<escript.Data>}
1888                              if not C.isEmpty() and not B.isEmpty():       @raise IllegalCoefficient: if name is not one of coefficients
1889                                 for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*(C[p,i,j]-B[p,j,i])*(C[p,k,l]-B[p,l,k])                    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                              elif C.isEmpty():       @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
1891                                 for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*B[p,j,i]*B[p,l,k]       """
1892                              else:       if name == "A" :
1893                                 for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*C[p,i,j]*C[p,k,l]           return escript.Data(util.kronecker(self.getDim()),escript.Function(self.getDomain()))
1894              else:       elif name == "B" :
1895                  for j in range(self.getDim()):           return escript.Data()
1896                     for l in range(self.getDim()):       elif name == "C" :
1897                        if not C.isEmpty() and not B.isEmpty():           return escript.Data()
1898                            Aout[j,l]+=Xi*(C[j]-B[j])*(C[l]-B[l])       elif name == "D" :
1899                        elif C.isEmpty():           return escript.Data()
1900                            Aout[j,l]+=Xi*B[j]*B[l]       elif name == "X" :
1901                        else:           return escript.Data()
1902                            Aout[j,l]+=Xi*C[j]*C[l]       elif name == "Y" :
1903           return Aout           return self.getCoefficient("f")
1904       elif name == "B" :       elif name == "d" :
1905           B=self.getCoefficient("B")           return escript.Data()
1906           C=self.getCoefficient("C")       elif name == "y" :
1907           D=self.getCoefficient("D")           return escript.Data()
1908           if C.isEmpty() or D.isEmpty():       elif name == "d_contact" :
1909              Bout=B           return escript.Data()
          else:  
             Xi=self.getXi()  
             if B.isEmpty():  
                 Bout=self.createNewCoefficient("B")  
             else:  
                 Bout=B[:]  
             if self.getNumEquations()>1:  
                for k in range(self.getNumSolutions()):  
                   for p in range(self.getNumEquations()):  
                      tmp=Xi*D[p,k]  
                      for i in range(self.getNumEquations()):  
                         for j in range(self.getDim()):  
                            Bout[i,j,k]+=tmp*C[p,i,j]  
             else:  
                tmp=Xi*D  
                for j in range(self.getDim()): Bout[j]+=tmp*C[j]  
          return Bout  
      elif name == "C" :  
          B=self.getCoefficient("B")  
          C=self.getCoefficient("C")  
          D=self.getCoefficient("D")  
          if B.isEmpty() or D.isEmpty():  
             Cout=C  
          else:  
             Xi=self.getXi()  
             if C.isEmpty():  
                 Cout=self.createNewCoefficient("C")  
             else:  
                 Cout=C[:]  
             if self.getNumEquations()>1:  
                for k in range(self.getNumSolutions()):  
                    for p in range(self.getNumEquations()):  
                       tmp=Xi*D[p,k]  
                       for i in range(self.getNumEquations()):  
                         for l in range(self.getDim()):  
                                  Cout[i,k,l]+=tmp*B[p,l,i]  
             else:  
                tmp=Xi*D  
                for j in range(self.getDim()): Cout[j]+=tmp*B[j]  
          return Cout  
      elif name == "D" :  
          return self.getCoefficient("D")  
      elif name == "X" :  
          X=self.getCoefficient("X")  
          Y=self.getCoefficient("Y")  
          B=self.getCoefficient("B")  
          C=self.getCoefficient("C")  
          if Y.isEmpty() or (B.isEmpty() and C.isEmpty()):  
             Xout=X  
          else:  
             if X.isEmpty():  
                 Xout=self.createNewCoefficient("X")  
             else:  
                 Xout=X[:]  
             Xi=self.getXi()  
             if self.getNumEquations()>1:  
                  for p in range(self.getNumEquations()):  
                     tmp=Xi*Y[p]  
                     for i in range(self.getNumEquations()):  
                        for j in range(self.getDim()):  
                           if not C.isEmpty() and not B.isEmpty():  
                              Xout[i,j]+=tmp*(C[p,i,j]-B[p,j,i])  
                           elif C.isEmpty():  
                              Xout[i,j]-=tmp*B[p,j,i]  
                           else:  
                              Xout[i,j]+=tmp*C[p,i,j]  
             else:  
                  tmp=Xi*Y  
                  for j in range(self.getDim()):  
                     if not C.isEmpty() and not B.isEmpty():  
                        Xout[j]+=tmp*(C[j]-B[j])  
                     elif C.isEmpty():  
                        Xout[j]-=tmp*B[j]  
                     else:  
                        Xout[j]+=tmp*C[j]  
          return Xout  
      elif name == "Y" :  
          return self.getCoefficient("Y")  
      elif name == "d" :  
          return self.getCoefficient("d")  
      elif name == "y" :  
          return self.getCoefficient("y")  
      elif name == "d_contact" :  
          return self.getCoefficient("d_contact")  
1910       elif name == "y_contact" :       elif name == "y_contact" :
1911           return self.getCoefficient("y_contact")           return escript.Data()
1912       elif name == "r" :       elif name == "A_reduced" :
1913           return self.getCoefficient("r")           return escript.Data()
1914       elif name == "q" :       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")           return self.getCoefficient("q")
1936       else:       else:
1937           raise SystemError,"unknown PDE coefficient %s",name          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
   
1938    
1939  class Poisson(LinearPDE):  class Helmholtz(LinearPDE):
1940     """     """
1941     Class to define a Poisson equation problem:     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     class to define a linear PDE of the form     with natural boundary conditons
1946     \f[  
1947     -u_{,jj} = f     M{k*n[j]*grad(u)[j] = g- S{alpha}u }
    \f]  
   
    with boundary conditons:  
   
    \f[  
    n_j*u_{,j} = 0  
    \f]  
1948    
1949     and constraints:     and constraints:
1950    
1951     \f[     M{u=r} where M{q>0}
1952     u=0 \quad \mathrm{where} \quad q>0  
    \f]  
1953     """     """
1954    
1955     def __init__(self,domain,f=escript.Data(),q=escript.Data()):     def __init__(self,domain,debug=False):
1956         LinearPDE.__init__(self,domain,1,1)       """
1957         self.COEFFICIENTS={       initializes a new Poisson equation
1958         "f"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),  
1959         "q"         : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.BOTH)}       @param domain: domain of the PDE
1960         self.setSymmetryOn()       @type domain: L{Domain<escript.Domain>}
1961         self.setValue(f,q)       @param debug: if True debug informations are printed.
1962    
1963     def setValue(self,f=escript.Data(),q=escript.Data()):       """
1964         """set value of PDE parameters f and q"""       super(Helmholtz, self).__init__(domain,1,1,debug)
1965         self._LinearPDE__setValue(f=f,q=q)       self.COEFFICIENTS={"omega": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
1966                            "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                            "g": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1971                            "g_reduced": PDECoefficient(PDECoefficient.BOUNDARY_REDUCED,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1972                            "r": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH),
1973                            "q": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
1974         self.setSymmetryOn()
1975    
1976       def setValue(self,**coefficients):
1977         """
1978         sets new values to coefficients
1979    
1980     def getCoefficientOfPDE(self,name):       @param coefficients: new values assigned to coefficients
1981         @keyword omega: value for coefficient M{S{omega}}
1982         @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         @type k: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1985         @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         @keyword alpha: value for right hand side M{S{alpha}}
1988         @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         @type g: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1991         @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                   depending of reduced order is used for the representation of the equation.
1994         @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                   depending of reduced order is used for the representation of the equation.
1997         @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1998       """       """
1999       return the value of the coefficient name of the general PDE       super(Helmholtz, self).setValue(**coefficients)
2000    
2001       @param name:     def getCoefficientOfGeneralPDE(self,name):
2002       """       """
2003       if name == "A" :       return the value of the coefficient name of the general PDE
2004           return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))  
2005       elif name == "B" :       @param name: name of the coefficient requested.
2006           return escript.Data()       @type name: C{string}
2007       elif name == "C" :       @return: the value of the coefficient  name
2008         @rtype: L{Data<escript.Data>}
2009         @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         @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
2012         """
2013         if name == "A" :
2014             return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))*self.getCoefficient("k")
2015         elif name == "B" :
2016           return escript.Data()           return escript.Data()
2017       elif name == "D" :       elif name == "C" :
2018           return escript.Data()           return escript.Data()
2019       elif name == "X" :       elif name == "D" :
2020             return self.getCoefficient("omega")
2021         elif name == "X" :
2022           return escript.Data()           return escript.Data()
2023       elif name == "Y" :       elif name == "Y" :
2024           return self.getCoefficient("f")           return self.getCoefficient("f")
2025       elif name == "d" :       elif name == "d" :
2026             return self.getCoefficient("alpha")
2027         elif name == "y" :
2028             return self.getCoefficient("g")
2029         elif name == "d_contact" :
2030           return escript.Data()           return escript.Data()
2031       elif name == "y" :       elif name == "y_contact" :
2032           return escript.Data()           return escript.Data()
2033       elif name == "d_contact" :       elif name == "A_reduced" :
2034           return escript.Data()           return escript.Data()
2035       elif name == "y_contact" :       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()           return escript.Data()
2041       elif name == "r" :       elif name == "X_reduced" :
2042           return escript.Data()           return escript.Data()
2043       elif name == "q" :       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")           return self.getCoefficient("q")
2057       else:       else:
2058           raise SystemError,"unknown PDE coefficient %s",name          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
2059    
2060  class LameEquation(LinearPDE):  class LameEquation(LinearPDE):
2061     """     """
2062     Class to define a Lame equation problem:     Class to define a Lame equation problem:
2063    
2064     class to define a linear PDE of the form     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     \f[  
2066     -(\mu (u_{i,j}+u_{j,i}))_{,j} - \lambda u_{j,ji}} = F_i -\sigma_{ij,j}     with natural boundary conditons:
2067     \f]  
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] }
    with boundary conditons:  
   
    \f[  
    n_j(\mu (u_{i,j}+u_{j,i})-sigma_{ij}) + n_i\lambda u_{j,j} = f_i  
    \f]  
2069    
2070     and constraints:     and constraints:
2071    
2072     \f[     M{u[i]=r[i]} where M{q[i]>0}
2073     u_i=r_i \quad \mathrm{where} \quad q_i>0  
    \f]  
2074     """     """
2075    
2076     def __init__(self,domain,f=escript.Data(),q=escript.Data()):     def __init__(self,domain,debug=False):
2077         LinearPDE.__init__(self,domain,domain.getDim(),domain.getDim())        super(LameEquation, self).__init__(domain,\
2078         self.COEFFICIENTS={                                           domain.getDim(),domain.getDim(),debug)
2079         "lame_lambda"  : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),        self.COEFFICIENTS={ "lame_lambda"  : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),
2080         "lame_mu"      : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),                            "lame_mu"      : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),
2081         "F"            : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),                            "F"            : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
2082         "sigma"        : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM),PDECoefficient.RIGHTHANDSIDE),                            "sigma"        : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM),PDECoefficient.RIGHTHANDSIDE),
2083         "f"            : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),                            "f"            : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
2084         "r"            : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.BOTH),                            "r"            : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH),
2085         "q"            : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.BOTH)}                            "q"            : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
2086         self.setSymmetryOn()        self.setSymmetryOn()
2087    
2088     def setValue(self,lame_lambda=escript.Data(),lame_mu=escript.Data(),F=escript.Data(),sigma=escript.Data(),f=escript.Data(),r=escript.Data(),q=escript.Data()):     def setValues(self,**coefficients):
2089         """set value of PDE parameters"""       """
2090         self._LinearPDE__setValue(lame_lambda=lame_lambda, \       sets new values to coefficients
2091                                   lame_mu=lame_mu, \  
2092                                   F=F, \       @param coefficients: new values assigned to coefficients
2093                                   sigma=sigma, \       @keyword lame_mu: value for coefficient M{S{mu}}
2094                                   f=f, \       @type lame_mu: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2095                                   r=r, \       @keyword lame_lambda: value for coefficient M{S{lambda}}
2096                                   q=q)       @type lame_lambda: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2097     def getCoefficientOfPDE(self,name):       @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       return the value of the coefficient name of the general PDE
2116    
2117       @param name:       @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" :       if name == "A" :
2126           out =self.createNewCoefficientOfPDE("A")           out =self.createCoefficientOfGeneralPDE("A")
2127           for i in range(self.getDim()):           for i in range(self.getDim()):
2128             for j in range(self.getDim()):             for j in range(self.getDim()):
2129               out[i,i,j,j] += self.getCoefficient("lame_lambda")               out[i,i,j,j] += self.getCoefficient("lame_lambda")
2130               out[i,j,j,i] += self.getCoefficient("lame_mu")               out[i,j,j,i] += self.getCoefficient("lame_mu")
2131               out[i,j,i,j] += self.getCoefficient("lame_mu")               out[i,j,i,j] += self.getCoefficient("lame_mu")
2132           return out           return out
2133       elif name == "B" :       elif name == "B" :
2134           return escript.Data()           return escript.Data()
2135       elif name == "C" :       elif name == "C" :
2136           return escript.Data()           return escript.Data()
2137       elif name == "D" :       elif name == "D" :
2138           return escript.Data()           return escript.Data()
2139       elif name == "X" :       elif name == "X" :
2140           return self.getCoefficient("sigma")           return self.getCoefficient("sigma")
2141       elif name == "Y" :       elif name == "Y" :
2142           return self.getCoefficient("F")           return self.getCoefficient("F")
2143       elif name == "d" :       elif name == "d" :
2144           return escript.Data()           return escript.Data()
2145       elif name == "y" :       elif name == "y" :
2146           return self.getCoefficient("f")           return self.getCoefficient("f")
2147       elif name == "d_contact" :       elif name == "d_contact" :
2148           return escript.Data()           return escript.Data()
2149       elif name == "y_contact" :       elif name == "y_contact" :
2150           return escript.Data()           return escript.Data()
2151       elif name == "r" :       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")           return self.getCoefficient("r")
2173       elif name == "q" :       elif name == "q" :
2174           return self.getCoefficient("q")           return self.getCoefficient("q")
2175       else:       else:
2176           raise SystemError,"unknown PDE coefficient %s",name          raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
2177    
 # $Log$  
 # Revision 1.10  2005/08/12 01:45:36  jgs  
 # erge of development branch dev-02 back to main trunk on 2005-08-12  
 #  
 # Revision 1.9.2.1  2005/07/29 07:10:27  gross  
 # new functions in util and a new pde type in linearPDEs  
 #  
 # Revision 1.1.2.25  2005/07/28 04:21:09  gross  
 # Lame equation: (linear elastic, isotropic) added  
 #  
 # Revision 1.1.2.24  2005/07/22 06:37:11  gross  
 # some extensions to modellib and linearPDEs  
 #  
 # Revision 1.1.2.23  2005/05/13 00:55:20  cochrane  
 # Fixed up some docstrings.  Moved module-level functions to top of file so  
 # that epydoc and doxygen can pick them up properly.  
 #  
 # Revision 1.1.2.22  2005/05/12 11:41:30  gross  
 # some basic Models have been added  
 #  
 # Revision 1.1.2.21  2005/05/12 07:16:12  cochrane  
 # Moved ELMAN_RAMAGE, SIMPLIFIED_BROOK_HUGHES, and HALF functions to bottom of  
 # file so that the AdvectivePDE class is picked up by doxygen.  Some  
 # reformatting of docstrings.  Addition of code to make equations come out  
 # as proper LaTeX.  
 #  
 # Revision 1.1.2.20  2005/04/15 07:09:08  gross  
 # some problems with functionspace and linearPDEs fixed.  
 #  
 # Revision 1.1.2.19  2005/03/04 05:27:07  gross  
 # bug in SystemPattern fixed.  
 #  
 # Revision 1.1.2.18  2005/02/08 06:16:45  gross  
 # Bugs in AdvectivePDE fixed, AdvectiveTest is stable but more testing is needed  
 #  
 # Revision 1.1.2.17  2005/02/08 05:56:19  gross  
 # Reference Number handling added  
 #  
 # Revision 1.1.2.16  2005/02/07 04:41:28  gross  
 # some function exposed to python to make mesh merging running  
 #  
 # Revision 1.1.2.15  2005/02/03 00:14:44  gross  
 # timeseries add and ESySParameter.py renames esysXML.py for consistence  
 #  
 # Revision 1.1.2.14  2005/02/01 06:44:10  gross  
 # new implementation of AdvectivePDE which now also updates right hand side. systems of PDEs are still not working  
 #  
 # Revision 1.1.2.13  2005/01/25 00:47:07  gross  
 # updates in the documentation  
 #  
 # Revision 1.1.2.12  2005/01/12 01:28:04  matt  
 # Added createCoefficient method for linearPDEs.  
 #  
 # Revision 1.1.2.11  2005/01/11 01:55:34  gross  
 # a problem in linearPDE class fixed  
 #  
 # Revision 1.1.2.10  2005/01/07 01:13:29  gross  
 # some bugs in linearPDE fixed  
 #  
 # Revision 1.1.2.9  2005/01/06 06:24:58  gross  
 # some bugs in slicing fixed  
 #  
 # Revision 1.1.2.8  2005/01/05 04:21:40  gross  
 # FunctionSpace checking/matchig in slicing added  
 #  
 # Revision 1.1.2.7  2004/12/29 10:03:41  gross  
 # bug in setValue fixed  
 #  
 # Revision 1.1.2.6  2004/12/29 05:29:59  gross  
 # AdvectivePDE successfully tested for Peclet number 1000000. there is still a problem with setValue and Data()  
 #  
 # Revision 1.1.2.5  2004/12/29 00:18:41  gross  
 # AdvectivePDE added  
 #  
 # Revision 1.1.2.4  2004/12/24 06:05:41  gross  
 # some changes in linearPDEs to add AdevectivePDE  
 #  
 # Revision 1.1.2.3  2004/12/16 00:12:34  gross  
 # __init__ of LinearPDE does not accept any coefficient anymore  
 #  
 # Revision 1.1.2.2  2004/12/14 03:55:01  jgs  
 # *** empty log message ***  
 #  
 # Revision 1.1.2.1  2004/12/12 22:53:47  gross  
 # linearPDE has been renamed LinearPDE  
 #  
 # Revision 1.1.1.1.2.7  2004/12/07 10:13:08  gross  
 # GMRES added  
 #  
 # Revision 1.1.1.1.2.6  2004/12/07 03:19:50  gross  
 # options for GMRES and PRES20 added  
 #  
 # Revision 1.1.1.1.2.5  2004/12/01 06:25:15  gross  
 # some small changes  
 #  
 # Revision 1.1.1.1.2.4  2004/11/24 01:50:21  gross  
 # Finley solves 4M unknowns now  
 #  
 # Revision 1.1.1.1.2.3  2004/11/15 06:05:26  gross  
 # poisson solver added  
 #  
 # Revision 1.1.1.1.2.2  2004/11/12 06:58:15  gross  
 # a lot of changes to get the linearPDE class running: most important change is that there is no matrix format exposed to the user anymore. the format is chosen by the Domain according to the solver and symmetry  
 #  
 # Revision 1.1.1.1.2.1  2004/10/28 22:59:22  gross  
 # finley's RecTest.py is running now: problem in SystemMatrixAdapater fixed  
 #  
 # Revision 1.1.1.1  2004/10/26 06:53:56  jgs  
 # initial import of project esys2  
 #  
 # Revision 1.3.2.3  2004/10/26 06:43:48  jgs  
 # committing Lutz's and Paul's changes to brach jgs  
 #  
 # Revision 1.3.4.1  2004/10/20 05:32:51  cochrane  
 # Added incomplete Doxygen comments to files, or merely put the docstrings that already exist into Doxygen form.  
 #  
 # Revision 1.3  2004/09/23 00:53:23  jgs  
 # minor fixes  
 #  
 # Revision 1.1  2004/08/28 12:58:06  gross  
 # SimpleSolve is not running yet: problem with == of functionsspace  
 #  
 #  

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

  ViewVC Help
Powered by ViewVC 1.1.26