/[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 971 - (show annotations)
Wed Feb 14 04:40:49 2007 UTC (12 years, 8 months ago) by ksteube
File MIME type: text/x-python
File size: 105977 byte(s)
Had to undo commit to new MPI branch. The changes went into the original and
not the branch. The files committed here are exactly the same as revision 969.


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26