/[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

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