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>
 
namespace FEAT
{
  namespace LAFEM
  {
    template<typename DT_, typename IT_>
    class MeanFilter
    {
    public:
      typedef DenseVector<DT_, IT_> VectorType;
      typedef typename VectorType::DataType DataType;
      typedef typename VectorType::IndexType IndexType;
 
      template <typename DT2_ = DT_, typename IT2_ = IT_>
      using FilterType = MeanFilter<DT2_, IT2_>;
 
      template <typename DataType2_, typename IndexType2_>
      using FilterTypeByDI = FilterType<DataType2_, IndexType2_>;
 
    protected:
      VectorType _vec_prim;
      VectorType _vec_dual;
      DataType _volume;
      DataType _sol_mean;
 
    public:
      // default CTOR
      MeanFilter() :
        _vec_prim(),
        _vec_dual(),
        _volume(),
        _sol_mean()
      {
      }
 
      explicit MeanFilter(VectorType && vec_prim, VectorType && vec_dual, DataType sol_mean = DataType(0)) :
        _vec_prim(std::forward<VectorType>(vec_prim)),
        _vec_dual(std::forward<VectorType>(vec_dual)),
        _volume(_vec_prim.dot(_vec_dual)),
        _sol_mean(sol_mean)
      {
        if(!_vec_prim.empty())
        {
          XASSERTM(_volume > Math::eps<DataType>(), "domain volume must not be zero");
        }
      }
 
      explicit MeanFilter(VectorType && vec_prim, VectorType && vec_dual, DataType sol_mean, DataType volume) :
        _vec_prim(std::forward<VectorType>(vec_prim)),
        _vec_dual(std::forward<VectorType>(vec_dual)),
        _volume(volume),
        _sol_mean(sol_mean)
      {
        if(!_vec_prim.empty())
        {
          XASSERTM(_volume > Math::eps<DataType>(), "domain volume must not be zero");
        }
      }
 
      MeanFilter(MeanFilter && other) :
        _vec_prim(std::move(other._vec_prim)),
        _vec_dual(std::move(other._vec_dual)),
        _volume(other._volume),
        _sol_mean(other._sol_mean)
      {
      }
 
      MeanFilter & operator=(MeanFilter && other)
      {
        if(this != &other)
        {
          _vec_prim = std::forward<VectorType>(other._vec_prim);
          _vec_dual = std::forward<VectorType>(other._vec_dual);
          _volume = other._volume;
          _sol_mean = other._sol_mean;
        }
        return *this;
      }
 
      virtual ~MeanFilter()
      {
      }
 
      MeanFilter clone(CloneMode clone_mode = CloneMode::Deep) const
      {
        return MeanFilter(_vec_prim.clone(clone_mode), _vec_dual.clone(clone_mode), _sol_mean, _volume);
      }
 
      void clone(const MeanFilter & other, CloneMode clone_mode = CloneMode::Deep)
      {
        _vec_prim.clone(other.get_vec_prim(), clone_mode);
        _vec_dual.clone(other.get_vec_dual(), clone_mode);
        _volume = other._volume;
        _sol_mean = other._sol_mean;
      }
 
      template<typename DT2_, typename IT2_>
      void convert(const MeanFilter<DT2_, IT2_>& other)
      {
        _vec_prim.convert(other.get_vec_prim());
        _vec_dual.convert(other.get_vec_dual());
        _volume = DataType(other.get_volume());
        _sol_mean = DataType(other.get_sol_mean());
        if(!_vec_prim.empty())
        {
          XASSERTM(_volume > Math::eps<DataType>(), "domain volume must not be zero");
        }
      }
 
      void clear()
      {
        _vec_prim.clear();
        _vec_dual.clear();
        _volume = _sol_mean = DataType(0);
      }
 
      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;
      }
 
      DataType get_volume() const
      {
        return _volume;
      }
 
      DataType get_sol_mean() const
      {
        return _sol_mean;
      }
 
      void filter_rhs(VectorType& vector) const
      {
        // subtract to dual integral mean zero
        if(!_vec_prim.empty())
          vector.axpy(_vec_dual, -vector.dot(_vec_prim) / _volume);
      }
 
      void filter_sol(VectorType& vector) const
      {
        // given: v=vec_prim, w=vec_dual, vol=v*w, input: x
        // note that int(v) = v*w and int(x) = x*w
        // we want:
        //      1/vol * int (x + alpha*v) = sol_mean
        // <==>     int(x) + alpha*int(v) = sol_mean*vol
        // <==>           x*w + alpha*v*w = sol_mean*v*w
        // <==>                 alpha*v*w = sol_mean*v*w - x*w
        // <==>                     alpha = sol_mean - (x*w/v*w)
        if(!_vec_prim.empty())
          vector.axpy(_vec_prim, _sol_mean - vector.dot(_vec_dual) / _volume);
      }
 
      void filter_def(VectorType& vector) const
      {
        // same as rhs
        filter_rhs(vector);
      }
 
      void filter_cor(VectorType& vector) const
      {
        // subtract to primal integral mean zero
        if(!_vec_prim.empty())
          vector.axpy(_vec_prim, -vector.dot(_vec_dual) / _volume);
      }
 
      template<typename MT_>
      void filter_mat(MT_& DOXY(matrix)) const
      {
        // nothing to do here
      }
    }; // class MeanFilter<...>
  } // namespace LAFEM
} // namespace FEAT