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