/[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 110 - (hide annotations)
Mon Feb 14 04:14:42 2005 UTC (14 years, 9 months ago) by jgs
Original Path: trunk/esys2/escript/py_src/linearPDEs.py
File MIME type: text/x-python
File size: 47700 byte(s)
*** empty log message ***

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    
176 jgs 104 def __init__(self,domain,numEquations=0,numSolutions=0):
177 jgs 102 """
178     @brief initializes a new linear PDE.
179    
180     @param args
181     """
182 jgs 108 # COEFFICIENTS can be overwritten by subclasses:
183     self.COEFFICIENTS={
184     "A" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM,PDECoefficient.SOLUTION,PDECoefficient.DIM),PDECoefficient.OPERATOR),
185     "B" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),
186     "C" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION,PDECoefficient.DIM),PDECoefficient.OPERATOR),
187     "D" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),
188     "X" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM),PDECoefficient.RIGHTHANDSIDE),
189     "Y" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
190     "d" : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),
191     "y" : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
192     "d_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),
193     "y_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
194     "r" : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
195     "q" : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.SOLUTION,),PDECoefficient.BOTH)}
196 jgs 102
197     # initialize attributes
198     self.__debug=None
199 jgs 104 self.__domain=domain
200     self.__numEquations=numEquations
201     self.__numSolutions=numSolutions
202 jgs 102 self.cleanCoefficients()
203    
204     self.__operator=escript.Operator()
205     self.__operator_isValid=False
206     self.__righthandside=escript.Data()
207     self.__righthandside_isValid=False
208     self.__solution=escript.Data()
209     self.__solution_isValid=False
210    
211     # set some default values:
212     self.__homogeneous_constraint=True
213     self.__row_function_space=escript.Solution(self.__domain)
214     self.__column_function_space=escript.Solution(self.__domain)
215     self.__tolerance=1.e-8
216     self.__solver_method=util.DEFAULT_METHOD
217     self.__matrix_type=self.__domain.getSystemMatrixTypeId(util.DEFAULT_METHOD,False)
218     self.__sym=False
219     self.__lumping=False
220    
221 jgs 108 def createCoefficient(self, name):
222     """
223     @brief create a data object corresponding to coefficient name
224     @param name
225     """
226     return escript.Data(shape = getShapeOfCoefficient(name), \
227 jgs 110 what = getFunctionSpaceForCoefficient(name))
228 jgs 108
229     def __del__(self):
230     pass
231    
232 jgs 102 def getCoefficient(self,name):
233     """
234 jgs 108 @brief return the value of the parameter name
235 jgs 102
236     @param name
237     """
238 jgs 108 return self.COEFFICIENTS[name].getValue()
239 jgs 102
240 jgs 108 def getCoefficientOfPDE(self,name):
241     """
242     @brief return the value of the coefficient name of the general PDE. This method is called by the assembling routine
243     it can be overwritten to map coefficients of a particualr PDE to the general PDE.
244     @param name
245     """
246     return self.getCoefficient(name)
247    
248     def hasCoefficient(self,name):
249 jgs 102 """
250 jgs 108 @brief return true if name is the name of a coefficient
251 jgs 102
252 jgs 108 @param name
253 jgs 102 """
254 jgs 108 return self.COEFFICIENTS.has_key(name)
255 jgs 102
256 jgs 108 def getFunctionSpaceForEquation(self):
257     """
258     @brief return true if the test functions should use reduced order
259     """
260     return self.__row_function_space
261    
262     def getFunctionSpaceForSolution(self):
263     """
264     @brief return true if the interpolation of the solution should use reduced order
265     """
266     return self.__column_function_space
267    
268     def setValue(self,**coefficients):
269 jgs 102 """
270     @brief sets new values to coefficients
271    
272     @param coefficients
273     """
274 jgs 108 self._setValue(**coefficients)
275 jgs 102
276    
277     def cleanCoefficients(self):
278     """
279     @brief resets all coefficients to default values.
280     """
281 jgs 108 for i in self.COEFFICIENTS.iterkeys():
282     self.COEFFICIENTS[i].resetValue()
283 jgs 102
284 jgs 108 def createNewCoefficient(self,name):
285     """
286     @brief returns a new coefficient appropriate for coefficient name:
287     """
288     return escript.Data(0,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))
289    
290    
291 jgs 102 def getShapeOfCoefficient(self,name):
292     """
293     @brief return the shape of the coefficient name
294    
295     @param name
296     """
297     if self.hasCoefficient(name):
298 jgs 108 return self.COEFFICIENTS[name].buildShape(self.getNumEquations(),self.getNumSolutions(),self.getDomain().getDim())
299 jgs 102 else:
300     raise ValueError,"Solution coefficient %s requested"%name
301    
302 jgs 108 def getFunctionSpaceForCoefficient(self,name):
303 jgs 102 """
304     @brief return the atoms of the coefficient name
305    
306     @param name
307     """
308     if self.hasCoefficient(name):
309 jgs 108 return self.COEFFICIENTS[name].getFunctionSpace(self.getDomain())
310 jgs 102 else:
311     raise ValueError,"Solution coefficient %s requested"%name
312    
313     def alteredCoefficient(self,name):
314     """
315     @brief annonced that coefficient name has been changed
316    
317     @param name
318     """
319     if self.hasCoefficient(name):
320 jgs 108 if self.COEFFICIENTS[name].isAlteringOperator(): self.__rebuildOperator()
321     if self.COEFFICIENTS[name].isAlteringRightHandSide(): self.__rebuildRightHandSide()
322 jgs 102 else:
323 jgs 108 raise ValueError,"unknown coefficient %s requested"%name
324 jgs 102
325     # ===== debug ==============================================================
326     def setDebugOn(self):
327     """
328     @brief
329     """
330     self.__debug=not None
331    
332     def setDebugOff(self):
333     """
334     @brief
335     """
336     self.__debug=None
337    
338     def debug(self):
339     """
340     @brief returns true if the PDE is in the debug mode
341     """
342     return self.__debug
343    
344     #===== Lumping ===========================
345     def setLumpingOn(self):
346     """
347     @brief indicates to use matrix lumping
348     """
349     if not self.isUsingLumping():
350     if self.debug() : print "PDE Debug: lumping is set on"
351     self.__rebuildOperator()
352     self.__lumping=True
353    
354     def setLumpingOff(self):
355     """
356     @brief switches off matrix lumping
357     """
358     if self.isUsingLumping():
359     if self.debug() : print "PDE Debug: lumping is set off"
360     self.__rebuildOperator()
361     self.__lumping=False
362    
363     def setLumping(self,flag=False):
364     """
365     @brief set the matrix lumping flag to flag
366     """
367     if flag:
368     self.setLumpingOn()
369     else:
370     self.setLumpingOff()
371    
372     def isUsingLumping(self):
373     """
374     @brief
375     """
376     return self.__lumping
377    
378     #============ method business =========================================================
379     def setSolverMethod(self,solver=util.DEFAULT_METHOD):
380     """
381     @brief sets a new solver
382     """
383     if not solver==self.getSolverMethod():
384     self.__solver_method=solver
385     if self.debug() : print "PDE Debug: New solver is %s"%solver
386     self.__checkMatrixType()
387    
388     def getSolverMethod(self):
389     """
390     @brief returns the solver method
391     """
392     return self.__solver_method
393    
394     #============ tolerance business =========================================================
395     def setTolerance(self,tol=1.e-8):
396     """
397     @brief resets the tolerance to tol.
398     """
399     if not tol>0:
400     raise ValueException,"Tolerance as to be positive"
401     if tol<self.getTolerance(): self.__rebuildSolution()
402     if self.debug() : print "PDE Debug: New tolerance %e",tol
403     self.__tolerance=tol
404     return
405     def getTolerance(self):
406     """
407     @brief returns the tolerance set for the solution
408     """
409     return self.__tolerance
410    
411     #===== symmetry flag ==========================
412     def isSymmetric(self):
413     """
414     @brief returns true is the operator is considered to be symmetric
415     """
416     return self.__sym
417    
418     def setSymmetryOn(self):
419     """
420     @brief sets the symmetry flag to true
421     """
422     if not self.isSymmetric():
423     if self.debug() : print "PDE Debug: Operator is set to be symmetric"
424     self.__sym=True
425     self.__checkMatrixType()
426    
427     def setSymmetryOff(self):
428     """
429     @brief sets the symmetry flag to false
430     """
431     if self.isSymmetric():
432     if self.debug() : print "PDE Debug: Operator is set to be unsymmetric"
433     self.__sym=False
434     self.__checkMatrixType()
435    
436     def setSymmetryTo(self,flag=False):
437     """
438     @brief sets the symmetry flag to flag
439    
440     @param flag
441     """
442     if flag:
443     self.setSymmetryOn()
444     else:
445     self.setSymmetryOff()
446    
447     #===== order reduction ==========================
448     def setReducedOrderOn(self):
449     """
450     @brief switches to on reduced order
451     """
452     self.setReducedOrderForSolutionOn()
453     self.setReducedOrderForEquationOn()
454    
455     def setReducedOrderOff(self):
456     """
457     @brief switches to full order
458     """
459     self.setReducedOrderForSolutionOff()
460     self.setReducedOrderForEquationOff()
461    
462     def setReducedOrderTo(self,flag=False):
463     """
464     @brief sets order according to flag
465    
466     @param flag
467     """
468     self.setReducedOrderForSolutionTo(flag)
469     self.setReducedOrderForEquationTo(flag)
470    
471    
472     #===== order reduction solution ==========================
473     def setReducedOrderForSolutionOn(self):
474     """
475     @brief switches to reduced order to interpolate solution
476     """
477     new_fs=escript.ReducedSolution(self.getDomain())
478     if self.getFunctionSpaceForSolution()!=new_fs:
479     if self.debug() : print "PDE Debug: Reduced order is used to interpolate solution."
480     self.__column_function_space=new_fs
481     self.__rebuildSystem(deep=True)
482    
483     def setReducedOrderForSolutionOff(self):
484     """
485     @brief switches to full order to interpolate solution
486     """
487     new_fs=escript.Solution(self.getDomain())
488     if self.getFunctionSpaceForSolution()!=new_fs:
489     if self.debug() : print "PDE Debug: Full order is used to interpolate solution."
490     self.__column_function_space=new_fs
491     self.__rebuildSystem(deep=True)
492    
493     def setReducedOrderForSolutionTo(self,flag=False):
494     """
495     @brief sets order for test functions according to flag
496    
497     @param flag
498     """
499     if flag:
500     self.setReducedOrderForSolutionOn()
501     else:
502     self.setReducedOrderForSolutionOff()
503    
504     #===== order reduction equation ==========================
505     def setReducedOrderForEquationOn(self):
506     """
507     @brief switches to reduced order for test functions
508     """
509     new_fs=escript.ReducedSolution(self.getDomain())
510     if self.getFunctionSpaceForEquation()!=new_fs:
511     if self.debug() : print "PDE Debug: Reduced order is used for test functions."
512     self.__row_function_space=new_fs
513     self.__rebuildSystem(deep=True)
514    
515     def setReducedOrderForEquationOff(self):
516     """
517     @brief switches to full order for test functions
518     """
519     new_fs=escript.Solution(self.getDomain())
520     if self.getFunctionSpaceForEquation()!=new_fs:
521     if self.debug() : print "PDE Debug: Full order is used for test functions."
522     self.__row_function_space=new_fs
523     self.__rebuildSystem(deep=True)
524    
525     def setReducedOrderForEquationTo(self,flag=False):
526     """
527     @brief sets order for test functions according to flag
528    
529     @param flag
530     """
531     if flag:
532     self.setReducedOrderForEquationOn()
533     else:
534     self.setReducedOrderForEquationOff()
535    
536     # ==== initialization =====================================================================
537     def __makeNewOperator(self):
538     """
539     @brief
540     """
541     return self.getDomain().newOperator( \
542     self.getNumEquations(), \
543     self.getFunctionSpaceForEquation(), \
544     self.getNumSolutions(), \
545     self.getFunctionSpaceForSolution(), \
546     self.__matrix_type)
547    
548     def __makeNewRightHandSide(self):
549     """
550     @brief
551     """
552     return escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)
553    
554     def __makeNewSolution(self):
555     """
556     @brief
557     """
558     return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)
559    
560     def __getFreshOperator(self):
561     """
562     @brief
563     """
564     if self.__operator.isEmpty():
565     self.__operator=self.__makeNewOperator()
566     if self.debug() : print "PDE Debug: New operator allocated"
567     else:
568     self.__operator.setValue(0.)
569 jgs 108 self.__operator.resetSolver()
570 jgs 102 if self.debug() : print "PDE Debug: Operator reset to zero"
571     return self.__operator
572    
573     def __getFreshRightHandSide(self):
574     """
575     @brief
576     """
577     if self.__righthandside.isEmpty():
578     self.__righthandside=self.__makeNewRightHandSide()
579     if self.debug() : print "PDE Debug: New right hand side allocated"
580     else:
581     print "fix self.__righthandside*=0"
582     self.__righthandside*=0.
583     if self.debug() : print "PDE Debug: Right hand side reset to zero"
584     return self.__righthandside
585    
586 jgs 108 #============ some serivice functions =====================================================
587     def getDomain(self):
588     """
589     @brief returns the domain of the PDE
590     """
591     return self.__domain
592    
593     def getDim(self):
594     """
595     @brief returns the spatial dimension of the PDE
596     """
597     return self.getDomain().getDim()
598    
599     def getNumEquations(self):
600     """
601     @brief returns the number of equations
602     """
603     if self.__numEquations>0:
604     return self.__numEquations
605     else:
606     raise ValueError,"Number of equations is undefined. Please specify argument numEquations."
607    
608     def getNumSolutions(self):
609     """
610     @brief returns the number of unknowns
611     """
612     if self.__numSolutions>0:
613     return self.__numSolutions
614     else:
615     raise ValueError,"Number of solution is undefined. Please specify argument numSolutions."
616    
617    
618     def checkSymmetry(self,verbose=True):
619     """
620     @brief returns if the Operator is symmetric. This is a very expensive operation!!! The symmetry flag is not altered.
621     """
622     verbose=verbose or self.debug()
623     out=True
624     if self.getNumSolutions()!=self.getNumEquations():
625     if verbose: print "non-symmetric PDE because of different number of equations and solutions"
626     out=False
627     else:
628     A=self.getCoefficientOfPDE("A")
629     if not A.isEmpty():
630     tol=util.Lsup(A)*self.TOL
631     if self.getNumSolutions()>1:
632     for i in range(self.getNumEquations()):
633     for j in range(self.getDim()):
634     for k in range(self.getNumSolutions()):
635     for l in range(self.getDim()):
636     if util.Lsup(A[i,j,k,l]-A[k,l,i,j])>tol:
637     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)
638     out=False
639     else:
640     for j in range(self.getDim()):
641     for l in range(self.getDim()):
642     if util.Lsup(A[j,l]-A[l,j])>tol:
643     if verbose: print "non-symmetric PDE because A[%d,%d]!=A[%d,%d]"%(j,l,l,j)
644     out=False
645     B=self.getCoefficientOfPDE("B")
646     C=self.getCoefficientOfPDE("C")
647     if B.isEmpty() and not C.isEmpty():
648     if verbose: print "non-symmetric PDE because B is not present but C is"
649     out=False
650     elif not B.isEmpty() and C.isEmpty():
651     if verbose: print "non-symmetric PDE because C is not present but B is"
652     out=False
653     elif not B.isEmpty() and not C.isEmpty():
654     tol=(util.Lsup(B)+util.Lsup(C))*self.TOL/2.
655     if self.getNumSolutions()>1:
656     for i in range(self.getNumEquations()):
657     for j in range(self.getDim()):
658     for k in range(self.getNumSolutions()):
659 jgs 110 if util.Lsup(B[i,j,k]-C[k,i,j])>tol:
660     if verbose: print "non-symmetric PDE because B[%d,%d,%d]!=C[%d,%d,%d]"%(i,j,k,k,i,j)
661 jgs 108 out=False
662     else:
663     for j in range(self.getDim()):
664     if util.Lsup(B[j]-C[j])>tol:
665     if verbose: print "non-symmetric PDE because B[%d]!=C[%d]"%(j,j)
666     out=False
667     if self.getNumSolutions()>1:
668     D=self.getCoefficientOfPDE("D")
669     if not D.isEmpty():
670     tol=util.Lsup(D)*self.TOL
671     for i in range(self.getNumEquations()):
672     for k in range(self.getNumSolutions()):
673     if util.Lsup(D[i,k]-D[k,i])>tol:
674     if verbose: print "non-symmetric PDE because D[%d,%d]!=D[%d,%d]"%(i,k,k,i)
675     out=False
676    
677     return out
678    
679     def getFlux(self,u):
680     """
681     @brief returns the flux J_ij for a given u
682    
683     J_ij=A_{ijkl}u_{k,l}+B_{ijk}u_k-X_{ij}
684    
685     @param u argument of the operator
686    
687     """
688     raise SystemError,"getFlux is not implemented yet"
689     return None
690    
691     def applyOperator(self,u):
692     """
693     @brief applies the operator of the PDE to a given solution u in weak from
694    
695     @param u argument of the operator
696    
697     """
698     return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())
699    
700     def getResidual(self,u):
701     """
702     @brief return the residual of u in the weak from
703    
704     @param u
705     """
706     return self.applyOperator(u)-self.getRightHandSide()
707    
708     def _setValue(self,**coefficients):
709     """
710     @brief sets new values to coefficient
711    
712     @param coefficients
713     """
714     # check if the coefficients are legal:
715     for i in coefficients.iterkeys():
716     if not self.hasCoefficient(i):
717     raise ValueError,"Attempt to set unknown coefficient %s"%i
718     # if the number of unknowns or equations is still unknown we try to estimate them:
719     if self.__numEquations<1 or self.__numSolutions<1:
720     for i,d in coefficients.iteritems():
721     if hasattr(d,"shape"):
722     s=d.shape
723     elif hasattr(d,"getShape"):
724     s=d.getShape()
725     else:
726     s=numarray.array(d).shape
727     if s!=None:
728     # get number of equations and number of unknowns:
729     res=self.COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(s,self.getDim())
730     if res==None:
731     raise ValueError,"Illegal shape %s of coefficient %s"%(s,i)
732     else:
733     if self.__numEquations<1: self.__numEquations=res[0]
734     if self.__numSolutions<1: self.__numSolutions=res[1]
735     if self.__numEquations<1: raise ValueError,"unidententified number of equations"
736     if self.__numSolutions<1: raise ValueError,"unidententified number of solutions"
737     # now we check the shape of the coefficient if numEquations and numSolutions are set:
738     for i,d in coefficients.iteritems():
739     if d==None:
740     d2=escript.Data()
741     elif isinstance(d,escript.Data):
742     if d.isEmpty():
743     d2=d
744     else:
745     d2=escript.Data(d,self.getFunctionSpaceForCoefficient(i))
746     else:
747     d2=escript.Data(d,self.getFunctionSpaceForCoefficient(i))
748     if not d2.isEmpty():
749     if not self.getShapeOfCoefficient(i)==d2.getShape():
750     raise ValueError,"Expected shape for coefficient %s is %s but actual shape is %s."%(i,self.getShapeOfCoefficient(i),d2.getShape())
751     # overwrite new values:
752     if self.debug(): print "PDE Debug: Coefficient %s has been altered."%i
753     self.COEFFICIENTS[i].setValue(d2)
754     self.alteredCoefficient(i)
755    
756     # reset the HomogeneousConstraintFlag:
757     self.__setHomogeneousConstraintFlag()
758     if len(coefficients)>0 and not self.isUsingLumping() and not self.__homogeneous_constraint: self.__rebuildSystem()
759    
760     def __setHomogeneousConstraintFlag(self):
761     """
762     @brief checks if the constraints are homogeneous and sets self.__homogeneous_constraint accordingly.
763     """
764     self.__homogeneous_constraint=True
765     q=self.getCoefficientOfPDE("q")
766     r=self.getCoefficientOfPDE("r")
767     if not q.isEmpty() and not r.isEmpty():
768     if (q*r).Lsup()>=1.e-13*r.Lsup(): self.__homogeneous_constraint=False
769     if self.debug():
770     if self.__homogeneous_constraint:
771     print "PDE Debug: Constraints are homogeneous."
772     else:
773     print "PDE Debug: Constraints are inhomogeneous."
774    
775    
776 jgs 102 # ==== rebuild switches =====================================================================
777     def __rebuildSolution(self,deep=False):
778     """
779     @brief indicates the PDE has to be reolved if the solution is requested
780     """
781     if self.__solution_isValid and self.debug() : print "PDE Debug: PDE has to be resolved."
782     self.__solution_isValid=False
783 jgs 108 if deep: self.__solution=escript.Data()
784 jgs 102
785    
786     def __rebuildOperator(self,deep=False):
787     """
788     @brief indicates the operator has to be rebuilt next time it is used
789     """
790     if self.__operator_isValid and self.debug() : print "PDE Debug: Operator has to be rebuilt."
791     self.__rebuildSolution(deep)
792     self.__operator_isValid=False
793     if deep: self.__operator=escript.Operator()
794    
795     def __rebuildRightHandSide(self,deep=False):
796     """
797     @brief indicates the right hand side has to be rebuild next time it is used
798     """
799     if self.__righthandside_isValid and self.debug() : print "PDE Debug: Right hand side has to be rebuilt."
800     self.__rebuildSolution(deep)
801     self.__righthandside_isValid=False
802     if deep: self.__righthandside=escript.Data()
803    
804     def __rebuildSystem(self,deep=False):
805     """
806     @brief annonced that all coefficient name has been changed
807     """
808     self.__rebuildSolution(deep)
809     self.__rebuildOperator(deep)
810     self.__rebuildRightHandSide(deep)
811    
812     def __checkMatrixType(self):
813     """
814     @brief reassess the matrix type and, if needed, initiates an operator rebuild
815     """
816     new_matrix_type=self.getDomain().getSystemMatrixTypeId(self.getSolverMethod(),self.isSymmetric())
817     if not new_matrix_type==self.__matrix_type:
818     if self.debug() : print "PDE Debug: Matrix type is now %d."%new_matrix_type
819     self.__matrix_type=new_matrix_type
820     self.__rebuildOperator(deep=True)
821    
822     #============ assembling =======================================================
823 jgs 108 def __copyConstraint(self):
824 jgs 102 """
825     @brief copies the constrint condition into u
826     """
827 jgs 108 if not self.__righthandside.isEmpty():
828     q=self.getCoefficientOfPDE("q")
829     r=self.getCoefficientOfPDE("r")
830     if not q.isEmpty():
831     if r.isEmpty():
832     r2=escript.Data(0,self.__righthandside.getShape(),self.__righthandside.getFunctionSpace())
833     else:
834     r2=escript.Data(r,self.__righthandside.getFunctionSpace())
835     self.__righthandside.copyWithMask(r2,escript.Data(q,self.__righthandside.getFunctionSpace()))
836 jgs 102
837 jgs 108 def __applyConstraint(self):
838 jgs 102 """
839 jgs 108 @brief applies the constraints defined by q and r to the system
840 jgs 102 """
841 jgs 108 q=self.getCoefficientOfPDE("q")
842     r=self.getCoefficientOfPDE("r")
843 jgs 102 if not q.isEmpty() and not self.__operator.isEmpty():
844     # q is the row and column mask to indicate where constraints are set:
845     row_q=escript.Data(q,self.getFunctionSpaceForEquation())
846     col_q=escript.Data(q,self.getFunctionSpaceForSolution())
847     u=self.__makeNewSolution()
848     if r.isEmpty():
849     r_s=self.__makeNewSolution()
850     else:
851     r_s=escript.Data(r,self.getFunctionSpaceForSolution())
852     u.copyWithMask(r_s,col_q)
853 jgs 108 if self.isUsingLumping():
854     self.__operator.copyWithMask(escript.Data(1,q.getShape(),self.getFunctionSpaceForEquation()),row_q)
855     else:
856     if not self.__righthandside.isEmpty(): self.__righthandside-=self.__operator*u
857     self.__operator.nullifyRowsAndCols(row_q,col_q,1.)
858 jgs 102
859 jgs 108 def getSystem(self):
860     """
861     @brief return the operator and right hand side of the PDE
862     """
863     if not self.__operator_isValid or not self.__righthandside_isValid:
864     if self.isUsingLumping():
865     if not self.__operator_isValid:
866     if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():
867     raise TypeError,"Lumped matrix requires same order for equations and unknowns"
868     if not self.getCoefficientOfPDE("A").isEmpty():
869     raise Warning,"Lumped matrix does not allow coefficient A"
870     if not self.getCoefficientOfPDE("B").isEmpty():
871     raise Warning,"Lumped matrix does not allow coefficient B"
872     if not self.getCoefficientOfPDE("C").isEmpty():
873     raise Warning,"Lumped matrix does not allow coefficient C"
874     if self.debug() : print "PDE Debug: New lumped operator is built."
875     mat=self.__makeNewOperator()
876     self.getDomain().addPDEToSystem(mat,escript.Data(), \
877     self.getCoefficientOfPDE("A"), \
878     self.getCoefficientOfPDE("B"), \
879     self.getCoefficientOfPDE("C"), \
880     self.getCoefficientOfPDE("D"), \
881     escript.Data(), \
882     escript.Data(), \
883     self.getCoefficientOfPDE("d"), \
884     escript.Data(),\
885     self.getCoefficientOfPDE("d_contact"), \
886     escript.Data())
887     self.__operator=mat*escript.Data(1,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)
888     self.__applyConstraint()
889     self.__operator_isValid=True
890     if not self.__righthandside_isValid:
891     if self.debug() : print "PDE Debug: New right hand side is built."
892     self.getDomain().addPDEToRHS(self.__getFreshRightHandSide(), \
893     self.getCoefficientOfPDE("X"), \
894     self.getCoefficientOfPDE("Y"),\
895     self.getCoefficientOfPDE("y"),\
896     self.getCoefficientOfPDE("y_contact"))
897     self.__copyConstraint()
898     self.__righthandside_isValid=True
899     else:
900     if not self.__operator_isValid and not self.__righthandside_isValid:
901     if self.debug() : print "PDE Debug: New system is built."
902     self.getDomain().addPDEToSystem(self.__getFreshOperator(),self.__getFreshRightHandSide(), \
903     self.getCoefficientOfPDE("A"), \
904     self.getCoefficientOfPDE("B"), \
905     self.getCoefficientOfPDE("C"), \
906     self.getCoefficientOfPDE("D"), \
907     self.getCoefficientOfPDE("X"), \
908     self.getCoefficientOfPDE("Y"), \
909     self.getCoefficientOfPDE("d"), \
910     self.getCoefficientOfPDE("y"), \
911     self.getCoefficientOfPDE("d_contact"), \
912     self.getCoefficientOfPDE("y_contact"))
913     self.__applyConstraint()
914     self.__copyConstraint()
915     self.__operator_isValid=True
916     self.__righthandside_isValid=True
917     elif not self.__righthandside_isValid:
918     if self.debug() : print "PDE Debug: New right hand side is built."
919     self.getDomain().addPDEToRHS(self.__getFreshRightHandSide(), \
920     self.getCoefficientOfPDE("X"), \
921     self.getCoefficientOfPDE("Y"),\
922     self.getCoefficientOfPDE("y"),\
923     self.getCoefficientOfPDE("y_contact"))
924     self.__copyConstraint()
925     self.__righthandside_isValid=True
926     elif not self.__operator_isValid:
927     if self.debug() : print "PDE Debug: New operator is built."
928     self.getDomain().addPDEToSystem(self.__getFreshOperator(),escript.Data(), \
929     self.getCoefficientOfPDE("A"), \
930     self.getCoefficientOfPDE("B"), \
931     self.getCoefficientOfPDE("C"), \
932     self.getCoefficientOfPDE("D"), \
933     escript.Data(), \
934     escript.Data(), \
935     self.getCoefficientOfPDE("d"), \
936     escript.Data(),\
937     self.getCoefficientOfPDE("d_contact"), \
938     escript.Data())
939     self.__applyConstraint()
940     self.__operator_isValid=True
941     return (self.__operator,self.__righthandside)
942 jgs 102 def getOperator(self):
943     """
944     @brief returns the operator of the PDE
945     """
946 jgs 108 return self.getSystem()[0]
947 jgs 102
948 jgs 108 def getRightHandSide(self):
949 jgs 102 """
950     @brief returns the right hand side of the PDE
951     """
952 jgs 108 return self.getSystem()[1]
953 jgs 102
954     def solve(self,**options):
955     """
956     @brief solve the PDE
957    
958     @param options
959     """
960     mat,f=self.getSystem()
961     if self.isUsingLumping():
962     out=f/mat
963     else:
964     options[util.TOLERANCE_KEY]=self.getTolerance()
965     options[util.METHOD_KEY]=self.getSolverMethod()
966     options[util.SYMMETRY_KEY]=self.isSymmetric()
967     if self.debug() : print "PDE Debug: solver options: ",options
968     out=mat.solve(f,options)
969     return out
970    
971     def getSolution(self,**options):
972     """
973     @brief returns the solution of the PDE
974    
975     @param options
976     """
977     if not self.__solution_isValid:
978     if self.debug() : print "PDE Debug: PDE is resolved."
979     self.__solution=self.solve(**options)
980     self.__solution_isValid=True
981     return self.__solution
982    
983 jgs 110
984    
985     def ELMAN_RAMAGE(P): return (P-1.).wherePositive()*0.5*(1.-1./(P+1.e-15))
986     def SIMPLIFIED_BROOK_HUGHES(P):
987     c=(P-3.).whereNegative()
988     return P/6.*c+1./2.*(1.-c)
989     def HALF(P): return escript.Scalar(0.5,P.getFunctionSpace())
990    
991    
992 jgs 108 class AdvectivePDE(LinearPDE):
993     """
994     @brief Class to handel a linear PDE domineated by advective terms:
995    
996     class to define a linear PDE of the form
997 jgs 104
998 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
999 jgs 102
1000 jgs 108 with boundary conditons:
1001 jgs 102
1002 jgs 108 n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_{ik}u_k = - n_j*X_{ij} + y_i
1003 jgs 102
1004 jgs 108 and contact conditions
1005 jgs 102
1006 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
1007 jgs 102
1008 jgs 108 and constraints:
1009 jgs 102
1010 jgs 108 u_i=r_i where q_i>0
1011 jgs 102
1012 jgs 110 """
1013     def __init__(self,domain,numEquations=0,numSolutions=0,xi=ELMAN_RAMAGE):
1014     LinearPDE.__init__(self,domain,numEquations,numSolutions)
1015     self.__xi=xi
1016     self.__Xi=escript.Data()
1017 jgs 102
1018 jgs 110 def __calculateXi(self,peclet_factor,Z,h):
1019     Z_max=util.Lsup(Z)
1020     if Z_max>0.:
1021     return h*self.__xi(Z*peclet_factor)/(Z+Z_max*self.TOL)
1022     else:
1023     return 0.
1024 jgs 102
1025 jgs 110 def setValue(self,**args):
1026     if "A" in args.keys() or "B" in args.keys() or "C" in args.keys(): self.__Xi=escript.Data()
1027     self._setValue(**args)
1028    
1029     def getXi(self):
1030     if self.__Xi.isEmpty():
1031     B=self.getCoefficient("B")
1032     C=self.getCoefficient("C")
1033     A=self.getCoefficient("A")
1034     h=self.getDomain().getSize()
1035     self.__Xi=escript.Scalar(0.,self.getFunctionSpaceForCoefficient("A"))
1036     if not C.isEmpty() or not B.isEmpty():
1037     if not C.isEmpty() and not B.isEmpty():
1038     Z2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))
1039     if self.getNumEquations()>1:
1040     if self.getNumSolutions()>1:
1041     for i in range(self.getNumEquations()):
1042     for k in range(self.getNumSolutions()):
1043     for l in range(self.getDim()): Z2+=(C[i,k,l]-B[i,l,k])**2
1044     else:
1045     for i in range(self.getNumEquations()):
1046     for l in range(self.getDim()): Z2+=(C[i,l]-B[i,l])**2
1047     else:
1048     if self.getNumSolutions()>1:
1049     for k in range(self.getNumSolutions()):
1050     for l in range(self.getDim()): Z2+=(C[k,l]-B[l,k])**2
1051     else:
1052     for l in range(self.getDim()): Z2+=(C[l]-B[l])**2
1053     length_of_Z=util.sqrt(Z2)
1054     elif C.isEmpty():
1055     length_of_Z=util.length(B)
1056     else:
1057     length_of_Z=util.length(C)
1058 jgs 102
1059 jgs 110 Z_max=util.Lsup(length_of_Z)
1060     if Z_max>0.:
1061     length_of_A=util.length(A)
1062     A_max=util.Lsup(length_of_A)
1063     if A_max>0:
1064     inv_A=1./(length_of_A+A_max*self.TOL)
1065     else:
1066     inv_A=1./self.TOL
1067     peclet_number=length_of_Z*h/2*inv_A
1068     xi=self.__xi(peclet_number)
1069     self.__Xi=h*xi/(length_of_Z+Z_max*self.TOL)
1070     print "@ preclet number = %e"%util.Lsup(peclet_number),util.Lsup(xi),util.Lsup(length_of_Z)
1071     return self.__Xi
1072    
1073 jgs 102
1074 jgs 108 def getCoefficientOfPDE(self,name):
1075     """
1076     @brief return the value of the coefficient name of the general PDE
1077     @param name
1078     """
1079 jgs 110 if not self.getNumEquations() == self.getNumSolutions():
1080     raise ValueError,"AdvectivePDE expects the number of solution componets and the number of equations to be equal."
1081    
1082 jgs 108 if name == "A" :
1083     A=self.getCoefficient("A")
1084     B=self.getCoefficient("B")
1085     C=self.getCoefficient("C")
1086 jgs 110 if B.isEmpty() and C.isEmpty():
1087     Aout=A
1088     else:
1089     if A.isEmpty():
1090     Aout=self.createNewCoefficient("A")
1091     else:
1092     Aout=A[:]
1093     Xi=self.getXi()
1094     if self.getNumEquations()>1:
1095     for i in range(self.getNumEquations()):
1096     for j in range(self.getDim()):
1097 jgs 108 for k in range(self.getNumSolutions()):
1098 jgs 110 for l in range(self.getDim()):
1099     if not C.isEmpty() and not B.isEmpty():
1100     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])
1101     elif C.isEmpty():
1102     for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*B[p,j,i]*B[p,l,k]
1103     else:
1104     for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*C[p,i,j]*C[p,k,l]
1105     else:
1106     for j in range(self.getDim()):
1107     for l in range(self.getDim()):
1108     if not C.isEmpty() and not B.isEmpty():
1109     Aout[j,l]+=Xi*(C[j]-B[j])*(C[l]-B[l])
1110     elif C.isEmpty():
1111     Aout[j,l]+=Xi*B[j]*B[l]
1112     else:
1113     Aout[j,l]+=Xi*C[j]*C[l]
1114     return Aout
1115 jgs 108 elif name == "B" :
1116 jgs 110 B=self.getCoefficient("B")
1117     C=self.getCoefficient("C")
1118     D=self.getCoefficient("D")
1119     if C.isEmpty() or D.isEmpty():
1120     Bout=B
1121     else:
1122     Xi=self.getXi()
1123     if B.isEmpty():
1124     Bout=self.createNewCoefficient("B")
1125     else:
1126     Bout=B[:]
1127     if self.getNumEquations()>1:
1128     for k in range(self.getNumSolutions()):
1129     for p in range(self.getNumEquations()):
1130     tmp=Xi*D[p,k]
1131     for i in range(self.getNumEquations()):
1132     for j in range(self.getDim()):
1133     Bout[i,j,k]+=tmp*C[p,i,j]
1134     else:
1135     tmp=Xi*D
1136     for j in range(self.getDim()): Bout[j]+=tmp*C[j]
1137     return Bout
1138 jgs 108 elif name == "C" :
1139 jgs 110 B=self.getCoefficient("B")
1140     C=self.getCoefficient("C")
1141     D=self.getCoefficient("D")
1142     if B.isEmpty() or D.isEmpty():
1143     Cout=C
1144     else:
1145     Xi=self.getXi()
1146     if C.isEmpty():
1147     Cout=self.createNewCoefficient("C")
1148     else:
1149     Cout=C[:]
1150     if self.getNumEquations()>1:
1151     for k in range(self.getNumSolutions()):
1152     for p in range(self.getNumEquations()):
1153     tmp=Xi*D[p,k]
1154     for i in range(self.getNumEquations()):
1155     for l in range(self.getDim()):
1156     Cout[i,k,l]+=tmp*B[p,l,i]
1157     else:
1158     tmp=Xi*D
1159     for j in range(self.getDim()): Cout[j]+=tmp*B[j]
1160     return Cout
1161 jgs 108 elif name == "D" :
1162     return self.getCoefficient("D")
1163     elif name == "X" :
1164 jgs 110 X=self.getCoefficient("X")
1165     Y=self.getCoefficient("Y")
1166     B=self.getCoefficient("B")
1167     C=self.getCoefficient("C")
1168     if Y.isEmpty() or (B.isEmpty() and C.isEmpty()):
1169     Xout=X
1170     else:
1171     if X.isEmpty():
1172     Xout=self.createNewCoefficient("X")
1173     else:
1174     Xout=X[:]
1175     Xi=self.getXi()
1176     if self.getNumEquations()>1:
1177     for p in range(self.getNumEquations()):
1178     tmp=Xi*Y[p]
1179     for i in range(self.getNumEquations()):
1180     for j in range(self.getDim()):
1181     if not C.isEmpty() and not B.isEmpty():
1182     Xout[i,j]+=tmp*(C[p,i,j]-B[p,j,i])
1183     elif C.isEmpty():
1184     Xout[i,j]-=tmp*B[p,j,i]
1185     else:
1186     Xout[i,j]+=tmp*C[p,i,j]
1187     else:
1188     tmp=Xi*Y
1189     for j in range(self.getDim()):
1190     if not C.isEmpty() and not B.isEmpty():
1191     Xout[j]+=tmp*(C[j]-B[j])
1192     elif C.isEmpty():
1193     Xout[j]-=tmp*B[j]
1194     else:
1195     Xout[j]+=tmp*C[j]
1196     return Xout
1197 jgs 108 elif name == "Y" :
1198     return self.getCoefficient("Y")
1199     elif name == "d" :
1200     return self.getCoefficient("d")
1201     elif name == "y" :
1202     return self.getCoefficient("y")
1203     elif name == "d_contact" :
1204     return self.getCoefficient("d_contact")
1205     elif name == "y_contact" :
1206     return self.getCoefficient("y_contact")
1207     elif name == "r" :
1208     return self.getCoefficient("r")
1209     elif name == "q" :
1210     return self.getCoefficient("q")
1211     else:
1212     raise SystemError,"unknown PDE coefficient %s",name
1213    
1214    
1215 jgs 102 class Poisson(LinearPDE):
1216     """
1217     @brief Class to define a Poisson equstion problem:
1218    
1219     class to define a linear PDE of the form
1220    
1221     -u_{,jj} = f
1222    
1223     with boundary conditons:
1224    
1225     n_j*u_{,j} = 0
1226    
1227     and constraints:
1228    
1229     u=0 where q>0
1230    
1231     """
1232    
1233 jgs 108 def __init__(self,domain,f=escript.Data(),q=escript.Data()):
1234     LinearPDE.__init__(self,domain,1,1)
1235     self.COEFFICIENTS={
1236     "f" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1237     "q" : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.BOTH)}
1238 jgs 102 self.setSymmetryOn()
1239     self.setValue(f,q)
1240    
1241     def setValue(self,f=escript.Data(),q=escript.Data()):
1242 jgs 108 self._setValue(f=f,q=q)
1243    
1244     def getCoefficientOfPDE(self,name):
1245     """
1246     @brief return the value of the coefficient name of the general PDE
1247     @param name
1248     """
1249     if name == "A" :
1250     return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))
1251     elif name == "B" :
1252     return escript.Data()
1253     elif name == "C" :
1254     return escript.Data()
1255     elif name == "D" :
1256     return escript.Data()
1257     elif name == "X" :
1258     return escript.Data()
1259     elif name == "Y" :
1260     return self.getCoefficient("f")
1261     elif name == "d" :
1262     return escript.Data()
1263     elif name == "y" :
1264     return escript.Data()
1265     elif name == "d_contact" :
1266     return escript.Data()
1267     elif name == "y_contact" :
1268     return escript.Data()
1269     elif name == "r" :
1270     return escript.Data()
1271     elif name == "q" :
1272     return self.getCoefficient("q")
1273     else:
1274     raise SystemError,"unknown PDE coefficient %s",name

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26