/[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 1223 - (show annotations)
Fri Aug 3 02:40:39 2007 UTC (11 years, 8 months ago) by gross
File MIME type: text/x-python
File size: 102166 byte(s)
first attemt towards an improved MPI version.  

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: solve=self.DEFAULT
940 if preconditioner==None: preconditioner=self.DEFAULT
941 if not (solver,preconditioner)==self.getSolverMethod():
942 self.__solver_method=solver
943 self.__preconditioner=preconditioner
944 self.__checkMatrixType()
945 self.trace("New solver is %s"%self.getSolverMethodName())
946
947 def getSolverMethodName(self):
948 """
949 returns the name of the solver currently used
950
951 @return: the name of the solver currently used.
952 @rtype: C{string}
953 """
954
955 m=self.getSolverMethod()
956 p=self.getSolverPackage()
957 method=""
958 if m[0]==self.DEFAULT: method="DEFAULT"
959 elif m[0]==self.DIRECT: method= "DIRECT"
960 elif m[0]==self.ITERATIVE: method= "ITERATIVE"
961 elif m[0]==self.CHOLEVSKY: method= "CHOLEVSKY"
962 elif m[0]==self.PCG: method= "PCG"
963 elif m[0]==self.CR: method= "CR"
964 elif m[0]==self.CGS: method= "CGS"
965 elif m[0]==self.BICGSTAB: method= "BICGSTAB"
966 elif m[0]==self.SSOR: method= "SSOR"
967 elif m[0]==self.GMRES: method= "GMRES"
968 elif m[0]==self.PRES20: method= "PRES20"
969 elif m[0]==self.LUMPING: method= "LUMPING"
970 elif m[0]==self.AMG: method= "AMG"
971 if m[1]==self.DEFAULT: method+="+DEFAULT"
972 elif m[1]==self.JACOBI: method+= "+JACOBI"
973 elif m[1]==self.ILU0: method+= "+ILU0"
974 elif m[1]==self.ILUT: method+= "+ILUT"
975 elif m[1]==self.SSOR: method+= "+SSOR"
976 elif m[1]==self.AMG: method+= "+AMG"
977 elif m[1]==self.RILU: method+= "+RILU"
978 if p==self.DEFAULT: package="DEFAULT"
979 elif p==self.PASO: package= "PASO"
980 elif p==self.MKL: package= "MKL"
981 elif p==self.SCSL: package= "SCSL"
982 elif p==self.UMFPACK: package= "UMFPACK"
983 elif p==self.TRILINOS: package= "TRILINOS"
984 else : method="unknown"
985 return "%s solver of %s package"%(method,package)
986
987
988 def getSolverMethod(self):
989 """
990 returns the solver method
991
992 @return: the solver method currently be used.
993 @rtype: C{int}
994 """
995 return self.__solver_method,self.__preconditioner
996
997 def setSolverPackage(self,package=None):
998 """
999 sets a new solver package
1000
1001 @param package: sets a new solver method.
1002 @type package: one of L{DEFAULT}, L{PASO} L{SCSL}, L{MKL}, L{UMFPACK}, L{TRILINOS}
1003 """
1004 if package==None: package=self.DEFAULT
1005 if not package==self.getSolverPackage():
1006 self.__solver_package=package
1007 self.__checkMatrixType()
1008 self.trace("New solver is %s"%self.getSolverMethodName())
1009
1010 def getSolverPackage(self):
1011 """
1012 returns the package of the solver
1013
1014 @return: the solver package currently being used.
1015 @rtype: C{int}
1016 """
1017 return self.__solver_package
1018
1019 def isUsingLumping(self):
1020 """
1021 checks if matrix lumping is used a solver method
1022
1023 @return: True is lumping is currently used a solver method.
1024 @rtype: C{bool}
1025 """
1026 return self.getSolverMethod()[0]==self.LUMPING
1027
1028 def setTolerance(self,tol=1.e-8):
1029 """
1030 resets the tolerance for the solver method to tol where for an appropriate norm M{|.|}
1031
1032 M{|L{getResidual}()|<tol*|L{getRightHandSide}()|}
1033
1034 defines the stopping criterion.
1035
1036 @param tol: new tolerance for the solver. If the tol is lower then the current tolerence
1037 the system will be resolved.
1038 @type tol: positive C{float}
1039 @raise ValueError: if tolerance is not positive.
1040 """
1041 if not tol>0:
1042 raise ValueError,"Tolerance as to be positive"
1043 if tol<self.getTolerance(): self.__invalidateSolution()
1044 self.trace("New tolerance %e"%tol)
1045 self.__tolerance=tol
1046 return
1047
1048 def getTolerance(self):
1049 """
1050 returns the tolerance set for the solution
1051
1052 @return: tolerance currently used.
1053 @rtype: C{float}
1054 """
1055 return self.__tolerance
1056
1057 # =============================================================================
1058 # symmetry flag:
1059 # =============================================================================
1060 def isSymmetric(self):
1061 """
1062 checks if symmetry is indicated.
1063
1064 @return: True is a symmetric PDE is indicated, otherwise False is returned
1065 @rtype: C{bool}
1066 """
1067 return self.__sym
1068
1069 def setSymmetryOn(self):
1070 """
1071 sets the symmetry flag.
1072 """
1073 if not self.isSymmetric():
1074 self.trace("PDE is set to be symmetric")
1075 self.__sym=True
1076 self.__checkMatrixType()
1077
1078 def setSymmetryOff(self):
1079 """
1080 removes the symmetry flag.
1081 """
1082 if self.isSymmetric():
1083 self.trace("PDE is set to be unsymmetric")
1084 self.__sym=False
1085 self.__checkMatrixType()
1086
1087 def setSymmetryTo(self,flag=False):
1088 """
1089 sets the symmetry flag to flag
1090
1091 @param flag: If flag, the symmetry flag is set otherwise the symmetry flag is released.
1092 @type flag: C{bool}
1093 """
1094 if flag:
1095 self.setSymmetryOn()
1096 else:
1097 self.setSymmetryOff()
1098
1099 # =============================================================================
1100 # function space handling for the equation as well as the solution
1101 # =============================================================================
1102 def setReducedOrderOn(self):
1103 """
1104 switches on reduced order for solution and equation representation
1105
1106 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1107 """
1108 self.setReducedOrderForSolutionOn()
1109 self.setReducedOrderForEquationOn()
1110
1111 def setReducedOrderOff(self):
1112 """
1113 switches off reduced order for solution and equation representation
1114
1115 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1116 """
1117 self.setReducedOrderForSolutionOff()
1118 self.setReducedOrderForEquationOff()
1119
1120 def setReducedOrderTo(self,flag=False):
1121 """
1122 sets order reduction for both solution and equation representation according to flag.
1123 @param flag: if flag is True, the order reduction is switched on for both solution and equation representation, otherwise or
1124 if flag is not present order reduction is switched off
1125 @type flag: C{bool}
1126 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1127 """
1128 self.setReducedOrderForSolutionTo(flag)
1129 self.setReducedOrderForEquationTo(flag)
1130
1131
1132 def setReducedOrderForSolutionOn(self):
1133 """
1134 switches on reduced order for solution representation
1135
1136 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1137 """
1138 if not self.__reduce_solution_order:
1139 if self.__altered_coefficients:
1140 raise RuntimeError,"order cannot be altered after coefficients have been defined."
1141 self.trace("Reduced order is used to solution representation.")
1142 self.__reduce_solution_order=True
1143 self.__resetSystem()
1144
1145 def setReducedOrderForSolutionOff(self):
1146 """
1147 switches off reduced order for solution representation
1148
1149 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1150 """
1151 if self.__reduce_solution_order:
1152 if self.__altered_coefficients:
1153 raise RuntimeError,"order cannot be altered after coefficients have been defined."
1154 self.trace("Full order is used to interpolate solution.")
1155 self.__reduce_solution_order=False
1156 self.__resetSystem()
1157
1158 def setReducedOrderForSolutionTo(self,flag=False):
1159 """
1160 sets order for test functions according to flag
1161
1162 @param flag: if flag is True, the order reduction is switched on for solution representation, otherwise or
1163 if flag is not present order reduction is switched off
1164 @type flag: C{bool}
1165 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1166 """
1167 if flag:
1168 self.setReducedOrderForSolutionOn()
1169 else:
1170 self.setReducedOrderForSolutionOff()
1171
1172 def setReducedOrderForEquationOn(self):
1173 """
1174 switches on reduced order for equation representation
1175
1176 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1177 """
1178 if not self.__reduce_equation_order:
1179 if self.__altered_coefficients:
1180 raise RuntimeError,"order cannot be altered after coefficients have been defined."
1181 self.trace("Reduced order is used for test functions.")
1182 self.__reduce_equation_order=True
1183 self.__resetSystem()
1184
1185 def setReducedOrderForEquationOff(self):
1186 """
1187 switches off reduced order for equation representation
1188
1189 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1190 """
1191 if self.__reduce_equation_order:
1192 if self.__altered_coefficients:
1193 raise RuntimeError,"order cannot be altered after coefficients have been defined."
1194 self.trace("Full order is used for test functions.")
1195 self.__reduce_equation_order=False
1196 self.__resetSystem()
1197
1198 def setReducedOrderForEquationTo(self,flag=False):
1199 """
1200 sets order for test functions according to flag
1201
1202 @param flag: if flag is True, the order reduction is switched on for equation representation, otherwise or
1203 if flag is not present order reduction is switched off
1204 @type flag: C{bool}
1205 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
1206 """
1207 if flag:
1208 self.setReducedOrderForEquationOn()
1209 else:
1210 self.setReducedOrderForEquationOff()
1211
1212 # =============================================================================
1213 # private method:
1214 # =============================================================================
1215 def __checkMatrixType(self):
1216 """
1217 reassess the matrix type and, if a new matrix is needed, resets the system.
1218 """
1219 new_matrix_type=self.getDomain().getSystemMatrixTypeId(self.getSolverMethod()[0],self.getSolverPackage(),self.isSymmetric())
1220 if not new_matrix_type==self.__matrix_type:
1221 self.trace("Matrix type is now %d."%new_matrix_type)
1222 self.__matrix_type=new_matrix_type
1223 self.__resetSystem()
1224 #
1225 # rebuild switches :
1226 #
1227 def __invalidateSolution(self):
1228 """
1229 indicates the PDE has to be resolved if the solution is requested
1230 """
1231 if self.__solution_isValid: self.trace("PDE has to be resolved.")
1232 self.__solution_isValid=False
1233
1234 def __invalidateOperator(self):
1235 """
1236 indicates the operator has to be rebuilt next time it is used
1237 """
1238 if self.__operator_is_Valid: self.trace("Operator has to be rebuilt.")
1239 self.__invalidateSolution()
1240 self.__operator_is_Valid=False
1241
1242 def __invalidateRightHandSide(self):
1243 """
1244 indicates the right hand side has to be rebuild next time it is used
1245 """
1246 if self.__righthandside_isValid: self.trace("Right hand side has to be rebuilt.")
1247 self.__invalidateSolution()
1248 self.__righthandside_isValid=False
1249
1250 def __invalidateSystem(self):
1251 """
1252 annonced that everthing has to be rebuild:
1253 """
1254 if self.__righthandside_isValid: self.trace("System has to be rebuilt.")
1255 self.__invalidateSolution()
1256 self.__invalidateOperator()
1257 self.__invalidateRightHandSide()
1258
1259 def __resetSystem(self):
1260 """
1261 annonced that everthing has to be rebuild:
1262 """
1263 self.trace("New System is built from scratch.")
1264 self.__operator=escript.Operator()
1265 self.__operator_is_Valid=False
1266 self.__righthandside=escript.Data()
1267 self.__righthandside_isValid=False
1268 self.__solution=escript.Data()
1269 self.__solution_isValid=False
1270 #
1271 # system initialization:
1272 #
1273 def __getNewOperator(self):
1274 """
1275 returns an instance of a new operator
1276 """
1277 self.trace("New operator is allocated.")
1278 return self.getDomain().newOperator( \
1279 self.getNumEquations(), \
1280 self.getFunctionSpaceForEquation(), \
1281 self.getNumSolutions(), \
1282 self.getFunctionSpaceForSolution(), \
1283 self.__matrix_type)
1284
1285 def __getNewRightHandSide(self):
1286 """
1287 returns an instance of a new right hand side
1288 """
1289 self.trace("New right hand side is allocated.")
1290 if self.getNumEquations()>1:
1291 return escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)
1292 else:
1293 return escript.Data(0.,(),self.getFunctionSpaceForEquation(),True)
1294
1295 def __getNewSolution(self):
1296 """
1297 returns an instance of a new solution
1298 """
1299 self.trace("New solution is allocated.")
1300 if self.getNumSolutions()>1:
1301 return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)
1302 else:
1303 return escript.Data(0.,(),self.getFunctionSpaceForSolution(),True)
1304
1305 def __makeFreshSolution(self):
1306 """
1307 makes sure that the solution is instantiated and returns it initialized by zeros
1308 """
1309 if self.__solution.isEmpty():
1310 self.__solution=self.__getNewSolution()
1311 else:
1312 self.__solution*=0
1313 self.trace("Solution is reset to zero.")
1314 return self.__solution
1315
1316 def __makeFreshRightHandSide(self):
1317 """
1318 makes sure that the right hand side is instantiated and returns it initialized by zeros
1319 """
1320 if self.__righthandside.isEmpty():
1321 self.__righthandside=self.__getNewRightHandSide()
1322 else:
1323 self.__righthandside.setToZero()
1324 self.trace("Right hand side is reset to zero.")
1325 return self.__righthandside
1326
1327 def __makeFreshOperator(self):
1328 """
1329 makes sure that the operator is instantiated and returns it initialized by zeros
1330 """
1331 if self.__operator.isEmpty():
1332 self.__operator=self.__getNewOperator()
1333 else:
1334 self.__operator.resetValues()
1335 self.trace("Operator reset to zero")
1336 return self.__operator
1337
1338 def __applyConstraint(self):
1339 """
1340 applies the constraints defined by q and r to the system
1341 """
1342 if not self.isUsingLumping():
1343 q=self.getCoefficientOfGeneralPDE("q")
1344 r=self.getCoefficientOfGeneralPDE("r")
1345 if not q.isEmpty() and not self.__operator.isEmpty():
1346 # q is the row and column mask to indicate where constraints are set:
1347 row_q=escript.Data(q,self.getFunctionSpaceForEquation())
1348 col_q=escript.Data(q,self.getFunctionSpaceForSolution())
1349 u=self.__getNewSolution()
1350 if r.isEmpty():
1351 r_s=self.__getNewSolution()
1352 else:
1353 r_s=escript.Data(r,self.getFunctionSpaceForSolution())
1354 u.copyWithMask(r_s,col_q)
1355 if not self.__righthandside.isEmpty():
1356 self.__righthandside-=self.__operator*u
1357 self.__righthandside=self.copyConstraint(self.__righthandside)
1358 self.__operator.nullifyRowsAndCols(row_q,col_q,1.)
1359 # =============================================================================
1360 # function giving access to coefficients of the general PDE:
1361 # =============================================================================
1362 def getCoefficientOfGeneralPDE(self,name):
1363 """
1364 return the value of the coefficient name of the general PDE.
1365
1366 @note: This method is called by the assembling routine it can be overwritten
1367 to map coefficients of a particular PDE to the general PDE.
1368 @param name: name of the coefficient requested.
1369 @type name: C{string}
1370 @return: the value of the coefficient name
1371 @rtype: L{Data<escript.Data>}
1372 @raise IllegalCoefficient: if name is not one of coefficients
1373 M{A}, M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact},
1374 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}.
1375 """
1376 if self.hasCoefficientOfGeneralPDE(name):
1377 return self.getCoefficient(name)
1378 else:
1379 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1380
1381 def hasCoefficientOfGeneralPDE(self,name):
1382 """
1383 checks if name is a the name of a coefficient of the general PDE.
1384
1385 @param name: name of the coefficient enquired.
1386 @type name: C{string}
1387 @return: True if name is the name of a coefficient of the general PDE. Otherwise False.
1388 @rtype: C{bool}
1389
1390 """
1391 return self.__COEFFICIENTS_OF_GENEARL_PDE.has_key(name)
1392
1393 def createCoefficientOfGeneralPDE(self,name):
1394 """
1395 returns a new instance of a coefficient for coefficient name of the general PDE
1396
1397 @param name: name of the coefficient requested.
1398 @type name: C{string}
1399 @return: a coefficient name initialized to 0.
1400 @rtype: L{Data<escript.Data>}
1401 @raise IllegalCoefficient: if name is not one of coefficients
1402 M{A}, M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact},
1403 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}.
1404 """
1405 if self.hasCoefficientOfGeneralPDE(name):
1406 return escript.Data(0,self.getShapeOfCoefficientOfGeneralPDE(name),self.getFunctionSpaceForCoefficientOfGeneralPDE(name))
1407 else:
1408 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1409
1410 def getFunctionSpaceForCoefficientOfGeneralPDE(self,name):
1411 """
1412 return the L{FunctionSpace<escript.FunctionSpace>} to be used for coefficient name of the general PDE
1413
1414 @param name: name of the coefficient enquired.
1415 @type name: C{string}
1416 @return: the function space to be used for coefficient name
1417 @rtype: L{FunctionSpace<escript.FunctionSpace>}
1418 @raise IllegalCoefficient: if name is not one of coefficients
1419 M{A}, M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact},
1420 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}.
1421 """
1422 if self.hasCoefficientOfGeneralPDE(name):
1423 return self.__COEFFICIENTS_OF_GENEARL_PDE[name].getFunctionSpace(self.getDomain())
1424 else:
1425 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1426
1427 def getShapeOfCoefficientOfGeneralPDE(self,name):
1428 """
1429 return the shape of the coefficient name of the general PDE
1430
1431 @param name: name of the coefficient enquired.
1432 @type name: C{string}
1433 @return: the shape of the coefficient name
1434 @rtype: C{tuple} of C{int}
1435 @raise IllegalCoefficient: if name is not one of coefficients
1436 M{A}, M{B}, M{C}, M{D}, M{X}, M{Y}, M{d}, M{y}, M{d_contact}, M{y_contact},
1437 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}.
1438 """
1439 if self.hasCoefficientOfGeneralPDE(name):
1440 return self.__COEFFICIENTS_OF_GENEARL_PDE[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())
1441 else:
1442 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1443
1444 # =============================================================================
1445 # functions giving access to coefficients of a particular PDE implementation:
1446 # =============================================================================
1447 def getCoefficient(self,name):
1448 """
1449 returns the value of the coefficient name
1450
1451 @param name: name of the coefficient requested.
1452 @type name: C{string}
1453 @return: the value of the coefficient name
1454 @rtype: L{Data<escript.Data>}
1455 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1456 """
1457 if self.hasCoefficient(name):
1458 return self.COEFFICIENTS[name].getValue()
1459 else:
1460 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1461
1462 def hasCoefficient(self,name):
1463 """
1464 return True if name is the name of a coefficient
1465
1466 @param name: name of the coefficient enquired.
1467 @type name: C{string}
1468 @return: True if name is the name of a coefficient of the general PDE. Otherwise False.
1469 @rtype: C{bool}
1470 """
1471 return self.COEFFICIENTS.has_key(name)
1472
1473 def createCoefficient(self, name):
1474 """
1475 create a L{Data<escript.Data>} object corresponding to coefficient name
1476
1477 @return: a coefficient name initialized to 0.
1478 @rtype: L{Data<escript.Data>}
1479 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1480 """
1481 if self.hasCoefficient(name):
1482 return escript.Data(0.,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))
1483 else:
1484 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1485
1486 def getFunctionSpaceForCoefficient(self,name):
1487 """
1488 return the L{FunctionSpace<escript.FunctionSpace>} to be used for coefficient name
1489
1490 @param name: name of the coefficient enquired.
1491 @type name: C{string}
1492 @return: the function space to be used for coefficient name
1493 @rtype: L{FunctionSpace<escript.FunctionSpace>}
1494 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1495 """
1496 if self.hasCoefficient(name):
1497 return self.COEFFICIENTS[name].getFunctionSpace(self.getDomain())
1498 else:
1499 raise ValueError,"unknown coefficient %s requested"%name
1500 def getShapeOfCoefficient(self,name):
1501 """
1502 return the shape of the coefficient name
1503
1504 @param name: name of the coefficient enquired.
1505 @type name: C{string}
1506 @return: the shape of the coefficient name
1507 @rtype: C{tuple} of C{int}
1508 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1509 """
1510 if self.hasCoefficient(name):
1511 return self.COEFFICIENTS[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())
1512 else:
1513 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1514
1515 def resetCoefficients(self):
1516 """
1517 resets all coefficients to there default values.
1518 """
1519 for i in self.COEFFICIENTS.iterkeys():
1520 self.COEFFICIENTS[i].resetValue()
1521
1522 def alteredCoefficient(self,name):
1523 """
1524 announce that coefficient name has been changed
1525
1526 @param name: name of the coefficient enquired.
1527 @type name: C{string}
1528 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1529 @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.
1530 """
1531 if self.hasCoefficient(name):
1532 self.trace("Coefficient %s has been altered."%name)
1533 if not ((name=="q" or name=="r") and self.isUsingLumping()):
1534 if self.COEFFICIENTS[name].isAlteringOperator(): self.__invalidateOperator()
1535 if self.COEFFICIENTS[name].isAlteringRightHandSide(): self.__invalidateRightHandSide()
1536 else:
1537 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1538
1539 def copyConstraint(self,u):
1540 """
1541 copies the constraint into u and returns u.
1542
1543 @param u: a function of rank 0 is a single PDE is solved and of shape (numSolution,) for a system of PDEs
1544 @type u: L{Data<escript.Data>}
1545 @return: the input u modified by the constraints.
1546 @rtype: L{Data<escript.Data>}
1547 @warning: u is altered if it has the appropriate L{FunctionSpace<escript.FunctionSpace>}
1548 """
1549 q=self.getCoefficientOfGeneralPDE("q")
1550 r=self.getCoefficientOfGeneralPDE("r")
1551 if not q.isEmpty():
1552 if u.isEmpty(): u=escript.Data(0.,q.getShape(),q.getFunctionSpace())
1553 if r.isEmpty():
1554 r=escript.Data(0,u.getShape(),u.getFunctionSpace())
1555 else:
1556 r=escript.Data(r,u.getFunctionSpace())
1557 u.copyWithMask(r,escript.Data(q,u.getFunctionSpace()))
1558 return u
1559
1560 def setValue(self,**coefficients):
1561 """
1562 sets new values to coefficients
1563
1564 @param coefficients: new values assigned to coefficients
1565 @keyword A: value for coefficient A.
1566 @type A: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1567 @keyword A_reduced: value for coefficient A_reduced.
1568 @type A_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1569 @keyword B: value for coefficient B
1570 @type B: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1571 @keyword B_reduced: value for coefficient B_reduced
1572 @type B_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1573 @keyword C: value for coefficient C
1574 @type C: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1575 @keyword C_reduced: value for coefficient C_reduced
1576 @type C_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1577 @keyword D: value for coefficient D
1578 @type D: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1579 @keyword D_reduced: value for coefficient D_reduced
1580 @type D_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1581 @keyword X: value for coefficient X
1582 @type X: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1583 @keyword X_reduced: value for coefficient X_reduced
1584 @type X_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1585 @keyword Y: value for coefficient Y
1586 @type Y: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1587 @keyword Y_reduced: value for coefficient Y_reduced
1588 @type Y_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.Function>}.
1589 @keyword d: value for coefficient d
1590 @type d: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1591 @keyword d_reduced: value for coefficient d_reduced
1592 @type d_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
1593 @keyword y: value for coefficient y
1594 @type y: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1595 @keyword d_contact: value for coefficient d_contact
1596 @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>}.
1597 @keyword d_contact_reduced: value for coefficient d_contact_reduced
1598 @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>}.
1599 @keyword y_contact: value for coefficient y_contact
1600 @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>}.
1601 @keyword y_contact_reduced: value for coefficient y_contact_reduced
1602 @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>}.
1603 @keyword r: values prescribed to the solution at the locations of constraints
1604 @type r: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1605 depending of reduced order is used for the solution.
1606 @keyword q: mask for location of constraints
1607 @type q: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1608 depending of reduced order is used for the representation of the equation.
1609 @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1610 """
1611 # check if the coefficients are legal:
1612 for i in coefficients.iterkeys():
1613 if not self.hasCoefficient(i):
1614 raise IllegalCoefficient,"Attempt to set unknown coefficient %s"%i
1615 # if the number of unknowns or equations is still unknown we try to estimate them:
1616 if self.__numEquations==None or self.__numSolutions==None:
1617 for i,d in coefficients.iteritems():
1618 if hasattr(d,"shape"):
1619 s=d.shape
1620 elif hasattr(d,"getShape"):
1621 s=d.getShape()
1622 else:
1623 s=numarray.array(d).shape
1624 if s!=None:
1625 # get number of equations and number of unknowns:
1626 res=self.COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(self.getDomain(),s)
1627 if res==None:
1628 raise IllegalCoefficientValue,"Illegal shape %s of coefficient %s"%(s,i)
1629 else:
1630 if self.__numEquations==None: self.__numEquations=res[0]
1631 if self.__numSolutions==None: self.__numSolutions=res[1]
1632 if self.__numEquations==None: raise UndefinedPDEError,"unidententified number of equations"
1633 if self.__numSolutions==None: raise UndefinedPDEError,"unidententified number of solutions"
1634 # now we check the shape of the coefficient if numEquations and numSolutions are set:
1635 for i,d in coefficients.iteritems():
1636 try:
1637 self.COEFFICIENTS[i].setValue(self.getDomain(),
1638 self.getNumEquations(),self.getNumSolutions(),
1639 self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
1640 self.alteredCoefficient(i)
1641 except IllegalCoefficientFunctionSpace,m:
1642 # if the function space is wrong then we try the reduced version:
1643 i_red=i+"_reduced"
1644 if (not i_red in coefficients.keys()) and i_red in self.COEFFICIENTS.keys():
1645 try:
1646 self.COEFFICIENTS[i_red].setValue(self.getDomain(),
1647 self.getNumEquations(),self.getNumSolutions(),
1648 self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
1649 self.alteredCoefficient(i_red)
1650 except IllegalCoefficientValue,m:
1651 raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
1652 except IllegalCoefficientFunctionSpace,m:
1653 raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m))
1654 else:
1655 raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m))
1656 except IllegalCoefficientValue,m:
1657 raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
1658 self.__altered_coefficients=True
1659 # check if the systrem is inhomogeneous:
1660 if len(coefficients)>0 and not self.isUsingLumping():
1661 q=self.getCoefficientOfGeneralPDE("q")
1662 r=self.getCoefficientOfGeneralPDE("r")
1663 homogeneous_constraint=True
1664 if not q.isEmpty() and not r.isEmpty():
1665 if util.Lsup(q*r)>0.:
1666 self.trace("Inhomogeneous constraint detected.")
1667 self.__invalidateSystem()
1668
1669 def getSystem(self):
1670 """
1671 return the operator and right hand side of the PDE
1672
1673 @return: the discrete version of the PDE
1674 @rtype: C{tuple} of L{Operator,<escript.Operator>} and L{Data<escript.Data>}.
1675 """
1676 if not self.__operator_is_Valid or not self.__righthandside_isValid:
1677 if self.isUsingLumping():
1678 if not self.__operator_is_Valid:
1679 if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():
1680 raise TypeError,"Lumped matrix requires same order for equations and unknowns"
1681 if not self.getCoefficientOfGeneralPDE("A").isEmpty():
1682 raise ValueError,"coefficient A in lumped matrix may not be present."
1683 if not self.getCoefficientOfGeneralPDE("B").isEmpty():
1684 raise ValueError,"coefficient B in lumped matrix may not be present."
1685 if not self.getCoefficientOfGeneralPDE("C").isEmpty():
1686 raise ValueError,"coefficient C in lumped matrix may not be present."
1687 if not self.getCoefficientOfGeneralPDE("d_contact").isEmpty():
1688 raise ValueError,"coefficient d_contact in lumped matrix may not be present."
1689 if not self.getCoefficientOfGeneralPDE("A_reduced").isEmpty():
1690 raise ValueError,"coefficient A_reduced in lumped matrix may not be present."
1691 if not self.getCoefficientOfGeneralPDE("B_reduced").isEmpty():
1692 raise ValueError,"coefficient B_reduced in lumped matrix may not be present."
1693 if not self.getCoefficientOfGeneralPDE("C_reduced").isEmpty():
1694 raise ValueError,"coefficient C_reduced in lumped matrix may not be present."
1695 if not self.getCoefficientOfGeneralPDE("d_contact_reduced").isEmpty():
1696 raise ValueError,"coefficient d_contact_reduced in lumped matrix may not be present."
1697 D=self.getCoefficientOfGeneralPDE("D")
1698 d=self.getCoefficientOfGeneralPDE("d")
1699 D_reduced=self.getCoefficientOfGeneralPDE("D_reduced")
1700 d_reduced=self.getCoefficientOfGeneralPDE("d_reduced")
1701 if not D.isEmpty():
1702 if self.getNumSolutions()>1:
1703 D_times_e=util.matrix_mult(D,numarray.ones((self.getNumSolutions(),)))
1704 else:
1705 D_times_e=D
1706 else:
1707 D_times_e=escript.Data()
1708 if not d.isEmpty():
1709 if self.getNumSolutions()>1:
1710 d_times_e=util.matrix_mult(d,numarray.ones((self.getNumSolutions(),)))
1711 else:
1712 d_times_e=d
1713 else:
1714 d_times_e=escript.Data()
1715
1716 if not D_reduced.isEmpty():
1717 if self.getNumSolutions()>1:
1718 D_reduced_times_e=util.matrix_mult(D_reduced,numarray.ones((self.getNumSolutions(),)))
1719 else:
1720 D_reduced_times_e=D_reduced
1721 else:
1722 D_reduced_times_e=escript.Data()
1723 if not d_reduced.isEmpty():
1724 if self.getNumSolutions()>1:
1725 d_reduced_times_e=util.matrix_mult(d_reduced,numarray.ones((self.getNumSolutions(),)))
1726 else:
1727 d_reduced_times_e=d_reduced
1728 else:
1729 d_reduced_times_e=escript.Data()
1730
1731 self.__operator=self.__getNewRightHandSide()
1732 if hasattr(self.getDomain(), "addPDEToLumpedSystem") :
1733 self.getDomain().addPDEToLumpedSystem(self.__operator, D_times_e, d_times_e)
1734 self.getDomain().addPDEToLumpedSystem(self.__operator, D_reduced_times_e, d_reduced_times_e)
1735 else:
1736 self.getDomain().addPDEToRHS(self.__operator, \
1737 escript.Data(), \
1738 D_times_e, \
1739 d_times_e,\
1740 escript.Data())
1741 self.getDomain().addPDEToRHS(self.__operator, \
1742 escript.Data(), \
1743 D_reduced_times_e, \
1744 d_reduced_times_e,\
1745 escript.Data())
1746 self.__operator=1./self.__operator
1747 self.trace("New lumped operator has been built.")
1748 self.__operator_is_Valid=True
1749 if not self.__righthandside_isValid:
1750 self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \
1751 self.getCoefficientOfGeneralPDE("X"), \
1752 self.getCoefficientOfGeneralPDE("Y"),\
1753 self.getCoefficientOfGeneralPDE("y"),\
1754 self.getCoefficientOfGeneralPDE("y_contact"))
1755 self.getDomain().addPDEToRHS(self.__righthandside, \
1756 self.getCoefficientOfGeneralPDE("X_reduced"), \
1757 self.getCoefficientOfGeneralPDE("Y_reduced"),\
1758 self.getCoefficientOfGeneralPDE("y_reduced"),\
1759 self.getCoefficientOfGeneralPDE("y_contact_reduced"))
1760 self.trace("New right hand side as been built.")
1761 self.__righthandside_isValid=True
1762 else:
1763 if not self.__operator_is_Valid and not self.__righthandside_isValid:
1764 self.getDomain().addPDEToSystem(self.__makeFreshOperator(),self.__makeFreshRightHandSide(), \
1765 self.getCoefficientOfGeneralPDE("A"), \
1766 self.getCoefficientOfGeneralPDE("B"), \
1767 self.getCoefficientOfGeneralPDE("C"), \
1768 self.getCoefficientOfGeneralPDE("D"), \
1769 self.getCoefficientOfGeneralPDE("X"), \
1770 self.getCoefficientOfGeneralPDE("Y"), \
1771 self.getCoefficientOfGeneralPDE("d"), \
1772 self.getCoefficientOfGeneralPDE("y"), \
1773 self.getCoefficientOfGeneralPDE("d_contact"), \
1774 self.getCoefficientOfGeneralPDE("y_contact"))
1775 self.getDomain().addPDEToSystem(self.__operator,self.__righthandside, \
1776 self.getCoefficientOfGeneralPDE("A_reduced"), \
1777 self.getCoefficientOfGeneralPDE("B_reduced"), \
1778 self.getCoefficientOfGeneralPDE("C_reduced"), \
1779 self.getCoefficientOfGeneralPDE("D_reduced"), \
1780 self.getCoefficientOfGeneralPDE("X_reduced"), \
1781 self.getCoefficientOfGeneralPDE("Y_reduced"), \
1782 self.getCoefficientOfGeneralPDE("d_reduced"), \
1783 self.getCoefficientOfGeneralPDE("y_reduced"), \
1784 self.getCoefficientOfGeneralPDE("d_contact_reduced"), \
1785 self.getCoefficientOfGeneralPDE("y_contact_reduced"))
1786 self.__applyConstraint()
1787 self.__righthandside=self.copyConstraint(self.__righthandside)
1788 self.trace("New system has been built.")
1789 self.__operator_is_Valid=True
1790 self.__righthandside_isValid=True
1791 elif not self.__righthandside_isValid:
1792 self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \
1793 self.getCoefficientOfGeneralPDE("X"), \
1794 self.getCoefficientOfGeneralPDE("Y"),\
1795 self.getCoefficientOfGeneralPDE("y"),\
1796 self.getCoefficientOfGeneralPDE("y_contact"))
1797 self.getDomain().addPDEToRHS(self.__righthandside, \
1798 self.getCoefficientOfGeneralPDE("X_reduced"), \
1799 self.getCoefficientOfGeneralPDE("Y_reduced"),\
1800 self.getCoefficientOfGeneralPDE("y_reduced"),\
1801 self.getCoefficientOfGeneralPDE("y_contact_reduced"))
1802 self.__righthandside=self.copyConstraint(self.__righthandside)
1803 self.trace("New right hand side has been built.")
1804 self.__righthandside_isValid=True
1805 elif not self.__operator_is_Valid:
1806 self.getDomain().addPDEToSystem(self.__makeFreshOperator(),escript.Data(), \
1807 self.getCoefficientOfGeneralPDE("A"), \
1808 self.getCoefficientOfGeneralPDE("B"), \
1809 self.getCoefficientOfGeneralPDE("C"), \
1810 self.getCoefficientOfGeneralPDE("D"), \
1811 escript.Data(), \
1812 escript.Data(), \
1813 self.getCoefficientOfGeneralPDE("d"), \
1814 escript.Data(),\
1815 self.getCoefficientOfGeneralPDE("d_contact"), \
1816 escript.Data())
1817 self.getDomain().addPDEToSystem(self.__operator,escript.Data(), \
1818 self.getCoefficientOfGeneralPDE("A_reduced"), \
1819 self.getCoefficientOfGeneralPDE("B_reduced"), \
1820 self.getCoefficientOfGeneralPDE("C_reduced"), \
1821 self.getCoefficientOfGeneralPDE("D_reduced"), \
1822 escript.Data(), \
1823 escript.Data(), \
1824 self.getCoefficientOfGeneralPDE("d_reduced"), \
1825 escript.Data(),\
1826 self.getCoefficientOfGeneralPDE("d_contact_reduced"), \
1827 escript.Data())
1828 self.__applyConstraint()
1829 self.trace("New operator has been built.")
1830 self.__operator_is_Valid=True
1831 return (self.__operator,self.__righthandside)
1832
1833
1834 class Poisson(LinearPDE):
1835 """
1836 Class to define a Poisson equation problem, which is genear L{LinearPDE} of the form
1837
1838 M{-grad(grad(u)[j])[j] = f}
1839
1840 with natural boundary conditons
1841
1842 M{n[j]*grad(u)[j] = 0 }
1843
1844 and constraints:
1845
1846 M{u=0} where M{q>0}
1847
1848 """
1849
1850 def __init__(self,domain,debug=False):
1851 """
1852 initializes a new Poisson equation
1853
1854 @param domain: domain of the PDE
1855 @type domain: L{Domain<escript.Domain>}
1856 @param debug: if True debug informations are printed.
1857
1858 """
1859 super(Poisson, self).__init__(domain,1,1,debug)
1860 self.COEFFICIENTS={"f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1861 "f_reduced": PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1862 "q": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
1863 self.setSymmetryOn()
1864
1865 def setValue(self,**coefficients):
1866 """
1867 sets new values to coefficients
1868
1869 @param coefficients: new values assigned to coefficients
1870 @keyword f: value for right hand side M{f}
1871 @type f: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1872 @keyword q: mask for location of constraints
1873 @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>}
1874 depending of reduced order is used for the representation of the equation.
1875 @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1876 """
1877 super(Poisson, self).setValue(**coefficients)
1878
1879 def getCoefficientOfGeneralPDE(self,name):
1880 """
1881 return the value of the coefficient name of the general PDE
1882 @param name: name of the coefficient requested.
1883 @type name: C{string}
1884 @return: the value of the coefficient name
1885 @rtype: L{Data<escript.Data>}
1886 @raise IllegalCoefficient: if name is not one of coefficients
1887 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}.
1888 @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
1889 """
1890 if name == "A" :
1891 return escript.Data(util.kronecker(self.getDim()),escript.Function(self.getDomain()))
1892 elif name == "B" :
1893 return escript.Data()
1894 elif name == "C" :
1895 return escript.Data()
1896 elif name == "D" :
1897 return escript.Data()
1898 elif name == "X" :
1899 return escript.Data()
1900 elif name == "Y" :
1901 return self.getCoefficient("f")
1902 elif name == "d" :
1903 return escript.Data()
1904 elif name == "y" :
1905 return escript.Data()
1906 elif name == "d_contact" :
1907 return escript.Data()
1908 elif name == "y_contact" :
1909 return escript.Data()
1910 elif name == "A_reduced" :
1911 return escript.Data()
1912 elif name == "B_reduced" :
1913 return escript.Data()
1914 elif name == "C_reduced" :
1915 return escript.Data()
1916 elif name == "D_reduced" :
1917 return escript.Data()
1918 elif name == "X_reduced" :
1919 return escript.Data()
1920 elif name == "Y_reduced" :
1921 return self.getCoefficient("f_reduced")
1922 elif name == "d_reduced" :
1923 return escript.Data()
1924 elif name == "y_reduced" :
1925 return escript.Data()
1926 elif name == "d_contact_reduced" :
1927 return escript.Data()
1928 elif name == "y_contact_reduced" :
1929 return escript.Data()
1930 elif name == "r" :
1931 return escript.Data()
1932 elif name == "q" :
1933 return self.getCoefficient("q")
1934 else:
1935 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1936
1937 class Helmholtz(LinearPDE):
1938 """
1939 Class to define a Helmhotz equation problem, which is genear L{LinearPDE} of the form
1940
1941 M{S{omega}*u - grad(k*grad(u)[j])[j] = f}
1942
1943 with natural boundary conditons
1944
1945 M{k*n[j]*grad(u)[j] = g- S{alpha}u }
1946
1947 and constraints:
1948
1949 M{u=r} where M{q>0}
1950
1951 """
1952
1953 def __init__(self,domain,debug=False):
1954 """
1955 initializes a new Poisson equation
1956
1957 @param domain: domain of the PDE
1958 @type domain: L{Domain<escript.Domain>}
1959 @param debug: if True debug informations are printed.
1960
1961 """
1962 super(Helmholtz, self).__init__(domain,1,1,debug)
1963 self.COEFFICIENTS={"omega": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
1964 "k": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
1965 "f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1966 "f_reduced": PDECoefficient(PDECoefficient.INTERIOR_REDUCED,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1967 "alpha": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.OPERATOR),
1968 "g": PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1969 "g_reduced": PDECoefficient(PDECoefficient.BOUNDARY_REDUCED,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1970 "r": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH),
1971 "q": PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
1972 self.setSymmetryOn()
1973
1974 def setValue(self,**coefficients):
1975 """
1976 sets new values to coefficients
1977
1978 @param coefficients: new values assigned to coefficients
1979 @keyword omega: value for coefficient M{S{omega}}
1980 @type omega: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1981 @keyword k: value for coefficeint M{k}
1982 @type k: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1983 @keyword f: value for right hand side M{f}
1984 @type f: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1985 @keyword alpha: value for right hand side M{S{alpha}}
1986 @type alpha: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1987 @keyword g: value for right hand side M{g}
1988 @type g: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1989 @keyword r: prescribed values M{r} for the solution in constraints.
1990 @type r: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1991 depending of reduced order is used for the representation of the equation.
1992 @keyword q: mask for location of constraints
1993 @type q: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1994 depending of reduced order is used for the representation of the equation.
1995 @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1996 """
1997 super(Helmholtz, self).setValue(**coefficients)
1998
1999 def getCoefficientOfGeneralPDE(self,name):
2000 """
2001 return the value of the coefficient name of the general PDE
2002
2003 @param name: name of the coefficient requested.
2004 @type name: C{string}
2005 @return: the value of the coefficient name
2006 @rtype: L{Data<escript.Data>}
2007 @raise IllegalCoefficient: if name is not one of coefficients
2008 "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}.
2009 @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
2010 """
2011 if name == "A" :
2012 return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))*self.getCoefficient("k")
2013 elif name == "B" :
2014 return escript.Data()
2015 elif name == "C" :
2016 return escript.Data()
2017 elif name == "D" :
2018 return self.getCoefficient("omega")
2019 elif name == "X" :
2020 return escript.Data()
2021 elif name == "Y" :
2022 return self.getCoefficient("f")
2023 elif name == "d" :
2024 return self.getCoefficient("alpha")
2025 elif name == "y" :
2026 return self.getCoefficient("g")
2027 elif name == "d_contact" :
2028 return escript.Data()
2029 elif name == "y_contact" :
2030 return escript.Data()
2031 elif name == "A_reduced" :
2032 return escript.Data()
2033 elif name == "B_reduced" :
2034 return escript.Data()
2035 elif name == "C_reduced" :
2036 return escript.Data()
2037 elif name == "D_reduced" :
2038 return escript.Data()
2039 elif name == "X_reduced" :
2040 return escript.Data()
2041 elif name == "Y_reduced" :
2042 return self.getCoefficient("f_reduced")
2043 elif name == "d_reduced" :
2044 return escript.Data()
2045 elif name == "y_reduced" :
2046 return self.getCoefficient("g_reduced")
2047 elif name == "d_contact_reduced" :
2048 return escript.Data()
2049 elif name == "y_contact_reduced" :
2050 return escript.Data()
2051 elif name == "r" :
2052 return self.getCoefficient("r")
2053 elif name == "q" :
2054 return self.getCoefficient("q")
2055 else:
2056 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
2057
2058 class LameEquation(LinearPDE):
2059 """
2060 Class to define a Lame equation problem:
2061
2062 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] }
2063
2064 with natural boundary conditons:
2065
2066 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] }
2067
2068 and constraints:
2069
2070 M{u[i]=r[i]} where M{q[i]>0}
2071
2072 """
2073
2074 def __init__(self,domain,debug=False):
2075 super(LameEquation, self).__init__(domain,\
2076 domain.getDim(),domain.getDim(),debug)
2077 self.COEFFICIENTS={ "lame_lambda" : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),
2078 "lame_mu" : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),
2079 "F" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
2080 "sigma" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.BY_EQUATION,PDECoefficient.BY_DIM),PDECoefficient.RIGHTHANDSIDE),
2081 "f" : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.BY_EQUATION,),PDECoefficient.RIGHTHANDSIDE),
2082 "r" : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH),
2083 "q" : PDECoefficient(PDECoefficient.SOLUTION,(PDECoefficient.BY_EQUATION,),PDECoefficient.BOTH)}
2084 self.setSymmetryOn()
2085
2086 def setValues(self,**coefficients):
2087 """
2088 sets new values to coefficients
2089
2090 @param coefficients: new values assigned to coefficients
2091 @keyword lame_mu: value for coefficient M{S{mu}}
2092 @type lame_mu: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2093 @keyword lame_lambda: value for coefficient M{S{lambda}}
2094 @type lame_lambda: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2095 @keyword F: value for internal force M{F}
2096 @type F: any type that can be casted to L{Vector<escript.Vector>} object on L{Function<escript.Function>}
2097 @keyword sigma: value for initial stress M{S{sigma}}
2098 @type sigma: any type that can be casted to L{Tensor<escript.Tensor>} object on L{Function<escript.Function>}
2099 @keyword f: value for extrenal force M{f}
2100 @type f: any type that can be casted to L{Vector<escript.Vector>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}
2101 @keyword r: prescribed values M{r} for the solution in constraints.
2102 @type r: any type that can be casted to L{Vector<escript.Vector>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2103 depending of reduced order is used for the representation of the equation.
2104 @keyword q: mask for location of constraints
2105 @type q: any type that can be casted to L{Vector<escript.Vector>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2106 depending of reduced order is used for the representation of the equation.
2107 @raise IllegalCoefficient: if an unknown coefficient keyword is used.
2108 """
2109 super(LameEquation, self).setValues(**coefficients)
2110
2111 def getCoefficientOfGeneralPDE(self,name):
2112 """
2113 return the value of the coefficient name of the general PDE
2114
2115 @param name: name of the coefficient requested.
2116 @type name: C{string}
2117 @return: the value of the coefficient name
2118 @rtype: L{Data<escript.Data>}
2119 @raise IllegalCoefficient: if name is not one of coefficients
2120 "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}.
2121 @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
2122 """
2123 if name == "A" :
2124 out =self.createCoefficientOfGeneralPDE("A")
2125 for i in range(self.getDim()):
2126 for j in range(self.getDim()):
2127 out[i,i,j,j] += self.getCoefficient("lame_lambda")
2128 out[i,j,j,i] += self.getCoefficient("lame_mu")
2129 out[i,j,i,j] += self.getCoefficient("lame_mu")
2130 return out
2131 elif name == "B" :
2132 return escript.Data()
2133 elif name == "C" :
2134 return escript.Data()
2135 elif name == "D" :
2136 return escript.Data()
2137 elif name == "X" :
2138 return self.getCoefficient("sigma")
2139 elif name == "Y" :
2140 return self.getCoefficient("F")
2141 elif name == "d" :
2142 return escript.Data()
2143 elif name == "y" :
2144 return self.getCoefficient("f")
2145 elif name == "d_contact" :
2146 return escript.Data()
2147 elif name == "y_contact" :
2148 return escript.Data()
2149 elif name == "A_reduced" :
2150 return escript.Data()
2151 elif name == "B_reduced" :
2152 return escript.Data()
2153 elif name == "C_reduced" :
2154 return escript.Data()
2155 elif name == "D_reduced" :
2156 return escript.Data()
2157 elif name == "X_reduced" :
2158 return escript.Data()
2159 elif name == "Y_reduced" :
2160 return escript.Data()
2161 elif name == "d_reduced" :
2162 return escript.Data()
2163 elif name == "y_reduced" :
2164 return escript.Data()
2165 elif name == "d_contact_reduced" :
2166 return escript.Data()
2167 elif name == "y_contact_reduced" :
2168 return escript.Data()
2169 elif name == "r" :
2170 return self.getCoefficient("r")
2171 elif name == "q" :
2172 return self.getCoefficient("q")
2173 else:
2174 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
2175

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26