/[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

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