/[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 614 - (hide annotations)
Wed Mar 22 01:37:07 2006 UTC (13 years, 8 months ago) by elspeth
File MIME type: text/x-python
File size: 105731 byte(s)
Corrected spelling of 'license' in url so that the link actually points to the license.

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