/[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 345 - (hide annotations)
Tue Dec 13 05:23:45 2005 UTC (13 years, 11 months ago) by gross
File MIME type: text/x-python
File size: 109712 byte(s)
advection diffusion class added
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 jgs 102
453 jgs 148
454     def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):
455 jgs 102 """
456 jgs 148 initializes a new linear PDE
457 jgs 102
458 jgs 148 @param domain: domain of the PDE
459 jgs 149 @type domain: L{Domain<escript.Domain>}
460 jgs 148 @param numEquations: number of equations. If numEquations==None the number of equations
461     is exracted from the PDE coefficients.
462     @param numSolutions: number of solution components. If numSolutions==None the number of solution components
463     is exracted from the PDE coefficients.
464     @param debug: if True debug informations are printed.
465    
466 jgs 102 """
467 jgs 151 super(LinearPDE, self).__init__()
468 jgs 148 #
469     # the coefficients of the general PDE:
470     #
471     self.__COEFFICIENTS_OF_GENEARL_PDE={
472 jgs 149 "A" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM,PDECoefficient.BY_SOLUTION,PDECoefficient.BY_DIM),PDECoefficient.OPERATOR),
473     "B" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
474     "C" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION,PDECoefficient.BY_DIM),PDECoefficient.OPERATOR),
475     "D" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
476     "X" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM),PDECoefficient.RIGHTHANDSIDE),
477     "Y" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
478     "d" : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
479     "y" : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
480     "d_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
481     "y_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
482     "r" : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_SOLUTION,),PDECoefficient.RIGHTHANDSIDE),
483     "q" : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_SOLUTION,),PDECoefficient.BOTH)}
484 jgs 102
485 jgs 148 # COEFFICIENTS can be overwritten by subclasses:
486     self.COEFFICIENTS=self.__COEFFICIENTS_OF_GENEARL_PDE
487 jgs 149 self.__altered_coefficients=False
488 jgs 102 # initialize attributes
489 jgs 148 self.__debug=debug
490 jgs 104 self.__domain=domain
491     self.__numEquations=numEquations
492     self.__numSolutions=numSolutions
493 jgs 148 self.__resetSystem()
494 jgs 102
495     # set some default values:
496 jgs 149 self.__reduce_equation_order=False
497     self.__reduce_solution_order=False
498 jgs 102 self.__tolerance=1.e-8
499 jgs 150 self.__solver_method=self.DEFAULT
500     self.__solver_package=self.DEFAULT
501     self.__matrix_type=self.__domain.getSystemMatrixTypeId(self.DEFAULT,self.DEFAULT,False)
502 jgs 102 self.__sym=False
503    
504 jgs 148 self.resetCoefficients()
505     self.trace("PDE Coeffients are %s"%str(self.COEFFICIENTS.keys()))
506     # =============================================================================
507     # general stuff:
508     # =============================================================================
509     def __str__(self):
510 jgs 149 """
511     returns string representation of the PDE
512    
513     @return: a simple representation of the PDE
514     @rtype: C{str}
515     """
516     return "<LinearPDE %d>"%id(self)
517 jgs 148 # =============================================================================
518     # debug :
519     # =============================================================================
520     def setDebugOn(self):
521 jgs 149 """
522 jgs 148 switches on debugging
523 jgs 108 """
524 jgs 148 self.__debug=not None
525    
526     def setDebugOff(self):
527 jgs 108 """
528 jgs 148 switches off debugging
529     """
530     self.__debug=None
531 jgs 108
532 jgs 148 def trace(self,text):
533     """
534 jgs 149 print the text message if debugging is swiched on.
535     @param text: message
536     @type text: C{string}
537 jgs 102 """
538 jgs 148 if self.__debug: print "%s: %s"%(str(self),text)
539 jgs 102
540 jgs 148 # =============================================================================
541     # some service functions:
542     # =============================================================================
543     def getDomain(self):
544 jgs 102 """
545 jgs 148 returns the domain of the PDE
546 jgs 102
547 jgs 149 @return: the domain of the PDE
548     @rtype: L{Domain<escript.Domain>}
549 jgs 108 """
550 jgs 148 return self.__domain
551 jgs 122
552 jgs 148 def getDim(self):
553 jgs 108 """
554 jgs 148 returns the spatial dimension of the PDE
555 jgs 108
556 jgs 149 @return: the spatial dimension of the PDE domain
557     @rtype: C{int}
558 jgs 148 """
559     return self.getDomain().getDim()
560 jgs 102
561 jgs 148 def getNumEquations(self):
562     """
563     returns the number of equations
564 jgs 102
565 jgs 149 @return: the number of equations
566     @rtype: C{int}
567 jgs 148 @raise UndefinedPDEError: if the number of equations is not be specified yet.
568     """
569     if self.__numEquations==None:
570     raise UndefinedPDEError,"Number of equations is undefined. Please specify argument numEquations."
571     else:
572     return self.__numEquations
573 jgs 147
574 jgs 148 def getNumSolutions(self):
575     """
576     returns the number of unknowns
577 jgs 147
578 jgs 149 @return: the number of unknowns
579     @rtype: C{int}
580 jgs 148 @raise UndefinedPDEError: if the number of unknowns is not be specified yet.
581     """
582     if self.__numSolutions==None:
583     raise UndefinedPDEError,"Number of solution is undefined. Please specify argument numSolutions."
584     else:
585     return self.__numSolutions
586    
587 jgs 149 def reduceEquationOrder(self):
588     """
589     return status for order reduction for equation
590    
591     @return: return True is reduced interpolation order is used for the represenation of the equation
592     @rtype: L{bool}
593     """
594     return self.__reduce_equation_order
595    
596     def reduceSolutionOrder(self):
597     """
598     return status for order reduction for the solution
599    
600     @return: return True is reduced interpolation order is used for the represenation of the solution
601     @rtype: L{bool}
602     """
603     return self.__reduce_solution_order
604 jgs 151
605 jgs 108 def getFunctionSpaceForEquation(self):
606     """
607 jgs 149 returns the L{FunctionSpace<escript.FunctionSpace>} used to discretize the equation
608 jgs 148
609 jgs 149 @return: representation space of equation
610     @rtype: L{FunctionSpace<escript.FunctionSpace>}
611 jgs 108 """
612 jgs 149 if self.reduceEquationOrder():
613     return escript.ReducedSolution(self.getDomain())
614     else:
615     return escript.Solution(self.getDomain())
616 jgs 108
617     def getFunctionSpaceForSolution(self):
618     """
619 jgs 149 returns the L{FunctionSpace<escript.FunctionSpace>} used to represent the solution
620 jgs 148
621 jgs 149 @return: representation space of solution
622     @rtype: L{FunctionSpace<escript.FunctionSpace>}
623 jgs 108 """
624 jgs 149 if self.reduceSolutionOrder():
625     return escript.ReducedSolution(self.getDomain())
626     else:
627     return escript.Solution(self.getDomain())
628 jgs 108
629 jgs 102
630 jgs 148 def getOperator(self):
631 jgs 102 """
632 jgs 148 provides access to the operator of the PDE
633 jgs 102
634 jgs 149 @return: the operator of the PDE
635     @rtype: L{Operator<escript.Operator>}
636 jgs 108 """
637 jgs 148 m=self.getSystem()[0]
638     if self.isUsingLumping():
639     return self.copyConstraint(1./m)
640 jgs 147 else:
641 jgs 148 return m
642 jgs 147
643 jgs 148 def getRightHandSide(self):
644 jgs 147 """
645 jgs 148 provides access to the right hand side of the PDE
646 jgs 149 @return: the right hand side of the PDE
647     @rtype: L{Data<escript.Data>}
648 jgs 147 """
649 jgs 148 r=self.getSystem()[1]
650     if self.isUsingLumping():
651     return self.copyConstraint(r)
652 jgs 147 else:
653 jgs 148 return r
654 jgs 147
655 jgs 148 def applyOperator(self,u=None):
656 jgs 102 """
657 jgs 148 applies the operator of the PDE to a given u or the solution of PDE if u is not present.
658 jgs 102
659 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}
660     the current solution is used.
661 jgs 149 @type u: L{Data<escript.Data>} or None
662     @return: image of u
663     @rtype: L{Data<escript.Data>}
664 jgs 102 """
665 jgs 149 if u==None:
666 jgs 148 return self.getOperator()*self.getSolution()
667 jgs 102 else:
668 jgs 148 self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())
669 jgs 102
670 jgs 148 def getResidual(self,u=None):
671 jgs 102 """
672 jgs 148 return the residual of u or the current solution if u is not present.
673 jgs 102
674 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}
675     the current solution is used.
676 jgs 149 @type u: L{Data<escript.Data>} or None
677     @return: residual of u
678     @rtype: L{Data<escript.Data>}
679 jgs 102 """
680 jgs 148 return self.applyOperator(u)-self.getRightHandSide()
681 jgs 102
682 jgs 148 def checkSymmetry(self,verbose=True):
683     """
684     test the PDE for symmetry.
685 jgs 102
686 jgs 149 @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.
687     @type verbose: C{bool}
688     @return: True if the PDE is symmetric.
689     @rtype: L{Data<escript.Data>}
690 jgs 148 @note: This is a very expensive operation. It should be used for degugging only! The symmetry flag is not altered.
691     """
692 jgs 149 verbose=verbose or self.__debug
693 jgs 148 out=True
694     if self.getNumSolutions()!=self.getNumEquations():
695     if verbose: print "non-symmetric PDE because of different number of equations and solutions"
696     out=False
697     else:
698     A=self.getCoefficientOfGeneralPDE("A")
699     if not A.isEmpty():
700 jgs 149 tol=util.Lsup(A)*self.__TOL
701 jgs 148 if self.getNumSolutions()>1:
702     for i in range(self.getNumEquations()):
703     for j in range(self.getDim()):
704     for k in range(self.getNumSolutions()):
705     for l in range(self.getDim()):
706     if util.Lsup(A[i,j,k,l]-A[k,l,i,j])>tol:
707     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)
708     out=False
709     else:
710     for j in range(self.getDim()):
711     for l in range(self.getDim()):
712     if util.Lsup(A[j,l]-A[l,j])>tol:
713     if verbose: print "non-symmetric PDE because A[%d,%d]!=A[%d,%d]"%(j,l,l,j)
714     out=False
715     B=self.getCoefficientOfGeneralPDE("B")
716     C=self.getCoefficientOfGeneralPDE("C")
717     if B.isEmpty() and not C.isEmpty():
718     if verbose: print "non-symmetric PDE because B is not present but C is"
719     out=False
720     elif not B.isEmpty() and C.isEmpty():
721     if verbose: print "non-symmetric PDE because C is not present but B is"
722     out=False
723     elif not B.isEmpty() and not C.isEmpty():
724 jgs 149 tol=(util.Lsup(B)+util.Lsup(C))*self.__TOL/2.
725 jgs 148 if self.getNumSolutions()>1:
726     for i in range(self.getNumEquations()):
727     for j in range(self.getDim()):
728     for k in range(self.getNumSolutions()):
729     if util.Lsup(B[i,j,k]-C[k,i,j])>tol:
730     if verbose: print "non-symmetric PDE because B[%d,%d,%d]!=C[%d,%d,%d]"%(i,j,k,k,i,j)
731     out=False
732     else:
733     for j in range(self.getDim()):
734     if util.Lsup(B[j]-C[j])>tol:
735     if verbose: print "non-symmetric PDE because B[%d]!=C[%d]"%(j,j)
736     out=False
737     if self.getNumSolutions()>1:
738     D=self.getCoefficientOfGeneralPDE("D")
739     if not D.isEmpty():
740 jgs 149 tol=util.Lsup(D)*self.__TOL
741 jgs 148 for i in range(self.getNumEquations()):
742     for k in range(self.getNumSolutions()):
743     if util.Lsup(D[i,k]-D[k,i])>tol:
744     if verbose: print "non-symmetric PDE because D[%d,%d]!=D[%d,%d]"%(i,k,k,i)
745     out=False
746 jgs 149 d=self.getCoefficientOfGeneralPDE("d")
747     if not d.isEmpty():
748     tol=util.Lsup(d)*self.__TOL
749     for i in range(self.getNumEquations()):
750     for k in range(self.getNumSolutions()):
751     if util.Lsup(d[i,k]-d[k,i])>tol:
752     if verbose: print "non-symmetric PDE because d[%d,%d]!=d[%d,%d]"%(i,k,k,i)
753     out=False
754     d_contact=self.getCoefficientOfGeneralPDE("d_contact")
755     if not d_contact.isEmpty():
756     tol=util.Lsup(d_contact)*self.__TOL
757     for i in range(self.getNumEquations()):
758     for k in range(self.getNumSolutions()):
759     if util.Lsup(d_contact[i,k]-d_contact[k,i])>tol:
760     if verbose: print "non-symmetric PDE because d_contact[%d,%d]!=d_contact[%d,%d]"%(i,k,k,i)
761     out=False
762 jgs 148 return out
763    
764     def getSolution(self,**options):
765 jgs 102 """
766 jgs 149 returns the solution of the PDE. If the solution is not valid the PDE is solved.
767 jgs 102
768 jgs 148 @return: the solution
769 jgs 149 @rtype: L{Data<escript.Data>}
770 jgs 148 @param options: solver options
771 jgs 149 @keyword verbose: True to get some information during PDE solution
772     @type verbose: C{bool}
773     @keyword reordering: reordering scheme to be used during elimination. Allowed values are
774     L{NO_REORDERING}, L{MINIMUM_FILL_IN}, L{NESTED_DISSECTION}
775     @keyword preconditioner: preconditioner method to be used. Allowed values are
776     L{SSOR}, L{ILU0}, L{ILUT}, L{JACOBI}
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     options[self.__METHOD_KEY]=self.getSolverMethod()
790 jgs 150 options[self.__PACKAGE_KEY]=self.getSolverPackage()
791 jgs 149 options[self.__SYMMETRY_KEY]=self.isSymmetric()
792 jgs 148 self.trace("PDE is resolved.")
793     self.trace("solver options: %s"%str(options))
794     self.__solution=mat.solve(f,options)
795     self.__solution_isValid=True
796     return self.__solution
797 jgs 102
798 jgs 148 def getFlux(self,u=None):
799     """
800 jgs 149 returns the flux M{J} for a given M{u}
801 jgs 102
802 jgs 149 M{J[i,j]=A[i,j,k,l]*grad(u[k])[l]+B[i,j,k]u[k]-X[i,j]}
803 jgs 102
804 jgs 151 or
805 jgs 149
806     M{J[j]=A[i,j]*grad(u)[l]+B[j]u-X[j]}
807    
808 jgs 148 @param u: argument in the flux. If u is not present or equals L{None} the current solution is used.
809 jgs 149 @type u: L{Data<escript.Data>} or None
810     @return: flux
811     @rtype: L{Data<escript.Data>}
812 jgs 148 """
813     if u==None: u=self.getSolution()
814     return util.tensormult(self.getCoefficientOfGeneralPDE("A"),util.grad(u))+util.matrixmult(self.getCoefficientOfGeneralPDE("B"),u)-util.self.getCoefficientOfGeneralPDE("X")
815     # =============================================================================
816     # solver settings:
817     # =============================================================================
818     def setSolverMethod(self,solver=None):
819 jgs 102 """
820 jgs 122 sets a new solver
821 jgs 148
822 jgs 149 @param solver: sets a new solver method.
823 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}.
824 jgs 102 """
825 jgs 150 if solver==None: solve=self.DEFAULT
826 jgs 102 if not solver==self.getSolverMethod():
827     self.__solver_method=solver
828     self.__checkMatrixType()
829 jgs 148 self.trace("New solver is %s"%self.getSolverMethodName())
830 jgs 102
831 jgs 148 def getSolverMethodName(self):
832     """
833     returns the name of the solver currently used
834    
835 jgs 149 @return: the name of the solver currently used.
836 jgs 148 @rtype: C{string}
837     """
838    
839     m=self.getSolverMethod()
840 jgs 150 p=self.getSolverPackage()
841     if m==self.DEFAULT: method="DEFAULT"
842     elif m==self.DIRECT: method= "DIRECT"
843     elif m==self.ITERATIVE: method= "ITERATIVE"
844     elif m==self.CHOLEVSKY: method= "CHOLEVSKY"
845     elif m==self.PCG: method= "PCG"
846     elif m==self.CR: method= "CR"
847     elif m==self.CGS: method= "CGS"
848     elif m==self.BICGSTAB: method= "BICGSTAB"
849     elif m==self.SSOR: method= "SSOR"
850     elif m==self.GMRES: method= "GMRES"
851     elif m==self.PRES20: method= "PRES20"
852     elif m==self.LUMPING: method= "LUMPING"
853     else : method="unknown"
854     if p==self.DEFAULT: package="DEFAULT"
855     elif p==self.PASO: package= "PASO"
856     elif p==self.MKL: package= "MKL"
857     elif p==self.SCSL: package= "SCSL"
858     elif p==self.UMFPACK: package= "UMFPACK"
859     else : method="unknown"
860     return "%s solver of %s package"%(method,package)
861 jgs 148
862 jgs 149
863 jgs 102 def getSolverMethod(self):
864     """
865 jgs 122 returns the solver method
866 jgs 149
867 jgs 151 @return: the solver method currently be used.
868 jgs 149 @rtype: C{int}
869 jgs 102 """
870     return self.__solver_method
871    
872 jgs 150 def setSolverPackage(self,package=None):
873     """
874     sets a new solver package
875    
876     @param solver: sets a new solver method.
877     @type solver: one of L{DEFAULT}, L{PASO} L{SCSL}, L{MKL}, L{UMLPACK}
878     """
879     if package==None: package=self.DEFAULT
880     if not package==self.getSolverPackage():
881     self.__solver_method=solver
882     self.__checkMatrixType()
883     self.trace("New solver is %s"%self.getSolverMethodName())
884    
885     def getSolverPackage(self):
886     """
887     returns the package of the solver
888    
889 jgs 151 @return: the solver package currently being used.
890 jgs 150 @rtype: C{int}
891     """
892     return self.__solver_package
893    
894 jgs 148 def isUsingLumping(self):
895     """
896     checks if matrix lumping is used a solver method
897    
898 jgs 149 @return: True is lumping is currently used a solver method.
899 jgs 148 @rtype: C{bool}
900     """
901     return self.getSolverMethod()==self.LUMPING
902    
903 jgs 102 def setTolerance(self,tol=1.e-8):
904     """
905 jgs 149 resets the tolerance for the solver method to tol where for an appropriate norm M{|.|}
906 jgs 148
907 jgs 149 M{|L{getResidual}()|<tol*|L{getRightHandSide}()|}
908 jgs 148
909 jgs 149 defines the stopping criterion.
910 jgs 148
911     @param tol: new tolerance for the solver. If the tol is lower then the current tolerence
912     the system will be resolved.
913 jgs 149 @type tol: positive C{float}
914 jgs 148 @raise ValueException: if tolerance is not positive.
915 jgs 102 """
916     if not tol>0:
917     raise ValueException,"Tolerance as to be positive"
918 jgs 148 if tol<self.getTolerance(): self.__invalidateSolution()
919     self.trace("New tolerance %e"%tol)
920 jgs 102 self.__tolerance=tol
921     return
922 jgs 148
923 jgs 102 def getTolerance(self):
924     """
925 jgs 122 returns the tolerance set for the solution
926 jgs 148
927     @return: tolerance currently used.
928     @rtype: C{float}
929 jgs 102 """
930     return self.__tolerance
931    
932 jgs 148 # =============================================================================
933     # symmetry flag:
934     # =============================================================================
935 jgs 102 def isSymmetric(self):
936     """
937 jgs 148 checks if symmetry is indicated.
938 jgs 149
939     @return: True is a symmetric PDE is indicated, otherwise False is returned
940     @rtype: C{bool}
941 jgs 102 """
942     return self.__sym
943    
944     def setSymmetryOn(self):
945     """
946 jgs 148 sets the symmetry flag.
947 jgs 102 """
948     if not self.isSymmetric():
949 jgs 148 self.trace("PDE is set to be symmetric")
950 jgs 102 self.__sym=True
951     self.__checkMatrixType()
952    
953     def setSymmetryOff(self):
954     """
955 jgs 148 removes the symmetry flag.
956 jgs 102 """
957     if self.isSymmetric():
958 jgs 148 self.trace("PDE is set to be unsymmetric")
959 jgs 102 self.__sym=False
960     self.__checkMatrixType()
961    
962     def setSymmetryTo(self,flag=False):
963 jgs 148 """
964     sets the symmetry flag to flag
965 jgs 149
966 jgs 148 @param flag: If flag, the symmetry flag is set otherwise the symmetry flag is released.
967     @type flag: C{bool}
968     """
969     if flag:
970     self.setSymmetryOn()
971     else:
972     self.setSymmetryOff()
973 jgs 102
974 jgs 148 # =============================================================================
975     # function space handling for the equation as well as the solution
976     # =============================================================================
977 jgs 102 def setReducedOrderOn(self):
978     """
979 jgs 148 switches on reduced order for solution and equation representation
980 jgs 149
981     @raise RuntimeError: if order reduction is altered after a coefficient has been set.
982 jgs 102 """
983     self.setReducedOrderForSolutionOn()
984     self.setReducedOrderForEquationOn()
985    
986     def setReducedOrderOff(self):
987     """
988 jgs 148 switches off reduced order for solution and equation representation
989 jgs 149
990     @raise RuntimeError: if order reduction is altered after a coefficient has been set.
991 jgs 102 """
992     self.setReducedOrderForSolutionOff()
993     self.setReducedOrderForEquationOff()
994    
995     def setReducedOrderTo(self,flag=False):
996     """
997 jgs 149 sets order reduction for both solution and equation representation according to flag.
998     @param flag: if flag is True, the order reduction is switched on for both solution and equation representation, otherwise or
999 jgs 148 if flag is not present order reduction is switched off
1000     @type flag: C{bool}
1001 jgs 149 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1002 jgs 102 """
1003     self.setReducedOrderForSolutionTo(flag)
1004     self.setReducedOrderForEquationTo(flag)
1005    
1006 jgs 148
1007 jgs 102 def setReducedOrderForSolutionOn(self):
1008     """
1009 jgs 148 switches on reduced order for solution representation
1010 jgs 149
1011     @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1012 jgs 102 """
1013 jgs 149 if not self.__reduce_solution_order:
1014     if self.__altered_coefficients:
1015     raise RuntimeError,"order cannot be altered after coefficients have been defined."
1016 jgs 148 self.trace("Reduced order is used to solution representation.")
1017 jgs 149 self.__reduce_solution_order=True
1018 jgs 148 self.__resetSystem()
1019 jgs 102
1020     def setReducedOrderForSolutionOff(self):
1021     """
1022 jgs 148 switches off reduced order for solution representation
1023 jgs 149
1024     @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1025 jgs 102 """
1026 jgs 149 if self.__reduce_solution_order:
1027     if self.__altered_coefficients:
1028     raise RuntimeError,"order cannot be altered after coefficients have been defined."
1029 jgs 148 self.trace("Full order is used to interpolate solution.")
1030 jgs 149 self.__reduce_solution_order=False
1031 jgs 148 self.__resetSystem()
1032 jgs 102
1033     def setReducedOrderForSolutionTo(self,flag=False):
1034     """
1035 jgs 122 sets order for test functions according to flag
1036 jgs 102
1037 jgs 149 @param flag: if flag is True, the order reduction is switched on for solution representation, otherwise or
1038 jgs 148 if flag is not present order reduction is switched off
1039     @type flag: C{bool}
1040 jgs 149 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1041 jgs 102 """
1042     if flag:
1043     self.setReducedOrderForSolutionOn()
1044     else:
1045     self.setReducedOrderForSolutionOff()
1046 jgs 148
1047 jgs 102 def setReducedOrderForEquationOn(self):
1048     """
1049 jgs 148 switches on reduced order for equation representation
1050 jgs 149
1051     @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1052 jgs 102 """
1053 jgs 149 if not self.__reduce_equation_order:
1054     if self.__altered_coefficients:
1055     raise RuntimeError,"order cannot be altered after coefficients have been defined."
1056 jgs 148 self.trace("Reduced order is used for test functions.")
1057 jgs 149 self.__reduce_equation_order=True
1058 jgs 148 self.__resetSystem()
1059 jgs 102
1060     def setReducedOrderForEquationOff(self):
1061     """
1062 jgs 148 switches off reduced order for equation representation
1063 jgs 149
1064     @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1065 jgs 102 """
1066 jgs 149 if self.__reduce_equation_order:
1067     if self.__altered_coefficients:
1068     raise RuntimeError,"order cannot be altered after coefficients have been defined."
1069 jgs 148 self.trace("Full order is used for test functions.")
1070 jgs 149 self.__reduce_equation_order=False
1071 jgs 148 self.__resetSystem()
1072 jgs 102
1073     def setReducedOrderForEquationTo(self,flag=False):
1074     """
1075 jgs 122 sets order for test functions according to flag
1076 jgs 102
1077 jgs 149 @param flag: if flag is True, the order reduction is switched on for equation representation, otherwise or
1078 jgs 148 if flag is not present order reduction is switched off
1079     @type flag: C{bool}
1080 jgs 149 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1081 jgs 102 """
1082     if flag:
1083     self.setReducedOrderForEquationOn()
1084     else:
1085     self.setReducedOrderForEquationOff()
1086 jgs 148
1087     # =============================================================================
1088     # private method:
1089     # =============================================================================
1090     def __checkMatrixType(self):
1091     """
1092     reassess the matrix type and, if a new matrix is needed, resets the system.
1093     """
1094 jgs 150 new_matrix_type=self.getDomain().getSystemMatrixTypeId(self.getSolverMethod(),self.getSolverPackage(),self.isSymmetric())
1095 jgs 148 if not new_matrix_type==self.__matrix_type:
1096     self.trace("Matrix type is now %d."%new_matrix_type)
1097     self.__matrix_type=new_matrix_type
1098     self.__resetSystem()
1099 jgs 149 #
1100 jgs 148 # rebuild switches :
1101 jgs 149 #
1102 jgs 148 def __invalidateSolution(self):
1103     """
1104     indicates the PDE has to be resolved if the solution is requested
1105     """
1106     if self.__solution_isValid: self.trace("PDE has to be resolved.")
1107     self.__solution_isValid=False
1108    
1109     def __invalidateOperator(self):
1110     """
1111     indicates the operator has to be rebuilt next time it is used
1112     """
1113 jgs 149 if self.__operator_is_Valid: self.trace("Operator has to be rebuilt.")
1114 jgs 148 self.__invalidateSolution()
1115 jgs 149 self.__operator_is_Valid=False
1116 jgs 148
1117     def __invalidateRightHandSide(self):
1118     """
1119     indicates the right hand side has to be rebuild next time it is used
1120     """
1121     if self.__righthandside_isValid: self.trace("Right hand side has to be rebuilt.")
1122     self.__invalidateSolution()
1123     self.__righthandside_isValid=False
1124    
1125     def __invalidateSystem(self):
1126     """
1127     annonced that everthing has to be rebuild:
1128     """
1129     if self.__righthandside_isValid: self.trace("System has to be rebuilt.")
1130     self.__invalidateSolution()
1131     self.__invalidateOperator()
1132     self.__invalidateRightHandSide()
1133    
1134     def __resetSystem(self):
1135     """
1136     annonced that everthing has to be rebuild:
1137     """
1138     self.trace("New System is built from scratch.")
1139     self.__operator=escript.Operator()
1140 jgs 149 self.__operator_is_Valid=False
1141 jgs 148 self.__righthandside=escript.Data()
1142     self.__righthandside_isValid=False
1143     self.__solution=escript.Data()
1144     self.__solution_isValid=False
1145 jgs 149 #
1146 jgs 148 # system initialization:
1147 jgs 149 #
1148 jgs 121 def __getNewOperator(self):
1149 jgs 102 """
1150 jgs 148 returns an instance of a new operator
1151 jgs 102 """
1152 jgs 148 self.trace("New operator is allocated.")
1153 jgs 102 return self.getDomain().newOperator( \
1154     self.getNumEquations(), \
1155     self.getFunctionSpaceForEquation(), \
1156     self.getNumSolutions(), \
1157     self.getFunctionSpaceForSolution(), \
1158     self.__matrix_type)
1159    
1160 jgs 148 def __getNewRightHandSide(self):
1161 jgs 102 """
1162 jgs 148 returns an instance of a new right hand side
1163 jgs 102 """
1164 jgs 148 self.trace("New right hand side is allocated.")
1165 jgs 121 if self.getNumEquations()>1:
1166 jgs 148 return escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)
1167 jgs 121 else:
1168 jgs 148 return escript.Data(0.,(),self.getFunctionSpaceForEquation(),True)
1169 jgs 102
1170 jgs 121 def __getNewSolution(self):
1171 jgs 102 """
1172 jgs 148 returns an instance of a new solution
1173 jgs 102 """
1174 jgs 148 self.trace("New solution is allocated.")
1175 jgs 121 if self.getNumSolutions()>1:
1176     return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)
1177     else:
1178     return escript.Data(0.,(),self.getFunctionSpaceForSolution(),True)
1179 jgs 102
1180 jgs 148 def __makeFreshSolution(self):
1181     """
1182     makes sure that the solution is instantiated and returns it initialized by zeros
1183     """
1184     if self.__solution.isEmpty():
1185     self.__solution=self.__getNewSolution()
1186     else:
1187     self.__solution*=0
1188     self.trace("Solution is reset to zero.")
1189     return self.__solution
1190    
1191     def __makeFreshRightHandSide(self):
1192     """
1193     makes sure that the right hand side is instantiated and returns it initialized by zeros
1194     """
1195     if self.__righthandside.isEmpty():
1196     self.__righthandside=self.__getNewRightHandSide()
1197     else:
1198     self.__righthandside*=0
1199     self.trace("Right hand side is reset to zero.")
1200     return self.__righthandside
1201    
1202 jgs 121 def __makeFreshOperator(self):
1203 jgs 102 """
1204 jgs 148 makes sure that the operator is instantiated and returns it initialized by zeros
1205 jgs 102 """
1206     if self.__operator.isEmpty():
1207 jgs 121 self.__operator=self.__getNewOperator()
1208 jgs 102 else:
1209 jgs 149 self.__operator.resetValues()
1210 jgs 148 self.trace("Operator reset to zero")
1211 jgs 102 return self.__operator
1212    
1213 jgs 148 def __applyConstraint(self):
1214     """
1215     applies the constraints defined by q and r to the system
1216     """
1217     if not self.isUsingLumping():
1218     q=self.getCoefficientOfGeneralPDE("q")
1219     r=self.getCoefficientOfGeneralPDE("r")
1220     if not q.isEmpty() and not self.__operator.isEmpty():
1221     # q is the row and column mask to indicate where constraints are set:
1222     row_q=escript.Data(q,self.getFunctionSpaceForEquation())
1223     col_q=escript.Data(q,self.getFunctionSpaceForSolution())
1224     u=self.__getNewSolution()
1225     if r.isEmpty():
1226     r_s=self.__getNewSolution()
1227     else:
1228     r_s=escript.Data(r,self.getFunctionSpaceForSolution())
1229     u.copyWithMask(r_s,col_q)
1230 jgs 149 if not self.__righthandside.isEmpty():
1231 jgs 148 self.__righthandside-=self.__operator*u
1232     self.__righthandside=self.copyConstraint(self.__righthandside)
1233     self.__operator.nullifyRowsAndCols(row_q,col_q,1.)
1234     # =============================================================================
1235     # function giving access to coefficients of the general PDE:
1236     # =============================================================================
1237     def getCoefficientOfGeneralPDE(self,name):
1238 jgs 108 """
1239 jgs 148 return the value of the coefficient name of the general PDE.
1240    
1241 jgs 149 @note: This method is called by the assembling routine it can be overwritten
1242 jgs 148 to map coefficients of a particular PDE to the general PDE.
1243 jgs 149 @param name: name of the coefficient requested.
1244 jgs 148 @type name: C{string}
1245 jgs 149 @return: the value of the coefficient name
1246     @rtype: L{Data<escript.Data>}
1247 jgs 148 @raise IllegalCoefficient: if name is not one of coefficients
1248 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}.
1249 jgs 108 """
1250 jgs 148 if self.hasCoefficientOfGeneralPDE(name):
1251     return self.getCoefficient(name)
1252     else:
1253     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1254 jgs 108
1255 jgs 148 def hasCoefficientOfGeneralPDE(self,name):
1256 jgs 108 """
1257 jgs 148 checks if name is a the name of a coefficient of the general PDE.
1258 jgs 149
1259 jgs 148 @param name: name of the coefficient enquired.
1260     @type name: C{string}
1261 jgs 149 @return: True if name is the name of a coefficient of the general PDE. Otherwise False.
1262     @rtype: C{bool}
1263    
1264 jgs 108 """
1265 jgs 148 return self.__COEFFICIENTS_OF_GENEARL_PDE.has_key(name)
1266 jgs 108
1267 jgs 148 def createCoefficientOfGeneralPDE(self,name):
1268 jgs 108 """
1269 jgs 148 returns a new instance of a coefficient for coefficient name of the general PDE
1270    
1271     @param name: name of the coefficient requested.
1272     @type name: C{string}
1273 jgs 149 @return: a coefficient name initialized to 0.
1274     @rtype: L{Data<escript.Data>}
1275 jgs 148 @raise IllegalCoefficient: if name is not one of coefficients
1276 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}.
1277 jgs 108 """
1278 jgs 148 if self.hasCoefficientOfGeneralPDE(name):
1279     return escript.Data(0,self.getShapeOfCoefficientOfGeneralPDE(name),self.getFunctionSpaceForCoefficientOfGeneralPDE(name))
1280 jgs 108 else:
1281 jgs 148 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1282 jgs 108
1283 jgs 148 def getFunctionSpaceForCoefficientOfGeneralPDE(self,name):
1284 jgs 108 """
1285 jgs 149 return the L{FunctionSpace<escript.FunctionSpace>} to be used for coefficient name of the general PDE
1286 jgs 148
1287     @param name: name of the coefficient enquired.
1288     @type name: C{string}
1289 jgs 149 @return: the function space to be used for coefficient name
1290     @rtype: L{FunctionSpace<escript.FunctionSpace>}
1291 jgs 148 @raise IllegalCoefficient: if name is not one of coefficients
1292 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}.
1293 jgs 108 """
1294 jgs 148 if self.hasCoefficientOfGeneralPDE(name):
1295     return self.__COEFFICIENTS_OF_GENEARL_PDE[name].getFunctionSpace(self.getDomain())
1296 jgs 108 else:
1297 jgs 148 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1298 jgs 108
1299 jgs 148 def getShapeOfCoefficientOfGeneralPDE(self,name):
1300     """
1301     return the shape of the coefficient name of the general PDE
1302 jgs 108
1303 jgs 148 @param name: name of the coefficient enquired.
1304     @type name: C{string}
1305 jgs 149 @return: the shape of the coefficient name
1306     @rtype: C{tuple} of C{int}
1307 jgs 148 @raise IllegalCoefficient: if name is not one of coefficients
1308 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}.
1309 jgs 148 """
1310     if self.hasCoefficientOfGeneralPDE(name):
1311     return self.__COEFFICIENTS_OF_GENEARL_PDE[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())
1312     else:
1313     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1314 jgs 108
1315 jgs 148 # =============================================================================
1316     # functions giving access to coefficients of a particular PDE implementation:
1317     # =============================================================================
1318     def getCoefficient(self,name):
1319     """
1320     returns the value of the coefficient name
1321 jgs 108
1322 jgs 149 @param name: name of the coefficient requested.
1323 jgs 148 @type name: C{string}
1324 jgs 149 @return: the value of the coefficient name
1325     @rtype: L{Data<escript.Data>}
1326 jgs 148 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1327     """
1328     if self.hasCoefficient(name):
1329     return self.COEFFICIENTS[name].getValue()
1330     else:
1331     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1332 jgs 108
1333 jgs 148 def hasCoefficient(self,name):
1334     """
1335     return True if name is the name of a coefficient
1336 jgs 108
1337 jgs 148 @param name: name of the coefficient enquired.
1338     @type name: C{string}
1339 jgs 149 @return: True if name is the name of a coefficient of the general PDE. Otherwise False.
1340     @rtype: C{bool}
1341 jgs 148 """
1342     return self.COEFFICIENTS.has_key(name)
1343 jgs 108
1344 jgs 148 def createCoefficient(self, name):
1345     """
1346 jgs 149 create a L{Data<escript.Data>} object corresponding to coefficient name
1347 jgs 108
1348 jgs 149 @return: a coefficient name initialized to 0.
1349     @rtype: L{Data<escript.Data>}
1350 jgs 148 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1351     """
1352     if self.hasCoefficient(name):
1353     return escript.Data(0.,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))
1354     else:
1355     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1356    
1357     def getFunctionSpaceForCoefficient(self,name):
1358     """
1359 jgs 149 return the L{FunctionSpace<escript.FunctionSpace>} to be used for coefficient name
1360 jgs 148
1361     @param name: name of the coefficient enquired.
1362     @type name: C{string}
1363 jgs 149 @return: the function space to be used for coefficient name
1364     @rtype: L{FunctionSpace<escript.FunctionSpace>}
1365 jgs 148 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1366     """
1367     if self.hasCoefficient(name):
1368     return self.COEFFICIENTS[name].getFunctionSpace(self.getDomain())
1369     else:
1370     raise ValueError,"unknown coefficient %s requested"%name
1371     def getShapeOfCoefficient(self,name):
1372     """
1373 jgs 149 return the shape of the coefficient name
1374 jgs 148
1375     @param name: name of the coefficient enquired.
1376     @type name: C{string}
1377 jgs 149 @return: the shape of the coefficient name
1378     @rtype: C{tuple} of C{int}
1379 jgs 148 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1380     """
1381     if self.hasCoefficient(name):
1382     return self.COEFFICIENTS[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())
1383     else:
1384     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1385    
1386     def resetCoefficients(self):
1387     """
1388     resets all coefficients to there default values.
1389     """
1390     for i in self.COEFFICIENTS.iterkeys():
1391     self.COEFFICIENTS[i].resetValue()
1392    
1393     def alteredCoefficient(self,name):
1394     """
1395     announce that coefficient name has been changed
1396    
1397     @param name: name of the coefficient enquired.
1398     @type name: C{string}
1399     @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1400 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.
1401 jgs 148 """
1402     if self.hasCoefficient(name):
1403     self.trace("Coefficient %s has been altered."%name)
1404 jgs 149 if not ((name=="q" or name=="r") and self.isUsingLumping()):
1405     if self.COEFFICIENTS[name].isAlteringOperator(): self.__invalidateOperator()
1406     if self.COEFFICIENTS[name].isAlteringRightHandSide(): self.__invalidateRightHandSide()
1407 jgs 148 else:
1408     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1409    
1410     def copyConstraint(self,u):
1411 jgs 108 """
1412 jgs 149 copies the constraint into u and returns u.
1413    
1414 jgs 148 @param u: a function of rank 0 is a single PDE is solved and of shape (numSolution,) for a system of PDEs
1415 jgs 149 @type u: L{Data<escript.Data>}
1416     @return: the input u modified by the constraints.
1417     @rtype: L{Data<escript.Data>}
1418     @warning: u is altered if it has the appropriate L{FunctionSpace<escript.FunctionSpace>}
1419 jgs 108 """
1420 jgs 148 q=self.getCoefficientOfGeneralPDE("q")
1421     r=self.getCoefficientOfGeneralPDE("r")
1422     if not q.isEmpty():
1423     if u.isEmpty(): u=escript.Data(0.,q.getShape(),q.getFunctionSpace())
1424     if r.isEmpty():
1425     r=escript.Data(0,u.getShape(),u.getFunctionSpace())
1426     else:
1427     r=escript.Data(r,u.getFunctionSpace())
1428     u.copyWithMask(r,escript.Data(q,u.getFunctionSpace()))
1429     return u
1430    
1431     def setValue(self,**coefficients):
1432     """
1433     sets new values to coefficients
1434    
1435 jgs 149 @param coefficients: new values assigned to coefficients
1436     @keyword A: value for coefficient A.
1437     @type A: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1438 jgs 148 @keyword B: value for coefficient B
1439 jgs 149 @type B: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1440 jgs 148 @keyword C: value for coefficient C
1441 jgs 149 @type C: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1442 jgs 148 @keyword D: value for coefficient D
1443 jgs 149 @type D: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1444 jgs 148 @keyword X: value for coefficient X
1445 jgs 149 @type X: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1446 jgs 148 @keyword Y: value for coefficient Y
1447 jgs 149 @type Y: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1448 jgs 148 @keyword d: value for coefficient d
1449 jgs 149 @type d: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1450 jgs 148 @keyword y: value for coefficient y
1451 jgs 149 @type y: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1452 jgs 148 @keyword d_contact: value for coefficient d_contact
1453 jgs 149 @type d_contact: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnContactOne<escript.FunctionOnContactOne>}.
1454     or L{FunctionOnContactZero<escript.FunctionOnContactZero>}.
1455 jgs 148 @keyword y_contact: value for coefficient y_contact
1456 jgs 149 @type y_contact: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnContactOne<escript.FunctionOnContactOne>}.
1457     or L{FunctionOnContactZero<escript.FunctionOnContactZero>}.
1458 jgs 148 @keyword r: values prescribed to the solution at the locations of constraints
1459 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>}
1460 jgs 148 depending of reduced order is used for the solution.
1461     @keyword q: mask for location of constraints
1462 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>}
1463 jgs 148 depending of reduced order is used for the representation of the equation.
1464     @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1465     """
1466 jgs 108 # check if the coefficients are legal:
1467     for i in coefficients.iterkeys():
1468     if not self.hasCoefficient(i):
1469 jgs 148 raise IllegalCoefficient,"Attempt to set unknown coefficient %s"%i
1470 jgs 108 # if the number of unknowns or equations is still unknown we try to estimate them:
1471 jgs 148 if self.__numEquations==None or self.__numSolutions==None:
1472 jgs 108 for i,d in coefficients.iteritems():
1473     if hasattr(d,"shape"):
1474     s=d.shape
1475     elif hasattr(d,"getShape"):
1476     s=d.getShape()
1477     else:
1478     s=numarray.array(d).shape
1479     if s!=None:
1480     # get number of equations and number of unknowns:
1481 jgs 148 res=self.COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(self.getDomain(),s)
1482 jgs 108 if res==None:
1483 jgs 148 raise IllegalCoefficientValue,"Illegal shape %s of coefficient %s"%(s,i)
1484 jgs 108 else:
1485 jgs 148 if self.__numEquations==None: self.__numEquations=res[0]
1486     if self.__numSolutions==None: self.__numSolutions=res[1]
1487     if self.__numEquations==None: raise UndefinedPDEError,"unidententified number of equations"
1488     if self.__numSolutions==None: raise UndefinedPDEError,"unidententified number of solutions"
1489 jgs 108 # now we check the shape of the coefficient if numEquations and numSolutions are set:
1490     for i,d in coefficients.iteritems():
1491 jgs 148 try:
1492 jgs 149 self.COEFFICIENTS[i].setValue(self.getDomain(),self.getNumEquations(),self.getNumSolutions(),self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
1493 jgs 148 except IllegalCoefficientValue,m:
1494     raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
1495 jgs 108 self.alteredCoefficient(i)
1496    
1497 jgs 149 self.__altered_coefficients=True
1498 jgs 148 # check if the systrem is inhomogeneous:
1499     if len(coefficients)>0 and not self.isUsingLumping():
1500     q=self.getCoefficientOfGeneralPDE("q")
1501     r=self.getCoefficientOfGeneralPDE("r")
1502     homogeneous_constraint=True
1503     if not q.isEmpty() and not r.isEmpty():
1504 jgs 149 if util.Lsup(q*r)>=1.e-13*util.Lsup(r):
1505 jgs 148 self.trace("Inhomogeneous constraint detected.")
1506     self.__invalidateSystem()
1507 jgs 108
1508     def getSystem(self):
1509     """
1510 jgs 122 return the operator and right hand side of the PDE
1511 jgs 149
1512     @return: the discrete version of the PDE
1513     @rtype: C{tuple} of L{Operator,<escript.Operator>} and L{Data<escript.Data>}.
1514 jgs 108 """
1515 jgs 149 if not self.__operator_is_Valid or not self.__righthandside_isValid:
1516 jgs 108 if self.isUsingLumping():
1517 jgs 149 if not self.__operator_is_Valid:
1518 jgs 148 if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution(): raise TypeError,"Lumped matrix requires same order for equations and unknowns"
1519 jgs 149 if not self.getCoefficientOfGeneralPDE("A").isEmpty(): raise Warning,"Using coefficient A in lumped matrix can produce wrong results"
1520     if not self.getCoefficientOfGeneralPDE("B").isEmpty(): raise Warning,"Using coefficient B in lumped matrix can produce wrong results"
1521     if not self.getCoefficientOfGeneralPDE("C").isEmpty(): raise Warning,"Using coefficient C in lumped matrix can produce wrong results"
1522 jgs 121 mat=self.__getNewOperator()
1523 jgs 108 self.getDomain().addPDEToSystem(mat,escript.Data(), \
1524 jgs 148 self.getCoefficientOfGeneralPDE("A"), \
1525     self.getCoefficientOfGeneralPDE("B"), \
1526     self.getCoefficientOfGeneralPDE("C"), \
1527     self.getCoefficientOfGeneralPDE("D"), \
1528 jgs 108 escript.Data(), \
1529     escript.Data(), \
1530 jgs 148 self.getCoefficientOfGeneralPDE("d"), \
1531 jgs 108 escript.Data(),\
1532 jgs 148 self.getCoefficientOfGeneralPDE("d_contact"), \
1533 jgs 108 escript.Data())
1534 jgs 148 self.__operator=1./(mat*escript.Data(1,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True))
1535     del mat
1536     self.trace("New lumped operator has been built.")
1537 jgs 149 self.__operator_is_Valid=True
1538 jgs 108 if not self.__righthandside_isValid:
1539 jgs 121 self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \
1540 jgs 148 self.getCoefficientOfGeneralPDE("X"), \
1541     self.getCoefficientOfGeneralPDE("Y"),\
1542     self.getCoefficientOfGeneralPDE("y"),\
1543     self.getCoefficientOfGeneralPDE("y_contact"))
1544     self.trace("New right hand side as been built.")
1545 jgs 108 self.__righthandside_isValid=True
1546     else:
1547 jgs 149 if not self.__operator_is_Valid and not self.__righthandside_isValid:
1548 jgs 121 self.getDomain().addPDEToSystem(self.__makeFreshOperator(),self.__makeFreshRightHandSide(), \
1549 jgs 148 self.getCoefficientOfGeneralPDE("A"), \
1550     self.getCoefficientOfGeneralPDE("B"), \
1551     self.getCoefficientOfGeneralPDE("C"), \
1552     self.getCoefficientOfGeneralPDE("D"), \
1553     self.getCoefficientOfGeneralPDE("X"), \
1554     self.getCoefficientOfGeneralPDE("Y"), \
1555     self.getCoefficientOfGeneralPDE("d"), \
1556     self.getCoefficientOfGeneralPDE("y"), \
1557     self.getCoefficientOfGeneralPDE("d_contact"), \
1558     self.getCoefficientOfGeneralPDE("y_contact"))
1559 jgs 108 self.__applyConstraint()
1560 jgs 148 self.__righthandside=self.copyConstraint(self.__righthandside)
1561     self.trace("New system has been built.")
1562 jgs 149 self.__operator_is_Valid=True
1563 jgs 108 self.__righthandside_isValid=True
1564     elif not self.__righthandside_isValid:
1565 jgs 121 self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \
1566 jgs 148 self.getCoefficientOfGeneralPDE("X"), \
1567     self.getCoefficientOfGeneralPDE("Y"),\
1568     self.getCoefficientOfGeneralPDE("y"),\
1569     self.getCoefficientOfGeneralPDE("y_contact"))
1570     self.__righthandside=self.copyConstraint(self.__righthandside)
1571     self.trace("New right hand side has been built.")
1572 jgs 108 self.__righthandside_isValid=True
1573 jgs 149 elif not self.__operator_is_Valid:
1574 jgs 121 self.getDomain().addPDEToSystem(self.__makeFreshOperator(),escript.Data(), \
1575 jgs 148 self.getCoefficientOfGeneralPDE("A"), \
1576     self.getCoefficientOfGeneralPDE("B"), \
1577     self.getCoefficientOfGeneralPDE("C"), \
1578     self.getCoefficientOfGeneralPDE("D"), \
1579 jgs 108 escript.Data(), \
1580     escript.Data(), \
1581 jgs 148 self.getCoefficientOfGeneralPDE("d"), \
1582 jgs 108 escript.Data(),\
1583 jgs 148 self.getCoefficientOfGeneralPDE("d_contact"), \
1584 jgs 108 escript.Data())
1585     self.__applyConstraint()
1586 jgs 148 self.trace("New operator has been built.")
1587 jgs 149 self.__operator_is_Valid=True
1588 jgs 108 return (self.__operator,self.__righthandside)
1589 jgs 102
1590    
1591 jgs 149 class Poisson(LinearPDE):
1592     """
1593     Class to define a Poisson equation problem, which is genear L{LinearPDE} of the form
1594 jgs 102
1595 jgs 149 M{-grad(grad(u)[j])[j] = f}
1596 jgs 102
1597 jgs 149 with natural boundary conditons
1598    
1599     M{n[j]*grad(u)[j] = 0 }
1600    
1601     and constraints:
1602    
1603     M{u=0} where M{q>0}
1604    
1605 jgs 108 """
1606 jgs 148
1607 gross 328 def __init__(self,domain,debug=False):
1608 jgs 149 """
1609     initializes a new Poisson equation
1610 jgs 104
1611 jgs 149 @param domain: domain of the PDE
1612     @type domain: L{Domain<escript.Domain>}
1613     @param debug: if True debug informations are printed.
1614 jgs 102
1615 jgs 149 """
1616 jgs 151 super(Poisson, self).__init__(domain,1,1,debug)
1617 jgs 149 self.COEFFICIENTS={"f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1618     "q": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
1619     self.setSymmetryOn()
1620 jgs 102
1621 jgs 149 def setValue(self,**coefficients):
1622     """
1623     sets new values to coefficients
1624 jgs 102
1625 jgs 149 @param coefficients: new values assigned to coefficients
1626     @keyword f: value for right hand side M{f}
1627     @type f: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1628     @keyword q: mask for location of constraints
1629     @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>}
1630     depending of reduced order is used for the representation of the equation.
1631     @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1632     """
1633 jgs 151 super(Poisson, self).setValue(**coefficients)
1634 jgs 102
1635 jgs 149 def getCoefficientOfGeneralPDE(self,name):
1636     """
1637     return the value of the coefficient name of the general PDE
1638     @param name: name of the coefficient requested.
1639     @type name: C{string}
1640     @return: the value of the coefficient name
1641     @rtype: L{Data<escript.Data>}
1642     @raise IllegalCoefficient: if name is not one of coefficients
1643     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}.
1644     @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
1645     """
1646     if name == "A" :
1647 gross 328 return escript.Data(util.kronecker(self.getDim()),escript.Function(self.getDomain()))
1648 jgs 149 elif name == "B" :
1649     return escript.Data()
1650     elif name == "C" :
1651     return escript.Data()
1652     elif name == "D" :
1653     return escript.Data()
1654     elif name == "X" :
1655     return escript.Data()
1656     elif name == "Y" :
1657     return self.getCoefficient("f")
1658     elif name == "d" :
1659     return escript.Data()
1660     elif name == "y" :
1661     return escript.Data()
1662     elif name == "d_contact" :
1663     return escript.Data()
1664     elif name == "y_contact" :
1665     return escript.Data()
1666     elif name == "r" :
1667     return escript.Data()
1668     elif name == "q" :
1669     return self.getCoefficient("q")
1670     else:
1671     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1672 jgs 102
1673 jgs 149 class Helmholtz(LinearPDE):
1674     """
1675     Class to define a Helmhotz equation problem, which is genear L{LinearPDE} of the form
1676    
1677     M{S{omega}*u - grad(k*grad(u)[j])[j] = f}
1678    
1679     with natural boundary conditons
1680    
1681     M{k*n[j]*grad(u)[j] = g- S{alpha}u }
1682    
1683 jgs 122 and constraints:
1684 jgs 102
1685 jgs 149 M{u=r} where M{q>0}
1686    
1687 jgs 110 """
1688 jgs 149
1689     def __init__(self,domain,debug=False):
1690     """
1691     initializes a new Poisson equation
1692    
1693     @param domain: domain of the PDE
1694     @type domain: L{Domain<escript.Domain>}
1695     @param debug: if True debug informations are printed.
1696    
1697     """
1698 jgs 151 super(Helmholtz, self).__init__(domain,1,1,debug)
1699 jgs 149 self.COEFFICIENTS={"omega": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
1700     "k": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
1701     "f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1702     "alpha": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
1703     "g": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1704     "r": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH),
1705     "q": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
1706     self.setSymmetryOn()
1707    
1708     def setValue(self,**coefficients):
1709     """
1710     sets new values to coefficients
1711    
1712     @param coefficients: new values assigned to coefficients
1713     @keyword omega: value for coefficient M{S{omega}}
1714     @type omega: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1715     @keyword k: value for coefficeint M{k}
1716     @type k: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1717     @keyword f: value for right hand side M{f}
1718     @type f: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1719     @keyword alpha: value for right hand side M{S{alpha}}
1720     @type alpha: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1721     @keyword g: value for right hand side M{g}
1722     @type g: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1723     @keyword r: prescribed values M{r} for the solution in constraints.
1724     @type r: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1725     depending of reduced order is used for the representation of the equation.
1726     @keyword q: mask for location of constraints
1727     @type q: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1728     depending of reduced order is used for the representation of the equation.
1729     @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1730     """
1731 jgs 151 super(Helmholtz, self).setValue(**coefficients)
1732 jgs 149
1733     def getCoefficientOfGeneralPDE(self,name):
1734     """
1735     return the value of the coefficient name of the general PDE
1736    
1737     @param name: name of the coefficient requested.
1738     @type name: C{string}
1739     @return: the value of the coefficient name
1740     @rtype: L{Data<escript.Data>}
1741     @raise IllegalCoefficient: if name is not one of coefficients
1742     "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}.
1743     @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
1744     """
1745     if name == "A" :
1746     return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))*self.getCoefficient("k")
1747     elif name == "B" :
1748     return escript.Data()
1749     elif name == "C" :
1750     return escript.Data()
1751     elif name == "D" :
1752     return self.getCoefficient("omega")
1753     elif name == "X" :
1754     return escript.Data()
1755     elif name == "Y" :
1756     return self.getCoefficient("f")
1757     elif name == "d" :
1758     return self.getCoefficient("alpha")
1759     elif name == "y" :
1760     return self.getCoefficient("g")
1761     elif name == "d_contact" :
1762     return escript.Data()
1763     elif name == "y_contact" :
1764     return escript.Data()
1765     elif name == "r" :
1766     return self.getCoefficient("r")
1767     elif name == "q" :
1768     return self.getCoefficient("q")
1769     else:
1770     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1771    
1772     class LameEquation(LinearPDE):
1773     """
1774     Class to define a Lame equation problem:
1775    
1776     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] }
1777    
1778     with natural boundary conditons:
1779    
1780     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] }
1781    
1782     and constraints:
1783    
1784     M{u[i]=r[i]} where M{q[i]>0}
1785    
1786     """
1787    
1788     def __init__(self,domain,debug=False):
1789 jgs 151 super(LameEquation, self).__init__(domain,\
1790     domain.getDim(),domain.getDim(),debug)
1791     self.COEFFICIENTS={ "lame_lambda" : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),
1792 jgs 149 "lame_mu" : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),
1793     "F" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1794     "sigma" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM),PDECoefficient.RIGHTHANDSIDE),
1795     "f" : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1796     "r" : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH),
1797     "q" : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
1798 jgs 151 self.setSymmetryOn()
1799 jgs 149
1800     def setValue(self,**coefficients):
1801     """
1802     sets new values to coefficients
1803    
1804     @param coefficients: new values assigned to coefficients
1805     @keyword lame_mu: value for coefficient M{S{mu}}
1806     @type lame_mu: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1807     @keyword lame_lambda: value for coefficient M{S{lambda}}
1808     @type lame_lambda: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1809     @keyword F: value for internal force M{F}
1810     @type F: any type that can be casted to L{Vector<escript.Vector>} object on L{Function<escript.Function>}
1811     @keyword sigma: value for initial stress M{S{sigma}}
1812     @type sigma: any type that can be casted to L{Tensor<escript.Tensor>} object on L{Function<escript.Function>}
1813     @keyword f: value for extrenal force M{f}
1814     @type f: any type that can be casted to L{Vector<escript.Vector>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}
1815     @keyword r: prescribed values M{r} for the solution in constraints.
1816     @type r: any type that can be casted to L{Vector<escript.Vector>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1817     depending of reduced order is used for the representation of the equation.
1818     @keyword q: mask for location of constraints
1819     @type q: any type that can be casted to L{Vector<escript.Vector>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1820     depending of reduced order is used for the representation of the equation.
1821     @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1822     """
1823 jgs 151 super(LameEquation, self).setValue(**coefficients)
1824 jgs 149
1825     def getCoefficientOfGeneralPDE(self,name):
1826     """
1827     return the value of the coefficient name of the general PDE
1828    
1829     @param name: name of the coefficient requested.
1830     @type name: C{string}
1831     @return: the value of the coefficient name
1832     @rtype: L{Data<escript.Data>}
1833     @raise IllegalCoefficient: if name is not one of coefficients
1834     "A", M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact}, M{r} or M{q}.
1835     @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
1836     """
1837     if name == "A" :
1838     out =self.createCoefficientOfGeneralPDE("A")
1839     for i in range(self.getDim()):
1840     for j in range(self.getDim()):
1841     out[i,i,j,j] += self.getCoefficient("lame_lambda")
1842     out[i,j,j,i] += self.getCoefficient("lame_mu")
1843     out[i,j,i,j] += self.getCoefficient("lame_mu")
1844     return out
1845     elif name == "B" :
1846     return escript.Data()
1847     elif name == "C" :
1848     return escript.Data()
1849     elif name == "D" :
1850     return escript.Data()
1851     elif name == "X" :
1852     return self.getCoefficient("sigma")
1853     elif name == "Y" :
1854     return self.getCoefficient("F")
1855     elif name == "d" :
1856     return escript.Data()
1857     elif name == "y" :
1858     return self.getCoefficient("f")
1859     elif name == "d_contact" :
1860     return escript.Data()
1861     elif name == "y_contact" :
1862     return escript.Data()
1863     elif name == "r" :
1864     return self.getCoefficient("r")
1865     elif name == "q" :
1866     return self.getCoefficient("q")
1867     else:
1868     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1869    
1870     class AdvectivePDE(LinearPDE):
1871     """
1872     In cases of PDEs dominated by the advection terms M{B} and M{C} against the adevctive terms M{A}
1873     up-winding has been used. The L{AdvectivePDE} class applies SUPG upwinding to the advective terms.
1874    
1875     In the following we set
1876    
1877     M{Z[j]=C[j]-B[j]}
1878    
1879 jgs 151 or
1880 jgs 149
1881     M{Z[i,k,l]=C[i,k,l]-B[i,l,k]}
1882    
1883 jgs 151 To measure the dominance of the advective terms over the diffusive term M{A} the
1884 jgs 149 X{Pelclet number} M{P} is used. It is defined as
1885 jgs 151
1886 jgs 149 M{P=h|Z|/(2|A|)}
1887    
1888 jgs 151 where M{|.|} denotes the L{length<util.length>} of the arument and M{h} is the local cell size
1889 jgs 149 from L{getSize<escript.Domain.getSize>}. Where M{|A|==0} M{P} is M{S{infinity}}.
1890    
1891     From the X{Pelclet number} the stabilization parameters M{S{Xi}} and M{S{Xi}} are calculated:
1892    
1893     M{S{Xi}=S{xi}(P) h/|Z|}
1894    
1895     where M{S{xi}} is a suitable function of the Peclet number.
1896    
1897     In the case of a single PDE the coefficient are up-dated in the following way:
1898     - M{A[i,j] S{<-} A[i,j] + S{Xi} * Z[j] * Z[l]}
1899     - M{B[j] S{<-} B[j] + S{Xi} * C[j] * D}
1900     - M{C[j] S{<-} C[j] + S{Xi} * B[j] * D}
1901     - M{X[j] S{<-} X[j] + S{Xi} * Z[j] * Y}
1902    
1903     Similar for the case of a systems of PDEs:
1904     - M{A[i,j,k,l] S{<-} A[i,j,k,l]+ S{delta}[p,m] * S{Xi} * Z[p,i,j] * Z[m,k,l]}
1905     - M{B[i,j,k] S{<-} B[i,j,k] + S{delta}[p,m] * S{Xi} * D[p,k] * C[m,i,j]}
1906     - M{C[i,k,l] S{<-} C[i,k,l] + S{delta}[p,m] * S{Xi} * D[p,k] * B[m,l,i]}
1907     - M{X[i,j] S{<-} X[i,j] + S{delta}[p,m] * S{Xi} * Y[p] * Z[m,i,j]}
1908    
1909     where M{S{delta}} is L{kronecker}.
1910     Using upwinding in this form, introduces an additonal error which is proprtional to the cell size M{h}
1911     but with the intension to stabilize the solution.
1912    
1913     """
1914     def __init__(self,domain,numEquations=None,numSolutions=None,xi=None,debug=False):
1915     """
1916 jgs 151 creates a linear, steady, second order PDE on a L{Domain<escript.Domain>}
1917 jgs 149
1918     @param domain: domain of the PDE
1919     @type domain: L{Domain<escript.Domain>}
1920     @param numEquations: number of equations. If numEquations==None the number of equations
1921     is exracted from the PDE coefficients.
1922     @param numSolutions: number of solution components. If numSolutions==None the number of solution components
1923     is exracted from the PDE coefficients.
1924 jgs 151 @param xi: defines a function which returns for any given Preclet number as L{Scalar<escript.Scalar>} object the
1925 jgs 149 M{S{xi}}-value used to define the stabilization parameters. If equal to None, L{ELMAN_RAMAGE} is used.
1926     @type xi: callable object which returns a L{Scalar<escript.Scalar>} object.
1927     @param debug: if True debug informations are printed.
1928     """
1929 jgs 151 super(AdvectivePDE, self).__init__(domain,\
1930     numEquations,numSolutions,debug)
1931 jgs 148 if xi==None:
1932     self.__xi=AdvectivePDE.ELMAN_RAMAGE
1933     else:
1934     self.__xi=xi
1935 jgs 110 self.__Xi=escript.Data()
1936 jgs 102
1937 jgs 151 def setValue(**coefficients):
1938 jgs 149 """
1939     sets new values to coefficients
1940 jgs 102
1941 jgs 149 @param coefficients: new values assigned to coefficients
1942     @keyword A: value for coefficient A.
1943     @type A: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1944     @keyword B: value for coefficient B
1945     @type B: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1946     @keyword C: value for coefficient C
1947     @type C: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1948     @keyword D: value for coefficient D
1949     @type D: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1950     @keyword X: value for coefficient X
1951     @type X: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1952     @keyword Y: value for coefficient Y
1953     @type Y: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1954     @keyword d: value for coefficient d
1955     @type d: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1956     @keyword y: value for coefficient y
1957     @type y: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1958     @keyword d_contact: value for coefficient d_contact
1959     @type d_contact: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnContactOne<escript.FunctionOnContactOne>}.
1960     or L{FunctionOnContactZero<escript.FunctionOnContactZero>}.
1961     @keyword y_contact: value for coefficient y_contact
1962     @type y_contact: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnContactOne<escript.FunctionOnContactOne>}.
1963     or L{FunctionOnContactZero<escript.FunctionOnContactZero>}.
1964     @keyword r: values prescribed to the solution at the locations of constraints
1965     @type r: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1966     depending of reduced order is used for the solution.
1967     @keyword q: mask for location of constraints
1968     @type q: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1969     depending of reduced order is used for the representation of the equation.
1970     @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1971 jgs 148
1972 jgs 149 """
1973     if "A" in coefficients.keys() or "B" in coefficients.keys() or "C" in coefficients.keys(): self.__Xi=escript.Data()
1974 jgs 151 super(AdvectivePDE, self).setValue(**coefficients)
1975 jgs 149
1976     def ELMAN_RAMAGE(self,P):
1977 jgs 151 """
1978     Predefined function to set a values for M{S{xi}} from a Preclet number M{P}.
1979 jgs 149 This function uses the method suggested by H.C. Elman and A. Ramage, I{SIAM J. Numer. Anal.}, B{40} (2002)
1980     - M{S{xi}(P)=0} for M{P<1}
1981     - M{S{xi}(P)=(1-1/P)/2} otherwise
1982    
1983 jgs 151 @param P: Preclet number
1984 jgs 149 @type P: L{Scalar<escript.Scalar>}
1985     @return: up-wind weightimg factor
1986     @rtype: L{Scalar<escript.Scalar>}
1987 jgs 151 """
1988 gross 304 return util.wherePositive(P-1.)*0.5*(1.-1./(P+1.e-15))
1989 jgs 149
1990     def SIMPLIFIED_BROOK_HUGHES(self,P):
1991 jgs 151 """
1992     Predefined function to set a values for M{S{xi}} from a Preclet number M{P}.
1993 jgs 149 The original methods is
1994 jgs 151
1995 jgs 149 M{S{xi}(P)=coth(P)-1/P}
1996    
1997     As the evaluation of M{coth} is expensive we are using the approximation:
1998 jgs 151
1999 jgs 149 - M{S{xi}(P)=P/3} where M{P<3}
2000     - M{S{xi}(P)=1/2} otherwise
2001    
2002 jgs 151 @param P: Preclet number
2003 jgs 149 @type P: L{Scalar<escript.Scalar>}
2004     @return: up-wind weightimg factor
2005     @rtype: L{Scalar<escript.Scalar>}
2006 jgs 151 """
2007 gross 304 c=util.whereNegative(P-3.)
2008 jgs 148 return P/6.*c+1./2.*(1.-c)
2009    
2010 jgs 149 def HALF(self,P):
2011 jgs 151 """
2012     Predefined function to set value M{1/2} for M{S{xi}}
2013    
2014     @param P: Preclet number
2015 jgs 149 @type P: L{Scalar<escript.Scalar>}
2016     @return: up-wind weightimg factor
2017     @rtype: L{Scalar<escript.Scalar>}
2018 jgs 151 """
2019 jgs 149 return escript.Scalar(0.5,P.getFunctionSpace())
2020 jgs 148
2021 gross 328 def __calculateXi(self,peclet_factor,flux,h):
2022     flux=util.Lsup(flux)
2023     if flux_max>0.:
2024     return h*self.__xi(flux*peclet_factor)/(flux+flux_max*self.__TOL)
2025 jgs 149 else:
2026     return 0.
2027    
2028     def __getXi(self):
2029 jgs 110 if self.__Xi.isEmpty():
2030     B=self.getCoefficient("B")
2031     C=self.getCoefficient("C")
2032     A=self.getCoefficient("A")
2033     h=self.getDomain().getSize()
2034     self.__Xi=escript.Scalar(0.,self.getFunctionSpaceForCoefficient("A"))
2035     if not C.isEmpty() or not B.isEmpty():
2036     if not C.isEmpty() and not B.isEmpty():
2037 gross 328 flux2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))
2038 jgs 110 if self.getNumEquations()>1:
2039     if self.getNumSolutions()>1:
2040     for i in range(self.getNumEquations()):
2041     for k in range(self.getNumSolutions()):
2042 gross 328 for l in range(self.getDim()): flux2+=(C[i,k,l]-B[i,l,k])**2
2043     # flux=C-util.reorderComponents(B,[0,2,1])
2044 jgs 110 else:
2045     for i in range(self.getNumEquations()):
2046 gross 328 for l in range(self.getDim()): flux2+=(C[i,l]-B[i,l])**2
2047     # flux=C-B
2048 jgs 110 else:
2049     if self.getNumSolutions()>1:
2050     for k in range(self.getNumSolutions()):
2051 gross 328 for l in range(self.getDim()): flux2+=(C[k,l]-B[l,k])**2
2052     # flux=C-util.reorderComponents(B,[1,0])
2053 jgs 110 else:
2054 gross 328 for l in range(self.getDim()): flux2+=(C[l]-B[l])**2
2055     #flux=C-B
2056     length_of_flux=util.sqrt(flux2)
2057 jgs 110 elif C.isEmpty():
2058 gross 328 length_of_flux=util.length(B)
2059     #flux=B
2060 jgs 110 else:
2061 gross 328 length_of_flux=util.length(C)
2062     #flux=C
2063 jgs 102
2064 gross 328 #length_of_flux=util.length(flux)
2065     flux_max=util.Lsup(length_of_flux)
2066     if flux_max>0.:
2067     # length_of_A=util.inner(flux,util.tensormutiply(A,flux))
2068 jgs 110 length_of_A=util.length(A)
2069     A_max=util.Lsup(length_of_A)
2070     if A_max>0:
2071 jgs 149 inv_A=1./(length_of_A+A_max*self.__TOL)
2072 jgs 110 else:
2073 jgs 149 inv_A=1./self.__TOL
2074 gross 328 peclet_number=length_of_flux*h/2*inv_A
2075 jgs 110 xi=self.__xi(peclet_number)
2076 gross 328 self.__Xi=h*xi/(length_of_flux+flux_max*self.__TOL)
2077 jgs 149 self.trace("preclet number = %e"%util.Lsup(peclet_number))
2078 jgs 110 return self.__Xi
2079 jgs 102
2080 jgs 148
2081     def getCoefficientOfGeneralPDE(self,name):
2082 jgs 108 """
2083 jgs 122 return the value of the coefficient name of the general PDE
2084    
2085 jgs 149 @param name: name of the coefficient requested.
2086     @type name: C{string}
2087     @return: the value of the coefficient name
2088     @rtype: L{Data<escript.Data>}
2089     @raise IllegalCoefficient: if name is not one of coefficients
2090     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}.
2091     @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
2092 jgs 108 """
2093 jgs 110 if not self.getNumEquations() == self.getNumSolutions():
2094     raise ValueError,"AdvectivePDE expects the number of solution componets and the number of equations to be equal."
2095    
2096 jgs 148 if name == "A" :
2097 jgs 108 A=self.getCoefficient("A")
2098     B=self.getCoefficient("B")
2099     C=self.getCoefficient("C")
2100 jgs 110 if B.isEmpty() and C.isEmpty():
2101     Aout=A
2102     else:
2103 jgs 148 if A.isEmpty():
2104 jgs 110 Aout=self.createNewCoefficient("A")
2105     else:
2106     Aout=A[:]
2107 jgs 149 Xi=self.__getXi()
2108 jgs 110 if self.getNumEquations()>1:
2109     for i in range(self.getNumEquations()):
2110     for j in range(self.getDim()):
2111 jgs 108 for k in range(self.getNumSolutions()):
2112 jgs 110 for l in range(self.getDim()):
2113     if not C.isEmpty() and not B.isEmpty():
2114 gross 328 # tmp=C-util.reorderComponents(B,[0,2,1])
2115     # Aout=Aout+Xi*util.generalTensorProduct(util.reorder(tmp,[1,2,0]),tmp,offset=1)
2116 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])
2117     elif C.isEmpty():
2118     for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*B[p,j,i]*B[p,l,k]
2119 gross 328 # Aout=Aout+Xi*util.generalTensorProduct(util.reorder(B,[2,1,0]),util.reorder(B,[0,2,1]),offset=1)
2120 jgs 110 else:
2121     for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*C[p,i,j]*C[p,k,l]
2122 gross 328 # Aout=Aout+Xi*util.generalTensorProduct(util.reorder(C,[1,2,0]),C,offset=1)
2123 jgs 110 else:
2124     for j in range(self.getDim()):
2125     for l in range(self.getDim()):
2126     if not C.isEmpty() and not B.isEmpty():
2127     Aout[j,l]+=Xi*(C[j]-B[j])*(C[l]-B[l])
2128     elif C.isEmpty():
2129     Aout[j,l]+=Xi*B[j]*B[l]
2130     else:
2131     Aout[j,l]+=Xi*C[j]*C[l]
2132 gross 328 # if not C.isEmpty() and not B.isEmpty():
2133     # tmp=C-B
2134     # Aout=Aout+Xi*util.outer(tmp,tmp)
2135     # elif C.isEmpty():
2136     # Aout=Aout+Xi*util.outer(B,B)
2137     # else:
2138     # Aout=Aout+Xi*util.outer(C,C)
2139 jgs 110 return Aout
2140 jgs 148 elif name == "B" :
2141 jgs 110 B=self.getCoefficient("B")
2142     C=self.getCoefficient("C")
2143     D=self.getCoefficient("D")
2144     if C.isEmpty() or D.isEmpty():
2145     Bout=B
2146     else:
2147 jgs 149 Xi=self.__getXi()
2148 jgs 148 if B.isEmpty():
2149 jgs 110 Bout=self.createNewCoefficient("B")
2150     else:
2151     Bout=B[:]
2152     if self.getNumEquations()>1:
2153     for k in range(self.getNumSolutions()):
2154 jgs 148 for p in range(self.getNumEquations()):
2155 jgs 110 tmp=Xi*D[p,k]
2156     for i in range(self.getNumEquations()):
2157     for j in range(self.getDim()):
2158     Bout[i,j,k]+=tmp*C[p,i,j]
2159 gross 328 # Bout=Bout+Xi*util.generalTensorProduct(util.reorder(C,[1,2,0]),D,offset=1)
2160 jgs 110 else:
2161     tmp=Xi*D
2162     for j in range(self.getDim()): Bout[j]+=tmp*C[j]
2163 gross 328 # Bout=Bout+Xi*D*C
2164 jgs 110 return Bout
2165 jgs 148 elif name == "C" :
2166 jgs 110 B=self.getCoefficient("B")
2167     C=self.getCoefficient("C")
2168     D=self.getCoefficient("D")
2169     if B.isEmpty() or D.isEmpty():
2170     Cout=C
2171     else:
2172 jgs 149 Xi=self.__getXi()
2173 jgs 148 if C.isEmpty():
2174 jgs 110 Cout=self.createNewCoefficient("C")
2175     else:
2176     Cout=C[:]
2177     if self.getNumEquations()>1:
2178     for k in range(self.getNumSolutions()):
2179     for p in range(self.getNumEquations()):
2180     tmp=Xi*D[p,k]
2181     for i in range(self.getNumEquations()):
2182     for l in range(self.getDim()):
2183     Cout[i,k,l]+=tmp*B[p,l,i]
2184 gross 328 # Cout=Cout+Xi*B[p,l,i]*D[p,k]
2185 jgs 110 else:
2186     tmp=Xi*D
2187     for j in range(self.getDim()): Cout[j]+=tmp*B[j]
2188 gross 328 # Cout=Cout+tmp*D*B
2189 jgs 110 return Cout
2190 jgs 148 elif name == "D" :
2191 jgs 108 return self.getCoefficient("D")
2192 jgs 148 elif name == "X" :
2193 jgs 110 X=self.getCoefficient("X")
2194     Y=self.getCoefficient("Y")
2195     B=self.getCoefficient("B")
2196     C=self.getCoefficient("C")
2197     if Y.isEmpty() or (B.isEmpty() and C.isEmpty()):
2198     Xout=X
2199     else:
2200     if X.isEmpty():
2201     Xout=self.createNewCoefficient("X")
2202     else:
2203     Xout=X[:]
2204 jgs 149 Xi=self.__getXi()
2205 jgs 110 if self.getNumEquations()>1:
2206 jgs 148 for p in range(self.getNumEquations()):
2207 jgs 110 tmp=Xi*Y[p]
2208     for i in range(self.getNumEquations()):
2209     for j in range(self.getDim()):
2210     if not C.isEmpty() and not B.isEmpty():
2211     Xout[i,j]+=tmp*(C[p,i,j]-B[p,j,i])
2212 gross 328 # Xout=X_out+Xi*util.inner(Y,C-util.reorderComponents(B,[0,2,1]),offset=1)
2213 jgs 110 elif C.isEmpty():
2214     Xout[i,j]-=tmp*B[p,j,i]
2215 gross 328 # Xout=X_out-Xi*util.inner(Y,util.reorderComponents(B,[0,2,1]),offset=1)
2216 jgs 110 else:
2217     Xout[i,j]+=tmp*C[p,i,j]
2218 gross 328 # Xout=X_out+Xi*util.inner(Y,C,offset=1)
2219 jgs 110 else:
2220     tmp=Xi*Y
2221     for j in range(self.getDim()):
2222     if not C.isEmpty() and not B.isEmpty():
2223     Xout[j]+=tmp*(C[j]-B[j])
2224 gross 328 # Xout=Xout+Xi*Y*(C-B)
2225 jgs 110 elif C.isEmpty():
2226     Xout[j]-=tmp*B[j]
2227 gross 328 # Xout=Xout-Xi*Y*B
2228 jgs 110 else:
2229     Xout[j]+=tmp*C[j]
2230 gross 328 # Xout=Xout+Xi*Y*C
2231 jgs 110 return Xout
2232 jgs 148 elif name == "Y" :
2233 jgs 108 return self.getCoefficient("Y")
2234 jgs 148 elif name == "d" :
2235 jgs 108 return self.getCoefficient("d")
2236 jgs 148 elif name == "y" :
2237 jgs 108 return self.getCoefficient("y")
2238 jgs 148 elif name == "d_contact" :
2239 jgs 108 return self.getCoefficient("d_contact")
2240     elif name == "y_contact" :
2241     return self.getCoefficient("y_contact")
2242 jgs 148 elif name == "r" :
2243 jgs 108 return self.getCoefficient("r")
2244 jgs 148 elif name == "q" :
2245 jgs 108 return self.getCoefficient("q")
2246     else:
2247 jgs 149 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
2248 jgs 108
2249 gross 345 class AdvectionDiffusion(LinearPDE):
2250     """
2251     Class to define PDE equation of the unisotropic advection-diffusion problem, which is genear L{LinearPDE} of the form
2252 jgs 108
2253 gross 345 M{S{omega}*u + inner(v,grad(u))- grad(matrixmult(k,grad(u))[j])[j] = f}
2254    
2255     with natural boundary conditons
2256    
2257     M{n[j]*matrixmult(k,grad(u))[j] = g- S{alpha}u }
2258    
2259     and constraints:
2260    
2261     M{u=r} where M{q>0}
2262    
2263     @remark: no upwinding is applied yet.
2264    
2265     """
2266    
2267     def __init__(self,domain,debug=False):
2268     """
2269     initializes a new Poisson equation
2270    
2271     @param domain: domain of the PDE
2272     @type domain: L{Domain<escript.Domain>}
2273     @param debug: if True debug informations are printed.
2274    
2275     """
2276     super(Helmholtz, self).__init__(domain,1,1,debug)
2277     self.COEFFICIENTS={"omega": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
2278     "k": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_DIM,PDECoefficient.BY_DIM),PDECoefficient.OPERATOR),
2279     "f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
2280     "v": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_DIM,),PDECoefficient.OPERATOR),
2281     "alpha": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
2282     "g": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
2283     "r": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH),
2284     "q": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
2285    
2286     def setValue(self,**coefficients):
2287     """
2288     sets new values to coefficients
2289    
2290     @param coefficients: new values assigned to coefficients
2291     @keyword omega: value for coefficient M{S{omega}}
2292     @type omega: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2293     @keyword k: value for coefficient M{k}
2294     @type k: any type that can be casted to L{Tensor<escript.Tensor>} object on L{Function<escript.Function>}.
2295     @keyword v: value for coefficient M{v}
2296     @type v: any type that can be casted to L{Vector<escript.Vector>} object on L{Function<escript.Function>}.
2297     @keyword f: value for right hand side M{f}
2298     @type f: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2299     @keyword alpha: value for right hand side M{S{alpha}}
2300     @type alpha: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
2301     @keyword g: value for right hand side M{g}
2302     @type g: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
2303     @keyword r: prescribed values M{r} for the solution in constraints.
2304     @type r: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2305     depending of reduced order is used for the representation of the equation.
2306     @keyword q: mask for location of constraints
2307     @type q: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2308     depending of reduced order is used for the representation of the equation.
2309     @raise IllegalCoefficient: if an unknown coefficient keyword is used.
2310     """
2311     super(Helmholtz, self).setValue(**coefficients)
2312    
2313     def getCoefficientOfGeneralPDE(self,name):
2314     """
2315     return the value of the coefficient name of the general PDE
2316    
2317     @param name: name of the coefficient requested.
2318     @type name: C{string}
2319     @return: the value of the coefficient name
2320     @rtype: L{Data<escript.Data>}
2321     @raise IllegalCoefficient: if name is not one of coefficients
2322     "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}.
2323     @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
2324     """
2325     if name == "A" :
2326     return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))*self.getCoefficient("k")
2327     elif name == "B" :
2328     return escript.Data()
2329     elif name == "C" :
2330     return escript.getCoefficient("v")
2331     elif name == "D" :
2332     return self.getCoefficient("omega")
2333     elif name == "X" :
2334     return escript.Data()
2335     elif name == "Y" :
2336     return self.getCoefficient("f")
2337     elif name == "d" :
2338     return self.getCoefficient("alpha")
2339     elif name == "y" :
2340     return self.getCoefficient("g")
2341