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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26