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

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

  ViewVC Help
Powered by ViewVC 1.1.26