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