/[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 1137 - (hide annotations)
Thu May 10 08:11:31 2007 UTC (12 years, 4 months ago) by gross
File MIME type: text/x-python
File size: 102616 byte(s)
This version passes the tests on windows except for 

   * vtk
   * netCDF

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

Properties

Name Value
svn:eol-style native
svn:keywords Author Date Id Revision

  ViewVC Help
Powered by ViewVC 1.1.26