/[escript]/trunk-mpi-branch/escript/py_src/linearPDEs.py
ViewVC logotype

Annotation of /trunk-mpi-branch/escript/py_src/linearPDEs.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1302 - (hide annotations)
Mon Sep 17 06:57:51 2007 UTC (12 years ago) by ksteube
File MIME type: text/x-python
File size: 102287 byte(s)
Some of the epydoc doesn't show up or doesn't get formatted, in this
case because of a missing } or because of improper indentation.

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26