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;
 
    enum {
      category=SolverCategory::overlapping
    };
 
    OverlappingSchwarzOperator (const 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_;
    }
 
  private:
    const matrix_type& _A_;
    const communication_type& communication;
  };
 
  template<class X, class C>
  class OverlappingSchwarzScalarProduct : public ScalarProduct<X>
  {
  public:
    typedef X domain_type;
    typedef typename X::field_type field_type;
    typedef typename FieldTraits<field_type>::real_type real_type;
    typedef C communication_type;
 
    enum {category=SolverCategory::overlapping};
 
    OverlappingSchwarzScalarProduct (const communication_type& com)
      : communication(com)
    {}
 
    virtual field_type dot (const X& x, const X& y)
    {
      field_type result;
      communication.dot(x,y,result);
      return result;
    }
 
    virtual real_type norm (const X& x)
    {
      return communication.norm(x);
    }
 
  private:
    const communication_type& communication;
  };
 
  template<class X, class C>
  struct ScalarProductChooser<X,C,SolverCategory::overlapping>
  {
    typedef OverlappingSchwarzScalarProduct<X,C> ScalarProduct;
    typedef C communication_type;
 
    enum {
      solverCategory=SolverCategory::overlapping
    };
 
    static ScalarProduct* construct(const communication_type& comm)
    {
      return new ScalarProduct(comm);
    }
  };
 
  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;
 
    // define the category
    enum {
      category=SolverCategory::overlapping
    };
 
    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) {}
 
  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 T=Preconditioner<X,Y> >
  class BlockPreconditioner : public Preconditioner<X,Y> {
    friend class Amg::ConstructionTraits<BlockPreconditioner<X,Y,C,T> >;
  public:
    typedef X domain_type;
    typedef Y range_type;
    typedef typename X::field_type field_type;
    typedef C communication_type;
 
    // define the category
    enum {
      category=SolverCategory::overlapping
    };
 
    BlockPreconditioner (T& 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);
    }
 
  private:
    T& preconditioner;
 
    const communication_type& communication;
  };
 
} // end namespace
 
#endif