/[escript]/branches/symbolic_from_3470/escript/py_src/linearPDEs.py
ViewVC logotype

Annotation of /branches/symbolic_from_3470/escript/py_src/linearPDEs.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 104 - (hide annotations)
Fri Dec 17 07:43:12 2004 UTC (14 years, 11 months ago) by jgs
Original Path: trunk/esys2/escript/py_src/linearPDEs.py
File MIME type: text/x-python
File size: 36737 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     def identifyDomain(domain=None,data={}):
14     """
15     @brief Return the Domain which is equal to the input domain (if not None)
16     and is the domain of all Data objects in the dictionary data.
17     An exception is raised if this is not possible
18    
19     @param domain
20     @param data
21     """
22     # get the domain used by any Data object in the list data:
23     data_domain=None
24     for d in data.itervalues():
25     if isinstance(d,escript.Data):
26     if not d.isEmpty(): data_domain=d.getDomain()
27     # check if domain and data_domain are identical?
28     if domain == None:
29     if data_domain == None:
30     raise ValueError,"Undefined PDE domain. Specify a domain or use a Data class object as coefficient"
31     else:
32     if data_domain == None:
33     data_domain=domain
34     else:
35     if not data_domain == domain:
36     raise ValueError,"Domain of coefficients doesnot match specified domain"
37     # now we check if all Data class object coefficients are defined on data_domain:
38     for i,d in data.iteritems():
39     if isinstance(d,escript.Data):
40     if not d.isEmpty():
41     if not data_domain==d.getDomain():
42     raise ValueError,"Illegal domain for coefficient %s."%i
43     # done:
44     return data_domain
45    
46     def identifyNumEquationsAndSolutions(dim,coef={}):
47     # get number of equations and number of unknowns:
48     numEquations=0
49     numSolutions=0
50     for i in coef.iterkeys():
51     if not coef[i].isEmpty():
52     res=_PDECoefficientTypes[i].estimateNumEquationsAndNumSolutions(coef[i].getShape(),dim)
53     if res==None:
54     raise ValueError,"Illegal shape %s of coefficient %s"%(coef[i].getShape().__str__(),i)
55     else:
56     numEquations=max(numEquations,res[0])
57     numSolutions=max(numSolutions,res[1])
58     return numEquations,numSolutions
59    
60    
61     def _CompTuple2(t1,t2):
62     """
63     @brief
64    
65     @param t1
66     @param t2
67     """
68     dif=t1[0]+t1[1]-(t2[0]+t2[1])
69     if dif<0: return 1
70     elif dif>0: return -1
71     else: return 0
72    
73     class PDECoefficientType:
74     """
75     @brief
76     """
77     # identifier for location of Data objects defining coefficients
78     INTERIOR=0
79     BOUNDARY=1
80     CONTACT=2
81     CONTINUOUS=3
82     # identifier in the pattern of coefficients:
83     # the pattern is a tuple of EQUATION,SOLUTION,DIM where DIM represents the spatial dimension, EQUATION the number of equations and SOLUTION the
84     # number of unknowns.
85     EQUATION=3
86     SOLUTION=4
87     DIM=5
88     # indicator for what is altered if the coefficient is altered:
89     OPERATOR=5
90     RIGHTHANDSIDE=6
91     BOTH=7
92     def __init__(self,where,pattern,altering):
93     """
94     @brief Initialise a PDE Coefficient type
95     """
96     self.what=where
97     self.pattern=pattern
98     self.altering=altering
99    
100     def getFunctionSpace(self,domain):
101     """
102     @brief defines the FunctionSpace of the coefficient on the domain
103    
104     @param domain
105     """
106     if self.what==self.INTERIOR: return escript.Function(domain)
107     elif self.what==self.BOUNDARY: return escript.FunctionOnBoundary(domain)
108     elif self.what==self.CONTACT: return escript.FunctionOnContactZero(domain)
109     elif self.what==self.CONTINUOUS: return escript.ContinuousFunction(domain)
110    
111     def isAlteringOperator(self):
112     """
113     @brief return true if the operator of the PDE is changed when the coefficient is changed
114     """
115     if self.altering==self.OPERATOR or self.altering==self.BOTH:
116     return not None
117     else:
118     return None
119    
120     def isAlteringRightHandSide(self):
121     """
122     @brief return true if the right hand side of the PDE is changed when the coefficient is changed
123     """
124     if self.altering==self.RIGHTHANDSIDE or self.altering==self.BOTH:
125     return not None
126     else:
127     return None
128    
129     def estimateNumEquationsAndNumSolutions(self,shape=(),dim=3):
130     """
131     @brief tries to estimate the number of equations in a given tensor shape for a given spatial dimension dim
132    
133     @param shape
134     @param dim
135     """
136     if len(shape)>0:
137     num=max(shape)+1
138     else:
139     num=1
140     search=[]
141     for u in range(num):
142     for e in range(num):
143     search.append((e,u))
144     search.sort(_CompTuple2)
145     for item in search:
146     s=self.buildShape(item[0],item[1],dim)
147     if len(s)==0 and len(shape)==0:
148     return (1,1)
149     else:
150     if s==shape: return item
151     return None
152    
153     def buildShape(self,e=1,u=1,dim=3):
154     """
155     @brief builds the required shape for a given number of equations e, number of unknowns u and spatial dimension dim
156    
157     @param e
158     @param u
159     @param dim
160     """
161     s=()
162     for i in self.pattern:
163     if i==self.EQUATION:
164     if e>1: s=s+(e,)
165     elif i==self.SOLUTION:
166     if u>1: s=s+(u,)
167     else:
168     s=s+(dim,)
169     return s
170    
171     _PDECoefficientTypes={
172     "A" : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,PDECoefficientType.DIM,PDECoefficientType.SOLUTION,PDECoefficientType.DIM),PDECoefficientType.OPERATOR),
173     "B" : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,PDECoefficientType.DIM,PDECoefficientType.SOLUTION),PDECoefficientType.OPERATOR),
174     "C" : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,PDECoefficientType.SOLUTION,PDECoefficientType.DIM),PDECoefficientType.OPERATOR),
175     "D" : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,PDECoefficientType.SOLUTION),PDECoefficientType.OPERATOR),
176     "X" : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,PDECoefficientType.DIM),PDECoefficientType.RIGHTHANDSIDE),
177     "Y" : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,),PDECoefficientType.RIGHTHANDSIDE),
178     "d" : PDECoefficientType(PDECoefficientType.BOUNDARY,(PDECoefficientType.EQUATION,PDECoefficientType.SOLUTION),PDECoefficientType.OPERATOR),
179     "y" : PDECoefficientType(PDECoefficientType.BOUNDARY,(PDECoefficientType.EQUATION,),PDECoefficientType.RIGHTHANDSIDE),
180     "d_contact" : PDECoefficientType(PDECoefficientType.CONTACT,(PDECoefficientType.EQUATION,PDECoefficientType.SOLUTION),PDECoefficientType.OPERATOR),
181     "y_contact" : PDECoefficientType(PDECoefficientType.CONTACT,(PDECoefficientType.EQUATION,),PDECoefficientType.RIGHTHANDSIDE),
182     "r" : PDECoefficientType(PDECoefficientType.CONTINUOUS,(PDECoefficientType.EQUATION,),PDECoefficientType.RIGHTHANDSIDE),
183     "q" : PDECoefficientType(PDECoefficientType.CONTINUOUS,(PDECoefficientType.SOLUTION,),PDECoefficientType.BOTH),
184     }
185    
186     class LinearPDE:
187     """
188     @brief Class to define a linear PDE
189    
190     class to define a linear PDE of the form
191    
192     -(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
193    
194     with boundary conditons:
195    
196     n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_{ik}u_k = - n_j*X_{ij} + y_i
197    
198     and contact conditions
199    
200     n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_contact_{ik}[u_k] = - n_j*X_{ij} + y_contact_i
201    
202     and constraints:
203    
204     u_i=r_i where q_i>0
205    
206     """
207     DEFAULT_METHOD=util.DEFAULT_METHOD
208     DIRECT=util.DIRECT
209     CHOLEVSKY=util.CHOLEVSKY
210     PCG=util.PCG
211     CR=util.CR
212     CGS=util.CGS
213     BICGSTAB=util.BICGSTAB
214     SSOR=util.SSOR
215     GMRES=util.GMRES
216     PRES20=util.PRES20
217    
218 jgs 104 def __init__(self,domain,numEquations=0,numSolutions=0):
219 jgs 102 """
220     @brief initializes a new linear PDE.
221    
222     @param args
223     """
224    
225     # initialize attributes
226     self.__debug=None
227 jgs 104 self.__domain=domain
228     self.__numEquations=numEquations
229     self.__numSolutions=numSolutions
230 jgs 102 self.cleanCoefficients()
231    
232     self.__operator=escript.Operator()
233     self.__operator_isValid=False
234     self.__righthandside=escript.Data()
235     self.__righthandside_isValid=False
236     self.__solution=escript.Data()
237     self.__solution_isValid=False
238    
239     # set some default values:
240     self.__homogeneous_constraint=True
241     self.__row_function_space=escript.Solution(self.__domain)
242     self.__column_function_space=escript.Solution(self.__domain)
243     self.__tolerance=1.e-8
244     self.__solver_method=util.DEFAULT_METHOD
245     self.__matrix_type=self.__domain.getSystemMatrixTypeId(util.DEFAULT_METHOD,False)
246     self.__sym=False
247     self.__lumping=False
248    
249     def getCoefficient(self,name):
250     """
251     @brief return the value of the coefficient name
252    
253     @param name
254     """
255     return self.__coefficient[name]
256    
257     def setValue(self,**coefficients):
258     """
259     @brief sets new values to coefficients
260    
261     @param coefficients
262     """
263 jgs 104 self._setValue(**coefficients)
264 jgs 102
265    
266 jgs 104 def _setValue(self,**coefficients):
267 jgs 102 """
268     @brief sets new values to coefficients
269    
270     @param coefficients
271     """
272    
273     # get the dictionary of the coefficinets been altered:
274     alteredCoefficients={}
275     for i,d in coefficients.iteritems():
276     if self.hasCoefficient(i):
277     if d == None:
278     alteredCoefficients[i]=escript.Data()
279     elif isinstance(d,escript.Data):
280     if d.isEmpty():
281     alteredCoefficients[i]=escript.Data()
282     else:
283     alteredCoefficients[i]=escript.Data(d,self.getFunctionSpaceOfCoefficient(i))
284     else:
285 jgs 104 alteredCoefficients[i]=escript.Data(d,self.getFunctionSpaceOfCoefficient(i))
286 jgs 102 else:
287     raise ValueError,"Attempt to set undefined coefficient %s"%i
288     # if numEquations and numSolutions is undefined we try identify their values based on the coefficients:
289     if self.__numEquations<1 or self.__numSolutions<1:
290     numEquations,numSolutions=identifyNumEquationsAndSolutions(self.getDomain().getDim(),alteredCoefficients)
291     if self.__numEquations<1 and numEquations>0: self.__numEquations=numEquations
292     if self.__numSolutions<1 and numSolutions>0: self.__numSolutions=numSolutions
293     if self.debug() and self.__numEquations>0: print "PDE Debug: identified number of equations is ",self.__numEquations
294     if self.debug() and self.__numSolutions>0: print "PDE Debug: identified number of solutions is ",self.__numSolutions
295    
296     # now we check the shape of the coefficient if numEquations and numSolutions are set:
297     if self.__numEquations>0 and self.__numSolutions>0:
298     for i in self.__coefficient.iterkeys():
299     if alteredCoefficients.has_key(i) and not alteredCoefficients[i].isEmpty():
300     if not self.getShapeOfCoefficient(i)==alteredCoefficients[i].getShape():
301     raise ValueError,"Expected shape for coefficient %s is %s but actual shape is %s."%(i,self.getShapeOfCoefficient(i),alteredCoefficients[i].getShape())
302     else:
303     if not self.__coefficient[i].isEmpty():
304     if not self.getShapeOfCoefficient(i)==self.__coefficient[i].getShape():
305     raise ValueError,"Expected shape for coefficient %s is %s but actual shape is %s."%(i,self.getShapeOfCoefficient(i),self.__coefficient[i].getShape())
306     # overwrite new values:
307     for i,d in alteredCoefficients.iteritems():
308     if self.debug(): print "PDE Debug: Coefficient %s has been altered."%i
309     self.__coefficient[i]=d
310     self.alteredCoefficient(i)
311    
312     # reset the HomogeneousConstraintFlag:
313     self.__setHomogeneousConstraintFlag()
314     if not "q" in alteredCoefficients and not self.__homogeneous_constraint: self.__rebuildSystem()
315    
316     def cleanCoefficients(self):
317     """
318     @brief resets all coefficients to default values.
319     """
320     self.__coefficient={}
321     for i in _PDECoefficientTypes.iterkeys():
322     self.__coefficient[i]=escript.Data()
323    
324     def getShapeOfCoefficient(self,name):
325     """
326     @brief return the shape of the coefficient name
327    
328     @param name
329     """
330     if self.hasCoefficient(name):
331     return _PDECoefficientTypes[name].buildShape(self.getNumEquations(),self.getNumSolutions(),self.getDomain().getDim())
332     else:
333     raise ValueError,"Solution coefficient %s requested"%name
334    
335     def getFunctionSpaceOfCoefficient(self,name):
336     """
337     @brief return the atoms of the coefficient name
338    
339     @param name
340     """
341     if self.hasCoefficient(name):
342     return _PDECoefficientTypes[name].getFunctionSpace(self.getDomain())
343     else:
344     raise ValueError,"Solution coefficient %s requested"%name
345    
346     def alteredCoefficient(self,name):
347     """
348     @brief annonced that coefficient name has been changed
349    
350     @param name
351     """
352     if self.hasCoefficient(name):
353     if _PDECoefficientTypes[name].isAlteringOperator(): self.__rebuildOperator()
354     if _PDECoefficientTypes[name].isAlteringRightHandSide(): self.__rebuildRightHandSide()
355     else:
356     raise ValueError,"Solution coefficient %s requested"%name
357    
358     def __setHomogeneousConstraintFlag(self):
359     """
360     @brief checks if the constraints are homogeneous and sets self.__homogeneous_constraint accordingly.
361     """
362     self.__homogeneous_constraint=True
363     q=self.getCoefficient("q")
364     r=self.getCoefficient("r")
365     if not q.isEmpty() and not r.isEmpty():
366     print (q*r).Lsup(), 1.e-13*r.Lsup()
367     if (q*r).Lsup()>=1.e-13*r.Lsup(): self.__homogeneous_constraint=False
368     if self.debug():
369     if self.__homogeneous_constraint:
370     print "PDE Debug: Constraints are homogeneous."
371     else:
372     print "PDE Debug: Constraints are inhomogeneous."
373    
374    
375     def hasCoefficient(self,name):
376     """
377     @brief return true if name is the name of a coefficient
378    
379     @param name
380     """
381     return self.__coefficient.has_key(name)
382    
383     def getFunctionSpaceForEquation(self):
384     """
385     @brief return true if the test functions should use reduced order
386     """
387     return self.__row_function_space
388    
389     def getFunctionSpaceForSolution(self):
390     """
391     @brief return true if the interpolation of the solution should use reduced order
392     """
393     return self.__column_function_space
394    
395     # ===== debug ==============================================================
396     def setDebugOn(self):
397     """
398     @brief
399     """
400     self.__debug=not None
401    
402     def setDebugOff(self):
403     """
404     @brief
405     """
406     self.__debug=None
407    
408     def debug(self):
409     """
410     @brief returns true if the PDE is in the debug mode
411     """
412     return self.__debug
413    
414     #===== Lumping ===========================
415     def setLumpingOn(self):
416     """
417     @brief indicates to use matrix lumping
418     """
419     if not self.isUsingLumping():
420     raise SystemError,"Lumping is not working yet! Talk to the experts"
421     if self.debug() : print "PDE Debug: lumping is set on"
422     self.__rebuildOperator()
423     self.__lumping=True
424    
425     def setLumpingOff(self):
426     """
427     @brief switches off matrix lumping
428     """
429     if self.isUsingLumping():
430     if self.debug() : print "PDE Debug: lumping is set off"
431     self.__rebuildOperator()
432     self.__lumping=False
433    
434     def setLumping(self,flag=False):
435     """
436     @brief set the matrix lumping flag to flag
437     """
438     if flag:
439     self.setLumpingOn()
440     else:
441     self.setLumpingOff()
442    
443     def isUsingLumping(self):
444     """
445     @brief
446     """
447     return self.__lumping
448    
449     #============ method business =========================================================
450     def setSolverMethod(self,solver=util.DEFAULT_METHOD):
451     """
452     @brief sets a new solver
453     """
454     if not solver==self.getSolverMethod():
455     self.__solver_method=solver
456     if self.debug() : print "PDE Debug: New solver is %s"%solver
457     self.__checkMatrixType()
458    
459     def getSolverMethod(self):
460     """
461     @brief returns the solver method
462     """
463     return self.__solver_method
464    
465     #============ tolerance business =========================================================
466     def setTolerance(self,tol=1.e-8):
467     """
468     @brief resets the tolerance to tol.
469     """
470     if not tol>0:
471     raise ValueException,"Tolerance as to be positive"
472     if tol<self.getTolerance(): self.__rebuildSolution()
473     if self.debug() : print "PDE Debug: New tolerance %e",tol
474     self.__tolerance=tol
475     return
476     def getTolerance(self):
477     """
478     @brief returns the tolerance set for the solution
479     """
480     return self.__tolerance
481    
482     #===== symmetry flag ==========================
483     def isSymmetric(self):
484     """
485     @brief returns true is the operator is considered to be symmetric
486     """
487     return self.__sym
488    
489     def setSymmetryOn(self):
490     """
491     @brief sets the symmetry flag to true
492     """
493     if not self.isSymmetric():
494     if self.debug() : print "PDE Debug: Operator is set to be symmetric"
495     self.__sym=True
496     self.__checkMatrixType()
497    
498     def setSymmetryOff(self):
499     """
500     @brief sets the symmetry flag to false
501     """
502     if self.isSymmetric():
503     if self.debug() : print "PDE Debug: Operator is set to be unsymmetric"
504     self.__sym=False
505     self.__checkMatrixType()
506    
507     def setSymmetryTo(self,flag=False):
508     """
509     @brief sets the symmetry flag to flag
510    
511     @param flag
512     """
513     if flag:
514     self.setSymmetryOn()
515     else:
516     self.setSymmetryOff()
517    
518     #===== order reduction ==========================
519     def setReducedOrderOn(self):
520     """
521     @brief switches to on reduced order
522     """
523     self.setReducedOrderForSolutionOn()
524     self.setReducedOrderForEquationOn()
525    
526     def setReducedOrderOff(self):
527     """
528     @brief switches to full order
529     """
530     self.setReducedOrderForSolutionOff()
531     self.setReducedOrderForEquationOff()
532    
533     def setReducedOrderTo(self,flag=False):
534     """
535     @brief sets order according to flag
536    
537     @param flag
538     """
539     self.setReducedOrderForSolutionTo(flag)
540     self.setReducedOrderForEquationTo(flag)
541    
542    
543     #===== order reduction solution ==========================
544     def setReducedOrderForSolutionOn(self):
545     """
546     @brief switches to reduced order to interpolate solution
547     """
548     new_fs=escript.ReducedSolution(self.getDomain())
549     if self.getFunctionSpaceForSolution()!=new_fs:
550     if self.debug() : print "PDE Debug: Reduced order is used to interpolate solution."
551     self.__column_function_space=new_fs
552     self.__rebuildSystem(deep=True)
553    
554     def setReducedOrderForSolutionOff(self):
555     """
556     @brief switches to full order to interpolate solution
557     """
558     new_fs=escript.Solution(self.getDomain())
559     if self.getFunctionSpaceForSolution()!=new_fs:
560     if self.debug() : print "PDE Debug: Full order is used to interpolate solution."
561     self.__column_function_space=new_fs
562     self.__rebuildSystem(deep=True)
563    
564     def setReducedOrderForSolutionTo(self,flag=False):
565     """
566     @brief sets order for test functions according to flag
567    
568     @param flag
569     """
570     if flag:
571     self.setReducedOrderForSolutionOn()
572     else:
573     self.setReducedOrderForSolutionOff()
574    
575     #===== order reduction equation ==========================
576     def setReducedOrderForEquationOn(self):
577     """
578     @brief switches to reduced order for test functions
579     """
580     new_fs=escript.ReducedSolution(self.getDomain())
581     if self.getFunctionSpaceForEquation()!=new_fs:
582     if self.debug() : print "PDE Debug: Reduced order is used for test functions."
583     self.__row_function_space=new_fs
584     self.__rebuildSystem(deep=True)
585    
586     def setReducedOrderForEquationOff(self):
587     """
588     @brief switches to full order for test functions
589     """
590     new_fs=escript.Solution(self.getDomain())
591     if self.getFunctionSpaceForEquation()!=new_fs:
592     if self.debug() : print "PDE Debug: Full order is used for test functions."
593     self.__row_function_space=new_fs
594     self.__rebuildSystem(deep=True)
595    
596     def setReducedOrderForEquationTo(self,flag=False):
597     """
598     @brief sets order for test functions according to flag
599    
600     @param flag
601     """
602     if flag:
603     self.setReducedOrderForEquationOn()
604     else:
605     self.setReducedOrderForEquationOff()
606    
607    
608     # ==== initialization =====================================================================
609     def __makeNewOperator(self):
610     """
611     @brief
612     """
613     return self.getDomain().newOperator( \
614     self.getNumEquations(), \
615     self.getFunctionSpaceForEquation(), \
616     self.getNumSolutions(), \
617     self.getFunctionSpaceForSolution(), \
618     self.__matrix_type)
619    
620     def __makeNewRightHandSide(self):
621     """
622     @brief
623     """
624     return escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)
625    
626     def __makeNewSolution(self):
627     """
628     @brief
629     """
630     return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)
631    
632     def __getFreshOperator(self):
633     """
634     @brief
635     """
636     if self.__operator.isEmpty():
637     self.__operator=self.__makeNewOperator()
638     if self.debug() : print "PDE Debug: New operator allocated"
639     else:
640     self.__operator.setValue(0.)
641     if self.debug() : print "PDE Debug: Operator reset to zero"
642     return self.__operator
643    
644     def __getFreshRightHandSide(self):
645     """
646     @brief
647     """
648     if self.__righthandside.isEmpty():
649     self.__righthandside=self.__makeNewRightHandSide()
650     if self.debug() : print "PDE Debug: New right hand side allocated"
651     else:
652     print "fix self.__righthandside*=0"
653     self.__righthandside*=0.
654     if self.debug() : print "PDE Debug: Right hand side reset to zero"
655     return self.__righthandside
656    
657     # ==== rebuild switches =====================================================================
658     def __rebuildSolution(self,deep=False):
659     """
660     @brief indicates the PDE has to be reolved if the solution is requested
661     """
662     if self.__solution_isValid and self.debug() : print "PDE Debug: PDE has to be resolved."
663     self.__solution_isValid=False
664     if deep: self.__solution=escript.Data(deep)
665    
666    
667     def __rebuildOperator(self,deep=False):
668     """
669     @brief indicates the operator has to be rebuilt next time it is used
670     """
671     if self.__operator_isValid and self.debug() : print "PDE Debug: Operator has to be rebuilt."
672     self.__rebuildSolution(deep)
673     self.__operator_isValid=False
674     if deep: self.__operator=escript.Operator()
675    
676     def __rebuildRightHandSide(self,deep=False):
677     """
678     @brief indicates the right hand side has to be rebuild next time it is used
679     """
680     if self.__righthandside_isValid and self.debug() : print "PDE Debug: Right hand side has to be rebuilt."
681     self.__rebuildSolution(deep)
682     self.__righthandside_isValid=False
683     if not self.__homogeneous_constraint: self.__rebuildOperator()
684     if deep: self.__righthandside=escript.Data()
685    
686     def __rebuildSystem(self,deep=False):
687     """
688     @brief annonced that all coefficient name has been changed
689     """
690     self.__rebuildSolution(deep)
691     self.__rebuildOperator(deep)
692     self.__rebuildRightHandSide(deep)
693    
694     def __checkMatrixType(self):
695     """
696     @brief reassess the matrix type and, if needed, initiates an operator rebuild
697     """
698     new_matrix_type=self.getDomain().getSystemMatrixTypeId(self.getSolverMethod(),self.isSymmetric())
699     if not new_matrix_type==self.__matrix_type:
700     if self.debug() : print "PDE Debug: Matrix type is now %d."%new_matrix_type
701     self.__matrix_type=new_matrix_type
702     self.__rebuildOperator(deep=True)
703    
704     #============ assembling =======================================================
705     def __copyConstraint(self,u):
706     """
707     @brief copies the constrint condition into u
708     """
709     q=self.getCoefficient("q")
710     r=self.getCoefficient("r")
711     if not q.isEmpty():
712     if r.isEmpty():
713     r2=escript.Data(0,u.getShape(),u.getFunctionSpace())
714     else:
715     r2=escript.Data(r,u.getFunctionSpace())
716     u.copyWithMask(r2,escript.Data(q,u.getFunctionSpace()))
717    
718     def __applyConstraint(self,rhs_update=True):
719     """
720     @brief applies the constraints defined by q and r to the system
721     """
722     q=self.getCoefficient("q")
723     r=self.getCoefficient("r")
724     if not q.isEmpty() and not self.__operator.isEmpty():
725     # q is the row and column mask to indicate where constraints are set:
726     row_q=escript.Data(q,self.getFunctionSpaceForEquation())
727     col_q=escript.Data(q,self.getFunctionSpaceForSolution())
728     u=self.__makeNewSolution()
729     if r.isEmpty():
730     r_s=self.__makeNewSolution()
731     else:
732     r_s=escript.Data(r,self.getFunctionSpaceForSolution())
733     u.copyWithMask(r_s,col_q)
734     if not self.__righthandside.isEmpty() and rhs_update: self.__righthandside-=self.__operator*u
735     self.__operator.nullifyRowsAndCols(row_q,col_q,1.)
736    
737     def getOperator(self):
738     """
739     @brief returns the operator of the PDE
740     """
741     if not self.__operator_isValid:
742     # some Constraints are applying for a lumpled stifness matrix:
743     if self.isUsingLumping():
744     if self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():
745     raise TypeError,"Lumped matrix requires same order for equations and unknowns"
746     if not self.getCoefficient("A").isEmpty():
747     raise Warning,"Lumped matrix does not allow coefficient A"
748     if not self.getCoefficient("B").isEmpty():
749     raise Warning,"Lumped matrix does not allow coefficient B"
750     if not self.getCoefficient("C").isEmpty():
751     raise Warning,"Lumped matrix does not allow coefficient C"
752     if not self.getCoefficient("D").isEmpty():
753     raise Warning,"Lumped matrix does not allow coefficient D"
754     if self.debug() : print "PDE Debug: New lumped operator is built."
755     mat=self.__makeNewOperator(self)
756     else:
757     if self.debug() : print "PDE Debug: New operator is built."
758     mat=self.__getFreshOperator()
759    
760     self.getDomain().addPDEToSystem(mat,escript.Data(), \
761     self.getCoefficient("A"), \
762     self.getCoefficient("B"), \
763     self.getCoefficient("C"), \
764     self.getCoefficient("D"), \
765     escript.Data(), \
766     escript.Data(), \
767     self.getCoefficient("d"), \
768     escript.Data(),\
769     self.getCoefficient("d_contact"), \
770     escript.Data())
771     if self.isUsingLumping():
772     self.__operator=mat*escript.Data(1,(self.getNumSolutions(),),self.getFunctionSpaceOfSolution(),True)
773     else:
774     self.__applyConstraint(rhs_update=False)
775     self.__operator_isValid=True
776     return self.__operator
777    
778     def getRightHandSide(self,ignoreConstraint=False):
779     """
780     @brief returns the right hand side of the PDE
781    
782     @param ignoreConstraint
783     """
784     if not self.__righthandside_isValid:
785     if self.debug() : print "PDE Debug: New right hand side is built."
786     self.getDomain().addPDEToRHS(self.__getFreshRightHandSide(), \
787     self.getCoefficient("X"), \
788     self.getCoefficient("Y"),\
789     self.getCoefficient("y"),\
790     self.getCoefficient("y_contact"))
791     self.__righthandside_isValid=True
792     if ignoreConstraint: self.__copyConstraint(self.__righthandside)
793     return self.__righthandside
794    
795     def getSystem(self):
796     """
797     @brief
798     """
799     if not self.__operator_isValid and not self.__righthandside_isValid:
800     if self.debug() : print "PDE Debug: New PDE is built."
801     if self.isUsingLumping():
802     self.getRightHandSide(ignoreConstraint=True)
803     self.getOperator()
804     else:
805     self.getDomain().addPDEToSystem(self.__getFreshOperator(),self.__getFreshRightHandSide(), \
806     self.getCoefficient("A"), \
807     self.getCoefficient("B"), \
808     self.getCoefficient("C"), \
809     self.getCoefficient("D"), \
810     self.getCoefficient("X"), \
811     self.getCoefficient("Y"), \
812     self.getCoefficient("d"), \
813     self.getCoefficient("y"), \
814     self.getCoefficient("d_contact"), \
815     self.getCoefficient("y_contact"))
816     self.__operator_isValid=True
817     self.__righthandside_isValid=True
818     self.__applyConstraint()
819     self.__copyConstraint(self.__righthandside)
820     elif not self.__operator_isValid:
821     self.getOperator()
822     elif not self.__righthandside_isValid:
823     self.getRightHandSide()
824     return (self.__operator,self.__righthandside)
825    
826     def solve(self,**options):
827     """
828     @brief solve the PDE
829    
830     @param options
831     """
832     mat,f=self.getSystem()
833     if self.isUsingLumping():
834     out=f/mat
835     self.__copyConstraint(out)
836     else:
837     options[util.TOLERANCE_KEY]=self.getTolerance()
838     options[util.METHOD_KEY]=self.getSolverMethod()
839     options[util.SYMMETRY_KEY]=self.isSymmetric()
840     if self.debug() : print "PDE Debug: solver options: ",options
841     out=mat.solve(f,options)
842     return out
843    
844     def getSolution(self,**options):
845     """
846     @brief returns the solution of the PDE
847    
848     @param options
849     """
850     if not self.__solution_isValid:
851     if self.debug() : print "PDE Debug: PDE is resolved."
852     self.__solution=self.solve(**options)
853     self.__solution_isValid=True
854     return self.__solution
855     #============ some serivice functions =====================================================
856     def getDomain(self):
857     """
858     @brief returns the domain of the PDE
859     """
860     return self.__domain
861    
862 jgs 104 def getDim(self):
863     """
864     @brief returns the spatial dimension of the PDE
865     """
866     return self.getDomain().getDim()
867    
868 jgs 102 def getNumEquations(self):
869     """
870     @brief returns the number of equations
871     """
872     if self.__numEquations>0:
873     return self.__numEquations
874     else:
875     raise ValueError,"Number of equations is undefined. Please specify argument numEquations."
876    
877     def getNumSolutions(self):
878     """
879     @brief returns the number of unknowns
880     """
881     if self.__numSolutions>0:
882     return self.__numSolutions
883     else:
884     raise ValueError,"Number of solution is undefined. Please specify argument numSolutions."
885    
886    
887     def checkSymmetry(self):
888     """
889     @brief returns if the Operator is symmetric. This is a very expensive operation!!! The symmetry flag is not altered.
890     """
891     raise SystemError,"checkSymmetry is not implemented yet"
892    
893     return None
894    
895     def getFlux(self,u):
896     """
897     @brief returns the flux J_ij for a given u
898    
899     J_ij=A_{ijkl}u_{k,l}+B_{ijk}u_k-X_{ij}
900    
901     @param u argument of the operator
902    
903     """
904     raise SystemError,"getFlux is not implemented yet"
905     return None
906    
907     def applyOperator(self,u):
908     """
909     @brief applies the operator of the PDE to a given solution u in weak from
910    
911     @param u argument of the operator
912    
913     """
914     return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())
915    
916     def getResidual(self,u):
917     """
918     @brief return the residual of u in the weak from
919    
920     @param u
921     """
922     return self.applyOperator(u)-self.getRightHandSide()
923    
924     class Poisson(LinearPDE):
925     """
926     @brief Class to define a Poisson equstion problem:
927    
928     class to define a linear PDE of the form
929    
930     -u_{,jj} = f
931    
932     with boundary conditons:
933    
934     n_j*u_{,j} = 0
935    
936     and constraints:
937    
938     u=0 where q>0
939    
940     """
941    
942     def __init__(self,domain=None,f=escript.Data(),q=escript.Data()):
943     LinearPDE.__init__(self,domain=identifyDomain(domain,{"f":f, "q":q}))
944 jgs 104 self._setValue(A=numarray.identity(self.getDomain().getDim()))
945 jgs 102 self.setSymmetryOn()
946     self.setValue(f,q)
947    
948     def setValue(self,f=escript.Data(),q=escript.Data()):
949 jgs 104 self._setValue(Y=f,q=q)
950 jgs 102
951    
952     # $Log$
953 jgs 104 # Revision 1.3 2004/12/17 07:43:10 jgs
954 jgs 102 # *** empty log message ***
955     #
956 jgs 104 # Revision 1.1.2.3 2004/12/16 00:12:34 gross
957     # __init__ of LinearPDE does not accept any coefficients anymore
958     #
959 jgs 102 # Revision 1.1.2.2 2004/12/14 03:55:01 jgs
960     # *** empty log message ***
961     #
962     # Revision 1.1.2.1 2004/12/12 22:53:47 gross
963     # linearPDE has been renamed LinearPDE
964     #
965     # Revision 1.1.1.1.2.7 2004/12/07 10:13:08 gross
966     # GMRES added
967     #
968     # Revision 1.1.1.1.2.6 2004/12/07 03:19:50 gross
969     # options for GMRES and PRES20 added
970     #
971     # Revision 1.1.1.1.2.5 2004/12/01 06:25:15 gross
972     # some small changes
973     #
974     # Revision 1.1.1.1.2.4 2004/11/24 01:50:21 gross
975     # Finley solves 4M unknowns now
976     #
977     # Revision 1.1.1.1.2.3 2004/11/15 06:05:26 gross
978     # poisson solver added
979     #
980     # Revision 1.1.1.1.2.2 2004/11/12 06:58:15 gross
981     # 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
982     #
983     # Revision 1.1.1.1.2.1 2004/10/28 22:59:22 gross
984     # finley's RecTest.py is running now: problem in SystemMatrixAdapater fixed
985     #
986     # Revision 1.1.1.1 2004/10/26 06:53:56 jgs
987     # initial import of project esys2
988     #
989     # Revision 1.3.2.3 2004/10/26 06:43:48 jgs
990     # committing Lutz's and Paul's changes to brach jgs
991     #
992     # Revision 1.3.4.1 2004/10/20 05:32:51 cochrane
993     # Added incomplete Doxygen comments to files, or merely put the docstrings that already exist into Doxygen form.
994     #
995     # Revision 1.3 2004/09/23 00:53:23 jgs
996     # minor fixes
997     #
998     # Revision 1.1 2004/08/28 12:58:06 gross
999     # SimpleSolve is not running yet: problem with == of functionsspace
1000     #
1001     #

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26