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

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26