/[escript]/trunk/escript/py_src/linearPDEs.py
ViewVC logotype

Contents of /trunk/escript/py_src/linearPDEs.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1072 - (show annotations)
Thu Mar 29 06:44:30 2007 UTC (12 years, 6 months ago) by gross
File MIME type: text/x-python
File size: 102599 byte(s)
PDE assemblage for reduced integration order + tests added.


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26