/[escript]/branches/symbolic_from_3470/escript/py_src/linearPDEs.py
ViewVC logotype

Annotation of /branches/symbolic_from_3470/escript/py_src/linearPDEs.py

Parent Directory Parent Directory | Revision Log Revision Log


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