/[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 3852 by jfenwick, Thu Mar 1 05:34:52 2012 UTC
# Line 1  Line 1 
1  # $Id$  
2    ########################################################
3    #
4    # Copyright (c) 2003-2010 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-2010 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             else:             r=data.getRank()
416               out=data.interpolate(self.getFunctionSpace()).convertToNumArrayFromDPNo(self.getId()[0],self.getId()[1])             if isinstance(id,list):
417             if data.getRank()==0:                 out=[]
418                return out[0]                 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                return out               out=numpy.array(dat.getTupleForGlobalDataPoint(*id))
427                 if data.getRank()==0:
428                    return out[0]
429                 else:
430                    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    #       data.expand()   # Need to ensure that this is done globally
439    #       if isinstance(data, escript.Data):
440    #    id=self.getId()
441    #    if isinstance(id, list):
442    #      for i in id:
443    #         data._setTupleForGlobalDataPoint(i[1], i[0], v)
444    #    else:
445    #      data._setTupleForGlobalDataPoint(id[1], id[0], v)
446    #       else:
447    #    raise TypeError, "setValue: Invalid argument type."
448    
449    
450    def getInfLocator(arg):
451        """
452        Return a Locator for a point with the inf value over all arg.
453        """
454        if not isinstance(arg, escript.Data):
455           raise TypeError("getInfLocator: Unknown argument type.")
456        a_inf=util.inf(arg)
457        loc=util.length(arg-a_inf).minGlobalDataPoint() # This gives us the location but not coords
458        x=arg.getFunctionSpace().getX()
459        x_min=x.getTupleForGlobalDataPoint(*loc)
460        return Locator(arg.getFunctionSpace(),x_min)
461    
462    def getSupLocator(arg):
463        """
464        Return a Locator for a point with the sup value over all arg.
465        """
466        if not isinstance(arg, escript.Data):
467           raise TypeError("getInfLocator: Unknown argument type.")
468        a_inf=util.sup(arg)
469        loc=util.length(arg-a_inf).minGlobalDataPoint() # This gives us the location but not coords
470        x=arg.getFunctionSpace().getX()
471        x_min=x.getTupleForGlobalDataPoint(*loc)
472        return Locator(arg.getFunctionSpace(),x_min)
473        
474    
475    class SolverSchemeException(Exception):
476       """
477       This is a generic exception thrown by solvers.
478       """
479       pass
480    
481    class IndefinitePreconditioner(SolverSchemeException):
482       """
483       Exception thrown if the preconditioner is not positive definite.
484       """
485       pass
486    
487    class MaxIterReached(SolverSchemeException):
488       """
489       Exception thrown if the maximum number of iteration steps is reached.
490       """
491       pass
492    
493    class CorrectionFailed(SolverSchemeException):
494       """
495       Exception thrown if no convergence has been achieved in the solution
496       correction scheme.
497       """
498       pass
499    
500    class IterationBreakDown(SolverSchemeException):
501       """
502       Exception thrown if the iteration scheme encountered an incurable breakdown.
503       """
504       pass
505    
506    class NegativeNorm(SolverSchemeException):
507       """
508       Exception thrown if a norm calculation returns a negative norm.
509       """
510       pass
511    
512    def PCG(r, Aprod, x, Msolve, bilinearform, atol=0, rtol=1.e-8, iter_max=100, initial_guess=True, verbose=False):
513       """
514       Solver for
515    
516       *Ax=b*
517    
518       with a symmetric and positive definite operator A (more details required!).
519       It uses the conjugate gradient method with preconditioner M providing an
520       approximation of A.
521    
522       The iteration is terminated if
523    
524       *|r| <= atol+rtol*|r0|*
525    
526       where *r0* is the initial residual and *|.|* is the energy norm. In fact
527    
528       *|r| = sqrt( bilinearform(Msolve(r),r))*
529    
530       For details on the preconditioned conjugate gradient method see the book:
531    
532       I{Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
533       T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
534       C. Romine, and H. van der Vorst}.
535    
536       :param r: initial residual *r=b-Ax*. ``r`` is altered.
537       :type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
538       :param x: an initial guess for the solution
539       :type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
540       :param Aprod: returns the value Ax
541       :type Aprod: function ``Aprod(x)`` where ``x`` is of the same object like
542                    argument ``x``. The returned object needs to be of the same type
543                    like argument ``r``.
544       :param Msolve: solves Mx=r
545       :type Msolve: function ``Msolve(r)`` where ``r`` is of the same type like
546                     argument ``r``. The returned object needs to be of the same
547                     type like argument ``x``.
548       :param bilinearform: inner product ``<x,r>``
549       :type bilinearform: function ``bilinearform(x,r)`` where ``x`` is of the same
550                           type like argument ``x`` and ``r`` is. The returned value
551                           is a ``float``.
552       :param atol: absolute tolerance
553       :type atol: non-negative ``float``
554       :param rtol: relative tolerance
555       :type rtol: non-negative ``float``
556       :param iter_max: maximum number of iteration steps
557       :type iter_max: ``int``
558       :return: the solution approximation and the corresponding residual
559       :rtype: ``tuple``
560       :warning: ``r`` and ``x`` are altered.
561       """
562       iter=0
563       rhat=Msolve(r)
564       d = rhat
565       rhat_dot_r = bilinearform(rhat, r)
566       if rhat_dot_r<0: raise NegativeNorm("negative norm.")
567       norm_r0=math.sqrt(rhat_dot_r)
568       atol2=atol+rtol*norm_r0
569       if atol2<=0:
570          raise ValueError("Non-positive tolarance.")
571       atol2=max(atol2, 100. * util.EPSILON * norm_r0)
572    
573       if verbose: print(("PCG: initial residual norm = %e (absolute tolerance = %e)"%(norm_r0, atol2)))
574    
575    
576       while not math.sqrt(rhat_dot_r) <= atol2:
577           iter+=1
578           if iter  >= iter_max: raise MaxIterReached("maximum number of %s steps reached."%iter_max)
579    
580           q=Aprod(d)
581           alpha = rhat_dot_r / bilinearform(d, q)
582           x += alpha * d
583           if isinstance(q,ArithmeticTuple):
584              r += q * (-alpha)      # Doing it the other way calls the float64.__mul__ not AT.__rmul__
585           else:
586               r += (-alpha) * q
587           rhat=Msolve(r)
588           rhat_dot_r_new = bilinearform(rhat, r)
589           beta = rhat_dot_r_new / rhat_dot_r
590           rhat+=beta * d
591           d=rhat
592    
593           rhat_dot_r = rhat_dot_r_new
594           if rhat_dot_r<0: raise NegativeNorm("negative norm.")
595           if verbose: print(("PCG: iteration step %s: residual norm = %e"%(iter, math.sqrt(rhat_dot_r))))
596       if verbose: print(("PCG: tolerance reached after %s steps."%iter))
597       return x,r,math.sqrt(rhat_dot_r)
598    
599    class Defect(object):
600        """
601        Defines a non-linear defect F(x) of a variable x.
602        """
603        def __init__(self):
604            """
605            Initializes defect.
606            """
607            self.setDerivativeIncrementLength()
608    
609        def bilinearform(self, x0, x1):
610            """
611            Returns the inner product of x0 and x1
612    
613            :param x0: value for x0
614            :param x1: value for x1
615            :return: the inner product of x0 and x1
616            :rtype: ``float``
617            """
618            return 0
619    
620        def norm(self,x):
621            """
622            Returns the norm of argument ``x``.
623    
624            :param x: a value
625            :return: norm of argument x
626            :rtype: ``float``
627            :note: by default ``sqrt(self.bilinearform(x,x)`` is returned.
628            """
629            s=self.bilinearform(x,x)
630            if s<0: raise NegativeNorm("negative norm.")
631            return math.sqrt(s)
632    
633        def eval(self,x):
634            """
635            Returns the value F of a given ``x``.
636    
637            :param x: value for which the defect ``F`` is evaluated
638            :return: value of the defect at ``x``
639            """
640            return 0
641    
642        def __call__(self,x):
643            return self.eval(x)
644    
645        def setDerivativeIncrementLength(self,inc=1000.*math.sqrt(util.EPSILON)):
646            """
647            Sets the relative length of the increment used to approximate the
648            derivative of the defect. The increment is inc*norm(x)/norm(v)*v in the
649            direction of v with x as a starting point.
650    
651            :param inc: relative increment length
652            :type inc: positive ``float``
653            """
654            if inc<=0: raise ValueError("positive increment required.")
655            self.__inc=inc
656    
657        def getDerivativeIncrementLength(self):
658            """
659            Returns the relative increment length used to approximate the
660            derivative of the defect.
661            :return: value of the defect at ``x``
662            :rtype: positive ``float``
663            """
664            return self.__inc
665    
666        def derivative(self, F0, x0, v, v_is_normalised=True):
667            """
668            Returns the directional derivative at ``x0`` in the direction of ``v``.
669    
670            :param F0: value of this defect at x0
671            :param x0: value at which derivative is calculated
672            :param v: direction
673            :param v_is_normalised: True to indicate that ``v`` is nomalized
674                                    (self.norm(v)=0)
675            :return: derivative of this defect at x0 in the direction of ``v``
676            :note: by default numerical evaluation (self.eval(x0+eps*v)-F0)/eps is
677                   used but this method maybe overwritten to use exact evaluation.
678            """
679            normx=self.norm(x0)
680            if normx>0:
681                 epsnew = self.getDerivativeIncrementLength() * normx
682            else:
683                 epsnew = self.getDerivativeIncrementLength()
684            if not v_is_normalised:
685                normv=self.norm(v)
686                if normv<=0:
687                   return F0*0
688                else:
689                   epsnew /= normv
690            F1=self.eval(x0 + epsnew * v)
691            return (F1-F0)/epsnew
692    
693    ######################################
694    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):
695       """
696       Solves a non-linear problem *F(x)=0* for unknown *x* using the stopping
697       criterion:
698    
699       *norm(F(x) <= atol + rtol * norm(F(x0)*
700    
701       where *x0* is the initial guess.
702    
703       :param defect: object defining the function *F*. ``defect.norm`` defines the
704                      *norm* used in the stopping criterion.
705       :type defect: `Defect`
706       :param x: initial guess for the solution, ``x`` is altered.
707       :type x: any object type allowing basic operations such as
708                ``numpy.ndarray``, `Data`
709       :param iter_max: maximum number of iteration steps
710       :type iter_max: positive ``int``
711       :param sub_iter_max: maximum number of inner iteration steps
712       :type sub_iter_max: positive ``int``
713       :param atol: absolute tolerance for the solution
714       :type atol: positive ``float``
715       :param rtol: relative tolerance for the solution
716       :type rtol: positive ``float``
717       :param gamma: tolerance safety factor for inner iteration
718       :type gamma: positive ``float``, less than 1
719       :param subtol_max: upper bound for inner tolerance
720       :type subtol_max: positive ``float``, less than 1
721       :return: an approximation of the solution with the desired accuracy
722       :rtype: same type as the initial guess
723       """
724       lmaxit=iter_max
725       if atol<0: raise ValueError("atol needs to be non-negative.")
726       if rtol<0: raise ValueError("rtol needs to be non-negative.")
727       if rtol+atol<=0: raise ValueError("rtol or atol needs to be non-negative.")
728       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)
729       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)
730    
731       F=defect(x)
732       fnrm=defect.norm(F)
733       stop_tol=atol + rtol*fnrm
734       subtol=subtol_max
735       if verbose: print(("NewtonGMRES: initial residual = %e."%fnrm))
736       if verbose: print(("             tolerance = %e."%subtol))
737       iter=1
738       #
739       # main iteration loop
740       #
741       while not fnrm<=stop_tol:
742                if iter  >= iter_max: raise MaxIterReached("maximum number of %s steps reached."%iter_max)
743                #
744            #   adjust subtol_
745            #
746                if iter > 1:
747                   rat=fnrm/fnrmo
748                   subtol_old=subtol
749                   subtol=gamma*rat**2
750                   if gamma*subtol_old**2 > .1: subtol=max(subtol,gamma*subtol_old**2)
751                   subtol=max(min(subtol,subtol_max), .5*stop_tol/fnrm)
752            #
753            # calculate newton increment xc
754                #     if iter_max in __FDGMRES is reached MaxIterReached is thrown
755                #     if iter_restart -1 is returned as sub_iter
756                #     if  atol is reached sub_iter returns the numer of steps performed to get there
757                #
758                #
759                if verbose: print(("             subiteration (GMRES) is called with relative tolerance %e."%subtol))
760                try:
761                   xc, sub_iter=__FDGMRES(F, defect, x, subtol*fnrm, iter_max=iter_max-iter, iter_restart=sub_iter_max)
762                except MaxIterReached:
763                   raise MaxIterReached("maximum number of %s steps reached."%iter_max)
764                if sub_iter<0:
765                   iter+=sub_iter_max
766                else:
767                   iter+=sub_iter
768                # ====
769                x+=xc
770                F=defect(x)
771                iter+=1
772                fnrmo, fnrm=fnrm, defect.norm(F)
773                if verbose: print(("             step %s: residual %e."%(iter,fnrm)))
774       if verbose: print(("NewtonGMRES: completed after %s steps."%iter))
775       return x
776    
777    def __givapp(c,s,vin):
778        """
779        Applies a sequence of Givens rotations (c,s) recursively to the vector
780        ``vin``
781    
782        :warning: ``vin`` is altered.
783        """
784        vrot=vin
785        if isinstance(c,float):
786            vrot=[c*vrot[0]-s*vrot[1],s*vrot[0]+c*vrot[1]]
787        else:
788            for i in range(len(c)):
789                w1=c[i]*vrot[i]-s[i]*vrot[i+1]
790                w2=s[i]*vrot[i]+c[i]*vrot[i+1]
791                vrot[i]=w1
792                vrot[i+1]=w2
793        return vrot
794    
795    def __FDGMRES(F0, defect, x0, atol, iter_max=100, iter_restart=20):
796       h=numpy.zeros((iter_restart,iter_restart),numpy.float64)
797       c=numpy.zeros(iter_restart,numpy.float64)
798       s=numpy.zeros(iter_restart,numpy.float64)
799       g=numpy.zeros(iter_restart,numpy.float64)
800       v=[]
801    
802       rho=defect.norm(F0)
803       if rho<=0.: return x0*0
804    
805       v.append(-F0/rho)
806       g[0]=rho
807       iter=0
808       while rho > atol and iter<iter_restart-1:
809            if iter  >= iter_max:
810                raise MaxIterReached("maximum number of %s steps reached."%iter_max)
811    
812            p=defect.derivative(F0,x0,v[iter], v_is_normalised=True)
813            v.append(p)
814    
815            v_norm1=defect.norm(v[iter+1])
816    
817            # Modified Gram-Schmidt
818            for j in range(iter+1):
819                h[j,iter]=defect.bilinearform(v[j],v[iter+1])
820                v[iter+1]-=h[j,iter]*v[j]
821    
822            h[iter+1,iter]=defect.norm(v[iter+1])
823            v_norm2=h[iter+1,iter]
824    
825            # Reorthogonalize if needed
826            if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
827                for j in range(iter+1):
828                    hr=defect.bilinearform(v[j],v[iter+1])
829                    h[j,iter]=h[j,iter]+hr
830                    v[iter+1] -= hr*v[j]
831    
832                v_norm2=defect.norm(v[iter+1])
833                h[iter+1,iter]=v_norm2
834            #   watch out for happy breakdown
835            if not v_norm2 == 0:
836                v[iter+1]=v[iter+1]/h[iter+1,iter]
837    
838            #   Form and store the information for the new Givens rotation
839            if iter > 0 :
840                hhat=numpy.zeros(iter+1,numpy.float64)
841                for i in range(iter+1) : hhat[i]=h[i,iter]
842                hhat=__givapp(c[0:iter],s[0:iter],hhat);
843                for i in range(iter+1) : h[i,iter]=hhat[i]
844    
845            mu=math.sqrt(h[iter,iter]*h[iter,iter]+h[iter+1,iter]*h[iter+1,iter])
846    
847            if mu!=0 :
848                c[iter]=h[iter,iter]/mu
849                s[iter]=-h[iter+1,iter]/mu
850                h[iter,iter]=c[iter]*h[iter,iter]-s[iter]*h[iter+1,iter]
851                h[iter+1,iter]=0.0
852                gg=__givapp(c[iter],s[iter],[g[iter],g[iter+1]])
853                g[iter]=gg[0]
854                g[iter+1]=gg[1]
855    
856            # Update the residual norm
857            rho=abs(g[iter+1])
858            iter+=1
859    
860       # At this point either iter > iter_max or rho < tol.
861       # It's time to compute x and leave.
862       if iter > 0 :
863         y=numpy.zeros(iter,numpy.float64)
864         y[iter-1] = g[iter-1] / h[iter-1,iter-1]
865         if iter > 1 :
866            i=iter-2
867            while i>=0 :
868              y[i] = ( g[i] - numpy.dot(h[i,i+1:iter], y[i+1:iter])) / h[i,i]
869              i=i-1
870         xhat=v[iter-1]*y[iter-1]
871         for i in range(iter-1):
872           xhat += v[i]*y[i]
873       else :
874          xhat=v[0] * 0
875    
876       if iter<iter_restart-1:
877          stopped=iter
878       else:
879          stopped=-1
880    
881       return xhat,stopped
882    
883    def GMRES(r, Aprod, x, bilinearform, atol=0, rtol=1.e-8, iter_max=100, iter_restart=20, verbose=False,P_R=None):
884       """
885       Solver for
886    
887       *Ax=b*
888    
889       with a general operator A (more details required!).
890       It uses the generalized minimum residual method (GMRES).
891    
892       The iteration is terminated if
893    
894       *|r| <= atol+rtol*|r0|*
895    
896       where *r0* is the initial residual and *|.|* is the energy norm. In fact
897    
898       *|r| = sqrt( bilinearform(r,r))*
899    
900       :param r: initial residual *r=b-Ax*. ``r`` is altered.
901       :type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
902       :param x: an initial guess for the solution
903       :type x: same like ``r``
904       :param Aprod: returns the value Ax
905       :type Aprod: function ``Aprod(x)`` where ``x`` is of the same object like
906                    argument ``x``. The returned object needs to be of the same
907                    type like argument ``r``.
908       :param bilinearform: inner product ``<x,r>``
909       :type bilinearform: function ``bilinearform(x,r)`` where ``x`` is of the same
910                           type like argument ``x`` and ``r``. The returned value is
911                           a ``float``.
912       :param atol: absolute tolerance
913       :type atol: non-negative ``float``
914       :param rtol: relative tolerance
915       :type rtol: non-negative ``float``
916       :param iter_max: maximum number of iteration steps
917       :type iter_max: ``int``
918       :param iter_restart: in order to save memory the orthogonalization process
919                            is terminated after ``iter_restart`` steps and the
920                            iteration is restarted.
921       :type iter_restart: ``int``
922       :return: the solution approximation and the corresponding residual
923       :rtype: ``tuple``
924       :warning: ``r`` and ``x`` are altered.
925       """
926       m=iter_restart
927       restarted=False
928       iter=0
929       if rtol>0:
930          r_dot_r = bilinearform(r, r)
931          if r_dot_r<0: raise NegativeNorm("negative norm.")
932          atol2=atol+rtol*math.sqrt(r_dot_r)
933          if verbose: print(("GMRES: norm of right hand side = %e (absolute tolerance = %e)"%(math.sqrt(r_dot_r), atol2)))
934       else:
935          atol2=atol
936          if verbose: print(("GMRES: absolute tolerance = %e"%atol2))
937       if atol2<=0:
938          raise ValueError("Non-positive tolarance.")
939    
940       while True:
941          if iter  >= iter_max: raise MaxIterReached("maximum number of %s steps reached"%iter_max)
942          if restarted:
943             r2 = r-Aprod(x-x2)
944          else:
945             r2=1*r
946          x2=x*1.
947          x,stopped=_GMRESm(r2, Aprod, x, bilinearform, atol2, iter_max=iter_max-iter, iter_restart=m, verbose=verbose,P_R=P_R)
948          iter+=iter_restart
949          if stopped: break
950          if verbose: print("GMRES: restart.")
951          restarted=True
952       if verbose: print("GMRES: tolerance has been reached.")
953       return x
954    
955    def _GMRESm(r, Aprod, x, bilinearform, atol, iter_max=100, iter_restart=20, verbose=False, P_R=None):
956       iter=0
957    
958       h=numpy.zeros((iter_restart+1,iter_restart),numpy.float64)
959       c=numpy.zeros(iter_restart,numpy.float64)
960       s=numpy.zeros(iter_restart,numpy.float64)
961       g=numpy.zeros(iter_restart+1,numpy.float64)
962       v=[]
963    
964       r_dot_r = bilinearform(r, r)
965       if r_dot_r<0: raise NegativeNorm("negative norm.")
966       rho=math.sqrt(r_dot_r)
967    
968       v.append(r/rho)
969       g[0]=rho
970    
971       if verbose: print(("GMRES: initial residual %e (absolute tolerance = %e)"%(rho,atol)))
972       while not (rho<=atol or iter==iter_restart):
973    
974            if iter  >= iter_max: raise MaxIterReached("maximum number of %s steps reached."%iter_max)
975    
976            if P_R!=None:
977                p=Aprod(P_R(v[iter]))
978            else:
979                p=Aprod(v[iter])
980            v.append(p)
981    
982            v_norm1=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
983    
984    # Modified Gram-Schmidt
985            for j in range(iter+1):
986              h[j,iter]=bilinearform(v[j],v[iter+1])
987              v[iter+1]-=h[j,iter]*v[j]
988    
989            h[iter+1,iter]=math.sqrt(bilinearform(v[iter+1],v[iter+1]))
990            v_norm2=h[iter+1,iter]
991    
992    # Reorthogonalize if needed
993            if v_norm1 + 0.001*v_norm2 == v_norm1:   #Brown/Hindmarsh condition (default)
994             for j in range(iter+1):
995                hr=bilinearform(v[j],v[iter+1])
996                h[j,iter]=h[j,iter]+hr
997                v[iter+1] -= hr*v[j]
998    
999             v_norm2=math.sqrt(bilinearform(v[iter+1], v[iter+1]))
1000             h[iter+1,iter]=v_norm2
1001    
1002    #   watch out for happy breakdown
1003            if not v_norm2 == 0:
1004             v[iter+1]=v[iter+1]/h[iter+1,iter]
1005    
1006    #   Form and store the information for the new Givens rotation
1007            if iter > 0: h[:iter+1,iter]=__givapp(c[:iter],s[:iter],h[:iter+1,iter])
1008            mu=math.sqrt(h[iter,iter]*h[iter,iter]+h[iter+1,iter]*h[iter+1,iter])
1009    
1010            if mu!=0 :
1011                    c[iter]=h[iter,iter]/mu
1012                    s[iter]=-h[iter+1,iter]/mu
1013                    h[iter,iter]=c[iter]*h[iter,iter]-s[iter]*h[iter+1,iter]
1014                    h[iter+1,iter]=0.0
1015                    gg=__givapp(c[iter],s[iter],[g[iter],g[iter+1]])
1016                    g[iter]=gg[0]
1017                    g[iter+1]=gg[1]
1018    # Update the residual norm
1019    
1020            rho=abs(g[iter+1])
1021            if verbose: print(("GMRES: iteration step %s: residual %e"%(iter,rho)))
1022            iter+=1
1023    
1024    # At this point either iter > iter_max or rho < tol.
1025    # It's time to compute x and leave.
1026    
1027       if verbose: print(("GMRES: iteration stopped after %s step."%iter))
1028       if iter > 0 :
1029         y=numpy.zeros(iter,numpy.float64)
1030         y[iter-1] = g[iter-1] / h[iter-1,iter-1]
1031         if iter > 1 :
1032            i=iter-2
1033            while i>=0 :
1034              y[i] = ( g[i] - numpy.dot(h[i,i+1:iter], y[i+1:iter])) / h[i,i]
1035              i=i-1
1036         xhat=v[iter-1]*y[iter-1]
1037         for i in range(iter-1):
1038           xhat += v[i]*y[i]
1039       else:
1040         xhat=v[0] * 0
1041       if P_R!=None:
1042          x += P_R(xhat)
1043       else:
1044          x += xhat
1045       if iter<iter_restart-1:
1046          stopped=True
1047       else:
1048          stopped=False
1049    
1050       return x,stopped
1051    
1052    def MINRES(r, Aprod, x, Msolve, bilinearform, atol=0, rtol=1.e-8, iter_max=100):
1053        """
1054        Solver for
1055    
1056        *Ax=b*
1057    
1058        with a symmetric and positive definite operator A (more details required!).
1059        It uses the minimum residual method (MINRES) with preconditioner M
1060        providing an approximation of A.
1061    
1062        The iteration is terminated if
1063    
1064        *|r| <= atol+rtol*|r0|*
1065    
1066        where *r0* is the initial residual and *|.|* is the energy norm. In fact
1067    
1068        *|r| = sqrt( bilinearform(Msolve(r),r))*
1069    
1070        For details on the preconditioned conjugate gradient method see the book:
1071    
1072        I{Templates for the Solution of Linear Systems by R. Barrett, M. Berry,
1073        T.F. Chan, J. Demmel, J. Donato, J. Dongarra, V. Eijkhout, R. Pozo,
1074        C. Romine, and H. van der Vorst}.
1075    
1076        :param r: initial residual *r=b-Ax*. ``r`` is altered.
1077        :type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1078        :param x: an initial guess for the solution
1079        :type x: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1080        :param Aprod: returns the value Ax
1081        :type Aprod: function ``Aprod(x)`` where ``x`` is of the same object like
1082                     argument ``x``. The returned object needs to be of the same
1083                     type like argument ``r``.
1084        :param Msolve: solves Mx=r
1085        :type Msolve: function ``Msolve(r)`` where ``r`` is of the same type like
1086                      argument ``r``. The returned object needs to be of the same
1087                      type like argument ``x``.
1088        :param bilinearform: inner product ``<x,r>``
1089        :type bilinearform: function ``bilinearform(x,r)`` where ``x`` is of the same
1090                            type like argument ``x`` and ``r`` is. The returned value
1091                            is a ``float``.
1092        :param atol: absolute tolerance
1093        :type atol: non-negative ``float``
1094        :param rtol: relative tolerance
1095        :type rtol: non-negative ``float``
1096        :param iter_max: maximum number of iteration steps
1097        :type iter_max: ``int``
1098        :return: the solution approximation and the corresponding residual
1099        :rtype: ``tuple``
1100        :warning: ``r`` and ``x`` are altered.
1101        """
1102        #------------------------------------------------------------------
1103        # Set up y and v for the first Lanczos vector v1.
1104        # y  =  beta1 P' v1,  where  P = C**(-1).
1105        # v is really P' v1.
1106        #------------------------------------------------------------------
1107        r1    = r
1108        y = Msolve(r)
1109        beta1 = bilinearform(y,r)
1110    
1111        if beta1< 0: raise NegativeNorm("negative norm.")
1112    
1113        #  If r = 0 exactly, stop with x
1114        if beta1==0: return x
1115    
1116        if beta1> 0: beta1  = math.sqrt(beta1)
1117    
1118        #------------------------------------------------------------------
1119        # Initialize quantities.
1120        # ------------------------------------------------------------------
1121        iter   = 0
1122        Anorm = 0
1123        ynorm = 0
1124        oldb   = 0
1125        beta   = beta1
1126        dbar   = 0
1127        epsln  = 0
1128        phibar = beta1
1129        rhs1   = beta1
1130        rhs2   = 0
1131        rnorm  = phibar
1132        tnorm2 = 0
1133        ynorm2 = 0
1134        cs     = -1
1135        sn     = 0
1136        w      = r*0.
1137        w2     = r*0.
1138        r2     = r1
1139        eps    = 0.0001
1140    
1141        #---------------------------------------------------------------------
1142        # Main iteration loop.
1143        # --------------------------------------------------------------------
1144        while not rnorm<=atol+rtol*Anorm*ynorm:    #  checks ||r|| < (||A|| ||x||) * TOL
1145    
1146            if iter  >= iter_max: raise MaxIterReached("maximum number of %s steps reached."%iter_max)
1147            iter    = iter  +  1
1148    
1149            #-----------------------------------------------------------------
1150            # Obtain quantities for the next Lanczos vector vk+1, k = 1, 2,...
1151            # The general iteration is similar to the case k = 1 with v0 = 0:
1152            #
1153            #   p1      = Operator * v1  -  beta1 * v0,
1154            #   alpha1  = v1'p1,
1155            #   q2      = p2  -  alpha1 * v1,
1156            #   beta2^2 = q2'q2,
1157            #   v2      = (1/beta2) q2.
1158            #
1159            # Again, y = betak P vk,  where  P = C**(-1).
1160            #-----------------------------------------------------------------
1161            s = 1/beta                 # Normalize previous vector (in y).
1162            v = s*y                    # v = vk if P = I
1163    
1164            y      = Aprod(v)
1165    
1166            if iter >= 2:
1167              y = y - (beta/oldb)*r1
1168    
1169            alfa   = bilinearform(v,y)              # alphak
1170            y      += (- alfa/beta)*r2
1171            r1     = r2
1172            r2     = y
1173            y = Msolve(r2)
1174            oldb   = beta                           # oldb = betak
1175            beta   = bilinearform(y,r2)             # beta = betak+1^2
1176            if beta < 0: raise NegativeNorm("negative norm.")
1177    
1178            beta   = math.sqrt( beta )
1179            tnorm2 = tnorm2 + alfa*alfa + oldb*oldb + beta*beta
1180    
1181            if iter==1:                 # Initialize a few things.
1182              gmax   = abs( alfa )      # alpha1
1183              gmin   = gmax             # alpha1
1184    
1185            # Apply previous rotation Qk-1 to get
1186            #   [deltak epslnk+1] = [cs  sn][dbark    0   ]
1187            #   [gbar k dbar k+1]   [sn -cs][alfak betak+1].
1188    
1189            oldeps = epsln
1190            delta  = cs * dbar  +  sn * alfa  # delta1 = 0         deltak
1191            gbar   = sn * dbar  -  cs * alfa  # gbar 1 = alfa1     gbar k
1192            epsln  =               sn * beta  # epsln2 = 0         epslnk+1
1193            dbar   =            -  cs * beta  # dbar 2 = beta2     dbar k+1
1194    
1195            # Compute the next plane rotation Qk
1196    
1197            gamma  = math.sqrt(gbar*gbar+beta*beta)  # gammak
1198            gamma  = max(gamma,eps)
1199            cs     = gbar / gamma             # ck
1200            sn     = beta / gamma             # sk
1201            phi    = cs * phibar              # phik
1202            phibar = sn * phibar              # phibark+1
1203    
1204            # Update  x.
1205    
1206            denom = 1/gamma
1207            w1    = w2
1208            w2    = w
1209            w     = (v - oldeps*w1 - delta*w2) * denom
1210            x     +=  phi*w
1211    
1212            # Go round again.
1213    
1214            gmax   = max(gmax,gamma)
1215            gmin   = min(gmin,gamma)
1216            z      = rhs1 / gamma
1217            ynorm2 = z*z  +  ynorm2
1218            rhs1   = rhs2 -  delta*z
1219            rhs2   =      -  epsln*z
1220    
1221            # Estimate various norms and test for convergence.
1222    
1223            Anorm  = math.sqrt( tnorm2 )
1224            ynorm  = math.sqrt( ynorm2 )
1225    
1226            rnorm  = phibar
1227    
1228        return x
1229    
1230    def TFQMR(r, Aprod, x, bilinearform, atol=0, rtol=1.e-8, iter_max=100):
1231      """
1232      Solver for
1233    
1234      *Ax=b*
1235    
1236      with a general operator A (more details required!).
1237      It uses the Transpose-Free Quasi-Minimal Residual method (TFQMR).
1238    
1239      The iteration is terminated if
1240    
1241      *|r| <= atol+rtol*|r0|*
1242    
1243      where *r0* is the initial residual and *|.|* is the energy norm. In fact
1244    
1245      *|r| = sqrt( bilinearform(r,r))*
1246    
1247      :param r: initial residual *r=b-Ax*. ``r`` is altered.
1248      :type r: any object supporting inplace add (x+=y) and scaling (x=scalar*y)
1249      :param x: an initial guess for the solution
1250      :type x: same like ``r``
1251      :param Aprod: returns the value Ax
1252      :type Aprod: function ``Aprod(x)`` where ``x`` is of the same object like
1253                   argument ``x``. The returned object needs to be of the same type
1254                   like argument ``r``.
1255      :param bilinearform: inner product ``<x,r>``
1256      :type bilinearform: function ``bilinearform(x,r)`` where ``x`` is of the same
1257                          type like argument ``x`` and ``r``. The returned value is
1258                          a ``float``.
1259      :param atol: absolute tolerance
1260      :type atol: non-negative ``float``
1261      :param rtol: relative tolerance
1262      :type rtol: non-negative ``float``
1263      :param iter_max: maximum number of iteration steps
1264      :type iter_max: ``int``
1265      :rtype: ``tuple``
1266      :warning: ``r`` and ``x`` are altered.
1267      """
1268      u1=0
1269      u2=0
1270      y1=0
1271      y2=0
1272    
1273      w = r
1274      y1 = r
1275      iter = 0
1276      d = 0
1277      v = Aprod(y1)
1278      u1 = v
1279    
1280      theta = 0.0;
1281      eta = 0.0;
1282      rho=bilinearform(r,r)
1283      if rho < 0: raise NegativeNorm("negative norm.")
1284      tau = math.sqrt(rho)
1285      norm_r0=tau
1286      while tau>atol+rtol*norm_r0:
1287        if iter  >= iter_max: raise MaxIterReached("maximum number of %s steps reached."%iter_max)
1288    
1289        sigma = bilinearform(r,v)
1290        if sigma == 0.0: raise IterationBreakDown('TFQMR breakdown, sigma=0')
1291    
1292        alpha = rho / sigma
1293    
1294        for j in range(2):
1295    #
1296    #   Compute y2 and u2 only if you have to
1297    #
1298          if ( j == 1 ):
1299            y2 = y1 - alpha * v
1300            u2 = Aprod(y2)
1301    
1302          m = 2 * (iter+1) - 2 + (j+1)
1303          if j==0:
1304             w = w - alpha * u1
1305             d = y1 + ( theta * theta * eta / alpha ) * d
1306          if j==1:
1307             w = w - alpha * u2
1308             d = y2 + ( theta * theta * eta / alpha ) * d
1309    
1310          theta = math.sqrt(bilinearform(w,w))/ tau
1311          c = 1.0 / math.sqrt ( 1.0 + theta * theta )
1312          tau = tau * theta * c
1313          eta = c * c * alpha
1314          x = x + eta * d
1315    #
1316    #  Try to terminate the iteration at each pass through the loop
1317    #
1318        if rho == 0.0: raise IterationBreakDown('TFQMR breakdown, rho=0')
1319    
1320        rhon = bilinearform(r,w)
1321        beta = rhon / rho;
1322        rho = rhon;
1323        y1 = w + beta * y2;
1324        u1 = Aprod(y1)
1325        v = u1 + beta * ( u2 + beta * v )
1326    
1327        iter += 1
1328    
1329      return x
1330    
1331    
1332    #############################################
1333    
1334    class ArithmeticTuple(object):
1335       """
1336       Tuple supporting inplace update x+=y and scaling x=a*y where ``x,y`` is an
1337       ArithmeticTuple and ``a`` is a float.
1338    
1339       Example of usage::
1340    
1341           from esys.escript import Data
1342           from numpy import array
1343           a=Data(...)
1344           b=array([1.,4.])
1345           x=ArithmeticTuple(a,b)
1346           y=5.*x
1347    
1348       """
1349       def __init__(self,*args):
1350           """
1351           Initializes object with elements ``args``.
1352    
1353           :param args: tuple of objects that support inplace add (x+=y) and
1354                        scaling (x=a*y)
1355           """
1356           self.__items=list(args)
1357    
1358       def __len__(self):
1359           """
1360           Returns the number of items.
1361    
1362           :return: number of items
1363           :rtype: ``int``
1364           """
1365           return len(self.__items)
1366    
1367       def __getitem__(self,index):
1368           """
1369           Returns item at specified position.
1370    
1371           :param index: index of item to be returned
1372           :type index: ``int``
1373           :return: item with index ``index``
1374           """
1375           return self.__items.__getitem__(index)
1376    
1377       def __mul__(self,other):
1378           """
1379           Scales by ``other`` from the right.
1380    
1381           :param other: scaling factor
1382           :type other: ``float``
1383           :return: itemwise self*other
1384           :rtype: `ArithmeticTuple`
1385           """
1386           out=[]
1387           try:
1388               l=len(other)
1389               if l!=len(self):
1390                   raise ValueError("length of arguments don't match.")
1391               for i in range(l): out.append(self[i]*other[i])
1392           except TypeError:
1393               for i in range(len(self)): out.append(self[i]*other)
1394           return ArithmeticTuple(*tuple(out))
1395    
1396       def __rmul__(self,other):
1397           """
1398           Scales by ``other`` from the left.
1399    
1400           :param other: scaling factor
1401           :type other: ``float``
1402           :return: itemwise other*self
1403           :rtype: `ArithmeticTuple`
1404           """
1405           out=[]
1406           try:
1407               l=len(other)
1408               if l!=len(self):
1409                   raise ValueError("length of arguments don't match.")
1410               for i in range(l): out.append(other[i]*self[i])
1411           except TypeError:
1412               for i in range(len(self)): out.append(other*self[i])
1413           return ArithmeticTuple(*tuple(out))
1414    
1415       def __div__(self,other):
1416           """
1417           Scales by (1/``other``) from the right.
1418    
1419           :param other: scaling factor
1420           :type other: ``float``
1421           :return: itemwise self/other
1422           :rtype: `ArithmeticTuple`
1423           """
1424           return self*(1/other)
1425    
1426       def __rdiv__(self,other):
1427           """
1428           Scales by (1/``other``) from the left.
1429    
1430           :param other: scaling factor
1431           :type other: ``float``
1432           :return: itemwise other/self
1433           :rtype: `ArithmeticTuple`
1434           """
1435           out=[]
1436           try:
1437               l=len(other)
1438               if l!=len(self):
1439                   raise ValueError("length of arguments don't match.")
1440               for i in range(l): out.append(other[i]/self[i])
1441           except TypeError:
1442               for i in range(len(self)): out.append(other/self[i])
1443           return ArithmeticTuple(*tuple(out))
1444    
1445       def __iadd__(self,other):
1446           """
1447           Inplace addition of ``other`` to self.
1448    
1449           :param other: increment
1450           :type other: ``ArithmeticTuple``
1451           """
1452           if len(self) != len(other):
1453               raise ValueError("tuple lengths must match.")
1454           for i in range(len(self)):
1455               self.__items[i]+=other[i]
1456           return self
1457    
1458       def __add__(self,other):
1459           """
1460           Adds ``other`` to self.
1461    
1462           :param other: increment
1463           :type other: ``ArithmeticTuple``
1464           """
1465           out=[]
1466           try:
1467               l=len(other)
1468               if l!=len(self):
1469                   raise ValueError("length of arguments don't match.")
1470               for i in range(l): out.append(self[i]+other[i])
1471           except TypeError:
1472               for i in range(len(self)): out.append(self[i]+other)
1473           return ArithmeticTuple(*tuple(out))
1474    
1475       def __sub__(self,other):
1476           """
1477           Subtracts ``other`` from self.
1478    
1479           :param other: decrement
1480           :type other: ``ArithmeticTuple``
1481           """
1482           out=[]
1483           try:
1484               l=len(other)
1485               if l!=len(self):
1486                   raise ValueError("length of arguments don't match.")
1487               for i in range(l): out.append(self[i]-other[i])
1488           except TypeError:
1489               for i in range(len(self)): out.append(self[i]-other)
1490           return ArithmeticTuple(*tuple(out))
1491    
1492       def __isub__(self,other):
1493           """
1494           Inplace subtraction of ``other`` from self.
1495    
1496           :param other: decrement
1497           :type other: ``ArithmeticTuple``
1498           """
1499           if len(self) != len(other):
1500               raise ValueError("tuple length must match.")
1501           for i in range(len(self)):
1502               self.__items[i]-=other[i]
1503           return self
1504    
1505       def __neg__(self):
1506           """
1507           Negates values.
1508           """
1509           out=[]
1510           for i in range(len(self)):
1511               out.append(-self[i])
1512           return ArithmeticTuple(*tuple(out))
1513    
1514    
1515    class HomogeneousSaddlePointProblem(object):
1516          """
1517          This class provides a framework for solving linear homogeneous saddle
1518          point problems of the form::
1519    
1520              *Av+B^*p=f*
1521              *Bv     =0*
1522    
1523          for the unknowns *v* and *p* and given operators *A* and *B* and
1524          given right hand side *f*. *B^** is the adjoint operator of *B*.
1525          *A* may depend weakly on *v* and *p*.
1526          """
1527          def __init__(self, **kwargs):
1528            """
1529            initializes the saddle point problem
1530            """
1531            self.resetControlParameters()
1532            self.setTolerance()
1533            self.setAbsoluteTolerance()
1534          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):
1535             """
1536             sets a control parameter
1537    
1538             :param K_p: initial value for constant to adjust pressure tolerance
1539             :type K_p: ``float``
1540             :param K_v: initial value for constant to adjust velocity tolerance
1541             :type K_v: ``float``
1542             :param rtol_max: maximuim relative tolerance used to calculate presssure and velocity increment.
1543             :type rtol_max: ``float``
1544             :param chi_max: maximum tolerable converegence rate.
1545             :type chi_max: ``float``
1546             :param reduction_factor: reduction factor for adjustment factors.
1547             :type reduction_factor: ``float``
1548             """
1549             self.setControlParameter(K_p, K_v, rtol_max, rtol_min, chi_max, reduction_factor, theta)
1550    
1551          def setControlParameter(self,K_p=None, K_v=None, rtol_max=None, rtol_min=None, chi_max=None, reduction_factor=None, theta=None):
1552             """
1553             sets a control parameter
1554    
1555    
1556             :param K_p: initial value for constant to adjust pressure tolerance
1557             :type K_p: ``float``
1558             :param K_v: initial value for constant to adjust velocity tolerance
1559             :type K_v: ``float``
1560             :param rtol_max: maximuim relative tolerance used to calculate presssure and velocity increment.
1561             :type rtol_max: ``float``
1562             :param chi_max: maximum tolerable converegence rate.
1563             :type chi_max: ``float``
1564             :type reduction_factor: ``float``
1565             """
1566             if not K_p == None:
1567                if K_p<1:
1568                   raise ValueError("K_p need to be greater or equal to 1.")
1569             else:
1570                K_p=self.__K_p
1571    
1572             if not K_v == None:
1573                if K_v<1:
1574                   raise ValueError("K_v need to be greater or equal to 1.")
1575             else:
1576                K_v=self.__K_v
1577    
1578             if not rtol_max == None:
1579                if rtol_max<=0 or rtol_max>=1:
1580                   raise ValueError("rtol_max needs to be positive and less than 1.")
1581             else:
1582                rtol_max=self.__rtol_max
1583    
1584             if not rtol_min == None:
1585                if rtol_min<=0 or rtol_min>=1:
1586                   raise ValueError("rtol_min needs to be positive and less than 1.")
1587             else:
1588                rtol_min=self.__rtol_min
1589    
1590             if not chi_max == None:
1591                if chi_max<=0 or chi_max>=1:
1592                   raise ValueError("chi_max needs to be positive and less than 1.")
1593             else:
1594                chi_max = self.__chi_max
1595    
1596             if not reduction_factor == None:
1597                if reduction_factor<=0 or reduction_factor>1:
1598                   raise ValueError("reduction_factor need to be between zero and one.")
1599             else:
1600                reduction_factor=self.__reduction_factor
1601    
1602             if not theta == None:
1603                if theta<=0 or theta>1:
1604                   raise ValueError("theta need to be between zero and one.")
1605             else:
1606                theta=self.__theta
1607    
1608             if rtol_min>=rtol_max:
1609                 raise ValueError("rtol_max = %e needs to be greater than rtol_min = %e"%(rtol_max,rtol_min))
1610             self.__chi_max = chi_max
1611             self.__rtol_max = rtol_max
1612             self.__K_p = K_p
1613             self.__K_v = K_v
1614             self.__reduction_factor = reduction_factor
1615             self.__theta = theta
1616             self.__rtol_min=rtol_min
1617    
1618          #=============================================================
1619          def inner_pBv(self,p,Bv):
1620             """
1621             Returns inner product of element p and Bv (overwrite).
1622    
1623             :param p: a pressure increment
1624             :param Bv: a residual
1625             :return: inner product of element p and Bv
1626             :rtype: ``float``
1627             :note: used if PCG is applied.
1628             """
1629             raise NotImplementedError("no inner product for p and Bv implemented.")
1630    
1631          def inner_p(self,p0,p1):
1632             """
1633             Returns inner product of p0 and p1 (overwrite).
1634    
1635             :param p0: a pressure
1636             :param p1: a pressure
1637             :return: inner product of p0 and p1
1638             :rtype: ``float``
1639             """
1640             raise NotImplementedError("no inner product for p implemented.")
1641      
1642          def norm_v(self,v):
1643             """
1644             Returns the norm of v (overwrite).
1645    
1646             :param v: a velovity
1647             :return: norm of v
1648             :rtype: non-negative ``float``
1649             """
1650             raise NotImplementedError("no norm of v implemented.")
1651          def getDV(self, p, v, tol):
1652             """
1653             return a correction to the value for a given v and a given p with accuracy `tol` (overwrite)
1654    
1655             :param p: pressure
1656             :param v: pressure
1657             :return: dv given as *dv= A^{-1} (f-A v-B^*p)*
1658             :note: Only *A* may depend on *v* and *p*
1659             """
1660             raise NotImplementedError("no dv calculation implemented.")
1661    
1662            
1663          def Bv(self,v, tol):
1664            """
1665            Returns Bv with accuracy `tol` (overwrite)
1666    
1667            :rtype: equal to the type of p
1668            :note: boundary conditions on p should be zero!
1669            """
1670            raise NotImplementedError("no operator B implemented.")
1671    
1672          def norm_Bv(self,Bv):
1673            """
1674            Returns the norm of Bv (overwrite).
1675    
1676            :rtype: equal to the type of p
1677            :note: boundary conditions on p should be zero!
1678            """
1679            raise NotImplementedError("no norm of Bv implemented.")
1680    
1681          def solve_AinvBt(self,dp, tol):
1682             """
1683             Solves *A dv=B^*dp* with accuracy `tol`
1684    
1685             :param dp: a pressure increment
1686             :return: the solution of *A dv=B^*dp*
1687             :note: boundary conditions on dv should be zero! *A* is the operator used in ``getDV`` and must not be altered.
1688             """
1689             raise NotImplementedError("no operator A implemented.")
1690    
1691          def solve_prec(self,Bv, tol):
1692             """
1693             Provides a preconditioner for *(BA^{-1}B^ * )* applied to Bv with accuracy `tol`
1694    
1695             :rtype: equal to the type of p
1696             :note: boundary conditions on p should be zero!
1697             """
1698             raise NotImplementedError("no preconditioner for Schur complement implemented.")
1699          #=============================================================
1700          def __Aprod_PCG(self,dp):
1701              dv=self.solve_AinvBt(dp, self.__subtol)
1702              return ArithmeticTuple(dv,self.Bv(dv, self.__subtol))
1703    
1704          def __inner_PCG(self,p,r):
1705             return self.inner_pBv(p,r[1])
1706    
1707          def __Msolve_PCG(self,r):
1708              return self.solve_prec(r[1], self.__subtol)
1709          #=============================================================
1710          def __Aprod_GMRES(self,p):
1711              return self.solve_prec(self.Bv(self.solve_AinvBt(p, self.__subtol), self.__subtol), self.__subtol)
1712          def __inner_GMRES(self,p0,p1):
1713             return self.inner_p(p0,p1)
1714    
1715          #=============================================================
1716          def norm_p(self,p):
1717              """
1718              calculates the norm of ``p``
1719              
1720              :param p: a pressure
1721              :return: the norm of ``p`` using the inner product for pressure
1722              :rtype: ``float``
1723              """
1724              f=self.inner_p(p,p)
1725              if f<0: raise ValueError("negative pressure norm.")
1726              return math.sqrt(f)
1727          
1728          def solve(self,v,p,max_iter=20, verbose=False, usePCG=True, iter_restart=20, max_correction_steps=10):
1729             """
1730             Solves the saddle point problem using initial guesses v and p.
1731    
1732             :param v: initial guess for velocity
1733             :param p: initial guess for pressure
1734             :type v: `Data`
1735             :type p: `Data`
1736             :param usePCG: indicates the usage of the PCG rather than GMRES scheme.
1737             :param max_iter: maximum number of iteration steps per correction
1738                              attempt
1739             :param verbose: if True, shows information on the progress of the
1740                             saddlepoint problem solver.
1741             :param iter_restart: restart the iteration after ``iter_restart`` steps
1742                                  (only used if useUzaw=False)
1743             :type usePCG: ``bool``
1744             :type max_iter: ``int``
1745             :type verbose: ``bool``
1746             :type iter_restart: ``int``
1747             :rtype: ``tuple`` of `Data` objects
1748             :note: typically this method is overwritten by a subclass. It provides a wrapper for the ``_solve`` method.
1749             """
1750             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)
1751    
1752          def _solve(self,v,p,max_iter=20, verbose=False, usePCG=True, iter_restart=20, max_correction_steps=10):
1753             """
1754             see `_solve` method.
1755             """
1756             self.verbose=verbose
1757             rtol=self.getTolerance()
1758             atol=self.getAbsoluteTolerance()
1759    
1760             K_p=self.__K_p
1761             K_v=self.__K_v
1762             correction_step=0
1763             converged=False
1764             chi=None
1765             eps=None
1766    
1767             if self.verbose: print(("HomogeneousSaddlePointProblem: start iteration: rtol= %e, atol=%e"%(rtol, atol)))
1768             while not converged:
1769    
1770                 # get tolerance for velecity increment:
1771                 if chi == None:
1772                    rtol_v=self.__rtol_max
1773                 else:
1774                    rtol_v=min(chi/K_v,self.__rtol_max)
1775                 rtol_v=max(rtol_v, self.__rtol_min)
1776                 if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: rtol_v= %e"%(correction_step,rtol_v)))
1777                 # get velocity increment:
1778                 dv1=self.getDV(p,v,rtol_v)
1779                 v1=v+dv1
1780                 Bv1=self.Bv(v1, rtol_v)
1781                 norm_Bv1=self.norm_Bv(Bv1)
1782                 norm_dv1=self.norm_v(dv1)
1783                 if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: norm_Bv1 = %e, norm_dv1 = %e"%(correction_step, norm_Bv1, norm_dv1)))
1784                 if norm_dv1*self.__theta < norm_Bv1:
1785                    # get tolerance for pressure increment:
1786                    large_Bv1=True
1787                    if chi == None or eps == None:
1788                       rtol_p=self.__rtol_max
1789                    else:
1790                       rtol_p=min(chi**2*eps/K_p/norm_Bv1, self.__rtol_max)
1791                    self.__subtol=max(rtol_p**2, self.__rtol_min)
1792                    if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: rtol_p= %e"%(correction_step,rtol_p)))
1793                    # now we solve for the pressure increment dp from B*A^{-1}B^* dp = Bv1
1794                    if usePCG:
1795                        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)
1796                        v2=r[0]
1797                        Bv2=r[1]
1798                    else:
1799                        # don't use!!!!
1800                        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)
1801                        dv2=self.solve_AinvBt(dp, self.__subtol)
1802                        v2=v1-dv2
1803                        Bv2=self.Bv(v2, self.__subtol)
1804                    p2=p+dp
1805                 else:
1806                    large_Bv1=False
1807                    v2=v1
1808                    p2=p
1809                 # update business:
1810                 norm_dv2=self.norm_v(v2-v)
1811                 norm_v2=self.norm_v(v2)
1812                 if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: v2 = %e, norm_dv2 = %e"%(correction_step, norm_v2, self.norm_v(v2-v))))
1813                 eps, eps_old = max(norm_Bv1, norm_dv2), eps
1814                 if eps_old == None:
1815                      chi, chi_old = None, chi
1816                 else:
1817                      chi, chi_old = min(eps/ eps_old, self.__chi_max), chi
1818                 if eps != None:
1819                     if chi !=None:
1820                        if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: convergence rate = %e, correction = %e"%(correction_step,chi, eps)))
1821                     else:
1822                        if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: correction = %e"%(correction_step, eps)))
1823                 if eps <= rtol*norm_v2+atol :
1824                     converged = True
1825                 else:
1826                     if correction_step>=max_correction_steps:
1827                          raise CorrectionFailed("Given up after %d correction steps."%correction_step)
1828                     if chi_old!=None:
1829                        K_p=max(1,self.__reduction_factor*K_p,(chi-chi_old)/chi_old**2*K_p)
1830                        K_v=max(1,self.__reduction_factor*K_v,(chi-chi_old)/chi_old**2*K_p)
1831                        if self.verbose: print(("HomogeneousSaddlePointProblem: step %s: new adjustment factor K = %e"%(correction_step,K_p)))
1832                 correction_step+=1
1833                 v,p =v2, p2
1834             if self.verbose: print(("HomogeneousSaddlePointProblem: tolerance reached after %s steps."%correction_step))
1835             return v,p
1836          #========================================================================
1837          def setTolerance(self,tolerance=1.e-4):
1838             """
1839             Sets the relative tolerance for (v,p).
1840    
1841             :param tolerance: tolerance to be used
1842             :type tolerance: non-negative ``float``
1843             """
1844             if tolerance<0:
1845                 raise ValueError("tolerance must be positive.")
1846             self.__rtol=tolerance
1847    
1848          def getTolerance(self):
1849             """
1850             Returns the relative tolerance.
1851    
1852             :return: relative tolerance
1853             :rtype: ``float``
1854             """
1855             return self.__rtol
1856    
1857          def setAbsoluteTolerance(self,tolerance=0.):
1858             """
1859             Sets the absolute tolerance.
1860    
1861             :param tolerance: tolerance to be used
1862             :type tolerance: non-negative ``float``
1863             """
1864             if tolerance<0:
1865                 raise ValueError("tolerance must be non-negative.")
1866             self.__atol=tolerance
1867    
1868          def getAbsoluteTolerance(self):
1869             """
1870             Returns the absolute tolerance.
1871    
1872             :return: absolute tolerance
1873             :rtype: ``float``
1874             """
1875             return self.__atol
1876    
1877    
1878    def MaskFromBoundaryTag(domain,*tags):
1879       """
1880       Creates a mask on the Solution(domain) function space where the value is
1881       one for samples that touch the boundary tagged by tags.
1882    
1883       Usage: m=MaskFromBoundaryTag(domain, "left", "right")
1884    
1885       :param domain: domain to be used
1886       :type domain: `escript.Domain`
1887       :param tags: boundary tags
1888       :type tags: ``str``
1889       :return: a mask which marks samples that are touching the boundary tagged
1890                by any of the given tags
1891       :rtype: `escript.Data` of rank 0
1892       """
1893       pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1894       d=escript.Scalar(0.,escript.FunctionOnBoundary(domain))
1895       for t in tags: d.setTaggedValue(t,1.)
1896       pde.setValue(y=d)
1897       return util.whereNonZero(pde.getRightHandSide())
1898    
1899    def MaskFromTag(domain,*tags):
1900       """
1901       Creates a mask on the Solution(domain) function space where the value is
1902       one for samples that touch regions tagged by tags.
1903    
1904       Usage: m=MaskFromTag(domain, "ham")
1905    
1906       :param domain: domain to be used
1907       :type domain: `escript.Domain`
1908       :param tags: boundary tags
1909       :type tags: ``str``
1910       :return: a mask which marks samples that are touching the boundary tagged
1911                by any of the given tags
1912       :rtype: `escript.Data` of rank 0
1913       """
1914       pde=linearPDEs.LinearPDE(domain,numEquations=1, numSolutions=1)
1915       d=escript.Scalar(0.,escript.Function(domain))
1916       for t in tags: d.setTaggedValue(t,1.)
1917       pde.setValue(Y=d)
1918       return util.whereNonZero(pde.getRightHandSide())
1919    
1920    
 # vim: expandtab shiftwidth=4:  

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

  ViewVC Help
Powered by ViewVC 1.1.26