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