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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 121 - (hide annotations)
Fri May 6 04:26:16 2005 UTC (14 years, 4 months ago) by jgs
Original Path: trunk/esys2/escript/py_src/linearPDEs.py
File MIME type: text/x-python
File size: 51196 byte(s)
Merge of development branch back to main trunk on 2005-05-06

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26