/[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 147 - (show annotations)
Fri Aug 12 01:45:47 2005 UTC (14 years, 4 months ago) by jgs
Original Path: trunk/esys2/escript/py_src/linearPDEs.py
File MIME type: text/x-python
File size: 55333 byte(s)
erge of development branch dev-02 back to main trunk on 2005-08-12

1 # $Id$
2
3 ## @file linearPDEs.py
4
5 """
6 Functions and classes for linear PDEs
7 """
8
9 import escript
10 import util
11 import numarray
12
13
14 def _CompTuple2(t1,t2):
15 """
16 Compare two tuples
17
18 \param t1 The first tuple
19 \param t2 The second tuple
20 """
21
22 dif=t1[0]+t1[1]-(t2[0]+t2[1])
23 if dif<0: return 1
24 elif dif>0: return -1
25 else: return 0
26
27 def ELMAN_RAMAGE(P):
28 return (P-1.).wherePositive()*0.5*(1.-1./(P+1.e-15))
29
30 def SIMPLIFIED_BROOK_HUGHES(P):
31 c=(P-3.).whereNegative()
32 return P/6.*c+1./2.*(1.-c)
33
34 def HALF(P):
35 return escript.Scalar(0.5,P.getFunctionSpace())
36
37 class PDECoefficient:
38 """
39 A class for PDE coefficients
40 """
41 # identifier for location of Data objects defining COEFFICIENTS
42 INTERIOR=0
43 BOUNDARY=1
44 CONTACT=2
45 CONTINUOUS=3
46 # identifier in the pattern of COEFFICIENTS:
47 # the pattern is a tuple of EQUATION,SOLUTION,DIM where DIM represents the spatial dimension, EQUATION the number of equations and SOLUTION the
48 # number of unknowns.
49 EQUATION=3
50 SOLUTION=4
51 DIM=5
52 # indicator for what is altered if the coefficient is altered:
53 OPERATOR=5
54 RIGHTHANDSIDE=6
55 BOTH=7
56 def __init__(self,where,pattern,altering):
57 """
58 Initialise a PDE Coefficient type
59 """
60 self.what=where
61 self.pattern=pattern
62 self.altering=altering
63 self.resetValue()
64
65 def resetValue(self):
66 """
67 resets coefficient value to default
68 """
69 self.value=escript.Data()
70
71 def getFunctionSpace(self,domain):
72 """
73 defines the FunctionSpace of the coefficient on the domain
74
75 @param domain:
76 """
77 if self.what==self.INTERIOR: return escript.Function(domain)
78 elif self.what==self.BOUNDARY: return escript.FunctionOnBoundary(domain)
79 elif self.what==self.CONTACT: return escript.FunctionOnContactZero(domain)
80 elif self.what==self.CONTINUOUS: return escript.ContinuousFunction(domain)
81
82 def getValue(self):
83 """
84 returns the value of the coefficient:
85 """
86 return self.value
87
88 def setValue(self,newValue):
89 """
90 set the value of the coefficient to new value
91 """
92 self.value=newValue
93
94 def isAlteringOperator(self):
95 """
96 return true if the operator of the PDE is changed when the coefficient is changed
97 """
98 if self.altering==self.OPERATOR or self.altering==self.BOTH:
99 return not None
100 else:
101 return None
102
103 def isAlteringRightHandSide(self):
104 """
105 return true if the right hand side of the PDE is changed when the coefficient is changed
106 """
107 if self.altering==self.RIGHTHANDSIDE or self.altering==self.BOTH:
108 return not None
109 else:
110 return None
111
112 def estimateNumEquationsAndNumSolutions(self,shape=(),dim=3):
113 """
114 tries to estimate the number of equations in a given tensor shape for a given spatial dimension dim
115
116 @param shape:
117 @param dim:
118 """
119 if len(shape)>0:
120 num=max(shape)+1
121 else:
122 num=1
123 search=[]
124 for u in range(num):
125 for e in range(num):
126 search.append((e,u))
127 search.sort(_CompTuple2)
128 for item in search:
129 s=self.buildShape(item[0],item[1],dim)
130 if len(s)==0 and len(shape)==0:
131 return (1,1)
132 else:
133 if s==shape: return item
134 return None
135
136 def buildShape(self,e=1,u=1,dim=3):
137 """
138 builds the required shape for a given number of equations e, number of unknowns u and spatial dimension dim
139
140 @param e:
141 @param u:
142 @param dim:
143 """
144 s=()
145 for i in self.pattern:
146 if i==self.EQUATION:
147 if e>1: s=s+(e,)
148 elif i==self.SOLUTION:
149 if u>1: s=s+(u,)
150 else:
151 s=s+(dim,)
152 return s
153
154 class LinearPDE:
155 """
156 Class to handle a linear PDE
157
158 class to define a linear PDE of the form
159
160 \f[
161 -(A_{ijkl}u_{k,l})_{,j} -(B_{ijk}u_k)_{,j} + C_{ikl}u_{k,l} +D_{ik}u_k = - (X_{ij})_{,j} + Y_i
162 \f]
163
164 with boundary conditons:
165
166 \f[
167 n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_{ik}u_k = - n_j*X_{ij} + y_i
168 \f]
169
170 and contact conditions
171
172 \f[
173 n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_contact_{ik}[u_k] = - n_j*X_{ij} + y_contact_i
174 \f]
175
176 and constraints:
177
178 \f[
179 u_i=r_i \quad \mathrm{where} \quad q_i>0
180 \f]
181
182 """
183 TOL=1.e-13
184 DEFAULT_METHOD=util.DEFAULT_METHOD
185 DIRECT=util.DIRECT
186 CHOLEVSKY=util.CHOLEVSKY
187 PCG=util.PCG
188 CR=util.CR
189 CGS=util.CGS
190 BICGSTAB=util.BICGSTAB
191 SSOR=util.SSOR
192 GMRES=util.GMRES
193 PRES20=util.PRES20
194 ILU0=util.ILU0
195 JACOBI=util.JACOBI
196
197 def __init__(self,domain,numEquations=0,numSolutions=0):
198 """
199 initializes a new linear PDE.
200
201 @param args:
202 """
203 # COEFFICIENTS can be overwritten by subclasses:
204 self.__COEFFICIENTS={
205 "A" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM,PDECoefficient.SOLUTION,PDECoefficient.DIM),PDECoefficient.OPERATOR),
206 "B" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),
207 "C" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION,PDECoefficient.DIM),PDECoefficient.OPERATOR),
208 "D" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),
209 "X" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM),PDECoefficient.RIGHTHANDSIDE),
210 "Y" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
211 "d" : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),
212 "y" : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
213 "d_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),
214 "y_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
215 "r" : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
216 "q" : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.SOLUTION,),PDECoefficient.BOTH)}
217
218 self.COEFFICIENTS=self.__COEFFICIENTS
219 # initialize attributes
220 self.__debug=None
221 self.__domain=domain
222 self.__numEquations=numEquations
223 self.__numSolutions=numSolutions
224 self.cleanCoefficients()
225
226 self.__operator=escript.Operator()
227 self.__operator_isValid=False
228 self.__righthandside=escript.Data()
229 self.__righthandside_isValid=False
230 self.__solution=escript.Data()
231 self.__solution_isValid=False
232
233 # set some default values:
234 self.__homogeneous_constraint=True
235 self.__row_function_space=escript.Solution(self.__domain)
236 self.__column_function_space=escript.Solution(self.__domain)
237 self.__tolerance=1.e-8
238 self.__solver_method=util.DEFAULT_METHOD
239 self.__matrix_type=self.__domain.getSystemMatrixTypeId(util.DEFAULT_METHOD,False)
240 self.__sym=False
241 self.__lumping=False
242
243 def createCoefficient(self, name):
244 """
245 create a data object corresponding to coefficient name
246 @param name:
247 """
248 return escript.Data(shape = getShapeOfCoefficient(name), \
249 what = getFunctionSpaceForCoefficient(name))
250
251 def __del__(self):
252 pass
253
254 def getCoefficient(self,name):
255 """
256 return the value of the parameter name
257
258 @param name:
259 """
260 return self.COEFFICIENTS[name].getValue()
261
262 def getCoefficientOfPDE(self,name):
263 """
264 return the value of the coefficient name of the general PDE.
265 This method is called by the assembling routine it can be
266 overwritten to map coefficients of a particualr PDE to the general PDE.
267
268 @param name:
269 """
270 return self.getCoefficient(name)
271
272 def hasCoefficient(self,name):
273 """
274 return true if name is the name of a coefficient
275
276 @param name:
277 """
278 return self.COEFFICIENTS.has_key(name)
279
280 def hasPDECoefficient(self,name):
281 """
282 return true if name is the name of a coefficient
283
284 @param name:
285 """
286 return self.__COEFFICIENTS.has_key(name)
287
288 def getFunctionSpaceForEquation(self):
289 """
290 return true if the test functions should use reduced order
291 """
292 return self.__row_function_space
293
294 def getFunctionSpaceForSolution(self):
295 """
296 return true if the interpolation of the solution should use reduced order
297 """
298 return self.__column_function_space
299
300 def setValue(self,**coefficients):
301 """
302 sets new values to coefficients
303
304 @param coefficients:
305 """
306 self.__setValue(**coefficients)
307
308
309 def cleanCoefficients(self):
310 """
311 resets all coefficients to default values.
312 """
313 for i in self.COEFFICIENTS.iterkeys():
314 self.COEFFICIENTS[i].resetValue()
315
316 def createNewCoefficient(self,name):
317 """
318 returns a new coefficient appropriate for coefficient name:
319 """
320 return escript.Data(0,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))
321
322 def createNewCoefficientOfPDE(self,name):
323 """
324 returns a new coefficient appropriate for coefficient name:
325 """
326 return escript.Data(0,self.getShapeOfCoefficientOfPDE(name),self.getFunctionSpaceForCoefficientOfPDE(name))
327
328 def getShapeOfCoefficientOfPDE(self,name):
329 """
330 return the shape of the coefficient name
331
332 @param name:
333 """
334 if self.hasPDECoefficient(name):
335 return self.__COEFFICIENTS[name].buildShape(self.getNumEquations(),self.getNumSolutions(),self.getDomain().getDim())
336 else:
337 raise ValueError,"Unknown coefficient %s requested"%name
338
339 def getFunctionSpaceForCoefficientOfPDE(self,name):
340 """
341 return the atoms of the coefficient name
342
343 @param name:
344 """
345 if self.hasPDECoefficient(name):
346 return self.__COEFFICIENTS[name].getFunctionSpace(self.getDomain())
347 else:
348 raise ValueError,"unknown coefficient %s requested"%name
349
350 def getShapeOfCoefficient(self,name):
351 """
352 return the shape of the coefficient name
353
354 @param name:
355 """
356 if self.hasCoefficient(name):
357 return self.COEFFICIENTS[name].buildShape(self.getNumEquations(),self.getNumSolutions(),self.getDomain().getDim())
358 else:
359 raise ValueError,"Unknown coefficient %s requested"%name
360
361 def getFunctionSpaceForCoefficient(self,name):
362 """
363 return the atoms of the coefficient name
364
365 @param name:
366 """
367 if self.hasCoefficient(name):
368 return self.COEFFICIENTS[name].getFunctionSpace(self.getDomain())
369 else:
370 raise ValueError,"unknown coefficient %s requested"%name
371
372 def alteredCoefficient(self,name):
373 """
374 announce that coefficient name has been changed
375
376 @param name:
377 """
378 if self.hasCoefficient(name):
379 if self.COEFFICIENTS[name].isAlteringOperator(): self.__rebuildOperator()
380 if self.COEFFICIENTS[name].isAlteringRightHandSide(): self.__rebuildRightHandSide()
381 else:
382 raise ValueError,"unknown coefficient %s requested"%name
383
384 # ===== debug ==============================================================
385 def setDebugOn(self):
386 """
387 """
388 self.__debug=not None
389
390 def setDebugOff(self):
391 """
392 """
393 self.__debug=None
394
395 def debug(self):
396 """
397 returns true if the PDE is in the debug mode
398 """
399 return self.__debug
400
401 #===== Lumping ===========================
402 def setLumpingOn(self):
403 """
404 indicates to use matrix lumping
405 """
406 if not self.isUsingLumping():
407 if self.debug() : print "PDE Debug: lumping is set on"
408 self.__rebuildOperator()
409 self.__lumping=True
410
411 def setLumpingOff(self):
412 """
413 switches off matrix lumping
414 """
415 if self.isUsingLumping():
416 if self.debug() : print "PDE Debug: lumping is set off"
417 self.__rebuildOperator()
418 self.__lumping=False
419
420 def setLumping(self,flag=False):
421 """
422 set the matrix lumping flag to flag
423 """
424 if flag:
425 self.setLumpingOn()
426 else:
427 self.setLumpingOff()
428
429 def isUsingLumping(self):
430 """
431
432 """
433 return self.__lumping
434
435 #============ method business =========================================================
436 def setSolverMethod(self,solver=util.DEFAULT_METHOD):
437 """
438 sets a new solver
439 """
440 if not solver==self.getSolverMethod():
441 self.__solver_method=solver
442 if self.debug() : print "PDE Debug: New solver is %s"%solver
443 self.__checkMatrixType()
444
445 def getSolverMethod(self):
446 """
447 returns the solver method
448 """
449 return self.__solver_method
450
451 #============ tolerance business =========================================================
452 def setTolerance(self,tol=1.e-8):
453 """
454 resets the tolerance to tol.
455 """
456 if not tol>0:
457 raise ValueException,"Tolerance as to be positive"
458 if tol<self.getTolerance(): self.__rebuildSolution()
459 if self.debug() : print "PDE Debug: New tolerance %e",tol
460 self.__tolerance=tol
461 return
462 def getTolerance(self):
463 """
464 returns the tolerance set for the solution
465 """
466 return self.__tolerance
467
468 #===== symmetry flag ==========================
469 def isSymmetric(self):
470 """
471 returns true is the operator is considered to be symmetric
472 """
473 return self.__sym
474
475 def setSymmetryOn(self):
476 """
477 sets the symmetry flag to true
478 """
479 if not self.isSymmetric():
480 if self.debug() : print "PDE Debug: Operator is set to be symmetric"
481 self.__sym=True
482 self.__checkMatrixType()
483
484 def setSymmetryOff(self):
485 """
486 sets the symmetry flag to false
487 """
488 if self.isSymmetric():
489 if self.debug() : print "PDE Debug: Operator is set to be unsymmetric"
490 self.__sym=False
491 self.__checkMatrixType()
492
493 def setSymmetryTo(self,flag=False):
494 """
495 sets the symmetry flag to flag
496
497 @param flag:
498 """
499 if flag:
500 self.setSymmetryOn()
501 else:
502 self.setSymmetryOff()
503
504 #===== order reduction ==========================
505 def setReducedOrderOn(self):
506 """
507 switches to on reduced order
508 """
509 self.setReducedOrderForSolutionOn()
510 self.setReducedOrderForEquationOn()
511
512 def setReducedOrderOff(self):
513 """
514 switches to full order
515 """
516 self.setReducedOrderForSolutionOff()
517 self.setReducedOrderForEquationOff()
518
519 def setReducedOrderTo(self,flag=False):
520 """
521 sets order according to flag
522
523 @param flag:
524 """
525 self.setReducedOrderForSolutionTo(flag)
526 self.setReducedOrderForEquationTo(flag)
527
528
529 #===== order reduction solution ==========================
530 def setReducedOrderForSolutionOn(self):
531 """
532 switches to reduced order to interpolate solution
533 """
534 new_fs=escript.ReducedSolution(self.getDomain())
535 if self.getFunctionSpaceForSolution()!=new_fs:
536 if self.debug() : print "PDE Debug: Reduced order is used to interpolate solution."
537 self.__column_function_space=new_fs
538 self.__rebuildSystem(deep=True)
539
540 def setReducedOrderForSolutionOff(self):
541 """
542 switches to full order to interpolate solution
543 """
544 new_fs=escript.Solution(self.getDomain())
545 if self.getFunctionSpaceForSolution()!=new_fs:
546 if self.debug() : print "PDE Debug: Full order is used to interpolate solution."
547 self.__column_function_space=new_fs
548 self.__rebuildSystem(deep=True)
549
550 def setReducedOrderForSolutionTo(self,flag=False):
551 """
552 sets order for test functions according to flag
553
554 @param flag:
555 """
556 if flag:
557 self.setReducedOrderForSolutionOn()
558 else:
559 self.setReducedOrderForSolutionOff()
560
561 #===== order reduction equation ==========================
562 def setReducedOrderForEquationOn(self):
563 """
564 switches to reduced order for test functions
565 """
566 new_fs=escript.ReducedSolution(self.getDomain())
567 if self.getFunctionSpaceForEquation()!=new_fs:
568 if self.debug() : print "PDE Debug: Reduced order is used for test functions."
569 self.__row_function_space=new_fs
570 self.__rebuildSystem(deep=True)
571
572 def setReducedOrderForEquationOff(self):
573 """
574 switches to full order for test functions
575 """
576 new_fs=escript.Solution(self.getDomain())
577 if self.getFunctionSpaceForEquation()!=new_fs:
578 if self.debug() : print "PDE Debug: Full order is used for test functions."
579 self.__row_function_space=new_fs
580 self.__rebuildSystem(deep=True)
581
582 def setReducedOrderForEquationTo(self,flag=False):
583 """
584 sets order for test functions according to flag
585
586 @param flag:
587 """
588 if flag:
589 self.setReducedOrderForEquationOn()
590 else:
591 self.setReducedOrderForEquationOff()
592
593 # ==== initialization =====================================================================
594 def __getNewOperator(self):
595 """
596 """
597 return self.getDomain().newOperator( \
598 self.getNumEquations(), \
599 self.getFunctionSpaceForEquation(), \
600 self.getNumSolutions(), \
601 self.getFunctionSpaceForSolution(), \
602 self.__matrix_type)
603
604 def __makeFreshRightHandSide(self):
605 """
606 """
607 if self.debug() : print "PDE Debug: New right hand side allocated"
608 if self.getNumEquations()>1:
609 self.__righthandside=escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)
610 else:
611 self.__righthandside=escript.Data(0.,(),self.getFunctionSpaceForEquation(),True)
612 return self.__righthandside
613
614 def __getNewSolution(self):
615 """
616 """
617 if self.debug() : print "PDE Debug: New right hand side allocated"
618 if self.getNumSolutions()>1:
619 return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)
620 else:
621 return escript.Data(0.,(),self.getFunctionSpaceForSolution(),True)
622
623 def __makeFreshOperator(self):
624 """
625 """
626 if self.__operator.isEmpty():
627 self.__operator=self.__getNewOperator()
628 if self.debug() : print "PDE Debug: New operator allocated"
629 else:
630 self.__operator.setValue(0.)
631 self.__operator.resetSolver()
632 if self.debug() : print "PDE Debug: Operator reset to zero"
633 return self.__operator
634
635 #============ some serivice functions =====================================================
636 def getDomain(self):
637 """
638 returns the domain of the PDE
639 """
640 return self.__domain
641
642 def getDim(self):
643 """
644 returns the spatial dimension of the PDE
645 """
646 return self.getDomain().getDim()
647
648 def getNumEquations(self):
649 """
650 returns the number of equations
651 """
652 if self.__numEquations>0:
653 return self.__numEquations
654 else:
655 raise ValueError,"Number of equations is undefined. Please specify argument numEquations."
656
657 def getNumSolutions(self):
658 """
659 returns the number of unknowns
660 """
661 if self.__numSolutions>0:
662 return self.__numSolutions
663 else:
664 raise ValueError,"Number of solution is undefined. Please specify argument numSolutions."
665
666
667 def checkSymmetry(self,verbose=True):
668 """
669 returns if the Operator is symmetric. This is a very expensive operation!!! The symmetry flag is not altered.
670 """
671 verbose=verbose or self.debug()
672 out=True
673 if self.getNumSolutions()!=self.getNumEquations():
674 if verbose: print "non-symmetric PDE because of different number of equations and solutions"
675 out=False
676 else:
677 A=self.getCoefficientOfPDE("A")
678 if not A.isEmpty():
679 tol=util.Lsup(A)*self.TOL
680 if self.getNumSolutions()>1:
681 for i in range(self.getNumEquations()):
682 for j in range(self.getDim()):
683 for k in range(self.getNumSolutions()):
684 for l in range(self.getDim()):
685 if util.Lsup(A[i,j,k,l]-A[k,l,i,j])>tol:
686 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)
687 out=False
688 else:
689 for j in range(self.getDim()):
690 for l in range(self.getDim()):
691 if util.Lsup(A[j,l]-A[l,j])>tol:
692 if verbose: print "non-symmetric PDE because A[%d,%d]!=A[%d,%d]"%(j,l,l,j)
693 out=False
694 B=self.getCoefficientOfPDE("B")
695 C=self.getCoefficientOfPDE("C")
696 if B.isEmpty() and not C.isEmpty():
697 if verbose: print "non-symmetric PDE because B is not present but C is"
698 out=False
699 elif not B.isEmpty() and C.isEmpty():
700 if verbose: print "non-symmetric PDE because C is not present but B is"
701 out=False
702 elif not B.isEmpty() and not C.isEmpty():
703 tol=(util.Lsup(B)+util.Lsup(C))*self.TOL/2.
704 if self.getNumSolutions()>1:
705 for i in range(self.getNumEquations()):
706 for j in range(self.getDim()):
707 for k in range(self.getNumSolutions()):
708 if util.Lsup(B[i,j,k]-C[k,i,j])>tol:
709 if verbose: print "non-symmetric PDE because B[%d,%d,%d]!=C[%d,%d,%d]"%(i,j,k,k,i,j)
710 out=False
711 else:
712 for j in range(self.getDim()):
713 if util.Lsup(B[j]-C[j])>tol:
714 if verbose: print "non-symmetric PDE because B[%d]!=C[%d]"%(j,j)
715 out=False
716 if self.getNumSolutions()>1:
717 D=self.getCoefficientOfPDE("D")
718 if not D.isEmpty():
719 tol=util.Lsup(D)*self.TOL
720 for i in range(self.getNumEquations()):
721 for k in range(self.getNumSolutions()):
722 if util.Lsup(D[i,k]-D[k,i])>tol:
723 if verbose: print "non-symmetric PDE because D[%d,%d]!=D[%d,%d]"%(i,k,k,i)
724 out=False
725
726 return out
727
728 def getFlux(self,u):
729 """
730 returns the flux J_ij for a given u
731
732 \f[
733 J_ij=A_{ijkl}u_{k,l}+B_{ijk}u_k-X_{ij}
734 \f]
735
736 @param u: argument of the operator
737 """
738 raise SystemError,"getFlux is not implemented yet"
739 return None
740
741 def applyOperator(self,u):
742 """
743 applies the operator of the PDE to a given solution u in weak from
744
745 @param u: argument of the operator
746 """
747 return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())
748
749 def getResidual(self,u):
750 """
751 return the residual of u in the weak from
752
753 @param u:
754 """
755 return self.applyOperator(u)-self.getRightHandSide()
756
757 def __setValue(self,**coefficients):
758 """
759 sets new values to coefficient
760
761 @param coefficients:
762 """
763 # check if the coefficients are legal:
764 for i in coefficients.iterkeys():
765 if not self.hasCoefficient(i):
766 raise ValueError,"Attempt to set unknown coefficient %s"%i
767 # if the number of unknowns or equations is still unknown we try to estimate them:
768 if self.__numEquations<1 or self.__numSolutions<1:
769 for i,d in coefficients.iteritems():
770 if hasattr(d,"shape"):
771 s=d.shape
772 elif hasattr(d,"getShape"):
773 s=d.getShape()
774 else:
775 s=numarray.array(d).shape
776 if s!=None:
777 # get number of equations and number of unknowns:
778 res=self.COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(s,self.getDim())
779 if res==None:
780 raise ValueError,"Illegal shape %s of coefficient %s"%(s,i)
781 else:
782 if self.__numEquations<1: self.__numEquations=res[0]
783 if self.__numSolutions<1: self.__numSolutions=res[1]
784 if self.__numEquations<1: raise ValueError,"unidententified number of equations"
785 if self.__numSolutions<1: raise ValueError,"unidententified number of solutions"
786 # now we check the shape of the coefficient if numEquations and numSolutions are set:
787 for i,d in coefficients.iteritems():
788 if d==None:
789 d2=escript.Data()
790 elif isinstance(d,escript.Data):
791 if d.isEmpty():
792 d2=d
793 else:
794 d2=escript.Data(d,self.getFunctionSpaceForCoefficient(i))
795 else:
796 d2=escript.Data(d,self.getFunctionSpaceForCoefficient(i))
797 if not d2.isEmpty():
798 if not self.getShapeOfCoefficient(i)==d2.getShape():
799 raise ValueError,"Expected shape for coefficient %s is %s but actual shape is %s."%(i,self.getShapeOfCoefficient(i),d2.getShape())
800 # overwrite new values:
801 if self.debug(): print "PDE Debug: Coefficient %s has been altered."%i
802 self.COEFFICIENTS[i].setValue(d2)
803 self.alteredCoefficient(i)
804
805 # reset the HomogeneousConstraintFlag:
806 self.__setHomogeneousConstraintFlag()
807 if len(coefficients)>0 and not self.isUsingLumping() and not self.__homogeneous_constraint: self.__rebuildSystem()
808
809 def __setHomogeneousConstraintFlag(self):
810 """
811 checks if the constraints are homogeneous and sets self.__homogeneous_constraint accordingly.
812 """
813 self.__homogeneous_constraint=True
814 q=self.getCoefficientOfPDE("q")
815 r=self.getCoefficientOfPDE("r")
816 if not q.isEmpty() and not r.isEmpty():
817 if (q*r).Lsup()>=1.e-13*r.Lsup(): self.__homogeneous_constraint=False
818 if self.debug():
819 if self.__homogeneous_constraint:
820 print "PDE Debug: Constraints are homogeneous."
821 else:
822 print "PDE Debug: Constraints are inhomogeneous."
823
824
825 # ==== rebuild switches =====================================================================
826 def __rebuildSolution(self,deep=False):
827 """
828 indicates the PDE has to be reolved if the solution is requested
829 """
830 if self.__solution_isValid and self.debug() : print "PDE Debug: PDE has to be resolved."
831 self.__solution_isValid=False
832 if deep: self.__solution=escript.Data()
833
834
835 def __rebuildOperator(self,deep=False):
836 """
837 indicates the operator has to be rebuilt next time it is used
838 """
839 if self.__operator_isValid and self.debug() : print "PDE Debug: Operator has to be rebuilt."
840 self.__rebuildSolution(deep)
841 self.__operator_isValid=False
842 if deep: self.__operator=escript.Operator()
843
844 def __rebuildRightHandSide(self,deep=False):
845 """
846 indicates the right hand side has to be rebuild next time it is used
847 """
848 if self.__righthandside_isValid and self.debug() : print "PDE Debug: Right hand side has to be rebuilt."
849 self.__rebuildSolution(deep)
850 self.__righthandside_isValid=False
851 if deep: self.__righthandside=escript.Data()
852
853 def __rebuildSystem(self,deep=False):
854 """
855 annonced that all coefficient name has been changed
856 """
857 self.__rebuildSolution(deep)
858 self.__rebuildOperator(deep)
859 self.__rebuildRightHandSide(deep)
860
861 def __checkMatrixType(self):
862 """
863 reassess the matrix type and, if needed, initiates an operator rebuild
864 """
865 new_matrix_type=self.getDomain().getSystemMatrixTypeId(self.getSolverMethod(),self.isSymmetric())
866 if not new_matrix_type==self.__matrix_type:
867 if self.debug() : print "PDE Debug: Matrix type is now %d."%new_matrix_type
868 self.__matrix_type=new_matrix_type
869 self.__rebuildOperator(deep=True)
870
871 #============ assembling =======================================================
872 def __copyConstraint(self):
873 """
874 copies the constrint condition into u
875 """
876 if not self.__righthandside.isEmpty():
877 q=self.getCoefficientOfPDE("q")
878 r=self.getCoefficientOfPDE("r")
879 if not q.isEmpty():
880 if r.isEmpty():
881 r2=escript.Data(0,self.__righthandside.getShape(),self.__righthandside.getFunctionSpace())
882 else:
883 r2=escript.Data(r,self.__righthandside.getFunctionSpace())
884 self.__righthandside.copyWithMask(r2,escript.Data(q,self.__righthandside.getFunctionSpace()))
885
886 def __applyConstraint(self):
887 """
888 applies the constraints defined by q and r to the system
889 """
890 q=self.getCoefficientOfPDE("q")
891 r=self.getCoefficientOfPDE("r")
892 if not q.isEmpty() and not self.__operator.isEmpty():
893 # q is the row and column mask to indicate where constraints are set:
894 row_q=escript.Data(q,self.getFunctionSpaceForEquation())
895 col_q=escript.Data(q,self.getFunctionSpaceForSolution())
896 u=self.__getNewSolution()
897 if r.isEmpty():
898 r_s=self.__getNewSolution()
899 else:
900 r_s=escript.Data(r,self.getFunctionSpaceForSolution())
901 u.copyWithMask(r_s,col_q)
902 if self.isUsingLumping():
903 self.__operator.copyWithMask(escript.Data(1,q.getShape(),self.getFunctionSpaceForEquation()),row_q)
904 else:
905 if not self.__righthandside.isEmpty(): self.__righthandside-=self.__operator*u
906 self.__operator.nullifyRowsAndCols(row_q,col_q,1.)
907
908 def getSystem(self):
909 """
910 return the operator and right hand side of the PDE
911 """
912 if not self.__operator_isValid or not self.__righthandside_isValid:
913 if self.isUsingLumping():
914 if not self.__operator_isValid:
915 if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():
916 raise TypeError,"Lumped matrix requires same order for equations and unknowns"
917 if not self.getCoefficientOfPDE("A").isEmpty():
918 raise Warning,"Lumped matrix does not allow coefficient A"
919 if not self.getCoefficientOfPDE("B").isEmpty():
920 raise Warning,"Lumped matrix does not allow coefficient B"
921 if not self.getCoefficientOfPDE("C").isEmpty():
922 raise Warning,"Lumped matrix does not allow coefficient C"
923 if self.debug() : print "PDE Debug: New lumped operator is built."
924 mat=self.__getNewOperator()
925 self.getDomain().addPDEToSystem(mat,escript.Data(), \
926 self.getCoefficientOfPDE("A"), \
927 self.getCoefficientOfPDE("B"), \
928 self.getCoefficientOfPDE("C"), \
929 self.getCoefficientOfPDE("D"), \
930 escript.Data(), \
931 escript.Data(), \
932 self.getCoefficientOfPDE("d"), \
933 escript.Data(),\
934 self.getCoefficientOfPDE("d_contact"), \
935 escript.Data())
936 self.__operator=mat*escript.Data(1,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)
937 self.__applyConstraint()
938 self.__operator_isValid=True
939 if not self.__righthandside_isValid:
940 if self.debug() : print "PDE Debug: New right hand side is built."
941 self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \
942 self.getCoefficientOfPDE("X"), \
943 self.getCoefficientOfPDE("Y"),\
944 self.getCoefficientOfPDE("y"),\
945 self.getCoefficientOfPDE("y_contact"))
946 self.__copyConstraint()
947 self.__righthandside_isValid=True
948 else:
949 if not self.__operator_isValid and not self.__righthandside_isValid:
950 if self.debug() : print "PDE Debug: New system is built."
951 self.getDomain().addPDEToSystem(self.__makeFreshOperator(),self.__makeFreshRightHandSide(), \
952 self.getCoefficientOfPDE("A"), \
953 self.getCoefficientOfPDE("B"), \
954 self.getCoefficientOfPDE("C"), \
955 self.getCoefficientOfPDE("D"), \
956 self.getCoefficientOfPDE("X"), \
957 self.getCoefficientOfPDE("Y"), \
958 self.getCoefficientOfPDE("d"), \
959 self.getCoefficientOfPDE("y"), \
960 self.getCoefficientOfPDE("d_contact"), \
961 self.getCoefficientOfPDE("y_contact"))
962 self.__applyConstraint()
963 self.__copyConstraint()
964 self.__operator_isValid=True
965 self.__righthandside_isValid=True
966 elif not self.__righthandside_isValid:
967 if self.debug() : print "PDE Debug: New right hand side is built."
968 self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \
969 self.getCoefficientOfPDE("X"), \
970 self.getCoefficientOfPDE("Y"),\
971 self.getCoefficientOfPDE("y"),\
972 self.getCoefficientOfPDE("y_contact"))
973 self.__copyConstraint()
974 self.__righthandside_isValid=True
975 elif not self.__operator_isValid:
976 if self.debug() : print "PDE Debug: New operator is built."
977 self.getDomain().addPDEToSystem(self.__makeFreshOperator(),escript.Data(), \
978 self.getCoefficientOfPDE("A"), \
979 self.getCoefficientOfPDE("B"), \
980 self.getCoefficientOfPDE("C"), \
981 self.getCoefficientOfPDE("D"), \
982 escript.Data(), \
983 escript.Data(), \
984 self.getCoefficientOfPDE("d"), \
985 escript.Data(),\
986 self.getCoefficientOfPDE("d_contact"), \
987 escript.Data())
988 self.__applyConstraint()
989 self.__operator_isValid=True
990 return (self.__operator,self.__righthandside)
991 def getOperator(self):
992 """
993 returns the operator of the PDE
994 """
995 return self.getSystem()[0]
996
997 def getRightHandSide(self):
998 """
999 returns the right hand side of the PDE
1000 """
1001 return self.getSystem()[1]
1002
1003 def solve(self,**options):
1004 """
1005 solve the PDE
1006
1007 @param options:
1008 """
1009 mat,f=self.getSystem()
1010 if self.isUsingLumping():
1011 out=f/mat
1012 else:
1013 options[util.TOLERANCE_KEY]=self.getTolerance()
1014 options[util.METHOD_KEY]=self.getSolverMethod()
1015 options[util.SYMMETRY_KEY]=self.isSymmetric()
1016 if self.debug() : print "PDE Debug: solver options: ",options
1017 out=mat.solve(f,options)
1018 return out
1019
1020 def getSolution(self,**options):
1021 """
1022 returns the solution of the PDE
1023
1024 @param options:
1025 """
1026 if not self.__solution_isValid:
1027 if self.debug() : print "PDE Debug: PDE is resolved."
1028 self.__solution=self.solve(**options)
1029 self.__solution_isValid=True
1030 return self.__solution
1031
1032
1033
1034 def ELMAN_RAMAGE(P):
1035 """ """
1036 return (P-1.).wherePositive()*0.5*(1.-1./(P+1.e-15))
1037 def SIMPLIFIED_BROOK_HUGHES(P):
1038 """ """
1039 c=(P-3.).whereNegative()
1040 return P/6.*c+1./2.*(1.-c)
1041
1042 def HALF(P):
1043 """ """
1044 return escript.Scalar(0.5,P.getFunctionSpace())
1045
1046 class AdvectivePDE(LinearPDE):
1047 """
1048 Class to handle a linear PDE dominated by advective terms:
1049
1050 class to define a linear PDE of the form
1051
1052 \f[
1053 -(A_{ijkl}u_{k,l})_{,j} -(B_{ijk}u_k)_{,j} + C_{ikl}u_{k,l} +D_{ik}u_k = - (X_{ij})_{,j} + Y_i
1054 \f]
1055
1056 with boundary conditons:
1057
1058 \f[
1059 n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_{ik}u_k = - n_j*X_{ij} + y_i
1060 \f]
1061
1062 and contact conditions
1063
1064 \f[
1065 n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d^{contact}_{ik}[u_k] = - n_j*X_{ij} + y^{contact}_{i}
1066 \f]
1067
1068 and constraints:
1069
1070 \f[
1071 u_i=r_i \quad \mathrm{where} \quad q_i>0
1072 \f]
1073 """
1074 def __init__(self,domain,numEquations=0,numSolutions=0,xi=ELMAN_RAMAGE):
1075 LinearPDE.__init__(self,domain,numEquations,numSolutions)
1076 self.__xi=xi
1077 self.__Xi=escript.Data()
1078
1079 def __calculateXi(self,peclet_factor,Z,h):
1080 Z_max=util.Lsup(Z)
1081 if Z_max>0.:
1082 return h*self.__xi(Z*peclet_factor)/(Z+Z_max*self.TOL)
1083 else:
1084 return 0.
1085
1086 def setValue(self,**args):
1087 if "A" in args.keys() or "B" in args.keys() or "C" in args.keys(): self.__Xi=escript.Data()
1088 self._LinearPDE__setValue(**args)
1089
1090 def getXi(self):
1091 if self.__Xi.isEmpty():
1092 B=self.getCoefficient("B")
1093 C=self.getCoefficient("C")
1094 A=self.getCoefficient("A")
1095 h=self.getDomain().getSize()
1096 self.__Xi=escript.Scalar(0.,self.getFunctionSpaceForCoefficient("A"))
1097 if not C.isEmpty() or not B.isEmpty():
1098 if not C.isEmpty() and not B.isEmpty():
1099 Z2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))
1100 if self.getNumEquations()>1:
1101 if self.getNumSolutions()>1:
1102 for i in range(self.getNumEquations()):
1103 for k in range(self.getNumSolutions()):
1104 for l in range(self.getDim()): Z2+=(C[i,k,l]-B[i,l,k])**2
1105 else:
1106 for i in range(self.getNumEquations()):
1107 for l in range(self.getDim()): Z2+=(C[i,l]-B[i,l])**2
1108 else:
1109 if self.getNumSolutions()>1:
1110 for k in range(self.getNumSolutions()):
1111 for l in range(self.getDim()): Z2+=(C[k,l]-B[l,k])**2
1112 else:
1113 for l in range(self.getDim()): Z2+=(C[l]-B[l])**2
1114 length_of_Z=util.sqrt(Z2)
1115 elif C.isEmpty():
1116 length_of_Z=util.length(B)
1117 else:
1118 length_of_Z=util.length(C)
1119
1120 Z_max=util.Lsup(length_of_Z)
1121 if Z_max>0.:
1122 length_of_A=util.length(A)
1123 A_max=util.Lsup(length_of_A)
1124 if A_max>0:
1125 inv_A=1./(length_of_A+A_max*self.TOL)
1126 else:
1127 inv_A=1./self.TOL
1128 peclet_number=length_of_Z*h/2*inv_A
1129 xi=self.__xi(peclet_number)
1130 self.__Xi=h*xi/(length_of_Z+Z_max*self.TOL)
1131 print "@ preclet number = %e"%util.Lsup(peclet_number),util.Lsup(xi),util.Lsup(length_of_Z)
1132 return self.__Xi
1133
1134
1135 def getCoefficientOfPDE(self,name):
1136 """
1137 return the value of the coefficient name of the general PDE
1138
1139 @param name:
1140 """
1141 if not self.getNumEquations() == self.getNumSolutions():
1142 raise ValueError,"AdvectivePDE expects the number of solution componets and the number of equations to be equal."
1143
1144 if name == "A" :
1145 A=self.getCoefficient("A")
1146 B=self.getCoefficient("B")
1147 C=self.getCoefficient("C")
1148 if B.isEmpty() and C.isEmpty():
1149 Aout=A
1150 else:
1151 if A.isEmpty():
1152 Aout=self.createNewCoefficient("A")
1153 else:
1154 Aout=A[:]
1155 Xi=self.getXi()
1156 if self.getNumEquations()>1:
1157 for i in range(self.getNumEquations()):
1158 for j in range(self.getDim()):
1159 for k in range(self.getNumSolutions()):
1160 for l in range(self.getDim()):
1161 if not C.isEmpty() and not B.isEmpty():
1162 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])
1163 elif C.isEmpty():
1164 for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*B[p,j,i]*B[p,l,k]
1165 else:
1166 for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*C[p,i,j]*C[p,k,l]
1167 else:
1168 for j in range(self.getDim()):
1169 for l in range(self.getDim()):
1170 if not C.isEmpty() and not B.isEmpty():
1171 Aout[j,l]+=Xi*(C[j]-B[j])*(C[l]-B[l])
1172 elif C.isEmpty():
1173 Aout[j,l]+=Xi*B[j]*B[l]
1174 else:
1175 Aout[j,l]+=Xi*C[j]*C[l]
1176 return Aout
1177 elif name == "B" :
1178 B=self.getCoefficient("B")
1179 C=self.getCoefficient("C")
1180 D=self.getCoefficient("D")
1181 if C.isEmpty() or D.isEmpty():
1182 Bout=B
1183 else:
1184 Xi=self.getXi()
1185 if B.isEmpty():
1186 Bout=self.createNewCoefficient("B")
1187 else:
1188 Bout=B[:]
1189 if self.getNumEquations()>1:
1190 for k in range(self.getNumSolutions()):
1191 for p in range(self.getNumEquations()):
1192 tmp=Xi*D[p,k]
1193 for i in range(self.getNumEquations()):
1194 for j in range(self.getDim()):
1195 Bout[i,j,k]+=tmp*C[p,i,j]
1196 else:
1197 tmp=Xi*D
1198 for j in range(self.getDim()): Bout[j]+=tmp*C[j]
1199 return Bout
1200 elif name == "C" :
1201 B=self.getCoefficient("B")
1202 C=self.getCoefficient("C")
1203 D=self.getCoefficient("D")
1204 if B.isEmpty() or D.isEmpty():
1205 Cout=C
1206 else:
1207 Xi=self.getXi()
1208 if C.isEmpty():
1209 Cout=self.createNewCoefficient("C")
1210 else:
1211 Cout=C[:]
1212 if self.getNumEquations()>1:
1213 for k in range(self.getNumSolutions()):
1214 for p in range(self.getNumEquations()):
1215 tmp=Xi*D[p,k]
1216 for i in range(self.getNumEquations()):
1217 for l in range(self.getDim()):
1218 Cout[i,k,l]+=tmp*B[p,l,i]
1219 else:
1220 tmp=Xi*D
1221 for j in range(self.getDim()): Cout[j]+=tmp*B[j]
1222 return Cout
1223 elif name == "D" :
1224 return self.getCoefficient("D")
1225 elif name == "X" :
1226 X=self.getCoefficient("X")
1227 Y=self.getCoefficient("Y")
1228 B=self.getCoefficient("B")
1229 C=self.getCoefficient("C")
1230 if Y.isEmpty() or (B.isEmpty() and C.isEmpty()):
1231 Xout=X
1232 else:
1233 if X.isEmpty():
1234 Xout=self.createNewCoefficient("X")
1235 else:
1236 Xout=X[:]
1237 Xi=self.getXi()
1238 if self.getNumEquations()>1:
1239 for p in range(self.getNumEquations()):
1240 tmp=Xi*Y[p]
1241 for i in range(self.getNumEquations()):
1242 for j in range(self.getDim()):
1243 if not C.isEmpty() and not B.isEmpty():
1244 Xout[i,j]+=tmp*(C[p,i,j]-B[p,j,i])
1245 elif C.isEmpty():
1246 Xout[i,j]-=tmp*B[p,j,i]
1247 else:
1248 Xout[i,j]+=tmp*C[p,i,j]
1249 else:
1250 tmp=Xi*Y
1251 for j in range(self.getDim()):
1252 if not C.isEmpty() and not B.isEmpty():
1253 Xout[j]+=tmp*(C[j]-B[j])
1254 elif C.isEmpty():
1255 Xout[j]-=tmp*B[j]
1256 else:
1257 Xout[j]+=tmp*C[j]
1258 return Xout
1259 elif name == "Y" :
1260 return self.getCoefficient("Y")
1261 elif name == "d" :
1262 return self.getCoefficient("d")
1263 elif name == "y" :
1264 return self.getCoefficient("y")
1265 elif name == "d_contact" :
1266 return self.getCoefficient("d_contact")
1267 elif name == "y_contact" :
1268 return self.getCoefficient("y_contact")
1269 elif name == "r" :
1270 return self.getCoefficient("r")
1271 elif name == "q" :
1272 return self.getCoefficient("q")
1273 else:
1274 raise SystemError,"unknown PDE coefficient %s",name
1275
1276
1277 class Poisson(LinearPDE):
1278 """
1279 Class to define a Poisson equation problem:
1280
1281 class to define a linear PDE of the form
1282 \f[
1283 -u_{,jj} = f
1284 \f]
1285
1286 with boundary conditons:
1287
1288 \f[
1289 n_j*u_{,j} = 0
1290 \f]
1291
1292 and constraints:
1293
1294 \f[
1295 u=0 \quad \mathrm{where} \quad q>0
1296 \f]
1297 """
1298
1299 def __init__(self,domain,f=escript.Data(),q=escript.Data()):
1300 LinearPDE.__init__(self,domain,1,1)
1301 self.COEFFICIENTS={
1302 "f" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1303 "q" : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.BOTH)}
1304 self.setSymmetryOn()
1305 self.setValue(f,q)
1306
1307 def setValue(self,f=escript.Data(),q=escript.Data()):
1308 """set value of PDE parameters f and q"""
1309 self._LinearPDE__setValue(f=f,q=q)
1310
1311 def getCoefficientOfPDE(self,name):
1312 """
1313 return the value of the coefficient name of the general PDE
1314
1315 @param name:
1316 """
1317 if name == "A" :
1318 return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))
1319 elif name == "B" :
1320 return escript.Data()
1321 elif name == "C" :
1322 return escript.Data()
1323 elif name == "D" :
1324 return escript.Data()
1325 elif name == "X" :
1326 return escript.Data()
1327 elif name == "Y" :
1328 return self.getCoefficient("f")
1329 elif name == "d" :
1330 return escript.Data()
1331 elif name == "y" :
1332 return escript.Data()
1333 elif name == "d_contact" :
1334 return escript.Data()
1335 elif name == "y_contact" :
1336 return escript.Data()
1337 elif name == "r" :
1338 return escript.Data()
1339 elif name == "q" :
1340 return self.getCoefficient("q")
1341 else:
1342 raise SystemError,"unknown PDE coefficient %s",name
1343
1344 class LameEquation(LinearPDE):
1345 """
1346 Class to define a Lame equation problem:
1347
1348 class to define a linear PDE of the form
1349 \f[
1350 -(\mu (u_{i,j}+u_{j,i}))_{,j} - \lambda u_{j,ji}} = F_i -\sigma_{ij,j}
1351 \f]
1352
1353 with boundary conditons:
1354
1355 \f[
1356 n_j(\mu (u_{i,j}+u_{j,i})-sigma_{ij}) + n_i\lambda u_{j,j} = f_i
1357 \f]
1358
1359 and constraints:
1360
1361 \f[
1362 u_i=r_i \quad \mathrm{where} \quad q_i>0
1363 \f]
1364 """
1365
1366 def __init__(self,domain,f=escript.Data(),q=escript.Data()):
1367 LinearPDE.__init__(self,domain,domain.getDim(),domain.getDim())
1368 self.COEFFICIENTS={
1369 "lame_lambda" : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),
1370 "lame_mu" : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),
1371 "F" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1372 "sigma" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM),PDECoefficient.RIGHTHANDSIDE),
1373 "f" : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1374 "r" : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.BOTH),
1375 "q" : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.BOTH)}
1376 self.setSymmetryOn()
1377
1378 def setValue(self,lame_lambda=escript.Data(),lame_mu=escript.Data(),F=escript.Data(),sigma=escript.Data(),f=escript.Data(),r=escript.Data(),q=escript.Data()):
1379 """set value of PDE parameters"""
1380 self._LinearPDE__setValue(lame_lambda=lame_lambda, \
1381 lame_mu=lame_mu, \
1382 F=F, \
1383 sigma=sigma, \
1384 f=f, \
1385 r=r, \
1386 q=q)
1387 def getCoefficientOfPDE(self,name):
1388 """
1389 return the value of the coefficient name of the general PDE
1390
1391 @param name:
1392 """
1393 if name == "A" :
1394 out =self.createNewCoefficientOfPDE("A")
1395 for i in range(self.getDim()):
1396 for j in range(self.getDim()):
1397 out[i,i,j,j] += self.getCoefficient("lame_lambda")
1398 out[i,j,j,i] += self.getCoefficient("lame_mu")
1399 out[i,j,i,j] += self.getCoefficient("lame_mu")
1400 return out
1401 elif name == "B" :
1402 return escript.Data()
1403 elif name == "C" :
1404 return escript.Data()
1405 elif name == "D" :
1406 return escript.Data()
1407 elif name == "X" :
1408 return self.getCoefficient("sigma")
1409 elif name == "Y" :
1410 return self.getCoefficient("F")
1411 elif name == "d" :
1412 return escript.Data()
1413 elif name == "y" :
1414 return self.getCoefficient("f")
1415 elif name == "d_contact" :
1416 return escript.Data()
1417 elif name == "y_contact" :
1418 return escript.Data()
1419 elif name == "r" :
1420 return self.getCoefficient("r")
1421 elif name == "q" :
1422 return self.getCoefficient("q")
1423 else:
1424 raise SystemError,"unknown PDE coefficient %s",name
1425
1426 # $Log$
1427 # Revision 1.10 2005/08/12 01:45:36 jgs
1428 # erge of development branch dev-02 back to main trunk on 2005-08-12
1429 #
1430 # Revision 1.9.2.1 2005/07/29 07:10:27 gross
1431 # new functions in util and a new pde type in linearPDEs
1432 #
1433 # Revision 1.1.2.25 2005/07/28 04:21:09 gross
1434 # Lame equation: (linear elastic, isotropic) added
1435 #
1436 # Revision 1.1.2.24 2005/07/22 06:37:11 gross
1437 # some extensions to modellib and linearPDEs
1438 #
1439 # Revision 1.1.2.23 2005/05/13 00:55:20 cochrane
1440 # Fixed up some docstrings. Moved module-level functions to top of file so
1441 # that epydoc and doxygen can pick them up properly.
1442 #
1443 # Revision 1.1.2.22 2005/05/12 11:41:30 gross
1444 # some basic Models have been added
1445 #
1446 # Revision 1.1.2.21 2005/05/12 07:16:12 cochrane
1447 # Moved ELMAN_RAMAGE, SIMPLIFIED_BROOK_HUGHES, and HALF functions to bottom of
1448 # file so that the AdvectivePDE class is picked up by doxygen. Some
1449 # reformatting of docstrings. Addition of code to make equations come out
1450 # as proper LaTeX.
1451 #
1452 # Revision 1.1.2.20 2005/04/15 07:09:08 gross
1453 # some problems with functionspace and linearPDEs fixed.
1454 #
1455 # Revision 1.1.2.19 2005/03/04 05:27:07 gross
1456 # bug in SystemPattern fixed.
1457 #
1458 # Revision 1.1.2.18 2005/02/08 06:16:45 gross
1459 # Bugs in AdvectivePDE fixed, AdvectiveTest is stable but more testing is needed
1460 #
1461 # Revision 1.1.2.17 2005/02/08 05:56:19 gross
1462 # Reference Number handling added
1463 #
1464 # Revision 1.1.2.16 2005/02/07 04:41:28 gross
1465 # some function exposed to python to make mesh merging running
1466 #
1467 # Revision 1.1.2.15 2005/02/03 00:14:44 gross
1468 # timeseries add and ESySParameter.py renames esysXML.py for consistence
1469 #
1470 # Revision 1.1.2.14 2005/02/01 06:44:10 gross
1471 # new implementation of AdvectivePDE which now also updates right hand side. systems of PDEs are still not working
1472 #
1473 # Revision 1.1.2.13 2005/01/25 00:47:07 gross
1474 # updates in the documentation
1475 #
1476 # Revision 1.1.2.12 2005/01/12 01:28:04 matt
1477 # Added createCoefficient method for linearPDEs.
1478 #
1479 # Revision 1.1.2.11 2005/01/11 01:55:34 gross
1480 # a problem in linearPDE class fixed
1481 #
1482 # Revision 1.1.2.10 2005/01/07 01:13:29 gross
1483 # some bugs in linearPDE fixed
1484 #
1485 # Revision 1.1.2.9 2005/01/06 06:24:58 gross
1486 # some bugs in slicing fixed
1487 #
1488 # Revision 1.1.2.8 2005/01/05 04:21:40 gross
1489 # FunctionSpace checking/matchig in slicing added
1490 #
1491 # Revision 1.1.2.7 2004/12/29 10:03:41 gross
1492 # bug in setValue fixed
1493 #
1494 # Revision 1.1.2.6 2004/12/29 05:29:59 gross
1495 # AdvectivePDE successfully tested for Peclet number 1000000. there is still a problem with setValue and Data()
1496 #
1497 # Revision 1.1.2.5 2004/12/29 00:18:41 gross
1498 # AdvectivePDE added
1499 #
1500 # Revision 1.1.2.4 2004/12/24 06:05:41 gross
1501 # some changes in linearPDEs to add AdevectivePDE
1502 #
1503 # Revision 1.1.2.3 2004/12/16 00:12:34 gross
1504 # __init__ of LinearPDE does not accept any coefficient anymore
1505 #
1506 # Revision 1.1.2.2 2004/12/14 03:55:01 jgs
1507 # *** empty log message ***
1508 #
1509 # Revision 1.1.2.1 2004/12/12 22:53:47 gross
1510 # linearPDE has been renamed LinearPDE
1511 #
1512 # Revision 1.1.1.1.2.7 2004/12/07 10:13:08 gross
1513 # GMRES added
1514 #
1515 # Revision 1.1.1.1.2.6 2004/12/07 03:19:50 gross
1516 # options for GMRES and PRES20 added
1517 #
1518 # Revision 1.1.1.1.2.5 2004/12/01 06:25:15 gross
1519 # some small changes
1520 #
1521 # Revision 1.1.1.1.2.4 2004/11/24 01:50:21 gross
1522 # Finley solves 4M unknowns now
1523 #
1524 # Revision 1.1.1.1.2.3 2004/11/15 06:05:26 gross
1525 # poisson solver added
1526 #
1527 # Revision 1.1.1.1.2.2 2004/11/12 06:58:15 gross
1528 # 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
1529 #
1530 # Revision 1.1.1.1.2.1 2004/10/28 22:59:22 gross
1531 # finley's RecTest.py is running now: problem in SystemMatrixAdapater fixed
1532 #
1533 # Revision 1.1.1.1 2004/10/26 06:53:56 jgs
1534 # initial import of project esys2
1535 #
1536 # Revision 1.3.2.3 2004/10/26 06:43:48 jgs
1537 # committing Lutz's and Paul's changes to brach jgs
1538 #
1539 # Revision 1.3.4.1 2004/10/20 05:32:51 cochrane
1540 # Added incomplete Doxygen comments to files, or merely put the docstrings that already exist into Doxygen form.
1541 #
1542 # Revision 1.3 2004/09/23 00:53:23 jgs
1543 # minor fixes
1544 #
1545 # Revision 1.1 2004/08/28 12:58:06 gross
1546 # SimpleSolve is not running yet: problem with == of functionsspace
1547 #
1548 #

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26