/[escript]/trunk/esys2/escript/py_src/linearPDEs.py
ViewVC logotype

Annotation of /trunk/esys2/escript/py_src/linearPDEs.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 122 - (hide annotations)
Thu Jun 9 05:38:05 2005 UTC (14 years, 4 months ago) by jgs
File MIME type: text/x-python
File size: 50734 byte(s)
Merge of development branch back to main trunk on 2005-06-09

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     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 jgs 102
218     # initialize attributes
219     self.__debug=None
220 jgs 104 self.__domain=domain
221     self.__numEquations=numEquations
222     self.__numSolutions=numSolutions
223 jgs 102 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 jgs 108 def createCoefficient(self, name):
243     """
244 jgs 122 create a data object corresponding to coefficient name
245     @param name:
246 jgs 108 """
247     return escript.Data(shape = getShapeOfCoefficient(name), \
248 jgs 110 what = getFunctionSpaceForCoefficient(name))
249 jgs 108
250     def __del__(self):
251     pass
252    
253 jgs 102 def getCoefficient(self,name):
254     """
255 jgs 122 return the value of the parameter name
256 jgs 102
257 jgs 122 @param name:
258 jgs 102 """
259 jgs 108 return self.COEFFICIENTS[name].getValue()
260 jgs 102
261 jgs 108 def getCoefficientOfPDE(self,name):
262     """
263 jgs 122 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 jgs 108 """
269     return self.getCoefficient(name)
270    
271     def hasCoefficient(self,name):
272 jgs 102 """
273 jgs 122 return true if name is the name of a coefficient
274 jgs 102
275 jgs 122 @param name:
276 jgs 102 """
277 jgs 108 return self.COEFFICIENTS.has_key(name)
278 jgs 102
279 jgs 108 def getFunctionSpaceForEquation(self):
280     """
281 jgs 122 return true if the test functions should use reduced order
282 jgs 108 """
283     return self.__row_function_space
284    
285     def getFunctionSpaceForSolution(self):
286     """
287 jgs 122 return true if the interpolation of the solution should use reduced order
288 jgs 108 """
289     return self.__column_function_space
290    
291     def setValue(self,**coefficients):
292 jgs 102 """
293 jgs 122 sets new values to coefficients
294 jgs 102
295 jgs 122 @param coefficients:
296 jgs 102 """
297 jgs 121 self.__setValue(**coefficients)
298 jgs 102
299    
300     def cleanCoefficients(self):
301     """
302 jgs 122 resets all coefficients to default values.
303 jgs 102 """
304 jgs 108 for i in self.COEFFICIENTS.iterkeys():
305     self.COEFFICIENTS[i].resetValue()
306 jgs 102
307 jgs 108 def createNewCoefficient(self,name):
308     """
309 jgs 122 returns a new coefficient appropriate for coefficient name:
310 jgs 108 """
311     return escript.Data(0,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))
312    
313    
314 jgs 102 def getShapeOfCoefficient(self,name):
315     """
316 jgs 122 return the shape of the coefficient name
317 jgs 102
318 jgs 122 @param name:
319 jgs 102 """
320     if self.hasCoefficient(name):
321 jgs 108 return self.COEFFICIENTS[name].buildShape(self.getNumEquations(),self.getNumSolutions(),self.getDomain().getDim())
322 jgs 102 else:
323     raise ValueError,"Solution coefficient %s requested"%name
324    
325 jgs 108 def getFunctionSpaceForCoefficient(self,name):
326 jgs 102 """
327 jgs 122 return the atoms of the coefficient name
328 jgs 102
329 jgs 122 @param name:
330 jgs 102 """
331     if self.hasCoefficient(name):
332 jgs 108 return self.COEFFICIENTS[name].getFunctionSpace(self.getDomain())
333 jgs 102 else:
334     raise ValueError,"Solution coefficient %s requested"%name
335    
336     def alteredCoefficient(self,name):
337     """
338 jgs 122 announce that coefficient name has been changed
339 jgs 102
340 jgs 122 @param name:
341 jgs 102 """
342     if self.hasCoefficient(name):
343 jgs 108 if self.COEFFICIENTS[name].isAlteringOperator(): self.__rebuildOperator()
344     if self.COEFFICIENTS[name].isAlteringRightHandSide(): self.__rebuildRightHandSide()
345 jgs 102 else:
346 jgs 108 raise ValueError,"unknown coefficient %s requested"%name
347 jgs 102
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 jgs 122 returns true if the PDE is in the debug mode
362 jgs 102 """
363     return self.__debug
364    
365     #===== Lumping ===========================
366     def setLumpingOn(self):
367     """
368 jgs 122 indicates to use matrix lumping
369 jgs 102 """
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 jgs 122 switches off matrix lumping
378 jgs 102 """
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 jgs 122 set the matrix lumping flag to flag
387 jgs 102 """
388     if flag:
389     self.setLumpingOn()
390     else:
391     self.setLumpingOff()
392    
393     def isUsingLumping(self):
394     """
395 jgs 122
396 jgs 102 """
397     return self.__lumping
398    
399     #============ method business =========================================================
400     def setSolverMethod(self,solver=util.DEFAULT_METHOD):
401     """
402 jgs 122 sets a new solver
403 jgs 102 """
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 jgs 122 returns the solver method
412 jgs 102 """
413     return self.__solver_method
414    
415     #============ tolerance business =========================================================
416     def setTolerance(self,tol=1.e-8):
417     """
418 jgs 122 resets the tolerance to tol.
419 jgs 102 """
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 jgs 122 returns the tolerance set for the solution
429 jgs 102 """
430     return self.__tolerance
431    
432     #===== symmetry flag ==========================
433     def isSymmetric(self):
434     """
435 jgs 122 returns true is the operator is considered to be symmetric
436 jgs 102 """
437     return self.__sym
438    
439     def setSymmetryOn(self):
440     """
441 jgs 122 sets the symmetry flag to true
442 jgs 102 """
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 jgs 122 sets the symmetry flag to false
451 jgs 102 """
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 jgs 122 sets the symmetry flag to flag
460 jgs 102
461 jgs 122 @param flag:
462 jgs 102 """
463     if flag:
464     self.setSymmetryOn()
465     else:
466     self.setSymmetryOff()
467    
468     #===== order reduction ==========================
469     def setReducedOrderOn(self):
470     """
471 jgs 122 switches to on reduced order
472 jgs 102 """
473     self.setReducedOrderForSolutionOn()
474     self.setReducedOrderForEquationOn()
475    
476     def setReducedOrderOff(self):
477     """
478 jgs 122 switches to full order
479 jgs 102 """
480     self.setReducedOrderForSolutionOff()
481     self.setReducedOrderForEquationOff()
482    
483     def setReducedOrderTo(self,flag=False):
484     """
485 jgs 122 sets order according to flag
486 jgs 102
487 jgs 122 @param flag:
488 jgs 102 """
489     self.setReducedOrderForSolutionTo(flag)
490     self.setReducedOrderForEquationTo(flag)
491    
492    
493     #===== order reduction solution ==========================
494     def setReducedOrderForSolutionOn(self):
495     """
496 jgs 122 switches to reduced order to interpolate solution
497 jgs 102 """
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 jgs 122 switches to full order to interpolate solution
507 jgs 102 """
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 jgs 122 sets order for test functions according to flag
517 jgs 102
518 jgs 122 @param flag:
519 jgs 102 """
520     if flag:
521     self.setReducedOrderForSolutionOn()
522     else:
523     self.setReducedOrderForSolutionOff()
524    
525     #===== order reduction equation ==========================
526     def setReducedOrderForEquationOn(self):
527     """
528 jgs 122 switches to reduced order for test functions
529 jgs 102 """
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 jgs 122 switches to full order for test functions
539 jgs 102 """
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 jgs 122 sets order for test functions according to flag
549 jgs 102
550 jgs 122 @param flag:
551 jgs 102 """
552     if flag:
553     self.setReducedOrderForEquationOn()
554     else:
555     self.setReducedOrderForEquationOff()
556    
557     # ==== initialization =====================================================================
558 jgs 121 def __getNewOperator(self):
559 jgs 102 """
560     """
561     return self.getDomain().newOperator( \
562     self.getNumEquations(), \
563     self.getFunctionSpaceForEquation(), \
564     self.getNumSolutions(), \
565     self.getFunctionSpaceForSolution(), \
566     self.__matrix_type)
567    
568 jgs 121 def __makeFreshRightHandSide(self):
569 jgs 102 """
570     """
571 jgs 121 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 jgs 102
578 jgs 121 def __getNewSolution(self):
579 jgs 102 """
580     """
581 jgs 121 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 jgs 102
587 jgs 121 def __makeFreshOperator(self):
588 jgs 102 """
589     """
590     if self.__operator.isEmpty():
591 jgs 121 self.__operator=self.__getNewOperator()
592 jgs 102 if self.debug() : print "PDE Debug: New operator allocated"
593     else:
594     self.__operator.setValue(0.)
595 jgs 108 self.__operator.resetSolver()
596 jgs 102 if self.debug() : print "PDE Debug: Operator reset to zero"
597     return self.__operator
598    
599 jgs 108 #============ some serivice functions =====================================================
600     def getDomain(self):
601     """
602 jgs 122 returns the domain of the PDE
603 jgs 108 """
604     return self.__domain
605    
606     def getDim(self):
607     """
608 jgs 122 returns the spatial dimension of the PDE
609 jgs 108 """
610     return self.getDomain().getDim()
611    
612     def getNumEquations(self):
613     """
614 jgs 122 returns the number of equations
615 jgs 108 """
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 jgs 122 returns the number of unknowns
624 jgs 108 """
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 jgs 122 returns if the Operator is symmetric. This is a very expensive operation!!! The symmetry flag is not altered.
634 jgs 108 """
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 jgs 110 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 jgs 108 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 jgs 122 returns the flux J_ij for a given u
695 jgs 108
696 jgs 122 \f[
697     J_ij=A_{ijkl}u_{k,l}+B_{ijk}u_k-X_{ij}
698     \f]
699 jgs 108
700 jgs 122 @param u: argument of the operator
701 jgs 108 """
702     raise SystemError,"getFlux is not implemented yet"
703     return None
704    
705     def applyOperator(self,u):
706     """
707 jgs 122 applies the operator of the PDE to a given solution u in weak from
708 jgs 108
709 jgs 122 @param u: argument of the operator
710 jgs 108 """
711     return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())
712    
713     def getResidual(self,u):
714     """
715 jgs 122 return the residual of u in the weak from
716 jgs 108
717 jgs 122 @param u:
718 jgs 108 """
719     return self.applyOperator(u)-self.getRightHandSide()
720    
721 jgs 121 def __setValue(self,**coefficients):
722 jgs 108 """
723 jgs 122 sets new values to coefficient
724 jgs 108
725 jgs 122 @param coefficients:
726 jgs 108 """
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 jgs 122 checks if the constraints are homogeneous and sets self.__homogeneous_constraint accordingly.
776 jgs 108 """
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 jgs 102 # ==== rebuild switches =====================================================================
790     def __rebuildSolution(self,deep=False):
791     """
792 jgs 122 indicates the PDE has to be reolved if the solution is requested
793 jgs 102 """
794     if self.__solution_isValid and self.debug() : print "PDE Debug: PDE has to be resolved."
795     self.__solution_isValid=False
796 jgs 108 if deep: self.__solution=escript.Data()
797 jgs 102
798    
799     def __rebuildOperator(self,deep=False):
800     """
801 jgs 122 indicates the operator has to be rebuilt next time it is used
802 jgs 102 """
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 jgs 122 indicates the right hand side has to be rebuild next time it is used
811 jgs 102 """
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 jgs 122 annonced that all coefficient name has been changed
820 jgs 102 """
821     self.__rebuildSolution(deep)
822     self.__rebuildOperator(deep)
823     self.__rebuildRightHandSide(deep)
824    
825     def __checkMatrixType(self):
826     """
827 jgs 122 reassess the matrix type and, if needed, initiates an operator rebuild
828 jgs 102 """
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 jgs 108 def __copyConstraint(self):
837 jgs 102 """
838 jgs 122 copies the constrint condition into u
839 jgs 102 """
840 jgs 108 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 jgs 102
850 jgs 108 def __applyConstraint(self):
851 jgs 102 """
852 jgs 122 applies the constraints defined by q and r to the system
853 jgs 102 """
854 jgs 108 q=self.getCoefficientOfPDE("q")
855     r=self.getCoefficientOfPDE("r")
856 jgs 102 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 jgs 121 u=self.__getNewSolution()
861 jgs 102 if r.isEmpty():
862 jgs 121 r_s=self.__getNewSolution()
863 jgs 102 else:
864     r_s=escript.Data(r,self.getFunctionSpaceForSolution())
865     u.copyWithMask(r_s,col_q)
866 jgs 108 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 jgs 102
872 jgs 108 def getSystem(self):
873     """
874 jgs 122 return the operator and right hand side of the PDE
875 jgs 108 """
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 jgs 121 mat=self.__getNewOperator()
889 jgs 108 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 jgs 121 self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \
906 jgs 108 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 jgs 121 self.getDomain().addPDEToSystem(self.__makeFreshOperator(),self.__makeFreshRightHandSide(), \
916 jgs 108 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 jgs 121 self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \
933 jgs 108 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 jgs 121 self.getDomain().addPDEToSystem(self.__makeFreshOperator(),escript.Data(), \
942 jgs 108 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 jgs 102 def getOperator(self):
956     """
957 jgs 122 returns the operator of the PDE
958 jgs 102 """
959 jgs 108 return self.getSystem()[0]
960 jgs 102
961 jgs 108 def getRightHandSide(self):
962 jgs 102 """
963 jgs 122 returns the right hand side of the PDE
964 jgs 102 """
965 jgs 108 return self.getSystem()[1]
966 jgs 102
967     def solve(self,**options):
968     """
969 jgs 122 solve the PDE
970 jgs 102
971 jgs 122 @param options:
972 jgs 102 """
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 jgs 122 returns the solution of the PDE
987 jgs 102
988 jgs 122 @param options:
989 jgs 102 """
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 jgs 110
997    
998 jgs 122 def ELMAN_RAMAGE(P):
999     """ """
1000     return (P-1.).wherePositive()*0.5*(1.-1./(P+1.e-15))
1001 jgs 110 def SIMPLIFIED_BROOK_HUGHES(P):
1002 jgs 122 """ """
1003     c=(P-3.).whereNegative()
1004     return P/6.*c+1./2.*(1.-c)
1005 jgs 110
1006 jgs 122 def HALF(P):
1007     """ """
1008     return escript.Scalar(0.5,P.getFunctionSpace())
1009 jgs 110
1010 jgs 108 class AdvectivePDE(LinearPDE):
1011     """
1012 jgs 122 Class to handle a linear PDE dominated by advective terms:
1013 jgs 108
1014     class to define a linear PDE of the form
1015 jgs 104
1016 jgs 122 \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 jgs 102
1020 jgs 122 with boundary conditons:
1021 jgs 102
1022 jgs 122 \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 jgs 102
1026 jgs 122 and contact conditions
1027 jgs 102
1028 jgs 122 \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 jgs 102
1032 jgs 122 and constraints:
1033 jgs 102
1034 jgs 122 \f[
1035     u_i=r_i \quad \mathrm{where} \quad q_i>0
1036     \f]
1037 jgs 110 """
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 jgs 102
1043 jgs 110 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 jgs 102
1050 jgs 110 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 jgs 121 self._LinearPDE__setValue(**args)
1053 jgs 110
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 jgs 102
1084 jgs 110 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 jgs 102
1099 jgs 108 def getCoefficientOfPDE(self,name):
1100     """
1101 jgs 122 return the value of the coefficient name of the general PDE
1102    
1103     @param name:
1104 jgs 108 """
1105 jgs 110 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 jgs 108 if name == "A" :
1109     A=self.getCoefficient("A")
1110     B=self.getCoefficient("B")
1111     C=self.getCoefficient("C")
1112 jgs 110 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 jgs 108 for k in range(self.getNumSolutions()):
1124 jgs 110 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 jgs 108 elif name == "B" :
1142 jgs 110 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 jgs 108 elif name == "C" :
1165 jgs 110 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 jgs 108 elif name == "D" :
1188     return self.getCoefficient("D")
1189     elif name == "X" :
1190 jgs 110 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 jgs 108 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 jgs 102 class Poisson(LinearPDE):
1242     """
1243 jgs 122 Class to define a Poisson equstion problem:
1244    
1245 jgs 102 class to define a linear PDE of the form
1246 jgs 122 \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 jgs 102 """
1262    
1263 jgs 108 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 jgs 102 self.setSymmetryOn()
1269     self.setValue(f,q)
1270    
1271     def setValue(self,f=escript.Data(),q=escript.Data()):
1272 jgs 121 self._LinearPDE__setValue(f=f,q=q)
1273 jgs 108
1274     def getCoefficientOfPDE(self,name):
1275     """
1276 jgs 122 return the value of the coefficient name of the general PDE
1277    
1278     @param name:
1279 jgs 108 """
1280     if name == "A" :
1281     return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))
1282     elif name == "B" :
1283     return escript.Data()
1284     elif name == "C" :
1285     return escript.Data()
1286     elif name == "D" :
1287     return escript.Data()
1288     elif name == "X" :
1289     return escript.Data()
1290     elif name == "Y" :
1291     return self.getCoefficient("f")
1292     elif name == "d" :
1293     return escript.Data()
1294     elif name == "y" :
1295     return escript.Data()
1296     elif name == "d_contact" :
1297     return escript.Data()
1298     elif name == "y_contact" :
1299     return escript.Data()
1300     elif name == "r" :
1301     return escript.Data()
1302     elif name == "q" :
1303     return self.getCoefficient("q")
1304     else:
1305     raise SystemError,"unknown PDE coefficient %s",name
1306 jgs 121
1307     # $Log$
1308 jgs 122 # Revision 1.8 2005/06/09 05:37:59 jgs
1309     # Merge of development branch back to main trunk on 2005-06-09
1310     #
1311 jgs 121 # Revision 1.7 2005/05/06 04:26:10 jgs
1312     # Merge of development branch back to main trunk on 2005-05-06
1313     #
1314 jgs 122 # Revision 1.1.2.23 2005/05/13 00:55:20 cochrane
1315     # Fixed up some docstrings. Moved module-level functions to top of file so
1316     # that epydoc and doxygen can pick them up properly.
1317     #
1318     # Revision 1.1.2.22 2005/05/12 11:41:30 gross
1319     # some basic Models have been added
1320     #
1321     # Revision 1.1.2.21 2005/05/12 07:16:12 cochrane
1322     # Moved ELMAN_RAMAGE, SIMPLIFIED_BROOK_HUGHES, and HALF functions to bottom of
1323     # file so that the AdvectivePDE class is picked up by doxygen. Some
1324     # reformatting of docstrings. Addition of code to make equations come out
1325     # as proper LaTeX.
1326     #
1327 jgs 121 # Revision 1.1.2.20 2005/04/15 07:09:08 gross
1328     # some problems with functionspace and linearPDEs fixed.
1329     #
1330     # Revision 1.1.2.19 2005/03/04 05:27:07 gross
1331     # bug in SystemPattern fixed.
1332     #
1333     # Revision 1.1.2.18 2005/02/08 06:16:45 gross
1334     # Bugs in AdvectivePDE fixed, AdvectiveTest is stable but more testing is needed
1335     #
1336     # Revision 1.1.2.17 2005/02/08 05:56:19 gross
1337     # Reference Number handling added
1338     #
1339     # Revision 1.1.2.16 2005/02/07 04:41:28 gross
1340     # some function exposed to python to make mesh merging running
1341     #
1342     # Revision 1.1.2.15 2005/02/03 00:14:44 gross
1343     # timeseries add and ESySParameter.py renames esysXML.py for consistence
1344     #
1345     # Revision 1.1.2.14 2005/02/01 06:44:10 gross
1346     # new implementation of AdvectivePDE which now also updates right hand side. systems of PDEs are still not working
1347     #
1348     # Revision 1.1.2.13 2005/01/25 00:47:07 gross
1349     # updates in the documentation
1350     #
1351     # Revision 1.1.2.12 2005/01/12 01:28:04 matt
1352     # Added createCoefficient method for linearPDEs.
1353     #
1354     # Revision 1.1.2.11 2005/01/11 01:55:34 gross
1355     # a problem in linearPDE class fixed
1356     #
1357     # Revision 1.1.2.10 2005/01/07 01:13:29 gross
1358     # some bugs in linearPDE fixed
1359     #
1360     # Revision 1.1.2.9 2005/01/06 06:24:58 gross
1361     # some bugs in slicing fixed
1362     #
1363     # Revision 1.1.2.8 2005/01/05 04:21:40 gross
1364     # FunctionSpace checking/matchig in slicing added
1365     #
1366     # Revision 1.1.2.7 2004/12/29 10:03:41 gross
1367     # bug in setValue fixed
1368     #
1369     # Revision 1.1.2.6 2004/12/29 05:29:59 gross
1370     # AdvectivePDE successfully tested for Peclet number 1000000. there is still a problem with setValue and Data()
1371     #
1372     # Revision 1.1.2.5 2004/12/29 00:18:41 gross
1373     # AdvectivePDE added
1374     #
1375     # Revision 1.1.2.4 2004/12/24 06:05:41 gross
1376     # some changes in linearPDEs to add AdevectivePDE
1377     #
1378     # Revision 1.1.2.3 2004/12/16 00:12:34 gross
1379     # __init__ of LinearPDE does not accept any coefficient anymore
1380     #
1381     # Revision 1.1.2.2 2004/12/14 03:55:01 jgs
1382     # *** empty log message ***
1383     #
1384     # Revision 1.1.2.1 2004/12/12 22:53:47 gross
1385     # linearPDE has been renamed LinearPDE
1386     #
1387     # Revision 1.1.1.1.2.7 2004/12/07 10:13:08 gross
1388     # GMRES added
1389     #
1390     # Revision 1.1.1.1.2.6 2004/12/07 03:19:50 gross
1391     # options for GMRES and PRES20 added
1392     #
1393     # Revision 1.1.1.1.2.5 2004/12/01 06:25:15 gross
1394     # some small changes
1395     #
1396     # Revision 1.1.1.1.2.4 2004/11/24 01:50:21 gross
1397     # Finley solves 4M unknowns now
1398     #
1399     # Revision 1.1.1.1.2.3 2004/11/15 06:05:26 gross
1400     # poisson solver added
1401     #
1402     # Revision 1.1.1.1.2.2 2004/11/12 06:58:15 gross
1403     # 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
1404     #
1405     # Revision 1.1.1.1.2.1 2004/10/28 22:59:22 gross
1406     # finley's RecTest.py is running now: problem in SystemMatrixAdapater fixed
1407     #
1408     # Revision 1.1.1.1 2004/10/26 06:53:56 jgs
1409     # initial import of project esys2
1410     #
1411     # Revision 1.3.2.3 2004/10/26 06:43:48 jgs
1412     # committing Lutz's and Paul's changes to brach jgs
1413     #
1414     # Revision 1.3.4.1 2004/10/20 05:32:51 cochrane
1415     # Added incomplete Doxygen comments to files, or merely put the docstrings that already exist into Doxygen form.
1416     #
1417     # Revision 1.3 2004/09/23 00:53:23 jgs
1418     # minor fixes
1419     #
1420     # Revision 1.1 2004/08/28 12:58:06 gross
1421     # SimpleSolve is not running yet: problem with == of functionsspace
1422     #
1423     #

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26