/[escript]/trunk/escript/py_src/linearPDEs.py
ViewVC logotype

Diff of /trunk/escript/py_src/linearPDEs.py

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

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