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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 82 - (hide annotations)
Tue Oct 26 06:53:54 2004 UTC (18 years, 3 months ago) by jgs
File MIME type: text/x-python
File size: 26373 byte(s)
Initial revision

1 jgs 82 # $Id$
2    
3     ## @file linearPDE.py
4    
5     """
6     @brief Functions and classes for linear PDEs
7     """
8    
9     import escript
10     import util
11    
12     def identifyDomain(domain=None,data={}):
13     """
14     @brief Return the Domain which is equal to the input domain (if not None)
15     and is the domain of all Data objects in the dictionary data.
16     An exception is raised if this is not possible
17    
18     @param domain
19     @param data
20     """
21     # get the domain used by any Data object in the list data:
22     data_domain=None
23     for d in data.itervalues():
24     if isinstance(d,escript.Data):
25     data_domain=d.getDomain()
26     print "identifyDomain: data_domain: ", data_domain
27     print "identifyDomain: domain: ", domain
28     # check if domain and data_domain are identical?
29     if domain == None:
30     if data_domain == None:
31     raise ValueError,"Undefined PDE domain. Specify a domain or use a Data class object as coefficient"
32     else:
33     if data_domain == None:
34     data_domain=domain
35     else:
36     if not data_domain==domain:
37     raise ValueError,"Domain of coefficients doesnot match specified domain"
38     # now we check if all Data class object coefficients are defined on data_domain:
39     for i,d in data.iteritems():
40     if isinstance(d,escript.Data):
41     if not data_domain==d.getDomain():
42     raise ValueError,"Illegal domain for coefficient %s."%i
43     # done:
44     return data_domain
45    
46    
47     def _CompTuple2(t1,t2):
48     """
49     @brief
50    
51     @param t1
52     @param t2
53     """
54     dif=t1[0]+t1[1]-(t2[0]+t2[1])
55     if dif<0: return 1
56     elif dif>0: return -1
57     else: return 0
58    
59     class PDECoefficientType:
60     """
61     @brief
62     """
63     # identifier for location of Data objects defining coefficients
64     INTERIOR=0
65     BOUNDARY=1
66     CONTACT=2
67     CONTINUOUS=3
68     # identifier in the pattern of coefficients:
69     # the pattern is a tuple of EQUATION,SOLUTION,DIM where DIM represents the spatial dimension, EQUATION the number of equations and SOLUTION the
70     # number of unknowns.
71     EQUATION=3
72     SOLUTION=4
73     DIM=5
74     # indicator for what is altered if the coefficient is altered:
75     OPERATOR=5
76     RIGHTHANDSIDE=6
77     BOTH=7
78     def __init__(self,where,pattern,altering):
79     """
80     @brief Initialise a PDE Coefficient type
81     """
82     self.what=where
83     self.pattern=pattern
84     self.altering=altering
85    
86     def getFunctionSpace(self,domain):
87     """
88     @brief defines the FunctionSpace of the coefficient on the domain
89    
90     @param domain
91     """
92     if self.what==self.INTERIOR: return escript.Function(domain)
93     elif self.what==self.BOUNDARY: return escript.FunctionOnBoundary(domain)
94     elif self.what==self.CONTACT: return escript.FunctionOnContactZero(domain)
95     elif self.what==self.CONTINUOUS: return escript.ContinuousFunction(domain)
96    
97     def isAlteringOperator(self):
98     """
99     @brief return true if the operator of the PDE is changed when the coefficient is changed
100     """
101     if self.altering==self.OPERATOR or self.altering==self.BOTH:
102     return not None
103     else:
104     return None
105    
106     def isAlteringRightHandSide(self):
107     """
108     @brief return true if the right hand side of the PDE is changed when the coefficient is changed
109     """
110     if self.altering==self.RIGHTHANDSIDE or self.altering==self.BOTH:
111     return not None
112     else:
113     return None
114    
115     def estimateNumEquationsAndNumSolutions(self,shape=(),dim=3):
116     """
117     @brief tries to estimate the number of equations in a given tensor shape for a given spatial dimension dim
118    
119     @param shape
120     @param dim
121     """
122     if len(shape)>0:
123     num=max(shape)+1
124     else:
125     num=1
126     search=[]
127     for u in range(num):
128     for e in range(num):
129     search.append((e,u))
130     search.sort(_CompTuple2)
131     for item in search:
132     s=self.buildShape(item[0],item[1],dim)
133     if s==shape: return item
134     return None
135    
136     def buildShape(self,e=1,u=1,dim=3):
137     """
138     @brief builds the required shape for a given number of equations e, number of unknowns u and spatial dimension dim
139    
140     @param e
141     @param u
142     @param dim
143     """
144     s=()
145     for i in self.pattern:
146     if i==self.EQUATION:
147     if e>1: s=s+(e,)
148     elif i==self.SOLUTION:
149     if u>1: s=s+(u,)
150     else:
151     s=s+(dim,)
152     return s
153    
154     _PDECoefficientTypes={
155     "A" : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,PDECoefficientType.DIM,PDECoefficientType.SOLUTION,PDECoefficientType.DIM),PDECoefficientType.OPERATOR),
156     "B" : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,PDECoefficientType.DIM,PDECoefficientType.SOLUTION),PDECoefficientType.OPERATOR),
157     "C" : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,PDECoefficientType.SOLUTION,PDECoefficientType.DIM),PDECoefficientType.OPERATOR),
158     "D" : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,PDECoefficientType.SOLUTION),PDECoefficientType.OPERATOR),
159     "X" : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,PDECoefficientType.DIM),PDECoefficientType.RIGHTHANDSIDE),
160     "Y" : PDECoefficientType(PDECoefficientType.INTERIOR,(PDECoefficientType.EQUATION,),PDECoefficientType.RIGHTHANDSIDE),
161     "d" : PDECoefficientType(PDECoefficientType.BOUNDARY,(PDECoefficientType.EQUATION,PDECoefficientType.SOLUTION),PDECoefficientType.OPERATOR),
162     "y" : PDECoefficientType(PDECoefficientType.BOUNDARY,(PDECoefficientType.EQUATION,),PDECoefficientType.RIGHTHANDSIDE),
163     "d_contact" : PDECoefficientType(PDECoefficientType.CONTACT,(PDECoefficientType.EQUATION,PDECoefficientType.SOLUTION),PDECoefficientType.OPERATOR),
164     "y_contact" : PDECoefficientType(PDECoefficientType.CONTACT,(PDECoefficientType.EQUATION,),PDECoefficientType.RIGHTHANDSIDE),
165     "r" : PDECoefficientType(PDECoefficientType.CONTINUOUS,(PDECoefficientType.EQUATION,),PDECoefficientType.RIGHTHANDSIDE),
166     "q" : PDECoefficientType(PDECoefficientType.CONTINUOUS,(PDECoefficientType.SOLUTION,),PDECoefficientType.BOTH),
167     }
168    
169     class linearPDE:
170     """
171     @brief Class to define a linear PDE
172    
173     class to define a linear PDE of the form
174    
175     -(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
176    
177     with boundary conditons:
178    
179     n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_{ik}u_k = - n_j*X_{ij} + y_i
180    
181     and contact conditions
182    
183     n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_contact_{ik}[u_k] = - n_j*X_{ij} + y_contact_i
184    
185     and constraints:
186    
187     u_i=r_i where q_i>0
188    
189     """
190    
191     def __init__(self,**args):
192     """
193     @brief initializes a new linear PDE.
194    
195     @param args
196     """
197    
198     print "Creating new linearPDE"
199    
200     # initialize attributes
201     self.__debug=None
202     self.__domain=None
203     self.__numEquations=0
204     self.__numSolutions=0
205     self.__coefficient={}
206     for i in _PDECoefficientTypes.iterkeys():
207     self.__coefficient[i]=escript.Data()
208     self.__operator=escript.Operator()
209     self.__righthandside=escript.Data()
210     self.__solution=escript.Data()
211     self.__solveroptions=None
212     self.__matrixtype=util.UNKNOWN
213     self.__sym=False
214     coef={}
215    
216     # check the arguments
217     for arg in args.iterkeys():
218     if arg=="domain":
219     self.__domain=args[arg]
220     elif arg=="numEquations":
221     self.__numEquations=args[arg]
222     elif arg=="numSolutions":
223     self.__numSolutions=args[arg]
224     elif _PDECoefficientTypes.has_key(arg):
225     coef[arg]=args[arg]
226     else:
227     raise ValueError,"Illegal argument %s"%arg
228    
229     # get the domain of the PDE
230     self.__domain=identifyDomain(self.__domain,coef)
231    
232     # now each coefficient is changed into a Data object if it is not already one or is None
233     for key in coef.iterkeys():
234     self.__coefficient[key]=escript.Data(coef[key],_PDECoefficientTypes[key].getFunctionSpace(self.__domain))
235    
236     # get number of equations and number of unknowns:
237     numEquations=0
238     numSolutions=0
239     numDim=self.__domain.getDim()
240     for i in self.__coefficient.iterkeys():
241     if not self.__coefficient[i].isEmpty():
242     res=_PDECoefficientTypes[i].estimateNumEquationsAndNumSolutions(self.__coefficient[i].getShape(),numDim)
243     if res==None:
244     raise ValueError,"Illegal shape %s of coefficient %s"%(self.__coefficient[i].getShape().__str__(),i)
245     else:
246     numEquations=max(numEquations,res[0])
247     numSolutions=max(numSolutions,res[1])
248    
249     # check the number of equations:
250     if numEquations>0:
251     if self.__numEquations>0:
252     if self.__numEquations!=numEquations:
253     raise ValueError,"Shape of coefficients is inconsistent with the specified number of equations (=%d)"%self.__numEquations
254     else:
255     self.__numEquations=numEquations
256     else:
257     if self.__numEquations<1:
258     raise ValueError,"Number of equations could not be identified. Please specify argument numEquations."
259    
260     # check the number of unknowns:
261     if numSolutions>0:
262     if self.__numSolutions>0:
263     if self.__numSolutions!=numSolutions:
264     raise ValueError,"Shape of coefficients is inconsistent with the specified number of unknowns (=%d)"%self.__numSolutions
265     else:
266     self.__numSolutions=numSolutions
267     else:
268     if self.__numSolutions<1: self.__numSolutions=self.__numEquations
269    
270     print "Identified number of equations and unknowns: (",self.__numEquations,self.__numSolutions,")"
271    
272     # check the shape of all coefficients:
273    
274     for i,d in self.__coefficient.iteritems():
275     if not d.isEmpty():
276     if not self.getShapeOfCoefficient(i)==d.getShape():
277     raise ValueError,"Expected shape for coefficient %s is %s but actual shape is %s."%(i,self.getShape(i).__str__(),d.getShape().__str__())
278     #
279     self.__row_function_space=escript.Solution(self.__domain)
280     self.__column_function_space=escript.Solution(self.__domain)
281    
282     def setDebugOn(self):
283     """
284     @brief
285     """
286     self.__debug=not None
287    
288     def setDebugOff(self):
289     """
290     @brief
291     """
292     self.__debug=None
293    
294     def debug(self):
295     """
296     @brief returns true if the PDE is in the debug mode
297     """
298     return self.__debug
299    
300     def getCoefficient(self,name):
301     """
302     @brief return the value of the coefficient name
303    
304     @param name
305     """
306     return self.__coefficient[name]
307    
308     def setCoefficient(self,**coefficients):
309     """
310     @brief sets new values to coefficients
311    
312     @param coefficients
313     """
314     alteredCoefficients=[]
315     for i,d in coefficients.iteritems():
316     if self.hasCoefficient(i):
317     if d == None:
318     if not self.__coefficient[i].isEmpty():
319     self.__coefficient[i]=escript.Data()
320     alteredCoefficients.append(i)
321     elif isinstance(d,escript.Data):
322     if d.isEmpty():
323     if not self.__coefficient[i].isEmpty():
324     self.__coefficient[i]=escript.Data()
325     alteredCoefficients.append(i)
326     else:
327     if not self.getShapeOfCoefficient(i)==d.getShape():
328     raise ValueError,"Illegal shape for coefficient %s"%i
329     if not self.getDomain()==d.getDomain():
330     raise ValueError,"Illegal domain for coefficient %s"%i
331     self.__coefficient[i]=escript.Data(d,self.getFunctionSpaceOfCoefficient(i))
332     alteredCoefficients.append(i)
333     else:
334     self.__coefficient[i]=escript.Data(d,self.getShapeOfCoefficient(i),self.getFunctionSpaceOfCoefficient(i))
335     alteredCoefficients.append(i)
336     if self.debug(): print "PDE Debug: Coefficient %s has been altered."%i
337     else:
338     raise ValueError,"Attempt to set unknown coefficient %s"%i
339    
340     if len(alteredCoefficients) > 0:
341     inhomogeneous=None
342     # even if q has not been altered, the system has to be rebuilt if the constraint is not homogeneous:
343     if not "q" in alteredCoefficients:
344     q=self.getCoefficient("q")
345     r=self.getCoefficient("r")
346     if not q.isEmpty() and not r.isEmpty():
347     if (q*r).Lsup()>1.e-13*r.Lsup(): inhomogeneous=not None
348    
349     if inhomogeneous:
350     if self.debug() : print "PDE Debug: System has to be rebuilt because of inhomogeneous constraints."
351     self.rebuildSystem()
352     else:
353     for i in alteredCoefficients:
354     self.alteredCoefficient(i)
355    
356     def hasCoefficient(self,name):
357     """
358     @brief return true if name is the name of a coefficient
359    
360     @param name
361     """
362     return self.__coefficient.has_key(name)
363    
364     def getFunctionSpaceForEquation(self):
365     """
366     @brief return true if the test functions should use reduced order
367     """
368     return self.__row_function_space
369    
370     def getFunctionSpaceForSolution(self):
371     """
372     @brief return true if the interpolation of the solution should use reduced order
373     """
374     return self.__column_function_space
375    
376     def getMatrixType(self):
377     """
378     @brief returns the matrix type to be used to store the operator
379     """
380     return self.__matrixtype
381    
382     def setMatrixType(self,type=util.UNKNOWN):
383     """
384     @brief sets the new matrix type
385    
386     @param type
387     """
388     if not type==self.__matrixtype:
389     if self.debug() : print "PDE Debug: Matrix type is now %d."%type
390     self.__matrixtype=type
391     self.rebuildOperator()
392    
393     def isSymmetric(self):
394     """
395     @brief returns true is the operator is considered to be symmetric
396     """
397     return self.__sym
398    
399     def setReducedOrderForSolutionsOn(self):
400     """
401     @brief switches to reduced order to interpolate solution
402     """
403     new_fs=escript.ReducedSolution(self.getDomain())
404     if self.getFunctionSpaceForSolution()!=new_fs:
405     if self.debug() : print "PDE Debug: Reduced order is used to interpolate solution."
406     self.__column_function_space=new_fs
407     self.rebuildOperator()
408    
409     def setReducedOrderForSolutionsOff(self):
410     """
411     @brief switches to full order to interpolate solution
412     """
413     new_fs=escript.Solution(self.getDomain())
414     if self.getFunctionSpaceForSolution()!=new_fs:
415     if self.debug() : print "PDE Debug: Full order is used to interpolate solution."
416     self.__column_function_space=new_fs
417     self.rebuildOperator()
418    
419     def setReducedOrderForEquationsOn(self):
420     """
421     @brief switches to reduced order for test functions
422     """
423     new_fs=escript.ReducedSolution(self.getDomain())
424     if self.getFunctionSpaceForEquation()!=new_fs:
425     if self.debug() : print "PDE Debug: Reduced order is used for test functions."
426     self.__row_function_space=new_fs
427     self.rebuildOperator()
428     self.rebuildRightHandSide()
429    
430     def setReducedOrderForSolutionsOff(self):
431     """
432     @brief switches to full order for test functions
433     """
434     new_fs=escript.Solution(self.getDomain())
435     if self.getFunctionSpaceForEquation()!=new_fs:
436     if self.debug() : print "PDE Debug: Full order is used for test functions."
437     self.__row_function_space=new_fs
438     self.rebuildOperator()
439     self.rebuildRightHandSide()
440    
441     def getDomain(self):
442     """
443     @brief returns the domain of the PDE
444     """
445     return self.__domain
446    
447     def getNumEquations(self):
448     """
449     @brief returns the number of equations
450     """
451     return self.__numEquations
452    
453     def getNumSolutions(self):
454     """
455     @brief returns the number of unknowns
456     """
457     return self.__numSolutions
458    
459     def getShapeOfCoefficient(self,name):
460     """
461     @brief return the shape of the coefficient name
462    
463     @param name
464     """
465     if self.hasCoefficient(name):
466     return _PDECoefficientTypes[name].buildShape(self.getNumEquations(),self.getNumSolutions(),self.getDomain().getDim())
467     else:
468     raise ValueError,"Solution coefficient %s requested"%name
469    
470     def getFunctionSpaceOfCoefficient(self,name):
471     """
472     @brief return the atoms of the coefficient name
473    
474     @param name
475     """
476     if self.hasCoefficient(name):
477     return _PDECoefficientTypes[name].getFunctionSpace(self.getDomain())
478     else:
479     raise ValueError,"Solution coefficient %s requested"%name
480    
481     def alteredCoefficient(self,name):
482     """
483     @brief annonced that coefficient name has been changed
484    
485     @param name
486     """
487     if self.hasCoefficient(name):
488     if _PDECoefficientTypes[name].isAlteringOperator(): self.rebuildOperator()
489     if _PDECoefficientTypes[name].isAlteringRightHandSide(): self.rebuildRightHandSide()
490     else:
491     raise ValueError,"Solution coefficient %s requested"%name
492    
493     def initialiseNewOperator(self):
494     """
495     @brief
496     """
497     return self.getDomain().newOperator( \
498     self.getNumEquations(), \
499     self.getFunctionSpaceForEquation(), \
500     self.getNumSolutions(), \
501     self.getFunctionSpaceForSolution(), \
502     self.getMatrixType(), \
503     self.isSymmetric())
504    
505     def initialiseNewRightHandSide(self,expanded=True):
506     """
507     @brief
508    
509     @param expanded
510     """
511     return escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),expanded)
512    
513     def initialiseNewSolution(self,expanded=True):
514     """
515     @brief
516    
517     @param expanded
518     """
519     return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),expanded)
520    
521     def rebuildSolution(self):
522     """
523     @brief indicates the PDE has to be reolved if the solution is requested
524     """
525     if not self.__solution.isEmpty() and self.debug() : print "PDE Debug: PDE has to be resolved."
526     self.__solution=escript.Data()
527    
528     def rebuildOperator(self):
529     """
530     @brief indicates the operator has to be rebuilt next time it is used
531     """
532     if not self.__operator.isEmpty() and self.debug() : print "PDE Debug: Operator has to be rebuilt."
533     self.rebuildSolution()
534     self.__operator=escript.Operator()
535    
536     def rebuildRightHandSide(self):
537     """
538     @brief indicates the right hand side has to be rebuild next time it is used
539     """
540     if not self.__righthandside.isEmpty() and self.debug() : print "PDE Debug: Right hand side has to be rebuilt."
541     self.rebuildSolution()
542     self.__righthandside=escript.Data()
543    
544     def rebuildSystem(self):
545     """
546     @brief annonced that all coefficient name has been changed
547     """
548     self.rebuildSolution()
549     self.rebuildOperator()
550     self.rebuildRightHandSide()
551    
552     def applyConstraint(self):
553     """
554     @brief applies the constraints defined by q and r to the system
555     """
556     q=self.getCoefficient("q")
557     r=self.getCoefficient("r")
558     rhs=self.getRightHandSide()
559     mat=self.getOperator()
560     if not q.isEmpty():
561     # q is the row and column mask to indicate where constraints are set:
562     row_q=escript.Data(q,self.getFunctionSpaceForEquation())
563     col_q=escript.Data(q,self.getFunctionSpaceForSolution())
564     if r.isEmpty():
565     r2=self.initialiseNewRightHandSide()
566     else:
567     u=self.initialiseNewSolution()
568     src=escript.Data(r,self.getFunctionSpaceForEquation())
569     u.copyWithMask(src,row_q)
570     if not rhs.isEmpty(): rhs-=mat*u
571     r2=escript.Data(r,self.getFunctionSpaceForEquation())
572     if not rhs.isEmpty(): rhs.copyWithMask(r2,row_q)
573     if not mat.isEmpty(): mat.nullifyRowsAndCols(row_q,col_q,1.)
574    
575     def getOperator(self):
576     """
577     @brief returns the operator of the PDE
578     """
579     if self.__operator.isEmpty():
580     if self.debug() : print "PDE Debug: New operator is built."
581     # some constrainst are applying for a lumpled stifness matrix:
582     if self.getMatrixType()==util.LUMPED:
583     if self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():
584     raise Warning,"Lumped matrix requires same order for equations and unknowns"
585     if not self.getCoefficient("A").isEmpty():
586     raise Warning,"Lumped matrix does not allow coefficient A"
587     if not self.getCoefficient("B").isEmpty():
588     raise Warning,"Lumped matrix does not allow coefficient B"
589     if not self.getCoefficient("C").isEmpty():
590     raise Warning,"Lumped matrix does not allow coefficient C"
591     if not self.getCoefficient("D").isEmpty():
592     raise Warning,"Lumped matrix does not allow coefficient D"
593     # get a new matrix:
594     mat=self.initialiseNewOperator()
595     #
596     # assemble stiffness matrix:
597     #
598     self.getDomain().addPDEToSystem(mat,escript.Data(), \
599     self.getCoefficient("A"), \
600     self.getCoefficient("B"), \
601     self.getCoefficient("C"), \
602     self.getCoefficient("D"), \
603     escript.Data(), \
604     escript.Data())
605     self.getDomain().addRobinConditionsToSystem(mat,escript.Data(), \
606     self.getCoefficient("d"), \
607     escript.Data())
608     self.getDomain().addContactToSystem(mat,escript.Data(), \
609     self.getCoefficient("d_contact"), \
610     escript.Data())
611     self.applyConstraint()
612     if self.getMatrixType()==util.LUMPED:
613     self.__operator=mat*escript.Data(1,(self.getNumSolutions(),),self.getFunctionSpaceOfSolution(),True)
614     else:
615     self.__operator=mat
616     return self.__operator
617    
618     def getRightHandSide(self,ignoreConstraint=None):
619     """
620     @brief returns the right hand side of the PDE
621    
622     @param ignoreConstraint
623     """
624     if self.__righthandside.isEmpty():
625     if self.debug() : print "PDE Debug: New right hand side is built."
626     self.__righthandside=self.initialiseNewRightHandSide()
627     self.getDomain().addPDEToSystem(escript.Operator(),self.__righthandside, \
628     escript.Data(), \
629     escript.Data(), \
630     escript.Data(), \
631     escript.Data(), \
632     self.getCoefficient("X"), \
633     self.getCoefficient("Y"))
634     self.getDomain().addRobinConditionsToSystem(escript.Operator(),self.__righthandside, \
635     escript.Data(), \
636     self.getCoefficient("y"))
637     self.getDomain().addContactToSystem(escript.Operator(),self.__righthandside, \
638     escript.Data(), \
639     self.getCoefficient("y_contact"))
640     if not ignoreConstraint: self.applyConstraint()
641     return self.__righthandside
642    
643     def getSystem(self):
644     """
645     @brief
646     """
647     if self.__operator.isEmpty() and self.__righthandside.isEmpty():
648     if self.debug() : print "PDE Debug: New PDE is built."
649     if self.getMatrixType()==util.LUMPED:
650     self.getRightHandSide(ignoreConstraint=not None)
651     self.getOperator()
652     else:
653     self.__righthandside=self.initialiseNewRightHandSide()
654     self.__operator=self.initialiseNewOperator()
655     self.getDomain().addPDEToSystem(self.__operator,self.__righthandside, \
656     self.getCoefficient("A"), \
657     self.getCoefficient("B"), \
658     self.getCoefficient("C"), \
659     self.getCoefficient("D"), \
660     self.getCoefficient("X"), \
661     self.getCoefficient("Y"))
662     self.getDomain().addRobinConditionsToSystem(self.__operator,self.__righthandside, \
663     self.getCoefficient("d"), \
664     self.getCoefficient("y"))
665     self.getDomain().addContactToSystem(self.__operator,self.__righthandside, \
666     self.getCoefficient("d_contact"), \
667     self.getCoefficient("y_contact"))
668     self.applyConstraint()
669     elif self.__operator.isEmpty():
670     self.getOperator()
671     elif self.__righthandside.isEmpty():
672     self.getRightHandSide()
673     return (self.__operator,self.__righthandside)
674    
675     def solve(self,**options):
676     """
677     @brief solve the PDE
678    
679     @param options
680     """
681     if not options.has_key("iterative"): options["iterative"]=True
682     mat,f=self.getSystem()
683     if isinstance(mat,escript.Data):
684     return f/mat
685     else:
686     return mat.solve(f,options)
687    
688     def getSolution(self,**options):
689     """
690     @brief returns the solution of the PDE
691    
692     @param options
693     """
694     if self.__solution.isEmpty() or not self.__solveroptions==options:
695     self.__solveroptions=options
696     self.__solution=self.solve(**options)
697     if self.debug() : print "PDE Debug: PDE is resolved."
698     return self.__solution
699    
700     # $Log$
701     # Revision 1.1 2004/10/26 06:53:56 jgs
702     # Initial revision
703     #
704     # Revision 1.3.2.3 2004/10/26 06:43:48 jgs
705     # committing Lutz's and Paul's changes to brach jgs
706     #
707     # Revision 1.3.4.1 2004/10/20 05:32:51 cochrane
708     # Added incomplete Doxygen comments to files, or merely put the docstrings that already exist into Doxygen form.
709     #
710     # Revision 1.3 2004/09/23 00:53:23 jgs
711     # minor fixes
712     #
713     # Revision 1.1 2004/08/28 12:58:06 gross
714     # SimpleSolve is not running yet: problem with == of functionsspace
715     #
716     #

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26