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