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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

1 # $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