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

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

Parent Directory Parent Directory | Revision Log Revision Log


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

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26