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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 148 - (hide annotations)
Tue Aug 23 01:24:31 2005 UTC (14 years ago) by jgs
Original Path: trunk/esys2/escript/py_src/linearPDEs.py
File MIME type: text/x-python
File size: 68376 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-08-23

1 jgs 102 # $Id$
2    
3     ## @file linearPDEs.py
4    
5     """
6 jgs 122 Functions and classes for linear PDEs
7 jgs 102 """
8    
9     import escript
10     import util
11     import numarray
12    
13    
14 jgs 148 class IllegalCoefficient(ValueError):
15 jgs 102 """
16 jgs 148 raised if an illegal coefficient of the general ar particular PDE is requested.
17     """
18 jgs 102
19 jgs 148 class IllegalCoefficientValue(ValueError):
20 jgs 102 """
21 jgs 148 raised if an incorrect value for a coefficient is used.
22     """
23 jgs 122
24 jgs 148 class UndefinedPDEError(ValueError):
25     """
26     raised if a PDE is not fully defined yet.
27     """
28 jgs 102
29 jgs 148 def _CompTuple2(t1,t2):
30     """
31     Compare two tuples
32    
33     @param t1 The first tuple
34     @param t2 The second tuple
35    
36     """
37    
38     dif=t1[0]+t1[1]-(t2[0]+t2[1])
39     if dif<0: return 1
40     elif dif>0: return -1
41     else: return 0
42    
43 jgs 108 class PDECoefficient:
44 jgs 102 """
45 jgs 122 A class for PDE coefficients
46 jgs 102 """
47 jgs 108 # identifier for location of Data objects defining COEFFICIENTS
48 jgs 102 INTERIOR=0
49     BOUNDARY=1
50     CONTACT=2
51     CONTINUOUS=3
52 jgs 108 # identifier in the pattern of COEFFICIENTS:
53 jgs 102 # the pattern is a tuple of EQUATION,SOLUTION,DIM where DIM represents the spatial dimension, EQUATION the number of equations and SOLUTION the
54     # number of unknowns.
55     EQUATION=3
56     SOLUTION=4
57     DIM=5
58     # indicator for what is altered if the coefficient is altered:
59     OPERATOR=5
60     RIGHTHANDSIDE=6
61     BOTH=7
62     def __init__(self,where,pattern,altering):
63     """
64 jgs 122 Initialise a PDE Coefficient type
65 jgs 102 """
66     self.what=where
67     self.pattern=pattern
68     self.altering=altering
69 jgs 108 self.resetValue()
70 jgs 102
71 jgs 108 def resetValue(self):
72     """
73 jgs 122 resets coefficient value to default
74 jgs 108 """
75     self.value=escript.Data()
76    
77 jgs 102 def getFunctionSpace(self,domain):
78     """
79 jgs 122 defines the FunctionSpace of the coefficient on the domain
80 jgs 102
81 jgs 122 @param domain:
82 jgs 102 """
83     if self.what==self.INTERIOR: return escript.Function(domain)
84     elif self.what==self.BOUNDARY: return escript.FunctionOnBoundary(domain)
85     elif self.what==self.CONTACT: return escript.FunctionOnContactZero(domain)
86     elif self.what==self.CONTINUOUS: return escript.ContinuousFunction(domain)
87    
88 jgs 108 def getValue(self):
89     """
90 jgs 122 returns the value of the coefficient:
91 jgs 108 """
92     return self.value
93 jgs 148
94     def setValue(self,domain,numEquations=1,numSolutions=1,newValue=None):
95 jgs 108 """
96 jgs 122 set the value of the coefficient to new value
97 jgs 108 """
98 jgs 148 if newValue==None:
99     newValue=escript.Data()
100     elif isinstance(newValue,escript.Data):
101     if not newValue.isEmpty():
102     newValue=escript.Data(newValue,self.getFunctionSpace(domain))
103     else:
104     newValue=escript.Data(newValue,self.getFunctionSpace(domain))
105     if not newValue.isEmpty():
106     if not self.getShape(domain,numEquations,numSolutions)==newValue.getShape():
107     raise IllegalCoefficientValue,"Expected shape for coefficient %s is %s but actual shape is %s."%(self.getShape(domain,numEquations,numSolutions),newValue.getShape())
108 jgs 108 self.value=newValue
109 jgs 148
110 jgs 102 def isAlteringOperator(self):
111     """
112 jgs 122 return true if the operator of the PDE is changed when the coefficient is changed
113 jgs 102 """
114     if self.altering==self.OPERATOR or self.altering==self.BOTH:
115     return not None
116     else:
117     return None
118    
119     def isAlteringRightHandSide(self):
120     """
121 jgs 122 return true if the right hand side of the PDE is changed when the coefficient is changed
122 jgs 102 """
123     if self.altering==self.RIGHTHANDSIDE or self.altering==self.BOTH:
124     return not None
125     else:
126     return None
127    
128 jgs 148 def estimateNumEquationsAndNumSolutions(self,domain,shape=()):
129 jgs 102 """
130 jgs 122 tries to estimate the number of equations in a given tensor shape for a given spatial dimension dim
131 jgs 102
132 jgs 122 @param shape:
133     @param dim:
134 jgs 102 """
135 jgs 148 dim=domain.getDim()
136 jgs 102 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 jgs 148 s=self.getShape(domain,item[0],item[1])
147 jgs 102 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 jgs 148 def getShape(self,domain,numEquations=1,numSolutions=1):
154 jgs 102 """
155 jgs 122 builds the required shape for a given number of equations e, number of unknowns u and spatial dimension dim
156 jgs 102
157 jgs 122 @param e:
158     @param u:
159     @param dim:
160 jgs 102 """
161 jgs 148 dim=domain.getDim()
162 jgs 102 s=()
163     for i in self.pattern:
164     if i==self.EQUATION:
165 jgs 148 if numEquations>1: s=s+(numEquations,)
166 jgs 102 elif i==self.SOLUTION:
167 jgs 148 if numSolutions>1: s=s+(numSolutions,)
168 jgs 102 else:
169     s=s+(dim,)
170     return s
171    
172     class LinearPDE:
173     """
174 jgs 148 Class to define a linear PDE of the form
175 jgs 102
176 jgs 122 \f[
177 jgs 102 -(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
178 jgs 122 \f]
179 jgs 102
180 jgs 122 with boundary conditons:
181 jgs 102
182 jgs 122 \f[
183     n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_{ik}u_k = - n_j*X_{ij} + y_i
184     \f]
185 jgs 102
186 jgs 122 and contact conditions
187 jgs 102
188 jgs 122 \f[
189     n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_contact_{ik}[u_k] = - n_j*X_{ij} + y_contact_i
190     \f]
191 jgs 102
192 jgs 122 and constraints:
193 jgs 102
194 jgs 122 \f[
195     u_i=r_i \quad \mathrm{where} \quad q_i>0
196     \f]
197 jgs 102
198     """
199 jgs 108 TOL=1.e-13
200 jgs 148 # solver methods
201     UNKNOWN=-1
202     DEFAULT_METHOD=0
203     DIRECT=1
204     CHOLEVSKY=2
205     PCG=3
206     CR=4
207     CGS=5
208     BICGSTAB=6
209     SSOR=7
210     ILU0=8
211     ILUT=9
212     JACOBI=10
213     GMRES=11
214     PRES20=12
215     LUMPING=13
216     # matrix reordering:
217     NO_REORDERING=0
218     MINIMUM_FILL_IN=1
219     NESTED_DISSECTION=2
220     # important keys in the dictonary used to hand over options to the solver library:
221     METHOD_KEY="method"
222     SYMMETRY_KEY="symmetric"
223     TOLERANCE_KEY="tolerance"
224 jgs 102
225 jgs 148
226     def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):
227 jgs 102 """
228 jgs 148 initializes a new linear PDE
229 jgs 102
230 jgs 148 @param domain: domain of the PDE
231     @type domain: L{Domain}
232     @param numEquations: number of equations. If numEquations==None the number of equations
233     is exracted from the PDE coefficients.
234     @param numSolutions: number of solution components. If numSolutions==None the number of solution components
235     is exracted from the PDE coefficients.
236     @param debug: if True debug informations are printed.
237    
238    
239 jgs 102 """
240 jgs 148 #
241     # the coefficients of the general PDE:
242     #
243     self.__COEFFICIENTS_OF_GENEARL_PDE={
244 jgs 108 "A" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM,PDECoefficient.SOLUTION,PDECoefficient.DIM),PDECoefficient.OPERATOR),
245     "B" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),
246     "C" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION,PDECoefficient.DIM),PDECoefficient.OPERATOR),
247     "D" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),
248     "X" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM),PDECoefficient.RIGHTHANDSIDE),
249     "Y" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
250     "d" : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),
251     "y" : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
252     "d_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),
253     "y_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
254     "r" : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
255     "q" : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.SOLUTION,),PDECoefficient.BOTH)}
256 jgs 102
257 jgs 148 # COEFFICIENTS can be overwritten by subclasses:
258     self.COEFFICIENTS=self.__COEFFICIENTS_OF_GENEARL_PDE
259 jgs 102 # initialize attributes
260 jgs 148 self.__debug=debug
261 jgs 104 self.__domain=domain
262     self.__numEquations=numEquations
263     self.__numSolutions=numSolutions
264 jgs 148 self.__resetSystem()
265 jgs 102
266     # set some default values:
267     self.__homogeneous_constraint=True
268     self.__row_function_space=escript.Solution(self.__domain)
269     self.__column_function_space=escript.Solution(self.__domain)
270     self.__tolerance=1.e-8
271 jgs 148 self.__solver_method=self.DEFAULT_METHOD
272     self.__matrix_type=self.__domain.getSystemMatrixTypeId(self.DEFAULT_METHOD,False)
273 jgs 102 self.__sym=False
274    
275 jgs 148 self.resetCoefficients()
276     self.trace("PDE Coeffients are %s"%str(self.COEFFICIENTS.keys()))
277     # =============================================================================
278     # general stuff:
279     # =============================================================================
280     def __str__(self):
281     return "<LinearPDE %d>"%id(self)
282     # =============================================================================
283     # debug :
284     # =============================================================================
285     def setDebugOn(self):
286     """
287     switches on debugging
288 jgs 108 """
289 jgs 148 self.__debug=not None
290    
291     def setDebugOff(self):
292 jgs 108 """
293 jgs 148 switches off debugging
294     """
295     self.__debug=None
296 jgs 108
297 jgs 148 def trace(self,text):
298     """
299     print the text message if debugging is swiched on.
300 jgs 108
301 jgs 148 @param name: name of the coefficient enquired.
302     @type name: C{string}
303 jgs 102 """
304 jgs 148 if self.__debug: print "%s: %s"%(str(self),text)
305 jgs 102
306 jgs 148 # =============================================================================
307     # some service functions:
308     # =============================================================================
309     def getDomain(self):
310 jgs 102 """
311 jgs 148 returns the domain of the PDE
312    
313     @return : the domain of the PDE
314     @rtype : C{Domain}
315 jgs 102
316 jgs 108 """
317 jgs 148 return self.__domain
318 jgs 122
319 jgs 148 def getDim(self):
320 jgs 108 """
321 jgs 148 returns the spatial dimension of the PDE
322 jgs 108
323 jgs 148 @return : the spatial dimension of the PDE domain
324     @rtype : C{int}
325     """
326     return self.getDomain().getDim()
327 jgs 102
328 jgs 148 def getNumEquations(self):
329     """
330     returns the number of equations
331 jgs 102
332 jgs 148 @return : the number of equations
333     @rtype : C{int}
334     @raise UndefinedPDEError: if the number of equations is not be specified yet.
335     """
336     if self.__numEquations==None:
337     raise UndefinedPDEError,"Number of equations is undefined. Please specify argument numEquations."
338     else:
339     return self.__numEquations
340 jgs 147
341 jgs 148 def getNumSolutions(self):
342     """
343     returns the number of unknowns
344 jgs 147
345 jgs 148 @return : the number of unknowns
346     @rtype : C{int}
347     @raise UndefinedPDEError: if the number of unknowns is not be specified yet.
348     """
349     if self.__numSolutions==None:
350     raise UndefinedPDEError,"Number of solution is undefined. Please specify argument numSolutions."
351     else:
352     return self.__numSolutions
353    
354 jgs 108 def getFunctionSpaceForEquation(self):
355     """
356 jgs 148 returns the L{escript.FunctionSpace} used to discretize the equation
357    
358     @return : representation space of equation
359     @rtype : L{escript.FunctionSpace}
360    
361 jgs 108 """
362     return self.__row_function_space
363    
364     def getFunctionSpaceForSolution(self):
365     """
366 jgs 148 returns the L{escript.FunctionSpace} used to represent the solution
367    
368     @return : representation space of solution
369     @rtype : L{escript.FunctionSpace}
370    
371 jgs 108 """
372     return self.__column_function_space
373    
374 jgs 102
375 jgs 148 def getOperator(self):
376 jgs 102 """
377 jgs 148 provides access to the operator of the PDE
378 jgs 102
379 jgs 148 @return : the operator of the PDE
380     @rtype : L{Operator}
381 jgs 108 """
382 jgs 148 m=self.getSystem()[0]
383     if self.isUsingLumping():
384     return self.copyConstraint(1./m)
385 jgs 147 else:
386 jgs 148 return m
387 jgs 147
388 jgs 148 def getRightHandSide(self):
389 jgs 147 """
390 jgs 148 provides access to the right hand side of the PDE
391 jgs 147
392 jgs 148 @return : the right hand side of the PDE
393     @rtype : L{escript.Data}
394 jgs 147 """
395 jgs 148 r=self.getSystem()[1]
396     if self.isUsingLumping():
397     return self.copyConstraint(r)
398 jgs 147 else:
399 jgs 148 return r
400 jgs 147
401 jgs 148 def applyOperator(self,u=None):
402 jgs 102 """
403 jgs 148 applies the operator of the PDE to a given u or the solution of PDE if u is not present.
404 jgs 102
405 jgs 148 @param u: argument of the operator. It must be representable in C{elf.getFunctionSpaceForSolution()}. If u is not present or equals L{None}
406     the current solution is used.
407     @type u: L{escript.Data} or None
408     @return : image of u
409     @rtype : L{escript.Data}
410 jgs 102 """
411 jgs 148 if u==None:
412     return self.getOperator()*self.getSolution()
413 jgs 102 else:
414 jgs 148 self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())
415 jgs 102
416 jgs 148 def getResidual(self,u=None):
417 jgs 102 """
418 jgs 148 return the residual of u or the current solution if u is not present.
419 jgs 102
420 jgs 148 @param u: argument in the residual calculation. It must be representable in C{elf.getFunctionSpaceForSolution()}. If u is not present or equals L{None}
421     the current solution is used.
422     @type u: L{escript.Data} or None
423     @return : residual of u
424     @rtype : L{escript.Data}
425 jgs 102 """
426 jgs 148 return self.applyOperator(u)-self.getRightHandSide()
427 jgs 102
428 jgs 148 def checkSymmetry(self,verbose=True):
429     """
430     test the PDE for symmetry.
431 jgs 102
432    
433 jgs 148 @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.
434     @type verbose: C{bool}
435     @return: True if the PDE is symmetric.
436     @rtype : C{escript.Data}
437 jgs 102
438 jgs 148 @note: This is a very expensive operation. It should be used for degugging only! The symmetry flag is not altered.
439     """
440     verbose=verbose or self.debug()
441     out=True
442     if self.getNumSolutions()!=self.getNumEquations():
443     if verbose: print "non-symmetric PDE because of different number of equations and solutions"
444     out=False
445     else:
446     A=self.getCoefficientOfGeneralPDE("A")
447     if not A.isEmpty():
448     tol=util.Lsup(A)*self.TOL
449     if self.getNumSolutions()>1:
450     for i in range(self.getNumEquations()):
451     for j in range(self.getDim()):
452     for k in range(self.getNumSolutions()):
453     for l in range(self.getDim()):
454     if util.Lsup(A[i,j,k,l]-A[k,l,i,j])>tol:
455     if verbose: print "non-symmetric PDE because A[%d,%d,%d,%d]!=A[%d,%d,%d,%d]"%(i,j,k,l,k,l,i,j)
456     out=False
457     else:
458     for j in range(self.getDim()):
459     for l in range(self.getDim()):
460     if util.Lsup(A[j,l]-A[l,j])>tol:
461     if verbose: print "non-symmetric PDE because A[%d,%d]!=A[%d,%d]"%(j,l,l,j)
462     out=False
463     B=self.getCoefficientOfGeneralPDE("B")
464     C=self.getCoefficientOfGeneralPDE("C")
465     if B.isEmpty() and not C.isEmpty():
466     if verbose: print "non-symmetric PDE because B is not present but C is"
467     out=False
468     elif not B.isEmpty() and C.isEmpty():
469     if verbose: print "non-symmetric PDE because C is not present but B is"
470     out=False
471     elif not B.isEmpty() and not C.isEmpty():
472     tol=(util.Lsup(B)+util.Lsup(C))*self.TOL/2.
473     if self.getNumSolutions()>1:
474     for i in range(self.getNumEquations()):
475     for j in range(self.getDim()):
476     for k in range(self.getNumSolutions()):
477     if util.Lsup(B[i,j,k]-C[k,i,j])>tol:
478     if verbose: print "non-symmetric PDE because B[%d,%d,%d]!=C[%d,%d,%d]"%(i,j,k,k,i,j)
479     out=False
480     else:
481     for j in range(self.getDim()):
482     if util.Lsup(B[j]-C[j])>tol:
483     if verbose: print "non-symmetric PDE because B[%d]!=C[%d]"%(j,j)
484     out=False
485     if self.getNumSolutions()>1:
486     D=self.getCoefficientOfGeneralPDE("D")
487     if not D.isEmpty():
488     tol=util.Lsup(D)*self.TOL
489     for i in range(self.getNumEquations()):
490     for k in range(self.getNumSolutions()):
491     if util.Lsup(D[i,k]-D[k,i])>tol:
492     if verbose: print "non-symmetric PDE because D[%d,%d]!=D[%d,%d]"%(i,k,k,i)
493     out=False
494    
495     return out
496    
497     def getSolution(self,**options):
498 jgs 102 """
499 jgs 148 returns the solution of the PDE. If the solution is not valid the PDE is solved.
500 jgs 102
501 jgs 148 @return: the solution
502     @rtype: L{escript.Data}
503     @param options: solver options
504     @keyword verbose:
505     @keyword reordering: reordering scheme to be used during elimination
506     @keyword preconditioner: preconditioner method to be used
507     @keyword iter_max: maximum number of iteration steps allowed.
508     @keyword drop_tolerance:
509     @keyword drop_storage:
510     @keyword truncation:
511     @keyword restart:
512 jgs 102 """
513 jgs 148 if not self.__solution_isValid:
514     mat,f=self.getSystem()
515     if self.isUsingLumping():
516     self.__solution=self.copyConstraint(f*mat)
517     else:
518     options[self.TOLERANCE_KEY]=self.getTolerance()
519     options[self.METHOD_KEY]=self.getSolverMethod()
520     options[self.SYMMETRY_KEY]=self.isSymmetric()
521     self.trace("PDE is resolved.")
522     self.trace("solver options: %s"%str(options))
523     self.__solution=mat.solve(f,options)
524     self.__solution_isValid=True
525     return self.__solution
526 jgs 102
527 jgs 148 def getFlux(self,u=None):
528     """
529     returns the flux J_ij for a given u
530 jgs 102
531 jgs 148 \f[
532     J_ij=A_{ijkl}u_{k,l}+B_{ijk}u_k-X_{ij}
533     \f]
534 jgs 102
535 jgs 148 @param u: argument in the flux. If u is not present or equals L{None} the current solution is used.
536     @type u: L{escript.Data} or None
537     @return : flux
538     @rtype : L{escript.Data}
539 jgs 102
540 jgs 148 """
541     if u==None: u=self.getSolution()
542     return util.tensormult(self.getCoefficientOfGeneralPDE("A"),util.grad(u))+util.matrixmult(self.getCoefficientOfGeneralPDE("B"),u)-util.self.getCoefficientOfGeneralPDE("X")
543 jgs 102
544 jgs 148 # =============================================================================
545     # solver settings:
546     # =============================================================================
547     def setSolverMethod(self,solver=None):
548 jgs 102 """
549 jgs 122 sets a new solver
550 jgs 148
551     @param solver: sets a new solver method.
552     @type solver: C{int}
553    
554 jgs 102 """
555 jgs 148 if solver==None: solve=self.DEFAULT_METHOD
556 jgs 102 if not solver==self.getSolverMethod():
557     self.__solver_method=solver
558     self.__checkMatrixType()
559 jgs 148 self.trace("New solver is %s"%self.getSolverMethodName())
560 jgs 102
561 jgs 148 def getSolverMethodName(self):
562     """
563     returns the name of the solver currently used
564    
565     @return : the name of the solver currently used.
566     @rtype: C{string}
567     """
568    
569     m=self.getSolverMethod()
570     if m==self.DEFAULT_METHOD: return "DEFAULT_METHOD"
571     elif m==self.DIRECT: return "DIRECT"
572     elif m==self.CHOLEVSKY: return "CHOLEVSKY"
573     elif m==self.PCG: return "PCG"
574     elif m==self.CR: return "CR"
575     elif m==self.CGS: return "CGS"
576     elif m==self.BICGSTAB: return "BICGSTAB"
577     elif m==self.SSOR: return "SSOR"
578     elif m==self.GMRES: return "GMRES"
579     elif m==self.PRES20: return "PRES20"
580     elif m==self.LUMPING: return "LUMPING"
581     return None
582    
583    
584 jgs 102 def getSolverMethod(self):
585     """
586 jgs 122 returns the solver method
587 jgs 148
588     @return : the solver method currently be used.
589     @rtype : C{int}
590 jgs 102 """
591     return self.__solver_method
592    
593 jgs 148 def isUsingLumping(self):
594     """
595     checks if matrix lumping is used a solver method
596    
597     @return : True is lumping is currently used a solver method.
598     @rtype: C{bool}
599     """
600     return self.getSolverMethod()==self.LUMPING
601    
602 jgs 102 def setTolerance(self,tol=1.e-8):
603     """
604 jgs 148 resets the tolerance for the solver method to tol where for an appropriate norm |.|
605    
606     |self.getResidual()|<tol*|self.getRightHandSide()|
607    
608     defines the stopping criterion.
609    
610     @param tol: new tolerance for the solver. If the tol is lower then the current tolerence
611     the system will be resolved.
612     @type solver: C{float}
613     @raise ValueException: if tolerance is not positive.
614 jgs 102 """
615     if not tol>0:
616     raise ValueException,"Tolerance as to be positive"
617 jgs 148 if tol<self.getTolerance(): self.__invalidateSolution()
618     self.trace("New tolerance %e"%tol)
619 jgs 102 self.__tolerance=tol
620     return
621 jgs 148
622 jgs 102 def getTolerance(self):
623     """
624 jgs 122 returns the tolerance set for the solution
625 jgs 148
626     @return: tolerance currently used.
627     @rtype: C{float}
628 jgs 102 """
629     return self.__tolerance
630    
631 jgs 148 # =============================================================================
632     # symmetry flag:
633     # =============================================================================
634 jgs 102 def isSymmetric(self):
635     """
636 jgs 148 checks if symmetry is indicated.
637    
638     @return : True is a symmetric PDE is indicated, otherwise False is returned
639     @rtype : C{bool}
640 jgs 102 """
641     return self.__sym
642    
643     def setSymmetryOn(self):
644     """
645 jgs 148 sets the symmetry flag.
646 jgs 102 """
647     if not self.isSymmetric():
648 jgs 148 self.trace("PDE is set to be symmetric")
649 jgs 102 self.__sym=True
650     self.__checkMatrixType()
651    
652     def setSymmetryOff(self):
653     """
654 jgs 148 removes the symmetry flag.
655 jgs 102 """
656     if self.isSymmetric():
657 jgs 148 self.trace("PDE is set to be unsymmetric")
658 jgs 102 self.__sym=False
659     self.__checkMatrixType()
660    
661     def setSymmetryTo(self,flag=False):
662 jgs 148 """
663     sets the symmetry flag to flag
664    
665     @param flag: If flag, the symmetry flag is set otherwise the symmetry flag is released.
666     @type flag: C{bool}
667     """
668     if flag:
669     self.setSymmetryOn()
670     else:
671     self.setSymmetryOff()
672 jgs 102
673 jgs 148
674     # =============================================================================
675     # function space handling for the equation as well as the solution
676     # =============================================================================
677 jgs 102 def setReducedOrderOn(self):
678     """
679 jgs 148 switches on reduced order for solution and equation representation
680 jgs 102 """
681     self.setReducedOrderForSolutionOn()
682     self.setReducedOrderForEquationOn()
683    
684     def setReducedOrderOff(self):
685     """
686 jgs 148 switches off reduced order for solution and equation representation
687 jgs 102 """
688     self.setReducedOrderForSolutionOff()
689     self.setReducedOrderForEquationOff()
690    
691     def setReducedOrderTo(self,flag=False):
692     """
693 jgs 148 sets order reduction for both solution and equation representation according to flag.
694 jgs 102
695 jgs 148 @param flag: if flag is True, the order reduction is switched on for both solution and equation representation, otherwise or
696     if flag is not present order reduction is switched off
697     @type flag: C{bool}
698 jgs 102 """
699     self.setReducedOrderForSolutionTo(flag)
700     self.setReducedOrderForEquationTo(flag)
701    
702 jgs 148
703 jgs 102 def setReducedOrderForSolutionOn(self):
704     """
705 jgs 148 switches on reduced order for solution representation
706 jgs 102 """
707     new_fs=escript.ReducedSolution(self.getDomain())
708     if self.getFunctionSpaceForSolution()!=new_fs:
709 jgs 148 self.trace("Reduced order is used to solution representation.")
710 jgs 102 self.__column_function_space=new_fs
711 jgs 148 self.__resetSystem()
712 jgs 102
713     def setReducedOrderForSolutionOff(self):
714     """
715 jgs 148 switches off reduced order for solution representation
716 jgs 102 """
717     new_fs=escript.Solution(self.getDomain())
718     if self.getFunctionSpaceForSolution()!=new_fs:
719 jgs 148 self.trace("Full order is used to interpolate solution.")
720 jgs 102 self.__column_function_space=new_fs
721 jgs 148 self.__resetSystem()
722 jgs 102
723     def setReducedOrderForSolutionTo(self,flag=False):
724     """
725 jgs 122 sets order for test functions according to flag
726 jgs 102
727 jgs 148 @param flag: if flag is True, the order reduction is switched on for solution representation, otherwise or
728     if flag is not present order reduction is switched off
729     @type flag: C{bool}
730 jgs 102 """
731     if flag:
732     self.setReducedOrderForSolutionOn()
733     else:
734     self.setReducedOrderForSolutionOff()
735 jgs 148
736 jgs 102 def setReducedOrderForEquationOn(self):
737     """
738 jgs 148 switches on reduced order for equation representation
739 jgs 102 """
740     new_fs=escript.ReducedSolution(self.getDomain())
741     if self.getFunctionSpaceForEquation()!=new_fs:
742 jgs 148 self.trace("Reduced order is used for test functions.")
743 jgs 102 self.__row_function_space=new_fs
744 jgs 148 self.__resetSystem()
745 jgs 102
746     def setReducedOrderForEquationOff(self):
747     """
748 jgs 148 switches off reduced order for equation representation
749 jgs 102 """
750     new_fs=escript.Solution(self.getDomain())
751     if self.getFunctionSpaceForEquation()!=new_fs:
752 jgs 148 self.trace("Full order is used for test functions.")
753 jgs 102 self.__row_function_space=new_fs
754 jgs 148 self.__resetSystem()
755 jgs 102
756     def setReducedOrderForEquationTo(self,flag=False):
757     """
758 jgs 122 sets order for test functions according to flag
759 jgs 102
760 jgs 148 @param flag: if flag is True, the order reduction is switched on for equation representation, otherwise or
761     if flag is not present order reduction is switched off
762     @type flag: C{bool}
763 jgs 102 """
764     if flag:
765     self.setReducedOrderForEquationOn()
766     else:
767     self.setReducedOrderForEquationOff()
768 jgs 148
769     # =============================================================================
770     # private method:
771     # =============================================================================
772     def __checkMatrixType(self):
773     """
774     reassess the matrix type and, if a new matrix is needed, resets the system.
775     """
776     new_matrix_type=self.getDomain().getSystemMatrixTypeId(self.getSolverMethod(),self.isSymmetric())
777     if not new_matrix_type==self.__matrix_type:
778     self.trace("Matrix type is now %d."%new_matrix_type)
779     self.__matrix_type=new_matrix_type
780     self.__resetSystem()
781     #
782     # rebuild switches :
783     #
784     def __invalidateSolution(self):
785     """
786     indicates the PDE has to be resolved if the solution is requested
787     """
788     if self.__solution_isValid: self.trace("PDE has to be resolved.")
789     self.__solution_isValid=False
790    
791     def __invalidateOperator(self):
792     """
793     indicates the operator has to be rebuilt next time it is used
794     """
795     if self.__operator_isValid: self.trace("Operator has to be rebuilt.")
796     self.__invalidateSolution()
797     self.__operator_isValid=False
798    
799     def __invalidateRightHandSide(self):
800     """
801     indicates the right hand side has to be rebuild next time it is used
802     """
803     if self.__righthandside_isValid: self.trace("Right hand side has to be rebuilt.")
804     self.__invalidateSolution()
805     self.__righthandside_isValid=False
806    
807     def __invalidateSystem(self):
808     """
809     annonced that everthing has to be rebuild:
810     """
811     if self.__righthandside_isValid: self.trace("System has to be rebuilt.")
812     self.__invalidateSolution()
813     self.__invalidateOperator()
814     self.__invalidateRightHandSide()
815    
816     def __resetSystem(self):
817     """
818     annonced that everthing has to be rebuild:
819     """
820     self.trace("New System is built from scratch.")
821     self.__operator=escript.Operator()
822     self.__operator_isValid=False
823     self.__righthandside=escript.Data()
824     self.__righthandside_isValid=False
825     self.__solution=escript.Data()
826     self.__solution_isValid=False
827     #
828     # system initialization:
829     #
830 jgs 121 def __getNewOperator(self):
831 jgs 102 """
832 jgs 148 returns an instance of a new operator
833 jgs 102 """
834 jgs 148 self.trace("New operator is allocated.")
835 jgs 102 return self.getDomain().newOperator( \
836     self.getNumEquations(), \
837     self.getFunctionSpaceForEquation(), \
838     self.getNumSolutions(), \
839     self.getFunctionSpaceForSolution(), \
840     self.__matrix_type)
841    
842 jgs 148 def __getNewRightHandSide(self):
843 jgs 102 """
844 jgs 148 returns an instance of a new right hand side
845 jgs 102 """
846 jgs 148 self.trace("New right hand side is allocated.")
847 jgs 121 if self.getNumEquations()>1:
848 jgs 148 return escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)
849 jgs 121 else:
850 jgs 148 return escript.Data(0.,(),self.getFunctionSpaceForEquation(),True)
851 jgs 102
852 jgs 121 def __getNewSolution(self):
853 jgs 102 """
854 jgs 148 returns an instance of a new solution
855 jgs 102 """
856 jgs 148 self.trace("New solution is allocated.")
857 jgs 121 if self.getNumSolutions()>1:
858     return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)
859     else:
860     return escript.Data(0.,(),self.getFunctionSpaceForSolution(),True)
861 jgs 102
862 jgs 148 def __makeFreshSolution(self):
863     """
864     makes sure that the solution is instantiated and returns it initialized by zeros
865     """
866     if self.__solution.isEmpty():
867     self.__solution=self.__getNewSolution()
868     else:
869     self.__solution*=0
870     self.trace("Solution is reset to zero.")
871     return self.__solution
872    
873     def __makeFreshRightHandSide(self):
874     """
875     makes sure that the right hand side is instantiated and returns it initialized by zeros
876     """
877     if self.__righthandside.isEmpty():
878     self.__righthandside=self.__getNewRightHandSide()
879     else:
880     self.__righthandside*=0
881     self.trace("Right hand side is reset to zero.")
882     return self.__righthandside
883    
884 jgs 121 def __makeFreshOperator(self):
885 jgs 102 """
886 jgs 148 makes sure that the operator is instantiated and returns it initialized by zeros
887 jgs 102 """
888     if self.__operator.isEmpty():
889 jgs 121 self.__operator=self.__getNewOperator()
890 jgs 102 else:
891     self.__operator.setValue(0.)
892 jgs 148 self.trace("Operator reset to zero")
893 jgs 102 return self.__operator
894    
895 jgs 148 def __applyConstraint(self):
896     """
897     applies the constraints defined by q and r to the system
898     """
899     if not self.isUsingLumping():
900     q=self.getCoefficientOfGeneralPDE("q")
901     r=self.getCoefficientOfGeneralPDE("r")
902     if not q.isEmpty() and not self.__operator.isEmpty():
903     # q is the row and column mask to indicate where constraints are set:
904     row_q=escript.Data(q,self.getFunctionSpaceForEquation())
905     col_q=escript.Data(q,self.getFunctionSpaceForSolution())
906     u=self.__getNewSolution()
907     if r.isEmpty():
908     r_s=self.__getNewSolution()
909     else:
910     r_s=escript.Data(r,self.getFunctionSpaceForSolution())
911     u.copyWithMask(r_s,col_q)
912     if not self.__righthandside.isEmpty():
913     self.__righthandside-=self.__operator*u
914     self.__righthandside=self.copyConstraint(self.__righthandside)
915     self.__operator.nullifyRowsAndCols(row_q,col_q,1.)
916     # =============================================================================
917     # function giving access to coefficients of the general PDE:
918     # =============================================================================
919     def getCoefficientOfGeneralPDE(self,name):
920 jgs 108 """
921 jgs 148 return the value of the coefficient name of the general PDE.
922    
923     @note This method is called by the assembling routine it can be overwritten
924     to map coefficients of a particular PDE to the general PDE.
925    
926     @param name: name of the coefficient requested.
927     @type name: C{string}
928     @return : the value of the coefficient name
929     @rtype : L{escript.Data}
930     @raise IllegalCoefficient: if name is not one of coefficients
931     "A", "B", "C", "D", "X", "Y", "d", "y", "d_contact", "y_contact", "r" or "q".
932 jgs 108 """
933 jgs 148 if self.hasCoefficientOfGeneralPDE(name):
934     return self.getCoefficient(name)
935     else:
936     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
937 jgs 108
938 jgs 148 def hasCoefficientOfGeneralPDE(self,name):
939 jgs 108 """
940 jgs 148 checks if name is a the name of a coefficient of the general PDE.
941    
942     @param name: name of the coefficient enquired.
943     @type name: C{string}
944     @return : True if name is the name of a coefficient of the general PDE. Otherwise False.
945     @rtype : C{bool}
946    
947 jgs 108 """
948 jgs 148 return self.__COEFFICIENTS_OF_GENEARL_PDE.has_key(name)
949 jgs 108
950 jgs 148 def createCoefficientOfGeneralPDE(self,name):
951 jgs 108 """
952 jgs 148 returns a new instance of a coefficient for coefficient name of the general PDE
953    
954     @param name: name of the coefficient requested.
955     @type name: C{string}
956     @return : a coefficient name initialized to 0.
957     @rtype : L{escript.Data}
958     @raise IllegalCoefficient: if name is not one of coefficients
959     "A", "B", "C", "D", "X", "Y", "d", "y", "d_contact", "y_contact", "r" or "q".
960 jgs 108 """
961 jgs 148 if self.hasCoefficientOfGeneralPDE(name):
962     return escript.Data(0,self.getShapeOfCoefficientOfGeneralPDE(name),self.getFunctionSpaceForCoefficientOfGeneralPDE(name))
963 jgs 108 else:
964 jgs 148 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
965 jgs 108
966 jgs 148 def getFunctionSpaceForCoefficientOfGeneralPDE(self,name):
967 jgs 108 """
968 jgs 148 return the L{escript.FunctionSpace} to be used for coefficient name of the general PDE
969    
970     @param name: name of the coefficient enquired.
971     @type name: C{string}
972     @return : the function space to be used for coefficient name
973     @rtype : L{escript.FunctionSpace}
974     @raise IllegalCoefficient: if name is not one of coefficients
975     "A", "B", "C", "D", "X", "Y", "d", "y", "d_contact", "y_contact", "r" or "q".
976 jgs 108 """
977 jgs 148 if self.hasCoefficientOfGeneralPDE(name):
978     return self.__COEFFICIENTS_OF_GENEARL_PDE[name].getFunctionSpace(self.getDomain())
979 jgs 108 else:
980 jgs 148 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
981 jgs 108
982 jgs 148 def getShapeOfCoefficientOfGeneralPDE(self,name):
983     """
984     return the shape of the coefficient name of the general PDE
985 jgs 108
986 jgs 148 @param name: name of the coefficient enquired.
987     @type name: C{string}
988     @return : the shape of the coefficient name
989     @rtype : C{tuple} of C{int}
990     @raise IllegalCoefficient: if name is not one of coefficients
991     "A", "B", "C", "D", "X", "Y", "d", "y", "d_contact", "y_contact", "r" or "q".
992 jgs 108
993 jgs 148 """
994     if self.hasCoefficientOfGeneralPDE(name):
995     return self.__COEFFICIENTS_OF_GENEARL_PDE[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())
996     else:
997     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
998 jgs 108
999 jgs 148 # =============================================================================
1000     # functions giving access to coefficients of a particular PDE implementation:
1001     # =============================================================================
1002     def getCoefficient(self,name):
1003     """
1004     returns the value of the coefficient name
1005 jgs 108
1006 jgs 148 @param name: name of the coefficient requested.
1007     @type name: C{string}
1008     @return : the value of the coefficient name
1009     @rtype : L{escript.Data}
1010     @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1011     """
1012     if self.hasCoefficient(name):
1013     return self.COEFFICIENTS[name].getValue()
1014     else:
1015     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1016 jgs 108
1017 jgs 148 def hasCoefficient(self,name):
1018     """
1019     return True if name is the name of a coefficient
1020 jgs 108
1021 jgs 148 @param name: name of the coefficient enquired.
1022     @type name: C{string}
1023     @return : True if name is the name of a coefficient of the general PDE. Otherwise False.
1024     @rtype : C{bool}
1025     """
1026     return self.COEFFICIENTS.has_key(name)
1027 jgs 108
1028 jgs 148 def createCoefficient(self, name):
1029     """
1030     create a L{escript.Data} object corresponding to coefficient name
1031 jgs 108
1032 jgs 148 @return : a coefficient name initialized to 0.
1033     @rtype : L{escript.Data}
1034     @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1035     """
1036     if self.hasCoefficient(name):
1037     return escript.Data(0.,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))
1038     else:
1039     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1040    
1041     def getFunctionSpaceForCoefficient(self,name):
1042     """
1043     return the L{escript.FunctionSpace} to be used for coefficient name
1044    
1045     @param name: name of the coefficient enquired.
1046     @type name: C{string}
1047     @return : the function space to be used for coefficient name
1048     @rtype : L{escript.FunctionSpace}
1049     @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1050     """
1051     if self.hasCoefficient(name):
1052     return self.COEFFICIENTS[name].getFunctionSpace(self.getDomain())
1053     else:
1054     raise ValueError,"unknown coefficient %s requested"%name
1055    
1056     def getShapeOfCoefficient(self,name):
1057     """
1058     return the shape of the coefficient name
1059    
1060     @param name: name of the coefficient enquired.
1061     @type name: C{string}
1062     @return : the shape of the coefficient name
1063     @rtype : C{tuple} of C{int}
1064     @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1065     """
1066     if self.hasCoefficient(name):
1067     return self.COEFFICIENTS[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())
1068     else:
1069     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1070    
1071     def resetCoefficients(self):
1072     """
1073     resets all coefficients to there default values.
1074     """
1075     for i in self.COEFFICIENTS.iterkeys():
1076     self.COEFFICIENTS[i].resetValue()
1077    
1078     def alteredCoefficient(self,name):
1079     """
1080     announce that coefficient name has been changed
1081    
1082     @param name: name of the coefficient enquired.
1083     @type name: C{string}
1084     @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1085     """
1086     if self.hasCoefficient(name):
1087     self.trace("Coefficient %s has been altered."%name)
1088     if self.COEFFICIENTS[name].isAlteringOperator(): self.__invalidateOperator()
1089     if self.COEFFICIENTS[name].isAlteringRightHandSide(): self.__invalidateRightHandSide()
1090     else:
1091     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1092    
1093     def copyConstraint(self,u):
1094 jgs 108 """
1095 jgs 148 copies the constraint into u and returns u.
1096    
1097     @param u: a function of rank 0 is a single PDE is solved and of shape (numSolution,) for a system of PDEs
1098     @type u: L{escript.Data}
1099     @return : the input u modified by the constraints.
1100     @rtype : L{escript.Data}
1101     @warning: u is altered if it has the appropriate L{escript.FunctionSpace}
1102 jgs 108
1103     """
1104 jgs 148 q=self.getCoefficientOfGeneralPDE("q")
1105     r=self.getCoefficientOfGeneralPDE("r")
1106     if not q.isEmpty():
1107     if u.isEmpty(): u=escript.Data(0.,q.getShape(),q.getFunctionSpace())
1108     if r.isEmpty():
1109     r=escript.Data(0,u.getShape(),u.getFunctionSpace())
1110     else:
1111     r=escript.Data(r,u.getFunctionSpace())
1112     u.copyWithMask(r,escript.Data(q,u.getFunctionSpace()))
1113     return u
1114    
1115     def setValue(self,**coefficients):
1116     """
1117     sets new values to coefficients
1118    
1119     @note This method is called by the assembling routine it can be overwritten
1120     to map coefficients of a particular PDE to the general PDE.
1121    
1122     @param name: name of the coefficient requested.
1123     @type name: C{string}
1124     @keyword A: value for coefficient A.
1125     @type A: any type that can be interpreted as L{escript.Data} object on L{escript.Function}.
1126     @keyword B: value for coefficient B
1127     @type B: any type that can be interpreted as L{escript.Data} object on L{escript.Function}.
1128     @keyword C: value for coefficient C
1129     @type C: any type that can be interpreted as L{escript.Data} object on L{escript.Function}.
1130     @keyword D: value for coefficient D
1131     @type D: any type that can be interpreted as L{escript.Data} object on L{escript.Function}.
1132     @keyword X: value for coefficient X
1133     @type X: any type that can be interpreted as L{escript.Data} object on L{escript.Function}.
1134     @keyword Y: value for coefficient Y
1135     @type Y: any type that can be interpreted as L{escript.Data} object on L{escript.Function}.
1136     @keyword d: value for coefficient d
1137     @type d: any type that can be interpreted as L{escript.Data} object on L{escript.FunctionOnBoundary}.
1138     @keyword y: value for coefficient y
1139     @type y: any type that can be interpreted as L{escript.Data} object on L{escript.FunctionOnBoundary}.
1140     @keyword d_contact: value for coefficient d_contact
1141     @type d_contact: any type that can be interpreted as L{escript.Data} object on L{escript.FunctionOnContactOne}.
1142     or L{escript.FunctionOnContactZero}.
1143     @keyword y_contact: value for coefficient y_contact
1144     @type y_contact: any type that can be interpreted as L{escript.Data} object on L{escript.FunctionOnContactOne}.
1145     or L{escript.FunctionOnContactZero}.
1146     @keyword r: values prescribed to the solution at the locations of constraints
1147     @type r: any type that can be interpreted as L{escript.Data} object on L{escript.Solution} or L{escript.ReducedSolution}
1148     depending of reduced order is used for the solution.
1149     @keyword q: mask for location of constraints
1150     @type q: any type that can be interpreted as L{escript.Data} object on L{escript.Solution} or L{escript.ReducedSolution}
1151     depending of reduced order is used for the representation of the equation.
1152     @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1153    
1154     """
1155 jgs 108 # check if the coefficients are legal:
1156     for i in coefficients.iterkeys():
1157     if not self.hasCoefficient(i):
1158 jgs 148 raise IllegalCoefficient,"Attempt to set unknown coefficient %s"%i
1159 jgs 108 # if the number of unknowns or equations is still unknown we try to estimate them:
1160 jgs 148 if self.__numEquations==None or self.__numSolutions==None:
1161 jgs 108 for i,d in coefficients.iteritems():
1162     if hasattr(d,"shape"):
1163     s=d.shape
1164     elif hasattr(d,"getShape"):
1165     s=d.getShape()
1166     else:
1167     s=numarray.array(d).shape
1168     if s!=None:
1169     # get number of equations and number of unknowns:
1170 jgs 148 res=self.COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(self.getDomain(),s)
1171 jgs 108 if res==None:
1172 jgs 148 raise IllegalCoefficientValue,"Illegal shape %s of coefficient %s"%(s,i)
1173 jgs 108 else:
1174 jgs 148 if self.__numEquations==None: self.__numEquations=res[0]
1175     if self.__numSolutions==None: self.__numSolutions=res[1]
1176     if self.__numEquations==None: raise UndefinedPDEError,"unidententified number of equations"
1177     if self.__numSolutions==None: raise UndefinedPDEError,"unidententified number of solutions"
1178 jgs 108 # now we check the shape of the coefficient if numEquations and numSolutions are set:
1179     for i,d in coefficients.iteritems():
1180 jgs 148 try:
1181     self.COEFFICIENTS[i].setValue(self.getDomain(),self.getNumEquations(),self.getNumSolutions(),d)
1182     except IllegalCoefficientValue,m:
1183     raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
1184 jgs 108 self.alteredCoefficient(i)
1185    
1186 jgs 148 # check if the systrem is inhomogeneous:
1187     if len(coefficients)>0 and not self.isUsingLumping():
1188     q=self.getCoefficientOfGeneralPDE("q")
1189     r=self.getCoefficientOfGeneralPDE("r")
1190     homogeneous_constraint=True
1191     if not q.isEmpty() and not r.isEmpty():
1192     if util.Lsup(q*r)>=1.e-13*util.Lsup(r):
1193     self.trace("Inhomogeneous constraint detected.")
1194     self.__invalidateSystem()
1195 jgs 108
1196 jgs 102
1197 jgs 108 def getSystem(self):
1198     """
1199 jgs 122 return the operator and right hand side of the PDE
1200 jgs 108 """
1201     if not self.__operator_isValid or not self.__righthandside_isValid:
1202     if self.isUsingLumping():
1203     if not self.__operator_isValid:
1204 jgs 148 if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution(): raise TypeError,"Lumped matrix requires same order for equations and unknowns"
1205     if not self.getCoefficientOfGeneralPDE("A").isEmpty(): raise Warning,"Lumped matrix does not allow coefficient A"
1206     if not self.getCoefficientOfGeneralPDE("B").isEmpty(): raise Warning,"Lumped matrix does not allow coefficient B"
1207     if not self.getCoefficientOfGeneralPDE("C").isEmpty(): raise Warning,"Lumped matrix does not allow coefficient C"
1208 jgs 121 mat=self.__getNewOperator()
1209 jgs 108 self.getDomain().addPDEToSystem(mat,escript.Data(), \
1210 jgs 148 self.getCoefficientOfGeneralPDE("A"), \
1211     self.getCoefficientOfGeneralPDE("B"), \
1212     self.getCoefficientOfGeneralPDE("C"), \
1213     self.getCoefficientOfGeneralPDE("D"), \
1214 jgs 108 escript.Data(), \
1215     escript.Data(), \
1216 jgs 148 self.getCoefficientOfGeneralPDE("d"), \
1217 jgs 108 escript.Data(),\
1218 jgs 148 self.getCoefficientOfGeneralPDE("d_contact"), \
1219 jgs 108 escript.Data())
1220 jgs 148 self.__operator=1./(mat*escript.Data(1,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True))
1221     del mat
1222     self.trace("New lumped operator has been built.")
1223 jgs 108 self.__operator_isValid=True
1224     if not self.__righthandside_isValid:
1225 jgs 121 self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \
1226 jgs 148 self.getCoefficientOfGeneralPDE("X"), \
1227     self.getCoefficientOfGeneralPDE("Y"),\
1228     self.getCoefficientOfGeneralPDE("y"),\
1229     self.getCoefficientOfGeneralPDE("y_contact"))
1230     self.trace("New right hand side as been built.")
1231 jgs 108 self.__righthandside_isValid=True
1232     else:
1233     if not self.__operator_isValid and not self.__righthandside_isValid:
1234 jgs 121 self.getDomain().addPDEToSystem(self.__makeFreshOperator(),self.__makeFreshRightHandSide(), \
1235 jgs 148 self.getCoefficientOfGeneralPDE("A"), \
1236     self.getCoefficientOfGeneralPDE("B"), \
1237     self.getCoefficientOfGeneralPDE("C"), \
1238     self.getCoefficientOfGeneralPDE("D"), \
1239     self.getCoefficientOfGeneralPDE("X"), \
1240     self.getCoefficientOfGeneralPDE("Y"), \
1241     self.getCoefficientOfGeneralPDE("d"), \
1242     self.getCoefficientOfGeneralPDE("y"), \
1243     self.getCoefficientOfGeneralPDE("d_contact"), \
1244     self.getCoefficientOfGeneralPDE("y_contact"))
1245 jgs 108 self.__applyConstraint()
1246 jgs 148 self.__righthandside=self.copyConstraint(self.__righthandside)
1247     self.trace("New system has been built.")
1248 jgs 108 self.__operator_isValid=True
1249     self.__righthandside_isValid=True
1250     elif not self.__righthandside_isValid:
1251 jgs 121 self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \
1252 jgs 148 self.getCoefficientOfGeneralPDE("X"), \
1253     self.getCoefficientOfGeneralPDE("Y"),\
1254     self.getCoefficientOfGeneralPDE("y"),\
1255     self.getCoefficientOfGeneralPDE("y_contact"))
1256     self.__righthandside=self.copyConstraint(self.__righthandside)
1257     self.trace("New right hand side has been built.")
1258 jgs 108 self.__righthandside_isValid=True
1259     elif not self.__operator_isValid:
1260 jgs 121 self.getDomain().addPDEToSystem(self.__makeFreshOperator(),escript.Data(), \
1261 jgs 148 self.getCoefficientOfGeneralPDE("A"), \
1262     self.getCoefficientOfGeneralPDE("B"), \
1263     self.getCoefficientOfGeneralPDE("C"), \
1264     self.getCoefficientOfGeneralPDE("D"), \
1265 jgs 108 escript.Data(), \
1266     escript.Data(), \
1267 jgs 148 self.getCoefficientOfGeneralPDE("d"), \
1268 jgs 108 escript.Data(),\
1269 jgs 148 self.getCoefficientOfGeneralPDE("d_contact"), \
1270 jgs 108 escript.Data())
1271     self.__applyConstraint()
1272 jgs 148 self.trace("New operator has been built.")
1273 jgs 108 self.__operator_isValid=True
1274     return (self.__operator,self.__righthandside)
1275 jgs 102
1276    
1277    
1278    
1279 jgs 108 class AdvectivePDE(LinearPDE):
1280     """
1281 jgs 122 Class to handle a linear PDE dominated by advective terms:
1282 jgs 148
1283 jgs 108 class to define a linear PDE of the form
1284 jgs 104
1285 jgs 122 \f[
1286     -(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
1287     \f]
1288 jgs 102
1289 jgs 122 with boundary conditons:
1290 jgs 102
1291 jgs 122 \f[
1292     n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_{ik}u_k = - n_j*X_{ij} + y_i
1293     \f]
1294 jgs 102
1295 jgs 122 and contact conditions
1296 jgs 102
1297 jgs 122 \f[
1298     n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d^{contact}_{ik}[u_k] = - n_j*X_{ij} + y^{contact}_{i}
1299     \f]
1300 jgs 102
1301 jgs 122 and constraints:
1302 jgs 102
1303 jgs 122 \f[
1304     u_i=r_i \quad \mathrm{where} \quad q_i>0
1305     \f]
1306 jgs 110 """
1307 jgs 148 def __init__(self,domain,numEquations=0,numSolutions=0,xi=None,debug=False):
1308     LinearPDE.__init__(self,domain,numEquations,numSolutions,debug)
1309     if xi==None:
1310     self.__xi=AdvectivePDE.ELMAN_RAMAGE
1311     else:
1312     self.__xi=xi
1313 jgs 110 self.__Xi=escript.Data()
1314 jgs 102
1315 jgs 110 def __calculateXi(self,peclet_factor,Z,h):
1316     Z_max=util.Lsup(Z)
1317     if Z_max>0.:
1318     return h*self.__xi(Z*peclet_factor)/(Z+Z_max*self.TOL)
1319     else:
1320     return 0.
1321 jgs 102
1322 jgs 110 def setValue(self,**args):
1323     if "A" in args.keys() or "B" in args.keys() or "C" in args.keys(): self.__Xi=escript.Data()
1324 jgs 148 LinearPDE.setValue(**args)
1325    
1326     def ELMAN_RAMAGE(P):
1327     """ """
1328     return (P-1.).wherePositive()*0.5*(1.-1./(P+1.e-15))
1329     def SIMPLIFIED_BROOK_HUGHES(P):
1330     """ """
1331     c=(P-3.).whereNegative()
1332     return P/6.*c+1./2.*(1.-c)
1333    
1334     def HALF(P):
1335     """ """
1336     return escript.Scalar(0.5,P.getFunctionSpace())
1337    
1338 jgs 110 def getXi(self):
1339     if self.__Xi.isEmpty():
1340     B=self.getCoefficient("B")
1341     C=self.getCoefficient("C")
1342     A=self.getCoefficient("A")
1343     h=self.getDomain().getSize()
1344     self.__Xi=escript.Scalar(0.,self.getFunctionSpaceForCoefficient("A"))
1345     if not C.isEmpty() or not B.isEmpty():
1346     if not C.isEmpty() and not B.isEmpty():
1347     Z2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))
1348     if self.getNumEquations()>1:
1349     if self.getNumSolutions()>1:
1350     for i in range(self.getNumEquations()):
1351     for k in range(self.getNumSolutions()):
1352     for l in range(self.getDim()): Z2+=(C[i,k,l]-B[i,l,k])**2
1353     else:
1354     for i in range(self.getNumEquations()):
1355     for l in range(self.getDim()): Z2+=(C[i,l]-B[i,l])**2
1356     else:
1357     if self.getNumSolutions()>1:
1358     for k in range(self.getNumSolutions()):
1359     for l in range(self.getDim()): Z2+=(C[k,l]-B[l,k])**2
1360     else:
1361     for l in range(self.getDim()): Z2+=(C[l]-B[l])**2
1362     length_of_Z=util.sqrt(Z2)
1363     elif C.isEmpty():
1364     length_of_Z=util.length(B)
1365     else:
1366     length_of_Z=util.length(C)
1367 jgs 102
1368 jgs 110 Z_max=util.Lsup(length_of_Z)
1369     if Z_max>0.:
1370     length_of_A=util.length(A)
1371     A_max=util.Lsup(length_of_A)
1372     if A_max>0:
1373     inv_A=1./(length_of_A+A_max*self.TOL)
1374     else:
1375     inv_A=1./self.TOL
1376     peclet_number=length_of_Z*h/2*inv_A
1377     xi=self.__xi(peclet_number)
1378     self.__Xi=h*xi/(length_of_Z+Z_max*self.TOL)
1379     print "@ preclet number = %e"%util.Lsup(peclet_number),util.Lsup(xi),util.Lsup(length_of_Z)
1380     return self.__Xi
1381 jgs 102
1382 jgs 148
1383     def getCoefficientOfGeneralPDE(self,name):
1384 jgs 108 """
1385 jgs 122 return the value of the coefficient name of the general PDE
1386    
1387     @param name:
1388 jgs 108 """
1389 jgs 110 if not self.getNumEquations() == self.getNumSolutions():
1390     raise ValueError,"AdvectivePDE expects the number of solution componets and the number of equations to be equal."
1391    
1392 jgs 148 if name == "A" :
1393 jgs 108 A=self.getCoefficient("A")
1394     B=self.getCoefficient("B")
1395     C=self.getCoefficient("C")
1396 jgs 110 if B.isEmpty() and C.isEmpty():
1397     Aout=A
1398     else:
1399 jgs 148 if A.isEmpty():
1400 jgs 110 Aout=self.createNewCoefficient("A")
1401     else:
1402     Aout=A[:]
1403     Xi=self.getXi()
1404     if self.getNumEquations()>1:
1405     for i in range(self.getNumEquations()):
1406     for j in range(self.getDim()):
1407 jgs 108 for k in range(self.getNumSolutions()):
1408 jgs 110 for l in range(self.getDim()):
1409     if not C.isEmpty() and not B.isEmpty():
1410     for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*(C[p,i,j]-B[p,j,i])*(C[p,k,l]-B[p,l,k])
1411     elif C.isEmpty():
1412     for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*B[p,j,i]*B[p,l,k]
1413     else:
1414     for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*C[p,i,j]*C[p,k,l]
1415     else:
1416     for j in range(self.getDim()):
1417     for l in range(self.getDim()):
1418     if not C.isEmpty() and not B.isEmpty():
1419     Aout[j,l]+=Xi*(C[j]-B[j])*(C[l]-B[l])
1420     elif C.isEmpty():
1421     Aout[j,l]+=Xi*B[j]*B[l]
1422     else:
1423     Aout[j,l]+=Xi*C[j]*C[l]
1424     return Aout
1425 jgs 148 elif name == "B" :
1426 jgs 110 B=self.getCoefficient("B")
1427     C=self.getCoefficient("C")
1428     D=self.getCoefficient("D")
1429     if C.isEmpty() or D.isEmpty():
1430     Bout=B
1431     else:
1432     Xi=self.getXi()
1433 jgs 148 if B.isEmpty():
1434 jgs 110 Bout=self.createNewCoefficient("B")
1435     else:
1436     Bout=B[:]
1437     if self.getNumEquations()>1:
1438     for k in range(self.getNumSolutions()):
1439 jgs 148 for p in range(self.getNumEquations()):
1440 jgs 110 tmp=Xi*D[p,k]
1441     for i in range(self.getNumEquations()):
1442     for j in range(self.getDim()):
1443     Bout[i,j,k]+=tmp*C[p,i,j]
1444     else:
1445     tmp=Xi*D
1446     for j in range(self.getDim()): Bout[j]+=tmp*C[j]
1447     return Bout
1448 jgs 148 elif name == "C" :
1449 jgs 110 B=self.getCoefficient("B")
1450     C=self.getCoefficient("C")
1451     D=self.getCoefficient("D")
1452     if B.isEmpty() or D.isEmpty():
1453     Cout=C
1454     else:
1455     Xi=self.getXi()
1456 jgs 148 if C.isEmpty():
1457 jgs 110 Cout=self.createNewCoefficient("C")
1458     else:
1459     Cout=C[:]
1460     if self.getNumEquations()>1:
1461     for k in range(self.getNumSolutions()):
1462     for p in range(self.getNumEquations()):
1463     tmp=Xi*D[p,k]
1464     for i in range(self.getNumEquations()):
1465     for l in range(self.getDim()):
1466     Cout[i,k,l]+=tmp*B[p,l,i]
1467     else:
1468     tmp=Xi*D
1469     for j in range(self.getDim()): Cout[j]+=tmp*B[j]
1470     return Cout
1471 jgs 148 elif name == "D" :
1472 jgs 108 return self.getCoefficient("D")
1473 jgs 148 elif name == "X" :
1474 jgs 110 X=self.getCoefficient("X")
1475     Y=self.getCoefficient("Y")
1476     B=self.getCoefficient("B")
1477     C=self.getCoefficient("C")
1478     if Y.isEmpty() or (B.isEmpty() and C.isEmpty()):
1479     Xout=X
1480     else:
1481     if X.isEmpty():
1482     Xout=self.createNewCoefficient("X")
1483     else:
1484     Xout=X[:]
1485     Xi=self.getXi()
1486     if self.getNumEquations()>1:
1487 jgs 148 for p in range(self.getNumEquations()):
1488 jgs 110 tmp=Xi*Y[p]
1489     for i in range(self.getNumEquations()):
1490     for j in range(self.getDim()):
1491     if not C.isEmpty() and not B.isEmpty():
1492     Xout[i,j]+=tmp*(C[p,i,j]-B[p,j,i])
1493     elif C.isEmpty():
1494     Xout[i,j]-=tmp*B[p,j,i]
1495     else:
1496     Xout[i,j]+=tmp*C[p,i,j]
1497     else:
1498     tmp=Xi*Y
1499     for j in range(self.getDim()):
1500     if not C.isEmpty() and not B.isEmpty():
1501     Xout[j]+=tmp*(C[j]-B[j])
1502     elif C.isEmpty():
1503     Xout[j]-=tmp*B[j]
1504     else:
1505     Xout[j]+=tmp*C[j]
1506     return Xout
1507 jgs 148 elif name == "Y" :
1508 jgs 108 return self.getCoefficient("Y")
1509 jgs 148 elif name == "d" :
1510 jgs 108 return self.getCoefficient("d")
1511 jgs 148 elif name == "y" :
1512 jgs 108 return self.getCoefficient("y")
1513 jgs 148 elif name == "d_contact" :
1514 jgs 108 return self.getCoefficient("d_contact")
1515     elif name == "y_contact" :
1516     return self.getCoefficient("y_contact")
1517 jgs 148 elif name == "r" :
1518 jgs 108 return self.getCoefficient("r")
1519 jgs 148 elif name == "q" :
1520 jgs 108 return self.getCoefficient("q")
1521     else:
1522     raise SystemError,"unknown PDE coefficient %s",name
1523    
1524    
1525 jgs 102 class Poisson(LinearPDE):
1526     """
1527 jgs 142 Class to define a Poisson equation problem:
1528 jgs 122
1529 jgs 102 class to define a linear PDE of the form
1530 jgs 122 \f[
1531     -u_{,jj} = f
1532     \f]
1533    
1534     with boundary conditons:
1535    
1536     \f[
1537     n_j*u_{,j} = 0
1538     \f]
1539    
1540     and constraints:
1541    
1542     \f[
1543     u=0 \quad \mathrm{where} \quad q>0
1544     \f]
1545 jgs 102 """
1546    
1547 jgs 148 def __init__(self,domain,f=escript.Data(),q=escript.Data(),debug=False):
1548     LinearPDE.__init__(self,domain,1,1,debug)
1549     self.COEFFICIENTS={"f": PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1550     "q": PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.BOTH)}
1551 jgs 102 self.setSymmetryOn()
1552     self.setValue(f,q)
1553    
1554     def setValue(self,f=escript.Data(),q=escript.Data()):
1555 jgs 142 """set value of PDE parameters f and q"""
1556 jgs 121 self._LinearPDE__setValue(f=f,q=q)
1557 jgs 108
1558 jgs 148 def getCoefficientOfGeneralPDE(self,name):
1559 jgs 108 """
1560 jgs 122 return the value of the coefficient name of the general PDE
1561    
1562     @param name:
1563 jgs 108 """
1564 jgs 148 if name == "A" :
1565 jgs 108 return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))
1566 jgs 148 elif name == "B" :
1567 jgs 108 return escript.Data()
1568 jgs 148 elif name == "C" :
1569 jgs 108 return escript.Data()
1570 jgs 148 elif name == "D" :
1571 jgs 108 return escript.Data()
1572 jgs 148 elif name == "X" :
1573 jgs 108 return escript.Data()
1574 jgs 148 elif name == "Y" :
1575 jgs 108 return self.getCoefficient("f")
1576 jgs 148 elif name == "d" :
1577 jgs 108 return escript.Data()
1578 jgs 148 elif name == "y" :
1579 jgs 108 return escript.Data()
1580 jgs 148 elif name == "d_contact" :
1581 jgs 108 return escript.Data()
1582     elif name == "y_contact" :
1583     return escript.Data()
1584 jgs 148 elif name == "r" :
1585 jgs 108 return escript.Data()
1586 jgs 148 elif name == "q" :
1587 jgs 108 return self.getCoefficient("q")
1588     else:
1589     raise SystemError,"unknown PDE coefficient %s",name
1590 jgs 121
1591 jgs 142 class LameEquation(LinearPDE):
1592     """
1593     Class to define a Lame equation problem:
1594    
1595     class to define a linear PDE of the form
1596     \f[
1597 jgs 147 -(\mu (u_{i,j}+u_{j,i}))_{,j} - \lambda u_{j,ji}} = F_i -\sigma_{ij,j}
1598 jgs 142 \f]
1599    
1600     with boundary conditons:
1601    
1602     \f[
1603 jgs 147 n_j(\mu (u_{i,j}+u_{j,i})-sigma_{ij}) + n_i\lambda u_{j,j} = f_i
1604 jgs 142 \f]
1605    
1606     and constraints:
1607    
1608     \f[
1609     u_i=r_i \quad \mathrm{where} \quad q_i>0
1610     \f]
1611     """
1612    
1613 jgs 148 def __init__(self,domain,f=escript.Data(),q=escript.Data(),debug=False):
1614     LinearPDE.__init__(self,domain,domain.getDim(),domain.getDim(),debug)
1615     self.COEFFICIENTS={ "lame_lambda" : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),
1616     "lame_mu" : PDECoefficient(PDECoefficient.INTERIOR,(),PDECoefficient.OPERATOR),
1617     "F" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1618     "sigma" : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM),PDECoefficient.RIGHTHANDSIDE),
1619     "f" : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),
1620     "r" : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.BOTH),
1621     "q" : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.BOTH)}
1622 jgs 142 self.setSymmetryOn()
1623    
1624     def setValue(self,lame_lambda=escript.Data(),lame_mu=escript.Data(),F=escript.Data(),sigma=escript.Data(),f=escript.Data(),r=escript.Data(),q=escript.Data()):
1625     """set value of PDE parameters"""
1626     self._LinearPDE__setValue(lame_lambda=lame_lambda, \
1627     lame_mu=lame_mu, \
1628     F=F, \
1629     sigma=sigma, \
1630     f=f, \
1631     r=r, \
1632 jgs 148 q=q)
1633     def getCoefficientOfGeneralPDE(self,name):
1634 jgs 142 """
1635     return the value of the coefficient name of the general PDE
1636    
1637     @param name:
1638     """
1639 jgs 148 if name == "A" :
1640     out =self.createCoefficientOfGeneralPDE("A")
1641 jgs 142 for i in range(self.getDim()):
1642     for j in range(self.getDim()):
1643 jgs 147 out[i,i,j,j] += self.getCoefficient("lame_lambda")
1644     out[i,j,j,i] += self.getCoefficient("lame_mu")
1645     out[i,j,i,j] += self.getCoefficient("lame_mu")
1646 jgs 142 return out
1647 jgs 148 elif name == "B" :
1648 jgs 142 return escript.Data()
1649 jgs 148 elif name == "C" :
1650 jgs 142 return escript.Data()
1651 jgs 148 elif name == "D" :
1652 jgs 142 return escript.Data()
1653 jgs 148 elif name == "X" :
1654 jgs 142 return self.getCoefficient("sigma")
1655 jgs 148 elif name == "Y" :
1656 jgs 142 return self.getCoefficient("F")
1657 jgs 148 elif name == "d" :
1658 jgs 142 return escript.Data()
1659 jgs 148 elif name == "y" :
1660 jgs 142 return self.getCoefficient("f")
1661 jgs 148 elif name == "d_contact" :
1662 jgs 142 return escript.Data()
1663     elif name == "y_contact" :
1664     return escript.Data()
1665 jgs 148 elif name == "r" :
1666 jgs 142 return self.getCoefficient("r")
1667 jgs 148 elif name == "q" :
1668 jgs 142 return self.getCoefficient("q")
1669     else:
1670     raise SystemError,"unknown PDE coefficient %s",name
1671    
1672 jgs 121 # $Log$
1673 jgs 148 # Revision 1.11 2005/08/23 01:24:28 jgs
1674     # Merge of development branch dev-02 back to main trunk on 2005-08-23
1675     #
1676 jgs 147 # Revision 1.10 2005/08/12 01:45:36 jgs
1677     # erge of development branch dev-02 back to main trunk on 2005-08-12
1678 jgs 142 #
1679 jgs 148 # Revision 1.9.2.4 2005/08/22 07:11:09 gross
1680     # some problems with LinearPDEs fixed.
1681     #
1682     # Revision 1.9.2.3 2005/08/18 04:48:48 gross
1683     # the methods SetLumping*() are removed. Lumping is set trough setSolverMethod(LinearPDE.LUMPING)
1684     #
1685     # Revision 1.9.2.2 2005/08/18 04:39:32 gross
1686     # the constants have been removed from util.py as they not needed anymore. PDE related constants are accessed through LinearPDE attributes now
1687     #
1688 jgs 147 # Revision 1.9.2.1 2005/07/29 07:10:27 gross
1689     # new functions in util and a new pde type in linearPDEs
1690 jgs 122 #
1691 jgs 147 # Revision 1.1.2.25 2005/07/28 04:21:09 gross
1692     # Lame equation: (linear elastic, isotropic) added
1693 jgs 121 #
1694 jgs 142 # Revision 1.1.2.24 2005/07/22 06:37:11 gross
1695     # some extensions to modellib and linearPDEs
1696     #
1697 jgs 122 # Revision 1.1.2.23 2005/05/13 00:55:20 cochrane
1698     # Fixed up some docstrings. Moved module-level functions to top of file so
1699     # that epydoc and doxygen can pick them up properly.
1700     #
1701     # Revision 1.1.2.22 2005/05/12 11:41:30 gross
1702     # some basic Models have been added
1703     #
1704     # Revision 1.1.2.21 2005/05/12 07:16:12 cochrane
1705     # Moved ELMAN_RAMAGE, SIMPLIFIED_BROOK_HUGHES, and HALF functions to bottom of
1706     # file so that the AdvectivePDE class is picked up by doxygen. Some
1707     # reformatting of docstrings. Addition of code to make equations come out
1708     # as proper LaTeX.
1709     #
1710 jgs 121 # Revision 1.1.2.20 2005/04/15 07:09:08 gross
1711     # some problems with functionspace and linearPDEs fixed.
1712     #
1713     # Revision 1.1.2.19 2005/03/04 05:27:07 gross
1714     # bug in SystemPattern fixed.
1715     #
1716     # Revision 1.1.2.18 2005/02/08 06:16:45 gross
1717     # Bugs in AdvectivePDE fixed, AdvectiveTest is stable but more testing is needed
1718     #
1719     # Revision 1.1.2.17 2005/02/08 05:56:19 gross
1720     # Reference Number handling added
1721     #
1722     # Revision 1.1.2.16 2005/02/07 04:41:28 gross
1723     # some function exposed to python to make mesh merging running
1724     #
1725     # Revision 1.1.2.15 2005/02/03 00:14:44 gross
1726     # timeseries add and ESySParameter.py renames esysXML.py for consistence
1727     #
1728     # Revision 1.1.2.14 2005/02/01 06:44:10 gross
1729     # new implementation of AdvectivePDE which now also updates right hand side. systems of PDEs are still not working
1730     #
1731     # Revision 1.1.2.13 2005/01/25 00:47:07 gross
1732     # updates in the documentation
1733     #
1734     # Revision 1.1.2.12 2005/01/12 01:28:04 matt
1735     # Added createCoefficient method for linearPDEs.
1736     #
1737     # Revision 1.1.2.11 2005/01/11 01:55:34 gross
1738     # a problem in linearPDE class fixed
1739     #
1740     # Revision 1.1.2.10 2005/01/07 01:13:29 gross
1741     # some bugs in linearPDE fixed
1742     #
1743     # Revision 1.1.2.9 2005/01/06 06:24:58 gross
1744     # some bugs in slicing fixed
1745     #
1746     # Revision 1.1.2.8 2005/01/05 04:21:40 gross
1747     # FunctionSpace checking/matchig in slicing added
1748     #
1749     # Revision 1.1.2.7 2004/12/29 10:03:41 gross
1750     # bug in setValue fixed
1751     #
1752     # Revision 1.1.2.6 2004/12/29 05:29:59 gross
1753     # AdvectivePDE successfully tested for Peclet number 1000000. there is still a problem with setValue and Data()
1754     #
1755     # Revision 1.1.2.5 2004/12/29 00:18:41 gross
1756     # AdvectivePDE added
1757     #
1758     # Revision 1.1.2.4 2004/12/24 06:05:41 gross
1759     # some changes in linearPDEs to add AdevectivePDE
1760     #
1761     # Revision 1.1.2.3 2004/12/16 00:12:34 gross
1762     # __init__ of LinearPDE does not accept any coefficient anymore
1763     #
1764     # Revision 1.1.2.2 2004/12/14 03:55:01 jgs
1765     # *** empty log message ***
1766     #
1767     # Revision 1.1.2.1 2004/12/12 22:53:47 gross
1768     # linearPDE has been renamed LinearPDE
1769     #
1770     # Revision 1.1.1.1.2.7 2004/12/07 10:13:08 gross
1771     # GMRES added
1772     #
1773     # Revision 1.1.1.1.2.6 2004/12/07 03:19:50 gross
1774     # options for GMRES and PRES20 added
1775     #
1776     # Revision 1.1.1.1.2.5 2004/12/01 06:25:15 gross
1777     # some small changes
1778     #
1779     # Revision 1.1.1.1.2.4 2004/11/24 01:50:21 gross
1780     # Finley solves 4M unknowns now
1781     #
1782     # Revision 1.1.1.1.2.3 2004/11/15 06:05:26 gross
1783     # poisson solver added
1784     #
1785     # Revision 1.1.1.1.2.2 2004/11/12 06:58:15 gross
1786     # 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
1787     #
1788     # Revision 1.1.1.1.2.1 2004/10/28 22:59:22 gross
1789     # finley's RecTest.py is running now: problem in SystemMatrixAdapater fixed
1790     #
1791     # Revision 1.1.1.1 2004/10/26 06:53:56 jgs
1792     # initial import of project esys2
1793     #
1794     # Revision 1.3.2.3 2004/10/26 06:43:48 jgs
1795     # committing Lutz's and Paul's changes to brach jgs
1796     #
1797     # Revision 1.3.4.1 2004/10/20 05:32:51 cochrane
1798     # Added incomplete Doxygen comments to files, or merely put the docstrings that already exist into Doxygen form.
1799     #
1800     # Revision 1.3 2004/09/23 00:53:23 jgs
1801     # minor fixes
1802     #
1803     # Revision 1.1 2004/08/28 12:58:06 gross
1804     # SimpleSolve is not running yet: problem with == of functionsspace
1805     #
1806     #

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26