/[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 2061 - (hide annotations)
Wed Nov 19 03:40:21 2008 UTC (10 years, 10 months ago) by jfenwick
File MIME type: text/x-python
File size: 117739 byte(s)
Fixing some warnings from epydoc.
Disabled c++ signatures in python docstrings.
Removed references to Bruce in epydoc and the users guide.
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 jgs 148 def trace(self,text):
482     """
483 jgs 149 print the text message if debugging is swiched on.
484     @param text: message
485     @type text: C{string}
486 jgs 102 """
487 jgs 148 if self.__debug: print "%s: %s"%(str(self),text)
488 jgs 102
489 jgs 148 # =============================================================================
490     # some service functions:
491     # =============================================================================
492 gross 1841 def introduceCoefficients(self,**coeff):
493     """
494     introduces a new coefficient into the problem.
495    
496     use
497    
498 jfenwick 2061 self.introduceCoefficients(A=PDECoef(...),B=PDECoef(...))
499 gross 1841
500     to introduce the coefficients M{A} ans M{B}.
501     """
502     for name, type in coeff.items():
503     if not isinstance(type,PDECoef):
504     raise ValueError,"coefficient %s has no type."%name
505     self.__COEFFICIENTS[name]=type
506     self.__COEFFICIENTS[name].resetValue()
507     self.trace("coefficient %s has been introduced."%name)
508    
509 jgs 148 def getDomain(self):
510 jgs 102 """
511 jgs 148 returns the domain of the PDE
512 jgs 102
513 jgs 149 @return: the domain of the PDE
514     @rtype: L{Domain<escript.Domain>}
515 jgs 108 """
516 jgs 148 return self.__domain
517 jgs 122
518 jgs 148 def getDim(self):
519 jgs 108 """
520 jgs 148 returns the spatial dimension of the PDE
521 jgs 108
522 jgs 149 @return: the spatial dimension of the PDE domain
523     @rtype: C{int}
524 jgs 148 """
525     return self.getDomain().getDim()
526 jgs 102
527 jgs 148 def getNumEquations(self):
528     """
529     returns the number of equations
530 jgs 102
531 jgs 149 @return: the number of equations
532     @rtype: C{int}
533 jgs 148 @raise UndefinedPDEError: if the number of equations is not be specified yet.
534     """
535     if self.__numEquations==None:
536 gross 1859 if self.__numSolutions==None:
537     raise UndefinedPDEError,"Number of equations is undefined. Please specify argument numEquations."
538     else:
539     self.__numEquations=self.__numSolutions
540     return self.__numEquations
541 jgs 147
542 jgs 148 def getNumSolutions(self):
543     """
544     returns the number of unknowns
545 jgs 147
546 jgs 149 @return: the number of unknowns
547     @rtype: C{int}
548 jgs 148 @raise UndefinedPDEError: if the number of unknowns is not be specified yet.
549     """
550     if self.__numSolutions==None:
551 gross 1859 if self.__numEquations==None:
552     raise UndefinedPDEError,"Number of solution is undefined. Please specify argument numSolutions."
553     else:
554     self.__numSolutions=self.__numEquations
555     return self.__numSolutions
556 jgs 148
557 jgs 149 def reduceEquationOrder(self):
558     """
559     return status for order reduction for equation
560    
561     @return: return True is reduced interpolation order is used for the represenation of the equation
562     @rtype: L{bool}
563     """
564     return self.__reduce_equation_order
565    
566     def reduceSolutionOrder(self):
567     """
568     return status for order reduction for the solution
569    
570     @return: return True is reduced interpolation order is used for the represenation of the solution
571     @rtype: L{bool}
572     """
573     return self.__reduce_solution_order
574 jgs 151
575 jgs 108 def getFunctionSpaceForEquation(self):
576     """
577 jgs 149 returns the L{FunctionSpace<escript.FunctionSpace>} used to discretize the equation
578 jgs 148
579 jgs 149 @return: representation space of equation
580     @rtype: L{FunctionSpace<escript.FunctionSpace>}
581 jgs 108 """
582 jgs 149 if self.reduceEquationOrder():
583     return escript.ReducedSolution(self.getDomain())
584     else:
585     return escript.Solution(self.getDomain())
586 jgs 108
587     def getFunctionSpaceForSolution(self):
588     """
589 jgs 149 returns the L{FunctionSpace<escript.FunctionSpace>} used to represent the solution
590 jgs 148
591 jgs 149 @return: representation space of solution
592     @rtype: L{FunctionSpace<escript.FunctionSpace>}
593 jgs 108 """
594 jgs 149 if self.reduceSolutionOrder():
595     return escript.ReducedSolution(self.getDomain())
596     else:
597     return escript.Solution(self.getDomain())
598 jgs 108
599 jgs 148 # =============================================================================
600     # solver settings:
601     # =============================================================================
602 gross 387 def setSolverMethod(self,solver=None,preconditioner=None):
603 jgs 102 """
604 jgs 122 sets a new solver
605 jgs 148
606 jgs 149 @param solver: sets a new solver method.
607 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}
608 gross 387 @param preconditioner: sets a new solver method.
609 artak 1819 @type preconditioner: one of L{DEFAULT}, L{JACOBI} L{ILU0}, L{ILUT},L{SSOR}, L{RILU}, L{GS}
610 gross 1841
611     @note: the solver method choosen may not be available by the selected package.
612 jgs 102 """
613 ksteube 1312 if solver==None: solver=self.__solver_method
614     if preconditioner==None: preconditioner=self.__preconditioner
615     if solver==None: solver=self.DEFAULT
616 gross 387 if preconditioner==None: preconditioner=self.DEFAULT
617     if not (solver,preconditioner)==self.getSolverMethod():
618 jgs 102 self.__solver_method=solver
619 gross 387 self.__preconditioner=preconditioner
620 gross 1841 self.updateSystemType()
621 jgs 148 self.trace("New solver is %s"%self.getSolverMethodName())
622 jgs 102
623 gross 1841 def getSolverMethodName(self):
624 jgs 148 """
625     returns the name of the solver currently used
626    
627 jgs 149 @return: the name of the solver currently used.
628 jgs 148 @rtype: C{string}
629     """
630    
631     m=self.getSolverMethod()
632 jgs 150 p=self.getSolverPackage()
633 gross 425 method=""
634 gross 387 if m[0]==self.DEFAULT: method="DEFAULT"
635     elif m[0]==self.DIRECT: method= "DIRECT"
636     elif m[0]==self.ITERATIVE: method= "ITERATIVE"
637     elif m[0]==self.CHOLEVSKY: method= "CHOLEVSKY"
638     elif m[0]==self.PCG: method= "PCG"
639 artak 1787 elif m[0]==self.TFQMR: method= "TFQMR"
640     elif m[0]==self.MINRES: method= "MINRES"
641 gross 387 elif m[0]==self.CR: method= "CR"
642     elif m[0]==self.CGS: method= "CGS"
643     elif m[0]==self.BICGSTAB: method= "BICGSTAB"
644     elif m[0]==self.SSOR: method= "SSOR"
645     elif m[0]==self.GMRES: method= "GMRES"
646     elif m[0]==self.PRES20: method= "PRES20"
647     elif m[0]==self.LUMPING: method= "LUMPING"
648 gross 431 elif m[0]==self.AMG: method= "AMG"
649 gross 425 if m[1]==self.DEFAULT: method+="+DEFAULT"
650     elif m[1]==self.JACOBI: method+= "+JACOBI"
651     elif m[1]==self.ILU0: method+= "+ILU0"
652     elif m[1]==self.ILUT: method+= "+ILUT"
653     elif m[1]==self.SSOR: method+= "+SSOR"
654 gross 431 elif m[1]==self.AMG: method+= "+AMG"
655     elif m[1]==self.RILU: method+= "+RILU"
656 artak 1819 elif m[1]==self.GS: method+= "+GS"
657 jgs 150 if p==self.DEFAULT: package="DEFAULT"
658     elif p==self.PASO: package= "PASO"
659     elif p==self.MKL: package= "MKL"
660     elif p==self.SCSL: package= "SCSL"
661     elif p==self.UMFPACK: package= "UMFPACK"
662 ksteube 1312 elif p==self.TRILINOS: package= "TRILINOS"
663 jgs 150 else : method="unknown"
664     return "%s solver of %s package"%(method,package)
665 jgs 148
666 gross 1841 def getSolverMethod(self):
667 jgs 102 """
668 gross 1841 returns the solver method and preconditioner used
669 jgs 149
670 jgs 151 @return: the solver method currently be used.
671 gross 1841 @rtype: C{tuple} of C{int}
672 jgs 102 """
673 gross 387 return self.__solver_method,self.__preconditioner
674 jgs 102
675 gross 1841 def setSolverPackage(self,package=None):
676 jgs 150 """
677     sets a new solver package
678    
679 gross 720 @param package: sets a new solver method.
680 ksteube 1312 @type package: one of L{DEFAULT}, L{PASO} L{SCSL}, L{MKL}, L{UMFPACK}, L{TRILINOS}
681 jgs 150 """
682     if package==None: package=self.DEFAULT
683     if not package==self.getSolverPackage():
684 gross 727 self.__solver_package=package
685 gross 1841 self.updateSystemType()
686 jgs 150 self.trace("New solver is %s"%self.getSolverMethodName())
687    
688 gross 1841 def getSolverPackage(self):
689 jgs 150 """
690     returns the package of the solver
691    
692 jgs 151 @return: the solver package currently being used.
693 jgs 150 @rtype: C{int}
694     """
695     return self.__solver_package
696    
697 jgs 148 def isUsingLumping(self):
698     """
699 gross 1841 checks if matrix lumping is used as a solver method
700 jgs 148
701 jgs 149 @return: True is lumping is currently used a solver method.
702 jgs 148 @rtype: C{bool}
703     """
704 gross 387 return self.getSolverMethod()[0]==self.LUMPING
705 jgs 148
706 gross 1841 def setTolerance(self,tol=1.e-8):
707 jgs 102 """
708 gross 1841 resets the tolerance for the solver method to tol.
709 jgs 148
710     @param tol: new tolerance for the solver. If the tol is lower then the current tolerence
711     the system will be resolved.
712 jgs 149 @type tol: positive C{float}
713 gross 877 @raise ValueError: if tolerance is not positive.
714 jgs 102 """
715     if not tol>0:
716 gross 877 raise ValueError,"Tolerance as to be positive"
717 gross 1841 if tol<self.getTolerance(): self.invalidateSolution()
718 jgs 148 self.trace("New tolerance %e"%tol)
719 jgs 102 self.__tolerance=tol
720     return
721 jgs 148
722 gross 1841 def getTolerance(self):
723 jgs 102 """
724 jgs 122 returns the tolerance set for the solution
725 jgs 148
726     @return: tolerance currently used.
727     @rtype: C{float}
728 jgs 102 """
729     return self.__tolerance
730    
731 jgs 148 # =============================================================================
732     # symmetry flag:
733     # =============================================================================
734 gross 1841 def isSymmetric(self):
735 jgs 102 """
736 jgs 148 checks if symmetry is indicated.
737 jgs 149
738     @return: True is a symmetric PDE is indicated, otherwise False is returned
739     @rtype: C{bool}
740 jgs 102 """
741     return self.__sym
742    
743 gross 1841 def setSymmetryOn(self):
744 jgs 102 """
745 jgs 148 sets the symmetry flag.
746 jgs 102 """
747     if not self.isSymmetric():
748 jgs 148 self.trace("PDE is set to be symmetric")
749 jgs 102 self.__sym=True
750 gross 1841 self.updateSystemType()
751 jgs 102
752 gross 1841 def setSymmetryOff(self):
753 jgs 102 """
754 jgs 148 removes the symmetry flag.
755 jgs 102 """
756     if self.isSymmetric():
757 jgs 148 self.trace("PDE is set to be unsymmetric")
758 jgs 102 self.__sym=False
759 gross 1841 self.updateSystemType()
760 jgs 102
761 gross 1841 def setSymmetryTo(self,flag=False):
762 jgs 148 """
763     sets the symmetry flag to flag
764 jgs 149
765 jgs 148 @param flag: If flag, the symmetry flag is set otherwise the symmetry flag is released.
766     @type flag: C{bool}
767     """
768     if flag:
769     self.setSymmetryOn()
770     else:
771     self.setSymmetryOff()
772 jgs 102
773 jgs 148 # =============================================================================
774     # function space handling for the equation as well as the solution
775     # =============================================================================
776 gross 1841 def setReducedOrderOn(self):
777 jgs 102 """
778 jgs 148 switches on reduced order for solution and equation representation
779 jgs 149
780     @raise RuntimeError: if order reduction is altered after a coefficient has been set.
781 jgs 102 """
782     self.setReducedOrderForSolutionOn()
783     self.setReducedOrderForEquationOn()
784    
785 gross 1841 def setReducedOrderOff(self):
786 jgs 102 """
787 jgs 148 switches off reduced order for solution and equation representation
788 jgs 149
789     @raise RuntimeError: if order reduction is altered after a coefficient has been set.
790 jgs 102 """
791     self.setReducedOrderForSolutionOff()
792     self.setReducedOrderForEquationOff()
793    
794 gross 1841 def setReducedOrderTo(self,flag=False):
795 jgs 102 """
796 jgs 149 sets order reduction for both solution and equation representation according to flag.
797     @param flag: if flag is True, the order reduction is switched on for both solution and equation representation, otherwise or
798 jgs 148 if flag is not present order reduction is switched off
799     @type flag: C{bool}
800 jgs 149 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
801 jgs 102 """
802     self.setReducedOrderForSolutionTo(flag)
803     self.setReducedOrderForEquationTo(flag)
804    
805 jgs 148
806 gross 1841 def setReducedOrderForSolutionOn(self):
807 jgs 102 """
808 jgs 148 switches on reduced order for solution representation
809 jgs 149
810     @raise RuntimeError: if order reduction is altered after a coefficient has been set.
811 jgs 102 """
812 jgs 149 if not self.__reduce_solution_order:
813     if self.__altered_coefficients:
814     raise RuntimeError,"order cannot be altered after coefficients have been defined."
815 jgs 148 self.trace("Reduced order is used to solution representation.")
816 jgs 149 self.__reduce_solution_order=True
817 gross 1841 self.initializeSystem()
818 jgs 102
819 gross 1841 def setReducedOrderForSolutionOff(self):
820 jgs 102 """
821 jgs 148 switches off reduced order for solution representation
822 jgs 149
823     @raise RuntimeError: if order reduction is altered after a coefficient has been set.
824 jgs 102 """
825 jgs 149 if self.__reduce_solution_order:
826     if self.__altered_coefficients:
827     raise RuntimeError,"order cannot be altered after coefficients have been defined."
828 jgs 148 self.trace("Full order is used to interpolate solution.")
829 jgs 149 self.__reduce_solution_order=False
830 gross 1841 self.initializeSystem()
831 jgs 102
832 gross 1841 def setReducedOrderForSolutionTo(self,flag=False):
833 jgs 102 """
834 jgs 122 sets order for test functions according to flag
835 jgs 102
836 jgs 149 @param flag: if flag is True, the order reduction is switched on for solution representation, otherwise or
837 jgs 148 if flag is not present order reduction is switched off
838     @type flag: C{bool}
839 jgs 149 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
840 jgs 102 """
841     if flag:
842     self.setReducedOrderForSolutionOn()
843     else:
844     self.setReducedOrderForSolutionOff()
845 jgs 148
846 gross 1841 def setReducedOrderForEquationOn(self):
847 jgs 102 """
848 jgs 148 switches on reduced order for equation representation
849 jgs 149
850     @raise RuntimeError: if order reduction is altered after a coefficient has been set.
851 jgs 102 """
852 jgs 149 if not self.__reduce_equation_order:
853     if self.__altered_coefficients:
854     raise RuntimeError,"order cannot be altered after coefficients have been defined."
855 jgs 148 self.trace("Reduced order is used for test functions.")
856 jgs 149 self.__reduce_equation_order=True
857 gross 1841 self.initializeSystem()
858 jgs 102
859 gross 1841 def setReducedOrderForEquationOff(self):
860 jgs 102 """
861 jgs 148 switches off reduced order for equation representation
862 jgs 149
863     @raise RuntimeError: if order reduction is altered after a coefficient has been set.
864 jgs 102 """
865 jgs 149 if self.__reduce_equation_order:
866     if self.__altered_coefficients:
867     raise RuntimeError,"order cannot be altered after coefficients have been defined."
868 jgs 148 self.trace("Full order is used for test functions.")
869 jgs 149 self.__reduce_equation_order=False
870 gross 1841 self.initializeSystem()
871 jgs 102
872 gross 1841 def setReducedOrderForEquationTo(self,flag=False):
873 jgs 102 """
874 jgs 122 sets order for test functions according to flag
875 jgs 102
876 jgs 149 @param flag: if flag is True, the order reduction is switched on for equation representation, otherwise or
877 jgs 148 if flag is not present order reduction is switched off
878     @type flag: C{bool}
879 jgs 149 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
880 jgs 102 """
881     if flag:
882     self.setReducedOrderForEquationOn()
883     else:
884     self.setReducedOrderForEquationOff()
885 jgs 148
886 gross 1841 def updateSystemType(self):
887 jgs 148 """
888 gross 1841 reasses the system type and, if a new matrix is needed, resets the system.
889 jgs 148 """
890 gross 1841 new_system_type=self.getRequiredSystemType()
891     if not new_system_type==self.__system_type:
892     self.trace("Matrix type is now %d."%new_system_type)
893     self.__system_type=new_system_type
894     self.initializeSystem()
895     def getSystemType(self):
896     """
897     return the current system type
898     """
899     return self.__system_type
900 jgs 148
901 gross 1841 def checkSymmetricTensor(self,name,verbose=True):
902     """
903     tests a coefficient for symmetry
904 jgs 148
905 gross 1841 @param name: name of the coefficient
906     @type name: C{str}
907     @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.
908     @type verbose: C{bool}
909     @return: True if coefficient C{name} is symmetric.
910     @rtype: C{bool}
911     """
912     A=self.getCoefficient(name)
913     verbose=verbose or self.__debug
914     out=True
915     if not A.isEmpty():
916     tol=util.Lsup(A)*self.SMALL_TOLERANCE
917     s=A.getShape()
918     if A.getRank() == 4:
919     if s[0]==s[2] and s[1] == s[3]:
920     for i in range(s[0]):
921     for j in range(s[1]):
922     for k in range(s[2]):
923     for l in range(s[3]):
924     if util.Lsup(A[i,j,k,l]-A[k,l,i,j])>tol:
925     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)
926     out=False
927     else:
928     if verbose: print "non-symmetric problem because of inappropriate shape %s of coefficient %s."%(s,name)
929     out=False
930     elif A.getRank() == 2:
931     if s[0]==s[1]:
932     for j in range(s[0]):
933     for l in range(s[1]):
934     if util.Lsup(A[j,l]-A[l,j])>tol:
935     if verbose: print "non-symmetric problem because %s[%d,%d]!=%s[%d,%d]"%(name,j,l,name,l,j)
936     out=False
937     else:
938     if verbose: print "non-symmetric problem because of inappropriate shape %s of coefficient %s."%(s,name)
939     out=False
940     elif A.getRank() == 0:
941     pass
942     else:
943     raise ValueError,"Cannot check rank %s of %s."%(A.getRank(),name)
944     return out
945 jgs 148
946 gross 1841 def checkReciprocalSymmetry(self,name0,name1,verbose=True):
947     """
948     test two coefficients for resciprocal symmetry
949 jgs 148
950 gross 1841 @param name0: name of the first coefficient
951 jfenwick 2061 @type name0: C{str}
952 gross 1841 @param name1: name of the second coefficient
953 jfenwick 2061 @type name1: C{str}
954 gross 1841 @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.
955     @type verbose: C{bool}
956     @return: True if coefficients C{name0} and C{name1} are resciprocally symmetric.
957     @rtype: C{bool}
958     """
959     B=self.getCoefficient(name0)
960     C=self.getCoefficient(name1)
961     verbose=verbose or self.__debug
962     out=True
963     if B.isEmpty() and not C.isEmpty():
964     if verbose: print "non-symmetric problem because %s is not present but %s is"%(name0,name1)
965     out=False
966     elif not B.isEmpty() and C.isEmpty():
967     if verbose: print "non-symmetric problem because %s is not present but %s is"%(name0,name1)
968     out=False
969     elif not B.isEmpty() and not C.isEmpty():
970     sB=B.getShape()
971     sC=C.getShape()
972     tol=(util.Lsup(B)+util.Lsup(C))*self.SMALL_TOLERANCE/2.
973     if len(sB) != len(sC):
974     if verbose: print "non-symmetric problem because ranks of %s (=%s) and %s (=%s) are different."%(name0,len(sB),name1,len(sC))
975     out=False
976     else:
977     if len(sB)==0:
978     if util.Lsup(B-C)>tol:
979     if verbose: print "non-symmetric problem because %s!=%s"%(name0,name1)
980     out=False
981     elif len(sB)==1:
982     if sB[0]==sC[0]:
983     for j in range(sB[0]):
984     if util.Lsup(B[j]-C[j])>tol:
985     if verbose: print "non-symmetric PDE because %s[%d]!=%s[%d]"%(name0,j,name1,j)
986     out=False
987     else:
988     if verbose: print "non-symmetric problem because of inappropriate shapes %s and %s of coefficients %s and %s, respectively."%(sB,sC,name0,name1)
989     elif len(sB)==3:
990     if sB[0]==sC[1] and sB[1]==sC[2] and sB[2]==sC[0]:
991     for i in range(sB[0]):
992     for j in range(sB[1]):
993     for k in range(sB[2]):
994     if util.Lsup(B[i,j,k]-C[k,i,j])>tol:
995     if verbose: print "non-symmetric problem because %s[%d,%d,%d]!=%s[%d,%d,%d]"%(name2,i,j,k,name1,k,i,j)
996     out=False
997     else:
998     if verbose: print "non-symmetric problem because of inappropriate shapes %s and %s of coefficients %s and %s, respectively."%(sB,sC,name0,name1)
999 jgs 148 else:
1000 gross 1841 raise ValueError,"Cannot check rank %s of %s and %s."%(len(sB),name0,name1)
1001     return out
1002 jgs 148
1003 gross 1841 def getCoefficient(self,name):
1004 jgs 108 """
1005 jgs 148 returns the value of the coefficient name
1006 jgs 108
1007 jgs 149 @param name: name of the coefficient requested.
1008 jgs 148 @type name: C{string}
1009 jgs 149 @return: the value of the coefficient name
1010     @rtype: L{Data<escript.Data>}
1011 jgs 148 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1012     """
1013     if self.hasCoefficient(name):
1014 gross 1841 return self.__COEFFICIENTS[name].getValue()
1015 jgs 148 else:
1016     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1017 jgs 108
1018 gross 1841 def hasCoefficient(self,name):
1019 jgs 148 """
1020     return True if name is the name of a coefficient
1021 jgs 108
1022 jgs 148 @param name: name of the coefficient enquired.
1023     @type name: C{string}
1024 jgs 149 @return: True if name is the name of a coefficient of the general PDE. Otherwise False.
1025     @rtype: C{bool}
1026 jgs 148 """
1027 gross 1841 return self.__COEFFICIENTS.has_key(name)
1028 jgs 108
1029 gross 1841 def createCoefficient(self, name):
1030 jgs 148 """
1031 jgs 149 create a L{Data<escript.Data>} object corresponding to coefficient name
1032 jgs 108
1033 jgs 149 @return: a coefficient name initialized to 0.
1034     @rtype: L{Data<escript.Data>}
1035 jgs 148 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1036     """
1037     if self.hasCoefficient(name):
1038     return escript.Data(0.,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))
1039     else:
1040     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1041    
1042 gross 1841 def getFunctionSpaceForCoefficient(self,name):
1043 jgs 148 """
1044 jgs 149 return the L{FunctionSpace<escript.FunctionSpace>} to be used for coefficient name
1045 jgs 148
1046     @param name: name of the coefficient enquired.
1047     @type name: C{string}
1048 jgs 149 @return: the function space to be used for coefficient name
1049     @rtype: L{FunctionSpace<escript.FunctionSpace>}
1050 jgs 148 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1051     """
1052     if self.hasCoefficient(name):
1053 gross 1841 return self.__COEFFICIENTS[name].getFunctionSpace(self.getDomain())
1054 jgs 148 else:
1055     raise ValueError,"unknown coefficient %s requested"%name
1056 gross 1841
1057     def getShapeOfCoefficient(self,name):
1058 jgs 148 """
1059 jgs 149 return the shape of the coefficient name
1060 jgs 148
1061     @param name: name of the coefficient enquired.
1062     @type name: C{string}
1063 jgs 149 @return: the shape of the coefficient name
1064     @rtype: C{tuple} of C{int}
1065 jgs 148 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1066     """
1067     if self.hasCoefficient(name):
1068 gross 1841 return self.__COEFFICIENTS[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())
1069 jgs 148 else:
1070     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1071    
1072 gross 1841 def resetAllCoefficients(self):
1073 jgs 148 """
1074     resets all coefficients to there default values.
1075     """
1076 gross 1841 for i in self.__COEFFICIENTS.iterkeys():
1077     self.__COEFFICIENTS[i].resetValue()
1078 jgs 148
1079 gross 1841 def alteredCoefficient(self,name):
1080 jgs 148 """
1081     announce that coefficient name has been changed
1082    
1083     @param name: name of the coefficient enquired.
1084     @type name: C{string}
1085     @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1086 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.
1087 jgs 148 """
1088     if self.hasCoefficient(name):
1089     self.trace("Coefficient %s has been altered."%name)
1090 jgs 149 if not ((name=="q" or name=="r") and self.isUsingLumping()):
1091 gross 1841 if self.__COEFFICIENTS[name].isAlteringOperator(): self.invalidateOperator()
1092     if self.__COEFFICIENTS[name].isAlteringRightHandSide(): self.invalidateRightHandSide()
1093 jgs 148 else:
1094     raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1095    
1096 gross 1841 def validSolution(self):
1097     """
1098     markes the solution as valid
1099     """
1100     self.__is_solution_valid=True
1101 jgs 149
1102 gross 1841 def invalidateSolution(self):
1103     """
1104     indicates the PDE has to be resolved if the solution is requested
1105     """
1106     self.trace("System will be resolved.")
1107     self.__is_solution_valid=False
1108 jgs 148
1109 gross 1841 def isSolutionValid(self):
1110     """
1111     returns True is the solution is still valid.
1112     """
1113     return self.__is_solution_valid
1114    
1115     def validOperator(self):
1116     """
1117     markes the operator as valid
1118     """
1119     self.__is_operator_valid=True
1120    
1121     def invalidateOperator(self):
1122     """
1123     indicates the operator has to be rebuilt next time it is used
1124     """
1125     self.trace("Operator will be rebuilt.")
1126     self.invalidateSolution()
1127     self.__is_operator_valid=False
1128    
1129     def isOperatorValid(self):
1130     """
1131     returns True is the operator is still valid.
1132     """
1133     return self.__is_operator_valid
1134    
1135     def validRightHandSide(self):
1136     """
1137     markes the right hand side as valid
1138     """
1139     self.__is_RHS_valid=True
1140    
1141     def invalidateRightHandSide(self):
1142     """
1143     indicates the right hand side has to be rebuild next time it is used
1144     """
1145     if self.isRightHandSideValid(): self.trace("Right hand side has to be rebuilt.")
1146     self.invalidateSolution()
1147     self.__is_RHS_valid=False
1148    
1149     def isRightHandSideValid(self):
1150     """
1151     returns True is the operator is still valid.
1152     """
1153     return self.__is_RHS_valid
1154    
1155     def invalidateSystem(self):
1156     """
1157     annonced that everthing has to be rebuild:
1158     """
1159     if self.isRightHandSideValid(): self.trace("System has to be rebuilt.")
1160     self.invalidateSolution()
1161     self.invalidateOperator()
1162     self.invalidateRightHandSide()
1163    
1164     def isSystemValid(self):
1165     """
1166     returns true if the system (including solution) is still vaild:
1167     """
1168     return self.isSolutionValid() and self.isOperatorValid() and self.isRightHandSideValid()
1169    
1170     def initializeSystem(self):
1171     """
1172     resets the system
1173     """
1174     self.trace("New System has been created.")
1175     self.__operator=escript.Operator()
1176     self.__righthandside=escript.Data()
1177     self.__solution=escript.Data()
1178     self.invalidateSystem()
1179    
1180     def getOperator(self):
1181     """
1182     returns the operator of the linear problem
1183    
1184     @return: the operator of the problem
1185     """
1186     return self.getSystem()[0]
1187    
1188     def getRightHandSide(self):
1189     """
1190     returns the right hand side of the linear problem
1191    
1192     @return: the right hand side of the problem
1193     @rtype: L{Data<escript.Data>}
1194     """
1195     return self.getSystem()[1]
1196    
1197     def createRightHandSide(self):
1198     """
1199     returns an instance of a new right hand side
1200     """
1201     self.trace("New right hand side is allocated.")
1202     if self.getNumEquations()>1:
1203     return escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)
1204     else:
1205     return escript.Data(0.,(),self.getFunctionSpaceForEquation(),True)
1206    
1207     def createSolution(self):
1208     """
1209     returns an instance of a new solution
1210     """
1211     self.trace("New solution is allocated.")
1212     if self.getNumSolutions()>1:
1213     return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)
1214     else:
1215     return escript.Data(0.,(),self.getFunctionSpaceForSolution(),True)
1216    
1217     def resetSolution(self):
1218     """
1219     sets the solution to zero
1220     """
1221     if self.__solution.isEmpty():
1222     self.__solution=self.createSolution()
1223     else:
1224     self.__solution.setToZero()
1225     self.trace("Solution is reset to zero.")
1226     def setSolution(self,u):
1227     """
1228     sets the solution assuming that makes the sytem valid.
1229     """
1230     self.__solution=u
1231     self.validSolution()
1232    
1233     def getCurrentSolution(self):
1234     """
1235     returns the solution in its current state
1236     """
1237     if self.__solution.isEmpty(): self.__solution=self.createSolution()
1238     return self.__solution
1239    
1240     def resetRightHandSide(self):
1241     """
1242     sets the right hand side to zero
1243     """
1244     if self.__righthandside.isEmpty():
1245     self.__righthandside=self.createRightHandSide()
1246     else:
1247     self.__righthandside.setToZero()
1248     self.trace("Right hand side is reset to zero.")
1249    
1250     def getCurrentRightHandSide(self):
1251     """
1252     returns the right hand side in its current state
1253     """
1254     if self.__righthandside.isEmpty(): self.__righthandside=self.createRightHandSide()
1255     return self.__righthandside
1256    
1257    
1258     def resetOperator(self):
1259     """
1260     makes sure that the operator is instantiated and returns it initialized by zeros
1261     """
1262     if self.__operator.isEmpty():
1263     if self.isUsingLumping():
1264     self.__operator=self.createSolution()
1265     else:
1266     self.__operator=self.createOperator()
1267     else:
1268     if self.isUsingLumping():
1269     self.__operator.setToZero()
1270     else:
1271     self.__operator.resetValues()
1272     self.trace("Operator reset to zero")
1273    
1274     def getCurrentOperator(self):
1275     """
1276     returns the operator in its current state
1277     """
1278     if self.__operator.isEmpty(): self.resetOperator()
1279     return self.__operator
1280    
1281     def setValue(self,**coefficients):
1282 jgs 148 """
1283     sets new values to coefficients
1284    
1285     @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1286     """
1287 jgs 108 # check if the coefficients are legal:
1288     for i in coefficients.iterkeys():
1289     if not self.hasCoefficient(i):
1290 jgs 148 raise IllegalCoefficient,"Attempt to set unknown coefficient %s"%i
1291 jgs 108 # if the number of unknowns or equations is still unknown we try to estimate them:
1292 jgs 148 if self.__numEquations==None or self.__numSolutions==None:
1293 jgs 108 for i,d in coefficients.iteritems():
1294     if hasattr(d,"shape"):
1295     s=d.shape
1296     elif hasattr(d,"getShape"):
1297     s=d.getShape()
1298     else:
1299     s=numarray.array(d).shape
1300     if s!=None:
1301     # get number of equations and number of unknowns:
1302 gross 1841 res=self.__COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(self.getDomain(),s)
1303 jgs 108 if res==None:
1304 jgs 148 raise IllegalCoefficientValue,"Illegal shape %s of coefficient %s"%(s,i)
1305 jgs 108 else:
1306 jgs 148 if self.__numEquations==None: self.__numEquations=res[0]
1307     if self.__numSolutions==None: self.__numSolutions=res[1]
1308     if self.__numEquations==None: raise UndefinedPDEError,"unidententified number of equations"
1309     if self.__numSolutions==None: raise UndefinedPDEError,"unidententified number of solutions"
1310 jgs 108 # now we check the shape of the coefficient if numEquations and numSolutions are set:
1311     for i,d in coefficients.iteritems():
1312 jgs 148 try:
1313 gross 1841 self.__COEFFICIENTS[i].setValue(self.getDomain(),
1314 gross 1072 self.getNumEquations(),self.getNumSolutions(),
1315     self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
1316     self.alteredCoefficient(i)
1317     except IllegalCoefficientFunctionSpace,m:
1318     # if the function space is wrong then we try the reduced version:
1319     i_red=i+"_reduced"
1320 gross 1841 if (not i_red in coefficients.keys()) and i_red in self.__COEFFICIENTS.keys():
1321 gross 1072 try:
1322 gross 1841 self.__COEFFICIENTS[i_red].setValue(self.getDomain(),
1323 gross 1072 self.getNumEquations(),self.getNumSolutions(),
1324     self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
1325     self.alteredCoefficient(i_red)
1326     except IllegalCoefficientValue,m:
1327     raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
1328     except IllegalCoefficientFunctionSpace,m:
1329     raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m))
1330     else:
1331     raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m))
1332 jgs 148 except IllegalCoefficientValue,m:
1333     raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
1334 jgs 149 self.__altered_coefficients=True
1335 jgs 108
1336 gross 1841 # ====================================================================================
1337     # methods that are typically overwritten when implementing a particular linear problem
1338     # ====================================================================================
1339     def getRequiredSystemType(self):
1340     """
1341     returns the system type which needs to be used by the current set up.
1342    
1343     @note: Typically this method is overwritten when implementing a particular linear problem.
1344     """
1345     return 0
1346    
1347     def createOperator(self):
1348 jgs 108 """
1349 gross 1841 returns an instance of a new operator
1350    
1351     @note: This method is overwritten when implementing a particular linear problem.
1352     """
1353     return escript.Operator()
1354    
1355     def checkSymmetry(self,verbose=True):
1356     """
1357     test the PDE for symmetry.
1358    
1359     @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.
1360     @type verbose: C{bool}
1361     @return: True if the problem is symmetric.
1362     @rtype: C{bool}
1363     @note: Typically this method is overwritten when implementing a particular linear problem.
1364     """
1365     out=True
1366     return out
1367    
1368     def getSolution(self,**options):
1369     """
1370     returns the solution of the problem
1371     @return: the solution
1372     @rtype: L{Data<escript.Data>}
1373    
1374     @note: This method is overwritten when implementing a particular linear problem.
1375     """
1376     return self.getCurrentSolution()
1377    
1378     def getSystem(self):
1379     """
1380 jgs 122 return the operator and right hand side of the PDE
1381 jgs 149
1382     @return: the discrete version of the PDE
1383     @rtype: C{tuple} of L{Operator,<escript.Operator>} and L{Data<escript.Data>}.
1384 gross 1841
1385     @note: This method is overwritten when implementing a particular linear problem.
1386 jgs 108 """
1387 gross 1841 return (self.getCurrentOperator(), self.getCurrentRightHandSide())
1388     #=====================
1389    
1390     class LinearPDE(LinearProblem):
1391     """
1392     This class is used to define a general linear, steady, second order PDE
1393     for an unknown function M{u} on a given domain defined through a L{Domain<escript.Domain>} object.
1394    
1395     For a single PDE with a solution with a single component the linear PDE is defined in the following form:
1396    
1397     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)}
1398    
1399    
1400     where M{grad(F)} denotes the spatial derivative of M{F}. Einstein's summation convention,
1401     ie. summation over indexes appearing twice in a term of a sum is performed, is used.
1402     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
1403     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
1404     L{ReducedFunction<escript.ReducedFunction>}. It is also allowd to use objects that can be converted into
1405     such L{Data<escript.Data>} objects. M{A} and M{A_reduced} are rank two, M{B}, M{C}, M{X},
1406     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.
1407    
1408     The following natural boundary conditions are considered:
1409    
1410     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}
1411    
1412     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>}.
1413    
1414    
1415     Constraints for the solution prescribing the value of the solution at certain locations in the domain. They have the form
1416    
1417     M{u=r} where M{q>0}
1418    
1419     M{r} and M{q} are each scalar where M{q} is the characteristic function defining where the constraint is applied.
1420     The constraints override any other condition set by the PDE or the boundary condition.
1421    
1422     The PDE is symmetrical if
1423    
1424     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]}
1425    
1426     For a system of PDEs and a solution with several components the PDE has the form
1427    
1428     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] }
1429    
1430     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.
1431     The natural boundary conditions take the form:
1432    
1433     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]}
1434    
1435    
1436     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>}.
1437    
1438     Constraints take the form
1439    
1440     M{u[i]=r[i]} where M{q[i]>0}
1441    
1442     M{r} and M{q} are each rank one. Notice that at some locations not necessarily all components must have a constraint.
1443    
1444     The system of PDEs is symmetrical if
1445    
1446     - M{A[i,j,k,l]=A[k,l,i,j]}
1447     - M{A_reduced[i,j,k,l]=A_reduced[k,l,i,j]}
1448     - M{B[i,j,k]=C[k,i,j]}
1449     - M{B_reduced[i,j,k]=C_reduced[k,i,j]}
1450     - M{D[i,k]=D[i,k]}
1451     - M{D_reduced[i,k]=D_reduced[i,k]}
1452     - M{d[i,k]=d[k,i]}
1453     - M{d_reduced[i,k]=d_reduced[k,i]}
1454    
1455     L{LinearPDE} also supports solution discontinuities over a contact region in the domain. To specify the conditions across the
1456     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
1457     defined as
1458    
1459     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]}
1460    
1461     For the case of single solution component and single PDE M{J} is defined
1462    
1463     M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[j]+(B[i]+B_reduced[i])*u-X[i]-X_reduced[i]}
1464    
1465     In the context of discontinuities M{n} denotes the normal on the discontinuity pointing from side 0 towards side 1
1466     calculated from L{getNormal<escript.FunctionSpace.getNormal>} of L{FunctionOnContactZero<escript.FunctionOnContactZero>}. For a system of PDEs
1467     the contact condition takes the form
1468    
1469     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]}
1470    
1471     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
1472     of the solution at side 1 and at side 0, denotes the jump of M{u} across discontinuity along the normal calculated by
1473     L{jump<util.jump>}.
1474     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>}.
1475     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>}.
1476     In case of a single PDE and a single component solution the contact condition takes the form
1477    
1478     M{n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)}
1479    
1480     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>}
1481    
1482     typical usage:
1483    
1484     p=LinearPDE(dom)
1485     p.setValue(A=kronecker(dom),D=1,Y=0.5)
1486     u=p.getSolution()
1487    
1488     """
1489    
1490    
1491     def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):
1492     """
1493     initializes a new linear PDE
1494    
1495     @param domain: domain of the PDE
1496     @type domain: L{Domain<escript.Domain>}
1497     @param numEquations: number of equations. If numEquations==None the number of equations
1498     is extracted from the PDE coefficients.
1499     @param numSolutions: number of solution components. If numSolutions==None the number of solution components
1500     is extracted from the PDE coefficients.
1501     @param debug: if True debug informations are printed.
1502    
1503     """
1504     super(LinearPDE, self).__init__(domain,numEquations,numSolutions,debug)
1505     #
1506     # the coefficients of the PDE:
1507     #
1508     self.introduceCoefficients(
1509     A=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
1510     B=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1511     C=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
1512     D=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1513     X=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
1514     Y=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1515     d=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1516     y=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1517     d_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1518     y_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1519     A_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
1520     B_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1521     C_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
1522     D_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1523     X_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
1524     Y_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1525     d_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1526     y_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1527     d_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1528     y_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1529     r=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.RIGHTHANDSIDE),
1530     q=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.BOTH) )
1531    
1532     def __str__(self):
1533     """
1534     returns string representation of the PDE
1535    
1536     @return: a simple representation of the PDE
1537     @rtype: C{str}
1538     """
1539     return "<LinearPDE %d>"%id(self)
1540    
1541     def getRequiredSystemType(self):
1542     """
1543     returns the system type which needs to be used by the current set up.
1544     """
1545 gross 1859 return self.getDomain().getSystemMatrixTypeId(self.getSolverMethod()[0],self.getSolverMethod()[1],self.getSolverPackage(),self.isSymmetric())
1546 gross 1841
1547     def checkSymmetry(self,verbose=True):
1548     """
1549     test the PDE for symmetry.
1550    
1551     @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.
1552     @type verbose: C{bool}
1553     @return: True if the PDE is symmetric.
1554     @rtype: L{bool}
1555     @note: This is a very expensive operation. It should be used for degugging only! The symmetry flag is not altered.
1556     """
1557     out=True
1558     out=out and self.checkSymmetricTensor("A", verbose)
1559     out=out and self.checkSymmetricTensor("A_reduced", verbose)
1560     out=out and self.checkReciprocalSymmetry("B","C", verbose)
1561     out=out and self.checkReciprocalSymmetry("B_reduced","C_reduced", verbose)
1562     out=out and self.checkSymmetricTensor("D", verbose)
1563     out=out and self.checkSymmetricTensor("D_reduced", verbose)
1564     out=out and self.checkSymmetricTensor("d", verbose)
1565     out=out and self.checkSymmetricTensor("d_reduced", verbose)
1566     out=out and self.checkSymmetricTensor("d_contact", verbose)
1567     out=out and self.checkSymmetricTensor("d_contact_reduced", verbose)
1568     return out
1569    
1570     def createOperator(self):
1571     """
1572     returns an instance of a new operator
1573     """
1574     self.trace("New operator is allocated.")
1575     return self.getDomain().newOperator( \
1576     self.getNumEquations(), \
1577     self.getFunctionSpaceForEquation(), \
1578     self.getNumSolutions(), \
1579     self.getFunctionSpaceForSolution(), \
1580     self.getSystemType())
1581    
1582     def getSolution(self,**options):
1583     """
1584     returns the solution of the PDE
1585    
1586     @return: the solution
1587     @rtype: L{Data<escript.Data>}
1588     @param options: solver options
1589     @keyword verbose: True to get some information during PDE solution
1590     @type verbose: C{bool}
1591     @keyword reordering: reordering scheme to be used during elimination. Allowed values are
1592     L{NO_REORDERING}, L{MINIMUM_FILL_IN}, L{NESTED_DISSECTION}
1593     @keyword iter_max: maximum number of iteration steps allowed.
1594     @keyword drop_tolerance: threshold for drupping in L{ILUT}
1595     @keyword drop_storage: maximum of allowed memory in L{ILUT}
1596     @keyword truncation: maximum number of residuals in L{GMRES}
1597     @keyword restart: restart cycle length in L{GMRES}
1598     """
1599     if not self.isSolutionValid():
1600     mat,f=self.getSystem()
1601 jgs 108 if self.isUsingLumping():
1602 gross 1841 self.setSolution(f*1/mat)
1603     else:
1604     options[self.TOLERANCE_KEY]=self.getTolerance()
1605     options[self.METHOD_KEY]=self.getSolverMethod()[0]
1606     options[self.PRECONDITIONER_KEY]=self.getSolverMethod()[1]
1607     options[self.PACKAGE_KEY]=self.getSolverPackage()
1608     options[self.SYMMETRY_KEY]=self.isSymmetric()
1609     self.trace("PDE is resolved.")
1610     self.trace("solver options: %s"%str(options))
1611     self.setSolution(mat.solve(f,options))
1612     return self.getCurrentSolution()
1613    
1614     def getSystem(self):
1615     """
1616     return the operator and right hand side of the PDE
1617    
1618     @return: the discrete version of the PDE
1619     @rtype: C{tuple} of L{Operator,<escript.Operator>} and L{Data<escript.Data>}.
1620     """
1621     if not self.isOperatorValid() or not self.isRightHandSideValid():
1622     if self.isUsingLumping():
1623     if not self.isOperatorValid():
1624 gross 531 if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():
1625     raise TypeError,"Lumped matrix requires same order for equations and unknowns"
1626 gross 1841 if not self.getCoefficient("A").isEmpty():
1627 gross 531 raise ValueError,"coefficient A in lumped matrix may not be present."
1628 gross 1841 if not self.getCoefficient("B").isEmpty():
1629 gross 1072 raise ValueError,"coefficient B in lumped matrix may not be present."
1630 gross 1841 if not self.getCoefficient("C").isEmpty():
1631 gross 1072 raise ValueError,"coefficient C in lumped matrix may not be present."
1632 gross 1841 if not self.getCoefficient("d_contact").isEmpty():
1633 gross 1204 raise ValueError,"coefficient d_contact in lumped matrix may not be present."
1634 gross 1841 if not self.getCoefficient("A_reduced").isEmpty():
1635 gross 1072 raise ValueError,"coefficient A_reduced in lumped matrix may not be present."
1636 gross 1841 if not self.getCoefficient("B_reduced").isEmpty():
1637 gross 1072 raise ValueError,"coefficient B_reduced in lumped matrix may not be present."
1638 gross 1841 if not self.getCoefficient("C_reduced").isEmpty():
1639 gross 1072 raise ValueError,"coefficient C_reduced in lumped matrix may not be present."
1640 gross 1841 if not self.getCoefficient("d_contact_reduced").isEmpty():
1641 gross 1204 raise ValueError,"coefficient d_contact_reduced in lumped matrix may not be present."
1642 gross 1841 D=self.getCoefficient("D")
1643     d=self.getCoefficient("d")
1644     D_reduced=self.getCoefficient("D_reduced")
1645     d_reduced=self.getCoefficient("d_reduced")
1646 gross 531 if not D.isEmpty():
1647     if self.getNumSolutions()>1:
1648 bcumming 791 D_times_e=util.matrix_mult(D,numarray.ones((self.getNumSolutions(),)))
1649 gross 531 else:
1650     D_times_e=D
1651     else:
1652     D_times_e=escript.Data()
1653     if not d.isEmpty():
1654     if self.getNumSolutions()>1:
1655 bcumming 791 d_times_e=util.matrix_mult(d,numarray.ones((self.getNumSolutions(),)))
1656 gross 531 else:
1657     d_times_e=d
1658     else:
1659     d_times_e=escript.Data()
1660 gross 1204
1661 gross 1072 if not D_reduced.isEmpty():
1662     if self.getNumSolutions()>1:
1663     D_reduced_times_e=util.matrix_mult(D_reduced,numarray.ones((self.getNumSolutions(),)))
1664     else:
1665     D_reduced_times_e=D_reduced
1666     else:
1667     D_reduced_times_e=escript.Data()
1668     if not d_reduced.isEmpty():
1669     if self.getNumSolutions()>1:
1670     d_reduced_times_e=util.matrix_mult(d_reduced,numarray.ones((self.getNumSolutions(),)))
1671     else:
1672     d_reduced_times_e=d_reduced
1673     else:
1674     d_reduced_times_e=escript.Data()
1675 gross 1204
1676 gross 1841 self.resetOperator()
1677     operator=self.getCurrentOperator()
1678 gross 1659 if False and hasattr(self.getDomain(), "addPDEToLumpedSystem") :
1679 gross 1841 self.getDomain().addPDEToLumpedSystem(operator, D_times_e, d_times_e)
1680     self.getDomain().addPDEToLumpedSystem(operator, D_reduced_times_e, d_reduced_times_e)
1681 gross 1072 else:
1682 gross 1841 self.getDomain().addPDEToRHS(operator, \
1683 gross 1204 escript.Data(), \
1684     D_times_e, \
1685     d_times_e,\
1686     escript.Data())
1687 gross 1841 self.getDomain().addPDEToRHS(operator, \
1688 gross 1204 escript.Data(), \
1689     D_reduced_times_e, \
1690     d_reduced_times_e,\
1691     escript.Data())
1692 jgs 148 self.trace("New lumped operator has been built.")
1693 gross 1841 if not self.isRightHandSideValid():
1694     self.resetRightHandSide()
1695     righthandside=self.getCurrentRightHandSide()
1696     self.getDomain().addPDEToRHS(righthandside, \
1697     self.getCoefficient("X"), \
1698     self.getCoefficient("Y"),\
1699     self.getCoefficient("y"),\
1700     self.getCoefficient("y_contact"))
1701     self.getDomain().addPDEToRHS(righthandside, \
1702     self.getCoefficient("X_reduced"), \
1703     self.getCoefficient("Y_reduced"),\
1704     self.getCoefficient("y_reduced"),\
1705     self.getCoefficient("y_contact_reduced"))
1706 jgs 148 self.trace("New right hand side as been built.")
1707 gross 1841 self.validRightHandSide()
1708     self.insertConstraint(rhs_only=False)
1709     self.validOperator()
1710 jgs 108 else:
1711 gross 1841 if not self.isOperatorValid() and not self.isRightHandSideValid():
1712     self.resetRightHandSide()
1713     righthandside=self.getCurrentRightHandSide()
1714     self.resetOperator()
1715     operator=self.getCurrentOperator()
1716     self.getDomain().addPDEToSystem(operator,righthandside, \
1717     self.getCoefficient("A"), \
1718     self.getCoefficient("B"), \
1719     self.getCoefficient("C"), \
1720     self.getCoefficient("D"), \
1721     self.getCoefficient("X"), \
1722     self.getCoefficient("Y"), \
1723     self.getCoefficient("d"), \
1724     self.getCoefficient("y"), \
1725     self.getCoefficient("d_contact"), \
1726     self.getCoefficient("y_contact"))
1727     self.getDomain().addPDEToSystem(operator,righthandside, \
1728     self.getCoefficient("A_reduced"), \
1729     self.getCoefficient("B_reduced"), \
1730     self.getCoefficient("C_reduced"), \
1731     self.getCoefficient("D_reduced"), \
1732     self.getCoefficient("X_reduced"), \
1733     self.getCoefficient("Y_reduced"), \
1734     self.getCoefficient("d_reduced"), \
1735     self.getCoefficient("y_reduced"), \
1736     self.getCoefficient("d_contact_reduced"), \
1737     self.getCoefficient("y_contact_reduced"))
1738     self.insertConstraint(rhs_only=False)
1739 jgs 148 self.trace("New system has been built.")
1740 gross 1841 self.validOperator()
1741     self.validRightHandSide()
1742     elif not self.isRightHandSideValid():
1743     self.resetRightHandSide()
1744     righthandside=self.getCurrentRightHandSide()
1745     self.getDomain().addPDEToRHS(righthandside,
1746     self.getCoefficient("X"), \
1747     self.getCoefficient("Y"),\
1748     self.getCoefficient("y"),\
1749     self.getCoefficient("y_contact"))
1750     self.getDomain().addPDEToRHS(righthandside,
1751     self.getCoefficient("X_reduced"), \
1752     self.getCoefficient("Y_reduced"),\
1753     self.getCoefficient("y_reduced"),\
1754     self.getCoefficient("y_contact_reduced"))
1755     self.insertConstraint(rhs_only=True)
1756 jgs 148 self.trace("New right hand side has been built.")
1757 gross 1841 self.validRightHandSide()
1758     elif not self.isOperatorValid():
1759     self.resetOperator()
1760     operator=self.getCurrentOperator()
1761     self.getDomain().addPDEToSystem(operator,escript.Data(), \
1762     self.getCoefficient("A"), \
1763     self.getCoefficient("B"), \
1764     self.getCoefficient("C"), \
1765     self.getCoefficient("D"), \
1766 jgs 108 escript.Data(), \
1767     escript.Data(), \
1768 gross 1841 self.getCoefficient("d"), \
1769 jgs 108 escript.Data(),\
1770 gross 1841 self.getCoefficient("d_contact"), \
1771 jgs 108 escript.Data())
1772 gross 1841 self.getDomain().addPDEToSystem(operator,escript.Data(), \
1773     self.getCoefficient("A_reduced"), \
1774     self.getCoefficient("B_reduced"), \
1775     self.getCoefficient("C_reduced"), \
1776     self.getCoefficient("D_reduced"), \
1777 gross 1072 escript.Data(), \
1778     escript.Data(), \
1779 gross 1841 self.getCoefficient("d_reduced"), \
1780 gross 1072 escript.Data(),\
1781 gross 1841 self.getCoefficient("d_contact_reduced"), \
1782 gross 1072 escript.Data())
1783 gross 1841 self.insertConstraint(rhs_only=False)
1784 jgs 148 self.trace("New operator has been built.")
1785 gross 1841 self.validOperator()
1786     return (self.getCurrentOperator(), self.getCurrentRightHandSide())
1787 jgs 102
1788 gross 1841 def insertConstraint(self, rhs_only=False):
1789     """
1790     applies the constraints defined by q and r to the PDE
1791    
1792     @param rhs_only: if True only the right hand side is altered by the constraint.
1793     @type rhs_only: C{bool}
1794     """
1795     q=self.getCoefficient("q")
1796     r=self.getCoefficient("r")
1797     righthandside=self.getCurrentRightHandSide()
1798     operator=self.getCurrentOperator()
1799 jgs 102
1800 gross 1841 if not q.isEmpty():
1801     if r.isEmpty():
1802     r_s=self.createSolution()
1803     else:
1804     r_s=r
1805     if not rhs_only and not operator.isEmpty():
1806     if self.isUsingLumping():
1807     operator.copyWithMask(escript.Data(1.,q.getShape(),q.getFunctionSpace()),q)
1808     else:
1809     row_q=escript.Data(q,self.getFunctionSpaceForEquation())
1810     col_q=escript.Data(q,self.getFunctionSpaceForSolution())
1811     u=self.createSolution()
1812     u.copyWithMask(r_s,col_q)
1813     righthandside-=operator*u
1814     operator.nullifyRowsAndCols(row_q,col_q,1.)
1815     righthandside.copyWithMask(r_s,q)
1816    
1817     def setValue(self,**coefficients):
1818     """
1819     sets new values to coefficients
1820    
1821     @param coefficients: new values assigned to coefficients
1822     @keyword A: value for coefficient A.
1823     @type A: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1824     @keyword A_reduced: value for coefficient A_reduced.
1825     @type A_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1826     @keyword B: value for coefficient B
1827     @type B: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1828     @keyword B_reduced: value for coefficient B_reduced
1829     @type B_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1830     @keyword C: value for coefficient C
1831     @type C: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1832     @keyword C_reduced: value for coefficient C_reduced
1833     @type C_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1834     @keyword D: value for coefficient D
1835     @type D: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1836     @keyword D_reduced: value for coefficient D_reduced
1837     @type D_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1838     @keyword X: value for coefficient X
1839     @type X: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1840     @keyword X_reduced: value for coefficient X_reduced
1841     @type X_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1842     @keyword Y: value for coefficient Y
1843     @type Y: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1844     @keyword Y_reduced: value for coefficient Y_reduced
1845     @type Y_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.Function>}.
1846     @keyword d: value for coefficient d
1847     @type d: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
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{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
1850     @keyword y: value for coefficient y
1851     @type y: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1852     @keyword d_contact: value for coefficient d_contact
1853     @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>}.
1854     @keyword d_contact_reduced: value for coefficient d_contact_reduced
1855     @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>}.
1856     @keyword y_contact: value for coefficient y_contact
1857     @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>}.
1858     @keyword y_contact_reduced: value for coefficient y_contact_reduced
1859     @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>}.
1860     @keyword r: values prescribed to the solution at the locations of constraints
1861     @type r: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1862     depending of reduced order is used for the solution.
1863     @keyword q: mask for location of constraints
1864     @type q: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1865     depending of reduced order is used for the representation of the equation.
1866     @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1867     """
1868     super(LinearPDE,self).setValue(**coefficients)
1869     # check if the systrem is inhomogeneous:
1870     if len(coefficients)>0 and not self.isUsingLumping():
1871     q=self.getCoefficient("q")
1872     r=self.getCoefficient("r")
1873     if not q.isEmpty() and not r.isEmpty():
1874     if util.Lsup(q*r)>0.:
1875     self.trace("Inhomogeneous constraint detected.")
1876     self.invalidateSystem()
1877    
1878    
1879     def getResidual(self,u=None):
1880     """
1881     return the residual of u or the current solution if u is not present.
1882    
1883     @param u: argument in the residual calculation. It must be representable in C{elf.getFunctionSpaceForSolution()}. If u is not present or equals L{None}
1884     the current solution is used.
1885     @type u: L{Data<escript.Data>} or None
1886     @return: residual of u
1887     @rtype: L{Data<escript.Data>}
1888     """
1889     if u==None:
1890     return self.getOperator()*self.getSolution()-self.getRightHandSide()
1891     else:
1892     return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())-self.getRightHandSide()
1893    
1894     def getFlux(self,u=None):
1895     """
1896     returns the flux M{J} for a given M{u}
1897    
1898     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]}
1899    
1900     or
1901    
1902     M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[l]+(B[j]+B_reduced[j])u-X[j]-X_reduced[j]}
1903    
1904     @param u: argument in the flux. If u is not present or equals L{None} the current solution is used.
1905     @type u: L{Data<escript.Data>} or None
1906     @return: flux
1907     @rtype: L{Data<escript.Data>}
1908     """
1909     if u==None: u=self.getSolution()
1910     return util.tensormult(self.getCoefficient("A"),util.grad(u,Funtion(self.getDomain))) \
1911     +util.matrixmult(self.getCoefficient("B"),u) \
1912     -util.self.getCoefficient("X") \
1913     +util.tensormult(self.getCoefficient("A_reduced"),util.grad(u,ReducedFuntion(self.getDomain))) \
1914     +util.matrixmult(self.getCoefficient("B_reduced"),u) \
1915     -util.self.getCoefficient("X_reduced")
1916    
1917    
1918 jgs 149 class Poisson(LinearPDE):
1919     """
1920     Class to define a Poisson equation problem, which is genear L{LinearPDE} of the form
1921 jgs 102
1922 jgs 149 M{-grad(grad(u)[j])[j] = f}
1923 jgs 102
1924 jgs 149 with natural boundary conditons
1925    
1926     M{n[j]*grad(u)[j] = 0 }
1927    
1928     and constraints:
1929    
1930     M{u=0} where M{q>0}
1931    
1932 jgs 108 """
1933 jgs 148
1934 gross 328 def __init__(self,domain,debug=False):
1935 jgs 149 """
1936     initializes a new Poisson equation
1937 jgs 104
1938 jgs 149 @param domain: domain of the PDE
1939     @type domain: L{Domain<escript.Domain>}
1940     @param debug: if True debug informations are printed.
1941 jgs 102
1942 jgs 149 """
1943 jgs 151 super(Poisson, self).__init__(domain,1,1,debug)
1944 gross 1841 self.introduceCoefficients(
1945     f=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1946     f_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE))
1947 jgs 149 self.setSymmetryOn()
1948 jgs 102
1949 jgs 149 def setValue(self,**coefficients):
1950     """
1951     sets new values to coefficients
1952 jgs 102
1953 jgs 149 @param coefficients: new values assigned to coefficients
1954     @keyword f: value for right hand side M{f}
1955     @type f: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1956     @keyword q: mask for location of constraints
1957     @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>}
1958     depending of reduced order is used for the representation of the equation.
1959     @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1960     """
1961 jgs 151 super(Poisson, self).setValue(**coefficients)
1962 jgs 102
1963 gross 1841
1964     def getCoefficient(self,name):
1965 jgs 149 """
1966     return the value of the coefficient name of the general PDE
1967     @param name: name of the coefficient requested.
1968     @type name: C{string}
1969     @return: the value of the coefficient name
1970     @rtype: L{Data<escript.Data>}
1971 gross 1841 @raise IllegalCoefficient: invalid coefficient name
1972 jgs 149 @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
1973     """
1974     if name == "A" :
1975 gross 328 return escript.Data(util.kronecker(self.getDim()),escript.Function(self.getDomain()))
1976 jgs 149 elif name == "Y" :
1977     return self.getCoefficient("f")
1978 gross 1072 elif name == "Y_reduced" :
1979     return self.getCoefficient("f_reduced")
1980 jgs 149 else:
1981 gross 1841 return super(Poisson, self).getCoefficient(name)
1982 jgs 102
1983 jgs 149 class Helmholtz(LinearPDE):
1984     """
1985     Class to define a Helmhotz equation problem, which is genear L{LinearPDE} of the form
1986    
1987     M{S{omega}*u - grad(k*grad(u)[j])[j] = f}
1988    
1989     with natural boundary conditons
1990    
1991     M{k*n[j]*grad(u)[j] = g- S{alpha}u }
1992    
1993 jgs 122 and constraints:
1994 jgs 102
1995 jgs 149 M{u=r} where M{q>0}
1996    
1997 jgs 110 """
1998 jgs 149
1999     def __init__(self,domain,debug=False):
2000     """
2001     initializes a new Poisson equation
2002    
2003     @param domain: domain of the PDE
2004     @type domain: L{Domain<escript.Domain>}
2005     @param debug: if True debug informations are printed.
2006    
2007     """
2008 jgs 151 super(Helmholtz, self).__init__(domain,1,1,debug)
2009 gross 1841 self.introduceCoefficients(
2010     omega=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.OPERATOR),
2011     k=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.OPERATOR),
2012     f=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2013     f_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2014     alpha=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.OPERATOR),
2015     g=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2016     g_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE))
2017 jgs 149 self.setSymmetryOn()
2018    
2019 gross 1841 def setValue(self,**coefficients):
2020 jgs 149 """
2021     sets new values to coefficients
2022    
2023     @param coefficients: new values assigned to coefficients
2024     @keyword omega: value for coefficient M{S{omega}}
2025     @type omega: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2026     @keyword k: value for coefficeint M{k}
2027     @type k: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2028     @keyword f: value for right hand side M{f}
2029     @type f: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2030     @keyword alpha: value for right hand side M{S{alpha}}
2031     @type alpha: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
2032     @keyword g: value for right hand side M{g}
2033     @type g: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
2034     @keyword r: prescribed values M{r} for the solution in constraints.
2035     @type r: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2036     depending of reduced order is used for the representation of the equation.
2037     @keyword q: mask for location of constraints
2038     @type q: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2039     depending of reduced order is used for the representation of the equation.
2040     @raise IllegalCoefficient: if an unknown coefficient keyword is used.
2041     """
2042 jgs 151 super(Helmholtz, self).setValue(**coefficients)
2043 jgs 149
2044 gross 1841 def getCoefficient(self,name):
2045 jgs 149 """
2046     return the value of the coefficient name of the general PDE
2047    
2048     @param name: name of the coefficient requested.
2049     @type name: C{string}
2050     @return: the value of the coefficient name
2051     @rtype: L{Data<escript.Data>}
2052 gross 1841 @raise IllegalCoefficient: invalid name
2053 jgs 149 """
2054     if name == "A" :
2055     return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))*self.getCoefficient("k")
2056     elif name == "D" :
2057     return self.getCoefficient("omega")
2058     elif name == "Y" :
2059     return self.getCoefficient("f")
2060     elif name == "d" :
2061     return self.getCoefficient("alpha")
2062     elif name == "y" :
2063     return self.getCoefficient("g")
2064 gross 1072 elif name == "Y_reduced" :
2065     return self.getCoefficient("f_reduced")
2066     elif name == "y_reduced" :
2067     return self.getCoefficient("g_reduced")
2068 jgs 149 else:
2069 gross 1841 return super(Helmholtz, self).getCoefficient(name)
2070 jgs 149
2071     class LameEquation(LinearPDE):
2072     """
2073     Class to define a Lame equation problem:
2074    
2075 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] }
2076 jgs 149
2077     with natural boundary conditons:
2078    
2079 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] }
2080 jgs 149
2081     and constraints:
2082    
2083     M{u[i]=r[i]} where M{q[i]>0}
2084    
2085     """
2086    
2087     def __init__(self,domain,debug=False):
2088 jgs 151 super(LameEquation, self).__init__(domain,\
2089     domain.getDim(),domain.getDim(),debug)
2090 gross 1841 self.introduceCoefficients(lame_lambda=PDECoef(PDECoef.INTERIOR,(),PDECoef.OPERATOR),
2091     lame_mu=PDECoef(PDECoef.INTERIOR,(),PDECoef.OPERATOR),
2092     F=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2093     sigma=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
2094     f=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE))
2095 jgs 151 self.setSymmetryOn()
2096 jgs 149
2097 gross 1841 def setValues(self,**coefficients):
2098 jgs 149 """
2099     sets new values to coefficients
2100    
2101     @param coefficients: new values assigned to coefficients
2102     @keyword lame_mu: value for coefficient M{S{mu}}
2103     @type lame_mu: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2104     @keyword lame_lambda: value for coefficient M{S{lambda}}
2105     @type lame_lambda: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2106     @keyword F: value for internal force M{F}
2107     @type F: any type that can be casted to L{Vector<escript.Vector>} object on L{Function<escript.Function>}
2108     @keyword sigma: value for initial stress M{S{sigma}}
2109     @type sigma: any type that can be casted to L{Tensor<escript.Tensor>} object on L{Function<escript.Function>}
2110     @keyword f: value for extrenal force M{f}
2111     @type f: any type that can be casted to L{Vector<escript.Vector>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}
2112     @keyword r: prescribed values M{r} for the solution in constraints.
2113     @type r: any type that can be casted to L{Vector<escript.Vector>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2114     depending of reduced order is used for the representation of the equation.
2115     @keyword q: mask for location of constraints
2116     @type q: any type that can be casted to L{Vector<escript.Vector>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2117     depending of reduced order is used for the representation of the equation.
2118     @raise IllegalCoefficient: if an unknown coefficient keyword is used.
2119     """
2120 gross 596 super(LameEquation, self).setValues(**coefficients)
2121 jgs 149
2122 gross 1841 def getCoefficient(self,name):
2123 jgs 149 """
2124     return the value of the coefficient name of the general PDE
2125    
2126     @param name: name of the coefficient requested.
2127     @type name: C{string}
2128     @return: the value of the coefficient name
2129     @rtype: L{Data<escript.Data>}
2130 gross 1841 @raise IllegalCoefficient: invalid coefficient name
2131 jgs 149 """
2132     if name == "A" :
2133 gross 1841 out =self.createCoefficient("A")
2134 jgs 149 for i in range(self.getDim()):
2135     for j in range(self.getDim()):
2136     out[i,i,j,j] += self.getCoefficient("lame_lambda")
2137     out[i,j,j,i] += self.getCoefficient("lame_mu")
2138     out[i,j,i,j] += self.getCoefficient("lame_mu")
2139     return out
2140     elif name == "X" :
2141     return self.getCoefficient("sigma")
2142     elif name == "Y" :
2143     return self.getCoefficient("F")
2144     elif name == "y" :
2145     return self.getCoefficient("f")
2146     else:
2147 gross 1841 return super(LameEquation, self).getCoefficient(name)
2148 jgs 149
2149 gross 1400 def LinearSinglePDE(domain,debug=False):
2150     """
2151     defines a single linear PDEs
2152    
2153     @param domain: domain of the PDE
2154     @type domain: L{Domain<escript.Domain>}
2155     @param debug: if True debug informations are printed.
2156     @rtype: L{LinearPDE}
2157     """
2158     return LinearPDE(domain,numEquations=1,numSolutions=1,debug=debug)
2159    
2160     def LinearPDESystem(domain,debug=False):
2161     """
2162     defines a system of linear PDEs
2163    
2164     @param domain: domain of the PDE
2165     @type domain: L{Domain<escript.Domain>}
2166     @param debug: if True debug informations are printed.
2167     @rtype: L{LinearPDE}
2168     """
2169     return LinearPDE(domain,numEquations=domain.getDim(),numSolutions=domain.getDim(),debug=debug)
2170 gross 1410
2171 gross 1841 class TransportPDE(LinearProblem):
2172     """
2173     This class is used to define a transport problem given by a general linear, time dependend, second order PDE
2174     for an unknown, non-negative, time-dependend function M{u} on a given domain defined through a L{Domain<escript.Domain>} object.
2175    
2176     For a single equation with a solution with a single component the transport problem is defined in the following form:
2177    
2178     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)}
2179    
2180    
2181     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,
2182     ie. summation over indexes appearing twice in a term of a sum is performed, is used.
2183     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
2184     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
2185     L{ReducedFunction<escript.ReducedFunction>}. It is also allowed to use objects that can be converted into
2186     such L{Data<escript.Data>} objects. M{M} and M{M_reduced} are scalar, M{A} and M{A_reduced} are rank two,
2187     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.
2188    
2189     The following natural boundary conditions are considered:
2190    
2191     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}
2192    
2193     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>}.
2194    
2195     Constraints for the solution prescribing the value of the solution at certain locations in the domain. They have the form
2196    
2197     M{u_t=r} where M{q>0}
2198    
2199     M{r} and M{q} are each scalar where M{q} is the characteristic function defining where the constraint is applied.
2200     The constraints override any other condition set by the transport problem or the boundary condition.
2201    
2202     The transport problem is symmetrical if
2203    
2204     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]}
2205    
2206     For a system and a solution with several components the transport problem has the form
2207    
2208     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] }
2209    
2210     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:
2211    
2212     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}
2213    
2214     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>}.
2215    
2216     Constraints take the form
2217    
2218     M{u[i]_t=r[i]} where M{q[i]>0}
2219    
2220     M{r} and M{q} are each rank one. Notice that at some locations not necessarily all components must have a constraint.
2221    
2222     The transport problem is symmetrical if
2223    
2224     - M{M[i,k]=M[i,k]}
2225     - M{M_reduced[i,k]=M_reduced[i,k]}
2226     - M{A[i,j,k,l]=A[k,l,i,j]}
2227     - M{A_reduced[i,j,k,l]=A_reduced[k,l,i,j]}
2228     - M{B[i,j,k]=C[k,i,j]}
2229     - M{B_reduced[i,j,k]=C_reduced[k,i,j]}
2230     - M{D[i,k]=D[i,k]}
2231     - M{D_reduced[i,k]=D_reduced[i,k]}
2232     - M{m[i,k]=m[k,i]}
2233     - M{m_reduced[i,k]=m_reduced[k,i]}
2234     - M{d[i,k]=d[k,i]}
2235     - M{d_reduced[i,k]=d_reduced[k,i]}
2236    
2237     L{TransportPDE} also supports solution discontinuities over a contact region in the domain. To specify the conditions across the
2238     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
2239     defined as
2240    
2241     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]}
2242    
2243     For the case of single solution component and single PDE M{J} is defined
2244    
2245     M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[j]+(B[i]+B_reduced[i])*u+X[i]+X_reduced[i]}
2246    
2247     In the context of discontinuities M{n} denotes the normal on the discontinuity pointing from side 0 towards side 1
2248     calculated from L{getNormal<escript.FunctionSpace.getNormal>} of L{FunctionOnContactZero<escript.FunctionOnContactZero>}. For a system of transport problems the contact condition takes the form
2249    
2250     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]}
2251    
2252     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
2253     of the solution at side 1 and at side 0, denotes the jump of M{u} across discontinuity along the normal calculated by
2254     L{jump<util.jump>}.
2255     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>}.
2256     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>}.
2257     In case of a single PDE and a single component solution the contact condition takes the form
2258    
2259     M{n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)}
2260    
2261     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>}
2262    
2263     typical usage:
2264    
2265     p=TransportPDE(dom)
2266     p.setValue(M=1.,C=[-1.,0.])
2267     p.setInitialSolution(u=exp(-length(dom.getX()-[0.1,0.1])**2)
2268     t=0
2269     dt=0.1
2270     while (t<1.):
2271     u=p.solve(dt)
2272    
2273     """
2274 gross 1859 def __init__(self,domain,numEquations=None,numSolutions=None, theta=0.5,debug=False):
2275 gross 1410 """
2276 gross 1841 initializes a linear problem
2277 gross 1410
2278 gross 1841 @param domain: domain of the PDE
2279     @type domain: L{Domain<escript.Domain>}
2280     @param numEquations: number of equations. If numEquations==None the number of equations
2281     is extracted from the coefficients.
2282     @param numSolutions: number of solution components. If numSolutions==None the number of solution components
2283     is extracted from the coefficients.
2284     @param debug: if True debug informations are printed.
2285     @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>}.
2286 gross 1410
2287 gross 1841 """
2288     if theta<0 or theta>1:
2289 gross 1859 raise ValueError,"theta (=%s) needs to be between 0 and 1 (inclusive)."%theta
2290     super(TransportPDE, self).__init__(domain,numEquations,numSolutions,debug)
2291 gross 1410
2292 gross 1841 self.__theta=theta
2293     #
2294     # the coefficients of the transport problem
2295     #
2296     self.introduceCoefficients(
2297     M=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2298     A=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2299     B=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2300     C=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2301     D=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2302     X=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
2303     Y=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2304     m=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2305     d=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2306     y=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2307     d_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2308     y_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2309     M_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2310     A_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2311     B_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2312     C_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2313     D_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2314     X_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
2315     Y_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2316     m_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2317     d_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2318     y_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2319     d_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2320     y_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2321     r=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.RIGHTHANDSIDE),
2322     q=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.BOTH) )
2323 gross 1410
2324 gross 1841 def __str__(self):
2325 gross 1410 """
2326 gross 1841 returns string representation of the transport problem
2327 gross 1410
2328 gross 1841 @return: a simple representation of the transport problem
2329     @rtype: C{str}
2330     """
2331     return "<TransportPDE %d>"%id(self)
2332 gross 1417
2333 gross 1841 def getTheta(self):
2334     """
2335     returns the time stepping control parameter
2336     @rtype: float
2337     """
2338     return self.__theta
2339 gross 1410
2340 gross 1859
2341 gross 1841 def checkSymmetry(self,verbose=True):
2342     """
2343     test the transport problem for symmetry.
2344    
2345     @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.
2346     @type verbose: C{bool}
2347     @return: True if the PDE is symmetric.
2348     @rtype: C{bool}
2349     @note: This is a very expensive operation. It should be used for degugging only! The symmetry flag is not altered.
2350     """
2351     out=True
2352     out=out and self.checkSymmetricTensor("M", verbose)
2353     out=out and self.checkSymmetricTensor("M_reduced", verbose)
2354     out=out and self.checkSymmetricTensor("A", verbose)
2355     out=out and self.checkSymmetricTensor("A_reduced", verbose)
2356     out=out and self.checkReciprocalSymmetry("B","C", verbose)
2357     out=out and self.checkReciprocalSymmetry("B_reduced","C_reduced", verbose)
2358     out=out and self.checkSymmetricTensor("D", verbose)
2359     out=out and self.checkSymmetricTensor("D_reduced", verbose)
2360     out=out and self.checkSymmetricTensor("m", verbose)
2361     out=out and self.checkSymmetricTensor("m_reduced", verbose)
2362     out=out and self.checkSymmetricTensor("d", verbose)
2363     out=out and self.checkSymmetricTensor("d_reduced", verbose)
2364     out=out and self.checkSymmetricTensor("d_contact", verbose)
2365     out=out and self.checkSymmetricTensor("d_contact_reduced", verbose)
2366     return out
2367    
2368     def setValue(self,**coefficients):
2369     """
2370     sets new values to coefficients
2371    
2372     @param coefficients: new values assigned to coefficients
2373     @keyword M: value for coefficient M.
2374     @type M: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
2375     @keyword M_reduced: value for coefficient M_reduced.
2376     @type M_reduced: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.ReducedFunction>}.
2377     @keyword A: value for coefficient A.
2378     @type A: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
2379     @keyword A_reduced: value for coefficient A_reduced.
2380     @type A_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
2381     @keyword B: value for coefficient B
2382     @type B: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
2383     @keyword B_reduced: value for coefficient B_reduced
2384     @type B_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
2385     @keyword C: value for coefficient C
2386     @type C: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
2387     @keyword C_reduced: value for coefficient C_reduced
2388     @type C_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
2389     @keyword D: value for coefficient D
2390     @type D: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
2391     @keyword D_reduced: value for coefficient D_reduced
2392     @type D_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
2393     @keyword X: value for coefficient X
2394     @type X: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
2395     @keyword X_reduced: value for coefficient X_reduced
2396     @type X_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
2397     @keyword Y: value for coefficient Y
2398     @type Y: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
2399     @keyword Y_reduced: value for coefficient Y_reduced
2400     @type Y_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.Function>}.
2401     @keyword m: value for coefficient m
2402     @type m: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
2403     @keyword m_reduced: value for coefficient m_reduced