stencil.tpp

/*---------------------------------------------------------------------------*\
 *
 *  bitpit
 *
 *  Copyright (C) 2015-2021 OPTIMAD engineering Srl
 *
 *  -------------------------------------------------------------------------
 *  License
 *  This file is part of bitpit.
 *
 *  bitpit is free software: you can redistribute it and/or modify it
 *  under the terms of the GNU Lesser General Public License v3 (LGPL)
 *  as published by the Free Software Foundation.
 *
 *  bitpit is distributed in the hope that it will be useful, but WITHOUT
 *  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 *  FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
 *  License for more details.
 *
 *  You should have received a copy of the GNU Lesser General Public License
 *  along with bitpit. If not, see <http://www.gnu.org/licenses/>.
 *
\*---------------------------------------------------------------------------*/
 
#ifndef __BITPIT_STENCIL_TPP__
#define __BITPIT_STENCIL_TPP__
 
namespace bitpit {

template<typename weight_t, typename value_t>
const typename DiscreteStencil<weight_t, value_t>::weight_manager_type DiscreteStencil<weight_t, value_t>::m_weightManager = typename DiscreteStencil<weight_t>::weight_manager_type();

template<typename weight_t, typename value_t>
const typename DiscreteStencil<weight_t, value_t>::weight_manager_type & DiscreteStencil<weight_t, value_t>::getWeightManager()
{
    return m_weightManager;
}

template<typename weight_t, typename value_t>
OBinaryStream& operator<<(OBinaryStream &buffer, const DiscreteStencil<weight_t, value_t> &stencil)
{
    buffer << stencil.m_zero;
 
    std::size_t nItems = stencil.size();
    buffer << nItems;
 
    const long *patternData = stencil.patternData();
    const weight_t *weightData = stencil.weightData();
    for (std::size_t n = 0; n < nItems; ++n) {
        buffer << patternData[n];
        buffer << weightData[n];
    }
 
    buffer << stencil.m_constant;
 
    return buffer;
}

/*!
* Input stream operator from class DiscreteStencil to communication buffer.
*
* \param[in] buffer is the input memory stream
* \param[in] stencil is the stencil to be streamed
* \result Returns the same input stream received in input.
*/
template<typename weight_t, typename value_t>
IBinaryStream& operator>>(IBinaryStream &buffer, DiscreteStencil<weight_t, value_t> &stencil)
{
    buffer >> stencil.m_zero;
 
    std::size_t nItems;
    buffer >> nItems;
 
    stencil.resize(nItems);
    long *patternData = stencil.patternData();
    weight_t *weightData = stencil.weightData();
    for (std::size_t n = 0; n < nItems; ++n) {
        buffer >> patternData[n];
        buffer >> weightData[n];
    }
 
    buffer >> stencil.m_constant;
 
    return buffer;
}

* Constructor
*
* \param zero is the value to be used as zero
*/
template<typename weight_t, typename value_t>
DiscreteStencil<weight_t, value_t>::DiscreteStencil(const weight_t &zero)
    : DiscreteStencil(0, zero)
{
}

/*!
* Constructor
*
* \param size is the stencil size, expressed in number of elements
* \param zero is the value to be used as zero
*/
template<typename weight_t, typename value_t>
DiscreteStencil<weight_t, value_t>::DiscreteStencil(std::size_t size, const weight_t &zero)
    : m_zero(zero),
      m_pattern(size, -1), m_weights(size, m_zero),
      m_constant(m_zero)
{
}

/*!
* Constructor
*
* \param size is the stencil size, expressed in number of elements
* \param pattern is the patterns of the stencil
* \param zero is the value to be used as zero
*/
template<typename weight_t, typename value_t>
DiscreteStencil<weight_t, value_t>::DiscreteStencil(std::size_t size, const long *pattern, const weight_t &zero)
    : m_zero(zero),
      m_pattern(pattern, pattern + size), m_weights(size, m_zero),
      m_constant(m_zero)
{
}

/*!
* Constructor
*
* \param size is the stencil size, expressed in number of elements
* \param pattern is the patterns of the stencil
* \param weights are the weights of the stencil
* \param zero is the value to be used as zero
*/
template<typename weight_t, typename value_t>
DiscreteStencil<weight_t, value_t>::DiscreteStencil(std::size_t size, const long *pattern, const weight_t *weights, const weight_t &zero)
    : m_zero(zero),
      m_pattern(pattern, pattern + size), m_weights(weights, weights + size),
      m_constant(m_zero)
{
}

* \param zero is the value to be used as zero
*/
template<typename weight_t, typename value_t>
void DiscreteStencil<weight_t, value_t>::initialize(std::size_t size, const weight_t &zero)
{
    m_weightManager.copy(zero, &m_zero);
 
    std::size_t previousSize = this->size();
    std::size_t commonSize   = std::min(previousSize, size);
    for (std::size_t n = 0; n < commonSize; ++n) {
        m_pattern[n] = -1;
        m_weightManager.copy(m_zero, m_weights.data() + n);
    }
    if (previousSize != size) {
        resize(size);
    }
 
    zeroConstant();
}

*/
template<typename weight_t, typename value_t>
void DiscreteStencil<weight_t, value_t>::initialize(std::size_t size, const long *pattern, const weight_t &zero)
{
    m_weightManager.copy(zero, &m_zero);
 
    std::size_t previousSize = this->size();
    std::size_t commonSize   = std::min(previousSize, size);
    for (std::size_t n = 0; n < commonSize; ++n) {
        m_pattern[n] = pattern[n];
        m_weightManager.copy(m_zero, m_weights.data() + n);
    }
    if (previousSize != size) {
        resize(size);
        for (std::size_t n = previousSize; n < size; ++n) {
            m_pattern[n] = pattern[n];
        }
    }
 
    zeroConstant();
}

template<typename weight_t, typename value_t>
void DiscreteStencil<weight_t, value_t>::initialize(std::size_t size, const long *pattern, const weight_t *weights, const weight_t &zero)
{
    m_weightManager.copy(zero, &m_zero);
 
    std::size_t previousSize = this->size();
    std::size_t commonSize   = std::min(previousSize, size);
    for (std::size_t n = 0; n < commonSize; ++n) {
        m_pattern[n] = pattern[n];
        m_weightManager.copy(weights[n], m_weights.data() + n);
    }
    if (size > previousSize) {
        reserve(size);
        for (std::size_t n = previousSize; n < size; ++n) {
            appendItem(pattern[n], weights[n]);
        }
    } else if (size < previousSize) {
        resize(size);
    }
 
    zeroConstant();
}

template<typename weight_t, typename value_t>
void DiscreteStencil<weight_t, value_t>::initialize(const DiscreteStencil<weight_t, value_t> &other)
{
    initialize(other.size(), other.m_pattern.data(), other.m_weights.data(), other.m_zero);
}

template<typename weight_t, typename value_t>
std::size_t DiscreteStencil<weight_t, value_t>::size() const
{
    return m_pattern.size();
}

template<typename weight_t, typename value_t>
void DiscreteStencil<weight_t, value_t>::resize(std::size_t size)
{
    m_pattern.resize(size, -1);
    m_weights.resize(size, m_zero);
}

template<typename weight_t, typename value_t>
void DiscreteStencil<weight_t, value_t>::reserve(std::size_t capacity)
{
    m_pattern.reserve(capacity);
    m_weights.reserve(capacity);
}

template<typename weight_t, typename value_t>
long & DiscreteStencil<weight_t, value_t>::getPattern(std::size_t pos)
{
    return m_pattern[pos];
}

template<typename weight_t, typename value_t>
const long & DiscreteStencil<weight_t, value_t>::getPattern(std::size_t pos) const
{
    return m_pattern[pos];
}

template<typename weight_t, typename value_t>
long * DiscreteStencil<weight_t, value_t>::patternData()
{
    return m_pattern.data();
}

template<typename weight_t, typename value_t>
const long * DiscreteStencil<weight_t, value_t>::patternData() const
{
    return m_pattern.data();
}

template<typename weight_t, typename value_t>
void DiscreteStencil<weight_t, value_t>::setPattern(std::size_t pos, long id)
{
    m_pattern[pos] = id;
}

template<typename weight_t, typename value_t>
weight_t & DiscreteStencil<weight_t, value_t>::getWeight(std::size_t pos)
{
    return m_weights[pos];
}

template<typename weight_t, typename value_t>
const weight_t & DiscreteStencil<weight_t, value_t>::getWeight(std::size_t pos) const
{
    return m_weights[pos];
}

template<typename weight_t, typename value_t>
weight_t * DiscreteStencil<weight_t, value_t>::weightData()
{
    return m_weights.data();
}

template<typename weight_t, typename value_t>
const weight_t * DiscreteStencil<weight_t, value_t>::weightData() const
{
    return m_weights.data();
}

template<typename weight_t, typename value_t>
void DiscreteStencil<weight_t, value_t>::setWeight(std::size_t pos, const weight_t &weight)
{
    m_weights[pos] = weight;
}

template<typename weight_t, typename value_t>
void DiscreteStencil<weight_t, value_t>::setWeight(std::size_t pos, weight_t &&weight)
{
    m_weights[pos] = std::move(weight);
}

template<typename weight_t, typename value_t>
void DiscreteStencil<weight_t, value_t>::sumWeight(std::size_t pos, const weight_t &weight, double factor)
{
    m_weightManager.sum(weight, factor, m_weights.data() + pos);
}

template<typename weight_t, typename value_t>
void DiscreteStencil<weight_t, value_t>::zeroWeight(std::size_t pos)
{
    m_weightManager.copy(m_zero, m_weights.data() + pos);
}

template<typename weight_t, typename value_t>
void DiscreteStencil<weight_t, value_t>::setItem(std::size_t pos, long id, const weight_t &weight)
{
    setPattern(pos, id);
    setWeight(pos, weight);
}

template<typename weight_t, typename value_t>
void DiscreteStencil<weight_t, value_t>::setItem(std::size_t pos, long id, weight_t &&weight)
{
    setPattern(pos, id);
    setWeight(pos, std::move(weight));
}

template<typename weight_t, typename value_t>
void DiscreteStencil<weight_t, value_t>::sumItem(long id, const weight_t &weight, double factor)
{
    weight_t *target = findWeight(id);
    if (target) {
        m_weightManager.sum(weight, factor, target);
    } else {
        appendItem(id, weight);
        if (factor != 1.) {
            m_weights.back() *= factor;
        }
    }
}

template<typename weight_t, typename value_t>
void DiscreteStencil<weight_t, value_t>::appendItem(long id, const weight_t &weight)
{
    m_pattern.push_back(id);
    appendWeight(weight);
}

template<typename weight_t, typename value_t>
void DiscreteStencil<weight_t, value_t>::appendItem(long id, weight_t &&weight)
{
    m_pattern.push_back(id);
    appendWeight(std::move(weight));
}

template<typename weight_t, typename value_t>
const weight_t & DiscreteStencil<weight_t, value_t>::getConstant() const
{
    return m_constant;
}

template<typename weight_t, typename value_t>
weight_t & DiscreteStencil<weight_t, value_t>::getConstant()
{
    return m_constant;
}

template<typename weight_t, typename value_t>
void DiscreteStencil<weight_t, value_t>::setConstant(const weight_t &constant)
{
    m_weightManager.copy(constant, &m_constant);
}

template<typename weight_t, typename value_t>
void DiscreteStencil<weight_t, value_t>::setConstant(weight_t &&constant)
{
    m_weightManager.move(std::move(constant), &m_constant);
}

template<typename weight_t, typename value_t>
void DiscreteStencil<weight_t, value_t>::sumConstant(const weight_t &constant, double factor)
{
    m_weightManager.sum(constant, factor, &m_constant);
}

template<typename weight_t, typename value_t>
void DiscreteStencil<weight_t, value_t>::zeroConstant()
{
    m_weightManager.copy(m_zero, &m_constant);
}

template<typename weight_t, typename value_t>
void DiscreteStencil<weight_t, value_t>::clear(bool release)
{
    m_pattern.clear();
    if (release) {
        m_pattern.shrink_to_fit();
    }
 
    clearWeights(release);
 
    zeroConstant();
}

template<typename weight_t, typename value_t>
void DiscreteStencil<weight_t, value_t>::sum(const DiscreteStencil<weight_t, value_t> &other, double factor)
{
    const std::size_t other_nItems = other.size();
    for (std::size_t n = 0; n < other_nItems; ++n) {
        sumItem(other.m_pattern[n], other.m_weights[n], factor);
    }
 
    sumConstant(other.m_constant, factor);
}

template<typename weight_t, typename value_t>
void DiscreteStencil<weight_t, value_t>::optimize(double tolerance)
{
    std::size_t nItems = size();
    for (std::size_t n = 0; n < nItems; ++n) {
        if (m_weightManager.isNegligible(m_weights[n], m_zero, tolerance)) {
            m_pattern.erase(m_pattern.begin() + n);
            m_weights.erase(m_weights.begin() + n);
            --n;
            --nItems;
        }
    }
}

template<typename weight_t, typename value_t>
void DiscreteStencil<weight_t, value_t>::renumber(const std::unordered_map<long, long> &map)
{
    renumber([&map] (long id) -> long { return map.at(id); });
}

template<typename weight_t, typename value_t>
template<typename Mapper>
void DiscreteStencil<weight_t, value_t>::renumber(const Mapper &map)
{
    for (long &id : m_pattern) {
        id = map(id);
    }
}

template<typename weight_t, typename value_t>
void DiscreteStencil<weight_t, value_t>::addComplementToZero(long id)
{
    const std::size_t nItems = size();
    if (nItems == 0) {
        return;
    }
 
    weight_t complement = m_zero;
    for (const weight_t &weight : m_weights) {
        m_weightManager.sum(weight, -1., &complement);
    }
 
    appendItem(id, complement);
}

template<typename weight_t, typename value_t>
void DiscreteStencil<weight_t, value_t>::zero()
{
    for (weight_t &weight : m_weights) {
        m_weightManager.copy(m_zero, &weight);
    }
 
    setConstant(m_zero);
}

template<typename weight_t, typename value_t>
weight_t * DiscreteStencil<weight_t, value_t>::findWeight(long id)
{
    auto patternBegin = m_pattern.cbegin();
    auto patternEnd   = m_pattern.cend();
    for (auto itr = patternBegin; itr != patternEnd; ++itr) {
        if (*itr == id) {
            return (m_weights.data() + std::distance(patternBegin, itr));
        }
    }
 
    return nullptr;
}

template<typename weight_t, typename value_t>
const weight_t * DiscreteStencil<weight_t, value_t>::findWeight(long id) const
{
    auto patternBegin = m_pattern.cbegin();
    auto patternEnd   = m_pattern.cend();
    for (auto itr = patternBegin; itr != patternEnd; ++itr) {
        if (*itr == id) {
            return (m_weights.data() + std::distance(patternBegin, itr));
        }
    }
 
    return nullptr;
}

template<typename weight_t, typename value_t>
void DiscreteStencil<weight_t, value_t>::appendWeight(weight_t &&weight)
{
    m_weights.push_back(std::move(weight));
}

template<typename weight_t, typename value_t>
void DiscreteStencil<weight_t, value_t>::appendWeight(const weight_t &weight)
{
    m_weights.push_back(weight);
}

template<typename weight_t, typename value_t>
void DiscreteStencil<weight_t, value_t>::clearWeights(bool release)
{
    m_weights.clear();
    if (release) {
        m_weights.shrink_to_fit();
    }
}

template<typename weight_t, typename value_t>
void DiscreteStencil<weight_t, value_t>::display(std::ostream &out, double factor) const
{
    const std::size_t nItems = size();
 
    weight_t sum = m_zero;
    for (std::size_t n = 0; n < nItems; ++n) {
        long id = m_pattern[n];
        weight_t weight = factor * m_weights[n];
        out << "   id: " << id << " weight: " << weight << std::endl;
        m_weightManager.sum(weight, 1., &sum);
    }
 
    out << " constant : " << (factor * m_constant) << std::endl;
    out << " sum      : " << sum << std::endl;
}

template<typename weight_t, typename value_t>
size_t DiscreteStencil<weight_t, value_t>::getBinarySize() const
{
    std::size_t nItems = size();
 
    return (sizeof(m_zero) + nItems * (sizeof(long) + sizeof(weight_t)) + sizeof(m_constant));
}

template<typename weight_t, typename value_t>
weight_t & DiscreteStencil<weight_t, value_t>::at(long id)
{
    return const_cast<weight_t &>(static_cast<const DiscreteStencil<weight_t, value_t> &>(*this).at(id));
}

template<typename weight_t, typename value_t>
const weight_t & DiscreteStencil<weight_t, value_t>::at(long id) const
{
    const weight_t *weight = findWeight(id);
    if (weight) {
        return *weight;
    }
 
    throw("The stencil doens't contain an item with the specified id");
}

template<typename weight_t, typename value_t>
weight_t & DiscreteStencil<weight_t, value_t>::operator[](long id)
{
    weight_t *weight = findWeight(id);
    if (weight) {
        return *weight;
    }
 
    appendItem(id, m_zero);
 
    return m_weights.back();
}

template<typename weight_t, typename value_t>
DiscreteStencil<weight_t, value_t> & DiscreteStencil<weight_t, value_t>::operator*=(double factor)
{
    for (weight_t &weight : m_weights) {
        weight *= factor;
    }
    m_constant *= factor;
 
    return *this;
}

template<typename weight_t, typename value_t>
DiscreteStencil<weight_t, value_t> & DiscreteStencil<weight_t, value_t>::operator/=(double factor)
{
    for (weight_t &weight : m_weights) {
        weight /= factor;
    }
    m_constant /= factor;
 
    return *this;
}

template<typename weight_t, typename value_t>
DiscreteStencil<weight_t, value_t> & DiscreteStencil<weight_t, value_t>::operator+=(const DiscreteStencil<weight_t, value_t> &other)
{
    sum(other, 1.);
 
    return *this;
}

template<typename weight_t, typename value_t>
DiscreteStencil<weight_t, value_t> & DiscreteStencil<weight_t, value_t>::operator-=(const DiscreteStencil<weight_t, value_t> &other)
{
    sum(other, - 1.);
 
    return *this;
}

template<typename weight_t, typename value_t>
MPDiscreteStencil<weight_t, value_t>::MPDiscreteStencil(const weight_t &zero)
    : DiscreteStencil<weight_t, value_t>(zero),
      m_weightPool(nullptr)
{
}

template<typename weight_t, typename value_t>
MPDiscreteStencil<weight_t, value_t>::MPDiscreteStencil(std::size_t size, const weight_t &zero)
    : DiscreteStencil<weight_t, value_t>(size, zero),
      m_weightPool(nullptr)
{
}

template<typename weight_t, typename value_t>
MPDiscreteStencil<weight_t, value_t>::MPDiscreteStencil(std::size_t size, const long *pattern, const weight_t &zero)
    : DiscreteStencil<weight_t, value_t>(size, pattern, zero),
      m_weightPool(nullptr)
{
}

template<typename weight_t, typename value_t>
MPDiscreteStencil<weight_t, value_t>::MPDiscreteStencil(std::size_t size, const long *pattern, const weight_t *weights, const weight_t &zero)
    : DiscreteStencil<weight_t, value_t>(size, pattern, weights, zero),
      m_weightPool(nullptr)
 
{
}

template<typename weight_t, typename value_t>
void MPDiscreteStencil<weight_t, value_t>::setWeightPool(MPDiscreteStencil<weight_t, value_t>::weight_pool_type *pool)
{
    m_weightPool = pool;
}

template<typename weight_t, typename value_t>
void MPDiscreteStencil<weight_t, value_t>::appendWeight(const weight_t &weight)
{
    if (m_weightPool && m_weightPool->size() > 0) {
        this->m_weights.emplace_back(m_weightPool->retrieve());
        this->m_weights.back() = weight;
    } else {
        DiscreteStencil<weight_t, value_t>::appendWeight(weight);
    }
}

template<typename weight_t, typename value_t>
void MPDiscreteStencil<weight_t, value_t>::clearWeights(bool release)
{
    if (m_weightPool) {
        m_weightPool->store(&(this->m_weights));
    }
 
    DiscreteStencil<weight_t, value_t>::clearWeights(release);
}
 
}

template<typename weight_t, typename value_t>
bitpit::DiscreteStencil<weight_t, value_t> operator*(const bitpit::DiscreteStencil<weight_t, value_t> &stencil, double factor)
{
    return (factor * stencil);
}

template<typename weight_t, typename value_t>
bitpit::DiscreteStencil<weight_t, value_t> operator*(double factor, const bitpit::DiscreteStencil<weight_t, value_t> &stencil)
{
    bitpit::DiscreteStencil<weight_t, value_t> stencil_result(stencil);
    stencil_result *= factor;
 
    return stencil_result;
}

template<typename weight_t, typename value_t>
bitpit::DiscreteStencil<weight_t, value_t> operator/(const bitpit::DiscreteStencil<weight_t, value_t> &stencil, double factor)
{
    bitpit::DiscreteStencil<weight_t, value_t> stencil_result(stencil);
    stencil_result /= factor;
 
    return stencil_result;
}

template<typename weight_t, typename value_t>
bitpit::DiscreteStencil<weight_t, value_t> operator+(const bitpit::DiscreteStencil<weight_t, value_t> &stencil_A, const bitpit::DiscreteStencil<weight_t, value_t> &stencil_B)
{
    bitpit::DiscreteStencil<weight_t, value_t> stencil_result(stencil_A);
    stencil_result += stencil_B;
 
    return stencil_result;
}

template<typename weight_t, typename value_t>
bitpit::DiscreteStencil<weight_t, value_t> operator-(const bitpit::DiscreteStencil<weight_t, value_t> &stencil_A, const bitpit::DiscreteStencil<weight_t, value_t> &stencil_B)
{
    bitpit::DiscreteStencil<weight_t, value_t> stencil_result(stencil_A);
    stencil_result -= stencil_B;
 
    return stencil_result;
}

template <typename V>
typename bitpit::DiscreteStencil<std::array<V, 3>> operator*(const typename bitpit::DiscreteStencil<V> &stencil, const std::array<V, 3> &vector)
{
    return (vector * stencil);
}

template <typename V>
typename bitpit::DiscreteStencil<std::array<V, 3>> operator*(const std::array<V, 3> &vector, const typename bitpit::DiscreteStencil<V> &stencil)
{
    const std::size_t nItems = stencil.size();
    bitpit::StencilVector stencil_B;
    stencil_B.resize(nItems);
    for (std::size_t n = 0; n < nItems; ++n) {
        stencil_B.setPattern(n, stencil.getPattern(n));
        stencil_B.setWeight(n, ::operator*(stencil.getWeight(n), vector));
    }
    stencil_B.setConstant(::operator*(stencil.getConstant(), vector));
 
    return stencil_B;
}

template <typename V>
typename bitpit::DiscreteStencil<V> dotProduct(const typename bitpit::DiscreteStencil<std::array<V, 3>> &stencil, const typename bitpit::DiscreteStencil<std::array<V, 3>>::weight_type &vector)
{
    typename bitpit::DiscreteStencil<V> stencil_dotProduct;
    dotProduct(stencil, vector, &stencil_dotProduct);
 
    return stencil_dotProduct;
}

template <typename V>
void dotProduct(const typename bitpit::DiscreteStencil<std::array<V, 3>> &stencil, const typename bitpit::DiscreteStencil<std::array<V, 3>>::weight_type &vector, typename bitpit::DiscreteStencil<V> *stencil_dotProduct)
{
    const std::size_t nItems = stencil.size();
    stencil_dotProduct->resize(nItems);
 
    const long *patternData = stencil.patternData();
    const bitpit::StencilVector::weight_type *weightData = stencil.weightData();
    long *patternData_dotProduct = stencil_dotProduct->patternData();
    typename bitpit::DiscreteStencil<V>::weight_type *weightData_dotProduct = stencil_dotProduct->weightData();
    for (std::size_t n = 0; n < nItems; ++n) {
        patternData_dotProduct[n] = patternData[n];
        weightData_dotProduct[n]  = ::dotProduct(weightData[n], vector);
    }
 
    stencil_dotProduct->setConstant(::dotProduct(stencil.getConstant(), vector));
}

template <typename V>
typename bitpit::DiscreteStencil<std::array<V, 3>> project(const typename bitpit::DiscreteStencil<std::array<V, 3>> &stencil, const std::array<V, 3> &direction)
{
    typename bitpit::DiscreteStencil<std::array<V, 3>> stencil_projection;
    project(stencil, direction, &stencil_projection);
 
    return stencil_projection;
}

template <typename V>
void project(const typename bitpit::DiscreteStencil<std::array<V, 3>> &stencil, const std::array<V, 3> &direction, typename bitpit::DiscreteStencil<std::array<V, 3>> *stencil_projection)
{
    stencil_projection->initialize(stencil);
    const std::size_t nItems = stencil_projection->size();
 
    typename bitpit::DiscreteStencil<std::array<V, 3>>::weight_type *weightData_projection = stencil_projection->weightData();
    for (std::size_t n = 0; n < nItems; ++n) {
        typename bitpit::DiscreteStencil<std::array<V, 3>>::weight_type &weight = weightData_projection[n];
        weight = ::dotProduct(weight, direction) * direction;
    }
 
    stencil_projection->setConstant(::dotProduct(stencil_projection->getConstant(), direction) * direction);
}
 
#endif