/[escript]/trunk-mpi-branch/escript/py_src/linearPDEs.py
ViewVC logotype

Contents of /trunk-mpi-branch/escript/py_src/linearPDEs.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1302 - (show annotations)
Mon Sep 17 06:57:51 2007 UTC (11 years, 11 months ago) by ksteube
File MIME type: text/x-python
File size: 102287 byte(s)
Some of the epydoc doesn't show up or doesn't get formatted, in this
case because of a missing } or because of improper indentation.

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 TRILINOS: the TRILINOS parallel solver class library from Sandia Natl Labs
447 @cvar ITERATIVE: The default iterative solver
448 @cvar AMG: algebraic multi grid
449 @cvar RILU: recursive ILU
450
451 """
452 DEFAULT= 0
453 DIRECT= 1
454 CHOLEVSKY= 2
455 PCG= 3
456 CR= 4
457 CGS= 5
458 BICGSTAB= 6
459 SSOR= 7
460 ILU0= 8
461 ILUT= 9
462 JACOBI= 10
463 GMRES= 11
464 PRES20= 12
465 LUMPING= 13
466 NO_REORDERING= 17
467 MINIMUM_FILL_IN= 18
468 NESTED_DISSECTION= 19
469 SCSL= 14
470 MKL= 15
471 UMFPACK= 16
472 ITERATIVE= 20
473 PASO= 21
474 AMG= 22
475 RILU = 23
476 TRILINOS = 24
477
478 SMALL_TOLERANCE=1.e-13
479 __PACKAGE_KEY="package"
480 __METHOD_KEY="method"
481 __SYMMETRY_KEY="symmetric"
482 __TOLERANCE_KEY="tolerance"
483 __PRECONDITIONER_KEY="preconditioner"
484
485
486 def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):
487 """
488 initializes a new linear PDE
489
490 @param domain: domain of the PDE
491 @type domain: L{Domain<escript.Domain>}
492 @param numEquations: number of equations. If numEquations==None the number of equations
493 is exracted from the PDE coefficients.
494 @param numSolutions: number of solution components. If numSolutions==None the number of solution components
495 is exracted from the PDE coefficients.
496 @param debug: if True debug informations are printed.
497
498 """
499 super(LinearPDE, self).__init__()
500 #
501 # the coefficients of the general PDE:
502 #
503 self.__COEFFICIENTS_OF_GENEARL_PDE={
504 "A" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM,PDECoefficient.BY_SOLUTION,PDECoefficient.BY_DIM),PDECoefficient.OPERATOR),
505 "B" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
506 "C" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION,PDECoefficient.BY_DIM),PDECoefficient.OPERATOR),
507 "D" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
508 "X" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM),PDECoefficient.RIGHTHANDSIDE),
509 "Y" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
510 "d" : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
511 "y" : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
512 "d_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
513 "y_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
514 "A_reduced" : PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM,PDECoefficient.BY_SOLUTION,PDECoefficient.BY_DIM),PDECoefficient.OPERATOR),
515 "B_reduced" : PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
516 "C_reduced" : PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION,PDECoefficient.BY_DIM),PDECoefficient.OPERATOR),
517 "D_reduced" : PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
518 "X_reduced" : PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM),PDECoefficient.RIGHTHANDSIDE),
519 "Y_reduced" : PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
520 "d_reduced" : PDECoefficient(PDECoefficient.BOUNDARY_REDUCED,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
521 "y_reduced" : PDECoefficient(PDECoefficient.BOUNDARY_REDUCED,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
522 "d_contact_reduced" : PDECoefficient(PDECoefficient.CONTACT_REDUCED,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_SOLUTION),PDECoefficient.OPERATOR),
523 "y_contact_reduced" : PDECoefficient(PDECoefficient.CONTACT_REDUCED,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
524 "r" : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_SOLUTION,),PDECoefficient.RIGHTHANDSIDE),
525 "q" : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_SOLUTION,),PDECoefficient.BOTH)}
526
527 # COEFFICIENTS can be overwritten by subclasses:
528 self.COEFFICIENTS=self.__COEFFICIENTS_OF_GENEARL_PDE
529 self.__altered_coefficients=False
530 # initialize attributes
531 self.__debug=debug
532 self.__domain=domain
533 self.__numEquations=numEquations
534 self.__numSolutions=numSolutions
535 self.__resetSystem()
536
537 # set some default values:
538 self.__reduce_equation_order=False
539 self.__reduce_solution_order=False
540 self.__tolerance=1.e-8
541 self.__solver_method=self.DEFAULT
542 self.__solver_package=self.DEFAULT
543 self.__preconditioner=self.DEFAULT
544 self.__matrix_type=self.__domain.getSystemMatrixTypeId(self.DEFAULT,self.DEFAULT,False)
545 self.__sym=False
546
547 self.resetCoefficients()
548 self.trace("PDE Coeffients are %s"%str(self.COEFFICIENTS.keys()))
549 # =============================================================================
550 # general stuff:
551 # =============================================================================
552 def __str__(self):
553 """
554 returns string representation of the PDE
555
556 @return: a simple representation of the PDE
557 @rtype: C{str}
558 """
559 return "<LinearPDE %d>"%id(self)
560 # =============================================================================
561 # debug :
562 # =============================================================================
563 def setDebugOn(self):
564 """
565 switches on debugging
566 """
567 self.__debug=not None
568
569 def setDebugOff(self):
570 """
571 switches off debugging
572 """
573 self.__debug=None
574
575 def trace(self,text):
576 """
577 print the text message if debugging is swiched on.
578 @param text: message
579 @type text: C{string}
580 """
581 if self.__debug: print "%s: %s"%(str(self),text)
582
583 # =============================================================================
584 # some service functions:
585 # =============================================================================
586 def getDomain(self):
587 """
588 returns the domain of the PDE
589
590 @return: the domain of the PDE
591 @rtype: L{Domain<escript.Domain>}
592 """
593 return self.__domain
594
595 def getDim(self):
596 """
597 returns the spatial dimension of the PDE
598
599 @return: the spatial dimension of the PDE domain
600 @rtype: C{int}
601 """
602 return self.getDomain().getDim()
603
604 def getNumEquations(self):
605 """
606 returns the number of equations
607
608 @return: the number of equations
609 @rtype: C{int}
610 @raise UndefinedPDEError: if the number of equations is not be specified yet.
611 """
612 if self.__numEquations==None:
613 raise UndefinedPDEError,"Number of equations is undefined. Please specify argument numEquations."
614 else:
615 return self.__numEquations
616
617 def getNumSolutions(self):
618 """
619 returns the number of unknowns
620
621 @return: the number of unknowns
622 @rtype: C{int}
623 @raise UndefinedPDEError: if the number of unknowns is not be specified yet.
624 """
625 if self.__numSolutions==None:
626 raise UndefinedPDEError,"Number of solution is undefined. Please specify argument numSolutions."
627 else:
628 return self.__numSolutions
629
630 def reduceEquationOrder(self):
631 """
632 return status for order reduction for equation
633
634 @return: return True is reduced interpolation order is used for the represenation of the equation
635 @rtype: L{bool}
636 """
637 return self.__reduce_equation_order
638
639 def reduceSolutionOrder(self):
640 """
641 return status for order reduction for the solution
642
643 @return: return True is reduced interpolation order is used for the represenation of the solution
644 @rtype: L{bool}
645 """
646 return self.__reduce_solution_order
647
648 def getFunctionSpaceForEquation(self):
649 """
650 returns the L{FunctionSpace<escript.FunctionSpace>} used to discretize the equation
651
652 @return: representation space of equation
653 @rtype: L{FunctionSpace<escript.FunctionSpace>}
654 """
655 if self.reduceEquationOrder():
656 return escript.ReducedSolution(self.getDomain())
657 else:
658 return escript.Solution(self.getDomain())
659
660 def getFunctionSpaceForSolution(self):
661 """
662 returns the L{FunctionSpace<escript.FunctionSpace>} used to represent the solution
663
664 @return: representation space of solution
665 @rtype: L{FunctionSpace<escript.FunctionSpace>}
666 """
667 if self.reduceSolutionOrder():
668 return escript.ReducedSolution(self.getDomain())
669 else:
670 return escript.Solution(self.getDomain())
671
672
673 def getOperator(self):
674 """
675 provides access to the operator of the PDE
676
677 @return: the operator of the PDE
678 @rtype: L{Operator<escript.Operator>}
679 """
680 m=self.getSystem()[0]
681 if self.isUsingLumping():
682 return self.copyConstraint(1./m)
683 else:
684 return m
685
686 def getRightHandSide(self):
687 """
688 provides access to the right hand side of the PDE
689 @return: the right hand side of the PDE
690 @rtype: L{Data<escript.Data>}
691 """
692 r=self.getSystem()[1]
693 if self.isUsingLumping():
694 return self.copyConstraint(r)
695 else:
696 return r
697
698 def applyOperator(self,u=None):
699 """
700 applies the operator of the PDE to a given u or the solution of PDE if u is not present.
701
702 @param u: argument of the operator. It must be representable in C{elf.getFunctionSpaceForSolution()}. If u is not present or equals L{None}
703 the current solution is used.
704 @type u: L{Data<escript.Data>} or None
705 @return: image of u
706 @rtype: L{Data<escript.Data>}
707 """
708 if u==None:
709 return self.getOperator()*self.getSolution()
710 else:
711 return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())
712
713 def getResidual(self,u=None):
714 """
715 return the residual of u or the current solution if u is not present.
716
717 @param u: argument in the residual calculation. It must be representable in C{elf.getFunctionSpaceForSolution()}. If u is not present or equals L{None}
718 the current solution is used.
719 @type u: L{Data<escript.Data>} or None
720 @return: residual of u
721 @rtype: L{Data<escript.Data>}
722 """
723 return self.applyOperator(u)-self.getRightHandSide()
724
725 def checkSymmetry(self,verbose=True):
726 """
727 test the PDE for symmetry.
728
729 @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.
730 @type verbose: C{bool}
731 @return: True if the PDE is symmetric.
732 @rtype: L{Data<escript.Data>}
733 @note: This is a very expensive operation. It should be used for degugging only! The symmetry flag is not altered.
734 """
735 verbose=verbose or self.__debug
736 out=True
737 if self.getNumSolutions()!=self.getNumEquations():
738 if verbose: print "non-symmetric PDE because of different number of equations and solutions"
739 out=False
740 else:
741 A=self.getCoefficientOfGeneralPDE("A")
742 if not A.isEmpty():
743 tol=util.Lsup(A)*self.SMALL_TOLERANCE
744 if self.getNumSolutions()>1:
745 for i in range(self.getNumEquations()):
746 for j in range(self.getDim()):
747 for k in range(self.getNumSolutions()):
748 for l in range(self.getDim()):
749 if util.Lsup(A[i,j,k,l]-A[k,l,i,j])>tol:
750 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)
751 out=False
752 else:
753 for j in range(self.getDim()):
754 for l in range(self.getDim()):
755 if util.Lsup(A[j,l]-A[l,j])>tol:
756 if verbose: print "non-symmetric PDE because A[%d,%d]!=A[%d,%d]"%(j,l,l,j)
757 out=False
758 B=self.getCoefficientOfGeneralPDE("B")
759 C=self.getCoefficientOfGeneralPDE("C")
760 if B.isEmpty() and not C.isEmpty():
761 if verbose: print "non-symmetric PDE because B is not present but C is"
762 out=False
763 elif not B.isEmpty() and C.isEmpty():
764 if verbose: print "non-symmetric PDE because C is not present but B is"
765 out=False
766 elif not B.isEmpty() and not C.isEmpty():
767 tol=(util.Lsup(B)+util.Lsup(C))*self.SMALL_TOLERANCE/2.
768 if self.getNumSolutions()>1:
769 for i in range(self.getNumEquations()):
770 for j in range(self.getDim()):
771 for k in range(self.getNumSolutions()):
772 if util.Lsup(B[i,j,k]-C[k,i,j])>tol:
773 if verbose: print "non-symmetric PDE because B[%d,%d,%d]!=C[%d,%d,%d]"%(i,j,k,k,i,j)
774 out=False
775 else:
776 for j in range(self.getDim()):
777 if util.Lsup(B[j]-C[j])>tol:
778 if verbose: print "non-symmetric PDE because B[%d]!=C[%d]"%(j,j)
779 out=False
780 if self.getNumSolutions()>1:
781 D=self.getCoefficientOfGeneralPDE("D")
782 if not D.isEmpty():
783 tol=util.Lsup(D)*self.SMALL_TOLERANCE
784 for i in range(self.getNumEquations()):
785 for k in range(self.getNumSolutions()):
786 if util.Lsup(D[i,k]-D[k,i])>tol:
787 if verbose: print "non-symmetric PDE because D[%d,%d]!=D[%d,%d]"%(i,k,k,i)
788 out=False
789 d=self.getCoefficientOfGeneralPDE("d")
790 if not d.isEmpty():
791 tol=util.Lsup(d)*self.SMALL_TOLERANCE
792 for i in range(self.getNumEquations()):
793 for k in range(self.getNumSolutions()):
794 if util.Lsup(d[i,k]-d[k,i])>tol:
795 if verbose: print "non-symmetric PDE because d[%d,%d]!=d[%d,%d]"%(i,k,k,i)
796 out=False
797 d_contact=self.getCoefficientOfGeneralPDE("d_contact")
798 if not d_contact.isEmpty():
799 tol=util.Lsup(d_contact)*self.SMALL_TOLERANCE
800 for i in range(self.getNumEquations()):
801 for k in range(self.getNumSolutions()):
802 if util.Lsup(d_contact[i,k]-d_contact[k,i])>tol:
803 if verbose: print "non-symmetric PDE because d_contact[%d,%d]!=d_contact[%d,%d]"%(i,k,k,i)
804 out=False
805 # and now the reduced coefficients
806 A_reduced=self.getCoefficientOfGeneralPDE("A_reduced")
807 if not A_reduced.isEmpty():
808 tol=util.Lsup(A_reduced)*self.SMALL_TOLERANCE
809 if self.getNumSolutions()>1:
810 for i in range(self.getNumEquations()):
811 for j in range(self.getDim()):
812 for k in range(self.getNumSolutions()):
813 for l in range(self.getDim()):
814 if util.Lsup(A_reduced[i,j,k,l]-A_reduced[k,l,i,j])>tol:
815 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)
816 out=False
817 else:
818 for j in range(self.getDim()):
819 for l in range(self.getDim()):
820 if util.Lsup(A_reduced[j,l]-A_reduced[l,j])>tol:
821 if verbose: print "non-symmetric PDE because A_reduced[%d,%d]!=A_reduced[%d,%d]"%(j,l,l,j)
822 out=False
823 B_reduced=self.getCoefficientOfGeneralPDE("B_reduced")
824 C_reduced=self.getCoefficientOfGeneralPDE("C_reduced")
825 if B_reduced.isEmpty() and not C_reduced.isEmpty():
826 if verbose: print "non-symmetric PDE because B_reduced is not present but C_reduced is"
827 out=False
828 elif not B_reduced.isEmpty() and C_reduced.isEmpty():
829 if verbose: print "non-symmetric PDE because C_reduced is not present but B_reduced is"
830 out=False
831 elif not B_reduced.isEmpty() and not C_reduced.isEmpty():
832 tol=(util.Lsup(B_reduced)+util.Lsup(C_reduced))*self.SMALL_TOLERANCE/2.
833 if self.getNumSolutions()>1:
834 for i in range(self.getNumEquations()):
835 for j in range(self.getDim()):
836 for k in range(self.getNumSolutions()):
837 if util.Lsup(B_reduced[i,j,k]-C_reduced[k,i,j])>tol:
838 if verbose: print "non-symmetric PDE because B_reduced[%d,%d,%d]!=C_reduced[%d,%d,%d]"%(i,j,k,k,i,j)
839 out=False
840 else:
841 for j in range(self.getDim()):
842 if util.Lsup(B_reduced[j]-C_reduced[j])>tol:
843 if verbose: print "non-symmetric PDE because B_reduced[%d]!=C_reduced[%d]"%(j,j)
844 out=False
845 if self.getNumSolutions()>1:
846 D_reduced=self.getCoefficientOfGeneralPDE("D_reduced")
847 if not D_reduced.isEmpty():
848 tol=util.Lsup(D_reduced)*self.SMALL_TOLERANCE
849 for i in range(self.getNumEquations()):
850 for k in range(self.getNumSolutions()):
851 if util.Lsup(D_reduced[i,k]-D_reduced[k,i])>tol:
852 if verbose: print "non-symmetric PDE because D_reduced[%d,%d]!=D_reduced[%d,%d]"%(i,k,k,i)
853 out=False
854 d_reduced=self.getCoefficientOfGeneralPDE("d_reduced")
855 if not d_reduced.isEmpty():
856 tol=util.Lsup(d_reduced)*self.SMALL_TOLERANCE
857 for i in range(self.getNumEquations()):
858 for k in range(self.getNumSolutions()):
859 if util.Lsup(d_reduced[i,k]-d_reduced[k,i])>tol:
860 if verbose: print "non-symmetric PDE because d_reduced[%d,%d]!=d_reduced[%d,%d]"%(i,k,k,i)
861 out=False
862 d_contact_reduced=self.getCoefficientOfGeneralPDE("d_contact_reduced")
863 if not d_contact_reduced.isEmpty():
864 tol=util.Lsup(d_contact_reduced)*self.SMALL_TOLERANCE
865 for i in range(self.getNumEquations()):
866 for k in range(self.getNumSolutions()):
867 if util.Lsup(d_contact_reduced[i,k]-d_contact_reduced[k,i])>tol:
868 if verbose: print "non-symmetric PDE because d_contact_reduced[%d,%d]!=d_contact_reduced[%d,%d]"%(i,k,k,i)
869 out=False
870 return out
871
872 def getSolution(self,**options):
873 """
874 returns the solution of the PDE. If the solution is not valid the PDE is solved.
875
876 @return: the solution
877 @rtype: L{Data<escript.Data>}
878 @param options: solver options
879 @keyword verbose: True to get some information during PDE solution
880 @type verbose: C{bool}
881 @keyword reordering: reordering scheme to be used during elimination. Allowed values are
882 L{NO_REORDERING}, L{MINIMUM_FILL_IN}, L{NESTED_DISSECTION}
883 @keyword iter_max: maximum number of iteration steps allowed.
884 @keyword drop_tolerance: threshold for drupping in L{ILUT}
885 @keyword drop_storage: maximum of allowed memory in L{ILUT}
886 @keyword truncation: maximum number of residuals in L{GMRES}
887 @keyword restart: restart cycle length in L{GMRES}
888 """
889 if not self.__solution_isValid:
890 mat,f=self.getSystem()
891 if self.isUsingLumping():
892 self.__solution=self.copyConstraint(f*mat)
893 else:
894 options[self.__TOLERANCE_KEY]=self.getTolerance()
895 options[self.__METHOD_KEY]=self.getSolverMethod()[0]
896 options[self.__PRECONDITIONER_KEY]=self.getSolverMethod()[1]
897 options[self.__PACKAGE_KEY]=self.getSolverPackage()
898 options[self.__SYMMETRY_KEY]=self.isSymmetric()
899 self.trace("PDE is resolved.")
900 self.trace("solver options: %s"%str(options))
901 self.__solution=mat.solve(f,options)
902 self.__solution_isValid=True
903 return self.__solution
904
905 def getFlux(self,u=None):
906 """
907 returns the flux M{J} for a given M{u}
908
909 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]}
910
911 or
912
913 M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[l]+(B[j]+B_reduced[j])u-X[j]-X_reduced[j]}
914
915 @param u: argument in the flux. If u is not present or equals L{None} the current solution is used.
916 @type u: L{Data<escript.Data>} or None
917 @return: flux
918 @rtype: L{Data<escript.Data>}
919 """
920 if u==None: u=self.getSolution()
921 return util.tensormult(self.getCoefficientOfGeneralPDE("A"),util.grad(u,Funtion(self.getDomain))) \
922 +util.matrixmult(self.getCoefficientOfGeneralPDE("B"),u) \
923 -util.self.getCoefficientOfGeneralPDE("X") \
924 +util.tensormult(self.getCoefficientOfGeneralPDE("A_reduced"),util.grad(u,ReducedFuntion(self.getDomain))) \
925 +util.matrixmult(self.getCoefficientOfGeneralPDE("B_reduced"),u) \
926 -util.self.getCoefficientOfGeneralPDE("X_reduced")
927 # =============================================================================
928 # solver settings:
929 # =============================================================================
930 def setSolverMethod(self,solver=None,preconditioner=None):
931 """
932 sets a new solver
933
934 @param solver: sets a new solver method.
935 @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}
936 @param preconditioner: sets a new solver method.
937 @type preconditioner: one of L{DEFAULT}, L{JACOBI} L{ILU0}, L{ILUT},L{SSOR}, L{RILU}
938 """
939 if solver==None: solver=self.__solver_method
940 if preconditioner==None: preconditioner=self.__preconditioner
941 if solver==None: solver=self.DEFAULT
942 if preconditioner==None: preconditioner=self.DEFAULT
943 if not (solver,preconditioner)==self.getSolverMethod():
944 self.__solver_method=solver
945 self.__preconditioner=preconditioner
946 self.__checkMatrixType()
947 self.trace("New solver is %s"%self.getSolverMethodName())
948
949 def getSolverMethodName(self):
950 """
951 returns the name of the solver currently used
952
953 @return: the name of the solver currently used.
954 @rtype: C{string}
955 """
956
957 m=self.getSolverMethod()
958 p=self.getSolverPackage()
959 method=""
960 if m[0]==self.DEFAULT: method="DEFAULT"
961 elif m[0]==self.DIRECT: method= "DIRECT"
962 elif m[0]==self.ITERATIVE: method= "ITERATIVE"
963 elif m[0]==self.CHOLEVSKY: method= "CHOLEVSKY"
964 elif m[0]==self.PCG: method= "PCG"
965 elif m[0]==self.CR: method= "CR"
966 elif m[0]==self.CGS: method= "CGS"
967 elif m[0]==self.BICGSTAB: method= "BICGSTAB"
968 elif m[0]==self.SSOR: method= "SSOR"
969 elif m[0]==self.GMRES: method= "GMRES"
970 elif m[0]==self.PRES20: method= "PRES20"
971 elif m[0]==self.LUMPING: method= "LUMPING"
972 elif m[0]==self.AMG: method= "AMG"
973 if m[1]==self.DEFAULT: method+="+DEFAULT"
974 elif m[1]==self.JACOBI: method+= "+JACOBI"
975 elif m[1]==self.ILU0: method+= "+ILU0"
976 elif m[1]==self.ILUT: method+= "+ILUT"
977 elif m[1]==self.SSOR: method+= "+SSOR"
978 elif m[1]==self.AMG: method+= "+AMG"
979 elif m[1]==self.RILU: method+= "+RILU"
980 if p==self.DEFAULT: package="DEFAULT"
981 elif p==self.PASO: package= "PASO"
982 elif p==self.MKL: package= "MKL"
983 elif p==self.SCSL: package= "SCSL"
984 elif p==self.UMFPACK: package= "UMFPACK"
985 elif p==self.TRILINOS: package= "TRILINOS"
986 else : method="unknown"
987 return "%s solver of %s package"%(method,package)
988
989
990 def getSolverMethod(self):
991 """
992 returns the solver method
993
994 @return: the solver method currently be used.
995 @rtype: C{int}
996 """
997 return self.__solver_method,self.__preconditioner
998
999 def setSolverPackage(self,package=None):
1000 """
1001 sets a new solver package
1002
1003 @param package: sets a new solver method.
1004 @type package: one of L{DEFAULT}, L{PASO} L{SCSL}, L{MKL}, L{UMFPACK}, L{TRILINOS}
1005 """
1006 if package==None: package=self.DEFAULT
1007 if not package==self.getSolverPackage():
1008 self.__solver_package=package
1009 self.__checkMatrixType()
1010 self.trace("New solver is %s"%self.getSolverMethodName())
1011
1012 def getSolverPackage(self):
1013 """
1014 returns the package of the solver
1015
1016 @return: the solver package currently being used.
1017 @rtype: C{int}
1018 """
1019 return self.__solver_package
1020
1021 def isUsingLumping(self):
1022 """
1023 checks if matrix lumping is used a solver method
1024
1025 @return: True is lumping is currently used a solver method.
1026 @rtype: C{bool}
1027 """
1028 return self.getSolverMethod()[0]==self.LUMPING
1029
1030 def setTolerance(self,tol=1.e-8):
1031 """
1032 resets the tolerance for the solver method to tol where for an appropriate norm M{|.|}
1033
1034 M{|L{getResidual}()|<tol*|L{getRightHandSide}()|}
1035
1036 defines the stopping criterion.
1037
1038 @param tol: new tolerance for the solver. If the tol is lower then the current tolerence
1039 the system will be resolved.
1040 @type tol: positive C{float}
1041 @raise ValueError: if tolerance is not positive.
1042 """
1043 if not tol>0:
1044 raise ValueError,"Tolerance as to be positive"
1045 if tol<self.getTolerance(): self.__invalidateSolution()
1046 self.trace("New tolerance %e"%tol)
1047 self.__tolerance=tol
1048 return
1049
1050 def getTolerance(self):
1051 """
1052 returns the tolerance set for the solution
1053
1054 @return: tolerance currently used.
1055 @rtype: C{float}
1056 """
1057 return self.__tolerance
1058
1059 # =============================================================================
1060 # symmetry flag:
1061 # =============================================================================
1062 def isSymmetric(self):
1063 """
1064 checks if symmetry is indicated.
1065
1066 @return: True is a symmetric PDE is indicated, otherwise False is returned
1067 @rtype: C{bool}
1068 """
1069 return self.__sym
1070
1071 def setSymmetryOn(self):
1072 """
1073 sets the symmetry flag.
1074 """
1075 if not self.isSymmetric():
1076 self.trace("PDE is set to be symmetric")
1077 self.__sym=True
1078 self.__checkMatrixType()
1079
1080 def setSymmetryOff(self):
1081 """
1082 removes the symmetry flag.
1083 """
1084 if self.isSymmetric():
1085 self.trace("PDE is set to be unsymmetric")
1086 self.__sym=False
1087 self.__checkMatrixType()
1088
1089 def setSymmetryTo(self,flag=False):
1090 """
1091 sets the symmetry flag to flag
1092
1093 @param flag: If flag, the symmetry flag is set otherwise the symmetry flag is released.
1094 @type flag: C{bool}
1095 """
1096 if flag:
1097 self.setSymmetryOn()
1098 else:
1099 self.setSymmetryOff()
1100
1101 # =============================================================================
1102 # function space handling for the equation as well as the solution
1103 # =============================================================================
1104 def setReducedOrderOn(self):
1105 """
1106 switches on reduced order for solution and equation representation
1107
1108 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1109 """
1110 self.setReducedOrderForSolutionOn()
1111 self.setReducedOrderForEquationOn()
1112
1113 def setReducedOrderOff(self):
1114 """
1115 switches off reduced order for solution and equation representation
1116
1117 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1118 """
1119 self.setReducedOrderForSolutionOff()
1120 self.setReducedOrderForEquationOff()
1121
1122 def setReducedOrderTo(self,flag=False):
1123 """
1124 sets order reduction for both solution and equation representation according to flag.
1125 @param flag: if flag is True, the order reduction is switched on for both solution and equation representation, otherwise or
1126 if flag is not present order reduction is switched off
1127 @type flag: C{bool}
1128 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1129 """
1130 self.setReducedOrderForSolutionTo(flag)
1131 self.setReducedOrderForEquationTo(flag)
1132
1133
1134 def setReducedOrderForSolutionOn(self):
1135 """
1136 switches on reduced order for solution representation
1137
1138 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1139 """
1140 if not self.__reduce_solution_order:
1141 if self.__altered_coefficients:
1142 raise RuntimeError,"order cannot be altered after coefficients have been defined."
1143 self.trace("Reduced order is used to solution representation.")
1144 self.__reduce_solution_order=True
1145 self.__resetSystem()
1146
1147 def setReducedOrderForSolutionOff(self):
1148 """
1149 switches off reduced order for solution representation
1150
1151 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1152 """
1153 if self.__reduce_solution_order:
1154 if self.__altered_coefficients:
1155 raise RuntimeError,"order cannot be altered after coefficients have been defined."
1156 self.trace("Full order is used to interpolate solution.")
1157 self.__reduce_solution_order=False
1158 self.__resetSystem()
1159
1160 def setReducedOrderForSolutionTo(self,flag=False):
1161 """
1162 sets order for test functions according to flag
1163
1164 @param flag: if flag is True, the order reduction is switched on for solution representation, otherwise or
1165 if flag is not present order reduction is switched off
1166 @type flag: C{bool}
1167 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1168 """
1169 if flag:
1170 self.setReducedOrderForSolutionOn()
1171 else:
1172 self.setReducedOrderForSolutionOff()
1173
1174 def setReducedOrderForEquationOn(self):
1175 """
1176 switches on reduced order for equation representation
1177
1178 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1179 """
1180 if not self.__reduce_equation_order:
1181 if self.__altered_coefficients:
1182 raise RuntimeError,"order cannot be altered after coefficients have been defined."
1183 self.trace("Reduced order is used for test functions.")
1184 self.__reduce_equation_order=True
1185 self.__resetSystem()
1186
1187 def setReducedOrderForEquationOff(self):
1188 """
1189 switches off reduced order for equation representation
1190
1191 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1192 """
1193 if self.__reduce_equation_order:
1194 if self.__altered_coefficients:
1195 raise RuntimeError,"order cannot be altered after coefficients have been defined."
1196 self.trace("Full order is used for test functions.")
1197 self.__reduce_equation_order=False
1198 self.__resetSystem()
1199
1200 def setReducedOrderForEquationTo(self,flag=False):
1201 """
1202 sets order for test functions according to flag
1203
1204 @param flag: if flag is True, the order reduction is switched on for equation representation, otherwise or
1205 if flag is not present order reduction is switched off
1206 @type flag: C{bool}
1207 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1208 """
1209 if flag:
1210 self.setReducedOrderForEquationOn()
1211 else:
1212 self.setReducedOrderForEquationOff()
1213
1214 # =============================================================================
1215 # private method:
1216 # =============================================================================
1217 def __checkMatrixType(self):
1218 """
1219 reassess the matrix type and, if a new matrix is needed, resets the system.
1220 """
1221 new_matrix_type=self.getDomain().getSystemMatrixTypeId(self.getSolverMethod()[0],self.getSolverPackage(),self.isSymmetric())
1222 if not new_matrix_type==self.__matrix_type:
1223 self.trace("Matrix type is now %d."%new_matrix_type)
1224 self.__matrix_type=new_matrix_type
1225 self.__resetSystem()
1226 #
1227 # rebuild switches :
1228 #
1229 def __invalidateSolution(self):
1230 """
1231 indicates the PDE has to be resolved if the solution is requested
1232 """
1233 if self.__solution_isValid: self.trace("PDE has to be resolved.")
1234 self.__solution_isValid=False
1235
1236 def __invalidateOperator(self):
1237 """
1238 indicates the operator has to be rebuilt next time it is used
1239 """
1240 if self.__operator_is_Valid: self.trace("Operator has to be rebuilt.")
1241 self.__invalidateSolution()
1242 self.__operator_is_Valid=False
1243
1244 def __invalidateRightHandSide(self):
1245 """
1246 indicates the right hand side has to be rebuild next time it is used
1247 """
1248 if self.__righthandside_isValid: self.trace("Right hand side has to be rebuilt.")
1249 self.__invalidateSolution()
1250 self.__righthandside_isValid=False
1251
1252 def __invalidateSystem(self):
1253 """
1254 annonced that everthing has to be rebuild:
1255 """
1256 if self.__righthandside_isValid: self.trace("System has to be rebuilt.")
1257 self.__invalidateSolution()
1258 self.__invalidateOperator()
1259 self.__invalidateRightHandSide()
1260
1261 def __resetSystem(self):
1262 """
1263 annonced that everthing has to be rebuild:
1264 """
1265 self.trace("New System is built from scratch.")
1266 self.__operator=escript.Operator()
1267 self.__operator_is_Valid=False
1268 self.__righthandside=escript.Data()
1269 self.__righthandside_isValid=False
1270 self.__solution=escript.Data()
1271 self.__solution_isValid=False
1272 #
1273 # system initialization:
1274 #
1275 def __getNewOperator(self):
1276 """
1277 returns an instance of a new operator
1278 """
1279 self.trace("New operator is allocated.")
1280 return self.getDomain().newOperator( \
1281 self.getNumEquations(), \
1282 self.getFunctionSpaceForEquation(), \
1283 self.getNumSolutions(), \
1284 self.getFunctionSpaceForSolution(), \
1285 self.__matrix_type)
1286
1287 def __getNewRightHandSide(self):
1288 """
1289 returns an instance of a new right hand side
1290 """
1291 self.trace("New right hand side is allocated.")
1292 if self.getNumEquations()>1:
1293 return escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)
1294 else:
1295 return escript.Data(0.,(),self.getFunctionSpaceForEquation(),True)
1296
1297 def __getNewSolution(self):
1298 """
1299 returns an instance of a new solution
1300 """
1301 self.trace("New solution is allocated.")
1302 if self.getNumSolutions()>1:
1303 return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)
1304 else:
1305 return escript.Data(0.,(),self.getFunctionSpaceForSolution(),True)
1306
1307 def __makeFreshSolution(self):
1308 """
1309 makes sure that the solution is instantiated and returns it initialized by zeros
1310 """
1311 if self.__solution.isEmpty():
1312 self.__solution=self.__getNewSolution()
1313 else:
1314 self.__solution*=0
1315 self.trace("Solution is reset to zero.")
1316 return self.__solution
1317
1318 def __makeFreshRightHandSide(self):
1319 """
1320 makes sure that the right hand side is instantiated and returns it initialized by zeros
1321 """
1322 if self.__righthandside.isEmpty():
1323 self.__righthandside=self.__getNewRightHandSide()
1324 else:
1325 self.__righthandside.setToZero()
1326 self.trace("Right hand side is reset to zero.")
1327 return self.__righthandside
1328
1329 def __makeFreshOperator(self):
1330 """
1331 makes sure that the operator is instantiated and returns it initialized by zeros
1332 """
1333 if self.__operator.isEmpty():
1334 self.__operator=self.__getNewOperator()
1335 else:
1336 self.__operator.resetValues()
1337 self.trace("Operator reset to zero")
1338 return self.__operator
1339
1340 def __applyConstraint(self):
1341 """
1342 applies the constraints defined by q and r to the system
1343 """
1344 if not self.isUsingLumping():
1345 q=self.getCoefficientOfGeneralPDE("q")
1346 r=self.getCoefficientOfGeneralPDE("r")
1347 if not q.isEmpty() and not self.__operator.isEmpty():
1348 # q is the row and column mask to indicate where constraints are set:
1349 row_q=escript.Data(q,self.getFunctionSpaceForEquation())
1350 col_q=escript.Data(q,self.getFunctionSpaceForSolution())
1351 u=self.__getNewSolution()
1352 if r.isEmpty():
1353 r_s=self.__getNewSolution()
1354 else:
1355 r_s=escript.Data(r,self.getFunctionSpaceForSolution())
1356 u.copyWithMask(r_s,col_q)
1357 if not self.__righthandside.isEmpty():
1358 self.__righthandside-=self.__operator*u
1359 self.__righthandside=self.copyConstraint(self.__righthandside)
1360 self.__operator.nullifyRowsAndCols(row_q,col_q,1.)
1361 # =============================================================================
1362 # function giving access to coefficients of the general PDE:
1363 # =============================================================================
1364 def getCoefficientOfGeneralPDE(self,name):
1365 """
1366 return the value of the coefficient name of the general PDE.
1367
1368 @note: This method is called by the assembling routine it can be overwritten
1369 to map coefficients of a particular PDE to the general PDE.
1370 @param name: name of the coefficient requested.
1371 @type name: C{string}
1372 @return: the value of the coefficient name
1373 @rtype: L{Data<escript.Data>}
1374 @raise IllegalCoefficient: if name is not one of coefficients
1375 M{A}, M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact},
1376 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}.
1377 """
1378 if self.hasCoefficientOfGeneralPDE(name):
1379 return self.getCoefficient(name)
1380 else:
1381 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1382
1383 def hasCoefficientOfGeneralPDE(self,name):
1384 """
1385 checks if name is a the name of a coefficient of the general PDE.
1386
1387 @param name: name of the coefficient enquired.
1388 @type name: C{string}
1389 @return: True if name is the name of a coefficient of the general PDE. Otherwise False.
1390 @rtype: C{bool}
1391
1392 """
1393 return self.__COEFFICIENTS_OF_GENEARL_PDE.has_key(name)
1394
1395 def createCoefficientOfGeneralPDE(self,name):
1396 """
1397 returns a new instance of a coefficient for coefficient name of the general PDE
1398
1399 @param name: name of the coefficient requested.
1400 @type name: C{string}
1401 @return: a coefficient name initialized to 0.
1402 @rtype: L{Data<escript.Data>}
1403 @raise IllegalCoefficient: if name is not one of coefficients
1404 M{A}, M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact},
1405 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}.
1406 """
1407 if self.hasCoefficientOfGeneralPDE(name):
1408 return escript.Data(0,self.getShapeOfCoefficientOfGeneralPDE(name),self.getFunctionSpaceForCoefficientOfGeneralPDE(name))
1409 else:
1410 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1411
1412 def getFunctionSpaceForCoefficientOfGeneralPDE(self,name):
1413 """
1414 return the L{FunctionSpace<escript.FunctionSpace>} to be used for coefficient name of the general PDE
1415
1416 @param name: name of the coefficient enquired.
1417 @type name: C{string}
1418 @return: the function space to be used for coefficient name
1419 @rtype: L{FunctionSpace<escript.FunctionSpace>}
1420 @raise IllegalCoefficient: if name is not one of coefficients
1421 M{A}, M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact},
1422 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}.
1423 """
1424 if self.hasCoefficientOfGeneralPDE(name):
1425 return self.__COEFFICIENTS_OF_GENEARL_PDE[name].getFunctionSpace(self.getDomain())
1426 else:
1427 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1428
1429 def getShapeOfCoefficientOfGeneralPDE(self,name):
1430 """
1431 return the shape of the coefficient name of the general PDE
1432
1433 @param name: name of the coefficient enquired.
1434 @type name: C{string}
1435 @return: the shape of the coefficient name
1436 @rtype: C{tuple} of C{int}
1437 @raise IllegalCoefficient: if name is not one of coefficients
1438 M{A}, M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact},
1439 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}.
1440 """
1441 if self.hasCoefficientOfGeneralPDE(name):
1442 return self.__COEFFICIENTS_OF_GENEARL_PDE[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())
1443 else:
1444 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1445
1446 # =============================================================================
1447 # functions giving access to coefficients of a particular PDE implementation:
1448 # =============================================================================
1449 def getCoefficient(self,name):
1450 """
1451 returns the value of the coefficient name
1452
1453 @param name: name of the coefficient requested.
1454 @type name: C{string}
1455 @return: the value of the coefficient name
1456 @rtype: L{Data<escript.Data>}
1457 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1458 """
1459 if self.hasCoefficient(name):
1460 return self.COEFFICIENTS[name].getValue()
1461 else:
1462 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1463
1464 def hasCoefficient(self,name):
1465 """
1466 return True if name is the name of a coefficient
1467
1468 @param name: name of the coefficient enquired.
1469 @type name: C{string}
1470 @return: True if name is the name of a coefficient of the general PDE. Otherwise False.
1471 @rtype: C{bool}
1472 """
1473 return self.COEFFICIENTS.has_key(name)
1474
1475 def createCoefficient(self, name):
1476 """
1477 create a L{Data<escript.Data>} object corresponding to coefficient name
1478
1479 @return: a coefficient name initialized to 0.
1480 @rtype: L{Data<escript.Data>}
1481 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1482 """
1483 if self.hasCoefficient(name):
1484 return escript.Data(0.,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))
1485 else:
1486 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1487
1488 def getFunctionSpaceForCoefficient(self,name):
1489 """
1490 return the L{FunctionSpace<escript.FunctionSpace>} to be used for coefficient name
1491
1492 @param name: name of the coefficient enquired.
1493 @type name: C{string}
1494 @return: the function space to be used for coefficient name
1495 @rtype: L{FunctionSpace<escript.FunctionSpace>}
1496 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1497 """
1498 if self.hasCoefficient(name):
1499 return self.COEFFICIENTS[name].getFunctionSpace(self.getDomain())
1500 else:
1501 raise ValueError,"unknown coefficient %s requested"%name
1502 def getShapeOfCoefficient(self,name):
1503 """
1504 return the shape of the coefficient name
1505
1506 @param name: name of the coefficient enquired.
1507 @type name: C{string}
1508 @return: the shape of the coefficient name
1509 @rtype: C{tuple} of C{int}
1510 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1511 """
1512 if self.hasCoefficient(name):
1513 return self.COEFFICIENTS[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())
1514 else:
1515 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1516
1517 def resetCoefficients(self):
1518 """
1519 resets all coefficients to there default values.
1520 """
1521 for i in self.COEFFICIENTS.iterkeys():
1522 self.COEFFICIENTS[i].resetValue()
1523
1524 def alteredCoefficient(self,name):
1525 """
1526 announce that coefficient name has been changed
1527
1528 @param name: name of the coefficient enquired.
1529 @type name: C{string}
1530 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1531 @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.
1532 """
1533 if self.hasCoefficient(name):
1534 self.trace("Coefficient %s has been altered."%name)
1535 if not ((name=="q" or name=="r") and self.isUsingLumping()):
1536 if self.COEFFICIENTS[name].isAlteringOperator(): self.__invalidateOperator()
1537 if self.COEFFICIENTS[name].isAlteringRightHandSide(): self.__invalidateRightHandSide()
1538 else:
1539 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1540
1541 def copyConstraint(self,u):
1542 """
1543 copies the constraint into u and returns u.
1544
1545 @param u: a function of rank 0 is a single PDE is solved and of shape (numSolution,) for a system of PDEs
1546 @type u: L{Data<escript.Data>}
1547 @return: the input u modified by the constraints.
1548 @rtype: L{Data<escript.Data>}
1549 @warning: u is altered if it has the appropriate L{FunctionSpace<escript.FunctionSpace>}
1550 """
1551 q=self.getCoefficientOfGeneralPDE("q")
1552 r=self.getCoefficientOfGeneralPDE("r")
1553 if not q.isEmpty():
1554 if u.isEmpty(): u=escript.Data(0.,q.getShape(),q.getFunctionSpace())
1555 if r.isEmpty():
1556 r=escript.Data(0,u.getShape(),u.getFunctionSpace())
1557 else:
1558 r=escript.Data(r,u.getFunctionSpace())
1559 u.copyWithMask(r,escript.Data(q,u.getFunctionSpace()))
1560 return u
1561
1562 def setValue(self,**coefficients):
1563 """
1564 sets new values to coefficients
1565
1566 @param coefficients: new values assigned to coefficients
1567 @keyword A: value for coefficient A.
1568 @type A: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1569 @keyword A_reduced: value for coefficient A_reduced.
1570 @type A_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1571 @keyword B: value for coefficient B
1572 @type B: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1573 @keyword B_reduced: value for coefficient B_reduced
1574 @type B_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1575 @keyword C: value for coefficient C
1576 @type C: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1577 @keyword C_reduced: value for coefficient C_reduced
1578 @type C_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1579 @keyword D: value for coefficient D
1580 @type D: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1581 @keyword D_reduced: value for coefficient D_reduced
1582 @type D_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1583 @keyword X: value for coefficient X
1584 @type X: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1585 @keyword X_reduced: value for coefficient X_reduced
1586 @type X_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1587 @keyword Y: value for coefficient Y
1588 @type Y: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1589 @keyword Y_reduced: value for coefficient Y_reduced
1590 @type Y_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.Function>}.
1591 @keyword d: value for coefficient d
1592 @type d: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1593 @keyword d_reduced: value for coefficient d_reduced
1594 @type d_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
1595 @keyword y: value for coefficient y
1596 @type y: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1597 @keyword d_contact: value for coefficient d_contact
1598 @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>}.
1599 @keyword d_contact_reduced: value for coefficient d_contact_reduced
1600 @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>}.
1601 @keyword y_contact: value for coefficient y_contact
1602 @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>}.
1603 @keyword y_contact_reduced: value for coefficient y_contact_reduced
1604 @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>}.
1605 @keyword r: values prescribed to the solution at the locations of constraints
1606 @type r: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1607 depending of reduced order is used for the solution.
1608 @keyword q: mask for location of constraints
1609 @type q: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1610 depending of reduced order is used for the representation of the equation.
1611 @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1612 """
1613 # check if the coefficients are legal:
1614 for i in coefficients.iterkeys():
1615 if not self.hasCoefficient(i):
1616 raise IllegalCoefficient,"Attempt to set unknown coefficient %s"%i
1617 # if the number of unknowns or equations is still unknown we try to estimate them:
1618 if self.__numEquations==None or self.__numSolutions==None:
1619 for i,d in coefficients.iteritems():
1620 if hasattr(d,"shape"):
1621 s=d.shape
1622 elif hasattr(d,"getShape"):
1623 s=d.getShape()
1624 else:
1625 s=numarray.array(d).shape
1626 if s!=None:
1627 # get number of equations and number of unknowns:
1628 res=self.COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(self.getDomain(),s)
1629 if res==None:
1630 raise IllegalCoefficientValue,"Illegal shape %s of coefficient %s"%(s,i)
1631 else:
1632 if self.__numEquations==None: self.__numEquations=res[0]
1633 if self.__numSolutions==None: self.__numSolutions=res[1]
1634 if self.__numEquations==None: raise UndefinedPDEError,"unidententified number of equations"
1635 if self.__numSolutions==None: raise UndefinedPDEError,"unidententified number of solutions"
1636 # now we check the shape of the coefficient if numEquations and numSolutions are set:
1637 for i,d in coefficients.iteritems():
1638 try:
1639 self.COEFFICIENTS[i].setValue(self.getDomain(),
1640 self.getNumEquations(),self.getNumSolutions(),
1641 self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
1642 self.alteredCoefficient(i)
1643 except IllegalCoefficientFunctionSpace,m:
1644 # if the function space is wrong then we try the reduced version:
1645 i_red=i+"_reduced"
1646 if (not i_red in coefficients.keys()) and i_red in self.COEFFICIENTS.keys():
1647 try:
1648 self.COEFFICIENTS[i_red].setValue(self.getDomain(),
1649 self.getNumEquations(),self.getNumSolutions(),
1650 self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
1651 self.alteredCoefficient(i_red)
1652 except IllegalCoefficientValue,m:
1653 raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
1654 except IllegalCoefficientFunctionSpace,m:
1655 raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m))
1656 else:
1657 raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m))
1658 except IllegalCoefficientValue,m:
1659 raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
1660 self.__altered_coefficients=True
1661 # check if the systrem is inhomogeneous:
1662 if len(coefficients)>0 and not self.isUsingLumping():
1663 q=self.getCoefficientOfGeneralPDE("q")
1664 r=self.getCoefficientOfGeneralPDE("r")
1665 homogeneous_constraint=True
1666 if not q.isEmpty() and not r.isEmpty():
1667 if util.Lsup(q*r)>0.:
1668 self.trace("Inhomogeneous constraint detected.")
1669 self.__invalidateSystem()
1670
1671 def getSystem(self):
1672 """
1673 return the operator and right hand side of the PDE
1674
1675 @return: the discrete version of the PDE
1676 @rtype: C{tuple} of L{Operator,<escript.Operator>} and L{Data<escript.Data>}.
1677 """
1678 if not self.__operator_is_Valid or not self.__righthandside_isValid:
1679 if self.isUsingLumping():
1680 if not self.__operator_is_Valid:
1681 if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():
1682 raise TypeError,"Lumped matrix requires same order for equations and unknowns"
1683 if not self.getCoefficientOfGeneralPDE("A").isEmpty():
1684 raise ValueError,"coefficient A in lumped matrix may not be present."
1685 if not self.getCoefficientOfGeneralPDE("B").isEmpty():
1686 raise ValueError,"coefficient B in lumped matrix may not be present."
1687 if not self.getCoefficientOfGeneralPDE("C").isEmpty():
1688 raise ValueError,"coefficient C in lumped matrix may not be present."
1689 if not self.getCoefficientOfGeneralPDE("d_contact").isEmpty():
1690 raise ValueError,"coefficient d_contact in lumped matrix may not be present."
1691 if not self.getCoefficientOfGeneralPDE("A_reduced").isEmpty():
1692 raise ValueError,"coefficient A_reduced in lumped matrix may not be present."
1693 if not self.getCoefficientOfGeneralPDE("B_reduced").isEmpty():
1694 raise ValueError,"coefficient B_reduced in lumped matrix may not be present."
1695 if not self.getCoefficientOfGeneralPDE("C_reduced").isEmpty():
1696 raise ValueError,"coefficient C_reduced in lumped matrix may not be present."
1697 if not self.getCoefficientOfGeneralPDE("d_contact_reduced").isEmpty():
1698 raise ValueError,"coefficient d_contact_reduced in lumped matrix may not be present."
1699 D=self.getCoefficientOfGeneralPDE("D")
1700 d=self.getCoefficientOfGeneralPDE("d")
1701 D_reduced=self.getCoefficientOfGeneralPDE("D_reduced")
1702 d_reduced=self.getCoefficientOfGeneralPDE("d_reduced")
1703 if not D.isEmpty():
1704 if self.getNumSolutions()>1:
1705 D_times_e=util.matrix_mult(D,numarray.ones((self.getNumSolutions(),)))
1706 else:
1707 D_times_e=D
1708 else:
1709 D_times_e=escript.Data()
1710 if not d.isEmpty():
1711 if self.getNumSolutions()>1:
1712 d_times_e=util.matrix_mult(d,numarray.ones((self.getNumSolutions(),)))
1713 else:
1714 d_times_e=d
1715 else:
1716 d_times_e=escript.Data()
1717
1718 if not D_reduced.isEmpty():
1719 if self.getNumSolutions()>1:
1720 D_reduced_times_e=util.matrix_mult(D_reduced,numarray.ones((self.getNumSolutions(),)))
1721 else:
1722 D_reduced_times_e=D_reduced
1723 else:
1724 D_reduced_times_e=escript.Data()
1725 if not d_reduced.isEmpty():
1726 if self.getNumSolutions()>1:
1727 d_reduced_times_e=util.matrix_mult(d_reduced,numarray.ones((self.getNumSolutions(),)))
1728 else:
1729 d_reduced_times_e=d_reduced
1730 else:
1731 d_reduced_times_e=escript.Data()
1732
1733 self.__operator=self.__getNewRightHandSide()
1734 if hasattr(self.getDomain(), "addPDEToLumpedSystem") :
1735 self.getDomain().addPDEToLumpedSystem(self.__operator, D_times_e, d_times_e)
1736 self.getDomain().addPDEToLumpedSystem(self.__operator, D_reduced_times_e, d_reduced_times_e)
1737 else:
1738 self.getDomain().addPDEToRHS(self.__operator, \
1739 escript.Data(), \
1740 D_times_e, \
1741 d_times_e,\
1742 escript.Data())
1743 self.getDomain().addPDEToRHS(self.__operator, \
1744 escript.Data(), \
1745 D_reduced_times_e, \
1746 d_reduced_times_e,\
1747 escript.Data())
1748 self.__operator=1./self.__operator
1749 self.trace("New lumped operator has been built.")
1750 self.__operator_is_Valid=True
1751 if not self.__righthandside_isValid:
1752 self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \
1753 self.getCoefficientOfGeneralPDE("X"), \
1754 self.getCoefficientOfGeneralPDE("Y"),\
1755 self.getCoefficientOfGeneralPDE("y"),\
1756 self.getCoefficientOfGeneralPDE("y_contact"))
1757 self.getDomain().addPDEToRHS(self.__righthandside, \
1758 self.getCoefficientOfGeneralPDE("X_reduced"), \
1759 self.getCoefficientOfGeneralPDE("Y_reduced"),\
1760 self.getCoefficientOfGeneralPDE("y_reduced"),\
1761 self.getCoefficientOfGeneralPDE("y_contact_reduced"))
1762 self.trace("New right hand side as been built.")
1763 self.__righthandside_isValid=True
1764 else:
1765 if not self.__operator_is_Valid and not self.__righthandside_isValid:
1766 self.getDomain().addPDEToSystem(self.__makeFreshOperator(),self.__makeFreshRightHandSide(), \
1767 self.getCoefficientOfGeneralPDE("A"), \
1768 self.getCoefficientOfGeneralPDE("B"), \
1769 self.getCoefficientOfGeneralPDE("C"), \
1770 self.getCoefficientOfGeneralPDE("D"), \
1771 self.getCoefficientOfGeneralPDE("X"), \
1772 self.getCoefficientOfGeneralPDE("Y"), \
1773 self.getCoefficientOfGeneralPDE("d"), \
1774 self.getCoefficientOfGeneralPDE("y"), \
1775 self.getCoefficientOfGeneralPDE("d_contact"), \
1776 self.getCoefficientOfGeneralPDE("y_contact"))
1777 self.getDomain().addPDEToSystem(self.__operator,self.__righthandside, \
1778 self.getCoefficientOfGeneralPDE("A_reduced"), \
1779 self.getCoefficientOfGeneralPDE("B_reduced"), \
1780 self.getCoefficientOfGeneralPDE("C_reduced"), \
1781 self.getCoefficientOfGeneralPDE("D_reduced"), \
1782 self.getCoefficientOfGeneralPDE("X_reduced"), \
1783 self.getCoefficientOfGeneralPDE("Y_reduced"), \
1784 self.getCoefficientOfGeneralPDE("d_reduced"), \
1785 self.getCoefficientOfGeneralPDE("y_reduced"), \
1786 self.getCoefficientOfGeneralPDE("d_contact_reduced"), \
1787 self.getCoefficientOfGeneralPDE("y_contact_reduced"))
1788 self.__applyConstraint()
1789 self.__righthandside=self.copyConstraint(self.__righthandside)
1790 self.trace("New system has been built.")
1791 self.__operator_is_Valid=True
1792 self.__righthandside_isValid=True
1793 elif not self.__righthandside_isValid:
1794 self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \
1795 self.getCoefficientOfGeneralPDE("X"), \
1796 self.getCoefficientOfGeneralPDE("Y"),\
1797 self.getCoefficientOfGeneralPDE("y"),\
1798 self.getCoefficientOfGeneralPDE("y_contact"))
1799 self.getDomain().addPDEToRHS(self.__righthandside, \
1800 self.getCoefficientOfGeneralPDE("X_reduced"), \
1801 self.getCoefficientOfGeneralPDE("Y_reduced"),\
1802 self.getCoefficientOfGeneralPDE("y_reduced"),\
1803 self.getCoefficientOfGeneralPDE("y_contact_reduced"))
1804 self.__righthandside=self.copyConstraint(self.__righthandside)
1805 self.trace("New right hand side has been built.")
1806 self.__righthandside_isValid=True
1807 elif not self.__operator_is_Valid:
1808 self.getDomain().addPDEToSystem(self.__makeFreshOperator(),escript.Data(), \
1809 self.getCoefficientOfGeneralPDE("A"), \
1810 self.getCoefficientOfGeneralPDE("B"), \
1811 self.getCoefficientOfGeneralPDE("C"), \
1812 self.getCoefficientOfGeneralPDE("D"), \
1813 escript.Data(), \
1814 escript.Data(), \
1815 self.getCoefficientOfGeneralPDE("d"), \
1816 escript.Data(),\
1817 self.getCoefficientOfGeneralPDE("d_contact"), \
1818 escript.Data())
1819 self.getDomain().addPDEToSystem(self.__operator,escript.Data(), \
1820 self.getCoefficientOfGeneralPDE("A_reduced"), \
1821 self.getCoefficientOfGeneralPDE("B_reduced"), \
1822 self.getCoefficientOfGeneralPDE("C_reduced"), \
1823 self.getCoefficientOfGeneralPDE("D_reduced"), \
1824 escript.Data(), \
1825 escript.Data(), \
1826 self.getCoefficientOfGeneralPDE("d_reduced"), \
1827 escript.Data(),\
1828 self.getCoefficientOfGeneralPDE("d_contact_reduced"), \
1829 escript.Data())
1830 self.__applyConstraint()
1831 self.trace("New operator has been built.")
1832 self.__operator_is_Valid=True
1833 return (self.__operator,self.__righthandside)
1834
1835
1836 class Poisson(LinearPDE):
1837 """
1838 Class to define a Poisson equation problem, which is genear L{LinearPDE} of the form
1839
1840 M{-grad(grad(u)[j])[j] = f}
1841
1842 with natural boundary conditons
1843
1844 M{n[j]*grad(u)[j] = 0 }
1845
1846 and constraints:
1847
1848 M{u=0} where M{q>0}
1849
1850 """
1851
1852 def __init__(self,domain,debug=False):
1853 """
1854 initializes a new Poisson equation
1855
1856 @param domain: domain of the PDE
1857 @type domain: L{Domain<escript.Domain>}
1858 @param debug: if True debug informations are printed.
1859
1860 """
1861 super(Poisson, self).__init__(domain,1,1,debug)
1862 self.COEFFICIENTS={"f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1863 "f_reduced": PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1864 "q": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
1865 self.setSymmetryOn()
1866
1867 def setValue(self,**coefficients):
1868 """
1869 sets new values to coefficients
1870
1871 @param coefficients: new values assigned to coefficients
1872 @keyword f: value for right hand side M{f}
1873 @type f: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1874 @keyword q: mask for location of constraints
1875 @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>}
1876 depending of reduced order is used for the representation of the equation.
1877 @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1878 """
1879 super(Poisson, self).setValue(**coefficients)
1880
1881 def getCoefficientOfGeneralPDE(self,name):
1882 """
1883 return the value of the coefficient name of the general PDE
1884 @param name: name of the coefficient requested.
1885 @type name: C{string}
1886 @return: the value of the coefficient name
1887 @rtype: L{Data<escript.Data>}
1888 @raise IllegalCoefficient: if name is not one of coefficients
1889 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}.
1890 @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
1891 """
1892 if name == "A" :
1893 return escript.Data(util.kronecker(self.getDim()),escript.Function(self.getDomain()))
1894 elif name == "B" :
1895 return escript.Data()
1896 elif name == "C" :
1897 return escript.Data()
1898 elif name == "D" :
1899 return escript.Data()
1900 elif name == "X" :
1901 return escript.Data()
1902 elif name == "Y" :
1903 return self.getCoefficient("f")
1904 elif name == "d" :
1905 return escript.Data()
1906 elif name == "y" :
1907 return escript.Data()
1908 elif name == "d_contact" :
1909 return escript.Data()
1910 elif name == "y_contact" :
1911 return escript.Data()
1912 elif name == "A_reduced" :
1913 return escript.Data()
1914 elif name == "B_reduced" :
1915 return escript.Data()
1916 elif name == "C_reduced" :
1917 return escript.Data()
1918 elif name == "D_reduced" :
1919 return escript.Data()
1920 elif name == "X_reduced" :
1921 return escript.Data()
1922 elif name == "Y_reduced" :
1923 return self.getCoefficient("f_reduced")
1924 elif name == "d_reduced" :
1925 return escript.Data()
1926 elif name == "y_reduced" :
1927 return escript.Data()
1928 elif name == "d_contact_reduced" :
1929 return escript.Data()
1930 elif name == "y_contact_reduced" :
1931 return escript.Data()
1932 elif name == "r" :
1933 return escript.Data()
1934 elif name == "q" :
1935 return self.getCoefficient("q")
1936 else:
1937 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1938
1939 class Helmholtz(LinearPDE):
1940 """
1941 Class to define a Helmhotz equation problem, which is genear L{LinearPDE} of the form
1942
1943 M{S{omega}*u - grad(k*grad(u)[j])[j] = f}
1944
1945 with natural boundary conditons
1946
1947 M{k*n[j]*grad(u)[j] = g- S{alpha}u }
1948
1949 and constraints:
1950
1951 M{u=r} where M{q>0}
1952
1953 """
1954
1955 def __init__(self,domain,debug=False):
1956 """
1957 initializes a new Poisson equation
1958
1959 @param domain: domain of the PDE
1960 @type domain: L{Domain<escript.Domain>}
1961 @param debug: if True debug informations are printed.
1962
1963 """
1964 super(Helmholtz, self).__init__(domain,1,1,debug)
1965 self.COEFFICIENTS={"omega": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
1966 "k": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
1967 "f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1968 "f_reduced": PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1969 "alpha": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
1970 "g": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1971 "g_reduced": PDECoefficient(PDECoefficient.BOUNDARY_REDUCED,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1972 "r": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH),
1973 "q": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
1974 self.setSymmetryOn()
1975
1976 def setValue(self,**coefficients):
1977 """
1978 sets new values to coefficients
1979
1980 @param coefficients: new values assigned to coefficients
1981 @keyword omega: value for coefficient M{S{omega}}
1982 @type omega: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1983 @keyword k: value for coefficeint M{k}
1984 @type k: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1985 @keyword f: value for right hand side M{f}
1986 @type f: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1987 @keyword alpha: value for right hand side M{S{alpha}}
1988 @type alpha: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1989 @keyword g: value for right hand side M{g}
1990 @type g: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1991 @keyword r: prescribed values M{r} for the solution in constraints.
1992 @type r: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1993 depending of reduced order is used for the representation of the equation.
1994 @keyword q: mask for location of constraints
1995 @type q: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1996 depending of reduced order is used for the representation of the equation.
1997 @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1998 """
1999 super(Helmholtz, self).setValue(**coefficients)
2000
2001 def getCoefficientOfGeneralPDE(self,name):
2002 """
2003 return the value of the coefficient name of the general PDE
2004
2005 @param name: name of the coefficient requested.
2006 @type name: C{string}
2007 @return: the value of the coefficient name
2008 @rtype: L{Data<escript.Data>}
2009 @raise IllegalCoefficient: if name is not one of coefficients
2010 "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}.
2011 @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
2012 """
2013 if name == "A" :
2014 return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))*self.getCoefficient("k")
2015 elif name == "B" :
2016 return escript.Data()
2017 elif name == "C" :
2018 return escript.Data()
2019 elif name == "D" :
2020 return self.getCoefficient("omega")
2021 elif name == "X" :
2022 return escript.Data()
2023 elif name == "Y" :
2024 return self.getCoefficient("f")
2025 elif name == "d" :
2026 return self.getCoefficient("alpha")
2027 elif name == "y" :
2028 return self.getCoefficient("g")
2029 elif name == "d_contact" :
2030 return escript.Data()
2031 elif name == "y_contact" :
2032 return escript.Data()
2033 elif name == "A_reduced" :
2034 return escript.Data()
2035 elif name == "B_reduced" :
2036 return escript.Data()
2037 elif name == "C_reduced" :
2038 return escript.Data()
2039 elif name == "D_reduced" :
2040 return escript.Data()
2041 elif name == "X_reduced" :
2042 return escript.Data()
2043 elif name == "Y_reduced" :
2044 return self.getCoefficient("f_reduced")
2045 elif name == "d_reduced" :
2046 return escript.Data()
2047 elif name == "y_reduced" :
2048 return self.getCoefficient("g_reduced")
2049 elif name == "d_contact_reduced" :
2050 return escript.Data()
2051 elif name == "y_contact_reduced" :
2052 return escript.Data()
2053 elif name == "r" :
2054 return self.getCoefficient("r")
2055 elif name == "q" :
2056 return self.getCoefficient("q")
2057 else:
2058 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
2059
2060 class LameEquation(LinearPDE):
2061 """
2062 Class to define a Lame equation problem:
2063
2064 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] }
2065
2066 with natural boundary conditons:
2067
2068 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] }
2069
2070 and constraints:
2071
2072 M{u[i]=r[i]} where M{q[i]>0}
2073
2074 """
2075
2076 def __init__(self,domain,debug=False):
2077 super(LameEquation, self).__init__(domain,\
2078 domain.getDim(),domain.getDim(),debug)
2079 self.COEFFICIENTS={ "lame_lambda" : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),
2080 "lame_mu" : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),
2081 "F" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
2082 "sigma" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM),PDECoefficient.RIGHTHANDSIDE),
2083 "f" : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
2084 "r" : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH),
2085 "q" : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
2086 self.setSymmetryOn()
2087
2088 def setValues(self,**coefficients):
2089 """
2090 sets new values to coefficients
2091
2092 @param coefficients: new values assigned to coefficients
2093 @keyword lame_mu: value for coefficient M{S{mu}}
2094 @type lame_mu: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2095 @keyword lame_lambda: value for coefficient M{S{lambda}}
2096 @type lame_lambda: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2097 @keyword F: value for internal force M{F}
2098 @type F: any type that can be casted to L{Vector<escript.Vector>} object on L{Function<escript.Function>}
2099 @keyword sigma: value for initial stress M{S{sigma}}
2100 @type sigma: any type that can be casted to L{Tensor<escript.Tensor>} object on L{Function<escript.Function>}
2101 @keyword f: value for extrenal force M{f}
2102 @type f: any type that can be casted to L{Vector<escript.Vector>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}
2103 @keyword r: prescribed values M{r} for the solution in constraints.
2104 @type r: any type that can be casted to L{Vector<escript.Vector>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2105 depending of reduced order is used for the representation of the equation.
2106 @keyword q: mask for location of constraints
2107 @type q: any type that can be casted to L{Vector<escript.Vector>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2108 depending of reduced order is used for the representation of the equation.
2109 @raise IllegalCoefficient: if an unknown coefficient keyword is used.
2110 """
2111 super(LameEquation, self).setValues(**coefficients)
2112
2113 def getCoefficientOfGeneralPDE(self,name):
2114 """
2115 return the value of the coefficient name of the general PDE
2116
2117 @param name: name of the coefficient requested.
2118 @type name: C{string}
2119 @return: the value of the coefficient name
2120 @rtype: L{Data<escript.Data>}
2121 @raise IllegalCoefficient: if name is not one of coefficients
2122 "A", M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact}, M{r} or M{q}.
2123 @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
2124 """
2125 if name == "A" :
2126 out =self.createCoefficientOfGeneralPDE("A")
2127 for i in range(self.getDim()):
2128 for j in range(self.getDim()):
2129 out[i,i,j,j] += self.getCoefficient("lame_lambda")
2130 out[i,j,j,i] += self.getCoefficient("lame_mu")
2131 out[i,j,i,j] += self.getCoefficient("lame_mu")
2132 return out
2133 elif name == "B" :
2134 return escript.Data()
2135 elif name == "C" :
2136 return escript.Data()
2137 elif name == "D" :
2138 return escript.Data()
2139 elif name == "X" :
2140 return self.getCoefficient("sigma")
2141 elif name == "Y" :
2142 return self.getCoefficient("F")
2143 elif name == "d" :
2144 return escript.Data()
2145 elif name == "y" :
2146 return self.getCoefficient("f")
2147 elif name == "d_contact" :
2148 return escript.Data()
2149 elif name == "y_contact" :
2150 return escript.Data()
2151 elif name == "A_reduced" :
2152 return escript.Data()
2153 elif name == "B_reduced" :
2154 return escript.Data()
2155 elif name == "C_reduced" :
2156 return escript.Data()
2157 elif name == "D_reduced" :
2158 return escript.Data()
2159 elif name == "X_reduced" :
2160 return escript.Data()
2161 elif name == "Y_reduced" :
2162 return escript.Data()
2163 elif name == "d_reduced" :
2164 return escript.Data()
2165 elif name == "y_reduced" :
2166 return escript.Data()
2167 elif name == "d_contact_reduced" :
2168 return escript.Data()
2169 elif name == "y_contact_reduced" :
2170 return escript.Data()
2171 elif name == "r" :
2172 return self.getCoefficient("r")
2173 elif name == "q" :
2174 return self.getCoefficient("q")
2175 else:
2176 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
2177

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26