/[escript]/trunk/escript/py_src/linearPDEs.py
ViewVC logotype

Annotation of /trunk/escript/py_src/linearPDEs.py

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26