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