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