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