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