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

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

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

revision 351 by gross, Tue Dec 13 09:12:15 2005 UTC revision 3928 by jfenwick, Mon Jul 9 23:23:43 2012 UTC
# Line 1  Line 1 
1  # $Id$  
2    ########################################################
3    #
4    # Copyright (c) 2003-2012 by University of Queensland
5    # Earth Systems Science Computational Center (ESSCC)
6    # http://www.uq.edu.au/esscc
7    #
8    # Primary Business: Queensland, Australia
9    # Licensed under the Open Software License version 3.0
10    # http://www.opensource.org/licenses/osl-3.0.php
11    #
12    ########################################################
13    
14    __copyright__="""Copyright (c) 2003-2012 by University of Queensland
15    Earth Systems Science Computational Center (ESSCC)
16    http://www.uq.edu.au/esscc
17    Primary Business: Queensland, Australia"""
18    __license__="""Licensed under the Open Software License version 3.0
19    http://www.opensource.org/licenses/osl-3.0.php"""
20    __url__="https://launchpad.net/escript-finley"
21    
22  """  """
23  Provides some tools related to PDEs.  Provides some tools related to PDEs.
24    
25  Currently includes:  Currently includes:
26      - Projector - to project a discontinuous      - Projector - to project a discontinuous function onto a continuous function
27      - Locator - to trace values in data objects at a certain location      - Locator - to trace values in data objects at a certain location
28      - TimeIntegrationManager - to handel extraplotion in time      - TimeIntegrationManager - to handle extrapolation in time
29        - SaddlePointProblem - solver for Saddle point problems using the inexact uszawa scheme
30    
31    :var __author__: name of author
32    :var __copyright__: copyrights
33    :var __license__: licence agreement
34    :var __url__: url entry point on documentation
35    :var __version__: version
36    :var __date__: date of the version
37  """  """
38    
39  import escript  __author__="Lutz Gross, l.gross@uq.edu.au"
40  import linearPDEs  
41  import numarray  
42  import util  from . import escript
43    from . import linearPDEs
44    import numpy
45    from . import util
46    import math
47    
48  class TimeIntegrationManager:  class TimeIntegrationManager:
49    """    """
50    a simple mechanism to manage time dependend values.    A simple mechanism to manage time dependend values.
51    
52    typical usage is:    Typical usage is::
53    
54    dt=0.1 # time increment       dt=0.1 # time increment
55    tm=TimeIntegrationManager(inital_value,p=1)       tm=TimeIntegrationManager(inital_value,p=1)
56    while t<1.       while t<1.
57        v_guess=tm.extrapolate(dt) # extrapolate to t+dt           v_guess=tm.extrapolate(dt) # extrapolate to t+dt
58        v=...           v=...
59        tm.checkin(dt,v)           tm.checkin(dt,v)
60        t+=dt           t+=dt
61    
62    @remark: currently only p=1 is supported.    :note: currently only p=1 is supported.
63    """    """
64    def __init__(self,*inital_values,**kwargs):    def __init__(self,*inital_values,**kwargs):
65       """       """
66       sets up the value manager where inital_value is the initial value and p is order used for extrapolation       Sets up the value manager where ``inital_values`` are the initial values
67         and p is the order used for extrapolation.
68       """       """
69       if kwargs.has_key("p"):       if "p" in kwargs:
70              self.__p=kwargs["p"]              self.__p=kwargs["p"]
71       else:       else:
72              self.__p=1              self.__p=1
73       if kwargs.has_key("time"):       if "time" in kwargs:
74              self.__t=kwargs["time"]              self.__t=kwargs["time"]
75       else:       else:
76              self.__t=0.              self.__t=0.
# Line 50  class TimeIntegrationManager: Line 82  class TimeIntegrationManager:
82    def getTime(self):    def getTime(self):
83        return self.__t        return self.__t
84    
85      def getValue(self):
86          out=self.__v_mem[0]
87          if len(out)==1:
88              return out[0]
89          else:
90              return out
91    
92    def checkin(self,dt,*values):    def checkin(self,dt,*values):
93        """        """
94        adds new values to the manager. the p+1 last value get lost        Adds new values to the manager. The p+1 last values are lost.
95        """        """
96        o=min(self.__order+1,self.__p)        o=min(self.__order+1,self.__p)
97        self.__order=min(self.__order+1,self.__p)        self.__order=min(self.__order+1,self.__p)
# Line 69  class TimeIntegrationManager: Line 108  class TimeIntegrationManager:
108    
109    def extrapolate(self,dt):    def extrapolate(self,dt):
110        """        """
111        extrapolates to dt forward in time.        Extrapolates to ``dt`` forward in time.
112        """        """
113        if self.__order==0:        if self.__order==0:
114           out=self.__v_mem[0]           out=self.__v_mem[0]
# Line 85  class TimeIntegrationManager: Line 124  class TimeIntegrationManager:
124        else:        else:
125           return out           return out
126    
127    
128  class Projector:  class Projector:
129    """    """
130    The Projector is a factory which projects a discontiuous function onto a    The Projector is a factory which projects a discontinuous function onto a
131    continuous function on the a given domain.    continuous function on a given domain.
132    """    """
133    def __init__(self, domain, reduce = True, fast=True):    def __init__(self, domain, reduce=True, fast=True):
134      """      """
135      Create a continuous function space projector for a domain.      Creates a continuous function space projector for a domain.
136    
137      @param domain: Domain of the projection.      :param domain: Domain of the projection.
138      @param reduce: Flag to reduce projection order (default is True)      :param reduce: Flag to reduce projection order
139      @param fast: Flag to use a fast method based on matrix lumping (default is true)      :param fast: Flag to use a fast method based on matrix lumping
140      """      """
141      self.__pde = linearPDEs.LinearPDE(domain)      self.__pde = linearPDEs.LinearPDE(domain)
142      if fast:      if fast:
143        self.__pde.setSolverMethod(linearPDEs.LinearPDE.LUMPING)          self.__pde.getSolverOptions().setSolverMethod(linearPDEs.SolverOptions.LUMPING)
144      self.__pde.setSymmetryOn()      self.__pde.setSymmetryOn()
145      self.__pde.setReducedOrderTo(reduce)      self.__pde.setReducedOrderTo(reduce)
146      self.__pde.setValue(D = 1.)      self.__pde.setValue(D = 1.)
147      return      return
148      def getSolverOptions(self):
149        """
150        Returns the solver options of the PDE solver.
151        
152        :rtype: `linearPDEs.SolverOptions`
153        """
154        return self.__pde.getSolverOptions()
155    
156    def __del__(self):    def getValue(self, input_data):
157      return      """
158        Projects ``input_data`` onto a continuous function.
159    
160        :param input_data: the data to be projected
161        """
162        return self(input_data)
163    
164    def __call__(self, input_data):    def __call__(self, input_data):
165      """      """
166      Projects input_data onto a continuous function      Projects ``input_data`` onto a continuous function.
167    
168      @param input_data: The input_data to be projected.      :param input_data: the data to be projected
169      """      """
170      out=escript.Data(0.,input_data.getShape(),what=escript.ContinuousFunction(self.__pde.getDomain()))      out=escript.Data(0.,input_data.getShape(),self.__pde.getFunctionSpaceForSolution())
171        self.__pde.setValue(Y = escript.Data(), Y_reduced = escript.Data())
172      if input_data.getRank()==0:      if input_data.getRank()==0:
173          self.__pde.setValue(Y = input_data)          self.__pde.setValue(Y = input_data)
174          out=self.__pde.getSolution()          out=self.__pde.getSolution()
# Line 143  class Projector: Line 196  class Projector:
196                      out[i0,i1,i2,i3]=self.__pde.getSolution()                      out[i0,i1,i2,i3]=self.__pde.getSolution()
197      return out      return out
198    
199    class NoPDE:
200         """
201         Solves the following problem for u:
202    
203         *kronecker[i,j]*D[j]*u[j]=Y[i]*
204    
205         with constraint
206    
207         *u[j]=r[j]*  where *q[j]>0*
208    
209         where *D*, *Y*, *r* and *q* are given functions of rank 1.
210    
211         In the case of scalars this takes the form
212    
213         *D*u=Y*
214    
215         with constraint
216    
217         *u=r* where *q>0*
218    
219         where *D*, *Y*, *r* and *q* are given scalar functions.
220    
221         The constraint overwrites any other condition.
222    
223         :note: This class is similar to the `linearPDEs.LinearPDE` class with
224                A=B=C=X=0 but has the intention that all input parameters are given
225                in `Solution` or `ReducedSolution`.
226         """
227         # The whole thing is a bit strange and I blame Rob Woodcock (CSIRO) for
228         # this.
229         def __init__(self,domain,D=None,Y=None,q=None,r=None):
230             """
231             Initializes the problem.
232    
233             :param domain: domain of the PDE
234             :type domain: `Domain`
235             :param D: coefficient of the solution
236             :type D: ``float``, ``int``, ``numpy.ndarray``, `Data`
237             :param Y: right hand side
238             :type Y: ``float``, ``int``, ``numpy.ndarray``, `Data`
239             :param q: location of constraints
240             :type q: ``float``, ``int``, ``numpy.ndarray``, `Data`
241             :param r: value of solution at locations of constraints
242             :type r: ``float``, ``int``, ``numpy.ndarray``, `Data`
243             """
244             self.__domain=domain
245             self.__D=D
246             self.__Y=Y
247             self.__q=q
248             self.__r=r
249             self.__u=None
250             self.__function_space=escript.Solution(self.__domain)
251    
252         def setReducedOn(self):
253             """
254             Sets the `FunctionSpace` of the solution to `ReducedSolution`.
255             """
256             self.__function_space=escript.ReducedSolution(self.__domain)
257             self.__u=None
258    
259         def setReducedOff(self):
260             """
261             Sets the `FunctionSpace` of the solution to `Solution`.
262             """
263             self.__function_space=escript.Solution(self.__domain)
264             self.__u=None
265    
266         def setValue(self,D=None,Y=None,q=None,r=None):
267             """
268             Assigns values to the parameters.
269    
270             :param D: coefficient of the solution
271             :type D: ``float``, ``int``, ``numpy.ndarray``, `Data`
272             :param Y: right hand side
273             :type Y: ``float``, ``int``, ``numpy.ndarray``, `Data`
274             :param q: location of constraints
275             :type q: ``float``, ``int``, ``numpy.ndarray``, `Data`
276             :param r: value of solution at locations of constraints
277             :type r: ``float``, ``int``, ``numpy.ndarray``, `Data`
278             """
279             if not D==None:
280                self.__D=D
281                self.__u=None
282             if not Y==None:
283                self.__Y=Y
284                self.__u=None
285             if not q==None:
286                self.__q=q
287                self.__u=None
288             if not r==None:
289                self.__r=r
290                self.__u=None
291    
292         def getSolution(self):
293             """
294             Returns the solution.
295    
296             :return: the solution of the problem
297             :rtype: `Data` object in the `FunctionSpace` `Solution` or
298                     `ReducedSolution`
299             """
300             if self.__u==None:
301                if self.__D==None:
302                   raise ValueError("coefficient D is undefined")
303                D=escript.Data(self.__D,self.__function_space)
304                if D.getRank()>1:
305                   raise ValueError("coefficient D must have rank 0 or 1")
306                if self.__Y==None:
307                   self.__u=escript.Data(0.,D.getShape(),self.__function_space)
308                else:
309                   self.__u=util.quotient(self.__Y,D)
310                if not self.__q==None:
311                    q=util.wherePositive(escript.Data(self.__q,self.__function_space))
312                    self.__u*=(1.-q)
313                    if not self.__r==None: self.__u+=q*self.__r
314             return self.__u
315    
316  class Locator:  class Locator:
317       """       """
318       Locator provides access to the values of data objects at a given       Locator provides access to the values of data objects at a given spatial
319       spatial coordinate x.         coordinate x.
320        
321       In fact, a Locator object finds the sample in the set of samples of a       In fact, a Locator object finds the sample in the set of samples of a
322       given function space or domain where which is closest to the given       given function space or domain which is closest to the given point x.
      point x.  
323       """       """
324    
325       def __init__(self,where,x=numarray.zeros((3,))):       def __init__(self,where,x=numpy.zeros((3,))):
326         """         """
327         Initializes a Locator to access values in Data objects on the Doamin         Initializes a Locator to access values in Data objects on the Doamin
328         or FunctionSpace where for the sample point which         or FunctionSpace for the sample point which is closest to the given
329         closest to the given point x.         point x.
330    
331           :param where: function space
332           :type where: `escript.FunctionSpace`
333           :param x: location(s) of the Locator
334           :type x: ``numpy.ndarray`` or ``list`` of ``numpy.ndarray``
335         """         """
336         if isinstance(where,escript.FunctionSpace):         if isinstance(where,escript.FunctionSpace):
337            self.__function_space=where            self.__function_space=where
338         else:         else:
339            self.__function_space=escript.ContinuousFunction(where)            self.__function_space=escript.ContinuousFunction(where)
340         self.__id=util.length(x[:self.__function_space.getDim()]-self.__function_space.getX()).mindp()         iterative=False
341           if isinstance(x, list):
342               if len(x)==0:
343                  raise ValueError("At least one point must be given.")
344               try:
345                 iter(x[0])
346                 iterative=True
347               except TypeError:
348                 iterative=False
349           xxx=self.__function_space.getX()
350           if iterative:
351               self.__id=[]
352               for p in x:
353                  self.__id.append(util.length(xxx-p[:self.__function_space.getDim()]).minGlobalDataPoint())
354           else:
355               self.__id=util.length(xxx-x[:self.__function_space.getDim()]).minGlobalDataPoint()
356    
357       def __str__(self):       def __str__(self):
358         """         """
359         Returns the coordinates of the Locator as a string.         Returns the coordinates of the Locator as a string.
360         """         """
361         return "<Locator %s>"%str(self.getX())         x=self.getX()
362           if isinstance(x,list):
363              out="["
364              first=True
365              for xx in x:
366                if not first:
367                    out+=","
368                else:
369                    first=False
370                out+=str(xx)
371              out+="]>"
372           else:
373              out=str(x)
374           return out
375    
376         def getX(self):
377            """
378            Returns the exact coordinates of the Locator.
379            """
380            return self(self.getFunctionSpace().getX())
381    
382       def getFunctionSpace(self):       def getFunctionSpace(self):
383          """          """
384      Returns the function space of the Locator.          Returns the function space of the Locator.
385      """          """
386          return self.__function_space          return self.__function_space
387    
388       def getId(self):       def getId(self,item=None):
389          """          """
390      Returns the identifier of the location.          Returns the identifier of the location.
     """  
         return self.__id  
   
      def getX(self):  
391          """          """
392      Returns the exact coordinates of the Locator.          if item == None:
393      """             return self.__id
394          return self(self.getFunctionSpace().getX())          else:
395               if isinstance(self.__id,list):
396                  return self.__id[item]
397               else:
398                  return self.__id
399    
400    
401       def __call__(self,data):       def __call__(self,data):
402          """          """
403      Returns the value of data at the Locator of a Data object otherwise          Returns the value of data at the Locator of a Data object.
404      the object is returned.          """
     """  
405          return self.getValue(data)          return self.getValue(data)
406    
407       def getValue(self,data):       def getValue(self,data):
408          """          """
409      Returns the value of data at the Locator if data is a Data object          Returns the value of ``data`` at the Locator if ``data`` is a `Data`
410      otherwise the object is returned.          object otherwise the object is returned.
411      """          """
412          if isinstance(data,escript.Data):          if isinstance(data,escript.Data):
413             if data.getFunctionSpace()==self.getFunctionSpace():             dat=util.interpolate(data,self.getFunctionSpace())
414               out=data.convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])             id=self.getId()
415               r=data.getRank()
416               if isinstance(id,list):
417                   out=[]
418                   for i in id:
419                      o=numpy.array(dat.getTupleForGlobalDataPoint(*i))
420                      if data.getRank()==0:
421                         out.append(o[0])
422                      else:
423                         out.append(o)
424                   return out
425             else:             else:
426               out=data.interpolate(self.getFunctionSpace()).convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])               out=numpy.array(dat.getTupleForGlobalDataPoint(*id))
427             if data.getRank()==0:               if data.getRank()==0:
428                return out[0]                  return out[0]
429             else:               else:
430                return out                  return out
431          else:          else:
432             return data             return data
433              
434         def setValue(self, data, v):
435          """
436          Sets the value of the ``data`` at the Locator.
437          """
438          if isinstance(data, escript.Data):
439             if data.getFunctionSpace()!=self.getFunctionSpace():
440               raise TypeError, "setValue: FunctionSpace of Locator and Data object must match."
441             data.expand()  
442             id=self.getId()
443             if isinstance(id, list):
444              for i in id:
445               data._setTupleForGlobalDataPoint(i[1], i[0], v)
446             else:
447               data._setTupleForGlobalDataPoint(id[1], id[0], v)
448          else:
449               raise TypeError, "setValue: Invalid argument type."
450    
451    
452    def getInfLocator(arg):
453        """
454        Return a Locator for a point with the inf value over all arg.
455        """
456        if not isinstance(arg, escript.Data):
457           raise TypeError("getInfLocator: Unknown argument type.")
458        a_inf=util.inf(arg)
459        loc=util.length(arg-a_inf).minGlobalDataPoint() # This gives us the location but not coords
460        x=arg.getFunctionSpace().getX()
461        x_min=x.getTupleForGlobalDataPoint(*loc)
462        return Locator(arg.getFunctionSpace(),x_min)
463    
464    def getSupLocator(arg):
465        """
466        Return a Locator for a point with the sup value over all arg.
467        """
468        if not isinstance(arg, escript.Data):
469           raise TypeError("getInfLocator: Unknown argument type.")
470        a_inf=util.sup(arg)
471        loc=util.length(arg-a_inf).minGlobalDataPoint() # This gives us the location but not coords
472        x=arg.getFunctionSpace().getX()
473        x_min=x.getTupleForGlobalDataPoint(*loc)
474        return Locator(arg.getFunctionSpace(),x_min)
475        
476    
477    class SolverSchemeException(Exception):
478       """
479       This is a generic exception thrown by solvers.
480       """
481       pass
482    
483    class IndefinitePreconditioner(SolverSchemeException):
484       """
485       Exception thrown if the preconditioner is not positive definite.
486       """
487       pass
488    
489    class MaxIterReached(SolverSchemeException):
490       """
491       Exception thrown if the maximum number of iteration steps is reached.
492       """
493       pass
494    
495    class CorrectionFailed(SolverSchemeException):
496       """
497       Exception thrown if no convergence has been achieved in the solution
498       correction scheme.
499       """
500       pass
501    
502    class IterationBreakDown(SolverSchemeException):
503       """
504       Exception thrown if the iteration scheme encountered an incurable breakdown.
505       """
506       pass
507    
508    class NegativeNorm(SolverSchemeException):
509       """
510       Exception thrown if a norm calculation returns a negative norm.
511       """
512       pass
513    
514    def PCG(r, Aprod, x, Msolve, bilinearform, atol=0, rtol=1.e-8, iter_max=100, initial_guess=True, verbose=False):
515       """
516       Solver for
517    
518       *Ax=b*
519    
520       with a symmetric and positive definite operator A (more details required!).
521       It uses the conjugate gradient method with preconditioner M providing an
522       approximation of A.
523    
524       The iteration is terminated if
525    
526       *|r| <= atol+rtol*|r0|*
527    
528       where *r0* is the initial residual and *|.|* is the energy norm. In fact
529    
530       *|r| = sqrt( bilinearform(Msolve(r),r))*
531    
532       For details on the preconditioned conjugate gradient method see the book:
533    
534       I{Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
535       T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
536       C. Romine, and H. van der Vorst}.
537    
538       :param r: initial residual *r=b-Ax*. ``r`` is altered.
539       :type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
540       :param x: an initial guess for the solution
541       :type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
542       :param Aprod: returns the value Ax
543       :type Aprod: function ``Aprod(x)`` where ``x`` is of the same object like
544                    argument ``x``. The returned object needs to be of the same type
545                    like argument ``r``.
546       :param Msolve: solves Mx=r
547       :type Msolve: function ``Msolve(r)`` where ``r`` is of the same type like
548                     argument ``r``. The returned object needs to be of the same
549                     type like argument ``x``.
550       :param bilinearform: inner product ``<x,r>``
551       :type bilinearform: function ``bilinearform(x,r)`` where ``x`` is of the same
552                           type like argument ``x`` and ``r`` is. The returned value
553                           is a ``float``.
554       :param atol: absolute tolerance
555       :type atol: non-negative ``float``
556       :param rtol: relative tolerance
557       :type rtol: non-negative ``float``
558       :param iter_max: maximum number of iteration steps
559       :type iter_max: ``int``
560       :return: the solution approximation and the corresponding residual
561       :rtype: ``tuple``
562       :warning: ``r`` and ``x`` are altered.
563       """
564       iter=0
565       rhat=Msolve(r)
566       d = rhat
567       rhat_dot_r = bilinearform(rhat, r)
568       if rhat_dot_r<0: raise NegativeNorm("negative norm.")
569       norm_r0=math.sqrt(rhat_dot_r)
570       atol2=atol+rtol*norm_r0
571       if atol2<=0:
572          raise ValueError("Non-positive tolarance.")
573       atol2=max(atol2, 100. * util.EPSILON * norm_r0)
574    
575       if verbose: print(("PCG: initial residual norm = %e (absolute tolerance = %e)"%(norm_r0, atol2)))
576    
577    
578       while not math.sqrt(rhat_dot_r) <= atol2:
579           iter+=1
580           if iter  >= iter_max: raise MaxIterReached("maximum number of %s steps reached."%iter_max)
581    
582           q=Aprod(d)
583           alpha = rhat_dot_r / bilinearform(d, q)
584           x += alpha * d
585           if isinstance(q,ArithmeticTuple):
586              r += q * (-alpha)      # Doing it the other way calls the float64.__mul__ not AT.__rmul__
587           else:
588               r += (-alpha) * q
589           rhat=Msolve(r)
590           rhat_dot_r_new = bilinearform(rhat, r)
591           beta = rhat_dot_r_new / rhat_dot_r
592           rhat+=beta * d
593           d=rhat
594    
595           rhat_dot_r = rhat_dot_r_new
596           if rhat_dot_r<0: raise NegativeNorm("negative norm.")
597           if verbose: print(("PCG: iteration step %s: residual norm = %e"%(iter, math.sqrt(rhat_dot_r))))
598       if verbose: print(("PCG: tolerance reached after %s steps."%iter))
599       return x,r,math.sqrt(rhat_dot_r)
600    
601    class Defect(object):
602        """
603        Defines a non-linear defect F(x) of a variable x.
604        """
605        def __init__(self):
606            """
607            Initializes defect.
608            """
609            self.setDerivativeIncrementLength()
610    
611        def bilinearform(self, x0, x1):
612            """
613            Returns the inner product of x0 and x1
614    
615            :param x0: value for x0
616            :param x1: value for x1
617            :return: the inner product of x0 and x1
618            :rtype: ``float``
619            """
620            return 0
621    
622        def norm(self,x):
623            """
624            Returns the norm of argument ``x``.
625    
626            :param x: a value
627            :return: norm of argument x
628            :rtype: ``float``
629            :note: by default ``sqrt(self.bilinearform(x,x)`` is returned.
630            """
631            s=self.bilinearform(x,x)
632            if s<0: raise NegativeNorm("negative norm.")
633            return math.sqrt(s)
634    
635        def eval(self,x):
636            """
637            Returns the value F of a given ``x``.
638    
639            :param x: value for which the defect ``F`` is evaluated
640            :return: value of the defect at ``x``
641            """
642            return 0
643    
644        def __call__(self,x):
645            return self.eval(x)
646    
647        def setDerivativeIncrementLength(self,inc=1000.*math.sqrt(util.EPSILON)):
648            """
649            Sets the relative length of the increment used to approximate the
650            derivative of the defect. The increment is inc*norm(x)/norm(v)*v in the
651            direction of v with x as a starting point.
652    
653            :param inc: relative increment length
654            :type inc: positive ``float``
655            """
656            if inc<=0: raise ValueError("positive increment required.")
657            self.__inc=inc
658    
659        def getDerivativeIncrementLength(self):
660            """
661            Returns the relative increment length used to approximate the
662            derivative of the defect.
663            :return: value of the defect at ``x``
664            :rtype: positive ``float``
665            """
666            return self.__inc
667    
668        def derivative(self, F0, x0, v, v_is_normalised=True):
669            """
670            Returns the directional derivative at ``x0`` in the direction of ``v``.
671    
672            :param F0: value of this defect at x0
673            :param x0: value at which derivative is calculated
674            :param v: direction
675            :param v_is_normalised: True to indicate that ``v`` is nomalized
676                                    (self.norm(v)=0)
677            :return: derivative of this defect at x0 in the direction of ``v``
678            :note: by default numerical evaluation (self.eval(x0+eps*v)-F0)/eps is
679                   used but this method maybe overwritten to use exact evaluation.
680            """
681            normx=self.norm(x0)
682            if normx>0:
683                 epsnew = self.getDerivativeIncrementLength() * normx
684            else:
685                 epsnew = self.getDerivativeIncrementLength()
686            if not v_is_normalised:
687                normv=self.norm(v)
688                if normv<=0:
689                   return F0*0
690                else:
691                   epsnew /= normv
692            F1=self.eval(x0 + epsnew * v)
693            return (F1-F0)/epsnew
694    
695    ######################################
696    def NewtonGMRES(defect, x, iter_max=100, sub_iter_max=20, atol=0,rtol=1.e-4, subtol_max=0.5, gamma=0.9, verbose=False):
697       """
698       Solves a non-linear problem *F(x)=0* for unknown *x* using the stopping
699       criterion:
700    
701       *norm(F(x) <= atol + rtol * norm(F(x0)*
702    
703       where *x0* is the initial guess.
704    
705       :param defect: object defining the function *F*. ``defect.norm`` defines the
706                      *norm* used in the stopping criterion.
707       :type defect: `Defect`
708       :param x: initial guess for the solution, ``x`` is altered.
709       :type x: any object type allowing basic operations such as
710                ``numpy.ndarray``, `Data`
711       :param iter_max: maximum number of iteration steps
712       :type iter_max: positive ``int``
713       :param sub_iter_max: maximum number of inner iteration steps
714       :type sub_iter_max: positive ``int``
715       :param atol: absolute tolerance for the solution
716       :type atol: positive ``float``
717       :param rtol: relative tolerance for the solution
718       :type rtol: positive ``float``
719       :param gamma: tolerance safety factor for inner iteration
720       :type gamma: positive ``float``, less than 1
721       :param subtol_max: upper bound for inner tolerance
722       :type subtol_max: positive ``float``, less than 1
723       :return: an approximation of the solution with the desired accuracy
724       :rtype: same type as the initial guess
725       """
726       lmaxit=iter_max
727       if atol<0: raise ValueError("atol needs to be non-negative.")
728       if rtol<0: raise ValueError("rtol needs to be non-negative.")
729       if rtol+atol<=0: raise ValueError("rtol or atol needs to be non-negative.")
730       if gamma<=0 or gamma>=1: raise ValueError("tolerance safety factor for inner iteration (gamma =%s) needs to be positive and less than 1."%gamma)
731       if subtol_max<=0 or subtol_max>=1: raise ValueError("upper bound for inner tolerance for inner iteration (subtol_max =%s) needs to be positive and less than 1."%subtol_max)
732    
733       F=defect(x)
734       fnrm=defect.norm(F)
735       stop_tol=atol + rtol*fnrm
736       subtol=subtol_max
737       if verbose: print(("NewtonGMRES: initial residual = %e."%fnrm))
738       if verbose: print(("             tolerance = %e."%subtol))
739       iter=1
740       #
741       # main iteration loop
742       #
743       while not fnrm<=stop_tol:
744                if iter  >= iter_max: raise MaxIterReached("maximum number of %s steps reached."%iter_max)
745                #
746            #   adjust subtol_
747            #
748                if iter > 1:
749                   rat=fnrm/fnrmo
750                   subtol_old=subtol
751                   subtol=gamma*rat**2
752                   if gamma*subtol_old**2 > .1: subtol=max(subtol,gamma*subtol_old**2)
753                   subtol=max(min(subtol,subtol_max), .5*stop_tol/fnrm)
754            #
755            # calculate newton increment xc
756                #     if iter_max in __FDGMRES is reached MaxIterReached is thrown
757                #     if iter_restart -1 is returned as sub_iter
758                #     if  atol is reached sub_iter returns the numer of steps performed to get there
759                #
760                #
761                if verbose: print(("             subiteration (GMRES) is called with relative tolerance %e."%subtol))
762                try:
763                   xc, sub_iter=__FDGMRES(F, defect, x, subtol*fnrm, iter_max=iter_max-iter, iter_restart=sub_iter_max)
764                except MaxIterReached:
765                   raise MaxIterReached("maximum number of %s steps reached."%iter_max)
766                if sub_iter<0:
767                   iter+=sub_iter_max
768                else:
769                   iter+=sub_iter
770                # ====
771                x+=xc
772                F=defect(x)
773                iter+=1
774                fnrmo, fnrm=fnrm, defect.norm(F)
775                if verbose: print(("             step %s: residual %e."%(iter,fnrm)))
776       if verbose: print(("NewtonGMRES: completed after %s steps."%iter))
777       return x
778    
779    def __givapp(c,s,vin):
780        """
781        Applies a sequence of Givens rotations (c,s) recursively to the vector
782        ``vin``
783    
784        :warning: ``vin`` is altered.
785        """
786        vrot=vin
787        if isinstance(c,float):
788            vrot=[c*vrot[0]-s*vrot[1],s*vrot[0]+c*vrot[1]]
789        else:
790            for i in range(len(c)):
791                w1=c[i]*vrot[i]-s[i]*vrot[i+1]
792                w2=s[i]*vrot[i]+c[i]*vrot[i+1]
793                vrot[i]=w1
794                vrot[i+1]=w2
795        return vrot
796    
797    def __FDGMRES(F0, defect, x0, atol, iter_max=100, iter_restart=20):
798       h=numpy.zeros((iter_restart,iter_restart),numpy.float64)
799       c=numpy.zeros(iter_restart,numpy.float64)
800       s=numpy.zeros(iter_restart,numpy.float64)
801       g=numpy.zeros(iter_restart,numpy.float64)
802       v=[]
803    
804       rho=defect.norm(F0)
805       if rho<=0.: return x0*0
806    
807       v.append(-F0/rho)
808       g[0]=rho
809       iter=0
810       while rho > atol and iter<iter_restart-1:
811            if iter  >= iter_max:
812                raise MaxIterReached("maximum number of %s steps reached."%iter_max)
813    
814            p=defect.derivative(F0,x0,v[iter], v_is_normalised=True)
815            v.append(p)
816    
817            v_norm1=defect.norm(v[iter+1])
818    
819            # Modified Gram-Schmidt
820            for j in range(iter+1):
821                h[j,iter]=defect.bilinearform(v[j],v[iter+1])
822                v[iter+1]-=h[j,iter]*v[j]
823    
824            h[iter+1,iter]=defect.norm(v[iter+1])
825            v_norm2=h[iter+1,iter]
826    
827            # Reorthogonalize if needed
828            if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
829                for j in range(iter+1):
830                    hr=defect.bilinearform(v[j],v[iter+1])
831                    h[j,iter]=h[j,iter]+hr
832                    v[iter+1] -= hr*v[j]
833    
834                v_norm2=defect.norm(v[iter+1])
835                h[iter+1,iter]=v_norm2
836            #   watch out for happy breakdown
837            if not v_norm2 == 0:
838                v[iter+1]=v[iter+1]/h[iter+1,iter]
839    
840            #   Form and store the information for the new Givens rotation
841            if iter > 0 :
842                hhat=numpy.zeros(iter+1,numpy.float64)
843                for i in range(iter+1) : hhat[i]=h[i,iter]
844                hhat=__givapp(c[0:iter],s[0:iter],hhat);
845                for i in range(iter+1) : h[i,iter]=hhat[i]
846    
847            mu=math.sqrt(h[iter,iter]*h[iter,iter]+h[iter+1,iter]*h[iter+1,iter])
848    
849            if mu!=0 :
850                c[iter]=h[iter,iter]/mu
851                s[iter]=-h[iter+1,iter]/mu
852                h[iter,iter]=c[iter]*h[iter,iter]-s[iter]*h[iter+1,iter]
853                h[iter+1,iter]=0.0
854                gg=__givapp(c[iter],s[iter],[g[iter],g[iter+1]])
855                g[iter]=gg[0]
856                g[iter+1]=gg[1]
857    
858            # Update the residual norm
859            rho=abs(g[iter+1])
860            iter+=1
861    
862       # At this point either iter > iter_max or rho < tol.
863       # It's time to compute x and leave.
864       if iter > 0 :
865         y=numpy.zeros(iter,numpy.float64)
866         y[iter-1] = g[iter-1] / h[iter-1,iter-1]
867         if iter > 1 :
868            i=iter-2
869            while i>=0 :
870              y[i] = ( g[i] - numpy.dot(h[i,i+1:iter], y[i+1:iter])) / h[i,i]
871              i=i-1
872         xhat=v[iter-1]*y[iter-1]
873         for i in range(iter-1):
874           xhat += v[i]*y[i]
875       else :
876          xhat=v[0] * 0
877    
878       if iter<iter_restart-1:
879          stopped=iter
880       else:
881          stopped=-1
882    
883       return xhat,stopped
884    
885    def GMRES(r, Aprod, x, bilinearform, atol=0, rtol=1.e-8, iter_max=100, iter_restart=20, verbose=False,P_R=None):
886       """
887       Solver for
888    
889       *Ax=b*
890    
891       with a general operator A (more details required!).
892       It uses the generalized minimum residual method (GMRES).
893    
894       The iteration is terminated if
895    
896       *|r| <= atol+rtol*|r0|*
897    
898       where *r0* is the initial residual and *|.|* is the energy norm. In fact
899    
900       *|r| = sqrt( bilinearform(r,r))*
901    
902       :param r: initial residual *r=b-Ax*. ``r`` is altered.
903       :type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
904       :param x: an initial guess for the solution
905       :type x: same like ``r``
906       :param Aprod: returns the value Ax
907       :type Aprod: function ``Aprod(x)`` where ``x`` is of the same object like
908                    argument ``x``. The returned object needs to be of the same
909                    type like argument ``r``.
910       :param bilinearform: inner product ``<x,r>``
911       :type bilinearform: function ``bilinearform(x,r)`` where ``x`` is of the same
912                           type like argument ``x`` and ``r``. The returned value is
913                           a ``float``.
914       :param atol: absolute tolerance
915       :type atol: non-negative ``float``
916       :param rtol: relative tolerance
917       :type rtol: non-negative ``float``
918       :param iter_max: maximum number of iteration steps
919       :type iter_max: ``int``
920       :param iter_restart: in order to save memory the orthogonalization process
921                            is terminated after ``iter_restart`` steps and the
922                            iteration is restarted.
923       :type iter_restart: ``int``
924       :return: the solution approximation and the corresponding residual
925       :rtype: ``tuple``
926       :warning: ``r`` and ``x`` are altered.
927       """
928       m=iter_restart
929       restarted=False
930       iter=0
931       if rtol>0:
932          r_dot_r = bilinearform(r, r)
933          if r_dot_r<0: raise NegativeNorm("negative norm.")
934          atol2=atol+rtol*math.sqrt(r_dot_r)
935          if verbose: print(("GMRES: norm of right hand side = %e (absolute tolerance = %e)"%(math.sqrt(r_dot_r), atol2)))
936       else:
937          atol2=atol
938          if verbose: print(("GMRES: absolute tolerance = %e"%atol2))
939       if atol2<=0:
940          raise ValueError("Non-positive tolarance.")
941    
942       while True:
943          if iter  >= iter_max: raise MaxIterReached("maximum number of %s steps reached"%iter_max)
944          if restarted:
945             r2 = r-Aprod(x-x2)
946          else:
947             r2=1*r
948          x2=x*1.
949          x,stopped=_GMRESm(r2, Aprod, x, bilinearform, atol2, iter_max=iter_max-iter, iter_restart=m, verbose=verbose,P_R=P_R)
950          iter+=iter_restart
951          if stopped: break
952          if verbose: print("GMRES: restart.")
953          restarted=True
954       if verbose: print("GMRES: tolerance has been reached.")
955       return x
956    
957    def _GMRESm(r, Aprod, x, bilinearform, atol, iter_max=100, iter_restart=20, verbose=False, P_R=None):
958       iter=0
959    
960       h=numpy.zeros((iter_restart+1,iter_restart),numpy.float64)
961       c=numpy.zeros(iter_restart,numpy.float64)
962       s=numpy.zeros(iter_restart,numpy.float64)
963       g=numpy.zeros(iter_restart+1,numpy.float64)
964       v=[]
965    
966       r_dot_r = bilinearform(r, r)
967       if r_dot_r<0: raise NegativeNorm("negative norm.")
968       rho=math.sqrt(r_dot_r)
969    
970       v.append(r/rho)
971       g[0]=rho
972    
973       if verbose: print(("GMRES: initial residual %e (absolute tolerance = %e)"%(rho,atol)))
974       while not (rho<=atol or iter==iter_restart):
975    
976            if iter  >= iter_max: raise MaxIterReached("maximum number of %s steps reached."%iter_max)
977    
978            if P_R!=None:
979                p=Aprod(P_R(v[iter]))
980            else:
981                p=Aprod(v[iter])
982            v.append(p)
983    
984            v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
985    
986    # Modified Gram-Schmidt
987            for j in range(iter+1):
988              h[j,iter]=bilinearform(v[j],v[iter+1])
989              v[iter+1]-=h[j,iter]*v[j]
990    
991            h[iter+1,iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
992            v_norm2=h[iter+1,iter]
993    
994    # Reorthogonalize if needed
995            if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
996             for j in range(iter+1):
997                hr=bilinearform(v[j],v[iter+1])
998                h[j,iter]=h[j,iter]+hr
999                v[iter+1] -= hr*v[j]
1000    
1001             v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
1002             h[iter+1,iter]=v_norm2
1003    
1004    #   watch out for happy breakdown
1005            if not v_norm2 == 0:
1006             v[iter+1]=v[iter+1]/h[iter+1,iter]
1007    
1008    #   Form and store the information for the new Givens rotation
1009            if iter > 0: h[:iter+1,iter]=__givapp(c[:iter],s[:iter],h[:iter+1,iter])
1010            mu=math.sqrt(h[iter,iter]*h[iter,iter]+h[iter+1,iter]*h[iter+1,iter])
1011    
1012            if mu!=0 :
1013                    c[iter]=h[iter,iter]/mu
1014                    s[iter]=-h[iter+1,iter]/mu
1015                    h[iter,iter]=c[iter]*h[iter,iter]-s[iter]*h[iter+1,iter]
1016                    h[iter+1,iter]=0.0
1017                    gg=__givapp(c[iter],s[iter],[g[iter],g[iter+1]])
1018                    g[iter]=gg[0]
1019                    g[iter+1]=gg[1]
1020    # Update the residual norm
1021    
1022            rho=abs(g[iter+1])
1023            if verbose: print(("GMRES: iteration step %s: residual %e"%(iter,rho)))
1024            iter+=1
1025    
1026    # At this point either iter > iter_max or rho < tol.
1027    # It's time to compute x and leave.
1028    
1029       if verbose: print(("GMRES: iteration stopped after %s step."%iter))
1030       if iter > 0 :
1031         y=numpy.zeros(iter,numpy.float64)
1032         y[iter-1] = g[iter-1] / h[iter-1,iter-1]
1033         if iter > 1 :
1034            i=iter-2
1035            while i>=0 :
1036              y[i] = ( g[i] - numpy.dot(h[i,i+1:iter], y[i+1:iter])) / h[i,i]
1037              i=i-1
1038         xhat=v[iter-1]*y[iter-1]
1039         for i in range(iter-1):
1040           xhat += v[i]*y[i]
1041       else:
1042         xhat=v[0] * 0
1043       if P_R!=None:
1044          x += P_R(xhat)
1045       else:
1046          x += xhat
1047       if iter<iter_restart-1:
1048          stopped=True
1049       else:
1050          stopped=False
1051    
1052       return x,stopped
1053    
1054    def MINRES(r, Aprod, x, Msolve, bilinearform, atol=0, rtol=1.e-8, iter_max=100):
1055        """
1056        Solver for
1057    
1058        *Ax=b*
1059    
1060        with a symmetric and positive definite operator A (more details required!).
1061        It uses the minimum residual method (MINRES) with preconditioner M
1062        providing an approximation of A.
1063    
1064        The iteration is terminated if
1065    
1066        *|r| <= atol+rtol*|r0|*
1067    
1068        where *r0* is the initial residual and *|.|* is the energy norm. In fact
1069    
1070        *|r| = sqrt( bilinearform(Msolve(r),r))*
1071    
1072        For details on the preconditioned conjugate gradient method see the book:
1073    
1074        I{Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
1075        T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
1076        C. Romine, and H. van der Vorst}.
1077    
1078        :param r: initial residual *r=b-Ax*. ``r`` is altered.
1079        :type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1080        :param x: an initial guess for the solution
1081        :type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1082        :param Aprod: returns the value Ax
1083        :type Aprod: function ``Aprod(x)`` where ``x`` is of the same object like
1084                     argument ``x``. The returned object needs to be of the same
1085                     type like argument ``r``.
1086        :param Msolve: solves Mx=r
1087        :type Msolve: function ``Msolve(r)`` where ``r`` is of the same type like
1088                      argument ``r``. The returned object needs to be of the same
1089                      type like argument ``x``.
1090        :param bilinearform: inner product ``<x,r>``
1091        :type bilinearform: function ``bilinearform(x,r)`` where ``x`` is of the same
1092                            type like argument ``x`` and ``r`` is. The returned value
1093                            is a ``float``.
1094        :param atol: absolute tolerance
1095        :type atol: non-negative ``float``
1096        :param rtol: relative tolerance
1097        :type rtol: non-negative ``float``
1098        :param iter_max: maximum number of iteration steps
1099        :type iter_max: ``int``
1100        :return: the solution approximation and the corresponding residual
1101        :rtype: ``tuple``
1102        :warning: ``r`` and ``x`` are altered.
1103        """
1104        #------------------------------------------------------------------
1105        # Set up y and v for the first Lanczos vector v1.
1106        # y  =  beta1 P' v1,  where  P = C**(-1).
1107        # v is really P' v1.
1108        #------------------------------------------------------------------
1109        r1    = r
1110        y = Msolve(r)
1111        beta1 = bilinearform(y,r)
1112    
1113        if beta1< 0: raise NegativeNorm("negative norm.")
1114    
1115        #  If r = 0 exactly, stop with x
1116        if beta1==0: return x
1117    
1118        if beta1> 0: beta1  = math.sqrt(beta1)
1119    
1120        #------------------------------------------------------------------
1121        # Initialize quantities.
1122        # ------------------------------------------------------------------
1123        iter   = 0
1124        Anorm = 0
1125        ynorm = 0
1126        oldb   = 0
1127        beta   = beta1
1128        dbar   = 0
1129        epsln  = 0
1130        phibar = beta1
1131        rhs1   = beta1
1132        rhs2   = 0
1133        rnorm  = phibar
1134        tnorm2 = 0
1135        ynorm2 = 0
1136        cs     = -1
1137        sn     = 0
1138        w      = r*0.
1139        w2     = r*0.
1140        r2     = r1
1141        eps    = 0.0001
1142    
1143        #---------------------------------------------------------------------
1144        # Main iteration loop.
1145        # --------------------------------------------------------------------
1146        while not rnorm<=atol+rtol*Anorm*ynorm:    #  checks ||r|| < (||A|| ||x||) * TOL
1147    
1148            if iter  >= iter_max: raise MaxIterReached("maximum number of %s steps reached."%iter_max)
1149            iter    = iter  +  1
1150    
1151            #-----------------------------------------------------------------
1152            # Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,...
1153            # The general iteration is similar to the case k = 1 with v0 = 0:
1154            #
1155            #   p1      = Operator * v1  -  beta1 * v0,
1156            #   alpha1  = v1'p1,
1157            #   q2      = p2  -  alpha1 * v1,
1158            #   beta2^2 = q2'q2,
1159            #   v2      = (1/beta2) q2.
1160            #
1161            # Again, y = betak P vk,  where  P = C**(-1).
1162            #-----------------------------------------------------------------
1163            s = 1/beta                 # Normalize previous vector (in y).
1164            v = s*y                    # v = vk if P = I
1165    
1166            y      = Aprod(v)
1167    
1168            if iter >= 2:
1169              y = y - (beta/oldb)*r1
1170    
1171            alfa   = bilinearform(v,y)              # alphak
1172            y      += (- alfa/beta)*r2
1173            r1     = r2
1174            r2     = y
1175            y = Msolve(r2)
1176            oldb   = beta                           # oldb = betak
1177            beta   = bilinearform(y,r2)             # beta = betak+1^2
1178            if beta < 0: raise NegativeNorm("negative norm.")
1179    
1180            beta   = math.sqrt( beta )
1181            tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta
1182    
1183            if iter==1:                 # Initialize a few things.
1184              gmax   = abs( alfa )      # alpha1
1185              gmin   = gmax             # alpha1
1186    
1187            # Apply previous rotation Qk-1 to get
1188            #   [deltak epslnk+1] = [cs  sn][dbark    0   ]
1189            #   [gbar k dbar k+1]   [sn -cs][alfak betak+1].
1190    
1191            oldeps = epsln
1192            delta  = cs * dbar  +  sn * alfa  # delta1 = 0         deltak
1193            gbar   = sn * dbar  -  cs * alfa  # gbar 1 = alfa1     gbar k
1194            epsln  =               sn * beta  # epsln2 = 0         epslnk+1
1195            dbar   =            -  cs * beta  # dbar 2 = beta2     dbar k+1
1196    
1197            # Compute the next plane rotation Qk
1198    
1199            gamma  = math.sqrt(gbar*gbar+beta*beta)  # gammak
1200            gamma  = max(gamma,eps)
1201            cs     = gbar / gamma             # ck
1202            sn     = beta / gamma             # sk
1203            phi    = cs * phibar              # phik
1204            phibar = sn * phibar              # phibark+1
1205    
1206            # Update  x.
1207    
1208            denom = 1/gamma
1209            w1    = w2
1210            w2    = w
1211            w     = (v - oldeps*w1 - delta*w2) * denom
1212            x     +=  phi*w
1213    
1214            # Go round again.
1215    
1216            gmax   = max(gmax,gamma)
1217            gmin   = min(gmin,gamma)
1218            z      = rhs1 / gamma
1219            ynorm2 = z*z  +  ynorm2
1220            rhs1   = rhs2 -  delta*z
1221            rhs2   =      -  epsln*z
1222    
1223            # Estimate various norms and test for convergence.
1224    
1225            Anorm  = math.sqrt( tnorm2 )
1226            ynorm  = math.sqrt( ynorm2 )
1227    
1228            rnorm  = phibar
1229    
1230        return x
1231    
1232    def TFQMR(r, Aprod, x, bilinearform, atol=0, rtol=1.e-8, iter_max=100):
1233      """
1234      Solver for
1235    
1236      *Ax=b*
1237    
1238      with a general operator A (more details required!).
1239      It uses the Transpose-Free Quasi-Minimal Residual method (TFQMR).
1240    
1241      The iteration is terminated if
1242    
1243      *|r| <= atol+rtol*|r0|*
1244    
1245      where *r0* is the initial residual and *|.|* is the energy norm. In fact
1246    
1247      *|r| = sqrt( bilinearform(r,r))*
1248    
1249      :param r: initial residual *r=b-Ax*. ``r`` is altered.
1250      :type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1251      :param x: an initial guess for the solution
1252      :type x: same like ``r``
1253      :param Aprod: returns the value Ax
1254      :type Aprod: function ``Aprod(x)`` where ``x`` is of the same object like
1255                   argument ``x``. The returned object needs to be of the same type
1256                   like argument ``r``.
1257      :param bilinearform: inner product ``<x,r>``
1258      :type bilinearform: function ``bilinearform(x,r)`` where ``x`` is of the same
1259                          type like argument ``x`` and ``r``. The returned value is
1260                          a ``float``.
1261      :param atol: absolute tolerance
1262      :type atol: non-negative ``float``
1263      :param rtol: relative tolerance
1264      :type rtol: non-negative ``float``
1265      :param iter_max: maximum number of iteration steps
1266      :type iter_max: ``int``
1267      :rtype: ``tuple``
1268      :warning: ``r`` and ``x`` are altered.
1269      """
1270      u1=0
1271      u2=0
1272      y1=0
1273      y2=0
1274    
1275      w = r
1276      y1 = r
1277      iter = 0
1278      d = 0
1279      v = Aprod(y1)
1280      u1 = v
1281    
1282      theta = 0.0;
1283      eta = 0.0;
1284      rho=bilinearform(r,r)
1285      if rho < 0: raise NegativeNorm("negative norm.")
1286      tau = math.sqrt(rho)
1287      norm_r0=tau
1288      while tau>atol+rtol*norm_r0:
1289        if iter  >= iter_max: raise MaxIterReached("maximum number of %s steps reached."%iter_max)
1290    
1291        sigma = bilinearform(r,v)
1292        if sigma == 0.0: raise IterationBreakDown('TFQMR breakdown, sigma=0')
1293    
1294        alpha = rho / sigma
1295    
1296        for j in range(2):
1297    #
1298    #   Compute y2 and u2 only if you have to
1299    #
1300          if ( j == 1 ):
1301            y2 = y1 - alpha * v
1302            u2 = Aprod(y2)
1303    
1304          m = 2 * (iter+1) - 2 + (j+1)
1305          if j==0:
1306             w = w - alpha * u1
1307             d = y1 + ( theta * theta * eta / alpha ) * d
1308          if j==1:
1309             w = w - alpha * u2
1310             d = y2 + ( theta * theta * eta / alpha ) * d
1311    
1312          theta = math.sqrt(bilinearform(w,w))/ tau
1313          c = 1.0 / math.sqrt ( 1.0 + theta * theta )
1314          tau = tau * theta * c
1315          eta = c * c * alpha
1316          x = x + eta * d
1317    #
1318    #  Try to terminate the iteration at each pass through the loop
1319    #
1320        if rho == 0.0: raise IterationBreakDown('TFQMR breakdown, rho=0')
1321    
1322        rhon = bilinearform(r,w)
1323        beta = rhon / rho;
1324        rho = rhon;
1325        y1 = w + beta * y2;
1326        u1 = Aprod(y1)
1327        v = u1 + beta * ( u2 + beta * v )
1328    
1329        iter += 1
1330    
1331      return x
1332    
1333    
1334    #############################################
1335    
1336    class ArithmeticTuple(object):
1337       """
1338       Tuple supporting inplace update x+=y and scaling x=a*y where ``x,y`` is an
1339       ArithmeticTuple and ``a`` is a float.
1340    
1341       Example of usage::
1342    
1343           from esys.escript import Data
1344           from numpy import array
1345           a=Data(...)
1346           b=array([1.,4.])
1347           x=ArithmeticTuple(a,b)
1348           y=5.*x
1349    
1350       """
1351       def __init__(self,*args):
1352           """
1353           Initializes object with elements ``args``.
1354    
1355           :param args: tuple of objects that support inplace add (x+=y) and
1356                        scaling (x=a*y)
1357           """
1358           self.__items=list(args)
1359    
1360       def __len__(self):
1361           """
1362           Returns the number of items.
1363    
1364           :return: number of items
1365           :rtype: ``int``
1366           """
1367           return len(self.__items)
1368    
1369       def __getitem__(self,index):
1370           """
1371           Returns item at specified position.
1372    
1373           :param index: index of item to be returned
1374           :type index: ``int``
1375           :return: item with index ``index``
1376           """
1377           return self.__items.__getitem__(index)
1378    
1379       def __mul__(self,other):
1380           """
1381           Scales by ``other`` from the right.
1382    
1383           :param other: scaling factor
1384           :type other: ``float``
1385           :return: itemwise self*other
1386           :rtype: `ArithmeticTuple`
1387           """
1388           out=[]
1389           try:
1390               l=len(other)
1391               if l!=len(self):
1392                   raise ValueError("length of arguments don't match.")
1393               for i in range(l): out.append(self[i]*other[i])
1394           except TypeError:
1395               for i in range(len(self)): out.append(self[i]*other)
1396           return ArithmeticTuple(*tuple(out))
1397    
1398       def __rmul__(self,other):
1399           """
1400           Scales by ``other`` from the left.
1401    
1402           :param other: scaling factor
1403           :type other: ``float``
1404           :return: itemwise other*self
1405           :rtype: `ArithmeticTuple`
1406           """
1407           out=[]
1408           try:
1409               l=len(other)
1410               if l!=len(self):
1411                   raise ValueError("length of arguments don't match.")
1412               for i in range(l): out.append(other[i]*self[i])
1413           except TypeError:
1414               for i in range(len(self)): out.append(other*self[i])
1415           return ArithmeticTuple(*tuple(out))
1416    
1417       def __div__(self,other):
1418           """
1419           Scales by (1/``other``) from the right.
1420    
1421           :param other: scaling factor
1422           :type other: ``float``
1423           :return: itemwise self/other
1424           :rtype: `ArithmeticTuple`
1425           """
1426           return self*(1/other)
1427    
1428       def __rdiv__(self,other):
1429           """
1430           Scales by (1/``other``) from the left.
1431    
1432           :param other: scaling factor
1433           :type other: ``float``
1434           :return: itemwise other/self
1435           :rtype: `ArithmeticTuple`
1436           """
1437           out=[]
1438           try:
1439               l=len(other)
1440               if l!=len(self):
1441                   raise ValueError("length of arguments don't match.")
1442               for i in range(l): out.append(other[i]/self[i])
1443           except TypeError:
1444               for i in range(len(self)): out.append(other/self[i])
1445           return ArithmeticTuple(*tuple(out))
1446    
1447       def __iadd__(self,other):
1448           """
1449           Inplace addition of ``other`` to self.
1450    
1451           :param other: increment
1452           :type other: ``ArithmeticTuple``
1453           """
1454           if len(self) != len(other):
1455               raise ValueError("tuple lengths must match.")
1456           for i in range(len(self)):
1457               self.__items[i]+=other[i]
1458           return self
1459    
1460       def __add__(self,other):
1461           """
1462           Adds ``other`` to self.
1463    
1464           :param other: increment
1465           :type other: ``ArithmeticTuple``
1466           """
1467           out=[]
1468           try:
1469               l=len(other)
1470               if l!=len(self):
1471                   raise ValueError("length of arguments don't match.")
1472               for i in range(l): out.append(self[i]+other[i])
1473           except TypeError:
1474               for i in range(len(self)): out.append(self[i]+other)
1475           return ArithmeticTuple(*tuple(out))
1476    
1477       def __sub__(self,other):
1478           """
1479           Subtracts ``other`` from self.
1480    
1481           :param other: decrement
1482           :type other: ``ArithmeticTuple``
1483           """
1484           out=[]
1485           try:
1486               l=len(other)
1487               if l!=len(self):
1488                   raise ValueError("length of arguments don't match.")
1489               for i in range(l): out.append(self[i]-other[i])
1490           except TypeError:
1491               for i in range(len(self)): out.append(self[i]-other)
1492           return ArithmeticTuple(*tuple(out))
1493    
1494       def __isub__(self,other):
1495           """
1496           Inplace subtraction of ``other`` from self.
1497    
1498           :param other: decrement
1499           :type other: ``ArithmeticTuple``
1500           """
1501           if len(self) != len(other):
1502               raise ValueError("tuple length must match.")
1503           for i in range(len(self)):
1504               self.__items[i]-=other[i]
1505           return self
1506    
1507       def __neg__(self):
1508           """
1509           Negates values.
1510           """
1511           out=[]
1512           for i in range(len(self)):
1513               out.append(-self[i])
1514           return ArithmeticTuple(*tuple(out))
1515    
1516    
1517    class HomogeneousSaddlePointProblem(object):
1518          """
1519          This class provides a framework for solving linear homogeneous saddle
1520          point problems of the form::
1521    
1522              *Av+B^*p=f*
1523              *Bv     =0*
1524    
1525          for the unknowns *v* and *p* and given operators *A* and *B* and
1526          given right hand side *f*. *B^** is the adjoint operator of *B*.
1527          *A* may depend weakly on *v* and *p*.
1528          """
1529          def __init__(self, **kwargs):
1530            """
1531            initializes the saddle point problem
1532            """
1533            self.resetControlParameters()
1534            self.setTolerance()
1535            self.setAbsoluteTolerance()
1536          def resetControlParameters(self, K_p=1., K_v=1., rtol_max=0.01, rtol_min = 1.e-7, chi_max=0.5, reduction_factor=0.3, theta = 0.1):
1537             """
1538             sets a control parameter
1539    
1540             :param K_p: initial value for constant to adjust pressure tolerance
1541             :type K_p: ``float``
1542             :param K_v: initial value for constant to adjust velocity tolerance
1543             :type K_v: ``float``
1544             :param rtol_max: maximuim relative tolerance used to calculate presssure and velocity increment.
1545             :type rtol_max: ``float``
1546             :param chi_max: maximum tolerable converegence rate.
1547             :type chi_max: ``float``
1548             :param reduction_factor: reduction factor for adjustment factors.
1549             :type reduction_factor: ``float``
1550             """
1551             self.setControlParameter(K_p, K_v, rtol_max, rtol_min, chi_max, reduction_factor, theta)
1552    
1553          def setControlParameter(self,K_p=None, K_v=None, rtol_max=None, rtol_min=None, chi_max=None, reduction_factor=None, theta=None):
1554             """
1555             sets a control parameter
1556    
1557    
1558             :param K_p: initial value for constant to adjust pressure tolerance
1559             :type K_p: ``float``
1560             :param K_v: initial value for constant to adjust velocity tolerance
1561             :type K_v: ``float``
1562             :param rtol_max: maximuim relative tolerance used to calculate presssure and velocity increment.
1563             :type rtol_max: ``float``
1564             :param chi_max: maximum tolerable converegence rate.
1565             :type chi_max: ``float``
1566             :type reduction_factor: ``float``
1567             """
1568             if not K_p == None:
1569                if K_p<1:
1570                   raise ValueError("K_p need to be greater or equal to 1.")
1571             else:
1572                K_p=self.__K_p
1573    
1574             if not K_v == None:
1575                if K_v<1:
1576                   raise ValueError("K_v need to be greater or equal to 1.")
1577             else:
1578                K_v=self.__K_v
1579    
1580             if not rtol_max == None:
1581                if rtol_max<=0 or rtol_max>=1:
1582                   raise ValueError("rtol_max needs to be positive and less than 1.")
1583             else:
1584                rtol_max=self.__rtol_max
1585    
1586             if not rtol_min == None:
1587                if rtol_min<=0 or rtol_min>=1:
1588                   raise ValueError("rtol_min needs to be positive and less than 1.")
1589             else:
1590                rtol_min=self.__rtol_min
1591    
1592             if not chi_max == None:
1593                if chi_max<=0 or chi_max>=1:
1594                   raise ValueError("chi_max needs to be positive and less than 1.")
1595             else:
1596                chi_max = self.__chi_max
1597    
1598             if not reduction_factor == None:
1599                if reduction_factor<=0 or reduction_factor>1:
1600                   raise ValueError("reduction_factor need to be between zero and one.")
1601             else:
1602                reduction_factor=self.__reduction_factor
1603    
1604             if not theta == None:
1605                if theta<=0 or theta>1:
1606                   raise ValueError("theta need to be between zero and one.")
1607             else:
1608                theta=self.__theta
1609    
1610             if rtol_min>=rtol_max:
1611                 raise ValueError("rtol_max = %e needs to be greater than rtol_min = %e"%(rtol_max,rtol_min))
1612             self.__chi_max = chi_max
1613             self.__rtol_max = rtol_max
1614             self.__K_p = K_p
1615             self.__K_v = K_v
1616             self.__reduction_factor = reduction_factor
1617             self.__theta = theta
1618             self.__rtol_min=rtol_min
1619    
1620          #=============================================================
1621          def inner_pBv(self,p,Bv):
1622             """
1623             Returns inner product of element p and Bv (overwrite).
1624    
1625             :param p: a pressure increment
1626             :param Bv: a residual
1627             :return: inner product of element p and Bv
1628             :rtype: ``float``
1629             :note: used if PCG is applied.
1630             """
1631             raise NotImplementedError("no inner product for p and Bv implemented.")
1632    
1633          def inner_p(self,p0,p1):
1634             """
1635             Returns inner product of p0 and p1 (overwrite).
1636    
1637             :param p0: a pressure
1638             :param p1: a pressure
1639             :return: inner product of p0 and p1
1640             :rtype: ``float``
1641             """
1642             raise NotImplementedError("no inner product for p implemented.")
1643      
1644          def norm_v(self,v):
1645             """
1646             Returns the norm of v (overwrite).
1647    
1648             :param v: a velovity
1649             :return: norm of v
1650             :rtype: non-negative ``float``
1651             """
1652             raise NotImplementedError("no norm of v implemented.")
1653          def getDV(self, p, v, tol):
1654             """
1655             return a correction to the value for a given v and a given p with accuracy `tol` (overwrite)
1656    
1657             :param p: pressure
1658             :param v: pressure
1659             :return: dv given as *dv= A^{-1} (f-A v-B^*p)*
1660             :note: Only *A* may depend on *v* and *p*
1661             """
1662             raise NotImplementedError("no dv calculation implemented.")
1663    
1664            
1665          def Bv(self,v, tol):
1666            """
1667            Returns Bv with accuracy `tol` (overwrite)
1668    
1669            :rtype: equal to the type of p
1670            :note: boundary conditions on p should be zero!
1671            """
1672            raise NotImplementedError("no operator B implemented.")
1673    
1674          def norm_Bv(self,Bv):
1675            """
1676            Returns the norm of Bv (overwrite).
1677    
1678            :rtype: equal to the type of p
1679            :note: boundary conditions on p should be zero!
1680            """
1681            raise NotImplementedError("no norm of Bv implemented.")
1682    
1683          def solve_AinvBt(self,dp, tol):
1684             """
1685             Solves *A dv=B^*dp* with accuracy `tol`
1686    
1687             :param dp: a pressure increment
1688             :return: the solution of *A dv=B^*dp*
1689             :note: boundary conditions on dv should be zero! *A* is the operator used in ``getDV`` and must not be altered.
1690             """
1691             raise NotImplementedError("no operator A implemented.")
1692    
1693          def solve_prec(self,Bv, tol):
1694             """
1695             Provides a preconditioner for *(BA^{-1}B^ * )* applied to Bv with accuracy `tol`
1696    
1697             :rtype: equal to the type of p
1698             :note: boundary conditions on p should be zero!
1699             """
1700             raise NotImplementedError("no preconditioner for Schur complement implemented.")
1701          #=============================================================
1702          def __Aprod_PCG(self,dp):
1703              dv=self.solve_AinvBt(dp, self.__subtol)
1704              return ArithmeticTuple(dv,self.Bv(dv, self.__subtol))
1705    
1706          def __inner_PCG(self,p,r):
1707             return self.inner_pBv(p,r[1])
1708    
1709          def __Msolve_PCG(self,r):
1710              return self.solve_prec(r[1], self.__subtol)
1711          #=============================================================
1712          def __Aprod_GMRES(self,p):
1713              return self.solve_prec(self.Bv(self.solve_AinvBt(p, self.__subtol), self.__subtol), self.__subtol)
1714          def __inner_GMRES(self,p0,p1):
1715             return self.inner_p(p0,p1)
1716    
1717          #=============================================================
1718          def norm_p(self,p):
1719              """
1720              calculates the norm of ``p``
1721              
1722              :param p: a pressure
1723              :return: the norm of ``p`` using the inner product for pressure
1724              :rtype: ``float``
1725              """
1726              f=self.inner_p(p,p)
1727              if f<0: raise ValueError("negative pressure norm.")
1728              return math.sqrt(f)
1729          
1730          def solve(self,v,p,max_iter=20, verbose=False, usePCG=True, iter_restart=20, max_correction_steps=10):
1731             """
1732             Solves the saddle point problem using initial guesses v and p.
1733    
1734             :param v: initial guess for velocity
1735             :param p: initial guess for pressure
1736             :type v: `Data`
1737             :type p: `Data`
1738             :param usePCG: indicates the usage of the PCG rather than GMRES scheme.
1739             :param max_iter: maximum number of iteration steps per correction
1740                              attempt
1741             :param verbose: if True, shows information on the progress of the
1742                             saddlepoint problem solver.
1743             :param iter_restart: restart the iteration after ``iter_restart`` steps
1744                                  (only used if useUzaw=False)
1745             :type usePCG: ``bool``
1746             :type max_iter: ``int``
1747             :type verbose: ``bool``
1748             :type iter_restart: ``int``
1749             :rtype: ``tuple`` of `Data` objects
1750             :note: typically this method is overwritten by a subclass. It provides a wrapper for the ``_solve`` method.
1751             """
1752             return self._solve(v=v,p=p,max_iter=max_iter,verbose=verbose, usePCG=usePCG, iter_restart=iter_restart, max_correction_steps=max_correction_steps)
1753    
1754          def _solve(self,v,p,max_iter=20, verbose=False, usePCG=True, iter_restart=20, max_correction_steps=10):
1755             """
1756             see `_solve` method.
1757             """
1758             self.verbose=verbose
1759             rtol=self.getTolerance()
1760             atol=self.getAbsoluteTolerance()
1761    
1762             K_p=self.__K_p
1763             K_v=self.__K_v
1764             correction_step=0
1765             converged=False
1766             chi=None
1767             eps=None
1768    
1769             if self.verbose: print(("HomogeneousSaddlePointProblem: start iteration: rtol= %e, atol=%e"%(rtol, atol)))
1770             while not converged:
1771    
1772                 # get tolerance for velecity increment:
1773                 if chi == None:
1774                    rtol_v=self.__rtol_max
1775                 else:
1776                    rtol_v=min(chi/K_v,self.__rtol_max)
1777                 rtol_v=max(rtol_v, self.__rtol_min)
1778                 if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: rtol_v= %e"%(correction_step,rtol_v)))
1779                 # get velocity increment:
1780                 dv1=self.getDV(p,v,rtol_v)
1781                 v1=v+dv1
1782                 Bv1=self.Bv(v1, rtol_v)
1783                 norm_Bv1=self.norm_Bv(Bv1)
1784                 norm_dv1=self.norm_v(dv1)
1785                 if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: norm_Bv1 = %e, norm_dv1 = %e"%(correction_step, norm_Bv1, norm_dv1)))
1786                 if norm_dv1*self.__theta < norm_Bv1:
1787                    # get tolerance for pressure increment:
1788                    large_Bv1=True
1789                    if chi == None or eps == None:
1790                       rtol_p=self.__rtol_max
1791                    else:
1792                       rtol_p=min(chi**2*eps/K_p/norm_Bv1, self.__rtol_max)
1793                    self.__subtol=max(rtol_p**2, self.__rtol_min)
1794                    if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: rtol_p= %e"%(correction_step,rtol_p)))
1795                    # now we solve for the pressure increment dp from B*A^{-1}B^* dp = Bv1
1796                    if usePCG:
1797                        dp,r,a_norm=PCG(ArithmeticTuple(v1,Bv1),self.__Aprod_PCG,0*p,self.__Msolve_PCG,self.__inner_PCG,atol=0, rtol=rtol_p,iter_max=max_iter, verbose=self.verbose)
1798                        v2=r[0]
1799                        Bv2=r[1]
1800                    else:
1801                        # don't use!!!!
1802                        dp=GMRES(self.solve_prec(Bv1,self.__subtol),self.__Aprod_GMRES, 0*p, self.__inner_GMRES,atol=0, rtol=rtol_p,iter_max=max_iter, iter_restart=iter_restart, verbose=self.verbose)
1803                        dv2=self.solve_AinvBt(dp, self.__subtol)
1804                        v2=v1-dv2
1805                        Bv2=self.Bv(v2, self.__subtol)
1806                    p2=p+dp
1807                 else:
1808                    large_Bv1=False
1809                    v2=v1
1810                    p2=p
1811                 # update business:
1812                 norm_dv2=self.norm_v(v2-v)
1813                 norm_v2=self.norm_v(v2)
1814                 if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: v2 = %e, norm_dv2 = %e"%(correction_step, norm_v2, self.norm_v(v2-v))))
1815                 eps, eps_old = max(norm_Bv1, norm_dv2), eps
1816                 if eps_old == None:
1817                      chi, chi_old = None, chi
1818                 else:
1819                      chi, chi_old = min(eps/ eps_old, self.__chi_max), chi
1820                 if eps != None:
1821                     if chi !=None:
1822                        if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: convergence rate = %e, correction = %e"%(correction_step,chi, eps)))
1823                     else:
1824                        if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: correction = %e"%(correction_step, eps)))
1825                 if eps <= rtol*norm_v2+atol :
1826                     converged = True
1827                 else:
1828                     if correction_step>=max_correction_steps:
1829                          raise CorrectionFailed("Given up after %d correction steps."%correction_step)
1830                     if chi_old!=None:
1831                        K_p=max(1,self.__reduction_factor*K_p,(chi-chi_old)/chi_old**2*K_p)
1832                        K_v=max(1,self.__reduction_factor*K_v,(chi-chi_old)/chi_old**2*K_p)
1833                        if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: new adjustment factor K = %e"%(correction_step,K_p)))
1834                 correction_step+=1
1835                 v,p =v2, p2
1836             if self.verbose: print(("HomogeneousSaddlePointProblem: tolerance reached after %s steps."%correction_step))
1837             return v,p
1838          #========================================================================
1839          def setTolerance(self,tolerance=1.e-4):
1840             """
1841             Sets the relative tolerance for (v,p).
1842    
1843             :param tolerance: tolerance to be used
1844             :type tolerance: non-negative ``float``
1845             """
1846             if tolerance<0:
1847                 raise ValueError("tolerance must be positive.")
1848             self.__rtol=tolerance
1849    
1850          def getTolerance(self):
1851             """
1852             Returns the relative tolerance.
1853    
1854             :return: relative tolerance
1855             :rtype: ``float``
1856             """
1857             return self.__rtol
1858    
1859          def setAbsoluteTolerance(self,tolerance=0.):
1860             """
1861             Sets the absolute tolerance.
1862    
1863             :param tolerance: tolerance to be used
1864             :type tolerance: non-negative ``float``
1865             """
1866             if tolerance<0:
1867                 raise ValueError("tolerance must be non-negative.")
1868             self.__atol=tolerance
1869    
1870          def getAbsoluteTolerance(self):
1871             """
1872             Returns the absolute tolerance.
1873    
1874             :return: absolute tolerance
1875             :rtype: ``float``
1876             """
1877             return self.__atol
1878    
1879    
1880    def MaskFromBoundaryTag(domain,*tags):
1881       """
1882       Creates a mask on the Solution(domain) function space where the value is
1883       one for samples that touch the boundary tagged by tags.
1884    
1885       Usage: m=MaskFromBoundaryTag(domain, "left", "right")
1886    
1887       :param domain: domain to be used
1888       :type domain: `escript.Domain`
1889       :param tags: boundary tags
1890       :type tags: ``str``
1891       :return: a mask which marks samples that are touching the boundary tagged
1892                by any of the given tags
1893       :rtype: `escript.Data` of rank 0
1894       """
1895       pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1896       d=escript.Scalar(0.,escript.FunctionOnBoundary(domain))
1897       for t in tags: d.setTaggedValue(t,1.)
1898       pde.setValue(y=d)
1899       return util.whereNonZero(pde.getRightHandSide())
1900    
1901    def MaskFromTag(domain,*tags):
1902       """
1903       Creates a mask on the Solution(domain) function space where the value is
1904       one for samples that touch regions tagged by tags.
1905    
1906       Usage: m=MaskFromTag(domain, "ham")
1907    
1908       :param domain: domain to be used
1909       :type domain: `escript.Domain`
1910       :param tags: boundary tags
1911       :type tags: ``str``
1912       :return: a mask which marks samples that are touching the boundary tagged
1913                by any of the given tags
1914       :rtype: `escript.Data` of rank 0
1915       """
1916       pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1917       d=escript.Scalar(0.,escript.Function(domain))
1918       for t in tags: d.setTaggedValue(t,1.)
1919       pde.setValue(Y=d)
1920       return util.whereNonZero(pde.getRightHandSide())
1921    
1922    
 # vim: expandtab shiftwidth=4:  

Legend:
Removed from v.351  
changed lines
  Added in v.3928

  ViewVC Help
Powered by ViewVC 1.1.26