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

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

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

trunk/esys2/escript/py_src/linearPDEs.py revision 147 by jgs, Fri Aug 12 01:45:47 2005 UTC trunk/escript/py_src/linearPDEs.py revision 3402 by gross, Tue Dec 7 07:36:12 2010 UTC
# Line 1  Line 1 
1  # $Id$  # -*- coding: utf-8 -*-
2    
3  ## @file linearPDEs.py  ########################################################
4    #
5    # Copyright (c) 2003-2010 by University of Queensland
6    # Earth Systems Science Computational Center (ESSCC)
7    # http://www.uq.edu.au/esscc
8    #
9    # Primary Business: Queensland, Australia
10    # Licensed under the Open Software License version 3.0
11    # http://www.opensource.org/licenses/osl-3.0.php
12    #
13    ########################################################
14    
15    __copyright__="""Copyright (c) 2003-2010 by University of Queensland
16    Earth Systems Science Computational Center (ESSCC)
17    http://www.uq.edu.au/esscc
18    Primary Business: Queensland, Australia"""
19    __license__="""Licensed under the Open Software License version 3.0
20    http://www.opensource.org/licenses/osl-3.0.php"""
21    __url__="https://launchpad.net/escript-finley"
22    
23  """  """
24  Functions and classes for linear PDEs  The module provides an interface to define and solve linear partial
25    differential equations (PDEs) and Transport problems within `escript`.
26    `linearPDEs` does not provide any solver capabilities in itself but hands the
27    PDE over to the PDE solver library defined through the `Domain`
28    of the PDE. The general interface is provided through the `LinearPDE` class.
29    `TransportProblem` provides an interface to initial value problems dominated
30    by its advective terms.
31    
32    :var __author__: name of author
33    :var __copyright__: copyrights
34    :var __license__: licence agreement
35    :var __url__: url entry point on documentation
36    :var __version__: version
37    :var __date__: date of the version
38  """  """
39    
40    import math
41  import escript  import escript
42  import util  import util
43  import numarray  import numpy
44    
45    __author__="Lutz Gross, l.gross@uq.edu.au"
46    
47    
48    class SolverOptions(object):
49        """
50        this class defines the solver options for a linear or non-linear solver.
51        
52        The option also supports the handling of diagnostic informations.
53        
54        Typical usage is
55        
56        ::
57        
58          opts=SolverOptions()
59          print opts
60          opts.resetDiagnostics()
61          u=solver(opts)
62          print "number of iteration steps: =",opts.getDiagnostics("num_iter")
63    
64        :cvar DEFAULT: The default method used to solve the system of linear equations
65        :cvar DIRECT: The direct solver based on LDU factorization
66        :cvar CHOLEVSKY: The direct solver based on LDLt factorization (can only be applied for symmetric PDEs)
67        :cvar PCG: The preconditioned conjugate gradient method (can only be applied for symmetric PDEs)
68        :cvar CR: The conjugate residual method
69        :cvar CGS: The conjugate gradient square method
70        :cvar BICGSTAB: The stabilized Bi-Conjugate Gradient method
71        :cvar TFQMR: Transpose Free Quasi Minimal Residual method
72        :cvar MINRES: Minimum residual method
73        :cvar ILU0: The incomplete LU factorization preconditioner with no fill-in
74        :cvar ILUT: The incomplete LU factorization preconditioner with fill-in
75        :cvar JACOBI: The Jacobi preconditioner
76        :cvar GMRES: The Gram-Schmidt minimum residual method
77        :cvar PRES20: Special GMRES with restart after 20 steps and truncation after 5 residuals
78        :cvar ROWSUM_LUMPING: Matrix lumping using row sum
79        :cvar HRZ_LUMPING: Matrix lumping using the HRZ approach
80        :cvar NO_REORDERING: No matrix reordering allowed
81        :cvar MINIMUM_FILL_IN: Reorder matrix to reduce fill-in during factorization
82        :cvar NESTED_DISSECTION: Reorder matrix to improve load balancing during factorization
83        :cvar PASO: PASO solver package
84        :cvar SCSL: SGI SCSL solver library
85        :cvar MKL: Intel's MKL solver library
86        :cvar UMFPACK: The UMFPACK library
87        :cvar TRILINOS: The TRILINOS parallel solver class library from Sandia National Labs
88        :cvar ITERATIVE: The default iterative solver
89        :cvar AMG: Algebraic Multi Grid
90        :cvar AMLI: Algebraic Multi Level Iteration
91        :cvar REC_ILU: recursive ILU0
92        :cvar RILU: relaxed ILU0
93        :cvar GAUSS_SEIDEL: Gauss-Seidel preconditioner
94        :cvar DEFAULT_REORDERING: the reordering method recommended by the solver
95        :cvar SUPER_LU: the Super_LU solver package
96        :cvar PASTIX: the Pastix direct solver_package
97        :cvar YAIR_SHAPIRA_COARSENING: AMG and AMLI coarsening method by Yair-Shapira
98        :cvar RUGE_STUEBEN_COARSENING: AMG and AMLI coarsening method by Ruge and Stueben
99        :cvar AGGREGATION_COARSENING: AMG and AMLI coarsening using (symmetric) aggregation
100        :cvar STANDARD_COARSENING: AMG and AMLI standard coarsening using mesure of importance of the unknowns
101        :cvar MIN_COARSE_MATRIX_SIZE: minimum size of the coarsest level matrix to use direct solver.
102        :cvar NO_PRECONDITIONER: no preconditioner is applied.
103        """
104        DEFAULT= 0
105        DIRECT= 1
106        CHOLEVSKY= 2
107        PCG= 3
108        CR= 4
109        CGS= 5
110        BICGSTAB= 6
111        ILU0= 8
112        ILUT= 9
113        JACOBI= 10
114        GMRES= 11
115        PRES20= 12
116        LUMPING= 13
117        ROWSUM_LUMPING= 13
118        HRZ_LUMPING= 14
119        NO_REORDERING= 17
120        MINIMUM_FILL_IN= 18
121        NESTED_DISSECTION= 19
122        MKL= 15
123        UMFPACK= 16
124        ITERATIVE= 20
125        PASO= 21
126        AMG= 22
127        REC_ILU = 23
128        TRILINOS = 24
129        NONLINEAR_GMRES = 25
130        TFQMR = 26
131        MINRES = 27
132        GAUSS_SEIDEL=28
133        RILU=29
134        DEFAULT_REORDERING=30
135        SUPER_LU=31
136        PASTIX=32
137        YAIR_SHAPIRA_COARSENING=33
138        RUGE_STUEBEN_COARSENING=34
139        AGGREGATION_COARSENING=35
140        NO_PRECONDITIONER=36
141        MIN_COARSE_MATRIX_SIZE=37
142        AMLI=38
143        STANDARD_COARSENING=39
144    
145        def __init__(self):
146            self.setLevelMax()
147            self.setCoarseningThreshold()
148            self.setSmoother()
149            self.setNumSweeps()
150            self.setNumPreSweeps()
151            self.setNumPostSweeps()
152            self.setTolerance()
153            self.setAbsoluteTolerance()
154            self.setInnerTolerance()
155            self.setDropTolerance()
156            self.setDropStorage()
157            self.setIterMax()
158            self.setInnerIterMax()
159            self.setTruncation()
160            self.setRestart()
161            self.setSymmetry()
162            self.setVerbosity()
163            self.setInnerToleranceAdaption()
164            self.setAcceptanceConvergenceFailure()
165            self.setReordering()
166            self.setPackage()
167            self.setSolverMethod()
168            self.setPreconditioner()
169            self.setCoarsening()
170            self.setMinCoarseMatrixSize()
171            self.setRelaxationFactor()        
172            self.setLocalPreconditionerOff()
173            self.resetDiagnostics(all=True)
174            self.setMinCoarseMatrixSparsity()
175            self.setNumRefinements()
176            self.setNumCoarseMatrixRefinements()
177            self.setUsePanel()
178            self.setUseDirectInterpolation()
179            self.setDiagonalDominanceThreshold()
180            
181    
182        def __str__(self):
183            return self.getSummary()
184        def getSummary(self):
185            """
186            Returns a string reporting the current settings
187            """
188            out="Solver Package: %s"%(self.getName(self.getPackage()))
189            out+="\nVerbosity = %s"%self.isVerbose()
190            out+="\nAccept failed convergence = %s"%self.acceptConvergenceFailure()
191            out+="\nRelative tolerance = %e"%self.getTolerance()
192            out+="\nAbsolute tolerance = %e"%self.getAbsoluteTolerance()
193            out+="\nSymmetric problem = %s"%self.isSymmetric()
194            out+="\nMaximum number of iteration steps = %s"%self.getIterMax()
195            # out+="\nInner tolerance = %e"%self.getInnerTolerance()
196            # out+="\nAdapt innner tolerance = %s"%self.adaptInnerTolerance()
197        
198            if self.getPackage() == self.PASO:
199                out+="\nSolver method = %s"%self.getName(self.getSolverMethod())
200                if self.getSolverMethod() == self.GMRES:
201                    out+="\nTruncation  = %s"%self.getTruncation()
202                    out+="\nRestart  = %s"%self.getRestart()
203                if self.getSolverMethod() == self.AMG:
204                    out+="\nNumber of pre / post sweeps = %s / %s, %s"%(self.getNumPreSweeps(), self.getNumPostSweeps(), self.getNumSweeps())
205                    out+="\nMaximum number of levels = %s"%self.getLevelMax()
206                    out+="\nCoarsening threshold = %e"%self.getCoarseningThreshold()
207                    out+="\nCoarsening method = %s"%self.getName(self.getCoarsening())
208                out+="\nPreconditioner = %s"%self.getName(self.getPreconditioner())
209                out+="\nApply preconditioner locally = %s"%self.useLocalPreconditioner()
210                if self.getPreconditioner() == self.AMG:
211                    out+="\nMaximum number of levels = %s"%self.getLevelMax()
212                    out+="\nCoarsening threshold = %e"%self.getCoarseningThreshold()
213                    out+="\nMinimal sparsity on coarsest level = %e"%(self.getMinCoarseMatrixSparsity(),)
214                    out+="\nSmoother = %s"%self.getName(self.getSmoother())
215                    out+="\nMinimum size of the coarsest level matrix = %e"%self.getCoarseningThreshold()
216                    out+="\nNumber of pre / post sweeps = %s / %s, %s"%(self.getNumPreSweeps(), self.getNumPostSweeps(), self.getNumSweeps())
217                    out+="\nNumber of refinement steps in coarsest level solver = %d"%(self.getNumCoarseMatrixRefinements(),)
218                    out+="\nUse node panel = %s"%(self.usePanel())
219                    out+="\nUse direct interpolation only = %s"%(self.useDirectInterpolation())
220                    out+="\nThreshold for diagonal dominant rows = %s"%(setDiagonalDominanceThreshold(),)
221    
222                if self.getPreconditioner() == self.AMLI:
223                    out+="\nMaximum number of levels = %s"%self.getLevelMax()
224                    out+="\nCoarsening method = %s"%self.getName(self.getCoarsening())
225                    out+="\nCoarsening threshold = %e"%self.getMinCoarseMatrixSize()
226                    out+="\nMinimum size of the coarsest level matrix = %e"%self.getCoarseningThreshold()
227                    out+="\nNumber of pre / post sweeps = %s / %s, %s"%(self.getNumPreSweeps(), self.getNumPostSweeps(), self.getNumSweeps())
228            if self.getPreconditioner() == self.GAUSS_SEIDEL:
229                    out+="\nNumber of sweeps = %s"%self.getNumSweeps()
230                if self.getPreconditioner() == self.ILUT:
231                    out+="\nDrop tolerance = %e"%self.getDropTolerance()
232                    out+="\nStorage increase = %e"%self.getDropStorage()
233                if self.getPreconditioner() == self.RILU:
234                    out+="\nRelaxation factor = %e"%self.getRelaxationFactor()
235            return out
236            
237        def getName(self,key):
238            """
239            returns the name of a given key
240            
241            :param key: a valid key
242            """
243            if key == self.DEFAULT: return "DEFAULT"
244            if key == self.DIRECT: return "DIRECT"
245            if key == self.CHOLEVSKY: return "CHOLEVSKY"
246            if key == self.PCG: return "PCG"
247            if key == self.CR: return "CR"
248            if key == self.CGS: return "CGS"
249            if key == self.BICGSTAB: return "BICGSTAB"
250            if key == self.ILU0: return "ILU0:"
251            if key == self.ILUT: return "ILUT"
252            if key == self.JACOBI: return "JACOBI"
253            if key == self.GMRES: return "GMRES"
254            if key == self.PRES20: return "PRES20"
255            if key == self.ROWSUM_LUMPING: return "ROWSUM_LUMPING"
256            if key == self.HRZ_LUMPING: return "HRZ_LUMPING"
257            if key == self.NO_REORDERING: return "NO_REORDERING"
258            if key == self.MINIMUM_FILL_IN: return "MINIMUM_FILL_IN"
259            if key == self.NESTED_DISSECTION: return "NESTED_DISSECTION"
260            if key == self.MKL: return "MKL"
261            if key == self.UMFPACK: return "UMFPACK"
262            if key == self.ITERATIVE: return "ITERATIVE"
263            if key == self.PASO: return "PASO"
264            if key == self.AMG: return "AMG"
265            if key == self.AMLI: return "AMLI"
266            if key == self.REC_ILU: return "REC_ILU"
267            if key == self.TRILINOS: return "TRILINOS"
268            if key == self.NONLINEAR_GMRES: return "NONLINEAR_GMRES"
269            if key == self.TFQMR: return "TFQMR"
270            if key == self.MINRES: return "MINRES"
271            if key == self.GAUSS_SEIDEL: return "GAUSS_SEIDEL"
272            if key == self.RILU: return "RILU"
273            if key == self.DEFAULT_REORDERING: return "DEFAULT_REORDERING"
274            if key == self.SUPER_LU: return "SUPER_LU"
275            if key == self.PASTIX: return "PASTIX"
276            if key == self.YAIR_SHAPIRA_COARSENING: return "YAIR_SHAPIRA_COARSENING"
277            if key == self.RUGE_STUEBEN_COARSENING: return "RUGE_STUEBEN_COARSENING"
278            if key == self.STANDARD_COARSENING: return "STANDARD_COARSENING"
279            if key == self.AGGREGATION_COARSENING: return "AGGREGATION_COARSENING"
280            if key == self.NO_PRECONDITIONER: return "NO_PRECONDITIONER"
281            
282        def resetDiagnostics(self,all=False):
283            """
284            resets the diagnostics
285            
286            :param all: if ``all`` is ``True`` all diagnostics including accumulative counters are reset.
287            :type all: ``bool``
288            """
289            self.__num_iter=None
290            self.__num_level=None
291            self.__num_inner_iter=None
292            self.__time=None
293            self.__set_up_time=None
294            self.__net_time=None
295            self.__residual_norm=None
296            self.__converged=None
297            self.__preconditioner_size=-1
298            self.__time_step_backtracking_used=None
299            if all:
300                self.__cum_num_inner_iter=0
301                self.__cum_num_iter=0
302                self.__cum_time=0
303                self.__cum_set_up_time=0
304                self.__cum_net_time=0
305    
306        def _updateDiagnostics(self, name, value):
307            """
308            Updates diagnostic information
309            
310            :param name: name of  diagnostic information
311            :type name: ``str`` in the list "num_iter", "num_level", "num_inner_iter", "time", "set_up_time", "net_time", "residual_norm", "converged".
312            :param value: new value of the diagnostic information
313            :note: this function is used by a solver to report diagnostics informations.
314            """
315            if name == "num_iter":
316                self.__num_iter=int(value)
317                self.__cum_num_iter+=self.__num_iter
318            if name == "num_level":
319                self.__num_level=int(value)
320            if name == "num_inner_iter":
321                self.__num_inner_iter=int(value)
322                self.__cum_num_inner_iter+=self.__num_inner_iter
323            if name == "time":
324                self.__time=float(value)
325                self.__cum_time+=self.__time
326            if name == "set_up_time":
327                self.__set_up_time=float(value)
328                self.__cum_set_up_time+=self.__set_up_time
329            if name == "net_time":
330                self.__net_time=float(value)
331                self.__cum_net_time+=self.__net_time
332            if name == "residual_norm":
333                self.__residual_norm=float(value)
334            if name == "converged":
335                self.__converged = (value == True)
336            if name == "time_step_backtracking_used":
337                self.__time_step_backtracking_used = (value == True)
338            if name == "coarse_level_sparsity":
339               self.__coarse_level_sparsity = value
340        if name == "num_coarse_unknowns":
341               self.__num_coarse_unknowns= value
342        def getDiagnostics(self, name):
343            """
344            Returns the diagnostic information ``name``. Possible values are:
345            
346        - "num_iter": the number of iteration steps
347            - "cum_num_iter": the cumulative number of iteration steps
348            - "num_level": the number of level in multi level solver
349            - "num_inner_iter": the number of inner iteration steps
350            - "cum_num_inner_iter": the cumulative number of inner iteration steps
351            - "time": execution time
352            - "cum_time": cumulative execution time
353            - "set_up_time": time to set up of the solver, typically this includes factorization and reordering
354            - "cum_set_up_time": cumulative time to set up of the solver
355            - "net_time": net execution time, excluding setup time for the solver and execution time for preconditioner
356            - "cum_net_time": cumulative net execution time
357            - "preconditioner_size": size of preconditioner [Bytes]
358            - "converged": return True if solution has converged.
359            - "time_step_backtracking_used": returns True if time step back tracking has been used.
360            - "coarse_level_sparsity": returns the sparsity of the matrix on the coarsest level
361            - "num_coarse_unknowns": returns the number of unknowns on the coarsest level
362        
363            
364            :param name: name of diagnostic information to return
365            :type name: ``str`` in the list above.
366            :return: requested value. ``None`` is returned if the value is undefined.
367            :note: If the solver has thrown an exception diagnostic values have an undefined status.
368            """
369            if name == "num_iter": return self.__num_iter
370            elif name == "cum_num_iter": return self.__cum_num_iter
371            elif name == "num_level": return self.__num_level
372            elif name == "num_inner_iter": return self.__num_inner_iter
373            elif name == "cum_num_inner_iter": return self.__cum_num_inner_iter
374            elif name == "time": return self.__time
375            elif name == "cum_time": return self.__cum_time
376            elif name == "set_up_time": return self.__set_up_time
377            elif name == "cum_set_up_time": return self.__cum_set_up_time
378            elif name == "net_time": return self.__net_time
379            elif name == "cum_net_time": return self.__cum_net_time
380            elif name == "residual_norm": return self.__residual_norm
381            elif name == "converged": return self.__converged      
382            elif name == "preconditioner_size": return  self.__preconditioner_size
383            elif name == "time_step_backtracking_used": return  self.__time_step_backtracking_used
384            elif name == "coarse_level_sparsity": return self.__coarse_level_sparsity
385            elif name == "num_coarse_unknowns": return self.__num_coarse_unknowns
386            else:
387                raise ValueError,"unknown diagnostic item %s"%name
388        def hasConverged(self):
389            """
390            Returns ``True`` if the last solver call has been finalized successfully.
391            :note: if an exception has been thrown by the solver the status of this flag is undefined.
392            """
393            return self.getDiagnostics("converged")
394        def setCoarsening(self,method=0):
395            """
396            Sets the key of the coarsening method to be applied in AMG or AMLI
397    
398            :param method: selects the coarsening method .
399            :type method: in {SolverOptions.DEFAULT}, `SolverOptions.YAIR_SHAPIRA_COARSENING`,  `SolverOptions.RUGE_STUEBEN_COARSENING`, `SolverOptions.AGGREGATION_COARSENING`
400            """
401        if method==None: method=0
402            if not method in [self.DEFAULT, self.YAIR_SHAPIRA_COARSENING, self.RUGE_STUEBEN_COARSENING, self.AGGREGATION_COARSENING, self.STANDARD_COARSENING,]:
403                 raise ValueError,"unknown coarsening method %s"%method
404            self.__coarsening=method
405        
406        def getCoarsening(self):
407            """
408            Returns the key of the coarsening algorithm to be applied AMG or AMLI
409    
410            :rtype: in the list `SolverOptions.DEFAULT`, `SolverOptions.YAIR_SHAPIRA_COARSENING`, `SolverOptions.RUGE_STUEBEN_COARSENING`, `SolverOptions.AGGREGATION_COARSENING`
411            """
412            return self.__coarsening
413          
414        def setMinCoarseMatrixSize(self,size=500):
415            """
416            Sets the minumum size of the coarsest level matrix in AMG or AMLI
417    
418            :param size: minumum size of the coarsest level matrix .
419            :type size: positive ``int`` or ``None``
420            """
421        if size==None: size=500
422            size=int(size)
423            if size<0:
424               raise ValueError,"minumum size of the coarsest level matrix must be non-negative."
425            self.__MinCoarseMatrixSize=size
426            
427        def getMinCoarseMatrixSize(self):
428            """
429            Returns the minumum size of the coarsest level matrix in AMG or AMLI
430    
431            :rtype: ``int``
432            """
433            return self.__MinCoarseMatrixSize
434          
435        def setPreconditioner(self, preconditioner=10):
436            """
437            Sets the preconditioner to be used.
438    
439            :param preconditioner: key of the preconditioner to be used.
440            :type preconditioner: in `SolverOptions.ILU0`, `SolverOptions.ILUT`, `SolverOptions.JACOBI`,
441                                        `SolverOptions.AMG`, `SolverOptions.AMLI`, `SolverOptions.REC_ILU`, `SolverOptions.GAUSS_SEIDEL`, `SolverOptions.RILU`,
442                                        `SolverOptions.NO_PRECONDITIONER`
443            :note: Not all packages support all preconditioner. It can be assumed that a package makes a reasonable choice if it encounters an unknown preconditioner.
444            """
445        if preconditioner==None: preconditioner=10
446            if not preconditioner in [  SolverOptions.ILU0, SolverOptions.ILUT, SolverOptions.JACOBI,
447                                        SolverOptions.AMG, SolverOptions.AMLI, SolverOptions.REC_ILU, SolverOptions.GAUSS_SEIDEL, SolverOptions.RILU,
448                                        SolverOptions.NO_PRECONDITIONER] :
449                 raise ValueError,"unknown preconditioner %s"%preconditioner
450            self.__preconditioner=preconditioner    
451        def getPreconditioner(self):
452            """
453            Returns key of the preconditioner to be used.
454    
455            :rtype: in the list `SolverOptions.ILU0`, `SolverOptions.ILUT`, `SolverOptions.JACOBI`, SolverOptions.AMLI,
456                                        `SolverOptions.AMG`, `SolverOptions.REC_ILU`, `SolverOptions.GAUSS_SEIDEL`, `SolverOptions.RILU`,
457                                        `SolverOptions.NO_PRECONDITIONER`
458            """
459            return self.__preconditioner
460        def setSmoother(self, smoother=None):
461            """
462            Sets the smoother to be used.
463    
464            :param smoother: key of the smoother to be used.
465            :type smoother: in `SolverOptions.JACOBI`, `SolverOptions.GAUSS_SEIDEL`
466            :note: Not all packages support all smoothers. It can be assumed that a package makes a reasonable choice if it encounters an unknown smoother.
467            """
468        if smoother==None: smoother=28
469            if not smoother in [ SolverOptions.JACOBI, SolverOptions.GAUSS_SEIDEL ] :
470                 raise ValueError,"unknown smoother %s"%smoother
471            self.__smoother=smoother    
472        def getSmoother(self):
473            """
474            Returns key of the smoother to be used.
475    
476            :rtype: in the list `SolverOptions.JACOBI`, `SolverOptions.GAUSS_SEIDEL`
477            """
478            return self.__smoother  
479        def setSolverMethod(self, method=0):
480            """
481            Sets the solver method to be used. Use ``method``=``SolverOptions.DIRECT`` to indicate that a direct rather than an iterative
482            solver should be used and Use ``method``=``SolverOptions.ITERATIVE`` to indicate that an iterative rather than a direct
483            solver should be used.
484    
485            :param method: key of the solver method to be used.
486            :type method: in `SolverOptions.DEFAULT`, `SolverOptions.DIRECT`, `SolverOptions.CHOLEVSKY`, `SolverOptions.PCG`,
487                            `SolverOptions.CR`, `SolverOptions.CGS`, `SolverOptions.BICGSTAB`,
488                            `SolverOptions.GMRES`, `SolverOptions.PRES20`, `SolverOptions.ROWSUM_LUMPING`, `SolverOptions.HRZ_LUMPING`, `SolverOptions.ITERATIVE`,
489                            `SolverOptions.NONLINEAR_GMRES`, `SolverOptions.TFQMR`, `SolverOptions.MINRES`
490            :note: Not all packages support all solvers. It can be assumed that a package makes a reasonable choice if it encounters an unknown solver method.
491            """
492        if method==None: method=0
493            if not method in [ SolverOptions.DEFAULT, SolverOptions.DIRECT, SolverOptions.CHOLEVSKY, SolverOptions.PCG,
494                               SolverOptions.CR, SolverOptions.CGS, SolverOptions.BICGSTAB,
495                               SolverOptions.GMRES, SolverOptions.PRES20, SolverOptions.ROWSUM_LUMPING, SolverOptions.HRZ_LUMPING,
496                               SolverOptions.ITERATIVE,
497                               SolverOptions.NONLINEAR_GMRES, SolverOptions.TFQMR, SolverOptions.MINRES ]:
498                 raise ValueError,"unknown solver method %s"%method
499            self.__method=method
500        def getSolverMethod(self):
501            """
502            Returns key of the solver method to be used.
503    
504            :rtype: in the list `SolverOptions.DEFAULT`, `SolverOptions.DIRECT`, `SolverOptions.CHOLEVSKY`, `SolverOptions.PCG`,
505                            `SolverOptions.CR`, `SolverOptions.CGS`, `SolverOptions.BICGSTAB`,
506                            `SolverOptions.GMRES`, `SolverOptions.PRES20`, `SolverOptions.ROWSUM_LUMPING`, `SolverOptions.HRZ_LUMPING`, `SolverOptions.ITERATIVE`,
507                            `SolverOptions.NONLINEAR_GMRES`, `SolverOptions.TFQMR`, `SolverOptions.MINRES`
508            """
509            return self.__method
510            
511        def setPackage(self, package=0):
512            """
513            Sets the solver package to be used as a solver.  
514    
515            :param package: key of the solver package to be used.
516            :type package: in `SolverOptions.DEFAULT`, `SolverOptions.PASO`, `SolverOptions.SUPER_LU`, `SolverOptions.PASTIX`, `SolverOptions.MKL`, `SolverOptions.UMFPACK`, `SolverOptions.TRILINOS`
517            :note: Not all packages are support on all implementation. An exception may be thrown on some platforms if a particular is requested.
518            """
519        if package==None: package=0
520            if not package in [SolverOptions.DEFAULT, SolverOptions.PASO, SolverOptions.SUPER_LU, SolverOptions.PASTIX, SolverOptions.MKL, SolverOptions.UMFPACK, SolverOptions.TRILINOS]:
521                 raise ValueError,"unknown solver package %s"%package
522            self.__package=package
523        def getPackage(self):
524            """
525            Returns the solver package key
526    
527            :rtype: in the list `SolverOptions.DEFAULT`, `SolverOptions.PASO`, `SolverOptions.SUPER_LU`, `SolverOptions.PASTIX`, `SolverOptions.MKL`, `SolverOptions.UMFPACK`, `SolverOptions.TRILINOS`
528            """
529            return self.__package
530        def setReordering(self,ordering=30):
531            """
532            Sets the key of the reordering method to be applied if supported by the solver. Some direct solvers support reordering
533            to optimize compute time and storage use during elimination.
534    
535            :param ordering: selects the reordering strategy.
536            :type ordering: in 'SolverOptions.NO_REORDERING', 'SolverOptions.MINIMUM_FILL_IN', 'SolverOptions.NESTED_DISSECTION', 'SolverOptions.DEFAULT_REORDERING'
537            """
538            if not ordering in [self.NO_REORDERING, self.MINIMUM_FILL_IN, self.NESTED_DISSECTION, self.DEFAULT_REORDERING]:
539                 raise ValueError,"unknown reordering strategy %s"%ordering
540            self.__reordering=ordering
541        def getReordering(self):
542            """
543            Returns the key of the reordering method to be applied if supported by the solver.
544    
545            :rtype: ordering in 'SolverOptions.NO_REORDERING', 'SolverOptions.MINIMUM_FILL_IN', 'SolverOptions.NESTED_DISSECTION', 'SolverOptions.DEFAULT_REORDERING'
546            """
547            return self.__reordering
548        def setRestart(self,restart=None):
549            """
550            Sets the number of iterations steps after which GMRES is performing a restart.
551    
552            :param restart: number of iteration steps after which to perform a restart. If equal to ``None`` no restart is performed.
553            :type restart: ``int`` or ``None``
554            """
555            if restart == None:
556                self.__restart=restart
557            else:
558                restart=int(restart)
559                if restart<1:
560                    raise ValueError,"restart must be positive."
561                self.__restart=restart
562            
563        def getRestart(self):
564            """
565            Returns the number of iterations steps after which GMRES is performing a restart.
566            If ``None`` is returned no restart is performed.
567    
568            :rtype: ``int`` or ``None``
569            """
570            if self.__restart < 0:
571                return None
572            else:
573                return self.__restart
574        def _getRestartForC(self):
575            r=self.getRestart()
576            if r == None:
577                return -1
578                else:
579                return r
580      
581        def setDiagonalDominanceThreshold(self,value=0.5):
582            """
583            Sets the threshold for diagonal dominant rows which are eliminated during AMG coarsening.
584    
585         :param value: threshold
586         :type value: ``float``
587            """
588            value=float(value)
589            if value<0 or value>1.:
590           raise ValueError,"Diagonal dominance threshold must be between 0 and 1."
591        self.__diagonal_dominance_threshold=value
592            
593        def getDiagonalDominanceThreshold(self):
594            """
595            Returns the threshold for diagonal dominant rows which are eliminated during AMG coarsening.
596    
597            :rtype: ``float``
598            """
599            return self.__diagonal_dominance_threshold
600            
601        def setTruncation(self,truncation=20):
602            """
603            Sets the number of residuals in GMRES to be stored for orthogonalization.  The more residuals are stored
604            the faster GMRES converged but
605    
606            :param truncation: truncation
607            :type truncation: ``int``
608            """
609            truncation=int(truncation)
610            if truncation<1:
611               raise ValueError,"truncation must be positive."
612            self.__truncation=truncation
613        def getTruncation(self):
614            """
615            Returns the number of residuals in GMRES to be stored for orthogonalization
616    
617            :rtype: ``int``
618            """
619            return self.__truncation
620        def setInnerIterMax(self,iter_max=10):
621            """
622            Sets the maximum number of iteration steps for the inner iteration.
623    
624            :param iter_max: maximum number of inner iterations
625            :type iter_max: ``int``
626            """
627            iter_max=int(iter_max)
628            if iter_max<1:
629               raise ValueError,"maximum number of inner iteration must be positive."
630            self.__inner_iter_max=iter_max
631        def getInnerIterMax(self):
632            """
633            Returns maximum number of inner iteration steps
634    
635            :rtype: ``int``
636            """
637            return self.__inner_iter_max
638        def setIterMax(self,iter_max=100000):
639            """
640            Sets the maximum number of iteration steps
641    
642            :param iter_max: maximum number of iteration steps
643            :type iter_max: ``int``
644            """
645            iter_max=int(iter_max)
646            if iter_max<1:
647               raise ValueError,"maximum number of iteration steps must be positive."
648            self.__iter_max=iter_max
649        def getIterMax(self):
650            """
651            Returns maximum number of iteration steps
652    
653            :rtype: ``int``
654            """
655            return self.__iter_max
656        def setLevelMax(self,level_max=100):
657            """
658            Sets the maximum number of coarsening levels to be used in an algebraic multi level solver or preconditioner
659    
660            :param level_max: maximum number of levels
661            :type level_max: ``int``
662            """
663            level_max=int(level_max)
664            if level_max<0:
665               raise ValueError,"maximum number of coarsening levels must be non-negative."
666            self.__level_max=level_max
667        def getLevelMax(self):
668            """
669            Returns the maximum number of coarsening levels to be used in an algebraic multi level solver or preconditioner
670    
671            :rtype: ``int``
672            """
673            return self.__level_max
674        def setCoarseningThreshold(self,theta=0.25):
675            """
676            Sets the threshold for coarsening in the algebraic multi level solver or preconditioner
677    
678            :param theta: threshold for coarsening
679            :type theta: positive ``float``
680            """
681            theta=float(theta)
682            if theta<0 or theta>1:
683               raise ValueError,"threshold must be non-negative and less or equal 1."
684            self.__coarsening_threshold=theta
685        def getCoarseningThreshold(self):
686            """
687            Returns the threshold for coarsening in the algebraic multi level solver or preconditioner
688    
689            :rtype: ``float``
690            """
691            return self.__coarsening_threshold
692        def setNumSweeps(self,sweeps=1):
693            """
694            Sets the number of sweeps in a Jacobi or Gauss-Seidel/SOR preconditioner.
695    
696            :param sweeps: number of sweeps
697            :type sweeps: positive ``int``
698            """
699            sweeps=int(sweeps)
700            if sweeps<1:
701               raise ValueError,"number of sweeps must be positive."
702            self.__sweeps=sweeps
703        def getNumSweeps(self):
704            """
705            Returns the number of sweeps in a Jacobi or Gauss-Seidel/SOR preconditioner.
706    
707            :rtype: ``int``
708            """
709            return self.__sweeps
710        def setNumPreSweeps(self,sweeps=2):
711            """
712            Sets the number of sweeps in the pre-smoothing step of a multi level solver or preconditioner
713    
714            :param sweeps: number of sweeps
715            :type sweeps: positive ``int``
716            """
717            sweeps=int(sweeps)
718            if sweeps<1:
719               raise ValueError,"number of sweeps must be positive."
720            self.__pre_sweeps=sweeps
721        def getNumPreSweeps(self):
722            """
723            Returns he number of sweeps in the pre-smoothing step of a multi level solver or preconditioner
724    
725            :rtype: ``int``
726            """
727            return self.__pre_sweeps
728        def setNumPostSweeps(self,sweeps=2):
729            """
730            Sets the number of sweeps in the post-smoothing step of a multi level solver or preconditioner
731    
732            :param sweeps: number of sweeps
733            :type sweeps: positive ``int``
734            """
735            sweeps=int(sweeps)
736            if sweeps<1:
737               raise ValueError,"number of sweeps must be positive."
738            self.__post_sweeps=sweeps
739        def getNumPostSweeps(self):
740            """
741            Returns he number of sweeps in the post-smoothing step of a multi level solver or preconditioner
742    
743            :rtype: ``int``
744            """
745            return self.__post_sweeps
746    
747        def setTolerance(self,rtol=1.e-8):
748            """
749            Sets the relative tolerance for the solver
750    
751            :param rtol: relative tolerance
752            :type rtol: non-negative ``float``
753            """
754            rtol=float(rtol)
755            if rtol<0 or rtol>1:
756               raise ValueError,"tolerance must be non-negative and less or equal 1."
757            self.__tolerance=rtol
758        def getTolerance(self):
759            """
760            Returns the relative tolerance for the solver
761    
762            :rtype: ``float``
763            """
764            return self.__tolerance
765        def setAbsoluteTolerance(self,atol=0.):
766            """
767            Sets the absolute tolerance for the solver
768    
769            :param atol:  absolute tolerance
770            :type atol: non-negative ``float``
771            """
772            atol=float(atol)
773            if atol<0:
774               raise ValueError,"tolerance must be non-negative."
775            self.__absolute_tolerance=atol
776        def getAbsoluteTolerance(self):
777            """
778            Returns the absolute tolerance for the solver
779    
780            :rtype: ``float``
781            """
782            return self.__absolute_tolerance
783    
784        def setInnerTolerance(self,rtol=0.9):
785            """
786             Sets the relative tolerance for an inner iteration scheme for instance on the coarsest level in a multi-level scheme.
787    
788            :param rtol: inner relative tolerance
789            :type rtol: positive ``float``
790            """
791            rtol=float(rtol)
792            if rtol<=0 or rtol>1:
793                raise ValueError,"tolerance must be positive and less or equal 1."
794            self.__inner_tolerance=rtol
795        def getInnerTolerance(self):
796            """
797            Returns the relative tolerance for an inner iteration scheme
798    
799            :rtype: ``float``
800            """
801            return self.__inner_tolerance
802        def setDropTolerance(self,drop_tol=0.01):
803            """
804            Sets the relative drop tolerance in ILUT
805    
806            :param drop_tol: drop tolerance
807            :type drop_tol: positive ``float``
808            """
809            drop_tol=float(drop_tol)
810            if drop_tol<=0 or drop_tol>1:
811                raise ValueError,"drop tolerance must be positive and less or equal 1."
812            self.__drop_tolerance=drop_tol
813        def getDropTolerance(self):
814            """
815            Returns the relative drop tolerance in ILUT
816    
817            :rtype: ``float``
818            """
819            return self.__drop_tolerance
820        def setDropStorage(self,storage=2.):
821            """
822            Sets the maximum allowed increase in storage for ILUT. ``storage`` =2 would mean that
823            a doubling of the storage needed for the coefficient matrix is allowed in the ILUT factorization.
824    
825            :param storage: allowed storage increase
826            :type storage: ``float``
827            """
828            storage=float(storage)
829            if storage<1:
830                raise ValueError,"allowed storage increase must be greater or equal to 1."
831            self.__drop_storage=storage
832        def getDropStorage(self):
833        
834            """
835            Returns the maximum allowed increase in storage for ILUT
836    
837            :rtype: ``float``
838            """
839            return self.__drop_storage
840        def setRelaxationFactor(self,factor=0.3):
841            """
842            Sets the relaxation factor used to add dropped elements in RILU to the main diagonal.
843    
844            :param factor: relaxation factor
845            :type factor: ``float``
846            :note: RILU with a relaxation factor 0 is identical to ILU0
847            """
848            factor=float(factor)
849            if factor<0:
850                raise ValueError,"relaxation factor must be non-negative."
851            self.__relaxation=factor
852        def getRelaxationFactor(self):
853        
854            """
855            Returns the relaxation factor used to add dropped elements in RILU to the main diagonal.
856    
857            :rtype: ``float``
858            """
859            return self.__relaxation
860        def isSymmetric(self):
861            """
862            Checks if symmetry of the  coefficient matrix is indicated.
863    
864            :return: True if a symmetric PDE is indicated, False otherwise
865            :rtype: ``bool``
866            """
867            return self.__symmetric
868        def setSymmetryOn(self):
869            """
870            Sets the symmetry flag to indicate that the coefficient matrix is symmetric.
871            """
872            self.__symmetric=True
873        def setSymmetryOff(self):
874            """
875            Clears the symmetry flag for the coefficient matrix.
876            """
877            self.__symmetric=False
878        def setSymmetry(self,flag=False):
879            """
880            Sets the symmetry flag for the coefficient matrix to ``flag``.
881    
882            :param flag: If True, the symmetry flag is set otherwise reset.
883            :type flag: ``bool``
884            """
885            if flag:
886                self.setSymmetryOn()
887            else:
888                self.setSymmetryOff()
889        def isVerbose(self):
890            """
891            Returns ``True`` if the solver is expected to be verbose.
892    
893            :return: True if verbosity of switched on.
894            :rtype: ``bool``
895            """
896            return self.__verbose
897    
898        def setVerbosityOn(self):
899            """
900            Switches the verbosity of the solver on.
901            """
902            self.__verbose=True
903        def setVerbosityOff(self):
904            """
905            Switches the verbosity of the solver off.
906            """
907            self.__verbose=False
908        def setVerbosity(self,verbose=False):
909            """
910            Sets the verbosity flag for the solver to ``flag``.
911    
912            :param verbose: If ``True``, the verbosity of the solver is switched on.
913            :type verbose: ``bool``
914            """
915            if verbose:
916                self.setVerbosityOn()
917            else:
918                self.setVerbosityOff()
919            
920        def adaptInnerTolerance(self):
921            """
922            Returns ``True`` if the tolerance of the inner solver is selected automatically.
923            Otherwise the inner tolerance set by `setInnerTolerance` is used.
924    
925            :return: ``True`` if inner tolerance adaption is chosen.
926            :rtype: ``bool``
927            """
928            return self.__adapt_inner_tolerance
929    
930        def setInnerToleranceAdaptionOn(self):
931            """
932            Switches the automatic selection of inner tolerance on
933            """
934            self.__adapt_inner_tolerance=True
935        def setInnerToleranceAdaptionOff(self):
936            """
937            Switches the automatic selection of inner tolerance off.
938            """
939            self.__adapt_inner_tolerance=False
940        def setInnerToleranceAdaption(self,adapt=True):
941            """
942            Sets the flag to indicate automatic selection of the inner tolerance.
943    
944            :param adapt: If ``True``, the inner tolerance is selected automatically.
945            :type adapt: ``bool``
946            """
947            if adapt:
948                self.setInnerToleranceAdaptionOn()
949            else:
950                self.setInnerToleranceAdaptionOff()
951    
952        def acceptConvergenceFailure(self):
953            """
954            Returns ``True`` if a failure to meet the stopping criteria within the
955            given number of iteration steps is not raising in exception. This is useful
956            if a solver is used in a non-linear context where the non-linear solver can
957            continue even if the returned the solution does not necessarily meet the
958            stopping criteria. One can use the `hasConverged` method to check if the
959            last call to the solver was successful.
960    
961            :return: ``True`` if a failure to achieve convergence is accepted.
962            :rtype: ``bool``
963            """
964            return self.__accept_convergence_failure
965    
966        def setAcceptanceConvergenceFailureOn(self):
967            """
968            Switches the acceptance of a failure of convergence on
969            """
970            self.__accept_convergence_failure=True
971        def setAcceptanceConvergenceFailureOff(self):
972            """
973            Switches the acceptance of a failure of convergence off.
974            """
975            self.__accept_convergence_failure=False
976        def setAcceptanceConvergenceFailure(self,accept=False):
977            """
978            Sets the flag to indicate the acceptance of a failure of convergence.
979    
980            :param accept: If ``True``, any failure to achieve convergence is accepted.
981            :type accept: ``bool``
982            """
983            if accept:
984                self.setAcceptanceConvergenceFailureOn()
985            else:
986                self.setAcceptanceConvergenceFailureOff()
987    
988        def useLocalPreconditioner(self):
989            """
990            Returns ``True`` if the preconditoner is applied locally on each MPI. This reducess communication costs
991            and speeds up the application of the preconditioner but at the costs of more iteration steps. This can be an
992            advantage on clusters with slower interconnects.
993    
994            :return: ``True`` if local preconditioning is applied
995            :rtype: ``bool``
996            """
997            return self.__use_local_preconditioner
998    
999        def setLocalPreconditionerOn(self):
1000            """
1001            Sets the flag to use  local preconditioning to on
1002            """
1003            self.__use_local_preconditioner=True
1004        def setLocalPreconditionerOff(self):
1005            """
1006            Sets the flag to use  local preconditioning to off
1007            """
1008            self.__use_local_preconditioner=False
1009            
1010        def setLocalPreconditioner(self,use=False):
1011            """
1012            Sets the flag to use  local preconditioning
1013    
1014            :param use: If ``True``, local proconditioning on each MPI rank is applied
1015            :type use: ``bool``
1016            """
1017            if use:
1018                self.setUseLocalPreconditionerOn()
1019            else:
1020                self.setUseLocalPreconditionerOff()
1021                
1022        def setMinCoarseMatrixSparsity(self,sparsity=0.05):
1023           """
1024           Sets the minimum sparsity on the coarsest level. Typically
1025           a direct solver is used when the sparsity becomes bigger than
1026           the set limit.
1027          
1028           :param sparsity: minimual sparsity
1029           :type sparsity: ``float``
1030           """
1031           sparsity=float(sparsity)
1032           if factor<0:
1033         raise ValueError,"sparsity must be non-negative."
1034           if factor>1:
1035               raise ValueError,"sparsity must be less than 1."
1036           self.__min_sparsity=factor
1037        def setMinCoarseMatrixSparsity(self,sparsity=0.05):
1038           """
1039           Sets the minimum sparsity on the coarsest level. Typically
1040           a direct solver is used when the sparsity becomes bigger than
1041           the set limit.
1042          
1043           :param sparsity: minimual sparsity
1044           :type sparsity: ``float``
1045           """
1046           sparsity=float(sparsity)
1047           if sparsity<0:
1048         raise ValueError,"sparsity must be non-negative."
1049           if sparsity>1:
1050              raise ValueError,"sparsity must be less than 1."
1051           self.__min_sparsity=sparsity
1052    
1053        def getMinCoarseMatrixSparsity(self):
1054          """
1055          Returns the minimum sparsity on the coarsest level. Typically
1056          a direct solver is used when the sparsity becomes bigger than
1057          the set limit.
1058      
1059          :return: minimual sparsity
1060          :rtype: ``float``
1061          """
1062          return self.__min_sparsity
1063    
1064        def setNumRefinements(self,refinements=2):
1065          """
1066          Sets the number of refinement steps to refine the solution when a direct solver is applied.
1067      
1068          :param refinements: number of refinements
1069          :type refinements: non-negative ``int``
1070          """
1071          refinements=int(refinements)
1072          if refinements<0:
1073         raise ValueError,"number of refinements must be non-negative."
1074          self.__refinements=refinements
1075      
1076        def getNumRefinements(self):
1077           """
1078           Returns the number of refinement steps to refine the solution when a direct solver is applied.
1079          
1080           :rtype: non-negative ``int``
1081           """
1082           return self.__refinements
1083    
1084        def setNumCoarseMatrixRefinements(self,refinements=2):
1085          """
1086          Sets the number of refinement steps to refine the solution on the coarset level when a direct solver is applied.
1087      
1088          :param refinements: number of refinements
1089          :type refinements: non-negative ``int``
1090          """
1091          refinements=int(refinements)
1092          if refinements<0:
1093         raise ValueError,"number of refinements must be non-negative."
1094          self.__coarse_refinements=refinements
1095      
1096        def getNumCoarseMatrixRefinements(self):
1097          """
1098          Returns the number of resfinement steps to refine the solution on the coarset level when a direct solver is applied.
1099          
1100          :rtype: non-negative ``int``
1101          """
1102          return self.__coarse_refinements
1103    
1104        def usePanel(self):
1105            """
1106            Returns ``True`` if a panel is used to serach for unknown in the AMG coarsening, The panel approach is normally faster
1107            but can lead to larger coarse level systems.
1108            
1109            :return: ``True`` if a panel is used to find unknowns in AMG coarsening
1110            :rtype: ``bool``
1111            """
1112            return self.__use_panel
1113    
1114        def setUsePanelOn(self):
1115            """
1116            Sets the flag to use a panel to find unknowns in AMG coarsening
1117            """
1118            self.__use_panel=True
1119            
1120        def setUsePanelOff(self):
1121            """
1122            Sets the flag to use a panel to find unknowns in AMG coarsening to off
1123            """
1124            self.__use_panel=False
1125            
1126        def setUsePanel(self,use=True):
1127            """
1128            Sets the flag to use  a panel to find unknowns in AMG coarsening
1129    
1130            :param use: If ``True``,a panel is used to find unknowns in AMG coarsening
1131            :type use: ``bool``
1132            """
1133            if use:
1134                self.setUsePanelOn()
1135            else:
1136                self.setUsePanelOff()
1137                
1138        def useDirectInterpolation(self):
1139            """
1140            Returns ``True`` if direct interpolation is used in AMG.
1141    
1142        :return: ``True`` if direct interpolation is used in AMG.
1143            :rtype: ``bool``
1144            """
1145            return self.__use_direct_interpolation
1146    
1147        def setUseDirectInterpolationOn(self):
1148            """
1149            Sets the flag to use direct interpolation in AMG
1150            """
1151            self.__use_direct_interpolation=True
1152            
1153        def setUseDirectInterpolationOff(self):
1154            """
1155            Sets the flag to use direct interpolation in AMG
1156            """
1157            self.__use_direct_interpolation=False
1158            
1159        def setUseDirectInterpolation(self,use=False):
1160            """
1161            Sets the flag to use direct interpolation in AMG
1162    
1163        :param use: If ``True``, direct interpolation in AMG
1164        :type use: ``bool``
1165            """
1166            if use:
1167                self.setUseDirectInterpolationOn()
1168            else:
1169                self.setUseDirectInterpolationOff()
1170    
1171  def _CompTuple2(t1,t2):  
1172    class IllegalCoefficient(ValueError):
1173       """
1174       Exception that is raised if an illegal coefficient of the general or
1175       particular PDE is requested.
1176     """     """
1177     Compare two tuples     pass
1178    
1179     \param t1 The first tuple  class IllegalCoefficientValue(ValueError):
    \param t2 The second tuple  
1180     """     """
1181       Exception that is raised if an incorrect value for a coefficient is used.
1182       """
1183       pass
1184    
1185     dif=t1[0]+t1[1]-(t2[0]+t2[1])  class IllegalCoefficientFunctionSpace(ValueError):
1186     if dif<0: return 1     """
1187     elif dif>0: return -1     Exception that is raised if an incorrect function space for a coefficient
1188     else: return 0     is used.
1189       """
1190    
1191  def ELMAN_RAMAGE(P):  class UndefinedPDEError(ValueError):
1192      return (P-1.).wherePositive()*0.5*(1.-1./(P+1.e-15))     """
1193       Exception that is raised if a PDE is not fully defined yet.
1194       """
1195       pass
1196    
1197  def SIMPLIFIED_BROOK_HUGHES(P):  class PDECoef(object):
1198      c=(P-3.).whereNegative()      """
1199      return P/6.*c+1./2.*(1.-c)      A class for describing a PDE coefficient.
1200    
1201  def HALF(P):      :cvar INTERIOR: indicator that coefficient is defined on the interior of
1202      return escript.Scalar(0.5,P.getFunctionSpace())                      the PDE domain
1203        :cvar BOUNDARY: indicator that coefficient is defined on the boundary of
1204                        the PDE domain
1205        :cvar CONTACT: indicator that coefficient is defined on the contact region
1206                       within the PDE domain
1207        :cvar INTERIOR_REDUCED: indicator that coefficient is defined on the
1208                                interior of the PDE domain using a reduced
1209                                integration order
1210        :cvar BOUNDARY_REDUCED: indicator that coefficient is defined on the
1211                                boundary of the PDE domain using a reduced
1212                                integration order
1213        :cvar CONTACT_REDUCED: indicator that coefficient is defined on the contact
1214                               region within the PDE domain using a reduced
1215                               integration order
1216        :cvar SOLUTION: indicator that coefficient is defined through a solution of
1217                        the PDE
1218        :cvar REDUCED: indicator that coefficient is defined through a reduced
1219                       solution of the PDE
1220        :cvar BY_EQUATION: indicator that the dimension of the coefficient shape is
1221                           defined by the number of PDE equations
1222        :cvar BY_SOLUTION: indicator that the dimension of the coefficient shape is
1223                           defined by the number of PDE solutions
1224        :cvar BY_DIM: indicator that the dimension of the coefficient shape is
1225                      defined by the spatial dimension
1226        :cvar OPERATOR: indicator that the the coefficient alters the operator of
1227                        the PDE
1228        :cvar RIGHTHANDSIDE: indicator that the the coefficient alters the right
1229                             hand side of the PDE
1230        :cvar BOTH: indicator that the the coefficient alters the operator as well
1231                    as the right hand side of the PDE
1232    
 class PDECoefficient:  
     """  
     A class for PDE coefficients  
1233      """      """
     # identifier for location of Data objects defining COEFFICIENTS  
1234      INTERIOR=0      INTERIOR=0
1235      BOUNDARY=1      BOUNDARY=1
1236      CONTACT=2      CONTACT=2
1237      CONTINUOUS=3      SOLUTION=3
1238      # identifier in the pattern of COEFFICIENTS:      REDUCED=4
1239      # the pattern is a tuple of EQUATION,SOLUTION,DIM where DIM represents the spatial dimension, EQUATION the number of equations and SOLUTION the      BY_EQUATION=5
1240      # number of unknowns.      BY_SOLUTION=6
1241      EQUATION=3      BY_DIM=7
1242      SOLUTION=4      OPERATOR=10
1243      DIM=5      RIGHTHANDSIDE=11
1244      # indicator for what is altered if the coefficient is altered:      BOTH=12
1245      OPERATOR=5      INTERIOR_REDUCED=13
1246      RIGHTHANDSIDE=6      BOUNDARY_REDUCED=14
1247      BOTH=7      CONTACT_REDUCED=15
1248      def __init__(self,where,pattern,altering):  
1249         """      def __init__(self, where, pattern, altering):
1250         Initialise a PDE Coefficient type         """
1251           Initialises a PDE coefficient type.
1252    
1253           :param where: describes where the coefficient lives
1254           :type where: one of `INTERIOR`, `BOUNDARY`, `CONTACT`, `SOLUTION`,
1255                        `REDUCED`, `INTERIOR_REDUCED`, `BOUNDARY_REDUCED`,
1256                        `CONTACT_REDUCED`
1257           :param pattern: describes the shape of the coefficient and how the shape
1258                           is built for a given spatial dimension and numbers of
1259                           equations and solutions in then PDE. For instance,
1260                           (`BY_EQUATION`,`BY_SOLUTION`,`BY_DIM`) describes a
1261                           rank 3 coefficient which is instantiated as shape (3,2,2)
1262                           in case of three equations and two solution components
1263                           on a 2-dimensional domain. In the case of single equation
1264                           and a single solution component the shape components
1265                           marked by `BY_EQUATION` or `BY_SOLUTION` are dropped.
1266                           In this case the example would be read as (2,).
1267           :type pattern: ``tuple`` of `BY_EQUATION`, `BY_SOLUTION`, `BY_DIM`
1268           :param altering: indicates what part of the PDE is altered if the
1269                            coefficient is altered
1270           :type altering: one of `OPERATOR`, `RIGHTHANDSIDE`, `BOTH`
1271         """         """
1272           super(PDECoef, self).__init__()
1273         self.what=where         self.what=where
1274         self.pattern=pattern         self.pattern=pattern
1275         self.altering=altering         self.altering=altering
# Line 64  class PDECoefficient: Line 1277  class PDECoefficient:
1277    
1278      def resetValue(self):      def resetValue(self):
1279         """         """
1280         resets coefficient value to default         Resets the coefficient value to the default.
1281         """         """
1282         self.value=escript.Data()         self.value=escript.Data()
1283    
1284      def getFunctionSpace(self,domain):      def getFunctionSpace(self,domain,reducedEquationOrder=False,reducedSolutionOrder=False):
1285         """         """
1286         defines the FunctionSpace of the coefficient on the domain         Returns the `FunctionSpace` of the coefficient.
1287    
1288         @param domain:         :param domain: domain on which the PDE uses the coefficient
1289         """         :type domain: `Domain`
1290         if self.what==self.INTERIOR: return escript.Function(domain)         :param reducedEquationOrder: True to indicate that reduced order is used
1291         elif self.what==self.BOUNDARY: return escript.FunctionOnBoundary(domain)                                      to represent the equation
1292         elif self.what==self.CONTACT: return escript.FunctionOnContactZero(domain)         :type reducedEquationOrder: ``bool``
1293         elif self.what==self.CONTINUOUS: return escript.ContinuousFunction(domain)         :param reducedSolutionOrder: True to indicate that reduced order is used
1294                                        to represent the solution
1295           :type reducedSolutionOrder: ``bool``
1296           :return: `FunctionSpace` of the coefficient
1297           :rtype: `FunctionSpace`
1298           """
1299           if self.what==self.INTERIOR:
1300                return escript.Function(domain)
1301           elif self.what==self.INTERIOR_REDUCED:
1302                return escript.ReducedFunction(domain)
1303           elif self.what==self.BOUNDARY:
1304                return escript.FunctionOnBoundary(domain)
1305           elif self.what==self.BOUNDARY_REDUCED:
1306                return escript.ReducedFunctionOnBoundary(domain)
1307           elif self.what==self.CONTACT:
1308                return escript.FunctionOnContactZero(domain)
1309           elif self.what==self.CONTACT_REDUCED:
1310                return escript.ReducedFunctionOnContactZero(domain)
1311           elif self.what==self.SOLUTION:
1312                if reducedEquationOrder and reducedSolutionOrder:
1313                    return escript.ReducedSolution(domain)
1314                else:
1315                    return escript.Solution(domain)
1316           elif self.what==self.REDUCED:
1317                return escript.ReducedSolution(domain)
1318    
1319      def getValue(self):      def getValue(self):
1320         """         """
1321         returns the value of the coefficient:         Returns the value of the coefficient.
1322    
1323           :return: value of the coefficient
1324           :rtype: `Data`
1325         """         """
1326         return self.value         return self.value
1327        
1328      def setValue(self,newValue):      def setValue(self,domain,numEquations=1,numSolutions=1,reducedEquationOrder=False,reducedSolutionOrder=False,newValue=None):
        """  
        set the value of the coefficient to new value  
1329         """         """
1330           Sets the value of the coefficient to a new value.
1331    
1332           :param domain: domain on which the PDE uses the coefficient
1333           :type domain: `Domain`
1334           :param numEquations: number of equations of the PDE
1335           :type numEquations: ``int``
1336           :param numSolutions: number of components of the PDE solution
1337           :type numSolutions: ``int``
1338           :param reducedEquationOrder: True to indicate that reduced order is used
1339                                        to represent the equation
1340           :type reducedEquationOrder: ``bool``
1341           :param reducedSolutionOrder: True to indicate that reduced order is used
1342                                        to represent the solution
1343           :type reducedSolutionOrder: ``bool``
1344           :param newValue: number of components of the PDE solution
1345           :type newValue: any object that can be converted into a
1346                           `Data` object with the appropriate shape
1347                           and `FunctionSpace`
1348           :raise IllegalCoefficientValue: if the shape of the assigned value does
1349                                           not match the shape of the coefficient
1350           :raise IllegalCoefficientFunctionSpace: if unable to interpolate value
1351                                                   to appropriate function space
1352           """
1353           if newValue==None:
1354               newValue=escript.Data()
1355           elif isinstance(newValue,escript.Data):
1356               if not newValue.isEmpty():
1357                  if not newValue.getFunctionSpace() == self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder):
1358                    try:
1359                      newValue=escript.Data(newValue,self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder))
1360                    except:
1361                      raise IllegalCoefficientFunctionSpace,"Unable to interpolate coefficient to function space %s"%self.getFunctionSpace(domain)
1362           else:
1363               newValue=escript.Data(newValue,self.getFunctionSpace(domain,reducedEquationOrder,reducedSolutionOrder))
1364           if not newValue.isEmpty():
1365               if not self.getShape(domain,numEquations,numSolutions)==newValue.getShape():
1366                   raise IllegalCoefficientValue,"Expected shape of coefficient is %s but actual shape is %s."%(self.getShape(domain,numEquations,numSolutions),newValue.getShape())
1367         self.value=newValue         self.value=newValue
1368        
1369      def isAlteringOperator(self):      def isAlteringOperator(self):
1370          """          """
1371      return true if the operator of the PDE is changed when the coefficient is changed          Checks if the coefficient alters the operator of the PDE.
1372      """  
1373            :return: True if the operator of the PDE is changed when the
1374                     coefficient is changed
1375            :rtype: ``bool``
1376            """
1377          if self.altering==self.OPERATOR or self.altering==self.BOTH:          if self.altering==self.OPERATOR or self.altering==self.BOTH:
1378              return not None              return not None
1379          else:          else:
# Line 102  class PDECoefficient: Line 1381  class PDECoefficient:
1381    
1382      def isAlteringRightHandSide(self):      def isAlteringRightHandSide(self):
1383          """          """
1384      return true if the right hand side of the PDE is changed when the coefficient is changed          Checks if the coefficient alters the right hand side of the PDE.
1385      """  
1386            :rtype: ``bool``
1387            :return: True if the right hand side of the PDE is changed when the
1388                     coefficient is changed, ``None`` otherwise.
1389            """
1390          if self.altering==self.RIGHTHANDSIDE or self.altering==self.BOTH:          if self.altering==self.RIGHTHANDSIDE or self.altering==self.BOTH:
1391              return not None              return not None
1392          else:          else:
1393              return None              return None
1394    
1395      def estimateNumEquationsAndNumSolutions(self,shape=(),dim=3):      def estimateNumEquationsAndNumSolutions(self,domain,shape=()):
1396         """         """
1397         tries to estimate the number of equations in a given tensor shape for a given spatial dimension dim         Tries to estimate the number of equations and number of solutions if
1398           the coefficient has the given shape.
1399    
1400         @param shape:         :param domain: domain on which the PDE uses the coefficient
1401         @param dim:         :type domain: `Domain`
1402           :param shape: suggested shape of the coefficient
1403           :type shape: ``tuple`` of ``int`` values
1404           :return: the number of equations and number of solutions of the PDE if
1405                    the coefficient has given shape. If no appropriate numbers
1406                    could be identified, ``None`` is returned
1407           :rtype: ``tuple`` of two ``int`` values or ``None``
1408         """         """
1409           dim=domain.getDim()
1410         if len(shape)>0:         if len(shape)>0:
1411             num=max(shape)+1             num=max(shape)+1
1412         else:         else:
1413             num=1             num=1
1414         search=[]         search=[]
1415         for u in range(num):         if self.definesNumEquation() and self.definesNumSolutions():
1416            for e in range(num):            for u in range(num):
1417               search.append((e,u))               for e in range(num):
1418         search.sort(_CompTuple2)                  search.append((e,u))
1419         for item in search:            search.sort(self.__CompTuple2)
1420               s=self.buildShape(item[0],item[1],dim)            for item in search:
1421                 s=self.getShape(domain,item[0],item[1])
1422               if len(s)==0 and len(shape)==0:               if len(s)==0 and len(shape)==0:
1423                   return (1,1)                   return (1,1)
1424               else:               else:
1425                   if s==shape: return item                   if s==shape: return item
1426           elif self.definesNumEquation():
1427              for e in range(num,0,-1):
1428                 s=self.getShape(domain,e,0)
1429                 if len(s)==0 and len(shape)==0:
1430                     return (1,None)
1431                 else:
1432                     if s==shape: return (e,None)
1433    
1434           elif self.definesNumSolutions():
1435              for u in range(num,0,-1):
1436                 s=self.getShape(domain,0,u)
1437                 if len(s)==0 and len(shape)==0:
1438                     return (None,1)
1439                 else:
1440                     if s==shape: return (None,u)
1441         return None         return None
1442    
1443      def buildShape(self,e=1,u=1,dim=3):      def definesNumSolutions(self):
1444          """         """
1445      builds the required shape for a given number of equations e, number of unknowns u and spatial dimension dim         Checks if the coefficient allows to estimate the number of solution
1446           components.
1447    
1448           :return: True if the coefficient allows an estimate of the number of
1449                    solution components, False otherwise
1450           :rtype: ``bool``
1451           """
1452           for i in self.pattern:
1453                 if i==self.BY_SOLUTION: return True
1454           return False
1455    
1456        def definesNumEquation(self):
1457           """
1458           Checks if the coefficient allows to estimate the number of equations.
1459    
1460      @param e:         :return: True if the coefficient allows an estimate of the number of
1461      @param u:                  equations, False otherwise
1462      @param dim:         :rtype: ``bool``
1463      """         """
1464          s=()         for i in self.pattern:
1465          for i in self.pattern:               if i==self.BY_EQUATION: return True
1466               if i==self.EQUATION:         return False
1467                  if e>1: s=s+(e,)  
1468               elif i==self.SOLUTION:      def __CompTuple2(self,t1,t2):
1469                  if u>1: s=s+(u,)        """
1470          Compares two tuples of possible number of equations and number of
1471          solutions.
1472    
1473          :param t1: the first tuple
1474          :param t2: the second tuple
1475          :return: 0, 1, or -1
1476          """
1477    
1478          dif=t1[0]+t1[1]-(t2[0]+t2[1])
1479          if dif<0: return 1
1480          elif dif>0: return -1
1481          else: return 0
1482    
1483        def getShape(self,domain,numEquations=1,numSolutions=1):
1484           """
1485           Builds the required shape of the coefficient.
1486    
1487           :param domain: domain on which the PDE uses the coefficient
1488           :type domain: `Domain`
1489           :param numEquations: number of equations of the PDE
1490           :type numEquations: ``int``
1491           :param numSolutions: number of components of the PDE solution
1492           :type numSolutions: ``int``
1493           :return: shape of the coefficient
1494           :rtype: ``tuple`` of ``int`` values
1495           """
1496           dim=domain.getDim()
1497           s=()
1498           for i in self.pattern:
1499                 if i==self.BY_EQUATION:
1500                    if numEquations>1: s=s+(numEquations,)
1501                 elif i==self.BY_SOLUTION:
1502                    if numSolutions>1: s=s+(numSolutions,)
1503               else:               else:
1504                  s=s+(dim,)                  s=s+(dim,)
1505          return s         return s
1506    
1507  class LinearPDE:  #====================================================================================================================
    """  
    Class to handle a linear PDE  
     
    class to define a linear PDE of the form  
1508    
1509     \f[  class LinearProblem(object):
1510       -(A_{ijkl}u_{k,l})_{,j} -(B_{ijk}u_k)_{,j} + C_{ikl}u_{k,l} +D_{ik}u_k = - (X_{ij})_{,j} + Y_i     """
1511     \f]     This is the base class to define a general linear PDE-type problem for
1512       for an unknown function *u* on a given domain defined through a
1513       `Domain` object. The problem can be given as a single
1514       equation or as a system of equations.
1515    
1516     with boundary conditons:     The class assumes that some sort of assembling process is required to form
1517       a problem of the form
1518    
1519     \f[     *L u=f*
    n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_{ik}u_k = - n_j*X_{ij} + y_i  
    \f]  
1520    
1521     and contact conditions     where *L* is an operator and *f* is the right hand side. This operator
1522       problem will be solved to get the unknown *u*.
1523    
1524     \f[     """
1525     n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_contact_{ik}[u_k] = - n_j*X_{ij} + y_contact_i     def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):
1526     \f]       """
1527         Initializes a linear problem.
1528    
1529     and constraints:       :param domain: domain of the PDE
1530         :type domain: `Domain`
1531         :param numEquations: number of equations. If ``None`` the number of
1532                              equations is extracted from the coefficients.
1533         :param numSolutions: number of solution components. If ``None`` the number
1534                              of solution components is extracted from the
1535                              coefficients.
1536         :param debug: if True debug information is printed
1537    
1538     \f[       """
1539     u_i=r_i \quad \mathrm{where} \quad q_i>0       super(LinearProblem, self).__init__()
    \f]  
   
    """  
    TOL=1.e-13  
    DEFAULT_METHOD=util.DEFAULT_METHOD  
    DIRECT=util.DIRECT  
    CHOLEVSKY=util.CHOLEVSKY  
    PCG=util.PCG  
    CR=util.CR  
    CGS=util.CGS  
    BICGSTAB=util.BICGSTAB  
    SSOR=util.SSOR  
    GMRES=util.GMRES  
    PRES20=util.PRES20  
    ILU0=util.ILU0  
    JACOBI=util.JACOBI  
   
    def __init__(self,domain,numEquations=0,numSolutions=0):  
      """  
      initializes a new linear PDE.  
   
      @param args:  
      """  
      # COEFFICIENTS can be overwritten by subclasses:  
      self.__COEFFICIENTS={  
        "A"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM,PDECoefficient.SOLUTION,PDECoefficient.DIM),PDECoefficient.OPERATOR),  
        "B"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),  
        "C"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION,PDECoefficient.DIM),PDECoefficient.OPERATOR),  
        "D"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),  
        "X"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM),PDECoefficient.RIGHTHANDSIDE),  
        "Y"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),  
        "d"         : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),  
        "y"         : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),  
        "d_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),  
        "y_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),  
        "r"         : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),  
        "q"         : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.SOLUTION,),PDECoefficient.BOTH)}  
1540    
1541       self.COEFFICIENTS=self.__COEFFICIENTS       self.__debug=debug
      # initialize attributes  
      self.__debug=None  
1542       self.__domain=domain       self.__domain=domain
1543       self.__numEquations=numEquations       self.__numEquations=numEquations
1544       self.__numSolutions=numSolutions       self.__numSolutions=numSolutions
1545       self.cleanCoefficients()       self.__altered_coefficients=False
1546         self.__reduce_equation_order=False
1547       self.__operator=escript.Operator()       self.__reduce_solution_order=False
      self.__operator_isValid=False  
      self.__righthandside=escript.Data()  
      self.__righthandside_isValid=False  
      self.__solution=escript.Data()  
      self.__solution_isValid=False  
   
      # set some default values:  
      self.__homogeneous_constraint=True  
      self.__row_function_space=escript.Solution(self.__domain)  
      self.__column_function_space=escript.Solution(self.__domain)  
      self.__tolerance=1.e-8  
      self.__solver_method=util.DEFAULT_METHOD  
      self.__matrix_type=self.__domain.getSystemMatrixTypeId(util.DEFAULT_METHOD,False)  
1548       self.__sym=False       self.__sym=False
1549       self.__lumping=False       self.setSolverOptions()
1550         self.__is_RHS_valid=False
1551     def createCoefficient(self, name):       self.__is_operator_valid=False
1552         self.__COEFFICIENTS={}
1553         self.__solution_rtol=1.e99
1554         self.__solution_atol=1.e99
1555         self.setSymmetryOff()
1556         # initialize things:
1557         self.resetAllCoefficients()
1558         self.initializeSystem()
1559       # ==========================================================================
1560       #    general stuff:
1561       # ==========================================================================
1562       def __str__(self):
1563         """
1564         Returns a string representation of the PDE.
1565    
1566         :return: a simple representation of the PDE
1567         :rtype: ``str``
1568         """
1569         return "<LinearProblem %d>"%id(self)
1570       # ==========================================================================
1571       #    debug :
1572       # ==========================================================================
1573       def setDebugOn(self):
1574       """       """
1575       create a data object corresponding to coefficient name       Switches debug output on.
      @param name:  
1576       """       """
1577       return escript.Data(shape = getShapeOfCoefficient(name), \       self.__debug=not None
                          what = getFunctionSpaceForCoefficient(name))  
1578    
1579     def __del__(self):     def setDebugOff(self):
1580       pass       """
1581         Switches debug output off.
1582         """
1583         self.__debug=None
1584    
1585     def getCoefficient(self,name):     def setDebug(self, flag):
1586       """       """
1587       return the value of the parameter name       Switches debug output on if ``flag`` is True otherwise it is switched off.
1588    
1589       @param name:       :param flag: desired debug status
1590         :type flag: ``bool``
1591       """       """
1592       return self.COEFFICIENTS[name].getValue()       if flag:
1593             self.setDebugOn()
1594         else:
1595             self.setDebugOff()
1596    
1597     def getCoefficientOfPDE(self,name):     def trace(self,text):
1598       """       """
1599       return the value of the coefficient name of the general PDE.       Prints the text message if debug mode is switched on.
      This method is called by the assembling routine it can be  
      overwritten to map coefficients of a particualr PDE to the general PDE.  
1600    
1601       @param name:       :param text: message to be printed
1602         :type text: ``string``
1603       """       """
1604       return self.getCoefficient(name)       if self.__debug: print "%s: %s"%(str(self),text)
1605    
1606     def hasCoefficient(self,name):     # ==========================================================================
1607        """     # some service functions:
1608        return true if name is the name of a coefficient     # ==========================================================================
1609       def introduceCoefficients(self,**coeff):
1610           """
1611           Introduces new coefficients into the problem.
1612    
1613        @param name:         Use:
       """  
       return self.COEFFICIENTS.has_key(name)  
1614    
1615     def hasPDECoefficient(self,name):         p.introduceCoefficients(A=PDECoef(...), B=PDECoef(...))
       """  
       return true if name is the name of a coefficient  
1616    
1617        @param name:         to introduce the coefficients *A* and *B*.
1618        """         """
1619        return self.__COEFFICIENTS.has_key(name)         for name, type in coeff.items():
1620               if not isinstance(type,PDECoef):
1621                  raise ValueError,"coefficient %s has no type."%name
1622               self.__COEFFICIENTS[name]=type
1623               self.__COEFFICIENTS[name].resetValue()
1624               self.trace("coefficient %s has been introduced."%name)
1625    
1626     def getFunctionSpaceForEquation(self):     def getDomain(self):
      """  
      return true if the test functions should use reduced order  
1627       """       """
1628       return self.__row_function_space       Returns the domain of the PDE.
1629    
1630     def getFunctionSpaceForSolution(self):       :return: the domain of the PDE
1631       """       :rtype: `Domain`
      return true if the interpolation of the solution should use reduced order  
1632       """       """
1633       return self.__column_function_space       return self.__domain
1634       def getDomainStatus(self):
    def setValue(self,**coefficients):  
       """  
       sets new values to coefficients  
   
       @param coefficients:  
       """  
       self.__setValue(**coefficients)  
         
   
    def cleanCoefficients(self):  
1635       """       """
1636       resets all coefficients to default values.       Return the status indicator of the domain
1637       """       """
1638       for i in self.COEFFICIENTS.iterkeys():       return self.getDomain().getStatus()
          self.COEFFICIENTS[i].resetValue()  
1639    
1640     def createNewCoefficient(self,name):     def getSystemStatus(self):
1641       """       """
1642       returns a new coefficient appropriate for coefficient name:       Return the domain status used to build the current system
1643       """       """
1644       return escript.Data(0,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))       return self.__system_status
1645             def setSystemStatus(self,status=None):
    def createNewCoefficientOfPDE(self,name):  
1646       """       """
1647       returns a new coefficient appropriate for coefficient name:       Sets the system status to ``status`` if ``status`` is not present the
1648         current status of the domain is used.
1649       """       """
1650       return escript.Data(0,self.getShapeOfCoefficientOfPDE(name),self.getFunctionSpaceForCoefficientOfPDE(name))       if status == None:
1651                   self.__system_status=self.getDomainStatus()
1652     def getShapeOfCoefficientOfPDE(self,name):       else:
1653             self.__system_status=status
1654    
1655       def getDim(self):
1656       """       """
1657       return the shape of the coefficient name       Returns the spatial dimension of the PDE.
1658    
1659       @param name:       :return: the spatial dimension of the PDE domain
1660         :rtype: ``int``
1661       """       """
1662       if self.hasPDECoefficient(name):       return self.getDomain().getDim()
         return self.__COEFFICIENTS[name].buildShape(self.getNumEquations(),self.getNumSolutions(),self.getDomain().getDim())  
      else:  
         raise ValueError,"Unknown coefficient %s requested"%name  
1663    
1664     def getFunctionSpaceForCoefficientOfPDE(self,name):     def getNumEquations(self):
1665       """       """
1666       return the atoms of the coefficient name       Returns the number of equations.
1667    
1668         :return: the number of equations
1669         :rtype: ``int``
1670         :raise UndefinedPDEError: if the number of equations is not specified yet
1671         """
1672         if self.__numEquations==None:
1673             if self.__numSolutions==None:
1674                raise UndefinedPDEError,"Number of equations is undefined. Please specify argument numEquations."
1675             else:
1676                self.__numEquations=self.__numSolutions
1677         return self.__numEquations
1678    
1679       @param name:     def getNumSolutions(self):
1680       """       """
1681       if self.hasPDECoefficient(name):       Returns the number of unknowns.
         return self.__COEFFICIENTS[name].getFunctionSpace(self.getDomain())  
      else:  
         raise ValueError,"unknown coefficient %s requested"%name  
1682    
1683     def getShapeOfCoefficient(self,name):       :return: the number of unknowns
1684         :rtype: ``int``
1685         :raise UndefinedPDEError: if the number of unknowns is not specified yet
1686         """
1687         if self.__numSolutions==None:
1688            if self.__numEquations==None:
1689                raise UndefinedPDEError,"Number of solution is undefined. Please specify argument numSolutions."
1690            else:
1691                self.__numSolutions=self.__numEquations
1692         return self.__numSolutions
1693    
1694       def reduceEquationOrder(self):
1695       """       """
1696       return the shape of the coefficient name       Returns the status of order reduction for the equation.
1697    
1698       @param name:       :return: True if reduced interpolation order is used for the
1699                  representation of the equation, False otherwise
1700         :rtype: `bool`
1701       """       """
1702       if self.hasCoefficient(name):       return self.__reduce_equation_order
         return self.COEFFICIENTS[name].buildShape(self.getNumEquations(),self.getNumSolutions(),self.getDomain().getDim())  
      else:  
         raise ValueError,"Unknown coefficient %s requested"%name  
1703    
1704     def getFunctionSpaceForCoefficient(self,name):     def reduceSolutionOrder(self):
1705       """       """
1706       return the atoms of the coefficient name       Returns the status of order reduction for the solution.
1707    
1708       @param name:       :return: True if reduced interpolation order is used for the
1709                  representation of the solution, False otherwise
1710         :rtype: `bool`
1711       """       """
1712       if self.hasCoefficient(name):       return self.__reduce_solution_order
         return self.COEFFICIENTS[name].getFunctionSpace(self.getDomain())  
      else:  
         raise ValueError,"unknown coefficient %s requested"%name  
1713    
1714     def alteredCoefficient(self,name):     def getFunctionSpaceForEquation(self):
1715       """       """
1716       announce that coefficient name has been changed       Returns the `FunctionSpace` used to discretize
1717         the equation.
1718    
1719       @param name:       :return: representation space of equation
1720         :rtype: `FunctionSpace`
1721       """       """
1722       if self.hasCoefficient(name):       if self.reduceEquationOrder():
1723          if self.COEFFICIENTS[name].isAlteringOperator(): self.__rebuildOperator()           return escript.ReducedSolution(self.getDomain())
         if self.COEFFICIENTS[name].isAlteringRightHandSide(): self.__rebuildRightHandSide()  
1724       else:       else:
1725          raise ValueError,"unknown coefficient %s requested"%name           return escript.Solution(self.getDomain())
1726    
1727     # ===== debug ==============================================================     def getFunctionSpaceForSolution(self):
1728     def setDebugOn(self):       """
1729         """       Returns the `FunctionSpace` used to represent
1730         """       the solution.
        self.__debug=not None  
1731    
1732     def setDebugOff(self):       :return: representation space of solution
1733         """       :rtype: `FunctionSpace`
1734         """       """
1735         self.__debug=None       if self.reduceSolutionOrder():
1736             return escript.ReducedSolution(self.getDomain())
1737         else:
1738             return escript.Solution(self.getDomain())
1739    
1740     def debug(self):     # ==========================================================================
1741       #   solver settings:
1742       # ==========================================================================
1743       def setSolverOptions(self,options=None):
1744           """
1745           Sets the solver options.
1746    
1747           :param options: the new solver options. If equal ``None``, the solver options are set to the default.
1748           :type options: `SolverOptions` or ``None``
1749           :note: The symmetry flag of options is overwritten by the symmetry flag of the `LinearProblem`.
1750           """
1751           if options==None:
1752              self.__solver_options=SolverOptions()
1753           elif isinstance(options, SolverOptions):
1754              self.__solver_options=options
1755           else:
1756              raise ValueError,"options must be a SolverOptions object."
1757           self.__solver_options.setSymmetry(self.__sym)
1758        
1759       def getSolverOptions(self):
1760         """         """
1761         returns true if the PDE is in the debug mode         Returns the solver options
1762      
1763           :rtype: `SolverOptions`
1764         """         """
1765         return self.__debug         self.__solver_options.setSymmetry(self.__sym)
1766           return self.__solver_options
1767     #===== Lumping ===========================        
    def setLumpingOn(self):  
       """  
       indicates to use matrix lumping  
       """  
       if not self.isUsingLumping():  
          if self.debug() : print "PDE Debug: lumping is set on"  
          self.__rebuildOperator()  
          self.__lumping=True  
   
    def setLumpingOff(self):  
       """  
       switches off matrix lumping  
       """  
       if self.isUsingLumping():  
          if self.debug() : print "PDE Debug: lumping is set off"  
          self.__rebuildOperator()  
          self.__lumping=False  
   
    def setLumping(self,flag=False):  
       """  
       set the matrix lumping flag to flag  
       """  
       if flag:  
          self.setLumpingOn()  
       else:  
          self.setLumpingOff()  
   
1768     def isUsingLumping(self):     def isUsingLumping(self):
1769        """        """
1770                Checks if matrix lumping is the current solver method.
       """  
       return self.__lumping  
   
    #============ method business =========================================================  
    def setSolverMethod(self,solver=util.DEFAULT_METHOD):  
        """  
        sets a new solver  
        """  
        if not solver==self.getSolverMethod():  
            self.__solver_method=solver  
            if self.debug() : print "PDE Debug: New solver is %s"%solver  
            self.__checkMatrixType()  
1771    
1772     def getSolverMethod(self):        :return: True if the current solver method is lumping
1773         """        :rtype: ``bool``
1774         returns the solver method        """
1775         """        return self.getSolverOptions().getSolverMethod() in [ SolverOptions.ROWSUM_LUMPING, SolverOptions.HRZ_LUMPING ]
1776         return self.__solver_method     # ==========================================================================
1777       #    symmetry  flag:
1778     #============ tolerance business =========================================================     # ==========================================================================
    def setTolerance(self,tol=1.e-8):  
        """  
        resets the tolerance to tol.  
        """  
        if not tol>0:  
            raise ValueException,"Tolerance as to be positive"  
        if tol<self.getTolerance(): self.__rebuildSolution()  
        if self.debug() : print "PDE Debug: New tolerance %e",tol  
        self.__tolerance=tol  
        return  
    def getTolerance(self):  
        """  
        returns the tolerance set for the solution  
        """  
        return self.__tolerance  
   
    #===== symmetry  flag ==========================  
1779     def isSymmetric(self):     def isSymmetric(self):
1780        """        """
1781        returns true is the operator is considered to be symmetric        Checks if symmetry is indicated.
1782    
1783          :return: True if a symmetric PDE is indicated, False otherwise
1784          :rtype: ``bool``
1785          :note: the method is equivalent to use getSolverOptions().isSymmetric()
1786        """        """
1787        return self.__sym        self.getSolverOptions().isSymmetric()
1788    
1789     def setSymmetryOn(self):     def setSymmetryOn(self):
1790        """        """
1791        sets the symmetry flag to true        Sets the symmetry flag.
1792          :note: The method overwrites the symmetry flag set by the solver options
1793        """        """
1794        if not self.isSymmetric():        self.__sym=True
1795           if self.debug() : print "PDE Debug: Operator is set to be symmetric"        self.getSolverOptions().setSymmetryOn()
          self.__sym=True  
          self.__checkMatrixType()  
1796    
1797     def setSymmetryOff(self):     def setSymmetryOff(self):
1798        """        """
1799        sets the symmetry flag to false        Clears the symmetry flag.
1800          :note: The method overwrites the symmetry flag set by the solver options
1801        """        """
1802        if self.isSymmetric():        self.__sym=False
1803           if self.debug() : print "PDE Debug: Operator is set to be unsymmetric"        self.getSolverOptions().setSymmetryOff()
          self.__sym=False  
          self.__checkMatrixType()  
   
    def setSymmetryTo(self,flag=False):  
      """  
      sets the symmetry flag to flag  
1804    
1805       @param flag:     def setSymmetry(self,flag=False):
1806       """        """
1807       if flag:        Sets the symmetry flag to ``flag``.
         self.setSymmetryOn()  
      else:  
         self.setSymmetryOff()  
1808    
1809     #===== order reduction ==========================        :param flag: If True, the symmetry flag is set otherwise reset.
1810          :type flag: ``bool``
1811          :note: The method overwrites the symmetry flag set by the solver options
1812          """
1813          self.getSolverOptions().setSymmetry(flag)
1814       # ==========================================================================
1815       # function space handling for the equation as well as the solution
1816       # ==========================================================================
1817     def setReducedOrderOn(self):     def setReducedOrderOn(self):
1818       """       """
1819       switches to on reduced order       Switches reduced order on for solution and equation representation.
1820    
1821         :raise RuntimeError: if order reduction is altered after a coefficient has
1822                              been set
1823       """       """
1824       self.setReducedOrderForSolutionOn()       self.setReducedOrderForSolutionOn()
1825       self.setReducedOrderForEquationOn()       self.setReducedOrderForEquationOn()
1826    
1827     def setReducedOrderOff(self):     def setReducedOrderOff(self):
1828       """       """
1829       switches to full order       Switches reduced order off for solution and equation representation
1830    
1831         :raise RuntimeError: if order reduction is altered after a coefficient has
1832                              been set
1833       """       """
1834       self.setReducedOrderForSolutionOff()       self.setReducedOrderForSolutionOff()
1835       self.setReducedOrderForEquationOff()       self.setReducedOrderForEquationOff()
1836    
1837     def setReducedOrderTo(self,flag=False):     def setReducedOrderTo(self,flag=False):
1838       """       """
1839       sets order according to flag       Sets order reduction state for both solution and equation representation
1840         according to flag.
1841    
1842       @param flag:       :param flag: if True, the order reduction is switched on for both solution
1843                      and equation representation, otherwise or if flag is not
1844                      present order reduction is switched off
1845         :type flag: ``bool``
1846         :raise RuntimeError: if order reduction is altered after a coefficient has
1847                              been set
1848       """       """
1849       self.setReducedOrderForSolutionTo(flag)       self.setReducedOrderForSolutionTo(flag)
1850       self.setReducedOrderForEquationTo(flag)       self.setReducedOrderForEquationTo(flag)
                                                                                                                                                             
1851    
1852     #===== order reduction solution ==========================  
1853     def setReducedOrderForSolutionOn(self):     def setReducedOrderForSolutionOn(self):
1854       """       """
1855       switches to reduced order to interpolate solution       Switches reduced order on for solution representation.
1856    
1857         :raise RuntimeError: if order reduction is altered after a coefficient has
1858                              been set
1859       """       """
1860       new_fs=escript.ReducedSolution(self.getDomain())       if not self.__reduce_solution_order:
1861       if self.getFunctionSpaceForSolution()!=new_fs:           if self.__altered_coefficients:
1862           if self.debug() : print "PDE Debug: Reduced order is used to interpolate solution."                raise RuntimeError,"order cannot be altered after coefficients have been defined."
1863           self.__column_function_space=new_fs           self.trace("Reduced order is used for solution representation.")
1864           self.__rebuildSystem(deep=True)           self.__reduce_solution_order=True
1865             self.initializeSystem()
1866    
1867     def setReducedOrderForSolutionOff(self):     def setReducedOrderForSolutionOff(self):
1868       """       """
1869       switches to full order to interpolate solution       Switches reduced order off for solution representation
1870    
1871         :raise RuntimeError: if order reduction is altered after a coefficient has
1872                              been set.
1873       """       """
1874       new_fs=escript.Solution(self.getDomain())       if self.__reduce_solution_order:
1875       if self.getFunctionSpaceForSolution()!=new_fs:           if self.__altered_coefficients:
1876           if self.debug() : print "PDE Debug: Full order is used to interpolate solution."                raise RuntimeError,"order cannot be altered after coefficients have been defined."
1877           self.__column_function_space=new_fs           self.trace("Full order is used to interpolate solution.")
1878           self.__rebuildSystem(deep=True)           self.__reduce_solution_order=False
1879             self.initializeSystem()
1880    
1881     def setReducedOrderForSolutionTo(self,flag=False):     def setReducedOrderForSolutionTo(self,flag=False):
1882       """       """
1883       sets order for test functions according to flag       Sets order reduction state for solution representation according to flag.
1884    
1885       @param flag:       :param flag: if flag is True, the order reduction is switched on for
1886                      solution representation, otherwise or if flag is not present
1887                      order reduction is switched off
1888         :type flag: ``bool``
1889         :raise RuntimeError: if order reduction is altered after a coefficient has
1890                              been set
1891       """       """
1892       if flag:       if flag:
1893          self.setReducedOrderForSolutionOn()          self.setReducedOrderForSolutionOn()
1894       else:       else:
1895          self.setReducedOrderForSolutionOff()          self.setReducedOrderForSolutionOff()
1896                                                                                                                                                              
    #===== order reduction equation ==========================  
1897     def setReducedOrderForEquationOn(self):     def setReducedOrderForEquationOn(self):
1898       """       """
1899       switches to reduced order for test functions       Switches reduced order on for equation representation.
1900    
1901         :raise RuntimeError: if order reduction is altered after a coefficient has
1902                              been set
1903       """       """
1904       new_fs=escript.ReducedSolution(self.getDomain())       if not self.__reduce_equation_order:
1905       if self.getFunctionSpaceForEquation()!=new_fs:           if self.__altered_coefficients:
1906           if self.debug() : print "PDE Debug: Reduced order is used for test functions."                raise RuntimeError,"order cannot be altered after coefficients have been defined."
1907           self.__row_function_space=new_fs           self.trace("Reduced order is used for test functions.")
1908           self.__rebuildSystem(deep=True)           self.__reduce_equation_order=True
1909             self.initializeSystem()
1910    
1911     def setReducedOrderForEquationOff(self):     def setReducedOrderForEquationOff(self):
1912       """       """
1913       switches to full order for test functions       Switches reduced order off for equation representation.
1914    
1915         :raise RuntimeError: if order reduction is altered after a coefficient has
1916                              been set
1917       """       """
1918       new_fs=escript.Solution(self.getDomain())       if self.__reduce_equation_order:
1919       if self.getFunctionSpaceForEquation()!=new_fs:           if self.__altered_coefficients:
1920           if self.debug() : print "PDE Debug: Full order is used for test functions."                raise RuntimeError,"order cannot be altered after coefficients have been defined."
1921           self.__row_function_space=new_fs           self.trace("Full order is used for test functions.")
1922           self.__rebuildSystem(deep=True)           self.__reduce_equation_order=False
1923             self.initializeSystem()
1924    
1925     def setReducedOrderForEquationTo(self,flag=False):     def setReducedOrderForEquationTo(self,flag=False):
1926       """       """
1927       sets order for test functions according to flag       Sets order reduction state for equation representation according to flag.
1928    
1929       @param flag:       :param flag: if flag is True, the order reduction is switched on for
1930                      equation representation, otherwise or if flag is not present
1931                      order reduction is switched off
1932         :type flag: ``bool``
1933         :raise RuntimeError: if order reduction is altered after a coefficient has
1934                              been set
1935       """       """
1936       if flag:       if flag:
1937          self.setReducedOrderForEquationOn()          self.setReducedOrderForEquationOn()
1938       else:       else:
1939          self.setReducedOrderForEquationOff()          self.setReducedOrderForEquationOff()
1940                                                                                                                                                                 def getOperatorType(self):
1941     # ==== initialization =====================================================================        """
1942     def __getNewOperator(self):        Returns the current system type.
1943          """
1944          return self.__operator_type
1945    
1946       def checkSymmetricTensor(self,name,verbose=True):
1947          """
1948          Tests a coefficient for symmetry.
1949    
1950          :param name: name of the coefficient
1951          :type name: ``str``
1952          :param verbose: if set to True or not present a report on coefficients
1953                          which break the symmetry is printed.
1954          :type verbose: ``bool``
1955          :return: True if coefficient ``name`` is symmetric
1956          :rtype: ``bool``
1957          """
1958          SMALL_TOLERANCE=util.EPSILON*10.
1959          A=self.getCoefficient(name)
1960          verbose=verbose or self.__debug
1961          out=True
1962          if not A.isEmpty():
1963             tol=util.Lsup(A)*SMALL_TOLERANCE
1964             s=A.getShape()
1965             if A.getRank() == 4:
1966                if s[0]==s[2] and s[1] == s[3]:
1967                   for i in range(s[0]):
1968                      for j in range(s[1]):
1969                         for k in range(s[2]):
1970                            for l in range(s[3]):
1971                                if util.Lsup(A[i,j,k,l]-A[k,l,i,j])>tol:
1972                                   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)
1973                                   out=False
1974                else:
1975                     if verbose: print "non-symmetric problem because of inappropriate shape %s of coefficient %s."%(s,name)
1976                     out=False
1977             elif A.getRank() == 2:
1978                if s[0]==s[1]:
1979                   for j in range(s[0]):
1980                      for l in range(s[1]):
1981                         if util.Lsup(A[j,l]-A[l,j])>tol:
1982                            if verbose: print "non-symmetric problem because %s[%d,%d]!=%s[%d,%d]"%(name,j,l,name,l,j)
1983                            out=False
1984                else:
1985                     if verbose: print "non-symmetric problem because of inappropriate shape %s of coefficient %s."%(s,name)
1986                     out=False
1987             elif A.getRank() == 0:
1988                pass
1989             else:
1990                 raise ValueError,"Cannot check rank %s of %s."%(A.getRank(),name)
1991          return out
1992    
1993       def checkReciprocalSymmetry(self,name0,name1,verbose=True):
1994          """
1995          Tests two coefficients for reciprocal symmetry.
1996    
1997          :param name0: name of the first coefficient
1998          :type name0: ``str``
1999          :param name1: name of the second coefficient
2000          :type name1: ``str``
2001          :param verbose: if set to True or not present a report on coefficients
2002                          which break the symmetry is printed
2003          :type verbose: ``bool``
2004          :return: True if coefficients ``name0`` and ``name1`` are reciprocally
2005                   symmetric.
2006          :rtype: ``bool``
2007          """
2008          SMALL_TOLERANCE=util.EPSILON*10.
2009          B=self.getCoefficient(name0)
2010          C=self.getCoefficient(name1)
2011          verbose=verbose or self.__debug
2012          out=True
2013          if B.isEmpty() and not C.isEmpty():
2014             if verbose: print "non-symmetric problem because %s is not present but %s is"%(name0,name1)
2015             out=False
2016          elif not B.isEmpty() and C.isEmpty():
2017             if verbose: print "non-symmetric problem because %s is not present but %s is"%(name0,name1)
2018             out=False
2019          elif not B.isEmpty() and not C.isEmpty():
2020             sB=B.getShape()
2021             sC=C.getShape()
2022             tol=(util.Lsup(B)+util.Lsup(C))*SMALL_TOLERANCE/2.
2023             if len(sB) != len(sC):
2024                 if verbose: print "non-symmetric problem because ranks of %s (=%s) and %s (=%s) are different."%(name0,len(sB),name1,len(sC))
2025                 out=False
2026             else:
2027                 if len(sB)==0:
2028                   if util.Lsup(B-C)>tol:
2029                      if verbose: print "non-symmetric problem because %s!=%s"%(name0,name1)
2030                      out=False
2031                 elif len(sB)==1:
2032                   if sB[0]==sC[0]:
2033                      for j in range(sB[0]):
2034                         if util.Lsup(B[j]-C[j])>tol:
2035                            if verbose: print "non-symmetric PDE because %s[%d]!=%s[%d]"%(name0,j,name1,j)
2036                            out=False
2037                   else:
2038                     if verbose: print "non-symmetric problem because of inappropriate shapes %s and %s of coefficients %s and %s, respectively."%(sB,sC,name0,name1)
2039                 elif len(sB)==3:
2040                   if sB[0]==sC[1] and sB[1]==sC[2] and sB[2]==sC[0]:
2041                       for i in range(sB[0]):
2042                          for j in range(sB[1]):
2043                             for k in range(sB[2]):
2044                                if util.Lsup(B[i,j,k]-C[k,i,j])>tol:
2045                                     if verbose: print "non-symmetric problem because %s[%d,%d,%d]!=%s[%d,%d,%d]"%(name0,i,j,k,name1,k,i,j)
2046                                     out=False
2047                   else:
2048                     if verbose: print "non-symmetric problem because of inappropriate shapes %s and %s of coefficients %s and %s, respectively."%(sB,sC,name0,name1)
2049                 else:
2050                     raise ValueError,"Cannot check rank %s of %s and %s."%(len(sB),name0,name1)
2051          return out
2052    
2053       def getCoefficient(self,name):
2054         """
2055         Returns the value of the coefficient ``name``.
2056    
2057         :param name: name of the coefficient requested
2058         :type name: ``string``
2059         :return: the value of the coefficient
2060         :rtype: `Data`
2061         :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
2062         """
2063         if self.hasCoefficient(name):
2064             return self.__COEFFICIENTS[name].getValue()
2065         else:
2066            raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
2067    
2068       def hasCoefficient(self,name):
2069         """
2070         Returns True if ``name`` is the name of a coefficient.
2071    
2072         :param name: name of the coefficient enquired
2073         :type name: ``string``
2074         :return: True if ``name`` is the name of a coefficient of the general PDE,
2075                  False otherwise
2076         :rtype: ``bool``
2077         """
2078         return self.__COEFFICIENTS.has_key(name)
2079    
2080       def createCoefficient(self, name):
2081         """
2082         Creates a `Data` object corresponding to coefficient
2083         ``name``.
2084    
2085         :return: the coefficient ``name`` initialized to 0
2086         :rtype: `Data`
2087         :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
2088         """
2089         if self.hasCoefficient(name):
2090            return escript.Data(0.,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))
2091         else:
2092            raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
2093    
2094       def getFunctionSpaceForCoefficient(self,name):
2095         """
2096         Returns the `FunctionSpace` to be used for
2097         coefficient ``name``.
2098    
2099         :param name: name of the coefficient enquired
2100         :type name: ``string``
2101         :return: the function space to be used for coefficient ``name``
2102         :rtype: `FunctionSpace`
2103         :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
2104         """
2105         if self.hasCoefficient(name):
2106            return self.__COEFFICIENTS[name].getFunctionSpace(self.getDomain())
2107         else:
2108            raise ValueError,"unknown coefficient %s requested"%name
2109    
2110       def getShapeOfCoefficient(self,name):
2111         """
2112         Returns the shape of the coefficient ``name``.
2113    
2114         :param name: name of the coefficient enquired
2115         :type name: ``string``
2116         :return: the shape of the coefficient ``name``
2117         :rtype: ``tuple`` of ``int``
2118         :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
2119         """
2120         if self.hasCoefficient(name):
2121            return self.__COEFFICIENTS[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())
2122         else:
2123            raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
2124    
2125       def resetAllCoefficients(self):
2126         """
2127         Resets all coefficients to their default values.
2128         """
2129         for i in self.__COEFFICIENTS.iterkeys():
2130             self.__COEFFICIENTS[i].resetValue()
2131    
2132       def alteredCoefficient(self,name):
2133         """
2134         Announces that coefficient ``name`` has been changed.
2135    
2136         :param name: name of the coefficient affected
2137         :type name: ``string``
2138         :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
2139         :note: if ``name`` is q or r, the method will not trigger a rebuild of the
2140                system as constraints are applied to the solved system.
2141         """
2142         if self.hasCoefficient(name):
2143            self.trace("Coefficient %s has been altered."%name)
2144            if not ((name=="q" or name=="r") and self.isUsingLumping()):
2145               if self.__COEFFICIENTS[name].isAlteringOperator(): self.invalidateOperator()
2146               if self.__COEFFICIENTS[name].isAlteringRightHandSide(): self.invalidateRightHandSide()
2147         else:
2148            raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
2149    
2150       def validSolution(self):
2151         """         """
2152           Marks the solution as valid.
2153         """         """
2154         return self.getDomain().newOperator( \         self.__is_solution_valid=True
                            self.getNumEquations(), \  
                            self.getFunctionSpaceForEquation(), \  
                            self.getNumSolutions(), \  
                            self.getFunctionSpaceForSolution(), \  
                            self.__matrix_type)  
2155    
2156     def __makeFreshRightHandSide(self):     def invalidateSolution(self):
2157         """         """
2158           Indicates the PDE has to be resolved if the solution is requested.
2159         """         """
2160         if self.debug() : print "PDE Debug: New right hand side allocated"         self.trace("System will be resolved.")
2161         if self.getNumEquations()>1:         self.__is_solution_valid=False
            self.__righthandside=escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)  
        else:  
            self.__righthandside=escript.Data(0.,(),self.getFunctionSpaceForEquation(),True)  
        return self.__righthandside  
2162    
2163     def __getNewSolution(self):     def isSolutionValid(self):
2164         """         """
2165           Returns True if the solution is still valid.
2166         """         """
2167         if self.debug() : print "PDE Debug: New right hand side allocated"         if not self.getDomainStatus()==self.getSystemStatus(): self.invalidateSolution()
2168         if self.getNumSolutions()>1:         if self.__solution_rtol>self.getSolverOptions().getTolerance() or \
2169             return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)            self.__solution_atol>self.getSolverOptions().getAbsoluteTolerance():
2170         else:           self.invalidateSolution()  
2171             return escript.Data(0.,(),self.getFunctionSpaceForSolution(),True)         return self.__is_solution_valid
2172    
2173     def __makeFreshOperator(self):     def validOperator(self):
2174         """         """
2175           Marks the operator as valid.
2176         """         """
2177         if self.__operator.isEmpty():         self.__is_operator_valid=True
            self.__operator=self.__getNewOperator()  
            if self.debug() : print "PDE Debug: New operator allocated"  
        else:  
            self.__operator.setValue(0.)  
            self.__operator.resetSolver()  
            if self.debug() : print "PDE Debug: Operator reset to zero"  
        return self.__operator  
2178    
2179     #============ some serivice functions  =====================================================     def invalidateOperator(self):
2180     def getDomain(self):         """
2181       """         Indicates the operator has to be rebuilt next time it is used.
2182       returns the domain of the PDE         """
2183       """         self.trace("Operator will be rebuilt.")
2184       return self.__domain         self.invalidateSolution()
2185           self.__is_operator_valid=False
2186    
2187     def getDim(self):     def isOperatorValid(self):
2188       """         """
2189       returns the spatial dimension of the PDE         Returns True if the operator is still valid.
2190       """         """
2191       return self.getDomain().getDim()         if not self.getDomainStatus()==self.getSystemStatus(): self.invalidateOperator()
2192           if not self.getRequiredOperatorType()==self.getOperatorType(): self.invalidateOperator()
2193           return self.__is_operator_valid
2194    
2195     def getNumEquations(self):     def validRightHandSide(self):
2196           """
2197           Marks the right hand side as valid.
2198           """
2199           self.__is_RHS_valid=True
2200    
2201       def invalidateRightHandSide(self):
2202           """
2203           Indicates the right hand side has to be rebuilt next time it is used.
2204           """
2205           self.trace("Right hand side has to be rebuilt.")
2206           self.invalidateSolution()
2207           self.__is_RHS_valid=False
2208    
2209       def isRightHandSideValid(self):
2210           """
2211           Returns True if the operator is still valid.
2212           """
2213           if not self.getDomainStatus()==self.getSystemStatus(): self.invalidateRightHandSide()
2214           return self.__is_RHS_valid
2215    
2216       def invalidateSystem(self):
2217           """
2218           Announces that everything has to be rebuilt.
2219           """
2220           self.invalidateSolution()
2221           self.invalidateOperator()
2222           self.invalidateRightHandSide()
2223    
2224       def isSystemValid(self):
2225           """
2226           Returns True if the system (including solution) is still vaild.
2227           """
2228           return self.isSolutionValid() and self.isOperatorValid() and self.isRightHandSideValid()
2229    
2230       def initializeSystem(self):
2231           """
2232           Resets the system clearing the operator, right hand side and solution.
2233           """
2234           self.trace("New System has been created.")
2235           self.__operator_type=None
2236           self.setSystemStatus()
2237           self.__operator=escript.Operator()
2238           self.__righthandside=escript.Data()
2239           self.__solution=escript.Data()
2240           self.invalidateSystem()
2241    
2242       def getOperator(self):
2243       """       """
2244       returns the number of equations       Returns the operator of the linear problem.
2245    
2246         :return: the operator of the problem
2247       """       """
2248       if self.__numEquations>0:       return self.getSystem()[0]
          return self.__numEquations  
      else:  
          raise ValueError,"Number of equations is undefined. Please specify argument numEquations."  
2249    
2250     def getNumSolutions(self):     def getRightHandSide(self):
2251       """       """
2252       returns the number of unknowns       Returns the right hand side of the linear problem.
2253    
2254         :return: the right hand side of the problem
2255         :rtype: `Data`
2256       """       """
2257       if self.__numSolutions>0:       return self.getSystem()[1]
         return self.__numSolutions  
      else:  
         raise ValueError,"Number of solution is undefined. Please specify argument numSolutions."  
2258    
2259       def createRightHandSide(self):
2260           """
2261           Returns an instance of a new right hand side.
2262           """
2263           self.trace("New right hand side is allocated.")
2264           if self.getNumEquations()>1:
2265               return escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)
2266           else:
2267               return escript.Data(0.,(),self.getFunctionSpaceForEquation(),True)
2268    
2269     def checkSymmetry(self,verbose=True):     def createSolution(self):
2270        """         """
2271        returns if the Operator is symmetric. This is a very expensive operation!!! The symmetry flag is not altered.         Returns an instance of a new solution.
2272        """         """
2273        verbose=verbose or self.debug()         self.trace("New solution is allocated.")
2274        out=True         if self.getNumSolutions()>1:
2275        if self.getNumSolutions()!=self.getNumEquations():             return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)
2276           if verbose: print "non-symmetric PDE because of different number of equations and solutions"         else:
2277           out=False             return escript.Data(0.,(),self.getFunctionSpaceForSolution(),True)
       else:  
          A=self.getCoefficientOfPDE("A")  
          if not A.isEmpty():  
             tol=util.Lsup(A)*self.TOL  
             if self.getNumSolutions()>1:  
                for i in range(self.getNumEquations()):  
                   for j in range(self.getDim()):  
                      for k in range(self.getNumSolutions()):  
                         for l in range(self.getDim()):  
                             if util.Lsup(A[i,j,k,l]-A[k,l,i,j])>tol:  
                                if verbose: print "non-symmetric PDE because A[%d,%d,%d,%d]!=A[%d,%d,%d,%d]"%(i,j,k,l,k,l,i,j)  
                                out=False  
             else:  
                for j in range(self.getDim()):  
                   for l in range(self.getDim()):  
                      if util.Lsup(A[j,l]-A[l,j])>tol:  
                         if verbose: print "non-symmetric PDE because A[%d,%d]!=A[%d,%d]"%(j,l,l,j)  
                         out=False  
          B=self.getCoefficientOfPDE("B")  
          C=self.getCoefficientOfPDE("C")  
          if B.isEmpty() and not C.isEmpty():  
             if verbose: print "non-symmetric PDE because B is not present but C is"  
             out=False  
          elif not B.isEmpty() and C.isEmpty():  
             if verbose: print "non-symmetric PDE because C is not present but B is"  
             out=False  
          elif not B.isEmpty() and not C.isEmpty():  
             tol=(util.Lsup(B)+util.Lsup(C))*self.TOL/2.  
             if self.getNumSolutions()>1:  
                for i in range(self.getNumEquations()):  
                    for j in range(self.getDim()):  
                       for k in range(self.getNumSolutions()):  
                          if util.Lsup(B[i,j,k]-C[k,i,j])>tol:  
                               if verbose: print "non-symmetric PDE because B[%d,%d,%d]!=C[%d,%d,%d]"%(i,j,k,k,i,j)  
                               out=False  
             else:  
                for j in range(self.getDim()):  
                   if util.Lsup(B[j]-C[j])>tol:  
                      if verbose: print "non-symmetric PDE because B[%d]!=C[%d]"%(j,j)  
                      out=False  
          if self.getNumSolutions()>1:  
            D=self.getCoefficientOfPDE("D")  
            if not D.isEmpty():  
              tol=util.Lsup(D)*self.TOL  
              for i in range(self.getNumEquations()):  
                 for k in range(self.getNumSolutions()):  
                   if util.Lsup(D[i,k]-D[k,i])>tol:  
                       if verbose: print "non-symmetric PDE because D[%d,%d]!=D[%d,%d]"%(i,k,k,i)  
                       out=False  
               
       return out  
2278    
2279     def getFlux(self,u):     def resetSolution(self):
2280         """         """
2281         returns the flux J_ij for a given u         Sets the solution to zero.
2282           """
2283           if self.__solution.isEmpty():
2284               self.__solution=self.createSolution()
2285           else:
2286               self.__solution.setToZero()
2287               self.trace("Solution is reset to zero.")
2288    
2289         \f[     def setSolution(self,u, validate=True):
2290         J_ij=A_{ijkl}u_{k,l}+B_{ijk}u_k-X_{ij}         """
2291         \f]         Sets the solution assuming that makes the system valid with the tolrance
2292           defined by the solver options
2293           """
2294           if validate:
2295          self.__solution_rtol=self.getSolverOptions().getTolerance()
2296          self.__solution_atol=self.getSolverOptions().getAbsoluteTolerance()
2297          self.validSolution()
2298           self.__solution=u
2299       def getCurrentSolution(self):
2300           """
2301           Returns the solution in its current state.
2302           """
2303           if self.__solution.isEmpty(): self.__solution=self.createSolution()
2304           return self.__solution
2305    
2306         @param u: argument of the operator     def resetRightHandSide(self):
2307         """         """
2308         raise SystemError,"getFlux is not implemented yet"         Sets the right hand side to zero.
2309         return None         """
2310           if self.__righthandside.isEmpty():
2311               self.__righthandside=self.createRightHandSide()
2312           else:
2313               self.__righthandside.setToZero()
2314               self.trace("Right hand side is reset to zero.")
2315    
2316     def applyOperator(self,u):     def getCurrentRightHandSide(self):
2317         """         """
2318         applies the operator of the PDE to a given solution u in weak from         Returns the right hand side in its current state.
2319           """
2320           if self.__righthandside.isEmpty(): self.__righthandside=self.createRightHandSide()
2321           return self.__righthandside
2322    
2323         @param u: argument of the operator     def resetOperator(self):
2324         """         """
2325         return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())         Makes sure that the operator is instantiated and returns it initialized
2326                                                                                                                                                                     with zeros.
    def getResidual(self,u):  
2327         """         """
2328         return the residual of u in the weak from         if self.getOperatorType() == None:
2329               if self.isUsingLumping():
2330                   self.__operator=self.createSolution()
2331               else:
2332                   self.__operator=self.createOperator()
2333           self.__operator_type=self.getRequiredOperatorType()
2334           else:
2335               if self.isUsingLumping():
2336                   self.__operator.setToZero()
2337               else:
2338                   if self.getOperatorType() == self.getRequiredOperatorType():
2339                       self.__operator.resetValues()
2340                   else:
2341                       self.__operator=self.createOperator()
2342                   self.__operator_type=self.getRequiredOperatorType()
2343               self.trace("Operator reset to zero")
2344    
2345         @param u:     def getCurrentOperator(self):
2346           """
2347           Returns the operator in its current state.
2348         """         """
2349         return self.applyOperator(u)-self.getRightHandSide()         return self.__operator
2350    
2351     def __setValue(self,**coefficients):     def setValue(self,**coefficients):
2352        """        """
2353        sets new values to coefficient        Sets new values to coefficients.
2354    
2355        @param coefficients:        :raise IllegalCoefficient: if an unknown coefficient keyword is used
2356        """        """
2357        # check if the coefficients are  legal:        # check if the coefficients are  legal:
2358        for i in coefficients.iterkeys():        for i in coefficients.iterkeys():
2359           if not self.hasCoefficient(i):           if not self.hasCoefficient(i):
2360              raise ValueError,"Attempt to set unknown coefficient %s"%i              raise IllegalCoefficient,"Attempt to set unknown coefficient %s"%i
2361        # if the number of unknowns or equations is still unknown we try to estimate them:        # if the number of unknowns or equations is still unknown we try to estimate them:
2362        if self.__numEquations<1 or self.__numSolutions<1:        if self.__numEquations==None or self.__numSolutions==None:
2363           for i,d in coefficients.iteritems():           for i,d in coefficients.iteritems():
2364              if hasattr(d,"shape"):              if hasattr(d,"shape"):
2365                  s=d.shape                  s=d.shape
2366              elif hasattr(d,"getShape"):              elif hasattr(d,"getShape"):
2367                  s=d.getShape()                  s=d.getShape()
2368              else:              else:
2369                  s=numarray.array(d).shape                  s=numpy.array(d).shape
2370              if s!=None:              if s!=None:
2371                  # get number of equations and number of unknowns:                  # get number of equations and number of unknowns:
2372                  res=self.COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(s,self.getDim())                  res=self.__COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(self.getDomain(),s)
2373                  if res==None:                  if res==None:
2374                      raise ValueError,"Illegal shape %s of coefficient %s"%(s,i)                      raise IllegalCoefficientValue,"Illegal shape %s of coefficient %s"%(s,i)
2375                  else:                  else:
2376                      if self.__numEquations<1: self.__numEquations=res[0]                      if self.__numEquations==None: self.__numEquations=res[0]
2377                      if self.__numSolutions<1: self.__numSolutions=res[1]                      if self.__numSolutions==None: self.__numSolutions=res[1]
2378        if self.__numEquations<1: raise ValueError,"unidententified number of equations"        if self.__numEquations==None: raise UndefinedPDEError,"unidentified number of equations"
2379        if self.__numSolutions<1: raise ValueError,"unidententified number of solutions"        if self.__numSolutions==None: raise UndefinedPDEError,"unidentified number of solutions"
2380        # now we check the shape of the coefficient if numEquations and numSolutions are set:        # now we check the shape of the coefficient if numEquations and numSolutions are set:
2381        for i,d in coefficients.iteritems():        for i,d in coefficients.iteritems():
2382          if d==None:          try:
2383               d2=escript.Data()             self.__COEFFICIENTS[i].setValue(self.getDomain(),
2384          elif isinstance(d,escript.Data):                       self.getNumEquations(),self.getNumSolutions(),
2385               if d.isEmpty():                       self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
2386                  d2=d             self.alteredCoefficient(i)
2387               else:          except IllegalCoefficientFunctionSpace,m:
2388                  d2=escript.Data(d,self.getFunctionSpaceForCoefficient(i))              # if the function space is wrong then we try the reduced version:
2389          else:              i_red=i+"_reduced"
2390                d2=escript.Data(d,self.getFunctionSpaceForCoefficient(i))              if (not i_red in coefficients.keys()) and i_red in self.__COEFFICIENTS.keys():
2391          if not d2.isEmpty():                  try:
2392             if not self.getShapeOfCoefficient(i)==d2.getShape():                      self.__COEFFICIENTS[i_red].setValue(self.getDomain(),
2393                 raise ValueError,"Expected shape for coefficient %s is %s but actual shape is %s."%(i,self.getShapeOfCoefficient(i),d2.getShape())                                                        self.getNumEquations(),self.getNumSolutions(),
2394          # overwrite new values:                                                        self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
2395          if self.debug(): print "PDE Debug: Coefficient %s has been altered."%i                      self.alteredCoefficient(i_red)
2396          self.COEFFICIENTS[i].setValue(d2)                  except IllegalCoefficientValue,m:
2397          self.alteredCoefficient(i)                      raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
2398                          except IllegalCoefficientFunctionSpace,m:
2399        # reset the HomogeneousConstraintFlag:                      raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m))
2400        self.__setHomogeneousConstraintFlag()              else:
2401        if len(coefficients)>0 and not self.isUsingLumping() and not self.__homogeneous_constraint: self.__rebuildSystem()                  raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m))
2402            except IllegalCoefficientValue,m:
2403     def __setHomogeneousConstraintFlag(self):             raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
2404        """        self.__altered_coefficients=True
       checks if the constraints are homogeneous and sets self.__homogeneous_constraint accordingly.  
       """  
       self.__homogeneous_constraint=True  
       q=self.getCoefficientOfPDE("q")  
       r=self.getCoefficientOfPDE("r")  
       if not q.isEmpty() and not r.isEmpty():  
          if (q*r).Lsup()>=1.e-13*r.Lsup(): self.__homogeneous_constraint=False  
       if self.debug():  
            if self.__homogeneous_constraint:  
                print "PDE Debug: Constraints are homogeneous."  
            else:  
                print "PDE Debug: Constraints are inhomogeneous."  
   
2405    
2406     # ==== rebuild switches =====================================================================     # ==========================================================================
2407     def __rebuildSolution(self,deep=False):     # methods that are typically overwritten when implementing a particular
2408         """     # linear problem
2409         indicates the PDE has to be reolved if the solution is requested     # ==========================================================================
2410         """     def getRequiredOperatorType(self):
2411         if self.__solution_isValid and self.debug() : print "PDE Debug: PDE has to be resolved."        """
2412         self.__solution_isValid=False        Returns the system type which needs to be used by the current set up.
        if deep: self.__solution=escript.Data()  
2413    
2414          :note: Typically this method is overwritten when implementing a
2415                 particular linear problem.
2416          """
2417          return None
2418    
2419     def __rebuildOperator(self,deep=False):     def createOperator(self):
2420         """         """
2421         indicates the operator has to be rebuilt next time it is used         Returns an instance of a new operator.
2422    
2423           :note: This method is overwritten when implementing a particular
2424                  linear problem.
2425         """         """
2426         if self.__operator_isValid and self.debug() : print "PDE Debug: Operator has to be rebuilt."         return escript.Operator()
2427         self.__rebuildSolution(deep)  
2428         self.__operator_isValid=False     def checkSymmetry(self,verbose=True):
2429         if deep: self.__operator=escript.Operator()        """
2430          Tests the PDE for symmetry.
2431    
2432          :param verbose: if set to True or not present a report on coefficients
2433                          which break the symmetry is printed
2434          :type verbose: ``bool``
2435          :return: True if the problem is symmetric
2436          :rtype: ``bool``
2437          :note: Typically this method is overwritten when implementing a
2438                 particular linear problem.
2439          """
2440          out=True
2441          return out
2442    
2443     def __rebuildRightHandSide(self,deep=False):     def getSolution(self,**options):
2444         """         """
2445         indicates the right hand side has to be rebuild next time it is used         Returns the solution of the problem.
2446    
2447           :return: the solution
2448           :rtype: `Data`
2449    
2450           :note: This method is overwritten when implementing a particular
2451                  linear problem.
2452         """         """
2453         if self.__righthandside_isValid and self.debug() : print "PDE Debug: Right hand side has to be rebuilt."         return self.getCurrentSolution()
        self.__rebuildSolution(deep)  
        self.__righthandside_isValid=False  
        if deep: self.__righthandside=escript.Data()  
2454    
2455     def __rebuildSystem(self,deep=False):     def getSystem(self):
2456         """         """
2457         annonced that all coefficient name has been changed         Returns the operator and right hand side of the PDE.
2458    
2459           :return: the discrete version of the PDE
2460           :rtype: ``tuple`` of `Operator` and `Data`.
2461    
2462           :note: This method is overwritten when implementing a particular
2463                  linear problem.
2464         """         """
2465         self.__rebuildSolution(deep)         return (self.getCurrentOperator(), self.getCurrentRightHandSide())
2466         self.__rebuildOperator(deep)  
2467         self.__rebuildRightHandSide(deep)  class LinearPDE(LinearProblem):
2468         """
2469     def __checkMatrixType(self):     This class is used to define a general linear, steady, second order PDE
2470       for an unknown function *u* on a given domain defined through a
2471       `Domain` object.
2472    
2473       For a single PDE having a solution with a single component the linear PDE
2474       is defined in the following form:
2475    
2476       *-(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)*
2477    
2478       where *grad(F)* denotes the spatial derivative of *F*. Einstein's
2479       summation convention, ie. summation over indexes appearing twice in a term
2480       of a sum performed, is used.
2481       The coefficients *A*, *B*, *C*, *D*, *X* and *Y* have to be specified
2482       through `Data` objects in `Function` and
2483       the coefficients *A_reduced*, *B_reduced*, *C_reduced*, *D_reduced*,
2484       *X_reduced* and *Y_reduced* have to be specified through
2485       `Data` objects in `ReducedFunction`.
2486       It is also allowed to use objects that can be converted into such
2487       `Data` objects. *A* and *A_reduced* are rank two, *B*,
2488       *C*, *X*, *B_reduced*, *C_reduced* and *X_reduced* are rank one and
2489       *D*, *D_reduced*, *Y* and *Y_reduced* are scalar.
2490    
2491       The following natural boundary conditions are considered:
2492    
2493       *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*
2494    
2495       where *n* is the outer normal field. Notice that the coefficients *A*,
2496       *A_reduced*, *B*, *B_reduced*, *X* and *X_reduced* are defined in the
2497       PDE. The coefficients *d* and *y* are each a scalar in
2498       `FunctionOnBoundary` and the coefficients
2499       *d_reduced* and *y_reduced* are each a scalar in
2500       `ReducedFunctionOnBoundary`.
2501    
2502       Constraints for the solution prescribe the value of the solution at certain
2503       locations in the domain. They have the form
2504    
2505       *u=r* where *q>0*
2506    
2507       *r* and *q* are each scalar where *q* is the characteristic function
2508       defining where the constraint is applied. The constraints override any
2509       other condition set by the PDE or the boundary condition.
2510    
2511       The PDE is symmetrical if
2512    
2513       *A[i,j]=A[j,i]*  and *B[j]=C[j]* and *A_reduced[i,j]=A_reduced[j,i]*
2514       and *B_reduced[j]=C_reduced[j]*
2515    
2516       For a system of PDEs and a solution with several components the PDE has the
2517       form
2518    
2519       *-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]*
2520    
2521       *A* and *A_reduced* are of rank four, *B*, *B_reduced*, *C* and
2522       *C_reduced* are each of rank three, *D*, *D_reduced*, *X_reduced* and
2523       *X* are each of rank two and *Y* and *Y_reduced* are of rank one.
2524       The natural boundary conditions take the form:
2525    
2526       *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]*
2527    
2528       The coefficient *d* is of rank two and *y* is of rank one both in
2529       `FunctionOnBoundary`. The coefficients
2530       *d_reduced* is of rank two and *y_reduced* is of rank one both in
2531       `ReducedFunctionOnBoundary`.
2532    
2533       Constraints take the form
2534    
2535       *u[i]=r[i]*  where  *q[i]>0*
2536    
2537       *r* and *q* are each rank one. Notice that at some locations not
2538       necessarily all components must have a constraint.
2539    
2540       The system of PDEs is symmetrical if
2541    
2542          - *A[i,j,k,l]=A[k,l,i,j]*
2543          - *A_reduced[i,j,k,l]=A_reduced[k,l,i,j]*
2544          - *B[i,j,k]=C[k,i,j]*
2545          - *B_reduced[i,j,k]=C_reduced[k,i,j]*
2546          - *D[i,k]=D[i,k]*
2547          - *D_reduced[i,k]=D_reduced[i,k]*
2548          - *d[i,k]=d[k,i]*
2549          - *d_reduced[i,k]=d_reduced[k,i]*
2550    
2551       `LinearPDE` also supports solution discontinuities over a contact region
2552       in the domain. To specify the conditions across the discontinuity we are
2553       using the generalised flux *J* which, in the case of a system of PDEs
2554       and several components of the solution, is defined as
2555    
2556       *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]*
2557    
2558       For the case of single solution component and single PDE *J* is defined as
2559    
2560       *J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[j]+(B[i]+B_reduced[i])*u-X[i]-X_reduced[i]*
2561    
2562       In the context of discontinuities *n* denotes the normal on the
2563       discontinuity pointing from side 0 towards side 1 calculated from
2564       `FunctionSpace.getNormal` of `FunctionOnContactZero`.
2565       For a system of PDEs the contact condition takes the form
2566    
2567       *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]*
2568    
2569       where *J0* and *J1* are the fluxes on side 0 and side 1 of the
2570       discontinuity, respectively. *jump(u)*, which is the difference of the
2571       solution at side 1 and at side 0, denotes the jump of *u* across
2572       discontinuity along the normal calculated by `jump`.
2573       The coefficient *d_contact* is of rank two and *y_contact* is of rank one
2574       both in `FunctionOnContactZero` or
2575       `FunctionOnContactOne`.
2576       The coefficient *d_contact_reduced* is of rank two and *y_contact_reduced*
2577       is of rank one both in `ReducedFunctionOnContactZero`
2578       or `ReducedFunctionOnContactOne`.
2579       In case of a single PDE and a single component solution the contact
2580       condition takes the form
2581    
2582       *n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)*
2583    
2584       In this case the coefficient *d_contact* and *y_contact* are each scalar
2585       both in `FunctionOnContactZero` or
2586       `FunctionOnContactOne` and the coefficient
2587       *d_contact_reduced* and *y_contact_reduced* are each scalar both in
2588       `ReducedFunctionOnContactZero` or
2589       `ReducedFunctionOnContactOne`.
2590    
2591       Typical usage::
2592    
2593           p = LinearPDE(dom)
2594           p.setValue(A=kronecker(dom), D=1, Y=0.5)
2595           u = p.getSolution()
2596    
2597       """
2598    
2599       def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):
2600       """       """
2601       reassess the matrix type and, if needed, initiates an operator rebuild       Initializes a new linear PDE.
2602    
2603         :param domain: domain of the PDE
2604         :type domain: `Domain`
2605         :param numEquations: number of equations. If ``None`` the number of
2606                              equations is extracted from the PDE coefficients.
2607         :param numSolutions: number of solution components. If ``None`` the number
2608                              of solution components is extracted from the PDE
2609                              coefficients.
2610         :param debug: if True debug information is printed
2611    
2612         """
2613         super(LinearPDE, self).__init__(domain,numEquations,numSolutions,debug)
2614         #
2615         #   the coefficients of the PDE:
2616         #
2617         self.introduceCoefficients(
2618           A=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2619           B=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2620           C=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2621           D=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2622           X=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
2623           Y=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2624           d=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2625           y=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2626           d_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2627           y_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2628           A_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2629           B_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2630           C_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2631           D_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2632           X_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
2633           Y_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2634           d_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2635           y_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2636           d_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2637           y_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2638           r=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.RIGHTHANDSIDE),
2639           q=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.BOTH) )
2640    
2641       def __str__(self):
2642         """
2643         Returns the string representation of the PDE.
2644    
2645         :return: a simple representation of the PDE
2646         :rtype: ``str``
2647       """       """
2648       new_matrix_type=self.getDomain().getSystemMatrixTypeId(self.getSolverMethod(),self.isSymmetric())       return "<LinearPDE %d>"%id(self)
2649       if not new_matrix_type==self.__matrix_type:  
2650           if self.debug() : print "PDE Debug: Matrix type is now %d."%new_matrix_type     def getRequiredOperatorType(self):
2651           self.__matrix_type=new_matrix_type        """
2652           self.__rebuildOperator(deep=True)        Returns the system type which needs to be used by the current set up.
2653          """
2654     #============ assembling =======================================================        if self.isUsingLumping():
2655     def __copyConstraint(self):       return "__ESCRIPT_DATA"
2656        """        else:
2657        copies the constrint condition into u       solver_options=self.getSolverOptions()
2658        """       return self.getDomain().getSystemMatrixTypeId(solver_options.getSolverMethod(), solver_options.getPreconditioner(),solver_options.getPackage(), solver_options.isSymmetric())
2659        if not self.__righthandside.isEmpty():  
2660           q=self.getCoefficientOfPDE("q")     def checkSymmetry(self,verbose=True):
2661           r=self.getCoefficientOfPDE("r")        """
2662           if not q.isEmpty():        Tests the PDE for symmetry.
              if r.isEmpty():  
                 r2=escript.Data(0,self.__righthandside.getShape(),self.__righthandside.getFunctionSpace())  
              else:  
                 r2=escript.Data(r,self.__righthandside.getFunctionSpace())  
              self.__righthandside.copyWithMask(r2,escript.Data(q,self.__righthandside.getFunctionSpace()))  
2663    
2664     def __applyConstraint(self):        :param verbose: if set to True or not present a report on coefficients
2665                          which break the symmetry is printed.
2666          :type verbose: ``bool``
2667          :return: True if the PDE is symmetric
2668          :rtype: `bool`
2669          :note: This is a very expensive operation. It should be used for
2670                 degugging only! The symmetry flag is not altered.
2671          """
2672          out=True
2673          out=out and self.checkSymmetricTensor("A", verbose)
2674          out=out and self.checkSymmetricTensor("A_reduced", verbose)
2675          out=out and self.checkReciprocalSymmetry("B","C", verbose)
2676          out=out and self.checkReciprocalSymmetry("B_reduced","C_reduced", verbose)
2677          out=out and self.checkSymmetricTensor("D", verbose)
2678          out=out and self.checkSymmetricTensor("D_reduced", verbose)
2679          out=out and self.checkSymmetricTensor("d", verbose)
2680          out=out and self.checkSymmetricTensor("d_reduced", verbose)
2681          out=out and self.checkSymmetricTensor("d_contact", verbose)
2682          out=out and self.checkSymmetricTensor("d_contact_reduced", verbose)
2683          return out
2684    
2685       def createOperator(self):
2686         """         """
2687         applies the constraints defined by q and r to the system         Returns an instance of a new operator.
2688         """         """
2689         q=self.getCoefficientOfPDE("q")         optype=self.getRequiredOperatorType()
2690         r=self.getCoefficientOfPDE("r")         self.trace("New operator of type %s is allocated."%optype)
2691         if not q.isEmpty() and not self.__operator.isEmpty():         return self.getDomain().newOperator( \
2692            # q is the row and column mask to indicate where constraints are set:                             self.getNumEquations(), \
2693            row_q=escript.Data(q,self.getFunctionSpaceForEquation())                             self.getFunctionSpaceForEquation(), \
2694            col_q=escript.Data(q,self.getFunctionSpaceForSolution())                             self.getNumSolutions(), \
2695            u=self.__getNewSolution()                             self.getFunctionSpaceForSolution(), \
2696            if r.isEmpty():                             optype)
2697               r_s=self.__getNewSolution()  
2698            else:     def getSolution(self):
2699               r_s=escript.Data(r,self.getFunctionSpaceForSolution())         """
2700            u.copyWithMask(r_s,col_q)         Returns the solution of the PDE.
2701    
2702           :return: the solution
2703           :rtype: `Data`
2704           """
2705           option_class=self.getSolverOptions()
2706           if not self.isSolutionValid():
2707              mat,f=self.getSystem()
2708            if self.isUsingLumping():            if self.isUsingLumping():
2709               self.__operator.copyWithMask(escript.Data(1,q.getShape(),self.getFunctionSpaceForEquation()),row_q)           if not util.inf(abs(mat)) > 0.:
2710             raise ZeroDivisionError,"Lumped mass matrix as zero entry (try order 1 elements or HRZ lumping)."
2711                 self.setSolution(f*1/mat)
2712            else:            else:
2713               if not self.__righthandside.isEmpty(): self.__righthandside-=self.__operator*u               self.trace("PDE is resolved.")
2714               self.__operator.nullifyRowsAndCols(row_q,col_q,1.)               self.trace("solver options: %s"%str(option_class))
2715                 self.setSolution(mat.solve(f,option_class))
2716           return self.getCurrentSolution()
2717    
2718     def getSystem(self):     def getSystem(self):
2719         """         """
2720         return the operator and right hand side of the PDE         Returns the operator and right hand side of the PDE.
2721    
2722           :return: the discrete version of the PDE
2723           :rtype: ``tuple`` of `Operator` and
2724                   `Data`
2725         """         """
2726         if not self.__operator_isValid or not self.__righthandside_isValid:         if not self.isOperatorValid() or not self.isRightHandSideValid():
2727            if self.isUsingLumping():            if self.isUsingLumping():
2728                if not self.__operator_isValid:                if not self.isOperatorValid():
2729                   if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():                   if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():
2730                         raise TypeError,"Lumped matrix requires same order for equations and unknowns"                        raise TypeError,"Lumped matrix requires same order for equations and unknowns"
2731                   if not self.getCoefficientOfPDE("A").isEmpty():                   if not self.getCoefficient("A").isEmpty():
2732                            raise Warning,"Lumped matrix does not allow coefficient A"                        raise ValueError,"coefficient A in lumped matrix may not be present."
2733                   if not self.getCoefficientOfPDE("B").isEmpty():                   if not self.getCoefficient("B").isEmpty():
2734                            raise Warning,"Lumped matrix does not allow coefficient B"                        raise ValueError,"coefficient B in lumped matrix may not be present."
2735                   if not self.getCoefficientOfPDE("C").isEmpty():                   if not self.getCoefficient("C").isEmpty():
2736                            raise Warning,"Lumped matrix does not allow coefficient C"                        raise ValueError,"coefficient C in lumped matrix may not be present."
2737                   if self.debug() : print "PDE Debug: New lumped operator is built."                   if not self.getCoefficient("d_contact").isEmpty():
2738                   mat=self.__getNewOperator()                        raise ValueError,"coefficient d_contact in lumped matrix may not be present."
2739                   self.getDomain().addPDEToSystem(mat,escript.Data(), \                   if not self.getCoefficient("A_reduced").isEmpty():
2740                             self.getCoefficientOfPDE("A"), \                        raise ValueError,"coefficient A_reduced in lumped matrix may not be present."
2741                             self.getCoefficientOfPDE("B"), \                   if not self.getCoefficient("B_reduced").isEmpty():
2742                             self.getCoefficientOfPDE("C"), \                        raise ValueError,"coefficient B_reduced in lumped matrix may not be present."
2743                             self.getCoefficientOfPDE("D"), \                   if not self.getCoefficient("C_reduced").isEmpty():
2744                             escript.Data(), \                        raise ValueError,"coefficient C_reduced in lumped matrix may not be present."
2745                             escript.Data(), \                   if not self.getCoefficient("d_contact_reduced").isEmpty():
2746                             self.getCoefficientOfPDE("d"), \                        raise ValueError,"coefficient d_contact_reduced in lumped matrix may not be present."
2747                             escript.Data(),\                   D=self.getCoefficient("D")
2748                             self.getCoefficientOfPDE("d_contact"), \                   d=self.getCoefficient("d")
2749                             escript.Data())                   D_reduced=self.getCoefficient("D_reduced")
2750                   self.__operator=mat*escript.Data(1,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)                   d_reduced=self.getCoefficient("d_reduced")
2751                   self.__applyConstraint()                   if not D.isEmpty():
2752                   self.__operator_isValid=True                       if self.getNumSolutions()>1:
2753                if not self.__righthandside_isValid:                          D_times_e=util.matrix_mult(D,numpy.ones((self.getNumSolutions(),)))
2754                   if self.debug() : print "PDE Debug: New right hand side is built."                       else:
2755                   self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \                          D_times_e=D
2756                                 self.getCoefficientOfPDE("X"), \                   else:
2757                                 self.getCoefficientOfPDE("Y"),\                      D_times_e=escript.Data()
2758                                 self.getCoefficientOfPDE("y"),\                   if not d.isEmpty():
2759                                 self.getCoefficientOfPDE("y_contact"))                       if self.getNumSolutions()>1:
2760                   self.__copyConstraint()                          d_times_e=util.matrix_mult(d,numpy.ones((self.getNumSolutions(),)))
2761                   self.__righthandside_isValid=True                       else:
2762                            d_times_e=d
2763                     else:
2764                        d_times_e=escript.Data()
2765    
2766                     if not D_reduced.isEmpty():
2767                         if self.getNumSolutions()>1:
2768                            D_reduced_times_e=util.matrix_mult(D_reduced,numpy.ones((self.getNumSolutions(),)))
2769                         else:
2770                            D_reduced_times_e=D_reduced
2771                     else:
2772                        D_reduced_times_e=escript.Data()
2773                     if not d_reduced.isEmpty():
2774                         if self.getNumSolutions()>1:
2775                            d_reduced_times_e=util.matrix_mult(d_reduced,numpy.ones((self.getNumSolutions(),)))
2776                         else:
2777                            d_reduced_times_e=d_reduced
2778                     else:
2779                        d_reduced_times_e=escript.Data()
2780    
2781                     self.resetOperator()
2782                     operator=self.getCurrentOperator()
2783                     if hasattr(self.getDomain(), "addPDEToLumpedSystem") :
2784                hrz_lumping=( self.getSolverOptions().getSolverMethod() ==  SolverOptions.HRZ_LUMPING )
2785                self.getDomain().addPDEToLumpedSystem(operator, D_times_e, d_times_e,  hrz_lumping )
2786                self.getDomain().addPDEToLumpedSystem(operator, D_reduced_times_e, d_reduced_times_e, hrz_lumping)
2787                     else:
2788                        self.getDomain().addPDEToRHS(operator, \
2789                                                     escript.Data(), \
2790                                                     D_times_e, \
2791                                                     d_times_e,\
2792                                                     escript.Data())
2793                        self.getDomain().addPDEToRHS(operator, \
2794                                                     escript.Data(), \
2795                                                     D_reduced_times_e, \
2796                                                     d_reduced_times_e,\
2797                                                     escript.Data())
2798                     self.trace("New lumped operator has been built.")
2799                  if not self.isRightHandSideValid():
2800                     self.resetRightHandSide()
2801                     righthandside=self.getCurrentRightHandSide()
2802                     self.getDomain().addPDEToRHS(righthandside, \
2803                                   self.getCoefficient("X"), \
2804                                   self.getCoefficient("Y"),\
2805                                   self.getCoefficient("y"),\
2806                                   self.getCoefficient("y_contact"))
2807                     self.getDomain().addPDEToRHS(righthandside, \
2808                                   self.getCoefficient("X_reduced"), \
2809                                   self.getCoefficient("Y_reduced"),\
2810                                   self.getCoefficient("y_reduced"),\
2811                                   self.getCoefficient("y_contact_reduced"))
2812                     self.trace("New right hand side has been built.")
2813                     self.validRightHandSide()
2814                  self.insertConstraint(rhs_only=False)
2815                  self.validOperator()
2816            else:            else:
2817               if not self.__operator_isValid and not self.__righthandside_isValid:               if not self.isOperatorValid() and not self.isRightHandSideValid():
2818                   if self.debug() : print "PDE Debug: New system is built."                   self.resetRightHandSide()
2819                   self.getDomain().addPDEToSystem(self.__makeFreshOperator(),self.__makeFreshRightHandSide(), \                   righthandside=self.getCurrentRightHandSide()
2820                                 self.getCoefficientOfPDE("A"), \                   self.resetOperator()
2821                                 self.getCoefficientOfPDE("B"), \                   operator=self.getCurrentOperator()
2822                                 self.getCoefficientOfPDE("C"), \                   self.getDomain().addPDEToSystem(operator,righthandside, \
2823                                 self.getCoefficientOfPDE("D"), \                                 self.getCoefficient("A"), \
2824                                 self.getCoefficientOfPDE("X"), \                                 self.getCoefficient("B"), \
2825                                 self.getCoefficientOfPDE("Y"), \                                 self.getCoefficient("C"), \
2826                                 self.getCoefficientOfPDE("d"), \                                 self.getCoefficient("D"), \
2827                                 self.getCoefficientOfPDE("y"), \                                 self.getCoefficient("X"), \
2828                                 self.getCoefficientOfPDE("d_contact"), \                                 self.getCoefficient("Y"), \
2829                                 self.getCoefficientOfPDE("y_contact"))                                 self.getCoefficient("d"), \
2830                   self.__applyConstraint()                                 self.getCoefficient("y"), \
2831                   self.__copyConstraint()                                 self.getCoefficient("d_contact"), \
2832                   self.__operator_isValid=True                                 self.getCoefficient("y_contact"))
2833                   self.__righthandside_isValid=True                   self.getDomain().addPDEToSystem(operator,righthandside, \
2834               elif not self.__righthandside_isValid:                                 self.getCoefficient("A_reduced"), \
2835                   if self.debug() : print "PDE Debug: New right hand side is built."                                 self.getCoefficient("B_reduced"), \
2836                   self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \                                 self.getCoefficient("C_reduced"), \
2837                                 self.getCoefficientOfPDE("X"), \                                 self.getCoefficient("D_reduced"), \
2838                                 self.getCoefficientOfPDE("Y"),\                                 self.getCoefficient("X_reduced"), \
2839                                 self.getCoefficientOfPDE("y"),\                                 self.getCoefficient("Y_reduced"), \
2840                                 self.getCoefficientOfPDE("y_contact"))                                 self.getCoefficient("d_reduced"), \
2841                   self.__copyConstraint()                                 self.getCoefficient("y_reduced"), \
2842                   self.__righthandside_isValid=True                                 self.getCoefficient("d_contact_reduced"), \
2843               elif not self.__operator_isValid:                                 self.getCoefficient("y_contact_reduced"))
2844                   if self.debug() : print "PDE Debug: New operator is built."                   self.insertConstraint(rhs_only=False)
2845                   self.getDomain().addPDEToSystem(self.__makeFreshOperator(),escript.Data(), \                   self.trace("New system has been built.")
2846                              self.getCoefficientOfPDE("A"), \                   self.validOperator()
2847                              self.getCoefficientOfPDE("B"), \                   self.validRightHandSide()
2848                              self.getCoefficientOfPDE("C"), \               elif not self.isRightHandSideValid():
2849                              self.getCoefficientOfPDE("D"), \                   self.resetRightHandSide()
2850                     righthandside=self.getCurrentRightHandSide()
2851                     self.getDomain().addPDEToRHS(righthandside,
2852                                   self.getCoefficient("X"), \
2853                                   self.getCoefficient("Y"),\
2854                                   self.getCoefficient("y"),\
2855                                   self.getCoefficient("y_contact"))
2856                     self.getDomain().addPDEToRHS(righthandside,
2857                                   self.getCoefficient("X_reduced"), \
2858                                   self.getCoefficient("Y_reduced"),\
2859                                   self.getCoefficient("y_reduced"),\
2860                                   self.getCoefficient("y_contact_reduced"))
2861                     self.insertConstraint(rhs_only=True)
2862                     self.trace("New right hand side has been built.")
2863                     self.validRightHandSide()
2864                 elif not self.isOperatorValid():
2865                     self.resetOperator()
2866                     operator=self.getCurrentOperator()
2867                     self.getDomain().addPDEToSystem(operator,escript.Data(), \
2868                                self.getCoefficient("A"), \
2869                                self.getCoefficient("B"), \
2870                                self.getCoefficient("C"), \
2871                                self.getCoefficient("D"), \
2872                              escript.Data(), \                              escript.Data(), \
2873                              escript.Data(), \                              escript.Data(), \
2874                              self.getCoefficientOfPDE("d"), \                              self.getCoefficient("d"), \
2875                              escript.Data(),\                              escript.Data(),\
2876                              self.getCoefficientOfPDE("d_contact"), \                              self.getCoefficient("d_contact"), \
2877                              escript.Data())                              escript.Data())
2878                   self.__applyConstraint()                   self.getDomain().addPDEToSystem(operator,escript.Data(), \
2879                   self.__operator_isValid=True                              self.getCoefficient("A_reduced"), \
2880         return (self.__operator,self.__righthandside)                              self.getCoefficient("B_reduced"), \
2881     def getOperator(self):                              self.getCoefficient("C_reduced"), \
2882         """                              self.getCoefficient("D_reduced"), \
2883         returns the operator of the PDE                              escript.Data(), \
2884         """                              escript.Data(), \
2885         return self.getSystem()[0]                              self.getCoefficient("d_reduced"), \
2886                                escript.Data(),\
2887     def getRightHandSide(self):                              self.getCoefficient("d_contact_reduced"), \