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