/[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 596 - (hide annotations)
Wed Mar 15 22:41:23 2006 UTC (13 years, 6 months ago) by gross
File MIME type: text/x-python
File size: 105871 byte(s)
some work on AdvectivePDE for the scalar case
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 gross 596 SMALL_TOLERANCE=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 gross 596 tol=util.Lsup(A)*self.SMALL_TOLERANCE
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 gross 596 tol=(util.Lsup(B)+util.Lsup(C))*self.SMALL_TOLERANCE/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 gross 596 tol=util.Lsup(D)*self.SMALL_TOLERANCE
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 gross 596 tol=util.Lsup(d)*self.SMALL_TOLERANCE
755 jgs 149 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 gross 596 tol=util.Lsup(d_contact)*self.SMALL_TOLERANCE
763 jgs 149 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 gross 431 @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}, L{AMG}
829 gross 387 @param preconditioner: sets a new solver method.
830 gross 431 @type preconditioner: one of L{DEFAULT}, L{JACOBI} L{ILU0}, L{ILUT},L{SSOR}, L{RILU}
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 431 elif m[0]==self.AMG: method= "AMG"
864 gross 425 if m[1]==self.DEFAULT: method+="+DEFAULT"
865     elif m[1]==self.JACOBI: method+= "+JACOBI"
866     elif m[1]==self.ILU0: method+= "+ILU0"
867     elif m[1]==self.ILUT: method+= "+ILUT"
868     elif m[1]==self.SSOR: method+= "+SSOR"
869 gross 431 elif m[1]==self.AMG: method+= "+AMG"
870     elif m[1]==self.RILU: method+= "+RILU"
871 jgs 150 if p==self.DEFAULT: package="DEFAULT"
872     elif p==self.PASO: package= "PASO"
873     elif p==self.MKL: package= "MKL"
874     elif p==self.SCSL: package= "SCSL"
875     elif p==self.UMFPACK: package= "UMFPACK"
876     else : method="unknown"
877     return "%s solver of %s package"%(method,package)
878 jgs 148
879 jgs 149
880 jgs 102 def getSolverMethod(self):
881     """
882 jgs 122 returns the solver method
883 jgs 149
884 jgs 151 @return: the solver method currently be used.
885 jgs 149 @rtype: C{int}
886 jgs 102 """
887 gross 387 return self.__solver_method,self.__preconditioner
888 jgs 102
889 jgs 150 def setSolverPackage(self,package=None):
890     """
891     sets a new solver package
892    
893     @param solver: sets a new solver method.
894     @type solver: one of L{DEFAULT}, L{PASO} L{SCSL}, L{MKL}, L{UMLPACK}
895     """
896     if package==None: package=self.DEFAULT
897     if not package==self.getSolverPackage():
898     self.__solver_method=solver
899     self.__checkMatrixType()
900     self.trace("New solver is %s"%self.getSolverMethodName())
901    
902     def getSolverPackage(self):
903     """
904     returns the package of the solver
905    
906 jgs 151 @return: the solver package currently being used.
907 jgs 150 @rtype: C{int}
908     """
909     return self.__solver_package
910    
911 jgs 148 def isUsingLumping(self):
912     """
913     checks if matrix lumping is used a solver method
914    
915 jgs 149 @return: True is lumping is currently used a solver method.
916 jgs 148 @rtype: C{bool}
917     """
918 gross 387 return self.getSolverMethod()[0]==self.LUMPING
919 jgs 148
920 jgs 102 def setTolerance(self,tol=1.e-8):
921     """
922 jgs 149 resets the tolerance for the solver method to tol where for an appropriate norm M{|.|}
923 jgs 148
924 jgs 149 M{|L{getResidual}()|<tol*|L{getRightHandSide}()|}
925 jgs 148
926 jgs 149 defines the stopping criterion.
927 jgs 148
928     @param tol: new tolerance for the solver. If the tol is lower then the current tolerence
929     the system will be resolved.
930 jgs 149 @type tol: positive C{float}
931 jgs 148 @raise ValueException: if tolerance is not positive.
932 jgs 102 """
933     if not tol>0:
934     raise ValueException,"Tolerance as to be positive"
935 jgs 148 if tol<self.getTolerance(): self.__invalidateSolution()
936     self.trace("New tolerance %e"%tol)
937 jgs 102 self.__tolerance=tol
938     return
939 jgs 148
940 jgs 102 def getTolerance(self):
941     """
942 jgs 122 returns the tolerance set for the solution
943 jgs 148
944     @return: tolerance currently used.
945     @rtype: C{float}
946 jgs 102 """
947     return self.__tolerance
948    
949 jgs 148 # =============================================================================
950     # symmetry flag:
951     # =============================================================================
952 jgs 102 def isSymmetric(self):
953     """
954 jgs 148 checks if symmetry is indicated.
955 jgs 149
956     @return: True is a symmetric PDE is indicated, otherwise False is returned
957     @rtype: C{bool}
958 jgs 102 """
959     return self.__sym
960    
961     def setSymmetryOn(self):
962     """
963 jgs 148 sets the symmetry flag.
964 jgs 102 """
965     if not self.isSymmetric():
966 jgs 148 self.trace("PDE is set to be symmetric")
967 jgs 102 self.__sym=True
968     self.__checkMatrixType()
969    
970     def setSymmetryOff(self):
971     """
972 jgs 148 removes the symmetry flag.
973 jgs 102 """
974     if self.isSymmetric():
975 jgs 148 self.trace("PDE is set to be unsymmetric")
976 jgs 102 self.__sym=False
977     self.__checkMatrixType()
978    
979     def setSymmetryTo(self,flag=False):
980 jgs 148 """
981     sets the symmetry flag to flag
982 jgs 149
983 jgs 148 @param flag: If flag, the symmetry flag is set otherwise the symmetry flag is released.
984     @type flag: C{bool}
985     """
986     if flag:
987     self.setSymmetryOn()
988     else:
989     self.setSymmetryOff()
990 jgs 102
991 jgs 148 # =============================================================================
992     # function space handling for the equation as well as the solution
993     # =============================================================================
994 jgs 102 def setReducedOrderOn(self):
995     """
996 jgs 148 switches on reduced order for solution and equation representation
997 jgs 149
998     @raise RuntimeError: if order reduction is altered after a coefficient has been set.
999 jgs 102 """
1000     self.setReducedOrderForSolutionOn()
1001     self.setReducedOrderForEquationOn()
1002    
1003     def setReducedOrderOff(self):
1004     """
1005 jgs 148 switches off reduced order for solution and equation representation
1006 jgs 149
1007     @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1008 jgs 102 """
1009     self.setReducedOrderForSolutionOff()
1010     self.setReducedOrderForEquationOff()
1011    
1012     def setReducedOrderTo(self,flag=False):
1013     """
1014 jgs 149 sets order reduction for both solution and equation representation according to flag.
1015     @param flag: if flag is True, the order reduction is switched on for both solution and equation representation, otherwise or
1016 jgs 148 if flag is not present order reduction is switched off
1017     @type flag: C{bool}
1018 jgs 149 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1019 jgs 102 """
1020     self.setReducedOrderForSolutionTo(flag)
1021     self.setReducedOrderForEquationTo(flag)
1022    
1023 jgs 148
1024 jgs 102 def setReducedOrderForSolutionOn(self):
1025     """
1026 jgs 148 switches on reduced order for solution representation
1027 jgs 149
1028     @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1029 jgs 102 """
1030 jgs 149 if not self.__reduce_solution_order:
1031     if self.__altered_coefficients:
1032     raise RuntimeError,"order cannot be altered after coefficients have been defined."
1033 jgs 148 self.trace("Reduced order is used to solution representation.")
1034 jgs 149 self.__reduce_solution_order=True
1035 jgs 148 self.__resetSystem()
1036 jgs 102
1037     def setReducedOrderForSolutionOff(self):
1038     """
1039 jgs 148 switches off reduced order for solution representation
1040 jgs 149
1041     @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1042 jgs 102 """
1043 jgs 149 if self.__reduce_solution_order:
1044     if self.__altered_coefficients:
1045     raise RuntimeError,"order cannot be altered after coefficients have been defined."
1046 jgs 148 self.trace("Full order is used to interpolate solution.")
1047 jgs 149 self.__reduce_solution_order=False
1048 jgs 148 self.__resetSystem()
1049 jgs 102
1050     def setReducedOrderForSolutionTo(self,flag=False):
1051     """
1052 jgs 122 sets order for test functions according to flag
1053 jgs 102
1054 jgs 149 @param flag: if flag is True, the order reduction is switched on for solution representation, otherwise or
1055 jgs 148 if flag is not present order reduction is switched off
1056     @type flag: C{bool}
1057 jgs 149 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1058 jgs 102 """
1059     if flag:
1060     self.setReducedOrderForSolutionOn()
1061     else:
1062     self.setReducedOrderForSolutionOff()
1063 jgs 148
1064 jgs 102 def setReducedOrderForEquationOn(self):
1065     """
1066 jgs 148 switches on reduced order for equation representation
1067 jgs 149
1068     @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1069 jgs 102 """
1070 jgs 149 if not self.__reduce_equation_order:
1071     if self.__altered_coefficients:
1072     raise RuntimeError,"order cannot be altered after coefficients have been defined."
1073 jgs 148 self.trace("Reduced order is used for test functions.")
1074 jgs 149 self.__reduce_equation_order=True
1075 jgs 148 self.__resetSystem()
1076 jgs 102
1077     def setReducedOrderForEquationOff(self):
1078     """
1079 jgs 148 switches off reduced order for equation representation
1080 jgs 149
1081     @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1082 jgs 102 """
1083 jgs 149 if self.__reduce_equation_order:
1084     if self.__altered_coefficients:
1085     raise RuntimeError,"order cannot be altered after coefficients have been defined."
1086 jgs 148 self.trace("Full order is used for test functions.")
1087 jgs 149 self.__reduce_equation_order=False
1088 jgs 148 self.__resetSystem()
1089 jgs 102
1090     def setReducedOrderForEquationTo(self,flag=False):
1091     """
1092 jgs 122 sets order for test functions according to flag
1093 jgs 102
1094 jgs 149 @param flag: if flag is True, the order reduction is switched on for equation representation, otherwise or
1095 jgs 148 if flag is not present order reduction is switched off
1096     @type flag: C{bool}
1097 jgs 149 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1098 jgs 102 """
1099     if flag:
1100     self.setReducedOrderForEquationOn()
1101     else:
1102     self.setReducedOrderForEquationOff()
1103 jgs 148
1104     # =============================================================================
1105     # private method:
1106     # =============================================================================
1107     def __checkMatrixType(self):
1108     """
1109     reassess the matrix type and, if a new matrix is needed, resets the system.
1110     """
1111 gross 387 new_matrix_type=self.getDomain().getSystemMatrixTypeId(self.getSolverMethod()[0],self.getSolverPackage(),self.isSymmetric())
1112 jgs 148 if not new_matrix_type==self.__matrix_type:
1113     self.trace("Matrix type is now %d."%new_matrix_type)
1114     self.__matrix_type=new_matrix_type
1115     self.__resetSystem()
1116 jgs 149 #
1117 jgs 148 # rebuild switches :
1118 jgs 149 #
1119 jgs 148 def __invalidateSolution(self):
1120     """
1121     indicates the PDE has to be resolved if the solution is requested
1122     """
1123     if self.__solution_isValid: self.trace("PDE has to be resolved.")
1124     self.__solution_isValid=False
1125    
1126     def __invalidateOperator(self):
1127     """
1128     indicates the operator has to be rebuilt next time it is used
1129     """
1130 jgs 149 if self.__operator_is_Valid: self.trace("Operator has to be rebuilt.")
1131 jgs 148 self.__invalidateSolution()
1132 jgs 149 self.__operator_is_Valid=False
1133 jgs 148
1134     def __invalidateRightHandSide(self):
1135     """
1136     indicates the right hand side has to be rebuild next time it is used
1137     """
1138     if self.__righthandside_isValid: self.trace("Right hand side has to be rebuilt.")
1139     self.__invalidateSolution()
1140     self.__righthandside_isValid=False
1141    
1142     def __invalidateSystem(self):
1143     """
1144     annonced that everthing has to be rebuild:
1145     """
1146     if self.__righthandside_isValid: self.trace("System has to be rebuilt.")
1147     self.__invalidateSolution()
1148     self.__invalidateOperator()
1149     self.__invalidateRightHandSide()
1150    
1151     def __resetSystem(self):
1152     """
1153     annonced that everthing has to be rebuild:
1154     """
1155     self.trace("New System is built from scratch.")
1156     self.__operator=escript.Operator()
1157 jgs 149 self.__operator_is_Valid=False
1158 jgs 148 self.__righthandside=escript.Data()
1159     self.__righthandside_isValid=False
1160     self.__solution=escript.Data()
1161     self.__solution_isValid=False
1162 jgs 149 #
1163 jgs 148 # system initialization:
1164 jgs 149 #
1165 jgs 121 def __getNewOperator(self):
1166 jgs 102 """
1167 jgs 148 returns an instance of a new operator
1168 jgs 102 """
1169 jgs 148 self.trace("New operator is allocated.")
1170 jgs 102 return self.getDomain().newOperator( \
1171     self.getNumEquations(), \
1172     self.getFunctionSpaceForEquation(), \
1173     self.getNumSolutions(), \
1174     self.getFunctionSpaceForSolution(), \
1175     self.__matrix_type)
1176    
1177 jgs 148 def __getNewRightHandSide(self):
1178 jgs 102 """
1179 jgs 148 returns an instance of a new right hand side
1180 jgs 102 """
1181 jgs 148 self.trace("New right hand side is allocated.")
1182 jgs 121 if self.getNumEquations()>1:
1183 jgs 148 return escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)
1184 jgs 121 else:
1185 jgs 148 return escript.Data(0.,(),self.getFunctionSpaceForEquation(),True)
1186 jgs 102
1187 jgs 121 def __getNewSolution(self):
1188 jgs 102 """
1189 jgs 148 returns an instance of a new solution
1190 jgs 102 """
1191 jgs 148 self.trace("New solution is allocated.")
1192 jgs 121 if self.getNumSolutions()>1:
1193     return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)
1194     else:
1195     return escript.Data(0.,(),self.getFunctionSpaceForSolution(),True)
1196 jgs 102
1197 jgs 148 def __makeFreshSolution(self):
1198     """
1199     makes sure that the solution is instantiated and returns it initialized by zeros
1200     """
1201     if self.__solution.isEmpty():
1202     self.__solution=self.__getNewSolution()
1203     else:
1204     self.__solution*=0
1205     self.trace("Solution is reset to zero.")
1206     return self.__solution
1207    
1208     def __makeFreshRightHandSide(self):
1209     """
1210     makes sure that the right hand side is instantiated and returns it initialized by zeros
1211     """
1212     if self.__righthandside.isEmpty():
1213     self.__righthandside=self.__getNewRightHandSide()
1214     else:
1215     self.__righthandside*=0
1216     self.trace("Right hand side is reset to zero.")
1217     return self.__righthandside
1218    
1219 jgs 121 def __makeFreshOperator(self):
1220 jgs 102 """
1221 jgs 148 makes sure that the operator is instantiated and returns it initialized by zeros
1222 jgs 102 """
1223     if self.__operator.isEmpty():
1224 jgs 121 self.__operator=self.__getNewOperator()
1225 jgs 102 else:
1226 jgs 149 self.__operator.resetValues()
1227 jgs 148 self.trace("Operator reset to zero")
1228 jgs 102 return self.__operator
1229    
1230 jgs 148 def __applyConstraint(self):
1231     """
1232     applies the constraints defined by q and r to the system
1233     """
1234     if not self.isUsingLumping():
1235     q=self.getCoefficientOfGeneralPDE("q")
1236     r=self.getCoefficientOfGeneralPDE("r")
1237     if not q.isEmpty() and not self.__operator.isEmpty():
1238     # q is the row and column mask to indicate where constraints are set:
1239     row_q=escript.Data(q,self.getFunctionSpaceForEquation())
1240     col_q=escript.Data(q,self.getFunctionSpaceForSolution())
1241     u=self.__getNewSolution()
1242     if r.isEmpty():
1243     r_s=self.__getNewSolution()
1244     else:
1245     r_s=escript.Data(r,self.getFunctionSpaceForSolution())
1246     u.copyWithMask(r_s,col_q)
1247 jgs 149 if not self.__righthandside.isEmpty():
1248 jgs 148 self.__righthandside-=self.__operator*u
1249     self.__righthandside=self.copyConstraint(self.__righthandside)
1250     self.__operator.nullifyRowsAndCols(row_q,col_q,1.)
1251     # =============================================================================
1252     # function giving access to coefficients of the general PDE:
1253     # =============================================================================
1254     def getCoefficientOfGeneralPDE(self,name):
1255 jgs 108 """
1256 jgs 148 return the value of the coefficient name of the general PDE.
1257    
1258 jgs 149 @note: This method is called by the assembling routine it can be overwritten
1259 jgs 148 to map coefficients of a particular PDE to the general PDE.
1260 jgs 149 @param name: name of the coefficient requested.
1261 jgs 148 @type name: C{string}
1262 jgs 149 @return: the value of the coefficient name
1263     @rtype: L{Data<escript.Data>}
1264 jgs 148 @raise IllegalCoefficient: if name is not one of coefficients
1265 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}.
1266 jgs 108 """
1267 jgs 148 if self.hasCoefficientOfGeneralPDE(name):
1268     return self.getCoefficient(name)
1269     else:
1270     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1271 jgs 108
1272 jgs 148 def hasCoefficientOfGeneralPDE(self,name):
1273 jgs 108 """
1274 jgs 148 checks if name is a the name of a coefficient of the general PDE.
1275 jgs 149
1276 jgs 148 @param name: name of the coefficient enquired.
1277     @type name: C{string}
1278 jgs 149 @return: True if name is the name of a coefficient of the general PDE. Otherwise False.
1279     @rtype: C{bool}
1280    
1281 jgs 108 """
1282 jgs 148 return self.__COEFFICIENTS_OF_GENEARL_PDE.has_key(name)
1283 jgs 108
1284 jgs 148 def createCoefficientOfGeneralPDE(self,name):
1285 jgs 108 """
1286 jgs 148 returns a new instance of a coefficient for coefficient name of the general PDE
1287    
1288     @param name: name of the coefficient requested.
1289     @type name: C{string}
1290 jgs 149 @return: a coefficient name initialized to 0.
1291     @rtype: L{Data<escript.Data>}
1292 jgs 148 @raise IllegalCoefficient: if name is not one of coefficients
1293 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}.
1294 jgs 108 """
1295 jgs 148 if self.hasCoefficientOfGeneralPDE(name):
1296     return escript.Data(0,self.getShapeOfCoefficientOfGeneralPDE(name),self.getFunctionSpaceForCoefficientOfGeneralPDE(name))
1297 jgs 108 else:
1298 jgs 148 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1299 jgs 108
1300 jgs 148 def getFunctionSpaceForCoefficientOfGeneralPDE(self,name):
1301 jgs 108 """
1302 jgs 149 return the L{FunctionSpace<escript.FunctionSpace>} to be used for coefficient name of the general PDE
1303 jgs 148
1304     @param name: name of the coefficient enquired.
1305     @type name: C{string}
1306 jgs 149 @return: the function space to be used for coefficient name
1307     @rtype: L{FunctionSpace<escript.FunctionSpace>}
1308 jgs 148 @raise IllegalCoefficient: if name is not one of coefficients
1309 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}.
1310 jgs 108 """
1311 jgs 148 if self.hasCoefficientOfGeneralPDE(name):
1312     return self.__COEFFICIENTS_OF_GENEARL_PDE[name].getFunctionSpace(self.getDomain())
1313 jgs 108 else:
1314 jgs 148 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1315 jgs 108
1316 jgs 148 def getShapeOfCoefficientOfGeneralPDE(self,name):
1317     """
1318     return the shape of the coefficient name of the general PDE
1319 jgs 108
1320 jgs 148 @param name: name of the coefficient enquired.
1321     @type name: C{string}
1322 jgs 149 @return: the shape of the coefficient name
1323     @rtype: C{tuple} of C{int}
1324 jgs 148 @raise IllegalCoefficient: if name is not one of coefficients
1325 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}.
1326 jgs 148 """
1327     if self.hasCoefficientOfGeneralPDE(name):
1328     return self.__COEFFICIENTS_OF_GENEARL_PDE[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())
1329     else:
1330     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1331 jgs 108
1332 jgs 148 # =============================================================================
1333     # functions giving access to coefficients of a particular PDE implementation:
1334     # =============================================================================
1335     def getCoefficient(self,name):
1336     """
1337     returns the value of the coefficient name
1338 jgs 108
1339 jgs 149 @param name: name of the coefficient requested.
1340 jgs 148 @type name: C{string}
1341 jgs 149 @return: the value of the coefficient name
1342     @rtype: L{Data<escript.Data>}
1343 jgs 148 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1344     """
1345     if self.hasCoefficient(name):
1346     return self.COEFFICIENTS[name].getValue()
1347     else:
1348     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1349 jgs 108
1350 jgs 148 def hasCoefficient(self,name):
1351     """
1352     return True if name is the name of a coefficient
1353 jgs 108
1354 jgs 148 @param name: name of the coefficient enquired.
1355     @type name: C{string}
1356 jgs 149 @return: True if name is the name of a coefficient of the general PDE. Otherwise False.
1357     @rtype: C{bool}
1358 jgs 148 """
1359     return self.COEFFICIENTS.has_key(name)
1360 jgs 108
1361 jgs 148 def createCoefficient(self, name):
1362     """
1363 jgs 149 create a L{Data<escript.Data>} object corresponding to coefficient name
1364 jgs 108
1365 jgs 149 @return: a coefficient name initialized to 0.
1366     @rtype: L{Data<escript.Data>}
1367 jgs 148 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1368     """
1369     if self.hasCoefficient(name):
1370     return escript.Data(0.,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))
1371     else:
1372     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1373    
1374     def getFunctionSpaceForCoefficient(self,name):
1375     """
1376 jgs 149 return the L{FunctionSpace<escript.FunctionSpace>} to be used for coefficient name
1377 jgs 148
1378     @param name: name of the coefficient enquired.
1379     @type name: C{string}
1380 jgs 149 @return: the function space to be used for coefficient name
1381     @rtype: L{FunctionSpace<escript.FunctionSpace>}
1382 jgs 148 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1383     """
1384     if self.hasCoefficient(name):
1385     return self.COEFFICIENTS[name].getFunctionSpace(self.getDomain())
1386     else:
1387     raise ValueError,"unknown coefficient %s requested"%name
1388     def getShapeOfCoefficient(self,name):
1389     """
1390 jgs 149 return the shape of the coefficient name
1391 jgs 148
1392     @param name: name of the coefficient enquired.
1393     @type name: C{string}
1394 jgs 149 @return: the shape of the coefficient name
1395     @rtype: C{tuple} of C{int}
1396 jgs 148 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1397     """
1398     if self.hasCoefficient(name):
1399     return self.COEFFICIENTS[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())
1400     else:
1401     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1402    
1403     def resetCoefficients(self):
1404     """
1405     resets all coefficients to there default values.
1406     """
1407     for i in self.COEFFICIENTS.iterkeys():
1408     self.COEFFICIENTS[i].resetValue()
1409    
1410     def alteredCoefficient(self,name):
1411     """
1412     announce that coefficient name has been changed
1413    
1414     @param name: name of the coefficient enquired.
1415     @type name: C{string}
1416     @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1417 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.
1418 jgs 148 """
1419     if self.hasCoefficient(name):
1420     self.trace("Coefficient %s has been altered."%name)
1421 jgs 149 if not ((name=="q" or name=="r") and self.isUsingLumping()):
1422     if self.COEFFICIENTS[name].isAlteringOperator(): self.__invalidateOperator()
1423     if self.COEFFICIENTS[name].isAlteringRightHandSide(): self.__invalidateRightHandSide()
1424 jgs 148 else:
1425     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1426    
1427     def copyConstraint(self,u):
1428 jgs 108 """
1429 jgs 149 copies the constraint into u and returns u.
1430    
1431 jgs 148 @param u: a function of rank 0 is a single PDE is solved and of shape (numSolution,) for a system of PDEs
1432 jgs 149 @type u: L{Data<escript.Data>}
1433     @return: the input u modified by the constraints.
1434     @rtype: L{Data<escript.Data>}
1435     @warning: u is altered if it has the appropriate L{FunctionSpace<escript.FunctionSpace>}
1436 jgs 108 """
1437 jgs 148 q=self.getCoefficientOfGeneralPDE("q")
1438     r=self.getCoefficientOfGeneralPDE("r")
1439     if not q.isEmpty():
1440     if u.isEmpty(): u=escript.Data(0.,q.getShape(),q.getFunctionSpace())
1441     if r.isEmpty():
1442     r=escript.Data(0,u.getShape(),u.getFunctionSpace())
1443     else:
1444     r=escript.Data(r,u.getFunctionSpace())
1445     u.copyWithMask(r,escript.Data(q,u.getFunctionSpace()))
1446     return u
1447    
1448     def setValue(self,**coefficients):
1449     """
1450     sets new values to coefficients
1451    
1452 jgs 149 @param coefficients: new values assigned to coefficients
1453     @keyword A: value for coefficient A.
1454     @type A: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1455 jgs 148 @keyword B: value for coefficient B
1456 jgs 149 @type B: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1457 jgs 148 @keyword C: value for coefficient C
1458 jgs 149 @type C: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1459 jgs 148 @keyword D: value for coefficient D
1460 jgs 149 @type D: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1461 jgs 148 @keyword X: value for coefficient X
1462 jgs 149 @type X: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1463 jgs 148 @keyword Y: value for coefficient Y
1464 jgs 149 @type Y: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1465 jgs 148 @keyword d: value for coefficient d
1466 jgs 149 @type d: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1467 jgs 148 @keyword y: value for coefficient y
1468 jgs 149 @type y: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1469 jgs 148 @keyword d_contact: value for coefficient d_contact
1470 jgs 149 @type d_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 y_contact: value for coefficient y_contact
1473 jgs 149 @type y_contact: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnContactOne<escript.FunctionOnContactOne>}.
1474     or L{FunctionOnContactZero<escript.FunctionOnContactZero>}.
1475 jgs 148 @keyword r: values prescribed to the solution at the locations of constraints
1476 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>}
1477 jgs 148 depending of reduced order is used for the solution.
1478     @keyword q: mask for location of constraints
1479 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>}
1480 jgs 148 depending of reduced order is used for the representation of the equation.
1481     @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1482     """
1483 jgs 108 # check if the coefficients are legal:
1484     for i in coefficients.iterkeys():
1485     if not self.hasCoefficient(i):
1486 jgs 148 raise IllegalCoefficient,"Attempt to set unknown coefficient %s"%i
1487 jgs 108 # if the number of unknowns or equations is still unknown we try to estimate them:
1488 jgs 148 if self.__numEquations==None or self.__numSolutions==None:
1489 jgs 108 for i,d in coefficients.iteritems():
1490     if hasattr(d,"shape"):
1491     s=d.shape
1492     elif hasattr(d,"getShape"):
1493     s=d.getShape()
1494     else:
1495     s=numarray.array(d).shape
1496     if s!=None:
1497     # get number of equations and number of unknowns:
1498 jgs 148 res=self.COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(self.getDomain(),s)
1499 jgs 108 if res==None:
1500 jgs 148 raise IllegalCoefficientValue,"Illegal shape %s of coefficient %s"%(s,i)
1501 jgs 108 else:
1502 jgs 148 if self.__numEquations==None: self.__numEquations=res[0]
1503     if self.__numSolutions==None: self.__numSolutions=res[1]
1504     if self.__numEquations==None: raise UndefinedPDEError,"unidententified number of equations"
1505     if self.__numSolutions==None: raise UndefinedPDEError,"unidententified number of solutions"
1506 jgs 108 # now we check the shape of the coefficient if numEquations and numSolutions are set:
1507     for i,d in coefficients.iteritems():
1508 jgs 148 try:
1509 jgs 149 self.COEFFICIENTS[i].setValue(self.getDomain(),self.getNumEquations(),self.getNumSolutions(),self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
1510 jgs 148 except IllegalCoefficientValue,m:
1511     raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
1512 jgs 108 self.alteredCoefficient(i)
1513    
1514 jgs 149 self.__altered_coefficients=True
1515 jgs 148 # check if the systrem is inhomogeneous:
1516     if len(coefficients)>0 and not self.isUsingLumping():
1517     q=self.getCoefficientOfGeneralPDE("q")
1518     r=self.getCoefficientOfGeneralPDE("r")
1519     homogeneous_constraint=True
1520     if not q.isEmpty() and not r.isEmpty():
1521 jgs 149 if util.Lsup(q*r)>=1.e-13*util.Lsup(r):
1522 jgs 148 self.trace("Inhomogeneous constraint detected.")
1523     self.__invalidateSystem()
1524 jgs 108
1525     def getSystem(self):
1526     """
1527 jgs 122 return the operator and right hand side of the PDE
1528 jgs 149
1529     @return: the discrete version of the PDE
1530     @rtype: C{tuple} of L{Operator,<escript.Operator>} and L{Data<escript.Data>}.
1531 jgs 108 """
1532 jgs 149 if not self.__operator_is_Valid or not self.__righthandside_isValid:
1533 jgs 108 if self.isUsingLumping():
1534 jgs 149 if not self.__operator_is_Valid:
1535 gross 531 if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():
1536     raise TypeError,"Lumped matrix requires same order for equations and unknowns"
1537     if not self.getCoefficientOfGeneralPDE("A").isEmpty():
1538     raise ValueError,"coefficient A in lumped matrix may not be present."
1539     if not self.getCoefficientOfGeneralPDE("B").isEmpty():
1540     raise ValueError,"coefficient A in lumped matrix may not be present."
1541     if not self.getCoefficientOfGeneralPDE("C").isEmpty():
1542     raise ValueError,"coefficient A in lumped matrix may not be present."
1543     D=self.getCoefficientOfGeneralPDE("D")
1544     if not D.isEmpty():
1545     if self.getNumSolutions()>1:
1546     D_times_e=util.matrixmult(D,numarray.ones((self.getNumSolutions(),)))
1547     else:
1548     D_times_e=D
1549     else:
1550     D_times_e=escript.Data()
1551     d=self.getCoefficientOfGeneralPDE("d")
1552     if not d.isEmpty():
1553     if self.getNumSolutions()>1:
1554     d_times_e=util.matrixmult(d,numarray.ones((self.getNumSolutions(),)))
1555     else:
1556     d_times_e=d
1557     else:
1558     d_times_e=escript.Data()
1559     d_contact=self.getCoefficientOfGeneralPDE("d_contact")
1560     if not d_contact.isEmpty():
1561     if self.getNumSolutions()>1:
1562     d_contact_times_e=util.matrixmult(d_contact,numarray.ones((self.getNumSolutions(),)))
1563     else:
1564     d_contact_times_e=d_contact
1565     else:
1566     d_contact_times_e=escript.Data()
1567    
1568     self.__operator=self.__getNewRightHandSide()
1569     self.getDomain().addPDEToRHS(self.__operator, \
1570     escript.Data(), \
1571     D_times_e, \
1572     d_times_e,\
1573     d_contact_times_e)
1574     self.__operator=1./self.__operator
1575 jgs 148 self.trace("New lumped operator has been built.")
1576 jgs 149 self.__operator_is_Valid=True
1577 jgs 108 if not self.__righthandside_isValid:
1578 jgs 121 self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \
1579 jgs 148 self.getCoefficientOfGeneralPDE("X"), \
1580     self.getCoefficientOfGeneralPDE("Y"),\
1581     self.getCoefficientOfGeneralPDE("y"),\
1582     self.getCoefficientOfGeneralPDE("y_contact"))
1583     self.trace("New right hand side as been built.")
1584 jgs 108 self.__righthandside_isValid=True
1585     else:
1586 jgs 149 if not self.__operator_is_Valid and not self.__righthandside_isValid:
1587 jgs 121 self.getDomain().addPDEToSystem(self.__makeFreshOperator(),self.__makeFreshRightHandSide(), \
1588 jgs 148 self.getCoefficientOfGeneralPDE("A"), \
1589     self.getCoefficientOfGeneralPDE("B"), \
1590     self.getCoefficientOfGeneralPDE("C"), \
1591     self.getCoefficientOfGeneralPDE("D"), \
1592     self.getCoefficientOfGeneralPDE("X"), \
1593     self.getCoefficientOfGeneralPDE("Y"), \
1594     self.getCoefficientOfGeneralPDE("d"), \
1595     self.getCoefficientOfGeneralPDE("y"), \
1596     self.getCoefficientOfGeneralPDE("d_contact"), \
1597     self.getCoefficientOfGeneralPDE("y_contact"))
1598 jgs 108 self.__applyConstraint()
1599 jgs 148 self.__righthandside=self.copyConstraint(self.__righthandside)
1600     self.trace("New system has been built.")
1601 jgs 149 self.__operator_is_Valid=True
1602 jgs 108 self.__righthandside_isValid=True
1603     elif not self.__righthandside_isValid:
1604 jgs 121 self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \
1605 jgs 148 self.getCoefficientOfGeneralPDE("X"), \
1606     self.getCoefficientOfGeneralPDE("Y"),\
1607     self.getCoefficientOfGeneralPDE("y"),\
1608     self.getCoefficientOfGeneralPDE("y_contact"))
1609     self.__righthandside=self.copyConstraint(self.__righthandside)
1610     self.trace("New right hand side has been built.")
1611 jgs 108 self.__righthandside_isValid=True
1612 jgs 149 elif not self.__operator_is_Valid:
1613 jgs 121 self.getDomain().addPDEToSystem(self.__makeFreshOperator(),escript.Data(), \
1614 jgs 148 self.getCoefficientOfGeneralPDE("A"), \
1615     self.getCoefficientOfGeneralPDE("B"), \
1616     self.getCoefficientOfGeneralPDE("C"), \
1617     self.getCoefficientOfGeneralPDE("D"), \
1618 jgs 108 escript.Data(), \
1619     escript.Data(), \
1620 jgs 148 self.getCoefficientOfGeneralPDE("d"), \
1621 jgs 108 escript.Data(),\
1622 jgs 148 self.getCoefficientOfGeneralPDE("d_contact"), \
1623 jgs 108 escript.Data())
1624     self.__applyConstraint()
1625 jgs 148 self.trace("New operator has been built.")
1626 jgs 149 self.__operator_is_Valid=True
1627 jgs 108 return (self.__operator,self.__righthandside)
1628 jgs 102
1629    
1630 jgs 149 class Poisson(LinearPDE):
1631     """
1632     Class to define a Poisson equation problem, which is genear L{LinearPDE} of the form
1633 jgs 102
1634 jgs 149 M{-grad(grad(u)[j])[j] = f}
1635 jgs 102
1636 jgs 149 with natural boundary conditons
1637    
1638     M{n[j]*grad(u)[j] = 0 }
1639    
1640     and constraints:
1641    
1642     M{u=0} where M{q>0}
1643    
1644 jgs 108 """
1645 jgs 148
1646 gross 328 def __init__(self,domain,debug=False):
1647 jgs 149 """
1648     initializes a new Poisson equation
1649 jgs 104
1650 jgs 149 @param domain: domain of the PDE
1651     @type domain: L{Domain<escript.Domain>}
1652     @param debug: if True debug informations are printed.
1653 jgs 102
1654 jgs 149 """
1655 jgs 151 super(Poisson, self).__init__(domain,1,1,debug)
1656 jgs 149 self.COEFFICIENTS={"f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1657     "q": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
1658     self.setSymmetryOn()
1659 jgs 102
1660 jgs 149 def setValue(self,**coefficients):
1661     """
1662     sets new values to coefficients
1663 jgs 102
1664 jgs 149 @param coefficients: new values assigned to coefficients
1665     @keyword f: value for right hand side M{f}
1666     @type f: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1667     @keyword q: mask for location of constraints
1668     @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>}
1669     depending of reduced order is used for the representation of the equation.
1670     @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1671     """
1672 jgs 151 super(Poisson, self).setValue(**coefficients)
1673 jgs 102
1674 jgs 149 def getCoefficientOfGeneralPDE(self,name):
1675     """
1676     return the value of the coefficient name of the general PDE
1677     @param name: name of the coefficient requested.
1678     @type name: C{string}
1679     @return: the value of the coefficient name
1680     @rtype: L{Data<escript.Data>}
1681     @raise IllegalCoefficient: if name is not one of coefficients
1682     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}.
1683     @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
1684     """
1685     if name == "A" :
1686 gross 328 return escript.Data(util.kronecker(self.getDim()),escript.Function(self.getDomain()))
1687 jgs 149 elif name == "B" :
1688     return escript.Data()
1689     elif name == "C" :
1690     return escript.Data()
1691     elif name == "D" :
1692     return escript.Data()
1693     elif name == "X" :
1694     return escript.Data()
1695     elif name == "Y" :
1696     return self.getCoefficient("f")
1697     elif name == "d" :
1698     return escript.Data()
1699     elif name == "y" :
1700     return escript.Data()
1701     elif name == "d_contact" :
1702     return escript.Data()
1703     elif name == "y_contact" :
1704     return escript.Data()
1705     elif name == "r" :
1706     return escript.Data()
1707     elif name == "q" :
1708     return self.getCoefficient("q")
1709     else:
1710     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1711 jgs 102
1712 jgs 149 class Helmholtz(LinearPDE):
1713     """
1714     Class to define a Helmhotz equation problem, which is genear L{LinearPDE} of the form
1715    
1716     M{S{omega}*u - grad(k*grad(u)[j])[j] = f}
1717    
1718     with natural boundary conditons
1719    
1720     M{k*n[j]*grad(u)[j] = g- S{alpha}u }
1721    
1722 jgs 122 and constraints:
1723 jgs 102
1724 jgs 149 M{u=r} where M{q>0}
1725    
1726 jgs 110 """
1727 jgs 149
1728     def __init__(self,domain,debug=False):
1729     """
1730     initializes a new Poisson equation
1731    
1732     @param domain: domain of the PDE
1733     @type domain: L{Domain<escript.Domain>}
1734     @param debug: if True debug informations are printed.
1735    
1736     """
1737 jgs 151 super(Helmholtz, self).__init__(domain,1,1,debug)
1738 jgs 149 self.COEFFICIENTS={"omega": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
1739     "k": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
1740     "f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1741     "alpha": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
1742     "g": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1743     "r": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH),
1744     "q": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
1745     self.setSymmetryOn()
1746    
1747     def setValue(self,**coefficients):
1748     """
1749     sets new values to coefficients
1750    
1751     @param coefficients: new values assigned to coefficients
1752     @keyword omega: value for coefficient M{S{omega}}
1753     @type omega: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1754     @keyword k: value for coefficeint M{k}
1755     @type k: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1756     @keyword f: value for right hand side M{f}
1757     @type f: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1758     @keyword alpha: value for right hand side M{S{alpha}}
1759     @type alpha: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1760     @keyword g: value for right hand side M{g}
1761     @type g: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1762     @keyword r: prescribed values M{r} for the solution in constraints.
1763     @type r: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1764     depending of reduced order is used for the representation of the equation.
1765     @keyword q: mask for location of constraints
1766     @type q: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1767     depending of reduced order is used for the representation of the equation.
1768     @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1769     """
1770 jgs 151 super(Helmholtz, self).setValue(**coefficients)
1771 jgs 149
1772     def getCoefficientOfGeneralPDE(self,name):
1773     """
1774     return the value of the coefficient name of the general PDE
1775    
1776     @param name: name of the coefficient requested.
1777     @type name: C{string}
1778     @return: the value of the coefficient name
1779     @rtype: L{Data<escript.Data>}
1780     @raise IllegalCoefficient: if name is not one of coefficients
1781     "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}.
1782     @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
1783     """
1784     if name == "A" :
1785     return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))*self.getCoefficient("k")
1786     elif name == "B" :
1787     return escript.Data()
1788     elif name == "C" :
1789     return escript.Data()
1790     elif name == "D" :
1791     return self.getCoefficient("omega")
1792     elif name == "X" :
1793     return escript.Data()
1794     elif name == "Y" :
1795     return self.getCoefficient("f")
1796     elif name == "d" :
1797     return self.getCoefficient("alpha")
1798     elif name == "y" :
1799     return self.getCoefficient("g")
1800     elif name == "d_contact" :
1801     return escript.Data()
1802     elif name == "y_contact" :
1803     return escript.Data()
1804     elif name == "r" :
1805     return self.getCoefficient("r")
1806     elif name == "q" :
1807     return self.getCoefficient("q")
1808     else:
1809     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1810    
1811     class LameEquation(LinearPDE):
1812     """
1813     Class to define a Lame equation problem:
1814    
1815     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] }
1816    
1817     with natural boundary conditons:
1818    
1819     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] }
1820    
1821     and constraints:
1822    
1823     M{u[i]=r[i]} where M{q[i]>0}
1824    
1825     """
1826    
1827     def __init__(self,domain,debug=False):
1828 jgs 151 super(LameEquation, self).__init__(domain,\
1829     domain.getDim(),domain.getDim(),debug)
1830     self.COEFFICIENTS={ "lame_lambda" : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),
1831 jgs 149 "lame_mu" : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),
1832     "F" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1833     "sigma" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM),PDECoefficient.RIGHTHANDSIDE),
1834     "f" : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1835     "r" : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH),
1836     "q" : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
1837 jgs 151 self.setSymmetryOn()
1838 jgs 149
1839 gross 596 def setValues(self,**coefficients):
1840 jgs 149 """
1841     sets new values to coefficients
1842    
1843     @param coefficients: new values assigned to coefficients
1844     @keyword lame_mu: value for coefficient M{S{mu}}
1845     @type lame_mu: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1846     @keyword lame_lambda: value for coefficient M{S{lambda}}
1847     @type lame_lambda: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1848     @keyword F: value for internal force M{F}
1849     @type F: any type that can be casted to L{Vector<escript.Vector>} object on L{Function<escript.Function>}
1850     @keyword sigma: value for initial stress M{S{sigma}}
1851     @type sigma: any type that can be casted to L{Tensor<escript.Tensor>} object on L{Function<escript.Function>}
1852     @keyword f: value for extrenal force M{f}
1853     @type f: any type that can be casted to L{Vector<escript.Vector>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}
1854     @keyword r: prescribed values M{r} for the solution in constraints.
1855     @type r: any type that can be casted to L{Vector<escript.Vector>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1856     depending of reduced order is used for the representation of the equation.
1857     @keyword q: mask for location of constraints
1858     @type q: any type that can be casted to L{Vector<escript.Vector>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1859     depending of reduced order is used for the representation of the equation.
1860     @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1861     """
1862 gross 596 super(LameEquation, self).setValues(**coefficients)
1863 jgs 149
1864     def getCoefficientOfGeneralPDE(self,name):
1865     """
1866     return the value of the coefficient name of the general PDE
1867    
1868     @param name: name of the coefficient requested.
1869     @type name: C{string}
1870     @return: the value of the coefficient name
1871     @rtype: L{Data<escript.Data>}
1872     @raise IllegalCoefficient: if name is not one of coefficients
1873     "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}.
1874     @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
1875     """
1876     if name == "A" :
1877     out =self.createCoefficientOfGeneralPDE("A")
1878     for i in range(self.getDim()):
1879     for j in range(self.getDim()):
1880     out[i,i,j,j] += self.getCoefficient("lame_lambda")
1881     out[i,j,j,i] += self.getCoefficient("lame_mu")
1882     out[i,j,i,j] += self.getCoefficient("lame_mu")
1883     return out
1884     elif name == "B" :
1885     return escript.Data()
1886     elif name == "C" :
1887     return escript.Data()
1888     elif name == "D" :
1889     return escript.Data()
1890     elif name == "X" :
1891     return self.getCoefficient("sigma")
1892     elif name == "Y" :
1893     return self.getCoefficient("F")
1894     elif name == "d" :
1895     return escript.Data()
1896     elif name == "y" :
1897     return self.getCoefficient("f")
1898     elif name == "d_contact" :
1899     return escript.Data()
1900     elif name == "y_contact" :
1901     return escript.Data()
1902     elif name == "r" :
1903     return self.getCoefficient("r")
1904     elif name == "q" :
1905     return self.getCoefficient("q")
1906     else:
1907     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1908    
1909     class AdvectivePDE(LinearPDE):
1910     """
1911     In cases of PDEs dominated by the advection terms M{B} and M{C} against the adevctive terms M{A}
1912     up-winding has been used. The L{AdvectivePDE} class applies SUPG upwinding to the advective terms.
1913    
1914     In the following we set
1915    
1916     M{Z[j]=C[j]-B[j]}
1917    
1918 jgs 151 or
1919 jgs 149
1920     M{Z[i,k,l]=C[i,k,l]-B[i,l,k]}
1921    
1922 jgs 151 To measure the dominance of the advective terms over the diffusive term M{A} the
1923 jgs 149 X{Pelclet number} M{P} is used. It is defined as
1924 jgs 151
1925 jgs 149 M{P=h|Z|/(2|A|)}
1926    
1927 jgs 151 where M{|.|} denotes the L{length<util.length>} of the arument and M{h} is the local cell size
1928 jgs 149 from L{getSize<escript.Domain.getSize>}. Where M{|A|==0} M{P} is M{S{infinity}}.
1929    
1930     From the X{Pelclet number} the stabilization parameters M{S{Xi}} and M{S{Xi}} are calculated:
1931    
1932     M{S{Xi}=S{xi}(P) h/|Z|}
1933    
1934     where M{S{xi}} is a suitable function of the Peclet number.
1935    
1936     In the case of a single PDE the coefficient are up-dated in the following way:
1937     - M{A[i,j] S{<-} A[i,j] + S{Xi} * Z[j] * Z[l]}
1938     - M{B[j] S{<-} B[j] + S{Xi} * C[j] * D}
1939     - M{C[j] S{<-} C[j] + S{Xi} * B[j] * D}
1940     - M{X[j] S{<-} X[j] + S{Xi} * Z[j] * Y}
1941    
1942     Similar for the case of a systems of PDEs:
1943     - 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]}
1944     - M{B[i,j,k] S{<-} B[i,j,k] + S{delta}[p,m] * S{Xi} * D[p,k] * C[m,i,j]}
1945     - M{C[i,k,l] S{<-} C[i,k,l] + S{delta}[p,m] * S{Xi} * D[p,k] * B[m,l,i]}
1946     - M{X[i,j] S{<-} X[i,j] + S{delta}[p,m] * S{Xi} * Y[p] * Z[m,i,j]}
1947    
1948     where M{S{delta}} is L{kronecker}.
1949     Using upwinding in this form, introduces an additonal error which is proprtional to the cell size M{h}
1950     but with the intension to stabilize the solution.
1951    
1952     """
1953     def __init__(self,domain,numEquations=None,numSolutions=None,xi=None,debug=False):
1954     """
1955 jgs 151 creates a linear, steady, second order PDE on a L{Domain<escript.Domain>}
1956 jgs 149
1957     @param domain: domain of the PDE
1958     @type domain: L{Domain<escript.Domain>}
1959     @param numEquations: number of equations. If numEquations==None the number of equations
1960     is exracted from the PDE coefficients.
1961     @param numSolutions: number of solution components. If numSolutions==None the number of solution components
1962     is exracted from the PDE coefficients.
1963 jgs 151 @param xi: defines a function which returns for any given Preclet number as L{Scalar<escript.Scalar>} object the
1964 jgs 149 M{S{xi}}-value used to define the stabilization parameters. If equal to None, L{ELMAN_RAMAGE} is used.
1965     @type xi: callable object which returns a L{Scalar<escript.Scalar>} object.
1966     @param debug: if True debug informations are printed.
1967     """
1968 jgs 151 super(AdvectivePDE, self).__init__(domain,\
1969     numEquations,numSolutions,debug)
1970 jgs 148 if xi==None:
1971     self.__xi=AdvectivePDE.ELMAN_RAMAGE
1972     else:
1973     self.__xi=xi
1974 jgs 110 self.__Xi=escript.Data()
1975 jgs 102
1976 gross 596 def setValue(self,**coefficients):
1977 jgs 149 """
1978     sets new values to coefficients
1979 jgs 102
1980 jgs 149 @param coefficients: new values assigned to coefficients
1981     @keyword A: value for coefficient A.
1982     @type A: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1983     @keyword B: value for coefficient B
1984     @type B: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1985     @keyword C: value for coefficient C
1986     @type C: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1987     @keyword D: value for coefficient D
1988     @type D: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1989     @keyword X: value for coefficient X
1990     @type X: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1991     @keyword Y: value for coefficient Y
1992     @type Y: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1993     @keyword d: value for coefficient d
1994     @type d: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1995     @keyword y: value for coefficient y
1996     @type y: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1997     @keyword d_contact: value for coefficient d_contact
1998     @type d_contact: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnContactOne<escript.FunctionOnContactOne>}.
1999     or L{FunctionOnContactZero<escript.FunctionOnContactZero>}.
2000     @keyword y_contact: value for coefficient y_contact
2001     @type y_contact: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnContactOne<escript.FunctionOnContactOne>}.
2002     or L{FunctionOnContactZero<escript.FunctionOnContactZero>}.
2003     @keyword r: values prescribed to the solution at the locations of constraints
2004     @type r: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2005     depending of reduced order is used for the solution.
2006     @keyword q: mask for location of constraints
2007     @type q: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2008     depending of reduced order is used for the representation of the equation.
2009     @raise IllegalCoefficient: if an unknown coefficient keyword is used.
2010 jgs 148
2011 jgs 149 """
2012     if "A" in coefficients.keys() or "B" in coefficients.keys() or "C" in coefficients.keys(): self.__Xi=escript.Data()
2013 jgs 151 super(AdvectivePDE, self).setValue(**coefficients)
2014 jgs 149
2015     def ELMAN_RAMAGE(self,P):
2016 jgs 151 """
2017     Predefined function to set a values for M{S{xi}} from a Preclet number M{P}.
2018 jgs 149 This function uses the method suggested by H.C. Elman and A. Ramage, I{SIAM J. Numer. Anal.}, B{40} (2002)
2019     - M{S{xi}(P)=0} for M{P<1}
2020     - M{S{xi}(P)=(1-1/P)/2} otherwise
2021    
2022 jgs 151 @param P: Preclet number
2023 jgs 149 @type P: L{Scalar<escript.Scalar>}
2024     @return: up-wind weightimg factor
2025     @rtype: L{Scalar<escript.Scalar>}
2026 jgs 151 """
2027 gross 304 return util.wherePositive(P-1.)*0.5*(1.-1./(P+1.e-15))
2028 jgs 149
2029     def SIMPLIFIED_BROOK_HUGHES(self,P):
2030 jgs 151 """
2031     Predefined function to set a values for M{S{xi}} from a Preclet number M{P}.
2032 jgs 149 The original methods is
2033 jgs 151
2034 jgs 149 M{S{xi}(P)=coth(P)-1/P}
2035    
2036     As the evaluation of M{coth} is expensive we are using the approximation:
2037 jgs 151
2038 jgs 149 - M{S{xi}(P)=P/3} where M{P<3}
2039     - M{S{xi}(P)=1/2} otherwise
2040    
2041 jgs 151 @param P: Preclet number
2042 jgs 149 @type P: L{Scalar<escript.Scalar>}
2043     @return: up-wind weightimg factor
2044     @rtype: L{Scalar<escript.Scalar>}
2045 jgs 151 """
2046 gross 304 c=util.whereNegative(P-3.)
2047 jgs 148 return P/6.*c+1./2.*(1.-c)
2048    
2049 jgs 149 def HALF(self,P):
2050 jgs 151 """
2051     Predefined function to set value M{1/2} for M{S{xi}}
2052    
2053     @param P: Preclet number
2054 jgs 149 @type P: L{Scalar<escript.Scalar>}
2055     @return: up-wind weightimg factor
2056     @rtype: L{Scalar<escript.Scalar>}
2057 jgs 151 """
2058 jgs 149 return escript.Scalar(0.5,P.getFunctionSpace())
2059 jgs 148
2060 jgs 149 def __getXi(self):
2061 jgs 110 if self.__Xi.isEmpty():
2062     B=self.getCoefficient("B")
2063     C=self.getCoefficient("C")
2064     A=self.getCoefficient("A")
2065     h=self.getDomain().getSize()
2066     self.__Xi=escript.Scalar(0.,self.getFunctionSpaceForCoefficient("A"))
2067     if not C.isEmpty() or not B.isEmpty():
2068     if not C.isEmpty() and not B.isEmpty():
2069     if self.getNumEquations()>1:
2070     if self.getNumSolutions()>1:
2071 gross 596 flux2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))
2072 jgs 110 for i in range(self.getNumEquations()):
2073     for k in range(self.getNumSolutions()):
2074 gross 328 for l in range(self.getDim()): flux2+=(C[i,k,l]-B[i,l,k])**2
2075 gross 596 length_of_flux=util.sqrt(flux2)
2076 gross 328 # flux=C-util.reorderComponents(B,[0,2,1])
2077 jgs 110 else:
2078 gross 596 flux2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))
2079 jgs 110 for i in range(self.getNumEquations()):
2080 gross 328 for l in range(self.getDim()): flux2+=(C[i,l]-B[i,l])**2
2081 gross 596 length_of_flux=util.sqrt(flux2)
2082 gross 328 # flux=C-B
2083 jgs 110 else:
2084     if self.getNumSolutions()>1:
2085 gross 596 flux2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))
2086 jgs 110 for k in range(self.getNumSolutions()):
2087 gross 328 for l in range(self.getDim()): flux2+=(C[k,l]-B[l,k])**2
2088     # flux=C-util.reorderComponents(B,[1,0])
2089 gross 596 length_of_flux=util.sqrt(flux2)
2090 jgs 110 else:
2091 gross 596 length_of_flux=util.length(C-B)
2092 jgs 110 elif C.isEmpty():
2093 gross 328 length_of_flux=util.length(B)
2094 jgs 110 else:
2095 gross 328 length_of_flux=util.length(C)
2096     flux_max=util.Lsup(length_of_flux)
2097     if flux_max>0.:
2098 gross 596 if A.isEmpty():
2099     inv_A=1./self.SMALL_TOLERANCE
2100     peclet_number=escript.Scalar(inv_A,length_of_flux.getFunctionSpace())
2101     xi=self.__xi(self,peclet_number)
2102     else:
2103     # length_of_A=util.inner(flux,util.tensormutiply(A,flux))
2104     length_of_A=util.length(A)
2105     A_max=util.Lsup(length_of_A)
2106     if A_max>0:
2107     inv_A=1./(length_of_A+A_max*self.SMALL_TOLERANCE)
2108     else:
2109     inv_A=1./self.SMALL_TOLERANCE
2110     peclet_number=length_of_flux*h/2*inv_A
2111     xi=self.__xi(self,peclet_number)
2112     self.__Xi=h*xi/(length_of_flux+flux_max*self.SMALL_TOLERANCE)
2113     self.trace("preclet number = %e"%util.Lsup(peclet_number))
2114     else:
2115     self.__Xi=escript.Scalar(0.,length_of_flux.getFunctionSpace())
2116 jgs 110 return self.__Xi
2117 jgs 102
2118 jgs 148
2119     def getCoefficientOfGeneralPDE(self,name):
2120 jgs 108 """
2121 jgs 122 return the value of the coefficient name of the general PDE
2122    
2123 jgs 149 @param name: name of the coefficient requested.
2124     @type name: C{string}
2125     @return: the value of the coefficient name
2126     @rtype: L{Data<escript.Data>}
2127     @raise IllegalCoefficient: if name is not one of coefficients
2128     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}.
2129     @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
2130 jgs 108 """
2131 jgs 110 if not self.getNumEquations() == self.getNumSolutions():
2132     raise ValueError,"AdvectivePDE expects the number of solution componets and the number of equations to be equal."
2133    
2134 jgs 148 if name == "A" :
2135 jgs 108 A=self.getCoefficient("A")
2136     B=self.getCoefficient("B")
2137     C=self.getCoefficient("C")
2138 jgs 110 if B.isEmpty() and C.isEmpty():
2139     Aout=A
2140     else:
2141 jgs 148 if A.isEmpty():
2142 gross 596 Aout=self.createCoefficientOfGeneralPDE("A")
2143 jgs 110 else:
2144     Aout=A[:]
2145 jgs 149 Xi=self.__getXi()
2146 jgs 110 if self.getNumEquations()>1:
2147     for i in range(self.getNumEquations()):
2148     for j in range(self.getDim()):
2149 jgs 108 for k in range(self.getNumSolutions()):
2150 jgs 110 for l in range(self.getDim()):
2151     if not C.isEmpty() and not B.isEmpty():
2152 gross 328 # tmp=C-util.reorderComponents(B,[0,2,1])
2153     # Aout=Aout+Xi*util.generalTensorProduct(util.reorder(tmp,[1,2,0]),tmp,offset=1)
2154 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])
2155     elif C.isEmpty():
2156     for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*B[p,j,i]*B[p,l,k]
2157 gross 328 # Aout=Aout+Xi*util.generalTensorProduct(util.reorder(B,[2,1,0]),util.reorder(B,[0,2,1]),offset=1)
2158 jgs 110 else:
2159     for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*C[p,i,j]*C[p,k,l]
2160 gross 328 # Aout=Aout+Xi*util.generalTensorProduct(util.reorder(C,[1,2,0]),C,offset=1)
2161 jgs 110 else:
2162 gross 596 if not C.isEmpty() and not B.isEmpty():
2163     delta=(C-B)
2164     Aout+=util.outer(Xi*delta,delta)
2165     elif not B.isEmpty():
2166     Aout+=util.outer(Xi*B,B)
2167     elif not C.isEmpty():
2168     Aout+=util.outer(Xi*C,C)
2169 jgs 110 return Aout
2170 jgs 148 elif name == "B" :
2171 gross 596 # return self.getCoefficient("B")
2172 jgs 110 B=self.getCoefficient("B")
2173     C=self.getCoefficient("C")
2174     D=self.getCoefficient("D")
2175     if C.isEmpty() or D.isEmpty():
2176     Bout=B
2177     else:
2178 jgs 149 Xi=self.__getXi()
2179 jgs 148 if B.isEmpty():
2180 gross 596 Bout=self.createCoefficientOfGeneralPDE("B")
2181 jgs 110 else:
2182     Bout=B[:]
2183     if self.getNumEquations()>1:
2184     for k in range(self.getNumSolutions()):
2185 jgs 148 for p in range(self.getNumEquations()):
2186 jgs 110 tmp=Xi*D[p,k]
2187     for i in range(self.getNumEquations()):
2188     for j in range(self.getDim()):
2189     Bout[i,j,k]+=tmp*C[p,i,j]
2190 gross 328 # Bout=Bout+Xi*util.generalTensorProduct(util.reorder(C,[1,2,0]),D,offset=1)
2191 jgs 110 else:
2192 gross 596 Bout+=(Xi*D)*C
2193 jgs 110 return Bout
2194 jgs 148 elif name == "C" :
2195 gross 596 # return self.getCoefficient("C")
2196 jgs 110 B=self.getCoefficient("B")
2197     C=self.getCoefficient("C")
2198     D=self.getCoefficient("D")
2199     if B.isEmpty() or D.isEmpty():
2200     Cout=C
2201     else:
2202 jgs 149 Xi=self.__getXi()
2203 jgs 148 if C.isEmpty():
2204 gross 596 Cout=self.createCoefficientOfGeneralPDE("C")
2205 jgs 110 else:
2206     Cout=C[:]
2207     if self.getNumEquations()>1:
2208     for k in range(self.getNumSolutions()):
2209     for p in range(self.getNumEquations()):
2210     tmp=Xi*D[p,k]
2211     for i in range(self.getNumEquations()):
2212     for l in range(self.getDim()):
2213     Cout[i,k,l]+=tmp*B[p,l,i]
2214 gross 328 # Cout=Cout+Xi*B[p,l,i]*D[p,k]
2215 jgs 110 else:
2216 gross 596 Cout+=(Xi*D)*B
2217 jgs 110 return Cout
2218 jgs 148 elif name == "D" :
2219 jgs 108 return self.getCoefficient("D")
2220 jgs 148 elif name == "X" :
2221 gross 596 # return self.getCoefficient("X")
2222 jgs 110 X=self.getCoefficient("X")
2223     Y=self.getCoefficient("Y")
2224     B=self.getCoefficient("B")
2225     C=self.getCoefficient("C")
2226     if Y.isEmpty() or (B.isEmpty() and C.isEmpty()):
2227     Xout=X
2228     else:
2229     if X.isEmpty():
2230 gross 596 Xout=self.createCoefficientOfGeneralPDE("X")
2231 jgs 110 else:
2232     Xout=X[:]
2233 jgs 149 Xi=self.__getXi()
2234 jgs 110 if self.getNumEquations()>1:
2235 jgs 148 for p in range(self.getNumEquations()):
2236 jgs 110 tmp=Xi*Y[p]
2237     for i in range(self.getNumEquations()):
2238     for j in range(self.getDim()):
2239     if not C.isEmpty() and not B.isEmpty():
2240     Xout[i,j]+=tmp*(C[p,i,j]-B[p,j,i])
2241 gross 328 # Xout=X_out+Xi*util.inner(Y,C-util.reorderComponents(B,[0,2,1]),offset=1)
2242 jgs 110 elif C.isEmpty():
2243     Xout[i,j]-=tmp*B[p,j,i]
2244 gross 328 # Xout=X_out-Xi*util.inner(Y,util.reorderComponents(B,[0,2,1]),offset=1)
2245 jgs 110 else:
2246     Xout[i,j]+=tmp*C[p,i,j]
2247 gross 328 # Xout=X_out+Xi*util.inner(Y,C,offset=1)
2248 jgs 110 else:
2249 gross 596 if not C.isEmpty() and not B.isEmpty():
2250     Xout+=(Xi*Y)*(C-B)
2251     elif C.isEmpty():
2252     Xout-=(Xi*Y)*B
2253     else:
2254     Xout+=(Xi*Y)*C
2255 jgs 110 return Xout
2256 jgs 148 elif name == "Y" :
2257 jgs 108 return self.getCoefficient("Y")
2258 jgs 148 elif name == "d" :
2259 jgs 108 return self.getCoefficient("d")
2260 jgs 148 elif name == "y" :
2261 jgs 108 return self.getCoefficient("y")
2262 jgs 148 elif name == "d_contact" :
2263 jgs 108 return self.getCoefficient("d_contact")
2264     elif name == "y_contact" :
2265     return self.getCoefficient("y_contact")
2266 jgs 148 elif name == "r" :
2267 jgs 108 return self.getCoefficient("r")
2268 jgs 148 elif name == "q" :
2269 jgs 108 return self.getCoefficient("q")
2270     else:
2271 jgs 149 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
2272 jgs 108
2273 jgs 121 # $Log$
2274 jgs 151 # Revision 1.14 2005/09/22 01:54:57 jgs
2275     # Merge of development branch dev-02 back to main trunk on 2005-09-22
2276     #
2277 jgs 150 # Revision 1.13 2005/09/15 03:44:19 jgs
2278     # Merge of development branch dev-02 back to main trunk on 2005-09-15
2279     #
2280 jgs 149 # Revision 1.12 2005/09/01 03:31:28 jgs
2281     # Merge of development branch dev-02 back to main trunk on 2005-09-01
2282     #
2283 jgs 148 # Revision 1.11 2005/08/23 01:24:28 jgs
2284     # Merge of development branch dev-02 back to main trunk on 2005-08-23
2285     #
2286 jgs 147 # Revision 1.10 2005/08/12 01:45:36 jgs
2287     # erge of development branch dev-02 back to main trunk on 2005-08-12
2288 jgs 142 #
2289 jgs 151 # Revision 1.9.2.17 2005/09/21 07:03:33 matt
2290     # PDECoefficient and LinearPDE are now new style classes (introduced in Python
2291     # 2.2). Classes Poisson, Helmholtz, LameEquation and AdvectivePDE have been
2292     # modified to instead use portable/cooperative "super" calls to extend base
2293     # class methods.
2294     #
2295     # Revision 1.9.2.16 2005/09/16 01:54:37 matt
2296     # Removed redundant if-loop.
2297     #
2298 jgs 150 # Revision 1.9.2.15 2005/09/14 08:09:18 matt
2299     # Added "REDUCED" solution PDECoefficient descriptors for LinearPDEs.
2300     #
2301     # Revision 1.9.2.14 2005/09/07 06:26:16 gross
2302     # the solver from finley are put into the standalone package paso now
2303     #
2304 jgs 149 # Revision 1.9.2.13 2005/08/31 08:45:03 gross
2305     # in the case of lumping no new system is allocated if the constraint is changed.
2306     #
2307     # Revision 1.9.2.12 2005/08/31 07:10:23 gross
2308     # test for Lumping added
2309     #
2310     # Revision 1.9.2.11 2005/08/30 01:53:45 gross
2311     # bug in format fixed.
2312     #
2313     # Revision 1.9.2.10 2005/08/26 07:14:17 gross
2314     # a few more bugs in linearPDE fixed. remaining problem are finley problems
2315     #
2316     # Revision 1.9.2.9 2005/08/26 06:30:45 gross
2317     # fix for reported bug 0000004. test_linearPDE passes a few more tests
2318     #
2319     # Revision 1.9.2.8 2005/08/26 04:30:13 gross
2320     # gneric unit testing for linearPDE
2321     #
2322     # Revision 1.9.2.7 2005/08/25 07:06:50 gross
2323     # linearPDE documentation is parsed now by epydoc. there is still a problem with links into escriptcpp.so
2324     #
2325     # Revision 1.9.2.6 2005/08/24 05:01:24 gross
2326     # problem with resetting the matrix in case of resetting its values to 0 fixed.
2327     #
2328     # Revision 1.9.2.5 2005/08/24 02:03:28 gross
2329     # epydoc mark up partially fixed
2330     #
2331 jgs 148 # Revision 1.9.2.4 2005/08/22 07:11:09 gross
2332     # some problems with LinearPDEs fixed.
2333     #
2334     # Revision 1.9.2.3 2005/08/18 04:48:48 gross
2335     # the methods SetLumping*() are removed. Lumping is set trough setSolverMethod(LinearPDE.LUMPING)
2336     #
2337     # Revision 1.9.2.2 2005/08/18 04:39:32 gross
2338     # the constants have been removed from util.py as they not needed anymore. PDE related constants are accessed through LinearPDE attributes now
2339     #
2340 jgs 147 # Revision 1.9.2.1 2005/07/29 07:10:27 gross
2341     # new functions in util and a new pde type in linearPDEs
2342 jgs 122 #
2343 jgs 147 # Revision 1.1.2.25 2005/07/28 04:21:09 gross
2344     # Lame equation: (linear elastic, isotropic) added
2345 jgs 121 #
2346 jgs 142 # Revision 1.1.2.24 2005/07/22 06:37:11 gross
2347     # some extensions to modellib and linearPDEs
2348     #
2349 jgs 122 # Revision 1.1.2.23 2005/05/13 00:55:20 cochrane
2350     # Fixed up some docstrings. Moved module-level functions to top of file so
2351     # that epydoc and doxygen can pick them up properly.
2352     #
2353     # Revision 1.1.2.22 2005/05/12 11:41:30 gross
2354     # some basic Models have been added
2355     #
2356     # Revision 1.1.2.21 2005/05/12 07:16:12 cochrane
2357     # Moved ELMAN_RAMAGE, SIMPLIFIED_BROOK_HUGHES, and HALF functions to bottom of
2358     # file so that the AdvectivePDE class is picked up by doxygen. Some
2359     # reformatting of docstrings. Addition of code to make equations come out
2360     # as proper LaTeX.
2361     #
2362 jgs 121 # Revision 1.1.2.20 2005/04/15 07:09:08 gross
2363     # some problems with functionspace and linearPDEs fixed.
2364     #
2365     # Revision 1.1.2.19 2005/03/04 05:27:07 gross
2366     # bug in SystemPattern fixed.
2367     #
2368     # Revision 1.1.2.18 2005/02/08 06:16:45 gross
2369     # Bugs in AdvectivePDE fixed, AdvectiveTest is stable but more testing is needed
2370     #
2371     # Revision 1.1.2.17 2005/02/08 05:56:19 gross
2372     # Reference Number handling added
2373     #
2374     # Revision 1.1.2.16 2005/02/07 04:41:28 gross
2375     # some function exposed to python to make mesh merging running
2376     #
2377     # Revision 1.1.2.15 2005/02/03 00:14:44 gross
2378     # timeseries add and ESySParameter.py renames esysXML.py for consistence
2379     #
2380     # Revision 1.1.2.14 2005/02/01 06:44:10 gross
2381     # new implementation of AdvectivePDE which now also updates right hand side. systems of PDE