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

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

  ViewVC Help
Powered by ViewVC 1.1.26