/[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

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