/[escript]/release/3.4.2/escriptcore/py_src/linearPDEs.py
ViewVC logotype

Contents of /release/3.4.2/escriptcore/py_src/linearPDEs.py

Parent Directory Parent Directory | Revision Log Revision Log


Revision 4924 - (show annotations)
Wed Apr 30 23:27:04 2014 UTC (5 years, 1 month ago) by sshaw
File MIME type: text/x-python
File size: 135355 byte(s)
disabling release use of fast assemblers by default
1 # -*- coding: utf-8 -*-
2
3 ##############################################################################
4 #
5 # Copyright (c) 2003-2014 by University of Queensland
6 # http://www.uq.edu.au
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 # Development until 2012 by Earth Systems Science Computational Center (ESSCC)
13 # Development 2012-2013 by School of Earth Sciences
14 # Development from 2014 by Centre for Geoscience Computing (GeoComp)
15 #
16 ##############################################################################
17
18 __copyright__="""Copyright (c) 2003-2014 by University of Queensland
19 http://www.uq.edu.au
20 Primary Business: Queensland, Australia"""
21 __license__="""Licensed under the Open Software License version 3.0
22 http://www.opensource.org/licenses/osl-3.0.php"""
23 __url__="https://launchpad.net/escript-finley"
24
25 """
26 The module provides an interface to define and solve linear partial
27 differential equations (PDEs) and Transport problems within `escript`.
28 `linearPDEs` does not provide any solver capabilities in itself but hands the
29 PDE over to the PDE solver library defined through the `Domain`
30 of the PDE. The general interface is provided through the `LinearPDE` class.
31 `TransportProblem` provides an interface to initial value problems dominated
32 by its advective terms.
33
34 :var __author__: name of author
35 :var __copyright__: copyrights
36 :var __license__: licence agreement
37 :var __url__: url entry point on documentation
38 :var __version__: version
39 :var __date__: date of the version
40 """
41
42 from . import escriptcpp as escore
43 from . import util
44 import math
45 import numpy
46
47 __author__="Lutz Gross, l.gross@uq.edu.au"
48
49 SolverOptions = escore.SolverOptions
50 SolverBuddy = escore.SolverBuddy
51
52 class IllegalCoefficient(ValueError):
53
54 """
55 Exception that is raised if an illegal coefficient of the general or
56 particular PDE is requested.
57 """
58 pass
59
60 class IllegalCoefficientValue(ValueError):
61 """
62 Exception that is raised if an incorrect value for a coefficient is used.
63 """
64 pass
65
66 class IllegalCoefficientFunctionSpace(ValueError):
67 """
68 Exception that is raised if an incorrect function space for a coefficient
69 is used.
70 """
71
72 class UndefinedPDEError(ValueError):
73 """
74 Exception that is raised if a PDE is not fully defined yet.
75 """
76 pass
77
78 class PDECoef(object):
79 """
80 A class for describing a PDE coefficient.
81
82 :cvar INTERIOR: indicator that coefficient is defined on the interior of
83 the PDE domain
84 :cvar BOUNDARY: indicator that coefficient is defined on the boundary of
85 the PDE domain
86 :cvar CONTACT: indicator that coefficient is defined on the contact region
87 within the PDE domain
88 :cvar INTERIOR_REDUCED: indicator that coefficient is defined on the
89 interior of the PDE domain using a reduced
90 integration order
91 :cvar BOUNDARY_REDUCED: indicator that coefficient is defined on the
92 boundary of the PDE domain using a reduced
93 integration order
94 :cvar CONTACT_REDUCED: indicator that coefficient is defined on the contact
95 region within the PDE domain using a reduced
96 integration order
97 :cvar SOLUTION: indicator that coefficient is defined through a solution of
98 the PDE
99 :cvar REDUCED: indicator that coefficient is defined through a reduced
100 solution of the PDE
101 :cvar DIRACDELTA: indicator that coefficient is defined as Dirac delta functions
102 :cvar BY_EQUATION: indicator that the dimension of the coefficient shape is
103 defined by the number of PDE equations
104 :cvar BY_SOLUTION: indicator that the dimension of the coefficient shape is
105 defined by the number of PDE solutions
106 :cvar BY_DIM: indicator that the dimension of the coefficient shape is
107 defined by the spatial dimension
108 :cvar OPERATOR: indicator that the coefficient alters the operator of
109 the PDE
110 :cvar RIGHTHANDSIDE: indicator that the coefficient alters the right
111 hand side of the PDE
112 :cvar BOTH: indicator that the coefficient alters the operator as well
113 as the right hand side of the PDE
114
115 """
116 INTERIOR=0
117 BOUNDARY=1
118 CONTACT=2
119 SOLUTION=3
120 REDUCED=4
121 BY_EQUATION=5
122 BY_SOLUTION=6
123 BY_DIM=7
124 OPERATOR=10
125 RIGHTHANDSIDE=11
126 BOTH=12
127 INTERIOR_REDUCED=13
128 BOUNDARY_REDUCED=14
129 CONTACT_REDUCED=15
130 DIRACDELTA=16
131
132 def __init__(self, where, pattern, altering):
133 """
134 Initialises a PDE coefficient type.
135
136 :param where: describes where the coefficient lives
137 :type where: one of `INTERIOR`, `BOUNDARY`, `CONTACT`, `SOLUTION`,
138 `REDUCED`, `INTERIOR_REDUCED`, `BOUNDARY_REDUCED`,
139 `CONTACT_REDUCED`, 'DIRACDELTA'
140 :param pattern: describes the shape of the coefficient and how the shape
141 is built for a given spatial dimension and numbers of
142 equations and solutions in then PDE. For instance,
143 (`BY_EQUATION`,`BY_SOLUTION`,`BY_DIM`) describes a
144 rank 3 coefficient which is instantiated as shape (3,2,2)
145 in case of three equations and two solution components
146 on a 2-dimensional domain. In the case of single equation
147 and a single solution component the shape components
148 marked by `BY_EQUATION` or `BY_SOLUTION` are dropped.
149 In this case the example would be read as (2,).
150 :type pattern: ``tuple`` of `BY_EQUATION`, `BY_SOLUTION`, `BY_DIM`
151 :param altering: indicates what part of the PDE is altered if the
152 coefficient is altered
153 :type altering: one of `OPERATOR`, `RIGHTHANDSIDE`, `BOTH`
154 """
155 super(PDECoef, self).__init__()
156 self.what=where
157 self.pattern=pattern
158 self.altering=altering
159 self.resetValue()
160
161 def resetValue(self):
162 """
163 Resets the coefficient value to the default.
164 """
165 self.value=escore.Data()
166
167 def getFunctionSpace(self,domain,reducedEquationOrder=False,reducedSolutionOrder=False):
168 """
169 Returns the `FunctionSpace` of the coefficient.
170
171 :param domain: domain on which the PDE uses the coefficient
172 :type domain: `Domain`
173 :param reducedEquationOrder: True to indicate that reduced order is used
174 to represent the equation
175 :type reducedEquationOrder: ``bool``
176 :param reducedSolutionOrder: True to indicate that reduced order is used
177 to represent the solution
178 :type reducedSolutionOrder: ``bool``
179 :return: `FunctionSpace` of the coefficient
180 :rtype: `FunctionSpace`
181 """
182 if self.what==self.INTERIOR:
183 return escore.Function(domain)
184 elif self.what==self.INTERIOR_REDUCED:
185 return escore.ReducedFunction(domain)
186 elif self.what==self.BOUNDARY:
187 return escore.FunctionOnBoundary(domain)
188 elif self.what==self.BOUNDARY_REDUCED:
189 return escore.ReducedFunctionOnBoundary(domain)
190 elif self.what==self.CONTACT:
191 return escore.FunctionOnContactZero(domain)
192 elif self.what==self.CONTACT_REDUCED:
193 return escore.ReducedFunctionOnContactZero(domain)
194 elif self.what==self.DIRACDELTA:
195 return escore.DiracDeltaFunctions(domain)
196 elif self.what==self.SOLUTION:
197 if reducedEquationOrder and reducedSolutionOrder:
198 return escore.ReducedSolution(domain)
199 else:
200 return escore.Solution(domain)
201 elif self.what==self.REDUCED:
202 return escore.ReducedSolution(domain)
203
204 def getValue(self):
205 """
206 Returns the value of the coefficient.
207
208 :return: value of the coefficient
209 :rtype: `Data`
210 """
211 return self.value
212
213 def setValue(self,domain,numEquations=1,numSolutions=1,reducedEquationOrder=False,reducedSolutionOrder=False,newValue=None):
214 """
215 Sets the value of the coefficient to a new value.
216
217 :param domain: domain on which the PDE uses the coefficient
218 :type domain: `Domain`
219 :param numEquations: number of equations of the PDE
220 :type numEquations: ``int``
221 :param numSolutions: number of components of the PDE solution
222 :type numSolutions: ``int``
223 :param reducedEquationOrder: True to indicate that reduced order is used
224 to represent the equation
225 :type reducedEquationOrder: ``bool``
226 :param reducedSolutionOrder: True to indicate that reduced order is used
227 to represent the solution
228 :type reducedSolutionOrder: ``bool``
229 :param newValue: number of components of the PDE solution
230 :type newValue: any object that can be converted into a
231 `Data` object with the appropriate shape
232 and `FunctionSpace`
233 :raise IllegalCoefficientValue: if the shape of the assigned value does
234 not match the shape of the coefficient
235 :raise IllegalCoefficientFunctionSpace: if unable to interpolate value
236 to appropriate function space
237 """
238 if newValue==None:
239 newValue=escore.Data()
240 elif isinstance(newValue,escore.Data):
241 if not newValue.isEmpty():
242 if not newValue.getFunctionSpace() == self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder):
243 try:
244 newValue=escore.Data(newValue,self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder))
245 except:
246 raise IllegalCoefficientFunctionSpace("Unable to interpolate coefficient to function space %s"%self.getFunctionSpace(domain))
247 else:
248 newValue=escore.Data(newValue,self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder))
249 if not newValue.isEmpty():
250 if not self.getShape(domain,numEquations,numSolutions)==newValue.getShape():
251 raise IllegalCoefficientValue("Expected shape of coefficient is %s but actual shape is %s."%(self.getShape(domain,numEquations,numSolutions),newValue.getShape()))
252 self.value=newValue
253
254 def isAlteringOperator(self):
255 """
256 Checks if the coefficient alters the operator of the PDE.
257
258 :return: True if the operator of the PDE is changed when the
259 coefficient is changed
260 :rtype: ``bool``
261 """
262 if self.altering==self.OPERATOR or self.altering==self.BOTH:
263 return not None
264 else:
265 return None
266
267 def isAlteringRightHandSide(self):
268 """
269 Checks if the coefficient alters the right hand side of the PDE.
270
271 :rtype: ``bool``
272 :return: True if the right hand side of the PDE is changed when the
273 coefficient is changed, ``None`` otherwise.
274 """
275 if self.altering==self.RIGHTHANDSIDE or self.altering==self.BOTH:
276 return not None
277 else:
278 return None
279
280 def estimateNumEquationsAndNumSolutions(self,domain,shape=()):
281 """
282 Tries to estimate the number of equations and number of solutions if
283 the coefficient has the given shape.
284
285 :param domain: domain on which the PDE uses the coefficient
286 :type domain: `Domain`
287 :param shape: suggested shape of the coefficient
288 :type shape: ``tuple`` of ``int`` values
289 :return: the number of equations and number of solutions of the PDE if
290 the coefficient has given shape. If no appropriate numbers
291 could be identified, ``None`` is returned
292 :rtype: ``tuple`` of two ``int`` values or ``None``
293 """
294 dim=domain.getDim()
295 if len(shape)>0:
296 num=max(shape)+1
297 else:
298 num=1
299 search=[]
300 if self.definesNumEquation() and self.definesNumSolutions():
301 for u in range(num):
302 for e in range(num):
303 search.append((e,u))
304 search.sort(key=lambda x: -(x[0]+x[1]))
305 for item in search:
306 s=self.getShape(domain,item[0],item[1])
307 if len(s)==0 and len(shape)==0:
308 return (1,1)
309 else:
310 if s==shape: return item
311 elif self.definesNumEquation():
312 for e in range(num,0,-1):
313 s=self.getShape(domain,e,0)
314 if len(s)==0 and len(shape)==0:
315 return (1,None)
316 else:
317 if s==shape: return (e,None)
318
319 elif self.definesNumSolutions():
320 for u in range(num,0,-1):
321 s=self.getShape(domain,0,u)
322 if len(s)==0 and len(shape)==0:
323 return (None,1)
324 else:
325 if s==shape: return (None,u)
326 return None
327
328 def definesNumSolutions(self):
329 """
330 Checks if the coefficient allows to estimate the number of solution
331 components.
332
333 :return: True if the coefficient allows an estimate of the number of
334 solution components, False otherwise
335 :rtype: ``bool``
336 """
337 for i in self.pattern:
338 if i==self.BY_SOLUTION: return True
339 return False
340
341 def definesNumEquation(self):
342 """
343 Checks if the coefficient allows to estimate the number of equations.
344
345 :return: True if the coefficient allows an estimate of the number of
346 equations, False otherwise
347 :rtype: ``bool``
348 """
349 for i in self.pattern:
350 if i==self.BY_EQUATION: return True
351 return False
352
353 def __CompTuple2(self,t1,t2):
354 """
355 Compares two tuples of possible number of equations and number of
356 solutions.
357
358 :param t1: the first tuple
359 :param t2: the second tuple
360 :return: 0, 1, or -1
361 """
362
363 dif=t1[0]+t1[1]-(t2[0]+t2[1])
364 if dif<0: return 1
365 elif dif>0: return -1
366 else: return 0
367
368 def getShape(self,domain,numEquations=1,numSolutions=1):
369 """
370 Builds the required shape of the coefficient.
371
372 :param domain: domain on which the PDE uses the coefficient
373 :type domain: `Domain`
374 :param numEquations: number of equations of the PDE
375 :type numEquations: ``int``
376 :param numSolutions: number of components of the PDE solution
377 :type numSolutions: ``int``
378 :return: shape of the coefficient
379 :rtype: ``tuple`` of ``int`` values
380 """
381 dim=domain.getDim()
382 s=()
383 for i in self.pattern:
384 if i==self.BY_EQUATION:
385 if numEquations>1: s=s+(numEquations,)
386 elif i==self.BY_SOLUTION:
387 if numSolutions>1: s=s+(numSolutions,)
388 else:
389 s=s+(dim,)
390 return s
391
392 #====================================================================================================================
393
394 class LinearProblem(object):
395 """
396 This is the base class to define a general linear PDE-type problem for
397 for an unknown function *u* on a given domain defined through a
398 `Domain` object. The problem can be given as a single
399 equation or as a system of equations.
400
401 The class assumes that some sort of assembling process is required to form
402 a problem of the form
403
404 *L u=f*
405
406 where *L* is an operator and *f* is the right hand side. This operator
407 problem will be solved to get the unknown *u*.
408
409 """
410 def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):
411 """
412 Initializes a linear problem.
413
414 :param domain: domain of the PDE
415 :type domain: `Domain`
416 :param numEquations: number of equations. If ``None`` the number of
417 equations is extracted from the coefficients.
418 :param numSolutions: number of solution components. If ``None`` the number
419 of solution components is extracted from the
420 coefficients.
421 :param debug: if True debug information is printed
422
423 """
424 super(LinearProblem, self).__init__()
425
426 self.__debug=debug
427 self.__domain=domain
428 self.__numEquations=numEquations
429 self.__numSolutions=numSolutions
430 self.__altered_coefficients=False
431 self.__reduce_equation_order=False
432 self.__reduce_solution_order=False
433 self.__sym=False
434 self.setSolverOptions()
435 self.__is_RHS_valid=False
436 self.__is_operator_valid=False
437 self.__COEFFICIENTS={}
438 self.__solution_rtol=1.e99
439 self.__solution_atol=1.e99
440 self.setSymmetryOff()
441 # initialize things:
442 self.resetAllCoefficients()
443 self.initializeSystem()
444 # ==========================================================================
445 # general stuff:
446 # ==========================================================================
447 def __str__(self):
448 """
449 Returns a string representation of the PDE.
450
451 :return: a simple representation of the PDE
452 :rtype: ``str``
453 """
454 return "<LinearProblem %d>"%id(self)
455 # ==========================================================================
456 # debug :
457 # ==========================================================================
458 def setDebugOn(self):
459 """
460 Switches debug output on.
461 """
462 self.__debug=not None
463
464 def setDebugOff(self):
465 """
466 Switches debug output off.
467 """
468 self.__debug=None
469
470 def setDebug(self, flag):
471 """
472 Switches debug output on if ``flag`` is True otherwise it is switched off.
473
474 :param flag: desired debug status
475 :type flag: ``bool``
476 """
477 if flag:
478 self.setDebugOn()
479 else:
480 self.setDebugOff()
481
482 def trace(self,text):
483 """
484 Prints the text message if debug mode is switched on.
485
486 :param text: message to be printed
487 :type text: ``string``
488 """
489 if self.__debug: print(("%s: %s"%(str(self),text)))
490
491 # ==========================================================================
492 # some service functions:
493 # ==========================================================================
494 def introduceCoefficients(self,**coeff):
495 """
496 Introduces new coefficients into the problem.
497
498 Use:
499
500 p.introduceCoefficients(A=PDECoef(...), B=PDECoef(...))
501
502 to introduce the coefficients *A* and *B*.
503 """
504 for name, type in list(coeff.items()):
505 if not isinstance(type,PDECoef):
506 raise ValueError("coefficient %s has no type."%name)
507 self.__COEFFICIENTS[name]=type
508 self.__COEFFICIENTS[name].resetValue()
509 self.trace("coefficient %s has been introduced."%name)
510 def resetRightHandSideCoefficients(self):
511 """
512 Resets all coefficients defining the right hand side
513 """
514 for name in self.__COEFFICIENTS:
515 if self.__COEFFICIENTS[name].altering == PDECoef.RIGHTHANDSIDE :
516 self.__COEFFICIENTS[name].resetValue()
517 self.trace("coefficient %s has been reset."%name)
518
519 def getDomain(self):
520 """
521 Returns the domain of the PDE.
522
523 :return: the domain of the PDE
524 :rtype: `Domain`
525 """
526 return self.__domain
527 def getDomainStatus(self):
528 """
529 Return the status indicator of the domain
530 """
531 return self.getDomain().getStatus()
532
533 def getSystemStatus(self):
534 """
535 Return the domain status used to build the current system
536 """
537 return self.__system_status
538 def setSystemStatus(self,status=None):
539 """
540 Sets the system status to ``status`` if ``status`` is not present the
541 current status of the domain is used.
542 """
543 if status == None:
544 self.__system_status=self.getDomainStatus()
545 else:
546 self.__system_status=status
547
548 def getDim(self):
549 """
550 Returns the spatial dimension of the PDE.
551
552 :return: the spatial dimension of the PDE domain
553 :rtype: ``int``
554 """
555 return self.getDomain().getDim()
556
557 def getNumEquations(self):
558 """
559 Returns the number of equations.
560
561 :return: the number of equations
562 :rtype: ``int``
563 :raise UndefinedPDEError: if the number of equations is not specified yet
564 """
565 if self.__numEquations==None:
566 if self.__numSolutions==None:
567 raise UndefinedPDEError("Number of equations is undefined. Please specify argument numEquations.")
568 else:
569 self.__numEquations=self.__numSolutions
570 return self.__numEquations
571
572 def getNumSolutions(self):
573 """
574 Returns the number of unknowns.
575
576 :return: the number of unknowns
577 :rtype: ``int``
578 :raise UndefinedPDEError: if the number of unknowns is not specified yet
579 """
580 if self.__numSolutions==None:
581 if self.__numEquations==None:
582 raise UndefinedPDEError("Number of solution is undefined. Please specify argument numSolutions.")
583 else:
584 self.__numSolutions=self.__numEquations
585 return self.__numSolutions
586
587 def reduceEquationOrder(self):
588 """
589 Returns the status of order reduction for the equation.
590
591 :return: True if reduced interpolation order is used for the
592 representation of the equation, False otherwise
593 :rtype: `bool`
594 """
595 return self.__reduce_equation_order
596
597 def reduceSolutionOrder(self):
598 """
599 Returns the status of order reduction for the solution.
600
601 :return: True if reduced interpolation order is used for the
602 representation of the solution, False otherwise
603 :rtype: `bool`
604 """
605 return self.__reduce_solution_order
606
607 def getFunctionSpaceForEquation(self):
608 """
609 Returns the `FunctionSpace` used to discretize
610 the equation.
611
612 :return: representation space of equation
613 :rtype: `FunctionSpace`
614 """
615 if self.reduceEquationOrder():
616 return escore.ReducedSolution(self.getDomain())
617 else:
618 return escore.Solution(self.getDomain())
619
620 def getFunctionSpaceForSolution(self):
621 """
622 Returns the `FunctionSpace` used to represent
623 the solution.
624
625 :return: representation space of solution
626 :rtype: `FunctionSpace`
627 """
628 if self.reduceSolutionOrder():
629 return escore.ReducedSolution(self.getDomain())
630 else:
631 return escore.Solution(self.getDomain())
632
633 # ==========================================================================
634 # solver settings:
635 # ==========================================================================
636 def setSolverOptions(self,options=None):
637 """
638 Sets the solver options.
639
640 :param options: the new solver options. If equal ``None``, the solver options are set to the default.
641 :type options: `SolverOptions` or ``None``
642 :note: The symmetry flag of options is overwritten by the symmetry flag of the `LinearProblem`.
643 """
644 if options==None:
645 self.__solver_options=SolverBuddy()
646 elif isinstance(options, SolverBuddy):
647 self.__solver_options=options
648 else:
649 raise ValueError("options must be a SolverOptions object.")
650 self.__solver_options.setSymmetry(self.__sym)
651
652 def getSolverOptions(self):
653 """
654 Returns the solver options
655
656 :rtype: `SolverOptions`
657 """
658 self.__solver_options.setSymmetry(self.__sym)
659 return self.__solver_options
660
661 def isUsingLumping(self):
662 """
663 Checks if matrix lumping is the current solver method.
664
665 :return: True if the current solver method is lumping
666 :rtype: ``bool``
667 """
668 return self.getSolverOptions().getSolverMethod() in [ SolverOptions.ROWSUM_LUMPING, SolverOptions.HRZ_LUMPING ]
669 # ==========================================================================
670 # symmetry flag:
671 # ==========================================================================
672 def isSymmetric(self):
673 """
674 Checks if symmetry is indicated.
675
676 :return: True if a symmetric PDE is indicated, False otherwise
677 :rtype: ``bool``
678 :note: the method is equivalent to use getSolverOptions().isSymmetric()
679 """
680 self.getSolverOptions().isSymmetric()
681
682 def setSymmetryOn(self):
683 """
684 Sets the symmetry flag.
685 :note: The method overwrites the symmetry flag set by the solver options
686 """
687 self.__sym=True
688 self.getSolverOptions().setSymmetryOn()
689
690 def setSymmetryOff(self):
691 """
692 Clears the symmetry flag.
693 :note: The method overwrites the symmetry flag set by the solver options
694 """
695 self.__sym=False
696 self.getSolverOptions().setSymmetryOff()
697
698 def setSymmetry(self,flag=False):
699 """
700 Sets the symmetry flag to ``flag``.
701
702 :param flag: If True, the symmetry flag is set otherwise reset.
703 :type flag: ``bool``
704 :note: The method overwrites the symmetry flag set by the solver options
705 """
706 self.getSolverOptions().setSymmetry(flag)
707 # ==========================================================================
708 # function space handling for the equation as well as the solution
709 # ==========================================================================
710 def setReducedOrderOn(self):
711 """
712 Switches reduced order on for solution and equation representation.
713
714 :raise RuntimeError: if order reduction is altered after a coefficient has
715 been set
716 """
717 self.setReducedOrderForSolutionOn()
718 self.setReducedOrderForEquationOn()
719
720 def setReducedOrderOff(self):
721 """
722 Switches reduced order off for solution and equation representation
723
724 :raise RuntimeError: if order reduction is altered after a coefficient has
725 been set
726 """
727 self.setReducedOrderForSolutionOff()
728 self.setReducedOrderForEquationOff()
729
730 def setReducedOrderTo(self,flag=False):
731 """
732 Sets order reduction state for both solution and equation representation
733 according to flag.
734
735 :param flag: if True, the order reduction is switched on for both solution
736 and equation representation, otherwise or if flag is not
737 present order reduction is switched off
738 :type flag: ``bool``
739 :raise RuntimeError: if order reduction is altered after a coefficient has
740 been set
741 """
742 self.setReducedOrderForSolutionTo(flag)
743 self.setReducedOrderForEquationTo(flag)
744
745
746 def setReducedOrderForSolutionOn(self):
747 """
748 Switches reduced order on for solution representation.
749
750 :raise RuntimeError: if order reduction is altered after a coefficient has
751 been set
752 """
753 if not self.__reduce_solution_order:
754 if self.__altered_coefficients:
755 raise RuntimeError("order cannot be altered after coefficients have been defined.")
756 self.trace("Reduced order is used for solution representation.")
757 self.__reduce_solution_order=True
758 self.initializeSystem()
759
760 def setReducedOrderForSolutionOff(self):
761 """
762 Switches reduced order off for solution representation
763
764 :raise RuntimeError: if order reduction is altered after a coefficient has
765 been set.
766 """
767 if self.__reduce_solution_order:
768 if self.__altered_coefficients:
769 raise RuntimeError("order cannot be altered after coefficients have been defined.")
770 self.trace("Full order is used to interpolate solution.")
771 self.__reduce_solution_order=False
772 self.initializeSystem()
773
774 def setReducedOrderForSolutionTo(self,flag=False):
775 """
776 Sets order reduction state for solution representation according to flag.
777
778 :param flag: if flag is True, the order reduction is switched on for
779 solution representation, otherwise or if flag is not present
780 order reduction is switched off
781 :type flag: ``bool``
782 :raise RuntimeError: if order reduction is altered after a coefficient has
783 been set
784 """
785 if flag:
786 self.setReducedOrderForSolutionOn()
787 else:
788 self.setReducedOrderForSolutionOff()
789
790 def setReducedOrderForEquationOn(self):
791 """
792 Switches reduced order on for equation representation.
793
794 :raise RuntimeError: if order reduction is altered after a coefficient has
795 been set
796 """
797 if not self.__reduce_equation_order:
798 if self.__altered_coefficients:
799 raise RuntimeError("order cannot be altered after coefficients have been defined.")
800 self.trace("Reduced order is used for test functions.")
801 self.__reduce_equation_order=True
802 self.initializeSystem()
803
804 def setReducedOrderForEquationOff(self):
805 """
806 Switches reduced order off for equation representation.
807
808 :raise RuntimeError: if order reduction is altered after a coefficient has
809 been set
810 """
811 if self.__reduce_equation_order:
812 if self.__altered_coefficients:
813 raise RuntimeError("order cannot be altered after coefficients have been defined.")
814 self.trace("Full order is used for test functions.")
815 self.__reduce_equation_order=False
816 self.initializeSystem()
817
818 def setReducedOrderForEquationTo(self,flag=False):
819 """
820 Sets order reduction state for equation representation according to flag.
821
822 :param flag: if flag is True, the order reduction is switched on for
823 equation representation, otherwise or if flag is not present
824 order reduction is switched off
825 :type flag: ``bool``
826 :raise RuntimeError: if order reduction is altered after a coefficient has
827 been set
828 """
829 if flag:
830 self.setReducedOrderForEquationOn()
831 else:
832 self.setReducedOrderForEquationOff()
833 def getOperatorType(self):
834 """
835 Returns the current system type.
836 """
837 return self.__operator_type
838
839 def checkSymmetricTensor(self,name,verbose=True):
840 """
841 Tests a coefficient for symmetry.
842
843 :param name: name of the coefficient
844 :type name: ``str``
845 :param verbose: if set to True or not present a report on coefficients
846 which break the symmetry is printed.
847 :type verbose: ``bool``
848 :return: True if coefficient ``name`` is symmetric
849 :rtype: ``bool``
850 """
851 SMALL_TOLERANCE=util.EPSILON*10.
852 A=self.getCoefficient(name)
853 verbose=verbose or self.__debug
854 out=True
855 if not A.isEmpty():
856 tol=util.Lsup(A)*SMALL_TOLERANCE
857 s=A.getShape()
858 if A.getRank() == 4:
859 if s[0]==s[2] and s[1] == s[3]:
860 for i in range(s[0]):
861 for j in range(s[1]):
862 for k in range(s[2]):
863 for l in range(s[3]):
864 if util.Lsup(A[i,j,k,l]-A[k,l,i,j])>tol:
865 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)))
866 out=False
867 else:
868 if verbose: print(("non-symmetric problem because of inappropriate shape %s of coefficient %s."%(s,name)))
869 out=False
870 elif A.getRank() == 2:
871 if s[0]==s[1]:
872 for j in range(s[0]):
873 for l in range(s[1]):
874 if util.Lsup(A[j,l]-A[l,j])>tol:
875 if verbose: print(("non-symmetric problem because %s[%d,%d]!=%s[%d,%d]"%(name,j,l,name,l,j)))
876 out=False
877 else:
878 if verbose: print(("non-symmetric problem because of inappropriate shape %s of coefficient %s."%(s,name)))
879 out=False
880 elif A.getRank() == 0:
881 pass
882 else:
883 raise ValueError("Cannot check rank %s of %s."%(A.getRank(),name))
884 return out
885
886 def checkReciprocalSymmetry(self,name0,name1,verbose=True):
887 """
888 Tests two coefficients for reciprocal symmetry.
889
890 :param name0: name of the first coefficient
891 :type name0: ``str``
892 :param name1: name of the second coefficient
893 :type name1: ``str``
894 :param verbose: if set to True or not present a report on coefficients
895 which break the symmetry is printed
896 :type verbose: ``bool``
897 :return: True if coefficients ``name0`` and ``name1`` are reciprocally
898 symmetric.
899 :rtype: ``bool``
900 """
901 SMALL_TOLERANCE=util.EPSILON*10.
902 B=self.getCoefficient(name0)
903 C=self.getCoefficient(name1)
904 verbose=verbose or self.__debug
905 out=True
906 if B.isEmpty() and not C.isEmpty():
907 if verbose: print(("non-symmetric problem because %s is not present but %s is"%(name0,name1)))
908 out=False
909 elif not B.isEmpty() and C.isEmpty():
910 if verbose: print(("non-symmetric problem because %s is not present but %s is"%(name0,name1)))
911 out=False
912 elif not B.isEmpty() and not C.isEmpty():
913 sB=B.getShape()
914 sC=C.getShape()
915 tol=(util.Lsup(B)+util.Lsup(C))*SMALL_TOLERANCE/2.
916 if len(sB) != len(sC):
917 if verbose: print(("non-symmetric problem because ranks of %s (=%s) and %s (=%s) are different."%(name0,len(sB),name1,len(sC))))
918 out=False
919 else:
920 if len(sB)==0:
921 if util.Lsup(B-C)>tol:
922 if verbose: print(("non-symmetric problem because %s!=%s"%(name0,name1)))
923 out=False
924 elif len(sB)==1:
925 if sB[0]==sC[0]:
926 for j in range(sB[0]):
927 if util.Lsup(B[j]-C[j])>tol:
928 if verbose: print(("non-symmetric PDE because %s[%d]!=%s[%d]"%(name0,j,name1,j)))
929 out=False
930 else:
931 if verbose: print(("non-symmetric problem because of inappropriate shapes %s and %s of coefficients %s and %s, respectively."%(sB,sC,name0,name1)))
932 elif len(sB)==3:
933 if sB[0]==sC[1] and sB[1]==sC[2] and sB[2]==sC[0]:
934 for i in range(sB[0]):
935 for j in range(sB[1]):
936 for k in range(sB[2]):
937 if util.Lsup(B[i,j,k]-C[k,i,j])>tol:
938 if verbose: print(("non-symmetric problem because %s[%d,%d,%d]!=%s[%d,%d,%d]"%(name0,i,j,k,name1,k,i,j)))
939 out=False
940 else:
941 if verbose: print(("non-symmetric problem because of inappropriate shapes %s and %s of coefficients %s and %s, respectively."%(sB,sC,name0,name1)))
942 else:
943 raise ValueError("Cannot check rank %s of %s and %s."%(len(sB),name0,name1))
944 return out
945
946 def getCoefficient(self,name):
947 """
948 Returns the value of the coefficient ``name``.
949
950 :param name: name of the coefficient requested
951 :type name: ``string``
952 :return: the value of the coefficient
953 :rtype: `Data`
954 :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
955 """
956 if self.hasCoefficient(name):
957 return self.__COEFFICIENTS[name].getValue()
958 else:
959 raise IllegalCoefficient("illegal coefficient %s requested for general PDE."%name)
960
961 def hasCoefficient(self,name):
962 """
963 Returns True if ``name`` is the name of a coefficient.
964
965 :param name: name of the coefficient enquired
966 :type name: ``string``
967 :return: True if ``name`` is the name of a coefficient of the general PDE,
968 False otherwise
969 :rtype: ``bool``
970 """
971 return name in self.__COEFFICIENTS
972
973 def createCoefficient(self, name):
974 """
975 Creates a `Data` object corresponding to coefficient
976 ``name``.
977
978 :return: the coefficient ``name`` initialized to 0
979 :rtype: `Data`
980 :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
981 """
982 if self.hasCoefficient(name):
983 return escore.Data(0.,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))
984 else:
985 raise IllegalCoefficient("illegal coefficient %s requested for general PDE."%name)
986
987 def getFunctionSpaceForCoefficient(self,name):
988 """
989 Returns the `FunctionSpace` to be used for
990 coefficient ``name``.
991
992 :param name: name of the coefficient enquired
993 :type name: ``string``
994 :return: the function space to be used for coefficient ``name``
995 :rtype: `FunctionSpace`
996 :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
997 """
998 if self.hasCoefficient(name):
999 return self.__COEFFICIENTS[name].getFunctionSpace(self.getDomain())
1000 else:
1001 raise ValueError("unknown coefficient %s requested"%name)
1002
1003 def getShapeOfCoefficient(self,name):
1004 """
1005 Returns the shape of the coefficient ``name``.
1006
1007 :param name: name of the coefficient enquired
1008 :type name: ``string``
1009 :return: the shape of the coefficient ``name``
1010 :rtype: ``tuple`` of ``int``
1011 :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
1012 """
1013 if self.hasCoefficient(name):
1014 return self.__COEFFICIENTS[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())
1015 else:
1016 raise IllegalCoefficient("illegal coefficient %s requested for general PDE."%name)
1017
1018 def resetAllCoefficients(self):
1019 """
1020 Resets all coefficients to their default values.
1021 """
1022 for i in list(self.__COEFFICIENTS.keys()):
1023 self.__COEFFICIENTS[i].resetValue()
1024
1025 def alteredCoefficient(self,name):
1026 """
1027 Announces that coefficient ``name`` has been changed.
1028
1029 :param name: name of the coefficient affected
1030 :type name: ``string``
1031 :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
1032 :note: if ``name`` is q or r, the method will not trigger a rebuild of the
1033 system as constraints are applied to the solved system.
1034 """
1035 if self.hasCoefficient(name):
1036 self.trace("Coefficient %s has been altered."%name)
1037 if not ((name=="q" or name=="r") and self.isUsingLumping()):
1038 if self.__COEFFICIENTS[name].isAlteringOperator(): self.invalidateOperator()
1039 if self.__COEFFICIENTS[name].isAlteringRightHandSide(): self.invalidateRightHandSide()
1040 else:
1041 raise IllegalCoefficient("illegal coefficient %s requested for general PDE."%name)
1042
1043 def validSolution(self):
1044 """
1045 Marks the solution as valid.
1046 """
1047 self.__is_solution_valid=True
1048
1049 def invalidateSolution(self):
1050 """
1051 Indicates the PDE has to be resolved if the solution is requested.
1052 """
1053 self.trace("System will be resolved.")
1054 self.__is_solution_valid=False
1055
1056 def isSolutionValid(self):
1057 """
1058 Returns True if the solution is still valid.
1059 """
1060 if not self.getDomainStatus()==self.getSystemStatus(): self.invalidateSolution()
1061 if self.__solution_rtol>self.getSolverOptions().getTolerance() or \
1062 self.__solution_atol>self.getSolverOptions().getAbsoluteTolerance():
1063 self.invalidateSolution()
1064 return self.__is_solution_valid
1065
1066 def validOperator(self):
1067 """
1068 Marks the operator as valid.
1069 """
1070 self.__is_operator_valid=True
1071
1072 def invalidateOperator(self):
1073 """
1074 Indicates the operator has to be rebuilt next time it is used.
1075 """
1076 self.trace("Operator will be rebuilt.")
1077 self.invalidateSolution()
1078 self.__is_operator_valid=False
1079
1080 def isOperatorValid(self):
1081 """
1082 Returns True if the operator is still valid.
1083 """
1084 if not self.getDomainStatus()==self.getSystemStatus(): self.invalidateOperator()
1085 if not self.getRequiredOperatorType()==self.getOperatorType(): self.invalidateOperator()
1086 return self.__is_operator_valid
1087
1088 def validRightHandSide(self):
1089 """
1090 Marks the right hand side as valid.
1091 """
1092 self.__is_RHS_valid=True
1093
1094 def invalidateRightHandSide(self):
1095 """
1096 Indicates the right hand side has to be rebuilt next time it is used.
1097 """
1098 self.trace("Right hand side has to be rebuilt.")
1099 self.invalidateSolution()
1100 self.__is_RHS_valid=False
1101
1102 def isRightHandSideValid(self):
1103 """
1104 Returns True if the operator is still valid.
1105 """
1106 if not self.getDomainStatus()==self.getSystemStatus(): self.invalidateRightHandSide()
1107 return self.__is_RHS_valid
1108
1109 def invalidateSystem(self):
1110 """
1111 Announces that everything has to be rebuilt.
1112 """
1113 self.invalidateSolution()
1114 self.invalidateOperator()
1115 self.invalidateRightHandSide()
1116
1117 def isSystemValid(self):
1118 """
1119 Returns True if the system (including solution) is still vaild.
1120 """
1121 return self.isSolutionValid() and self.isOperatorValid() and self.isRightHandSideValid()
1122
1123 def initializeSystem(self):
1124 """
1125 Resets the system clearing the operator, right hand side and solution.
1126 """
1127 self.trace("New System has been created.")
1128 self.__operator_type=None
1129 self.setSystemStatus()
1130 self.__operator=escore.Operator()
1131 self.__righthandside=escore.Data()
1132 self.__solution=escore.Data()
1133 self.invalidateSystem()
1134
1135 def getOperator(self):
1136 """
1137 Returns the operator of the linear problem.
1138
1139 :return: the operator of the problem
1140 """
1141 return self.getSystem()[0]
1142
1143 def getRightHandSide(self):
1144 """
1145 Returns the right hand side of the linear problem.
1146
1147 :return: the right hand side of the problem
1148 :rtype: `Data`
1149 """
1150 return self.getSystem()[1]
1151
1152 def createRightHandSide(self):
1153 """
1154 Returns an instance of a new right hand side.
1155 """
1156 self.trace("New right hand side is allocated.")
1157 if self.getNumEquations()>1:
1158 return escore.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)
1159 else:
1160 return escore.Data(0.,(),self.getFunctionSpaceForEquation(),True)
1161
1162 def createSolution(self):
1163 """
1164 Returns an instance of a new solution.
1165 """
1166 self.trace("New solution is allocated.")
1167 if self.getNumSolutions()>1:
1168 return escore.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)
1169 else:
1170 return escore.Data(0.,(),self.getFunctionSpaceForSolution(),True)
1171
1172 def resetSolution(self):
1173 """
1174 Sets the solution to zero.
1175 """
1176 if self.__solution.isEmpty():
1177 self.__solution=self.createSolution()
1178 else:
1179 self.__solution.setToZero()
1180 self.trace("Solution is reset to zero.")
1181
1182 def setSolution(self,u, validate=True):
1183 """
1184 Sets the solution assuming that makes the system valid with the tolrance
1185 defined by the solver options
1186 """
1187 if validate:
1188 self.__solution_rtol=self.getSolverOptions().getTolerance()
1189 self.__solution_atol=self.getSolverOptions().getAbsoluteTolerance()
1190 self.validSolution()
1191 self.__solution=u
1192 def getCurrentSolution(self):
1193 """
1194 Returns the solution in its current state.
1195 """
1196 if self.__solution.isEmpty(): self.__solution=self.createSolution()
1197 return self.__solution
1198
1199 def resetRightHandSide(self):
1200 """
1201 Sets the right hand side to zero.
1202 """
1203 if self.__righthandside.isEmpty():
1204 self.__righthandside=self.createRightHandSide()
1205 else:
1206 self.__righthandside.setToZero()
1207 self.trace("Right hand side is reset to zero.")
1208
1209 def getCurrentRightHandSide(self):
1210 """
1211 Returns the right hand side in its current state.
1212 """
1213 if self.__righthandside.isEmpty(): self.__righthandside=self.createRightHandSide()
1214 return self.__righthandside
1215
1216 def resetOperator(self):
1217 """
1218 Makes sure that the operator is instantiated and returns it initialized
1219 with zeros.
1220 """
1221 if self.getOperatorType() == None:
1222 if self.isUsingLumping():
1223 self.__operator=self.createSolution()
1224 else:
1225 self.__operator=self.createOperator()
1226 self.__operator_type=self.getRequiredOperatorType()
1227 else:
1228 if self.isUsingLumping():
1229 self.__operator.setToZero()
1230 else:
1231 if self.getOperatorType() == self.getRequiredOperatorType():
1232 self.__operator.resetValues()
1233 else:
1234 self.__operator=self.createOperator()
1235 self.__operator_type=self.getRequiredOperatorType()
1236 self.trace("Operator reset to zero")
1237
1238 def getCurrentOperator(self):
1239 """
1240 Returns the operator in its current state.
1241 """
1242 return self.__operator
1243
1244 def setValue(self,**coefficients):
1245 """
1246 Sets new values to coefficients.
1247
1248 :raise IllegalCoefficient: if an unknown coefficient keyword is used
1249 """
1250 # check if the coefficients are legal:
1251 for i in list(coefficients.keys()):
1252 if not self.hasCoefficient(i):
1253 raise IllegalCoefficient("Attempt to set unknown coefficient %s"%i)
1254 # if the number of unknowns or equations is still unknown we try to estimate them:
1255 if self.__numEquations==None or self.__numSolutions==None:
1256 for i,d in list(coefficients.items()):
1257 if hasattr(d,"shape"):
1258 s=d.shape
1259 elif isinstance(d, escore.Data) and not d.isEmpty():
1260 s=d.getShape()
1261 else:
1262 s=numpy.array(d).shape
1263 if s!=None:
1264 # get number of equations and number of unknowns:
1265 res=self.__COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(self.getDomain(),s)
1266 if res==None:
1267 raise IllegalCoefficientValue("Illegal shape %s of coefficient %s"%(s,i))
1268 else:
1269 if self.__numEquations==None: self.__numEquations=res[0]
1270 if self.__numSolutions==None: self.__numSolutions=res[1]
1271 if self.__numEquations==None: raise UndefinedPDEError("unidentified number of equations")
1272 if self.__numSolutions==None: raise UndefinedPDEError("unidentified number of solutions")
1273 # now we check the shape of the coefficient if numEquations and numSolutions are set:
1274 for i,d in list(coefficients.items()):
1275 try:
1276 self.__COEFFICIENTS[i].setValue(self.getDomain(),
1277 self.getNumEquations(),self.getNumSolutions(),
1278 self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
1279 self.alteredCoefficient(i)
1280 except IllegalCoefficientFunctionSpace as m:
1281 # if the function space is wrong then we try the reduced version:
1282 i_red=i+"_reduced"
1283 if (not i_red in list(coefficients.keys())) and i_red in list(self.__COEFFICIENTS.keys()):
1284 try:
1285 self.__COEFFICIENTS[i_red].setValue(self.getDomain(),
1286 self.getNumEquations(),self.getNumSolutions(),
1287 self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
1288 self.alteredCoefficient(i_red)
1289 except IllegalCoefficientValue as m:
1290 raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
1291 except IllegalCoefficientFunctionSpace as m:
1292 raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m))
1293 else:
1294 raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m))
1295 except IllegalCoefficientValue as m:
1296 raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
1297 self.__altered_coefficients=True
1298
1299 # ==========================================================================
1300 # methods that are typically overwritten when implementing a particular
1301 # linear problem
1302 # ==========================================================================
1303 def getRequiredOperatorType(self):
1304 """
1305 Returns the system type which needs to be used by the current set up.
1306
1307 :note: Typically this method is overwritten when implementing a
1308 particular linear problem.
1309 """
1310 return None
1311
1312 def createOperator(self):
1313 """
1314 Returns an instance of a new operator.
1315
1316 :note: This method is overwritten when implementing a particular
1317 linear problem.
1318 """
1319 return escore.Operator()
1320
1321 def checkSymmetry(self,verbose=True):
1322 """
1323 Tests the PDE for symmetry.
1324
1325 :param verbose: if set to True or not present a report on coefficients
1326 which break the symmetry is printed
1327 :type verbose: ``bool``
1328 :return: True if the problem is symmetric
1329 :rtype: ``bool``
1330 :note: Typically this method is overwritten when implementing a
1331 particular linear problem.
1332 """
1333 out=True
1334 return out
1335
1336 def getSolution(self,**options):
1337 """
1338 Returns the solution of the problem.
1339
1340 :return: the solution
1341 :rtype: `Data`
1342
1343 :note: This method is overwritten when implementing a particular
1344 linear problem.
1345 """
1346 return self.getCurrentSolution()
1347
1348 def getSystem(self):
1349 """
1350 Returns the operator and right hand side of the PDE.
1351
1352 :return: the discrete version of the PDE
1353 :rtype: ``tuple`` of `Operator` and `Data`.
1354
1355 :note: This method is overwritten when implementing a particular
1356 linear problem.
1357 """
1358 return (self.getCurrentOperator(), self.getCurrentRightHandSide())
1359
1360 class LinearPDE(LinearProblem):
1361 """
1362 This class is used to define a general linear, steady, second order PDE
1363 for an unknown function *u* on a given domain defined through a
1364 `Domain` object.
1365
1366 For a single PDE having a solution with a single component the linear PDE
1367 is defined in the following form:
1368
1369 *-(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)*
1370
1371 where *grad(F)* denotes the spatial derivative of *F*. Einstein's
1372 summation convention, ie. summation over indexes appearing twice in a term
1373 of a sum performed, is used.
1374 The coefficients *A*, *B*, *C*, *D*, *X* and *Y* have to be specified
1375 through `Data` objects in `Function` and
1376 the coefficients *A_reduced*, *B_reduced*, *C_reduced*, *D_reduced*,
1377 *X_reduced* and *Y_reduced* have to be specified through
1378 `Data` objects in `ReducedFunction`.
1379 It is also allowed to use objects that can be converted into such
1380 `Data` objects. *A* and *A_reduced* are rank two, *B*,
1381 *C*, *X*, *B_reduced*, *C_reduced* and *X_reduced* are rank one and
1382 *D*, *D_reduced*, *Y* and *Y_reduced* are scalar.
1383
1384 The following natural boundary conditions are considered:
1385
1386 *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*
1387
1388 where *n* is the outer normal field. Notice that the coefficients *A*,
1389 *A_reduced*, *B*, *B_reduced*, *X* and *X_reduced* are defined in the
1390 PDE. The coefficients *d* and *y* are each a scalar in
1391 `FunctionOnBoundary` and the coefficients
1392 *d_reduced* and *y_reduced* are each a scalar in
1393 `ReducedFunctionOnBoundary`.
1394
1395 Constraints for the solution prescribe the value of the solution at certain
1396 locations in the domain. They have the form
1397
1398 *u=r* where *q>0*
1399
1400 *r* and *q* are each scalar where *q* is the characteristic function
1401 defining where the constraint is applied. The constraints override any
1402 other condition set by the PDE or the boundary condition.
1403
1404 The PDE is symmetrical if
1405
1406 *A[i,j]=A[j,i]* and *B[j]=C[j]* and *A_reduced[i,j]=A_reduced[j,i]*
1407 and *B_reduced[j]=C_reduced[j]*
1408
1409 For a system of PDEs and a solution with several components the PDE has the
1410 form
1411
1412 *-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]*
1413
1414 *A* and *A_reduced* are of rank four, *B*, *B_reduced*, *C* and
1415 *C_reduced* are each of rank three, *D*, *D_reduced*, *X_reduced* and
1416 *X* are each of rank two and *Y* and *Y_reduced* are of rank one.
1417 The natural boundary conditions take the form:
1418
1419 *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]*
1420
1421 The coefficient *d* is of rank two and *y* is of rank one both in
1422 `FunctionOnBoundary`. The coefficients
1423 *d_reduced* is of rank two and *y_reduced* is of rank one both in
1424 `ReducedFunctionOnBoundary`.
1425
1426 Constraints take the form
1427
1428 *u[i]=r[i]* where *q[i]>0*
1429
1430 *r* and *q* are each rank one. Notice that at some locations not
1431 necessarily all components must have a constraint.
1432
1433 The system of PDEs is symmetrical if
1434
1435 - *A[i,j,k,l]=A[k,l,i,j]*
1436 - *A_reduced[i,j,k,l]=A_reduced[k,l,i,j]*
1437 - *B[i,j,k]=C[k,i,j]*
1438 - *B_reduced[i,j,k]=C_reduced[k,i,j]*
1439 - *D[i,k]=D[i,k]*
1440 - *D_reduced[i,k]=D_reduced[i,k]*
1441 - *d[i,k]=d[k,i]*
1442 - *d_reduced[i,k]=d_reduced[k,i]*
1443
1444 `LinearPDE` also supports solution discontinuities over a contact region
1445 in the domain. To specify the conditions across the discontinuity we are
1446 using the generalised flux *J* which, in the case of a system of PDEs
1447 and several components of the solution, is defined as
1448
1449 *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]*
1450
1451 For the case of single solution component and single PDE *J* is defined as
1452
1453 *J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[j]+(B[i]+B_reduced[i])*u-X[i]-X_reduced[i]*
1454
1455 In the context of discontinuities *n* denotes the normal on the
1456 discontinuity pointing from side 0 towards side 1 calculated from
1457 `FunctionSpace.getNormal` of `FunctionOnContactZero`.
1458 For a system of PDEs the contact condition takes the form
1459
1460 *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]*
1461
1462 where *J0* and *J1* are the fluxes on side 0 and side 1 of the
1463 discontinuity, respectively. *jump(u)*, which is the difference of the
1464 solution at side 1 and at side 0, denotes the jump of *u* across
1465 discontinuity along the normal calculated by `jump`.
1466 The coefficient *d_contact* is of rank two and *y_contact* is of rank one
1467 both in `FunctionOnContactZero` or
1468 `FunctionOnContactOne`.
1469 The coefficient *d_contact_reduced* is of rank two and *y_contact_reduced*
1470 is of rank one both in `ReducedFunctionOnContactZero`
1471 or `ReducedFunctionOnContactOne`.
1472 In case of a single PDE and a single component solution the contact
1473 condition takes the form
1474
1475 *n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)*
1476
1477 In this case the coefficient *d_contact* and *y_contact* are each scalar
1478 both in `FunctionOnContactZero` or
1479 `FunctionOnContactOne` and the coefficient
1480 *d_contact_reduced* and *y_contact_reduced* are each scalar both in
1481 `ReducedFunctionOnContactZero` or
1482 `ReducedFunctionOnContactOne`.
1483
1484 Typical usage::
1485
1486 p = LinearPDE(dom)
1487 p.setValue(A=kronecker(dom), D=1, Y=0.5)
1488 u = p.getSolution()
1489
1490 """
1491
1492 def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):
1493 """
1494 Initializes a new linear PDE.
1495
1496 :param domain: domain of the PDE
1497 :type domain: `Domain`
1498 :param numEquations: number of equations. If ``None`` the number of
1499 equations is extracted from the PDE coefficients.
1500 :param numSolutions: number of solution components. If ``None`` the number
1501 of solution components is extracted from the PDE
1502 coefficients.
1503 :param debug: if True debug information is printed
1504
1505 """
1506 super(LinearPDE, self).__init__(domain,numEquations,numSolutions,debug)
1507 #
1508 # the coefficients of the PDE:
1509 #
1510 self.introduceCoefficients(
1511 A=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
1512 B=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1513 C=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
1514 D=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1515 X=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
1516 Y=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1517 d=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1518 y=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1519 d_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1520 y_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1521 A_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
1522 B_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1523 C_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
1524 D_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1525 X_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
1526 Y_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1527 d_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1528 y_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1529 d_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1530 y_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1531 d_dirac=PDECoef(PDECoef.DIRACDELTA,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
1532 y_dirac=PDECoef(PDECoef.DIRACDELTA,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
1533 r=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.RIGHTHANDSIDE),
1534 q=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.BOTH) )
1535
1536 def __str__(self):
1537 """
1538 Returns the string representation of the PDE.
1539
1540 :return: a simple representation of the PDE
1541 :rtype: ``str``
1542 """
1543 return "<LinearPDE %d>"%id(self)
1544
1545 def getRequiredOperatorType(self):
1546 """
1547 Returns the system type which needs to be used by the current set up.
1548 """
1549 if self.isUsingLumping():
1550 return "__ESCRIPT_DATA"
1551 else:
1552 solver_options=self.getSolverOptions()
1553 return self.getDomain().getSystemMatrixTypeId(solver_options.getSolverMethod(), solver_options.getPreconditioner(),solver_options.getPackage(), solver_options.isSymmetric())
1554
1555 def checkSymmetry(self,verbose=True):
1556 """
1557 Tests the PDE for symmetry.
1558
1559 :param verbose: if set to True or not present a report on coefficients
1560 which break the symmetry is printed.
1561 :type verbose: ``bool``
1562 :return: True if the PDE is symmetric
1563 :rtype: `bool`
1564 :note: This is a very expensive operation. It should be used for
1565 degugging only! The symmetry flag is not altered.
1566 """
1567 out=True
1568 out=out and self.checkSymmetricTensor("A", verbose)
1569 out=out and self.checkSymmetricTensor("A_reduced", verbose)
1570 out=out and self.checkReciprocalSymmetry("B","C", verbose)
1571 out=out and self.checkReciprocalSymmetry("B_reduced","C_reduced", verbose)
1572 out=out and self.checkSymmetricTensor("D", verbose)
1573 out=out and self.checkSymmetricTensor("D_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_contact", verbose)
1577 out=out and self.checkSymmetricTensor("d_contact_reduced", verbose)
1578 out=out and self.checkSymmetricTensor("d_dirac", verbose)
1579 return out
1580
1581 def createOperator(self):
1582 """
1583 Returns an instance of a new operator.
1584 """
1585 optype=self.getRequiredOperatorType()
1586 self.trace("New operator of type %s is allocated."%optype)
1587 return self.getDomain().newOperator( \
1588 self.getNumEquations(), \
1589 self.getFunctionSpaceForEquation(), \
1590 self.getNumSolutions(), \
1591 self.getFunctionSpaceForSolution(), \
1592 optype)
1593
1594 def getSolution(self):
1595 """
1596 Returns the solution of the PDE.
1597
1598 :return: the solution
1599 :rtype: `Data`
1600 """
1601 option_class=self.getSolverOptions()
1602 if not self.isSolutionValid():
1603 mat,f=self.getSystem()
1604 if self.isUsingLumping():
1605 if not util.inf(abs(mat)) > 0.:
1606 raise ZeroDivisionError("Lumped mass matrix has zero entry (try order 1 elements or HRZ lumping).")
1607 self.setSolution(f*1/mat)
1608 else:
1609 self.trace("PDE is resolved.")
1610 self.trace("solver options: %s"%str(option_class))
1611 self.setSolution(mat.solve(f,option_class))
1612 return self.getCurrentSolution()
1613
1614 def getSystem(self):
1615 """
1616 Returns the operator and right hand side of the PDE.
1617
1618 :return: the discrete version of the PDE
1619 :rtype: ``tuple`` of `Operator` and
1620 `Data`
1621 """
1622 if not self.isOperatorValid() or not self.isRightHandSideValid():
1623 if self.isUsingLumping():
1624 if not self.isOperatorValid():
1625 if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():
1626 raise TypeError("Lumped matrix requires same order for equations and unknowns")
1627 if not self.getCoefficient("A").isEmpty():
1628 raise ValueError("coefficient A in lumped matrix may not be present.")
1629 if not self.getCoefficient("B").isEmpty():
1630 raise ValueError("coefficient B in lumped matrix may not be present.")
1631 if not self.getCoefficient("C").isEmpty():
1632 raise ValueError("coefficient C in lumped matrix may not be present.")
1633 if not self.getCoefficient("d_contact").isEmpty():
1634 raise ValueError("coefficient d_contact in lumped matrix may not be present.")
1635 if not self.getCoefficient("A_reduced").isEmpty():
1636 raise ValueError("coefficient A_reduced in lumped matrix may not be present.")
1637 if not self.getCoefficient("B_reduced").isEmpty():
1638 raise ValueError("coefficient B_reduced in lumped matrix may not be present.")
1639 if not self.getCoefficient("C_reduced").isEmpty():
1640 raise ValueError("coefficient C_reduced in lumped matrix may not be present.")
1641 if not self.getCoefficient("d_contact_reduced").isEmpty():
1642 raise ValueError("coefficient d_contact_reduced in lumped matrix may not be present.")
1643 D=self.getCoefficient("D")
1644 d=self.getCoefficient("d")
1645 D_reduced=self.getCoefficient("D_reduced")
1646 d_reduced=self.getCoefficient("d_reduced")
1647 d_dirac=self.getCoefficient("d_dirac")
1648
1649 if not D.isEmpty():
1650 if self.getNumSolutions()>1:
1651 D_times_e=util.matrix_mult(D,numpy.ones((self.getNumSolutions(),)))
1652 else:
1653 D_times_e=D
1654 else:
1655 D_times_e=escore.Data()
1656 if not d.isEmpty():
1657 if self.getNumSolutions()>1:
1658 d_times_e=util.matrix_mult(d,numpy.ones((self.getNumSolutions(),)))
1659 else:
1660 d_times_e=d
1661 else:
1662 d_times_e=escore.Data()
1663
1664 if not D_reduced.isEmpty():
1665 if self.getNumSolutions()>1:
1666 D_reduced_times_e=util.matrix_mult(D_reduced,numpy.ones((self.getNumSolutions(),)))
1667 else:
1668 D_reduced_times_e=D_reduced
1669 else:
1670 D_reduced_times_e=escore.Data()
1671
1672 if not d_reduced.isEmpty():
1673 if self.getNumSolutions()>1:
1674 d_reduced_times_e=util.matrix_mult(d_reduced,numpy.ones((self.getNumSolutions(),)))
1675 else:
1676 d_reduced_times_e=d_reduced
1677 else:
1678 d_reduced_times_e=escore.Data()
1679
1680 if not d_dirac.isEmpty():
1681 if self.getNumSolutions()>1:
1682 d_dirac_times_e=util.matrix_mult(d_dirac,numpy.ones((self.getNumSolutions(),)))
1683 else:
1684 d_reduced_dirac_e=d_dirac
1685 else:
1686 d_dirac_times_e=escore.Data()
1687
1688 self.resetOperator()
1689 operator=self.getCurrentOperator()
1690 if hasattr(self.getDomain(), "addPDEToLumpedSystem") :
1691 hrz_lumping=( self.getSolverOptions().getSolverMethod() == SolverOptions.HRZ_LUMPING )
1692 self.getDomain().addPDEToLumpedSystem(operator, D_times_e, d_times_e, d_dirac_times_e, hrz_lumping )
1693 self.getDomain().addPDEToLumpedSystem(operator, D_reduced_times_e, d_reduced_times_e, escore.Data(), hrz_lumping)
1694 else:
1695 self.getDomain().addPDEToRHS(operator, \
1696 escore.Data(), \
1697 D_times_e, \
1698 d_times_e,\
1699 escore.Data(),\
1700 d_dirac_times_e)
1701 self.getDomain().addPDEToRHS(operator, \
1702 escore.Data(), \
1703 D_reduced_times_e, \
1704 d_reduced_times_e,\
1705 escore.Data(), \
1706 escore.Data())
1707 self.trace("New lumped operator has been built.")
1708 if not self.isRightHandSideValid():
1709 self.resetRightHandSide()
1710 righthandside=self.getCurrentRightHandSide()
1711 self.getDomain().addPDEToRHS(righthandside, \
1712 self.getCoefficient("X"), \
1713 self.getCoefficient("Y"),\
1714 self.getCoefficient("y"),\
1715 self.getCoefficient("y_contact"), \
1716 self.getCoefficient("y_dirac"))
1717 self.getDomain().addPDEToRHS(righthandside, \
1718 self.getCoefficient("X_reduced"), \
1719 self.getCoefficient("Y_reduced"),\
1720 self.getCoefficient("y_reduced"),\
1721 self.getCoefficient("y_contact_reduced"), \
1722 escore.Data())
1723 self.trace("New right hand side has been built.")
1724 self.validRightHandSide()
1725 self.insertConstraint(rhs_only=False)
1726 self.validOperator()
1727 else:
1728 if not self.isOperatorValid() and not self.isRightHandSideValid():
1729 self.resetRightHandSide()
1730 righthandside=self.getCurrentRightHandSide()
1731 self.resetOperator()
1732 operator=self.getCurrentOperator()
1733 self.getDomain().addPDEToSystem(operator,righthandside, \
1734 self.getCoefficient("A"), \
1735 self.getCoefficient("B"), \
1736 self.getCoefficient("C"), \
1737 self.getCoefficient("D"), \
1738 self.getCoefficient("X"), \
1739 self.getCoefficient("Y"), \
1740 self.getCoefficient("d"), \
1741 self.getCoefficient("y"), \
1742 self.getCoefficient("d_contact"), \
1743 self.getCoefficient("y_contact"), \
1744 self.getCoefficient("d_dirac"), \
1745 self.getCoefficient("y_dirac"))
1746 self.getDomain().addPDEToSystem(operator,righthandside, \
1747 self.getCoefficient("A_reduced"), \
1748 self.getCoefficient("B_reduced"), \
1749 self.getCoefficient("C_reduced"), \
1750 self.getCoefficient("D_reduced"), \
1751 self.getCoefficient("X_reduced"), \
1752 self.getCoefficient("Y_reduced"), \
1753 self.getCoefficient("d_reduced"), \
1754 self.getCoefficient("y_reduced"), \
1755 self.getCoefficient("d_contact_reduced"), \
1756 self.getCoefficient("y_contact_reduced"), \
1757 escore.Data(), \
1758 escore.Data())
1759 self.insertConstraint(rhs_only=False)
1760 self.trace("New system has been built.")
1761 self.validOperator()
1762 self.validRightHandSide()
1763 elif not self.isRightHandSideValid():
1764 self.resetRightHandSide()
1765 righthandside=self.getCurrentRightHandSide()
1766 self.getDomain().addPDEToRHS(righthandside,
1767 self.getCoefficient("X"), \
1768 self.getCoefficient("Y"),\
1769 self.getCoefficient("y"),\
1770 self.getCoefficient("y_contact"), \
1771 self.getCoefficient("y_dirac") )
1772 self.getDomain().addPDEToRHS(righthandside,
1773 self.getCoefficient("X_reduced"), \
1774 self.getCoefficient("Y_reduced"),\
1775 self.getCoefficient("y_reduced"),\
1776 self.getCoefficient("y_contact_reduced"), \
1777 escore.Data())
1778 self.insertConstraint(rhs_only=True)
1779 self.trace("New right hand side has been built.")
1780 self.validRightHandSide()
1781 elif not self.isOperatorValid():
1782 self.resetOperator()
1783 operator=self.getCurrentOperator()
1784 self.getDomain().addPDEToSystem(operator,escore.Data(), \
1785 self.getCoefficient("A"), \
1786 self.getCoefficient("B"), \
1787 self.getCoefficient("C"), \
1788 self.getCoefficient("D"), \
1789 escore.Data(), \
1790 escore.Data(), \
1791 self.getCoefficient("d"), \
1792 escore.Data(),\
1793 self.getCoefficient("d_contact"), \
1794 escore.Data(), \
1795 self.getCoefficient("d_dirac"), \
1796 escore.Data())
1797 self.getDomain().addPDEToSystem(operator,escore.Data(), \
1798 self.getCoefficient("A_reduced"), \
1799 self.getCoefficient("B_reduced"), \
1800 self.getCoefficient("C_reduced"), \
1801 self.getCoefficient("D_reduced"), \
1802 escore.Data(), \
1803 escore.Data(), \
1804 self.getCoefficient("d_reduced"), \
1805 escore.Data(),\
1806 self.getCoefficient("d_contact_reduced"), \
1807 escore.Data(), \
1808 escore.Data(), \
1809 escore.Data())
1810 self.insertConstraint(rhs_only=False)
1811 self.trace("New operator has been built.")
1812 self.validOperator()
1813 self.setSystemStatus()
1814 self.trace("System status is %s."%self.getSystemStatus())
1815 return (self.getCurrentOperator(), self.getCurrentRightHandSide())
1816
1817 def insertConstraint(self, rhs_only=False):
1818 """
1819 Applies the constraints defined by q and r to the PDE.
1820
1821 :param rhs_only: if True only the right hand side is altered by the
1822 constraint
1823 :type rhs_only: ``bool``
1824 """
1825 q=self.getCoefficient("q")
1826 r=self.getCoefficient("r")
1827 righthandside=self.getCurrentRightHandSide()
1828 operator=self.getCurrentOperator()
1829
1830 if not q.isEmpty():
1831 if r.isEmpty():
1832 r_s=self.createSolution()
1833 else:
1834 r_s=r
1835 if not rhs_only and not operator.isEmpty():
1836 if self.isUsingLumping():
1837 operator.copyWithMask(escore.Data(1.,q.getShape(),q.getFunctionSpace()),q)
1838 else:
1839 row_q=escore.Data(q,self.getFunctionSpaceForEquation())
1840 col_q=escore.Data(q,self.getFunctionSpaceForSolution())
1841 u=self.createSolution()
1842 u.copyWithMask(r_s,col_q)
1843 righthandside-=operator*u
1844 operator.nullifyRowsAndCols(row_q,col_q,1.)
1845 righthandside.copyWithMask(r_s,q)
1846
1847 def setValue(self,**coefficients):
1848 """
1849 Sets new values to coefficients.
1850
1851 :param coefficients: new values assigned to coefficients
1852 :keyword A: value for coefficient ``A``
1853 :type A: any type that can be cast to a `Data` object on
1854 `Function`
1855 :keyword A_reduced: value for coefficient ``A_reduced``
1856 :type A_reduced: any type that can be cast to a `Data`
1857 object on `ReducedFunction`
1858 :keyword B: value for coefficient ``B``
1859 :type B: any type that can be cast to a `Data` object on
1860 `Function`
1861 :keyword B_reduced: value for coefficient ``B_reduced``
1862 :type B_reduced: any type that can be cast to a `Data`
1863 object on `ReducedFunction`
1864 :keyword C: value for coefficient ``C``
1865 :type C: any type that can be cast to a `Data` object on
1866 `Function`
1867 :keyword C_reduced: value for coefficient ``C_reduced``
1868 :type C_reduced: any type that can be cast to a `Data`
1869 object on `ReducedFunction`
1870 :keyword D: value for coefficient ``D``
1871 :type D: any type that can be cast to a `Data` object on
1872 `Function`
1873 :keyword D_reduced: value for coefficient ``D_reduced``
1874 :type D_reduced: any type that can be cast to a `Data`
1875 object on `ReducedFunction`
1876 :keyword X: value for coefficient ``X``
1877 :type X: any type that can be cast to a `Data` object on
1878 `Function`
1879 :keyword X_reduced: value for coefficient ``X_reduced``
1880 :type X_reduced: any type that can be cast to a `Data`
1881 object on `ReducedFunction`
1882 :keyword Y: value for coefficient ``Y``
1883 :type Y: any type that can be cast to a `Data` object on
1884 `Function`
1885 :keyword Y_reduced: value for coefficient ``Y_reduced``
1886 :type Y_reduced: any type that can be cast to a `Data`
1887 object on `ReducedFunction`
1888 :keyword d: value for coefficient ``d``
1889 :type d: any type that can be cast to a `Data` object on
1890 `FunctionOnBoundary`
1891 :keyword d_reduced: value for coefficient ``d_reduced``
1892 :type d_reduced: any type that can be cast to a `Data`
1893 object on `ReducedFunctionOnBoundary`
1894 :keyword y: value for coefficient ``y``
1895 :type y: any type that can be cast to a `Data` object on
1896 `FunctionOnBoundary`
1897 :keyword d_contact: value for coefficient ``d_contact``
1898 :type d_contact: any type that can be cast to a `Data`
1899 object on `FunctionOnContactOne`
1900 or `FunctionOnContactZero`
1901 :keyword d_contact_reduced: value for coefficient ``d_contact_reduced``
1902 :type d_contact_reduced: any type that can be cast to a `Data`
1903 object on `ReducedFunctionOnContactOne`
1904 or `ReducedFunctionOnContactZero`
1905 :keyword y_contact: value for coefficient ``y_contact``
1906 :type y_contact: any type that can be cast to a `Data`
1907 object on `FunctionOnContactOne`
1908 or `FunctionOnContactZero`
1909 :keyword y_contact_reduced: value for coefficient ``y_contact_reduced``
1910 :type y_contact_reduced: any type that can be cast to a `Data`
1911 object on `ReducedFunctionOnContactOne`
1912 or `ReducedFunctionOnContactZero`
1913 :keyword d_dirac: value for coefficient ``d_dirac``
1914 :type d_dirac: any type that can be cast to a `Data` object on `DiracDeltaFunctions`
1915 :keyword y_dirac: value for coefficient ``y_dirac``
1916 :type y_dirac: any type that can be cast to a `Data` object on `DiracDeltaFunctions`
1917 :keyword r: values prescribed to the solution at the locations of
1918 constraints
1919 :type r: any type that can be cast to a `Data` object on
1920 `Solution` or `ReducedSolution`
1921 depending on whether reduced order is used for the solution
1922 :keyword q: mask for location of constraints
1923 :type q: any type that can be cast to a `Data` object on
1924 `Solution` or `ReducedSolution`
1925 depending on whether reduced order is used for the
1926 representation of the equation
1927 :raise IllegalCoefficient: if an unknown coefficient keyword is used
1928 """
1929 super(LinearPDE,self).setValue(**coefficients)
1930 # check if the systrem is inhomogeneous:
1931 if len(coefficients)>0 and not self.isUsingLumping():
1932 q=self.getCoefficient("q")
1933 r=self.getCoefficient("r")
1934 if not q.isEmpty() and not r.isEmpty():
1935 if util.Lsup(q*r)>0.:
1936 self.trace("Inhomogeneous constraint detected.")
1937 self.invalidateSystem()
1938
1939
1940 def getResidual(self,u=None):
1941 """
1942 Returns the residual of u or the current solution if u is not present.
1943
1944 :param u: argument in the residual calculation. It must be representable
1945 in `self.getFunctionSpaceForSolution()`. If u is not present
1946 or equals ``None`` the current solution is used.
1947 :type u: `Data` or None
1948 :return: residual of u
1949 :rtype: `Data`
1950 """
1951 if u==None:
1952 return self.getOperator()*self.getSolution()-self.getRightHandSide()
1953 else:
1954 return self.getOperator()*escore.Data(u,self.getFunctionSpaceForSolution())-self.getRightHandSide()
1955
1956 def getFlux(self,u=None):
1957 """
1958 Returns the flux *J* for a given *u*.
1959
1960 *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]*
1961
1962 or
1963
1964 *J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[l]+(B[j]+B_reduced[j])u-X[j]-X_reduced[j]*
1965
1966 :param u: argument in the flux. If u is not present or equals `None` the
1967 current solution is used.
1968 :type u: `Data` or None
1969 :return: flux
1970 :rtype: `Data`
1971 """
1972 if u==None: u=self.getSolution()
1973 if self.getNumEquations()>1:
1974 out = escore.Data(0.,(self.getNumEquations(),self.getDim()),self.getFunctionSpaceForCoefficient("X"))
1975 else:
1976 out = escore.Data(0.,(self.getDim(), ),self.getFunctionSpaceForCoefficient("X"))
1977
1978 A=self.getCoefficient("A")
1979 if not A.isEmpty():
1980 out+=util.tensormult(A,util.grad(u,self.getFunctionSpaceForCoefficient("A")))
1981
1982 B=self.getCoefficient("B")
1983 if not B.isEmpty():
1984 if B.getRank() == 1:
1985 out+=B * u
1986 else:
1987 out+=util.generalTensorProduct(B,u,axis_offset=1)
1988
1989 X=self.getCoefficient("X")
1990 if not X.isEmpty():
1991 out-=X
1992
1993 A_reduced=self.getCoefficient("A_reduced")
1994 if not A_reduced.isEmpty():
1995 out+=util.tensormult(A_reduced, util.grad(u,self.getFunctionSpaceForCoefficient("A_reduced"))) \
1996
1997 B_reduced=self.getCoefficient("B_reduced")
1998 if not B_reduced.isEmpty():
1999 if B_reduced.getRank() == 1:
2000 out+=B_reduced*u
2001 else:
2002 out+=util.generalTensorProduct(B_reduced,u,axis_offset=1)
2003
2004 X_reduced=self.getCoefficient("X_reduced")
2005 if not X_reduced.isEmpty():
2006 out-=X_reduced
2007 return out
2008
2009 class Poisson(LinearPDE):
2010 """
2011 Class to define a Poisson equation problem. This is generally a
2012 `LinearPDE` of the form
2013
2014 *-grad(grad(u)[j])[j] = f*
2015
2016 with natural boundary conditions
2017
2018 *n[j]*grad(u)[j] = 0*
2019
2020 and constraints:
2021
2022 *u=0* where *q>0*
2023
2024 """
2025
2026 def __init__(self,domain,debug=False):
2027 """
2028 Initializes a new Poisson equation.
2029
2030 :param domain: domain of the PDE
2031 :type domain: `Domain`
2032 :param debug: if True debug information is printed
2033
2034 """
2035 super(Poisson, self).__init__(domain,1,1,debug)
2036 self.introduceCoefficients(
2037 f=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2038 f_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE))
2039 self.setSymmetryOn()
2040
2041 def setValue(self,**coefficients):
2042 """
2043 Sets new values to coefficients.
2044
2045 :param coefficients: new values assigned to coefficients
2046 :keyword f: value for right hand side *f*
2047 :type f: any type that can be cast to a `Scalar` object
2048 on `Function`
2049 :keyword q: mask for location of constraints
2050 :type q: any type that can be cast to a rank zero `Data`
2051 object on `Solution` or
2052 `ReducedSolution` depending on whether
2053 reduced order is used for the representation of the equation
2054 :raise IllegalCoefficient: if an unknown coefficient keyword is used
2055 """
2056 super(Poisson, self).setValue(**coefficients)
2057
2058
2059 def getCoefficient(self,name):
2060 """
2061 Returns the value of the coefficient ``name`` of the general PDE.
2062
2063 :param name: name of the coefficient requested
2064 :type name: ``string``
2065 :return: the value of the coefficient ``name``
2066 :rtype: `Data`
2067 :raise IllegalCoefficient: invalid coefficient name
2068 :note: This method is called by the assembling routine to map the Poisson
2069 equation onto the general PDE.
2070 """
2071 if name == "A" :
2072 return escore.Data(util.kronecker(self.getDim()),escore.Function(self.getDomain()))
2073 elif name == "Y" :
2074 return self.getCoefficient("f")
2075 elif name == "Y_reduced" :
2076 return self.getCoefficient("f_reduced")
2077 else:
2078 return super(Poisson, self).getCoefficient(name)
2079
2080 class Helmholtz(LinearPDE):
2081 """
2082 Class to define a Helmholtz equation problem. This is generally a
2083 `LinearPDE` of the form
2084
2085 *omega*u - grad(k*grad(u)[j])[j] = f*
2086
2087 with natural boundary conditions
2088
2089 *k*n[j]*grad(u)[j] = g- alphau*
2090
2091 and constraints:
2092
2093 *u=r* where *q>0*
2094
2095 """
2096
2097 def __init__(self,domain,debug=False):
2098 """
2099 Initializes a new Helmholtz equation.
2100
2101 :param domain: domain of the PDE
2102 :type domain: `Domain`
2103 :param debug: if True debug information is printed
2104
2105 """
2106 super(Helmholtz, self).__init__(domain,1,1,debug)
2107 self.introduceCoefficients(
2108 omega=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.OPERATOR),
2109 k=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.OPERATOR),
2110 f=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2111 f_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2112 alpha=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.OPERATOR),
2113 g=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2114 g_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE))
2115 self.setSymmetryOn()
2116
2117 def setValue(self,**coefficients):
2118 """
2119 Sets new values to coefficients.
2120
2121 :param coefficients: new values assigned to coefficients
2122 :keyword omega: value for coefficient *omega*
2123 :type omega: any type that can be cast to a `Scalar`
2124 object on `Function`
2125 :keyword k: value for coefficient *k*
2126 :type k: any type that can be cast to a `Scalar` object
2127 on `Function`
2128 :keyword f: value for right hand side *f*
2129 :type f: any type that can be cast to a `Scalar` object
2130 on `Function`
2131 :keyword alpha: value for right hand side *alpha*
2132 :type alpha: any type that can be cast to a `Scalar`
2133 object on `FunctionOnBoundary`
2134 :keyword g: value for right hand side *g*
2135 :type g: any type that can be cast to a `Scalar` object
2136 on `FunctionOnBoundary`
2137 :keyword r: prescribed values *r* for the solution in constraints
2138 :type r: any type that can be cast to a `Scalar` object
2139 on `Solution` or
2140 `ReducedSolution` depending on whether
2141 reduced order is used for the representation of the equation
2142 :keyword q: mask for the location of constraints
2143 :type q: any type that can be cast to a `Scalar` object
2144 on `Solution` or
2145 `ReducedSolution` depending on whether
2146 reduced order is used for the representation of the equation
2147 :raise IllegalCoefficient: if an unknown coefficient keyword is used
2148 """
2149 super(Helmholtz, self).setValue(**coefficients)
2150
2151 def getCoefficient(self,name):
2152 """
2153 Returns the value of the coefficient ``name`` of the general PDE.
2154
2155 :param name: name of the coefficient requested
2156 :type name: ``string``
2157 :return: the value of the coefficient ``name``
2158 :rtype: `Data`
2159 :raise IllegalCoefficient: invalid name
2160 """
2161 if name == "A" :
2162 if self.getCoefficient("k").isEmpty():
2163 return escore.Data(numpy.identity(self.getDim()),escore.Function(self.getDomain()))
2164 else:
2165 return escore.Data(numpy.identity(self.getDim()),escore.Function(self.getDomain()))*self.getCoefficient("k")
2166 elif name == "D" :
2167 return self.getCoefficient("omega")
2168 elif name == "Y" :
2169 return self.getCoefficient("f")
2170 elif name == "d" :
2171 return self.getCoefficient("alpha")
2172 elif name == "y" :
2173 return self.getCoefficient("g")
2174 elif name == "Y_reduced" :
2175 return self.getCoefficient("f_reduced")
2176 elif name == "y_reduced" :
2177 return self.getCoefficient("g_reduced")
2178 else:
2179 return super(Helmholtz, self).getCoefficient(name)
2180
2181 class WavePDE(LinearPDE):
2182 """
2183 A class specifically for waves, passes along values to native implementation
2184 to save computational time.
2185 """
2186 def __init__(self,domain,c,numEquations=None,numSolutions=None,debug=False):
2187 """
2188 Initializes a new linear PDE.
2189
2190 :param domain: domain of the PDE
2191 :type domain: `Domain`
2192 :param numEquations: number of equations. If ``None`` the number of
2193 equations is extracted from the PDE coefficients.
2194 :param numSolutions: number of solution components. If ``None`` the number
2195 of solution components is extracted from the PDE
2196 coefficients.
2197 :param debug: if True debug information is printed
2198
2199 """
2200 super(WavePDE, self).__init__(domain,numEquations,numSolutions,debug)
2201 #
2202 # the coefficients of the PDE:
2203 #
2204 self.introduceCoefficients(
2205 A=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2206 B=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2207 C=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2208 D=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2209 du=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
2210 Y=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2211 d=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2212 y=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2213 d_dirac=PDECoef(PDECoef.DIRACDELTA,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2214 y_dirac=PDECoef(PDECoef.DIRACDELTA,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2215 r=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.RIGHTHANDSIDE),
2216 q=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.BOTH))
2217 domain.setAssembler("WaveAssembler", c)
2218
2219
2220 def getSystem(self):
2221 """
2222 Returns the operator and right hand side of the PDE.
2223
2224 :return: the discrete version of the PDE
2225 :rtype: ``tuple`` of `Operator` and
2226 `Data`
2227 """
2228 if not self.isOperatorValid() or not self.isRightHandSideValid():
2229 if self.isUsingLumping():
2230 if not self.isOperatorValid():
2231 if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():
2232 raise TypeError("Lumped matrix requires same order for equations and unknowns")
2233 if not self.getCoefficient("A").isEmpty():
2234 raise ValueError("coefficient A in lumped matrix may not be present.")
2235 if not self.getCoefficient("B").isEmpty():
2236 raise ValueError("coefficient B in lumped matrix may not be present.")
2237 if not self.getCoefficient("C").isEmpty():
2238 raise ValueError("coefficient C in lumped matrix may not be present.")
2239
2240 D=self.getCoefficient("D")
2241 d=self.getCoefficient("d")
2242 d_dirac=self.getCoefficient("d_dirac")
2243
2244 if not D.isEmpty():
2245 if self.getNumSolutions()>1:
2246 D_times_e=util.matrix_mult(D,numpy.ones((self.getNumSolutions(),)))
2247 else:
2248 D_times_e=D
2249 else:
2250 D_times_e=escore.Data()
2251 if not d.isEmpty():
2252 if self.getNumSolutions()>1:
2253 d_times_e=util.matrix_mult(d,numpy.ones((self.getNumSolutions(),)))
2254 else:
2255 d_times_e=d
2256 else:
2257 d_times_e=escore.Data()
2258
2259 if not d_dirac.isEmpty():
2260 if self.getNumSolutions()>1:
2261 d_dirac_times_e=util.matrix_mult(d_dirac,numpy.ones((self.getNumSolutions(),)))
2262 else:
2263 d_dirac_times_e=escore.Data()
2264
2265 self.resetOperator()
2266 operator=self.getCurrentOperator()
2267 if hasattr(self.getDomain(), "addPDEToLumpedSystem") :
2268 hrz_lumping=( self.getSolverOptions().getSolverMethod() == SolverOptions.HRZ_LUMPING )
2269 self.getDomain().addPDEToLumpedSystem(operator, D_times_e, d_times_e, d_dirac_times_e, hrz_lumping )
2270 else:
2271 self.getDomain().addToRHS(operator,
2272 [("Y", D_times_e), ("y", d_times_e),
2273 ("y_dirac", d_dirac_times_e)])
2274 self.trace("New lumped operator has been built.")
2275 if not self.isRightHandSideValid():
2276 self.resetRightHandSide()
2277 righthandside=self.getCurrentRightHandSide()
2278 self.getDomain().addToRHS(righthandside,
2279 [(i, self.getCoefficient(i)) for i in
2280 ["du", "Y", "y", "y_dirac"]
2281 ])
2282 self.trace("New right hand side has been built.")
2283 self.validRightHandSide()
2284 self.insertConstraint(rhs_only=False)
2285 self.validOperator()
2286 else:
2287 if not self.isOperatorValid() and not self.isRightHandSideValid():
2288 self.resetRightHandSide()
2289 righthandside=self.getCurrentRightHandSide()
2290 self.resetOperator()
2291 operator=self.getCurrentOperator()
2292 data = [(i, self.getCoefficient(i)) for i in ["A", "B", "C",
2293 "D", "Y", "d", "y", "d_contact",
2294 "y_contact", "d_dirac", "y_dirac", "du"]
2295 ]
2296 self.getDomain().addToSystem(operator, righthandside, data)
2297 self.insertConstraint(rhs_only=False)
2298 self.trace("New system has been built.")
2299 self.validOperator()
2300 self.validRightHandSide()
2301 elif not self.isRightHandSideValid():
2302 self.resetRightHandSide()
2303 righthandside=self.getCurrentRightHandSide()
2304 self.getDomain().addToRHS(righthandside,
2305 [(i, self.getCoefficient(i)) for i in
2306 ["du", "Y", "y", "y_contact", "y_dirac"]
2307 ])
2308 self.insertConstraint(rhs_only=True)
2309 self.trace("New right hand side has been built.")
2310 self.validRightHandSide()
2311 elif not self.isOperatorValid():
2312 self.resetOperator()
2313 operator=self.getCurrentOperator()
2314 data = [(i, self.getCoefficient(i)) for i in ["A", "B", "C",
2315 "D", "d", "d_contact", "d_dirac", "du"]]
2316 self.getDomain().addToSystem(operator, escore.Data(), data)
2317 self.insertConstraint(rhs_only=False)
2318 self.trace("New operator has been built.")
2319 self.validOperator()
2320 self.setSystemStatus()
2321 self.trace("System status is %s."%self.getSystemStatus())
2322 return (self.getCurrentOperator(), self.getCurrentRightHandSide())
2323
2324
2325 class LameEquation(LinearPDE):
2326 """
2327 Class to define a Lame equation problem. This problem is defined as:
2328
2329 *-grad(mu*(grad(u[i])[j]+grad(u[j])[i]))[j] - grad(lambda*grad(u[k])[k])[j] = F_i -grad(sigma[ij])[j]*
2330
2331 with natural boundary conditions:
2332
2333 *n[j]*(mu*(grad(u[i])[j]+grad(u[j])[i]) + lambda*grad(u[k])[k]) = f_i +n[j]*sigma[ij]*
2334
2335 and constraints:
2336
2337 *u[i]=r[i]* where *q[i]>0*
2338
2339 """
2340
2341 def __init__(self,domain,debug=False,useFast=False):
2342 """
2343 Initializes a new Lame equation.
2344
2345 :param domain: domain of the PDE
2346 :type domain: `Domain`
2347 :param debug: if True debug information is printed
2348
2349 """
2350 self.fastAssembler = False
2351 if useFast and hasattr(domain, "setAssembler"):
2352 self.fastAssembler = True
2353 domain.setAssembler("LameAssembler", [])
2354 super(LameEquation, self).__init__(domain,\
2355 domain.getDim(),domain.getDim(),debug)
2356 self.introduceCoefficients(lame_lambda=PDECoef(PDECoef.INTERIOR,(),PDECoef.OPERATOR),
2357 lame_mu=PDECoef(PDECoef.INTERIOR,(),PDECoef.OPERATOR),
2358 F=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2359 sigma=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
2360 f=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE))
2361 self.setSymmetryOn()
2362
2363 def setValues(self,**coefficients):
2364 """
2365 Sets new values to coefficients.
2366
2367 :param coefficients: new values assigned to coefficients
2368 :keyword lame_mu: value for coefficient *mu*
2369 :type lame_mu: any type that can be cast to a `Scalar`
2370 object on `Function`
2371 :keyword lame_lambda: value for coefficient *lambda*
2372 :type lame_lambda: any type that can be cast to a `Scalar`
2373 object on `Function`
2374 :keyword F: value for internal force *F*
2375 :type F: any type that can be cast to a `Vector` object
2376 on `Function`
2377 :keyword sigma: value for initial stress *sigma*
2378 :type sigma: any type that can be cast to a `Tensor`
2379 object on `Function`
2380 :keyword f: value for external force *f*
2381 :type f: any type that can be cast to a `Vector` object
2382 on `FunctionOnBoundary`
2383 :keyword r: prescribed values *r* for the solution in constraints
2384 :type r: any type that can be cast to a `Vector` object
2385 on `Solution` or
2386 `ReducedSolution` depending on whether
2387 reduced order is used for the representation of the equation
2388 :keyword q: mask for the location of constraints
2389 :type q: any type that can be cast to a `Vector` object
2390 on `Solution` or
2391 `ReducedSolution` depending on whether
2392 reduced order is used for the representation of the equation
2393 :raise IllegalCoefficient: if an unknown coefficient keyword is used
2394 """
2395 super(LameEquation, self).setValues(**coefficients)
2396
2397 def getCoefficient(self,name):
2398 """
2399 Returns the value of the coefficient ``name`` of the general PDE.
2400
2401 :param name: name of the coefficient requested
2402 :type name: ``string``
2403 :return: the value of the coefficient ``name``
2404 :rtype: `Data`
2405 :raise IllegalCoefficient: invalid coefficient name
2406 """
2407
2408 if name == "A" :
2409 out = self.createCoefficient("A")
2410 if self.getCoefficient("lame_lambda").isEmpty():
2411 if self.getCoefficient("lame_mu").isEmpty():
2412 pass
2413 else:
2414 for i in range(self.getDim()):
2415 for j in range(self.getDim()):
2416 out[i,j,j,i] += self.getCoefficient("lame_mu")
2417 out[i,j,i,j] += self.getCoefficient("lame_mu")
2418 else:
2419 if self.getCoefficient("lame_mu").isEmpty():
2420 for i in range(self.getDim()):
2421 for j in range(self.getDim()):
2422 out[i,i,j,j] += self.getCoefficient("lame_lambda")
2423 else:
2424 for i in range(self.getDim()):
2425 for j in range(self.getDim()):
2426 out[i,i,j,j] += self.getCoefficient("lame_lambda")
2427 out[i,j,j,i] += self.getCoefficient("lame_mu")
2428 out[i,j,i,j] += self.getCoefficient("lame_mu")
2429 return out
2430 elif name == "X" :
2431 return self.getCoefficient("sigma")
2432 elif name == "Y" :
2433 return self.getCoefficient("F")
2434 elif name == "y" :
2435 return self.getCoefficient("f")
2436 return super(LameEquation, self).getCoefficient(name)
2437
2438 def getSystem(self):
2439 """
2440 Returns the operator and right hand side of the PDE.
2441
2442 :return: the discrete version of the PDE
2443 :rtype: ``tuple`` of `Operator` and
2444 `Data`
2445 """
2446
2447 if not self.fastAssembler:
2448 return super(LameEquation, self).getSystem()
2449
2450 if not self.isOperatorValid() or not self.isRightHandSideValid():
2451 if self.isUsingLumping():
2452 if not self.isOperatorValid():
2453 if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():
2454 raise TypeError("Lumped matrix requires same order for equations and unknowns")
2455 if not self.getCoefficient("lame_mu").isEmpty() and not self.getCoefficient("lame_mu").isEmpty():
2456 raise ValueError("coefficient A in lumped matrix may not be present.")
2457 if not self.getCoefficient("B").isEmpty():
2458 raise ValueError("coefficient B in lumped matrix may not be present.")
2459 if not self.getCoefficient("C").isEmpty():
2460 raise ValueError("coefficient C in lumped matrix may not be present.")
2461 if not self.getCoefficient("d_contact").isEmpty():
2462 raise ValueError("coefficient d_contact in lumped matrix may not be present.")
2463 if not self.getCoefficient("A_reduced").isEmpty():
2464 raise ValueError("coefficient A_reduced in lumped matrix may not be present.")
2465 if not self.getCoefficient("B_reduced").isEmpty():
2466 raise ValueError("coefficient B_reduced in lumped matrix may not be present.")
2467 if not self.getCoefficient("C_reduced").isEmpty():
2468 raise ValueError("coefficient C_reduced in lumped matrix may not be present.")
2469 if not self.getCoefficient("d_contact_reduced").isEmpty():
2470 raise ValueError("coefficient d_contact_reduced in lumped matrix may not be present.")
2471 D=self.getCoefficient("D")
2472 d=self.getCoefficient("d")
2473 D_reduced=self.getCoefficient("D_reduced")
2474 d_reduced=self.getCoefficient("d_reduced")
2475 d_dirac=self.getCoefficient("d_dirac")
2476
2477 if not D.isEmpty():
2478 if self.getNumSolutions()>1:
2479 D_times_e=util.matrix_mult(D,numpy.ones((self.getNumSolutions(),)))
2480 else:
2481 D_times_e=D
2482 else:
2483 D_times_e=escore.Data()
2484 if not d.isEmpty():
2485 if self.getNumSolutions()>1:
2486 d_times_e=util.matrix_mult(d,numpy.ones((self.getNumSolutions(),)))
2487 else:
2488 d_times_e=d
2489 else:
2490 d_times_e=escore.Data()
2491
2492 if not D_reduced.isEmpty():
2493 if self.getNumSolutions()>1:
2494 D_reduced_times_e=util.matrix_mult(D_reduced,numpy.ones((self.getNumSolutions(),)))
2495 else:
2496 D_reduced_times_e=D_reduced
2497 else:
2498 D_reduced_times_e=escore.Data()
2499
2500 if not d_reduced.isEmpty():
2501 if self.getNumSolutions()>1:
2502 d_reduced_times_e=util.matrix_mult(d_reduced,numpy.ones((self.getNumSolutions(),)))
2503 else:
2504 d_reduced_times_e=d_reduced
2505 else:
2506 d_reduced_times_e=escore.Data()
2507
2508 if not d_dirac.isEmpty():
2509 if self.getNumSolutions()>1:
2510 d_dirac_times_e=util.matrix_mult(d_dirac,numpy.ones((self.getNumSolutions(),)))
2511 else:
2512 d_reduced_dirac_e=d_dirac
2513 else:
2514 d_dirac_times_e=escore.Data()
2515
2516 self.resetOperator()
2517 operator=self.getCurrentOperator()
2518 if hasattr(self.getDomain(), "addPDEToLumpedSystem") :
2519 hrz_lumping=( self.getSolverOptions().getSolverMethod() == SolverOptions.HRZ_LUMPING )
2520 self.getDomain().addPDEToLumpedSystem(operator, D_times_e, d_times_e, d_dirac_times_e, hrz_lumping )
2521 self.getDomain().addPDEToLumpedSystem(operator, D_reduced_times_e, d_reduced_times_e, escore.Data(), hrz_lumping)
2522 else:
2523 self.getDomain().addToRHS(operator, [
2524 ("Y", D_times_e),
2525 ("y", d_times_e),
2526 ("y_dirac", d_dirac_times_e)])
2527 self.getDomain().addToRHS(operator, [
2528 ("Y",D_reduced_times_e),
2529 ("y",d_reduced_times_e)])
2530 self.trace("New lumped operator has been built.")
2531 if not self.isRightHandSideValid():
2532 self.resetRightHandSide()
2533 righthandside=self.getCurrentRightHandSide()
2534 data = [(i, self.getCoefficient(i)) for i in ["X", "Y", "y",
2535 "y_contact", "y_dirac"]]
2536 self.getDomain().addToRHS(righthandside, data)
2537 data = [(i, self.getCoefficient(i+"_reduced")) for i in ["X",
2538 "Y", "y", "y_contact"]]
2539 self.getDomain().addToRHS(righthandside, data)
2540 self.trace("New right hand side has been built.")
2541 self.validRightHandSide()
2542 self.insertConstraint(rhs_only=False)
2543 self.validOperator()
2544 else:
2545 if not self.isOperatorValid() and not self.isRightHandSideValid():
2546 self.resetRightHandSide()
2547 righthandside=self.getCurrentRightHandSide()
2548 self.resetOperator()
2549 operator=self.getCurrentOperator()
2550 data = [(i, self.getCoefficient(i)) for i in ["lame_mu",
2551 "lame_lambda", "B", "C", "D",
2552 "X", "Y", "d", "y", "d_contact", "y_contact",
2553 "d_dirac", "y_dirac"]]
2554 self.getDomain().addToSystem(operator,righthandside, data)
2555 data = [(i, self.getCoefficient(i+"_reduced")) for i in [
2556 "A", "B", "C", "D",
2557 "X", "Y", "d", "y", "d_contact", "y_contact"]
2558 ]
2559 self.getDomain().addToSystem(operator,righthandside, data)
2560 self.insertConstraint(rhs_only=False)
2561 self.trace("New system has been built.")
2562 self.validOperator()
2563 self.validRightHandSide()
2564 elif not self.isRightHandSideValid():
2565 self.resetRightHandSide()
2566 righthandside=self.getCurrentRightHandSide()
2567 data = [(i, self.getCoefficient(i)) for i in ["X", "Y", "y",
2568 "y_contact", "y_dirac"]]
2569 self.getDomain().addToRHS(righthandside, data)
2570 data = [(i, self.getCoefficient(i+"_reduced")) for i in [
2571 "X", "Y", "y","y_contact"]]
2572 self.getDomain().addToRHS(righthandside, data)
2573 self.insertConstraint(rhs_only=True)
2574 self.trace("New right hand side has been built.")
2575 self.validRightHandSide()
2576 elif not self.isOperatorValid():
2577 self.resetOperator()
2578 operator=self.getCurrentOperator()
2579 data = [(i, self.getCoefficient(i)) for i in ["lame_mu",
2580 "lame_lambda", "B","C","D","d","d_contact","d_dirac"]]
2581 self.getDomain().addToSystem(operator, data)
2582 data = [(i, self.getCoefficient(i)) for i in [
2583 "lame_mu","lame_lambda"]
2584 ] + [(i, self.getCoefficient(i+"_reduced")) for i in [
2585 "B","C","D","d","d_contact"]
2586 ]
2587 self.getDomain().addToSystem(operator,data)
2588 self.insertConstraint(rhs_only=False)
2589 self.trace("New operator has been built.")
2590 self.validOperator()
2591 self.setSystemStatus()
2592 self.trace("System status is %s."%self.getSystemStatus())
2593 return (self.getCurrentOperator(), self.getCurrentRightHandSide())
2594
2595 def LinearSinglePDE(domain,debug=False):
2596 """
2597 Defines a single linear PDE.
2598
2599 :param domain: domain of the PDE
2600 :type domain: `Domain`
2601 :param debug: if True debug information is printed
2602 :rtype: `LinearPDE`
2603 """
2604 return LinearPDE(domain,numEquations=1,numSolutions=1,debug=debug)
2605
2606 def LinearPDESystem(domain,debug=False):
2607 """
2608 Defines a system of linear PDEs.
2609
2610 :param domain: domain of the PDEs
2611 :type domain: `Domain`
2612 :param debug: if True debug information is printed
2613 :rtype: `LinearPDE`
2614 """
2615 return LinearPDE(domain,numEquations=domain.getDim(),numSolutions=domain.getDim(),debug=debug)
2616
2617
2618 class TransportPDE(LinearProblem):
2619 """
2620 This class is used to define a transport problem given by a general linear,
2621 time dependent, second order PDE for an unknown, non-negative,
2622 time-dependent function *u* on a given domain defined through a
2623 `Domain` object.
2624
2625 For a single equation with a solution with a single component the transport
2626 problem is defined in the following form:
2627
2628 *(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)*
2629
2630 where *u_t* denotes the time derivative of *u* and *grad(F)* denotes the
2631 spatial derivative of *F*. Einstein's summation convention, ie. summation
2632 over indexes appearing twice in a term of a sum performed, is used.
2633 The coefficients *M*, *A*, *B*, *C*, *D*, *X* and *Y* have to be
2634 specified through `Data` objects in `Function`
2635 and the coefficients *M_reduced*, *A_reduced*, *B_reduced*, *C_reduced*,
2636 *D_reduced*, *X_reduced* and *Y_reduced* have to be specified through
2637 `Data` objects in `ReducedFunction`.
2638 It is also allowed to use objects that can be converted into such
2639 `Data` objects. *M* and *M_reduced* are scalar, *A* and
2640 *A_reduced* are rank two, *B*, *C*, *X*, *B_reduced*, *C_reduced* and
2641 *X_reduced* are rank one and *D*, *D_reduced*, *Y* and *Y_reduced* are
2642 scalar.
2643
2644 The following natural boundary conditions are considered:
2645
2646 *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*
2647
2648 where *n* is the outer normal field. Notice that the coefficients *A*,
2649 *A_reduced*, *B*, *B_reduced*, *X* and *X_reduced* are defined in the
2650 transport problem. The coefficients *m*, *d* and *y* are each a scalar in
2651 `FunctionOnBoundary` and the coefficients
2652 *m_reduced*, *d_reduced* and *y_reduced* are each a scalar in
2653 `ReducedFunctionOnBoundary`.
2654
2655 Constraints for the solution prescribing the value of the solution at
2656 certain locations in the domain have the form
2657
2658 *u_t=r* where *q>0*
2659
2660 *r* and *q* are each scalar where *q* is the characteristic function
2661 defining where the constraint is applied. The constraints override any other
2662 condition set by the transport problem or the boundary condition.
2663
2664 The transport problem is symmetrical if
2665
2666 *A[i,j]=A[j,i]* and *B[j]=C[j]* and *A_reduced[i,j]=A_reduced[j,i]* and
2667 *B_reduced[j]=C_reduced[j]*
2668
2669 For a system and a solution with several components the transport problem
2670 has the form
2671
2672 *(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]*
2673
2674 *A* and *A_reduced* are of rank four, *B*, *B_reduced*, *C* and
2675 *C_reduced* are each of rank three, *M*, *M_reduced*, *D*, *D_reduced*,
2676 *X_reduced* and *X* are each of rank two and *Y* and *Y_reduced* are of
2677 rank one. The natural boundary conditions take the form:
2678
2679 *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*
2680
2681 The coefficient *d* and *m* are of rank two and *y* is of rank one with
2682 all in `FunctionOnBoundary`. The coefficients
2683 *d_reduced* and *m_reduced* are of rank two and *y_reduced* is of rank
2684 one all in `ReducedFunctionOnBoundary`.
2685
2686 Constraints take the form
2687
2688 *u[i]_t=r[i]* where *q[i]>0*
2689
2690 *r* and *q* are each rank one. Notice that at some locations not
2691 necessarily all components must have a constraint.
2692
2693 The transport problem is symmetrical if
2694
2695 - *M[i,k]=M[i,k]*
2696 - *M_reduced[i,k]=M_reduced[i,k]*
2697 - *A[i,j,k,l]=A[k,l,i,j]*
2698 - *A_reduced[i,j,k,l]=A_reduced[k,l,i,j]*
2699 - *B[i,j,k]=C[k,i,j]*
2700 - *B_reduced[i,j,k]=C_reduced[k,i,j]*
2701 - *D[i,k]=D[i,k]*
2702 - *D_reduced[i,k]=D_reduced[i,k]*
2703 - *m[i,k]=m[k,i]*
2704 - *m_reduced[i,k]=m_reduced[k,i]*
2705 - *d[i,k]=d[k,i]*
2706 - *d_reduced[i,k]=d_reduced[k,i]*
2707 - *d_dirac[i,k]=d_dirac[k,i]*
2708
2709 `TransportPDE` also supports solution discontinuities over a contact region
2710 in the domain. To specify the conditions across the discontinuity we are
2711 using the generalised flux *J* which, in the case of a system of PDEs and
2712 several components of the solution, is defined as
2713
2714 *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]*
2715
2716 For the case of single solution component and single PDE *J* is defined as
2717
2718 *J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[j]+(B[i]+B_reduced[i])*u+X[i]+X_reduced[i]*
2719
2720 In the context of discontinuities *n* denotes the normal on the
2721 discontinuity pointing from side 0 towards side 1 calculated from
2722 `FunctionSpace.getNormal` of `FunctionOnContactZero`.
2723 For a system of transport problems the contact condition takes the form
2724
2725 *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]*
2726
2727 where *J0* and *J1* are the fluxes on side 0 and side 1 of the
2728 discontinuity, respectively. *jump(u)*, which is the difference of the
2729 solution at side 1 and at side 0, denotes the jump of *u* across
2730 discontinuity along the normal calculated by `jump`.
2731 The coefficient *d_contact* is of rank two and *y_contact* is of rank one
2732 both in `FunctionOnContactZero` or `FunctionOnContactOne`.
2733 The coefficient *d_contact_reduced* is of rank two and *y_contact_reduced*
2734 is of rank one both in `ReducedFunctionOnContactZero` or `ReducedFunctionOnContactOne`.
2735 In case of a single PDE and a single component solution the contact
2736 condition takes the form
2737
2738 *n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)*
2739
2740 In this case the coefficient *d_contact* and *y_contact* are each scalar
2741 both in `FunctionOnContactZero` or
2742 `FunctionOnContactOne` and the coefficient
2743 *d_contact_reduced* and *y_contact_reduced* are each scalar both in
2744 `ReducedFunctionOnContactZero` or
2745 `ReducedFunctionOnContactOne`.
2746
2747 Typical usage::
2748
2749 p = TransportPDE(dom)
2750 p.setValue(M=1., C=[-1.,0.])
2751 p.setInitialSolution(u=exp(-length(dom.getX()-[0.1,0.1])**2)
2752 t = 0
2753 dt = 0.1
2754 while (t < 1.):
2755 u = p.solve(dt)
2756
2757 """
2758 def __init__(self,domain,numEquations=None,numSolutions=None, useBackwardEuler=None, debug=False):
2759 """
2760 Initializes a transport problem.
2761
2762 :param domain: domain of the PDE
2763 :type domain: `Domain`
2764 :param numEquations: number of equations. If ``None`` the number of
2765 equations is extracted from the coefficients.
2766 :param numSolutions: number of solution components. If ``None`` the number
2767 of solution components is extracted from the
2768 coefficients.
2769 :param debug: if True debug information is printed
2770 """
2771 super(TransportPDE, self).__init__(domain,numEquations,numSolutions,debug)
2772 #
2773 # the coefficients of the transport problem
2774 #
2775 self.introduceCoefficients(
2776 M=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2777 A=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2778 B=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2779 C=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2780 D=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2781 X=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
2782 Y=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2783 m=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2784 d=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2785 y=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2786 d_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2787 y_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2788 M_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2789 A_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2790 B_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2791 C_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2792 D_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2793 X_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
2794 Y_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2795 m_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2796 d_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2797 y_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2798 d_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2799 y_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2800 d_dirac=PDECoef(PDECoef.DIRACDELTA,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2801 y_dirac=PDECoef(PDECoef.DIRACDELTA,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2802 r=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.RIGHTHANDSIDE),
2803 q=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.BOTH) )
2804 if not useBackwardEuler == None:
2805 import warnings
2806 warnings.warn("Argument useBackwardEuler has expired and will be removed in a later release. Please use SolverOptions.setODESolver() instead.", PendingDeprecationWarning, stacklevel=2)
2807 if useBackwardEuler: self.getSolverOptions().setODESolver(SolverOptions.BACKWARD_EULER)
2808
2809 def __str__(self):
2810 """
2811 Returns the string representation of the transport problem.
2812
2813 :return: a simple representation of the transport problem
2814 :rtype: ``str``
2815 """
2816 return "<TransportPDE %d>"%id(self)
2817
2818 def checkSymmetry(self,verbose=True):
2819 """
2820 Tests the transport problem for symmetry.
2821
2822 :param verbose: if set to True or not present a report on coefficients
2823 which break the symmetry is printed.
2824 :type verbose: ``bool``
2825 :return: True if the PDE is symmetric
2826 :rtype: ``bool``
2827 :note: This is a very expensive operation. It should be used for
2828 degugging only! The symmetry flag is not altered.
2829 """
2830 out=True
2831 out=out and self.checkSymmetricTensor("M", verbose)
2832 out=out and self.checkSymmetricTensor("M_reduced", verbose)
2833 out=out and self.checkSymmetricTensor("A", verbose)
2834 out=out and self.checkSymmetricTensor("A_reduced", verbose)
2835 out=out and self.checkReciprocalSymmetry("B","C", verbose)
2836 out=out and self.checkReciprocalSymmetry("B_reduced","C_reduced", verbose)
2837 out=out and self.checkSymmetricTensor("D", verbose)
2838 out=out and self.checkSymmetricTensor("D_reduced", verbose)
2839 out=out and self.checkSymmetricTensor("m", verbose)
2840 out=out and self.checkSymmetricTensor("m_reduced", verbose)
2841 out=out and self.checkSymmetricTensor("d", verbose)
2842 out=out and self.checkSymmetricTensor("d_reduced", verbose)
2843 out=out and self.checkSymmetricTensor("d_contact", verbose)
2844 out=out and self.checkSymmetricTensor("d_contact_reduced", verbose)
2845 out=out and self.checkSymmetricTensor("d_dirac", verbose)
2846 return out
2847
2848 def setValue(self,**coefficients):
2849 """
2850 Sets new values to coefficients.
2851
2852 :param coefficients: new values assigned to coefficients
2853 :keyword M: value for coefficient ``M``
2854 :type M: any type that can be cast to a `Data` object on
2855 `Function`
2856 :keyword M_reduced: value for coefficient ``M_reduced``
2857 :type M_reduced: any type that can be cast to a `Data`
2858 object on `Function`
2859 :keyword A: value for coefficient ``A``
2860 :type A: any type that can be cast to a `Data` object on
2861 `Function`
2862 :keyword A_reduced: value for coefficient ``A_reduced``
2863 :type A_reduced: any type that can be cast to a `Data`
2864 object on `ReducedFunction`
2865 :keyword B: value for coefficient ``B``
2866 :type B: any type that can be cast to a `Data` object on
2867 `Function`
2868 :keyword B_reduced: value for coefficient ``B_reduced``
2869 :type B_reduced: any type that can be cast to a `Data`
2870 object on `ReducedFunction`
2871 :keyword C: value for coefficient ``C``
2872 :type C: any type that can be cast to a `Data` object on
2873 `Function`
2874 :keyword C_reduced: value for coefficient ``C_reduced``
2875 :type C_reduced: any type that can be cast to a `Data`
2876 object on `ReducedFunction`
2877 :keyword D: value for coefficient ``D``
2878 :type D: any type that can be cast to a `Data` object on
2879 `Function`
2880 :keyword D_reduced: value for coefficient ``D_reduced``
2881 :type D_reduced: any type that can be cast to a `Data`
2882 object on `ReducedFunction`
2883 :keyword X: value for coefficient ``X``
2884 :type X: any type that can be cast to a `Data` object on
2885 `Function`
2886 :keyword X_reduced: value for coefficient ``X_reduced``
2887 :type X_reduced: any type that can be cast to a `Data`
2888 object on `ReducedFunction`
2889 :keyword Y: value for coefficient ``Y``
2890 :type Y: any type that can be cast to a `Data` object on
2891 `Function`
2892 :keyword Y_reduced: value for coefficient ``Y_reduced``
2893 :type Y_reduced: any type that can be cast to a `Data`
2894 object on `ReducedFunction`
2895 :keyword m: value for coefficient ``m``
2896 :type m: any type that can be cast to a `Data` object on
2897 `FunctionOnBoundary`
2898 :keyword m_reduced: value for coefficient ``m_reduced``
2899 :type m_reduced: any type that can be cast to a `Data`
2900 object on `FunctionOnBoundary`
2901 :keyword d: value for coefficient ``d``
2902 :type d: any type that can be cast to a `Data` object on
2903 `FunctionOnBoundary`
2904 :keyword d_reduced: value for coefficient ``d_reduced``
2905 :type d_reduced: any type that can be cast to a `Data`
2906 object on `ReducedFunctionOnBoundary`
2907 :keyword y: value for coefficient ``y``
2908 :type y: any type that can be cast to a `Data` object on
2909 `FunctionOnBoundary`
2910 :keyword d_contact: value for coefficient ``d_contact``
2911 :type d_contact: any type that can be cast to a `Data`
2912 object on `FunctionOnContactOne` or `FunctionOnContactZero`
2913 :keyword d_contact_reduced: value for coefficient ``d_contact_reduced``
2914 :type d_contact_reduced: any type that can be cast to a `Data` object on `ReducedFunctionOnContactOne` or `ReducedFunctionOnContactZero`
2915 :keyword y_contact: value for coefficient ``y_contact``
2916 :type y_contact: any type that can be cast to a `Data`
2917 object on `FunctionOnContactOne` or `FunctionOnContactZero`
2918 :keyword y_contact_reduced: value for coefficient ``y_contact_reduced``
2919 :type y_contact_reduced: any type that can be cast to a `Data` object on `ReducedFunctionOnContactOne` or `ReducedFunctionOnContactZero`
2920
2921 :keyword d_dirac: value for coefficient ``d_dirac``
2922 :type d_dirac: any type that can be cast to a `Data` object on `DiracDeltaFunctions`
2923 :keyword y_dirac: value for coefficient ``y_dirac``
2924 :type y_dirac: any type that can be cast to a `Data` object on `DiracDeltaFunctions`
2925
2926 :keyword r: values prescribed to the solution at the locations of constraints
2927 :type r: any type that can be cast to a `Data` object on
2928 `Solution` or `ReducedSolution`
2929 depending on whether reduced order is used for the solution
2930 :keyword q: mask for the location of constraints
2931 :type q: any type that can be cast to a `Data` object on
2932 `Solution` or
2933 `ReducedSolution` depending on whether
2934 reduced order is used for the representation of the equation
2935 :raise IllegalCoefficient: if an unknown coefficient keyword is used
2936 """
2937 super(TransportPDE,self).setValue(**coefficients)
2938
2939 def createOperator(self):
2940 """
2941 Returns an instance of a new transport operator.
2942 """
2943 optype=self.getRequiredOperatorType()
2944 self.trace("New Transport problem pf type %s is allocated."%optype)
2945 return self.getDomain().newTransportProblem( \
2946 self.getNumEquations(), \
2947 self.getFunctionSpaceForSolution(), \
2948 optype)
2949
2950
2951 def getRequiredOperatorType(self):
2952 """
2953 Returns the system type which needs to be used by the current set up.
2954
2955 :return: a code to indicate the type of transport problem scheme used
2956 :rtype: ``float``
2957 """
2958 solver_options=self.getSolverOptions()
2959 return self.getDomain().getTransportTypeId(solver_options.getSolverMethod(), solver_options.getPreconditioner(),solver_options.getPackage(), solver_options.isSymmetric())
2960
2961 def getUnlimitedTimeStepSize(self):
2962 """
2963 Returns the value returned by the ``getSafeTimeStepSize`` method to
2964 indicate no limit on the safe time step size.
2965
2966 :return: the value used to indicate that no limit is set to the time
2967 step size
2968 :rtype: ``float``
2969 :note: Typically the biggest positive float is returned
2970 """
2971 return self.getOperator().getUnlimitedTimeStepSize()
2972
2973 def getSafeTimeStepSize(self):
2974 """
2975 Returns a safe time step size to do the next time step.
2976
2977 :return: safe time step size
2978 :rtype: ``float``
2979 :note: If not ``getSafeTimeStepSize()`` < ``getUnlimitedTimeStepSize()``
2980 any time step size can be used.
2981 """
2982 return self.getOperator().getSafeTimeStepSize()
2983
2984 #====================================================================
2985 def getSolution(self, dt=None, u0=None):
2986 """
2987 Returns the solution by marching forward by time step dt. if ''u0'' is present,
2988 ''u0'' is used as the initial value otherwise the solution from the last call is used.
2989
2990 :param dt: time step size. If ``None`` the last solution is returned.
2991 :type dt: positive ``float`` or ``None``
2992 :param u0: new initial solution or ``None``
2993 :type u0: any object that can be interpolated to a `Data`
2994 object on `Solution` or `ReducedSolution`
2995 :return: the solution
2996 :rtype: `Data`
2997 """
2998 if not dt == None:
2999 option_class=self.getSolverOptions()
3000 if dt<=0:
3001 raise ValueError("step size needs to be positive.")
3002 if u0 == None:
3003 u0=self.getCurrentSolution()
3004 else:
3005 u0=util.interpolate(u0,self.getFunctionSpaceForSolution())
3006 if self.getNumSolutions() == 1:
3007 if u0.getShape()!=():
3008 raise ValueError("Illegal shape %s of initial solution."%(u0.getShape(),))
3009 else:
3010 if u0.getShape()!=(self.getNumSolutions(),):
3011 raise ValueError("Illegal shape %s of initial solution."%(u0.getShape(),))
3012 self.setSolution(self.getOperator().solve(u0, self.getRightHandSide(),dt,option_class))
3013 self.validSolution()
3014 return self.getCurrentSolution()
3015
3016 def setInitialSolution(self,u):
3017 """
3018 Sets the initial solution.
3019
3020 :param u: initial solution
3021 :type u: any object that can be interpolated to a `Data`
3022 object on `Solution` or `ReducedSolution`
3023 """
3024 u2=util.interpolate(u,self.getFunctionSpaceForSolution())
3025 if self.getNumSolutions() == 1:
3026 if u2.getShape()!=():
3027 raise ValueError("Illegal shape %s of initial solution."%(u2.getShape(),))
3028 else:
3029 if u2.getShape()!=(self.getNumSolutions(),):
3030 raise ValueError("Illegal shape %s of initial solution."%(u2.getShape(),))
3031 self.setSolution(u2,validate=False)
3032
3033
3034 def getSystem(self):
3035 """
3036 Returns the operator and right hand side of the PDE.
3037
3038 :return: the discrete version of the PDE
3039 :rtype: ``tuple`` of `Operator` and
3040 `Data`
3041
3042 """
3043 if not self.isOperatorValid() or not self.isRightHandSideValid():
3044 self.resetRightHandSide()
3045 righthandside=self.getCurrentRightHandSide()
3046 self.resetOperator()
3047 operator=self.getCurrentOperator()
3048 self.getDomain().addPDEToTransportProblem(
3049 operator,
3050 righthandside,
3051 self.getCoefficient("M"),
3052 self.getCoefficient("A"),
3053 self.getCoefficient("B"),
3054 self.getCoefficient("C"),
3055 self.getCoefficient("D"),
3056 self.getCoefficient("X"),
3057 self.getCoefficient("Y"),
3058 self.getCoefficient("d"),
3059 self.getCoefficient("y"),
3060 self.getCoefficient("d_contact"),
3061 self.getCoefficient("y_contact"),
3062 self.getCoefficient("d_dirac"),
3063 self.getCoefficient("y_dirac") )
3064 self.getDomain().addPDEToTransportProblem(
3065 operator,
3066 righthandside,
3067 self.getCoefficient("M_reduced"),
3068 self.getCoefficient("A_reduced"),
3069 self.getCoefficient("B_reduced"),
3070 self.getCoefficient("C_reduced"),
3071 self.getCoefficient("D_reduced"),
3072 self.getCoefficient("X_reduced"),
3073 self.getCoefficient("Y_reduced"),
3074 self.getCoefficient("d_reduced"),
3075 self.getCoefficient("y_reduced"),
3076 self.getCoefficient("d_contact_reduced"),
3077 self.getCoefficient("y_contact_reduced"),
3078 escore.Data(),
3079 escore.Data() )
3080 operator.insertConstraint(righthandside,self.getCoefficient("q"),self.getCoefficient("r"))
3081 self.trace("New system has been built.")
3082 self.validOperator()
3083 self.validRightHandSide()
3084 self.setSystemStatus()
3085 self.trace("System status is %s."%self.getSystemStatus())
3086 return (self.getCurrentOperator(), self.getCurrentRightHandSide())
3087
3088 def setDebug(self, flag):
3089 """
3090 Switches debug output on if ``flag`` is True,
3091 otherwise it is switched off.
3092
3093 :param flag: desired debug status
3094 :type flag: ``bool``
3095 """
3096 if flag:
3097 self.setDebugOn()
3098 else:
3099 self.setDebugOff()
3100
3101 def setDebugOn(self):
3102 """
3103 Switches debug output on.
3104 """
3105 super(TransportPDE,self).setDebugOn()
3106
3107 def setDebugOff(self):
3108 """
3109 Switches debug output off.
3110 """
3111 super(TransportPDE,self).setDebugOff()
3112
3113 def SingleTransportPDE(domain, debug=False):
3114 """
3115 Defines a single transport problem
3116
3117 :param domain: domain of the PDE
3118 :type domain: `Domain`
3119 :param debug: if True debug information is printed
3120 :rtype: `TransportPDE`
3121 """
3122 return TransportPDE(domain,numEquations=1,numSolutions=1, debug=debug)

Properties

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

  ViewVC Help
Powered by ViewVC 1.1.26