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