/[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 1639 - (hide annotations)
Mon Jul 14 08:55:25 2008 UTC (11 years, 3 months ago) by gross
File MIME type: text/x-python
File size: 112746 byte(s)


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