/[escript]/trunk-mpi-branch/escript/py_src/linearPDEs.py
ViewVC logotype

Annotation of /trunk-mpi-branch/escript/py_src/linearPDEs.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 147 - (hide annotations)
Fri Aug 12 01:45:47 2005 UTC (14 years, 1 month 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 jgs 102 # $Id$
2    
3     ## @file linearPDEs.py
4    
5     """
6 jgs 122 Functions and classes for linear PDEs
7 jgs 102 """
8    
9     import escript
10     import util
11     import numarray
12    
13    
14     def _CompTuple2(t1,t2):
15     """
16 jgs 122 Compare two tuples
17 jgs 102
18 jgs 122 \param t1 The first tuple
19     \param t2 The second tuple
20 jgs 102 """
21 jgs 122
22 jgs 102 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 jgs 122 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 jgs 108 class PDECoefficient:
38 jgs 102 """
39 jgs 122 A class for PDE coefficients
40 jgs 102 """
41 jgs 108 # identifier for location of Data objects defining COEFFICIENTS
42 jgs 102 INTERIOR=0
43     BOUNDARY=1
44     CONTACT=2
45     CONTINUOUS=3
46 jgs 108 # identifier in the pattern of COEFFICIENTS:
47 jgs 102 # 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 jgs 122 Initialise a PDE Coefficient type
59 jgs 102 """
60     self.what=where
61     self.pattern=pattern
62     self.altering=altering
63 jgs 108 self.resetValue()
64 jgs 102
65 jgs 108 def resetValue(self):
66     """
67 jgs 122 resets coefficient value to default
68 jgs 108 """
69     self.value=escript.Data()
70    
71 jgs 102 def getFunctionSpace(self,domain):
72     """
73 jgs 122 defines the FunctionSpace of the coefficient on the domain
74 jgs 102
75 jgs 122 @param domain:
76 jgs 102 """
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 jgs 108 def getValue(self):
83     """
84 jgs 122 returns the value of the coefficient:
85 jgs 108 """
86     return self.value
87    
88     def setValue(self,newValue):
89     """
90 jgs 122 set the value of the coefficient to new value
91 jgs 108 """
92     self.value=newValue
93    
94 jgs 102 def isAlteringOperator(self):
95     """
96 jgs 122 return true if the operator of the PDE is changed when the coefficient is changed
97 jgs 102 """
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 jgs 122 return true if the right hand side of the PDE is changed when the coefficient is changed
106 jgs 102 """
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 jgs 122 tries to estimate the number of equations in a given tensor shape for a given spatial dimension dim
115 jgs 102
116 jgs 122 @param shape:
117     @param dim:
118 jgs 102 """
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 jgs 122 builds the required shape for a given number of equations e, number of unknowns u and spatial dimension dim
139 jgs 102
140 jgs 122 @param e:
141     @param u:
142     @param dim:
143 jgs 102 """
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 jgs 122 Class to handle a linear PDE
157 jgs 102
158     class to define a linear PDE of the form
159    
160 jgs 122 \f[
161 jgs 102 -(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 jgs 122 \f]
163 jgs 102
164 jgs 122 with boundary conditons:
165 jgs 102
166 jgs 122 \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 jgs 102
170 jgs 122 and contact conditions
171 jgs 102
172 jgs 122 \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 jgs 102
176 jgs 122 and constraints:
177 jgs 102
178 jgs 122 \f[
179     u_i=r_i \quad \mathrm{where} \quad q_i>0
180     \f]
181 jgs 102
182     """
183 jgs 108 TOL=1.e-13
184 jgs 102 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 jgs 114 ILU0=util.ILU0
195     JACOBI=util.JACOBI
196 jgs 102
197 jgs 104 def __init__(self,domain,numEquations=0,numSolutions=0):
198 jgs 102 """
199 jgs 122 initializes a new linear PDE.
200 jgs 102
201 jgs 122 @param args:
202 jgs 102 """
203 jgs 108 # COEFFICIENTS can be overwritten by subclasses:
204 jgs 147 self.__COEFFICIENTS={
205 jgs 108 "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 jgs 102
218 jgs 147 self.COEFFICIENTS=self.__COEFFICIENTS
219 jgs 102 # initialize attributes
220     self.__debug=None
221 jgs 104 self.__domain=domain
222     self.__numEquations=numEquations
223     self.__numSolutions=numSolutions
224 jgs 102 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 jgs 108 def createCoefficient(self, name):
244     """
245 jgs 122 create a data object corresponding to coefficient name
246     @param name:
247 jgs 108 """
248     return escript.Data(shape = getShapeOfCoefficient(name), \
249 jgs 110 what = getFunctionSpaceForCoefficient(name))
250 jgs 108
251     def __del__(self):
252     pass
253    
254 jgs 102 def getCoefficient(self,name):
255     """
256 jgs 122 return the value of the parameter name
257 jgs 102
258 jgs 122 @param name:
259 jgs 102 """
260 jgs 108 return self.COEFFICIENTS[name].getValue()
261 jgs 102
262 jgs 108 def getCoefficientOfPDE(self,name):
263     """
264 jgs 122 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 jgs 108 """
270     return self.getCoefficient(name)
271    
272     def hasCoefficient(self,name):
273 jgs 102 """
274 jgs 122 return true if name is the name of a coefficient
275 jgs 102
276 jgs 122 @param name:
277 jgs 102 """
278 jgs 108 return self.COEFFICIENTS.has_key(name)
279 jgs 102
280 jgs 147 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 jgs 108 def getFunctionSpaceForEquation(self):
289     """
290 jgs 122 return true if the test functions should use reduced order
291 jgs 108 """
292     return self.__row_function_space
293    
294     def getFunctionSpaceForSolution(self):
295     """
296 jgs 122 return true if the interpolation of the solution should use reduced order
297 jgs 108 """
298     return self.__column_function_space
299    
300     def setValue(self,**coefficients):
301 jgs 102 """
302 jgs 122 sets new values to coefficients
303 jgs 102
304 jgs 122 @param coefficients:
305 jgs 102 """
306 jgs 121 self.__setValue(**coefficients)
307 jgs 102
308    
309     def cleanCoefficients(self):
310     """
311 jgs 122 resets all coefficients to default values.
312 jgs 102 """
313 jgs 108 for i in self.COEFFICIENTS.iterkeys():
314     self.COEFFICIENTS[i].resetValue()
315 jgs 102
316 jgs 108 def createNewCoefficient(self,name):
317     """
318 jgs 122 returns a new coefficient appropriate for coefficient name:
319 jgs 108 """
320     return escript.Data(0,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))
321    
322 jgs 147 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 jgs 108
332 jgs 147 @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 jgs 102 def getShapeOfCoefficient(self,name):
351     """
352 jgs 122 return the shape of the coefficient name
353 jgs 102
354 jgs 122 @param name:
355 jgs 102 """
356     if self.hasCoefficient(name):
357 jgs 108 return self.COEFFICIENTS[name].buildShape(self.getNumEquations(),self.getNumSolutions(),self.getDomain().getDim())
358 jgs 102 else:
359 jgs 147 raise ValueError,"Unknown coefficient %s requested"%name
360 jgs 102
361 jgs 108 def getFunctionSpaceForCoefficient(self,name):
362 jgs 102 """
363 jgs 122 return the atoms of the coefficient name
364 jgs 102
365 jgs 122 @param name:
366 jgs 102 """
367     if self.hasCoefficient(name):
368 jgs 108 return self.COEFFICIENTS[name].getFunctionSpace(self.getDomain())
369 jgs 102 else:
370 jgs 147 raise ValueError,"unknown coefficient %s requested"%name
371 jgs 102
372     def alteredCoefficient(self,name):
373     """
374 jgs 122 announce that coefficient name has been changed
375 jgs 102
376 jgs 122 @param name:
377 jgs 102 """
378     if self.hasCoefficient(name):
379 jgs 108 if self.COEFFICIENTS[name].isAlteringOperator(): self.__rebuildOperator()
380     if self.COEFFICIENTS[name].isAlteringRightHandSide(): self.__rebuildRightHandSide()
381 jgs 102 else:
382 jgs 108 raise ValueError,"unknown coefficient %s requested"%name
383 jgs 102
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 jgs 122 returns true if the PDE is in the debug mode
398 jgs 102 """
399     return self.__debug
400    
401     #===== Lumping ===========================
402     def setLumpingOn(self):
403     """
404 jgs 122 indicates to use matrix lumping
405 jgs 102 """
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 jgs 122 switches off matrix lumping
414 jgs 102 """
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 jgs 122 set the matrix lumping flag to flag
423 jgs 102 """
424     if flag:
425     self.setLumpingOn()
426     else:
427     self.setLumpingOff()
428    
429     def isUsingLumping(self):
430     """
431 jgs 122
432 jgs 102 """
433     return self.__lumping
434    
435     #============ method business =========================================================
436     def setSolverMethod(self,solver=util.DEFAULT_METHOD):
437     """
438 jgs 122 sets a new solver
439 jgs 102 """
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 jgs 122 returns the solver method
448 jgs 102 """
449     return self.__solver_method
450    
451     #============ tolerance business =========================================================
452     def setTolerance(self,tol=1.e-8):
453     """
454 jgs 122 resets the tolerance to tol.
455 jgs 102 """
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 jgs 122 returns the tolerance set for the solution
465 jgs 102 """
466     return self.__tolerance
467    
468     #===== symmetry flag ==========================
469     def isSymmetric(self):
470     """
471 jgs 122 returns true is the operator is considered to be symmetric
472 jgs 102 """
473     return self.__sym
474    
475     def setSymmetryOn(self):
476     """
477 jgs 122 sets the symmetry flag to true
478 jgs 102 """
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 jgs 122 sets the symmetry flag to false
487 jgs 102 """
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 jgs 122 sets the symmetry flag to flag
496 jgs 102
497 jgs 122 @param flag:
498 jgs 102 """
499     if flag:
500     self.setSymmetryOn()
501     else:
502     self.setSymmetryOff()
503    
504     #===== order reduction ==========================
505     def setReducedOrderOn(self):
506     """
507 jgs 122 switches to on reduced order
508 jgs 102 """
509     self.setReducedOrderForSolutionOn()
510     self.setReducedOrderForEquationOn()
511    
512     def setReducedOrderOff(self):
513     """
514 jgs 122 switches to full order
515 jgs 102 """
516     self.setReducedOrderForSolutionOff()
517     self.setReducedOrderForEquationOff()
518    
519     def setReducedOrderTo(self,flag=False):
520     """
521 jgs 122 sets order according to flag
522 jgs 102
523 jgs 122 @param flag:
524 jgs 102 """
525     self.setReducedOrderForSolutionTo(flag)
526     self.setReducedOrderForEquationTo(flag)
527    
528    
529     #===== order reduction solution ==========================
530     def setReducedOrderForSolutionOn(self):
531     """
532 jgs 122 switches to reduced order to interpolate solution
533 jgs 102 """
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 jgs 122 switches to full order to interpolate solution
543 jgs 102 """
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 jgs 122 sets order for test functions according to flag
553 jgs 102
554 jgs 122 @param flag:
555 jgs 102 """
556     if flag:
557     self.setReducedOrderForSolutionOn()
558     else:
559     self.setReducedOrderForSolutionOff()
560    
561     #===== order reduction equation ==========================
562     def setReducedOrderForEquationOn(self):
563     """
564 jgs 122 switches to reduced order for test functions
565 jgs 102 """
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 jgs 122 switches to full order for test functions
575 jgs 102 """
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 jgs 122 sets order for test functions according to flag
585 jgs 102
586 jgs 122 @param flag:
587 jgs 102 """
588     if flag:
589     self.setReducedOrderForEquationOn()
590     else:
591     self.setReducedOrderForEquationOff()
592    
593     # ==== initialization =====================================================================
594 jgs 121 def __getNewOperator(self):
595 jgs 102 """
596     """
597     return self.getDomain().newOperator( \
598     self.getNumEquations(), \
599     self.getFunctionSpaceForEquation(), \
600     self.getNumSolutions(), \
601     self.getFunctionSpaceForSolution(), \
602     self.__matrix_type)
603    
604 jgs 121 def __makeFreshRightHandSide(self):
605 jgs 102 """
606     """
607 jgs 121 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 jgs 102
614 jgs 121 def __getNewSolution(self):
615 jgs 102 """
616     """
617 jgs 121 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 jgs 102
623 jgs 121 def __makeFreshOperator(self):
624 jgs 102 """
625     """
626     if self.__operator.isEmpty():
627 jgs 121 self.__operator=self.__getNewOperator()
628 jgs 102 if self.debug() : print "PDE Debug: New operator allocated"
629     else:
630     self.__operator.setValue(0.)
631 jgs 108 self.__operator.resetSolver()
632 jgs 102 if self.debug() : print "PDE Debug: Operator reset to zero"
633     return self.__operator
634    
635 jgs 108 #============ some serivice functions =====================================================
636     def getDomain(self):
637     """
638 jgs 122 returns the domain of the PDE
639 jgs 108 """
640     return self.__domain
641    
642     def getDim(self):
643     """
644 jgs 122 returns the spatial dimension of the PDE
645 jgs 108 """
646     return self.getDomain().getDim()
647    
648     def getNumEquations(self):
649     """
650 jgs 122 returns the number of equations
651 jgs 108 """
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 jgs 122 returns the number of unknowns
660 jgs 108 """
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 jgs 122 returns if the Operator is symmetric. This is a very expensive operation!!! The symmetry flag is not altered.
670 jgs 108 """
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 jgs 110 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 jgs 108 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 jgs 122 returns the flux J_ij for a given u
731 jgs 108
732 jgs 122 \f[
733     J_ij=A_{ijkl}u_{k,l}+B_{ijk}u_k-X_{ij}
734     \f]
735 jgs 108
736 jgs 122 @param u: argument of the operator
737 jgs 108 """
738     raise SystemError,"getFlux is not implemented yet"
739     return None
740    
741     def applyOperator(self,u):
742     """
743 jgs 122 applies the operator of the PDE to a given solution u in weak from
744 jgs 108
745 jgs 122 @param u: argument of the operator
746 jgs 108 """
747     return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())
748    
749     def getResidual(self,u):
750     """
751 jgs 122 return the residual of u in the weak from
752 jgs 108
753 jgs 122 @param u:
754 jgs 108 """
755     return self.applyOperator(u)-self.getRightHandSide()
756    
757 jgs 121 def __setValue(self,**coefficients):
758 jgs 108 """
759 jgs 122 sets new values to coefficient
760 jgs 108
761 jgs 122 @param coefficients:
762 jgs 108 """
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 jgs 122 checks if the constraints are homogeneous and sets self.__homogeneous_constraint accordingly.
812 jgs 108 """
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 jgs 102 # ==== rebuild switches =====================================================================
826     def __rebuildSolution(self,deep=False):
827     """
828 jgs 122 indicates the PDE has to be reolved if the solution is requested
829 jgs 102 """
830     if self.__solution_isValid and self.debug() : print "PDE Debug: PDE has to be resolved."
831     self.__solution_isValid=False
832 jgs 108 if deep: self.__solution=escript.Data()
833 jgs 102
834    
835     def __rebuildOperator(self,deep=False):
836     """
837 jgs 122 indicates the operator has to be rebuilt next time it is used
838 jgs 102 """
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 jgs 122 indicates the right hand side has to be rebuild next time it is used
847 jgs 102 """
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 jgs 122 annonced that all coefficient name has been changed
856 jgs 102 """
857     self.__rebuildSolution(deep)
858     self.__rebuildOperator(deep)
859     self.__rebuildRightHandSide(deep)
860    
861     def __checkMatrixType(self):
862     """
863 jgs 122 reassess the matrix type and, if needed, initiates an operator rebuild
864 jgs 102 """
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 jgs 108 def __copyConstraint(self):
873 jgs 102 """
874 jgs 122 copies the constrint condition into u
875 jgs 102 """
876 jgs 108 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 jgs 102
886 jgs 108 def __applyConstraint(self):
887 jgs 102 """
888 jgs 122 applies the constraints defined by q and r to the system
889 jgs 102 """
890 jgs 108 q=self.getCoefficientOfPDE("q")
891     r=self.getCoefficientOfPDE("r")
892 jgs 102 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 jgs 121 u=self.__getNewSolution()
897 jgs 102 if r.isEmpty():
898 jgs 121 r_s=self.__getNewSolution()
899 jgs 102 else:
900     r_s=escript.Data(r,self.getFunctionSpaceForSolution())
901     u.copyWithMask(r_s,col_q)
902 jgs 108 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 jgs 102
908 jgs 108 def getSystem(self):
909     """
910 jgs 122 return the operator and right hand side of the PDE
911 jgs 108 """
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 jgs 121 mat=self.__getNewOperator()
925 jgs 108 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 jgs 121 self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \
942 jgs 108 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 jgs 121 self.getDomain().addPDEToSystem(self.__makeFreshOperator(),self.__makeFreshRightHandSide(), \
952 jgs 108 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 jgs 121 self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \
969 jgs 108 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 jgs 121 self.getDomain().addPDEToSystem(self.__makeFreshOperator(),escript.Data(), \
978 jgs 108 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 jgs 102 def getOperator(self):
992     """
993 jgs 122 returns the operator of the PDE
994 jgs 102 """
995 jgs 108 return self.getSystem()[0]
996 jgs 102
997 jgs 108 def getRightHandSide(self):
998 jgs 102 """
999 jgs 122 returns the right hand side of the PDE
1000 jgs 102 """
1001 jgs 108 return self.getSystem()[1]
1002 jgs 102
1003     def solve(self,**options):
1004     """
1005 jgs 122 solve the PDE
1006 jgs 102
1007 jgs 122 @param options:
1008 jgs 102 """
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 jgs 122 returns the solution of the PDE
1023 jgs 102
1024 jgs 122 @param options:
1025 jgs 102 """
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 jgs 110
1033    
1034 jgs 122 def ELMAN_RAMAGE(P):
1035     """ """
1036     return (P-1.).wherePositive()*0.5*(1.-1./(P+1.e-15))
1037 jgs 110 def SIMPLIFIED_BROOK_HUGHES(P):
1038 jgs 122 """ """
1039     c=(P-3.).whereNegative()
1040     return P/6.*c+1./2.*(1.-c)
1041 jgs 110
1042 jgs 122 def HALF(P):
1043     """ """
1044     return escript.Scalar(0.5,P.getFunctionSpace())
1045 jgs 110
1046 jgs 108 class AdvectivePDE(LinearPDE):
1047     """
1048 jgs 122 Class to handle a linear PDE dominated by advective terms:
1049 jgs 108
1050     class to define a linear PDE of the form
1051 jgs 104
1052 jgs 122 \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 jgs 102
1056 jgs 122 with boundary conditons:
1057 jgs 102
1058 jgs 122 \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 jgs 102
1062 jgs 122 and contact conditions
1063 jgs 102
1064 jgs 122 \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 jgs 102
1068 jgs 122 and constraints:
1069 jgs 102
1070 jgs 122 \f[
1071     u_i=r_i \quad \mathrm{where} \quad q_i>0
1072     \f]
1073 jgs 110 """
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 jgs 102
1079 jgs 110 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 jgs 102
1086 jgs 110 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 jgs 121 self._LinearPDE__setValue(**args)
1089 jgs 110
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 jgs 102
1120 jgs 110 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 jgs 102
1135 jgs 108 def getCoefficientOfPDE(self,name):
1136     """
1137 jgs 122 return the value of the coefficient name of the general PDE
1138    
1139     @param name:
1140 jgs 108 """
1141 jgs 110 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 jgs 108 if name == "A" :
1145     A=self.getCoefficient("A")
1146     B=self.getCoefficient("B")
1147     C=self.getCoefficient("C")
1148 jgs 110 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 jgs 108 for k in range(self.getNumSolutions()):
1160 jgs 110 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 jgs 108 elif name == "B" :
1178 jgs 110 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 jgs 108 elif name == "C" :
1201 jgs 110 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 jgs 108 elif name == "D" :
1224     return self.getCoefficient("D")
1225     elif name == "X" :
1226 jgs 110 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 jgs 108 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 jgs 102 class Poisson(LinearPDE):
1278     """
1279 jgs 142 Class to define a Poisson equation problem:
1280 jgs 122
1281 jgs 102 class to define a linear PDE of the form
1282 jgs 122 \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 jgs 102 """
1298    
1299 jgs 108 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 jgs 102 self.setSymmetryOn()
1305     self.setValue(f,q)
1306    
1307     def setValue(self,f=escript.Data(),q=escript.Data()):
1308 jgs 142 """set value of PDE parameters f and q"""
1309 jgs 121 self._LinearPDE__setValue(f=f,q=q)
1310 jgs 108
1311     def getCoefficientOfPDE(self,name):
1312     """
1313 jgs 122 return the value of the coefficient name of the general PDE
1314    
1315     @param name:
1316 jgs 108 """
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 jgs 121
1344 jgs 142 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 jgs 147 -(\mu (u_{i,j}+u_{j,i}))_{,j} - \lambda u_{j,ji}} = F_i -\sigma_{ij,j}
1351 jgs 142 \f]
1352    
1353     with boundary conditons:
1354    
1355     \f[
1356 jgs 147 n_j(\mu (u_{i,j}+u_{j,i})-sigma_{ij}) + n_i\lambda u_{j,j} = f_i
1357 jgs 142 \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 jgs 147 "F" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1372     "sigma" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM),PDECoefficient.RIGHTHANDSIDE),
1373 jgs 142 "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 jgs 147 out =self.createNewCoefficientOfPDE("A")
1395 jgs 142 for i in range(self.getDim()):
1396     for j in range(self.getDim()):
1397 jgs 147 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 jgs 142 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 jgs 121 # $Log$
1427 jgs 147 # 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 jgs 142 #
1430 jgs 147 # 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 jgs 122 #
1433 jgs 147 # Revision 1.1.2.25 2005/07/28 04:21:09 gross
1434     # Lame equation: (linear elastic, isotropic) added
1435 jgs 121 #
1436 jgs 142 # Revision 1.1.2.24 2005/07/22 06:37:11 gross
1437     # some extensions to modellib and linearPDEs
1438     #
1439 jgs 122 # 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 jgs 121 # 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