mean_filter.hpp

// FEAT3: Finite Element Analysis Toolbox, Version 3
// Copyright (C) 2010 by Stefan Turek & the FEAT group
// FEAT3 is released under the GNU General Public License version 3,
// see the file 'copyright.txt' in the top level directory for details.
 
#pragma once
 
// includes, FEAT
#include <kernel/lafem/dense_vector.hpp>
#include <kernel/util/dist.hpp>
 
namespace FEAT
{
  namespace Global
  {
    template<typename DT_, typename IT_>
    class MeanFilter
    {
    public:
      typedef LAFEM::DenseVector<DT_, IT_> VectorType;
      typedef typename VectorType::DataType DataType;
      typedef typename VectorType::IndexType IndexType;
 
      template<typename DT2_, typename IT2_>
      using FilterType = MeanFilter<DT2_, IT2_>;
 
      template <typename DataType2_, typename IndexType2_>
      using FilterTypeByDI = FilterType<DataType2_, IndexType2_>;
 
    protected:
      const Dist::Comm* _comm;
      VectorType _vec_prim;
      VectorType _vec_dual;
      VectorType _vec_freq;
      DataType _volume;
 
    public:
      // default CTOR
      MeanFilter() :
        _comm(nullptr),
        _vec_prim(),
        _vec_dual(),
        _vec_freq(),
        _volume(DataType(0))
      {
      }
 
      explicit MeanFilter(VectorType&& vec_prim, VectorType&& vec_dual, VectorType&& vec_freq, const Dist::Comm* comm) :
        _comm(comm),
        _vec_prim(std::forward<VectorType>(vec_prim)),
        _vec_dual(std::forward<VectorType>(vec_dual)),
        _vec_freq(std::forward<VectorType>(vec_freq)),
        _volume()
      {
        if((!_vec_freq.empty()) && (_comm != nullptr))
        {
          _volume = _vec_freq.triple_dot(_vec_prim, _vec_dual);
          _comm->allreduce(&_volume, &_volume, std::size_t(1), Dist::op_sum);
        }
        else
          _volume = _vec_prim.dot(_vec_dual);
 
        //XASSERTM(_volume > Math::eps<DataType>(), "domain volume must not be zero");
      }
 
      MeanFilter(MeanFilter && other) :
        _comm(other._comm),
        _vec_prim(std::move(other._vec_prim)),
        _vec_dual(std::move(other._vec_dual)),
        _vec_freq(std::move(other._vec_freq)),
        _volume(other._volume)
      {
      }
 
      MeanFilter & operator=(MeanFilter && other)
      {
        if(this != &other)
        {
          _comm = other._comm;
          _vec_prim = std::move(other._vec_prim);
          _vec_dual = std::move(other._vec_dual);
          _vec_freq = std::move(other._vec_freq);
          _volume = other._volume;
        }
        return *this;
      }
 
      virtual ~MeanFilter()
      {
      }
 
      MeanFilter clone(LAFEM::CloneMode clone_mode = LAFEM::CloneMode::Deep) const
      {
        return MeanFilter(_vec_prim.clone(clone_mode), _vec_dual.clone(clone_mode), _vec_freq.clone(clone_mode), _comm);
      }
 
      void clone(const MeanFilter & other, LAFEM::CloneMode clone_mode = LAFEM::CloneMode::Deep)
      {
        if(this == &other)
          return;
 
        _comm = other.get_comm();
        _vec_prim.clone(other.get_vec_prim(), clone_mode);
        _vec_dual.clone(other.get_vec_dual(), clone_mode);
        _vec_freq.clone(other.get_vec_freq(), clone_mode);
        _volume = other.get_volume();
      }
 
      template<typename DT2_, typename IT2_>
      void convert(const MeanFilter<DT2_, IT2_>& other)
      {
        if((void*)this == (void*)&other)
          return;
 
        _comm = other.get_comm();
        _vec_prim.convert(other.get_vec_prim());
        _vec_dual.convert(other.get_vec_dual());
        _vec_freq.convert(other.get_vec_freq());
        _volume = DataType(other.get_volume());
      }
 
      VectorType & get_vec_prim()
      {
        return _vec_prim;
      }
 
      const VectorType & get_vec_prim() const
      {
        return _vec_prim;
      }
 
      VectorType & get_vec_dual()
      {
        return _vec_dual;
      }
 
      const VectorType & get_vec_dual() const
      {
        return _vec_dual;
      }
 
      VectorType & get_vec_freq()
      {
        return _vec_freq;
      }
 
      const VectorType & get_vec_freq() const
      {
        return _vec_freq;
      }
 
      DataType get_volume() const
      {
        return _volume;
      }
 
      const Dist::Comm* get_comm() const
      {
        return _comm;
      }
 
      std::size_t bytes() const
      {
        return _vec_prim.bytes() + _vec_dual.bytes();
      }
 
      bool empty() const
      {
        return _vec_prim.empty();
      }
 
      void filter_rhs(VectorType& vector) const
      {
        if(empty())
          return;
 
        // compute dual integral
        DataType integ = DataType(0);
        if(_vec_freq.empty() || (_comm == nullptr))
          integ = vector.dot(_vec_prim);
        else
        {
          integ = _vec_freq.triple_dot(vector, _vec_prim);
          _comm->allreduce(&integ, &integ, std::size_t(1), Dist::op_sum);
        }
        // subtract mean
        vector.axpy(_vec_dual, -integ / _volume);
      }
 
      void filter_sol(VectorType& vector) const
      {
        if(empty())
          return;
        // compute primal integral
        DataType integ = DataType(0);
        if(_vec_freq.empty() || (_comm == nullptr))
          integ = vector.dot(_vec_dual);
        else
        {
          integ = _vec_freq.triple_dot(vector, _vec_dual);
          _comm->allreduce(&integ, &integ, std::size_t(1), Dist::op_sum);
        }
        // subtract mean
        vector.axpy(_vec_prim, -integ / _volume);
      }
 
      void filter_def(VectorType& vector) const
      {
        // same as rhs
        filter_rhs(vector);
      }
 
      void filter_cor(VectorType& vector) const
      {
        // same as sol
        filter_sol(vector);
      }
    }; // class MeanFilter<...>
  } // namespace Global
} // namespace FEAT