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