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