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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26