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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26