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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26