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