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