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

revision 2624 by jfenwick, Mon Jul 20 06:43:47 2009 UTC revision 2625 by jfenwick, Fri Aug 21 06:30:25 2009 UTC
# Line 28  Currently includes: Line 28  Currently includes:
28      - TimeIntegrationManager - to handle extrapolation in time      - TimeIntegrationManager - to handle extrapolation in time
29      - SaddlePointProblem - solver for Saddle point problems using the inexact uszawa scheme      - SaddlePointProblem - solver for Saddle point problems using the inexact uszawa scheme
30
31  @var __author__: name of author  :var __author__: name of author
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  __author__="Lutz Gross, l.gross@uq.edu.au"  __author__="Lutz Gross, l.gross@uq.edu.au"
# Line 64  class TimeIntegrationManager: Line 64  class TimeIntegrationManager:
64           tm.checkin(dt,v)           tm.checkin(dt,v)
65           t+=dt           t+=dt
66
67    @note: currently only p=1 is supported.    :note: currently only p=1 is supported.
68    """    """
69    def __init__(self,*inital_values,**kwargs):    def __init__(self,*inital_values,**kwargs):
70       """       """
71       Sets up the value manager where C{inital_values} are the initial values       Sets up the value manager where ``inital_values`` are the initial values
72       and p is the order used for extrapolation.       and p is the order used for extrapolation.
73       """       """
74       if kwargs.has_key("p"):       if kwargs.has_key("p"):
# Line 113  class TimeIntegrationManager: Line 113  class TimeIntegrationManager:
113
114    def extrapolate(self,dt):    def extrapolate(self,dt):
115        """        """
116        Extrapolates to C{dt} forward in time.        Extrapolates to ``dt`` forward in time.
117        """        """
118        if self.__order==0:        if self.__order==0:
119           out=self.__v_mem[0]           out=self.__v_mem[0]
# Line 139  class Projector: Line 139  class Projector:
139      """      """
140      Creates a continuous function space projector for a domain.      Creates a continuous function space projector for a domain.
141
142      @param domain: Domain of the projection.      :param domain: Domain of the projection.
143      @param reduce: Flag to reduce projection order      :param reduce: Flag to reduce projection order
144      @param fast: Flag to use a fast method based on matrix lumping      :param fast: Flag to use a fast method based on matrix lumping
145      """      """
146      self.__pde = linearPDEs.LinearPDE(domain)      self.__pde = linearPDEs.LinearPDE(domain)
147      if fast:      if fast:
# Line 154  class Projector: Line 154  class Projector:
154      """      """
155      Returns the solver options of the PDE solver.      Returns the solver options of the PDE solver.
156
157      @rtype: L{linearPDEs.SolverOptions}      :rtype: `linearPDEs.SolverOptions`
158      """      """
159
160    def __call__(self, input_data):    def __call__(self, input_data):
161      """      """
162      Projects C{input_data} onto a continuous function.      Projects ``input_data`` onto a continuous function.
163
164      @param input_data: the data to be projected      :param input_data: the data to be projected
165      """      """
166      out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())      out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
167      self.__pde.setValue(Y = escript.Data(), Y_reduced = escript.Data())      self.__pde.setValue(Y = escript.Data(), Y_reduced = escript.Data())
# Line 196  class NoPDE: Line 196  class NoPDE:
196       """       """
197       Solves the following problem for u:       Solves the following problem for u:
198
199       M{kronecker[i,j]*D[j]*u[j]=Y[i]}       *kronecker[i,j]*D[j]*u[j]=Y[i]*
200
201       with constraint       with constraint
202
203       M{u[j]=r[j]}  where M{q[j]>0}       *u[j]=r[j]*  where *q[j]>0*
204
205       where M{D}, M{Y}, M{r} and M{q} are given functions of rank 1.       where *D*, *Y*, *r* and *q* are given functions of rank 1.
206
207       In the case of scalars this takes the form       In the case of scalars this takes the form
208
209       M{D*u=Y}       *D*u=Y*
210
211       with constraint       with constraint
212
213       M{u=r} where M{q>0}       *u=r* where *q>0*
214
215       where M{D}, M{Y}, M{r} and M{q} are given scalar functions.       where *D*, *Y*, *r* and *q* are given scalar functions.
216
217       The constraint overwrites any other condition.       The constraint overwrites any other condition.
218
219       @note: This class is similar to the L{linearPDEs.LinearPDE} class with       :note: This class is similar to the `linearPDEs.LinearPDE` class with
220              A=B=C=X=0 but has the intention that all input parameters are given              A=B=C=X=0 but has the intention that all input parameters are given
221              in L{Solution} or L{ReducedSolution}.              in `Solution` or `ReducedSolution`.
222       """       """
223       # The whole thing is a bit strange and I blame Rob Woodcock (CSIRO) for       # The whole thing is a bit strange and I blame Rob Woodcock (CSIRO) for
224       # this.       # this.
# Line 226  class NoPDE: Line 226  class NoPDE:
226           """           """
227           Initializes the problem.           Initializes the problem.
228
229           @param domain: domain of the PDE           :param domain: domain of the PDE
230           @type domain: L{Domain}           :type domain: `Domain`
231           @param D: coefficient of the solution           :param D: coefficient of the solution
232           @type D: C{float}, C{int}, C{numpy.ndarray}, L{Data}           :type D: ``float``, ``int``, ``numpy.ndarray``, `Data`
233           @param Y: right hand side           :param Y: right hand side
234           @type Y: C{float}, C{int}, C{numpy.ndarray}, L{Data}           :type Y: ``float``, ``int``, ``numpy.ndarray``, `Data`
235           @param q: location of constraints           :param q: location of constraints
236           @type q: C{float}, C{int}, C{numpy.ndarray}, L{Data}           :type q: ``float``, ``int``, ``numpy.ndarray``, `Data`
237           @param r: value of solution at locations of constraints           :param r: value of solution at locations of constraints
238           @type r: C{float}, C{int}, C{numpy.ndarray}, L{Data}           :type r: ``float``, ``int``, ``numpy.ndarray``, `Data`
239           """           """
240           self.__domain=domain           self.__domain=domain
241           self.__D=D           self.__D=D
# Line 247  class NoPDE: Line 247  class NoPDE:
247
248       def setReducedOn(self):       def setReducedOn(self):
249           """           """
250           Sets the L{FunctionSpace} of the solution to L{ReducedSolution}.           Sets the `FunctionSpace` of the solution to `ReducedSolution`.
251           """           """
252           self.__function_space=escript.ReducedSolution(self.__domain)           self.__function_space=escript.ReducedSolution(self.__domain)
253           self.__u=None           self.__u=None
254
255       def setReducedOff(self):       def setReducedOff(self):
256           """           """
257           Sets the L{FunctionSpace} of the solution to L{Solution}.           Sets the `FunctionSpace` of the solution to `Solution`.
258           """           """
259           self.__function_space=escript.Solution(self.__domain)           self.__function_space=escript.Solution(self.__domain)
260           self.__u=None           self.__u=None
# Line 263  class NoPDE: Line 263  class NoPDE:
263           """           """
264           Assigns values to the parameters.           Assigns values to the parameters.
265
266           @param D: coefficient of the solution           :param D: coefficient of the solution
267           @type D: C{float}, C{int}, C{numpy.ndarray}, L{Data}           :type D: ``float``, ``int``, ``numpy.ndarray``, `Data`
268           @param Y: right hand side           :param Y: right hand side
269           @type Y: C{float}, C{int}, C{numpy.ndarray}, L{Data}           :type Y: ``float``, ``int``, ``numpy.ndarray``, `Data`
270           @param q: location of constraints           :param q: location of constraints
271           @type q: C{float}, C{int}, C{numpy.ndarray}, L{Data}           :type q: ``float``, ``int``, ``numpy.ndarray``, `Data`
272           @param r: value of solution at locations of constraints           :param r: value of solution at locations of constraints
273           @type r: C{float}, C{int}, C{numpy.ndarray}, L{Data}           :type r: ``float``, ``int``, ``numpy.ndarray``, `Data`
274           """           """
275           if not D==None:           if not D==None:
276              self.__D=D              self.__D=D
# Line 289  class NoPDE: Line 289  class NoPDE:
289           """           """
290           Returns the solution.           Returns the solution.
291
292           @return: the solution of the problem           :return: the solution of the problem
293           @rtype: L{Data} object in the L{FunctionSpace} L{Solution} or           :rtype: `Data` object in the `FunctionSpace` `Solution` or
294                   L{ReducedSolution}                   `ReducedSolution`
295           """           """
296           if self.__u==None:           if self.__u==None:
297              if self.__D==None:              if self.__D==None:
# Line 324  class Locator: Line 324  class Locator:
324         or FunctionSpace for the sample point which is closest to the given         or FunctionSpace for the sample point which is closest to the given
325         point x.         point x.
326
327         @param where: function space         :param where: function space
328         @type where: L{escript.FunctionSpace}         :type where: `escript.FunctionSpace`
329         @param x: location(s) of the Locator         :param x: location(s) of the Locator
330         @type x: C{numpy.ndarray} or C{list} of C{numpy.ndarray}         :type x: ``numpy.ndarray`` or ``list`` of ``numpy.ndarray``
331         """         """
332         if isinstance(where,escript.FunctionSpace):         if isinstance(where,escript.FunctionSpace):
333            self.__function_space=where            self.__function_space=where
# Line 401  class Locator: Line 401  class Locator:
401
402       def getValue(self,data):       def getValue(self,data):
403          """          """
404          Returns the value of C{data} at the Locator if C{data} is a L{Data}          Returns the value of ``data`` at the Locator if ``data`` is a `Data`
405          object otherwise the object is returned.          object otherwise the object is returned.
406          """          """
407          if isinstance(data,escript.Data):          if isinstance(data,escript.Data):
# Line 467  def PCG(r, Aprod, x, Msolve, bilinearfor Line 467  def PCG(r, Aprod, x, Msolve, bilinearfor
467     """     """
468     Solver for     Solver for
469
470     M{Ax=b}     *Ax=b*
471
472     with a symmetric and positive definite operator A (more details required!).     with a symmetric and positive definite operator A (more details required!).
473     It uses the conjugate gradient method with preconditioner M providing an     It uses the conjugate gradient method with preconditioner M providing an
# Line 475  def PCG(r, Aprod, x, Msolve, bilinearfor Line 475  def PCG(r, Aprod, x, Msolve, bilinearfor
475
476     The iteration is terminated if     The iteration is terminated if
477
478     M{|r| <= atol+rtol*|r0|}     *|r| <= atol+rtol*|r0|*
479
480     where M{r0} is the initial residual and M{|.|} is the energy norm. In fact     where *r0* is the initial residual and *|.|* is the energy norm. In fact
481
482     M{|r| = sqrt( bilinearform(Msolve(r),r))}     *|r| = sqrt( bilinearform(Msolve(r),r))*
483
484     For details on the preconditioned conjugate gradient method see the book:     For details on the preconditioned conjugate gradient method see the book:
485
# Line 487  def PCG(r, Aprod, x, Msolve, bilinearfor Line 487  def PCG(r, Aprod, x, Msolve, bilinearfor
487     T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,     T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
488     C. Romine, and H. van der Vorst}.     C. Romine, and H. van der Vorst}.
489
490     @param r: initial residual M{r=b-Ax}. C{r} is altered.     :param r: initial residual *r=b-Ax*. ``r`` is altered.
491     @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)     :type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
492     @param x: an initial guess for the solution     :param x: an initial guess for the solution
493     @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)     :type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
494     @param Aprod: returns the value Ax     :param Aprod: returns the value Ax
495     @type Aprod: function C{Aprod(x)} where C{x} is of the same object like     :type Aprod: function ``Aprod(x)`` where ``x`` is of the same object like
496                  argument C{x}. The returned object needs to be of the same type                  argument ``x``. The returned object needs to be of the same type
497                  like argument C{r}.                  like argument ``r``.
498     @param Msolve: solves Mx=r     :param Msolve: solves Mx=r
499     @type Msolve: function C{Msolve(r)} where C{r} is of the same type like     :type Msolve: function ``Msolve(r)`` where ``r`` is of the same type like
500                   argument C{r}. The returned object needs to be of the same                   argument ``r``. The returned object needs to be of the same
501                   type like argument C{x}.                   type like argument ``x``.
502     @param bilinearform: inner product C{<x,r>}     :param bilinearform: inner product ``<x,r>``
503     @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same     :type bilinearform: function ``bilinearform(x,r)`` where ``x`` is of the same
504                         type like argument C{x} and C{r} is. The returned value                         type like argument ``x`` and ``r`` is. The returned value
505                         is a C{float}.                         is a ``float``.
506     @param atol: absolute tolerance     :param atol: absolute tolerance
507     @type atol: non-negative C{float}     :type atol: non-negative ``float``
508     @param rtol: relative tolerance     :param rtol: relative tolerance
509     @type rtol: non-negative C{float}     :type rtol: non-negative ``float``
510     @param iter_max: maximum number of iteration steps     :param iter_max: maximum number of iteration steps
511     @type iter_max: C{int}     :type iter_max: ``int``
512     @return: the solution approximation and the corresponding residual     :return: the solution approximation and the corresponding residual
513     @rtype: C{tuple}     :rtype: ``tuple``
514     @warning: C{r} and C{x} are altered.     :warning: ``r`` and ``x`` are altered.
515     """     """
516     iter=0     iter=0
517     rhat=Msolve(r)     rhat=Msolve(r)
# Line 564  class Defect(object): Line 564  class Defect(object):
564          """          """
565          Returns the inner product of x0 and x1          Returns the inner product of x0 and x1
566
567          @param x0: value for x0          :param x0: value for x0
568          @param x1: value for x1          :param x1: value for x1
569          @return: the inner product of x0 and x1          :return: the inner product of x0 and x1
570          @rtype: C{float}          :rtype: ``float``
571          """          """
572          return 0          return 0
573
574      def norm(self,x):      def norm(self,x):
575          """          """
576          Returns the norm of argument C{x}.          Returns the norm of argument ``x``.
577
578          @param x: a value          :param x: a value
579          @return: norm of argument x          :return: norm of argument x
580          @rtype: C{float}          :rtype: ``float``
581          @note: by default C{sqrt(self.bilinearform(x,x)} is returned.          :note: by default ``sqrt(self.bilinearform(x,x)`` is returned.
582          """          """
583          s=self.bilinearform(x,x)          s=self.bilinearform(x,x)
584          if s<0: raise NegativeNorm,"negative norm."          if s<0: raise NegativeNorm,"negative norm."
# Line 586  class Defect(object): Line 586  class Defect(object):
586
587      def eval(self,x):      def eval(self,x):
588          """          """
589          Returns the value F of a given C{x}.          Returns the value F of a given ``x``.
590
591          @param x: value for which the defect C{F} is evaluated          :param x: value for which the defect ``F`` is evaluated
592          @return: value of the defect at C{x}          :return: value of the defect at ``x``
593          """          """
594          return 0          return 0
595
# Line 602  class Defect(object): Line 602  class Defect(object):
602          derivative of the defect. The increment is inc*norm(x)/norm(v)*v in the          derivative of the defect. The increment is inc*norm(x)/norm(v)*v in the
603          direction of v with x as a starting point.          direction of v with x as a starting point.
604
605          @param inc: relative increment length          :param inc: relative increment length
606          @type inc: positive C{float}          :type inc: positive ``float``
607          """          """
608          if inc<=0: raise ValueError,"positive increment required."          if inc<=0: raise ValueError,"positive increment required."
609          self.__inc=inc          self.__inc=inc
# Line 612  class Defect(object): Line 612  class Defect(object):
612          """          """
613          Returns the relative increment length used to approximate the          Returns the relative increment length used to approximate the
614          derivative of the defect.          derivative of the defect.
615          @return: value of the defect at C{x}          :return: value of the defect at ``x``
616          @rtype: positive C{float}          :rtype: positive ``float``
617          """          """
618          return self.__inc          return self.__inc
619
620      def derivative(self, F0, x0, v, v_is_normalised=True):      def derivative(self, F0, x0, v, v_is_normalised=True):
621          """          """
622          Returns the directional derivative at C{x0} in the direction of C{v}.          Returns the directional derivative at ``x0`` in the direction of ``v``.
623
624          @param F0: value of this defect at x0          :param F0: value of this defect at x0
625          @param x0: value at which derivative is calculated          :param x0: value at which derivative is calculated
626          @param v: direction          :param v: direction
627          @param v_is_normalised: True to indicate that C{v} is nomalized          :param v_is_normalised: True to indicate that ``v`` is nomalized
628                                  (self.norm(v)=0)                                  (self.norm(v)=0)
629          @return: derivative of this defect at x0 in the direction of C{v}          :return: derivative of this defect at x0 in the direction of ``v``
630          @note: by default numerical evaluation (self.eval(x0+eps*v)-F0)/eps is          :note: by default numerical evaluation (self.eval(x0+eps*v)-F0)/eps is
631                 used but this method maybe overwritten to use exact evaluation.                 used but this method maybe overwritten to use exact evaluation.
632          """          """
633          normx=self.norm(x0)          normx=self.norm(x0)
# Line 647  class Defect(object): Line 647  class Defect(object):
647  ######################################  ######################################
648  def NewtonGMRES(defect, x, iter_max=100, sub_iter_max=20, atol=0,rtol=1.e-4, sub_tol_max=0.5, gamma=0.9, verbose=False):  def NewtonGMRES(defect, x, iter_max=100, sub_iter_max=20, atol=0,rtol=1.e-4, sub_tol_max=0.5, gamma=0.9, verbose=False):
649     """     """
650     Solves a non-linear problem M{F(x)=0} for unknown M{x} using the stopping     Solves a non-linear problem *F(x)=0* for unknown *x* using the stopping
651     criterion:     criterion:
652
653     M{norm(F(x) <= atol + rtol * norm(F(x0)}     *norm(F(x) <= atol + rtol * norm(F(x0)*
654
655     where M{x0} is the initial guess.     where *x0* is the initial guess.
656
657     @param defect: object defining the function M{F}. C{defect.norm} defines the     :param defect: object defining the function *F*. ``defect.norm`` defines the
658                    M{norm} used in the stopping criterion.                    *norm* used in the stopping criterion.
659     @type defect: L{Defect}     :type defect: `Defect`
660     @param x: initial guess for the solution, C{x} is altered.     :param x: initial guess for the solution, ``x`` is altered.
661     @type x: any object type allowing basic operations such as     :type x: any object type allowing basic operations such as
662              C{numpy.ndarray}, L{Data}              ``numpy.ndarray``, `Data`
663     @param iter_max: maximum number of iteration steps     :param iter_max: maximum number of iteration steps
664     @type iter_max: positive C{int}     :type iter_max: positive ``int``
665     @param sub_iter_max: maximum number of inner iteration steps     :param sub_iter_max: maximum number of inner iteration steps
666     @type sub_iter_max: positive C{int}     :type sub_iter_max: positive ``int``
667     @param atol: absolute tolerance for the solution     :param atol: absolute tolerance for the solution
668     @type atol: positive C{float}     :type atol: positive ``float``
669     @param rtol: relative tolerance for the solution     :param rtol: relative tolerance for the solution
670     @type rtol: positive C{float}     :type rtol: positive ``float``
671     @param gamma: tolerance safety factor for inner iteration     :param gamma: tolerance safety factor for inner iteration
672     @type gamma: positive C{float}, less than 1     :type gamma: positive ``float``, less than 1
673     @param sub_tol_max: upper bound for inner tolerance     :param sub_tol_max: upper bound for inner tolerance
674     @type sub_tol_max: positive C{float}, less than 1     :type sub_tol_max: positive ``float``, less than 1
675     @return: an approximation of the solution with the desired accuracy     :return: an approximation of the solution with the desired accuracy
676     @rtype: same type as the initial guess     :rtype: same type as the initial guess
677     """     """
678     lmaxit=iter_max     lmaxit=iter_max
679     if atol<0: raise ValueError,"atol needs to be non-negative."     if atol<0: raise ValueError,"atol needs to be non-negative."
# Line 731  def NewtonGMRES(defect, x, iter_max=100, Line 731  def NewtonGMRES(defect, x, iter_max=100,
731  def __givapp(c,s,vin):  def __givapp(c,s,vin):
732      """      """
733      Applies a sequence of Givens rotations (c,s) recursively to the vector      Applies a sequence of Givens rotations (c,s) recursively to the vector
734      C{vin}      ``vin``
735
736      @warning: C{vin} is altered.      :warning: ``vin`` is altered.
737      """      """
738      vrot=vin      vrot=vin
739      if isinstance(c,float):      if isinstance(c,float):
# Line 838  def GMRES(r, Aprod, x, bilinearform, ato Line 838  def GMRES(r, Aprod, x, bilinearform, ato
838     """     """
839     Solver for     Solver for
840
841     M{Ax=b}     *Ax=b*
842
843     with a general operator A (more details required!).     with a general operator A (more details required!).
844     It uses the generalized minimum residual method (GMRES).     It uses the generalized minimum residual method (GMRES).
845
846     The iteration is terminated if     The iteration is terminated if
847
848     M{|r| <= atol+rtol*|r0|}     *|r| <= atol+rtol*|r0|*
849
850     where M{r0} is the initial residual and M{|.|} is the energy norm. In fact     where *r0* is the initial residual and *|.|* is the energy norm. In fact
851
852     M{|r| = sqrt( bilinearform(r,r))}     *|r| = sqrt( bilinearform(r,r))*
853
854     @param r: initial residual M{r=b-Ax}. C{r} is altered.     :param r: initial residual *r=b-Ax*. ``r`` is altered.
855     @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)     :type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
856     @param x: an initial guess for the solution     :param x: an initial guess for the solution
857     @type x: same like C{r}     :type x: same like ``r``
858     @param Aprod: returns the value Ax     :param Aprod: returns the value Ax
859     @type Aprod: function C{Aprod(x)} where C{x} is of the same object like     :type Aprod: function ``Aprod(x)`` where ``x`` is of the same object like
860                  argument C{x}. The returned object needs to be of the same                  argument ``x``. The returned object needs to be of the same
861                  type like argument C{r}.                  type like argument ``r``.
862     @param bilinearform: inner product C{<x,r>}     :param bilinearform: inner product ``<x,r>``
863     @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same     :type bilinearform: function ``bilinearform(x,r)`` where ``x`` is of the same
864                         type like argument C{x} and C{r}. The returned value is                         type like argument ``x`` and ``r``. The returned value is
865                         a C{float}.                         a ``float``.
866     @param atol: absolute tolerance     :param atol: absolute tolerance
867     @type atol: non-negative C{float}     :type atol: non-negative ``float``
868     @param rtol: relative tolerance     :param rtol: relative tolerance
869     @type rtol: non-negative C{float}     :type rtol: non-negative ``float``
870     @param iter_max: maximum number of iteration steps     :param iter_max: maximum number of iteration steps
871     @type iter_max: C{int}     :type iter_max: ``int``
872     @param iter_restart: in order to save memory the orthogonalization process     :param iter_restart: in order to save memory the orthogonalization process
873                          is terminated after C{iter_restart} steps and the                          is terminated after ``iter_restart`` steps and the
874                          iteration is restarted.                          iteration is restarted.
875     @type iter_restart: C{int}     :type iter_restart: ``int``
876     @return: the solution approximation and the corresponding residual     :return: the solution approximation and the corresponding residual
877     @rtype: C{tuple}     :rtype: ``tuple``
878     @warning: C{r} and C{x} are altered.     :warning: ``r`` and ``x`` are altered.
879     """     """
880     m=iter_restart     m=iter_restart
881     restarted=False     restarted=False
# Line 1007  def MINRES(r, Aprod, x, Msolve, bilinear Line 1007  def MINRES(r, Aprod, x, Msolve, bilinear
1007      """      """
1008      Solver for      Solver for
1009
1010      M{Ax=b}      *Ax=b*
1011
1012      with a symmetric and positive definite operator A (more details required!).      with a symmetric and positive definite operator A (more details required!).
1013      It uses the minimum residual method (MINRES) with preconditioner M      It uses the minimum residual method (MINRES) with preconditioner M
# Line 1015  def MINRES(r, Aprod, x, Msolve, bilinear Line 1015  def MINRES(r, Aprod, x, Msolve, bilinear
1015
1016      The iteration is terminated if      The iteration is terminated if
1017
1018      M{|r| <= atol+rtol*|r0|}      *|r| <= atol+rtol*|r0|*
1019
1020      where M{r0} is the initial residual and M{|.|} is the energy norm. In fact      where *r0* is the initial residual and *|.|* is the energy norm. In fact
1021
1022      M{|r| = sqrt( bilinearform(Msolve(r),r))}      *|r| = sqrt( bilinearform(Msolve(r),r))*
1023
1024      For details on the preconditioned conjugate gradient method see the book:      For details on the preconditioned conjugate gradient method see the book:
1025
# Line 1027  def MINRES(r, Aprod, x, Msolve, bilinear Line 1027  def MINRES(r, Aprod, x, Msolve, bilinear
1027      T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,      T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
1028      C. Romine, and H. van der Vorst}.      C. Romine, and H. van der Vorst}.
1029
1030      @param r: initial residual M{r=b-Ax}. C{r} is altered.      :param r: initial residual *r=b-Ax*. ``r`` is altered.
1031      @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)      :type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1032      @param x: an initial guess for the solution      :param x: an initial guess for the solution
1033      @type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)      :type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1034      @param Aprod: returns the value Ax      :param Aprod: returns the value Ax
1035      @type Aprod: function C{Aprod(x)} where C{x} is of the same object like      :type Aprod: function ``Aprod(x)`` where ``x`` is of the same object like
1036                   argument C{x}. The returned object needs to be of the same                   argument ``x``. The returned object needs to be of the same
1037                   type like argument C{r}.                   type like argument ``r``.
1038      @param Msolve: solves Mx=r      :param Msolve: solves Mx=r
1039      @type Msolve: function C{Msolve(r)} where C{r} is of the same type like      :type Msolve: function ``Msolve(r)`` where ``r`` is of the same type like
1040                    argument C{r}. The returned object needs to be of the same                    argument ``r``. The returned object needs to be of the same
1041                    type like argument C{x}.                    type like argument ``x``.
1042      @param bilinearform: inner product C{<x,r>}      :param bilinearform: inner product ``<x,r>``
1043      @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same      :type bilinearform: function ``bilinearform(x,r)`` where ``x`` is of the same
1044                          type like argument C{x} and C{r} is. The returned value                          type like argument ``x`` and ``r`` is. The returned value
1045                          is a C{float}.                          is a ``float``.
1046      @param atol: absolute tolerance      :param atol: absolute tolerance
1047      @type atol: non-negative C{float}      :type atol: non-negative ``float``
1048      @param rtol: relative tolerance      :param rtol: relative tolerance
1049      @type rtol: non-negative C{float}      :type rtol: non-negative ``float``
1050      @param iter_max: maximum number of iteration steps      :param iter_max: maximum number of iteration steps
1051      @type iter_max: C{int}      :type iter_max: ``int``
1052      @return: the solution approximation and the corresponding residual      :return: the solution approximation and the corresponding residual
1053      @rtype: C{tuple}      :rtype: ``tuple``
1054      @warning: C{r} and C{x} are altered.      :warning: ``r`` and ``x`` are altered.
1055      """      """
1056      #------------------------------------------------------------------      #------------------------------------------------------------------
1057      # Set up y and v for the first Lanczos vector v1.      # Set up y and v for the first Lanczos vector v1.
# Line 1185  def TFQMR(r, Aprod, x, bilinearform, ato Line 1185  def TFQMR(r, Aprod, x, bilinearform, ato
1185    """    """
1186    Solver for    Solver for
1187
1188    M{Ax=b}    *Ax=b*
1189
1190    with a general operator A (more details required!).    with a general operator A (more details required!).
1191    It uses the Transpose-Free Quasi-Minimal Residual method (TFQMR).    It uses the Transpose-Free Quasi-Minimal Residual method (TFQMR).
1192
1193    The iteration is terminated if    The iteration is terminated if
1194
1195    M{|r| <= atol+rtol*|r0|}    *|r| <= atol+rtol*|r0|*
1196
1197    where M{r0} is the initial residual and M{|.|} is the energy norm. In fact    where *r0* is the initial residual and *|.|* is the energy norm. In fact
1198
1199    M{|r| = sqrt( bilinearform(r,r))}    *|r| = sqrt( bilinearform(r,r))*
1200
1201    @param r: initial residual M{r=b-Ax}. C{r} is altered.    :param r: initial residual *r=b-Ax*. ``r`` is altered.
1202    @type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)    :type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1203    @param x: an initial guess for the solution    :param x: an initial guess for the solution
1204    @type x: same like C{r}    :type x: same like ``r``
1205    @param Aprod: returns the value Ax    :param Aprod: returns the value Ax
1206    @type Aprod: function C{Aprod(x)} where C{x} is of the same object like    :type Aprod: function ``Aprod(x)`` where ``x`` is of the same object like
1207                 argument C{x}. The returned object needs to be of the same type                 argument ``x``. The returned object needs to be of the same type
1208                 like argument C{r}.                 like argument ``r``.
1209    @param bilinearform: inner product C{<x,r>}    :param bilinearform: inner product ``<x,r>``
1210    @type bilinearform: function C{bilinearform(x,r)} where C{x} is of the same    :type bilinearform: function ``bilinearform(x,r)`` where ``x`` is of the same
1211                        type like argument C{x} and C{r}. The returned value is                        type like argument ``x`` and ``r``. The returned value is
1212                        a C{float}.                        a ``float``.
1213    @param atol: absolute tolerance    :param atol: absolute tolerance
1214    @type atol: non-negative C{float}    :type atol: non-negative ``float``
1215    @param rtol: relative tolerance    :param rtol: relative tolerance
1216    @type rtol: non-negative C{float}    :type rtol: non-negative ``float``
1217    @param iter_max: maximum number of iteration steps    :param iter_max: maximum number of iteration steps
1218    @type iter_max: C{int}    :type iter_max: ``int``
1219    @rtype: C{tuple}    :rtype: ``tuple``
1220    @warning: C{r} and C{x} are altered.    :warning: ``r`` and ``x`` are altered.
1221    """    """
1222    u1=0    u1=0
1223    u2=0    u2=0
# Line 1287  def TFQMR(r, Aprod, x, bilinearform, ato Line 1287  def TFQMR(r, Aprod, x, bilinearform, ato
1287
1288  class ArithmeticTuple(object):  class ArithmeticTuple(object):
1289     """     """
1290     Tuple supporting inplace update x+=y and scaling x=a*y where C{x,y} is an     Tuple supporting inplace update x+=y and scaling x=a*y where ``x,y`` is an
1291     ArithmeticTuple and C{a} is a float.     ArithmeticTuple and ``a`` is a float.
1292
1293     Example of usage::     Example of usage::
1294
# Line 1302  class ArithmeticTuple(object): Line 1302  class ArithmeticTuple(object):
1302     """     """
1303     def __init__(self,*args):     def __init__(self,*args):
1304         """         """
1305         Initializes object with elements C{args}.         Initializes object with elements ``args``.
1306
1307         @param args: tuple of objects that support inplace add (x+=y) and         :param args: tuple of objects that support inplace add (x+=y) and
1308                      scaling (x=a*y)                      scaling (x=a*y)
1309         """         """
1310         self.__items=list(args)         self.__items=list(args)
# Line 1313  class ArithmeticTuple(object): Line 1313  class ArithmeticTuple(object):
1313         """         """
1314         Returns the number of items.         Returns the number of items.
1315
1316         @return: number of items         :return: number of items
1317         @rtype: C{int}         :rtype: ``int``
1318         """         """
1319         return len(self.__items)         return len(self.__items)
1320
# Line 1322  class ArithmeticTuple(object): Line 1322  class ArithmeticTuple(object):
1322         """         """
1323         Returns item at specified position.         Returns item at specified position.
1324
1325         @param index: index of item to be returned         :param index: index of item to be returned
1326         @type index: C{int}         :type index: ``int``
1327         @return: item with index C{index}         :return: item with index ``index``
1328         """         """
1329         return self.__items.__getitem__(index)         return self.__items.__getitem__(index)
1330
1331     def __mul__(self,other):     def __mul__(self,other):
1332         """         """
1333         Scales by C{other} from the right.         Scales by ``other`` from the right.
1334
1335         @param other: scaling factor         :param other: scaling factor
1336         @type other: C{float}         :type other: ``float``
1337         @return: itemwise self*other         :return: itemwise self*other
1338         @rtype: L{ArithmeticTuple}         :rtype: `ArithmeticTuple`
1339         """         """
1340         out=[]         out=[]
1341         try:         try:
# Line 1349  class ArithmeticTuple(object): Line 1349  class ArithmeticTuple(object):
1349
1350     def __rmul__(self,other):     def __rmul__(self,other):
1351         """         """
1352         Scales by C{other} from the left.         Scales by ``other`` from the left.
1353
1354         @param other: scaling factor         :param other: scaling factor
1355         @type other: C{float}         :type other: ``float``
1356         @return: itemwise other*self         :return: itemwise other*self
1357         @rtype: L{ArithmeticTuple}         :rtype: `ArithmeticTuple`
1358         """         """
1359         out=[]         out=[]
1360         try:         try:
# Line 1368  class ArithmeticTuple(object): Line 1368  class ArithmeticTuple(object):
1368
1369     def __div__(self,other):     def __div__(self,other):
1370         """         """
1371         Scales by (1/C{other}) from the right.         Scales by (1/``other``) from the right.
1372
1373         @param other: scaling factor         :param other: scaling factor
1374         @type other: C{float}         :type other: ``float``
1375         @return: itemwise self/other         :return: itemwise self/other
1376         @rtype: L{ArithmeticTuple}         :rtype: `ArithmeticTuple`
1377         """         """
1378         return self*(1/other)         return self*(1/other)
1379
1380     def __rdiv__(self,other):     def __rdiv__(self,other):
1381         """         """
1382         Scales by (1/C{other}) from the left.         Scales by (1/``other``) from the left.
1383
1384         @param other: scaling factor         :param other: scaling factor
1385         @type other: C{float}         :type other: ``float``
1386         @return: itemwise other/self         :return: itemwise other/self
1387         @rtype: L{ArithmeticTuple}         :rtype: `ArithmeticTuple`
1388         """         """
1389         out=[]         out=[]
1390         try:         try:
# Line 1398  class ArithmeticTuple(object): Line 1398  class ArithmeticTuple(object):
1398
1400         """         """
1401         Inplace addition of C{other} to self.         Inplace addition of ``other`` to self.
1402
1403         @param other: increment         :param other: increment
1404         @type other: C{ArithmeticTuple}         :type other: ``ArithmeticTuple``
1405         """         """
1406         if len(self) != len(other):         if len(self) != len(other):
1407             raise ValueError,"tuple lengths must match."             raise ValueError,"tuple lengths must match."
# Line 1411  class ArithmeticTuple(object): Line 1411  class ArithmeticTuple(object):
1411
1413         """         """
1414         Adds C{other} to self.         Adds ``other`` to self.
1415
1416         @param other: increment         :param other: increment
1417         @type other: C{ArithmeticTuple}         :type other: ``ArithmeticTuple``
1418         """         """
1419         out=[]         out=[]
1420         try:         try:
# Line 1428  class ArithmeticTuple(object): Line 1428  class ArithmeticTuple(object):
1428
1429     def __sub__(self,other):     def __sub__(self,other):
1430         """         """
1431         Subtracts C{other} from self.         Subtracts ``other`` from self.
1432
1433         @param other: decrement         :param other: decrement
1434         @type other: C{ArithmeticTuple}         :type other: ``ArithmeticTuple``
1435         """         """
1436         out=[]         out=[]
1437         try:         try:
# Line 1445  class ArithmeticTuple(object): Line 1445  class ArithmeticTuple(object):
1445
1446     def __isub__(self,other):     def __isub__(self,other):
1447         """         """
1448         Inplace subtraction of C{other} from self.         Inplace subtraction of ``other`` from self.
1449
1450         @param other: decrement         :param other: decrement
1451         @type other: C{ArithmeticTuple}         :type other: ``ArithmeticTuple``
1452         """         """
1453         if len(self) != len(other):         if len(self) != len(other):
1454             raise ValueError,"tuple length must match."             raise ValueError,"tuple length must match."
# Line 1471  class HomogeneousSaddlePointProblem(obje Line 1471  class HomogeneousSaddlePointProblem(obje
1471        This class provides a framework for solving linear homogeneous saddle        This class provides a framework for solving linear homogeneous saddle
1472        point problems of the form::        point problems of the form::
1473
1474            M{Av+B^*p=f}            *Av+B^*p=f*
1475            M{Bv     =0}            *Bv     =0*
1476
1477        for the unknowns M{v} and M{p} and given operators M{A} and M{B} and        for the unknowns *v* and *p* and given operators *A* and *B* and
1478        given right hand side M{f}. M{B^*} is the adjoint operator of M{B}.        given right hand side *f*. *B^** is the adjoint operator of *B*.
1479        """        """
1480        def __init__(self, adaptSubTolerance=True, **kwargs):        def __init__(self, adaptSubTolerance=True, **kwargs):
1481      """      """
1482      initializes the saddle point problem      initializes the saddle point problem
1483
1484      @param adaptSubTolerance: If True the tolerance for subproblem is set automatically.      :param adaptSubTolerance: If True the tolerance for subproblem is set automatically.
1485      @type adaptSubTolerance: C{bool}      :type adaptSubTolerance: ``bool``
1486      """      """
1487          self.setTolerance()          self.setTolerance()
1488          self.setAbsoluteTolerance()          self.setAbsoluteTolerance()
# Line 1498  class HomogeneousSaddlePointProblem(obje Line 1498  class HomogeneousSaddlePointProblem(obje
1498           """           """
1499           Returns inner product of element p and Bv (overwrite).           Returns inner product of element p and Bv (overwrite).
1500
1501           @param p: a pressure increment           :param p: a pressure increment
1502           @param v: a residual           :param Bv: a residual
1503           @return: inner product of element p and Bv           :return: inner product of element p and Bv
1504           @rtype: C{float}           :rtype: ``float``
1505           @note: used if PCG is applied.           :note: used if PCG is applied.
1506           """           """
1507           raise NotImplementedError,"no inner product for p and Bv implemented."           raise NotImplementedError,"no inner product for p and Bv implemented."
1508
# Line 1510  class HomogeneousSaddlePointProblem(obje Line 1510  class HomogeneousSaddlePointProblem(obje
1510           """           """
1511           Returns inner product of p0 and p1 (overwrite).           Returns inner product of p0 and p1 (overwrite).
1512
1513           @param p0: a pressure           :param p0: a pressure
1514           @param p1: a pressure           :param p1: a pressure
1515           @return: inner product of p0 and p1           :return: inner product of p0 and p1
1516           @rtype: C{float}           :rtype: ``float``
1517           """           """
1518           raise NotImplementedError,"no inner product for p implemented."           raise NotImplementedError,"no inner product for p implemented."
1519
# Line 1521  class HomogeneousSaddlePointProblem(obje Line 1521  class HomogeneousSaddlePointProblem(obje
1521           """           """
1522           Returns the norm of v (overwrite).           Returns the norm of v (overwrite).
1523
1524           @param v: a velovity           :param v: a velovity
1525           @return: norm of v           :return: norm of v
1526           @rtype: non-negative C{float}           :rtype: non-negative ``float``
1527           """           """
1528           raise NotImplementedError,"no norm of v implemented."           raise NotImplementedError,"no norm of v implemented."
1529        def getV(self, p, v0):        def getV(self, p, v0):
1530           """           """
1531           return the value for v for a given p (overwrite)           return the value for v for a given p (overwrite)
1532
1533           @param p: a pressure           :param p: a pressure
1534           @param v0: a initial guess for the value v to return.           :param v0: a initial guess for the value v to return.
1535           @return: v given as M{v= A^{-1} (f-B^*p)}           :return: v given as *v= A^{-1} (f-B^*p)*
1536           """           """
1537           raise NotImplementedError,"no v calculation implemented."           raise NotImplementedError,"no v calculation implemented."
1538
# Line 1541  class HomogeneousSaddlePointProblem(obje Line 1541  class HomogeneousSaddlePointProblem(obje
1541          """          """
1542          Returns Bv (overwrite).          Returns Bv (overwrite).
1543
1544          @rtype: equal to the type of p          :rtype: equal to the type of p
1545          @note: boundary conditions on p should be zero!          :note: boundary conditions on p should be zero!
1546          """          """
1547          raise NotImplementedError, "no operator B implemented."          raise NotImplementedError, "no operator B implemented."
1548
# Line 1550  class HomogeneousSaddlePointProblem(obje Line 1550  class HomogeneousSaddlePointProblem(obje
1550          """          """
1551          Returns the norm of Bv (overwrite).          Returns the norm of Bv (overwrite).
1552
1553          @rtype: equal to the type of p          :rtype: equal to the type of p
1554          @note: boundary conditions on p should be zero!          :note: boundary conditions on p should be zero!
1555          """          """
1556          raise NotImplementedError, "no norm of Bv implemented."          raise NotImplementedError, "no norm of Bv implemented."
1557
1558        def solve_AinvBt(self,p):        def solve_AinvBt(self,p):
1559           """           """
1560           Solves M{Av=B^*p} with accuracy L{self.getSubProblemTolerance()}           Solves *Av=B^*p* with accuracy `self.getSubProblemTolerance()`
1561           (overwrite).           (overwrite).
1562
1563           @param p: a pressure increment           :param p: a pressure increment
1564           @return: the solution of M{Av=B^*p}           :return: the solution of *Av=B^*p*
1565           @note: boundary conditions on v should be zero!           :note: boundary conditions on v should be zero!
1566           """           """
1567           raise NotImplementedError,"no operator A implemented."           raise NotImplementedError,"no operator A implemented."
1568
1569        def solve_prec(self,Bv):        def solve_prec(self,Bv):
1570           """           """
1571           Provides a preconditioner for M{BA^{-1}B^*} applied to Bv with accuracy           Provides a preconditioner for *(BA^{-1}B^ * )* applied to Bv with accuracy
1572           L{self.getSubProblemTolerance()} (overwrite).           `self.getSubProblemTolerance()` (overwrite).
1573
1574           @rtype: equal to the type of p           :rtype: equal to the type of p
1575           @note: boundary conditions on p should be zero!           :note: boundary conditions on p should be zero!
1576           """           """
1577           raise NotImplementedError,"no preconditioner for Schur complement implemented."           raise NotImplementedError,"no preconditioner for Schur complement implemented."
1578        def setSubProblemTolerance(self):        def setSubProblemTolerance(self):
1579           """           """
1580       Updates the tolerance for subproblems       Updates the tolerance for subproblems
1581       @note: method is typically the method is overwritten.       :note: method is typically the method is overwritten.
1582           """           """
1583           pass           pass
1584        #=============================================================        #=============================================================
# Line 1600  class HomogeneousSaddlePointProblem(obje Line 1600  class HomogeneousSaddlePointProblem(obje
1600        #=============================================================        #=============================================================
1601        def norm_p(self,p):        def norm_p(self,p):
1602            """            """
1603            calculates the norm of C{p}            calculates the norm of ``p``
1604
1605            @param p: a pressure            :param p: a pressure
1606            @return: the norm of C{p} using the inner product for pressure            :return: the norm of ``p`` using the inner product for pressure
1607            @rtype: C{float}            :rtype: ``float``
1608            """            """
1609            f=self.inner_p(p,p)            f=self.inner_p(p,p)
1610            if f<0: raise ValueError,"negative pressure norm."            if f<0: raise ValueError,"negative pressure norm."
# Line 1619  class HomogeneousSaddlePointProblem(obje Line 1619  class HomogeneousSaddlePointProblem(obje
1619           """           """
1620           Solves the saddle point problem using initial guesses v and p.           Solves the saddle point problem using initial guesses v and p.
1621
1622           @param v: initial guess for velocity           :param v: initial guess for velocity
1623           @param p: initial guess for pressure           :param p: initial guess for pressure
1624           @type v: L{Data}           :type v: `Data`
1625           @type p: L{Data}           :type p: `Data`
1626           @param usePCG: indicates the usage of the PCG rather than GMRES scheme.           :param usePCG: indicates the usage of the PCG rather than GMRES scheme.
1627           @param max_iter: maximum number of iteration steps per correction           :param max_iter: maximum number of iteration steps per correction
1628                            attempt                            attempt
1629           @param verbose: if True, shows information on the progress of the           :param verbose: if True, shows information on the progress of the
1630                           saddlepoint problem solver.                           saddlepoint problem solver.
1631           @param iter_restart: restart the iteration after C{iter_restart} steps           :param iter_restart: restart the iteration after ``iter_restart`` steps
1632                                (only used if useUzaw=False)                                (only used if useUzaw=False)
1633           @type usePCG: C{bool}           :type usePCG: ``bool``
1634           @type max_iter: C{int}           :type max_iter: ``int``
1635           @type verbose: C{bool}           :type verbose: ``bool``
1636           @type iter_restart: C{int}           :type iter_restart: ``int``
1637           @rtype: C{tuple} of L{Data} objects           :rtype: ``tuple`` of `Data` objects
1638           """           """
1639           self.verbose=verbose           self.verbose=verbose
1640           rtol=self.getTolerance()           rtol=self.getTolerance()
# Line 1678  class HomogeneousSaddlePointProblem(obje Line 1678  class HomogeneousSaddlePointProblem(obje
1678           """           """
1679           Sets the relative tolerance for (v,p).           Sets the relative tolerance for (v,p).
1680
1681           @param tolerance: tolerance to be used           :param tolerance: tolerance to be used
1682           @type tolerance: non-negative C{float}           :type tolerance: non-negative ``float``
1683           """           """
1684           if tolerance<0:           if tolerance<0:
1685               raise ValueError,"tolerance must be positive."               raise ValueError,"tolerance must be positive."
# Line 1689  class HomogeneousSaddlePointProblem(obje Line 1689  class HomogeneousSaddlePointProblem(obje
1689           """           """
1690           Returns the relative tolerance.           Returns the relative tolerance.
1691
1692           @return: relative tolerance           :return: relative tolerance
1693           @rtype: C{float}           :rtype: ``float``
1694           """           """
1695           return self.__rtol           return self.__rtol
1696
# Line 1698  class HomogeneousSaddlePointProblem(obje Line 1698  class HomogeneousSaddlePointProblem(obje
1698           """           """
1699           Sets the absolute tolerance.           Sets the absolute tolerance.
1700
1701           @param tolerance: tolerance to be used           :param tolerance: tolerance to be used
1702           @type tolerance: non-negative C{float}           :type tolerance: non-negative ``float``
1703           """           """
1704           if tolerance<0:           if tolerance<0:
1705               raise ValueError,"tolerance must be non-negative."               raise ValueError,"tolerance must be non-negative."
# Line 1709  class HomogeneousSaddlePointProblem(obje Line 1709  class HomogeneousSaddlePointProblem(obje
1709           """           """
1710           Returns the absolute tolerance.           Returns the absolute tolerance.
1711
1712           @return: absolute tolerance           :return: absolute tolerance
1713           @rtype: C{float}           :rtype: ``float``
1714           """           """
1715           return self.__atol           return self.__atol
1716
1717        def getSubProblemTolerance(self):        def getSubProblemTolerance(self):
1718           """           """
1719           Sets the relative tolerance to solve the subproblem(s).           Sets the relative tolerance to solve the subproblem(s).

@param rtol: relative tolerence
@type rtol: positive C{float}
1720           """           """
1721           return max(200.*util.EPSILON,self.getTolerance()**2)           return max(200.*util.EPSILON,self.getTolerance()**2)
1722
# Line 1730  def MaskFromBoundaryTag(domain,*tags): Line 1727  def MaskFromBoundaryTag(domain,*tags):
1727
1728     Usage: m=MaskFromBoundaryTag(domain, "left", "right")     Usage: m=MaskFromBoundaryTag(domain, "left", "right")
1729
1730     @param domain: domain to be used     :param domain: domain to be used
1731     @type domain: L{escript.Domain}     :type domain: `escript.Domain`
1732     @param tags: boundary tags     :param tags: boundary tags
1733     @type tags: C{str}     :type tags: ``str``
1734     @return: a mask which marks samples that are touching the boundary tagged     :return: a mask which marks samples that are touching the boundary tagged
1735              by any of the given tags              by any of the given tags
1736     @rtype: L{escript.Data} of rank 0     :rtype: `escript.Data` of rank 0
1737     """     """
1738     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1739     d=escript.Scalar(0.,escript.FunctionOnBoundary(domain))     d=escript.Scalar(0.,escript.FunctionOnBoundary(domain))
# Line 1751  def MaskFromTag(domain,*tags): Line 1748  def MaskFromTag(domain,*tags):
1748
1749     Usage: m=MaskFromTag(domain, "ham")     Usage: m=MaskFromTag(domain, "ham")
1750
1751     @param domain: domain to be used     :param domain: domain to be used
1752     @type domain: L{escript.Domain}     :type domain: `escript.Domain`
1753     @param tags: boundary tags     :param tags: boundary tags
1754     @type tags: C{str}     :type tags: ``str``
1755     @return: a mask which marks samples that are touching the boundary tagged     :return: a mask which marks samples that are touching the boundary tagged
1756              by any of the given tags              by any of the given tags
1757     @rtype: L{escript.Data} of rank 0     :rtype: `escript.Data` of rank 0
1758     """     """
1759     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)     pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1760     d=escript.Scalar(0.,escript.Function(domain))     d=escript.Scalar(0.,escript.Function(domain))

Legend:
 Removed from v.2624 changed lines Added in v.2625