/[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 969 - (hide annotations)
Tue Feb 13 23:02:23 2007 UTC (12 years, 7 months ago) by ksteube
File MIME type: text/x-python
File size: 106141 byte(s)
Parallelization using MPI for solution of implicit problems.

Parallelization for explicit problems has already been accomplished in
the main SVN branch.

This is incomplete and is not ready for use.


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