/[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 142 - (show annotations)
Mon Jul 25 05:28:20 2005 UTC (14 years, 2 months ago) by jgs
Original Path: trunk/esys2/escript/py_src/linearPDEs.py
File MIME type: text/x-python
File size: 54182 byte(s)
Merge of development branch back to main trunk on 2005-07-25

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26