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