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