/[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 969 - (show annotations)
Tue Feb 13 23:02:23 2007 UTC (12 years, 10 months ago) by ksteube
File MIME type: text/x-python
File size: 106141 byte(s)
Parallelization using MPI for solution of implicit problems.

Parallelization for explicit problems has already been accomplished in
the main SVN branch.

This is incomplete and is not ready for use.


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26