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