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