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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1137 - (show annotations)
Thu May 10 08:11:31 2007 UTC (12 years, 7 months ago) by gross
File MIME type: text/x-python
File size: 102616 byte(s)
This version passes the tests on windows except for 

   * vtk
   * netCDF

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26