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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1072 - (hide annotations)
Thu Mar 29 06:44:30 2007 UTC (12 years, 5 months ago) by gross
File MIME type: text/x-python
File size: 102599 byte(s)
PDE assemblage for reduced integration order + tests added.


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26