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