/[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 142 by jgs, Mon Jul 25 05:28:20 2005 UTC trunk/escript/py_src/linearPDEs.py revision 2524 by artak, Wed Jul 8 00:39:26 2009 UTC
# Line 1  Line 1 
 # $Id$  
1    
2  ## @file linearPDEs.py  ########################################################
3    #
4    # Copyright (c) 2003-2008 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-2008 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  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    __author__="Lutz Gross, l.gross@uq.edu.au"
45    
46    
47    class SolverOptions(object):
48        """
49        this class defines the solver options for a linear or non-linear solver.
50        
51        The option also supports the handling of diagnostic informations.
52        
53        Typical usage is
54        
55        opts=SolverOptions()
56        print opts
57        opts.resetDiagnostics()
58        u=solver(opts)
59        print "number of iteration steps: =",opts.getDiagnostics("num_iter")
60        
61    
62        @cvar DEFAULT: The default method used to solve the system of linear equations
63        @cvar DIRECT: The direct solver based on LDU factorization
64        @cvar CHOLEVSKY: The direct solver based on LDLt factorization (can only be applied for symmetric PDEs)
65        @cvar PCG: The preconditioned conjugate gradient method (can only be applied for symmetric PDEs)
66        @cvar CR: The conjugate residual method
67        @cvar CGS: The conjugate gradient square method
68        @cvar BICGSTAB: The stabilized Bi-Conjugate Gradient method
69        @cvar TFQMR: Transport Free Quasi Minimal Residual method
70        @cvar MINRES: Minimum residual method
71        @cvar SSOR: The symmetric over-relaxation method
72        @cvar ILU0: The incomplete LU factorization preconditioner with no fill-in
73        @cvar ILUT: The incomplete LU factorization preconditioner with fill-in
74        @cvar JACOBI: The Jacobi preconditioner
75        @cvar GMRES: The Gram-Schmidt minimum residual method
76        @cvar PRES20: Special GMRES with restart after 20 steps and truncation after 5 residuals
77        @cvar LUMPING: Matrix lumping
78        @cvar NO_REORDERING: No matrix reordering allowed
79        @cvar MINIMUM_FILL_IN: Reorder matrix to reduce fill-in during factorization
80        @cvar NESTED_DISSECTION: Reorder matrix to improve load balancing during factorization
81        @cvar PASO: PASO solver package
82        @cvar SCSL: SGI SCSL solver library
83        @cvar MKL: Intel's MKL solver library
84        @cvar UMFPACK: The UMFPACK library
85        @cvar TRILINOS: The TRILINOS parallel solver class library from Sandia National Labs
86        @cvar ITERATIVE: The default iterative solver
87        @cvar AMG: Algebraic Multi Grid
88        @cvar REC_ILU: recursive ILU0
89        @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  def _CompTuple2(t1,t2):          @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    class IllegalCoefficient(ValueError):
902       """
903       Exception that is raised if an illegal coefficient of the general or
904       particular PDE is requested.
905     """     """
906     Compare two tuples     pass
907    
908     \param t1 The first tuple  class IllegalCoefficientValue(ValueError):
    \param t2 The second tuple  
909     """     """
910       Exception that is raised if an incorrect value for a coefficient is used.
911       """
912       pass
913    
914     dif=t1[0]+t1[1]-(t2[0]+t2[1])  class IllegalCoefficientFunctionSpace(ValueError):
915     if dif<0: return 1     """
916     elif dif>0: return -1     Exception that is raised if an incorrect function space for a coefficient
917     else: return 0     is used.
918       """
919    
920  def ELMAN_RAMAGE(P):  class UndefinedPDEError(ValueError):
921      return (P-1.).wherePositive()*0.5*(1.-1./(P+1.e-15))     """
922       Exception that is raised if a PDE is not fully defined yet.
923       """
924       pass
925    
926  def SIMPLIFIED_BROOK_HUGHES(P):  class PDECoef(object):
927      c=(P-3.).whereNegative()      """
928      return P/6.*c+1./2.*(1.-c)      A class for describing a PDE coefficient.
929    
930  def HALF(P):      @cvar INTERIOR: indicator that coefficient is defined on the interior of
931      return escript.Scalar(0.5,P.getFunctionSpace())                      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    
 class PDECoefficient:  
     """  
     A class for PDE coefficients  
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         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
# Line 64  class PDECoefficient: Line 1006  class PDECoefficient:
1006    
1007      def resetValue(self):      def resetValue(self):
1008         """         """
1009         resets coefficient value to default         Resets the coefficient value to the default.
1010         """         """
1011         self.value=escript.Data()         self.value=escript.Data()
1012    
1013      def getFunctionSpace(self,domain):      def getFunctionSpace(self,domain,reducedEquationOrder=False,reducedSolutionOrder=False):
1014         """         """
1015         defines the FunctionSpace of the coefficient on the domain         Returns the L{FunctionSpace<escript.FunctionSpace>} of the coefficient.
1016    
1017         @param domain:         @param domain: domain on which the PDE uses the coefficient
1018         """         @type domain: L{Domain<escript.Domain>}
1019         if self.what==self.INTERIOR: return escript.Function(domain)         @param reducedEquationOrder: True to indicate that reduced order is used
1020         elif self.what==self.BOUNDARY: return escript.FunctionOnBoundary(domain)                                      to represent the equation
1021         elif self.what==self.CONTACT: return escript.FunctionOnContactZero(domain)         @type reducedEquationOrder: C{bool}
1022         elif self.what==self.CONTINUOUS: return escript.ContinuousFunction(domain)         @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):      def getValue(self):
1049         """         """
1050         returns the value of the coefficient:         Returns the value of the coefficient.
1051    
1052           @return: value of the coefficient
1053           @rtype: L{Data<escript.Data>}
1054         """         """
1055         return self.value         return self.value
1056        
1057      def setValue(self,newValue):      def setValue(self,domain,numEquations=1,numSolutions=1,reducedEquationOrder=False,reducedSolutionOrder=False,newValue=None):
        """  
        set the value of the coefficient to new value  
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         self.value=newValue
1097        
1098      def isAlteringOperator(self):      def isAlteringOperator(self):
1099          """          """
1100      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 102  class PDECoefficient: Line 1110  class PDECoefficient:
1110    
1111      def isAlteringRightHandSide(self):      def isAlteringRightHandSide(self):
1112          """          """
1113      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         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           elif self.definesNumEquation():
1156              for e in range(num,0,-1):
1157                 s=self.getShape(domain,e,0)
1158                 if len(s)==0 and len(shape)==0:
1159                     return (1,None)
1160                 else:
1161                     if s==shape: return (e,None)
1162    
1163           elif self.definesNumSolutions():
1164              for u in range(num,0,-1):
1165                 s=self.getShape(domain,0,u)
1166                 if len(s)==0 and len(shape)==0:
1167                     return (None,1)
1168                 else:
1169                     if s==shape: return (None,u)
1170         return None         return None
1171    
1172      def buildShape(self,e=1,u=1,dim=3):      def definesNumSolutions(self):
1173          """         """
1174      builds the required shape for a given number of equations e, number of unknowns u and spatial dimension dim         Checks if the coefficient allows to estimate the number of solution
1175           components.
1176    
1177      @param e:         @return: True if the coefficient allows an estimate of the number of
1178      @param u:                  solution components, False otherwise
1179      @param dim:         @rtype: C{bool}
1180      """         """
1181          s=()         for i in self.pattern:
1182          for i in self.pattern:               if i==self.BY_SOLUTION: return True
1183               if i==self.EQUATION:         return False
1184                  if e>1: s=s+(e,)  
1185               elif i==self.SOLUTION:      def definesNumEquation(self):
1186                  if u>1: s=s+(u,)         """
1187           Checks if the coefficient allows to estimate the number of equations.
1188    
1189           @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        def __CompTuple2(self,t1,t2):
1198          """
1199          Compares two tuples of possible number of equations and number of
1200          solutions.
1201    
1202          @param t1: the first tuple
1203          @param t2: the second tuple
1204          @return: 0, 1, or -1
1205          """
1206    
1207          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        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:               else:
1233                  s=s+(dim,)                  s=s+(dim,)
1234          return s         return s
1235    
1236  class LinearPDE:  #====================================================================================================================
    """  
    Class to handle a linear PDE  
     
    class to define a linear PDE of the form  
1237    
1238     \f[  class LinearProblem(object):
1239       -(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     """
1240     \f]     This is the base class to define a general linear PDE-type problem for
1241       for an unknown function M{u} on a given domain defined through a
1242       L{Domain<escript.Domain>} object. The problem can be given as a single
1243       equation or as a system of equations.
1244    
1245     with boundary conditons:     The class assumes that some sort of assembling process is required to form
1246       a problem of the form
1247    
1248     \f[     M{L u=f}
    n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_{ik}u_k = - n_j*X_{ij} + y_i  
    \f]  
1249    
1250     and contact conditions     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     \f[     """
1254     n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_contact_{ik}[u_k] = - n_j*X_{ij} + y_contact_i     def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):
1255     \f]       """
1256         Initializes a linear problem.
1257    
1258     and constraints:       @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     \f[       """
1268     u_i=r_i \quad \mathrm{where} \quad q_i>0       super(LinearProblem, self).__init__()
    \f]  
   
    """  
    TOL=1.e-13  
    DEFAULT_METHOD=util.DEFAULT_METHOD  
    DIRECT=util.DIRECT  
    CHOLEVSKY=util.CHOLEVSKY  
    PCG=util.PCG  
    CR=util.CR  
    CGS=util.CGS  
    BICGSTAB=util.BICGSTAB  
    SSOR=util.SSOR  
    GMRES=util.GMRES  
    PRES20=util.PRES20  
    ILU0=util.ILU0  
    JACOBI=util.JACOBI  
   
    def __init__(self,domain,numEquations=0,numSolutions=0):  
      """  
      initializes a new linear PDE.  
   
      @param args:  
      """  
      # COEFFICIENTS can be overwritten by subclasses:  
      self.COEFFICIENTS={  
        "A"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM,PDECoefficient.SOLUTION,PDECoefficient.DIM),PDECoefficient.OPERATOR),  
        "B"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),  
        "C"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION,PDECoefficient.DIM),PDECoefficient.OPERATOR),  
        "D"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),  
        "X"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,PDECoefficient.DIM),PDECoefficient.RIGHTHANDSIDE),  
        "Y"         : PDECoefficient(PDECoefficient.INTERIOR,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),  
        "d"         : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),  
        "y"         : PDECoefficient(PDECoefficient.BOUNDARY,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),  
        "d_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.EQUATION,PDECoefficient.SOLUTION),PDECoefficient.OPERATOR),  
        "y_contact" : PDECoefficient(PDECoefficient.CONTACT,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),  
        "r"         : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.EQUATION,),PDECoefficient.RIGHTHANDSIDE),  
        "q"         : PDECoefficient(PDECoefficient.CONTINUOUS,(PDECoefficient.SOLUTION,),PDECoefficient.BOTH)}  
1269    
1270       # initialize attributes       self.__debug=debug
      self.__debug=None  
1271       self.__domain=domain       self.__domain=domain
1272       self.__numEquations=numEquations       self.__numEquations=numEquations
1273       self.__numSolutions=numSolutions       self.__numSolutions=numSolutions
1274       self.cleanCoefficients()       self.__altered_coefficients=False
1275         self.__reduce_equation_order=False
1276       self.__operator=escript.Operator()       self.__reduce_solution_order=False
      self.__operator_isValid=False  
      self.__righthandside=escript.Data()  
      self.__righthandside_isValid=False  
      self.__solution=escript.Data()  
      self.__solution_isValid=False  
   
      # set some default values:  
      self.__homogeneous_constraint=True  
      self.__row_function_space=escript.Solution(self.__domain)  
      self.__column_function_space=escript.Solution(self.__domain)  
      self.__tolerance=1.e-8  
      self.__solver_method=util.DEFAULT_METHOD  
      self.__matrix_type=self.__domain.getSystemMatrixTypeId(util.DEFAULT_METHOD,False)  
1277       self.__sym=False       self.__sym=False
1278       self.__lumping=False       self.setSolverOptions()
1279         self.__is_RHS_valid=False
1280     def createCoefficient(self, name):       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       create a data object corresponding to coefficient name       Switches debug output on.
      @param name:  
1305       """       """
1306       return escript.Data(shape = getShapeOfCoefficient(name), \       self.__debug=not None
                          what = getFunctionSpaceForCoefficient(name))  
1307    
1308     def __del__(self):     def setDebugOff(self):
1309       pass       """
1310         Switches debug output off.
1311         """
1312         self.__debug=None
1313    
1314     def getCoefficient(self,name):     def setDebug(self, flag):
1315       """       """
1316       return the value of the parameter name       Switches debug output on if C{flag} is True otherwise it is switched off.
1317    
1318       @param name:       @param flag: desired debug status
1319         @type flag: C{bool}
1320       """       """
1321       return self.COEFFICIENTS[name].getValue()       if flag:
1322             self.setDebugOn()
1323         else:
1324             self.setDebugOff()
1325    
1326     def getCoefficientOfPDE(self,name):     def trace(self,text):
1327       """       """
1328       return the value of the coefficient name of the general PDE.       Prints the text message if debug mode is switched on.
      This method is called by the assembling routine it can be  
      overwritten to map coefficients of a particualr PDE to the general PDE.  
1329    
1330       @param name:       @param text: message to be printed
1331         @type text: C{string}
1332       """       """
1333       return self.getCoefficient(name)       if self.__debug: print "%s: %s"%(str(self),text)
1334    
1335     def hasCoefficient(self,name):     # ==========================================================================
1336        """     # some service functions:
1337        return true if name is the name of a coefficient     # ==========================================================================
1338       def introduceCoefficients(self,**coeff):
1339           """
1340           Introduces new coefficients into the problem.
1341    
1342        @param name:         Use::
       """  
       return self.COEFFICIENTS.has_key(name)  
1343    
1344     def getFunctionSpaceForEquation(self):         p.introduceCoefficients(A=PDECoef(...), B=PDECoef(...))
1345    
1346           to introduce the coefficients M{A} ans M{B}.
1347           """
1348           for name, type in coeff.items():
1349               if not isinstance(type,PDECoef):
1350                  raise ValueError,"coefficient %s has no type."%name
1351               self.__COEFFICIENTS[name]=type
1352               self.__COEFFICIENTS[name].resetValue()
1353               self.trace("coefficient %s has been introduced."%name)
1354    
1355       def getDomain(self):
1356       """       """
1357       return true if the test functions should use reduced order       Returns the domain of the PDE.
1358    
1359         @return: the domain of the PDE
1360         @rtype: L{Domain<escript.Domain>}
1361       """       """
1362       return self.__row_function_space       return self.__domain
1363    
1364     def getFunctionSpaceForSolution(self):     def getDim(self):
1365       """       """
1366       return true if the interpolation of the solution should use reduced order       Returns the spatial dimension of the PDE.
1367    
1368         @return: the spatial dimension of the PDE domain
1369         @rtype: C{int}
1370       """       """
1371       return self.__column_function_space       return self.getDomain().getDim()
1372    
1373     def setValue(self,**coefficients):     def getNumEquations(self):
1374        """       """
1375        sets new values to coefficients       Returns the number of equations.
1376    
1377        @param coefficients:       @return: the number of equations
1378        """       @rtype: C{int}
1379        self.__setValue(**coefficients)       @raise UndefinedPDEError: if the number of equations is not specified yet
1380               """
1381         if self.__numEquations==None:
1382             if self.__numSolutions==None:
1383                raise UndefinedPDEError,"Number of equations is undefined. Please specify argument numEquations."
1384             else:
1385                self.__numEquations=self.__numSolutions
1386         return self.__numEquations
1387    
1388     def cleanCoefficients(self):     def getNumSolutions(self):
      """  
      resets all coefficients to default values.  
1389       """       """
1390       for i in self.COEFFICIENTS.iterkeys():       Returns the number of unknowns.
          self.COEFFICIENTS[i].resetValue()  
1391    
1392     def createNewCoefficient(self,name):       @return: the number of unknowns
1393         @rtype: C{int}
1394         @raise UndefinedPDEError: if the number of unknowns is not specified yet
1395         """
1396         if self.__numSolutions==None:
1397            if self.__numEquations==None:
1398                raise UndefinedPDEError,"Number of solution is undefined. Please specify argument numSolutions."
1399            else:
1400                self.__numSolutions=self.__numEquations
1401         return self.__numSolutions
1402    
1403       def reduceEquationOrder(self):
1404       """       """
1405       returns a new coefficient appropriate for coefficient name:       Returns the status of order reduction for the equation.
1406    
1407         @return: True if reduced interpolation order is used for the
1408                  representation of the equation, False otherwise
1409         @rtype: L{bool}
1410       """       """
1411       return escript.Data(0,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))       return self.__reduce_equation_order
         
1412    
1413     def getShapeOfCoefficient(self,name):     def reduceSolutionOrder(self):
1414       """       """
1415       return the shape of the coefficient name       Returns the status of order reduction for the solution.
1416    
1417       @param name:       @return: True if reduced interpolation order is used for the
1418                  representation of the solution, False otherwise
1419         @rtype: L{bool}
1420       """       """
1421       if self.hasCoefficient(name):       return self.__reduce_solution_order
         return self.COEFFICIENTS[name].buildShape(self.getNumEquations(),self.getNumSolutions(),self.getDomain().getDim())  
      else:  
         raise ValueError,"Solution coefficient %s requested"%name  
1422    
1423     def getFunctionSpaceForCoefficient(self,name):     def getFunctionSpaceForEquation(self):
1424       """       """
1425       return the atoms of the coefficient name       Returns the L{FunctionSpace<escript.FunctionSpace>} used to discretize
1426         the equation.
1427    
1428       @param name:       @return: representation space of equation
1429         @rtype: L{FunctionSpace<escript.FunctionSpace>}
1430       """       """
1431       if self.hasCoefficient(name):       if self.reduceEquationOrder():
1432          return self.COEFFICIENTS[name].getFunctionSpace(self.getDomain())           return escript.ReducedSolution(self.getDomain())
1433       else:       else:
1434          raise ValueError,"Solution coefficient %s requested"%name           return escript.Solution(self.getDomain())
1435    
1436     def alteredCoefficient(self,name):     def getFunctionSpaceForSolution(self):
1437       """       """
1438       announce that coefficient name has been changed       Returns the L{FunctionSpace<escript.FunctionSpace>} used to represent
1439         the solution.
1440    
1441       @param name:       @return: representation space of solution
1442         @rtype: L{FunctionSpace<escript.FunctionSpace>}
1443       """       """
1444       if self.hasCoefficient(name):       if self.reduceSolutionOrder():
1445          if self.COEFFICIENTS[name].isAlteringOperator(): self.__rebuildOperator()           return escript.ReducedSolution(self.getDomain())
         if self.COEFFICIENTS[name].isAlteringRightHandSide(): self.__rebuildRightHandSide()  
1446       else:       else:
1447          raise ValueError,"unknown coefficient %s requested"%name           return escript.Solution(self.getDomain())
   
    # ===== debug ==============================================================  
    def setDebugOn(self):  
        """  
        """  
        self.__debug=not None  
1448    
1449     def setDebugOff(self):     # ==========================================================================
1450         """     #   solver settings:
1451         """     # ==========================================================================
1452         self.__debug=None     def setSolverOptions(self,options=None):
1453           """
1454     def debug(self):         Sets the solver options.
1455    
1456           @param options: the new solver options. If equal C{None}, the solver options are set to the default.
1457           @type options: L{SolverOptions} or C{None}
1458           @note: The symmetry flag of options is overwritten by the symmetry flag of the L{LinearProblem}.
1459           """
1460           if options==None:
1461              self.__solver_options=SolverOptions()
1462           elif isinstance(options, SolverOptions):
1463              self.__solver_options=options
1464           else:
1465              raise ValueError,"options must be a SolverOptions object."
1466           self.__solver_options.setSymmetry(self.__sym)
1467        
1468       def getSolverOptions(self):
1469         """         """
1470         returns true if the PDE is in the debug mode         Returns the solver options
1471      
1472           @rtype: L{SolverOptions}
1473         """         """
1474         return self.__debug         self.__solver_options.setSymmetry(self.__sym)
1475           return self.__solver_options
1476     #===== Lumping ===========================        
    def setLumpingOn(self):  
       """  
       indicates to use matrix lumping  
       """  
       if not self.isUsingLumping():  
          if self.debug() : print "PDE Debug: lumping is set on"  
          self.__rebuildOperator()  
          self.__lumping=True  
   
    def setLumpingOff(self):  
       """  
       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):  
       """  
       set the matrix lumping flag to flag  
       """  
       if flag:  
          self.setLumpingOn()  
       else:  
          self.setLumpingOff()  
   
1477     def isUsingLumping(self):     def isUsingLumping(self):
1478        """        """
1479                Checks if matrix lumping is the current solver method.
       """  
       return self.__lumping  
1480    
1481     #============ method business =========================================================        @return: True if the current solver method is lumping
1482     def setSolverMethod(self,solver=util.DEFAULT_METHOD):        @rtype: C{bool}
1483         """        """
1484         sets a new solver        return self.getSolverOptions().getSolverMethod()==self.getSolverOptions().LUMPING
1485         """     # ==========================================================================
1486         if not solver==self.getSolverMethod():     #    symmetry  flag:
1487             self.__solver_method=solver     # ==========================================================================
            if self.debug() : print "PDE Debug: New solver is %s"%solver  
            self.__checkMatrixType()  
   
    def getSolverMethod(self):  
        """  
        returns the solver method  
        """  
        return self.__solver_method  
   
    #============ tolerance business =========================================================  
    def setTolerance(self,tol=1.e-8):  
        """  
        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):  
        """  
        returns the tolerance set for the solution  
        """  
        return self.__tolerance  
   
    #===== symmetry  flag ==========================  
1488     def isSymmetric(self):     def isSymmetric(self):
1489        """        """
1490        returns true is the operator is considered to be symmetric        Checks if symmetry is indicated.
1491    
1492          @return: True if a symmetric PDE is indicated, False otherwise
1493          @rtype: C{bool}
1494          @note: the method is equivalent to use getSolverOptions().isSymmetric()
1495        """        """
1496        return self.__sym        self.getSolverOptions().isSymmetric()
1497    
1498     def setSymmetryOn(self):     def setSymmetryOn(self):
1499        """        """
1500        sets the symmetry flag to true        Sets the symmetry flag.
1501          @note: The method overwrites the symmetry flag set by the solver options
1502        """        """
1503        if not self.isSymmetric():        self.__sym=True
1504           if self.debug() : print "PDE Debug: Operator is set to be symmetric"        self.getSolverOptions().setSymmetryOn()
          self.__sym=True  
          self.__checkMatrixType()  
1505    
1506     def setSymmetryOff(self):     def setSymmetryOff(self):
1507        """        """
1508        sets the symmetry flag to false        Clears the symmetry flag.
1509          @note: The method overwrites the symmetry flag set by the solver options
1510        """        """
1511        if self.isSymmetric():        self.__sym=False
1512           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):  
      """  
      sets the symmetry flag to flag  
1513    
1514       @param flag:     def setSymmetry(self,flag=False):
1515       """        """
1516       if flag:        Sets the symmetry flag to C{flag}.
         self.setSymmetryOn()  
      else:  
         self.setSymmetryOff()  
1517    
1518     #===== order reduction ==========================        @param flag: If True, the symmetry flag is set otherwise reset.
1519          @type flag: C{bool}
1520          @note: The method overwrites the symmetry flag set by the solver options
1521          """
1522          self.getSolverOptions().setSymmetry(flag)
1523       # ==========================================================================
1524       # function space handling for the equation as well as the solution
1525       # ==========================================================================
1526     def setReducedOrderOn(self):     def setReducedOrderOn(self):
1527       """       """
1528       switches to on reduced order       Switches reduced order on for solution and equation representation.
1529    
1530         @raise RuntimeError: if order reduction is altered after a coefficient has
1531                              been set
1532       """       """
1533       self.setReducedOrderForSolutionOn()       self.setReducedOrderForSolutionOn()
1534       self.setReducedOrderForEquationOn()       self.setReducedOrderForEquationOn()
1535    
1536     def setReducedOrderOff(self):     def setReducedOrderOff(self):
1537       """       """
1538       switches to full order       Switches reduced order off for solution and equation representation
1539    
1540         @raise RuntimeError: if order reduction is altered after a coefficient has
1541                              been set
1542       """       """
1543       self.setReducedOrderForSolutionOff()       self.setReducedOrderForSolutionOff()
1544       self.setReducedOrderForEquationOff()       self.setReducedOrderForEquationOff()
1545    
1546     def setReducedOrderTo(self,flag=False):     def setReducedOrderTo(self,flag=False):
1547       """       """
1548       sets order according to flag       Sets order reduction state for both solution and equation representation
1549         according to flag.
1550    
1551       @param flag:       @param flag: if True, the order reduction is switched on for both solution
1552                      and equation representation, otherwise or if flag is not
1553                      present order reduction is switched off
1554         @type flag: C{bool}
1555         @raise RuntimeError: if order reduction is altered after a coefficient has
1556                              been set
1557       """       """
1558       self.setReducedOrderForSolutionTo(flag)       self.setReducedOrderForSolutionTo(flag)
1559       self.setReducedOrderForEquationTo(flag)       self.setReducedOrderForEquationTo(flag)
                                                                                                                                                             
1560    
1561     #===== order reduction solution ==========================  
1562     def setReducedOrderForSolutionOn(self):     def setReducedOrderForSolutionOn(self):
1563       """       """
1564       switches to reduced order to interpolate solution       Switches reduced order on for solution representation.
1565    
1566         @raise RuntimeError: if order reduction is altered after a coefficient has
1567                              been set
1568       """       """
1569       new_fs=escript.ReducedSolution(self.getDomain())       if not self.__reduce_solution_order:
1570       if self.getFunctionSpaceForSolution()!=new_fs:           if self.__altered_coefficients:
1571           if self.debug() : print "PDE Debug: Reduced order is used to interpolate solution."                raise RuntimeError,"order cannot be altered after coefficients have been defined."
1572           self.__column_function_space=new_fs           self.trace("Reduced order is used for solution representation.")
1573           self.__rebuildSystem(deep=True)           self.__reduce_solution_order=True
1574             self.initializeSystem()
1575    
1576     def setReducedOrderForSolutionOff(self):     def setReducedOrderForSolutionOff(self):
1577       """       """
1578       switches to full order to interpolate solution       Switches reduced order off for solution representation
1579    
1580         @raise RuntimeError: if order reduction is altered after a coefficient has
1581                              been set.
1582       """       """
1583       new_fs=escript.Solution(self.getDomain())       if self.__reduce_solution_order:
1584       if self.getFunctionSpaceForSolution()!=new_fs:           if self.__altered_coefficients:
1585           if self.debug() : print "PDE Debug: Full order is used to interpolate solution."                raise RuntimeError,"order cannot be altered after coefficients have been defined."
1586           self.__column_function_space=new_fs           self.trace("Full order is used to interpolate solution.")
1587           self.__rebuildSystem(deep=True)           self.__reduce_solution_order=False
1588             self.initializeSystem()
1589    
1590     def setReducedOrderForSolutionTo(self,flag=False):     def setReducedOrderForSolutionTo(self,flag=False):
1591       """       """
1592       sets order for test functions according to flag       Sets order reduction state for solution representation according to flag.
1593    
1594       @param flag:       @param flag: if flag is True, the order reduction is switched on for
1595                      solution representation, otherwise or if flag is not present
1596                      order reduction is switched off
1597         @type flag: C{bool}
1598         @raise RuntimeError: if order reduction is altered after a coefficient has
1599                              been set
1600       """       """
1601       if flag:       if flag:
1602          self.setReducedOrderForSolutionOn()          self.setReducedOrderForSolutionOn()
1603       else:       else:
1604          self.setReducedOrderForSolutionOff()          self.setReducedOrderForSolutionOff()
1605                                                                                                                                                              
    #===== order reduction equation ==========================  
1606     def setReducedOrderForEquationOn(self):     def setReducedOrderForEquationOn(self):
1607       """       """
1608       switches to reduced order for test functions       Switches reduced order on for equation representation.
1609    
1610         @raise RuntimeError: if order reduction is altered after a coefficient has
1611                              been set
1612       """       """
1613       new_fs=escript.ReducedSolution(self.getDomain())       if not self.__reduce_equation_order:
1614       if self.getFunctionSpaceForEquation()!=new_fs:           if self.__altered_coefficients:
1615           if self.debug() : print "PDE Debug: Reduced order is used for test functions."                raise RuntimeError,"order cannot be altered after coefficients have been defined."
1616           self.__row_function_space=new_fs           self.trace("Reduced order is used for test functions.")
1617           self.__rebuildSystem(deep=True)           self.__reduce_equation_order=True
1618             self.initializeSystem()
1619    
1620     def setReducedOrderForEquationOff(self):     def setReducedOrderForEquationOff(self):
1621       """       """
1622       switches to full order for test functions       Switches reduced order off for equation representation.
1623    
1624         @raise RuntimeError: if order reduction is altered after a coefficient has
1625                              been set
1626       """       """
1627       new_fs=escript.Solution(self.getDomain())       if self.__reduce_equation_order:
1628       if self.getFunctionSpaceForEquation()!=new_fs:           if self.__altered_coefficients:
1629           if self.debug() : print "PDE Debug: Full order is used for test functions."                raise RuntimeError,"order cannot be altered after coefficients have been defined."
1630           self.__row_function_space=new_fs           self.trace("Full order is used for test functions.")
1631           self.__rebuildSystem(deep=True)           self.__reduce_equation_order=False
1632             self.initializeSystem()
1633    
1634     def setReducedOrderForEquationTo(self,flag=False):     def setReducedOrderForEquationTo(self,flag=False):
1635       """       """
1636       sets order for test functions according to flag       Sets order reduction state for equation representation according to flag.
1637    
1638       @param flag:       @param flag: if flag is True, the order reduction is switched on for
1639                      equation representation, otherwise or if flag is not present
1640                      order reduction is switched off
1641         @type flag: C{bool}
1642         @raise RuntimeError: if order reduction is altered after a coefficient has
1643                              been set
1644       """       """
1645       if flag:       if flag:
1646          self.setReducedOrderForEquationOn()          self.setReducedOrderForEquationOn()
1647       else:       else:
1648          self.setReducedOrderForEquationOff()          self.setReducedOrderForEquationOff()
1649                                                                                                                                                                 def getOperatorType(self):
1650     # ==== initialization =====================================================================        """
1651     def __getNewOperator(self):        Returns the current system type.
1652          """
1653          return self.__operator_type
1654    
1655       def checkSymmetricTensor(self,name,verbose=True):
1656          """
1657          Tests a coefficient for symmetry.
1658    
1659          @param name: name of the coefficient
1660          @type name: C{str}
1661          @param verbose: if set to True or not present a report on coefficients
1662                          which break the symmetry is printed.
1663          @type verbose: C{bool}
1664          @return: True if coefficient C{name} is symmetric
1665          @rtype: C{bool}
1666          """
1667          SMALL_TOLERANCE=util.EPSILON*10.
1668          A=self.getCoefficient(name)
1669          verbose=verbose or self.__debug
1670          out=True
1671          if not A.isEmpty():
1672             tol=util.Lsup(A)*SMALL_TOLERANCE
1673             s=A.getShape()
1674             if A.getRank() == 4:
1675                if s[0]==s[2] and s[1] == s[3]:
1676                   for i in range(s[0]):
1677                      for j in range(s[1]):
1678                         for k in range(s[2]):
1679                            for l in range(s[3]):
1680                                if util.Lsup(A[i,j,k,l]-A[k,l,i,j])>tol:
1681                                   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)
1682                                   out=False
1683                else:
1684                     if verbose: print "non-symmetric problem because of inappropriate shape %s of coefficient %s."%(s,name)
1685                     out=False
1686             elif A.getRank() == 2:
1687                if s[0]==s[1]:
1688                   for j in range(s[0]):
1689                      for l in range(s[1]):
1690                         if util.Lsup(A[j,l]-A[l,j])>tol:
1691                            if verbose: print "non-symmetric problem because %s[%d,%d]!=%s[%d,%d]"%(name,j,l,name,l,j)
1692                            out=False
1693                else:
1694                     if verbose: print "non-symmetric problem because of inappropriate shape %s of coefficient %s."%(s,name)
1695                     out=False
1696             elif A.getRank() == 0:
1697                pass
1698             else:
1699                 raise ValueError,"Cannot check rank %s of %s."%(A.getRank(),name)
1700          return out
1701    
1702       def checkReciprocalSymmetry(self,name0,name1,verbose=True):
1703          """
1704          Tests two coefficients for reciprocal symmetry.
1705    
1706          @param name0: name of the first coefficient
1707          @type name0: C{str}
1708          @param name1: name of the second coefficient
1709          @type name1: C{str}
1710          @param verbose: if set to True or not present a report on coefficients
1711                          which break the symmetry is printed
1712          @type verbose: C{bool}
1713          @return: True if coefficients C{name0} and C{name1} are reciprocally
1714                   symmetric.
1715          @rtype: C{bool}
1716          """
1717          SMALL_TOLERANCE=util.EPSILON*10.
1718          B=self.getCoefficient(name0)
1719          C=self.getCoefficient(name1)
1720          verbose=verbose or self.__debug
1721          out=True
1722          if B.isEmpty() and not C.isEmpty():
1723             if verbose: print "non-symmetric problem because %s is not present but %s is"%(name0,name1)
1724             out=False
1725          elif not B.isEmpty() and C.isEmpty():
1726             if verbose: print "non-symmetric problem because %s is not present but %s is"%(name0,name1)
1727             out=False
1728          elif not B.isEmpty() and not C.isEmpty():
1729             sB=B.getShape()
1730             sC=C.getShape()
1731             tol=(util.Lsup(B)+util.Lsup(C))*SMALL_TOLERANCE/2.
1732             if len(sB) != len(sC):
1733                 if verbose: print "non-symmetric problem because ranks of %s (=%s) and %s (=%s) are different."%(name0,len(sB),name1,len(sC))
1734                 out=False
1735             else:
1736                 if len(sB)==0:
1737                   if util.Lsup(B-C)>tol:
1738                      if verbose: print "non-symmetric problem because %s!=%s"%(name0,name1)
1739                      out=False
1740                 elif len(sB)==1:
1741                   if sB[0]==sC[0]:
1742                      for j in range(sB[0]):
1743                         if util.Lsup(B[j]-C[j])>tol:
1744                            if verbose: print "non-symmetric PDE because %s[%d]!=%s[%d]"%(name0,j,name1,j)
1745                            out=False
1746                   else:
1747                     if verbose: print "non-symmetric problem because of inappropriate shapes %s and %s of coefficients %s and %s, respectively."%(sB,sC,name0,name1)
1748                 elif len(sB)==3:
1749                   if sB[0]==sC[1] and sB[1]==sC[2] and sB[2]==sC[0]:
1750                       for i in range(sB[0]):
1751                          for j in range(sB[1]):
1752                             for k in range(sB[2]):
1753                                if util.Lsup(B[i,j,k]-C[k,i,j])>tol:
1754                                     if verbose: print "non-symmetric problem because %s[%d,%d,%d]!=%s[%d,%d,%d]"%(name0,i,j,k,name1,k,i,j)
1755                                     out=False
1756                   else:
1757                     if verbose: print "non-symmetric problem because of inappropriate shapes %s and %s of coefficients %s and %s, respectively."%(sB,sC,name0,name1)
1758                 else:
1759                     raise ValueError,"Cannot check rank %s of %s and %s."%(len(sB),name0,name1)
1760          return out
1761    
1762       def getCoefficient(self,name):
1763         """
1764         Returns the value of the coefficient C{name}.
1765    
1766         @param name: name of the coefficient requested
1767         @type name: C{string}
1768         @return: the value of the coefficient
1769         @rtype: L{Data<escript.Data>}
1770         @raise IllegalCoefficient: if C{name} is not a coefficient of the PDE
1771         """
1772         if self.hasCoefficient(name):
1773             return self.__COEFFICIENTS[name].getValue()
1774         else:
1775            raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1776    
1777       def hasCoefficient(self,name):
1778         """
1779         Returns True if C{name} is the name of a coefficient.
1780    
1781         @param name: name of the coefficient enquired
1782         @type name: C{string}
1783         @return: True if C{name} is the name of a coefficient of the general PDE,
1784                  False otherwise
1785         @rtype: C{bool}
1786         """
1787         return self.__COEFFICIENTS.has_key(name)
1788    
1789       def createCoefficient(self, name):
1790         """
1791         Creates a L{Data<escript.Data>} object corresponding to coefficient
1792         C{name}.
1793    
1794         @return: the coefficient C{name} initialized to 0
1795         @rtype: L{Data<escript.Data>}
1796         @raise IllegalCoefficient: if C{name} is not a coefficient of the PDE
1797         """
1798         if self.hasCoefficient(name):
1799            return escript.Data(0.,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))
1800         else:
1801            raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1802    
1803       def getFunctionSpaceForCoefficient(self,name):
1804         """
1805         Returns the L{FunctionSpace<escript.FunctionSpace>} to be used for
1806         coefficient C{name}.
1807    
1808         @param name: name of the coefficient enquired
1809         @type name: C{string}
1810         @return: the function space to be used for coefficient C{name}
1811         @rtype: L{FunctionSpace<escript.FunctionSpace>}
1812         @raise IllegalCoefficient: if C{name} is not a coefficient of the PDE
1813         """
1814         if self.hasCoefficient(name):
1815            return self.__COEFFICIENTS[name].getFunctionSpace(self.getDomain())
1816         else:
1817            raise ValueError,"unknown coefficient %s requested"%name
1818    
1819       def getShapeOfCoefficient(self,name):
1820         """
1821         Returns the shape of the coefficient C{name}.
1822    
1823         @param name: name of the coefficient enquired
1824         @type name: C{string}
1825         @return: the shape of the coefficient C{name}
1826         @rtype: C{tuple} of C{int}
1827         @raise IllegalCoefficient: if C{name} is not a coefficient of the PDE
1828         """
1829         if self.hasCoefficient(name):
1830            return self.__COEFFICIENTS[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())
1831         else:
1832            raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1833    
1834       def resetAllCoefficients(self):
1835         """
1836         Resets all coefficients to their default values.
1837         """
1838         for i in self.__COEFFICIENTS.iterkeys():
1839             self.__COEFFICIENTS[i].resetValue()
1840    
1841       def alteredCoefficient(self,name):
1842         """
1843         Announces that coefficient C{name} has been changed.
1844    
1845         @param name: name of the coefficient affected
1846         @type name: C{string}
1847         @raise IllegalCoefficient: if C{name} is not a coefficient of the PDE
1848         @note: if C{name} is q or r, the method will not trigger a rebuild of the
1849                system as constraints are applied to the solved system.
1850         """
1851         if self.hasCoefficient(name):
1852            self.trace("Coefficient %s has been altered."%name)
1853            if not ((name=="q" or name=="r") and self.isUsingLumping()):
1854               if self.__COEFFICIENTS[name].isAlteringOperator(): self.invalidateOperator()
1855               if self.__COEFFICIENTS[name].isAlteringRightHandSide(): self.invalidateRightHandSide()
1856         else:
1857            raise IllegalCoefficient,"illegal coefficient %s requested for general PDE."%name
1858    
1859       def validSolution(self):
1860         """         """
1861           Marks the solution as valid.
1862         """         """
1863         return self.getDomain().newOperator( \         self.__is_solution_valid=True
                            self.getNumEquations(), \  
                            self.getFunctionSpaceForEquation(), \  
                            self.getNumSolutions(), \  
                            self.getFunctionSpaceForSolution(), \  
                            self.__matrix_type)  
1864    
1865     def __makeFreshRightHandSide(self):     def invalidateSolution(self):
1866         """         """
1867           Indicates the PDE has to be resolved if the solution is requested.
1868         """         """
1869         if self.debug() : print "PDE Debug: New right hand side allocated"         self.trace("System will be resolved.")
1870         if self.getNumEquations()>1:         self.__is_solution_valid=False
            self.__righthandside=escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)  
        else:  
            self.__righthandside=escript.Data(0.,(),self.getFunctionSpaceForEquation(),True)  
        return self.__righthandside  
1871    
1872     def __getNewSolution(self):     def isSolutionValid(self):
1873         """         """
1874           Returns True if the solution is still valid.
1875         """         """
1876         if self.debug() : print "PDE Debug: New right hand side allocated"         if self.__solution_rtol>self.getSolverOptions().getTolerance() or \
1877         if self.getNumSolutions()>1:            self.__solution_atol>self.getSolverOptions().getAbsoluteTolerance():
1878             return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)           self.invalidateSolution()  
1879         else:         return self.__is_solution_valid
            return escript.Data(0.,(),self.getFunctionSpaceForSolution(),True)  
1880    
1881     def __makeFreshOperator(self):     def validOperator(self):
1882         """         """
1883           Marks the operator as valid.
1884         """         """
1885         if self.__operator.isEmpty():         self.__is_operator_valid=True
            self.__operator=self.__getNewOperator()  
            if self.debug() : print "PDE Debug: New operator allocated"  
        else:  
            self.__operator.setValue(0.)  
            self.__operator.resetSolver()  
            if self.debug() : print "PDE Debug: Operator reset to zero"  
        return self.__operator  
1886    
1887     #============ some serivice functions  =====================================================     def invalidateOperator(self):
1888     def getDomain(self):         """
1889       """         Indicates the operator has to be rebuilt next time it is used.
1890       returns the domain of the PDE         """
1891       """         self.trace("Operator will be rebuilt.")
1892       return self.__domain         self.invalidateSolution()
1893           self.__is_operator_valid=False
1894    
1895     def getDim(self):     def isOperatorValid(self):
1896       """         """
1897       returns the spatial dimension of the PDE         Returns True if the operator is still valid.
1898       """         """
1899       return self.getDomain().getDim()         if self.getRequiredOperatorType()==self.getOperatorType(): self.invalidateOperator()
1900           return self.__is_operator_valid
1901    
1902     def getNumEquations(self):     def validRightHandSide(self):
1903           """
1904           Marks the right hand side as valid.
1905           """
1906           self.__is_RHS_valid=True
1907    
1908       def invalidateRightHandSide(self):
1909           """
1910           Indicates the right hand side has to be rebuilt next time it is used.
1911           """
1912           if self.isRightHandSideValid(): self.trace("Right hand side has to be rebuilt.")
1913           self.invalidateSolution()
1914           self.__is_RHS_valid=False
1915    
1916       def isRightHandSideValid(self):
1917           """
1918           Returns True if the operator is still valid.
1919           """
1920           return self.__is_RHS_valid
1921    
1922       def invalidateSystem(self):
1923           """
1924           Announces that everything has to be rebuilt.
1925           """
1926           self.invalidateSolution()
1927           self.invalidateOperator()
1928           self.invalidateRightHandSide()
1929    
1930       def isSystemValid(self):
1931           """
1932           Returns True if the system (including solution) is still vaild.
1933           """
1934           return self.isSolutionValid() and self.isOperatorValid() and self.isRightHandSideValid()
1935    
1936       def initializeSystem(self):
1937           """
1938           Resets the system clearing the operator, right hand side and solution.
1939           """
1940           self.trace("New System has been created.")
1941           self.__operator_type=None
1942           self.__operator=escript.Operator()
1943           self.__righthandside=escript.Data()
1944           self.__solution=escript.Data()
1945           self.invalidateSystem()
1946    
1947       def getOperator(self):
1948       """       """
1949       returns the number of equations       Returns the operator of the linear problem.
1950    
1951         @return: the operator of the problem
1952       """       """
1953       if self.__numEquations>0:       return self.getSystem()[0]
          return self.__numEquations  
      else:  
          raise ValueError,"Number of equations is undefined. Please specify argument numEquations."  
1954    
1955     def getNumSolutions(self):     def getRightHandSide(self):
1956       """       """
1957       returns the number of unknowns       Returns the right hand side of the linear problem.
1958    
1959         @return: the right hand side of the problem
1960         @rtype: L{Data<escript.Data>}
1961       """       """
1962       if self.__numSolutions>0:       return self.getSystem()[1]
         return self.__numSolutions  
      else:  
         raise ValueError,"Number of solution is undefined. Please specify argument numSolutions."  
1963    
1964       def createRightHandSide(self):
1965           """
1966           Returns an instance of a new right hand side.
1967           """
1968           self.trace("New right hand side is allocated.")
1969           if self.getNumEquations()>1:
1970               return escript.Data(0.,(self.getNumEquations(),),self.getFunctionSpaceForEquation(),True)
1971           else:
1972               return escript.Data(0.,(),self.getFunctionSpaceForEquation(),True)
1973    
1974     def checkSymmetry(self,verbose=True):     def createSolution(self):
1975        """         """
1976        returns if the Operator is symmetric. This is a very expensive operation!!! The symmetry flag is not altered.         Returns an instance of a new solution.
1977        """         """
1978        verbose=verbose or self.debug()         self.trace("New solution is allocated.")
1979        out=True         if self.getNumSolutions()>1:
1980        if self.getNumSolutions()!=self.getNumEquations():             return escript.Data(0.,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)
1981           if verbose: print "non-symmetric PDE because of different number of equations and solutions"         else:
1982           out=False             return escript.Data(0.,(),self.getFunctionSpaceForSolution(),True)
       else:  
          A=self.getCoefficientOfPDE("A")  
          if not A.isEmpty():  
             tol=util.Lsup(A)*self.TOL  
             if self.getNumSolutions()>1:  
                for i in range(self.getNumEquations()):  
                   for j in range(self.getDim()):  
                      for k in range(self.getNumSolutions()):  
                         for l in range(self.getDim()):  
                             if util.Lsup(A[i,j,k,l]-A[k,l,i,j])>tol:  
                                if verbose: print "non-symmetric PDE because A[%d,%d,%d,%d]!=A[%d,%d,%d,%d]"%(i,j,k,l,k,l,i,j)  
                                out=False  
             else:  
                for j in range(self.getDim()):  
                   for l in range(self.getDim()):  
                      if util.Lsup(A[j,l]-A[l,j])>tol:  
                         if verbose: print "non-symmetric PDE because A[%d,%d]!=A[%d,%d]"%(j,l,l,j)  
                         out=False  
          B=self.getCoefficientOfPDE("B")  
          C=self.getCoefficientOfPDE("C")  
          if B.isEmpty() and not C.isEmpty():  
             if verbose: print "non-symmetric PDE because B is not present but C is"  
             out=False  
          elif not B.isEmpty() and C.isEmpty():  
             if verbose: print "non-symmetric PDE because C is not present but B is"  
             out=False  
          elif not B.isEmpty() and not C.isEmpty():  
             tol=(util.Lsup(B)+util.Lsup(C))*self.TOL/2.  
             if self.getNumSolutions()>1:  
                for i in range(self.getNumEquations()):  
                    for j in range(self.getDim()):  
                       for k in range(self.getNumSolutions()):  
                          if util.Lsup(B[i,j,k]-C[k,i,j])>tol:  
                               if verbose: print "non-symmetric PDE because B[%d,%d,%d]!=C[%d,%d,%d]"%(i,j,k,k,i,j)  
                               out=False  
             else:  
                for j in range(self.getDim()):  
                   if util.Lsup(B[j]-C[j])>tol:  
                      if verbose: print "non-symmetric PDE because B[%d]!=C[%d]"%(j,j)  
                      out=False  
          if self.getNumSolutions()>1:  
            D=self.getCoefficientOfPDE("D")  
            if not D.isEmpty():  
              tol=util.Lsup(D)*self.TOL  
              for i in range(self.getNumEquations()):  
                 for k in range(self.getNumSolutions()):  
                   if util.Lsup(D[i,k]-D[k,i])>tol:  
                       if verbose: print "non-symmetric PDE because D[%d,%d]!=D[%d,%d]"%(i,k,k,i)  
                       out=False  
               
       return out  
1983    
1984     def getFlux(self,u):     def resetSolution(self):
1985           """
1986           Sets the solution to zero.
1987         """         """
1988         returns the flux J_ij for a given u         if self.__solution.isEmpty():
1989               self.__solution=self.createSolution()
1990           else:
1991               self.__solution.setToZero()
1992               self.trace("Solution is reset to zero.")
1993    
1994         \f[     def setSolution(self,u):
1995         J_ij=A_{ijkl}u_{k,l}+B_{ijk}u_k-X_{ij}         """
1996         \f]         Sets the solution assuming that makes the system valid with the tolrance
1997           defined by the solver options
1998           """
1999           self.__solution_rtol=self.getSolverOptions().getTolerance()
2000           self.__solution_atol=self.getSolverOptions().getAbsoluteTolerance()
2001           self.__solution=u
2002           self.validSolution()
2003    
2004         @param u: argument of the operator     def getCurrentSolution(self):
2005         """         """
2006         raise SystemError,"getFlux is not implemented yet"         Returns the solution in its current state.
2007         return None         """
2008           if self.__solution.isEmpty(): self.__solution=self.createSolution()
2009           return self.__solution
2010    
2011     def applyOperator(self,u):     def resetRightHandSide(self):
2012         """         """
2013         applies the operator of the PDE to a given solution u in weak from         Sets the right hand side to zero.
2014           """
2015           if self.__righthandside.isEmpty():
2016               self.__righthandside=self.createRightHandSide()
2017           else:
2018               self.__righthandside.setToZero()
2019               self.trace("Right hand side is reset to zero.")
2020    
2021         @param u: argument of the operator     def getCurrentRightHandSide(self):
2022         """         """
2023         return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())         Returns the right hand side in its current state.
                                                                                                                                                             
    def getResidual(self,u):  
2024         """         """
2025         return the residual of u in the weak from         if self.__righthandside.isEmpty(): self.__righthandside=self.createRightHandSide()
2026           return self.__righthandside
2027    
2028         @param u:     def resetOperator(self):
2029         """         """
2030         return self.applyOperator(u)-self.getRightHandSide()         Makes sure that the operator is instantiated and returns it initialized
2031           with zeros.
2032           """
2033           if self.getOperatorType() == None:
2034               if self.isUsingLumping():
2035                   self.__operator=self.createSolution()
2036               else:
2037                   self.__operator=self.createOperator()
2038           self.__operator_type=self.getRequiredOperatorType()
2039           else:
2040               if self.isUsingLumping():
2041                   self.__operator.setToZero()
2042               else:
2043                   self.__operator.resetValues()
2044               self.trace("Operator reset to zero")
2045    
2046     def __setValue(self,**coefficients):     def getCurrentOperator(self):
2047           """
2048           Returns the operator in its current state.
2049           """
2050           return self.__operator
2051    
2052       def setValue(self,**coefficients):
2053        """        """
2054        sets new values to coefficient        Sets new values to coefficients.
2055    
2056        @param coefficients:        @raise IllegalCoefficient: if an unknown coefficient keyword is used
2057        """        """
2058        # check if the coefficients are  legal:        # check if the coefficients are  legal:
2059        for i in coefficients.iterkeys():        for i in coefficients.iterkeys():
2060           if not self.hasCoefficient(i):           if not self.hasCoefficient(i):
2061              raise ValueError,"Attempt to set unknown coefficient %s"%i              raise IllegalCoefficient,"Attempt to set unknown coefficient %s"%i
2062        # if the number of unknowns or equations is still unknown we try to estimate them:        # if the number of unknowns or equations is still unknown we try to estimate them:
2063        if self.__numEquations<1 or self.__numSolutions<1:        if self.__numEquations==None or self.__numSolutions==None:
2064           for i,d in coefficients.iteritems():           for i,d in coefficients.iteritems():
2065              if hasattr(d,"shape"):              if hasattr(d,"shape"):
2066                  s=d.shape                  s=d.shape
2067              elif hasattr(d,"getShape"):              elif hasattr(d,"getShape"):
2068                  s=d.getShape()                  s=d.getShape()
2069              else:              else:
2070                  s=numarray.array(d).shape                  s=numpy.array(d).shape
2071              if s!=None:              if s!=None:
2072                  # get number of equations and number of unknowns:                  # get number of equations and number of unknowns:
2073                  res=self.COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(s,self.getDim())                  res=self.__COEFFICIENTS[i].estimateNumEquationsAndNumSolutions(self.getDomain(),s)
2074                  if res==None:                  if res==None:
2075                      raise ValueError,"Illegal shape %s of coefficient %s"%(s,i)                      raise IllegalCoefficientValue,"Illegal shape %s of coefficient %s"%(s,i)
2076                  else:                  else:
2077                      if self.__numEquations<1: self.__numEquations=res[0]                      if self.__numEquations==None: self.__numEquations=res[0]
2078                      if self.__numSolutions<1: self.__numSolutions=res[1]                      if self.__numSolutions==None: self.__numSolutions=res[1]
2079        if self.__numEquations<1: raise ValueError,"unidententified number of equations"        if self.__numEquations==None: raise UndefinedPDEError,"unidentified number of equations"
2080        if self.__numSolutions<1: raise ValueError,"unidententified number of solutions"        if self.__numSolutions==None: raise UndefinedPDEError,"unidentified number of solutions"
2081        # now we check the shape of the coefficient if numEquations and numSolutions are set:        # now we check the shape of the coefficient if numEquations and numSolutions are set:
2082        for i,d in coefficients.iteritems():        for i,d in coefficients.iteritems():
2083          if d==None:          try:
2084               d2=escript.Data()             self.__COEFFICIENTS[i].setValue(self.getDomain(),
2085          elif isinstance(d,escript.Data):                       self.getNumEquations(),self.getNumSolutions(),
2086               if d.isEmpty():                       self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
2087                  d2=d             self.alteredCoefficient(i)
2088               else:          except IllegalCoefficientFunctionSpace,m:
2089                  d2=escript.Data(d,self.getFunctionSpaceForCoefficient(i))              # if the function space is wrong then we try the reduced version:
2090          else:              i_red=i+"_reduced"
2091                d2=escript.Data(d,self.getFunctionSpaceForCoefficient(i))              if (not i_red in coefficients.keys()) and i_red in self.__COEFFICIENTS.keys():
2092          if not d2.isEmpty():                  try:
2093             if not self.getShapeOfCoefficient(i)==d2.getShape():                      self.__COEFFICIENTS[i_red].setValue(self.getDomain(),
2094                 raise ValueError,"Expected shape for coefficient %s is %s but actual shape is %s."%(i,self.getShapeOfCoefficient(i),d2.getShape())                                                        self.getNumEquations(),self.getNumSolutions(),
2095          # overwrite new values:                                                        self.reduceEquationOrder(),self.reduceSolutionOrder(),d)
2096          if self.debug(): print "PDE Debug: Coefficient %s has been altered."%i                      self.alteredCoefficient(i_red)
2097          self.COEFFICIENTS[i].setValue(d2)                  except IllegalCoefficientValue,m:
2098          self.alteredCoefficient(i)                      raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
2099                          except IllegalCoefficientFunctionSpace,m:
2100        # reset the HomogeneousConstraintFlag:                      raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m))
2101        self.__setHomogeneousConstraintFlag()              else:
2102        if len(coefficients)>0 and not self.isUsingLumping() and not self.__homogeneous_constraint: self.__rebuildSystem()                  raise IllegalCoefficientFunctionSpace("Coefficient %s:%s"%(i,m))
2103            except IllegalCoefficientValue,m:
2104     def __setHomogeneousConstraintFlag(self):             raise IllegalCoefficientValue("Coefficient %s:%s"%(i,m))
2105        """        self.__altered_coefficients=True
       checks if the constraints are homogeneous and sets self.__homogeneous_constraint accordingly.  
       """  
       self.__homogeneous_constraint=True  
       q=self.getCoefficientOfPDE("q")  
       r=self.getCoefficientOfPDE("r")  
       if not q.isEmpty() and not r.isEmpty():  
          if (q*r).Lsup()>=1.e-13*r.Lsup(): self.__homogeneous_constraint=False  
       if self.debug():  
            if self.__homogeneous_constraint:  
                print "PDE Debug: Constraints are homogeneous."  
            else:  
                print "PDE Debug: Constraints are inhomogeneous."  
   
2106    
2107     # ==== rebuild switches =====================================================================     # ==========================================================================
2108     def __rebuildSolution(self,deep=False):     # methods that are typically overwritten when implementing a particular
2109         """     # linear problem
2110         indicates the PDE has to be reolved if the solution is requested     # ==========================================================================
2111         """     def getRequiredOperatorType(self):
2112         if self.__solution_isValid and self.debug() : print "PDE Debug: PDE has to be resolved."        """
2113         self.__solution_isValid=False        Returns the system type which needs to be used by the current set up.
        if deep: self.__solution=escript.Data()  
2114    
2115          @note: Typically this method is overwritten when implementing a
2116                 particular linear problem.
2117          """
2118          return None
2119    
2120     def __rebuildOperator(self,deep=False):     def createOperator(self):
2121         """         """
2122         indicates the operator has to be rebuilt next time it is used         Returns an instance of a new operator.
2123    
2124           @note: This method is overwritten when implementing a particular
2125                  linear problem.
2126         """         """
2127         if self.__operator_isValid and self.debug() : print "PDE Debug: Operator has to be rebuilt."         return escript.Operator()
2128         self.__rebuildSolution(deep)  
2129         self.__operator_isValid=False     def checkSymmetry(self,verbose=True):
2130         if deep: self.__operator=escript.Operator()        """
2131          Tests the PDE for symmetry.
2132    
2133          @param verbose: if set to True or not present a report on coefficients
2134                          which break the symmetry is printed
2135          @type verbose: C{bool}
2136          @return: True if the problem is symmetric
2137          @rtype: C{bool}
2138          @note: Typically this method is overwritten when implementing a
2139                 particular linear problem.
2140          """
2141          out=True
2142          return out
2143    
2144     def __rebuildRightHandSide(self,deep=False):     def getSolution(self,**options):
2145         """         """
2146         indicates the right hand side has to be rebuild next time it is used         Returns the solution of the problem.
2147    
2148           @return: the solution
2149           @rtype: L{Data<escript.Data>}
2150    
2151           @note: This method is overwritten when implementing a particular
2152                  linear problem.
2153         """         """
2154         if self.__righthandside_isValid and self.debug() : print "PDE Debug: Right hand side has to be rebuilt."         return self.getCurrentSolution()
        self.__rebuildSolution(deep)  
        self.__righthandside_isValid=False  
        if deep: self.__righthandside=escript.Data()  
2155    
2156     def __rebuildSystem(self,deep=False):     def getSystem(self):
2157         """         """
2158         annonced that all coefficient name has been changed         Returns the operator and right hand side of the PDE.
2159    
2160           @return: the discrete version of the PDE
2161           @rtype: C{tuple} of L{Operator,<escript.Operator>} and L{Data<escript.Data>}.
2162    
2163           @note: This method is overwritten when implementing a particular
2164                  linear problem.
2165         """         """
2166         self.__rebuildSolution(deep)         return (self.getCurrentOperator(), self.getCurrentRightHandSide())
2167         self.__rebuildOperator(deep)  
2168         self.__rebuildRightHandSide(deep)  class LinearPDE(LinearProblem):
2169         """
2170     def __checkMatrixType(self):     This class is used to define a general linear, steady, second order PDE
2171       for an unknown function M{u} on a given domain defined through a
2172       L{Domain<escript.Domain>} object.
2173    
2174       For a single PDE having a solution with a single component the linear PDE
2175       is defined in the following form:
2176    
2177       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)}
2178    
2179       where M{grad(F)} denotes the spatial derivative of M{F}. Einstein's
2180       summation convention, ie. summation over indexes appearing twice in a term
2181       of a sum performed, is used.
2182       The coefficients M{A}, M{B}, M{C}, M{D}, M{X} and M{Y} have to be specified
2183       through L{Data<escript.Data>} objects in L{Function<escript.Function>} and
2184       the coefficients M{A_reduced}, M{B_reduced}, M{C_reduced}, M{D_reduced},
2185       M{X_reduced} and M{Y_reduced} have to be specified through
2186       L{Data<escript.Data>} objects in L{ReducedFunction<escript.ReducedFunction>}.
2187       It is also allowed to use objects that can be converted into such
2188       L{Data<escript.Data>} objects. M{A} and M{A_reduced} are rank two, M{B},
2189       M{C}, M{X}, M{B_reduced}, M{C_reduced} and M{X_reduced} are rank one and
2190       M{D}, M{D_reduced}, M{Y} and M{Y_reduced} are scalar.
2191    
2192       The following natural boundary conditions are considered:
2193    
2194       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}
2195    
2196       where M{n} is the outer normal field. Notice that the coefficients M{A},
2197       M{A_reduced}, M{B}, M{B_reduced}, M{X} and M{X_reduced} are defined in the
2198       PDE. The coefficients M{d} and M{y} are each a scalar in
2199       L{FunctionOnBoundary<escript.FunctionOnBoundary>} and the coefficients
2200       M{d_reduced} and M{y_reduced} are each a scalar in
2201       L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
2202    
2203       Constraints for the solution prescribe the value of the solution at certain
2204       locations in the domain. They have the form
2205    
2206       M{u=r} where M{q>0}
2207    
2208       M{r} and M{q} are each scalar where M{q} is the characteristic function
2209       defining where the constraint is applied. The constraints override any
2210       other condition set by the PDE or the boundary condition.
2211    
2212       The PDE is symmetrical if
2213    
2214       M{A[i,j]=A[j,i]}  and M{B[j]=C[j]} and M{A_reduced[i,j]=A_reduced[j,i]}
2215       and M{B_reduced[j]=C_reduced[j]}
2216    
2217       For a system of PDEs and a solution with several components the PDE has the
2218       form
2219    
2220       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] }
2221    
2222       M{A} and M{A_reduced} are of rank four, M{B}, M{B_reduced}, M{C} and
2223       M{C_reduced} are each of rank three, M{D}, M{D_reduced}, M{X_reduced} and
2224       M{X} are each of rank two and M{Y} and M{Y_reduced} are of rank one.
2225       The natural boundary conditions take the form:
2226    
2227       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]}
2228    
2229       The coefficient M{d} is of rank two and M{y} is of rank one both in
2230       L{FunctionOnBoundary<escript.FunctionOnBoundary>}. The coefficients
2231       M{d_reduced} is of rank two and M{y_reduced} is of rank one both in
2232       L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.
2233    
2234       Constraints take the form
2235    
2236       M{u[i]=r[i]}  where  M{q[i]>0}
2237    
2238       M{r} and M{q} are each rank one. Notice that at some locations not
2239       necessarily all components must have a constraint.
2240    
2241       The system of PDEs is symmetrical if
2242    
2243          - M{A[i,j,k,l]=A[k,l,i,j]}
2244          - M{A_reduced[i,j,k,l]=A_reduced[k,l,i,j]}
2245          - M{B[i,j,k]=C[k,i,j]}
2246          - M{B_reduced[i,j,k]=C_reduced[k,i,j]}
2247          - M{D[i,k]=D[i,k]}
2248          - M{D_reduced[i,k]=D_reduced[i,k]}
2249          - M{d[i,k]=d[k,i]}
2250          - M{d_reduced[i,k]=d_reduced[k,i]}
2251    
2252       L{LinearPDE} also supports solution discontinuities over a contact region
2253       in the domain. To specify the conditions across the discontinuity we are
2254       using the generalised flux M{J} which, in the case of a system of PDEs
2255       and several components of the solution, is defined as
2256    
2257       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]}
2258    
2259       For the case of single solution component and single PDE M{J} is defined as
2260    
2261       M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[j]+(B[i]+B_reduced[i])*u-X[i]-X_reduced[i]}
2262    
2263       In the context of discontinuities M{n} denotes the normal on the
2264       discontinuity pointing from side 0 towards side 1 calculated from
2265       L{getNormal<escript.FunctionSpace.getNormal>} of L{FunctionOnContactZero<escript.FunctionOnContactZero>}.
2266       For a system of PDEs the contact condition takes the form
2267    
2268       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]}
2269    
2270       where M{J0} and M{J1} are the fluxes on side 0 and side 1 of the
2271       discontinuity, respectively. M{jump(u)}, which is the difference of the
2272       solution at side 1 and at side 0, denotes the jump of M{u} across
2273       discontinuity along the normal calculated by L{jump<util.jump>}.
2274       The coefficient M{d_contact} is of rank two and M{y_contact} is of rank one
2275       both in L{FunctionOnContactZero<escript.FunctionOnContactZero>} or
2276       L{FunctionOnContactOne<escript.FunctionOnContactOne>}.
2277       The coefficient M{d_contact_reduced} is of rank two and M{y_contact_reduced}
2278       is of rank one both in L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>}
2279       or L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.
2280       In case of a single PDE and a single component solution the contact
2281       condition takes the form
2282    
2283       M{n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)}
2284    
2285       In this case the coefficient M{d_contact} and M{y_contact} are each scalar
2286       both in L{FunctionOnContactZero<escript.FunctionOnContactZero>} or
2287       L{FunctionOnContactOne<escript.FunctionOnContactOne>} and the coefficient
2288       M{d_contact_reduced} and M{y_contact_reduced} are each scalar both in
2289       L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or
2290       L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.
2291    
2292       Typical usage::
2293    
2294           p = LinearPDE(dom)
2295           p.setValue(A=kronecker(dom), D=1, Y=0.5)
2296           u = p.getSolution()
2297    
2298       """
2299    
2300       def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):
2301       """       """
2302       reassess the matrix type and, if needed, initiates an operator rebuild       Initializes a new linear PDE.
2303    
2304         @param domain: domain of the PDE
2305         @type domain: L{Domain<escript.Domain>}
2306         @param numEquations: number of equations. If C{None} the number of
2307                              equations is extracted from the PDE coefficients.
2308         @param numSolutions: number of solution components. If C{None} the number
2309                              of solution components is extracted from the PDE
2310                              coefficients.
2311         @param debug: if True debug information is printed
2312    
2313         """
2314         super(LinearPDE, self).__init__(domain,numEquations,numSolutions,debug)
2315         #
2316         #   the coefficients of the PDE:
2317         #
2318         self.introduceCoefficients(
2319           A=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2320           B=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2321           C=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2322           D=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2323           X=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
2324           Y=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2325           d=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2326           y=PDECoef(PDECoef.BOUNDARY,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2327           d_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2328           y_contact=PDECoef(PDECoef.CONTACT,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2329           A_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2330           B_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2331           C_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION,PDECoef.BY_DIM),PDECoef.OPERATOR),
2332           D_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2333           X_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_DIM),PDECoef.RIGHTHANDSIDE),
2334           Y_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2335           d_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2336           y_reduced=PDECoef(PDECoef.BOUNDARY_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2337           d_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,PDECoef.BY_SOLUTION),PDECoef.OPERATOR),
2338           y_contact_reduced=PDECoef(PDECoef.CONTACT_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2339           r=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.RIGHTHANDSIDE),
2340           q=PDECoef(PDECoef.SOLUTION,(PDECoef.BY_SOLUTION,),PDECoef.BOTH) )
2341    
2342       def __str__(self):
2343         """
2344         Returns the string representation of the PDE.
2345    
2346         @return: a simple representation of the PDE
2347         @rtype: C{str}
2348       """       """
2349       new_matrix_type=self.getDomain().getSystemMatrixTypeId(self.getSolverMethod(),self.isSymmetric())       return "<LinearPDE %d>"%id(self)
2350       if not new_matrix_type==self.__matrix_type:  
2351           if self.debug() : print "PDE Debug: Matrix type is now %d."%new_matrix_type     def getRequiredOperatorType(self):
2352           self.__matrix_type=new_matrix_type        """
2353           self.__rebuildOperator(deep=True)        Returns the system type which needs to be used by the current set up.
2354          """
2355     #============ assembling =======================================================        solver_options=self.getSolverOptions()
2356     def __copyConstraint(self):        return self.getDomain().getSystemMatrixTypeId(solver_options.getSolverMethod(), solver_options.getPreconditioner(),solver_options.getPackage(), solver_options.isSymmetric())
2357        """  
2358        copies the constrint condition into u     def checkSymmetry(self,verbose=True):
2359        """        """
2360        if not self.__righthandside.isEmpty():        Tests the PDE for symmetry.
          q=self.getCoefficientOfPDE("q")  
          r=self.getCoefficientOfPDE("r")  
          if not q.isEmpty():  
              if r.isEmpty():  
                 r2=escript.Data(0,self.__righthandside.getShape(),self.__righthandside.getFunctionSpace())  
              else:  
                 r2=escript.Data(r,self.__righthandside.getFunctionSpace())  
              self.__righthandside.copyWithMask(r2,escript.Data(q,self.__righthandside.getFunctionSpace()))  
2361    
2362     def __applyConstraint(self):        @param verbose: if set to True or not present a report on coefficients
2363                          which break the symmetry is printed.
2364          @type verbose: C{bool}
2365          @return: True if the PDE is symmetric
2366          @rtype: L{bool}
2367          @note: This is a very expensive operation. It should be used for
2368                 degugging only! The symmetry flag is not altered.
2369          """
2370          out=True
2371          out=out and self.checkSymmetricTensor("A", verbose)
2372          out=out and self.checkSymmetricTensor("A_reduced", verbose)
2373          out=out and self.checkReciprocalSymmetry("B","C", verbose)
2374          out=out and self.checkReciprocalSymmetry("B_reduced","C_reduced", verbose)
2375          out=out and self.checkSymmetricTensor("D", verbose)
2376          out=out and self.checkSymmetricTensor("D_reduced", verbose)
2377          out=out and self.checkSymmetricTensor("d", verbose)
2378          out=out and self.checkSymmetricTensor("d_reduced", verbose)
2379          out=out and self.checkSymmetricTensor("d_contact", verbose)
2380          out=out and self.checkSymmetricTensor("d_contact_reduced", verbose)
2381          return out
2382    
2383       def createOperator(self):
2384         """         """
2385         applies the constraints defined by q and r to the system         Returns an instance of a new operator.
2386         """         """
2387         q=self.getCoefficientOfPDE("q")         optype=self.getRequiredOperatorType()
2388         r=self.getCoefficientOfPDE("r")         self.trace("New operator of type %s is allocated."%optype)
2389         if not q.isEmpty() and not self.__operator.isEmpty():         return self.getDomain().newOperator( \
2390            # q is the row and column mask to indicate where constraints are set:                             self.getNumEquations(), \
2391            row_q=escript.Data(q,self.getFunctionSpaceForEquation())                             self.getFunctionSpaceForEquation(), \
2392            col_q=escript.Data(q,self.getFunctionSpaceForSolution())                             self.getNumSolutions(), \
2393            u=self.__getNewSolution()                             self.getFunctionSpaceForSolution(), \
2394            if r.isEmpty():                             optype)
2395               r_s=self.__getNewSolution()  
2396            else:     def getSolution(self):
2397               r_s=escript.Data(r,self.getFunctionSpaceForSolution())         """
2398            u.copyWithMask(r_s,col_q)         Returns the solution of the PDE.
2399    
2400           @return: the solution
2401           @rtype: L{Data<escript.Data>}
2402           """
2403           option_class=self.getSolverOptions()
2404           if not self.isSolutionValid():
2405              mat,f=self.getSystem()
2406            if self.isUsingLumping():            if self.isUsingLumping():
2407               self.__operator.copyWithMask(escript.Data(1,q.getShape(),self.getFunctionSpaceForEquation()),row_q)               self.setSolution(f*1/mat)
2408            else:            else:
2409               if not self.__righthandside.isEmpty(): self.__righthandside-=self.__operator*u               self.trace("PDE is resolved.")
2410               self.__operator.nullifyRowsAndCols(row_q,col_q,1.)               self.trace("solver options: %s"%str(option_class))
2411                 self.setSolution(mat.solve(f,option_class))
2412           return self.getCurrentSolution()
2413    
2414     def getSystem(self):     def getSystem(self):
2415         """         """
2416         return the operator and right hand side of the PDE         Returns the operator and right hand side of the PDE.
2417    
2418           @return: the discrete version of the PDE
2419           @rtype: C{tuple} of L{Operator,<escript.Operator>} and
2420                   L{Data<escript.Data>}
2421         """         """
2422         if not self.__operator_isValid or not self.__righthandside_isValid:         if not self.isOperatorValid() or not self.isRightHandSideValid():
2423            if self.isUsingLumping():            if self.isUsingLumping():
2424                if not self.__operator_isValid:                if not self.isOperatorValid():
2425                   if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():                   if not self.getFunctionSpaceForEquation()==self.getFunctionSpaceForSolution():
2426                         raise TypeError,"Lumped matrix requires same order for equations and unknowns"                        raise TypeError,"Lumped matrix requires same order for equations and unknowns"
2427                   if not self.getCoefficientOfPDE("A").isEmpty():                   if not self.getCoefficient("A").isEmpty():
2428                            raise Warning,"Lumped matrix does not allow coefficient A"                        raise ValueError,"coefficient A in lumped matrix may not be present."
2429                   if not self.getCoefficientOfPDE("B").isEmpty():                   if not self.getCoefficient("B").isEmpty():
2430                            raise Warning,"Lumped matrix does not allow coefficient B"                        raise ValueError,"coefficient B in lumped matrix may not be present."
2431                   if not self.getCoefficientOfPDE("C").isEmpty():                   if not self.getCoefficient("C").isEmpty():
2432                            raise Warning,"Lumped matrix does not allow coefficient C"                        raise ValueError,"coefficient C in lumped matrix may not be present."
2433                   if self.debug() : print "PDE Debug: New lumped operator is built."                   if not self.getCoefficient("d_contact").isEmpty():
2434                   mat=self.__getNewOperator()                        raise ValueError,"coefficient d_contact in lumped matrix may not be present."
2435                   self.getDomain().addPDEToSystem(mat,escript.Data(), \                   if not self.getCoefficient("A_reduced").isEmpty():
2436                             self.getCoefficientOfPDE("A"), \                        raise ValueError,"coefficient A_reduced in lumped matrix may not be present."
2437                             self.getCoefficientOfPDE("B"), \                   if not self.getCoefficient("B_reduced").isEmpty():
2438                             self.getCoefficientOfPDE("C"), \                        raise ValueError,"coefficient B_reduced in lumped matrix may not be present."
2439                             self.getCoefficientOfPDE("D"), \                   if not self.getCoefficient("C_reduced").isEmpty():
2440                             escript.Data(), \                        raise ValueError,"coefficient C_reduced in lumped matrix may not be present."
2441                             escript.Data(), \                   if not self.getCoefficient("d_contact_reduced").isEmpty():
2442                             self.getCoefficientOfPDE("d"), \                        raise ValueError,"coefficient d_contact_reduced in lumped matrix may not be present."
2443                             escript.Data(),\                   D=self.getCoefficient("D")
2444                             self.getCoefficientOfPDE("d_contact"), \                   d=self.getCoefficient("d")
2445                             escript.Data())                   D_reduced=self.getCoefficient("D_reduced")
2446                   self.__operator=mat*escript.Data(1,(self.getNumSolutions(),),self.getFunctionSpaceForSolution(),True)                   d_reduced=self.getCoefficient("d_reduced")
2447                   self.__applyConstraint()                   if not D.isEmpty():
2448                   self.__operator_isValid=True                       if self.getNumSolutions()>1:
2449                if not self.__righthandside_isValid:                          D_times_e=util.matrix_mult(D,numpy.ones((self.getNumSolutions(),)))
2450                   if self.debug() : print "PDE Debug: New right hand side is built."                       else:
2451                   self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \                          D_times_e=D
2452                                 self.getCoefficientOfPDE("X"), \                   else:
2453                                 self.getCoefficientOfPDE("Y"),\                      D_times_e=escript.Data()
2454                                 self.getCoefficientOfPDE("y"),\                   if not d.isEmpty():
2455                                 self.getCoefficientOfPDE("y_contact"))                       if self.getNumSolutions()>1:
2456                   self.__copyConstraint()                          d_times_e=util.matrix_mult(d,numpy.ones((self.getNumSolutions(),)))
2457                   self.__righthandside_isValid=True                       else:
2458                            d_times_e=d
2459                     else:
2460                        d_times_e=escript.Data()
2461    
2462                     if not D_reduced.isEmpty():
2463                         if self.getNumSolutions()>1:
2464                            D_reduced_times_e=util.matrix_mult(D_reduced,numpy.ones((self.getNumSolutions(),)))
2465                         else:
2466                            D_reduced_times_e=D_reduced
2467                     else:
2468                        D_reduced_times_e=escript.Data()
2469                     if not d_reduced.isEmpty():
2470                         if self.getNumSolutions()>1:
2471                            d_reduced_times_e=util.matrix_mult(d_reduced,numpy.ones((self.getNumSolutions(),)))
2472                         else:
2473                            d_reduced_times_e=d_reduced
2474                     else:
2475                        d_reduced_times_e=escript.Data()
2476    
2477                     self.resetOperator()
2478                     operator=self.getCurrentOperator()
2479                     if False and hasattr(self.getDomain(), "addPDEToLumpedSystem") :
2480                        self.getDomain().addPDEToLumpedSystem(operator, D_times_e, d_times_e)
2481                        self.getDomain().addPDEToLumpedSystem(operator, D_reduced_times_e, d_reduced_times_e)
2482                     else:
2483                        self.getDomain().addPDEToRHS(operator, \
2484                                                     escript.Data(), \
2485                                                     D_times_e, \
2486                                                     d_times_e,\
2487                                                     escript.Data())
2488                        self.getDomain().addPDEToRHS(operator, \
2489                                                     escript.Data(), \
2490                                                     D_reduced_times_e, \
2491                                                     d_reduced_times_e,\
2492                                                     escript.Data())
2493                     self.trace("New lumped operator has been built.")
2494                  if not self.isRightHandSideValid():
2495                     self.resetRightHandSide()
2496                     righthandside=self.getCurrentRightHandSide()
2497                     self.getDomain().addPDEToRHS(righthandside, \
2498                                   self.getCoefficient("X"), \
2499                                   self.getCoefficient("Y"),\
2500                                   self.getCoefficient("y"),\
2501                                   self.getCoefficient("y_contact"))
2502                     self.getDomain().addPDEToRHS(righthandside, \
2503                                   self.getCoefficient("X_reduced"), \
2504                                   self.getCoefficient("Y_reduced"),\
2505                                   self.getCoefficient("y_reduced"),\
2506                                   self.getCoefficient("y_contact_reduced"))
2507                     self.trace("New right hand side has been built.")
2508                     self.validRightHandSide()
2509                  self.insertConstraint(rhs_only=False)
2510                  self.validOperator()
2511            else:            else:
2512               if not self.__operator_isValid and not self.__righthandside_isValid:               if not self.isOperatorValid() and not self.isRightHandSideValid():
2513                   if self.debug() : print "PDE Debug: New system is built."                   self.resetRightHandSide()
2514                   self.getDomain().addPDEToSystem(self.__makeFreshOperator(),self.__makeFreshRightHandSide(), \                   righthandside=self.getCurrentRightHandSide()
2515                                 self.getCoefficientOfPDE("A"), \                   self.resetOperator()
2516                                 self.getCoefficientOfPDE("B"), \                   operator=self.getCurrentOperator()
2517                                 self.getCoefficientOfPDE("C"), \                   self.getDomain().addPDEToSystem(operator,righthandside, \
2518                                 self.getCoefficientOfPDE("D"), \                                 self.getCoefficient("A"), \
2519                                 self.getCoefficientOfPDE("X"), \                                 self.getCoefficient("B"), \
2520                                 self.getCoefficientOfPDE("Y"), \                                 self.getCoefficient("C"), \
2521                                 self.getCoefficientOfPDE("d"), \                                 self.getCoefficient("D"), \
2522                                 self.getCoefficientOfPDE("y"), \                                 self.getCoefficient("X"), \
2523                                 self.getCoefficientOfPDE("d_contact"), \                                 self.getCoefficient("Y"), \
2524                                 self.getCoefficientOfPDE("y_contact"))                                 self.getCoefficient("d"), \
2525                   self.__applyConstraint()                                 self.getCoefficient("y"), \
2526                   self.__copyConstraint()                                 self.getCoefficient("d_contact"), \
2527                   self.__operator_isValid=True                                 self.getCoefficient("y_contact"))
2528                   self.__righthandside_isValid=True                   self.getDomain().addPDEToSystem(operator,righthandside, \
2529               elif not self.__righthandside_isValid:                                 self.getCoefficient("A_reduced"), \
2530                   if self.debug() : print "PDE Debug: New right hand side is built."                                 self.getCoefficient("B_reduced"), \
2531                   self.getDomain().addPDEToRHS(self.__makeFreshRightHandSide(), \                                 self.getCoefficient("C_reduced"), \
2532                                 self.getCoefficientOfPDE("X"), \                                 self.getCoefficient("D_reduced"), \
2533                                 self.getCoefficientOfPDE("Y"),\                                 self.getCoefficient("X_reduced"), \
2534                                 self.getCoefficientOfPDE("y"),\                                 self.getCoefficient("Y_reduced"), \
2535                                 self.getCoefficientOfPDE("y_contact"))                                 self.getCoefficient("d_reduced"), \
2536                   self.__copyConstraint()                                 self.getCoefficient("y_reduced"), \
2537                   self.__righthandside_isValid=True                                 self.getCoefficient("d_contact_reduced"), \
2538               elif not self.__operator_isValid:                                 self.getCoefficient("y_contact_reduced"))
2539                   if self.debug() : print "PDE Debug: New operator is built."                   self.insertConstraint(rhs_only=False)
2540                   self.getDomain().addPDEToSystem(self.__makeFreshOperator(),escript.Data(), \                   self.trace("New system has been built.")
2541                              self.getCoefficientOfPDE("A"), \                   self.validOperator()
2542                              self.getCoefficientOfPDE("B"), \                   self.validRightHandSide()
2543                              self.getCoefficientOfPDE("C"), \               elif not self.isRightHandSideValid():
2544                              self.getCoefficientOfPDE("D"), \                   self.resetRightHandSide()
2545                     righthandside=self.getCurrentRightHandSide()
2546                     self.getDomain().addPDEToRHS(righthandside,
2547                                   self.getCoefficient("X"), \
2548                                   self.getCoefficient("Y"),\
2549                                   self.getCoefficient("y"),\
2550                                   self.getCoefficient("y_contact"))
2551                     self.getDomain().addPDEToRHS(righthandside,
2552                                   self.getCoefficient("X_reduced"), \
2553                                   self.getCoefficient("Y_reduced"),\
2554                                   self.getCoefficient("y_reduced"),\
2555                                   self.getCoefficient("y_contact_reduced"))
2556                     self.insertConstraint(rhs_only=True)
2557                     self.trace("New right hand side has been built.")
2558                     self.validRightHandSide()
2559                 elif not self.isOperatorValid():
2560                     self.resetOperator()
2561                     operator=self.getCurrentOperator()
2562                     self.getDomain().addPDEToSystem(operator,escript.Data(), \
2563                                self.getCoefficient("A"), \
2564                                self.getCoefficient("B"), \
2565                                self.getCoefficient("C"), \
2566                                self.getCoefficient("D"), \
2567                              escript.Data(), \                              escript.Data(), \
2568                              escript.Data(), \                              escript.Data(), \
2569                              self.getCoefficientOfPDE("d"), \                              self.getCoefficient("d"), \
2570                              escript.Data(),\                              escript.Data(),\
2571                              self.getCoefficientOfPDE("d_contact"), \                              self.getCoefficient("d_contact"), \
2572                              escript.Data())                              escript.Data())
2573                   self.__applyConstraint()                   self.getDomain().addPDEToSystem(operator,escript.Data(), \
2574                   self.__operator_isValid=True                              self.getCoefficient("A_reduced"), \
2575         return (self.__operator,self.__righthandside)                              self.getCoefficient("B_reduced"), \
2576     def getOperator(self):                              self.getCoefficient("C_reduced"), \
2577         """                              self.getCoefficient("D_reduced"), \
2578         returns the operator of the PDE                              escript.Data(), \
2579         """                              escript.Data(), \
2580         return self.getSystem()[0]                              self.getCoefficient("d_reduced"), \
2581                                escript.Data(),\
2582     def getRightHandSide(self):                              self.getCoefficient("d_contact_reduced"), \
2583         """                              escript.Data())
2584         returns the right hand side of the PDE                   self.insertConstraint(rhs_only=False)
2585         """                   self.trace("New operator has been built.")
2586         return self.getSystem()[1]                   self.validOperator()
2587           return (self.getCurrentOperator(), self.getCurrentRightHandSide())
2588    
2589       def insertConstraint(self, rhs_only=False):
2590          """
2591          Applies the constraints defined by q and r to the PDE.
2592    
2593          @param rhs_only: if True only the right hand side is altered by the
2594                           constraint
2595          @type rhs_only: C{bool}
2596          """
2597          q=self.getCoefficient("q")
2598          r=self.getCoefficient("r")
2599          righthandside=self.getCurrentRightHandSide()
2600          operator=self.getCurrentOperator()
2601    
2602          if not q.isEmpty():
2603             if r.isEmpty():
2604                r_s=self.createSolution()
2605             else:
2606                r_s=r
2607             if not rhs_only and not operator.isEmpty():
2608                 if self.isUsingLumping():
2609                     operator.copyWithMask(escript.Data(1.,q.getShape(),q.getFunctionSpace()),q)
2610                 else:
2611                     row_q=escript.Data(q,self.getFunctionSpaceForEquation())
2612                     col_q=escript.Data(q,self.getFunctionSpaceForSolution())
2613                     u=self.createSolution()
2614                     u.copyWithMask(r_s,col_q)
2615                     righthandside-=operator*u
2616                     operator.nullifyRowsAndCols(row_q,col_q,1.)
2617             righthandside.copyWithMask(r_s,q)
2618    
2619     def solve(self,**options):     def setValue(self,**coefficients):
2620        """        """
2621        solve the PDE        Sets new values to coefficients.
2622    
2623        @param options:        @param coefficients: new values assigned to coefficients
2624        """        @keyword A: value for coefficient C{A}
2625        mat,f=self.getSystem()        @type A: any type that can be cast to a L{Data<escript.Data>} object on
2626        if self.isUsingLumping():                 L{Function<escript.Function>}
2627           out=f/mat        @keyword A_reduced: value for coefficient C{A_reduced}
2628        else:        @type A_reduced: any type that can be cast to a L{Data<escript.Data>}
2629           options[util.TOLERANCE_KEY]=self.getTolerance()                         object on L{ReducedFunction<escript.ReducedFunction>}
2630           options[util.METHOD_KEY]=self.getSolverMethod()        @keyword B: value for coefficient C{B}
2631           options[util.SYMMETRY_KEY]=self.isSymmetric()        @type B: any type that can be cast to a L{Data<escript.Data>} object on
2632           if self.debug() : print "PDE Debug: solver options: ",options                 L{Function<escript.Function>}
2633           out=mat.solve(f,options)        @keyword B_reduced: value for coefficient C{B_reduced}
2634        return out        @type B_reduced: any type that can be cast to a L{Data<escript.Data>}
2635                           object on L{ReducedFunction<escript.ReducedFunction>}
2636          @keyword C: value for coefficient C{C}
2637          @type C: any type that can be cast to a L{Data<escript.Data>} object on
2638                   L{Function<escript.Function>}
2639          @keyword C_reduced: value for coefficient C{C_reduced}
2640          @type C_reduced: any type that can be cast to a L{Data<escript.Data>}
2641                           object on L{ReducedFunction<escript.ReducedFunction>}
2642          @keyword D: value for coefficient C{D}
2643          @type D: any type that can be cast to a L{Data<escript.Data>} object on
2644                   L{Function<escript.Function>}
2645          @keyword D_reduced: value for coefficient C{D_reduced}
2646          @type D_reduced: any type that can be cast to a L{Data<escript.Data>}
2647                           object on L{ReducedFunction<escript.ReducedFunction>}
2648          @keyword X: value for coefficient C{X}
2649          @type X: any type that can be cast to a L{Data<escript.Data>} object on
2650                   L{Function<escript.Function>}
2651          @keyword X_reduced: value for coefficient C{X_reduced}
2652          @type X_reduced: any type that can be cast to a L{Data<escript.Data>}
2653                           object on L{ReducedFunction<escript.ReducedFunction>}
2654          @keyword Y: value for coefficient C{Y}
2655          @type Y: any type that can be cast to a L{Data<escript.Data>} object on
2656                   L{Function<escript.Function>}
2657          @keyword Y_reduced: value for coefficient C{Y_reduced}
2658          @type Y_reduced: any type that can be cast to a L{Data<escript.Data>}
2659                           object on L{ReducedFunction<escript.Function>}
2660          @keyword d: value for coefficient C{d}
2661          @type d: any type that can be cast to a L{Data<escript.Data>} object on
2662                   L{FunctionOnBoundary<escript.FunctionOnBoundary>}
2663          @keyword d_reduced: value for coefficient C{d_reduced}
2664          @type d_reduced: any type that can be cast to a L{Data<escript.Data>}
2665                           object on L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}
2666          @keyword y: value for coefficient C{y}
2667          @type y: any type that can be cast to a L{Data<escript.Data>} object on
2668                   L{FunctionOnBoundary<escript.FunctionOnBoundary>}
2669          @keyword d_contact: value for coefficient C{d_contact}
2670          @type d_contact: any type that can be cast to a L{Data<escript.Data>}
2671                           object on L{FunctionOnContactOne<escript.FunctionOnContactOne>}
2672                           or L{FunctionOnContactZero<escript.FunctionOnContactZero>}
2673          @keyword d_contact_reduced: value for coefficient C{d_contact_reduced}
2674          @type d_contact_reduced: any type that can be cast to a L{Data<escript.Data>}
2675                                   object on L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}
2676                                   or L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>}
2677          @keyword y_contact: value for coefficient C{y_contact}
2678          @type y_contact: any type that can be cast to a L{Data<escript.Data>}
2679                           object on L{FunctionOnContactOne<escript.FunctionOnContactOne>}
2680                           or L{FunctionOnContactZero<escript.FunctionOnContactZero>}
2681          @keyword y_contact_reduced: value for coefficient C{y_contact_reduced}
2682          @type y_contact_reduced: any type that can be cast to a L{Data<escript.Data>}
2683                                   object on L{ReducedFunctionOnContactOne<escript.FunctionOnContactOne>}
2684                                   or L{ReducedFunctionOnContactZero<escript.FunctionOnContactZero>}
2685          @keyword r: values prescribed to the solution at the locations of
2686                      constraints
2687          @type r: any type that can be cast to a L{Data<escript.Data>} object on
2688                   L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2689                   depending on whether reduced order is used for the solution
2690          @keyword q: mask for location of constraints
2691          @type q: any type that can be cast to a L{Data<escript.Data>} object on
2692                   L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}
2693                   depending on whether reduced order is used for the
2694                   representation of the equation
2695          @raise IllegalCoefficient: if an unknown coefficient keyword is used
2696          """
2697          super(LinearPDE,self).setValue(**coefficients)
2698          # check if the systrem is inhomogeneous:
2699          if len(coefficients)>0 and not self.isUsingLumping():
2700             q=self.getCoefficient("q")
2701             r=self.getCoefficient("r")
2702             if not q.isEmpty() and not r.isEmpty():
2703                 if util.Lsup(q*r)>0.:
2704                   self.trace("Inhomogeneous constraint detected.")
2705                   self.invalidateSystem()
2706    
2707    
2708       def getResidual(self,u=None):
2709         """
2710         Returns the residual of u or the current solution if u is not present.
2711    
2712         @param u: argument in the residual calculation. It must be representable
2713                   in L{self.getFunctionSpaceForSolution()}. If u is not present
2714                   or equals C{None} the current solution is used.
2715         @type u: L{Data<escript.Data>} or None
2716         @return: residual of u
2717         @rtype: L{Data<escript.Data>}
2718         """
2719         if u==None:
2720            return self.getOperator()*self.getSolution()-self.getRightHandSide()
2721         else:
2722            return self.getOperator()*escript.Data(u,self.getFunctionSpaceForSolution())-self.getRightHandSide()
2723    
2724     def getSolution(self,**options):     def getFlux(self,u=None):
2725         """       """
2726         returns the solution of the PDE       Returns the flux M{J} for a given M{u}.
2727    
2728         @param options:       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]}
        """  
        if not self.__solution_isValid:  
            if self.debug() : print "PDE Debug: PDE is resolved."  
            self.__solution=self.solve(**options)  
            self.__solution_isValid=True  
        return self.__solution  
2729    
2730         or
2731    
2732         M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[l]+(B[j]+B_reduced[j])u-X[j]-X_reduced[j]}
2733    
2734         @param u: argument in the flux. If u is not present or equals L{None} the
2735                   current solution is used.
2736         @type u: L{Data<escript.Data>} or None
2737         @return: flux
2738         @rtype: L{Data<escript.Data>}
2739         """
2740         if u==None: u=self.getSolution()
2741         return util.tensormult(self.getCoefficient("A"),util.grad(u,Funtion(self.getDomain))) \
2742               +util.matrixmult(self.getCoefficient("B"),u) \
2743               -util.self.getCoefficient("X") \
2744               +util.tensormult(self.getCoefficient("A_reduced"),util.grad(u,ReducedFuntion(self.getDomain))) \
2745               +util.matrixmult(self.getCoefficient("B_reduced"),u) \
2746               -util.self.getCoefficient("X_reduced")
2747    
 def ELMAN_RAMAGE(P):  
      """   """  
      return (P-1.).wherePositive()*0.5*(1.-1./(P+1.e-15))  
 def SIMPLIFIED_BROOK_HUGHES(P):  
      """   """  
      c=(P-3.).whereNegative()  
      return P/6.*c+1./2.*(1.-c)  
   
 def HALF(P):  
     """ """  
     return escript.Scalar(0.5,P.getFunctionSpace())  
2748    
2749  class AdvectivePDE(LinearPDE):  class Poisson(LinearPDE):
2750     """     """
2751     Class to handle a linear PDE dominated by advective terms:     Class to define a Poisson equation problem. This is generally a
2752         L{LinearPDE} of the form
    class to define a linear PDE of the form  
2753    
2754     \f[     M{-grad(grad(u)[j])[j] = f}
    -(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  
    \f]  
2755    
2756     with boundary conditons:     with natural boundary conditions
2757    
2758     \f[     M{n[j]*grad(u)[j] = 0 }
    n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d_{ik}u_k = - n_j*X_{ij} + y_i  
    \f]  
2759    
2760     and contact conditions     and constraints:
2761    
2762     \f[     M{u=0} where M{q>0}
    n_j*(A_{ijkl}u_{k,l}+B_{ijk}u_k)_{,j} + d^{contact}_{ik}[u_k] = - n_j*X_{ij} + y^{contact}_{i}  
    \f]  
2763    
2764     and constraints:     """
2765    
2766     \f[     def __init__(self,domain,debug=False):
2767     u_i=r_i \quad \mathrm{where} \quad q_i>0       """
2768     \f]       Initializes a new Poisson equation.
    """  
    def __init__(self,domain,numEquations=0,numSolutions=0,xi=ELMAN_RAMAGE):  
       LinearPDE.__init__(self,domain,numEquations,numSolutions)  
       self.__xi=xi  
       self.__Xi=escript.Data()  
   
    def __calculateXi(self,peclet_factor,Z,h):  
        Z_max=util.Lsup(Z)  
        if Z_max>0.:  
           return h*self.__xi(Z*peclet_factor)/(Z+Z_max*self.TOL)  
        else:  
           return 0.  
2769    
2770     def setValue(self,**args):       @param domain: domain of the PDE
2771         if "A" in args.keys()   or "B" in args.keys() or "C" in args.keys(): self.__Xi=escript.Data()       @type domain: L{Domain<escript.Domain>}
2772         self._LinearPDE__setValue(**args)       @param debug: if True debug information is printed
             
    def getXi(self):  
       if self.__Xi.isEmpty():  
          B=self.getCoefficient("B")  
          C=self.getCoefficient("C")  
          A=self.getCoefficient("A")  
          h=self.getDomain().getSize()  
          self.__Xi=escript.Scalar(0.,self.getFunctionSpaceForCoefficient("A"))  
          if not C.isEmpty() or not B.isEmpty():  
             if not C.isEmpty() and not B.isEmpty():  
                 Z2=escript.Scalar(0,self.getFunctionSpaceForCoefficient("A"))  
                 if self.getNumEquations()>1:  
                    if self.getNumSolutions()>1:  
                       for i in range(self.getNumEquations()):  
                          for k in range(self.getNumSolutions()):  
                             for l in range(self.getDim()): Z2+=(C[i,k,l]-B[i,l,k])**2  
                    else:  
                       for i in range(self.getNumEquations()):  
                          for l in range(self.getDim()): Z2+=(C[i,l]-B[i,l])**2  
                 else:  
                    if self.getNumSolutions()>1:  
                       for k in range(self.getNumSolutions()):  
                          for l in range(self.getDim()): Z2+=(C[k,l]-B[l,k])**2  
                    else:  
                       for l in range(self.getDim()): Z2+=(C[l]-B[l])**2  
                 length_of_Z=util.sqrt(Z2)  
             elif C.isEmpty():  
               length_of_Z=util.length(B)  
             else:  
               length_of_Z=util.length(C)  
2773    
2774              Z_max=util.Lsup(length_of_Z)       """
2775              if Z_max>0.:       super(Poisson, self).__init__(domain,1,1,debug)
2776                 length_of_A=util.length(A)       self.introduceCoefficients(
2777                 A_max=util.Lsup(length_of_A)                          f=PDECoef(PDECoef.INTERIOR,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE),
2778                 if A_max>0:                          f_reduced=PDECoef(PDECoef.INTERIOR_REDUCED,(PDECoef.BY_EQUATION,),PDECoef.RIGHTHANDSIDE))
2779                      inv_A=1./(length_of_A+A_max*self.TOL)       self.setSymmetryOn()
                else:  
                     inv_A=1./self.TOL  
                peclet_number=length_of_Z*h/2*inv_A  
                xi=self.__xi(peclet_number)  
                self.__Xi=h*xi/(length_of_Z+Z_max*self.TOL)  
                print "@ preclet number = %e"%util.Lsup(peclet_number),util.Lsup(xi),util.Lsup(length_of_Z)  
       return self.__Xi  
         
2780    
2781     def getCoefficientOfPDE(self,name):     def setValue(self,**coefficients):
2782       """       """
2783       return the value of the coefficient name of the general PDE       Sets new values to coefficients.
2784    
2785       @param name:       @param coefficients: new values assigned to coefficients
2786         @keyword f: value for right hand side M{f}
2787         @type f: any type that can be cast to a L{Scalar<escript.Scalar>} object
2788                  on L{Function<escript.Function>}
2789         @keyword q: mask for location of constraints
2790         @type q: any type that can be cast to a rank zero L{Data<escript.Data>}
2791                  object on L{Solution<escript.Solution>} or
2792                  L{ReducedSolution<escript.ReducedSolution>} depending on whether
2793                  reduced order is used for the representation of the equation
2794         @raise IllegalCoefficient: if an unknown coefficient keyword is used
2795       """       """
2796       if not self.getNumEquations() == self.getNumSolutions():       super(Poisson, self).setValue(**coefficients)
           raise ValueError,"AdvectivePDE expects the number of solution componets and the number of equations to be equal."  
2797    
      if name == "A" :  
          A=self.getCoefficient("A")  
          B=self.getCoefficient("B")  
          C=self.getCoefficient("C")  
          if B.isEmpty() and C.isEmpty():  
             Aout=A  
          else:  
             if A.isEmpty():  
                Aout=self.createNewCoefficient("A")  
             else:  
                Aout=A[:]  
             Xi=self.getXi()  
             if self.getNumEquations()>1:  
                 for i in range(self.getNumEquations()):  
                    for j in range(self.getDim()):  
                       for k in range(self.getNumSolutions()):  
                          for l in range(self.getDim()):  
                             if not C.isEmpty() and not B.isEmpty():  
                                for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*(C[p,i,j]-B[p,j,i])*(C[p,k,l]-B[p,l,k])  
                             elif C.isEmpty():  
                                for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*B[p,j,i]*B[p,l,k]  
                             else:  
                                for p in range(self.getNumEquations()): Aout[i,j,k,l]+=Xi*C[p,i,j]*C[p,k,l]  
             else:  
                 for j in range(self.getDim()):  
                    for l in range(self.getDim()):  
                       if not C.isEmpty() and not B.isEmpty():  
                           Aout[j,l]+=Xi*(C[j]-B[j])*(C[l]-B[l])  
                       elif C.isEmpty():  
                           Aout[j,l]+=Xi*B[j]*B[l]  
                       else:  
                           Aout[j,l]+=Xi*C[j]*C[l]  
          return Aout  
      elif name == "B" :  
          B=self.getCoefficient("B")  
          C=self.getCoefficient("C")  
          D=self.getCoefficient("D")  
          if C.isEmpty() or D.isEmpty():  
             Bout=B  
          else:  
             Xi=self.getXi()  
             if B.isEmpty():  
                 Bout=self.createNewCoefficient("B")  
             else:  
                 Bout=B[:]  
             if self.getNumEquations()>1:  
                for k in range(self.getNumSolutions()):  
                   for p in range(self.getNumEquations()):  
                      tmp=Xi*D[p,k]  
                      for i in range(self.getNumEquations()):  
                         for j in range(self.getDim()):  
                            Bout[i,j,k]+=tmp*C[p,i,j]  
             else:  
                tmp=Xi*D  
                for j in range(self.getDim()): Bout[j]+=tmp*C[j]  
          return Bout  
      elif name == "C" :  
          B=self.getCoefficient("B")  
          C=self.getCoefficient("C")  
          D=self.getCoefficient("D")  
          if B.isEmpty() or D.isEmpty():  
             Cout=C  
          else:  
             Xi=self.getXi()  
             if C.isEmpty():  
                 Cout=self.createNewCoefficient("C")  
             else:  
                 Cout=C[:]  
             if self.getNumEquations()>1:  
                for k in range(self.getNumSolutions()):  
                    for p in range(self.getNumEquations()):  
                       tmp=Xi*D[p,k]  
                       for i in range(self.getNumEquations()):  
                         for l in range(self.getDim()):  
                                  Cout[i,k,l]+=tmp*B[p,l,i]  
             else:  
                tmp=Xi*D  
                for j in range(self.getDim()): Cout[j]+=tmp*B[j]