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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 425 - (hide annotations)
Tue Jan 10 04:10:39 2006 UTC (13 years, 9 months ago) by gross
File MIME type: text/x-python
File size: 110594 byte(s)
The sparse solver can be called by paso now. 

the building has been change to reduce some code redundancy:
now all scons SCscripts are importing scons/esys_options.py which
imports platform specific settings. 



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