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