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