/[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 150 - (show annotations)
Thu Sep 15 03:44:45 2005 UTC (14 years, 2 months ago) by jgs
Original Path: trunk/esys2/escript/py_src/linearPDEs.py
File MIME type: text/x-python
File size: 102221 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-09-15

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26