/[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 2987 by gross, Tue Mar 16 01:32:43 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: Transport 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 DEFAULT_REORDERING: the reordering method recommended by the solver
94      @cvar YAIR_SHAPIRA_COARSENING: AMG coarsening method by Yair-Shapira      :cvar SUPER_LU: the Super_LU solver package
95      @cvar RUGE_STUEBEN_COARSENING: AMG coarsening method by Ruge and Stueben      :cvar PASTIX: the Pastix direct solver_package
96      @cvar AGGREGATION_COARSENING: AMG coarsening using (symmetric) aggregation      :cvar YAIR_SHAPIRA_COARSENING: AMG and AMLI coarsening method by Yair-Shapira
97      @cvar MIN_COARSE_MATRIX_SIZE: minimum size of the coarsest level matrix to use direct solver.      :cvar RUGE_STUEBEN_COARSENING: AMG and AMLI coarsening method by Ruge and Stueben
98      @cvar NO_PRECONDITIONER: no preconditioner is applied.      :cvar AGGREGATION_COARSENING: AMG and AMLI coarsening using (symmetric) aggregation
99        :cvar STANDARD_COARSENING: AMG and AMLI standard coarsening using mesure of importance of the unknowns
100        :cvar MIN_COARSE_MATRIX_SIZE: minimum size of the coarsest level matrix to use direct solver.
101        :cvar NO_PRECONDITIONER: no preconditioner is applied.
102      """      """
103      DEFAULT= 0      DEFAULT= 0
104      DIRECT= 1      DIRECT= 1
# Line 134  class SolverOptions(object): Line 137  class SolverOptions(object):
137      AGGREGATION_COARSENING=35      AGGREGATION_COARSENING=35
138      NO_PRECONDITIONER=36      NO_PRECONDITIONER=36
139      MIN_COARSE_MATRIX_SIZE=37      MIN_COARSE_MATRIX_SIZE=37
140            AMLI=38
141        STANDARD_COARSENING=39
142    
143      def __init__(self):      def __init__(self):
144          self.setLevelMax()          self.setLevelMax()
145          self.setCoarseningThreshold()          self.setCoarseningThreshold()
146            self.setSmoother()
147          self.setNumSweeps()          self.setNumSweeps()
148          self.setNumPreSweeps()          self.setNumPreSweeps()
149          self.setNumPostSweeps()          self.setNumPostSweeps()
# Line 194  class SolverOptions(object): Line 200  class SolverOptions(object):
200                  out+="\nMaximum number of levels = %s"%self.LevelMax()                  out+="\nMaximum number of levels = %s"%self.LevelMax()
201                  out+="\nCoarsening method = %s"%self.getName(self.getCoarsening())                  out+="\nCoarsening method = %s"%self.getName(self.getCoarsening())
202                  out+="\nCoarsening threshold = %e"%self.getMinCoarseMatrixSize()                  out+="\nCoarsening threshold = %e"%self.getMinCoarseMatrixSize()
203                    out+="\nSmoother = %e"%self.getName(self.getSmoother())
204                    out+="\nMinimum size of the coarsest level matrix = %e"%self.getCoarseningThreshold()
205                    out+="\nNumber of pre / post sweeps = %s / %s, %s"%(self.getNumPreSweeps(), self.getNumPostSweeps(), self.getNumSweeps())
206                if self.getPreconditioner() == self.AMLI:
207                    out+="\nMaximum number of levels = %s"%self.LevelMax()
208                    out+="\nCoarsening method = %s"%self.getName(self.getCoarsening())
209                    out+="\nCoarsening threshold = %e"%self.getMinCoarseMatrixSize()
210                  out+="\nMinimum size of the coarsest level matrix = %e"%self.getCoarseningThreshold()                  out+="\nMinimum size of the coarsest level matrix = %e"%self.getCoarseningThreshold()
211                  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())
212              if self.getPreconditioner() == self.GAUSS_SEIDEL:          if self.getPreconditioner() == self.GAUSS_SEIDEL:
213                  out+="\nNumber of sweeps = %s"%self.getNumSweeps()                  out+="\nNumber of sweeps = %s"%self.getNumSweeps()
214              if self.getPreconditioner() == self.ILUT:              if self.getPreconditioner() == self.ILUT:
215                  out+="\nDrop tolerance = %e"%self.getDropTolerance()                  out+="\nDrop tolerance = %e"%self.getDropTolerance()
# Line 209  class SolverOptions(object): Line 222  class SolverOptions(object):
222          """          """
223          returns the name of a given key          returns the name of a given key
224                    
225          @param key: a valid key          :param key: a valid key
226          """          """
227          if key == self.DEFAULT: return "DEFAULT"          if key == self.DEFAULT: return "DEFAULT"
228          if key == self.DIRECT: return "DIRECT"          if key == self.DIRECT: return "DIRECT"
# Line 233  class SolverOptions(object): Line 246  class SolverOptions(object):
246          if key == self.ITERATIVE: return "ITERATIVE"          if key == self.ITERATIVE: return "ITERATIVE"
247          if key == self.PASO: return "PASO"          if key == self.PASO: return "PASO"
248          if key == self.AMG: return "AMG"          if key == self.AMG: return "AMG"
249            if key == self.AMLI: return "AMLI"
250          if key == self.REC_ILU: return "REC_ILU"          if key == self.REC_ILU: return "REC_ILU"
251          if key == self.TRILINOS: return "TRILINOS"          if key == self.TRILINOS: return "TRILINOS"
252          if key == self.NONLINEAR_GMRES: return "NONLINEAR_GMRES"          if key == self.NONLINEAR_GMRES: return "NONLINEAR_GMRES"
# Line 245  class SolverOptions(object): Line 259  class SolverOptions(object):
259          if key == self.PASTIX: return "PASTIX"          if key == self.PASTIX: return "PASTIX"
260          if key == self.YAIR_SHAPIRA_COARSENING: return "YAIR_SHAPIRA_COARSENING"          if key == self.YAIR_SHAPIRA_COARSENING: return "YAIR_SHAPIRA_COARSENING"
261          if key == self.RUGE_STUEBEN_COARSENING: return "RUGE_STUEBEN_COARSENING"          if key == self.RUGE_STUEBEN_COARSENING: return "RUGE_STUEBEN_COARSENING"
262            if key == self.STANDARD_COARSENING: return "STANDARD_COARSENING"
263          if key == self.AGGREGATION_COARSENING: return "AGGREGATION_COARSENING"          if key == self.AGGREGATION_COARSENING: return "AGGREGATION_COARSENING"
264          if key == self.NO_PRECONDITIONER: return "NO_PRECONDITIONER"          if key == self.NO_PRECONDITIONER: return "NO_PRECONDITIONER"
265          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 268  class SolverOptions(object):
268          """          """
269          resets the diagnostics          resets the diagnostics
270                    
271          @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.
272          @type all: C{bool}          :type all: ``bool``
273          """          """
274          self.__num_iter=None          self.__num_iter=None
275          self.__num_level=None          self.__num_level=None
276          self.__num_inner_iter=None          self.__num_inner_iter=None
277          self.__time=None          self.__time=None
278          self.__set_up_time=None          self.__set_up_time=None
279            self.__net_time=None
280          self.__residual_norm=None          self.__residual_norm=None
281          self.__converged=None          self.__converged=None
282          if all:          if all:
# Line 268  class SolverOptions(object): Line 284  class SolverOptions(object):
284              self.__cum_num_iter=0              self.__cum_num_iter=0
285              self.__cum_time=0              self.__cum_time=0
286              self.__cum_set_up_time=0              self.__cum_set_up_time=0
287                self.__cum_net_time=0
288    
289      def _updateDiagnostics(self, name, value):      def _updateDiagnostics(self, name, value):
290          """          """
291          Updates diagnostic information          Updates diagnostic information
292                    
293          @param name: name of  diagnostic information          :param name: name of  diagnostic information
294          @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".
295          @param vale: new value of the diagnostic information          :param value: new value of the diagnostic information
296          @note: this function is used by a solver to report diagnostics informations.          :note: this function is used by a solver to report diagnostics informations.
297          """          """
298          if name == "num_iter":          if name == "num_iter":
299              self.__num_iter=int(value)              self.__num_iter=int(value)
300              self.__cum_num_iter+=self.__num_iter              self.__cum_num_iter+=self.__num_iter
301          if name == "num_level":          if name == "num_level":
302              self.__num_iter=int(value)              self.__num_level=int(value)
303          if name == "num_inner_iter":          if name == "num_inner_iter":
304              self.__num_inner_iter=int(value)              self.__num_inner_iter=int(value)
305              self.__cum_num_inner_iter+=self.__num_inner_iter              self.__cum_num_inner_iter+=self.__num_inner_iter
# Line 292  class SolverOptions(object): Line 309  class SolverOptions(object):
309          if name == "set_up_time":          if name == "set_up_time":
310              self.__set_up_time=float(value)              self.__set_up_time=float(value)
311              self.__cum_set_up_time+=self.__set_up_time              self.__cum_set_up_time+=self.__set_up_time
312            if name == "net_time":
313                self.__net_time=float(value)
314                self.__cum_net_time+=self.__net_time
315          if name == "residual_norm":          if name == "residual_norm":
316              self.__residual_norm=float(value)              self.__residual_norm=float(value)
317          if name == "converged":          if name == "converged":
318              self.__converged = (value == True)              self.__converged = (value == True)
319      def getDiagnostics(self, name):      def getDiagnostics(self, name):
320          """          """
321          Returns the diagnostic information C{name}          Returns the diagnostic information ``name``. Possible values are:
322                    
323          @param name: name of diagnostic information where      - "num_iter": the number of iteration steps
         - "num_iter": the number of iteration steps  
324          - "cum_num_iter": the cumulative number of iteration steps          - "cum_num_iter": the cumulative number of iteration steps
325          - "num_level": the number of level in multi level solver          - "num_level": the number of level in multi level solver
326          - "num_inner_iter": the number of inner iteration steps          - "num_inner_iter": the number of inner iteration steps
# Line 310  class SolverOptions(object): Line 329  class SolverOptions(object):
329          - "cum_time": cumulative execution time          - "cum_time": cumulative execution time
330          - "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
331          - "cum_set_up_time": cumulative time to set up of the solver          - "cum_set_up_time": cumulative time to set up of the solver
332            - "net_time": net execution time, excluding setup time for the solver and execution time for preconditioner
333            - "cum_net_time": cumulative net execution time
334          - "residual_norm": norm of the final residual          - "residual_norm": norm of the final residual
335          - "converged": return self.__converged              - "converged": return self.__converged
336          @type name: C{str} in the list "num_iter", "num_level", "num_inner_iter", "time", "set_up_time", "residual_norm", "converged".      
337          @return: requested value. C{None} is returned if the value is undefined.      
338          @note: If the solver has thrown an exception diagnostic values have an undefined status.          
339            :param name: name of diagnostic information to return
340            :type name: ``str`` in the list above.
341            :return: requested value. ``None`` is returned if the value is undefined.
342            :note: If the solver has thrown an exception diagnostic values have an undefined status.
343          """          """
344          if name == "num_iter": return self.__num_iter          if name == "num_iter": return self.__num_iter
345          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 350  class SolverOptions(object):
350          elif name == "cum_time": return self.__cum_time          elif name == "cum_time": return self.__cum_time
351          elif name == "set_up_time": return self.__set_up_time          elif name == "set_up_time": return self.__set_up_time
352          elif name == "cum_set_up_time": return self.__cum_set_up_time          elif name == "cum_set_up_time": return self.__cum_set_up_time
353            elif name == "net_time": return self.__net_time
354            elif name == "cum_net_time": return self.__cum_net_time
355          elif name == "residual_norm": return self.__residual_norm          elif name == "residual_norm": return self.__residual_norm
356          elif name == "converged": return self.__converged                elif name == "converged": return self.__converged      
357          else:          else:
358              raise ValueError,"unknown diagnostic item %s"%name              raise ValueError,"unknown diagnostic item %s"%name
359      def hasConverged(self):      def hasConverged(self):
360          """          """
361          Returns C{True} if the last solver call has been finalized successfully.          Returns ``True`` if the last solver call has been finalized successfully.
362          @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.
363          """          """
364          return self.getDiagnostics("converged")          return self.getDiagnostics("converged")
365      def setCoarsening(self,method=0):      def setCoarsening(self,method=0):
366          """          """
367          Sets the key of the coarsening method to be applied in AMG.          Sets the key of the coarsening method to be applied in AMG.
368    
369          @param method: selects the coarsening method .          :param method: selects the coarsening method .
370          @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}  
371          """          """
372      if method==None: method=0      if method==None: method=0
373          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,]:
374               raise ValueError,"unknown coarsening method %s"%method               raise ValueError,"unknown coarsening method %s"%method
375          self.__coarsening=method          self.__coarsening=method
376            
# Line 352  class SolverOptions(object): Line 378  class SolverOptions(object):
378          """          """
379          Returns the key of the coarsening algorithm to be applied AMG.          Returns the key of the coarsening algorithm to be applied AMG.
380    
381          @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}  
382          """          """
383          return self.__coarsening          return self.__coarsening
384                
# Line 361  class SolverOptions(object): Line 386  class SolverOptions(object):
386          """          """
387          Sets the minumum size of the coarsest level matrix in AMG.          Sets the minumum size of the coarsest level matrix in AMG.
388    
389          @param size: minumum size of the coarsest level matrix .          :param size: minumum size of the coarsest level matrix .
390          @type size: positive C{int} or C{None}          :type size: positive ``int`` or ``None``
391          """          """
392          size=int(size)          size=int(size)
393          if size<0:          if size<0:
# Line 374  class SolverOptions(object): Line 399  class SolverOptions(object):
399          """          """
400          Returns the minumum size of the coarsest level matrix in AMG.          Returns the minumum size of the coarsest level matrix in AMG.
401    
402          @rtype: C{int}          :rtype: ``int``
403          """          """
404          return self.__MinCoarseMatrixSize          return self.__MinCoarseMatrixSize
405                
# Line 382  class SolverOptions(object): Line 407  class SolverOptions(object):
407          """          """
408          Sets the preconditioner to be used.          Sets the preconditioner to be used.
409    
410          @param preconditioner: key of the preconditioner to be used.          :param preconditioner: key of the preconditioner to be used.
411          @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`,
412                                      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.RILU`,
413                                      L{SolverOptions.NO_PRECONDITIONER}                                      `SolverOptions.NO_PRECONDITIONER`
414          @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.  
415          """          """
416      if preconditioner==None: preconditioner=10      if preconditioner==None: preconditioner=10
417          if not preconditioner in [ SolverOptions.SSOR, SolverOptions.ILU0, SolverOptions.ILUT, SolverOptions.JACOBI,          if not preconditioner in [ SolverOptions.SSOR, SolverOptions.ILU0, SolverOptions.ILUT, SolverOptions.JACOBI,
418                                      SolverOptions.AMG, SolverOptions.REC_ILU, SolverOptions.GAUSS_SEIDEL, SolverOptions.RILU,                                      SolverOptions.AMG, SolverOptions.AMLI, SolverOptions.REC_ILU, SolverOptions.GAUSS_SEIDEL, SolverOptions.RILU,
419                                      SolverOptions.NO_PRECONDITIONER] :                                      SolverOptions.NO_PRECONDITIONER] :
420               raise ValueError,"unknown preconditioner %s"%preconditioner               raise ValueError,"unknown preconditioner %s"%preconditioner
421          self.__preconditioner=preconditioner              self.__preconditioner=preconditioner    
# Line 399  class SolverOptions(object): Line 423  class SolverOptions(object):
423          """          """
424          Returns key of the preconditioner to be used.          Returns key of the preconditioner to be used.
425    
426          @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`,
427                                      L{SolverOptions.AMG}, L{SolverOptions.REC_ILU}, L{SolverOptions.GAUSS_SEIDEL}, L{SolverOptions.RILU},                                      `SolverOptions.AMG`, `SolverOptions.REC_ILU`, `SolverOptions.GAUSS_SEIDEL`, `SolverOptions.RILU`,
428                                      L{SolverOptions.NO_PRECONDITIONER}                                      `SolverOptions.NO_PRECONDITIONER`
429          """          """
430          return self.__preconditioner          return self.__preconditioner
431        def setSmoother(self, smoother=28):
432            """
433            Sets the smoother to be used.
434    
435            :param smoother: key of the smoother to be used.
436            :type smoother: in `SolverOptions.JACOBI`, `SolverOptions.GAUSS_SEIDEL`
437            :note: Not all packages support all smoothers. It can be assumed that a package makes a reasonable choice if it encounters an unknown smoother.
438            """
439        if smoother==None: smoother=28
440            if not smoother in [ SolverOptions.JACOBI, SolverOptions.GAUSS_SEIDEL ] :
441                 raise ValueError,"unknown smoother %s"%smoother
442            self.__smoother=smoother    
443        def getSmoother(self):
444            """
445            Returns key of the smoother to be used.
446    
447            :rtype: in the list `SolverOptions.JACOBI`, `SolverOptions.GAUSS_SEIDEL`
448            """
449            return self.__smoother  
450      def setSolverMethod(self, method=0):      def setSolverMethod(self, method=0):
451          """          """
452          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
453          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
454          solver should be used.          solver should be used.
455    
456          @param method: key of the solver method to be used.          :param method: key of the solver method to be used.
457          @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`,
458                          L{SolverOptions.CR}, L{SolverOptions.CGS}, L{SolverOptions.BICGSTAB}, L{SolverOptions.SSOR},                          `SolverOptions.CR`, `SolverOptions.CGS`, `SolverOptions.BICGSTAB`, `SolverOptions.SSOR`,
459                          L{SolverOptions.GMRES}, L{SolverOptions.PRES20}, L{SolverOptions.LUMPING}, L{SolverOptions.ITERATIVE},                          `SolverOptions.GMRES`, `SolverOptions.PRES20`, `SolverOptions.LUMPING`, `SolverOptions.ITERATIVE`,
460                          L{SolverOptions.AMG}, L{SolverOptions.NONLINEAR_GMRES}, L{SolverOptions.TFQMR}, L{SolverOptions.MINRES},                          `SolverOptions.AMG`, `SolverOptions.NONLINEAR_GMRES`, `SolverOptions.TFQMR`, `SolverOptions.MINRES`,
461                          L{SolverOptions.GAUSS_SEIDEL}                          `SolverOptions.GAUSS_SEIDEL`
462          @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.  
463          """          """
464      if method==None: method=0      if method==None: method=0
465          if not method in [ SolverOptions.DEFAULT, SolverOptions.DIRECT, SolverOptions.CHOLEVSKY, SolverOptions.PCG,          if not method in [ SolverOptions.DEFAULT, SolverOptions.DIRECT, SolverOptions.CHOLEVSKY, SolverOptions.PCG,
# Line 430  class SolverOptions(object): Line 472  class SolverOptions(object):
472          """          """
473          Returns key of the solver method to be used.          Returns key of the solver method to be used.
474    
475          @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`,
476                          L{SolverOptions.CR}, L{SolverOptions.CGS}, L{SolverOptions.BICGSTAB}, L{SolverOptions.SSOR},                          `SolverOptions.CR`, `SolverOptions.CGS`, `SolverOptions.BICGSTAB`, `SolverOptions.SSOR`,
477                          L{SolverOptions.GMRES}, L{SolverOptions.PRES20}, L{SolverOptions.LUMPING}, L{SolverOptions.ITERATIVE},                          `SolverOptions.GMRES`, `SolverOptions.PRES20`, `SolverOptions.LUMPING`, `SolverOptions.ITERATIVE`,
478                          L{SolverOptions.AMG}, L{SolverOptions.NONLINEAR_GMRES}, L{SolverOptions.TFQMR}, L{SolverOptions.MINRES},                          `SolverOptions.AMG`, `SolverOptions.NONLINEAR_GMRES`, `SolverOptions.TFQMR`, `SolverOptions.MINRES`,
479                          L{SolverOptions.GAUSS_SEIDEL}                          `SolverOptions.GAUSS_SEIDEL`
480          """          """
481          return self.__method          return self.__method
482                    
# Line 442  class SolverOptions(object): Line 484  class SolverOptions(object):
484          """          """
485          Sets the solver package to be used as a solver.            Sets the solver package to be used as a solver.  
486    
487          @param package: key of the solver package to be used.          :param package: key of the solver package to be used.
488          @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`
489          @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.
490          """          """
491      if package==None: package=0      if package==None: package=0
492          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 496  class SolverOptions(object):
496          """          """
497          Returns the solver package key          Returns the solver package key
498    
499          @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`
500          """          """
501          return self.__package          return self.__package
502      def setReordering(self,ordering=30):      def setReordering(self,ordering=30):
# Line 462  class SolverOptions(object): Line 504  class SolverOptions(object):
504          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
505          to optimize compute time and storage use during elimination.          to optimize compute time and storage use during elimination.
506    
507          @param ordering: selects the reordering strategy.          :param ordering: selects the reordering strategy.
508          @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}  
509          """          """
510          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]:
511               raise ValueError,"unknown reordering strategy %s"%ordering               raise ValueError,"unknown reordering strategy %s"%ordering
# Line 473  class SolverOptions(object): Line 514  class SolverOptions(object):
514          """          """
515          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.
516    
517          @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}  
518          """          """
519          return self.__reordering          return self.__reordering
520      def setRestart(self,restart=None):      def setRestart(self,restart=None):
521          """          """
522          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.
523    
524          @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.
525                          restart is performed.          :type restart: ``int`` or ``None``
         @type restart: C{int} or C{None}  
526          """          """
527          if restart == None:          if restart == None:
528              self.__restart=restart              self.__restart=restart
# Line 496  class SolverOptions(object): Line 535  class SolverOptions(object):
535      def getRestart(self):      def getRestart(self):
536          """          """
537          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.
538          If C{None} is returned no restart is performed.          If ``None`` is returned no restart is performed.
539    
540          @rtype: C{int} or C{None}          :rtype: ``int`` or ``None``
541          """          """
542          if self.__restart < 0:          if self.__restart < 0:
543              return None              return None
# Line 515  class SolverOptions(object): Line 554  class SolverOptions(object):
554          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
555          the faster GMRES converged but          the faster GMRES converged but
556    
557          @param truncation: truncation          :param truncation: truncation
558          @type truncation: C{int}          :type truncation: ``int``
559          """          """
560          truncation=int(truncation)          truncation=int(truncation)
561          if truncation<1:          if truncation<1:
# Line 526  class SolverOptions(object): Line 565  class SolverOptions(object):
565          """          """
566          Returns the number of residuals in GMRES to be stored for orthogonalization          Returns the number of residuals in GMRES to be stored for orthogonalization
567    
568          @rtype: C{int}          :rtype: ``int``
569          """          """
570          return self.__truncation          return self.__truncation
571      def setInnerIterMax(self,iter_max=10):      def setInnerIterMax(self,iter_max=10):
572          """          """
573          Sets the maximum number of iteration steps for the inner iteration.          Sets the maximum number of iteration steps for the inner iteration.
574    
575          @param iter_max: maximum number of inner iterations          :param iter_max: maximum number of inner iterations
576          @type iter_max: C{int}          :type iter_max: ``int``
577          """          """
578          iter_max=int(iter_max)          iter_max=int(iter_max)
579          if iter_max<1:          if iter_max<1:
# Line 544  class SolverOptions(object): Line 583  class SolverOptions(object):
583          """          """
584          Returns maximum number of inner iteration steps          Returns maximum number of inner iteration steps
585    
586          @rtype: C{int}          :rtype: ``int``
587          """          """
588          return self.__inner_iter_max          return self.__inner_iter_max
589      def setIterMax(self,iter_max=100000):      def setIterMax(self,iter_max=100000):
590          """          """
591          Sets the maximum number of iteration steps          Sets the maximum number of iteration steps
592    
593          @param iter_max: maximum number of iteration steps          :param iter_max: maximum number of iteration steps
594          @type iter_max: C{int}          :type iter_max: ``int``
595          """          """
596          iter_max=int(iter_max)          iter_max=int(iter_max)
597          if iter_max<1:          if iter_max<1:
# Line 562  class SolverOptions(object): Line 601  class SolverOptions(object):
601          """          """
602          Returns maximum number of iteration steps          Returns maximum number of iteration steps
603    
604          @rtype: C{int}          :rtype: ``int``
605          """          """
606          return self.__iter_max          return self.__iter_max
607      def setLevelMax(self,level_max=10):      def setLevelMax(self,level_max=3):
608          """          """
609          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
610    
611          @param level_max: maximum number of levels          :param level_max: maximum number of levels
612          @type level_max: C{int}          :type level_max: ``int``
613          """          """
614          level_max=int(level_max)          level_max=int(level_max)
615          if level_max<0:          if level_max<0:
# Line 580  class SolverOptions(object): Line 619  class SolverOptions(object):
619          """          """
620          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
621    
622          @rtype: C{int}          :rtype: ``int``
623          """          """
624          return self.__level_max          return self.__level_max
625      def setCoarseningThreshold(self,theta=0.05):      def setCoarseningThreshold(self,theta=0.25):
626          """          """
627          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
628    
629          @param theta: threshold for coarsening          :param theta: threshold for coarsening
630          @type theta: positive C{float}          :type theta: positive ``float``
631          """          """
632          theta=float(theta)          theta=float(theta)
633          if theta<0 or theta>1:          if theta<0 or theta>1:
# Line 598  class SolverOptions(object): Line 637  class SolverOptions(object):
637          """          """
638          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
639    
640          @rtype: C{float}          :rtype: ``float``
641          """          """
642          return self.__coarsening_threshold          return self.__coarsening_threshold
643      def setNumSweeps(self,sweeps=2):      def setNumSweeps(self,sweeps=2):
644          """          """
645          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.
646    
647          @param sweeps: number of sweeps          :param sweeps: number of sweeps
648          @type theta: positive C{int}          :type sweeps: positive ``int``
649          """          """
650          sweeps=int(sweeps)          sweeps=int(sweeps)
651          if sweeps<1:          if sweeps<1:
# Line 616  class SolverOptions(object): Line 655  class SolverOptions(object):
655          """          """
656          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.
657    
658          @rtype: C{int}          :rtype: ``int``
659          """          """
660          return self.__sweeps          return self.__sweeps
661      def setNumPreSweeps(self,sweeps=2):      def setNumPreSweeps(self,sweeps=2):
662          """          """
663          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
664    
665          @param sweeps: number of sweeps          :param sweeps: number of sweeps
666          @type theta: positive C{int}          :type sweeps: positive ``int``
667          """          """
668          sweeps=int(sweeps)          sweeps=int(sweeps)
669          if sweeps<1:          if sweeps<1:
# Line 634  class SolverOptions(object): Line 673  class SolverOptions(object):
673          """          """
674          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
675    
676          @rtype: C{int}          :rtype: ``int``
677          """          """
678          return self.__pre_sweeps          return self.__pre_sweeps
679      def setNumPostSweeps(self,sweeps=2):      def setNumPostSweeps(self,sweeps=2):
680          """          """
681          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
682    
683          @param sweeps: number of sweeps          :param sweeps: number of sweeps
684          @type theta: positive C{int}          :type sweeps: positive ``int``
685          """          """
686          sweeps=int(sweeps)          sweeps=int(sweeps)
687          if sweeps<1:          if sweeps<1:
# Line 652  class SolverOptions(object): Line 691  class SolverOptions(object):
691          """          """
692          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
693    
694          @rtype: C{int}          :rtype: ``int``
695          """          """
696          return self.__post_sweeps          return self.__post_sweeps
697    
# Line 660  class SolverOptions(object): Line 699  class SolverOptions(object):
699          """          """
700          Sets the relative tolerance for the solver          Sets the relative tolerance for the solver
701    
702          @param rtol: relative tolerance          :param rtol: relative tolerance
703          @type rtol: non-negative C{float}          :type rtol: non-negative ``float``
704          """          """
705          rtol=float(rtol)          rtol=float(rtol)
706          if rtol<0 or rtol>1:          if rtol<0 or rtol>1:
# Line 671  class SolverOptions(object): Line 710  class SolverOptions(object):
710          """          """
711          Returns the relative tolerance for the solver          Returns the relative tolerance for the solver
712    
713          @rtype: C{float}          :rtype: ``float``
714          """          """
715          return self.__tolerance          return self.__tolerance
716      def setAbsoluteTolerance(self,atol=0.):      def setAbsoluteTolerance(self,atol=0.):
717          """          """
718          Sets the absolute tolerance for the solver          Sets the absolute tolerance for the solver
719    
720          @param atol:  absolute tolerance          :param atol:  absolute tolerance
721          @type atol: non-negative C{float}          :type atol: non-negative ``float``
722          """          """
723          atol=float(atol)          atol=float(atol)
724          if atol<0:          if atol<0:
# Line 689  class SolverOptions(object): Line 728  class SolverOptions(object):
728          """          """
729          Returns the absolute tolerance for the solver          Returns the absolute tolerance for the solver
730    
731          @rtype: C{float}          :rtype: ``float``
732          """          """
733          return self.__absolute_tolerance          return self.__absolute_tolerance
734    
735      def setInnerTolerance(self,rtol=0.9):      def setInnerTolerance(self,rtol=0.9):
736          """          """
737           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.  
738    
739          @param rtol: inner relative tolerance          :param rtol: inner relative tolerance
740          @type rtol: positive C{float}          :type rtol: positive ``float``
741          """          """
742          rtol=float(rtol)          rtol=float(rtol)
743          if rtol<=0 or rtol>1:          if rtol<=0 or rtol>1:
# Line 709  class SolverOptions(object): Line 747  class SolverOptions(object):
747          """          """
748          Returns the relative tolerance for an inner iteration scheme          Returns the relative tolerance for an inner iteration scheme
749    
750          @rtype: C{float}          :rtype: ``float``
751          """          """
752          return self.__inner_tolerance          return self.__inner_tolerance
753      def setDropTolerance(self,drop_tol=0.01):      def setDropTolerance(self,drop_tol=0.01):
754          """          """
755          Sets the relative drop tolerance in ILUT          Sets the relative drop tolerance in ILUT
756    
757          @param drop_tol: drop tolerance          :param drop_tol: drop tolerance
758          @type drop_tol: positive C{float}          :type drop_tol: positive ``float``
759          """          """
760          drop_tol=float(drop_tol)          drop_tol=float(drop_tol)
761          if drop_tol<=0 or drop_tol>1:          if drop_tol<=0 or drop_tol>1:
# Line 727  class SolverOptions(object): Line 765  class SolverOptions(object):
765          """          """
766          Returns the relative drop tolerance in ILUT          Returns the relative drop tolerance in ILUT
767    
768          @rtype: C{float}          :rtype: ``float``
769          """          """
770          return self.__drop_tolerance          return self.__drop_tolerance
771      def setDropStorage(self,storage=2.):      def setDropStorage(self,storage=2.):
772          """          """
773          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
774          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.
775    
776          @param storage: allowed storage increase          :param storage: allowed storage increase
777          @type storage: C{float}          :type storage: ``float``
778          """          """
779          storage=float(storage)          storage=float(storage)
780          if storage<1:          if storage<1:
# Line 747  class SolverOptions(object): Line 785  class SolverOptions(object):
785          """          """
786          Returns the maximum allowed increase in storage for ILUT          Returns the maximum allowed increase in storage for ILUT
787    
788          @rtype: C{float}          :rtype: ``float``
789          """          """
790          return self.__drop_storage          return self.__drop_storage
791      def setRelaxationFactor(self,factor=0.3):      def setRelaxationFactor(self,factor=0.3):
792          """          """
793          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.
794    
795          @param factor: relaxation factor          :param factor: relaxation factor
796          @type factor: C{float}          :type factor: ``float``
797          @note: RILU with a relaxation factor 0 is identical to ILU0          :note: RILU with a relaxation factor 0 is identical to ILU0
798          """          """
799          factor=float(factor)          factor=float(factor)
800          if factor<0:          if factor<0:
# Line 767  class SolverOptions(object): Line 805  class SolverOptions(object):
805          """          """
806          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.
807    
808          @rtype: C{float}          :rtype: ``float``
809          """          """
810          return self.__relaxation          return self.__relaxation
811      def isSymmetric(self):      def isSymmetric(self):
812          """          """
813          Checks if symmetry of the  coefficient matrix is indicated.          Checks if symmetry of the  coefficient matrix is indicated.
814    
815          @return: True if a symmetric PDE is indicated, False otherwise          :return: True if a symmetric PDE is indicated, False otherwise
816          @rtype: C{bool}          :rtype: ``bool``
817          """          """
818          return self.__symmetric          return self.__symmetric
819      def setSymmetryOn(self):      def setSymmetryOn(self):
# Line 790  class SolverOptions(object): Line 828  class SolverOptions(object):
828          self.__symmetric=False          self.__symmetric=False
829      def setSymmetry(self,flag=False):      def setSymmetry(self,flag=False):
830          """          """
831          Sets the symmetry flag for the coefficient matrix to C{flag}.          Sets the symmetry flag for the coefficient matrix to ``flag``.
832    
833          @param flag: If True, the symmetry flag is set otherwise reset.          :param flag: If True, the symmetry flag is set otherwise reset.
834          @type flag: C{bool}          :type flag: ``bool``
835          """          """
836          if flag:          if flag:
837              self.setSymmetryOn()              self.setSymmetryOn()
# Line 801  class SolverOptions(object): Line 839  class SolverOptions(object):
839              self.setSymmetryOff()              self.setSymmetryOff()
840      def isVerbose(self):      def isVerbose(self):
841          """          """
842          Returns C{True} if the solver is expected to be verbose.          Returns ``True`` if the solver is expected to be verbose.
843    
844          @return: True if verbosity of switched on.          :return: True if verbosity of switched on.
845          @rtype: C{bool}          :rtype: ``bool``
846          """          """
847          return self.__verbose          return self.__verbose
848    
# Line 820  class SolverOptions(object): Line 858  class SolverOptions(object):
858          self.__verbose=False          self.__verbose=False
859      def setVerbosity(self,verbose=False):      def setVerbosity(self,verbose=False):
860          """          """
861          Sets the verbosity flag for the solver to C{flag}.          Sets the verbosity flag for the solver to ``flag``.
862    
863          @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.
864          @type flag: C{bool}          :type verbose: ``bool``
865          """          """
866          if verbose:          if verbose:
867              self.setVerbosityOn()              self.setVerbosityOn()
# Line 832  class SolverOptions(object): Line 870  class SolverOptions(object):
870                    
871      def adaptInnerTolerance(self):      def adaptInnerTolerance(self):
872          """          """
873          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.
874          Otherwise the inner tolerance set by L{setInnerTolerance} is used.          Otherwise the inner tolerance set by `setInnerTolerance` is used.
875    
876          @return: C{True} if inner tolerance adaption is chosen.          :return: ``True`` if inner tolerance adaption is chosen.
877          @rtype: C{bool}          :rtype: ``bool``
878          """          """
879          return self.__adapt_inner_tolerance          return self.__adapt_inner_tolerance
880    
# Line 854  class SolverOptions(object): Line 892  class SolverOptions(object):
892          """          """
893          Sets a flag to indicate automatic selection of the inner tolerance.          Sets a flag to indicate automatic selection of the inner tolerance.
894    
895          @param adapt: If C{True}, the inner tolerance is selected automatically.          :param adapt: If ``True``, the inner tolerance is selected automatically.
896          @type adapt: C{bool}          :type adapt: ``bool``
897          """          """
898          if adapt:          if adapt:
899              self.setInnerToleranceAdaptionOn()              self.setInnerToleranceAdaptionOn()
# Line 864  class SolverOptions(object): Line 902  class SolverOptions(object):
902    
903      def acceptConvergenceFailure(self):      def acceptConvergenceFailure(self):
904          """          """
905          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
906          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
907          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
908          continue even if the returned the solution does not necessarily meet the          continue even if the returned the solution does not necessarily meet the
909          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
910          last call to the solver was successful.          last call to the solver was successful.
911    
912          @return: C{True} if a failure to achieve convergence is accepted.          :return: ``True`` if a failure to achieve convergence is accepted.
913          @rtype: C{bool}          :rtype: ``bool``
914          """          """
915          return self.__accept_convergence_failure          return self.__accept_convergence_failure
916    
# Line 890  class SolverOptions(object): Line 928  class SolverOptions(object):
928          """          """
929          Sets a flag to indicate the acceptance of a failure of convergence.          Sets a flag to indicate the acceptance of a failure of convergence.
930    
931          @param accept: If C{True}, any failure to achieve convergence is accepted.          :param accept: If ``True``, any failure to achieve convergence is accepted.
932          @type accept: C{bool}          :type accept: ``bool``
933          """          """
934          if accept:          if accept:
935              self.setAcceptanceConvergenceFailureOn()              self.setAcceptanceConvergenceFailureOn()
# Line 927  class PDECoef(object): Line 965  class PDECoef(object):
965      """      """
966      A class for describing a PDE coefficient.      A class for describing a PDE coefficient.
967    
968      @cvar INTERIOR: indicator that coefficient is defined on the interior of      :cvar INTERIOR: indicator that coefficient is defined on the interior of
969                      the PDE domain                      the PDE domain
970      @cvar BOUNDARY: indicator that coefficient is defined on the boundary of      :cvar BOUNDARY: indicator that coefficient is defined on the boundary of
971                      the PDE domain                      the PDE domain
972      @cvar CONTACT: indicator that coefficient is defined on the contact region      :cvar CONTACT: indicator that coefficient is defined on the contact region
973                     within the PDE domain                     within the PDE domain
974      @cvar INTERIOR_REDUCED: indicator that coefficient is defined on the      :cvar INTERIOR_REDUCED: indicator that coefficient is defined on the
975                              interior of the PDE domain using a reduced                              interior of the PDE domain using a reduced
976                              integration order                              integration order
977      @cvar BOUNDARY_REDUCED: indicator that coefficient is defined on the      :cvar BOUNDARY_REDUCED: indicator that coefficient is defined on the
978                              boundary of the PDE domain using a reduced                              boundary of the PDE domain using a reduced
979                              integration order                              integration order
980      @cvar CONTACT_REDUCED: indicator that coefficient is defined on the contact      :cvar CONTACT_REDUCED: indicator that coefficient is defined on the contact
981                             region within the PDE domain using a reduced                             region within the PDE domain using a reduced
982                             integration order                             integration order
983      @cvar SOLUTION: indicator that coefficient is defined through a solution of      :cvar SOLUTION: indicator that coefficient is defined through a solution of
984                      the PDE                      the PDE
985      @cvar REDUCED: indicator that coefficient is defined through a reduced      :cvar REDUCED: indicator that coefficient is defined through a reduced
986                     solution of the PDE                     solution of the PDE
987      @cvar BY_EQUATION: indicator that the dimension of the coefficient shape is      :cvar BY_EQUATION: indicator that the dimension of the coefficient shape is
988                         defined by the number of PDE equations                         defined by the number of PDE equations
989      @cvar BY_SOLUTION: indicator that the dimension of the coefficient shape is      :cvar BY_SOLUTION: indicator that the dimension of the coefficient shape is
990                         defined by the number of PDE solutions                         defined by the number of PDE solutions
991      @cvar BY_DIM: indicator that the dimension of the coefficient shape is      :cvar BY_DIM: indicator that the dimension of the coefficient shape is
992                    defined by the spatial dimension                    defined by the spatial dimension
993      @cvar OPERATOR: indicator that the the coefficient alters the operator of      :cvar OPERATOR: indicator that the the coefficient alters the operator of
994                      the PDE                      the PDE
995      @cvar RIGHTHANDSIDE: indicator that the the coefficient alters the right      :cvar RIGHTHANDSIDE: indicator that the the coefficient alters the right
996                           hand side of the PDE                           hand side of the PDE
997      @cvar BOTH: indicator that the the coefficient alters the operator as well      :cvar BOTH: indicator that the the coefficient alters the operator as well
998                  as the right hand side of the PDE                  as the right hand side of the PDE
999    
1000      """      """
# Line 979  class PDECoef(object): Line 1017  class PDECoef(object):
1017         """         """
1018         Initialises a PDE coefficient type.         Initialises a PDE coefficient type.
1019    
1020         @param where: describes where the coefficient lives         :param where: describes where the coefficient lives
1021         @type where: one of L{INTERIOR}, L{BOUNDARY}, L{CONTACT}, L{SOLUTION},         :type where: one of `INTERIOR`, `BOUNDARY`, `CONTACT`, `SOLUTION`,
1022                      L{REDUCED}, L{INTERIOR_REDUCED}, L{BOUNDARY_REDUCED},                      `REDUCED`, `INTERIOR_REDUCED`, `BOUNDARY_REDUCED`,
1023                      L{CONTACT_REDUCED}                      `CONTACT_REDUCED`
1024         @param pattern: describes the shape of the coefficient and how the shape         :param pattern: describes the shape of the coefficient and how the shape
1025                         is built for a given spatial dimension and numbers of                         is built for a given spatial dimension and numbers of
1026                         equations and solutions in then PDE. For instance,                         equations and solutions in then PDE. For instance,
1027                         (L{BY_EQUATION},L{BY_SOLUTION},L{BY_DIM}) describes a                         (`BY_EQUATION`,`BY_SOLUTION`,`BY_DIM`) describes a
1028                         rank 3 coefficient which is instantiated as shape (3,2,2)                         rank 3 coefficient which is instantiated as shape (3,2,2)
1029                         in case of three equations and two solution components                         in case of three equations and two solution components
1030                         on a 2-dimensional domain. In the case of single equation                         on a 2-dimensional domain. In the case of single equation
1031                         and a single solution component the shape components                         and a single solution component the shape components
1032                         marked by L{BY_EQUATION} or L{BY_SOLUTION} are dropped.                         marked by `BY_EQUATION` or `BY_SOLUTION` are dropped.
1033                         In this case the example would be read as (2,).                         In this case the example would be read as (2,).
1034         @type pattern: C{tuple} of L{BY_EQUATION}, L{BY_SOLUTION}, L{BY_DIM}         :type pattern: ``tuple`` of `BY_EQUATION`, `BY_SOLUTION`, `BY_DIM`
1035         @param altering: indicates what part of the PDE is altered if the         :param altering: indicates what part of the PDE is altered if the
1036                          coefficient is altered                          coefficient is altered
1037         @type altering: one of L{OPERATOR}, L{RIGHTHANDSIDE}, L{BOTH}         :type altering: one of `OPERATOR`, `RIGHTHANDSIDE`, `BOTH`
1038         """         """
1039         super(PDECoef, self).__init__()         super(PDECoef, self).__init__()
1040         self.what=where         self.what=where
# Line 1012  class PDECoef(object): Line 1050  class PDECoef(object):
1050    
1051      def getFunctionSpace(self,domain,reducedEquationOrder=False,reducedSolutionOrder=False):      def getFunctionSpace(self,domain,reducedEquationOrder=False,reducedSolutionOrder=False):
1052         """         """
1053         Returns the L{FunctionSpace<escript.FunctionSpace>} of the coefficient.         Returns the `FunctionSpace` of the coefficient.
1054    
1055         @param domain: domain on which the PDE uses the coefficient         :param domain: domain on which the PDE uses the coefficient
1056         @type domain: L{Domain<escript.Domain>}         :type domain: `Domain`
1057         @param reducedEquationOrder: True to indicate that reduced order is used         :param reducedEquationOrder: True to indicate that reduced order is used
1058                                      to represent the equation                                      to represent the equation
1059         @type reducedEquationOrder: C{bool}         :type reducedEquationOrder: ``bool``
1060         @param reducedSolutionOrder: True to indicate that reduced order is used         :param reducedSolutionOrder: True to indicate that reduced order is used
1061                                      to represent the solution                                      to represent the solution
1062         @type reducedSolutionOrder: C{bool}         :type reducedSolutionOrder: ``bool``
1063         @return: L{FunctionSpace<escript.FunctionSpace>} of the coefficient         :return: `FunctionSpace` of the coefficient
1064         @rtype: L{FunctionSpace<escript.FunctionSpace>}         :rtype: `FunctionSpace`
1065         """         """
1066         if self.what==self.INTERIOR:         if self.what==self.INTERIOR:
1067              return escript.Function(domain)              return escript.Function(domain)
# Line 1049  class PDECoef(object): Line 1087  class PDECoef(object):
1087         """         """
1088         Returns the value of the coefficient.         Returns the value of the coefficient.
1089    
1090         @return: value of the coefficient         :return: value of the coefficient
1091         @rtype: L{Data<escript.Data>}         :rtype: `Data`
1092         """         """
1093         return self.value         return self.value
1094    
# Line 1058  class PDECoef(object): Line 1096  class PDECoef(object):
1096         """         """
1097         Sets the value of the coefficient to a new value.         Sets the value of the coefficient to a new value.
1098    
1099         @param domain: domain on which the PDE uses the coefficient         :param domain: domain on which the PDE uses the coefficient
1100         @type domain: L{Domain<escript.Domain>}         :type domain: `Domain`
1101         @param numEquations: number of equations of the PDE         :param numEquations: number of equations of the PDE
1102         @type numEquations: C{int}         :type numEquations: ``int``
1103         @param numSolutions: number of components of the PDE solution         :param numSolutions: number of components of the PDE solution
1104         @type numSolutions: C{int}         :type numSolutions: ``int``
1105         @param reducedEquationOrder: True to indicate that reduced order is used         :param reducedEquationOrder: True to indicate that reduced order is used
1106                                      to represent the equation                                      to represent the equation
1107         @type reducedEquationOrder: C{bool}         :type reducedEquationOrder: ``bool``
1108         @param reducedSolutionOrder: True to indicate that reduced order is used         :param reducedSolutionOrder: True to indicate that reduced order is used
1109                                      to represent the solution                                      to represent the solution
1110         @type reducedSolutionOrder: C{bool}         :type reducedSolutionOrder: ``bool``
1111         @param newValue: number of components of the PDE solution         :param newValue: number of components of the PDE solution
1112         @type newValue: any object that can be converted into a         :type newValue: any object that can be converted into a
1113                         L{Data<escript.Data>} object with the appropriate shape                         `Data` object with the appropriate shape
1114                         and L{FunctionSpace<escript.FunctionSpace>}                         and `FunctionSpace`
1115         @raise IllegalCoefficientValue: if the shape of the assigned value does         :raise IllegalCoefficientValue: if the shape of the assigned value does
1116                                         not match the shape of the coefficient                                         not match the shape of the coefficient
1117         @raise IllegalCoefficientFunctionSpace: if unable to interpolate value         :raise IllegalCoefficientFunctionSpace: if unable to interpolate value
1118                                                 to appropriate function space                                                 to appropriate function space
1119         """         """
1120         if newValue==None:         if newValue==None:
# Line 1099  class PDECoef(object): Line 1137  class PDECoef(object):
1137          """          """
1138          Checks if the coefficient alters the operator of the PDE.          Checks if the coefficient alters the operator of the PDE.
1139    
1140          @return: True if the operator of the PDE is changed when the          :return: True if the operator of the PDE is changed when the
1141                   coefficient is changed                   coefficient is changed
1142          @rtype: C{bool}          :rtype: ``bool``
1143          """          """
1144          if self.altering==self.OPERATOR or self.altering==self.BOTH:          if self.altering==self.OPERATOR or self.altering==self.BOTH:
1145              return not None              return not None
# Line 1112  class PDECoef(object): Line 1150  class PDECoef(object):
1150          """          """
1151          Checks if the coefficient alters the right hand side of the PDE.          Checks if the coefficient alters the right hand side of the PDE.
1152    
1153          @rtype: C{bool}          :rtype: ``bool``
1154          @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
1155                   coefficient is changed, C{None} otherwise.                   coefficient is changed, ``None`` otherwise.
1156          """          """
1157          if self.altering==self.RIGHTHANDSIDE or self.altering==self.BOTH:          if self.altering==self.RIGHTHANDSIDE or self.altering==self.BOTH:
1158              return not None              return not None
# Line 1126  class PDECoef(object): Line 1164  class PDECoef(object):
1164         Tries to estimate the number of equations and number of solutions if         Tries to estimate the number of equations and number of solutions if
1165         the coefficient has the given shape.         the coefficient has the given shape.
1166    
1167         @param domain: domain on which the PDE uses the coefficient         :param domain: domain on which the PDE uses the coefficient
1168         @type domain: L{Domain<escript.Domain>}         :type domain: `Domain`
1169         @param shape: suggested shape of the coefficient         :param shape: suggested shape of the coefficient
1170         @type shape: C{tuple} of C{int} values         :type shape: ``tuple`` of ``int`` values
1171         @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
1172                  the coefficient has given shape. If no appropriate numbers                  the coefficient has given shape. If no appropriate numbers
1173                  could be identified, C{None} is returned                  could be identified, ``None`` is returned
1174         @rtype: C{tuple} of two C{int} values or C{None}         :rtype: ``tuple`` of two ``int`` values or ``None``
1175         """         """
1176         dim=domain.getDim()         dim=domain.getDim()
1177         if len(shape)>0:         if len(shape)>0:
# Line 1174  class PDECoef(object): Line 1212  class PDECoef(object):
1212         Checks if the coefficient allows to estimate the number of solution         Checks if the coefficient allows to estimate the number of solution
1213         components.         components.
1214    
1215         @return: True if the coefficient allows an estimate of the number of         :return: True if the coefficient allows an estimate of the number of
1216                  solution components, False otherwise                  solution components, False otherwise
1217         @rtype: C{bool}         :rtype: ``bool``
1218         """         """
1219         for i in self.pattern:         for i in self.pattern:
1220               if i==self.BY_SOLUTION: return True               if i==self.BY_SOLUTION: return True
# Line 1186  class PDECoef(object): Line 1224  class PDECoef(object):
1224         """         """
1225         Checks if the coefficient allows to estimate the number of equations.         Checks if the coefficient allows to estimate the number of equations.
1226    
1227         @return: True if the coefficient allows an estimate of the number of         :return: True if the coefficient allows an estimate of the number of
1228                  equations, False otherwise                  equations, False otherwise
1229         @rtype: C{bool}         :rtype: ``bool``
1230         """         """
1231         for i in self.pattern:         for i in self.pattern:
1232               if i==self.BY_EQUATION: return True               if i==self.BY_EQUATION: return True
# Line 1199  class PDECoef(object): Line 1237  class PDECoef(object):
1237        Compares two tuples of possible number of equations and number of        Compares two tuples of possible number of equations and number of
1238        solutions.        solutions.
1239    
1240        @param t1: the first tuple        :param t1: the first tuple
1241        @param t2: the second tuple        :param t2: the second tuple
1242        @return: 0, 1, or -1        :return: 0, 1, or -1
1243        """        """
1244    
1245        dif=t1[0]+t1[1]-(t2[0]+t2[1])        dif=t1[0]+t1[1]-(t2[0]+t2[1])
# Line 1213  class PDECoef(object): Line 1251  class PDECoef(object):
1251         """         """
1252         Builds the required shape of the coefficient.         Builds the required shape of the coefficient.
1253    
1254         @param domain: domain on which the PDE uses the coefficient         :param domain: domain on which the PDE uses the coefficient
1255         @type domain: L{Domain<escript.Domain>}         :type domain: `Domain`
1256         @param numEquations: number of equations of the PDE         :param numEquations: number of equations of the PDE
1257         @type numEquations: C{int}         :type numEquations: ``int``
1258         @param numSolutions: number of components of the PDE solution         :param numSolutions: number of components of the PDE solution
1259         @type numSolutions: C{int}         :type numSolutions: ``int``
1260         @return: shape of the coefficient         :return: shape of the coefficient
1261         @rtype: C{tuple} of C{int} values         :rtype: ``tuple`` of ``int`` values
1262         """         """
1263         dim=domain.getDim()         dim=domain.getDim()
1264         s=()         s=()
# Line 1238  class PDECoef(object): Line 1276  class PDECoef(object):
1276  class LinearProblem(object):  class LinearProblem(object):
1277     """     """
1278     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
1279     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
1280     L{Domain<escript.Domain>} object. The problem can be given as a single     `Domain` object. The problem can be given as a single
1281     equation or as a system of equations.     equation or as a system of equations.
1282    
1283     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
1284     a problem of the form     a problem of the form
1285    
1286     M{L u=f}     *L u=f*
1287    
1288     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
1289     problem will be solved to get the unknown M{u}.     problem will be solved to get the unknown *u*.
1290    
1291     """     """
1292     def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):     def __init__(self,domain,numEquations=None,numSolutions=None,debug=False):
1293       """       """
1294       Initializes a linear problem.       Initializes a linear problem.
1295    
1296       @param domain: domain of the PDE       :param domain: domain of the PDE
1297       @type domain: L{Domain<escript.Domain>}       :type domain: `Domain`
1298       @param numEquations: number of equations. If C{None} the number of       :param numEquations: number of equations. If ``None`` the number of
1299                            equations is extracted from the coefficients.                            equations is extracted from the coefficients.
1300       @param numSolutions: number of solution components. If C{None} the number       :param numSolutions: number of solution components. If ``None`` the number
1301                            of solution components is extracted from the                            of solution components is extracted from the
1302                            coefficients.                            coefficients.
1303       @param debug: if True debug information is printed       :param debug: if True debug information is printed
1304    
1305       """       """
1306       super(LinearProblem, self).__init__()       super(LinearProblem, self).__init__()
# Line 1292  class LinearProblem(object): Line 1330  class LinearProblem(object):
1330       """       """
1331       Returns a string representation of the PDE.       Returns a string representation of the PDE.
1332    
1333       @return: a simple representation of the PDE       :return: a simple representation of the PDE
1334       @rtype: C{str}       :rtype: ``str``
1335       """       """
1336       return "<LinearProblem %d>"%id(self)       return "<LinearProblem %d>"%id(self)
1337     # ==========================================================================     # ==========================================================================
# Line 1313  class LinearProblem(object): Line 1351  class LinearProblem(object):
1351    
1352     def setDebug(self, flag):     def setDebug(self, flag):
1353       """       """
1354       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.
1355    
1356       @param flag: desired debug status       :param flag: desired debug status
1357       @type flag: C{bool}       :type flag: ``bool``
1358       """       """
1359       if flag:       if flag:
1360           self.setDebugOn()           self.setDebugOn()
# Line 1327  class LinearProblem(object): Line 1365  class LinearProblem(object):
1365       """       """
1366       Prints the text message if debug mode is switched on.       Prints the text message if debug mode is switched on.
1367    
1368       @param text: message to be printed       :param text: message to be printed
1369       @type text: C{string}       :type text: ``string``
1370       """       """
1371       if self.__debug: print "%s: %s"%(str(self),text)       if self.__debug: print "%s: %s"%(str(self),text)
1372    
# Line 1339  class LinearProblem(object): Line 1377  class LinearProblem(object):
1377         """         """
1378         Introduces new coefficients into the problem.         Introduces new coefficients into the problem.
1379    
1380         Use::         Use:
1381    
1382         p.introduceCoefficients(A=PDECoef(...), B=PDECoef(...))         p.introduceCoefficients(A=PDECoef(...), B=PDECoef(...))
1383    
1384         to introduce the coefficients M{A} ans M{B}.         to introduce the coefficients *A* and *B*.
1385         """         """
1386         for name, type in coeff.items():         for name, type in coeff.items():
1387             if not isinstance(type,PDECoef):             if not isinstance(type,PDECoef):
# Line 1356  class LinearProblem(object): Line 1394  class LinearProblem(object):
1394       """       """
1395       Returns the domain of the PDE.       Returns the domain of the PDE.
1396    
1397       @return: the domain of the PDE       :return: the domain of the PDE
1398       @rtype: L{Domain<escript.Domain>}       :rtype: `Domain`
1399       """       """
1400       return self.__domain       return self.__domain
1401     def getDomainStatus(self):     def getDomainStatus(self):
# Line 1373  class LinearProblem(object): Line 1411  class LinearProblem(object):
1411       return self.__system_status       return self.__system_status
1412     def setSystemStatus(self,status=None):     def setSystemStatus(self,status=None):
1413       """       """
1414       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
1415       current status of the domain is used.       current status of the domain is used.
1416       """       """
1417       if status == None:       if status == None:
# Line 1385  class LinearProblem(object): Line 1423  class LinearProblem(object):
1423       """       """
1424       Returns the spatial dimension of the PDE.       Returns the spatial dimension of the PDE.
1425    
1426       @return: the spatial dimension of the PDE domain       :return: the spatial dimension of the PDE domain
1427       @rtype: C{int}       :rtype: ``int``
1428       """       """
1429       return self.getDomain().getDim()       return self.getDomain().getDim()
1430    
# Line 1394  class LinearProblem(object): Line 1432  class LinearProblem(object):
1432       """       """
1433       Returns the number of equations.       Returns the number of equations.
1434    
1435       @return: the number of equations       :return: the number of equations
1436       @rtype: C{int}       :rtype: ``int``
1437       @raise UndefinedPDEError: if the number of equations is not specified yet       :raise UndefinedPDEError: if the number of equations is not specified yet
1438       """       """
1439       if self.__numEquations==None:       if self.__numEquations==None:
1440           if self.__numSolutions==None:           if self.__numSolutions==None:
# Line 1409  class LinearProblem(object): Line 1447  class LinearProblem(object):
1447       """       """
1448       Returns the number of unknowns.       Returns the number of unknowns.
1449    
1450       @return: the number of unknowns       :return: the number of unknowns
1451       @rtype: C{int}       :rtype: ``int``
1452       @raise UndefinedPDEError: if the number of unknowns is not specified yet       :raise UndefinedPDEError: if the number of unknowns is not specified yet
1453       """       """
1454       if self.__numSolutions==None:       if self.__numSolutions==None:
1455          if self.__numEquations==None:          if self.__numEquations==None:
# Line 1424  class LinearProblem(object): Line 1462  class LinearProblem(object):
1462       """       """
1463       Returns the status of order reduction for the equation.       Returns the status of order reduction for the equation.
1464    
1465       @return: True if reduced interpolation order is used for the       :return: True if reduced interpolation order is used for the
1466                representation of the equation, False otherwise                representation of the equation, False otherwise
1467       @rtype: L{bool}       :rtype: `bool`
1468       """       """
1469       return self.__reduce_equation_order       return self.__reduce_equation_order
1470    
# Line 1434  class LinearProblem(object): Line 1472  class LinearProblem(object):
1472       """       """
1473       Returns the status of order reduction for the solution.       Returns the status of order reduction for the solution.
1474    
1475       @return: True if reduced interpolation order is used for the       :return: True if reduced interpolation order is used for the
1476                representation of the solution, False otherwise                representation of the solution, False otherwise
1477       @rtype: L{bool}       :rtype: `bool`
1478       """       """
1479       return self.__reduce_solution_order       return self.__reduce_solution_order
1480    
1481     def getFunctionSpaceForEquation(self):     def getFunctionSpaceForEquation(self):
1482       """       """
1483       Returns the L{FunctionSpace<escript.FunctionSpace>} used to discretize       Returns the `FunctionSpace` used to discretize
1484       the equation.       the equation.
1485    
1486       @return: representation space of equation       :return: representation space of equation
1487       @rtype: L{FunctionSpace<escript.FunctionSpace>}       :rtype: `FunctionSpace`
1488       """       """
1489       if self.reduceEquationOrder():       if self.reduceEquationOrder():
1490           return escript.ReducedSolution(self.getDomain())           return escript.ReducedSolution(self.getDomain())
# Line 1455  class LinearProblem(object): Line 1493  class LinearProblem(object):
1493    
1494     def getFunctionSpaceForSolution(self):     def getFunctionSpaceForSolution(self):
1495       """       """
1496       Returns the L{FunctionSpace<escript.FunctionSpace>} used to represent       Returns the `FunctionSpace` used to represent
1497       the solution.       the solution.
1498    
1499       @return: representation space of solution       :return: representation space of solution
1500       @rtype: L{FunctionSpace<escript.FunctionSpace>}       :rtype: `FunctionSpace`
1501       """       """
1502       if self.reduceSolutionOrder():       if self.reduceSolutionOrder():
1503           return escript.ReducedSolution(self.getDomain())           return escript.ReducedSolution(self.getDomain())
# Line 1473  class LinearProblem(object): Line 1511  class LinearProblem(object):
1511         """         """
1512         Sets the solver options.         Sets the solver options.
1513    
1514         @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.
1515         @type options: L{SolverOptions} or C{None}         :type options: `SolverOptions` or ``None``
1516         @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`.
1517         """         """
1518         if options==None:         if options==None:
1519            self.__solver_options=SolverOptions()            self.__solver_options=SolverOptions()
# Line 1489  class LinearProblem(object): Line 1527  class LinearProblem(object):
1527         """         """
1528         Returns the solver options         Returns the solver options
1529        
1530         @rtype: L{SolverOptions}         :rtype: `SolverOptions`
1531         """         """
1532         self.__solver_options.setSymmetry(self.__sym)         self.__solver_options.setSymmetry(self.__sym)
1533         return self.__solver_options         return self.__solver_options
# Line 1498  class LinearProblem(object): Line 1536  class LinearProblem(object):
1536        """        """
1537        Checks if matrix lumping is the current solver method.        Checks if matrix lumping is the current solver method.
1538    
1539        @return: True if the current solver method is lumping        :return: True if the current solver method is lumping
1540        @rtype: C{bool}        :rtype: ``bool``
1541        """        """
1542        return self.getSolverOptions().getSolverMethod()==self.getSolverOptions().LUMPING        return self.getSolverOptions().getSolverMethod()==self.getSolverOptions().LUMPING
1543     # ==========================================================================     # ==========================================================================
# Line 1509  class LinearProblem(object): Line 1547  class LinearProblem(object):
1547        """        """
1548        Checks if symmetry is indicated.        Checks if symmetry is indicated.
1549    
1550        @return: True if a symmetric PDE is indicated, False otherwise        :return: True if a symmetric PDE is indicated, False otherwise
1551        @rtype: C{bool}        :rtype: ``bool``
1552        @note: the method is equivalent to use getSolverOptions().isSymmetric()        :note: the method is equivalent to use getSolverOptions().isSymmetric()
1553        """        """
1554        self.getSolverOptions().isSymmetric()        self.getSolverOptions().isSymmetric()
1555    
1556     def setSymmetryOn(self):     def setSymmetryOn(self):
1557        """        """
1558        Sets the symmetry flag.        Sets the symmetry flag.
1559        @note: The method overwrites the symmetry flag set by the solver options        :note: The method overwrites the symmetry flag set by the solver options
1560        """        """
1561        self.__sym=True        self.__sym=True
1562        self.getSolverOptions().setSymmetryOn()        self.getSolverOptions().setSymmetryOn()
# Line 1526  class LinearProblem(object): Line 1564  class LinearProblem(object):
1564     def setSymmetryOff(self):     def setSymmetryOff(self):
1565        """        """
1566        Clears the symmetry flag.        Clears the symmetry flag.
1567        @note: The method overwrites the symmetry flag set by the solver options        :note: The method overwrites the symmetry flag set by the solver options
1568        """        """
1569        self.__sym=False        self.__sym=False
1570        self.getSolverOptions().setSymmetryOff()        self.getSolverOptions().setSymmetryOff()
1571    
1572     def setSymmetry(self,flag=False):     def setSymmetry(self,flag=False):
1573        """        """
1574        Sets the symmetry flag to C{flag}.        Sets the symmetry flag to ``flag``.
1575    
1576        @param flag: If True, the symmetry flag is set otherwise reset.        :param flag: If True, the symmetry flag is set otherwise reset.
1577        @type flag: C{bool}        :type flag: ``bool``
1578        @note: The method overwrites the symmetry flag set by the solver options        :note: The method overwrites the symmetry flag set by the solver options
1579        """        """
1580        self.getSolverOptions().setSymmetry(flag)        self.getSolverOptions().setSymmetry(flag)
1581     # ==========================================================================     # ==========================================================================
# Line 1547  class LinearProblem(object): Line 1585  class LinearProblem(object):
1585       """       """
1586       Switches reduced order on for solution and equation representation.       Switches reduced order on for solution and equation representation.
1587    
1588       @raise RuntimeError: if order reduction is altered after a coefficient has       :raise RuntimeError: if order reduction is altered after a coefficient has
1589                            been set                            been set
1590       """       """
1591       self.setReducedOrderForSolutionOn()       self.setReducedOrderForSolutionOn()
# Line 1557  class LinearProblem(object): Line 1595  class LinearProblem(object):
1595       """       """
1596       Switches reduced order off for solution and equation representation       Switches reduced order off for solution and equation representation
1597    
1598       @raise RuntimeError: if order reduction is altered after a coefficient has       :raise RuntimeError: if order reduction is altered after a coefficient has
1599                            been set                            been set
1600       """       """
1601       self.setReducedOrderForSolutionOff()       self.setReducedOrderForSolutionOff()
# Line 1568  class LinearProblem(object): Line 1606  class LinearProblem(object):
1606       Sets order reduction state for both solution and equation representation       Sets order reduction state for both solution and equation representation
1607       according to flag.       according to flag.
1608    
1609       @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
1610                    and equation representation, otherwise or if flag is not                    and equation representation, otherwise or if flag is not
1611                    present order reduction is switched off                    present order reduction is switched off
1612       @type flag: C{bool}       :type flag: ``bool``
1613       @raise RuntimeError: if order reduction is altered after a coefficient has       :raise RuntimeError: if order reduction is altered after a coefficient has
1614                            been set                            been set
1615       """       """
1616       self.setReducedOrderForSolutionTo(flag)       self.setReducedOrderForSolutionTo(flag)
# Line 1583  class LinearProblem(object): Line 1621  class LinearProblem(object):
1621       """       """
1622       Switches reduced order on for solution representation.       Switches reduced order on for solution representation.
1623    
1624       @raise RuntimeError: if order reduction is altered after a coefficient has       :raise RuntimeError: if order reduction is altered after a coefficient has
1625                            been set                            been set
1626       """       """
1627       if not self.__reduce_solution_order:       if not self.__reduce_solution_order:
# Line 1597  class LinearProblem(object): Line 1635  class LinearProblem(object):
1635       """       """
1636       Switches reduced order off for solution representation       Switches reduced order off for solution representation
1637    
1638       @raise RuntimeError: if order reduction is altered after a coefficient has       :raise RuntimeError: if order reduction is altered after a coefficient has
1639                            been set.                            been set.
1640       """       """
1641       if self.__reduce_solution_order:       if self.__reduce_solution_order:
# Line 1611  class LinearProblem(object): Line 1649  class LinearProblem(object):
1649       """       """
1650       Sets order reduction state for solution representation according to flag.       Sets order reduction state for solution representation according to flag.
1651    
1652       @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
1653                    solution representation, otherwise or if flag is not present                    solution representation, otherwise or if flag is not present
1654                    order reduction is switched off                    order reduction is switched off
1655       @type flag: C{bool}       :type flag: ``bool``
1656       @raise RuntimeError: if order reduction is altered after a coefficient has       :raise RuntimeError: if order reduction is altered after a coefficient has
1657                            been set                            been set
1658       """       """
1659       if flag:       if flag:
# Line 1627  class LinearProblem(object): Line 1665  class LinearProblem(object):
1665       """       """
1666       Switches reduced order on for equation representation.       Switches reduced order on for equation representation.
1667    
1668       @raise RuntimeError: if order reduction is altered after a coefficient has       :raise RuntimeError: if order reduction is altered after a coefficient has
1669                            been set                            been set
1670       """       """
1671       if not self.__reduce_equation_order:       if not self.__reduce_equation_order:
# Line 1641  class LinearProblem(object): Line 1679  class LinearProblem(object):
1679       """       """
1680       Switches reduced order off for equation representation.       Switches reduced order off for equation representation.
1681    
1682       @raise RuntimeError: if order reduction is altered after a coefficient has       :raise RuntimeError: if order reduction is altered after a coefficient has
1683                            been set                            been set
1684       """       """
1685       if self.__reduce_equation_order:       if self.__reduce_equation_order:
# Line 1655  class LinearProblem(object): Line 1693  class LinearProblem(object):
1693       """       """
1694       Sets order reduction state for equation representation according to flag.       Sets order reduction state for equation representation according to flag.
1695    
1696       @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
1697                    equation representation, otherwise or if flag is not present                    equation representation, otherwise or if flag is not present
1698                    order reduction is switched off                    order reduction is switched off
1699       @type flag: C{bool}       :type flag: ``bool``
1700       @raise RuntimeError: if order reduction is altered after a coefficient has       :raise RuntimeError: if order reduction is altered after a coefficient has
1701                            been set                            been set
1702       """       """
1703       if flag:       if flag:
# Line 1676  class LinearProblem(object): Line 1714  class LinearProblem(object):
1714        """        """
1715        Tests a coefficient for symmetry.        Tests a coefficient for symmetry.
1716    
1717        @param name: name of the coefficient        :param name: name of the coefficient
1718        @type name: C{str}        :type name: ``str``
1719        @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
1720                        which break the symmetry is printed.                        which break the symmetry is printed.
1721        @type verbose: C{bool}        :type verbose: ``bool``
1722        @return: True if coefficient C{name} is symmetric        :return: True if coefficient ``name`` is symmetric
1723        @rtype: C{bool}        :rtype: ``bool``
1724        """        """
1725        SMALL_TOLERANCE=util.EPSILON*10.        SMALL_TOLERANCE=util.EPSILON*10.
1726        A=self.getCoefficient(name)        A=self.getCoefficient(name)
# Line 1723  class LinearProblem(object): Line 1761  class LinearProblem(object):
1761        """        """
1762        Tests two coefficients for reciprocal symmetry.        Tests two coefficients for reciprocal symmetry.
1763    
1764        @param name0: name of the first coefficient        :param name0: name of the first coefficient
1765        @type name0: C{str}        :type name0: ``str``
1766        @param name1: name of the second coefficient        :param name1: name of the second coefficient
1767        @type name1: C{str}        :type name1: ``str``
1768        @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
1769                        which break the symmetry is printed                        which break the symmetry is printed
1770        @type verbose: C{bool}        :type verbose: ``bool``
1771        @return: True if coefficients C{name0} and C{name1} are reciprocally        :return: True if coefficients ``name0`` and ``name1`` are reciprocally
1772                 symmetric.                 symmetric.
1773        @rtype: C{bool}        :rtype: ``bool``
1774        """        """
1775        SMALL_TOLERANCE=util.EPSILON*10.        SMALL_TOLERANCE=util.EPSILON*10.
1776        B=self.getCoefficient(name0)        B=self.getCoefficient(name0)
# Line 1781  class LinearProblem(object): Line 1819  class LinearProblem(object):
1819    
1820     def getCoefficient(self,name):     def getCoefficient(self,name):
1821       """       """
1822       Returns the value of the coefficient C{name}.       Returns the value of the coefficient ``name``.
1823    
1824       @param name: name of the coefficient requested       :param name: name of the coefficient requested
1825       @type name: C{string}       :type name: ``string``
1826       @return: the value of the coefficient       :return: the value of the coefficient
1827       @rtype: L{Data<escript.Data>}       :rtype: `Data`
1828       @raise IllegalCoefficient: if C{name} is not a coefficient of the PDE       :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
1829       """       """
1830       if self.hasCoefficient(name):       if self.hasCoefficient(name):
1831           return self.__COEFFICIENTS[name].getValue()           return self.__COEFFICIENTS[name].getValue()
# Line 1796  class LinearProblem(object): Line 1834  class LinearProblem(object):
1834    
1835     def hasCoefficient(self,name):     def hasCoefficient(self,name):
1836       """       """
1837       Returns True if C{name} is the name of a coefficient.       Returns True if ``name`` is the name of a coefficient.
1838    
1839       @param name: name of the coefficient enquired       :param name: name of the coefficient enquired
1840       @type name: C{string}       :type name: ``string``
1841       @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,
1842                False otherwise                False otherwise
1843       @rtype: C{bool}       :rtype: ``bool``
1844       """       """
1845       return self.__COEFFICIENTS.has_key(name)       return self.__COEFFICIENTS.has_key(name)
1846    
1847     def createCoefficient(self, name):     def createCoefficient(self, name):
1848       """       """
1849       Creates a L{Data<escript.Data>} object corresponding to coefficient       Creates a `Data` object corresponding to coefficient
1850       C{name}.       ``name``.
1851    
1852       @return: the coefficient C{name} initialized to 0       :return: the coefficient ``name`` initialized to 0
1853       @rtype: L{Data<escript.Data>}       :rtype: `Data`
1854       @raise IllegalCoefficient: if C{name} is not a coefficient of the PDE       :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
1855       """       """
1856       if self.hasCoefficient(name):       if self.hasCoefficient(name):
1857          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 1860  class LinearProblem(object):
1860    
1861     def getFunctionSpaceForCoefficient(self,name):     def getFunctionSpaceForCoefficient(self,name):
1862       """       """
1863       Returns the L{FunctionSpace<escript.FunctionSpace>} to be used for       Returns the `FunctionSpace` to be used for
1864       coefficient C{name}.       coefficient ``name``.
1865    
1866       @param name: name of the coefficient enquired       :param name: name of the coefficient enquired
1867       @type name: C{string}       :type name: ``string``
1868       @return: the function space to be used for coefficient C{name}       :return: the function space to be used for coefficient ``name``
1869       @rtype: L{FunctionSpace<escript.FunctionSpace>}       :rtype: `FunctionSpace`
1870       @raise IllegalCoefficient: if C{name} is not a coefficient of the PDE       :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
1871       """       """
1872       if self.hasCoefficient(name):       if self.hasCoefficient(name):
1873          return self.__COEFFICIENTS[name].getFunctionSpace(self.getDomain())          return self.__COEFFICIENTS[name].getFunctionSpace(self.getDomain())
# Line 1838  class LinearProblem(object): Line 1876  class LinearProblem(object):
1876    
1877     def getShapeOfCoefficient(self,name):     def getShapeOfCoefficient(self,name):
1878       """       """
1879       Returns the shape of the coefficient C{name}.       Returns the shape of the coefficient ``name``.
1880    
1881       @param name: name of the coefficient enquired       :param name: name of the coefficient enquired
1882       @type name: C{string}       :type name: ``string``
1883       @return: the shape of the coefficient C{name}       :return: the shape of the coefficient ``name``
1884       @rtype: C{tuple} of C{int}       :rtype: ``tuple`` of ``int``
1885       @raise IllegalCoefficient: if C{name} is not a coefficient of the PDE       :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
1886       """       """
1887       if self.hasCoefficient(name):       if self.hasCoefficient(name):
1888          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 1898  class LinearProblem(object):
1898    
1899     def alteredCoefficient(self,name):     def alteredCoefficient(self,name):
1900       """       """
1901       Announces that coefficient C{name} has been changed.       Announces that coefficient ``name`` has been changed.
1902    
1903       @param name: name of the coefficient affected       :param name: name of the coefficient affected
1904       @type name: C{string}       :type name: ``string``
1905       @raise IllegalCoefficient: if C{name} is not a coefficient of the PDE       :raise IllegalCoefficient: if ``name`` is not a coefficient of the PDE
1906       @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
1907              system as constraints are applied to the solved system.              system as constraints are applied to the solved system.
1908       """       """
1909       if self.hasCoefficient(name):       if self.hasCoefficient(name):
# Line 1972  class LinearProblem(object): Line 2010  class LinearProblem(object):
2010       """       """
2011       Returns the operator of the linear problem.       Returns the operator of the linear problem.
2012    
2013       @return: the operator of the problem       :return: the operator of the problem
2014       """       """
2015       return self.getSystem()[0]       return self.getSystem()[0]
2016    
# Line 1980  class LinearProblem(object): Line 2018  class LinearProblem(object):
2018       """       """
2019       Returns the right hand side of the linear problem.       Returns the right hand side of the linear problem.
2020    
2021       @return: the right hand side of the problem       :return: the right hand side of the problem
2022       @rtype: L{Data<escript.Data>}       :rtype: `Data`
2023       """       """
2024       return self.getSystem()[1]       return self.getSystem()[1]
2025    
# Line 2015  class LinearProblem(object): Line 2053  class LinearProblem(object):
2053             self.__solution.setToZero()             self.__solution.setToZero()
2054             self.trace("Solution is reset to zero.")             self.trace("Solution is reset to zero.")
2055    
2056     def setSolution(self,u):     def setSolution(self,u, validate=True):
2057         """         """
2058         Sets the solution assuming that makes the system valid with the tolrance         Sets the solution assuming that makes the system valid with the tolrance
2059         defined by the solver options         defined by the solver options
2060         """         """
2061         self.__solution_rtol=self.getSolverOptions().getTolerance()         if validate:
2062         self.__solution_atol=self.getSolverOptions().getAbsoluteTolerance()        self.__solution_rtol=self.getSolverOptions().getTolerance()
2063          self.__solution_atol=self.getSolverOptions().getAbsoluteTolerance()
2064          self.validSolution()
2065         self.__solution=u         self.__solution=u
        self.validSolution()  
   
2066     def getCurrentSolution(self):     def getCurrentSolution(self):
2067         """         """
2068         Returns the solution in its current state.         Returns the solution in its current state.
# Line 2064  class LinearProblem(object): Line 2102  class LinearProblem(object):
2102             if self.isUsingLumping():             if self.isUsingLumping():
2103                 self.__operator.setToZero()                 self.__operator.setToZero()
2104             else:             else:
2105                 self.__operator.resetValues()                 if self.getOperatorType() == self.getRequiredOperatorType():
2106                       self.__operator.resetValues()
2107                   else:
2108                       self.__operator=self.createOperator()
2109                   self.__operator_type=self.getRequiredOperatorType()
2110             self.trace("Operator reset to zero")             self.trace("Operator reset to zero")
2111    
2112     def getCurrentOperator(self):     def getCurrentOperator(self):
# Line 2077  class LinearProblem(object): Line 2119  class LinearProblem(object):
2119        """        """
2120        Sets new values to coefficients.        Sets new values to coefficients.
2121    
2122        @raise IllegalCoefficient: if an unknown coefficient keyword is used        :raise IllegalCoefficient: if an unknown coefficient keyword is used
2123        """        """
2124        # check if the coefficients are  legal:        # check if the coefficients are  legal:
2125        for i in coefficients.iterkeys():        for i in coefficients.iterkeys():
# Line 2136  class LinearProblem(object): Line 2178  class LinearProblem(object):
2178        """        """
2179        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.
2180    
2181        @note: Typically this method is overwritten when implementing a        :note: Typically this method is overwritten when implementing a
2182               particular linear problem.               particular linear problem.
2183        """        """
2184        return None        return None
# Line 2145  class LinearProblem(object): Line 2187  class LinearProblem(object):
2187         """         """
2188         Returns an instance of a new operator.         Returns an instance of a new operator.
2189    
2190         @note: This method is overwritten when implementing a particular         :note: This method is overwritten when implementing a particular
2191                linear problem.                linear problem.
2192         """         """
2193         return escript.Operator()         return escript.Operator()
# Line 2154  class LinearProblem(object): Line 2196  class LinearProblem(object):
2196        """        """
2197        Tests the PDE for symmetry.        Tests the PDE for symmetry.
2198    
2199        @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
2200                        which break the symmetry is printed                        which break the symmetry is printed
2201        @type verbose: C{bool}        :type verbose: ``bool``
2202        @return: True if the problem is symmetric        :return: True if the problem is symmetric
2203        @rtype: C{bool}        :rtype: ``bool``
2204        @note: Typically this method is overwritten when implementing a        :note: Typically this method is overwritten when implementing a
2205               particular linear problem.               particular linear problem.
2206        """        """
2207        out=True        out=True
# Line 2169  class LinearProblem(object): Line 2211  class LinearProblem(object):
2211         """         """
2212         Returns the solution of the problem.         Returns the solution of the problem.
2213    
2214         @return: the solution         :return: the solution
2215         @rtype: L{Data<escript.Data>}         :rtype: `Data`
2216    
2217         @note: This method is overwritten when implementing a particular         :note: This method is overwritten when implementing a particular
2218                linear problem.                linear problem.
2219         """         """
2220         return self.getCurrentSolution()         return self.getCurrentSolution()
# Line 2181  class LinearProblem(object): Line 2223  class LinearProblem(object):
2223         """         """
2224         Returns the operator and right hand side of the PDE.         Returns the operator and right hand side of the PDE.
2225    
2226         @return: the discrete version of the PDE         :return: the discrete version of the PDE
2227         @rtype: C{tuple} of L{Operator,<escript.Operator>} and L{Data<escript.Data>}.         :rtype: ``tuple`` of `Operator` and `Data`.
2228    
2229         @note: This method is overwritten when implementing a particular         :note: This method is overwritten when implementing a particular
2230                linear problem.                linear problem.
2231         """         """
2232         return (self.getCurrentOperator(), self.getCurrentRightHandSide())         return (self.getCurrentOperator(), self.getCurrentRightHandSide())
# Line 2192  class LinearProblem(object): Line 2234  class LinearProblem(object):
2234  class LinearPDE(LinearProblem):  class LinearPDE(LinearProblem):
2235     """     """
2236     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
2237     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
2238     L{Domain<escript.Domain>} object.     `Domain` object.
2239    
2240     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
2241     is defined in the following form:     is defined in the following form:
2242    
2243     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)*
2244    
2245     where M{grad(F)} denotes the spatial derivative of M{F}. Einstein's     where *grad(F)* denotes the spatial derivative of *F*. Einstein's
2246     summation convention, ie. summation over indexes appearing twice in a term     summation convention, ie. summation over indexes appearing twice in a term
2247     of a sum performed, is used.     of a sum performed, is used.
2248     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
2249     through L{Data<escript.Data>} objects in L{Function<escript.Function>} and     through `Data` objects in `Function` and
2250     the coefficients M{A_reduced}, M{B_reduced}, M{C_reduced}, M{D_reduced},     the coefficients *A_reduced*, *B_reduced*, *C_reduced*, *D_reduced*,
2251     M{X_reduced} and M{Y_reduced} have to be specified through     *X_reduced* and *Y_reduced* have to be specified through
2252     L{Data<escript.Data>} objects in L{ReducedFunction<escript.ReducedFunction>}.     `Data` objects in `ReducedFunction`.
2253     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
2254     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*,
2255     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
2256     M{D}, M{D_reduced}, M{Y} and M{Y_reduced} are scalar.     *D*, *D_reduced*, *Y* and *Y_reduced* are scalar.
2257    
2258     The following natural boundary conditions are considered:     The following natural boundary conditions are considered:
2259    
2260     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*
2261    
2262     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*,
2263     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
2264     PDE. The coefficients M{d} and M{y} are each a scalar in     PDE. The coefficients *d* and *y* are each a scalar in
2265     L{FunctionOnBoundary<escript.FunctionOnBoundary>} and the coefficients     `FunctionOnBoundary` and the coefficients
2266     M{d_reduced} and M{y_reduced} are each a scalar in     *d_reduced* and *y_reduced* are each a scalar in
2267     L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.     `ReducedFunctionOnBoundary`.
2268    
2269     Constraints for the solution prescribe the value of the solution at certain     Constraints for the solution prescribe the value of the solution at certain
2270     locations in the domain. They have the form     locations in the domain. They have the form
2271    
2272     M{u=r} where M{q>0}     *u=r* where *q>0*
2273    
2274     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
2275     defining where the constraint is applied. The constraints override any     defining where the constraint is applied. The constraints override any
2276     other condition set by the PDE or the boundary condition.     other condition set by the PDE or the boundary condition.
2277    
2278     The PDE is symmetrical if     The PDE is symmetrical if
2279    
2280     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]*
2281     and M{B_reduced[j]=C_reduced[j]}     and *B_reduced[j]=C_reduced[j]*
2282    
2283     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
2284     form     form
2285    
2286     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]*
2287    
2288     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
2289     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
2290     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.
2291     The natural boundary conditions take the form:     The natural boundary conditions take the form:
2292    
2293     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]*
2294    
2295     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
2296     L{FunctionOnBoundary<escript.FunctionOnBoundary>}. The coefficients     `FunctionOnBoundary`. The coefficients
2297     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
2298     L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.     `ReducedFunctionOnBoundary`.
2299    
2300     Constraints take the form     Constraints take the form
2301    
2302     M{u[i]=r[i]}  where  M{q[i]>0}     *u[i]=r[i]*  where  *q[i]>0*
2303    
2304     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
2305     necessarily all components must have a constraint.     necessarily all components must have a constraint.
2306    
2307     The system of PDEs is symmetrical if     The system of PDEs is symmetrical if
2308    
2309        - M{A[i,j,k,l]=A[k,l,i,j]}        - *A[i,j,k,l]=A[k,l,i,j]*
2310        - 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]*
2311        - M{B[i,j,k]=C[k,i,j]}        - *B[i,j,k]=C[k,i,j]*
2312        - M{B_reduced[i,j,k]=C_reduced[k,i,j]}        - *B_reduced[i,j,k]=C_reduced[k,i,j]*
2313        - M{D[i,k]=D[i,k]}        - *D[i,k]=D[i,k]*
2314        - M{D_reduced[i,k]=D_reduced[i,k]}        - *D_reduced[i,k]=D_reduced[i,k]*
2315        - M{d[i,k]=d[k,i]}        - *d[i,k]=d[k,i]*
2316        - M{d_reduced[i,k]=d_reduced[k,i]}        - *d_reduced[i,k]=d_reduced[k,i]*
2317    
2318     L{LinearPDE} also supports solution discontinuities over a contact region     `LinearPDE` also supports solution discontinuities over a contact region
2319     in the domain. To specify the conditions across the discontinuity we are     in the domain. To specify the conditions across the discontinuity we are
2320     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
2321     and several components of the solution, is defined as     and several components of the solution, is defined as
2322    
2323     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]*
2324    
2325     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
2326    
2327     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]*
2328    
2329     In the context of discontinuities M{n} denotes the normal on the     In the context of discontinuities *n* denotes the normal on the
2330     discontinuity pointing from side 0 towards side 1 calculated from     discontinuity pointing from side 0 towards side 1 calculated from
2331     L{getNormal<escript.FunctionSpace.getNormal>} of L{FunctionOnContactZero<escript.FunctionOnContactZero>}.     `FunctionSpace.getNormal` of `FunctionOnContactZero`.
2332     For a system of PDEs the contact condition takes the form     For a system of PDEs the contact condition takes the form
2333    
2334     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]*
2335    
2336     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
2337     discontinuity, respectively. M{jump(u)}, which is the difference of the     discontinuity, respectively. *jump(u)*, which is the difference of the
2338     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
2339     discontinuity along the normal calculated by L{jump<util.jump>}.     discontinuity along the normal calculated by `jump`.
2340     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
2341     both in L{FunctionOnContactZero<escript.FunctionOnContactZero>} or     both in `FunctionOnContactZero` or
2342     L{FunctionOnContactOne<escript.FunctionOnContactOne>}.     `FunctionOnContactOne`.
2343     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*
2344     is of rank one both in L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>}     is of rank one both in `ReducedFunctionOnContactZero`
2345     or L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.     or `ReducedFunctionOnContactOne`.
2346     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
2347     condition takes the form     condition takes the form
2348    
2349     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)*
2350    
2351     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
2352     both in L{FunctionOnContactZero<escript.FunctionOnContactZero>} or     both in `FunctionOnContactZero` or
2353     L{FunctionOnContactOne<escript.FunctionOnContactOne>} and the coefficient     `FunctionOnContactOne` and the coefficient
2354     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
2355     L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or     `ReducedFunctionOnContactZero` or
2356     L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.     `ReducedFunctionOnContactOne`.
2357    
2358     Typical usage::     Typical usage::
2359    
# Line 2325  class LinearPDE(LinearProblem): Line 2367  class LinearPDE(LinearProblem):
2367       """       """
2368       Initializes a new linear PDE.       Initializes a new linear PDE.
2369    
2370       @param domain: domain of the PDE       :param domain: domain of the PDE
2371       @type domain: L{Domain<escript.Domain>}       :type domain: `Domain`
2372       @param numEquations: number of equations. If C{None} the number of       :param numEquations: number of equations. If ``None`` the number of
2373                            equations is extracted from the PDE coefficients.                            equations is extracted from the PDE coefficients.
2374       @param numSolutions: number of solution components. If C{None} the number       :param numSolutions: number of solution components. If ``None`` the number
2375                            of solution components is extracted from the PDE                            of solution components is extracted from the PDE
2376                            coefficients.                            coefficients.
2377       @param debug: if True debug information is printed       :param debug: if True debug information is printed
2378    
2379       """       """
2380       super(LinearPDE, self).__init__(domain,numEquations,numSolutions,debug)       super(LinearPDE, self).__init__(domain,numEquations,numSolutions,debug)
# Line 2367  class LinearPDE(LinearProblem): Line 2409  class LinearPDE(LinearProblem):
2409       """       """
2410       Returns the string representation of the PDE.       Returns the string representation of the PDE.
2411    
2412       @return: a simple representation of the PDE       :return: a simple representation of the PDE
2413       @rtype: C{str}       :rtype: ``str``
2414       """       """
2415       return "<LinearPDE %d>"%id(self)       return "<LinearPDE %d>"%id(self)
2416    
# Line 2383  class LinearPDE(LinearProblem): Line 2425  class LinearPDE(LinearProblem):
2425        """        """
2426        Tests the PDE for symmetry.        Tests the PDE for symmetry.
2427    
2428        @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
2429                        which break the symmetry is printed.                        which break the symmetry is printed.
2430        @type verbose: C{bool}        :type verbose: ``bool``
2431        @return: True if the PDE is symmetric        :return: True if the PDE is symmetric
2432        @rtype: L{bool}        :rtype: `bool`
2433        @note: This is a very expensive operation. It should be used for        :note: This is a very expensive operation. It should be used for
2434               degugging only! The symmetry flag is not altered.               degugging only! The symmetry flag is not altered.
2435        """        """
2436        out=True        out=True
# Line 2421  class LinearPDE(LinearProblem): Line 2463  class LinearPDE(LinearProblem):
2463         """         """
2464         Returns the solution of the PDE.         Returns the solution of the PDE.
2465    
2466         @return: the solution         :return: the solution
2467         @rtype: L{Data<escript.Data>}         :rtype: `Data`
2468         """         """
2469         option_class=self.getSolverOptions()         option_class=self.getSolverOptions()
2470         if not self.isSolutionValid():         if not self.isSolutionValid():
# Line 2439  class LinearPDE(LinearProblem): Line 2481  class LinearPDE(LinearProblem):
2481         """         """
2482         Returns the operator and right hand side of the PDE.         Returns the operator and right hand side of the PDE.
2483    
2484         @return: the discrete version of the PDE         :return: the discrete version of the PDE
2485         @rtype: C{tuple} of L{Operator,<escript.Operator>} and         :rtype: ``tuple`` of `Operator` and
2486                 L{Data<escript.Data>}                 `Data`
2487         """         """
2488         if not self.isOperatorValid() or not self.isRightHandSideValid():         if not self.isOperatorValid() or not self.isRightHandSideValid():
2489            if self.isUsingLumping():            if self.isUsingLumping():
# Line 2616  class LinearPDE(LinearProblem): Line 2658  class LinearPDE(LinearProblem):
2658        """        """
2659        Applies the constraints defined by q and r to the PDE.        Applies the constraints defined by q and r to the PDE.
2660    
2661        @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
2662                         constraint                         constraint
2663        @type rhs_only: C{bool}        :type rhs_only: ``bool``
2664        """        """
2665        q=self.getCoefficient("q")        q=self.getCoefficient("q")
2666        r=self.getCoefficient("r")        r=self.getCoefficient("r")
# Line 2646  class LinearPDE(LinearProblem): Line 2688  class LinearPDE(LinearProblem):
2688        """        """
2689        Sets new values to coefficients.        Sets new values to coefficients.
2690    
2691        @param coefficients: new values assigned to coefficients        :param coefficients: new values assigned to coefficients
2692        @keyword A: value for coefficient C{A}        :keyword A: value for coefficient ``A``
2693        @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
2694                 L{Function<escript.Function>}                 `Function`
2695        @keyword A_reduced: value for coefficient C{A_reduced}        :keyword A_reduced: value for coefficient ``A_reduced``
2696        @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`
2697                         object on L{ReducedFunction<escript.ReducedFunction>}                         object on `ReducedFunction`
2698        @keyword B: value for coefficient C{B}        :keyword B: value for coefficient ``B``
2699        @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
2700                 L{Function<escript.Function>}                 `Function`
2701        @keyword B_reduced: value for coefficient C{B_reduced}        :keyword B_reduced: value for coefficient ``B_reduced``
2702        @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`
2703                         object on L{ReducedFunction<escript.ReducedFunction>}                         object on `ReducedFunction`
2704        @keyword C: value for coefficient C{C}        :keyword C: value for coefficient ``C``
2705        @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
2706                 L{Function<escript.Function>}                 `Function`
2707        @keyword C_reduced: value for coefficient C{C_reduced}        :keyword C_reduced: value for coefficient ``C_reduced``
2708        @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`
2709                         object on L{ReducedFunction<escript.ReducedFunction>}                         object on `ReducedFunction`
2710        @keyword D: value for coefficient C{D}        :keyword D: value for coefficient ``D``
2711        @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
2712                 L{Function<escript.Function>}                 `Function`
2713        @keyword D_reduced: value for coefficient C{D_reduced}        :keyword D_reduced: value for coefficient ``D_reduced``
2714        @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`
2715                         object on L{ReducedFunction<escript.ReducedFunction>}                         object on `ReducedFunction`
2716        @keyword X: value for coefficient C{X}        :keyword X: value for coefficient ``X``
2717        @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
2718                 L{Function<escript.Function>}                 `Function`
2719        @keyword X_reduced: value for coefficient C{X_reduced}        :keyword X_reduced: value for coefficient ``X_reduced``
2720        @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`
2721                         object on L{ReducedFunction<escript.ReducedFunction>}                         object on `ReducedFunction`
2722        @keyword Y: value for coefficient C{Y}        :keyword Y: value for coefficient ``Y``
2723        @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
2724                 L{Function<escript.Function>}                 `Function`
2725        @keyword Y_reduced: value for coefficient C{Y_reduced}        :keyword Y_reduced: value for coefficient ``Y_reduced``
2726        @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`
2727                         object on L{ReducedFunction<escript.Function>}                         object on `ReducedFunction`
2728        @keyword d: value for coefficient C{d}        :keyword d: value for coefficient ``d``
2729        @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
2730                 L{FunctionOnBoundary<escript.FunctionOnBoundary>}                 `FunctionOnBoundary`
2731        @keyword d_reduced: value for coefficient C{d_reduced}        :keyword d_reduced: value for coefficient ``d_reduced``
2732        @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`
2733                         object on L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}                         object on `ReducedFunctionOnBoundary`
2734        @keyword y: value for coefficient C{y}        :keyword y: value for coefficient ``y``
2735        @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
2736                 L{FunctionOnBoundary<escript.FunctionOnBoundary>}                 `FunctionOnBoundary`
2737        @keyword d_contact: value for coefficient C{d_contact}        :keyword d_contact: value for coefficient ``d_contact``
2738        @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`
2739                         object on L{FunctionOnContactOne<escript.FunctionOnContactOne>}                         object on `FunctionOnContactOne`
2740                         or L{FunctionOnContactZero<escript.FunctionOnContactZero>}                         or `FunctionOnContactZero`
2741        @keyword d_contact_reduced: value for coefficient C{d_contact_reduced}        :keyword d_contact_reduced: value for coefficient ``d_contact_reduced``
2742        @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`
2743                                 object on L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}                                 object on `ReducedFunctionOnContactOne`
2744                                 or L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>}                                 or `ReducedFunctionOnContactZero`
2745        @keyword y_contact: value for coefficient C{y_contact}        :keyword y_contact: value for coefficient ``y_contact``
2746        @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`
2747                         object on L{FunctionOnContactOne<escript.FunctionOnContactOne>}                         object on `FunctionOnContactOne`
2748                         or L{FunctionOnContactZero<escript.FunctionOnContactZero>}                         or `FunctionOnContactZero`
2749        @keyword y_contact_reduced: value for coefficient C{y_contact_reduced}        :keyword y_contact_reduced: value for coefficient ``y_contact_reduced``
2750        @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`
2751                                 object on L{ReducedFunctionOnContactOne<escript.FunctionOnContactOne>}                                 object on `ReducedFunctionOnContactOne`
2752                                 or L{ReducedFunctionOnContactZero<escript.FunctionOnContactZero>}                                 or `ReducedFunctionOnContactZero`
2753        @keyword r: values prescribed to the solution at the locations of        :keyword r: values prescribed to the solution at the locations of
2754                    constraints                    constraints
2755        @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
2756                 L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}                 `Solution` or `ReducedSolution`
2757                 depending on whether reduced order is used for the solution                 depending on whether reduced order is used for the solution
2758        @keyword q: mask for location of constraints        :keyword q: mask for location of constraints
2759        @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
2760                 L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}                 `Solution` or `ReducedSolution`
2761                 depending on whether reduced order is used for the                 depending on whether reduced order is used for the
2762                 representation of the equation                 representation of the equation
2763        @raise IllegalCoefficient: if an unknown coefficient keyword is used        :raise IllegalCoefficient: if an unknown coefficient keyword is used
2764        """        """
2765        super(LinearPDE,self).setValue(**coefficients)        super(LinearPDE,self).setValue(**coefficients)
2766        # check if the systrem is inhomogeneous:        # check if the systrem is inhomogeneous:
# Line 2735  class LinearPDE(LinearProblem): Line 2777  class LinearPDE(LinearProblem):
2777       """       """
2778       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.
2779    
2780       @param u: argument in the residual calculation. It must be representable       :param u: argument in the residual calculation. It must be representable
2781                 in L{self.getFunctionSpaceForSolution()}. If u is not present                 in `self.getFunctionSpaceForSolution()`. If u is not present
2782                 or equals C{None} the current solution is used.                 or equals ``None`` the current solution is used.
2783       @type u: L{Data<escript.Data>} or None       :type u: `Data` or None
2784       @return: residual of u       :return: residual of u
2785       @rtype: L{Data<escript.Data>}       :rtype: `Data`
2786       """       """
2787       if u==None:       if u==None:
2788          return self.getOperator()*self.getSolution()-self.getRightHandSide()          return self.getOperator()*self.getSolution()-self.getRightHandSide()
# Line 2749  class LinearPDE(LinearProblem): Line 2791  class LinearPDE(LinearProblem):
2791    
2792     def getFlux(self,u=None):     def getFlux(self,u=None):
2793       """       """
2794       Returns the flux M{J} for a given M{u}.       Returns the flux *J* for a given *u*.
2795    
2796       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]*
2797    
2798       or       or
2799    
2800       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]*
2801    
2802       @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
2803                 current solution is used.                 current solution is used.
2804       @type u: L{Data<escript.Data>} or None       :type u: `Data` or None
2805       @return: flux       :return: flux
2806       @rtype: L{Data<escript.Data>}       :rtype: `Data`
2807       """       """
2808       if u==None: u=self.getSolution()       if u==None: u=self.getSolution()
2809       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 2817  class LinearPDE(LinearProblem):
2817  class Poisson(LinearPDE):  class Poisson(LinearPDE):
2818     """     """
2819     Class to define a Poisson equation problem. This is generally a     Class to define a Poisson equation problem. This is generally a
2820     L{LinearPDE} of the form     `LinearPDE` of the form
2821    
2822     M{-grad(grad(u)[j])[j] = f}     *-grad(grad(u)[j])[j] = f*
2823    
2824     with natural boundary conditions     with natural boundary conditions
2825    
2826     M{n[j]*grad(u)[j] = 0 }     *n[j]*grad(u)[j] = 0*
2827    
2828     and constraints:     and constraints:
2829    
2830     M{u=0} where M{q>0}     *u=0* where *q>0*
2831    
2832     """     """
2833    
# Line 2793  class Poisson(LinearPDE): Line 2835  class Poisson(LinearPDE):
2835       """       """
2836       Initializes a new Poisson equation.       Initializes a new Poisson equation.
2837    
2838       @param domain: domain of the PDE       :param domain: domain of the PDE
2839       @type domain: L{Domain<escript.Domain>}       :type domain: `Domain`
2840       @param debug: if True debug information is printed       :param debug: if True debug information is printed
2841    
2842       """       """
2843       super(Poisson, self).__init__(domain,1,1,debug)       super(Poisson, self).__init__(domain,1,1,debug)
# Line 2808  class Poisson(LinearPDE): Line 2850  class Poisson(LinearPDE):
2850       """       """
2851       Sets new values to coefficients.       Sets new values to coefficients.
2852    
2853       @param coefficients: new values assigned to coefficients       :param coefficients: new values assigned to coefficients
2854       @keyword f: value for right hand side M{f}       :keyword f: value for right hand side *f*
2855       @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
2856                on L{Function<escript.Function>}                on `Function`
2857       @keyword q: mask for location of constraints       :keyword q: mask for location of constraints
2858       @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`
2859                object on L{Solution<escript.Solution>} or                object on `Solution` or
2860                L{ReducedSolution<escript.ReducedSolution>} depending on whether                `ReducedSolution` depending on whether
2861                reduced order is used for the representation of the equation                reduced order is used for the representation of the equation
2862       @raise IllegalCoefficient: if an unknown coefficient keyword is used       :raise IllegalCoefficient: if an unknown coefficient keyword is used
2863       """       """
2864       super(Poisson, self).setValue(**coefficients)       super(Poisson, self).setValue(**coefficients)
2865    
2866    
2867     def getCoefficient(self,name):     def getCoefficient(self,name):
2868       """       """
2869       Returns the value of the coefficient C{name} of the general PDE.       Returns the value of the coefficient ``name`` of the general PDE.
2870    
2871       @param name: name of the coefficient requested       :param name: name of the coefficient requested
2872       @type name: C{string}       :type name: ``string``
2873       @return: the value of the coefficient C{name}       :return: the value of the coefficient ``name``
2874       @rtype: L{Data<escript.Data>}       :rtype: `Data`
2875       @raise IllegalCoefficient: invalid coefficient name       :raise IllegalCoefficient: invalid coefficient name
2876       @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
2877              equation onto the general PDE.              equation onto the general PDE.
2878       """       """
2879       if name == "A" :       if name == "A" :
# Line 2846  class Poisson(LinearPDE): Line 2888  class Poisson(LinearPDE):
2888  class Helmholtz(LinearPDE):  class Helmholtz(LinearPDE):
2889     """     """
2890     Class to define a Helmholtz equation problem. This is generally a     Class to define a Helmholtz equation problem. This is generally a
2891     L{LinearPDE} of the form     `LinearPDE` of the form
2892    
2893     M{S{omega}*u - grad(k*grad(u)[j])[j] = f}     *omega*u - grad(k*grad(u)[j])[j] = f*
2894    
2895     with natural boundary conditions     with natural boundary conditions
2896    
2897     M{k*n[j]*grad(u)[j] = g- S{alpha}u }     *k*n[j]*grad(u)[j] = g- alphau*
2898    
2899     and constraints:     and constraints:
2900    
2901     M{u=r} where M{q>0}     *u=r* where *q>0*
2902    
2903     """     """
2904    
# Line 2864  class Helmholtz(LinearPDE): Line 2906  class Helmholtz(LinearPDE):
2906       """       """
2907       Initializes a new Helmholtz equation.       Initializes a new Helmholtz equation.
2908    
2909       @param domain: domain of the PDE       :param domain: domain of the PDE
2910       @type domain: L{Domain<escript.Domain>}       :type domain: `Domain`
2911       @param debug: if True debug information is printed       :param debug: if True debug information is printed
2912    
2913       """       """
2914       super(Helmholtz, self).__init__(domain,1,1,debug)       super(Helmholtz, self).__init__(domain,1,1,debug)
# Line 2884  class Helmholtz(LinearPDE): Line 2926  class Helmholtz(LinearPDE):
2926       """       """
2927       Sets new values to coefficients.       Sets new values to coefficients.
2928    
2929       @param coefficients: new values assigned to coefficients       :param coefficients: new values assigned to coefficients
2930       @keyword omega: value for coefficient M{S{omega}}       :keyword omega: value for coefficient *omega*
2931       @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`
2932                    object on L{Function<escript.Function>}                    object on `Function`
2933       @keyword k: value for coefficient M{k}       :keyword k: value for coefficient *k*
2934       @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
2935                on L{Function<escript.Function>}                on `Function`
2936       @keyword f: value for right hand side M{f}       :keyword f: value for right hand side *f*
2937       @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
2938                on L{Function<escript.Function>}                on `Function`
2939       @keyword alpha: value for right hand side M{S{alpha}}       :keyword alpha: value for right hand side *alpha*
2940       @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`
2941                    object on L{FunctionOnBoundary<escript.FunctionOnBoundary>}                    object on `FunctionOnBoundary`
2942       @keyword g: value for right hand side M{g}       :keyword g: value for right hand side *g*
2943       @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
2944                on L{FunctionOnBoundary<escript.FunctionOnBoundary>}                on `FunctionOnBoundary`
2945       @keyword r: prescribed values M{r} for the solution in constraints       :keyword r: prescribed values *r* for the solution in constraints
2946       @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
2947                on L{Solution<escript.Solution>} or                on `Solution` or
2948                L{ReducedSolution<escript.ReducedSolution>} depending on whether                `ReducedSolution` depending on whether
2949                reduced order is used for the representation of the equation                reduced order is used for the representation of the equation
2950       @keyword q: mask for the location of constraints       :keyword q: mask for the location of constraints
2951       @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
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       @raise IllegalCoefficient: if an unknown coefficient keyword is used       :raise IllegalCoefficient: if an unknown coefficient keyword is used
2956       """       """
2957       super(Helmholtz, self).setValue(**coefficients)       super(Helmholtz, self).setValue(**coefficients)
2958    
2959     def getCoefficient(self,name):     def getCoefficient(self,name):
2960       """       """
2961       Returns the value of the coefficient C{name} of the general PDE.       Returns the value of the coefficient ``name`` of the general PDE.
2962    
2963       @param name: name of the coefficient requested       :param name: name of the coefficient requested
2964       @type name: C{string}       :type name: ``string``
2965       @return: the value of the coefficient C{name}       :return: the value of the coefficient ``name``
2966       @rtype: L{Data<escript.Data>}       :rtype: `Data`
2967       @raise IllegalCoefficient: invalid name       :raise IllegalCoefficient: invalid name
2968       """       """
2969       if name == "A" :       if name == "A" :
2970           if self.getCoefficient("k").isEmpty():           if self.getCoefficient("k").isEmpty():
# Line 2948  class LameEquation(LinearPDE): Line 2990  class LameEquation(LinearPDE):
2990     """     """
2991     Class to define a Lame equation problem. This problem is defined as:     Class to define a Lame equation problem. This problem is defined as:
2992    
2993     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]*
2994    
2995     with natural boundary conditions:     with natural boundary conditions:
2996    
2997     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]*
2998    
2999     and constraints:     and constraints:
3000    
3001     M{u[i]=r[i]} where M{q[i]>0}     *u[i]=r[i]* where *q[i]>0*
3002    
3003     """     """
3004    
# Line 2964  class LameEquation(LinearPDE): Line 3006  class LameEquation(LinearPDE):
3006        """        """
3007        Initializes a new Lame equation.        Initializes a new Lame equation.
3008    
3009        @param domain: domain of the PDE        :param domain: domain of the PDE
3010        @type domain: L{Domain<escript.Domain>}        :type domain: `Domain`
3011        @param debug: if True debug information is printed        :param debug: if True debug information is printed
3012    
3013        """        """
3014        super(LameEquation, self).__init__(domain,\        super(LameEquation, self).__init__(domain,\
# Line 2982  class LameEquation(LinearPDE): Line 3024  class LameEquation(LinearPDE):
3024       """       """
3025       Sets new values to coefficients.       Sets new values to coefficients.
3026    
3027       @param coefficients: new values assigned to coefficients       :param coefficients: new values assigned to coefficients
3028       @keyword lame_mu: value for coefficient M{S{mu}}       :keyword lame_mu: value for coefficient *mu*
3029       @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`
3030                      object on L{Function<escript.Function>}                      object on `Function`
3031       @keyword lame_lambda: value for coefficient M{S{lambda}}       :keyword lame_lambda: value for coefficient *lambda*
3032       @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`
3033                          object on L{Function<escript.Function>}                          object on `Function`
3034       @keyword F: value for internal force M{F}       :keyword F: value for internal force *F*
3035       @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
3036                on L{Function<escript.Function>}                on `Function`
3037       @keyword sigma: value for initial stress M{S{sigma}}       :keyword sigma: value for initial stress *sigma*
3038       @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`
3039                    object on L{Function<escript.Function>}                    object on `Function`
3040       @keyword f: value for external force M{f}       :keyword f: value for external force *f*
3041       @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
3042                on L{FunctionOnBoundary<escript.FunctionOnBoundary>}                on `FunctionOnBoundary`
3043       @keyword r: prescribed values M{r} for the solution in constraints       :keyword r: prescribed values *r* for the solution in constraints
3044       @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
3045                on L{Solution<escript.Solution>} or                on `Solution` or
3046                L{ReducedSolution<escript.ReducedSolution>} depending on whether                `ReducedSolution` depending on whether
3047                reduced order is used for the representation of the equation                reduced order is used for the representation of the equation
3048       @keyword q: mask for the location of constraints       :keyword q: mask for the location of constraints
3049       @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
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       @raise IllegalCoefficient: if an unknown coefficient keyword is used       :raise IllegalCoefficient: if an unknown coefficient keyword is used
3054       """       """
3055       super(LameEquation, self).setValues(**coefficients)       super(LameEquation, self).setValues(**coefficients)
3056    
3057     def getCoefficient(self,name):     def getCoefficient(self,name):
3058       """       """
3059       Returns the value of the coefficient C{name} of the general PDE.       Returns the value of the coefficient ``name`` of the general PDE.
3060    
3061       @param name: name of the coefficient requested       :param name: name of the coefficient requested
3062       @type name: C{string}       :type name: ``string``
3063       @return: the value of the coefficient C{name}       :return: the value of the coefficient ``name``
3064       @rtype: L{Data<escript.Data>}       :rtype: `Data`
3065       @raise IllegalCoefficient: invalid coefficient name       :raise IllegalCoefficient: invalid coefficient name
3066       """       """
3067       out =self.createCoefficient("A")       out =self.createCoefficient("A")
3068       if name == "A" :       if name == "A" :
# Line 3057  def LinearSinglePDE(domain,debug=False): Line 3099  def LinearSinglePDE(domain,debug=False):
3099     """     """
3100     Defines a single linear PDE.     Defines a single linear PDE.
3101    
3102     @param domain: domain of the PDE     :param domain: domain of the PDE
3103     @type domain: L{Domain<escript.Domain>}     :type domain: `Domain`
3104     @param debug: if True debug information is printed     :param debug: if True debug information is printed
3105     @rtype: L{LinearPDE}     :rtype: `LinearPDE`
3106     """     """
3107     return LinearPDE(domain,numEquations=1,numSolutions=1,debug=debug)     return LinearPDE(domain,numEquations=1,numSolutions=1,debug=debug)
3108    
# Line 3068  def LinearPDESystem(domain,debug=False): Line 3110  def LinearPDESystem(domain,debug=False):
3110     """     """
3111     Defines a system of linear PDEs.     Defines a system of linear PDEs.
3112    
3113     @param domain: domain of the PDEs     :param domain: domain of the PDEs
3114     @type domain: L{Domain<escript.Domain>}     :type domain: `Domain`
3115     @param debug: if True debug information is printed     :param debug: if True debug information is printed
3116     @rtype: L{LinearPDE}     :rtype: `LinearPDE`
3117     """     """
3118     return LinearPDE(domain,numEquations=domain.getDim(),numSolutions=domain.getDim(),debug=debug)     return LinearPDE(domain,numEquations=domain.getDim(),numSolutions=domain.getDim(),debug=debug)
3119    
# Line 3080  class TransportPDE(LinearProblem): Line 3122  class TransportPDE(LinearProblem):
3122     """     """
3123     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,
3124     time dependent, second order PDE for an unknown, non-negative,     time dependent, second order PDE for an unknown, non-negative,
3125     time-dependent function M{u} on a given domain defined through a     time-dependent function *u* on a given domain defined through a
3126     L{Domain<escript.Domain>} object.     `Domain` object.
3127    
3128     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
3129     problem is defined in the following form:     problem is defined in the following form:
3130    
3131     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)*
3132    
3133     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
3134     spatial derivative of M{F}. Einstein's summation convention,  ie. summation     spatial derivative of *F*. Einstein's summation convention,  ie. summation
3135     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.
3136     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
3137     specified through L{Data<escript.Data>} objects in L{Function<escript.Function>}     specified through `Data` objects in `Function`
3138     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*,
3139     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
3140     L{Data<escript.Data>} objects in L{ReducedFunction<escript.ReducedFunction>}.     `Data` objects in `ReducedFunction`.
3141     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
3142     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
3143     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
3144     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
3145     scalar.     scalar.
3146    
3147     The following natural boundary conditions are considered:     The following natural boundary conditions are considered:
3148    
3149     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*
3150    
3151     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*,
3152     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
3153     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
3154     L{FunctionOnBoundary<escript.FunctionOnBoundary>} and the coefficients     `FunctionOnBoundary` and the coefficients
3155     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
3156     L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.     `ReducedFunctionOnBoundary`.
3157    
3158     Constraints for the solution prescribing the value of the solution at     Constraints for the solution prescribing the value of the solution at
3159     certain locations in the domain have the form     certain locations in the domain have the form
3160    
3161     M{u_t=r} where M{q>0}     *u_t=r* where *q>0*
3162    
3163     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
3164     defining where the constraint is applied. The constraints override any other     defining where the constraint is applied. The constraints override any other
3165     condition set by the transport problem or the boundary condition.     condition set by the transport problem or the boundary condition.
3166    
3167     The transport problem is symmetrical if     The transport problem is symmetrical if
3168    
3169     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
3170     M{B_reduced[j]=C_reduced[j]}     *B_reduced[j]=C_reduced[j]*
3171    
3172     For a system and a solution with several components the transport problem     For a system and a solution with several components the transport problem
3173     has the form     has the form
3174    
3175     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]*
3176    
3177     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
3178     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*,
3179     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
3180     rank one. The natural boundary conditions take the form:     rank one. The natural boundary conditions take the form:
3181    
3182     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*
3183    
3184     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
3185     all in L{FunctionOnBoundary<escript.FunctionOnBoundary>}. The coefficients     all in `FunctionOnBoundary`. The coefficients
3186     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
3187     one all in L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}.     one all in `ReducedFunctionOnBoundary`.
3188    
3189     Constraints take the form     Constraints take the form
3190    
3191     M{u[i]_t=r[i]} where M{q[i]>0}     *u[i]_t=r[i]* where *q[i]>0*
3192    
3193     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
3194     necessarily all components must have a constraint.     necessarily all components must have a constraint.
3195    
3196     The transport problem is symmetrical if     The transport problem is symmetrical if
3197    
3198        - M{M[i,k]=M[i,k]}        - *M[i,k]=M[i,k]*
3199        - M{M_reduced[i,k]=M_reduced[i,k]}        - *M_reduced[i,k]=M_reduced[i,k]*
3200        - M{A[i,j,k,l]=A[k,l,i,j]}        - *A[i,j,k,l]=A[k,l,i,j]*
3201        - 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]*
3202        - M{B[i,j,k]=C[k,i,j]}        - *B[i,j,k]=C[k,i,j]*
3203        - M{B_reduced[i,j,k]=C_reduced[k,i,j]}        - *B_reduced[i,j,k]=C_reduced[k,i,j]*
3204        - M{D[i,k]=D[i,k]}        - *D[i,k]=D[i,k]*
3205        - M{D_reduced[i,k]=D_reduced[i,k]}        - *D_reduced[i,k]=D_reduced[i,k]*
3206        - M{m[i,k]=m[k,i]}        - *m[i,k]=m[k,i]*
3207        - M{m_reduced[i,k]=m_reduced[k,i]}        - *m_reduced[i,k]=m_reduced[k,i]*
3208        - M{d[i,k]=d[k,i]}        - *d[i,k]=d[k,i]*
3209        - M{d_reduced[i,k]=d_reduced[k,i]}        - *d_reduced[i,k]=d_reduced[k,i]*
3210    
3211     L{TransportPDE} also supports solution discontinuities over a contact region     `TransportPDE` also supports solution discontinuities over a contact region
3212     in the domain. To specify the conditions across the discontinuity we are     in the domain. To specify the conditions across the discontinuity we are
3213     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
3214     several components of the solution, is defined as     several components of the solution, is defined as
3215    
3216     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]*
3217    
3218     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
3219    
3220     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]*
3221    
3222     In the context of discontinuities M{n} denotes the normal on the     In the context of discontinuities *n* denotes the normal on the
3223     discontinuity pointing from side 0 towards side 1 calculated from     discontinuity pointing from side 0 towards side 1 calculated from
3224     L{getNormal<escript.FunctionSpace.getNormal>} of L{FunctionOnContactZero<escript.FunctionOnContactZero>}.     `FunctionSpace.getNormal` of `FunctionOnContactZero`.
3225     For a system of transport problems the contact condition takes the form     For a system of transport problems the contact condition takes the form
3226    
3227     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]*
3228    
3229     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
3230     discontinuity, respectively. M{jump(u)}, which is the difference of the     discontinuity, respectively. *jump(u)*, which is the difference of the
3231     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
3232     discontinuity along the normal calculated by L{jump<util.jump>}.     discontinuity along the normal calculated by `jump`.
3233     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
3234     both in L{FunctionOnContactZero<escript.FunctionOnContactZero>} or L{FunctionOnContactOne<escript.FunctionOnContactOne>}.     both in `FunctionOnContactZero` or `FunctionOnContactOne`.
3235     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*
3236     is of rank one both in L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.     is of rank one both in `ReducedFunctionOnContactZero` or `ReducedFunctionOnContactOne`.
3237     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
3238     condition takes the form     condition takes the form
3239    
3240     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)*
3241    
3242     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
3243     both in L{FunctionOnContactZero<escript.FunctionOnContactZero>} or     both in `FunctionOnContactZero` or
3244     L{FunctionOnContactOne<escript.FunctionOnContactOne>} and the coefficient     `FunctionOnContactOne` and the coefficient
3245     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
3246     L{ReducedFunctionOnContactZero<escript.ReducedFunctionOnContactZero>} or     `ReducedFunctionOnContactZero` or
3247     L{ReducedFunctionOnContactOne<escript.ReducedFunctionOnContactOne>}.     `ReducedFunctionOnContactOne`.
3248    
3249     Typical usage::     Typical usage::
3250    
# Line 3219  class TransportPDE(LinearProblem): Line 3261  class TransportPDE(LinearProblem):
3261       """       """
3262       Initializes a transport problem.       Initializes a transport problem.
3263    
3264       @param domain: domain of the PDE       :param domain: domain of the PDE
3265       @type domain: L{Domain<escript.Domain>}       :type domain: `Domain`
3266       @param numEquations: number of equations. If C{None} the number of       :param numEquations: number of equations. If ``None`` the number of
3267                            equations is extracted from the coefficients.                            equations is extracted from the coefficients.
3268       @param numSolutions: number of solution components. If C{None} the number       :param numSolutions: number of solution components. If ``None`` the number
3269                            of solution components is extracted from the                            of solution components is extracted from the
3270                            coefficients.                            coefficients.
3271       @param debug: if True debug information is printed       :param debug: if True debug information is printed
3272       @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-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.
3273       @type useBackwardEuler: C{bool}       :type useBackwardEuler: ``bool``
3274       """       """
3275       if useBackwardEuler:       if useBackwardEuler:
3276           self.__useBackwardEuler=True           self.__useBackwardEuler=True
# Line 3272  class TransportPDE(LinearProblem): Line 3314  class TransportPDE(LinearProblem):
3314       """       """
3315       Returns the string representation of the transport problem.       Returns the string representation of the transport problem.
3316    
3317       @return: a simple representation of the transport problem       :return: a simple representation of the transport problem
3318       @rtype: C{str}       :rtype: ``str``
3319       """       """
3320       return "<TransportPDE %d>"%id(self)       return "<TransportPDE %d>"%id(self)
3321    
3322     def useBackwardEuler(self):     def useBackwardEuler(self):
3323        """        """
3324        Returns true if backward Euler is used. Otherwise false is returned.        Returns true if backward Euler is used. Otherwise false is returned.
3325        @rtype: bool        :rtype: bool
3326        """        """
3327        return self.__useBackwardEuler        return self.__useBackwardEuler
3328    
# Line 3289  class TransportPDE(LinearProblem): Line 3331  class TransportPDE(LinearProblem):
3331        """        """
3332        Tests the transport problem for symmetry.        Tests the transport problem for symmetry.
3333    
3334        @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
3335                        which break the symmetry is printed.                        which break the symmetry is printed.
3336        @type verbose: C{bool}        :type verbose: ``bool``
3337        @return:  True if the PDE is symmetric        :return:  True if the PDE is symmetric
3338        @rtype: C{bool}        :rtype: ``bool``
3339        @note: This is a very expensive operation. It should be used for        :note: This is a very expensive operation. It should be used for
3340               degugging only! The symmetry flag is not altered.               degugging only! The symmetry flag is not altered.
3341        """        """
3342        out=True        out=True
# Line 3318  class TransportPDE(LinearProblem): Line 3360  class TransportPDE(LinearProblem):
3360        """        """
3361        Sets new values to coefficients.        Sets new values to coefficients.
3362    
3363        @param coefficients: new values assigned to coefficients        :param coefficients: new values assigned to coefficients
3364        @keyword M: value for coefficient C{M}        :keyword M: value for coefficient ``M``
3365        @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
3366                 L{Function<escript.Function>}                 `Function`
3367        @keyword M_reduced: value for coefficient C{M_reduced}        :keyword M_reduced: value for coefficient ``M_reduced``
3368        @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`
3369                         object on L{Function<escript.ReducedFunction>}                         object on `Function`
3370        @keyword A: value for coefficient C{A}        :keyword A: value for coefficient ``A``
3371        @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
3372                 L{Function<escript.Function>}                 `Function`
3373        @keyword A_reduced: value for coefficient C{A_reduced}        :keyword A_reduced: value for coefficient ``A_reduced``
3374        @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`
3375                         object on L{ReducedFunction<escript.ReducedFunction>}                         object on `ReducedFunction`
3376        @keyword B: value for coefficient C{B}        :keyword B: value for coefficient ``B``
3377        @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
3378                 L{Function<escript.Function>}                 `Function`
3379        @keyword B_reduced: value for coefficient C{B_reduced}        :keyword B_reduced: value for coefficient ``B_reduced``
3380        @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`
3381                         object on L{ReducedFunction<escript.ReducedFunction>}                         object on `ReducedFunction`
3382        @keyword C: value for coefficient C{C}        :keyword C: value for coefficient ``C``
3383        @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
3384                 L{Function<escript.Function>}                 `Function`
3385        @keyword C_reduced: value for coefficient C{C_reduced}        :keyword C_reduced: value for coefficient ``C_reduced``
3386        @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`
3387                         object on L{ReducedFunction<escript.ReducedFunction>}                         object on `ReducedFunction`
3388        @keyword D: value for coefficient C{D}        :keyword D: value for coefficient ``D``
3389        @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
3390                 L{Function<escript.Function>}                 `Function`
3391        @keyword D_reduced: value for coefficient C{D_reduced}        :keyword D_reduced: value for coefficient ``D_reduced``
3392        @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`
3393                         object on L{ReducedFunction<escript.ReducedFunction>}                         object on `ReducedFunction`
3394        @keyword X: value for coefficient C{X}        :keyword X: value for coefficient ``X``
3395        @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
3396                 L{Function<escript.Function>}                 `Function`
3397        @keyword X_reduced: value for coefficient C{X_reduced}        :keyword X_reduced: value for coefficient ``X_reduced``
3398        @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`
3399                         object on L{ReducedFunction<escript.ReducedFunction>}                         object on `ReducedFunction`
3400        @keyword Y: value for coefficient C{Y}        :keyword Y: value for coefficient ``Y``
3401        @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
3402                 L{Function<escript.Function>}                 `Function`
3403        @keyword Y_reduced: value for coefficient C{Y_reduced}        :keyword Y_reduced: value for coefficient ``Y_reduced``
3404        @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`
3405                         object on L{ReducedFunction<escript.Function>}                         object on `ReducedFunction`
3406        @keyword m: value for coefficient C{m}        :keyword m: value for coefficient ``m``
3407        @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
3408                 L{FunctionOnBoundary<escript.FunctionOnBoundary>}                 `FunctionOnBoundary`
3409        @keyword m_reduced: value for coefficient C{m_reduced}        :keyword m_reduced: value for coefficient ``m_reduced``
3410        @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`
3411                         object on L{FunctionOnBoundary<escript.ReducedFunctionOnBoundary>}                         object on `FunctionOnBoundary`
3412        @keyword d: value for coefficient C{d}        :keyword d: value for coefficient ``d``
3413        @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
3414                 L{FunctionOnBoundary<escript.FunctionOnBoundary>}                 `FunctionOnBoundary`
3415        @keyword d_reduced: value for coefficient C{d_reduced}        :keyword d_reduced: value for coefficient ``d_reduced``
3416        @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`
3417                         object on L{ReducedFunctionOnBoundary<escript.ReducedFunctionOnBoundary>}                         object on `ReducedFunctionOnBoundary`
3418        @keyword y: value for coefficient C{y}        :keyword y: value for coefficient ``y``
3419        @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
3420                 L{FunctionOnBoundary<escript.FunctionOnBoundary>}                 `FunctionOnBoundary`
3421        @keyword d_contact: value for coefficient C{d_contact}        :keyword d_contact: value for coefficient ``d_contact``
3422        @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`
3423                         object on L{FunctionOnContactOne<escript.FunctionOnContactOne>} or L{FunctionOnContactZero<escript.FunctionOnContactZero>}                         object on `FunctionOnContactOne` or `FunctionOnContactZero`
3424        @keyword d_contact_reduced: value for coefficient C{d_contact_reduced}        :keyword d_contact_reduced: value for coefficient ``d_contact_reduced``
3425        @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`
3426        @keyword y_contact: value for coefficient C{y_contact}        :keyword y_contact: value for coefficient ``y_contact``
3427        @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`
3428                         object on L{FunctionOnContactOne<escript.FunctionOnContactOne>} or L{FunctionOnContactZero<escript.FunctionOnContactZero>}                         object on `FunctionOnContactOne` or `FunctionOnContactZero`
3429        @keyword y_contact_reduced: value for coefficient C{y_contact_reduced}        :keyword y_contact_reduced: value for coefficient ``y_contact_reduced``
3430        @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`
3431        @keyword r: values prescribed to the solution at the locations of constraints        :keyword r: values prescribed to the solution at the locations of constraints
3432        @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
3433                 L{Solution<escript.Solution>} or L{ReducedSolution<escript.ReducedSolution>}                 `Solution` or `ReducedSolution`
3434                 depending on whether reduced order is used for the solution                 depending on whether reduced order is used for the solution
3435        @keyword q: mask for the location of constraints        :keyword q: mask for the location of constraints
3436        @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
3437                 L{Solution<escript.Solution>} or                 `Solution` or
3438                 L{ReducedSolution<escript.ReducedSolution>} depending on whether                 `ReducedSolution` depending on whether
3439                 reduced order is used for the representation of the equation                 reduced order is used for the representation of the equation
3440        @raise IllegalCoefficient: if an unknown coefficient keyword is used        :raise IllegalCoefficient: if an unknown coefficient keyword is used
3441        """        """
3442        super(TransportPDE,self).setValue(**coefficients)        super(TransportPDE,self).setValue(**coefficients)
3443    
# Line 3403  class TransportPDE(LinearProblem): Line 3445  class TransportPDE(LinearProblem):
3445         """         """
3446         Returns an instance of a new transport operator.         Returns an instance of a new transport operator.
3447         """         """
        if self.useBackwardEuler():  
          theta=1.  
        else:  
          theta=0.5  
3448         optype=self.getRequiredOperatorType()         optype=self.getRequiredOperatorType()
3449         self.trace("New Transport problem pf type %s is allocated."%optype)         self.trace("New Transport problem pf type %s is allocated."%optype)
3450         return self.getDomain().newTransportProblem( \         return self.getDomain().newTransportProblem( \
3451                                 theta,                                 self.useBackwardEuler(),
3452                                 self.getNumEquations(), \                                 self.getNumEquations(), \
3453                                 self.getFunctionSpaceForSolution(), \                                 self.getFunctionSpaceForSolution(), \
3454                                 optype)                                 optype)
3455    
    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)  
3456    
3457     def getRequiredOperatorType(self):     def getRequiredOperatorType(self):
3458        """        """
3459        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.
3460    
3461        @return: a code to indicate the type of transport problem scheme used        :return: a code to indicate the type of transport problem scheme used
3462        @rtype: C{float}        :rtype: ``float``
3463        """        """
3464        solver_options=self.getSolverOptions()        solver_options=self.getSolverOptions()
3465        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())
3466    
3467     def getUnlimitedTimeStepSize(self):     def getUnlimitedTimeStepSize(self):
3468        """        """
3469        Returns the value returned by the C{getSafeTimeStepSize} method to        Returns the value returned by the ``getSafeTimeStepSize`` method to
3470        indicate no limit on the safe time step size.        indicate no limit on the safe time step size.
3471    
3472         @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
3473                  step size                  step size
3474         @rtype: C{float}         :rtype: ``float``
3475         @note: Typically the biggest positive float is returned         :note: Typically the biggest positive float is returned
3476        """        """
3477        return self.getOperator().getUnlimitedTimeStepSize()        return self.getOperator().getUnlimitedTimeStepSize()
3478    
# Line 3459  class TransportPDE(LinearProblem): Line 3480  class TransportPDE(LinearProblem):
3480         """         """
3481         Returns a safe time step size to do the next time step.         Returns a safe time step size to do the next time step.
3482    
3483         @return: safe time step size         :return: safe time step size
3484         @rtype: C{float}         :rtype: ``float``
3485         @note: If not C{getSafeTimeStepSize()} < C{getUnlimitedTimeStepSize()}         :note: If not ``getSafeTimeStepSize()`` < ``getUnlimitedTimeStepSize()``
3486                any time step size can be used.                any time step size can be used.
3487         """         """
3488         return self.getOperator().getSafeTimeStepSize()         return self.getOperator().getSafeTimeStepSize()
# Line 3470  class TransportPDE(LinearProblem): Line 3491  class TransportPDE(LinearProblem):
3491         """         """
3492         Sets the weighting factor used to insert the constraints into the problem         Sets the weighting factor used to insert the constraints into the problem
3493    
3494         @param value: value for the weighting factor         :param value: value for the weighting factor
3495         @type value: large positive C{float}         :type value: large positive ``float``
3496         """         """
3497         if not value>0:         if not value>0:
3498           raise ValueError,"weighting factor needs to be positive."           raise ValueError,"weighting factor needs to be positive."
# Line 3481  class TransportPDE(LinearProblem): Line 3502  class TransportPDE(LinearProblem):
3502     def getConstraintWeightingFactor(self):     def getConstraintWeightingFactor(self):
3503         """         """
3504         returns the weighting factor used to insert the constraints into the problem         returns the weighting factor used to insert the constraints into the problem
3505         @return: value for the weighting factor         :return: value for the weighting factor
3506         @rtype: C{float}         :rtype: ``float``
3507         """         """
3508         return self.__constraint_factor         return self.__constraint_factor
3509     #====================================================================     #====================================================================
3510     def getSolution(self,dt):     def getSolution(self, dt=None, u0=None):
3511         """         """
3512         Returns the solution of the problem.         Returns the solution by marching forward by time step dt. if ''u0'' is present,
3513           ''u0'' is used as the initial value otherwise the solution from the last call is used.
3514    
3515         @return: the solution         :param dt: time step size. If ``None`` the last solution is returned.
3516         @rtype: L{Data<escript.Data>}         :type dt: positive ``float`` or ``None``
3517         """         :param u0: new initial solution or ``None``
3518         option_class=self.getSolverOptions()         :type u0: any object that can be interpolated to a `Data`
3519         if dt<=0:                  object on `Solution` or `ReducedSolution`
3520             raise ValueError,"step size needs to be positive."         :return: the solution
3521         self.setSolution(self.getOperator().solve(self.getRightHandSide(),dt,option_class))         :rtype: `Data`
3522         self.validSolution()         """
3523           if not dt == None:
3524          option_class=self.getSolverOptions()
3525          if dt<=0:
3526              raise ValueError,"step size needs to be positive."
3527          if u0 == None:
3528              u0=self.getCurrentSolution()
3529          else:
3530              u0=util.interpolate(u0,self.getFunctionSpaceForSolution())
3531              if self.getNumSolutions() == 1:
3532            if u0.getShape()!=():
3533                raise ValueError,"Illegal shape %s of initial solution."%(u0.getShape(),)
3534            else:
3535                if u0.getShape()!=(self.getNumSolutions(),):
3536                  raise ValueError,"Illegal shape %s of initial solution."%(u0.getShape(),)
3537          self.setSolution(self.getOperator().solve(u0, self.getRightHandSide(),dt,option_class))
3538          self.validSolution()
3539         return self.getCurrentSolution()         return self.getCurrentSolution()
3540    
3541       def setInitialSolution(self,u):
3542           """
3543           Sets the initial solution.
3544    
3545           :param u: initial solution
3546           :type u: any object that can be interpolated to a `Data`
3547                    object on `Solution` or `ReducedSolution`
3548           """
3549           u2=util.interpolate(u,self.getFunctionSpaceForSolution())
3550           if self.getNumSolutions() == 1:
3551              if u2.getShape()!=():
3552                  raise ValueError,"Illegal shape %s of initial solution."%(u2.getShape(),)
3553           else:
3554              if u2.getShape()!=(self.getNumSolutions(),):
3555                  raise ValueError,"Illegal shape %s of initial solution."%(u2.getShape(),)
3556           self.setSolution(u2,validate=False)
3557    
3558    
3559     def getSystem(self):     def getSystem(self):
3560         """         """
3561         Returns the operator and right hand side of the PDE.         Returns the operator and right hand side of the PDE.
3562    
3563         @return: the discrete version of the PDE         :return: the discrete version of the PDE
3564         @rtype: C{tuple} of L{Operator,<escript.Operator>} and         :rtype: ``tuple`` of `Operator` and
3565                 L{Data<escript.Data>}                 `Data`
3566    
3567         """         """
3568         if not self.isOperatorValid() or not self.isRightHandSideValid():         if not self.isOperatorValid() or not self.isRightHandSideValid():
# Line 3552  class TransportPDE(LinearProblem): Line 3608  class TransportPDE(LinearProblem):
3608    
3609     def setDebug(self, flag):     def setDebug(self, flag):
3610       """       """
3611       Switches debug output on if C{flag} is True,       Switches debug output on if ``flag`` is True,
3612       otherwise it is switched off.       otherwise it is switched off.
3613    
3614       @param flag: desired debug status       :param flag: desired debug status
3615       @type flag: C{bool}       :type flag: ``bool``
3616       """       """
3617       if flag:       if flag:
3618           self.setDebugOn()           self.setDebugOn()

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

  ViewVC Help
Powered by ViewVC 1.1.26