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