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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 614 - (show annotations)
Wed Mar 22 01:37:07 2006 UTC (13 years, 6 months ago) by elspeth
Original Path: trunk/escript/py_src/linearPDEs.py
File MIME type: text/x-python
File size: 105731 byte(s)
Corrected spelling of 'license' in url so that the link actually points to the license.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26