/[escript]/trunk/escript/py_src/linearPDEs.py
ViewVC logotype

Diff of /trunk/escript/py_src/linearPDEs.py

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

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

Legend:
Removed from v.122  
changed lines
  Added in v.149

  ViewVC Help
Powered by ViewVC 1.1.26