/[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 2431 - (hide annotations)
Wed May 20 03:39:57 2009 UTC (10 years, 4 months ago) by jfenwick
File MIME type: text/x-python
File size: 124242 byte(s)
Added at the request of Dian.

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