/[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 1809 - (hide annotations)
Thu Sep 25 06:43:44 2008 UTC (11 years, 2 months ago) by ksteube
File MIME type: text/x-python
File size: 112970 byte(s)
Copyright updated in all python files

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