/[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 2100 - (show annotations)
Wed Nov 26 08:13:00 2008 UTC (10 years, 9 months ago) by gross
File MIME type: text/x-python
File size: 118017 byte(s)
This commit cleans up the incompressible solver and adds a DarcyFlux solver in model module. 
Some documentation for both classes has been added.
The convection code is only linear at the moment.



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 setDebug(self, flag):
482 """
483 switches debugging to on if C{flag} is true outherwise debugging is set to off
484
485 @param flag: desired debug status
486 @type flag: C{bool}
487 """
488 if flag:
489 self.setDebugOn()
490 else:
491 self.setDebugOff()
492
493 def trace(self,text):
494 """
495 print the text message if debugging is swiched on.
496 @param text: message
497 @type text: C{string}
498 """
499 if self.__debug: print "%s: %s"%(str(self),text)
500
501 # =============================================================================
502 # some service functions:
503 # =============================================================================
504 def introduceCoefficients(self,**coeff):
505 """
506 introduces a new coefficient into the problem.
507
508 use
509
510 self.introduceCoefficients(A=PDECoef(...),B=PDECoef(...))
511
512 to introduce the coefficients M{A} ans M{B}.
513 """
514 for name, type in coeff.items():
515 if not isinstance(type,PDECoef):
516 raise ValueError,"coefficient %s has no type."%name
517 self.__COEFFICIENTS[name]=type
518 self.__COEFFICIENTS[name].resetValue()
519 self.trace("coefficient %s has been introduced."%name)
520
521 def getDomain(self):
522 """
523 returns the domain of the PDE
524
525 @return: the domain of the PDE
526 @rtype: L{Domain<escript.Domain>}
527 """
528 return self.__domain
529
530 def getDim(self):
531 """
532 returns the spatial dimension of the PDE
533
534 @return: the spatial dimension of the PDE domain
535 @rtype: C{int}
536 """
537 return self.getDomain().getDim()
538
539 def getNumEquations(self):
540 """
541 returns the number of equations
542
543 @return: the number of equations
544 @rtype: C{int}
545 @raise UndefinedPDEError: if the number of equations is not be specified yet.
546 """
547 if self.__numEquations==None:
548 if self.__numSolutions==None:
549 raise UndefinedPDEError,"Number of equations is undefined. Please specify argument numEquations."
550 else:
551 self.__numEquations=self.__numSolutions
552 return self.__numEquations
553
554 def getNumSolutions(self):
555 """
556 returns the number of unknowns
557
558 @return: the number of unknowns
559 @rtype: C{int}
560 @raise UndefinedPDEError: if the number of unknowns is not be specified yet.
561 """
562 if self.__numSolutions==None:
563 if self.__numEquations==None:
564 raise UndefinedPDEError,"Number of solution is undefined. Please specify argument numSolutions."
565 else:
566 self.__numSolutions=self.__numEquations
567 return self.__numSolutions
568
569 def reduceEquationOrder(self):
570 """
571 return status for order reduction for equation
572
573 @return: return True is reduced interpolation order is used for the represenation of the equation
574 @rtype: L{bool}
575 """
576 return self.__reduce_equation_order
577
578 def reduceSolutionOrder(self):
579 """
580 return status for order reduction for the solution
581
582 @return: return True is reduced interpolation order is used for the represenation of the solution
583 @rtype: L{bool}
584 """
585 return self.__reduce_solution_order
586
587 def getFunctionSpaceForEquation(self):
588 """
589 returns the L{FunctionSpace<escript.FunctionSpace>} used to discretize the equation
590
591 @return: representation space of equation
592 @rtype: L{FunctionSpace<escript.FunctionSpace>}
593 """
594 if self.reduceEquationOrder():
595 return escript.ReducedSolution(self.getDomain())
596 else:
597 return escript.Solution(self.getDomain())
598
599 def getFunctionSpaceForSolution(self):
600 """
601 returns the L{FunctionSpace<escript.FunctionSpace>} used to represent the solution
602
603 @return: representation space of solution
604 @rtype: L{FunctionSpace<escript.FunctionSpace>}
605 """
606 if self.reduceSolutionOrder():
607 return escript.ReducedSolution(self.getDomain())
608 else:
609 return escript.Solution(self.getDomain())
610
611 # =============================================================================
612 # solver settings:
613 # =============================================================================
614 def setSolverMethod(self,solver=None,preconditioner=None):
615 """
616 sets a new solver
617
618 @param solver: sets a new solver method.
619 @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}
620 @param preconditioner: sets a new solver method.
621 @type preconditioner: one of L{DEFAULT}, L{JACOBI} L{ILU0}, L{ILUT},L{SSOR}, L{RILU}, L{GS}
622
623 @note: the solver method choosen may not be available by the selected package.
624 """
625 if solver==None: solver=self.__solver_method
626 if preconditioner==None: preconditioner=self.__preconditioner
627 if solver==None: solver=self.DEFAULT
628 if preconditioner==None: preconditioner=self.DEFAULT
629 if not (solver,preconditioner)==self.getSolverMethod():
630 self.__solver_method=solver
631 self.__preconditioner=preconditioner
632 self.updateSystemType()
633 self.trace("New solver is %s"%self.getSolverMethodName())
634
635 def getSolverMethodName(self):
636 """
637 returns the name of the solver currently used
638
639 @return: the name of the solver currently used.
640 @rtype: C{string}
641 """
642
643 m=self.getSolverMethod()
644 p=self.getSolverPackage()
645 method=""
646 if m[0]==self.DEFAULT: method="DEFAULT"
647 elif m[0]==self.DIRECT: method= "DIRECT"
648 elif m[0]==self.ITERATIVE: method= "ITERATIVE"
649 elif m[0]==self.CHOLEVSKY: method= "CHOLEVSKY"
650 elif m[0]==self.PCG: method= "PCG"
651 elif m[0]==self.TFQMR: method= "TFQMR"
652 elif m[0]==self.MINRES: method= "MINRES"
653 elif m[0]==self.CR: method= "CR"
654 elif m[0]==self.CGS: method= "CGS"
655 elif m[0]==self.BICGSTAB: method= "BICGSTAB"
656 elif m[0]==self.SSOR: method= "SSOR"
657 elif m[0]==self.GMRES: method= "GMRES"
658 elif m[0]==self.PRES20: method= "PRES20"
659 elif m[0]==self.LUMPING: method= "LUMPING"
660 elif m[0]==self.AMG: method= "AMG"
661 if m[1]==self.DEFAULT: method+="+DEFAULT"
662 elif m[1]==self.JACOBI: method+= "+JACOBI"
663 elif m[1]==self.ILU0: method+= "+ILU0"
664 elif m[1]==self.ILUT: method+= "+ILUT"
665 elif m[1]==self.SSOR: method+= "+SSOR"
666 elif m[1]==self.AMG: method+= "+AMG"
667 elif m[1]==self.RILU: method+= "+RILU"
668 elif m[1]==self.GS: method+= "+GS"
669 if p==self.DEFAULT: package="DEFAULT"
670 elif p==self.PASO: package= "PASO"
671 elif p==self.MKL: package= "MKL"
672 elif p==self.SCSL: package= "SCSL"
673 elif p==self.UMFPACK: package= "UMFPACK"
674 elif p==self.TRILINOS: package= "TRILINOS"
675 else : method="unknown"
676 return "%s solver of %s package"%(method,package)
677
678 def getSolverMethod(self):
679 """
680 returns the solver method and preconditioner used
681
682 @return: the solver method currently be used.
683 @rtype: C{tuple} of C{int}
684 """
685 return self.__solver_method,self.__preconditioner
686
687 def setSolverPackage(self,package=None):
688 """
689 sets a new solver package
690
691 @param package: sets a new solver method.
692 @type package: one of L{DEFAULT}, L{PASO} L{SCSL}, L{MKL}, L{UMFPACK}, L{TRILINOS}
693 """
694 if package==None: package=self.DEFAULT
695 if not package==self.getSolverPackage():
696 self.__solver_package=package
697 self.updateSystemType()
698 self.trace("New solver is %s"%self.getSolverMethodName())
699
700 def getSolverPackage(self):
701 """
702 returns the package of the solver
703
704 @return: the solver package currently being used.
705 @rtype: C{int}
706 """
707 return self.__solver_package
708
709 def isUsingLumping(self):
710 """
711 checks if matrix lumping is used as a solver method
712
713 @return: True is lumping is currently used a solver method.
714 @rtype: C{bool}
715 """
716 return self.getSolverMethod()[0]==self.LUMPING
717
718 def setTolerance(self,tol=1.e-8):
719 """
720 resets the tolerance for the solver method to tol.
721
722 @param tol: new tolerance for the solver. If the tol is lower then the current tolerence
723 the system will be resolved.
724 @type tol: positive C{float}
725 @raise ValueError: if tolerance is not positive.
726 """
727 if not tol>0:
728 raise ValueError,"Tolerance as to be positive"
729 if tol<self.getTolerance(): self.invalidateSolution()
730 self.trace("New tolerance %e"%tol)
731 self.__tolerance=tol
732 return
733
734 def getTolerance(self):
735 """
736 returns the tolerance set for the solution
737
738 @return: tolerance currently used.
739 @rtype: C{float}
740 """
741 return self.__tolerance
742
743 # =============================================================================
744 # symmetry flag:
745 # =============================================================================
746 def isSymmetric(self):
747 """
748 checks if symmetry is indicated.
749
750 @return: True is a symmetric PDE is indicated, otherwise False is returned
751 @rtype: C{bool}
752 """
753 return self.__sym
754
755 def setSymmetryOn(self):
756 """
757 sets the symmetry flag.
758 """
759 if not self.isSymmetric():
760 self.trace("PDE is set to be symmetric")
761 self.__sym=True
762 self.updateSystemType()
763
764 def setSymmetryOff(self):
765 """
766 removes the symmetry flag.
767 """
768 if self.isSymmetric():
769 self.trace("PDE is set to be unsymmetric")
770 self.__sym=False
771 self.updateSystemType()
772
773 def setSymmetryTo(self,flag=False):
774 """
775 sets the symmetry flag to flag
776
777 @param flag: If flag, the symmetry flag is set otherwise the symmetry flag is released.
778 @type flag: C{bool}
779 """
780 if flag:
781 self.setSymmetryOn()
782 else:
783 self.setSymmetryOff()
784
785 # =============================================================================
786 # function space handling for the equation as well as the solution
787 # =============================================================================
788 def setReducedOrderOn(self):
789 """
790 switches on reduced order for solution and equation representation
791
792 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
793 """
794 self.setReducedOrderForSolutionOn()
795 self.setReducedOrderForEquationOn()
796
797 def setReducedOrderOff(self):
798 """
799 switches off reduced order for solution and equation representation
800
801 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
802 """
803 self.setReducedOrderForSolutionOff()
804 self.setReducedOrderForEquationOff()
805
806 def setReducedOrderTo(self,flag=False):
807 """
808 sets order reduction for both solution and equation representation according to flag.
809 @param flag: if flag is True, the order reduction is switched on for both solution and equation representation, otherwise or
810 if flag is not present order reduction is switched off
811 @type flag: C{bool}
812 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
813 """
814 self.setReducedOrderForSolutionTo(flag)
815 self.setReducedOrderForEquationTo(flag)
816
817
818 def setReducedOrderForSolutionOn(self):
819 """
820 switches on reduced order for solution representation
821
822 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
823 """
824 if not self.__reduce_solution_order:
825 if self.__altered_coefficients:
826 raise RuntimeError,"order cannot be altered after coefficients have been defined."
827 self.trace("Reduced order is used to solution representation.")
828 self.__reduce_solution_order=True
829 self.initializeSystem()
830
831 def setReducedOrderForSolutionOff(self):
832 """
833 switches off reduced order for solution representation
834
835 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
836 """
837 if self.__reduce_solution_order:
838 if self.__altered_coefficients:
839 raise RuntimeError,"order cannot be altered after coefficients have been defined."
840 self.trace("Full order is used to interpolate solution.")
841 self.__reduce_solution_order=False
842 self.initializeSystem()
843
844 def setReducedOrderForSolutionTo(self,flag=False):
845 """
846 sets order for test functions according to flag
847
848 @param flag: if flag is True, the order reduction is switched on for solution representation, otherwise or
849 if flag is not present order reduction is switched off
850 @type flag: C{bool}
851 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
852 """
853 if flag:
854 self.setReducedOrderForSolutionOn()
855 else:
856 self.setReducedOrderForSolutionOff()
857
858 def setReducedOrderForEquationOn(self):
859 """
860 switches on reduced order for equation representation
861
862 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
863 """
864 if not self.__reduce_equation_order:
865 if self.__altered_coefficients:
866 raise RuntimeError,"order cannot be altered after coefficients have been defined."
867 self.trace("Reduced order is used for test functions.")
868 self.__reduce_equation_order=True
869 self.initializeSystem()
870
871 def setReducedOrderForEquationOff(self):
872 """
873 switches off reduced order for equation representation
874
875 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
876 """
877 if self.__reduce_equation_order:
878 if self.__altered_coefficients:
879 raise RuntimeError,"order cannot be altered after coefficients have been defined."
880 self.trace("Full order is used for test functions.")
881 self.__reduce_equation_order=False
882 self.initializeSystem()
883
884 def setReducedOrderForEquationTo(self,flag=False):
885 """
886 sets order for test functions according to flag
887
888 @param flag: if flag is True, the order reduction is switched on for equation representation, otherwise or
889 if flag is not present order reduction is switched off
890 @type flag: C{bool}
891 @raise RuntimeError: if order reduction is altered after a coefficient has been set.
892 """
893 if flag:
894 self.setReducedOrderForEquationOn()
895 else:
896 self.setReducedOrderForEquationOff()
897
898 def updateSystemType(self):
899 """
900 reasses the system type and, if a new matrix is needed, resets the system.
901 """
902 new_system_type=self.getRequiredSystemType()
903 if not new_system_type==self.__system_type:
904 self.trace("Matrix type is now %d."%new_system_type)
905 self.__system_type=new_system_type
906 self.initializeSystem()
907 def getSystemType(self):
908 """
909 return the current system type
910 """
911 return self.__system_type
912
913 def checkSymmetricTensor(self,name,verbose=True):
914 """
915 tests a coefficient for symmetry
916
917 @param name: name of the coefficient
918 @type name: C{str}
919 @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.
920 @type verbose: C{bool}
921 @return: True if coefficient C{name} is symmetric.
922 @rtype: C{bool}
923 """
924 A=self.getCoefficient(name)
925 verbose=verbose or self.__debug
926 out=True
927 if not A.isEmpty():
928 tol=util.Lsup(A)*self.SMALL_TOLERANCE
929 s=A.getShape()
930 if A.getRank() == 4:
931 if s[0]==s[2] and s[1] == s[3]:
932 for i in range(s[0]):
933 for j in range(s[1]):
934 for k in range(s[2]):
935 for l in range(s[3]):
936 if util.Lsup(A[i,j,k,l]-A[k,l,i,j])>tol:
937 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)
938 out=False
939 else:
940 if verbose: print "non-symmetric problem because of inappropriate shape %s of coefficient %s."%(s,name)
941 out=False
942 elif A.getRank() == 2:
943 if s[0]==s[1]:
944 for j in range(s[0]):
945 for l in range(s[1]):
946 if util.Lsup(A[j,l]-A[l,j])>tol:
947 if verbose: print "non-symmetric problem because %s[%d,%d]!=%s[%d,%d]"%(name,j,l,name,l,j)
948 out=False
949 else:
950 if verbose: print "non-symmetric problem because of inappropriate shape %s of coefficient %s."%(s,name)
951 out=False
952 elif A.getRank() == 0:
953 pass
954 else:
955 raise ValueError,"Cannot check rank %s of %s."%(A.getRank(),name)
956 return out
957
958 def checkReciprocalSymmetry(self,name0,name1,verbose=True):
959 """
960 test two coefficients for resciprocal symmetry
961
962 @param name0: name of the first coefficient
963 @type name0: C{str}
964 @param name1: name of the second coefficient
965 @type name1: C{str}
966 @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.
967 @type verbose: C{bool}
968 @return: True if coefficients C{name0} and C{name1} are resciprocally symmetric.
969 @rtype: C{bool}
970 """
971 B=self.getCoefficient(name0)
972 C=self.getCoefficient(name1)
973 verbose=verbose or self.__debug
974 out=True
975 if B.isEmpty() and not C.isEmpty():
976 if verbose: print "non-symmetric problem because %s is not present but %s is"%(name0,name1)
977 out=False
978 elif not B.isEmpty() and C.isEmpty():
979 if verbose: print "non-symmetric problem because %s is not present but %s is"%(name0,name1)
980 out=False
981 elif not B.isEmpty() and not C.isEmpty():
982 sB=B.getShape()
983 sC=C.getShape()
984 tol=(util.Lsup(B)+util.Lsup(C))*self.SMALL_TOLERANCE/2.
985 if len(sB) != len(sC):
986 if verbose: print "non-symmetric problem because ranks of %s (=%s) and %s (=%s) are different."%(name0,len(sB),name1,len(sC))
987 out=False
988 else:
989 if len(sB)==0:
990 if util.Lsup(B-C)>tol:
991 if verbose: print "non-symmetric problem because %s!=%s"%(name0,name1)
992 out=False
993 elif len(sB)==1:
994 if sB[0]==sC[0]:
995 for j in range(sB[0]):
996 if util.Lsup(B[j]-C[j])>tol:
997 if verbose: print "non-symmetric PDE because %s[%d]!=%s[%d]"%(name0,j,name1,j)
998 out=False
999 else:
1000 if verbose: print "non-symmetric problem because of inappropriate shapes %s and %s of coefficients %s and %s, respectively."%(sB,sC,name0,name1)
1001 elif len(sB)==3:
1002 if sB[0]==sC[1] and sB[1]==sC[2] and sB[2]==sC[0]:
1003 for i in range(sB[0]):
1004 for j in range(sB[1]):
1005 for k in range(sB[2]):
1006 if util.Lsup(B[i,j,k]-C[k,i,j])>tol:
1007 if verbose: print "non-symmetric problem because %s[%d,%d,%d]!=%s[%d,%d,%d]"%(name0,i,j,k,name1,k,i,j)
1008 out=False
1009 else:
1010 if verbose: print "non-symmetric problem because of inappropriate shapes %s and %s of coefficients %s and %s, respectively."%(sB,sC,name0,name1)
1011 else:
1012 raise ValueError,"Cannot check rank %s of %s and %s."%(len(sB),name0,name1)
1013 return out
1014
1015 def getCoefficient(self,name):
1016 """
1017 returns the value of the coefficient name
1018
1019 @param name: name of the coefficient requested.
1020 @type name: C{string}
1021 @return: the value of the coefficient name
1022 @rtype: L{Data<escript.Data>}
1023 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1024 """
1025 if self.hasCoefficient(name):
1026 return self.__COEFFICIENTS[name].getValue()
1027 else:
1028 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1029
1030 def hasCoefficient(self,name):
1031 """
1032 return True if name is the name of a coefficient
1033
1034 @param name: name of the coefficient enquired.
1035 @type name: C{string}
1036 @return: True if name is the name of a coefficient of the general PDE. Otherwise False.
1037 @rtype: C{bool}
1038 """
1039 return self.__COEFFICIENTS.has_key(name)
1040
1041 def createCoefficient(self, name):
1042 """
1043 create a L{Data<escript.Data>} object corresponding to coefficient name
1044
1045 @return: a coefficient name initialized to 0.
1046 @rtype: L{Data<escript.Data>}
1047 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1048 """
1049 if self.hasCoefficient(name):
1050 return escript.Data(0.,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))
1051 else:
1052 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1053
1054 def getFunctionSpaceForCoefficient(self,name):
1055 """
1056 return the L{FunctionSpace<escript.FunctionSpace>} to be used for coefficient name
1057
1058 @param name: name of the coefficient enquired.
1059 @type name: C{string}
1060 @return: the function space to be used for coefficient name
1061 @rtype: L{FunctionSpace<escript.FunctionSpace>}
1062 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1063 """
1064 if self.hasCoefficient(name):
1065 return self.__COEFFICIENTS[name].getFunctionSpace(self.getDomain())
1066 else:
1067 raise ValueError,"unknown coefficient %s requested"%name
1068
1069 def getShapeOfCoefficient(self,name):
1070 """
1071 return the shape of the coefficient name
1072
1073 @param name: name of the coefficient enquired.
1074 @type name: C{string}
1075 @return: the shape of the coefficient name
1076 @rtype: C{tuple} of C{int}
1077 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1078 """
1079 if self.hasCoefficient(name):
1080 return self.__COEFFICIENTS[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())
1081 else:
1082 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1083
1084 def resetAllCoefficients(self):
1085 """
1086 resets all coefficients to there default values.
1087 """
1088 for i in self.__COEFFICIENTS.iterkeys():
1089 self.__COEFFICIENTS[i].resetValue()
1090
1091 def alteredCoefficient(self,name):
1092 """
1093 announce that coefficient name has been changed
1094
1095 @param name: name of the coefficient enquired.
1096 @type name: C{string}
1097 @raise IllegalCoefficient: if name is not a coefficient of the PDE.
1098 @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.
1099 """
1100 if self.hasCoefficient(name):
1101 self.trace("Coefficient %s has been altered."%name)
1102 if not ((name=="q" or name=="r") and self.isUsingLumping()):
1103 if self.__COEFFICIENTS[name].isAlteringOperator(): self.invalidateOperator()
1104 if self.__COEFFICIENTS[name].isAlteringRightHandSide(): self.invalidateRightHandSide()
1105 else:
1106 raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1107
1108 def validSolution(self):
1109 """
1110 markes the solution as valid
1111 """
1112 self.__is_solution_valid=True
1113
1114 def invalidateSolution(self):
1115 """
1116 indicates the PDE has to be resolved if the solution is requested
1117 """
1118 self.trace("System will be resolved.")
1119 self.__is_solution_valid=False
1120
1121 def isSolutionValid(self):
1122 """
1123 returns True is the solution is still valid.
1124 """
1125 return self.__is_solution_valid
1126
1127 def validOperator(self):
1128 """
1129 markes the operator as valid
1130 """
1131 self.__is_operator_valid=True
1132
1133 def invalidateOperator(self):
1134 """
1135 indicates the operator has to be rebuilt next time it is used
1136 """
1137 self.trace("Operator will be rebuilt.")
1138 self.invalidateSolution()
1139 self.__is_operator_valid=False
1140
1141 def isOperatorValid(self):
1142 """
1143 returns True is the operator is still valid.
1144 """
1145 return self.__is_operator_valid
1146
1147 def validRightHandSide(self):
1148 """
1149 markes the right hand side as valid
1150 """
1151 self.__is_RHS_valid=True
1152
1153 def invalidateRightHandSide(self):
1154 """
1155 indicates the right hand side has to be rebuild next time it is used
1156 """
1157 if self.isRightHandSideValid(): self.trace("Right hand side has to be rebuilt.")
1158 self.invalidateSolution()
1159 self.__is_RHS_valid=False
1160
1161 def isRightHandSideValid(self):
1162 """
1163 returns True is the operator is still valid.
1164 """
1165 return self.__is_RHS_valid
1166
1167 def invalidateSystem(self):
1168 """
1169 annonced that everthing has to be rebuild:
1170 """
1171 if self.isRightHandSideValid(): self.trace("System has to be rebuilt.")
1172 self.invalidateSolution()
1173 self.invalidateOperator()
1174 self.invalidateRightHandSide()
1175
1176 def isSystemValid(self):
1177 """
1178 returns true if the system (including solution) is still vaild:
1179 """
1180 return self.isSolutionValid() and self.isOperatorValid() and self.isRightHandSideValid()
1181
1182 def initializeSystem(self):
1183 """
1184 resets the system
1185 """
1186 self.trace("New System has been created.")
1187 self.__operator=escript.Operator()
1188 self.__righthandside=escript.Data()
1189 self.__solution=escript.Data()
1190 self.invalidateSystem()
1191
1192 def getOperator(self):
1193 """
1194 returns the operator of the linear problem
1195
1196 @return: the operator of the problem
1197 """
1198 return self.getSystem()[0]
1199
1200 def getRightHandSide(self):
1201 """
1202 returns the right hand side of the linear problem
1203
1204 @return: the right hand side of the problem
1205 @rtype: L{Data<escript.Data>}
1206 """
1207 return self.getSystem()[1]
1208
1209 def createRightHandSide(self):
1210 """
1211 returns an instance of a new right hand side
1212 """
1213 self.trace("New right hand side is allocated.")
1214 if self.getNumEquations()>1:
1215 return escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)
1216 else:
1217 return escript.Data(0.,(),self.getFunctionSpaceForEquation(),True)
1218
1219 def createSolution(self):
1220 """
1221 returns an instance of a new solution
1222 """
1223 self.trace("New solution is allocated.")
1224 if self.getNumSolutions()>1:
1225 return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)
1226 else:
1227 return escript.Data(0.,(),self.getFunctionSpaceForSolution(),True)
1228
1229 def resetSolution(self):
1230 """
1231 sets the solution to zero
1232 """
1233 if self.__solution.isEmpty():
1234 self.__solution=self.createSolution()
1235 else:
1236 self.__solution.setToZero()
1237 self.trace("Solution is reset to zero.")
1238 def setSolution(self,u):
1239 """
1240 sets the solution assuming that makes the sytem valid.
1241 """
1242 self.__solution=u
1243 self.validSolution()
1244
1245 def getCurrentSolution(self):
1246 """
1247 returns the solution in its current state
1248 """
1249 if self.__solution.isEmpty(): self.__solution=self.createSolution()
1250 return self.__solution
1251
1252 def resetRightHandSide(self):
1253 """
1254 sets the right hand side to zero
1255 """
1256 if self.__righthandside.isEmpty():
1257 self.__righthandside=self.createRightHandSide()
1258 else:
1259 self.__righthandside.setToZero()
1260 self.trace("Right hand side is reset to zero.")
1261
1262 def getCurrentRightHandSide(self):
1263 """
1264 returns the right hand side in its current state
1265 """
1266 if self.__righthandside.isEmpty(): self.__righthandside=self.createRightHandSide()
1267 return self.__righthandside
1268
1269
1270 def resetOperator(self):
1271 """
1272 makes sure that the operator is instantiated and returns it initialized by zeros
1273 """
1274 if self.__operator.isEmpty():
1275 if self.isUsingLumping():
1276 self.__operator=self.createSolution()
1277 else:
1278 self.__operator=self.createOperator()
1279 else:
1280 if self.isUsingLumping():
1281 self.__operator.setToZero()
1282 else:
1283 self.__operator.resetValues()
1284 self.trace("Operator reset to zero")
1285
1286 def getCurrentOperator(self):
1287 """
1288 returns the operator in its current state
1289 """
1290 if self.__operator.isEmpty(): self.resetOperator()
1291 return self.__operator
1292
1293 def setValue(self,**coefficients):
1294 """
1295 sets new values to coefficients
1296
1297 @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1298 """
1299 # check if the coefficients are legal:
1300 for i in coefficients.iterkeys():
1301 if not self.hasCoefficient(i):
1302 raise IllegalCoefficient,"Attempt to set unknown coefficient %s"%i
1303 # if the number of unknowns or equations is still unknown we try to estimate them:
1304 if self.__numEquations==None or self.__numSolutions==None:
1305 for i,d in coefficients.iteritems():
1306 if hasattr(d,"shape"):
1307 s=d.shape
1308 elif hasattr(d,"getShape"):
1309 s=d.getShape()
1310 else:
1311 s=numarray.array(d).shape
1312 if s!=None:
1313 # get number of equations and number of unknowns:
1314 res=self.__COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(self.getDomain(),s)
1315 if res==None:
1316 raise IllegalCoefficientValue,"Illegal shape %s of coefficient %s"%(s,i)
1317 else:
1318 if self.__numEquations==None: self.__numEquations=res[0]
1319 if self.__numSolutions==None: self.__numSolutions=res[1]
1320 if self.__numEquations==None: raise UndefinedPDEError,"unidententified number of equations"
1321 if self.__numSolutions==None: raise UndefinedPDEError,"unidententified number of solutions"
1322 # now we check the shape of the coefficient if numEquations and numSolutions are set:
1323 for i,d in coefficients.iteritems():
1324 try:
1325 self.__COEFFICIENTS[i].setValue(self.getDomain(),
1326 self.getNumEquations(),self.getNumSolutions(),
1327 self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
1328 self.alteredCoefficient(i)
1329 except IllegalCoefficientFunctionSpace,m:
1330 # if the function space is wrong then we try the reduced version:
1331 i_red=i+"_reduced"
1332 if (not i_red in coefficients.keys()) and i_red in self.__COEFFICIENTS.keys():
1333 try:
1334 self.__COEFFICIENTS[i_red].setValue(self.getDomain(),
1335 self.getNumEquations(),self.getNumSolutions(),
1336 self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
1337 self.alteredCoefficient(i_red)
1338 except IllegalCoefficientValue,m:
1339 raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
1340 except IllegalCoefficientFunctionSpace,m:
1341 raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m))
1342 else:
1343 raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m))
1344 except IllegalCoefficientValue,m:
1345 raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
1346 self.__altered_coefficients=True
1347
1348 # ====================================================================================
1349 # methods that are typically overwritten when implementing a particular linear problem
1350 # ====================================================================================
1351 def getRequiredSystemType(self):
1352 """
1353 returns the system type which needs to be used by the current set up.
1354
1355 @note: Typically this method is overwritten when implementing a particular linear problem.
1356 """
1357 return 0
1358
1359 def createOperator(self):
1360 """
1361 returns an instance of a new operator
1362
1363 @note: This method is overwritten when implementing a particular linear problem.
1364 """
1365 return escript.Operator()
1366
1367 def checkSymmetry(self,verbose=True):
1368 """
1369 test the PDE for symmetry.
1370
1371 @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.
1372 @type verbose: C{bool}
1373 @return: True if the problem is symmetric.
1374 @rtype: C{bool}
1375 @note: Typically this method is overwritten when implementing a particular linear problem.
1376 """
1377 out=True
1378 return out
1379
1380 def getSolution(self,**options):
1381 """
1382 returns the solution of the problem
1383 @return: the solution
1384 @rtype: L{Data<escript.Data>}
1385
1386 @note: This method is overwritten when implementing a particular linear problem.
1387 """
1388 return self.getCurrentSolution()
1389
1390 def getSystem(self):
1391 """
1392 return the operator and right hand side of the PDE
1393
1394 @return: the discrete version of the PDE
1395 @rtype: C{tuple} of L{Operator,<escript.Operator>} and L{Data<escript.Data>}.
1396
1397 @note: This method is overwritten when implementing a particular linear problem.
1398 """
1399 return (self.getCurrentOperator(), self.getCurrentRightHandSide())
1400 #=====================
1401
1402 class LinearPDE(LinearProblem):
1403 """
1404 This class is used to define a general linear, steady, second order PDE
1405 for an unknown function M{u} on a given domain defined through a L{Domain<escript.Domain>} object.
1406
1407 For a single PDE with a solution with a single component the linear PDE is defined in the following form:
1408
1409 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)}
1410
1411
1412 where M{grad(F)} denotes the spatial derivative of M{F}. Einstein's summation convention,
1413 ie. summation over indexes appearing twice in a term of a sum is performed, is used.
1414 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
1415 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
1416 L{ReducedFunction<escript.ReducedFunction>}. It is also allowd to use objects that can be converted into
1417 such L{Data<escript.Data>} objects. M{A} and M{A_reduced} are rank two, M{B}, M{C}, M{X},
1418 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.
1419
1420 The following natural boundary conditions are considered:
1421
1422 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}
1423
1424 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>}.
1425
1426
1427 Constraints for the solution prescribing the value of the solution at certain locations in the domain. They have the form
1428
1429 M{u=r} where M{q>0}
1430
1431 M{r} and M{q} are each scalar where M{q} is the characteristic function defining where the constraint is applied.
1432 The constraints override any other condition set by the PDE or the boundary condition.
1433
1434 The PDE is symmetrical if
1435
1436 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]}
1437
1438 For a system of PDEs and a solution with several components the PDE has the form
1439
1440 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] }
1441
1442 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.
1443 The natural boundary conditions take the form:
1444
1445 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]}
1446
1447
1448 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>}.
1449
1450 Constraints take the form
1451
1452 M{u[i]=r[i]} where M{q[i]>0}
1453
1454 M{r} and M{q} are each rank one. Notice that at some locations not necessarily all components must have a constraint.
1455
1456 The system of PDEs is symmetrical if
1457
1458 - M{A[i,j,k,l]=A[k,l,i,j]}
1459 - M{A_reduced[i,j,k,l]=A_reduced[k,l,i,j]}
1460 - M{B[i,j,k]=C[k,i,j]}
1461 - M{B_reduced[i,j,k]=C_reduced[k,i,j]}
1462 - M{D[i,k]=D[i,k]}
1463 - M{D_reduced[i,k]=D_reduced[i,k]}
1464 - M{d[i,k]=d[k,i]}
1465 - M{d_reduced[i,k]=d_reduced[k,i]}
1466
1467 L{LinearPDE} also supports solution discontinuities over a contact region in the domain. To specify the conditions across the
1468 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
1469 defined as
1470
1471 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]}
1472
1473 For the case of single solution component and single PDE M{J} is defined
1474
1475 M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[j]+(B[i]+B_reduced[i])*u-X[i]-X_reduced[i]}
1476
1477 In the context of discontinuities M{n} denotes the normal on the discontinuity pointing from side 0 towards side 1
1478 calculated from L{getNormal<escript.FunctionSpace.getNormal>} of L{FunctionOnContactZero<escript.FunctionOnContactZero>}. For a system of PDEs
1479 the contact condition takes the form
1480
1481 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]}
1482
1483 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
1484 of the solution at side 1 and at side 0, denotes the jump of M{u} across discontinuity along the normal calculated by
1485 L{jump<util.jump>}.
1486 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>}.
1487 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>}.
1488 In case of a single PDE and a single component solution the contact condition takes the form
1489
1490 M{n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)}
1491
1492 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>}
1493
1494 typical usage:
1495
1496 p=LinearPDE(dom)
1497 p.setValue(A=kronecker(dom),D=1,Y=0.5)
1498 u=p.getSolution()
1499
1500 """
1501
1502
1503 def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):
1504 """
1505 initializes a new linear PDE
1506
1507 @param domain: domain of the PDE
1508 @type domain: L{Domain<escript.Domain>}
1509 @param numEquations: number of equations. If numEquations==None the number of equations
1510 is extracted from the PDE coefficients.
1511 @param numSolutions: number of solution components. If numSolutions==None the number of solution components
1512 is extracted from the PDE coefficients.
1513 @param debug: if True debug informations are printed.
1514
1515 """
1516 super(LinearPDE, self).__init__(domain,numEquations,numSolutions,debug)
1517 #
1518 # the coefficients of the PDE:
1519 #
1520 self.introduceCoefficients(
1521 A=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
1522 B=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1523 C=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
1524 D=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1525 X=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
1526 Y=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1527 d=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1528 y=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1529 d_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1530 y_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1531 A_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
1532 B_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1533 C_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
1534 D_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1535 X_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
1536 Y_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1537 d_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1538 y_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1539 d_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1540 y_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1541 r=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.RIGHTHANDSIDE),
1542 q=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.BOTH) )
1543
1544 def __str__(self):
1545 """
1546 returns string representation of the PDE
1547
1548 @return: a simple representation of the PDE
1549 @rtype: C{str}
1550 """
1551 return "<LinearPDE %d>"%id(self)
1552
1553 def getRequiredSystemType(self):
1554 """
1555 returns the system type which needs to be used by the current set up.
1556 """
1557 return self.getDomain().getSystemMatrixTypeId(self.getSolverMethod()[0],self.getSolverMethod()[1],self.getSolverPackage(),self.isSymmetric())
1558
1559 def checkSymmetry(self,verbose=True):
1560 """
1561 test the PDE for symmetry.
1562
1563 @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.
1564 @type verbose: C{bool}
1565 @return: True if the PDE is symmetric.
1566 @rtype: L{bool}
1567 @note: This is a very expensive operation. It should be used for degugging only! The symmetry flag is not altered.
1568 """
1569 out=True
1570 out=out and self.checkSymmetricTensor("A", verbose)
1571 out=out and self.checkSymmetricTensor("A_reduced", verbose)
1572 out=out and self.checkReciprocalSymmetry("B","C", verbose)
1573 out=out and self.checkReciprocalSymmetry("B_reduced","C_reduced", verbose)
1574 out=out and self.checkSymmetricTensor("D", verbose)
1575 out=out and self.checkSymmetricTensor("D_reduced", verbose)
1576 out=out and self.checkSymmetricTensor("d", verbose)
1577 out=out and self.checkSymmetricTensor("d_reduced", verbose)
1578 out=out and self.checkSymmetricTensor("d_contact", verbose)
1579 out=out and self.checkSymmetricTensor("d_contact_reduced", verbose)
1580 return out
1581
1582 def createOperator(self):
1583 """
1584 returns an instance of a new operator
1585 """
1586 self.trace("New operator is allocated.")
1587 return self.getDomain().newOperator( \
1588 self.getNumEquations(), \
1589 self.getFunctionSpaceForEquation(), \
1590 self.getNumSolutions(), \
1591 self.getFunctionSpaceForSolution(), \
1592 self.getSystemType())
1593
1594 def getSolution(self,**options):
1595 """
1596 returns the solution of the PDE
1597
1598 @return: the solution
1599 @rtype: L{Data<escript.Data>}
1600 @param options: solver options
1601 @keyword verbose: True to get some information during PDE solution
1602 @type verbose: C{bool}
1603 @keyword reordering: reordering scheme to be used during elimination. Allowed values are
1604 L{NO_REORDERING}, L{MINIMUM_FILL_IN}, L{NESTED_DISSECTION}
1605 @keyword iter_max: maximum number of iteration steps allowed.
1606 @keyword drop_tolerance: threshold for drupping in L{ILUT}
1607 @keyword drop_storage: maximum of allowed memory in L{ILUT}
1608 @keyword truncation: maximum number of residuals in L{GMRES}
1609 @keyword restart: restart cycle length in L{GMRES}
1610 """
1611 if not self.isSolutionValid():
1612 mat,f=self.getSystem()
1613 if self.isUsingLumping():
1614 self.setSolution(f*1/mat)
1615 else:
1616 options[self.TOLERANCE_KEY]=self.getTolerance()
1617 options[self.METHOD_KEY]=self.getSolverMethod()[0]
1618 options[self.PRECONDITIONER_KEY]=self.getSolverMethod()[1]
1619 options[self.PACKAGE_KEY]=self.getSolverPackage()
1620 options[self.SYMMETRY_KEY]=self.isSymmetric()
1621 self.trace("PDE is resolved.")
1622 self.trace("solver options: %s"%str(options))
1623 self.setSolution(mat.solve(f,options))
1624 return self.getCurrentSolution()
1625
1626 def getSystem(self):
1627 """
1628 return the operator and right hand side of the PDE
1629
1630 @return: the discrete version of the PDE
1631 @rtype: C{tuple} of L{Operator,<escript.Operator>} and L{Data<escript.Data>}.
1632 """
1633 if not self.isOperatorValid() or not self.isRightHandSideValid():
1634 if self.isUsingLumping():
1635 if not self.isOperatorValid():
1636 if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():
1637 raise TypeError,"Lumped matrix requires same order for equations and unknowns"
1638 if not self.getCoefficient("A").isEmpty():
1639 raise ValueError,"coefficient A in lumped matrix may not be present."
1640 if not self.getCoefficient("B").isEmpty():
1641 raise ValueError,"coefficient B in lumped matrix may not be present."
1642 if not self.getCoefficient("C").isEmpty():
1643 raise ValueError,"coefficient C in lumped matrix may not be present."
1644 if not self.getCoefficient("d_contact").isEmpty():
1645 raise ValueError,"coefficient d_contact in lumped matrix may not be present."
1646 if not self.getCoefficient("A_reduced").isEmpty():
1647 raise ValueError,"coefficient A_reduced in lumped matrix may not be present."
1648 if not self.getCoefficient("B_reduced").isEmpty():
1649 raise ValueError,"coefficient B_reduced in lumped matrix may not be present."
1650 if not self.getCoefficient("C_reduced").isEmpty():
1651 raise ValueError,"coefficient C_reduced in lumped matrix may not be present."
1652 if not self.getCoefficient("d_contact_reduced").isEmpty():
1653 raise ValueError,"coefficient d_contact_reduced in lumped matrix may not be present."
1654 D=self.getCoefficient("D")
1655 d=self.getCoefficient("d")
1656 D_reduced=self.getCoefficient("D_reduced")
1657 d_reduced=self.getCoefficient("d_reduced")
1658 if not D.isEmpty():
1659 if self.getNumSolutions()>1:
1660 D_times_e=util.matrix_mult(D,numarray.ones((self.getNumSolutions(),)))
1661 else:
1662 D_times_e=D
1663 else:
1664 D_times_e=escript.Data()
1665 if not d.isEmpty():
1666 if self.getNumSolutions()>1:
1667 d_times_e=util.matrix_mult(d,numarray.ones((self.getNumSolutions(),)))
1668 else:
1669 d_times_e=d
1670 else:
1671 d_times_e=escript.Data()
1672
1673 if not D_reduced.isEmpty():
1674 if self.getNumSolutions()>1:
1675 D_reduced_times_e=util.matrix_mult(D_reduced,numarray.ones((self.getNumSolutions(),)))
1676 else:
1677 D_reduced_times_e=D_reduced
1678 else:
1679 D_reduced_times_e=escript.Data()
1680 if not d_reduced.isEmpty():
1681 if self.getNumSolutions()>1:
1682 d_reduced_times_e=util.matrix_mult(d_reduced,numarray.ones((self.getNumSolutions(),)))
1683 else:
1684 d_reduced_times_e=d_reduced
1685 else:
1686 d_reduced_times_e=escript.Data()
1687
1688 self.resetOperator()
1689 operator=self.getCurrentOperator()
1690 if False and hasattr(self.getDomain(), "addPDEToLumpedSystem") :
1691 self.getDomain().addPDEToLumpedSystem(operator, D_times_e, d_times_e)
1692 self.getDomain().addPDEToLumpedSystem(operator, D_reduced_times_e, d_reduced_times_e)
1693 else:
1694 self.getDomain().addPDEToRHS(operator, \
1695 escript.Data(), \
1696 D_times_e, \
1697 d_times_e,\
1698 escript.Data())
1699 self.getDomain().addPDEToRHS(operator, \
1700 escript.Data(), \
1701 D_reduced_times_e, \
1702 d_reduced_times_e,\
1703 escript.Data())
1704 self.trace("New lumped operator has been built.")
1705 if not self.isRightHandSideValid():
1706 self.resetRightHandSide()
1707 righthandside=self.getCurrentRightHandSide()
1708 self.getDomain().addPDEToRHS(righthandside, \
1709 self.getCoefficient("X"), \
1710 self.getCoefficient("Y"),\
1711 self.getCoefficient("y"),\
1712 self.getCoefficient("y_contact"))
1713 self.getDomain().addPDEToRHS(righthandside, \
1714 self.getCoefficient("X_reduced"), \
1715 self.getCoefficient("Y_reduced"),\
1716 self.getCoefficient("y_reduced"),\
1717 self.getCoefficient("y_contact_reduced"))
1718 self.trace("New right hand side as been built.")
1719 self.validRightHandSide()
1720 self.insertConstraint(rhs_only=False)
1721 self.validOperator()
1722 else:
1723 if not self.isOperatorValid() and not self.isRightHandSideValid():
1724 self.resetRightHandSide()
1725 righthandside=self.getCurrentRightHandSide()
1726 self.resetOperator()
1727 operator=self.getCurrentOperator()
1728 self.getDomain().addPDEToSystem(operator,righthandside, \
1729 self.getCoefficient("A"), \
1730 self.getCoefficient("B"), \
1731 self.getCoefficient("C"), \
1732 self.getCoefficient("D"), \
1733 self.getCoefficient("X"), \
1734 self.getCoefficient("Y"), \
1735 self.getCoefficient("d"), \
1736 self.getCoefficient("y"), \
1737 self.getCoefficient("d_contact"), \
1738 self.getCoefficient("y_contact"))
1739 self.getDomain().addPDEToSystem(operator,righthandside, \
1740 self.getCoefficient("A_reduced"), \
1741 self.getCoefficient("B_reduced"), \
1742 self.getCoefficient("C_reduced"), \
1743 self.getCoefficient("D_reduced"), \
1744 self.getCoefficient("X_reduced"), \
1745 self.getCoefficient("Y_reduced"), \
1746 self.getCoefficient("d_reduced"), \
1747 self.getCoefficient("y_reduced"), \
1748 self.getCoefficient("d_contact_reduced"), \
1749 self.getCoefficient("y_contact_reduced"))
1750 self.insertConstraint(rhs_only=False)
1751 self.trace("New system has been built.")
1752 self.validOperator()
1753 self.validRightHandSide()
1754 elif not self.isRightHandSideValid():
1755 self.resetRightHandSide()
1756 righthandside=self.getCurrentRightHandSide()
1757 self.getDomain().addPDEToRHS(righthandside,
1758 self.getCoefficient("X"), \
1759 self.getCoefficient("Y"),\
1760 self.getCoefficient("y"),\
1761 self.getCoefficient("y_contact"))
1762 self.getDomain().addPDEToRHS(righthandside,
1763 self.getCoefficient("X_reduced"), \
1764 self.getCoefficient("Y_reduced"),\
1765 self.getCoefficient("y_reduced"),\
1766 self.getCoefficient("y_contact_reduced"))
1767 self.insertConstraint(rhs_only=True)
1768 self.trace("New right hand side has been built.")
1769 self.validRightHandSide()
1770 elif not self.isOperatorValid():
1771 self.resetOperator()
1772 operator=self.getCurrentOperator()
1773 self.getDomain().addPDEToSystem(operator,escript.Data(), \
1774 self.getCoefficient("A"), \
1775 self.getCoefficient("B"), \
1776 self.getCoefficient("C"), \
1777 self.getCoefficient("D"), \
1778 escript.Data(), \
1779 escript.Data(), \
1780 self.getCoefficient("d"), \
1781 escript.Data(),\
1782 self.getCoefficient("d_contact"), \
1783 escript.Data())
1784 self.getDomain().addPDEToSystem(operator,escript.Data(), \
1785 self.getCoefficient("A_reduced"), \
1786 self.getCoefficient("B_reduced"), \
1787 self.getCoefficient("C_reduced"), \
1788 self.getCoefficient("D_reduced"), \
1789 escript.Data(), \
1790 escript.Data(), \
1791 self.getCoefficient("d_reduced"), \
1792 escript.Data(),\
1793 self.getCoefficient("d_contact_reduced"), \
1794 escript.Data())
1795 self.insertConstraint(rhs_only=False)
1796 self.trace("New operator has been built.")
1797 self.validOperator()
1798 return (self.getCurrentOperator(), self.getCurrentRightHandSide())
1799
1800 def insertConstraint(self, rhs_only=False):
1801 """
1802 applies the constraints defined by q and r to the PDE
1803
1804 @param rhs_only: if True only the right hand side is altered by the constraint.
1805 @type rhs_only: C{bool}
1806 """
1807 q=self.getCoefficient("q")
1808 r=self.getCoefficient("r")
1809 righthandside=self.getCurrentRightHandSide()
1810 operator=self.getCurrentOperator()
1811
1812 if not q.isEmpty():
1813 if r.isEmpty():
1814 r_s=self.createSolution()
1815 else:
1816 r_s=r
1817 if not rhs_only and not operator.isEmpty():
1818 if self.isUsingLumping():
1819 operator.copyWithMask(escript.Data(1.,q.getShape(),q.getFunctionSpace()),q)
1820 else:
1821 row_q=escript.Data(q,self.getFunctionSpaceForEquation())
1822 col_q=escript.Data(q,self.getFunctionSpaceForSolution())
1823 u=self.createSolution()
1824 u.copyWithMask(r_s,col_q)
1825 righthandside-=operator*u
1826 operator.nullifyRowsAndCols(row_q,col_q,1.)
1827 righthandside.copyWithMask(r_s,q)
1828
1829 def setValue(self,**coefficients):
1830 """
1831 sets new values to coefficients
1832
1833 @param coefficients: new values assigned to coefficients
1834 @keyword A: value for coefficient A.
1835 @type A: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1836 @keyword A_reduced: value for coefficient A_reduced.
1837 @type A_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1838 @keyword B: value for coefficient B
1839 @type B: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1840 @keyword B_reduced: value for coefficient B_reduced
1841 @type B_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1842 @keyword C: value for coefficient C
1843 @type C: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1844 @keyword C_reduced: value for coefficient C_reduced
1845 @type C_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1846 @keyword D: value for coefficient D
1847 @type D: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
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{ReducedFunction<escript.ReducedFunction>}.
1850 @keyword X: value for coefficient X
1851 @type X: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1852 @keyword X_reduced: value for coefficient X_reduced
1853 @type X_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
1854 @keyword Y: value for coefficient Y
1855 @type Y: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
1856 @keyword Y_reduced: value for coefficient Y_reduced
1857 @type Y_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.Function>}.
1858 @keyword d: value for coefficient d
1859 @type d: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1860 @keyword d_reduced: value for coefficient d_reduced
1861 @type d_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
1862 @keyword y: value for coefficient y
1863 @type y: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
1864 @keyword d_contact: value for coefficient d_contact
1865 @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>}.
1866 @keyword d_contact_reduced: value for coefficient d_contact_reduced
1867 @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>}.
1868 @keyword y_contact: value for coefficient y_contact
1869 @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>}.
1870 @keyword y_contact_reduced: value for coefficient y_contact_reduced
1871 @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>}.
1872 @keyword r: values prescribed to the solution at the locations of constraints
1873 @type r: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1874 depending of reduced order is used for the solution.
1875 @keyword q: mask for location of constraints
1876 @type q: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
1877 depending of reduced order is used for the representation of the equation.
1878 @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1879 """
1880 super(LinearPDE,self).setValue(**coefficients)
1881 # check if the systrem is inhomogeneous:
1882 if len(coefficients)>0 and not self.isUsingLumping():
1883 q=self.getCoefficient("q")
1884 r=self.getCoefficient("r")
1885 if not q.isEmpty() and not r.isEmpty():
1886 if util.Lsup(q*r)>0.:
1887 self.trace("Inhomogeneous constraint detected.")
1888 self.invalidateSystem()
1889
1890
1891 def getResidual(self,u=None):
1892 """
1893 return the residual of u or the current solution if u is not present.
1894
1895 @param u: argument in the residual calculation. It must be representable in C{elf.getFunctionSpaceForSolution()}. If u is not present or equals L{None}
1896 the current solution is used.
1897 @type u: L{Data<escript.Data>} or None
1898 @return: residual of u
1899 @rtype: L{Data<escript.Data>}
1900 """
1901 if u==None:
1902 return self.getOperator()*self.getSolution()-self.getRightHandSide()
1903 else:
1904 return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())-self.getRightHandSide()
1905
1906 def getFlux(self,u=None):
1907 """
1908 returns the flux M{J} for a given M{u}
1909
1910 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]}
1911
1912 or
1913
1914 M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[l]+(B[j]+B_reduced[j])u-X[j]-X_reduced[j]}
1915
1916 @param u: argument in the flux. If u is not present or equals L{None} the current solution is used.
1917 @type u: L{Data<escript.Data>} or None
1918 @return: flux
1919 @rtype: L{Data<escript.Data>}
1920 """
1921 if u==None: u=self.getSolution()
1922 return util.tensormult(self.getCoefficient("A"),util.grad(u,Funtion(self.getDomain))) \
1923 +util.matrixmult(self.getCoefficient("B"),u) \
1924 -util.self.getCoefficient("X") \
1925 +util.tensormult(self.getCoefficient("A_reduced"),util.grad(u,ReducedFuntion(self.getDomain))) \
1926 +util.matrixmult(self.getCoefficient("B_reduced"),u) \
1927 -util.self.getCoefficient("X_reduced")
1928
1929
1930 class Poisson(LinearPDE):
1931 """
1932 Class to define a Poisson equation problem, which is genear L{LinearPDE} of the form
1933
1934 M{-grad(grad(u)[j])[j] = f}
1935
1936 with natural boundary conditons
1937
1938 M{n[j]*grad(u)[j] = 0 }
1939
1940 and constraints:
1941
1942 M{u=0} where M{q>0}
1943
1944 """
1945
1946 def __init__(self,domain,debug=False):
1947 """
1948 initializes a new Poisson equation
1949
1950 @param domain: domain of the PDE
1951 @type domain: L{Domain<escript.Domain>}
1952 @param debug: if True debug informations are printed.
1953
1954 """
1955 super(Poisson, self).__init__(domain,1,1,debug)
1956 self.introduceCoefficients(
1957 f=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1958 f_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE))
1959 self.setSymmetryOn()
1960
1961 def setValue(self,**coefficients):
1962 """
1963 sets new values to coefficients
1964
1965 @param coefficients: new values assigned to coefficients
1966 @keyword f: value for right hand side M{f}
1967 @type f: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
1968 @keyword q: mask for location of constraints
1969 @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>}
1970 depending of reduced order is used for the representation of the equation.
1971 @raise IllegalCoefficient: if an unknown coefficient keyword is used.
1972 """
1973 super(Poisson, self).setValue(**coefficients)
1974
1975
1976 def getCoefficient(self,name):
1977 """
1978 return the value of the coefficient name of the general PDE
1979 @param name: name of the coefficient requested.
1980 @type name: C{string}
1981 @return: the value of the coefficient name
1982 @rtype: L{Data<escript.Data>}
1983 @raise IllegalCoefficient: invalid coefficient name
1984 @note: This method is called by the assembling routine to map the Possion equation onto the general PDE.
1985 """
1986 if name == "A" :
1987 return escript.Data(util.kronecker(self.getDim()),escript.Function(self.getDomain()))
1988 elif name == "Y" :
1989 return self.getCoefficient("f")
1990 elif name == "Y_reduced" :
1991 return self.getCoefficient("f_reduced")
1992 else:
1993 return super(Poisson, self).getCoefficient(name)
1994
1995 class Helmholtz(LinearPDE):
1996 """
1997 Class to define a Helmhotz equation problem, which is genear L{LinearPDE} of the form
1998
1999 M{S{omega}*u - grad(k*grad(u)[j])[j] = f}
2000
2001 with natural boundary conditons
2002
2003 M{k*n[j]*grad(u)[j] = g- S{alpha}u }
2004
2005 and constraints:
2006
2007 M{u=r} where M{q>0}
2008
2009 """
2010
2011 def __init__(self,domain,debug=False):
2012 """
2013 initializes a new Poisson equation
2014
2015 @param domain: domain of the PDE
2016 @type domain: L{Domain<escript.Domain>}
2017 @param debug: if True debug informations are printed.
2018
2019 """
2020 super(Helmholtz, self).__init__(domain,1,1,debug)
2021 self.introduceCoefficients(
2022 omega=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.OPERATOR),
2023 k=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.OPERATOR),
2024 f=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2025 f_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2026 alpha=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.OPERATOR),
2027 g=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2028 g_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE))
2029 self.setSymmetryOn()
2030
2031 def setValue(self,**coefficients):
2032 """
2033 sets new values to coefficients
2034
2035 @param coefficients: new values assigned to coefficients
2036 @keyword omega: value for coefficient M{S{omega}}
2037 @type omega: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2038 @keyword k: value for coefficeint M{k}
2039 @type k: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2040 @keyword f: value for right hand side M{f}
2041 @type f: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2042 @keyword alpha: value for right hand side M{S{alpha}}
2043 @type alpha: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
2044 @keyword g: value for right hand side M{g}
2045 @type g: any type that can be casted to L{Scalar<escript.Scalar>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
2046 @keyword r: prescribed values M{r} for the solution in constraints.
2047 @type r: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2048 depending of reduced order is used for the representation of the equation.
2049 @keyword q: mask for location of constraints
2050 @type q: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2051 depending of reduced order is used for the representation of the equation.
2052 @raise IllegalCoefficient: if an unknown coefficient keyword is used.
2053 """
2054 super(Helmholtz, self).setValue(**coefficients)
2055
2056 def getCoefficient(self,name):
2057 """
2058 return the value of the coefficient name of the general PDE
2059
2060 @param name: name of the coefficient requested.
2061 @type name: C{string}
2062 @return: the value of the coefficient name
2063 @rtype: L{Data<escript.Data>}
2064 @raise IllegalCoefficient: invalid name
2065 """
2066 if name == "A" :
2067 return escript.Data(numarray.identity(self.getDim()),escript.Function(self.getDomain()))*self.getCoefficient("k")
2068 elif name == "D" :
2069 return self.getCoefficient("omega")
2070 elif name == "Y" :
2071 return self.getCoefficient("f")
2072 elif name == "d" :
2073 return self.getCoefficient("alpha")
2074 elif name == "y" :
2075 return self.getCoefficient("g")
2076 elif name == "Y_reduced" :
2077 return self.getCoefficient("f_reduced")
2078 elif name == "y_reduced" :
2079 return self.getCoefficient("g_reduced")
2080 else:
2081 return super(Helmholtz, self).getCoefficient(name)
2082
2083 class LameEquation(LinearPDE):
2084 """
2085 Class to define a Lame equation problem:
2086
2087 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] }
2088
2089 with natural boundary conditons:
2090
2091 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] }
2092
2093 and constraints:
2094
2095 M{u[i]=r[i]} where M{q[i]>0}
2096
2097 """
2098
2099 def __init__(self,domain,debug=False):
2100 super(LameEquation, self).__init__(domain,\
2101 domain.getDim(),domain.getDim(),debug)
2102 self.introduceCoefficients(lame_lambda=PDECoef(PDECoef.INTERIOR,(),PDECoef.OPERATOR),
2103 lame_mu=PDECoef(PDECoef.INTERIOR,(),PDECoef.OPERATOR),
2104 F=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2105 sigma=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
2106 f=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE))
2107 self.setSymmetryOn()
2108
2109 def setValues(self,**coefficients):
2110 """
2111 sets new values to coefficients
2112
2113 @param coefficients: new values assigned to coefficients
2114 @keyword lame_mu: value for coefficient M{S{mu}}
2115 @type lame_mu: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2116 @keyword lame_lambda: value for coefficient M{S{lambda}}
2117 @type lame_lambda: any type that can be casted to L{Scalar<escript.Scalar>} object on L{Function<escript.Function>}.
2118 @keyword F: value for internal force M{F}
2119 @type F: any type that can be casted to L{Vector<escript.Vector>} object on L{Function<escript.Function>}
2120 @keyword sigma: value for initial stress M{S{sigma}}
2121 @type sigma: any type that can be casted to L{Tensor<escript.Tensor>} object on L{Function<escript.Function>}
2122 @keyword f: value for extrenal force M{f}
2123 @type f: any type that can be casted to L{Vector<escript.Vector>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}
2124 @keyword r: prescribed values M{r} for the solution in constraints.
2125 @type r: any type that can be casted to L{Vector<escript.Vector>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2126 depending of reduced order is used for the representation of the equation.
2127 @keyword q: mask for location of constraints
2128 @type q: any type that can be casted to L{Vector<escript.Vector>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2129 depending of reduced order is used for the representation of the equation.
2130 @raise IllegalCoefficient: if an unknown coefficient keyword is used.
2131 """
2132 super(LameEquation, self).setValues(**coefficients)
2133
2134 def getCoefficient(self,name):
2135 """
2136 return the value of the coefficient name of the general PDE
2137
2138 @param name: name of the coefficient requested.
2139 @type name: C{string}
2140 @return: the value of the coefficient name
2141 @rtype: L{Data<escript.Data>}
2142 @raise IllegalCoefficient: invalid coefficient name
2143 """
2144 if name == "A" :
2145 out =self.createCoefficient("A")
2146 for i in range(self.getDim()):
2147 for j in range(self.getDim()):
2148 out[i,i,j,j] += self.getCoefficient("lame_lambda")
2149 out[i,j,j,i] += self.getCoefficient("lame_mu")
2150 out[i,j,i,j] += self.getCoefficient("lame_mu")
2151 return out
2152 elif name == "X" :
2153 return self.getCoefficient("sigma")
2154 elif name == "Y" :
2155 return self.getCoefficient("F")
2156 elif name == "y" :
2157 return self.getCoefficient("f")
2158 else:
2159 return super(LameEquation, self).getCoefficient(name)
2160
2161 def LinearSinglePDE(domain,debug=False):
2162 """
2163 defines a single linear PDEs
2164
2165 @param domain: domain of the PDE
2166 @type domain: L{Domain<escript.Domain>}
2167 @param debug: if True debug informations are printed.
2168 @rtype: L{LinearPDE}
2169 """
2170 return LinearPDE(domain,numEquations=1,numSolutions=1,debug=debug)
2171
2172 def LinearPDESystem(domain,debug=False):
2173 """
2174 defines a system of linear PDEs
2175
2176 @param domain: domain of the PDE
2177 @type domain: L{Domain<escript.Domain>}
2178 @param debug: if True debug informations are printed.
2179 @rtype: L{LinearPDE}
2180 """
2181 return LinearPDE(domain,numEquations=domain.getDim(),numSolutions=domain.getDim(),debug=debug)
2182
2183 class TransportPDE(LinearProblem):
2184 """
2185 This class is used to define a transport problem given by a general linear, time dependend, second order PDE
2186 for an unknown, non-negative, time-dependend function M{u} on a given domain defined through a L{Domain<escript.Domain>} object.
2187
2188 For a single equation with a solution with a single component the transport problem is defined in the following form:
2189
2190 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)}
2191
2192
2193 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,
2194 ie. summation over indexes appearing twice in a term of a sum is performed, is used.
2195 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
2196 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
2197 L{ReducedFunction<escript.ReducedFunction>}. It is also allowed to use objects that can be converted into
2198 such L{Data<escript.Data>} objects. M{M} and M{M_reduced} are scalar, M{A} and M{A_reduced} are rank two,
2199 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.
2200
2201 The following natural boundary conditions are considered:
2202
2203 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}
2204
2205 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>}.
2206
2207 Constraints for the solution prescribing the value of the solution at certain locations in the domain. They have the form
2208
2209 M{u_t=r} where M{q>0}
2210
2211 M{r} and M{q} are each scalar where M{q} is the characteristic function defining where the constraint is applied.
2212 The constraints override any other condition set by the transport problem or the boundary condition.
2213
2214 The transport problem is symmetrical if
2215
2216 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]}
2217
2218 For a system and a solution with several components the transport problem has the form
2219
2220 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] }
2221
2222 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:
2223
2224 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}
2225
2226 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>}.
2227
2228 Constraints take the form
2229
2230 M{u[i]_t=r[i]} where M{q[i]>0}
2231
2232 M{r} and M{q} are each rank one. Notice that at some locations not necessarily all components must have a constraint.
2233
2234 The transport problem is symmetrical if
2235
2236 - M{M[i,k]=M[i,k]}
2237 - M{M_reduced[i,k]=M_reduced[i,k]}
2238 - M{A[i,j,k,l]=A[k,l,i,j]}
2239 - M{A_reduced[i,j,k,l]=A_reduced[k,l,i,j]}
2240 - M{B[i,j,k]=C[k,i,j]}
2241 - M{B_reduced[i,j,k]=C_reduced[k,i,j]}
2242 - M{D[i,k]=D[i,k]}
2243 - M{D_reduced[i,k]=D_reduced[i,k]}
2244 - M{m[i,k]=m[k,i]}
2245 - M{m_reduced[i,k]=m_reduced[k,i]}
2246 - M{d[i,k]=d[k,i]}
2247 - M{d_reduced[i,k]=d_reduced[k,i]}
2248
2249 L{TransportPDE} also supports solution discontinuities over a contact region in the domain. To specify the conditions across the
2250 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
2251 defined as
2252
2253 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]}
2254
2255 For the case of single solution component and single PDE M{J} is defined
2256
2257 M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[j]+(B[i]+B_reduced[i])*u+X[i]+X_reduced[i]}
2258
2259 In the context of discontinuities M{n} denotes the normal on the discontinuity pointing from side 0 towards side 1
2260 calculated from L{getNormal<escript.FunctionSpace.getNormal>} of L{FunctionOnContactZero<escript.FunctionOnContactZero>}. For a system of transport problems the contact condition takes the form
2261
2262 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]}
2263
2264 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
2265 of the solution at side 1 and at side 0, denotes the jump of M{u} across discontinuity along the normal calculated by
2266 L{jump<util.jump>}.
2267 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>}.
2268 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>}.
2269 In case of a single PDE and a single component solution the contact condition takes the form
2270
2271 M{n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)}
2272
2273 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>}
2274
2275 typical usage:
2276
2277 p=TransportPDE(dom)
2278 p.setValue(M=1.,C=[-1.,0.])
2279 p.setInitialSolution(u=exp(-length(dom.getX()-[0.1,0.1])**2)
2280 t=0
2281 dt=0.1
2282 while (t<1.):
2283 u=p.solve(dt)
2284
2285 """
2286 def __init__(self,domain,numEquations=None,numSolutions=None, theta=0.5,debug=False):
2287 """
2288 initializes a linear problem
2289
2290 @param domain: domain of the PDE
2291 @type domain: L{Domain<escript.Domain>}
2292 @param numEquations: number of equations. If numEquations==None the number of equations
2293 is extracted from the coefficients.
2294 @param numSolutions: number of solution components. If numSolutions==None the number of solution components
2295 is extracted from the coefficients.
2296 @param debug: if True debug informations are printed.
2297 @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>}.
2298
2299 """
2300 if theta<0 or theta>1:
2301 raise ValueError,"theta (=%s) needs to be between 0 and 1 (inclusive)."%theta
2302 super(TransportPDE, self).__init__(domain,numEquations,numSolutions,debug)
2303
2304 self.__theta=theta
2305 #
2306 # the coefficients of the transport problem
2307 #
2308 self.introduceCoefficients(
2309 M=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2310 A=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2311 B=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2312 C=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2313 D=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2314 X=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
2315 Y=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2316 m=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2317 d=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2318 y=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2319 d_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2320 y_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2321 M_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2322 A_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2323 B_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2324 C_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2325 D_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2326 X_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
2327 Y_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2328 m_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2329 d_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2330 y_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2331 d_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2332 y_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2333 r=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.RIGHTHANDSIDE),
2334 q=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.BOTH) )
2335
2336 def __str__(self):
2337 """
2338 returns string representation of the transport problem
2339
2340 @return: a simple representation of the transport problem
2341 @rtype: C{str}
2342 """
2343 return "<TransportPDE %d>"%id(self)
2344
2345 def getTheta(self):
2346 """
2347 returns the time stepping control parameter
2348 @rtype: float
2349 """
2350 return self.__theta
2351
2352
2353 def checkSymmetry(self,verbose=True):
2354 """
2355 test the transport problem for symmetry.
2356
2357 @param verbose: if equal to True or not present a report on coefficients which are breaking the symmetry is printed.
2358 @type verbose: C{bool}
2359 @return: True if the PDE is symmetric.
2360 @rtype: C{bool}
2361 @note: This is a very expensive operation. It should be used for degugging only! The symmetry flag is not altered.
2362 """
2363 out=True
2364 out=out and self.checkSymmetricTensor("M", verbose)
2365 out=out and self.checkSymmetricTensor("M_reduced", verbose)
2366 out=out and self.checkSymmetricTensor("A", verbose)
2367 out=out and self.checkSymmetricTensor("A_reduced", verbose)
2368 out=out and self.checkReciprocalSymmetry("B","C", verbose)
2369 out=out and self.checkReciprocalSymmetry("B_reduced","C_reduced", verbose)
2370 out=out and self.checkSymmetricTensor("D", verbose)
2371 out=out and self.checkSymmetricTensor("D_reduced", verbose)
2372 out=out and self.checkSymmetricTensor("m", verbose)
2373 out=out and self.checkSymmetricTensor("m_reduced", verbose)
2374 out=out and self.checkSymmetricTensor("d", verbose)
2375 out=out and self.checkSymmetricTensor("d_reduced", verbose)
2376 out=out and self.checkSymmetricTensor("d_contact", verbose)
2377 out=out and self.checkSymmetricTensor("d_contact_reduced", verbose)
2378 return out
2379
2380 def setValue(self,**coefficients):
2381 """
2382 sets new values to coefficients
2383
2384 @param coefficients: new values assigned to coefficients
2385 @keyword M: value for coefficient M.
2386 @type M: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
2387 @keyword M_reduced: value for coefficient M_reduced.
2388 @type M_reduced: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.ReducedFunction>}.
2389 @keyword A: value for coefficient A.
2390 @type A: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
2391 @keyword A_reduced: value for coefficient A_reduced.
2392 @type A_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
2393 @keyword B: value for coefficient B
2394 @type B: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
2395 @keyword B_reduced: value for coefficient B_reduced
2396 @type B_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
2397 @keyword C: value for coefficient C
2398 @type C: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
2399 @keyword C_reduced: value for coefficient C_reduced
2400 @type C_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
2401 @keyword D: value for coefficient D
2402 @type D: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
2403 @keyword D_reduced: value for coefficient D_reduced
2404 @type D_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
2405 @keyword X: value for coefficient X
2406 @type X: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
2407 @keyword X_reduced: value for coefficient X_reduced
2408 @type X_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.ReducedFunction>}.
2409 @keyword Y: value for coefficient Y
2410 @type Y: any type that can be casted to L{Data<escript.Data>} object on L{Function<escript.Function>}.
2411 @keyword Y_reduced: value for coefficient Y_reduced
2412 @type Y_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunction<escript.Function>}.
2413 @keyword m: value for coefficient m
2414 @type m: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
2415 @keyword m_reduced: value for coefficient m_reduced
2416 @type m_reduced: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
2417 @keyword d: value for coefficient d
2418 @type d: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
2419 @keyword d_reduced: value for coefficient d_reduced
2420 @type d_reduced: any type that can be casted to L{Data<escript.Data>} object on L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
2421 @keyword y: value for coefficient y
2422 @type y: any type that can be casted to L{Data<escript.Data>} object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}.
2423 @keyword d_contact: value for coefficient d_contact
2424 @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>}.
2425 @keyword d_contact_reduced: value for coefficient d_contact_reduced
2426 @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>}.
2427 @keyword y_contact: value for coefficient y_contact
2428 @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>}.
2429 @keyword y_contact_reduced: value for coefficient y_contact_reduced
2430 @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>}.
2431 @keyword r: values prescribed to the solution at the locations of constraints
2432 @type r: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2433 depending of reduced order is used for the solution.
2434 @keyword q: mask for location of constraints
2435 @type q: any type that can be casted to L{Data<escript.Data>} object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2436 depending of reduced order is used for the representation of the equation.
2437 @raise IllegalCoefficient: if an unknown coefficient keyword is used.
2438 """
2439 super(TransportPDE,self).setValue(**coefficients)
2440
2441
2442 def createOperator(self):
2443 """
2444 returns an instance of a new transport operator
2445 """
2446 self.trace("New Transport problem is allocated.")
2447 return self.getDomain().newTransportProblem( \
2448 self.getTheta(),
2449 self.getNumEquations(), \
2450 self.getFunctionSpaceForSolution(), \
2451 self.getSystemType())
2452
2453 def setInitialSolution(self,u):
2454 """
2455 sets the initial solution
2456
2457 @param u: new initial solution
2458 @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>}.
2459 @note: C{u} must be non negative.
2460 """
2461 u2=util.interpolate(u,self.getFunctionSpaceForSolution())
2462 if self.getNumSolutions() == 1:
2463 if u2.getShape()!=():
2464 raise ValueError,"Illegal shape %s of initial solution."%(u2.getShape(),)
2465 else:
2466 if u2.getShape()!=(self.getNumSolutions(),):
2467 raise ValueError,"Illegal shape %s of initial solution."%(u2.getShape(),)
2468 self.getOperator().setInitialValue(u2)
2469 def getRequiredSystemType(self):
2470 """
2471 returns the system type which needs to be used by the current set up.
2472
2473 @return: a code to indicate the type of transport problem scheme used.
2474 @rtype: C{float}
2475 """
2476 return self.getDomain().getTransportTypeId(self.getSolverMethod()[0],self.getSolverMethod()[1],self.getSolverPackage(),self.isSymmetric())
2477
2478 def getUnlimitedTimeStepSize(self):
2479 """
2480 returns the value returned by the C{getSafeTimeStepSize} method to indicate no limit on the save time step size.
2481
2482 @return: the value used to indicate that no limit is set to the time step size
2483 @rtype: C{float}
2484 @note: Typically the biggest positive float is returned.
2485 """
2486 return self.getOperator().getUnlimitedTimeStepSize()
2487
2488 def getSafeTimeStepSize(self):
2489 """
2490 returns a safe time step size to do the next time step.
2491
2492 @return: safe time step size
2493 @rtype: C{float}
2494 @note: If not C{getSafeTimeStepSize()} < C{getUnlimitedTimeStepSize()} any time step size can be used.
2495 """
2496 return self.getOperator().getSafeTimeStepSize()
2497 #====================================================================
2498 def getSolution(self,dt,**options):
2499 """
2500 returns the solution of the problem
2501 @return: the solution
2502 @rtype: L{Data<escript.Data>}
2503
2504 @note: This method is overwritten when implementing a particular linear problem.
2505 """
2506 if dt<=0:
2507 raise ValueError,"step size needs to be positive."
2508 options[self.TOLERANCE_KEY]=self.getTolerance()
2509 self.setSolution(self.getOperator().solve(self.getRightHandSide(),dt,options))
2510 self.validSolution()
2511 return self.getCurrentSolution()
2512
2513 def getSystem(self):
2514 """
2515 return the operator and right hand side of the PDE
2516
2517 @return: the discrete version of the PDE
2518 @rtype: C{tuple} of L{Operator,<escript.Operator>} and L{Data<escript.Data>}.
2519
2520 @note: This method is overwritten when implementing a particular linear problem.
2521 """
2522 if not self.isOperatorValid() or not self.isRightHandSideValid():
2523 self.resetRightHandSide()
2524 righthandside=self.getCurrentRightHandSide()
2525 self.resetOperator()
2526 operator=self.getCurrentOperator()
2527 self.getDomain().addPDEToTransportProblem(
2528 operator,
2529 righthandside,
2530 self.getCoefficient("M"),
2531 self.getCoefficient("A"),
2532 self.getCoefficient("B"),
2533 self.getCoefficient("C"),
2534 self.getCoefficient("D"),
2535 self.getCoefficient("X"),
2536 self.getCoefficient("Y"),
2537 self.getCoefficient("d"),
2538 self.getCoefficient("y"),
2539 self.getCoefficient("d_contact"),
2540 self.getCoefficient("y_contact"))
2541 self.getDomain().addPDEToTransportProblem(
2542 operator,
2543 righthandside,
2544 self.getCoefficient("M_reduced"),
2545 self.getCoefficient("A_reduced"),
2546 self.getCoefficient("B_reduced"),
2547 self.getCoefficient("C_reduced"),
2548 self.getCoefficient("D_reduced"),
2549 self.getCoefficient("X_reduced"),
2550 self.getCoefficient("Y_reduced"),
2551 self.getCoefficient("d_reduced"),
2552 self.getCoefficient("y_reduced"),
2553 self.getCoefficient("d_contact_reduced"),
2554 self.getCoefficient("y_contact_reduced"))
2555 operator.insertConstraint(righthandside,self.getCoefficient("q"),self.getCoefficient("r"))
2556 self.trace("New system has been built.")
2557 self.validOperator()
2558 self.validRightHandSide()
2559 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