/[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 150 - (hide annotations)
Thu Sep 15 03:44:45 2005 UTC (14 years, 1 month ago) by jgs
Original Path: trunk/esys2/escript/py_src/linearPDEs.py
File MIME type: text/x-python
File size: 102221 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-09-15

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

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26