/[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

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