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