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

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

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

revision 2549 by jfenwick, Mon Jul 20 06:43:47 2009 UTC revision 3071 by gross, Wed Jul 21 05:37:30 2010 UTC
# Line 1  Line 1 
1    
2  ########################################################  ########################################################
3  #  #
4  # Copyright (c) 2003-2009 by University of Queensland  # Copyright (c) 2003-2010 by University of Queensland
5  # Earth Systems Science Computational Center (ESSCC)  # Earth Systems Science Computational Center (ESSCC)
6  # http://www.uq.edu.au/esscc  # http://www.uq.edu.au/esscc
7  #  #
# Line 11  Line 11 
11  #  #
12  ########################################################  ########################################################
13    
14  __copyright__="""Copyright (c) 2003-2009 by University of Queensland  __copyright__="""Copyright (c) 2003-2010 by University of Queensland
15  Earth Systems Science Computational Center (ESSCC)  Earth Systems Science Computational Center (ESSCC)
16  http://www.uq.edu.au/esscc  http://www.uq.edu.au/esscc
17  Primary Business: Queensland, Australia"""  Primary Business: Queensland, Australia"""
# Line 21  __url__="https://launchpad.net/escript-f Line 21  __url__="https://launchpad.net/escript-f
21    
22  """  """
23  The module provides an interface to define and solve linear partial  The module provides an interface to define and solve linear partial
24  differential equations (PDEs) and Transport problems within L{escript}.  differential equations (PDEs) and Transport problems within `escript`.
25  L{linearPDEs} does not provide any solver capabilities in itself but hands the  `linearPDEs` does not provide any solver capabilities in itself but hands the
26  PDE over to the PDE solver library defined through the L{Domain<escript.Domain>}  PDE over to the PDE solver library defined through the `Domain`
27  of the PDE. The general interface is provided through the L{LinearPDE} class.  of the PDE. The general interface is provided through the `LinearPDE` class.
28  L{TransportProblem} provides an interface to initial value problems dominated  `TransportProblem` provides an interface to initial value problems dominated
29  by its advective terms.  by its advective terms.
30    
31  @var __author__: name of author  :var __author__: name of author
32  @var __copyright__: copyrights  :var __copyright__: copyrights
33  @var __license__: licence agreement  :var __license__: licence agreement
34  @var __url__: url entry point on documentation  :var __url__: url entry point on documentation
35  @var __version__: version  :var __version__: version
36  @var __date__: date of the version  :var __date__: date of the version
37  """  """
38    
39  import math  import math
# Line 52  class SolverOptions(object): Line 52  class SolverOptions(object):
52            
53      Typical usage is      Typical usage is
54            
55      opts=SolverOptions()      ::
56      print opts      
57      opts.resetDiagnostics()        opts=SolverOptions()
58      u=solver(opts)        print opts
59      print "number of iteration steps: =",opts.getDiagnostics("num_iter")        opts.resetDiagnostics()
60              u=solver(opts)
61          print "number of iteration steps: =",opts.getDiagnostics("num_iter")
62      @cvar DEFAULT: The default method used to solve the system of linear equations  
63      @cvar DIRECT: The direct solver based on LDU factorization      :cvar DEFAULT: The default method used to solve the system of linear equations
64      @cvar CHOLEVSKY: The direct solver based on LDLt factorization (can only be applied for symmetric PDEs)      :cvar DIRECT: The direct solver based on LDU factorization
65      @cvar PCG: The preconditioned conjugate gradient method (can only be applied for symmetric PDEs)      :cvar CHOLEVSKY: The direct solver based on LDLt factorization (can only be applied for symmetric PDEs)
66      @cvar CR: The conjugate residual method      :cvar PCG: The preconditioned conjugate gradient method (can only be applied for symmetric PDEs)
67      @cvar CGS: The conjugate gradient square method      :cvar CR: The conjugate residual method
68      @cvar BICGSTAB: The stabilized Bi-Conjugate Gradient method      :cvar CGS: The conjugate gradient square method
69      @cvar TFQMR: Transport Free Quasi Minimal Residual method      :cvar BICGSTAB: The stabilized Bi-Conjugate Gradient method
70      @cvar MINRES: Minimum residual method      :cvar TFQMR: Transpose Free Quasi Minimal Residual method
71      @cvar SSOR: The symmetric over-relaxation method      :cvar MINRES: Minimum residual method
72      @cvar ILU0: The incomplete LU factorization preconditioner with no fill-in      :cvar SSOR: The symmetric over-relaxation method
73      @cvar ILUT: The incomplete LU factorization preconditioner with fill-in      :cvar ILU0: The incomplete LU factorization preconditioner with no fill-in
74      @cvar JACOBI: The Jacobi preconditioner      :cvar ILUT: The incomplete LU factorization preconditioner with fill-in
75      @cvar GMRES: The Gram-Schmidt minimum residual method      :cvar JACOBI: The Jacobi preconditioner
76      @cvar PRES20: Special GMRES with restart after 20 steps and truncation after 5 residuals      :cvar GMRES: The Gram-Schmidt minimum residual method
77      @cvar LUMPING: Matrix lumping      :cvar PRES20: Special GMRES with restart after 20 steps and truncation after 5 residuals
78      @cvar NO_REORDERING: No matrix reordering allowed      :cvar LUMPING: Matrix lumping
79      @cvar MINIMUM_FILL_IN: Reorder matrix to reduce fill-in during factorization      :cvar NO_REORDERING: No matrix reordering allowed
80      @cvar NESTED_DISSECTION: Reorder matrix to improve load balancing during factorization      :cvar MINIMUM_FILL_IN: Reorder matrix to reduce fill-in during factorization
81      @cvar PASO: PASO solver package      :cvar NESTED_DISSECTION: Reorder matrix to improve load balancing during factorization
82      @cvar SCSL: SGI SCSL solver library      :cvar PASO: PASO solver package
83      @cvar MKL: Intel's MKL solver library      :cvar SCSL: SGI SCSL solver library
84      @cvar UMFPACK: The UMFPACK library      :cvar MKL: Intel's MKL solver library
85      @cvar TRILINOS: The TRILINOS parallel solver class library from Sandia National Labs      :cvar UMFPACK: The UMFPACK library
86      @cvar ITERATIVE: The default iterative solver      :cvar TRILINOS: The TRILINOS parallel solver class library from Sandia National Labs
87      @cvar AMG: Algebraic Multi Grid      :cvar ITERATIVE: The default iterative solver
88      @cvar REC_ILU: recursive ILU0      :cvar AMG: Algebraic Multi Grid
89      @cvar RILU: relaxed ILU0      :cvar AMLI: Algebraic Multi Level Iteration
90      @cvar GAUSS_SEIDEL: Gauss-Seidel solver      :cvar REC_ILU: recursive ILU0
91      @cvar DEFAULT_REORDERING: the reordering method recommended by the solver      :cvar RILU: relaxed ILU0
92      @cvar SUPER_LU: the Super_LU solver package      :cvar GAUSS_SEIDEL: Gauss-Seidel solver
93      @cvar PASTIX: the Pastix direct solver_package      :cvar GAUSS_SEIDEL_MPI: MPI versioned Gauss-Seidel solver
94      @cvar YAIR_SHAPIRA_COARSENING: AMG coarsening method by Yair-Shapira      :cvar DEFAULT_REORDERING: the reordering method recommended by the solver
95      @cvar RUGE_STUEBEN_COARSENING: AMG coarsening method by Ruge and Stueben      :cvar SUPER_LU: the Super_LU solver package
96      @cvar AGGREGATION_COARSENING: AMG coarsening using (symmetric) aggregation      :cvar PASTIX: the Pastix direct solver_package
97      @cvar MIN_COARSE_MATRIX_SIZE: minimum size of the coarsest level matrix to use direct solver.      :cvar YAIR_SHAPIRA_COARSENING: AMG and AMLI coarsening method by Yair-Shapira
98      @cvar NO_PRECONDITIONER: no preconditioner is applied.      :cvar RUGE_STUEBEN_COARSENING: AMG and AMLI coarsening method by Ruge and Stueben
99        :cvar AGGREGATION_COARSENING: AMG and AMLI coarsening using (symmetric) aggregation
100        :cvar STANDARD_COARSENING: AMG and AMLI standard coarsening using mesure of importance of the unknowns
101        :cvar MIN_COARSE_MATRIX_SIZE: minimum size of the coarsest level matrix to use direct solver.
102        :cvar NO_PRECONDITIONER: no preconditioner is applied.
103      """      """
104      DEFAULT= 0      DEFAULT= 0
105      DIRECT= 1      DIRECT= 1
# Line 134  class SolverOptions(object): Line 138  class SolverOptions(object):
138      AGGREGATION_COARSENING=35      AGGREGATION_COARSENING=35
139      NO_PRECONDITIONER=36      NO_PRECONDITIONER=36
140      MIN_COARSE_MATRIX_SIZE=37      MIN_COARSE_MATRIX_SIZE=37
141            AMLI=38
142        STANDARD_COARSENING=39
143        GAUSS_SEIDEL_MPI=40
144    
145      def __init__(self):      def __init__(self):
146          self.setLevelMax()          self.setLevelMax()
147          self.setCoarseningThreshold()          self.setCoarseningThreshold()
148            self.setSmoother()
149          self.setNumSweeps()          self.setNumSweeps()
150          self.setNumPreSweeps()          self.setNumPreSweeps()
151          self.setNumPostSweeps()          self.setNumPostSweeps()
# Line 194  class SolverOptions(object): Line 202  class SolverOptions(object):
202                  out+="\nMaximum number of levels = %s"%self.LevelMax()                  out+="\nMaximum number of levels = %s"%self.LevelMax()
203                  out+="\nCoarsening method = %s"%self.getName(self.getCoarsening())                  out+="\nCoarsening method = %s"%self.getName(self.getCoarsening())
204                  out+="\nCoarsening threshold = %e"%self.getMinCoarseMatrixSize()                  out+="\nCoarsening threshold = %e"%self.getMinCoarseMatrixSize()
205                    out+="\nSmoother = %e"%self.getName(self.getSmoother())
206                  out+="\nMinimum size of the coarsest level matrix = %e"%self.getCoarseningThreshold()                  out+="\nMinimum size of the coarsest level matrix = %e"%self.getCoarseningThreshold()
207                  out+="\nNumber of pre / post sweeps = %s / %s, %s"%(self.getNumPreSweeps(), self.getNumPostSweeps(), self.getNumSweeps())                  out+="\nNumber of pre / post sweeps = %s / %s, %s"%(self.getNumPreSweeps(), self.getNumPostSweeps(), self.getNumSweeps())
208              if self.getPreconditioner() == self.GAUSS_SEIDEL:              if self.getPreconditioner() == self.AMLI:
209                    out+="\nMaximum number of levels = %s"%self.LevelMax()
210                    out+="\nCoarsening method = %s"%self.getName(self.getCoarsening())
211                    out+="\nCoarsening threshold = %e"%self.getMinCoarseMatrixSize()
212                    out+="\nMinimum size of the coarsest level matrix = %e"%self.getCoarseningThreshold()
213                    out+="\nNumber of pre / post sweeps = %s / %s, %s"%(self.getNumPreSweeps(), self.getNumPostSweeps(), self.getNumSweeps())
214            if self.getPreconditioner() == self.GAUSS_SEIDEL:
215                    out+="\nNumber of sweeps = %s"%self.getNumSweeps()
216                if self.getPreconditioner() == self.GAUSS_SEIDEL_MPI:
217                  out+="\nNumber of sweeps = %s"%self.getNumSweeps()                  out+="\nNumber of sweeps = %s"%self.getNumSweeps()
218              if self.getPreconditioner() == self.ILUT:              if self.getPreconditioner() == self.ILUT:
219                  out+="\nDrop tolerance = %e"%self.getDropTolerance()                  out+="\nDrop tolerance = %e"%self.getDropTolerance()
# Line 209  class SolverOptions(object): Line 226  class SolverOptions(object):
226          """          """
227          returns the name of a given key          returns the name of a given key
228                    
229          @param key: a valid key          :param key: a valid key
230          """          """
231          if key == self.DEFAULT: return "DEFAULT"          if key == self.DEFAULT: return "DEFAULT"
232          if key == self.DIRECT: return "DIRECT"          if key == self.DIRECT: return "DIRECT"
# Line 233  class SolverOptions(object): Line 250  class SolverOptions(object):
250          if key == self.ITERATIVE: return "ITERATIVE"          if key == self.ITERATIVE: return "ITERATIVE"
251          if key == self.PASO: return "PASO"          if key == self.PASO: return "PASO"
252          if key == self.AMG: return "AMG"          if key == self.AMG: return "AMG"
253            if key == self.AMLI: return "AMLI"
254          if key == self.REC_ILU: return "REC_ILU"          if key == self.REC_ILU: return "REC_ILU"
255          if key == self.TRILINOS: return "TRILINOS"          if key == self.TRILINOS: return "TRILINOS"
256          if key == self.NONLINEAR_GMRES: return "NONLINEAR_GMRES"          if key == self.NONLINEAR_GMRES: return "NONLINEAR_GMRES"
257          if key == self.TFQMR: return "TFQMR"          if key == self.TFQMR: return "TFQMR"
258          if key == self.MINRES: return "MINRES"          if key == self.MINRES: return "MINRES"
259          if key == self.GAUSS_SEIDEL: return "GAUSS_SEIDEL"          if key == self.GAUSS_SEIDEL: return "GAUSS_SEIDEL"
260            if key == self.GAUSS_SEIDEL_MPI: return "GAUSS_SEIDEL_MPI"
261          if key == self.RILU: return "RILU"          if key == self.RILU: return "RILU"
262          if key == self.DEFAULT_REORDERING: return "DEFAULT_REORDERING"          if key == self.DEFAULT_REORDERING: return "DEFAULT_REORDERING"
263          if key == self.SUPER_LU: return "SUPER_LU"          if key == self.SUPER_LU: return "SUPER_LU"
264          if key == self.PASTIX: return "PASTIX"          if key == self.PASTIX: return "PASTIX"
265          if key == self.YAIR_SHAPIRA_COARSENING: return "YAIR_SHAPIRA_COARSENING"          if key == self.YAIR_SHAPIRA_COARSENING: return "YAIR_SHAPIRA_COARSENING"
266          if key == self.RUGE_STUEBEN_COARSENING: return "RUGE_STUEBEN_COARSENING"          if key == self.RUGE_STUEBEN_COARSENING: return "RUGE_STUEBEN_COARSENING"
267            if key == self.STANDARD_COARSENING: return "STANDARD_COARSENING"
268          if key == self.AGGREGATION_COARSENING: return "AGGREGATION_COARSENING"          if key == self.AGGREGATION_COARSENING: return "AGGREGATION_COARSENING"
269          if key == self.NO_PRECONDITIONER: return "NO_PRECONDITIONER"          if key == self.NO_PRECONDITIONER: return "NO_PRECONDITIONER"
270          if key == self.MIN_COARSE_MATRIX_SIZE: return "MIN_COARSE_MATRIX_SIZE"          if key == self.MIN_COARSE_MATRIX_SIZE: return "MIN_COARSE_MATRIX_SIZE"
# Line 253  class SolverOptions(object): Line 273  class SolverOptions(object):
273          """          """
274          resets the diagnostics          resets the diagnostics
275                    
276          @param all: if C{all} is C{True} all diagnostics including accumulative counters are reset.          :param all: if ``all`` is ``True`` all diagnostics including accumulative counters are reset.
277          @type all: C{bool}          :type all: ``bool``
278          """          """
279          self.__num_iter=None          self.__num_iter=None
280          self.__num_level=None          self.__num_level=None
281          self.__num_inner_iter=None          self.__num_inner_iter=None
282          self.__time=None          self.__time=None
283          self.__set_up_time=None          self.__set_up_time=None
284            self.__net_time=None
285          self.__residual_norm=None          self.__residual_norm=None
286          self.__converged=None          self.__converged=None
287          if all:          if all:
# Line 268  class SolverOptions(object): Line 289  class SolverOptions(object):
289              self.__cum_num_iter=0              self.__cum_num_iter=0
290              self.__cum_time=0              self.__cum_time=0
291              self.__cum_set_up_time=0              self.__cum_set_up_time=0
292                self.__cum_net_time=0
293    
294      def _updateDiagnostics(self, name, value):      def _updateDiagnostics(self, name, value):
295          """          """
296          Updates diagnostic information          Updates diagnostic information
297                    
298          @param name: name of  diagnostic information          :param name: name of  diagnostic information
299          @type name: C{str} in the list "num_iter", "num_level", "num_inner_iter", "time", "set_up_time", "residual_norm", "converged".          :type name: ``str`` in the list "num_iter", "num_level", "num_inner_iter", "time", "set_up_time", "net_time", "residual_norm", "converged".
300          @param vale: new value of the diagnostic information          :param value: new value of the diagnostic information
301          @note: this function is used by a solver to report diagnostics informations.          :note: this function is used by a solver to report diagnostics informations.
302          """          """
303          if name == "num_iter":          if name == "num_iter":
304              self.__num_iter=int(value)              self.__num_iter=int(value)
305              self.__cum_num_iter+=self.__num_iter              self.__cum_num_iter+=self.__num_iter
306          if name == "num_level":          if name == "num_level":
307              self.__num_iter=int(value)              self.__num_level=int(value)
308          if name == "num_inner_iter":          if name == "num_inner_iter":
309              self.__num_inner_iter=int(value)              self.__num_inner_iter=int(value)
310              self.__cum_num_inner_iter+=self.__num_inner_iter              self.__cum_num_inner_iter+=self.__num_inner_iter
# Line 292  class SolverOptions(object): Line 314  class SolverOptions(object):
314          if name == "set_up_time":          if name == "set_up_time":
315              self.__set_up_time=float(value)              self.__set_up_time=float(value)
316              self.__cum_set_up_time+=self.__set_up_time              self.__cum_set_up_time+=self.__set_up_time
317            if name == "net_time":
318                self.__net_time=float(value)
319                self.__cum_net_time+=self.__net_time
320          if name == "residual_norm":          if name == "residual_norm":
321              self.__residual_norm=float(value)              self.__residual_norm=float(value)
322          if name == "converged":          if name == "converged":
323              self.__converged = (value == True)              self.__converged = (value == True)
324      def getDiagnostics(self, name):      def getDiagnostics(self, name):
325          """          """
326          Returns the diagnostic information C{name}          Returns the diagnostic information ``name``. Possible values are:
327                    
328          @param name: name of diagnostic information where      - "num_iter": the number of iteration steps
         - "num_iter": the number of iteration steps  
329          - "cum_num_iter": the cumulative number of iteration steps          - "cum_num_iter": the cumulative number of iteration steps
330          - "num_level": the number of level in multi level solver          - "num_level": the number of level in multi level solver
331          - "num_inner_iter": the number of inner iteration steps          - "num_inner_iter": the number of inner iteration steps
# Line 310  class SolverOptions(object): Line 334  class SolverOptions(object):
334          - "cum_time": cumulative execution time          - "cum_time": cumulative execution time
335          - "set_up_time": time to set up of the solver, typically this includes factorization and reordering          - "set_up_time": time to set up of the solver, typically this includes factorization and reordering
336          - "cum_set_up_time": cumulative time to set up of the solver          - "cum_set_up_time": cumulative time to set up of the solver
337            - "net_time": net execution time, excluding setup time for the solver and execution time for preconditioner
338            - "cum_net_time": cumulative net execution time
339          - "residual_norm": norm of the final residual          - "residual_norm": norm of the final residual
340          - "converged": return self.__converged              - "converged": return self.__converged
341          @type name: C{str} in the list "num_iter", "num_level", "num_inner_iter", "time", "set_up_time", "residual_norm", "converged".      
342          @return: requested value. C{None} is returned if the value is undefined.      
343          @note: If the solver has thrown an exception diagnostic values have an undefined status.          
344            :param name: name of diagnostic information to return
345            :type name: ``str`` in the list above.
346            :return: requested value. ``None`` is returned if the value is undefined.
347            :note: If the solver has thrown an exception diagnostic values have an undefined status.
348          """          """
349          if name == "num_iter": return self.__num_iter          if name == "num_iter": return self.__num_iter
350          elif name == "cum_num_iter": return self.__cum_num_iter          elif name == "cum_num_iter": return self.__cum_num_iter
# Line 325  class SolverOptions(object): Line 355  class SolverOptions(object):
355          elif name == "cum_time": return self.__cum_time          elif name == "cum_time": return self.__cum_time
356          elif name == "set_up_time": return self.__set_up_time          elif name == "set_up_time": return self.__set_up_time
357          elif name == "cum_set_up_time": return self.__cum_set_up_time          elif name == "cum_set_up_time": return self.__cum_set_up_time
358            elif name == "net_time": return self.__net_time
359            elif name == "cum_net_time": return self.__cum_net_time
360          elif name == "residual_norm": return self.__residual_norm          elif name == "residual_norm": return self.__residual_norm
361          elif name == "converged": return self.__converged                elif name == "converged": return self.__converged      
362          else:          else:
363              raise ValueError,"unknown diagnostic item %s"%name              raise ValueError,"unknown diagnostic item %s"%name
364      def hasConverged(self):      def hasConverged(self):
365          """          """
366          Returns C{True} if the last solver call has been finalized successfully.          Returns ``True`` if the last solver call has been finalized successfully.
367          @note: if an exception has been thrown by the solver the status of this flag is undefined.          :note: if an exception has been thrown by the solver the status of this flag is undefined.
368          """          """
369          return self.getDiagnostics("converged")          return self.getDiagnostics("converged")
370      def setCoarsening(self,method=0):      def setCoarsening(self,method=0):
371          """          """
372          Sets the key of the coarsening method to be applied in AMG.          Sets the key of the coarsening method to be applied in AMG.
373    
374          @param method: selects the coarsening method .          :param method: selects the coarsening method .
375          @type method: in {SolverOptions.DEFAULT}, L{SolverOptions.YAIR_SHAPIRA_COARSENING},          :type method: in {SolverOptions.DEFAULT}, `SolverOptions.YAIR_SHAPIRA_COARSENING`,  `SolverOptions.RUGE_STUEBEN_COARSENING`, `SolverOptions.AGGREGATION_COARSENING`
         L{SolverOptions.RUGE_STUEBEN_COARSENING}, L{SolverOptions.AGGREGATION_COARSENING}  
376          """          """
377      if method==None: method=0      if method==None: method=0
378          if not method in [self.DEFAULT, self.YAIR_SHAPIRA_COARSENING, self.RUGE_STUEBEN_COARSENING, self.AGGREGATION_COARSENING]:          if not method in [self.DEFAULT, self.YAIR_SHAPIRA_COARSENING, self.RUGE_STUEBEN_COARSENING, self.AGGREGATION_COARSENING, self.STANDARD_COARSENING,]:
379               raise ValueError,"unknown coarsening method %s"%method               raise ValueError,"unknown coarsening method %s"%method
380          self.__coarsening=method          self.__coarsening=method
381            
# Line 352  class SolverOptions(object): Line 383  class SolverOptions(object):
383          """          """
384          Returns the key of the coarsening algorithm to be applied AMG.          Returns the key of the coarsening algorithm to be applied AMG.
385    
386          @rtype: in the list L{SolverOptions.DEFAULT}, L{SolverOptions.YAIR_SHAPIRA_COARSENING},          :rtype: in the list `SolverOptions.DEFAULT`, `SolverOptions.YAIR_SHAPIRA_COARSENING`, `SolverOptions.RUGE_STUEBEN_COARSENING`, `SolverOptions.AGGREGATION_COARSENING`
         L{SolverOptions.RUGE_STUEBEN_COARSENING}, L{SolverOptions.AGGREGATION_COARSENING}  
387          """          """
388          return self.__coarsening          return self.__coarsening
389                
# Line 361  class SolverOptions(object): Line 391  class SolverOptions(object):
391          """          """
392          Sets the minumum size of the coarsest level matrix in AMG.          Sets the minumum size of the coarsest level matrix in AMG.
393    
394          @param size: minumum size of the coarsest level matrix .          :param size: minumum size of the coarsest level matrix .
395          @type size: positive C{int} or C{None}          :type size: positive ``int`` or ``None``
396          """          """
397          size=int(size)          size=int(size)
398          if size<0:          if size<0:
# Line 374  class SolverOptions(object): Line 404  class SolverOptions(object):
404          """          """
405          Returns the minumum size of the coarsest level matrix in AMG.          Returns the minumum size of the coarsest level matrix in AMG.
406    
407          @rtype: C{int}          :rtype: ``int``
408          """          """
409          return self.__MinCoarseMatrixSize          return self.__MinCoarseMatrixSize
410                
# Line 382  class SolverOptions(object): Line 412  class SolverOptions(object):
412          """          """
413          Sets the preconditioner to be used.          Sets the preconditioner to be used.
414    
415          @param preconditioner: key of the preconditioner to be used.          :param preconditioner: key of the preconditioner to be used.
416          @type preconditioner: in L{SolverOptions.SSOR}, L{SolverOptions.ILU0}, L{SolverOptions.ILUT}, L{SolverOptions.JACOBI},          :type preconditioner: in `SolverOptions.SSOR`, `SolverOptions.ILU0`, `SolverOptions.ILUT`, `SolverOptions.JACOBI`,
417                                      L{SolverOptions.AMG}, L{SolverOptions.REC_ILU}, L{SolverOptions.GAUSS_SEIDEL}, L{SolverOptions.RILU},                                      `SolverOptions.AMG`, `SolverOptions.AMLI`, `SolverOptions.REC_ILU`, `SolverOptions.GAUSS_SEIDEL`, `SolverOptions.GAUSS_SEIDEL_MPI`, `SolverOptions.RILU`,
418                                      L{SolverOptions.NO_PRECONDITIONER}                                      `SolverOptions.NO_PRECONDITIONER`
419          @note: Not all packages support all preconditioner. It can be assumed that a package makes a reasonable choice if it encounters          :note: Not all packages support all preconditioner. It can be assumed that a package makes a reasonable choice if it encounters an unknown preconditioner.
         an unknown preconditioner.  
420          """          """
421      if preconditioner==None: preconditioner=10      if preconditioner==None: preconditioner=10
422          if not preconditioner in [ SolverOptions.SSOR, SolverOptions.ILU0, SolverOptions.ILUT, SolverOptions.JACOBI,          if not preconditioner in [ SolverOptions.SSOR, SolverOptions.ILU0, SolverOptions.ILUT, SolverOptions.JACOBI,
423                                      SolverOptions.AMG, SolverOptions.REC_ILU, SolverOptions.GAUSS_SEIDEL, SolverOptions.RILU,                                      SolverOptions.AMG, SolverOptions.AMLI, SolverOptions.REC_ILU, SolverOptions.GAUSS_SEIDEL, SolverOptions.GAUSS_SEIDEL_MPI, SolverOptions.RILU,
424                                      SolverOptions.NO_PRECONDITIONER] :                                      SolverOptions.NO_PRECONDITIONER] :
425               raise ValueError,"unknown preconditioner %s"%preconditioner               raise ValueError,"unknown preconditioner %s"%preconditioner
426          self.__preconditioner=preconditioner              self.__preconditioner=preconditioner    
# Line 399  class SolverOptions(object): Line 428  class SolverOptions(object):
428          """          """
429          Returns key of the preconditioner to be used.          Returns key of the preconditioner to be used.
430    
431          @rtype: in the list L{SolverOptions.SSOR}, L{SolverOptions.ILU0}, L{SolverOptions.ILUT}, L{SolverOptions.JACOBI},          :rtype: in the list `SolverOptions.SSOR`, `SolverOptions.ILU0`, `SolverOptions.ILUT`, `SolverOptions.JACOBI`,
432                                      L{SolverOptions.AMG}, L{SolverOptions.REC_ILU}, L{SolverOptions.GAUSS_SEIDEL}, L{SolverOptions.RILU},                                      `SolverOptions.AMG`, `SolverOptions.REC_ILU`, `SolverOptions.GAUSS_SEIDEL`, `SolverOptions.GAUSS_SEIDEL_MPI`, `SolverOptions.RILU`,
433                                      L{SolverOptions.NO_PRECONDITIONER}                                      `SolverOptions.NO_PRECONDITIONER`
434          """          """
435          return self.__preconditioner          return self.__preconditioner
436        def setSmoother(self, smoother=28):
437            """
438            Sets the smoother to be used.
439    
440            :param smoother: key of the smoother to be used.
441            :type smoother: in `SolverOptions.JACOBI`, `SolverOptions.GAUSS_SEIDEL`
442            :note: Not all packages support all smoothers. It can be assumed that a package makes a reasonable choice if it encounters an unknown smoother.
443            """
444        if smoother==None: smoother=28
445            if not smoother in [ SolverOptions.JACOBI, SolverOptions.GAUSS_SEIDEL ] :
446                 raise ValueError,"unknown smoother %s"%smoother
447            self.__smoother=smoother    
448        def getSmoother(self):
449            """
450            Returns key of the smoother to be used.
451    
452            :rtype: in the list `SolverOptions.JACOBI`, `SolverOptions.GAUSS_SEIDEL`
453            """
454            return self.__smoother  
455      def setSolverMethod(self, method=0):      def setSolverMethod(self, method=0):
456          """          """
457          Sets the solver method to be used. Use C{method}=C{SolverOptions.DIRECT} to indicate that a direct rather than an iterative          Sets the solver method to be used. Use ``method``=``SolverOptions.DIRECT`` to indicate that a direct rather than an iterative
458          solver should be used and Use C{method}=C{SolverOptions.ITERATIVE} to indicate that an iterative rather than a direct          solver should be used and Use ``method``=``SolverOptions.ITERATIVE`` to indicate that an iterative rather than a direct
459          solver should be used.          solver should be used.
460    
461          @param method: key of the solver method to be used.          :param method: key of the solver method to be used.
462          @type method: in L{SolverOptions.DEFAULT}, L{SolverOptions.DIRECT}, L{SolverOptions.CHOLEVSKY}, L{SolverOptions.PCG},          :type method: in `SolverOptions.DEFAULT`, `SolverOptions.DIRECT`, `SolverOptions.CHOLEVSKY`, `SolverOptions.PCG`,
463                          L{SolverOptions.CR}, L{SolverOptions.CGS}, L{SolverOptions.BICGSTAB}, L{SolverOptions.SSOR},                          `SolverOptions.CR`, `SolverOptions.CGS`, `SolverOptions.BICGSTAB`, `SolverOptions.SSOR`,
464                          L{SolverOptions.GMRES}, L{SolverOptions.PRES20}, L{SolverOptions.LUMPING}, L{SolverOptions.ITERATIVE},                          `SolverOptions.GMRES`, `SolverOptions.PRES20`, `SolverOptions.LUMPING`, `SolverOptions.ITERATIVE`,
465                          L{SolverOptions.AMG}, L{SolverOptions.NONLINEAR_GMRES}, L{SolverOptions.TFQMR}, L{SolverOptions.MINRES},                          `SolverOptions.AMG`, `SolverOptions.NONLINEAR_GMRES`, `SolverOptions.TFQMR`, `SolverOptions.MINRES`,
466                          L{SolverOptions.GAUSS_SEIDEL}                          `SolverOptions.GAUSS_SEIDEL`, `SolverOptions.GAUSS_SEIDEL_MPI`
467          @note: Not all packages support all solvers. It can be assumed that a package makes a reasonable choice if it encounters          :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.
         an unknown solver method.  
468          """          """
469      if method==None: method=0      if method==None: method=0
470          if not method in [ SolverOptions.DEFAULT, SolverOptions.DIRECT, SolverOptions.CHOLEVSKY, SolverOptions.PCG,          if not method in [ SolverOptions.DEFAULT, SolverOptions.DIRECT, SolverOptions.CHOLEVSKY, SolverOptions.PCG,
471                             SolverOptions.CR, SolverOptions.CGS, SolverOptions.BICGSTAB, SolverOptions.SSOR,                             SolverOptions.CR, SolverOptions.CGS, SolverOptions.BICGSTAB, SolverOptions.SSOR,
472                             SolverOptions.GMRES, SolverOptions.PRES20, SolverOptions.LUMPING, SolverOptions.ITERATIVE, SolverOptions.AMG,                             SolverOptions.GMRES, SolverOptions.PRES20, SolverOptions.LUMPING, SolverOptions.ITERATIVE, SolverOptions.AMG,
473                             SolverOptions.NONLINEAR_GMRES, SolverOptions.TFQMR, SolverOptions.MINRES, SolverOptions.GAUSS_SEIDEL]:                             SolverOptions.NONLINEAR_GMRES, SolverOptions.TFQMR, SolverOptions.MINRES, SolverOptions.GAUSS_SEIDEL, SolverOptions.GAUSS_SEIDEL_MPI]:
474               raise ValueError,"unknown solver method %s"%method               raise ValueError,"unknown solver method %s"%method
475          self.__method=method          self.__method=method
476      def getSolverMethod(self):      def getSolverMethod(self):
477          """          """
478          Returns key of the solver method to be used.          Returns key of the solver method to be used.
479    
480          @rtype: in the list L{SolverOptions.DEFAULT}, L{SolverOptions.DIRECT}, L{SolverOptions.CHOLEVSKY}, L{SolverOptions.PCG},          :rtype: in the list `SolverOptions.DEFAULT`, `SolverOptions.DIRECT`, `SolverOptions.CHOLEVSKY`, `SolverOptions.PCG`,
481                          L{SolverOptions.CR}, L{SolverOptions.CGS}, L{SolverOptions.BICGSTAB}, L{SolverOptions.SSOR},                          `SolverOptions.CR`, `SolverOptions.CGS`, `SolverOptions.BICGSTAB`, `SolverOptions.SSOR`,
482                          L{SolverOptions.GMRES}, L{SolverOptions.PRES20}, L{SolverOptions.LUMPING}, L{SolverOptions.ITERATIVE},                          `SolverOptions.GMRES`, `SolverOptions.PRES20`, `SolverOptions.LUMPING`, `SolverOptions.ITERATIVE`,
483                          L{SolverOptions.AMG}, L{SolverOptions.NONLINEAR_GMRES}, L{SolverOptions.TFQMR}, L{SolverOptions.MINRES},                          `SolverOptions.AMG`, `SolverOptions.NONLINEAR_GMRES`, `SolverOptions.TFQMR`, `SolverOptions.MINRES`,
484                          L{SolverOptions.GAUSS_SEIDEL}                          `SolverOptions.GAUSS_SEIDEL`, `SolverOptions.GAUSS_SEIDEL_MPI`
485          """          """
486          return self.__method          return self.__method
487                    
# Line 442  class SolverOptions(object): Line 489  class SolverOptions(object):
489          """          """
490          Sets the solver package to be used as a solver.            Sets the solver package to be used as a solver.  
491    
492          @param package: key of the solver package to be used.          :param package: key of the solver package to be used.
493          @type package: in L{SolverOptions.DEFAULT}, L{SolverOptions.PASO}, L{SolverOptions.SUPER_LU}, L{SolverOptions.PASTIX}, L{SolverOptions.MKL}, L{SolverOptions.UMFPACK}, L{SolverOptions.TRILINOS}          :type package: in `SolverOptions.DEFAULT`, `SolverOptions.PASO`, `SolverOptions.SUPER_LU`, `SolverOptions.PASTIX`, `SolverOptions.MKL`, `SolverOptions.UMFPACK`, `SolverOptions.TRILINOS`
494          @note: Not all packages are support on all implementation. An exception may be thrown on some platforms if a particular is requested.          :note: Not all packages are support on all implementation. An exception may be thrown on some platforms if a particular is requested.
495          """          """
496      if package==None: package=0      if package==None: package=0
497          if not package in [SolverOptions.DEFAULT, SolverOptions.PASO, SolverOptions.SUPER_LU, SolverOptions.PASTIX, SolverOptions.MKL, SolverOptions.UMFPACK, SolverOptions.TRILINOS]:          if not package in [SolverOptions.DEFAULT, SolverOptions.PASO, SolverOptions.SUPER_LU, SolverOptions.PASTIX, SolverOptions.MKL, SolverOptions.UMFPACK, SolverOptions.TRILINOS]:
# Line 454  class SolverOptions(object): Line 501  class SolverOptions(object):
501          """          """
502          Returns the solver package key          Returns the solver package key
503    
504          @rtype: in the list L{SolverOptions.DEFAULT}, L{SolverOptions.PASO}, L{SolverOptions.SUPER_LU}, L{SolverOptions.PASTIX}, L{SolverOptions.MKL}, L{SolverOptions.UMFPACK}, L{SolverOptions.TRILINOS}          :rtype: in the list `SolverOptions.DEFAULT`, `SolverOptions.PASO`, `SolverOptions.SUPER_LU`, `SolverOptions.PASTIX`, `SolverOptions.MKL`, `SolverOptions.UMFPACK`, `SolverOptions.TRILINOS`
505          """          """
506          return self.__package          return self.__package
507      def setReordering(self,ordering=30):      def setReordering(self,ordering=30):
# Line 462  class SolverOptions(object): Line 509  class SolverOptions(object):
509          Sets the key of the reordering method to be applied if supported by the solver. Some direct solvers support reordering          Sets the key of the reordering method to be applied if supported by the solver. Some direct solvers support reordering
510          to optimize compute time and storage use during elimination.          to optimize compute time and storage use during elimination.
511    
512          @param ordering: selects the reordering strategy.          :param ordering: selects the reordering strategy.
513          @type ordering: in L{SolverOptions.NO_REORDERING}, L{SolverOptions.NO_REORDERING},          :type ordering: in `SolverOptions.NO_REORDERING`, `SolverOptions.NO_REORDERING`, `SolverOptions.NO_REORDERING`, `SolverOptions.DEFAULT_REORDERING`
         L{SolverOptions.NO_REORDERING}, L{SolverOptions.DEFAULT_REORDERING}  
514          """          """
515          if not ordering in [self.NO_REORDERING, self.MINIMUM_FILL_IN, self.NESTED_DISSECTION, self.DEFAULT_REORDERING]:          if not ordering in [self.NO_REORDERING, self.MINIMUM_FILL_IN, self.NESTED_DISSECTION, self.DEFAULT_REORDERING]:
516               raise ValueError,"unknown reordering strategy %s"%ordering               raise ValueError,"unknown reordering strategy %s"%ordering
# Line 473  class SolverOptions(object): Line 519  class SolverOptions(object):
519          """          """
520          Returns the key of the reordering method to be applied if supported by the solver.          Returns the key of the reordering method to be applied if supported by the solver.
521    
522          @rtype: in the list L{SolverOptions.NO_REORDERING}, L{SolverOptions.NO_REORDERING},          :rtype: in the list `SolverOptions.NO_REORDERING`, `SolverOptions.NO_REORDERING`,  `SolverOptions.NO_REORDERING`, `SolverOptions.DEFAULT_REORDERING`
         L{SolverOptions.NO_REORDERING}, L{SolverOptions.DEFAULT_REORDERING}  
523          """          """
524          return self.__reordering          return self.__reordering
525      def setRestart(self,restart=None):      def setRestart(self,restart=None):
526          """          """
527          Sets the number of iterations steps after which GMRES is performing a restart.          Sets the number of iterations steps after which GMRES is performing a restart.
528    
529          @param restart: number of iteration steps after which to perform a restart. If equal to C{None} no          :param restart: number of iteration steps after which to perform a restart. If equal to ``None`` no restart is performed.
530                          restart is performed.          :type restart: ``int`` or ``None``
         @type restart: C{int} or C{None}  
531          """          """
532          if restart == None:          if restart == None:
533              self.__restart=restart              self.__restart=restart
# Line 496  class SolverOptions(object): Line 540  class SolverOptions(object):
540      def getRestart(self):      def getRestart(self):
541          """          """
542          Returns the number of iterations steps after which GMRES is performing a restart.          Returns the number of iterations steps after which GMRES is performing a restart.
543          If C{None} is returned no restart is performed.          If ``None`` is returned no restart is performed.
544    
545          @rtype: C{int} or C{None}          :rtype: ``int`` or ``None``
546          """          """
547          if self.__restart < 0:          if self.__restart < 0:
548              return None              return None
# Line 515  class SolverOptions(object): Line 559  class SolverOptions(object):
559          Sets the number of residuals in GMRES to be stored for orthogonalization.  The more residuals are stored          Sets the number of residuals in GMRES to be stored for orthogonalization.  The more residuals are stored
560          the faster GMRES converged but          the faster GMRES converged but
561    
562          @param truncation: truncation          :param truncation: truncation
563          @type truncation: C{int}          :type truncation: ``int``
564          """          """
565          truncation=int(truncation)          truncation=int(truncation)
566          if truncation<1:          if truncation<1:
# Line 526  class SolverOptions(object): Line 570  class SolverOptions(object):
570          """          """
571          Returns the number of residuals in GMRES to be stored for orthogonalization          Returns the number of residuals in GMRES to be stored for orthogonalization
572    
573          @rtype: C{int}          :rtype: ``int``
574          """          """
575          return self.__truncation          return self.__truncation
576      def setInnerIterMax(self,iter_max=10):      def setInnerIterMax(self,iter_max=10):
577          """          """
578          Sets the maximum number of iteration steps for the inner iteration.          Sets the maximum number of iteration steps for the inner iteration.
579    
580          @param iter_max: maximum number of inner iterations          :param iter_max: maximum number of inner iterations
581          @type iter_max: C{int}          :type iter_max: ``int``
582          """          """
583          iter_max=int(iter_max)          iter_max=int(iter_max)
584          if iter_max<1:          if iter_max<1:
# Line 544  class SolverOptions(object): Line 588  class SolverOptions(object):
588          """          """
589          Returns maximum number of inner iteration steps          Returns maximum number of inner iteration steps
590    
591          @rtype: C{int}          :rtype: ``int``
592          """          """
593          return self.__inner_iter_max          return self.__inner_iter_max
594      def setIterMax(self,iter_max=100000):      def setIterMax(self,iter_max=100000):
595          """          """
596          Sets the maximum number of iteration steps          Sets the maximum number of iteration steps
597    
598          @param iter_max: maximum number of iteration steps          :param iter_max: maximum number of iteration steps
599          @type iter_max: C{int}          :type iter_max: ``int``
600          """          """
601          iter_max=int(iter_max)          iter_max=int(iter_max)
602          if iter_max<1:          if iter_max<1:
# Line 562  class SolverOptions(object): Line 606  class SolverOptions(object):
606          """          """
607          Returns maximum number of iteration steps          Returns maximum number of iteration steps
608    
609          @rtype: C{int}          :rtype: ``int``
610          """          """
611          return self.__iter_max          return self.__iter_max
612      def setLevelMax(self,level_max=10):      def setLevelMax(self,level_max=5):
613          """          """
614          Sets the maximum number of coarsening levels to be used in an algebraic multi level solver or preconditioner          Sets the maximum number of coarsening levels to be used in an algebraic multi level solver or preconditioner
615    
616          @param level_max: maximum number of levels          :param level_max: maximum number of levels
617          @type level_max: C{int}          :type level_max: ``int``
618          """          """
619          level_max=int(level_max)          level_max=int(level_max)
620          if level_max<0:          if level_max<0:
# Line 580  class SolverOptions(object): Line 624  class SolverOptions(object):
624          """          """
625          Returns the maximum number of coarsening levels to be used in an algebraic multi level solver or preconditioner          Returns the maximum number of coarsening levels to be used in an algebraic multi level solver or preconditioner
626    
627          @rtype: C{int}          :rtype: ``int``
628          """          """
629          return self.__level_max          return self.__level_max
630      def setCoarseningThreshold(self,theta=0.05):      def setCoarseningThreshold(self,theta=0.25):
631          """          """
632          Sets the threshold for coarsening in the algebraic multi level solver or preconditioner          Sets the threshold for coarsening in the algebraic multi level solver or preconditioner
633    
634          @param theta: threshold for coarsening          :param theta: threshold for coarsening
635          @type theta: positive C{float}          :type theta: positive ``float``
636          """          """
637          theta=float(theta)          theta=float(theta)
638          if theta<0 or theta>1:          if theta<0 or theta>1:
# Line 598  class SolverOptions(object): Line 642  class SolverOptions(object):
642          """          """
643          Returns the threshold for coarsening in the algebraic multi level solver or preconditioner          Returns the threshold for coarsening in the algebraic multi level solver or preconditioner
644    
645          @rtype: C{float}          :rtype: ``float``
646          """          """
647          return self.__coarsening_threshold          return self.__coarsening_threshold
648      def setNumSweeps(self,sweeps=2):      def setNumSweeps(self,sweeps=2):
649          """          """
650          Sets the number of sweeps in a Jacobi or Gauss-Seidel/SOR preconditioner.          Sets the number of sweeps in a Jacobi or Gauss-Seidel/SOR preconditioner.
651    
652          @param sweeps: number of sweeps          :param sweeps: number of sweeps
653          @type theta: positive C{int}          :type sweeps: positive ``int``
654          """          """
655          sweeps=int(sweeps)          sweeps=int(sweeps)
656          if sweeps<1:          if sweeps<1:
# Line 616  class SolverOptions(object): Line 660  class SolverOptions(object):
660          """          """
661          Returns the number of sweeps in a Jacobi or Gauss-Seidel/SOR preconditioner.          Returns the number of sweeps in a Jacobi or Gauss-Seidel/SOR preconditioner.
662    
663          @rtype: C{int}          :rtype: ``int``
664          """          """
665          return self.__sweeps          return self.__sweeps
666      def setNumPreSweeps(self,sweeps=2):      def setNumPreSweeps(self,sweeps=2):
667          """          """
668          Sets the number of sweeps in the pre-smoothing step of a multi level solver or preconditioner          Sets the number of sweeps in the pre-smoothing step of a multi level solver or preconditioner
669    
670          @param sweeps: number of sweeps          :param sweeps: number of sweeps
671          @type theta: positive C{int}          :type sweeps: positive ``int``
672          """          """
673          sweeps=int(sweeps)          sweeps=int(sweeps)
674          if sweeps<1:          if sweeps<1:
# Line 634  class SolverOptions(object): Line 678  class SolverOptions(object):
678          """          """
679          Returns he number of sweeps in the pre-smoothing step of a multi level solver or preconditioner          Returns he number of sweeps in the pre-smoothing step of a multi level solver or preconditioner
680    
681          @rtype: C{int}          :rtype: ``int``
682          """          """
683          return self.__pre_sweeps          return self.__pre_sweeps
684      def setNumPostSweeps(self,sweeps=2):      def setNumPostSweeps(self,sweeps=2):
685          """          """
686          Sets the number of sweeps in the post-smoothing step of a multi level solver or preconditioner          Sets the number of sweeps in the post-smoothing step of a multi level solver or preconditioner
687    
688          @param sweeps: number of sweeps          :param sweeps: number of sweeps
689          @type theta: positive C{int}          :type sweeps: positive ``int``
690          """          """
691          sweeps=int(sweeps)          sweeps=int(sweeps)
692          if sweeps<1:          if sweeps<1:
# Line 652  class SolverOptions(object): Line 696  class SolverOptions(object):
696          """          """
697          Returns he number of sweeps in the post-smoothing step of a multi level solver or preconditioner          Returns he number of sweeps in the post-smoothing step of a multi level solver or preconditioner
698    
699          @rtype: C{int}          :rtype: ``int``
700          """          """
701          return self.__post_sweeps          return self.__post_sweeps
702    
# Line 660  class SolverOptions(object): Line 704  class SolverOptions(object):
704          """          """
705          Sets the relative tolerance for the solver          Sets the relative tolerance for the solver
706    
707          @param rtol: relative tolerance          :param rtol: relative tolerance
708          @type rtol: non-negative C{float}          :type rtol: non-negative ``float``
709          """          """
710          rtol=float(rtol)          rtol=float(rtol)
711          if rtol<0 or rtol>1:          if rtol<0 or rtol>1:
# Line 671  class SolverOptions(object): Line 715  class SolverOptions(object):
715          """          """
716          Returns the relative tolerance for the solver          Returns the relative tolerance for the solver
717    
718          @rtype: C{float}          :rtype: ``float``
719          """          """
720          return self.__tolerance          return self.__tolerance
721      def setAbsoluteTolerance(self,atol=0.):      def setAbsoluteTolerance(self,atol=0.):
722          """          """
723          Sets the absolute tolerance for the solver          Sets the absolute tolerance for the solver
724    
725          @param atol:  absolute tolerance          :param atol:  absolute tolerance
726          @type atol: non-negative C{float}          :type atol: non-negative ``float``
727          """          """
728          atol=float(atol)          atol=float(atol)
729          if atol<0:          if atol<0:
# Line 689  class SolverOptions(object): Line 733  class SolverOptions(object):
733          """          """
734          Returns the absolute tolerance for the solver          Returns the absolute tolerance for the solver
735    
736          @rtype: C{float}          :rtype: ``float``
737          """          """
738          return self.__absolute_tolerance          return self.__absolute_tolerance
739    
740      def setInnerTolerance(self,rtol=0.9):      def setInnerTolerance(self,rtol=0.9):
741          """          """
742           Sets the relative tolerance for an inner iteration scheme for instance           Sets the relative tolerance for an inner iteration scheme for instance on the coarsest level in a multi-level scheme.
         on the coarsest level in a multi-level scheme.  
743    
744          @param rtol: inner relative tolerance          :param rtol: inner relative tolerance
745          @type rtol: positive C{float}          :type rtol: positive ``float``
746          """          """
747          rtol=float(rtol)          rtol=float(rtol)
748          if rtol<=0 or rtol>1:          if rtol<=0 or rtol>1:
# Line 709  class SolverOptions(object): Line 752  class SolverOptions(object):
752          """          """
753          Returns the relative tolerance for an inner iteration scheme          Returns the relative tolerance for an inner iteration scheme
754    
755          @rtype: C{float}          :rtype: ``float``
756          """          """
757          return self.__inner_tolerance          return self.__inner_tolerance
758      def setDropTolerance(self,drop_tol=0.01):      def setDropTolerance(self,drop_tol=0.01):
759          """          """
760          Sets the relative drop tolerance in ILUT          Sets the relative drop tolerance in ILUT
761    
762          @param drop_tol: drop tolerance          :param drop_tol: drop tolerance
763          @type drop_tol: positive C{float}          :type drop_tol: positive ``float``
764          """          """
765          drop_tol=float(drop_tol)          drop_tol=float(drop_tol)
766          if drop_tol<=0 or drop_tol>1:          if drop_tol<=0 or drop_tol>1:
# Line 727  class SolverOptions(object): Line 770  class SolverOptions(object):
770          """          """
771          Returns the relative drop tolerance in ILUT          Returns the relative drop tolerance in ILUT
772    
773          @rtype: C{float}          :rtype: ``float``
774          """          """
775          return self.__drop_tolerance          return self.__drop_tolerance
776      def setDropStorage(self,storage=2.):      def setDropStorage(self,storage=2.):
777          """          """
778          Sets the maximum allowed increase in storage for ILUT. C{storage}=2 would mean that          Sets the maximum allowed increase in storage for ILUT. ``storage`` =2 would mean that
779          a doubling of the storage needed for the coefficient matrix is allowed in the ILUT factorization.          a doubling of the storage needed for the coefficient matrix is allowed in the ILUT factorization.
780    
781          @param storage: allowed storage increase          :param storage: allowed storage increase
782          @type storage: C{float}          :type storage: ``float``
783          """          """
784          storage=float(storage)          storage=float(storage)
785          if storage<1:          if storage<1:
# Line 747  class SolverOptions(object): Line 790  class SolverOptions(object):
790          """          """
791          Returns the maximum allowed increase in storage for ILUT          Returns the maximum allowed increase in storage for ILUT
792    
793          @rtype: C{float}          :rtype: ``float``
794          """          """
795          return self.__drop_storage          return self.__drop_storage
796      def setRelaxationFactor(self,factor=0.3):      def setRelaxationFactor(self,factor=0.3):
797          """          """
798          Sets the relaxation factor used to add dropped elements in RILU to the main diagonal.          Sets the relaxation factor used to add dropped elements in RILU to the main diagonal.
799    
800          @param factor: relaxation factor          :param factor: relaxation factor
801          @type factor: C{float}          :type factor: ``float``
802          @note: RILU with a relaxation factor 0 is identical to ILU0          :note: RILU with a relaxation factor 0 is identical to ILU0
803          """          """
804          factor=float(factor)          factor=float(factor)
805          if factor<0:          if factor<0:
# Line 767  class SolverOptions(object): Line 810  class SolverOptions(object):
810          """          """
811          Returns the relaxation factor used to add dropped elements in RILU to the main diagonal.          Returns the relaxation factor used to add dropped elements in RILU to the main diagonal.
812    
813          @rtype: C{float}          :rtype: ``float``
814          """          """
815          return self.__relaxation          return self.__relaxation
816      def isSymmetric(self):      def isSymmetric(self):
817          """          """
818          Checks if symmetry of the  coefficient matrix is indicated.          Checks if symmetry of the  coefficient matrix is indicated.
819    
820          @return: True if a symmetric PDE is indicated, False otherwise          :return: True if a symmetric PDE is indicated, False otherwise
821          @rtype: C{bool}          :rtype: ``bool``
822          """          """
823          return self.__symmetric          return self.__symmetric
824      def setSymmetryOn(self):      def setSymmetryOn(self):
# Line 790  class SolverOptions(object): Line 833  class SolverOptions(object):
833          self.__symmetric=False          self.__symmetric=False
834      def setSymmetry(self,flag=False):      def setSymmetry(self,flag=False):
835          """          """
836          Sets the symmetry flag for the coefficient matrix to C{flag}.          Sets the symmetry flag for the coefficient matrix to ``flag``.
837    
838          @param flag: If True, the symmetry flag is set otherwise reset.          :param flag: If True, the symmetry flag is set otherwise reset.
839          @type flag: C{bool}          :type flag: ``bool``
840          """          """
841          if flag:          if flag:
842              self.setSymmetryOn()              self.setSymmetryOn()
# Line 801  class SolverOptions(object): Line 844  class SolverOptions(object):
844              self.setSymmetryOff()              self.setSymmetryOff()
845      def isVerbose(self):      def isVerbose(self):
846          """          """
847          Returns C{True} if the solver is expected to be verbose.          Returns ``True`` if the solver is expected to be verbose.
848    
849          @return: True if verbosity of switched on.          :return: True if verbosity of switched on.
850          @rtype: C{bool}          :rtype: ``bool``
851          """          """
852          return self.__verbose          return self.__verbose
853    
# Line 820  class SolverOptions(object): Line 863  class SolverOptions(object):
863          self.__verbose=False          self.__verbose=False
864      def setVerbosity(self,verbose=False):      def setVerbosity(self,verbose=False):
865          """          """
866          Sets the verbosity flag for the solver to C{flag}.          Sets the verbosity flag for the solver to ``flag``.
867    
868          @param flag: If C{True}, the verbosity of the solver is switched on.          :param verbose: If ``True``, the verbosity of the solver is switched on.
869          @type flag: C{bool}          :type verbose: ``bool``
870          """          """
871          if verbose:          if verbose:
872              self.setVerbosityOn()              self.setVerbosityOn()
# Line 832  class SolverOptions(object): Line 875  class SolverOptions(object):
875                    
876      def adaptInnerTolerance(self):      def adaptInnerTolerance(self):
877          """          """
878          Returns C{True} if the tolerance of the inner solver is selected automatically.          Returns ``True`` if the tolerance of the inner solver is selected automatically.
879          Otherwise the inner tolerance set by L{setInnerTolerance} is used.          Otherwise the inner tolerance set by `setInnerTolerance` is used.
880    
881          @return: C{True} if inner tolerance adaption is chosen.          :return: ``True`` if inner tolerance adaption is chosen.
882          @rtype: C{bool}          :rtype: ``bool``
883          """          """
884          return self.__adapt_inner_tolerance          return self.__adapt_inner_tolerance
885    
# Line 854  class SolverOptions(object): Line 897  class SolverOptions(object):
897          """          """
898          Sets a flag to indicate automatic selection of the inner tolerance.          Sets a flag to indicate automatic selection of the inner tolerance.
899    
900          @param adapt: If C{True}, the inner tolerance is selected automatically.          :param adapt: If ``True``, the inner tolerance is selected automatically.
901          @type adapt: C{bool}          :type adapt: ``bool``
902          """          """
903          if adapt:          if adapt:
904              self.setInnerToleranceAdaptionOn()              self.setInnerToleranceAdaptionOn()
# Line 864  class SolverOptions(object): Line 907  class SolverOptions(object):
907    
908      def acceptConvergenceFailure(self):      def acceptConvergenceFailure(self):
909          """          """
910          Returns C{True} if a failure to meet the stopping criteria within the          Returns ``True`` if a failure to meet the stopping criteria within the
911          given number of iteration steps is not raising in exception. This is useful          given number of iteration steps is not raising in exception. This is useful
912          if a solver is used in a non-linear context where the non-linear solver can          if a solver is used in a non-linear context where the non-linear solver can
913          continue even if the returned the solution does not necessarily meet the          continue even if the returned the solution does not necessarily meet the
914          stopping criteria. One can use the L{hasConverged} method to check if the          stopping criteria. One can use the `hasConverged` method to check if the
915          last call to the solver was successful.          last call to the solver was successful.
916    
917          @return: C{True} if a failure to achieve convergence is accepted.          :return: ``True`` if a failure to achieve convergence is accepted.
918          @rtype: C{bool}          :rtype: ``bool``
919          """          """
920          return self.__accept_convergence_failure          return self.__accept_convergence_failure
921    
# Line 890  class SolverOptions(object): Line 933  class SolverOptions(object):
933          """          """
934          Sets a flag to indicate the acceptance of a failure of convergence.          Sets a flag to indicate the acceptance of a failure of convergence.
935    
936          @param accept: If C{True}, any failure to achieve convergence is accepted.          :param accept: If ``True``, any failure to achieve convergence is accepted.
937          @type accept: C{bool}          :type accept: ``bool``
938          """          """
939          if accept:          if accept:
940              self.setAcceptanceConvergenceFailureOn()              self.setAcceptanceConvergenceFailureOn()
# Line 927  class PDECoef(object): Line 970  class PDECoef(object):
970      """      """
971      A class for describing a PDE coefficient.      A class for describing a PDE coefficient.
972    
973      @cvar INTERIOR: indicator that coefficient is defined on the interior of      :cvar INTERIOR: indicator that coefficient is defined on the interior of
974                      the PDE domain                      the PDE domain
975      @cvar BOUNDARY: indicator that coefficient is defined on the boundary of      :cvar BOUNDARY: indicator that coefficient is defined on the boundary of
976                      the PDE domain                      the PDE domain
977      @cvar CONTACT: indicator that coefficient is defined on the contact region      :cvar CONTACT: indicator that coefficient is defined on the contact region
978                     within the PDE domain                     within the PDE domain
979      @cvar INTERIOR_REDUCED: indicator that coefficient is defined on the      :cvar INTERIOR_REDUCED: indicator that coefficient is defined on the
980                              interior of the PDE domain using a reduced                              interior of the PDE domain using a reduced
981                              integration order                              integration order
982      @cvar BOUNDARY_REDUCED: indicator that coefficient is defined on the      :cvar BOUNDARY_REDUCED: indicator that coefficient is defined on the
983                              boundary of the PDE domain using a reduced                              boundary of the PDE domain using a reduced
984                              integration order                              integration order
985      @cvar CONTACT_REDUCED: indicator that coefficient is defined on the contact      :cvar CONTACT_REDUCED: indicator that coefficient is defined on the contact
986                             region within the PDE domain using a reduced                             region within the PDE domain using a reduced
987                             integration order                             integration order
988      @cvar SOLUTION: indicator that coefficient is defined through a solution of      :cvar SOLUTION: indicator that coefficient is defined through a solution of
989                      the PDE                      the PDE
990      @cvar REDUCED: indicator that coefficient is defined through a reduced      :cvar REDUCED: indicator that coefficient is defined through a reduced
991                     solution of the PDE                     solution of the PDE
992      @cvar BY_EQUATION: indicator that the dimension of the coefficient shape is      :cvar BY_EQUATION: indicator that the dimension of the coefficient shape is
993                         defined by the number of PDE equations                         defined by the number of PDE equations
994      @cvar BY_SOLUTION: indicator that the dimension of the coefficient shape is      :cvar BY_SOLUTION: indicator that the dimension of the coefficient shape is
995                         defined by the number of PDE solutions                         defined by the number of PDE solutions
996      @cvar BY_DIM: indicator that the dimension of the coefficient shape is      :cvar BY_DIM: indicator that the dimension of the coefficient shape is
997                    defined by the spatial dimension                    defined by the spatial dimension
998      @cvar OPERATOR: indicator that the the coefficient alters the operator of      :cvar OPERATOR: indicator that the the coefficient alters the operator of
999                      the PDE                      the PDE
1000      @cvar RIGHTHANDSIDE: indicator that the the coefficient alters the right      :cvar RIGHTHANDSIDE: indicator that the the coefficient alters the right
1001                           hand side of the PDE                           hand side of the PDE
1002      @cvar BOTH: indicator that the the coefficient alters the operator as well      :cvar BOTH: indicator that the the coefficient alters the operator as well
1003                  as the right hand side of the PDE                  as the right hand side of the PDE
1004    
1005      """      """
# Line 979  class PDECoef(object): Line 1022  class PDECoef(object):
1022         """         """
1023         Initialises a PDE coefficient type.         Initialises a PDE coefficient type.
1024    
1025         @param where: describes where the coefficient lives         :param where: describes where the coefficient lives
1026         @type where: one of L{INTERIOR}, L{BOUNDARY}, L{CONTACT}, L{SOLUTION},         :type where: one of `INTERIOR`, `BOUNDARY`, `CONTACT`, `SOLUTION`,
1027                      L{REDUCED}, L{INTERIOR_REDUCED}, L{BOUNDARY_REDUCED},                      `REDUCED`, `INTERIOR_REDUCED`, `BOUNDARY_REDUCED`,
1028                      L{CONTACT_REDUCED}                      `CONTACT_REDUCED`
1029         @param pattern: describes the shape of the coefficient and how the shape         :param pattern: describes the shape of the coefficient and how the shape
1030                         is built for a given spatial dimension and numbers of                         is built for a given spatial dimension and numbers of
1031                         equations and solutions in then PDE. For instance,                         equations and solutions in then PDE. For instance,
1032                         (L{BY_EQUATION},L{BY_SOLUTION},L{BY_DIM}) describes a                         (`BY_EQUATION`,`BY_SOLUTION`,`BY_DIM`) describes a
1033                         rank 3 coefficient which is instantiated as shape (3,2,2)                         rank 3 coefficient which is instantiated as shape (3,2,2)
1034                         in case of three equations and two solution components                         in case of three equations and two solution components
1035                         on a 2-dimensional domain. In the case of single equation                         on a 2-dimensional domain. In the case of single equation
1036                         and a single solution component the shape components                         and a single solution component the shape components
1037                         marked by L{BY_EQUATION} or L{BY_SOLUTION} are dropped.                         marked by `BY_EQUATION` or `BY_SOLUTION` are dropped.
1038                         In this case the example would be read as (2,).                         In this case the example would be read as (2,).
1039         @type pattern: C{tuple} of L{BY_EQUATION}, L{BY_SOLUTION}, L{BY_DIM}         :type pattern: ``tuple`` of `BY_EQUATION`, `BY_SOLUTION`, `BY_DIM`
1040         @param altering: indicates what part of the PDE is altered if the         :param altering: indicates what part of the PDE is altered if the
1041                          coefficient is altered                          coefficient is altered
1042         @type altering: one of L{OPERATOR}, L{RIGHTHANDSIDE}, L{BOTH}         :type altering: one of `OPERATOR`, `RIGHTHANDSIDE`, `BOTH`
1043         """         """
1044         super(PDECoef, self).__init__()         super(PDECoef, self).__init__()
1045         self.what=where         self.what=where
# Line 1012  class PDECoef(object): Line 1055  class PDECoef(object):
1055    
1056      def getFunctionSpace(self,domain,reducedEquationOrder=False,reducedSolutionOrder=False):      def getFunctionSpace(self,domain,reducedEquationOrder=False,reducedSolutionOrder=False):
1057         """         """
1058         Returns the L{FunctionSpace<escript.FunctionSpace>} of the coefficient.         Returns the `FunctionSpace` of the coefficient.
1059    
1060         @param domain: domain on which the PDE uses the coefficient         :param domain: domain on which the PDE uses the coefficient
1061         @type domain: L{Domain<escript.Domain>}         :type domain: `Domain`
1062         @param reducedEquationOrder: True to indicate that reduced order is used         :param reducedEquationOrder: True to indicate that reduced order is used
1063                                      to represent the equation                                      to represent the equation
1064         @type reducedEquationOrder: C{bool}         :type reducedEquationOrder: ``bool``
1065         @param reducedSolutionOrder: True to indicate that reduced order is used         :param reducedSolutionOrder: True to indicate that reduced order is used
1066                                      to represent the solution                                      to represent the solution
1067         @type reducedSolutionOrder: C{bool}         :type reducedSolutionOrder: ``bool``
1068         @return: L{FunctionSpace<escript.FunctionSpace>} of the coefficient         :return: `FunctionSpace` of the coefficient
1069         @rtype: L{FunctionSpace<escript.FunctionSpace>}         :rtype: `FunctionSpace`
1070         """         """
1071         if self.what==self.INTERIOR:         if self.what==self.INTERIOR:
1072              return escript.Function(domain)              return escript.Function(domain)
# Line 1049  class PDECoef(object): Line 1092  class PDECoef(object):
1092         """         """
1093         Returns the value of the coefficient.         Returns the value of the coefficient.
1094    
1095         @return: value of the coefficient         :return: value of the coefficient
1096         @rtype: L{Data<escript.Data>}         :rtype: `Data`
1097         """         """
1098         return self.value         return self.value
1099    
# Line 1058  class PDECoef(object): Line 1101  class PDECoef(object):
1101         """         """
1102         Sets the value of the coefficient to a new value.         Sets the value of the coefficient to a new value.
1103    
1104         @param domain: domain on which the PDE uses the coefficient         :param domain: domain on which the PDE uses the coefficient
1105         @type domain: L{Domain<escript.Domain>}         :type domain: `Domain`
1106         @param numEquations: number of equations of the PDE         :param numEquations: number of equations of the PDE
1107         @type numEquations: C{int}         :type numEquations: ``int``
1108         @param numSolutions: number of components of the PDE solution         :param numSolutions: number of components of the PDE solution
1109         @type numSolutions: C{int}         :type numSolutions: ``int``
1110         @param reducedEquationOrder: True to indicate that reduced order is used         :param reducedEquationOrder: True to indicate that reduced order is used
1111                                      to represent the equation                                      to represent the equation
1112         @type reducedEquationOrder: C{bool}         :type reducedEquationOrder: ``bool``
1113         @param reducedSolutionOrder: True to indicate that reduced order is used         :param reducedSolutionOrder: True to indicate that reduced order is used
1114                                      to represent the solution                                      to represent the solution
1115         @type reducedSolutionOrder: C{bool}         :type reducedSolutionOrder: ``bool``
1116         @param newValue: number of components of the PDE solution         :param newValue: number of components of the PDE solution
1117         @type newValue: any object that can be converted into a         :type newValue: any object that can be converted into a
1118                         L{Data<escript.Data>} object with the appropriate shape                         `Data` object with the appropriate shape
1119                         and L{FunctionSpace<escript.FunctionSpace>}                         and `FunctionSpace`
1120         @raise IllegalCoefficientValue: if the shape of the assigned value does         :raise IllegalCoefficientValue: if the shape of the assigned value does
1121                                         not match the shape of the coefficient                                         not match the shape of the coefficient
1122         @raise IllegalCoefficientFunctionSpace: if unable to interpolate value         :raise IllegalCoefficientFunctionSpace: if unable to interpolate value
1123                                                 to appropriate function space                                                 to appropriate function space
1124         """         """
1125         if newValue==None:         if newValue==None:
# Line 1099  class PDECoef(object): Line 1142  class PDECoef(object):
1142          """          """
1143          Checks if the coefficient alters the operator of the PDE.          Checks if the coefficient alters the operator of the PDE.
1144    
1145          @return: True if the operator of the PDE is changed when the          :return: True if the operator of the PDE is changed when the
1146                   coefficient is changed                   coefficient is changed
1147          @rtype: C{bool}          :rtype: ``bool``
1148          """          """
1149          if self.altering==self.OPERATOR or self.altering==self.BOTH:          if self.altering==self.OPERATOR or self.altering==self.BOTH:
1150              return not None              return not None
# Line 1112  class PDECoef(object): Line 1155  class PDECoef(object):
1155          """          """
1156          Checks if the coefficient alters the right hand side of the PDE.          Checks if the coefficient alters the right hand side of the PDE.
1157    
1158          @rtype: C{bool}          :rtype: ``bool``
1159          @return: True if the right hand side of the PDE is changed when the          :return: True if the right hand side of the PDE is changed when the
1160                   coefficient is changed, C{None} otherwise.                   coefficient is changed, ``None`` otherwise.
1161          """          """
1162          if self.altering==self.RIGHTHANDSIDE or self.altering==self.BOTH:          if self.altering==self.RIGHTHANDSIDE or self.altering==self.BOTH:
1163              return not None              return not None
# Line 1126  class PDECoef(object): Line 1169  class PDECoef(object):
1169         Tries to estimate the number of equations and number of solutions if         Tries to estimate the number of equations and number of solutions if
1170         the coefficient has the given shape.         the coefficient has the given shape.
1171    
1172         @param domain: domain on which the PDE uses the coefficient         :param domain: domain on which the PDE uses the coefficient
1173         @type domain: L{Domain<escript.Domain>}         :type domain: `Domain`
1174         @param shape: suggested shape of the coefficient         :param shape: suggested shape of the coefficient
1175         @type shape: C{tuple} of C{int} values         :type shape: ``tuple`` of ``int`` values
1176         @return: the number of equations and number of solutions of the PDE if         :return: the number of equations and number of solutions of the PDE if
1177                  the coefficient has given shape. If no appropriate numbers                  the coefficient has given shape. If no appropriate numbers
1178                  could be identified, C{None} is returned                  could be identified, ``None`` is returned
1179         @rtype: C{tuple} of two C{int} values or C{None}         :rtype: ``tuple`` of two ``int`` values or ``None``
1180         """         """
1181         dim=domain.getDim()         dim=domain.getDim()
1182         if len(shape)>0:         if len(shape)>0:
# Line 1174  class PDECoef(object): Line 1217  class PDECoef(object):
1217         Checks if the coefficient allows to estimate the number of solution         Checks if the coefficient allows to estimate the number of solution
1218         components.         components.
1219    
1220         @return: True if the coefficient allows an estimate of the number of         :return: True if the coefficient allows an estimate of the number of
1221                  solution components, False otherwise                  solution components, False otherwise
1222         @rtype: C{bool}         :rtype: ``bool``
1223         """         """
1224         for i in self.pattern:         for i in self.pattern:
1225               if i==self.BY_SOLUTION: return True               if i==self.BY_SOLUTION: return True
# Line 1186  class PDECoef(object): Line 1229  class PDECoef(object):
1229         """         """
1230         Checks if the coefficient allows to estimate the number of equations.         Checks if the coefficient allows to estimate the number of equations.
1231    
1232         @return: True if the coefficient allows an estimate of the number of         :return: True if the coefficient allows an estimate of the number of
1233                  equations, False otherwise                  equations, False otherwise
1234         @rtype: C{bool}         :rtype: ``bool``
1235         """         """
1236         for i in self.pattern:         for i in self.pattern:
1237               if i==self.BY_EQUATION: return True               if i==self.BY_EQUATION: return True
# Line 1199  class PDECoef(object): Line 1242  class PDECoef(object):
1242        Compares two tuples of possible number of equations and number of        Compares two tuples of possible number of equations and number of
1243        solutions.        solutions.
1244    
1245        @param t1: the first tuple        :param t1: the first tuple
1246        @param t2: the second tuple        :param t2: the second tuple
1247        @return: 0, 1, or -1        :return: 0, 1, or -1
1248        """        """
1249    
1250        dif=t1[0]+t1[1]-(t2[0]+t2[1])        dif=t1[0]+t1[1]-(t2[0]+t2[1])
# Line 1213  class PDECoef(object): Line 1256  class PDECoef(object):
1256         """         """
1257         Builds the required shape of the coefficient.         Builds the required shape of the coefficient.
1258    
1259         @param domain: domain on which the PDE uses the coefficient         :param domain: domain on which the PDE uses the coefficient
1260         @type domain: L{Domain<escript.Domain>}         :type domain: `Domain`
1261         @param numEquations: number of equations of the PDE         :param numEquations: number of equations of the PDE
1262         @type numEquations: C{int}         :type numEquations: ``int``
1263         @param numSolutions: number of components of the PDE solution         :param numSolutions: number of components of the PDE solution
1264         @type numSolutions: C{int}         :type numSolutions: ``int``
1265         @return: shape of the coefficient         :return: shape of the coefficient
1266         @rtype: C{tuple} of C{int} values         :rtype: ``tuple`` of ``int`` values
1267         """         """
1268         dim=domain.getDim()         dim=domain.getDim()
1269         s=()         s=()
# Line 1238  class PDECoef(object): Line 1281  class PDECoef(object):
1281  class LinearProblem(object):  class LinearProblem(object):
1282     """     """
1283     This is the base class to define a general linear PDE-type problem for     This is the base class to define a general linear PDE-type problem for
1284     for an unknown function M{u} on a given domain defined through a     for an unknown function *u* on a given domain defined through a
1285     L{Domain<escript.Domain>} object. The problem can be given as a single     `Domain` object. The problem can be given as a single
1286     equation or as a system of equations.     equation or as a system of equations.
1287    
1288     The class assumes that some sort of assembling process is required to form     The class assumes that some sort of assembling process is required to form
1289     a problem of the form     a problem of the form
1290    
1291     M{L u=f}     *L u=f*
1292    
1293     where M{L} is an operator and M{f} is the right hand side. This operator     where *L* is an operator and *f* is the right hand side. This operator
1294     problem will be solved to get the unknown M{u}.     problem will be solved to get the unknown *u*.
1295    
1296     """     """
1297     def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):     def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):
1298       """       """
1299       Initializes a linear problem.       Initializes a linear problem.
1300    
1301       @param domain: domain of the PDE       :param domain: domain of the PDE
1302       @type domain: L{Domain<escript.Domain>}       :type domain: `Domain`
1303       @param numEquations: number of equations. If C{None} the number of       :param numEquations: number of equations. If ``None`` the number of
1304                            equations is extracted from the coefficients.                            equations is extracted from the coefficients.
1305       @param numSolutions: number of solution components. If C{None} the number       :param numSolutions: number of solution components. If ``None`` the number
1306                            of solution components is extracted from the                            of solution components is extracted from the
1307                            coefficients.                            coefficients.
1308       @param debug: if True debug information is printed       :param debug: if True debug information is printed
1309    
1310       """       """
1311       super(LinearProblem, self).__init__()       super(LinearProblem, self).__init__()
# Line 1292  class LinearProblem(object): Line 1335  class LinearProblem(object):
1335       """       """
1336       Returns a string representation of the PDE.       Returns a string representation of the PDE.
1337    
1338       @return: a simple representation of the PDE       :return: a simple representation of the PDE
1339       @rtype: C{str}       :rtype: ``str``
1340       """       """
1341       return "<LinearProblem %d>"%id(self)       return "<LinearProblem %d>"%id(self)
1342     # ==========================================================================     # ==========================================================================
# Line 1313  class LinearProblem(object): Line 1356  class LinearProblem(object):
1356    
1357     def setDebug(self, flag):     def setDebug(self, flag):
1358       """       """
1359       Switches debug output on if C{flag} is True otherwise it is switched off.       Switches debug output on if ``flag`` is True otherwise it is switched off.
1360    
1361       @param flag: desired debug status       :param flag: desired debug status
1362       @type flag: C{bool}       :type flag: ``bool``
1363       """       """
1364       if flag:       if flag:
1365           self.setDebugOn()           self.setDebugOn()
# Line 1327  class LinearProblem(object): Line 1370  class LinearProblem(object):
1370       """       """
1371       Prints the text message if debug mode is switched on.       Prints the text message if debug mode is switched on.
1372    
1373       @param text: message to be printed       :param text: message to be printed
1374       @type text: C{string}       :type text: ``string``
1375       """       """
1376       if self.__debug: print "%s: %s"%(str(self),text)       if self.__debug: print "%s: %s"%(str(self),text)
1377    
# Line 1339  class LinearProblem(object): Line 1382  class LinearProblem(object):
1382         """         """
1383         Introduces new coefficients into the problem.         Introduces new coefficients into the problem.
1384    
1385         Use::         Use:
1386    
1387         p.introduceCoefficients(A=PDECoef(...), B=PDECoef(...))         p.introduceCoefficients(A=PDECoef(...), B=PDECoef(...))
1388    
1389         to introduce the coefficients M{A} ans M{B}.         to introduce the coefficients *A* and *B*.
1390         """         """
1391         for name, type in coeff.items():         for name, type in coeff.items():
1392             if not isinstance(type,PDECoef):             if not isinstance(type,PDECoef):
# Line 1356  class LinearProblem(object): Line 1399  class LinearProblem(object):
1399       """       """
1400       Returns the domain of the PDE.       Returns the domain of the PDE.
1401    
1402       @return: the domain of the PDE       :return: the domain of the PDE
1403       @rtype: L{Domain<escript.Domain>}       :rtype: `Domain`
1404       """       """
1405       return self.__domain       return self.__domain
1406     def getDomainStatus(self):     def getDomainStatus(self):
# Line 1373  class LinearProblem(object): Line 1416  class LinearProblem(object):
1416       return self.__system_status       return self.__system_status
1417     def setSystemStatus(self,status=None):     def setSystemStatus(self,status=None):
1418       """       """
1419       Sets the system status to C{status} if C{status} is not present the       Sets the system status to ``status`` if ``status`` is not present the
1420       current status of the domain is used.       current status of the domain is used.
1421       """       """
1422       if status == None:       if status == None:
# Line 1385  class LinearProblem(object): Line 1428  class LinearProblem(object):
1428       """       """
1429       Returns the spatial dimension of the PDE.       Returns the spatial dimension of the PDE.
1430    
1431       @return: the spatial dimension of the PDE domain       :return: the spatial dimension of the PDE domain
1432       @rtype: C{int}       :rtype: ``int``
1433       """       """
1434       return self.getDomain().getDim()       return self.getDomain().getDim()
1435    
# Line 1394  class LinearProblem(object): Line 1437  class LinearProblem(object):
1437       """       """
1438       Returns the number of equations.       Returns the number of equations.
1439    
1440       @return: the number of equations       :return: the number of equations
1441       @rtype: C{int}       :rtype: ``int``
1442       @raise UndefinedPDEError: if the number of equations is not specified yet       :raise UndefinedPDEError: if the number of equations is not specified yet
1443       """       """
1444       if self.__numEquations==None:       if self.__numEquations==None:
1445           if self.__numSolutions==None:           if self.__numSolutions==None:
# Line 1409  class LinearProblem(object): Line 1452  class LinearProblem(object):
1452       """       """
1453       Returns the number of unknowns.       Returns the number of unknowns.
1454    
1455       @return: the number of unknowns       :return: the number of unknowns
1456       @rtype: C{int}       :rtype: ``int``
1457       @raise UndefinedPDEError: if the number of unknowns is not specified yet       :raise UndefinedPDEError: if the number of unknowns is not specified yet
1458       """       """
1459       if self.__numSolutions==None:       if self.__numSolutions==None:
1460          if self.__numEquations==None:          if self.__numEquations==None:
# Line 1424  class LinearProblem(object): Line 1467  class LinearProblem(object):
1467       """       """
1468       Returns the status of order reduction for the equation.       Returns the status of order reduction for the equation.
1469    
1470       @return: True if reduced interpolation order is used for the       :return: True if reduced interpolation order is used for the
1471                representation of the equation, False otherwise                representation of the equation, False otherwise
1472       @rtype: L{bool}       :rtype: `bool`
1473       """       """
1474       return self.__reduce_equation_order       return self.__reduce_equation_order
1475    
# Line 1434  class LinearProblem(object): Line 1477  class LinearProblem(object):
1477       """       """
1478       Returns the status of order reduction for the solution.       Returns the status of order reduction for the solution.
1479    
1480       @return: True if reduced interpolation order is used for the       :return: True if reduced interpolation order is used for the
1481                representation of the solution, False otherwise                representation of the solution, False otherwise
1482       @rtype: L{bool}       :rtype: `bool`
1483       """       """
1484       return self.__reduce_solution_order       return self.__reduce_solution_order
1485    
1486     def getFunctionSpaceForEquation(self):     def getFunctionSpaceForEquation(self):
1487       """       """
1488       Returns the L{FunctionSpace<escript.FunctionSpace>} used to discretize       Returns the `FunctionSpace` used to discretize
1489       the equation.       the equation.
1490    
1491       @return: representation space of equation       :return: representation space of equation
1492       @rtype: L{FunctionSpace<escript.FunctionSpace>}       :rtype: `FunctionSpace`
1493       """       """
1494       if self.reduceEquationOrder():       if self.reduceEquationOrder():
1495           return escript.ReducedSolution(self.getDomain())           return escript.ReducedSolution(self.getDomain())
# Line 1455  class LinearProblem(object): Line 1498  class LinearProblem(object):
1498    
1499     def getFunctionSpaceForSolution(self):     def getFunctionSpaceForSolution(self):
1500       """       """
1501       Returns the L{FunctionSpace<escript.FunctionSpace>} used to represent       Returns the `FunctionSpace` used to represent
1502       the solution.       the solution.
1503    
1504       @return: representation space of solution       :return: representation space of solution
1505       @rtype: L{FunctionSpace<escript.FunctionSpace>}       :rtype: `FunctionSpace`
1506       """       """
1507       if self.reduceSolutionOrder():       if self.reduceSolutionOrder():
1508           return escript.ReducedSolution(self.getDomain())           return escript.ReducedSolution(self.getDomain())
# Line 1473  class LinearProblem(object): Line 1516  class LinearProblem(object):
1516         """         """
1517         Sets the solver options.         Sets the solver options.
1518    
1519         @param options: the new solver options. If equal C{None}, the solver options are set to the default.         :param options: the new solver options. If equal ``None``, the solver options are set to the default.
1520         @type options: L{SolverOptions} or C{None}         :type options: `SolverOptions` or ``None``
1521         @note: The symmetry flag of options is overwritten by the symmetry flag of the L{LinearProblem}.         :note: The symmetry flag of options is overwritten by the symmetry flag of the `LinearProblem`.
1522         """         """
1523         if options==None:         if options==None:
1524            self.__solver_options=SolverOptions()            self.__solver_options=SolverOptions()
# Line 1489  class LinearProblem(object): Line 1532  class LinearProblem(object):
1532         """         """
1533         Returns the solver options         Returns the solver options
1534        
1535         @rtype: L{SolverOptions}         :rtype: `SolverOptions`
1536         """         """
1537         self.__solver_options.setSymmetry(self.__sym)         self.__solver_options.setSymmetry(self.__sym)
1538         return self.__solver_options         return self.__solver_options
# Line 1498  class LinearProblem(object): Line 1541  class LinearProblem(object):
1541        """        """
1542        Checks if matrix lumping is the current solver method.        Checks if matrix lumping is the current solver method.
1543    
1544        @return: True if the current solver method is lumping        :return: True if the current solver method is lumping
1545        @rtype: C{bool}        :rtype: ``bool``
1546        """        """
1547        return self.getSolverOptions().getSolverMethod()==self.getSolverOptions().LUMPING        return self.getSolverOptions().getSolverMethod()==self.getSolverOptions().LUMPING
1548     # ==========================================================================     # ==========================================================================
# Line 1509  class LinearProblem(object): Line 1552  class LinearProblem(object):
1552        """        """
1553        Checks if symmetry is indicated.        Checks if symmetry is indicated.
1554    
1555        @return: True if a symmetric PDE is indicated, False otherwise        :return: True if a symmetric PDE is indicated, False otherwise
1556        @rtype: C{bool}        :rtype: ``bool``
1557        @note: the method is equivalent to use getSolverOptions().isSymmetric()        :note: the method is equivalent to use getSolverOptions().isSymmetric()
1558        """        """
1559        self.getSolverOptions().isSymmetric()        self.getSolverOptions().isSymmetric()
1560    
1561     def setSymmetryOn(self):     def setSymmetryOn(self):
1562        """        """
1563        Sets the symmetry flag.        Sets the symmetry flag.
1564        @note: The method overwrites the symmetry flag set by the solver options        :note: The method overwrites the symmetry flag set by the solver options
1565        """        """
1566        self.__sym=True        self.__sym=True
1567        self.getSolverOptions().setSymmetryOn()        self.getSolverOptions().setSymmetryOn()
# Line 1526  class LinearProblem(object): Line 1569  class LinearProblem(object):
1569     def setSymmetryOff(self):     def setSymmetryOff(self):
1570        """        """
1571        Clears the symmetry flag.        Clears the symmetry flag.
1572        @note: The method overwrites the symmetry flag set by the solver options        :note: The method overwrites the symmetry flag set by the solver options
1573        """        """
1574        self.__sym=False        self.__sym=False
1575        self.getSolverOptions().setSymmetryOff()        self.getSolverOptions().setSymmetryOff()
1576    
1577     def setSymmetry(self,flag=False):     def setSymmetry(self,flag=False):
1578        """        """
1579        Sets the symmetry flag to C{flag}.        Sets the symmetry flag to ``flag``.
1580    
1581        @param flag: If True, the symmetry flag is set otherwise reset.        :param flag: If True, the symmetry flag is set otherwise reset.
1582        @type flag: C{bool}        :type flag: ``bool``
1583        @note: The method overwrites the symmetry flag set by the solver options        :note: The method overwrites the symmetry flag set by the solver options
1584        """        """
1585        self.getSolverOptions().setSymmetry(flag)        self.getSolverOptions().setSymmetry(flag)
1586     # ==========================================================================     # ==========================================================================
# Line 1547  class LinearProblem(object): Line 1590  class LinearProblem(object):
1590       """       """
1591       Switches reduced order on for solution and equation representation.       Switches reduced order on for solution and equation representation.
1592    
1593       @raise RuntimeError: if order reduction is altered after a coefficient has       :raise RuntimeError: if order reduction is altered after a coefficient has
1594                            been set                            been set
1595       """       """
1596       self.setReducedOrderForSolutionOn()       self.setReducedOrderForSolutionOn()
# Line 1557  class LinearProblem(object): Line 1600  class LinearProblem(object):
1600       """       """
1601       Switches reduced order off for solution and equation representation       Switches reduced order off for solution and equation representation
1602    
1603       @raise RuntimeError: if order reduction is altered after a coefficient has       :raise RuntimeError: if order reduction is altered after a coefficient has
1604                            been set                            been set
1605       """       """
1606       self.setReducedOrderForSolutionOff()       self.setReducedOrderForSolutionOff()
# Line 1568  class LinearProblem(object): Line 1611  class LinearProblem(object):
1611       Sets order reduction state for both solution and equation representation       Sets order reduction state for both solution and equation representation
1612       according to flag.       according to flag.
1613    
1614       @param flag: if True, the order reduction is switched on for both solution       :param flag: if True, the order reduction is switched on for both solution
1615                    and equation representation, otherwise or if flag is not                    and equation representation, otherwise or if flag is not
1616                    present order reduction is switched off                    present order reduction is switched off
1617       @type flag: C{bool}       :type flag: ``bool``
1618       @raise RuntimeError: if order reduction is altered after a coefficient has       :raise RuntimeError: if order reduction is altered after a coefficient has
1619                            been set                            been set
1620       """       """
1621       self.setReducedOrderForSolutionTo(flag)       self.setReducedOrderForSolutionTo(flag)
# Line 1583  class LinearProblem(object): Line 1626  class LinearProblem(object):
1626       """       """
1627       Switches reduced order on for solution representation.       Switches reduced order on for solution representation.
1628    
1629       @raise RuntimeError: if order reduction is altered after a coefficient has       :raise RuntimeError: if order reduction is altered after a coefficient has
1630                            been set                            been set
1631       """       """
1632       if not self.__reduce_solution_order:       if not self.__reduce_solution_order:
# Line 1597  class LinearProblem(object): Line 1640  class LinearProblem(object):
1640       """       """
1641       Switches reduced order off for solution representation       Switches reduced order off for solution representation
1642    
1643       @raise RuntimeError: if order reduction is altered after a coefficient has       :raise RuntimeError: if order reduction is altered after a coefficient has
1644                            been set.                            been set.
1645       """       """
1646       if self.__reduce_solution_order:       if self.__reduce_solution_order:
# Line 1611  class LinearProblem(object): Line 1654  class LinearProblem(object):
1654       """       """
1655       Sets order reduction state for solution representation according to flag.       Sets order reduction state for solution representation according to flag.
1656    
1657       @param flag: if flag is True, the order reduction is switched on for       :param flag: if flag is True, the order reduction is switched on for
1658                    solution representation, otherwise or if flag is not present                    solution representation, otherwise or if flag is not present
1659                    order reduction is switched off                    order reduction is switched off
1660       @type flag: C{bool}       :type flag: ``bool``
1661       @raise RuntimeError: if order reduction is altered after a coefficient has       :raise RuntimeError: if order reduction is altered after a coefficient has
1662                            been set                            been set
1663       """       """
1664       if flag:       if flag:
# Line 1627  class LinearProblem(object): Line 1670  class LinearProblem(object):
1670       """       """
1671       Switches reduced order on for equation representation.       Switches reduced order on for equation representation.
1672    
1673       @raise RuntimeError: if order reduction is altered after a coefficient has       :raise RuntimeError: if order reduction is altered after a coefficient has
1674                            been set                            been set
1675       """       """
1676       if not self.__reduce_equation_order:       if not self.__reduce_equation_order:
# Line 1641  class LinearProblem(object): Line 1684  class LinearProblem(object):
1684       """       """
1685       Switches reduced order off for equation representation.       Switches reduced order off for equation representation.
1686    
1687       @raise RuntimeError: if order reduction is altered after a coefficient has       :raise RuntimeError: if order reduction is altered after a coefficient has
1688                            been set                            been set
1689       """       """
1690       if self.__reduce_equation_order:       if self.__reduce_equation_order:
# Line 1655  class LinearProblem(object): Line 1698  class LinearProblem(object):
1698       """       """
1699       Sets order reduction state for equation representation according to flag.       Sets order reduction state for equation representation according to flag.
1700    
1701       @param flag: if flag is True, the order reduction is switched on for       :param flag: if flag is True, the order reduction is switched on for
1702                    equation representation, otherwise or if flag is not present                    equation representation, otherwise or if flag is not present
1703                    order reduction is switched off                    order reduction is switched off
1704       @type flag: C{bool}       :type flag: ``bool``
1705       @raise RuntimeError: if order reduction is altered after a coefficient has       :raise RuntimeError: if order reduction is altered after a coefficient has
1706                            been set                            been set
1707       """       """
1708       if flag:       if flag:
# Line 1676  class LinearProblem(object): Line 1719  class LinearProblem(object):
1719        """        """
1720        Tests a coefficient for symmetry.        Tests a coefficient for symmetry.
1721    
1722        @param name: name of the coefficient        :param name: name of the coefficient
1723        @type name: C{str}        :type name: ``str``
1724        @param verbose: if set to True or not present a report on coefficients        :param verbose: if set to True or not present a report on coefficients
1725                        which break the symmetry is printed.                        which break the symmetry is printed.
1726        @type verbose: C{bool}        :type verbose: ``bool``
1727        @return: True if coefficient C{name} is symmetric        :return: True if coefficient ``name`` is symmetric
1728        @rtype: C{bool}        :rtype: ``bool``
1729        """        """
1730        SMALL_TOLERANCE=util.EPSILON*10.        SMALL_TOLERANCE=util.EPSILON*10.
1731        A=self.getCoefficient(name)        A=self.getCoefficient(name)
# Line 1723  class LinearProblem(object): Line 1766  class LinearProblem(object):
1766        """        """
1767        Tests two coefficients for reciprocal symmetry.        Tests two coefficients for reciprocal symmetry.
1768    
1769        @param name0: name of the first coefficient        :param name0: name of the first coefficient
1770        @type name0: C{str}        :type name0: ``str``
1771        @param name1: name of the second coefficient        :param name1: name of the second coefficient
1772        @type name1: C{str}        :type name1: ``str``
1773        @param verbose: if set to True or not present a report on coefficients        :param verbose: if set to True or not present a report on coefficients
1774                        which break the symmetry is printed                        which break the symmetry is printed
1775        @type verbose: C{bool}        :type verbose: ``bool``
1776        @return: True if coefficients C{name0} and C{name1} are reciprocally        :return: True if coefficients ``name0`` and ``name1`` are reciprocally
1777                 symmetric.                 symmetric.
1778        @rtype: C{bool}        :rtype: ``bool``
1779        """        """
1780        SMALL_TOLERANCE=util.EPSILON*10.        SMALL_TOLERANCE=util.EPSILON*10.
1781        B=self.getCoefficient(name0)        B=self.getCoefficient(name0)
# Line 1781  class LinearProblem(object): Line 1824  class LinearProblem(object):
1824    
1825     def getCoefficient(self,name):     def getCoefficient(self,name):
1826       """       """
1827       Returns the value of the coefficient C{name}.       Returns the value of the coefficient ``name``.
1828    
1829       @param name: name of the coefficient requested       :param name: name of the coefficient requested
1830       @type name: C{string}       :type name: ``string``
1831       @return: the value of the coefficient       :return: the value of the coefficient
1832       @rtype: L{Data<escript.Data>}       :rtype: `Data`
1833       @raise IllegalCoefficient: if C{name} is not a coefficient of the PDE       :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
1834       """       """
1835       if self.hasCoefficient(name):       if self.hasCoefficient(name):
1836           return self.__COEFFICIENTS[name].getValue()           return self.__COEFFICIENTS[name].getValue()
# Line 1796  class LinearProblem(object): Line 1839  class LinearProblem(object):
1839    
1840     def hasCoefficient(self,name):     def hasCoefficient(self,name):
1841       """       """
1842       Returns True if C{name} is the name of a coefficient.       Returns True if ``name`` is the name of a coefficient.
1843    
1844       @param name: name of the coefficient enquired       :param name: name of the coefficient enquired
1845       @type name: C{string}       :type name: ``string``
1846       @return: True if C{name} is the name of a coefficient of the general PDE,       :return: True if ``name`` is the name of a coefficient of the general PDE,
1847                False otherwise                False otherwise
1848       @rtype: C{bool}       :rtype: ``bool``
1849       """       """
1850       return self.__COEFFICIENTS.has_key(name)       return self.__COEFFICIENTS.has_key(name)
1851    
1852     def createCoefficient(self, name):     def createCoefficient(self, name):
1853       """       """
1854       Creates a L{Data<escript.Data>} object corresponding to coefficient       Creates a `Data` object corresponding to coefficient
1855       C{name}.       ``name``.
1856    
1857       @return: the coefficient C{name} initialized to 0       :return: the coefficient ``name`` initialized to 0
1858       @rtype: L{Data<escript.Data>}       :rtype: `Data`
1859       @raise IllegalCoefficient: if C{name} is not a coefficient of the PDE       :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
1860       """       """
1861       if self.hasCoefficient(name):       if self.hasCoefficient(name):
1862          return escript.Data(0.,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))          return escript.Data(0.,self.getShapeOfCoefficient(name),self.getFunctionSpaceForCoefficient(name))
# Line 1822  class LinearProblem(object): Line 1865  class LinearProblem(object):
1865    
1866     def getFunctionSpaceForCoefficient(self,name):     def getFunctionSpaceForCoefficient(self,name):
1867       """       """
1868       Returns the L{FunctionSpace<escript.FunctionSpace>} to be used for       Returns the `FunctionSpace` to be used for
1869       coefficient C{name}.       coefficient ``name``.
1870    
1871       @param name: name of the coefficient enquired       :param name: name of the coefficient enquired
1872       @type name: C{string}       :type name: ``string``
1873       @return: the function space to be used for coefficient C{name}       :return: the function space to be used for coefficient ``name``
1874       @rtype: L{FunctionSpace<escript.FunctionSpace>}       :rtype: `FunctionSpace`
1875       @raise IllegalCoefficient: if C{name} is not a coefficient of the PDE       :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
1876       """       """
1877       if self.hasCoefficient(name):       if self.hasCoefficient(name):
1878          return self.__COEFFICIENTS[name].getFunctionSpace(self.getDomain())          return self.__COEFFICIENTS[name].getFunctionSpace(self.getDomain())
# Line 1838  class LinearProblem(object): Line 1881  class LinearProblem(object):
1881    
1882     def getShapeOfCoefficient(self,name):     def getShapeOfCoefficient(self,name):
1883       """       """
1884       Returns the shape of the coefficient C{name}.       Returns the shape of the coefficient ``name``.
1885    
1886       @param name: name of the coefficient enquired       :param name: name of the coefficient enquired
1887       @type name: C{string}       :type name: ``string``
1888       @return: the shape of the coefficient C{name}       :return: the shape of the coefficient ``name``
1889       @rtype: C{tuple} of C{int}       :rtype: ``tuple`` of ``int``
1890       @raise IllegalCoefficient: if C{name} is not a coefficient of the PDE       :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
1891       """       """
1892       if self.hasCoefficient(name):       if self.hasCoefficient(name):
1893          return self.__COEFFICIENTS[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())          return self.__COEFFICIENTS[name].getShape(self.getDomain(),self.getNumEquations(),self.getNumSolutions())
# Line 1860  class LinearProblem(object): Line 1903  class LinearProblem(object):
1903    
1904     def alteredCoefficient(self,name):     def alteredCoefficient(self,name):
1905       """       """
1906       Announces that coefficient C{name} has been changed.       Announces that coefficient ``name`` has been changed.
1907    
1908       @param name: name of the coefficient affected       :param name: name of the coefficient affected
1909       @type name: C{string}       :type name: ``string``
1910       @raise IllegalCoefficient: if C{name} is not a coefficient of the PDE       :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
1911       @note: if C{name} is q or r, the method will not trigger a rebuild of the       :note: if ``name`` is q or r, the method will not trigger a rebuild of the
1912              system as constraints are applied to the solved system.              system as constraints are applied to the solved system.
1913       """       """
1914       if self.hasCoefficient(name):       if self.hasCoefficient(name):
# Line 1972  class LinearProblem(object): Line 2015  class LinearProblem(object):
2015       """       """
2016       Returns the operator of the linear problem.       Returns the operator of the linear problem.
2017    
2018       @return: the operator of the problem       :return: the operator of the problem
2019       """       """
2020       return self.getSystem()[0]       return self.getSystem()[0]
2021    
# Line 1980  class LinearProblem(object): Line 2023  class LinearProblem(object):
2023       """       """
2024       Returns the right hand side of the linear problem.       Returns the right hand side of the linear problem.
2025    
2026       @return: the right hand side of the problem       :return: the right hand side of the problem
2027       @rtype: L{Data<escript.Data>}       :rtype: `Data`
2028       """       """
2029       return self.getSystem()[1]       return self.getSystem()[1]
2030    
# Line 2015  class LinearProblem(object): Line 2058  class LinearProblem(object):
2058             self.__solution.setToZero()             self.__solution.setToZero()
2059             self.trace("Solution is reset to zero.")             self.trace("Solution is reset to zero.")
2060    
2061     def setSolution(self,u):     def setSolution(self,u, validate=True):
2062         """         """
2063         Sets the solution assuming that makes the system valid with the tolrance         Sets the solution assuming that makes the system valid with the tolrance
2064         defined by the solver options         defined by the solver options
2065         """         """
2066         self.__solution_rtol=self.getSolverOptions().getTolerance()         if validate:
2067         self.__solution_atol=self.getSolverOptions().getAbsoluteTolerance()        self.__solution_rtol=self.getSolverOptions().getTolerance()
2068          self.__solution_atol=self.getSolverOptions().getAbsoluteTolerance()
2069          self.validSolution()
2070         self.__solution=u         self.__solution=u
        self.validSolution()  
   
2071     def getCurrentSolution(self):     def getCurrentSolution(self):
2072         """         """
2073         Returns the solution in its current state.         Returns the solution in its current state.
# Line 2064  class LinearProblem(object): Line 2107  class LinearProblem(object):
2107             if self.isUsingLumping():             if self.isUsingLumping():
2108                 self.__operator.setToZero()                 self.__operator.setToZero()
2109             else:             else:
2110                 self.__operator.resetValues()                 if self.getOperatorType() == self.getRequiredOperatorType():
2111                       self.__operator.resetValues()
2112                   else:
2113                       self.__operator=self.createOperator()
2114                   self.__operator_type=self.getRequiredOperatorType()
2115             self.trace("Operator reset to zero")             self.trace("Operator reset to zero")
2116    
2117     def getCurrentOperator(self):     def getCurrentOperator(self):
# Line 2077  class LinearProblem(object): Line 2124  class LinearProblem(object):
2124        """        """
2125        Sets new values to coefficients.        Sets new values to coefficients.
2126    
2127        @raise IllegalCoefficient: if an unknown coefficient keyword is used        :raise IllegalCoefficient: if an unknown coefficient keyword is used
2128        """        """
2129        # check if the coefficients are  legal:        # check if the coefficients are  legal:
2130        for i in coefficients.iterkeys():        for i in coefficients.iterkeys():
# Line 2136  class LinearProblem(object): Line 2183  class LinearProblem(object):
2183        """        """
2184        Returns the system type which needs to be used by the current set up.        Returns the system type which needs to be used by the current set up.
2185    
2186        @note: Typically this method is overwritten when implementing a        :note: Typically this method is overwritten when implementing a
2187               particular linear problem.               particular linear problem.
2188        """        """
2189        return None        return None
# Line 2145  class LinearProblem(object): Line 2192  class LinearProblem(object):
2192         """         """
2193         Returns an instance of a new operator.         Returns an instance of a new operator.
2194    
2195         @note: This method is overwritten when implementing a particular         :note: This method is overwritten when implementing a particular
2196                linear problem.                linear problem.
2197         """         """
2198         return escript.Operator()         return escript.Operator()
# Line 2154  class LinearProblem(object): Line 2201  class LinearProblem(object):
2201        """        """
2202        Tests the PDE for symmetry.        Tests the PDE for symmetry.
2203    
2204        @param verbose: if set to True or not present a report on coefficients        :param verbose: if set to True or not present a report on coefficients
2205                        which break the symmetry is printed                        which break the symmetry is printed
2206        @type verbose: C{bool}        :type verbose: ``bool``
2207        @return: True if the problem is symmetric        :return: True if the problem is symmetric
2208        @rtype: C{bool}        :rtype: ``bool``
2209        @note: Typically this method is overwritten when implementing a        :note: Typically this method is overwritten when implementing a
2210               particular linear problem.               particular linear problem.
2211        """        """
2212        out=True        out=True
# Line 2169  class LinearProblem(object): Line 2216  class LinearProblem(object):
2216         """         """
2217         Returns the solution of the problem.         Returns the solution of the problem.
2218    
2219         @return: the solution         :return: the solution
2220         @rtype: L{Data<escript.Data>}         :rtype: `Data`
2221    
2222         @note: This method is overwritten when implementing a particular         :note: This method is overwritten when implementing a particular
2223                linear problem.                linear problem.
2224         """         """
2225         return self.getCurrentSolution()         return self.getCurrentSolution()
# Line 2181  class LinearProblem(object): Line 2228  class LinearProblem(object):
2228         """         """
2229         Returns the operator and right hand side of the PDE.         Returns the operator and right hand side of the PDE.
2230    
2231         @return: the discrete version of the PDE         :return: the discrete version of the PDE
2232         @rtype: C{tuple} of L{Operator,<escript.Operator>} and L{Data<escript.Data>}.         :rtype: ``tuple`` of `Operator` and `Data`.
2233    
2234         @note: This method is overwritten when implementing a particular         :note: This method is overwritten when implementing a particular
2235                linear problem.                linear problem.
2236         """         """
2237         return (self.getCurrentOperator(), self.getCurrentRightHandSide())         return (self.getCurrentOperator(), self.getCurrentRightHandSide())
# Line 2192  class LinearProblem(object): Line 2239  class LinearProblem(object):
2239  class LinearPDE(LinearProblem):  class LinearPDE(LinearProblem):
2240     """     """
2241     This class is used to define a general linear, steady, second order PDE     This class is used to define a general linear, steady, second order PDE
2242     for an unknown function M{u} on a given domain defined through a     for an unknown function *u* on a given domain defined through a
2243     L{Domain<escript.Domain>} object.     `Domain` object.
2244    
2245     For a single PDE having a solution with a single component the linear PDE     For a single PDE having a solution with a single component the linear PDE
2246     is defined in the following form:     is defined in the following form:
2247    
2248     M{-(grad(A[j,l]+A_reduced[j,l])*grad(u)[l]+(B[j]+B_reduced[j])u)[j]+(C[l]+C_reduced[l])*grad(u)[l]+(D+D_reduced)=-grad(X+X_reduced)[j,j]+(Y+Y_reduced)}     *-(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)*
2249    
2250     where M{grad(F)} denotes the spatial derivative of M{F}. Einstein's     where *grad(F)* denotes the spatial derivative of *F*. Einstein's
2251     summation convention, ie. summation over indexes appearing twice in a term     summation convention, ie. summation over indexes appearing twice in a term
2252     of a sum performed, is used.     of a sum performed, is used.
2253     The coefficients M{A}, M{B}, M{C}, M{D}, M{X} and M{Y} have to be specified     The coefficients *A*, *B*, *C*, *D*, *X* and *Y* have to be specified
2254     through L{Data<escript.Data>} objects in L{Function<escript.Function>} and     through `Data` objects in `Function` and
2255     the coefficients M{A_reduced}, M{B_reduced}, M{C_reduced}, M{D_reduced},     the coefficients *A_reduced*, *B_reduced*, *C_reduced*, *D_reduced*,
2256     M{X_reduced} and M{Y_reduced} have to be specified through     *X_reduced* and *Y_reduced* have to be specified through
2257     L{Data<escript.Data>} objects in L{ReducedFunction<escript.ReducedFunction>}.     `Data` objects in `ReducedFunction`.
2258     It is also allowed to use objects that can be converted into such     It is also allowed to use objects that can be converted into such
2259     L{Data<escript.Data>} objects. M{A} and M{A_reduced} are rank two, M{B},     `Data` objects. *A* and *A_reduced* are rank two, *B*,
2260     M{C}, M{X}, M{B_reduced}, M{C_reduced} and M{X_reduced} are rank one and     *C*, *X*, *B_reduced*, *C_reduced* and *X_reduced* are rank one and
2261     M{D}, M{D_reduced}, M{Y} and M{Y_reduced} are scalar.     *D*, *D_reduced*, *Y* and *Y_reduced* are scalar.
2262    
2263     The following natural boundary conditions are considered:     The following natural boundary conditions are considered:
2264    
2265     M{n[j]*((A[i,j]+A_reduced[i,j])*grad(u)[l]+(B+B_reduced)[j]*u)+(d+d_reduced)*u=n[j]*(X[j]+X_reduced[j])+y}     *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*
2266    
2267     where M{n} is the outer normal field. Notice that the coefficients M{A},     where *n* is the outer normal field. Notice that the coefficients *A*,
2268     M{A_reduced}, M{B}, M{B_reduced}, M{X} and M{X_reduced} are defined in the     *A_reduced*, *B*, *B_reduced*, *X* and *X_reduced* are defined in the
2269     PDE. The coefficients M{d} and M{y} are each a scalar in     PDE. The coefficients *d* and *y* are each a scalar in
2270     L{FunctionOnBoundary<escript.FunctionOnBoundary>} and the coefficients     `FunctionOnBoundary` and the coefficients
2271     M{d_reduced} and M{y_reduced} are each a scalar in     *d_reduced* and *y_reduced* are each a scalar in
2272     L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.     `ReducedFunctionOnBoundary`.
2273    
2274     Constraints for the solution prescribe the value of the solution at certain     Constraints for the solution prescribe the value of the solution at certain
2275     locations in the domain. They have the form     locations in the domain. They have the form
2276    
2277     M{u=r} where M{q>0}     *u=r* where *q>0*
2278    
2279     M{r} and M{q} are each scalar where M{q} is the characteristic function     *r* and *q* are each scalar where *q* is the characteristic function
2280     defining where the constraint is applied. The constraints override any     defining where the constraint is applied. The constraints override any
2281     other condition set by the PDE or the boundary condition.     other condition set by the PDE or the boundary condition.
2282    
2283     The PDE is symmetrical if     The PDE is symmetrical if
2284    
2285     M{A[i,j]=A[j,i]}  and M{B[j]=C[j]} and M{A_reduced[i,j]=A_reduced[j,i]}     *A[i,j]=A[j,i]*  and *B[j]=C[j]* and *A_reduced[i,j]=A_reduced[j,i]*
2286     and M{B_reduced[j]=C_reduced[j]}     and *B_reduced[j]=C_reduced[j]*
2287    
2288     For a system of PDEs and a solution with several components the PDE has the     For a system of PDEs and a solution with several components the PDE has the
2289     form     form
2290    
2291     M{-grad((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k])[j]+(C[i,k,l]+C_reduced[i,k,l])*grad(u[k])[l]+(D[i,k]+D_reduced[i,k]*u[k] =-grad(X[i,j]+X_reduced[i,j])[j]+Y[i]+Y_reduced[i] }     *-grad((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k])[j]+(C[i,k,l]+C_reduced[i,k,l])*grad(u[k])[l]+(D[i,k]+D_reduced[i,k]*u[k] =-grad(X[i,j]+X_reduced[i,j])[j]+Y[i]+Y_reduced[i]*
2292    
2293     M{A} and M{A_reduced} are of rank four, M{B}, M{B_reduced}, M{C} and     *A* and *A_reduced* are of rank four, *B*, *B_reduced*, *C* and
2294     M{C_reduced} are each of rank three, M{D}, M{D_reduced}, M{X_reduced} and     *C_reduced* are each of rank three, *D*, *D_reduced*, *X_reduced* and
2295     M{X} are each of rank two and M{Y} and M{Y_reduced} are of rank one.     *X* are each of rank two and *Y* and *Y_reduced* are of rank one.
2296     The natural boundary conditions take the form:     The natural boundary conditions take the form:
2297    
2298     M{n[j]*((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k])+(d[i,k]+d_reduced[i,k])*u[k]=n[j]*(X[i,j]+X_reduced[i,j])+y[i]+y_reduced[i]}     *n[j]*((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k])+(d[i,k]+d_reduced[i,k])*u[k]=n[j]*(X[i,j]+X_reduced[i,j])+y[i]+y_reduced[i]*
2299    
2300     The coefficient M{d} is of rank two and M{y} is of rank one both in     The coefficient *d* is of rank two and *y* is of rank one both in
2301     L{FunctionOnBoundary<escript.FunctionOnBoundary>}. The coefficients     `FunctionOnBoundary`. The coefficients
2302     M{d_reduced} is of rank two and M{y_reduced} is of rank one both in     *d_reduced* is of rank two and *y_reduced* is of rank one both in
2303     L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.     `ReducedFunctionOnBoundary`.
2304    
2305     Constraints take the form     Constraints take the form
2306    
2307     M{u[i]=r[i]}  where  M{q[i]>0}     *u[i]=r[i]*  where  *q[i]>0*
2308    
2309     M{r} and M{q} are each rank one. Notice that at some locations not     *r* and *q* are each rank one. Notice that at some locations not
2310     necessarily all components must have a constraint.     necessarily all components must have a constraint.
2311    
2312     The system of PDEs is symmetrical if     The system of PDEs is symmetrical if
2313    
2314        - M{A[i,j,k,l]=A[k,l,i,j]}        - *A[i,j,k,l]=A[k,l,i,j]*
2315        - M{A_reduced[i,j,k,l]=A_reduced[k,l,i,j]}        - *A_reduced[i,j,k,l]=A_reduced[k,l,i,j]*
2316        - M{B[i,j,k]=C[k,i,j]}        - *B[i,j,k]=C[k,i,j]*
2317        - M{B_reduced[i,j,k]=C_reduced[k,i,j]}        - *B_reduced[i,j,k]=C_reduced[k,i,j]*
2318        - M{D[i,k]=D[i,k]}        - *D[i,k]=D[i,k]*
2319        - M{D_reduced[i,k]=D_reduced[i,k]}        - *D_reduced[i,k]=D_reduced[i,k]*
2320        - M{d[i,k]=d[k,i]}        - *d[i,k]=d[k,i]*
2321        - M{d_reduced[i,k]=d_reduced[k,i]}        - *d_reduced[i,k]=d_reduced[k,i]*
2322    
2323     L{LinearPDE} also supports solution discontinuities over a contact region     `LinearPDE` also supports solution discontinuities over a contact region
2324     in the domain. To specify the conditions across the discontinuity we are     in the domain. To specify the conditions across the discontinuity we are
2325     using the generalised flux M{J} which, in the case of a system of PDEs     using the generalised flux *J* which, in the case of a system of PDEs
2326     and several components of the solution, is defined as     and several components of the solution, is defined as
2327    
2328     M{J[i,j]=(A[i,j,k,l]+A_reduced[[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k]-X[i,j]-X_reduced[i,j]}     *J[i,j]=(A[i,j,k,l]+A_reduced[[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k]-X[i,j]-X_reduced[i,j]*
2329    
2330     For the case of single solution component and single PDE M{J} is defined as     For the case of single solution component and single PDE *J* is defined as
2331    
2332     M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[j]+(B[i]+B_reduced[i])*u-X[i]-X_reduced[i]}     *J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[j]+(B[i]+B_reduced[i])*u-X[i]-X_reduced[i]*
2333    
2334     In the context of discontinuities M{n} denotes the normal on the     In the context of discontinuities *n* denotes the normal on the
2335     discontinuity pointing from side 0 towards side 1 calculated from     discontinuity pointing from side 0 towards side 1 calculated from
2336     L{getNormal<escript.FunctionSpace.getNormal>} of L{FunctionOnContactZero<escript.FunctionOnContactZero>}.     `FunctionSpace.getNormal` of `FunctionOnContactZero`.
2337     For a system of PDEs the contact condition takes the form     For a system of PDEs the contact condition takes the form
2338    
2339     M{n[j]*J0[i,j]=n[j]*J1[i,j]=(y_contact[i]+y_contact_reduced[i])- (d_contact[i,k]+d_contact_reduced[i,k])*jump(u)[k]}     *n[j]*J0[i,j]=n[j]*J1[i,j]=(y_contact[i]+y_contact_reduced[i])- (d_contact[i,k]+d_contact_reduced[i,k])*jump(u)[k]*
2340    
2341     where M{J0} and M{J1} are the fluxes on side 0 and side 1 of the     where *J0* and *J1* are the fluxes on side 0 and side 1 of the
2342     discontinuity, respectively. M{jump(u)}, which is the difference of the     discontinuity, respectively. *jump(u)*, which is the difference of the
2343     solution at side 1 and at side 0, denotes the jump of M{u} across     solution at side 1 and at side 0, denotes the jump of *u* across
2344     discontinuity along the normal calculated by L{jump<util.jump>}.     discontinuity along the normal calculated by `jump`.
2345     The coefficient M{d_contact} is of rank two and M{y_contact} is of rank one     The coefficient *d_contact* is of rank two and *y_contact* is of rank one
2346     both in L{FunctionOnContactZero<escript.FunctionOnContactZero>} or     both in `FunctionOnContactZero` or
2347     L{FunctionOnContactOne<escript.FunctionOnContactOne>}.     `FunctionOnContactOne`.
2348     The coefficient M{d_contact_reduced} is of rank two and M{y_contact_reduced}     The coefficient *d_contact_reduced* is of rank two and *y_contact_reduced*
2349     is of rank one both in L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>}     is of rank one both in `ReducedFunctionOnContactZero`
2350     or L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.     or `ReducedFunctionOnContactOne`.
2351     In case of a single PDE and a single component solution the contact     In case of a single PDE and a single component solution the contact
2352     condition takes the form     condition takes the form
2353    
2354     M{n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)}     *n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)*
2355    
2356     In this case the coefficient M{d_contact} and M{y_contact} are each scalar     In this case the coefficient *d_contact* and *y_contact* are each scalar
2357     both in L{FunctionOnContactZero<escript.FunctionOnContactZero>} or     both in `FunctionOnContactZero` or
2358     L{FunctionOnContactOne<escript.FunctionOnContactOne>} and the coefficient     `FunctionOnContactOne` and the coefficient
2359     M{d_contact_reduced} and M{y_contact_reduced} are each scalar both in     *d_contact_reduced* and *y_contact_reduced* are each scalar both in
2360     L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or     `ReducedFunctionOnContactZero` or
2361     L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.     `ReducedFunctionOnContactOne`.
2362    
2363     Typical usage::     Typical usage::
2364    
# Line 2325  class LinearPDE(LinearProblem): Line 2372  class LinearPDE(LinearProblem):
2372       """       """
2373       Initializes a new linear PDE.       Initializes a new linear PDE.
2374    
2375       @param domain: domain of the PDE       :param domain: domain of the PDE
2376       @type domain: L{Domain<escript.Domain>}       :type domain: `Domain`
2377       @param numEquations: number of equations. If C{None} the number of       :param numEquations: number of equations. If ``None`` the number of
2378                            equations is extracted from the PDE coefficients.                            equations is extracted from the PDE coefficients.
2379       @param numSolutions: number of solution components. If C{None} the number       :param numSolutions: number of solution components. If ``None`` the number
2380                            of solution components is extracted from the PDE                            of solution components is extracted from the PDE
2381                            coefficients.                            coefficients.
2382       @param debug: if True debug information is printed       :param debug: if True debug information is printed
2383    
2384       """       """
2385       super(LinearPDE, self).__init__(domain,numEquations,numSolutions,debug)       super(LinearPDE, self).__init__(domain,numEquations,numSolutions,debug)
# Line 2367  class LinearPDE(LinearProblem): Line 2414  class LinearPDE(LinearProblem):
2414       """       """
2415       Returns the string representation of the PDE.       Returns the string representation of the PDE.
2416    
2417       @return: a simple representation of the PDE       :return: a simple representation of the PDE
2418       @rtype: C{str}       :rtype: ``str``
2419       """       """
2420       return "<LinearPDE %d>"%id(self)       return "<LinearPDE %d>"%id(self)
2421    
# Line 2383  class LinearPDE(LinearProblem): Line 2430  class LinearPDE(LinearProblem):
2430        """        """
2431        Tests the PDE for symmetry.        Tests the PDE for symmetry.
2432    
2433        @param verbose: if set to True or not present a report on coefficients        :param verbose: if set to True or not present a report on coefficients
2434                        which break the symmetry is printed.                        which break the symmetry is printed.
2435        @type verbose: C{bool}        :type verbose: ``bool``
2436        @return: True if the PDE is symmetric        :return: True if the PDE is symmetric
2437        @rtype: L{bool}        :rtype: `bool`
2438        @note: This is a very expensive operation. It should be used for        :note: This is a very expensive operation. It should be used for
2439               degugging only! The symmetry flag is not altered.               degugging only! The symmetry flag is not altered.
2440        """        """
2441        out=True        out=True
# Line 2421  class LinearPDE(LinearProblem): Line 2468  class LinearPDE(LinearProblem):
2468         """         """
2469         Returns the solution of the PDE.         Returns the solution of the PDE.
2470    
2471         @return: the solution         :return: the solution
2472         @rtype: L{Data<escript.Data>}         :rtype: `Data`
2473         """         """
2474         option_class=self.getSolverOptions()         option_class=self.getSolverOptions()
2475         if not self.isSolutionValid():         if not self.isSolutionValid():
# Line 2439  class LinearPDE(LinearProblem): Line 2486  class LinearPDE(LinearProblem):
2486         """         """
2487         Returns the operator and right hand side of the PDE.         Returns the operator and right hand side of the PDE.
2488    
2489         @return: the discrete version of the PDE         :return: the discrete version of the PDE
2490         @rtype: C{tuple} of L{Operator,<escript.Operator>} and         :rtype: ``tuple`` of `Operator` and
2491                 L{Data<escript.Data>}                 `Data`
2492         """         """
2493         if not self.isOperatorValid() or not self.isRightHandSideValid():         if not self.isOperatorValid() or not self.isRightHandSideValid():
2494            if self.isUsingLumping():            if self.isUsingLumping():
# Line 2500  class LinearPDE(LinearProblem): Line 2547  class LinearPDE(LinearProblem):
2547    
2548                   self.resetOperator()                   self.resetOperator()
2549                   operator=self.getCurrentOperator()                   operator=self.getCurrentOperator()
2550                   if False and hasattr(self.getDomain(), "addPDEToLumpedSystem") :                   if hasattr(self.getDomain(), "addPDEToLumpedSystem") :
2551                      self.getDomain().addPDEToLumpedSystem(operator, D_times_e, d_times_e)                      self.getDomain().addPDEToLumpedSystem(operator, D_times_e, d_times_e)
2552                      self.getDomain().addPDEToLumpedSystem(operator, D_reduced_times_e, d_reduced_times_e)                      self.getDomain().addPDEToLumpedSystem(operator, D_reduced_times_e, d_reduced_times_e)
2553                   else:                   else:
# Line 2616  class LinearPDE(LinearProblem): Line 2663  class LinearPDE(LinearProblem):
2663        """        """
2664        Applies the constraints defined by q and r to the PDE.        Applies the constraints defined by q and r to the PDE.
2665    
2666        @param rhs_only: if True only the right hand side is altered by the        :param rhs_only: if True only the right hand side is altered by the
2667                         constraint                         constraint
2668        @type rhs_only: C{bool}        :type rhs_only: ``bool``
2669        """        """
2670        q=self.getCoefficient("q")        q=self.getCoefficient("q")
2671        r=self.getCoefficient("r")        r=self.getCoefficient("r")
# Line 2646  class LinearPDE(LinearProblem): Line 2693  class LinearPDE(LinearProblem):
2693        """        """
2694        Sets new values to coefficients.        Sets new values to coefficients.
2695    
2696        @param coefficients: new values assigned to coefficients        :param coefficients: new values assigned to coefficients
2697        @keyword A: value for coefficient C{A}        :keyword A: value for coefficient ``A``
2698        @type A: any type that can be cast to a L{Data<escript.Data>} object on        :type A: any type that can be cast to a `Data` object on
2699                 L{Function<escript.Function>}                 `Function`
2700        @keyword A_reduced: value for coefficient C{A_reduced}        :keyword A_reduced: value for coefficient ``A_reduced``
2701        @type A_reduced: any type that can be cast to a L{Data<escript.Data>}        :type A_reduced: any type that can be cast to a `Data`
2702                         object on L{ReducedFunction<escript.ReducedFunction>}                         object on `ReducedFunction`
2703        @keyword B: value for coefficient C{B}        :keyword B: value for coefficient ``B``
2704        @type B: any type that can be cast to a L{Data<escript.Data>} object on        :type B: any type that can be cast to a `Data` object on
2705                 L{Function<escript.Function>}                 `Function`
2706        @keyword B_reduced: value for coefficient C{B_reduced}        :keyword B_reduced: value for coefficient ``B_reduced``
2707        @type B_reduced: any type that can be cast to a L{Data<escript.Data>}        :type B_reduced: any type that can be cast to a `Data`
2708                         object on L{ReducedFunction<escript.ReducedFunction>}                         object on `ReducedFunction`
2709        @keyword C: value for coefficient C{C}        :keyword C: value for coefficient ``C``
2710        @type C: any type that can be cast to a L{Data<escript.Data>} object on        :type C: any type that can be cast to a `Data` object on
2711                 L{Function<escript.Function>}                 `Function`
2712        @keyword C_reduced: value for coefficient C{C_reduced}        :keyword C_reduced: value for coefficient ``C_reduced``
2713        @type C_reduced: any type that can be cast to a L{Data<escript.Data>}        :type C_reduced: any type that can be cast to a `Data`
2714                         object on L{ReducedFunction<escript.ReducedFunction>}                         object on `ReducedFunction`
2715        @keyword D: value for coefficient C{D}        :keyword D: value for coefficient ``D``
2716        @type D: any type that can be cast to a L{Data<escript.Data>} object on        :type D: any type that can be cast to a `Data` object on
2717                 L{Function<escript.Function>}                 `Function`
2718        @keyword D_reduced: value for coefficient C{D_reduced}        :keyword D_reduced: value for coefficient ``D_reduced``
2719        @type D_reduced: any type that can be cast to a L{Data<escript.Data>}        :type D_reduced: any type that can be cast to a `Data`
2720                         object on L{ReducedFunction<escript.ReducedFunction>}                         object on `ReducedFunction`
2721        @keyword X: value for coefficient C{X}        :keyword X: value for coefficient ``X``
2722        @type X: any type that can be cast to a L{Data<escript.Data>} object on        :type X: any type that can be cast to a `Data` object on
2723                 L{Function<escript.Function>}                 `Function`
2724        @keyword X_reduced: value for coefficient C{X_reduced}        :keyword X_reduced: value for coefficient ``X_reduced``
2725        @type X_reduced: any type that can be cast to a L{Data<escript.Data>}        :type X_reduced: any type that can be cast to a `Data`
2726                         object on L{ReducedFunction<escript.ReducedFunction>}                         object on `ReducedFunction`
2727        @keyword Y: value for coefficient C{Y}        :keyword Y: value for coefficient ``Y``
2728        @type Y: any type that can be cast to a L{Data<escript.Data>} object on        :type Y: any type that can be cast to a `Data` object on
2729                 L{Function<escript.Function>}                 `Function`
2730        @keyword Y_reduced: value for coefficient C{Y_reduced}        :keyword Y_reduced: value for coefficient ``Y_reduced``
2731        @type Y_reduced: any type that can be cast to a L{Data<escript.Data>}        :type Y_reduced: any type that can be cast to a `Data`
2732                         object on L{ReducedFunction<escript.Function>}                         object on `ReducedFunction`
2733        @keyword d: value for coefficient C{d}        :keyword d: value for coefficient ``d``
2734        @type d: any type that can be cast to a L{Data<escript.Data>} object on        :type d: any type that can be cast to a `Data` object on
2735                 L{FunctionOnBoundary<escript.FunctionOnBoundary>}                 `FunctionOnBoundary`
2736        @keyword d_reduced: value for coefficient C{d_reduced}        :keyword d_reduced: value for coefficient ``d_reduced``
2737        @type d_reduced: any type that can be cast to a L{Data<escript.Data>}        :type d_reduced: any type that can be cast to a `Data`
2738                         object on L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}                         object on `ReducedFunctionOnBoundary`
2739        @keyword y: value for coefficient C{y}        :keyword y: value for coefficient ``y``
2740        @type y: any type that can be cast to a L{Data<escript.Data>} object on        :type y: any type that can be cast to a `Data` object on
2741                 L{FunctionOnBoundary<escript.FunctionOnBoundary>}                 `FunctionOnBoundary`
2742        @keyword d_contact: value for coefficient C{d_contact}        :keyword d_contact: value for coefficient ``d_contact``
2743        @type d_contact: any type that can be cast to a L{Data<escript.Data>}        :type d_contact: any type that can be cast to a `Data`
2744                         object on L{FunctionOnContactOne<escript.FunctionOnContactOne>}                         object on `FunctionOnContactOne`
2745                         or L{FunctionOnContactZero<escript.FunctionOnContactZero>}                         or `FunctionOnContactZero`
2746        @keyword d_contact_reduced: value for coefficient C{d_contact_reduced}        :keyword d_contact_reduced: value for coefficient ``d_contact_reduced``
2747        @type d_contact_reduced: any type that can be cast to a L{Data<escript.Data>}        :type d_contact_reduced: any type that can be cast to a `Data`
2748                                 object on L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}                                 object on `ReducedFunctionOnContactOne`
2749                                 or L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>}                                 or `ReducedFunctionOnContactZero`
2750        @keyword y_contact: value for coefficient C{y_contact}        :keyword y_contact: value for coefficient ``y_contact``
2751        @type y_contact: any type that can be cast to a L{Data<escript.Data>}        :type y_contact: any type that can be cast to a `Data`
2752                         object on L{FunctionOnContactOne<escript.FunctionOnContactOne>}                         object on `FunctionOnContactOne`
2753                         or L{FunctionOnContactZero<escript.FunctionOnContactZero>}                         or `FunctionOnContactZero`
2754        @keyword y_contact_reduced: value for coefficient C{y_contact_reduced}        :keyword y_contact_reduced: value for coefficient ``y_contact_reduced``
2755        @type y_contact_reduced: any type that can be cast to a L{Data<escript.Data>}        :type y_contact_reduced: any type that can be cast to a `Data`
2756                                 object on L{ReducedFunctionOnContactOne<escript.FunctionOnContactOne>}                                 object on `ReducedFunctionOnContactOne`
2757                                 or L{ReducedFunctionOnContactZero<escript.FunctionOnContactZero>}                                 or `ReducedFunctionOnContactZero`
2758        @keyword r: values prescribed to the solution at the locations of        :keyword r: values prescribed to the solution at the locations of
2759                    constraints                    constraints
2760        @type r: any type that can be cast to a L{Data<escript.Data>} object on        :type r: any type that can be cast to a `Data` object on
2761                 L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}                 `Solution` or `ReducedSolution`
2762                 depending on whether reduced order is used for the solution                 depending on whether reduced order is used for the solution
2763        @keyword q: mask for location of constraints        :keyword q: mask for location of constraints
2764        @type q: any type that can be cast to a L{Data<escript.Data>} object on        :type q: any type that can be cast to a `Data` object on
2765                 L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}                 `Solution` or `ReducedSolution`
2766                 depending on whether reduced order is used for the                 depending on whether reduced order is used for the
2767                 representation of the equation                 representation of the equation
2768        @raise IllegalCoefficient: if an unknown coefficient keyword is used        :raise IllegalCoefficient: if an unknown coefficient keyword is used
2769        """        """
2770        super(LinearPDE,self).setValue(**coefficients)        super(LinearPDE,self).setValue(**coefficients)
2771        # check if the systrem is inhomogeneous:        # check if the systrem is inhomogeneous:
# Line 2735  class LinearPDE(LinearProblem): Line 2782  class LinearPDE(LinearProblem):
2782       """       """
2783       Returns the residual of u or the current solution if u is not present.       Returns the residual of u or the current solution if u is not present.
2784    
2785       @param u: argument in the residual calculation. It must be representable       :param u: argument in the residual calculation. It must be representable
2786                 in L{self.getFunctionSpaceForSolution()}. If u is not present                 in `self.getFunctionSpaceForSolution()`. If u is not present
2787                 or equals C{None} the current solution is used.                 or equals ``None`` the current solution is used.
2788       @type u: L{Data<escript.Data>} or None       :type u: `Data` or None
2789       @return: residual of u       :return: residual of u
2790       @rtype: L{Data<escript.Data>}       :rtype: `Data`
2791       """       """
2792       if u==None:       if u==None:
2793          return self.getOperator()*self.getSolution()-self.getRightHandSide()          return self.getOperator()*self.getSolution()-self.getRightHandSide()
# Line 2749  class LinearPDE(LinearProblem): Line 2796  class LinearPDE(LinearProblem):
2796    
2797     def getFlux(self,u=None):     def getFlux(self,u=None):
2798       """       """
2799       Returns the flux M{J} for a given M{u}.       Returns the flux *J* for a given *u*.
2800    
2801       M{J[i,j]=(A[i,j,k,l]+A_reduced[A[i,j,k,l]]*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])u[k]-X[i,j]-X_reduced[i,j]}       *J[i,j]=(A[i,j,k,l]+A_reduced[A[i,j,k,l]]*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])u[k]-X[i,j]-X_reduced[i,j]*
2802    
2803       or       or
2804    
2805       M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[l]+(B[j]+B_reduced[j])u-X[j]-X_reduced[j]}       *J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[l]+(B[j]+B_reduced[j])u-X[j]-X_reduced[j]*
2806    
2807       @param u: argument in the flux. If u is not present or equals L{None} the       :param u: argument in the flux. If u is not present or equals `None` the
2808                 current solution is used.                 current solution is used.
2809       @type u: L{Data<escript.Data>} or None       :type u: `Data` or None
2810       @return: flux       :return: flux
2811       @rtype: L{Data<escript.Data>}       :rtype: `Data`
2812       """       """
2813       if u==None: u=self.getSolution()       if u==None: u=self.getSolution()
2814       return util.tensormult(self.getCoefficient("A"),util.grad(u,Funtion(self.getDomain))) \       return util.tensormult(self.getCoefficient("A"),util.grad(u,Funtion(self.getDomain))) \
# Line 2775  class LinearPDE(LinearProblem): Line 2822  class LinearPDE(LinearProblem):
2822  class Poisson(LinearPDE):  class Poisson(LinearPDE):
2823     """     """
2824     Class to define a Poisson equation problem. This is generally a     Class to define a Poisson equation problem. This is generally a
2825     L{LinearPDE} of the form     `LinearPDE` of the form
2826    
2827     M{-grad(grad(u)[j])[j] = f}     *-grad(grad(u)[j])[j] = f*
2828    
2829     with natural boundary conditions     with natural boundary conditions
2830    
2831     M{n[j]*grad(u)[j] = 0 }     *n[j]*grad(u)[j] = 0*
2832    
2833     and constraints:     and constraints:
2834    
2835     M{u=0} where M{q>0}     *u=0* where *q>0*
2836    
2837     """     """
2838    
# Line 2793  class Poisson(LinearPDE): Line 2840  class Poisson(LinearPDE):
2840       """       """
2841       Initializes a new Poisson equation.       Initializes a new Poisson equation.
2842    
2843       @param domain: domain of the PDE       :param domain: domain of the PDE
2844       @type domain: L{Domain<escript.Domain>}       :type domain: `Domain`
2845       @param debug: if True debug information is printed       :param debug: if True debug information is printed
2846    
2847       """       """
2848       super(Poisson, self).__init__(domain,1,1,debug)       super(Poisson, self).__init__(domain,1,1,debug)
# Line 2808  class Poisson(LinearPDE): Line 2855  class Poisson(LinearPDE):
2855       """       """
2856       Sets new values to coefficients.       Sets new values to coefficients.
2857    
2858       @param coefficients: new values assigned to coefficients       :param coefficients: new values assigned to coefficients
2859       @keyword f: value for right hand side M{f}       :keyword f: value for right hand side *f*
2860       @type f: any type that can be cast to a L{Scalar<escript.Scalar>} object       :type f: any type that can be cast to a `Scalar` object
2861                on L{Function<escript.Function>}                on `Function`
2862       @keyword q: mask for location of constraints       :keyword q: mask for location of constraints
2863       @type q: any type that can be cast to a rank zero L{Data<escript.Data>}       :type q: any type that can be cast to a rank zero `Data`
2864                object on L{Solution<escript.Solution>} or                object on `Solution` or
2865                L{ReducedSolution<escript.ReducedSolution>} depending on whether                `ReducedSolution` depending on whether
2866                reduced order is used for the representation of the equation                reduced order is used for the representation of the equation
2867       @raise IllegalCoefficient: if an unknown coefficient keyword is used       :raise IllegalCoefficient: if an unknown coefficient keyword is used
2868       """       """
2869       super(Poisson, self).setValue(**coefficients)       super(Poisson, self).setValue(**coefficients)
2870    
2871    
2872     def getCoefficient(self,name):     def getCoefficient(self,name):
2873       """       """
2874       Returns the value of the coefficient C{name} of the general PDE.       Returns the value of the coefficient ``name`` of the general PDE.
2875    
2876       @param name: name of the coefficient requested       :param name: name of the coefficient requested
2877       @type name: C{string}       :type name: ``string``
2878       @return: the value of the coefficient C{name}       :return: the value of the coefficient ``name``
2879       @rtype: L{Data<escript.Data>}       :rtype: `Data`
2880       @raise IllegalCoefficient: invalid coefficient name       :raise IllegalCoefficient: invalid coefficient name
2881       @note: This method is called by the assembling routine to map the Poisson       :note: This method is called by the assembling routine to map the Poisson
2882              equation onto the general PDE.              equation onto the general PDE.
2883       """       """
2884       if name == "A" :       if name == "A" :
# Line 2846  class Poisson(LinearPDE): Line 2893  class Poisson(LinearPDE):
2893  class Helmholtz(LinearPDE):  class Helmholtz(LinearPDE):
2894     """     """
2895     Class to define a Helmholtz equation problem. This is generally a     Class to define a Helmholtz equation problem. This is generally a
2896     L{LinearPDE} of the form     `LinearPDE` of the form
2897    
2898     M{S{omega}*u - grad(k*grad(u)[j])[j] = f}     *omega*u - grad(k*grad(u)[j])[j] = f*
2899    
2900     with natural boundary conditions     with natural boundary conditions
2901    
2902     M{k*n[j]*grad(u)[j] = g- S{alpha}u }     *k*n[j]*grad(u)[j] = g- alphau*
2903    
2904     and constraints:     and constraints:
2905    
2906     M{u=r} where M{q>0}     *u=r* where *q>0*
2907    
2908     """     """
2909    
# Line 2864  class Helmholtz(LinearPDE): Line 2911  class Helmholtz(LinearPDE):
2911       """       """
2912       Initializes a new Helmholtz equation.       Initializes a new Helmholtz equation.
2913    
2914       @param domain: domain of the PDE       :param domain: domain of the PDE
2915       @type domain: L{Domain<escript.Domain>}       :type domain: `Domain`
2916       @param debug: if True debug information is printed       :param debug: if True debug information is printed
2917    
2918       """       """
2919       super(Helmholtz, self).__init__(domain,1,1,debug)       super(Helmholtz, self).__init__(domain,1,1,debug)
# Line 2884  class Helmholtz(LinearPDE): Line 2931  class Helmholtz(LinearPDE):
2931       """       """
2932       Sets new values to coefficients.       Sets new values to coefficients.
2933    
2934       @param coefficients: new values assigned to coefficients       :param coefficients: new values assigned to coefficients
2935       @keyword omega: value for coefficient M{S{omega}}       :keyword omega: value for coefficient *omega*
2936       @type omega: any type that can be cast to a L{Scalar<escript.Scalar>}       :type omega: any type that can be cast to a `Scalar`
2937                    object on L{Function<escript.Function>}                    object on `Function`
2938       @keyword k: value for coefficient M{k}       :keyword k: value for coefficient *k*
2939       @type k: any type that can be cast to a L{Scalar<escript.Scalar>} object       :type k: any type that can be cast to a `Scalar` object
2940                on L{Function<escript.Function>}                on `Function`
2941       @keyword f: value for right hand side M{f}       :keyword f: value for right hand side *f*
2942       @type f: any type that can be cast to a L{Scalar<escript.Scalar>} object       :type f: any type that can be cast to a `Scalar` object
2943                on L{Function<escript.Function>}                on `Function`
2944       @keyword alpha: value for right hand side M{S{alpha}}       :keyword alpha: value for right hand side *alpha*
2945       @type alpha: any type that can be cast to a L{Scalar<escript.Scalar>}       :type alpha: any type that can be cast to a `Scalar`
2946                    object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}                    object on `FunctionOnBoundary`
2947       @keyword g: value for right hand side M{g}       :keyword g: value for right hand side *g*
2948       @type g: any type that can be cast to a L{Scalar<escript.Scalar>} object       :type g: any type that can be cast to a `Scalar` object
2949                on L{FunctionOnBoundary<escript.FunctionOnBoundary>}                on `FunctionOnBoundary`
2950       @keyword r: prescribed values M{r} for the solution in constraints       :keyword r: prescribed values *r* for the solution in constraints
2951       @type r: any type that can be cast to a L{Scalar<escript.Scalar>} object       :type r: any type that can be cast to a `Scalar` object
2952                on L{Solution<escript.Solution>} or                on `Solution` or
2953                L{ReducedSolution<escript.ReducedSolution>} depending on whether                `ReducedSolution` depending on whether
2954                reduced order is used for the representation of the equation                reduced order is used for the representation of the equation
2955       @keyword q: mask for the location of constraints       :keyword q: mask for the location of constraints
2956       @type q: any type that can be cast to a L{Scalar<escript.Scalar>} object       :type q: any type that can be cast to a `Scalar` object
2957                on L{Solution<escript.Solution>} or                on `Solution` or
2958                L{ReducedSolution<escript.ReducedSolution>} depending on whether                `ReducedSolution` depending on whether
2959                reduced order is used for the representation of the equation                reduced order is used for the representation of the equation
2960       @raise IllegalCoefficient: if an unknown coefficient keyword is used       :raise IllegalCoefficient: if an unknown coefficient keyword is used
2961       """       """
2962       super(Helmholtz, self).setValue(**coefficients)       super(Helmholtz, self).setValue(**coefficients)
2963    
2964     def getCoefficient(self,name):     def getCoefficient(self,name):
2965       """       """
2966       Returns the value of the coefficient C{name} of the general PDE.       Returns the value of the coefficient ``name`` of the general PDE.
2967    
2968       @param name: name of the coefficient requested       :param name: name of the coefficient requested
2969       @type name: C{string}       :type name: ``string``
2970       @return: the value of the coefficient C{name}       :return: the value of the coefficient ``name``
2971       @rtype: L{Data<escript.Data>}       :rtype: `Data`
2972       @raise IllegalCoefficient: invalid name       :raise IllegalCoefficient: invalid name
2973       """       """
2974       if name == "A" :       if name == "A" :
2975           if self.getCoefficient("k").isEmpty():           if self.getCoefficient("k").isEmpty():
# Line 2948  class LameEquation(LinearPDE): Line 2995  class LameEquation(LinearPDE):
2995     """     """
2996     Class to define a Lame equation problem. This problem is defined as:     Class to define a Lame equation problem. This problem is defined as:
2997    
2998     M{-grad(S{mu}*(grad(u[i])[j]+grad(u[j])[i]))[j] - grad(S{lambda}*grad(u[k])[k])[j] = F_i -grad(S{sigma}[ij])[j] }     *-grad(mu*(grad(u[i])[j]+grad(u[j])[i]))[j] - grad(lambda*grad(u[k])[k])[j] = F_i -grad(sigma[ij])[j]*
2999    
3000     with natural boundary conditions:     with natural boundary conditions:
3001    
3002     M{n[j]*(S{mu}*(grad(u[i])[j]+grad(u[j])[i]) + S{lambda}*grad(u[k])[k]) = f_i +n[j]*S{sigma}[ij] }     *n[j]*(mu*(grad(u[i])[j]+grad(u[j])[i]) + lambda*grad(u[k])[k]) = f_i +n[j]*sigma[ij]*
3003    
3004     and constraints:     and constraints:
3005    
3006     M{u[i]=r[i]} where M{q[i]>0}     *u[i]=r[i]* where *q[i]>0*
3007    
3008     """     """
3009    
# Line 2964  class LameEquation(LinearPDE): Line 3011  class LameEquation(LinearPDE):
3011        """        """
3012        Initializes a new Lame equation.        Initializes a new Lame equation.
3013    
3014        @param domain: domain of the PDE        :param domain: domain of the PDE
3015        @type domain: L{Domain<escript.Domain>}        :type domain: `Domain`
3016        @param debug: if True debug information is printed        :param debug: if True debug information is printed
3017    
3018        """        """
3019        super(LameEquation, self).__init__(domain,\        super(LameEquation, self).__init__(domain,\
# Line 2982  class LameEquation(LinearPDE): Line 3029  class LameEquation(LinearPDE):
3029       """       """
3030       Sets new values to coefficients.       Sets new values to coefficients.
3031    
3032       @param coefficients: new values assigned to coefficients       :param coefficients: new values assigned to coefficients
3033       @keyword lame_mu: value for coefficient M{S{mu}}       :keyword lame_mu: value for coefficient *mu*
3034       @type lame_mu: any type that can be cast to a L{Scalar<escript.Scalar>}       :type lame_mu: any type that can be cast to a `Scalar`
3035                      object on L{Function<escript.Function>}                      object on `Function`
3036       @keyword lame_lambda: value for coefficient M{S{lambda}}       :keyword lame_lambda: value for coefficient *lambda*
3037       @type lame_lambda: any type that can be cast to a L{Scalar<escript.Scalar>}       :type lame_lambda: any type that can be cast to a `Scalar`
3038                          object on L{Function<escript.Function>}                          object on `Function`
3039       @keyword F: value for internal force M{F}       :keyword F: value for internal force *F*
3040       @type F: any type that can be cast to a L{Vector<escript.Vector>} object       :type F: any type that can be cast to a `Vector` object
3041                on L{Function<escript.Function>}                on `Function`
3042       @keyword sigma: value for initial stress M{S{sigma}}       :keyword sigma: value for initial stress *sigma*
3043       @type sigma: any type that can be cast to a L{Tensor<escript.Tensor>}       :type sigma: any type that can be cast to a `Tensor`
3044                    object on L{Function<escript.Function>}                    object on `Function`
3045       @keyword f: value for external force M{f}       :keyword f: value for external force *f*
3046       @type f: any type that can be cast to a L{Vector<escript.Vector>} object       :type f: any type that can be cast to a `Vector` object
3047                on L{FunctionOnBoundary<escript.FunctionOnBoundary>}                on `FunctionOnBoundary`
3048       @keyword r: prescribed values M{r} for the solution in constraints       :keyword r: prescribed values *r* for the solution in constraints
3049       @type r: any type that can be cast to a L{Vector<escript.Vector>} object       :type r: any type that can be cast to a `Vector` object
3050                on L{Solution<escript.Solution>} or                on `Solution` or
3051                L{ReducedSolution<escript.ReducedSolution>} depending on whether                `ReducedSolution` depending on whether
3052                reduced order is used for the representation of the equation                reduced order is used for the representation of the equation
3053       @keyword q: mask for the location of constraints       :keyword q: mask for the location of constraints
3054       @type q: any type that can be cast to a L{Vector<escript.Vector>} object       :type q: any type that can be cast to a `Vector` object
3055                on L{Solution<escript.Solution>} or                on `Solution` or
3056                L{ReducedSolution<escript.ReducedSolution>} depending on whether                `ReducedSolution` depending on whether
3057                reduced order is used for the representation of the equation                reduced order is used for the representation of the equation
3058       @raise IllegalCoefficient: if an unknown coefficient keyword is used       :raise IllegalCoefficient: if an unknown coefficient keyword is used
3059       """       """
3060       super(LameEquation, self).setValues(**coefficients)       super(LameEquation, self).setValues(**coefficients)
3061    
3062     def getCoefficient(self,name):     def getCoefficient(self,name):
3063       """       """
3064       Returns the value of the coefficient C{name} of the general PDE.       Returns the value of the coefficient ``name`` of the general PDE.
3065    
3066       @param name: name of the coefficient requested       :param name: name of the coefficient requested
3067       @type name: C{string}       :type name: ``string``
3068       @return: the value of the coefficient C{name}       :return: the value of the coefficient ``name``
3069       @rtype: L{Data<escript.Data>}       :rtype: `Data`
3070       @raise IllegalCoefficient: invalid coefficient name       :raise IllegalCoefficient: invalid coefficient name
3071       """       """
3072       out =self.createCoefficient("A")       out =self.createCoefficient("A")
3073       if name == "A" :       if name == "A" :
# Line 3053  class LameEquation(LinearPDE): Line 3100  class LameEquation(LinearPDE):
3100       else:       else:
3101          return super(LameEquation, self).getCoefficient(name)          return super(LameEquation, self).getCoefficient(name)
3102    
3103    
3104  def LinearSinglePDE(domain,debug=False):  def LinearSinglePDE(domain,debug=False):
3105     """     """
3106     Defines a single linear PDE.     Defines a single linear PDE.
3107    
3108     @param domain: domain of the PDE     :param domain: domain of the PDE
3109     @type domain: L{Domain<escript.Domain>}     :type domain: `Domain`
3110     @param debug: if True debug information is printed     :param debug: if True debug information is printed
3111     @rtype: L{LinearPDE}     :rtype: `LinearPDE`
3112     """     """
3113     return LinearPDE(domain,numEquations=1,numSolutions=1,debug=debug)     return LinearPDE(domain,numEquations=1,numSolutions=1,debug=debug)
3114    
# Line 3068  def LinearPDESystem(domain,debug=False): Line 3116  def LinearPDESystem(domain,debug=False):
3116     """     """
3117     Defines a system of linear PDEs.     Defines a system of linear PDEs.
3118    
3119     @param domain: domain of the PDEs     :param domain: domain of the PDEs
3120     @type domain: L{Domain<escript.Domain>}     :type domain: `Domain`
3121     @param debug: if True debug information is printed     :param debug: if True debug information is printed
3122     @rtype: L{LinearPDE}     :rtype: `LinearPDE`
3123     """     """
3124     return LinearPDE(domain,numEquations=domain.getDim(),numSolutions=domain.getDim(),debug=debug)     return LinearPDE(domain,numEquations=domain.getDim(),numSolutions=domain.getDim(),debug=debug)
3125    
# Line 3080  class TransportPDE(LinearProblem): Line 3128  class TransportPDE(LinearProblem):
3128     """     """
3129     This class is used to define a transport problem given by a general linear,     This class is used to define a transport problem given by a general linear,
3130     time dependent, second order PDE for an unknown, non-negative,     time dependent, second order PDE for an unknown, non-negative,
3131     time-dependent function M{u} on a given domain defined through a     time-dependent function *u* on a given domain defined through a
3132     L{Domain<escript.Domain>} object.     `Domain` object.
3133    
3134     For a single equation with a solution with a single component the transport     For a single equation with a solution with a single component the transport
3135     problem is defined in the following form:     problem is defined in the following form:
3136    
3137     M{(M+M_reduced)*u_t=-(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)}     *(M+M_reduced)*u_t=-(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)*
3138    
3139     where M{u_t} denotes the time derivative of M{u} and M{grad(F)} denotes the     where *u_t* denotes the time derivative of *u* and *grad(F)* denotes the
3140     spatial derivative of M{F}. Einstein's summation convention,  ie. summation     spatial derivative of *F*. Einstein's summation convention,  ie. summation
3141     over indexes appearing twice in a term of a sum performed, is used.     over indexes appearing twice in a term of a sum performed, is used.
3142     The coefficients M{M}, M{A}, M{B}, M{C}, M{D}, M{X} and M{Y} have to be     The coefficients *M*, *A*, *B*, *C*, *D*, *X* and *Y* have to be
3143     specified through L{Data<escript.Data>} objects in L{Function<escript.Function>}     specified through `Data` objects in `Function`
3144     and the coefficients M{M_reduced}, M{A_reduced}, M{B_reduced}, M{C_reduced},     and the coefficients *M_reduced*, *A_reduced*, *B_reduced*, *C_reduced*,
3145     M{D_reduced}, M{X_reduced} and M{Y_reduced} have to be specified through     *D_reduced*, *X_reduced* and *Y_reduced* have to be specified through
3146     L{Data<escript.Data>} objects in L{ReducedFunction<escript.ReducedFunction>}.     `Data` objects in `ReducedFunction`.
3147     It is also allowed to use objects that can be converted into such     It is also allowed to use objects that can be converted into such
3148     L{Data<escript.Data>} objects. M{M} and M{M_reduced} are scalar, M{A} and     `Data` objects. *M* and *M_reduced* are scalar, *A* and
3149     M{A_reduced} are rank two, M{B}, M{C}, M{X}, M{B_reduced}, M{C_reduced} and     *A_reduced* are rank two, *B*, *C*, *X*, *B_reduced*, *C_reduced* and
3150     M{X_reduced} are rank one and M{D}, M{D_reduced}, M{Y} and M{Y_reduced} are     *X_reduced* are rank one and *D*, *D_reduced*, *Y* and *Y_reduced* are
3151     scalar.     scalar.
3152    
3153     The following natural boundary conditions are considered:     The following natural boundary conditions are considered:
3154    
3155     M{n[j]*((A[i,j]+A_reduced[i,j])*grad(u)[l]+(B+B_reduced)[j]*u+X[j]+X_reduced[j])+(d+d_reduced)*u+y+y_reduced=(m+m_reduced)*u_t}     *n[j]*((A[i,j]+A_reduced[i,j])*grad(u)[l]+(B+B_reduced)[j]*u+X[j]+X_reduced[j])+(d+d_reduced)*u+y+y_reduced=(m+m_reduced)*u_t*
3156    
3157     where M{n} is the outer normal field. Notice that the coefficients M{A},     where *n* is the outer normal field. Notice that the coefficients *A*,
3158     M{A_reduced}, M{B}, M{B_reduced}, M{X} and M{X_reduced} are defined in the     *A_reduced*, *B*, *B_reduced*, *X* and *X_reduced* are defined in the
3159     transport problem. The coefficients M{m}, M{d} and M{y} are each a scalar in     transport problem. The coefficients *m*, *d* and *y* are each a scalar in
3160     L{FunctionOnBoundary<escript.FunctionOnBoundary>} and the coefficients     `FunctionOnBoundary` and the coefficients
3161     M{m_reduced}, M{d_reduced} and M{y_reduced} are each a scalar in     *m_reduced*, *d_reduced* and *y_reduced* are each a scalar in
3162     L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.     `ReducedFunctionOnBoundary`.
3163    
3164     Constraints for the solution prescribing the value of the solution at     Constraints for the solution prescribing the value of the solution at
3165     certain locations in the domain have the form     certain locations in the domain have the form
3166    
3167     M{u_t=r} where M{q>0}     *u_t=r* where *q>0*
3168    
3169     M{r} and M{q} are each scalar where M{q} is the characteristic function     *r* and *q* are each scalar where *q* is the characteristic function
3170     defining where the constraint is applied. The constraints override any other     defining where the constraint is applied. The constraints override any other
3171     condition set by the transport problem or the boundary condition.     condition set by the transport problem or the boundary condition.
3172    
3173     The transport problem is symmetrical if     The transport problem is symmetrical if
3174    
3175     M{A[i,j]=A[j,i]} and M{B[j]=C[j]} and M{A_reduced[i,j]=A_reduced[j,i]} and     *A[i,j]=A[j,i]* and *B[j]=C[j]* and *A_reduced[i,j]=A_reduced[j,i]* and
3176     M{B_reduced[j]=C_reduced[j]}     *B_reduced[j]=C_reduced[j]*
3177    
3178     For a system and a solution with several components the transport problem     For a system and a solution with several components the transport problem
3179     has the form     has the form
3180    
3181     M{(M[i,k]+M_reduced[i,k])*u[k]_t=-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] }     *(M[i,k]+M_reduced[i,k]) * u[k]_t=-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]*
3182    
3183     M{A} and M{A_reduced} are of rank four, M{B}, M{B_reduced}, M{C} and     *A* and *A_reduced* are of rank four, *B*, *B_reduced*, *C* and
3184     M{C_reduced} are each of rank three, M{M}, M{M_reduced}, M{D}, M{D_reduced},     *C_reduced* are each of rank three, *M*, *M_reduced*, *D*, *D_reduced*,
3185     M{X_reduced} and M{X} are each of rank two and M{Y} and M{Y_reduced} are of     *X_reduced* and *X* are each of rank two and *Y* and *Y_reduced* are of
3186     rank one. The natural boundary conditions take the form:     rank one. The natural boundary conditions take the form:
3187    
3188     M{n[j]*((A[i,j,k,l]+A_reduced[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k]+X[i,j]+X_reduced[i,j])+(d[i,k]+d_reduced[i,k])*u[k]+y[i]+y_reduced[i]= (m[i,k]+m_reduced[i,k])*u[k]_t}     *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]+X[i,j]+X_reduced[i,j])+(d[i,k]+d_reduced[i,k])*u[k]+y[i]+y_reduced[i]= (m[i,k]+m_reduced[i,k])*u[k]_t*
3189    
3190     The coefficient M{d} and M{m} are of rank two and M{y} is of rank one with     The coefficient *d* and *m* are of rank two and *y* is of rank one with
3191     all in L{FunctionOnBoundary<escript.FunctionOnBoundary>}. The coefficients     all in `FunctionOnBoundary`. The coefficients
3192     M{d_reduced} and M{m_reduced} are of rank two and M{y_reduced} is of rank     *d_reduced* and *m_reduced* are of rank two and *y_reduced* is of rank
3193     one all in L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.     one all in `ReducedFunctionOnBoundary`.
3194    
3195     Constraints take the form     Constraints take the form
3196    
3197     M{u[i]_t=r[i]} where M{q[i]>0}     *u[i]_t=r[i]* where *q[i]>0*
3198    
3199     M{r} and M{q} are each rank one. Notice that at some locations not     *r* and *q* are each rank one. Notice that at some locations not
3200     necessarily all components must have a constraint.     necessarily all components must have a constraint.
3201    
3202     The transport problem is symmetrical if     The transport problem is symmetrical if
3203    
3204        - M{M[i,k]=M[i,k]}        - *M[i,k]=M[i,k]*
3205        - M{M_reduced[i,k]=M_reduced[i,k]}        - *M_reduced[i,k]=M_reduced[i,k]*
3206        - M{A[i,j,k,l]=A[k,l,i,j]}        - *A[i,j,k,l]=A[k,l,i,j]*
3207        - M{A_reduced[i,j,k,l]=A_reduced[k,l,i,j]}        - *A_reduced[i,j,k,l]=A_reduced[k,l,i,j]*
3208        - M{B[i,j,k]=C[k,i,j]}        - *B[i,j,k]=C[k,i,j]*
3209        - M{B_reduced[i,j,k]=C_reduced[k,i,j]}        - *B_reduced[i,j,k]=C_reduced[k,i,j]*
3210        - M{D[i,k]=D[i,k]}        - *D[i,k]=D[i,k]*
3211        - M{D_reduced[i,k]=D_reduced[i,k]}        - *D_reduced[i,k]=D_reduced[i,k]*
3212        - M{m[i,k]=m[k,i]}        - *m[i,k]=m[k,i]*
3213        - M{m_reduced[i,k]=m_reduced[k,i]}        - *m_reduced[i,k]=m_reduced[k,i]*
3214        - M{d[i,k]=d[k,i]}        - *d[i,k]=d[k,i]*
3215        - M{d_reduced[i,k]=d_reduced[k,i]}        - *d_reduced[i,k]=d_reduced[k,i]*
3216    
3217     L{TransportPDE} also supports solution discontinuities over a contact region     `TransportPDE` also supports solution discontinuities over a contact region
3218     in the domain. To specify the conditions across the discontinuity we are     in the domain. To specify the conditions across the discontinuity we are
3219     using the generalised flux M{J} which, in the case of a system of PDEs and     using the generalised flux *J* which, in the case of a system of PDEs and
3220     several components of the solution, is defined as     several components of the solution, is defined as
3221    
3222     M{J[i,j]=(A[i,j,k,l]+A_reduced[[i,j,k,l])*grad(u[k])[l]+(B[i,j,k]+B_reduced[i,j,k])*u[k]+X[i,j]+X_reduced[i,j]}     *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]*
3223    
3224     For the case of single solution component and single PDE M{J} is defined as     For the case of single solution component and single PDE *J* is defined as
3225    
3226     M{J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[j]+(B[i]+B_reduced[i])*u+X[i]+X_reduced[i]}     *J[j]=(A[i,j]+A_reduced[i,j])*grad(u)[j]+(B[i]+B_reduced[i])*u+X[i]+X_reduced[i]*
3227    
3228     In the context of discontinuities M{n} denotes the normal on the     In the context of discontinuities *n* denotes the normal on the
3229     discontinuity pointing from side 0 towards side 1 calculated from     discontinuity pointing from side 0 towards side 1 calculated from
3230     L{getNormal<escript.FunctionSpace.getNormal>} of L{FunctionOnContactZero<escript.FunctionOnContactZero>}.     `FunctionSpace.getNormal` of `FunctionOnContactZero`.
3231     For a system of transport problems the contact condition takes the form     For a system of transport problems the contact condition takes the form
3232    
3233     M{n[j]*J0[i,j]=n[j]*J1[i,j]=(y_contact[i]+y_contact_reduced[i])- (d_contact[i,k]+d_contact_reduced[i,k])*jump(u)[k]}     *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]*
3234    
3235     where M{J0} and M{J1} are the fluxes on side 0 and side 1 of the     where *J0* and *J1* are the fluxes on side 0 and side 1 of the
3236     discontinuity, respectively. M{jump(u)}, which is the difference of the     discontinuity, respectively. *jump(u)*, which is the difference of the
3237     solution at side 1 and at side 0, denotes the jump of M{u} across     solution at side 1 and at side 0, denotes the jump of *u* across
3238     discontinuity along the normal calculated by L{jump<util.jump>}.     discontinuity along the normal calculated by `jump`.
3239     The coefficient M{d_contact} is of rank two and M{y_contact} is of rank one     The coefficient *d_contact* is of rank two and *y_contact* is of rank one
3240     both in L{FunctionOnContactZero<escript.FunctionOnContactZero>} or L{FunctionOnContactOne<escript.FunctionOnContactOne>}.     both in `FunctionOnContactZero` or `FunctionOnContactOne`.
3241     The coefficient M{d_contact_reduced} is of rank two and M{y_contact_reduced}     The coefficient *d_contact_reduced* is of rank two and *y_contact_reduced*
3242     is of rank one both in L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.     is of rank one both in `ReducedFunctionOnContactZero` or `ReducedFunctionOnContactOne`.
3243     In case of a single PDE and a single component solution the contact     In case of a single PDE and a single component solution the contact
3244     condition takes the form     condition takes the form
3245    
3246     M{n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)}     *n[j]*J0_{j}=n[j]*J1_{j}=(y_contact+y_contact_reduced)-(d_contact+y_contact_reduced)*jump(u)*
3247    
3248     In this case the coefficient M{d_contact} and M{y_contact} are each scalar     In this case the coefficient *d_contact* and *y_contact* are each scalar
3249     both in L{FunctionOnContactZero<escript.FunctionOnContactZero>} or     both in `FunctionOnContactZero` or
3250     L{FunctionOnContactOne<escript.FunctionOnContactOne>} and the coefficient     `FunctionOnContactOne` and the coefficient
3251     M{d_contact_reduced} and M{y_contact_reduced} are each scalar both in     *d_contact_reduced* and *y_contact_reduced* are each scalar both in
3252     L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or     `ReducedFunctionOnContactZero` or
3253     L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.     `ReducedFunctionOnContactOne`.
3254    
3255     Typical usage::     Typical usage::
3256    
# Line 3219  class TransportPDE(LinearProblem): Line 3267  class TransportPDE(LinearProblem):
3267       """       """
3268       Initializes a transport problem.       Initializes a transport problem.
3269    
3270       @param domain: domain of the PDE       :param domain: domain of the PDE
3271       @type domain: L{Domain<escript.Domain>}       :type domain: `Domain`
3272       @param numEquations: number of equations. If C{None} the number of       :param numEquations: number of equations. If ``None`` the number of
3273                            equations is extracted from the coefficients.                            equations is extracted from the coefficients.
3274       @param numSolutions: number of solution components. If C{None} the number       :param numSolutions: number of solution components. If ``None`` the number
3275                            of solution components is extracted from the                            of solution components is extracted from the
3276                            coefficients.                            coefficients.
3277       @param debug: if True debug information is printed       :param debug: if True debug information is printed
3278       @param useBackwardEuler: if set the backward Euler scheme is used. Otherwise the Crank-Nicholson scheme is applied. Note that backward Euler scheme will return a safe time step size which is practically infinity as the scheme is unconditional unstable. The Crank-Nicholson scheme provides a higher accuracy but requires to limit the time step size to be stable.       :param useBackwardEuler: if set the backward Euler scheme is used. Otherwise the Crank-Nicolson scheme is applied. Note that backward Euler scheme will return a safe time step size which is practically infinity as the scheme is unconditional unstable. The Crank-Nicolson scheme provides a higher accuracy but requires to limit the time step size to be stable.
3279       @type useBackwardEuler: C{bool}       :type useBackwardEuler: ``bool``
3280       """       """
3281       if useBackwardEuler:       if useBackwardEuler:
3282           self.__useBackwardEuler=True           self.__useBackwardEuler=True
# Line 3272  class TransportPDE(LinearProblem): Line 3320  class TransportPDE(LinearProblem):
3320       """       """
3321       Returns the string representation of the transport problem.       Returns the string representation of the transport problem.
3322    
3323       @return: a simple representation of the transport problem       :return: a simple representation of the transport problem
3324       @rtype: C{str}       :rtype: ``str``
3325       """       """
3326       return "<TransportPDE %d>"%id(self)       return "<TransportPDE %d>"%id(self)
3327    
3328     def useBackwardEuler(self):     def useBackwardEuler(self):
3329        """        """
3330        Returns true if backward Euler is used. Otherwise false is returned.        Returns true if backward Euler is used. Otherwise false is returned.
3331        @rtype: bool        :rtype: bool
3332        """        """
3333        return self.__useBackwardEuler        return self.__useBackwardEuler
3334    
# Line 3289  class TransportPDE(LinearProblem): Line 3337  class TransportPDE(LinearProblem):
3337        """        """
3338        Tests the transport problem for symmetry.        Tests the transport problem for symmetry.
3339    
3340        @param verbose: if set to True or not present a report on coefficients        :param verbose: if set to True or not present a report on coefficients
3341                        which break the symmetry is printed.                        which break the symmetry is printed.
3342        @type verbose: C{bool}        :type verbose: ``bool``
3343        @return:  True if the PDE is symmetric        :return:  True if the PDE is symmetric
3344        @rtype: C{bool}        :rtype: ``bool``
3345        @note: This is a very expensive operation. It should be used for        :note: This is a very expensive operation. It should be used for
3346               degugging only! The symmetry flag is not altered.               degugging only! The symmetry flag is not altered.
3347        """        """
3348        out=True        out=True
# Line 3318  class TransportPDE(LinearProblem): Line 3366  class TransportPDE(LinearProblem):
3366        """        """
3367        Sets new values to coefficients.        Sets new values to coefficients.
3368    
3369        @param coefficients: new values assigned to coefficients        :param coefficients: new values assigned to coefficients
3370        @keyword M: value for coefficient C{M}        :keyword M: value for coefficient ``M``
3371        @type M: any type that can be cast to a L{Data<escript.Data>} object on        :type M: any type that can be cast to a `Data` object on
3372                 L{Function<escript.Function>}                 `Function`
3373        @keyword M_reduced: value for coefficient C{M_reduced}        :keyword M_reduced: value for coefficient ``M_reduced``
3374        @type M_reduced: any type that can be cast to a L{Data<escript.Data>}        :type M_reduced: any type that can be cast to a `Data`
3375                         object on L{Function<escript.ReducedFunction>}                         object on `Function`
3376        @keyword A: value for coefficient C{A}        :keyword A: value for coefficient ``A``
3377        @type A: any type that can be cast to a L{Data<escript.Data>} object on        :type A: any type that can be cast to a `Data` object on
3378                 L{Function<escript.Function>}                 `Function`
3379        @keyword A_reduced: value for coefficient C{A_reduced}        :keyword A_reduced: value for coefficient ``A_reduced``
3380        @type A_reduced: any type that can be cast to a L{Data<escript.Data>}        :type A_reduced: any type that can be cast to a `Data`
3381                         object on L{ReducedFunction<escript.ReducedFunction>}                         object on `ReducedFunction`
3382        @keyword B: value for coefficient C{B}        :keyword B: value for coefficient ``B``
3383        @type B: any type that can be cast to a L{Data<escript.Data>} object on        :type B: any type that can be cast to a `Data` object on
3384                 L{Function<escript.Function>}                 `Function`
3385        @keyword B_reduced: value for coefficient C{B_reduced}        :keyword B_reduced: value for coefficient ``B_reduced``
3386        @type B_reduced: any type that can be cast to a L{Data<escript.Data>}        :type B_reduced: any type that can be cast to a `Data`
3387                         object on L{ReducedFunction<escript.ReducedFunction>}                         object on `ReducedFunction`
3388        @keyword C: value for coefficient C{C}        :keyword C: value for coefficient ``C``
3389        @type C: any type that can be cast to a L{Data<escript.Data>} object on        :type C: any type that can be cast to a `Data` object on
3390                 L{Function<escript.Function>}                 `Function`
3391        @keyword C_reduced: value for coefficient C{C_reduced}        :keyword C_reduced: value for coefficient ``C_reduced``
3392        @type C_reduced: any type that can be cast to a L{Data<escript.Data>}        :type C_reduced: any type that can be cast to a `Data`
3393                         object on L{ReducedFunction<escript.ReducedFunction>}                         object on `ReducedFunction`
3394        @keyword D: value for coefficient C{D}        :keyword D: value for coefficient ``D``
3395        @type D: any type that can be cast to a L{Data<escript.Data>} object on        :type D: any type that can be cast to a `Data` object on
3396                 L{Function<escript.Function>}                 `Function`
3397        @keyword D_reduced: value for coefficient C{D_reduced}        :keyword D_reduced: value for coefficient ``D_reduced``
3398        @type D_reduced: any type that can be cast to a L{Data<escript.Data>}        :type D_reduced: any type that can be cast to a `Data`
3399                         object on L{ReducedFunction<escript.ReducedFunction>}                         object on `ReducedFunction`
3400        @keyword X: value for coefficient C{X}        :keyword X: value for coefficient ``X``
3401        @type X: any type that can be cast to a L{Data<escript.Data>} object on        :type X: any type that can be cast to a `Data` object on
3402                 L{Function<escript.Function>}                 `Function`
3403        @keyword X_reduced: value for coefficient C{X_reduced}        :keyword X_reduced: value for coefficient ``X_reduced``
3404        @type X_reduced: any type that can be cast to a L{Data<escript.Data>}        :type X_reduced: any type that can be cast to a `Data`
3405                         object on L{ReducedFunction<escript.ReducedFunction>}                         object on `ReducedFunction`
3406        @keyword Y: value for coefficient C{Y}        :keyword Y: value for coefficient ``Y``
3407        @type Y: any type that can be cast to a L{Data<escript.Data>} object on        :type Y: any type that can be cast to a `Data` object on
3408                 L{Function<escript.Function>}                 `Function`
3409        @keyword Y_reduced: value for coefficient C{Y_reduced}        :keyword Y_reduced: value for coefficient ``Y_reduced``
3410        @type Y_reduced: any type that can be cast to a L{Data<escript.Data>}        :type Y_reduced: any type that can be cast to a `Data`
3411                         object on L{ReducedFunction<escript.Function>}                         object on `ReducedFunction`
3412        @keyword m: value for coefficient C{m}        :keyword m: value for coefficient ``m``
3413        @type m: any type that can be cast to a L{Data<escript.Data>} object on        :type m: any type that can be cast to a `Data` object on
3414                 L{FunctionOnBoundary<escript.FunctionOnBoundary>}                 `FunctionOnBoundary`
3415        @keyword m_reduced: value for coefficient C{m_reduced}        :keyword m_reduced: value for coefficient ``m_reduced``
3416        @type m_reduced: any type that can be cast to a L{Data<escript.Data>}        :type m_reduced: any type that can be cast to a `Data`
3417                         object on L{FunctionOnBoundary<escript.ReducedFunctionOnBoundary>}                         object on `FunctionOnBoundary`
3418        @keyword d: value for coefficient C{d}        :keyword d: value for coefficient ``d``
3419        @type d: any type that can be cast to a L{Data<escript.Data>} object on        :type d: any type that can be cast to a `Data` object on
3420                 L{FunctionOnBoundary<escript.FunctionOnBoundary>}                 `FunctionOnBoundary`
3421        @keyword d_reduced: value for coefficient C{d_reduced}        :keyword d_reduced: value for coefficient ``d_reduced``
3422        @type d_reduced: any type that can be cast to a L{Data<escript.Data>}        :type d_reduced: any type that can be cast to a `Data`
3423                         object on L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}                         object on `ReducedFunctionOnBoundary`
3424        @keyword y: value for coefficient C{y}        :keyword y: value for coefficient ``y``
3425        @type y: any type that can be cast to a L{Data<escript.Data>} object on        :type y: any type that can be cast to a `Data` object on
3426                 L{FunctionOnBoundary<escript.FunctionOnBoundary>}                 `FunctionOnBoundary`
3427        @keyword d_contact: value for coefficient C{d_contact}        :keyword d_contact: value for coefficient ``d_contact``
3428        @type d_contact: any type that can be cast to a L{Data<escript.Data>}        :type d_contact: any type that can be cast to a `Data`
3429                         object on L{FunctionOnContactOne<escript.FunctionOnContactOne>} or L{FunctionOnContactZero<escript.FunctionOnContactZero>}                         object on `FunctionOnContactOne` or `FunctionOnContactZero`
3430        @keyword d_contact_reduced: value for coefficient C{d_contact_reduced}        :keyword d_contact_reduced: value for coefficient ``d_contact_reduced``
3431        @type d_contact_reduced: any type that can be cast to a L{Data<escript.Data>} object on L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>} or L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>}        :type d_contact_reduced: any type that can be cast to a `Data` object on `ReducedFunctionOnContactOne` or `ReducedFunctionOnContactZero`
3432        @keyword y_contact: value for coefficient C{y_contact}        :keyword y_contact: value for coefficient ``y_contact``
3433        @type y_contact: any type that can be cast to a L{Data<escript.Data>}        :type y_contact: any type that can be cast to a `Data`
3434                         object on L{FunctionOnContactOne<escript.FunctionOnContactOne>} or L{FunctionOnContactZero<escript.FunctionOnContactZero>}                         object on `FunctionOnContactOne` or `FunctionOnContactZero`
3435        @keyword y_contact_reduced: value for coefficient C{y_contact_reduced}        :keyword y_contact_reduced: value for coefficient ``y_contact_reduced``
3436        @type y_contact_reduced: any type that can be cast to a L{Data<escript.Data>} object on L{ReducedFunctionOnContactOne<escript.FunctionOnContactOne>} or L{ReducedFunctionOnContactZero<escript.FunctionOnContactZero>}        :type y_contact_reduced: any type that can be cast to a `Data` object on `ReducedFunctionOnContactOne` or `ReducedFunctionOnContactZero`
3437        @keyword r: values prescribed to the solution at the locations of constraints        :keyword r: values prescribed to the solution at the locations of constraints
3438        @type r: any type that can be cast to a L{Data<escript.Data>} object on        :type r: any type that can be cast to a `Data` object on
3439                 L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}                 `Solution` or `ReducedSolution`
3440                 depending on whether reduced order is used for the solution                 depending on whether reduced order is used for the solution
3441        @keyword q: mask for the location of constraints        :keyword q: mask for the location of constraints
3442        @type q: any type that can be cast to a L{Data<escript.Data>} object on        :type q: any type that can be cast to a `Data` object on
3443                 L{Solution<escript.Solution>} or                 `Solution` or
3444                 L{ReducedSolution<escript.ReducedSolution>} depending on whether                 `ReducedSolution` depending on whether
3445                 reduced order is used for the representation of the equation                 reduced order is used for the representation of the equation
3446        @raise IllegalCoefficient: if an unknown coefficient keyword is used        :raise IllegalCoefficient: if an unknown coefficient keyword is used
3447        """        """
3448        super(TransportPDE,self).setValue(**coefficients)        super(TransportPDE,self).setValue(**coefficients)
3449    
# Line 3403  class TransportPDE(LinearProblem): Line 3451  class TransportPDE(LinearProblem):
3451         """         """
3452         Returns an instance of a new transport operator.         Returns an instance of a new transport operator.
3453         """         """
        if self.useBackwardEuler():  
          theta=1.  
        else:  
          theta=0.5  
3454         optype=self.getRequiredOperatorType()         optype=self.getRequiredOperatorType()
3455         self.trace("New Transport problem pf type %s is allocated."%optype)         self.trace("New Transport problem pf type %s is allocated."%optype)
3456         return self.getDomain().newTransportProblem( \         return self.getDomain().newTransportProblem( \
3457                                 theta,                                 self.useBackwardEuler(),
3458                                 self.getNumEquations(), \                                 self.getNumEquations(), \
3459                                 self.getFunctionSpaceForSolution(), \                                 self.getFunctionSpaceForSolution(), \
3460                                 optype)                                 optype)
3461    
    def setInitialSolution(self,u):  
        """  
        Sets the initial solution.  
   
        @param u: new initial solution  
        @type u: any object that can be interpolated to a L{Data<escript.Data>}  
                 object on L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}  
        @note: C{u} must be non-negative  
        """  
        u2=util.interpolate(u,self.getFunctionSpaceForSolution())  
        if self.getNumSolutions() == 1:  
           if u2.getShape()!=():  
               raise ValueError,"Illegal shape %s of initial solution."%(u2.getShape(),)  
        else:  
           if u2.getShape()!=(self.getNumSolutions(),):  
               raise ValueError,"Illegal shape %s of initial solution."%(u2.getShape(),)  
        self.getOperator().setInitialValue(u2)  
3462    
3463     def getRequiredOperatorType(self):     def getRequiredOperatorType(self):
3464        """        """
3465        Returns the system type which needs to be used by the current set up.        Returns the system type which needs to be used by the current set up.
3466    
3467        @return: a code to indicate the type of transport problem scheme used        :return: a code to indicate the type of transport problem scheme used
3468        @rtype: C{float}        :rtype: ``float``
3469        """        """
3470        solver_options=self.getSolverOptions()        solver_options=self.getSolverOptions()
3471        return self.getDomain().getTransportTypeId(solver_options.getSolverMethod(), solver_options.getPreconditioner(),solver_options.getPackage(), solver_options.isSymmetric())        return self.getDomain().getTransportTypeId(solver_options.getSolverMethod(), solver_options.getPreconditioner(),solver_options.getPackage(), solver_options.isSymmetric())
3472    
3473     def getUnlimitedTimeStepSize(self):     def getUnlimitedTimeStepSize(self):
3474        """        """
3475        Returns the value returned by the C{getSafeTimeStepSize} method to        Returns the value returned by the ``getSafeTimeStepSize`` method to
3476        indicate no limit on the safe time step size.        indicate no limit on the safe time step size.
3477    
3478         @return: the value used to indicate that no limit is set to the time         :return: the value used to indicate that no limit is set to the time
3479                  step size                  step size
3480         @rtype: C{float}         :rtype: ``float``
3481         @note: Typically the biggest positive float is returned         :note: Typically the biggest positive float is returned
3482        """        """
3483        return self.getOperator().getUnlimitedTimeStepSize()        return self.getOperator().getUnlimitedTimeStepSize()
3484    
# Line 3459  class TransportPDE(LinearProblem): Line 3486  class TransportPDE(LinearProblem):
3486         """         """
3487         Returns a safe time step size to do the next time step.         Returns a safe time step size to do the next time step.
3488    
3489         @return: safe time step size         :return: safe time step size
3490         @rtype: C{float}         :rtype: ``float``
3491         @note: If not C{getSafeTimeStepSize()} < C{getUnlimitedTimeStepSize()}         :note: If not ``getSafeTimeStepSize()`` < ``getUnlimitedTimeStepSize()``
3492                any time step size can be used.                any time step size can be used.
3493         """         """
3494         return self.getOperator().getSafeTimeStepSize()         return self.getOperator().getSafeTimeStepSize()
# Line 3470  class TransportPDE(LinearProblem): Line 3497  class TransportPDE(LinearProblem):
3497         """         """
3498         Sets the weighting factor used to insert the constraints into the problem         Sets the weighting factor used to insert the constraints into the problem
3499    
3500         @param value: value for the weighting factor         :param value: value for the weighting factor
3501         @type value: large positive C{float}         :type value: large positive ``float``
3502         """         """
3503         if not value>0:         if not value>0:
3504           raise ValueError,"weighting factor needs to be positive."           raise ValueError,"weighting factor needs to be positive."
# Line 3481  class TransportPDE(LinearProblem): Line 3508  class TransportPDE(LinearProblem):
3508     def getConstraintWeightingFactor(self):     def getConstraintWeightingFactor(self):
3509         """         """
3510         returns the weighting factor used to insert the constraints into the problem         returns the weighting factor used to insert the constraints into the problem
3511         @return: value for the weighting factor         :return: value for the weighting factor
3512         @rtype: C{float}         :rtype: ``float``
3513         """         """
3514         return self.__constraint_factor         return self.__constraint_factor
3515     #====================================================================     #====================================================================
3516     def getSolution(self,dt):     def getSolution(self, dt=None, u0=None):
3517         """         """
3518         Returns the solution of the problem.         Returns the solution by marching forward by time step dt. if ''u0'' is present,
3519           ''u0'' is used as the initial value otherwise the solution from the last call is used.
3520    
3521         @return: the solution         :param dt: time step size. If ``None`` the last solution is returned.
3522         @rtype: L{Data<escript.Data>}         :type dt: positive ``float`` or ``None``
3523         """         :param u0: new initial solution or ``None``
3524         option_class=self.getSolverOptions()         :type u0: any object that can be interpolated to a `Data`
3525         if dt<=0:                  object on `Solution` or `ReducedSolution`
3526             raise ValueError,"step size needs to be positive."         :return: the solution
3527         self.setSolution(self.getOperator().solve(self.getRightHandSide(),dt,option_class))         :rtype: `Data`
3528         self.validSolution()         """
3529           if not dt == None:
3530          option_class=self.getSolverOptions()
3531          if dt<=0:
3532              raise ValueError,"step size needs to be positive."
3533          if u0 == None:
3534              u0=self.getCurrentSolution()
3535          else:
3536              u0=util.interpolate(u0,self.getFunctionSpaceForSolution())
3537              if self.getNumSolutions() == 1:
3538            if u0.getShape()!=():
3539                raise ValueError,"Illegal shape %s of initial solution."%(u0.getShape(),)
3540              else:
3541                if u0.getShape()!=(self.getNumSolutions(),):
3542                  raise ValueError,"Illegal shape %s of initial solution."%(u0.getShape(),)
3543          self.setSolution(self.getOperator().solve(u0, self.getRightHandSide(),dt,option_class))
3544          self.validSolution()
3545         return self.getCurrentSolution()         return self.getCurrentSolution()
3546    
3547       def setInitialSolution(self,u):
3548           """
3549           Sets the initial solution.
3550    
3551           :param u: initial solution
3552           :type u: any object that can be interpolated to a `Data`
3553                    object on `Solution` or `ReducedSolution`
3554           """
3555           u2=util.interpolate(u,self.getFunctionSpaceForSolution())
3556           if self.getNumSolutions() == 1:
3557              if u2.getShape()!=():
3558                  raise ValueError,"Illegal shape %s of initial solution."%(u2.getShape(),)
3559           else:
3560              if u2.getShape()!=(self.getNumSolutions(),):
3561                  raise ValueError,"Illegal shape %s of initial solution."%(u2.getShape(),)
3562           self.setSolution(u2,validate=False)
3563    
3564    
3565     def getSystem(self):     def getSystem(self):
3566         """         """
3567         Returns the operator and right hand side of the PDE.         Returns the operator and right hand side of the PDE.
3568    
3569         @return: the discrete version of the PDE         :return: the discrete version of the PDE
3570         @rtype: C{tuple} of L{Operator,<escript.Operator>} and         :rtype: ``tuple`` of `Operator` and
3571                 L{Data<escript.Data>}                 `Data`
3572    
3573         """         """
3574         if not self.isOperatorValid() or not self.isRightHandSideValid():         if not self.isOperatorValid() or not self.isRightHandSideValid():
# Line 3552  class TransportPDE(LinearProblem): Line 3614  class TransportPDE(LinearProblem):
3614    
3615     def setDebug(self, flag):     def setDebug(self, flag):
3616       """       """
3617       Switches debug output on if C{flag} is True,       Switches debug output on if ``flag`` is True,
3618       otherwise it is switched off.       otherwise it is switched off.
3619    
3620       @param flag: desired debug status       :param flag: desired debug status
3621       @type flag: C{bool}       :type flag: ``bool``
3622       """       """
3623       if flag:       if flag:
3624           self.setDebugOn()           self.setDebugOn()
# Line 3575  class TransportPDE(LinearProblem): Line 3637  class TransportPDE(LinearProblem):
3637       """       """
3638       super(TransportPDE,self).setDebugOff()       super(TransportPDE,self).setDebugOff()
3639            
3640    def SingleTransportPDE(domain,useBackwardEuler=False, debug=False):
3641       """
3642       Defines a single transport problem
3643    
3644       :param domain: domain of the PDE
3645       :type domain: `Domain`
3646       :param debug: if True debug information is printed
3647       :param useBackwardEuler: if set the backward Euler scheme is used. Otherwise the Crank-Nicolson scheme is applied. Note that backward Euler scheme will return a safe time step size which is practically infinity as the scheme is unconditional unstable. The Crank-Nicolson scheme provides a higher accuracy but requires to limit the time step size to be stable.
3648       :rtype: `TransportPDE`
3649       """
3650       return TransportPDE(domain,numEquations=1,numSolutions=1, useBackwardEuler=useBackwardEuler, debug=debug)

Legend:
Removed from v.2549  
changed lines
  Added in v.3071

  ViewVC Help
Powered by ViewVC 1.1.26