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