linearspace.hh

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
/*                                                                           */
/*  This file is part of the library KASKADE 7                               */
/*  https://www.zib.de/research/projects/kaskade7-finite-element-toolbox     */
/*                                                                           */
/*  Copyright (C) 2002-2023 Zuse Institute Berlin                            */
/*                                                                           */
/*  KASKADE 7 is distributed under the terms of the ZIB Academic License.    */
/*    see $KASKADE/academic.txt                                              */
/*                                                                           */
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#ifndef LINEARSPACE_HH
#define LINEARSPACE_HH
 
#include <cmath>
#include <type_traits>
 
#include <boost/version.hpp>
#include <boost/fusion/algorithm.hpp>
#include <boost/fusion/sequence.hpp>
 
#include "fem/functionspace.hh"
 
#include "linalg/crsutil.hh"
 
namespace Kaskade
{
  // forward declaration
  template <class> class VariableSet;
 
  namespace LinearSpace_Detail
  {
 
    template <class Space0, class Space2, int id> struct Copy;
 
    template <class Space, int id, class VSDescriptions>
    struct Copy<VariableSet<VSDescriptions>,Space,id>
    {
      static void apply(const VariableSet<VSDescriptions>& from, Space& to)
      {
        boost::fusion::at_c<id>(to.data) = boost::fusion::at_c<id>(from.data).coefficients();
        Copy<VariableSet<VSDescriptions>,Space,id-1>::apply(from,to);
      }
    };
 
    template <class Space, class VSDescriptions>
    struct Copy<VariableSet<VSDescriptions>,Space,0>
    {
      static void apply(const VariableSet<VSDescriptions>& from, Space& to)
      {
        boost::fusion::at_c<0>(to.data) = boost::fusion::at_c<0>(from.data).coefficients();
      }
    };
  }
 
  template <class Scalar_, class Seq>
  class LinearProductSpace
  {
  private:
    typedef LinearProductSpace<Scalar_,Seq> Self;
 
    // forward declarations (definitions at end of file)
    struct Add;
    struct Sub;
    struct Assign;
    struct ScalarMult;
    struct Axpy;
    struct ScalarProduct;
    template <typename> struct ReadBlock;
    template <typename> struct WriteBlock;
 
  public:
    typedef Scalar_ Scalar;
    typedef Scalar field_type; // for compatibility with Dune conventions
    
    typedef Seq Sequence;
    typedef Seq Functions;
    
    template <int m>
    using Component = std::remove_reference_t<typename boost::fusion::result_of::at_c<Sequence const,m>::type>;
 
 
    // use with care: default constructed subvectors may not be properly initialized (e.g. size?)
    LinearProductSpace() = delete;
 
    LinearProductSpace(LinearProductSpace<Scalar,Sequence> const& y): data(y.data) {}
    
    LinearProductSpace(LinearProductSpace<Scalar,Sequence>&& y): data(std::move(y.data)) {}
 
    template <class VSDescriptions>
    explicit LinearProductSpace(const VariableSet<VSDescriptions>& y)
    {
      *this = y;
    }
 
    template <class S>
    explicit LinearProductSpace(S const& init): data(init)
    {}
 
    template <typename... Args>
    explicit LinearProductSpace(boost::fusion::transform_view<Args...> const& data_) : data(data_)
    {
      applyUnaryOp(Assign(0));
    }
    
    size_t dim() const
    {
      return boost::fusion::accumulate(data,0,[](auto dim, auto const& v) { return dim+v.dim(); });
    }
 
    Self& operator=(Self const& y) { if (this!=&y) data = y.data; return *this; }
 
    template <class VSDescriptions>
    Self& operator=(VariableSet<VSDescriptions> const& y)
    {
      using namespace boost::fusion;
      static_assert(result_of::size<Sequence>::type::value == VSDescriptions::noOfVariables,"Numbers of variables do not match.");
      LinearSpace_Detail::Copy<VariableSet<VSDescriptions>,LinearProductSpace,
                               result_of::size<typename VSDescriptions::RepresentationData>::type::value-1>::apply(y,*this);
      return *this;
    } // need not check for self assignment here
 
    template <class OtherSeq>
    Self& operator=(LinearProductSpace<Scalar,OtherSeq> const& y) 
    { 
      using namespace boost::fusion;
      static_assert(result_of::size<Sequence>::type::value == result_of::size<OtherSeq>::type::value, "Numbers of variables do not match.");
      data = y.data;
      return *this; 
    } // need not check for self assignment here
 
    template <class FunctionSpace, int m>
    Self& operator=(FunctionSpaceElement<FunctionSpace,m> const& fse) 
    { 
      data = fse.coefficients(); 
      return *this; 
    }
 
    Self& operator=(Scalar a) 
    { 
      return applyUnaryOp(Assign(a));
    }
    
    template <class First, class Last>
    Self& operator=(boost::fusion::iterator_range<First,Last> range)
    {
      static_assert(boost::fusion::result_of::size<Sequence>::type::value == boost::fusion::result_of::distance<First,Last>::type::value,
                    "Numbers of variables do not match.");
      data = range;
      return *this;
    }
    
    Self& operator*=(Scalar a) { return applyUnaryOp(ScalarMult(a)); }
 
    template <class SequenceY>
    Self& operator += (LinearProductSpace<Scalar,SequenceY> const& y)
    {
      return applyBinaryOp(Add(),y);
    }
 
    Self& operator-=(Self const& y) 
    { 
      return applyBinaryOp(Sub(),y); 
    }
 
    Self& axpy(Scalar a, Self const& y) 
    { 
      return applyBinaryOp(Axpy(a),y); 
    }
 
    template <class InIterator>
    void read(InIterator i) [[deprecated]] 
    { 
      vectorFromSequence(*this,i);
    }
 
    template <class InIterator>
    void read(int rbegin, int rend, InIterator i) { boost::fusion::for_each(data,ReadBlock<InIterator>(rbegin,rend,i)); }
 
    template <class OutIterator>
    void write(OutIterator i) const  [[deprecated]]
    { 
      vectorToSequence(*this,i);
    }
 
    template <class DataOutIter>
    void write(int rbegin, int rend, DataOutIter i) const
    {
      boost::fusion::for_each(data,WriteBlock<DataOutIter>(rbegin,rend,i));
    }
 
    field_type operator*(Self const& y) const 
    {
      using namespace boost::fusion;
      zip_view<vector<Sequence const&,Sequence const&> > zipped(vector<Sequence const&,Sequence const&>(data,y.data));
      field_type s = 0;
      for_each(zipped,[&](auto const& p) { s += at_c<0>(p) * at_c<1>(p); });
      return s;
    }
    
    field_type dot(Self const& y) const
    {
      return *this * y;
    }
    
    double two_norm2() const
    {
      return (*this)*(*this);
    }
    
    double two_norm() const 
    { 
      return std::sqrt(two_norm2()); 
    }
    
    Sequence data;
 
  private:
 
    // applies a binary operator to the elements of a pair
    template <class BinaryOp>
    struct PairOp
    {
      PairOp(BinaryOp const& op_): op(op_) {}
 
      template <class Pair> void operator()(Pair p) const {
        op(boost::fusion::at_c<0>(p),boost::fusion::at_c<1>(p));
      }
 
    private:
      BinaryOp const& op;
    };
 
    template <class UnaryOp>
    Self& applyUnaryOp(UnaryOp const& op) {
      boost::fusion::for_each(data,op);
      return *this;
    }
 
    template <class BinaryOp, class SpaceY>
    Self& applyBinaryOp(BinaryOp const& op, SpaceY const& y) {
      using namespace boost::fusion;
      typedef typename SpaceY::Sequence SequenceY;
      zip_view<vector<Sequence&,SequenceY const&> > zipped(vector<Sequence&,SequenceY const&>(data,y.data));
      for_each(zipped,PairOp<BinaryOp>(op));
      return *this;
    }
 
    struct Add { template <class T, class S> void operator()(T& a, S const& b) const { a += b; } };
    struct Sub { template <class T> void operator()(T& a, T const& b) const { a -= b; } };
    struct Assign {
      Assign(Scalar s_): s(s_) {}
      template <class T> void operator()(T& a) const { a = s; }
    private: Scalar s;
    };
 
 
    struct ScalarMult {
      ScalarMult(Scalar s_): s(s_) {}
      template <class Element> void operator()(Element& e) const { e *= s; }
    private: Scalar s;
    };
 
    struct Axpy
    {
      Axpy(Scalar s_): s(s_) {}
      template <class T> void operator()(T& a, T const& b) const {
 
        // Funny runtime error could only be avoided by this funny code.
        // Probably a compiler bug.
        // assert(a.size()==b.size());
        std::cout << "";
 
 
        a.axpy(s,b);
      }
 
    private: Scalar s;
    };
 
    template <class DataOutIter>
    struct WriteBlock
    {
      WriteBlock(int& rbegin_, int& rend_, DataOutIter& out_): rbegin(rbegin_), rend(rend_), out(out_) {}
      template <class VectorBlock> void operator()(VectorBlock const& v) const
      {
        if (rbegin<=0 && rend>0)
          vectorToSequence(v,out);
        --rbegin;
        --rend;
      }
    private:
      int& rbegin;
      int& rend;
      DataOutIter& out;
 
    };
 
    template <class DataInIter>
    struct ReadBlock
    {
      ReadBlock(int& rbegin_, int& rend_, DataInIter& in_): rbegin(rbegin_), rend(rend_), in(in_) {}
      
      template <class VectorBlock> void operator()(VectorBlock& v) const
      {
        if (rbegin<=0 && rend>0)
          vectorFromSequence(v,in);
        --rbegin;
        --rend;
      }
    private:
      int& rbegin;
      int& rend;
      DataInIter& in;
 
    };
 
  };
  
  // ----------------------------------------------------------------------------------------------
 
  template <int m, class Scalar, class Sequence>
  auto const& component(LinearProductSpace<Scalar,Sequence> const& x)
  {
    return boost::fusion::at_c<m>(x.data);
  }
 
  template <int m, class Scalar, class Sequence>
  auto& component(LinearProductSpace<Scalar,Sequence>& x)
  {
    return boost::fusion::at_c<m>(x.data);
  }
  
  // ----------------------------------------------------------------------------------------------
  
 
  template <class X0, class ...Xs>
  auto cartesianProduct(X0 const& x0, Xs const&... xs)
  {
    using Sequence = boost::fusion::vector<X0,Xs...>;
    return LinearProductSpace<typename X0::field_type,Sequence>(Sequence(x0,xs...));
  }
  
  // ----------------------------------------------------------------------------------------------
  
 
  template <class Scalar, class Seq, class OutIter>
  OutIter vectorToSequence(LinearProductSpace<Scalar,Seq> const& v, OutIter i)
  {
    boost::fusion::for_each(v.data,[&](auto& x) { i = vectorToSequence(x,i); }); 
    return i;
  }
  
  template <class Scalar, class Seq, class InIter>
  InIter vectorFromSequence(LinearProductSpace<Scalar,Seq>& v, InIter i)
  {
    boost::fusion::for_each(v.data,[&](auto& x) { i = vectorFromSequence(x,i); }); 
    return i;
  }
  
} // end of namespace Kaskade
 
 
#endif