/[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 2100 - (hide annotations)
Wed Nov 26 08:13:00 2008 UTC (10 years, 9 months ago) by gross
File MIME type: text/x-python
File size: 118017 byte(s)
This commit cleans up the incompressible solver and adds a DarcyFlux solver in model module. 
Some documentation for both classes has been added.
The convection code is only linear at the moment.



1 ksteube 1809
2     ########################################################
3 ksteube 1312 #
4 ksteube 1809 # Copyright (c) 2003-2008 by University of Queensland
5     # Earth Systems Science Computational Center (ESSCC)
6     # http://www.uq.edu.au/esscc
7 ksteube 1312 #
8 ksteube 1809 # Primary Business: Queensland, Australia
9     # Licensed under the Open Software License version 3.0
10     # http://www.opensource.org/licenses/osl-3.0.php
11 ksteube 1312 #
12 ksteube 1809 ########################################################
13 ksteube 1312
14 ksteube 1809 __copyright__="""Copyright (c) 2003-2008 by University of Queensland
15     Earth Systems Science Computational Center (ESSCC)
16     http://www.uq.edu.au/esscc
17     Primary Business: Queensland, Australia"""
18     __license__="""Licensed under the Open Software License version 3.0
19     http://www.opensource.org/licenses/osl-3.0.php"""
20     __url__="http://www.uq.edu.au/esscc/escript-finley"
21    
22 jgs 149 """
23 jgs 151 The module provides an interface to define and solve linear partial
24 gross 1841 differential equations (PDEs) and Transport problems within L{escript}. L{linearPDEs} does not provide any
25 jgs 151 solver capabilities in itself but hands the PDE over to
26 jgs 149 the PDE solver library defined through the L{Domain<escript.Domain>} of the PDE.
27 gross 1841 The general interface is provided through the L{LinearPDE} class.
28     L{TransportProblem} provides an interface to initial value problems dominated by its advective terms.
29 jgs 102
30 jgs 149 @var __author__: name of author
31 gross 637 @var __copyright__: copyrights
32 elspeth 614 @var __license__: licence agreement
33 jgs 149 @var __url__: url entry point on documentation
34     @var __version__: version
35     @var __date__: date of the version
36 jgs 102 """
37    
38 gross 1417 import math
39 jgs 102 import escript
40     import util
41     import numarray
42    
43 jgs 149 __author__="Lutz Gross, l.gross@uq.edu.au"
44 jgs 102
45 jgs 149
46 jgs 148 class IllegalCoefficient(ValueError):
47 jgs 102 """
48 jgs 148 raised if an illegal coefficient of the general ar particular PDE is requested.
49     """
50 gross 1072 pass
51 jgs 102
52 jgs 148 class IllegalCoefficientValue(ValueError):
53 jgs 102 """
54 jgs 148 raised if an incorrect value for a coefficient is used.
55     """
56 gross 1072 pass
57 jgs 122
58 gross 1072 class IllegalCoefficientFunctionSpace(ValueError):
59     """
60     raised if an incorrect function space for a coefficient is used.
61     """
62    
63 jgs 148 class UndefinedPDEError(ValueError):
64     """
65     raised if a PDE is not fully defined yet.
66     """
67 gross 1072 pass
68 jgs 102
69 gross 1841 class PDECoef(object):
70 jgs 102 """
71 jgs 149 A class for describing a PDE coefficient
72    
73     @cvar INTERIOR: indicator that coefficient is defined on the interior of the PDE domain
74     @cvar BOUNDARY: indicator that coefficient is defined on the boundary of the PDE domain
75     @cvar CONTACT: indicator that coefficient is defined on the contact region within the PDE domain
76 gross 1072 @cvar INTERIOR_REDUCED: indicator that coefficient is defined on the interior of the PDE domain using a reduced integration order
77     @cvar BOUNDARY_REDUCED: indicator that coefficient is defined on the boundary of the PDE domain using a reduced integration order
78     @cvar CONTACT_REDUCED: indicator that coefficient is defined on the contact region within the PDE domain using a reduced integration order
79 jgs 149 @cvar SOLUTION: indicator that coefficient is defined trough a solution of the PDE
80 jgs 150 @cvar REDUCED: indicator that coefficient is defined trough a reduced solution of the PDE
81 jgs 149 @cvar BY_EQUATION: indicator that the dimension of the coefficient shape is defined by the number PDE equations
82     @cvar BY_SOLUTION: indicator that the dimension of the coefficient shape is defined by the number PDE solutions
83     @cvar BY_DIM: indicator that the dimension of the coefficient shape is defined by the spatial dimension
84     @cvar OPERATOR: indicator that the the coefficient alters the operator of the PDE
85     @cvar RIGHTHANDSIDE: indicator that the the coefficient alters the right hand side of the PDE
86     @cvar BOTH: indicator that the the coefficient alters the operator as well as the right hand side of the PDE
87    
88 jgs 102 """
89     INTERIOR=0
90     BOUNDARY=1
91     CONTACT=2
92 jgs 149 SOLUTION=3
93 jgs 150 REDUCED=4
94 jgs 149 BY_EQUATION=5
95     BY_SOLUTION=6
96     BY_DIM=7
97     OPERATOR=10
98     RIGHTHANDSIDE=11
99     BOTH=12
100 gross 1072 INTERIOR_REDUCED=13
101     BOUNDARY_REDUCED=14
102     CONTACT_REDUCED=15
103 jgs 149
104 gross 1072 def __init__(self, where, pattern, altering):
105 jgs 102 """
106 jgs 122 Initialise a PDE Coefficient type
107 jgs 151
108 jgs 149 @param where: describes where the coefficient lives
109 gross 1072 @type where: one of L{INTERIOR}, L{BOUNDARY}, L{CONTACT}, L{SOLUTION}, L{REDUCED},
110     L{INTERIOR_REDUCED}, L{BOUNDARY_REDUCED}, L{CONTACT_REDUCED}.
111 jgs 151 @param pattern: describes the shape of the coefficient and how the shape is build for a given
112 jgs 149 spatial dimension and numbers of equation and solution in then PDE. For instance,
113     (L{BY_EQUATION},L{BY_SOLUTION},L{BY_DIM}) descrbes a rank 3 coefficient which
114     is instanciated as shape (3,2,2) in case of a three equations and two solution components
115     on a 2-dimensional domain. In the case of single equation and a single solution component
116     the shape compoments marked by L{BY_EQUATION} or L{BY_SOLUTION} are dropped. In this case
117     the example would be read as (2,).
118     @type pattern: C{tuple} of L{BY_EQUATION}, L{BY_SOLUTION}, L{BY_DIM}
119     @param altering: indicates what part of the PDE is altered if the coefficiennt is altered
120     @type altering: one of L{OPERATOR}, L{RIGHTHANDSIDE}, L{BOTH}
121 jgs 102 """
122 gross 1841 super(PDECoef, self).__init__()
123 jgs 102 self.what=where
124     self.pattern=pattern
125     self.altering=altering
126 jgs 108 self.resetValue()
127 jgs 102
128 jgs 108 def resetValue(self):
129     """
130 jgs 122 resets coefficient value to default
131 jgs 108 """
132     self.value=escript.Data()
133    
134 jgs 149 def getFunctionSpace(self,domain,reducedEquationOrder=False,reducedSolutionOrder=False):
135 jgs 102 """
136 jgs 149 defines the L{FunctionSpace<escript.FunctionSpace>} of the coefficient
137 jgs 102
138 jgs 149 @param domain: domain on which the PDE uses the coefficient
139     @type domain: L{Domain<escript.Domain>}
140     @param reducedEquationOrder: True to indicate that reduced order is used to represent the equation
141 gross 720 @type reducedEquationOrder: C{bool}
142 jgs 149 @param reducedSolutionOrder: True to indicate that reduced order is used to represent the solution
143 gross 720 @type reducedSolutionOrder: C{bool}
144 jgs 149 @return: L{FunctionSpace<escript.FunctionSpace>} of the coefficient
145     @rtype: L{FunctionSpace<escript.FunctionSpace>}
146 jgs 102 """
147 jgs 151 if self.what==self.INTERIOR:
148 jgs 149 return escript.Function(domain)
149 gross 1072 elif self.what==self.INTERIOR_REDUCED:
150     return escript.ReducedFunction(domain)
151 jgs 151 elif self.what==self.BOUNDARY:
152 jgs 149 return escript.FunctionOnBoundary(domain)
153 gross 1072 elif self.what==self.BOUNDARY_REDUCED:
154     return escript.ReducedFunctionOnBoundary(domain)
155 jgs 151 elif self.what==self.CONTACT:
156 jgs 149 return escript.FunctionOnContactZero(domain)
157 gross 1072 elif self.what==self.CONTACT_REDUCED:
158     return escript.ReducedFunctionOnContactZero(domain)
159 jgs 151 elif self.what==self.SOLUTION:
160 jgs 149 if reducedEquationOrder and reducedSolutionOrder:
161     return escript.ReducedSolution(domain)
162     else:
163     return escript.Solution(domain)
164 jgs 151 elif self.what==self.REDUCED:
165     return escript.ReducedSolution(domain)
166 jgs 102
167 jgs 108 def getValue(self):
168     """
169 jgs 149 returns the value of the coefficient
170    
171     @return: value of the coefficient
172     @rtype: L{Data<escript.Data>}
173 jgs 108 """
174     return self.value
175 jgs 148
176 jgs 149 def setValue(self,domain,numEquations=1,numSolutions=1,reducedEquationOrder=False,reducedSolutionOrder=False,newValue=None):
177 jgs 108 """
178 jgs 149 set the value of the coefficient to a new value
179    
180     @param domain: domain on which the PDE uses the coefficient
181     @type domain: L{Domain<escript.Domain>}
182     @param numEquations: number of equations of the PDE
183     @type numEquations: C{int}
184     @param numSolutions: number of components of the PDE solution
185     @type numSolutions: C{int}
186     @param reducedEquationOrder: True to indicate that reduced order is used to represent the equation
187 gross 720 @type reducedEquationOrder: C{bool}
188 jgs 149 @param reducedSolutionOrder: True to indicate that reduced order is used to represent the solution
189 gross 720 @type reducedSolutionOrder: C{bool}
190 jgs 149 @param newValue: number of components of the PDE solution
191     @type newValue: any object that can be converted into a L{Data<escript.Data>} object with the appropriate shape and L{FunctionSpace<escript.FunctionSpace>}
192     @raise IllegalCoefficientValue: if the shape of the assigned value does not match the shape of the coefficient
193 gross 1072 @raise IllegalCoefficientFunctionSpace: if unable to interploate value to appropriate function space
194 jgs 108 """
195 jgs 148 if newValue==None:
196     newValue=escript.Data()
197     elif isinstance(newValue,escript.Data):
198     if not newValue.isEmpty():
199 gross 1072 if not newValue.getFunctionSpace() == self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder):
200     try:
201     newValue=escript.Data(newValue,self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder))
202     except:
203     raise IllegalCoefficientFunctionSpace,"Unable to interpolate coefficient to function space %s"%self.getFunctionSpace(domain)
204 jgs 148 else:
205 jgs 149 newValue=escript.Data(newValue,self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder))
206 jgs 148 if not newValue.isEmpty():
207     if not self.getShape(domain,numEquations,numSolutions)==newValue.getShape():
208 jgs 149 raise IllegalCoefficientValue,"Expected shape of coefficient is %s but actual shape is %s."%(self.getShape(domain,numEquations,numSolutions),newValue.getShape())
209 jgs 108 self.value=newValue
210 jgs 148
211 jgs 102 def isAlteringOperator(self):
212     """
213 jgs 149 checks if the coefficient alters the operator of the PDE
214    
215     @return: True if the operator of the PDE is changed when the coefficient is changed
216     @rtype: C{bool}
217 jgs 102 """
218     if self.altering==self.OPERATOR or self.altering==self.BOTH:
219     return not None
220     else:
221     return None
222    
223     def isAlteringRightHandSide(self):
224     """
225 jgs 149 checks if the coefficeint alters the right hand side of the PDE
226    
227     @rtype: C{bool}
228     @return: True if the right hand side of the PDE is changed when the coefficient is changed
229 jgs 102 """
230     if self.altering==self.RIGHTHANDSIDE or self.altering==self.BOTH:
231     return not None
232     else:
233     return None
234    
235 jgs 148 def estimateNumEquationsAndNumSolutions(self,domain,shape=()):
236 jgs 102 """
237 jgs 149 tries to estimate the number of equations and number of solutions if the coefficient has the given shape
238 jgs 102
239 jgs 149 @param domain: domain on which the PDE uses the coefficient
240     @type domain: L{Domain<escript.Domain>}
241     @param shape: suggested shape of the coefficient
242     @type shape: C{tuple} of C{int} values
243     @return: the number of equations and number of solutions of the PDE is the coefficient has shape s.
244     If no appropriate numbers could be identified, C{None} is returned
245     @rtype: C{tuple} of two C{int} values or C{None}
246 jgs 102 """
247 jgs 148 dim=domain.getDim()
248 jgs 102 if len(shape)>0:
249     num=max(shape)+1
250     else:
251     num=1
252     search=[]
253 jgs 149 if self.definesNumEquation() and self.definesNumSolutions():
254     for u in range(num):
255     for e in range(num):
256     search.append((e,u))
257     search.sort(self.__CompTuple2)
258     for item in search:
259 jgs 148 s=self.getShape(domain,item[0],item[1])
260 jgs 102 if len(s)==0 and len(shape)==0:
261     return (1,1)
262     else:
263     if s==shape: return item
264 jgs 149 elif self.definesNumEquation():
265     for e in range(num,0,-1):
266     s=self.getShape(domain,e,0)
267     if len(s)==0 and len(shape)==0:
268     return (1,None)
269     else:
270     if s==shape: return (e,None)
271    
272     elif self.definesNumSolutions():
273     for u in range(num,0,-1):
274     s=self.getShape(domain,0,u)
275     if len(s)==0 and len(shape)==0:
276     return (None,1)
277     else:
278     if s==shape: return (None,u)
279 jgs 102 return None
280 jgs 149 def definesNumSolutions(self):
281     """
282     checks if the coefficient allows to estimate the number of solution components
283 jgs 102
284 jgs 149 @return: True if the coefficient allows an estimate of the number of solution components
285     @rtype: C{bool}
286     """
287     for i in self.pattern:
288     if i==self.BY_SOLUTION: return True
289     return False
290    
291     def definesNumEquation(self):
292     """
293     checks if the coefficient allows to estimate the number of equations
294    
295     @return: True if the coefficient allows an estimate of the number of equations
296     @rtype: C{bool}
297     """
298     for i in self.pattern:
299     if i==self.BY_EQUATION: return True
300     return False
301    
302     def __CompTuple2(self,t1,t2):
303     """
304     Compare two tuples of possible number of equations and number of solutions
305    
306     @param t1: The first tuple
307     @param t2: The second tuple
308    
309     """
310    
311     dif=t1[0]+t1[1]-(t2[0]+t2[1])
312     if dif<0: return 1
313     elif dif>0: return -1
314     else: return 0
315    
316 jgs 148 def getShape(self,domain,numEquations=1,numSolutions=1):
317 jgs 149 """
318     builds the required shape of the coefficient
319 jgs 102
320 jgs 149 @param domain: domain on which the PDE uses the coefficient
321     @type domain: L{Domain<escript.Domain>}
322     @param numEquations: number of equations of the PDE
323     @type numEquations: C{int}
324     @param numSolutions: number of components of the PDE solution
325     @type numSolutions: C{int}
326     @return: shape of the coefficient
327     @rtype: C{tuple} of C{int} values
328     """
329     dim=domain.getDim()
330     s=()
331     for i in self.pattern:
332     if i==self.BY_EQUATION:
333 jgs 148 if numEquations>1: s=s+(numEquations,)
334 jgs 149 elif i==self.BY_SOLUTION:
335 jgs 148 if numSolutions>1: s=s+(numSolutions,)
336 jgs 102 else:
337     s=s+(dim,)
338 jgs 149 return s
339 jgs 102
340 gross 1841 class LinearProblem(object):
341 jgs 102 """
342 gross 1841 This is the base class is to define a general linear PDE-type problem for an u
343     for an unknown function M{u} on a given domain defined through a L{Domain<escript.Domain>} object.
344     The problem can be given as a single equations or as a system of equations.
345 jgs 102
346 gross 1841 The class assumes that some sort of assembling process is required to form a problem of the form
347 jgs 151
348 gross 1841 M{L u=f}
349    
350     where M{L} is an operator and M{f} is the right hand side. This operator problem will be solved to get the
351     unknown M{u}.
352 jgs 102
353 jgs 150 @cvar DEFAULT: The default method used to solve the system of linear equations
354 jgs 149 @cvar DIRECT: The direct solver based on LDU factorization
355     @cvar CHOLEVSKY: The direct solver based on LDLt factorization (can only be applied for symmetric PDEs)
356     @cvar PCG: The preconditioned conjugate gradient method (can only be applied for symmetric PDEs)
357 jgs 151 @cvar CR: The conjugate residual method
358 jgs 149 @cvar CGS: The conjugate gardient square method
359     @cvar BICGSTAB: The stabilized BiConjugate Gradient method.
360 artak 1703 @cvar TFQMR: Transport Free Quasi Minimal Residual method.
361 artak 1787 @cvar MINRES: Minimum residual method.
362 jgs 149 @cvar SSOR: The symmetric overrealaxtion method
363     @cvar ILU0: The incomplete LU factorization preconditioner with no fill in
364     @cvar ILUT: The incomplete LU factorization preconditioner with will in
365     @cvar JACOBI: The Jacobi preconditioner
366     @cvar GMRES: The Gram-Schmidt minimum residual method
367     @cvar PRES20: Special GMRES with restart after 20 steps and truncation after 5 residuals
368     @cvar LUMPING: Matrix lumping.
369     @cvar NO_REORDERING: No matrix reordering allowed
370     @cvar MINIMUM_FILL_IN: Reorder matrix to reduce fill-in during factorization
371     @cvar NESTED_DISSECTION: Reorder matrix to improve load balancing during factorization
372 jgs 150 @cvar PASO: PASO solver package
373     @cvar SCSL: SGI SCSL solver library
374     @cvar MKL: Intel's MKL solver library
375 jgs 151 @cvar UMFPACK: the UMFPACK library
376 ksteube 1312 @cvar TRILINOS: the TRILINOS parallel solver class library from Sandia Natl Labs
377 jgs 150 @cvar ITERATIVE: The default iterative solver
378 gross 430 @cvar AMG: algebraic multi grid
379     @cvar RILU: recursive ILU
380 artak 1819 @cvar GS: Gauss-Seidel solver
381 jgs 149
382 jgs 102 """
383 jgs 150 DEFAULT= 0
384     DIRECT= 1
385     CHOLEVSKY= 2
386     PCG= 3
387     CR= 4
388     CGS= 5
389     BICGSTAB= 6
390     SSOR= 7
391     ILU0= 8
392     ILUT= 9
393     JACOBI= 10
394     GMRES= 11
395     PRES20= 12
396     LUMPING= 13
397     NO_REORDERING= 17
398     MINIMUM_FILL_IN= 18
399     NESTED_DISSECTION= 19
400     SCSL= 14
401     MKL= 15
402     UMFPACK= 16
403     ITERATIVE= 20
404     PASO= 21
405 gross 430 AMG= 22
406     RILU = 23
407 ksteube 1312 TRILINOS = 24
408 gross 1639 NONLINEAR_GMRES = 25
409 artak 1703 TFQMR = 26
410 artak 1787 MINRES = 27
411 artak 1819 GS=28
412 jgs 150
413 gross 596 SMALL_TOLERANCE=1.e-13
414 gross 1841 PACKAGE_KEY="package"
415     METHOD_KEY="method"
416     SYMMETRY_KEY="symmetric"
417     TOLERANCE_KEY="tolerance"
418     PRECONDITIONER_KEY="preconditioner"
419 jgs 102
420 jgs 148
421     def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):
422 jgs 102 """
423 gross 1841 initializes a linear problem
424 jgs 102
425 jgs 148 @param domain: domain of the PDE
426 jgs 149 @type domain: L{Domain<escript.Domain>}
427 jgs 148 @param numEquations: number of equations. If numEquations==None the number of equations
428 gross 1841 is extracted from the coefficients.
429 jgs 148 @param numSolutions: number of solution components. If numSolutions==None the number of solution components
430 gross 1841 is extracted from the coefficients.
431 jgs 148 @param debug: if True debug informations are printed.
432    
433 jgs 102 """
434 gross 1841 super(LinearProblem, self).__init__()
435 jgs 102
436 jgs 148 self.__debug=debug
437 jgs 104 self.__domain=domain
438     self.__numEquations=numEquations
439     self.__numSolutions=numSolutions
440 gross 1841 self.__altered_coefficients=False
441 jgs 149 self.__reduce_equation_order=False
442     self.__reduce_solution_order=False
443 jgs 102 self.__tolerance=1.e-8
444 jgs 150 self.__solver_method=self.DEFAULT
445     self.__solver_package=self.DEFAULT
446 gross 387 self.__preconditioner=self.DEFAULT
447 gross 1841 self.__is_RHS_valid=False
448     self.__is_operator_valid=False
449 jgs 102 self.__sym=False
450 gross 1841 self.__COEFFICIENTS={}
451     # initialize things:
452     self.resetAllCoefficients()
453     self.__system_type=None
454     self.updateSystemType()
455 jgs 148 # =============================================================================
456     # general stuff:
457     # =============================================================================
458     def __str__(self):
459 jgs 149 """
460     returns string representation of the PDE
461    
462     @return: a simple representation of the PDE
463     @rtype: C{str}
464     """
465 gross 1841 return "<LinearProblem %d>"%id(self)
466 jgs 148 # =============================================================================
467     # debug :
468     # =============================================================================
469     def setDebugOn(self):
470 jgs 149 """
471 jgs 148 switches on debugging
472 jgs 108 """
473 jgs 148 self.__debug=not None
474    
475     def setDebugOff(self):
476 jgs 108 """
477 jgs 148 switches off debugging
478     """
479     self.__debug=None
480 jgs 108
481 gross 2100 def setDebug(self, flag):
482     """
483     switches debugging to on if C{flag} is true outherwise debugging is set to off
484    
485     @param flag: desired debug status
486     @type flag: C{bool}
487     """
488     if flag:
489     self.setDebugOn()
490     else:
491     self.setDebugOff()
492    
493 jgs 148 def trace(self,text):
494     """
495 jgs 149 print the text message if debugging is swiched on.
496     @param text: message
497     @type text: C{string}
498 jgs 102 """
499 jgs 148 if self.__debug: print "%s: %s"%(str(self),text)
500 jgs 102
501 jgs 148 # =============================================================================
502     # some service functions:
503     # =============================================================================
504 gross 1841 def introduceCoefficients(self,**coeff):
505     """
506     introduces a new coefficient into the problem.
507    
508     use
509    
510 jfenwick 2061 self.introduceCoefficients(A=PDECoef(...),B=PDECoef(...))
511 gross 1841
512     to introduce the coefficients M{A} ans M{B}.
513     """
514     for name, type in coeff.items():
515     if not isinstance(type,PDECoef):
516     raise ValueError,"coefficient %s has no type."%name
517     self.__COEFFICIENTS[name]=type
518     self.__COEFFICIENTS[name].resetValue()
519     self.trace("coefficient %s has been introduced."%name)
520    
521 jgs 148 def getDomain(self):
522 jgs 102 """
523 jgs 148 returns the domain of the PDE
524 jgs 102
525 jgs 149 @return: the domain of the PDE
526     @rtype: L{Domain<escript.Domain>}
527 jgs 108 """
528 jgs 148 return self.__domain
529 jgs 122
530 jgs 148 def getDim(self):
531 jgs 108 """
532 jgs 148 returns the spatial dimension of the PDE
533 jgs 108
534 jgs 149 @return: the spatial dimension of the PDE domain
535     @rtype: C{int}
536 jgs 148 """
537     return self.getDomain().getDim()
538 jgs 102
539 jgs 148 def getNumEquations(self):
540     """
541     returns the number of equations
542 jgs 102
543 jgs 149 @return: the number of equations
544     @rtype: C{int}
545 jgs 148 @raise UndefinedPDEError: if the number of equations is not be specified yet.
546     """
547     if self.__numEquations==None:
548 gross 1859 if self.__numSolutions==None:
549     raise UndefinedPDEError,"Number of equations is undefined. Please specify argument numEquations."
550     else:
551     self.__numEquations=self.__numSolutions
552     return self.__numEquations
553 jgs 147
554 jgs 148 def getNumSolutions(self):
555     """
556     returns the number of unknowns
557 jgs 147
558 jgs 149 @return: the number of unknowns
559     @rtype: C{int}
560 jgs 148 @raise UndefinedPDEError: if the number of unknowns is not be specified yet.
561     """
562     if self.__numSolutions==None:
563 gross 1859 if self.__numEquations==None:
564     raise UndefinedPDEError,"Number of solution is undefined. Please specify argument numSolutions."
565     else:
566     self.__numSolutions=self.__numEquations
567     return self.__numSolutions
568 jgs 148
569 jgs 149 def reduceEquationOrder(self):
570     """
571     return status for order reduction for equation
572    
573     @return: return True is reduced interpolation order is used for the represenation of the equation
574     @rtype: L{bool}
575     """
576     return self.__reduce_equation_order
577    
578     def reduceSolutionOrder(self):
579     """
580     return status for order reduction for the solution
581    
582     @return: return True is reduced interpolation order is used for the represenation of the solution
583     @rtype: L{bool}
584     """
585     return self.__reduce_solution_order
586 jgs 151
587 jgs 108 def getFunctionSpaceForEquation(self):
588     """
589 jgs 149 returns the L{FunctionSpace<escript.FunctionSpace>} used to discretize the equation
590 jgs 148
591 jgs 149 @return: representation space of equation
592     @rtype: L{FunctionSpace<escript.FunctionSpace>}
593 jgs 108 """
594 jgs 149 if self.reduceEquationOrder():
595     return escript.ReducedSolution(self.getDomain())
596     else:
597     return escript.Solution(self.getDomain())
598 jgs 108
599     def getFunctionSpaceForSolution(self):
600     """
601 jgs 149 returns the L{FunctionSpace<escript.FunctionSpace>} used to represent the solution
602 jgs 148
603 jgs 149 @return: representation space of solution
604     @rtype: L{FunctionSpace<escript.FunctionSpace>}
605 jgs 108 """
606 jgs 149 if self.reduceSolutionOrder():
607     return escript.ReducedSolution(self.getDomain())
608     else:
609     return escript.Solution(self.getDomain())
610 jgs 108
611 jgs 148 # =============================================================================
612     # solver settings:
613     # =============================================================================
614 gross 387 def setSolverMethod(self,solver=None,preconditioner=None):
615 jgs 102 """
616 jgs 122 sets a new solver
617 jgs 148
618 jgs 149 @param solver: sets a new solver method.
619 artak 1787 @type solver: one of L{DEFAULT}, L{ITERATIVE} L{DIRECT}, L{CHOLEVSKY}, L{PCG}, L{CR}, L{CGS}, L{BICGSTAB}, L{SSOR}, L{GMRES}, L{TFQMR}, L{MINRES}, L{PRES20}, L{LUMPING}, L{AMG}
620 gross 387 @param preconditioner: sets a new solver method.
621 artak 1819 @type preconditioner: one of L{DEFAULT}, L{JACOBI} L{ILU0}, L{ILUT},L{SSOR}, L{RILU}, L{GS}
622 gross 1841
623     @note: the solver method choosen may not be available by the selected package.
624 jgs 102 """
625 ksteube 1312 if solver==None: solver=self.__solver_method
626     if preconditioner==None: preconditioner=self.__preconditioner
627     if solver==None: solver=self.DEFAULT
628 gross 387 if preconditioner==None: preconditioner=self.DEFAULT
629     if not (solver,preconditioner)==self.getSolverMethod():
630 jgs 102 self.__solver_method=solver
631 gross 387 self.__preconditioner=preconditioner
632 gross 1841 self.updateSystemType()
633 jgs 148 self.trace("New solver is %s"%self.getSolverMethodName())
634 jgs 102
635 gross 1841 def getSolverMethodName(self):
636 jgs 148 """
637     returns the name of the solver currently used
638    
639 jgs 149 @return: the name of the solver currently used.
640 jgs 148 @rtype: C{string}
641     """
642    
643     m=self.getSolverMethod()
644 jgs 150 p=self.getSolverPackage()
645 gross 425 method=""
646 gross 387 if m[0]==self.DEFAULT: method="DEFAULT"
647     elif m[0]==self.DIRECT: method= "DIRECT"
648     elif m[0]==self.ITERATIVE: method= "ITERATIVE"
649     elif m[0]==self.CHOLEVSKY: method= "CHOLEVSKY"
650     elif m[0]==self.PCG: method= "PCG"
651 artak 1787 elif m[0]==self.TFQMR: method= "TFQMR"
652     elif m[0]==self.MINRES: method= "MINRES"
653 gross 387 elif m[0]==self.CR: method= "CR"
654     elif m[0]==self.CGS: method= "CGS"
655     elif m[0]==self.BICGSTAB: method= "BICGSTAB"
656     elif m[0]==self.SSOR: method= "SSOR"
657     elif m[0]==self.GMRES: method= "GMRES"
658     elif m[0]==self.PRES20: method= "PRES20"
659     elif m[0]==self.LUMPING: method= "LUMPING"
660 gross 431 elif m[0]==self.AMG: method= "AMG"
661 gross 425 if m[1]==self.DEFAULT: method+="+DEFAULT"
662     elif m[1]==self.JACOBI: method+= "+JACOBI"
663     elif m[1]==self.ILU0: method+= "+ILU0"
664     elif m[1]==self.ILUT: method+= "+ILUT"
665     elif m[1]==self.SSOR: method+= "+SSOR"
666 gross 431 elif m[1]==self.AMG: method+= "+AMG"
667     elif m[1]==self.RILU: method+= "+RILU"
668 artak 1819 elif m[1]==self.GS: method+= "+GS"
669 jgs 150 if p==self.DEFAULT: package="DEFAULT"
670     elif p==self.PASO: package= "PASO"
671     elif p==self.MKL: package= "MKL"
672     elif p==self.SCSL: package= "SCSL"
673     elif p==self.UMFPACK: package= "UMFPACK"
674 ksteube 1312 elif p==self.TRILINOS: package= "TRILINOS"
675 jgs 150 else : method="unknown"
676     return "%s solver of %s package"%(method,package)
677 jgs 148
678 gross 1841 def getSolverMethod(self):
679 jgs 102 """
680 gross 1841 returns the solver method and preconditioner used
681 jgs 149
682 jgs 151 @return: the solver method currently be used.
683 gross 1841 @rtype: C{tuple} of C{int}
684 jgs 102 """
685 gross 387 return self.__solver_method,self.__preconditioner
686 jgs 102
687 gross 1841 def setSolverPackage(self,package=None):
688 jgs 150 """
689     sets a new solver package
690    
691 gross 720 @param package: sets a new solver method.
692 ksteube 1312 @type package: one of L{DEFAULT}, L{PASO} L{SCSL}, L{MKL}, L{UMFPACK}, L{TRILINOS}
693 jgs 150 """
694     if package==None: package=self.DEFAULT
695     if not package==self.getSolverPackage():
696 gross 727 self.__solver_package=package
697 gross 1841 self.updateSystemType()
698 jgs 150 self.trace("New solver is %s"%self.getSolverMethodName())
699    
700 gross 1841 def getSolverPackage(self):
701 jgs 150 """
702     returns the package of the solver
703    
704 jgs 151 @return: the solver package currently being used.
705 jgs 150 @rtype: C{int}
706     """
707     return self.__solver_package
708    
709 jgs 148 def isUsingLumping(self):
710     """
711 gross 1841 checks if matrix lumping is used as a solver method
712 jgs 148
713 jgs 149 @return: True is lumping is currently used a solver method.
714 jgs 148 @rtype: C{bool}
715     """
716 gross 387 return self.getSolverMethod()[0]==self.LUMPING
717 jgs 148
718 gross 1841 def setTolerance(self,tol=1.e-8):
719 jgs 102 """
720 gross 1841 resets the tolerance for the solver method to tol.
721 jgs 148
722     @param tol: new tolerance for the solver. If the tol is lower then the current tolerence
723     the system will be resolved.
724 jgs 149 @type tol: positive C{float}
725 gross 877 @raise ValueError: if tolerance is not positive.
726 jgs 102 """
727     if not tol>0:
728 gross 877 raise ValueError,"Tolerance as to be positive"
729 gross 1841 if tol<self.getTolerance(): self.invalidateSolution()
730 jgs 148 self.trace("New tolerance %e"%tol)
731 jgs 102 self.__tolerance=tol
732     return
733 jgs 148
734 gross 1841 def getTolerance(self):
735 jgs 102 """
736 jgs 122 returns the tolerance set for the solution
737 jgs 148
738     @return: tolerance currently used.
739     @rtype: C{float}
740 jgs 102 """
741     return self.__tolerance
742    
743 jgs 148 # =============================================================================
744     # symmetry flag:
745     # =============================================================================
746 gross 1841 def isSymmetric(self):
747 jgs 102 """
748 jgs 148 checks if symmetry is indicated.
749 jgs 149
750     @return: True is a symmetric PDE is indicated, otherwise False is returned
751     @rtype: C{bool}
752 jgs 102 """
753     return self.__sym
754    
755 gross 1841 def setSymmetryOn(self):
756 jgs 102 """
757 jgs 148 sets the symmetry flag.
758 jgs 102 """
759     if not self.isSymmetric():
760 jgs 148 self.trace("PDE is set to be symmetric")
761 jgs 102 self.__sym=True
762 gross 1841 self.updateSystemType()
763 jgs 102
764 gross 1841 def setSymmetryOff(self):
765 jgs 102 """
766 jgs 148 removes the symmetry flag.
767 jgs 102 """
768     if self.isSymmetric():
769 jgs 148 self.trace("PDE is set to be unsymmetric")
770 jgs 102 self.__sym=False
771 gross 1841 self.updateSystemType()
772 jgs 102
773 gross 1841 def setSymmetryTo(self,flag=False):
774 jgs 148 """
775     sets the symmetry flag to flag
776 jgs 149
777 jgs 148 @param flag: If flag, the symmetry flag is set otherwise the symmetry flag is released.
778     @type flag: C{bool}
779     """
780     if flag:
781     self.setSymmetryOn()
782     else:
783     self.setSymmetryOff()
784 jgs 102
785 jgs 148 # =============================================================================
786     # function space handling for the equation as well as the solution
787     # =============================================================================
788 gross 1841 def setReducedOrderOn(self):
789 jgs 102 """
790 jgs 148 switches on reduced order for solution and equation representation
791 jgs 149
792     @raise RuntimeError: if order reduction is altered after a coefficient has been set.
793 jgs 102 """
794     self.setReducedOrderForSolutionOn()
795     self.setReducedOrderForEquationOn()
796    
797 gross 1841 def setReducedOrderOff(self):
798 jgs 102 """
799 jgs 148 switches off reduced order for solution and equation representation
800 jgs 149
801     @raise RuntimeError: if order reduction is altered after a coefficient has been set.
802 jgs 102 """
803     self.setReducedOrderForSolutionOff()
804     self.setReducedOrderForEquationOff()
805    
806 gross 1841 def setReducedOrderTo(self,flag=False):
807 jgs 102 """
808 jgs 149 sets order reduction for both solution and equation representation according to flag.
809     @param flag: if flag is True, the order reduction is switched on for both solution and equation representation, otherwise or
810 jgs 148 if flag is not present order reduction is switched off
811     @type flag: C{bool}
812 jgs 149 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
813 jgs 102 """
814     self.setReducedOrderForSolutionTo(flag)
815     self.setReducedOrderForEquationTo(flag)
816    
817 jgs 148
818 gross 1841 def setReducedOrderForSolutionOn(self):
819 jgs 102 """
820 jgs 148 switches on reduced order for solution representation
821 jgs 149
822     @raise RuntimeError: if order reduction is altered after a coefficient has been set.
823 jgs 102 """
824 jgs 149 if not self.__reduce_solution_order:
825     if self.__altered_coefficients:
826     raise RuntimeError,"order cannot be altered after coefficients have been defined."
827 jgs 148 self.trace("Reduced order is used to solution representation.")
828 jgs 149 self.__reduce_solution_order=True
829 gross 1841 self.initializeSystem()
830 jgs 102
831 gross 1841 def setReducedOrderForSolutionOff(self):
832 jgs 102 """
833 jgs 148 switches off reduced order for solution representation
834 jgs 149
835     @raise RuntimeError: if order reduction is altered after a coefficient has been set.
836 jgs 102 """
837 jgs 149 if self.__reduce_solution_order:
838     if self.__altered_coefficients:
839     raise RuntimeError,"order cannot be altered after coefficients have been defined."
840 jgs 148 self.trace("Full order is used to interpolate solution.")
841 jgs 149 self.__reduce_solution_order=False
842 gross 1841 self.initializeSystem()
843 jgs 102
844 gross 1841 def setReducedOrderForSolutionTo(self,flag=False):
845 jgs 102 """
846 jgs 122 sets order for test functions according to flag
847 jgs 102
848 jgs 149 @param flag: if flag is True, the order reduction is switched on for solution representation, otherwise or
849 jgs 148 if flag is not present order reduction is switched off
850     @type flag: C{bool}
851 jgs 149 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
852 jgs 102 """
853     if flag:
854     self.setReducedOrderForSolutionOn()
855     else:
856     self.setReducedOrderForSolutionOff()
857 jgs 148
858 gross 1841 def setReducedOrderForEquationOn(self):
859 jgs 102 """
860 jgs 148 switches on reduced order for equation representation
861 jgs 149
862     @raise RuntimeError: if order reduction is altered after a coefficient has been set.
863 jgs 102 """
864 jgs 149 if not self.__reduce_equation_order:
865     if self.__altered_coefficients:
866     raise RuntimeError,"order cannot be altered after coefficients have been defined."
867 jgs 148 self.trace("Reduced order is used for test functions.")
868 jgs 149 self.__reduce_equation_order=True
869 gross 1841 self.initializeSystem()
870 jgs 102
871 gross 1841 def setReducedOrderForEquationOff(self):
872 jgs 102 """
873 jgs 148 switches off reduced order for equation representation
874 jgs 149
875     @raise RuntimeError: if order reduction is altered after a coefficient has been set.
876 jgs 102 """
877 jgs 149 if self.__reduce_equation_order:
878     if self.__altered_coefficients:
879     raise RuntimeError,"order cannot be altered after coefficients have been defined."
880 jgs 148 self.trace("Full order is used for test functions.")
881 jgs 149 self.__reduce_equation_order=False
882 gross 1841 self.initializeSystem()
883 jgs 102
884 gross 1841 def setReducedOrderForEquationTo(self,flag=False):
885 jgs 102 """
886 jgs 122 sets order for test functions according to flag
887 jgs 102
888 jgs 149 @param flag: if flag is True, the order reduction is switched on for equation representation, otherwise or
889 jgs 148 if flag is not present order reduction is switched off
890     @type flag: C{bool}
891 jgs 149 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
892 jgs 102 """
893     if flag:
894     self.setReducedOrderForEquationOn()
895     else:
896     self.setReducedOrderForEquationOff()
897 jgs 148
898 gross 1841 def updateSystemType(self):
899 jgs 148 """
900 gross 1841 reasses the system type and, if a new matrix is needed, resets the system.
901 jgs 148 """
902 gross 1841 new_system_type=self.getRequiredSystemType()
903     if not new_system_type==self.__system_type:
904     self.trace("Matrix type is now %d."%new_system_type)
905     self.__system_type=new_system_type
906     self.initializeSystem()
907     def getSystemType(self):
908     """
909     return the current system type
910     """
911     return self.__system_type
912 jgs 148
913 gross 1841 def checkSymmetricTensor(self,name,verbose=True):
914     """
915     tests a coefficient for symmetry
916 jgs 148
917 gross 1841 @param name: name of the coefficient
918     @type name: C{str}
919     @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.
920     @type verbose: C{bool}
921     @return: True if coefficient C{name} is symmetric.
922     @rtype: C{bool}
923     """
924     A=self.getCoefficient(name)
925     verbose=verbose or self.__debug
926     out=True
927     if not A.isEmpty():
928     tol=util.Lsup(A)*self.SMALL_TOLERANCE
929     s=A.getShape()
930     if A.getRank() == 4:
931     if s[0]==s[2] and s[1] == s[3]:
932     for i in range(s[0]):
933     for j in range(s[1]):
934     for k in range(s[2]):
935     for l in range(s[3]):
936     if util.Lsup(A[i,j,k,l]-A[k,l,i,j])>tol:
937     if verbose: print "non-symmetric problem as %s[%d,%d,%d,%d]!=%s[%d,%d,%d,%d]"%(name,i,j,k,l,name,k,l,i,j)
938     out=False
939     else:
940     if verbose: print "non-symmetric problem because of inappropriate shape %s of coefficient %s."%(s,name)
941     out=False
942     elif A.getRank() == 2:
943     if s[0]==s[1]:
944     for j in range(s[0]):
945     for l in range(s[1]):
946     if util.Lsup(A[j,l]-A[l,j])>tol:
947     if verbose: print "non-symmetric problem because %s[%d,%d]!=%s[%d,%d]"%(name,j,l,name,l,j)
948     out=False
949     else:
950     if verbose: print "non-symmetric problem because of inappropriate shape %s of coefficient %s."%(s,name)
951     out=False
952     elif A.getRank() == 0:
953     pass
954     else:
955     raise ValueError,"Cannot check rank %s of %s."%(A.getRank(),name)
956     return out
957 jgs 148
958 gross 1841 def checkReciprocalSymmetry(self,name0,name1,verbose=True):
959     """
960     test two coefficients for resciprocal symmetry
961 jgs 148
962 gross 1841 @param name0: name of the first coefficient
963 jfenwick 2061 @type name0: C{str}
964 gross 1841 @param name1: name of the second coefficient
965 jfenwick 2061 @type name1: C{str}
966 gross 1841 @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.
967     @type verbose: C{bool}
968     @return: True if coefficients C{name0} and C{name1} are resciprocally symmetric.
969     @rtype: C{bool}
970     """
971     B=self.getCoefficient(name0)
972     C=self.getCoefficient(name1)
973     verbose=verbose or self.__debug
974     out=True
975     if B.isEmpty() and not C.isEmpty():
976     if verbose: print "non-symmetric problem because %s is not present but %s is"%(name0,name1)
977     out=False
978     elif not B.isEmpty() and C.isEmpty():
979     if verbose: print "non-symmetric problem because %s is not present but %s is"%(name0,name1)
980     out=False
981     elif not B.isEmpty() and not C.isEmpty():
982     sB=B.getShape()
983     sC=C.getShape()
984     tol=(util.Lsup(B)+util.Lsup(C))*self.SMALL_TOLERANCE/2.
985     if len(sB) != len(sC):
986     if verbose: print "non-symmetric problem because ranks of %s (=%s) and %s (=%s) are different."%(name0,len(sB),name1,len(sC))
987     out=False
988     else:
989     if len(sB)==0:
990     if util.Lsup(B-C)>tol:
991     if verbose: print "non-symmetric problem because %s!=%s"%(name0,name1)
992     out=False
993     elif len(sB)==1:
994     if sB[0]==sC[0]:
995     for j in range(sB[0]):
996     if util.Lsup(B[j]-C[j])>tol:
997     if verbose: print "non-symmetric PDE because %s[%d]!=%s[%d]"%(name0,j,name1,j)
998     out=False
999     else:
1000     if verbose: print "non-symmetric problem because of inappropriate shapes %s and %s of coefficients %s and %s, respectively."%(sB,sC,name0,name1)
1001     elif len(sB)==3:
1002     if sB[0]==sC[1] and sB[1]==sC[2] and sB[2]==sC[0]:
1003     for i in range(sB[0]):
1004     for j in range(sB[1]):
1005     for k in range(sB[2]):
1006     if util.Lsup(B[i,j,k]-C[k,i,j])>tol:
1007 gross 2100 if verbose: print "non-symmetric problem because %s[%d,%d,%d]!=%s[%d,%d,%d]"%(name0,i,j,k,name1,k,i,j)
1008 gross 1841 out=False
1009     else:
1010     if verbose: print "non-symmetric problem because of inappropriate shapes %s and %s of coefficients %s and %s, respectively."%(sB,sC,name0,name1)
1011 jgs 148 else:
1012 gross 1841 raise ValueError,"Cannot check rank %s of %s and %s."%(len(sB),name0,name1)
1013     return out
1014 jgs 148
1015 gross 1841 def getCoefficient(self,name):
1016 jgs 108 """
1017 jgs 148 returns the value of the coefficient name
1018 jgs 108
1019 jgs 149 @param name: name of the coefficient requested.
1020 jgs 148 @type name: C{string}
1021 jgs 149 @return: the value of the coefficient name
1022     @rtype: L{Data<escript.Data>}
1023 jgs 148 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1024     """
1025     if self.hasCoefficient(name):
1026 gross 1841 return self.__COEFFICIENTS[name].getValue()
1027 jgs 148 else:
1028     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1029 jgs 108
1030 gross 1841 def hasCoefficient(self,name):
1031 jgs 148 """
1032     return True if name is the name of a coefficient
1033 jgs 108
1034 jgs 148 @param name: name of the coefficient enquired.
1035     @type name: C{string}
1036 jgs 149 @return: True if name is the name of a coefficient of the general PDE. Otherwise False.
1037     @rtype: C{bool}
1038 jgs 148 """
1039 gross 1841 return self.__COEFFICIENTS.has_key(name)
1040 jgs 108
1041 gross 1841 def createCoefficient(self, name):
1042 jgs 148 """
1043 jgs 149 create a L{Data<escript.Data>} object corresponding to coefficient name
1044 jgs 108
1045 jgs 149 @return: a coefficient name initialized to 0.
1046     @rtype: L{Data<escript.Data>}
1047 jgs 148 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1048     """
1049     if self.hasCoefficient(name):
1050     return escript.Data(0.,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))
1051     else:
1052     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1053    
1054 gross 1841 def getFunctionSpaceForCoefficient(self,name):
1055 jgs 148 """
1056 jgs 149 return the L{FunctionSpace<escript.FunctionSpace>} to be used for coefficient name
1057 jgs 148
1058     @param name: name of the coefficient enquired.
1059     @type name: C{string}
1060 jgs 149 @return: the function space to be used for coefficient name
1061     @rtype: L{FunctionSpace<escript.FunctionSpace>}
1062 jgs 148 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1063     """
1064     if self.hasCoefficient(name):
1065 gross 1841 return self.__COEFFICIENTS[name].getFunctionSpace(self.getDomain())
1066 jgs 148 else:
1067     raise ValueError,"unknown coefficient %s requested"%name
1068 gross 1841
1069     def getShapeOfCoefficient(self,name):
1070 jgs 148 """
1071 jgs 149 return the shape of the coefficient name
1072 jgs 148
1073     @param name: name of the coefficient enquired.
1074     @type name: C{string}
1075 jgs 149 @return: the shape of the coefficient name
1076     @rtype: C{tuple} of C{int}
1077 jgs 148 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1078     """
1079     if self.hasCoefficient(name):
1080 gross 1841 return self.__COEFFICIENTS[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())
1081 jgs 148 else:
1082     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1083    
1084 gross 1841 def resetAllCoefficients(self):
1085 jgs 148 """
1086     resets all coefficients to there default values.
1087     """
1088 gross 1841 for i in self.__COEFFICIENTS.iterkeys():
1089     self.__COEFFICIENTS[i].resetValue()
1090 jgs 148
1091 gross 1841 def alteredCoefficient(self,name):
1092 jgs 148 """
1093     announce that coefficient name has been changed
1094    
1095     @param name: name of the coefficient enquired.
1096     @type name: C{string}
1097     @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1098 jgs 149 @note: if name is q or r, the method will not trigger a rebuilt of the system as constraints are applied to the solved system.
1099 jgs 148 """
1100     if self.hasCoefficient(name):
1101     self.trace("Coefficient %s has been altered."%name)
1102 jgs 149 if not ((name=="q" or name=="r") and self.isUsingLumping()):
1103 gross 1841 if self.__COEFFICIENTS[name].isAlteringOperator(): self.invalidateOperator()
1104     if self.__COEFFICIENTS[name].isAlteringRightHandSide(): self.invalidateRightHandSide()
1105 jgs 148 else:
1106     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1107    
1108 gross 1841 def validSolution(self):
1109     """
1110     markes the solution as valid
1111     """
1112     self.__is_solution_valid=True
1113 jgs 149
1114 gross 1841 def invalidateSolution(self):
1115     """
1116     indicates the PDE has to be resolved if the solution is requested
1117     """
1118     self.trace("System will be resolved.")
1119     self.__is_solution_valid=False
1120 jgs 148
1121 gross 1841 def isSolutionValid(self):
1122     """
1123     returns True is the solution is still valid.
1124     """
1125     return self.__is_solution_valid
1126    
1127     def validOperator(self):
1128     """
1129     markes the operator as valid
1130     """
1131     self.__is_operator_valid=True
1132    
1133     def invalidateOperator(self):
1134     """
1135     indicates the operator has to be rebuilt next time it is used
1136     """
1137     self.trace("Operator will be rebuilt.")
1138     self.invalidateSolution()
1139     self.__is_operator_valid=False
1140    
1141     def isOperatorValid(self):
1142     """
1143     returns True is the operator is still valid.
1144     """
1145     return self.__is_operator_valid
1146    
1147     def validRightHandSide(self):
1148     """
1149     markes the right hand side as valid
1150     """
1151     self.__is_RHS_valid=True
1152    
1153     def invalidateRightHandSide(self):
1154     """
1155     indicates the right hand side has to be rebuild next time it is used
1156     """
1157     if self.isRightHandSideValid(): self.trace("Right hand side has to be rebuilt.")
1158     self.invalidateSolution()
1159     self.__is_RHS_valid=False
1160    
1161     def isRightHandSideValid(self):
1162     """
1163     returns True is the operator is still valid.
1164     """
1165     return self.__is_RHS_valid
1166    
1167     def invalidateSystem(self):
1168     """
1169     annonced that everthing has to be rebuild:
1170     """
1171     if self.isRightHandSideValid(): self.trace("System has to be rebuilt.")
1172     self.invalidateSolution()
1173     self.invalidateOperator()
1174     self.invalidateRightHandSide()
1175    
1176     def isSystemValid(self):
1177     """
1178     returns true if the system (including solution) is still vaild:
1179     """
1180     return self.isSolutionValid() and self.isOperatorValid() and self.isRightHandSideValid()
1181    
1182     def initializeSystem(self):
1183     """
1184     resets the system
1185     """
1186     self.trace("New System has been created.")
1187     self.__operator=escript.Operator()
1188     self.__righthandside=escript.Data()
1189     self.__solution=escript.Data()
1190     self.invalidateSystem()
1191    
1192     def getOperator(self):
1193     """
1194     returns the operator of the linear problem
1195    
1196     @return: the operator of the problem
1197     """
1198     return self.getSystem()[0]
1199    
1200     def getRightHandSide(self):
1201     """
1202     returns the right hand side of the linear problem
1203    
1204     @return: the right hand side of the problem
1205     @rtype: L{Data<escript.Data>}
1206     """
1207     return self.getSystem()[1]
1208    
1209     def createRightHandSide(self):
1210     """
1211     returns an instance of a new right hand side
1212     """
1213     self.trace("New right hand side is allocated.")
1214     if self.getNumEquations()>1:
1215     return escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)
1216     else:
1217     return escript.Data(0.,(),self.getFunctionSpaceForEquation(),True)
1218    
1219     def createSolution(self):
1220     """
1221     returns an instance of a new solution
1222     """
1223     self.trace("New solution is allocated.")
1224     if self.getNumSolutions()>1:
1225     return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)
1226     else:
1227     return escript.Data(0.,(),self.getFunctionSpaceForSolution(),True)
1228    
1229     def resetSolution(self):
1230     """
1231     sets the solution to zero
1232     """
1233     if self.__solution.isEmpty():
1234     self.__solution=self.createSolution()
1235     else:
1236     self.__solution.setToZero()
1237     self.trace("Solution is reset to zero.")
1238     def setSolution(self,u):
1239     """
1240     sets the solution assuming that makes the sytem valid.
1241     """
1242     self.__solution=u
1243     self.validSolution()
1244    
1245     def getCurrentSolution(self):
1246     """
1247     returns the solution in its current state
1248     """
1249     if self.__solution.isEmpty(): self.__solution=self.createSolution()
1250     return self.__solution
1251    
1252     def resetRightHandSide(self):
1253     """
1254     sets the right hand side to zero
1255     """
1256     if self.__righthandside.isEmpty():
1257     self.__righthandside=self.createRightHandSide()
1258     else:
1259     self.__righthandside.setToZero()
1260     self.trace("Right hand side is reset to zero.")
1261    
1262     def getCurrentRightHandSide(self):
1263     """
1264     returns the right hand side in its current state
1265     """
1266     if self.__righthandside.isEmpty(): self.__righthandside=self.createRightHandSide()
1267     return self.__righthandside
1268    
1269    
1270     def resetOperator(self):
1271     """
1272     makes sure that the operator is instantiated and returns it initialized by zeros
1273     """
1274     if self.__operator.isEmpty():
1275     if self.isUsingLumping():
1276     self.__operator=self.createSolution()
1277     else:
1278     self.__operator=self.createOperator()
1279     else:
1280     if self.isUsingLumping():
1281     self.__operator.setToZero()
1282     else:
1283     self.__operator.resetValues()
1284     self.trace("Operator reset to zero")
1285    
1286     def getCurrentOperator(self):
1287     """
1288     returns the operator in its current state
1289     """
1290     if self.__operator.isEmpty(): self.resetOperator()
1291     return self.__operator
1292    
1293     def setValue(self,**coefficients):
1294 jgs 148 """
1295     sets new values to coefficients
1296    
1297     @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1298     """
1299 jgs 108 # check if the coefficients are legal:
1300     for i in coefficients.iterkeys():
1301     if not self.hasCoefficient(i):
1302 jgs 148 raise IllegalCoefficient,"Attempt to set unknown coefficient %s"%i
1303 jgs 108 # if the number of unknowns or equations is still unknown we try to estimate them:
1304 jgs 148 if self.__numEquations==None or self.__numSolutions==None:
1305 jgs 108 for i,d in coefficients.iteritems():
1306     if hasattr(d,"shape"):
1307     s=d.shape
1308     elif hasattr(d,"getShape"):
1309     s=d.getShape()
1310     else:
1311     s=numarray.array(d).shape
1312     if s!=None:
1313     # get number of equations and number of unknowns:
1314 gross 1841 res=self.__COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(self.getDomain(),s)
1315 jgs 108 if res==None:
1316 jgs 148 raise IllegalCoefficientValue,"Illegal shape %s of coefficient %s"%(s,i)
1317 jgs 108 else:
1318 jgs 148 if self.__numEquations==None: self.__numEquations=res[0]
1319     if self.__numSolutions==None: self.__numSolutions=res[1]
1320     if self.__numEquations==None: raise UndefinedPDEError,"unidententified number of equations"
1321     if self.__numSolutions==None: raise UndefinedPDEError,"unidententified number of solutions"
1322 jgs 108 # now we check the shape of the coefficient if numEquations and numSolutions are set:
1323     for i,d in coefficients.iteritems():
1324 jgs 148 try:
1325 gross 1841 self.__COEFFICIENTS[i].setValue(self.getDomain(),
1326 gross 1072 self.getNumEquations(),self.getNumSolutions(),
1327     self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
1328     self.alteredCoefficient(i)
1329     except IllegalCoefficientFunctionSpace,m:
1330     # if the function space is wrong then we try the reduced version:
1331     i_red=i+"_reduced"
1332 gross 1841 if (not i_red in coefficients.keys()) and i_red in self.__COEFFICIENTS.keys():
1333 gross 1072 try:
1334 gross 1841 self.__COEFFICIENTS[i_red].setValue(self.getDomain(),
1335 gross 1072 self.getNumEquations(),self.getNumSolutions(),
1336     self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
1337     self.alteredCoefficient(i_red)
1338     except IllegalCoefficientValue,m:
1339     raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
1340     except IllegalCoefficientFunctionSpace,m:
1341     raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m))
1342     else:
1343     raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m))
1344 jgs 148 except IllegalCoefficientValue,m:
1345     raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
1346 jgs 149 self.__altered_coefficients=True
1347 jgs 108
1348 gross 1841 # ====================================================================================
1349     # methods that are typically overwritten when implementing a particular linear problem
1350     # ====================================================================================
1351     def getRequiredSystemType(self):
1352     """
1353     returns the system type which needs to be used by the current set up.
1354    
1355     @note: Typically this method is overwritten when implementing a particular linear problem.
1356     """
1357     return 0
1358    
1359     def createOperator(self):
1360 jgs 108 """
1361 gross 1841 returns an instance of a new operator
1362    
1363     @note: This method is overwritten when implementing a particular linear problem.
1364     """
1365     return escript.Operator()
1366    
1367     def checkSymmetry(self,verbose=True):
1368     """
1369     test the PDE for symmetry.
1370    
1371     @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.
1372     @type verbose: C{bool}
1373     @return: True if the problem is symmetric.
1374     @rtype: C{bool}
1375     @note: Typically this method is overwritten when implementing a particular linear problem.
1376     """
1377     out=True
1378     return out
1379    
1380     def getSolution(self,**options):
1381     """
1382     returns the solution of the problem
1383     @return: the solution
1384     @rtype: L{Data<escript.Data>}
1385    
1386     @note: This method is overwritten when implementing a particular linear problem.
1387     """
1388     return self.getCurrentSolution()
1389    
1390     def getSystem(self):
1391     """
1392 jgs 122 return the operator and right hand side of the PDE
1393 jgs 149
1394     @return: the discrete version of the PDE
1395     @rtype: C{tuple} of L{Operator,<escript.Operator>} and L{Data<escript.Data>}.
1396 gross 1841
1397     @note: This method is overwritten when implementing a particular linear problem.
1398 jgs 108 """
1399 gross 1841 return (self.getCurrentOperator(), self.getCurrentRightHandSide())
1400     #=====================
1401    
1402     class LinearPDE(LinearProblem):
1403     """
1404     This class is used to define a general linear, steady, second order PDE
1405     for an unknown function M{u} on a given domain defined through a L{Domain<escript.Domain>} object.
1406    
1407     For a single PDE with a solution with a single component the linear PDE is defined in the following form:
1408    
1409     M{-(grad(A[j,l]+A_reduced[j,l])*grad(u)[l]+(B[j]+B_reduced[j])u)[j]+(C[l]+C_reduced[l])*grad(u)[l]+(D+D_reduced)=-grad(X+X_reduced)[j,j]+(Y+Y_reduced)}
1410    
1411    
1412     where M{grad(F)} denotes the spatial derivative of M{F}. Einstein's summation convention,
1413     ie. summation over indexes appearing twice in a term of a sum is performed, is used.
1414     The coefficients M{A}, M{B}, M{C}, M{D}, M{X} and M{Y} have to be specified through L{Data<escript.Data>} objects in the
1415     L{Function<escript.Function>} and the coefficients M{A_reduced}, M{B_reduced}, M{C_reduced}, M{D_reduced}, M{X_reduced} and M{Y_reduced} have to be specified through L{Data<escript.Data>} objects in the
1416     L{ReducedFunction<escript.ReducedFunction>}. It is also allowd to use objects that can be converted into
1417     such L{Data<escript.Data>} objects. M{A} and M{A_reduced} are rank two, M{B}, M{C}, M{X},
1418     M{B_reduced}, M{C_reduced} and M{X_reduced} are rank one and M{D}, M{D_reduced} M{Y} and M{Y_reduced} are scalar.
1419    
1420     The following natural boundary conditions are considered:
1421    
1422     M{n[j]*((A[i,j]+A_reduced[i,j])*grad(u)[l]+(B+B_reduced)[j]*u)+(d+d_reduced)*u=n[j]*(X[j]+X_reduced[j])+y}
1423    
1424     where M{n} is the outer normal field. Notice that the coefficients M{A}, M{A_reduced}, M{B}, M{B_reduced}, M{X} and M{X_reduced} are defined in the PDE. The coefficients M{d} and M{y} and are each a scalar in the L{FunctionOnBoundary<escript.FunctionOnBoundary>} and the coefficients M{d_reduced} and M{y_reduced} and are each a scalar in the L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
1425    
1426    
1427     Constraints for the solution prescribing the value of the solution at certain locations in the domain. They have the form
1428    
1429     M{u=r} where M{q>0}
1430    
1431     M{r} and M{q} are each scalar where M{q} is the characteristic function defining where the constraint is applied.
1432     The constraints override any other condition set by the PDE or the boundary condition.
1433    
1434     The PDE is symmetrical if
1435    
1436     M{A[i,j]=A[j,i]} and M{B[j]=C[j]} and M{A_reduced[i,j]=A_reduced[j,i]} and M{B_reduced[j]=C_reduced[j]}
1437    
1438     For a system of PDEs and a solution with several components the PDE has the form
1439    
1440     M{-grad((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k])[j]+(C[i,k,l]+C_reduced[i,k,l])*grad(u[k])[l]+(D[i,k]+D_reduced[i,k]*u[k] =-grad(X[i,j]+X_reduced[i,j])[j]+Y[i]+Y_reduced[i] }
1441    
1442     M{A} and M{A_reduced} are of rank four, M{B}, M{B_reduced}, M{C} and M{C_reduced} are each of rank three, M{D}, M{D_reduced}, M{X_reduced} and M{X} are each a rank two and M{Y} and M{Y_reduced} are of rank one.
1443     The natural boundary conditions take the form:
1444    
1445     M{n[j]*((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k])+(d[i,k]+d_reduced[i,k])*u[k]=n[j]*(X[i,j]+X_reduced[i,j])+y[i]+y_reduced[i]}
1446    
1447    
1448     The coefficient M{d} is a rank two and M{y} is a rank one both in the L{FunctionOnBoundary<escript.FunctionOnBoundary>}. Constraints take the form and the coefficients M{d_reduced} is a rank two and M{y_reduced} is a rank one both in the L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
1449    
1450     Constraints take the form
1451    
1452     M{u[i]=r[i]} where M{q[i]>0}
1453    
1454     M{r} and M{q} are each rank one. Notice that at some locations not necessarily all components must have a constraint.
1455    
1456     The system of PDEs is symmetrical if
1457    
1458     - M{A[i,j,k,l]=A[k,l,i,j]}
1459     - M{A_reduced[i,j,k,l]=A_reduced[k,l,i,j]}
1460     - M{B[i,j,k]=C[k,i,j]}
1461     - M{B_reduced[i,j,k]=C_reduced[k,i,j]}
1462     - M{D[i,k]=D[i,k]}
1463     - M{D_reduced[i,k]=D_reduced[i,k]}
1464     - M{d[i,k]=d[k,i]}
1465     - M{d_reduced[i,k]=d_reduced[k,i]}
1466    
1467     L{LinearPDE} also supports solution discontinuities over a contact region in the domain. To specify the conditions across the
1468     discontinuity we are using the generalised flux M{J} which is in the case of a systems of PDEs and several components of the solution
1469     defined as
1470    
1471     M{J[i,j]=(A[i,j,k,l]+A_reduced[[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k]-X[i,j]-X_reduced[i,j]}
1472    
1473     For the case of single solution component and single PDE M{J} is defined
1474    
1475     M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[j]+(B[i]+B_reduced[i])*u-X[i]-X_reduced[i]}
1476    
1477     In the context of discontinuities M{n} denotes the normal on the discontinuity pointing from side 0 towards side 1
1478     calculated from L{getNormal<escript.FunctionSpace.getNormal>} of L{FunctionOnContactZero<escript.FunctionOnContactZero>}. For a system of PDEs
1479     the contact condition takes the form
1480    
1481     M{n[j]*J0[i,j]=n[j]*J1[i,j]=(y_contact[i]+y_contact_reduced[i])- (d_contact[i,k]+d_contact_reduced[i,k])*jump(u)[k]}
1482    
1483     where M{J0} and M{J1} are the fluxes on side 0 and side 1 of the discontinuity, respectively. M{jump(u)}, which is the difference
1484     of the solution at side 1 and at side 0, denotes the jump of M{u} across discontinuity along the normal calculated by
1485     L{jump<util.jump>}.
1486     The coefficient M{d_contact} is a rank two and M{y_contact} is a rank one both in the L{FunctionOnContactZero<escript.FunctionOnContactZero>} or L{FunctionOnContactOne<escript.FunctionOnContactOne>}.
1487     The coefficient M{d_contact_reduced} is a rank two and M{y_contact_reduced} is a rank one both in the L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.
1488     In case of a single PDE and a single component solution the contact condition takes the form
1489    
1490     M{n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)}
1491    
1492     In this case the coefficient M{d_contact} and M{y_contact} are each scalar both in the L{FunctionOnContactZero<escript.FunctionOnContactZero>} or L{FunctionOnContactOne<escript.FunctionOnContactOne>} and the coefficient M{d_contact_reduced} and M{y_contact_reduced} are each scalar both in the L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}
1493    
1494     typical usage:
1495    
1496     p=LinearPDE(dom)
1497     p.setValue(A=kronecker(dom),D=1,Y=0.5)
1498     u=p.getSolution()
1499    
1500     """
1501    
1502    
1503     def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):
1504     """
1505     initializes a new linear PDE
1506    
1507     @param domain: domain of the PDE
1508     @type domain: L{Domain<escript.Domain>}
1509     @param numEquations: number of equations. If numEquations==None the number of equations
1510     is extracted from the PDE coefficients.
1511     @param numSolutions: number of solution components. If numSolutions==None the number of solution components
1512     is extracted from the PDE coefficients.
1513     @param debug: if True debug informations are printed.
1514    
1515     """
1516     super(LinearPDE, self).__init__(domain,numEquations,numSolutions,debug)
1517     #
1518     # the coefficients of the PDE:
1519     #
1520     self.introduceCoefficients(
1521     A=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
1522     B=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1523     C=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
1524     D=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1525     X=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
1526     Y=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1527     d=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1528     y=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1529     d_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1530     y_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1531     A_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
1532     B_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1533     C_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
1534     D_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1535     X_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
1536     Y_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1537     d_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1538     y_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1539     d_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1540     y_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1541     r=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.RIGHTHANDSIDE),
1542     q=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.BOTH) )
1543    
1544     def __str__(self):
1545     """
1546     returns string representation of the PDE
1547    
1548     @return: a simple representation of the PDE
1549     @rtype: C{str}
1550     """
1551     return "<LinearPDE %d>"%id(self)
1552    
1553     def getRequiredSystemType(self):
1554     """
1555     returns the system type which needs to be used by the current set up.
1556     """
1557 gross 1859 return self.getDomain().getSystemMatrixTypeId(self.getSolverMethod()[0],self.getSolverMethod()[1],self.getSolverPackage(),self.isSymmetric())
1558 gross 1841
1559     def checkSymmetry(self,verbose=True):
1560     """
1561     test the PDE for symmetry.
1562    
1563     @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.
1564     @type verbose: C{bool}
1565     @return: True if the PDE is symmetric.
1566     @rtype: L{bool}
1567     @note: This is a very expensive operation. It should be used for degugging only! The symmetry flag is not altered.
1568     """
1569     out=True
1570     out=out and self.checkSymmetricTensor("A", verbose)
1571     out=out and self.checkSymmetricTensor("A_reduced", verbose)
1572     out=out and self.checkReciprocalSymmetry("B","C", verbose)
1573     out=out and self.checkReciprocalSymmetry("B_reduced","C_reduced", verbose)
1574     out=out and self.checkSymmetricTensor("D", verbose)
1575     out=out and self.checkSymmetricTensor("D_reduced", verbose)
1576     out=out and self.checkSymmetricTensor("d", verbose)
1577     out=out and self.checkSymmetricTensor("d_reduced", verbose)
1578     out=out and self.checkSymmetricTensor("d_contact", verbose)
1579     out=out and self.checkSymmetricTensor("d_contact_reduced", verbose)
1580     return out
1581    
1582     def createOperator(self):
1583     """
1584     returns an instance of a new operator
1585     """
1586     self.trace("New operator is allocated.")
1587     return self.getDomain().newOperator( \
1588     self.getNumEquations(), \
1589     self.getFunctionSpaceForEquation(), \
1590     self.getNumSolutions(), \
1591     self.getFunctionSpaceForSolution(), \
1592     self.getSystemType())
1593    
1594     def getSolution(self,**options):
1595     """
1596     returns the solution of the PDE
1597    
1598     @return: the solution
1599     @rtype: L{Data<escript.Data>}
1600     @param options: solver options
1601     @keyword verbose: True to get some information during PDE solution
1602     @type verbose: C{bool}
1603     @keyword reordering: reordering scheme to be used during elimination. Allowed values are
1604     L{NO_REORDERING}, L{MINIMUM_FILL_IN}, L{NESTED_DISSECTION}
1605     @keyword iter_max: maximum number of iteration steps allowed.
1606     @keyword drop_tolerance: threshold for drupping in L{ILUT}
1607     @keyword drop_storage: maximum of allowed memory in L{ILUT}
1608     @keyword truncation: maximum number of residuals in L{GMRES}
1609     @keyword restart: restart cycle length in L{GMRES}
1610     """
1611     if not self.isSolutionValid():
1612     mat,f=self.getSystem()
1613 jgs 108 if self.isUsingLumping():
1614 gross 1841 self.setSolution(f*1/mat)
1615     else:
1616     options[self.TOLERANCE_KEY]=self.getTolerance()
1617     options[self.METHOD_KEY]=self.getSolverMethod()[0]
1618     options[self.PRECONDITIONER_KEY]=self.getSolverMethod()[1]
1619     options[self.PACKAGE_KEY]=self.getSolverPackage()
1620     options[self.SYMMETRY_KEY]=self.isSymmetric()
1621     self.trace("PDE is resolved.")
1622     self.trace("solver options: %s"%str(options))
1623     self.setSolution(mat.solve(f,options))
1624     return self.getCurrentSolution()
1625    
1626     def getSystem(self):
1627     """
1628     return the operator and right hand side of the PDE
1629    
1630     @return: the discrete version of the PDE
1631     @rtype: C{tuple} of L{Operator,<escript.Operator>} and L{Data<escript.Data>}.
1632     """
1633     if not self.isOperatorValid() or not self.isRightHandSideValid():
1634     if self.isUsingLumping():
1635     if not self.isOperatorValid():
1636 gross 531 if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():
1637     raise TypeError,"Lumped matrix requires same order for equations and unknowns"
1638 gross 1841 if not self.getCoefficient("A").isEmpty():
1639 gross 531 raise ValueError,"coefficient A in lumped matrix may not be present."
1640 gross 1841 if not self.getCoefficient("B").isEmpty():
1641 gross 1072 raise ValueError,"coefficient B in lumped matrix may not be present."
1642 gross 1841 if not self.getCoefficient("C").isEmpty():
1643 gross 1072 raise ValueError,"coefficient C in lumped matrix may not be present."
1644 gross 1841 if not self.getCoefficient("d_contact").isEmpty():
1645 gross 1204 raise ValueError,"coefficient d_contact in lumped matrix may not be present."
1646 gross 1841 if not self.getCoefficient("A_reduced").isEmpty():
1647 gross 1072 raise ValueError,"coefficient A_reduced in lumped matrix may not be present."
1648 gross 1841 if not self.getCoefficient("B_reduced").isEmpty():
1649 gross 1072 raise ValueError,"coefficient B_reduced in lumped matrix may not be present."
1650 gross 1841 if not self.getCoefficient("C_reduced").isEmpty():
1651 gross 1072 raise ValueError,"coefficient C_reduced in lumped matrix may not be present."
1652 gross 1841 if not self.getCoefficient("d_contact_reduced").isEmpty():
1653 gross 1204 raise ValueError,"coefficient d_contact_reduced in lumped matrix may not be present."
1654 gross 1841 D=self.getCoefficient("D")
1655     d=self.getCoefficient("d")
1656     D_reduced=self.getCoefficient("D_reduced")
1657     d_reduced=self.getCoefficient("d_reduced")
1658 gross 531 if not D.isEmpty():
1659     if self.getNumSolutions()>1:
1660 bcumming 791 D_times_e=util.matrix_mult(D,numarray.ones((self.getNumSolutions(),)))
1661 gross 531 else:
1662     D_times_e=D
1663     else:
1664     D_times_e=escript.Data()
1665     if not d.isEmpty():
1666     if self.getNumSolutions()>1:
1667 bcumming 791 d_times_e=util.matrix_mult(d,numarray.ones((self.getNumSolutions(),)))
1668 gross 531 else:
1669     d_times_e=d
1670     else:
1671     d_times_e=escript.Data()
1672 gross 1204
1673 gross 1072 if not D_reduced.isEmpty():
1674     if self.getNumSolutions()>1:
1675     D_reduced_times_e=util.matrix_mult(D_reduced,numarray.ones((self.getNumSolutions(),)))
1676     else:
1677     D_reduced_times_e=D_reduced
1678     else:
1679     D_reduced_times_e=escript.Data()
1680     if not d_reduced.isEmpty():
1681     if self.getNumSolutions()>1:
1682     d_reduced_times_e=util.matrix_mult(d_reduced,numarray.ones((self.getNumSolutions(),)))
1683     else:
1684     d_reduced_times_e=d_reduced
1685     else:
1686     d_reduced_times_e=escript.Data()
1687 gross 1204
1688 gross 1841 self.resetOperator()
1689     operator=self.getCurrentOperator()
1690 gross 1659 if False and hasattr(self.getDomain(), "addPDEToLumpedSystem") :
1691 gross 1841 self.getDomain().addPDEToLumpedSystem(operator, D_times_e, d_times_e)
1692     self.getDomain().addPDEToLumpedSystem(operator, D_reduced_times_e, d_reduced_times_e)
1693 gross 1072 else:
1694 gross 1841 self.getDomain().addPDEToRHS(operator, \
1695 gross 1204 escript.Data(), \
1696     D_times_e, \
1697     d_times_e,\
1698     escript.Data())
1699 gross 1841 self.getDomain().addPDEToRHS(operator, \
1700 gross 1204 escript.Data(), \
1701     D_reduced_times_e, \
1702     d_reduced_times_e,\
1703     escript.Data())
1704 jgs 148 self.trace("New lumped operator has been built.")
1705 gross 1841 if not self.isRightHandSideValid():
1706     self.resetRightHandSide()
1707     righthandside=self.getCurrentRightHandSide()
1708     self.getDomain().addPDEToRHS(righthandside, \
1709     self.getCoefficient("X"), \
1710     self.getCoefficient("Y"),\
1711     self.getCoefficient("y"),\
1712     self.getCoefficient("y_contact"))
1713     self.getDomain().addPDEToRHS(righthandside, \
1714     self.getCoefficient("X_reduced"), \
1715     self.getCoefficient("Y_reduced"),\
1716     self.getCoefficient("y_reduced"),\
1717     self.getCoefficient("y_contact_reduced"))
1718 jgs 148 self.trace("New right hand side as been built.")
1719 gross 1841 self.validRightHandSide()
1720     self.insertConstraint(rhs_only=False)
1721     self.validOperator()
1722 jgs 108 else:
1723 gross 1841 if not self.isOperatorValid() and not self.isRightHandSideValid():
1724     self.resetRightHandSide()
1725     righthandside=self.getCurrentRightHandSide()
1726     self.resetOperator()
1727     operator=self.getCurrentOperator()
1728     self.getDomain().addPDEToSystem(operator,righthandside, \
1729     self.getCoefficient("A"), \
1730     self.getCoefficient("B"), \
1731     self.getCoefficient("C"), \
1732     self.getCoefficient("D"), \
1733     self.getCoefficient("X"), \
1734     self.getCoefficient("Y"), \
1735     self.getCoefficient("d"), \
1736     self.getCoefficient("y"), \
1737     self.getCoefficient("d_contact"), \
1738     self.getCoefficient("y_contact"))
1739     self.getDomain().addPDEToSystem(operator,righthandside, \
1740     self.getCoefficient("A_reduced"), \
1741     self.getCoefficient("B_reduced"), \
1742     self.getCoefficient("C_reduced"), \
1743     self.getCoefficient("D_reduced"), \
1744     self.getCoefficient("X_reduced"), \
1745     self.getCoefficient("Y_reduced"), \
1746     self.getCoefficient("d_reduced"), \
1747     self.getCoefficient("y_reduced"), \
1748     self.getCoefficient("d_contact_reduced"), \
1749     self.getCoefficient("y_contact_reduced"))
1750     self.insertConstraint(rhs_only=False)
1751 jgs 148 self.trace("New system has been built.")
1752 gross 1841 self.validOperator()
1753     self.validRightHandSide()
1754     elif not self.isRightHandSideValid():
1755     self.resetRightHandSide()
1756     righthandside=self.getCurrentRightHandSide()
1757     self.getDomain().addPDEToRHS(righthandside,
1758     self.getCoefficient("X"), \
1759     self.getCoefficient("Y"),\
1760     self.getCoefficient("y"),\
1761     self.getCoefficient("y_contact"))
1762     self.getDomain().addPDEToRHS(righthandside,
1763     self.getCoefficient("X_reduced"), \
1764     self.getCoefficient("Y_reduced"),\
1765     self.getCoefficient("y_reduced"),\
1766     self.getCoefficient("y_contact_reduced"))
1767     self.insertConstraint(rhs_only=True)
1768 jgs 148 self.trace("New right hand side has been built.")
1769 gross 1841 self.validRightHandSide()
1770     elif not self.isOperatorValid():
1771     self.resetOperator()
1772     operator=self.getCurrentOperator()
1773     self.getDomain().addPDEToSystem(operator,escript.Data(), \
1774     self.getCoefficient("A"), \
1775     self.getCoefficient("B"), \
1776     self.getCoefficient("C"), \
1777     self.getCoefficient("D"), \
1778 jgs 108 escript.Data(), \
1779     escript.Data(), \
1780 gross 1841 self.getCoefficient("d"), \
1781 jgs 108 escript.Data(),\
1782 gross 1841 self.getCoefficient("d_contact"), \
1783 jgs 108 escript.Data())
1784 gross 1841 self.getDomain().addPDEToSystem(operator,escript.Data(), \
1785     self.getCoefficient("A_reduced"), \
1786     self.getCoefficient("B_reduced"), \
1787     self.getCoefficient("C_reduced"), \
1788     self.getCoefficient("D_reduced"), \
1789 gross 1072 escript.Data(), \
1790     escript.Data(), \
1791 gross 1841 self.getCoefficient("d_reduced"), \
1792 gross 1072 escript.Data(),\
1793 gross 1841 self.getCoefficient("d_contact_reduced"), \
1794 gross 1072 escript.Data())
1795 gross 1841 self.insertConstraint(rhs_only=False)
1796 jgs 148 self.trace("New operator has been built.")
1797 gross 1841 self.validOperator()
1798     return (self.getCurrentOperator(), self.getCurrentRightHandSide())
1799 jgs 102
1800 gross 1841 def insertConstraint(self, rhs_only=False):
1801     """
1802     applies the constraints defined by q and r to the PDE
1803    
1804     @param rhs_only: if True only the right hand side is altered by the constraint.
1805     @type rhs_only: C{bool}
1806     """
1807     q=self.getCoefficient("q")
1808     r=self.getCoefficient("r")
1809     righthandside=self.getCurrentRightHandSide()
1810     operator=self.getCurrentOperator()
1811 jgs 102
1812 gross 1841 if not q.isEmpty():
1813     if r.isEmpty():
1814     r_s=self.createSolution()
1815     else:
1816     r_s=r
1817     if not rhs_only and not operator.isEmpty():
1818     if self.isUsingLumping():
1819     operator.copyWithMask(escript.Data(1.,q.getShape(),q.getFunctionSpace()),q)
1820     else:
1821     row_q=escript.Data(q,self.getFunctionSpaceForEquation())
1822     col_q=escript.Data(q,self.getFunctionSpaceForSolution())
1823     u=self.createSolution()
1824     u.copyWithMask(r_s,col_q)
1825     righthandside-=operator*u
1826     operator.nullifyRowsAndCols(row_q,col_q,1.)
1827     righthandside.copyWithMask(r_s,q)
1828    
1829     def setValue(self,**coefficients):
1830     """
1831     sets new values to coefficients
1832    
1833     @param coefficients: new values assigned to coefficients
1834     @keyword A: value for coefficient A.
1835     @type A: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1836     @keyword A_reduced: value for coefficient A_reduced.
1837     @type A_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1838     @keyword B: value for coefficient B
1839     @type B: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1840     @keyword B_reduced: value for coefficient B_reduced
1841     @type B_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1842     @keyword C: value for coefficient C
1843     @type C: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1844     @keyword C_reduced: value for coefficient C_reduced
1845     @type C_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1846     @keyword D: value for coefficient D
1847     @type D: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1848     @keyword D_reduced: value for coefficient D_reduced
1849     @type D_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1850     @keyword X: value for coefficient X
1851     @type X: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1852     @keyword X_reduced: value for coefficient X_reduced
1853     @type X_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1854     @keyword Y: value for coefficient Y
1855     @type Y: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1856     @keyword Y_reduced: value for coefficient Y_reduced
1857     @type Y_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.Function>}.
1858     @keyword d: value for coefficient d
1859     @type d: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1860     @keyword d_reduced: value for coefficient d_reduced
1861     @type d_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
1862     @keyword y: value for coefficient y
1863     @type y: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1864     @keyword d_contact: value for coefficient d_contact
1865     @type d_contact: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnContactOne<escript.FunctionOnContactOne>} or L{FunctionOnContactZero<escript.FunctionOnContactZero>}.
1866     @keyword d_contact_reduced: value for coefficient d_contact_reduced
1867     @type d_contact_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>} or L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>}.
1868     @keyword y_contact: value for coefficient y_contact
1869     @type y_contact: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnContactOne<escript.FunctionOnContactOne>} or L{FunctionOnContactZero<escript.FunctionOnContactZero>}.
1870     @keyword y_contact_reduced: value for coefficient y_contact_reduced
1871     @type y_contact_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunctionOnContactOne<escript.FunctionOnContactOne>} or L{ReducedFunctionOnContactZero<escript.FunctionOnContactZero>}.
1872     @keyword r: values prescribed to the solution at the locations of constraints
1873     @type r: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1874     depending of reduced order is used for the solution.
1875     @keyword q: mask for location of constraints
1876     @type q: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1877     depending of reduced order is used for the representation of the equation.
1878     @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1879     """
1880     super(LinearPDE,self).setValue(**coefficients)
1881     # check if the systrem is inhomogeneous:
1882     if len(coefficients)>0 and not self.isUsingLumping():
1883     q=self.getCoefficient("q")
1884     r=self.getCoefficient("r")
1885     if not q.isEmpty() and not r.isEmpty():
1886     if util.Lsup(q*r)>0.:
1887     self.trace("Inhomogeneous constraint detected.")
1888     self.invalidateSystem()
1889    
1890    
1891     def getResidual(self,u=None):
1892     """
1893     return the residual of u or the current solution if u is not present.
1894    
1895     @param u: argument in the residual calculation. It must be representable in C{elf.getFunctionSpaceForSolution()}. If u is not present or equals L{None}
1896     the current solution is used.
1897     @type u: L{Data<escript.Data>} or None
1898     @return: residual of u
1899     @rtype: L{Data<escript.Data>}
1900     """
1901     if u==None:
1902     return self.getOperator()*self.getSolution()-self.getRightHandSide()
1903     else:
1904     return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())-self.getRightHandSide()
1905    
1906     def getFlux(self,u=None):
1907     """
1908     returns the flux M{J} for a given M{u}
1909    
1910     M{J[i,j]=(A[i,j,k,l]+A_reduced[A[i,j,k,l]]*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])u[k]-X[i,j]-X_reduced[i,j]}
1911    
1912     or
1913    
1914     M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[l]+(B[j]+B_reduced[j])u-X[j]-X_reduced[j]}
1915    
1916     @param u: argument in the flux. If u is not present or equals L{None} the current solution is used.
1917     @type u: L{Data<escript.Data>} or None
1918     @return: flux
1919     @rtype: L{Data<escript.Data>}
1920     """
1921     if u==None: u=self.getSolution()
1922     return util.tensormult(self.getCoefficient("A"),util.grad(u,Funtion(self.getDomain))) \
1923     +util.matrixmult(self.getCoefficient("B"),u) \
1924     -util.self.getCoefficient("X") \
1925     +util.tensormult(self.getCoefficient("A_reduced"),util.grad(u,ReducedFuntion(self.getDomain))) \
1926     +util.matrixmult(self.getCoefficient("B_reduced"),u) \
1927     -util.self.getCoefficient("X_reduced")
1928    
1929    
1930 jgs 149 class Poisson(LinearPDE):
1931     """
1932     Class to define a Poisson equation problem, which is genear L{LinearPDE} of the form
1933 jgs 102
1934 jgs 149 M{-grad(grad(u)[j])[j] = f}
1935 jgs 102
1936 jgs 149 with natural boundary conditons
1937    
1938     M{n[j]*grad(u)[j] = 0 }
1939    
1940     and constraints:
1941    
1942     M{u=0} where M{q>0}
1943    
1944 jgs 108 """
1945 jgs 148
1946 gross 328 def __init__(self,domain,debug=False):
1947 jgs 149 """
1948     initializes a new Poisson equation
1949 jgs 104
1950 jgs 149 @param domain: domain of the PDE
1951     @type domain: L{Domain<escript.Domain>}
1952     @param debug: if True debug informations are printed.
1953 jgs 102
1954 jgs 149 """
1955 jgs 151 super(Poisson, self).__init__(domain,1,1,debug)
1956 gross 1841 self.introduceCoefficients(
1957     f=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1958     f_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE))
1959 jgs 149 self.setSymmetryOn()
1960 jgs 102
1961 jgs 149 def setValue(self,**coefficients):
1962     """
1963     sets new values to coefficients
1964 jgs 102
1965 jgs 149 @param coefficients: new values assigned to coefficients
1966     @keyword f: value for right hand side M{f}
1967     @type f: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1968     @keyword q: mask for location of constraints
1969     @type q: any type that can be casted to rank zeo L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1970     depending of reduced order is used for the representation of the equation.
1971     @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1972     """
1973 jgs 151 super(Poisson, self).setValue(**coefficients)
1974 jgs 102
1975 gross 1841
1976     def getCoefficient(self,name):
1977 jgs 149 """
1978     return the value of the coefficient name of the general PDE
1979     @param name: name of the coefficient requested.
1980     @type name: C{string}
1981     @return: the value of the coefficient name
1982     @rtype: L{Data<escript.Data>}
1983 gross 1841 @raise IllegalCoefficient: invalid coefficient name
1984 jgs 149 @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
1985     """
1986     if name == "A" :
1987 gross 328 return escript.Data(util.kronecker(self.getDim()),escript.Function(self.getDomain()))
1988 jgs 149 elif name == "Y" :
1989     return self.getCoefficient("f")
1990 gross 1072 elif name == "Y_reduced" :
1991     return self.getCoefficient("f_reduced")
1992 jgs 149 else:
1993 gross 1841 return super(Poisson, self).getCoefficient(name)
1994 jgs 102
1995 jgs 149 class Helmholtz(LinearPDE):
1996     """
1997     Class to define a Helmhotz equation problem, which is genear L{LinearPDE} of the form
1998    
1999     M{S{omega}*u - grad(k*grad(u)[j])[j] = f}
2000    
2001     with natural boundary conditons
2002    
2003     M{k*n[j]*grad(u)[j] = g- S{alpha}u }
2004    
2005 jgs 122 and constraints:
2006 jgs 102
2007 jgs 149 M{u=r} where M{q>0}
2008    
2009 jgs 110 """
2010 jgs 149
2011     def __init__(self,domain,debug=False):
2012     """
2013     initializes a new Poisson equation
2014    
2015     @param domain: domain of the PDE
2016     @type domain: L{Domain<escript.Domain>}
2017     @param debug: if True debug informations are printed.
2018    
2019     """
2020 jgs 151 super(Helmholtz, self).__init__(domain,1,1,debug)
2021 gross 1841 self.introduceCoefficients(
2022     omega=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.OPERATOR),
2023     k=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.OPERATOR),
2024     f=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2025     f_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2026     alpha=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.OPERATOR),
2027     g=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2028     g_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE))
2029 jgs 149 self.setSymmetryOn()
2030    
2031 gross 1841 def setValue(self,**coefficients):
2032 jgs 149 """
2033     sets new values to coefficients
2034    
2035     @param coefficients: new values assigned to coefficients
2036     @keyword omega: value for coefficient M{S{omega}}
2037     @type omega: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2038     @keyword k: value for coefficeint M{k}
2039     @type k: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2040     @keyword f: value for right hand side M{f}
2041     @type f: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2042     @keyword alpha: value for right hand side M{S{alpha}}
2043     @type alpha: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
2044     @keyword g: value for right hand side M{g}
2045     @type g: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
2046     @keyword r: prescribed values M{r} for the solution in constraints.
2047     @type r: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2048     depending of reduced order is used for the representation of the equation.
2049     @keyword q: mask for location of constraints
2050     @type q: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2051     depending of reduced order is used for the representation of the equation.
2052     @raise IllegalCoefficient: if an unknown coefficient keyword is used.
2053     """
2054 jgs 151 super(Helmholtz, self).setValue(**coefficients)
2055 jgs 149
2056 gross 1841 def getCoefficient(self,name):
2057 jgs 149 """
2058     return the value of the coefficient name of the general PDE
2059    
2060     @param name: name of the coefficient requested.
2061     @type name: C{string}
2062     @return: the value of the coefficient name
2063     @rtype: L{Data<escript.Data>}
2064 gross 1841 @raise IllegalCoefficient: invalid name
2065 jgs 149 """
2066     if name == "A" :
2067     return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))*self.getCoefficient("k")
2068     elif name == "D" :
2069     return self.getCoefficient("omega")
2070     elif name == "Y" :
2071     return self.getCoefficient("f")
2072     elif name == "d" :
2073     return self.getCoefficient("alpha")
2074     elif name == "y" :
2075     return self.getCoefficient("g")
2076 gross 1072 elif name == "Y_reduced" :
2077     return self.getCoefficient("f_reduced")
2078     elif name == "y_reduced" :
2079     return self.getCoefficient("g_reduced")
2080 jgs 149 else:
2081 gross 1841 return super(Helmholtz, self).getCoefficient(name)
2082 jgs 149
2083     class LameEquation(LinearPDE):
2084     """
2085     Class to define a Lame equation problem:
2086    
2087 gross 657 M{-grad(S{mu}*(grad(u[i])[j]+grad(u[j])[i]))[j] - grad(S{lambda}*grad(u[k])[k])[j] = F_i -grad(S{sigma}[ij])[j] }
2088 jgs 149
2089     with natural boundary conditons:
2090    
2091 gross 657 M{n[j]*(S{mu}*(grad(u[i])[j]+grad(u[j])[i]) + S{lambda}*grad(u[k])[k]) = f_i +n[j]*S{sigma}[ij] }
2092 jgs 149
2093     and constraints:
2094    
2095     M{u[i]=r[i]} where M{q[i]>0}
2096    
2097     """
2098    
2099     def __init__(self,domain,debug=False):
2100 jgs 151 super(LameEquation, self).__init__(domain,\
2101     domain.getDim(),domain.getDim(),debug)
2102 gross 1841 self.introduceCoefficients(lame_lambda=PDECoef(PDECoef.INTERIOR,(),PDECoef.OPERATOR),
2103     lame_mu=PDECoef(PDECoef.INTERIOR,(),PDECoef.OPERATOR),
2104     F=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2105     sigma=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
2106     f=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE))
2107 jgs 151 self.setSymmetryOn()
2108 jgs 149
2109 gross 1841 def setValues(self,**coefficients):
2110 jgs 149 """
2111     sets new values to coefficients
2112    
2113     @param coefficients: new values assigned to coefficients
2114     @keyword lame_mu: value for coefficient M{S{mu}}
2115     @type lame_mu: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2116     @keyword lame_lambda: value for coefficient M{S{lambda}}
2117     @type lame_lambda: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2118     @keyword F: value for internal force M{F}
2119     @type F: any type that can be casted to L{Vector<escript.Vector>} object on L{Function<escript.Function>}
2120     @keyword sigma: value for initial stress M{S{sigma}}
2121     @type sigma: any type that can be casted to L{Tensor<escript.Tensor>} object on L{Function<escript.Function>}
2122     @keyword f: value for extrenal force M{f}
2123     @type f: any type that can be casted to L{Vector<escript.Vector>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}
2124     @keyword r: prescribed values M{r} for the solution in constraints.
2125     @type r: any type that can be casted to L{Vector<escript.Vector>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2126     depending of reduced order is used for the representation of the equation.
2127     @keyword q: mask for location of constraints
2128     @type q: any type that can be casted to L{Vector<escript.Vector>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2129     depending of reduced order is used for the representation of the equation.
2130     @raise IllegalCoefficient: if an unknown coefficient keyword is used.
2131     """
2132 gross 596 super(LameEquation, self).setValues(**coefficients)
2133 jgs 149
2134 gross 1841 def getCoefficient(self,name):
2135 jgs 149 """
2136     return the value of the coefficient name of the general PDE
2137    
2138     @param name: name of the coefficient requested.
2139     @type name: C{string}
2140     @return: the value of the coefficient name
2141     @rtype: L{Data<escript.Data>}
2142 gross 1841 @raise IllegalCoefficient: invalid coefficient name
2143 jgs 149 """
2144     if name == "A" :
2145 gross 1841 out =self.createCoefficient("A")
2146 jgs 149 for i in range(self.getDim()):
2147     for j in range(self.getDim()):
2148     out[i,i,j,j] += self.getCoefficient("lame_lambda")
2149     out[i,j,j,i] += self.getCoefficient("lame_mu")
2150     out[i,j,i,j] += self.getCoefficient("lame_mu")
2151     return out
2152     elif name == "X" :
2153     return self.getCoefficient("sigma")
2154     elif name == "Y" :
2155     return self.getCoefficient("F")
2156     elif name == "y" :
2157     return self.getCoefficient("f")
2158     else:
2159 gross 1841 return super(LameEquation, self).getCoefficient(name)
2160 jgs 149
2161 gross 1400 def LinearSinglePDE(domain,debug=False):
2162     """
2163     defines a single linear PDEs
2164    
2165     @param domain: domain of the PDE
2166     @type domain: L{Domain<escript.Domain>}
2167     @param debug: if True debug informations are printed.
2168     @rtype: L{LinearPDE}
2169     """
2170     return LinearPDE(domain,numEquations=1,numSolutions=1,debug=debug)
2171    
2172     def LinearPDESystem(domain,debug=False):
2173     """
2174     defines a system of linear PDEs
2175    
2176     @param domain: domain of the PDE
2177     @type domain: L{Domain<escript.Domain>}
2178     @param debug: if True debug informations are printed.
2179     @rtype: L{LinearPDE}
2180     """
2181     return LinearPDE(domain,numEquations=domain.getDim(),numSolutions=domain.getDim(),debug=debug)
2182 gross 1410
2183 gross 1841 class TransportPDE(LinearProblem):
2184     """
2185     This class is used to define a transport problem given by a general linear, time dependend, second order PDE
2186     for an unknown, non-negative, time-dependend function M{u} on a given domain defined through a L{Domain<escript.Domain>} object.
2187    
2188     For a single equation with a solution with a single component the transport problem is defined in the following form:
2189    
2190     M{(M+M_reduced)*u_t=-(grad(A[j,l]+A_reduced[j,l])*grad(u)[l]+(B[j]+B_reduced[j])u)[j]+(C[l]+C_reduced[l])*grad(u)[l]+(D+D_reduced)-grad(X+X_reduced)[j,j]+(Y+Y_reduced)}
2191    
2192    
2193     where M{u_t} denotes the time derivative of M{u} and M{grad(F)} denotes the spatial derivative of M{F}. Einstein's summation convention,
2194     ie. summation over indexes appearing twice in a term of a sum is performed, is used.
2195     The coefficients M{M}, M{A}, M{B}, M{C}, M{D}, M{X} and M{Y} have to be specified through L{Data<escript.Data>} objects in the
2196     L{Function<escript.Function>} and the coefficients M{M_reduced}, M{A_reduced}, M{B_reduced}, M{C_reduced}, M{D_reduced}, M{X_reduced} and M{Y_reduced} have to be specified through L{Data<escript.Data>} objects in the
2197     L{ReducedFunction<escript.ReducedFunction>}. It is also allowed to use objects that can be converted into
2198     such L{Data<escript.Data>} objects. M{M} and M{M_reduced} are scalar, M{A} and M{A_reduced} are rank two,
2199     M{B}, M{C}, M{X}, M{B_reduced}, M{C_reduced} and M{X_reduced} are rank one and M{D}, M{D_reduced}, M{Y} and M{Y_reduced} are scalar.
2200    
2201     The following natural boundary conditions are considered:
2202    
2203     M{n[j]*((A[i,j]+A_reduced[i,j])*grad(u)[l]+(B+B_reduced)[j]*u+X[j]+X_reduced[j])+(d+d_reduced)*u+y+y_reduced=(m+m_reduced)*u_t}
2204    
2205     where M{n} is the outer normal field. Notice that the coefficients M{A}, M{A_reduced}, M{B}, M{B_reduced}, M{X} and M{X_reduced} are defined in the transport problem. The coefficients M{m}, M{d} and M{y} and are each a scalar in the L{FunctionOnBoundary<escript.FunctionOnBoundary>} and the coefficients M{m_reduced}, M{d_reduced} and M{y_reduced} and are each a scalar in the L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
2206    
2207     Constraints for the solution prescribing the value of the solution at certain locations in the domain. They have the form
2208    
2209     M{u_t=r} where M{q>0}
2210    
2211     M{r} and M{q} are each scalar where M{q} is the characteristic function defining where the constraint is applied.
2212     The constraints override any other condition set by the transport problem or the boundary condition.
2213    
2214     The transport problem is symmetrical if
2215    
2216     M{A[i,j]=A[j,i]} and M{B[j]=C[j]} and M{A_reduced[i,j]=A_reduced[j,i]} and M{B_reduced[j]=C_reduced[j]}
2217    
2218     For a system and a solution with several components the transport problem has the form
2219    
2220     M{(M[i,k]+M_reduced[i,k])*u[k]_t=-grad((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k])[j]+(C[i,k,l]+C_reduced[i,k,l])*grad(u[k])[l]+(D[i,k]+D_reduced[i,k]*u[k]-grad(X[i,j]+X_reduced[i,j])[j]+Y[i]+Y_reduced[i] }
2221    
2222     M{A} and M{A_reduced} are of rank four, M{B}, M{B_reduced}, M{C} and M{C_reduced} are each of rank three, M{M}, M{M_reduced}, M{D}, M{D_reduced}, M{X_reduced} and M{X} are each a rank two and M{Y} and M{Y_reduced} are of rank one. The natural boundary conditions take the form:
2223    
2224     M{n[j]*((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k]+X[i,j]+X_reduced[i,j])+(d[i,k]+d_reduced[i,k])*u[k]+y[i]+y_reduced[i]= (m[i,k]+m_reduced[i,k])*u[k]_t}
2225    
2226     The coefficient M{d} and M{m} are of rank two and M{y} are of rank one with all in the L{FunctionOnBoundary<escript.FunctionOnBoundary>}. Constraints take the form and the coefficients M{d_reduced} and M{m_reduced} are of rank two and M{y_reduced} is of rank one with all in the L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
2227    
2228     Constraints take the form
2229    
2230     M{u[i]_t=r[i]} where M{q[i]>0}
2231    
2232     M{r} and M{q} are each rank one. Notice that at some locations not necessarily all components must have a constraint.
2233    
2234     The transport problem is symmetrical if
2235    
2236     - M{M[i,k]=M[i,k]}
2237     - M{M_reduced[i,k]=M_reduced[i,k]}
2238     - M{A[i,j,k,l]=A[k,l,i,j]}
2239     - M{A_reduced[i,j,k,l]=A_reduced[k,l,i,j]}
2240     - M{B[i,j,k]=C[k,i,j]}
2241     - M{B_reduced[i,j,k]=C_reduced[k,i,j]}
2242     - M{D[i,k]=D[i,k]}
2243     - M{D_reduced[i,k]=D_reduced[i,k]}
2244     - M{m[i,k]=m[k,i]}
2245     - M{m_reduced[i,k]=m_reduced[k,i]}
2246     - M{d[i,k]=d[k,i]}
2247     - M{d_reduced[i,k]=d_reduced[k,i]}
2248    
2249     L{TransportPDE} also supports solution discontinuities over a contact region in the domain. To specify the conditions across the
2250     discontinuity we are using the generalised flux M{J} which is in the case of a systems of PDEs and several components of the solution
2251     defined as
2252    
2253     M{J[i,j]=(A[i,j,k,l]+A_reduced[[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k]+X[i,j]+X_reduced[i,j]}
2254    
2255     For the case of single solution component and single PDE M{J} is defined
2256    
2257     M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[j]+(B[i]+B_reduced[i])*u+X[i]+X_reduced[i]}
2258    
2259     In the context of discontinuities M{n} denotes the normal on the discontinuity pointing from side 0 towards side 1
2260     calculated from L{getNormal<escript.FunctionSpace.getNormal>} of L{FunctionOnContactZero<escript.FunctionOnContactZero>}. For a system of transport problems the contact condition takes the form
2261    
2262     M{n[j]*J0[i,j]=n[j]*J1[i,j]=(y_contact[i]+y_contact_reduced[i])- (d_contact[i,k]+d_contact_reduced[i,k])*jump(u)[k]}
2263    
2264     where M{J0} and M{J1} are the fluxes on side 0 and side 1 of the discontinuity, respectively. M{jump(u)}, which is the difference
2265     of the solution at side 1 and at side 0, denotes the jump of M{u} across discontinuity along the normal calculated by
2266     L{jump<util.jump>}.
2267     The coefficient M{d_contact} is a rank two and M{y_contact} is a rank one both in the L{FunctionOnContactZero<escript.FunctionOnContactZero>} or L{FunctionOnContactOne<escript.FunctionOnContactOne>}.
2268     The coefficient M{d_contact_reduced} is a rank two and M{y_contact_reduced} is a rank one both in the L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.
2269     In case of a single PDE and a single component solution the contact condition takes the form
2270    
2271     M{n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)}
2272    
2273     In this case the coefficient M{d_contact} and M{y_contact} are each scalar both in the L{FunctionOnContactZero<escript.FunctionOnContactZero>} or L{FunctionOnContactOne<escript.FunctionOnContactOne>} and the coefficient M{d_contact_reduced} and M{y_contact_reduced} are each scalar both in the L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}
2274    
2275     typical usage:
2276    
2277     p=TransportPDE(dom)
2278     p.setValue(M=1.,C=[-1.,0.])
2279     p.setInitialSolution(u=exp(-length(dom.getX()-[0.1,0.1])**2)
2280     t=0
2281     dt=0.1
2282     while (t<1.):
2283     u=p.solve(dt)
2284    
2285     """
2286 gross 1859 def __init__(self,domain,numEquations=None,numSolutions=None, theta=0.5,debug=False):
2287 gross 1410 """
2288 gross 1841 initializes a linear problem
2289 gross 1410
2290 gross 1841 @param domain: domain of the PDE
2291     @type domain: L{Domain<escript.Domain>}
2292     @param numEquations: number of equations. If numEquations==None the number of equations
2293     is extracted from the coefficients.
2294     @param numSolutions: number of solution components. If numSolutions==None the number of solution components
2295     is extracted from the coefficients.
2296     @param debug: if True debug informations are printed.
2297     @param theta: time stepping control: =1 backward Euler, =0 forward Euler, =0.5 Crank-Nicholson scheme, U{see here<http://en.wikipedia.org/wiki/Crank-Nicolson_method>}.
2298 gross 1410
2299 gross 1841 """
2300     if theta<0 or theta>1:
2301 gross 1859 raise ValueError,"theta (=%s) needs to be between 0 and 1 (inclusive)."%theta
2302     super(TransportPDE, self).__init__(domain,numEquations,numSolutions,debug)
2303 gross 1410
2304 gross 1841 self.__theta=theta
2305     #
2306     # the coefficients of the transport problem
2307     #
2308     self.introduceCoefficients(
2309     M=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2310     A=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2311     B=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2312     C=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2313     D=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2314     X=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
2315     Y=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2316     m=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2317     d=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2318     y=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2319     d_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2320     y_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2321     M_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2322     A_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2323     B_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2324     C_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2325     D_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2326     X_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
2327     Y_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2328     m_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2329     d_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2330     y_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2331     d_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2332     y_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2333     r=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.RIGHTHANDSIDE),
2334     q=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.BOTH) )
2335 gross 1410
2336 gross 1841 def __str__(self):
2337 gross 1410 """
2338 gross 1841 returns string representation of the transport problem
2339 gross 1410
2340 gross 1841 @return: a simple representation of the transport problem
2341     @rtype: C{str}
2342     """
2343     return "<TransportPDE %d>"%id(self)
2344 gross 1417
2345 gross 1841 def getTheta(self):
2346     """
2347     returns the time stepping control parameter
2348     @rtype: float
2349     """
2350     return self.__theta
2351 gross 1410
2352 gross 1859
2353 gross 1841 def checkSymmetry(self,verbose=True):
2354     """
2355     test the transport problem for symmetry.
2356    
2357     @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.
2358     @type verbose: C{bool}
2359     @return: True if the PDE is symmetric.
2360     @rtype: C{bool}
2361     @note: This is a very expensive operation. It should be used for degugging only! The symmetry flag is not altered.
2362     """
2363     out=True
2364     out=out and self.checkSymmetricTensor("M", verbose)
2365     out=out and self.checkSymmetricTensor("M_reduced", verbose)
2366     out=out and self.checkSymmetricTensor("A", verbose)
2367     out=out and self.checkSymmetricTensor("A_reduced", verbose)
2368     out=out and self.checkReciprocalSymmetry("B","C", verbose)
2369     out=out and self.checkReciprocalSymmetry("B_reduced","C_reduced", verbose)
2370     out=out and self.checkSymmetricTensor("D", verbose)
2371     out=out and self.checkSymmetricTensor("D_reduced", verbose)
2372     out=out and self.checkSymmetricTensor("m", verbose)
2373     out=out and self.checkSymmetricTensor("m_reduced", verbose)
2374     out=out and self.checkSymmetricTensor("d", verbose)
2375     out=out and self.checkSymmetricTensor("d_reduced", verbose)
2376     out=out and self.checkSymmetricTensor("d_contact", verbose)
2377     out=out and self.checkSymmetricTensor("d_contact_reduced", verbose)
2378     return out
2379    
2380     def setValue(self,**coefficients):
2381     """
2382     sets new values to coefficients
2383    
2384     @param coefficients: new values assigned to coefficients
2385     @keyword M: value for coefficient M.
2386     @type M: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
2387     @keyword M_reduced: value for coefficient M_reduced.
2388     @type M_reduced: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.ReducedFunction>}.
2389     @keyword A: value for coefficient A.
2390     @type A: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
2391     @keyword A_reduced: value for coefficient A_reduced.
2392     @type A_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
2393     @keyword B: value for coefficient B
2394     @type B: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
2395     @keyword B_reduced: value for coefficient B_reduced
2396     @type B_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
2397     @keyword C: value for coefficient C
2398     @type C: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
2399     @keyword C_reduced: value for coefficient C_reduced
2400     @type C_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
2401     @keyword D: value for coefficient D
2402     @type D: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
2403