schwarz.hh

// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_ISTL_SCHWARZ_HH
#define DUNE_ISTL_SCHWARZ_HH
 
#include <iostream>              // for input/output to shell
#include <fstream>               // for input/output to files
#include <vector>                // STL vector class
#include <sstream>
 
#include <cmath>                // Yes, we do some math here
 
#include <dune/common/timer.hh>
 
#include "io.hh"
#include "bvector.hh"
#include "vbvector.hh"
#include "bcrsmatrix.hh"
#include "io.hh"
#include "gsetc.hh"
#include "ilu.hh"
#include "operators.hh"
#include "solvers.hh"
#include "preconditioners.hh"
#include "scalarproducts.hh"
#include "owneroverlapcopy.hh"
 
namespace Dune {
 
  template<class M, class X, class Y, class C>
  class OverlappingSchwarzOperator : public AssembledLinearOperator<M,X,Y>
  {
  public:
    typedef M matrix_type;
    typedef X domain_type;
    typedef Y range_type;
    typedef typename X::field_type field_type;
    typedef C communication_type;
 
    OverlappingSchwarzOperator (const matrix_type& A, const communication_type& com)
      : _A_(stackobject_to_shared_ptr(A)), communication(com)
    {}
 
    OverlappingSchwarzOperator (const std::shared_ptr<matrix_type> A, const communication_type& com)
      : _A_(A), communication(com)
    {}
 
    virtual void apply (const X& x, Y& y) const
    {
      y = 0;
      _A_->umv(x,y);     // result is consistent on interior+border
      communication.project(y);     // we want this here to avoid it before the preconditioner
                                    // since there d is const!
    }
 
    virtual void applyscaleadd (field_type alpha, const X& x, Y& y) const
    {
      _A_->usmv(alpha,x,y);     // result is consistent on interior+border
      communication.project(y);     // we want this here to avoid it before the preconditioner
                                    // since there d is const!
    }
 
    virtual const matrix_type& getmat () const
    {
      return *_A_;
    }
 
    virtual SolverCategory::Category category() const
    {
      return SolverCategory::overlapping;
    }
 
  private:
    const std::shared_ptr<const matrix_type>_A_;
    const communication_type& communication;
  };
 
  /*
   * @addtogroup ISTL_Prec
   * @{
   */
  template<class M, class X, class Y, class C>
  class ParSSOR : public Preconditioner<X,Y> {
  public:
    typedef M matrix_type;
    typedef X domain_type;
    typedef Y range_type;
    typedef typename X::field_type field_type;
    typedef C communication_type;
 
    ParSSOR (const matrix_type& A, int n, field_type w, const communication_type& c)
      : _A_(A), _n(n), _w(w), communication(c)
    {   }
 
    virtual void pre (X& x, Y& b)
    {
      communication.copyOwnerToAll(x,x);     // make dirichlet values consistent
    }
 
    virtual void apply (X& v, const Y& d)
    {
      for (int i=0; i<_n; i++) {
        bsorf(_A_,v,d,_w);
        bsorb(_A_,v,d,_w);
      }
      communication.copyOwnerToAll(v,v);
    }
 
    virtual void post (X& x) {}
 
    virtual SolverCategory::Category category() const
    {
      return SolverCategory::overlapping;
    }
 
  private:
    const matrix_type& _A_;
    int _n;
    field_type _w;
    const communication_type& communication;
  };
 
  namespace Amg
  {
    template<class T> class ConstructionTraits;
  }
 
  template<class X, class Y, class C, class P=Preconditioner<X,Y> >
  class BlockPreconditioner : public Preconditioner<X,Y> {
    friend struct Amg::ConstructionTraits<BlockPreconditioner<X,Y,C,P> >;
  public:
    typedef X domain_type;
    typedef Y range_type;
    typedef typename X::field_type field_type;
    typedef C communication_type;
 
    BlockPreconditioner (P& p, const communication_type& c)
      : _preconditioner(stackobject_to_shared_ptr(p)), _communication(c)
    {   }
 
    BlockPreconditioner (const std::shared_ptr<P>& p, const communication_type& c)
      : _preconditioner(p), _communication(c)
    {   }
 
    virtual void pre (X& x, Y& b)
    {
      _communication.copyOwnerToAll(x,x);     // make dirichlet values consistent
      _preconditioner->pre(x,b);
    }
 
    virtual void apply (X& v, const Y& d)
    {
      _preconditioner->apply(v,d);
      _communication.copyOwnerToAll(v,v);
    }
 
    template<bool forward>
    void apply (X& v, const Y& d)
    {
      _preconditioner->template apply<forward>(v,d);
      _communication.copyOwnerToAll(v,v);
    }
 
    virtual void post (X& x)
    {
      _preconditioner->post(x);
    }
 
    virtual SolverCategory::Category category() const
    {
      return SolverCategory::overlapping;
    }
 
  private:
    std::shared_ptr<P> _preconditioner;
 
    const communication_type& _communication;
  };
 
} // end namespace
 
#endif