/[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 149 - (hide annotations)
Thu Sep 1 03:31:39 2005 UTC (14 years, 3 months ago) by jgs
Original Path: trunk/esys2/escript/py_src/linearPDEs.py
File MIME type: text/x-python
File size: 100020 byte(s)
Merge of development branch dev-02 back to main trunk on 2005-09-01

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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26