DOFVectorView.inc.hpp

#pragma once

#include <amdis/common/FieldMatVec.hpp>
#include <amdis/LocalBasisCache.hpp>

namespace AMDiS {

template <class GB, class RT, class TP>
typename DOFVectorConstView<GB,RT,TP>::Range DOFVectorConstView<GB,RT,TP>::
LocalFunction::operator()(LocalDomain const& x) const
{
  assert( bound_ );
  auto y = Range(0);

  auto&& coefficients = *globalFunction_->dofVector_;
  auto&& nodeToRangeEntry = globalFunction_->nodeToRangeEntry_;
  forEachLeafNode_(*subTree_, [&,this](auto const& node, auto const& tp)
  {
    auto&& fe = node.finiteElement();
    auto&& localBasis = fe.localBasis();

    LocalBasisCache<std::decay_t<decltype(localBasis)>> localBasisCache(localBasis);
    auto const& shapeFunctionValues = localBasisCache.evaluateFunction(localView_.element(),x);

    // Get range entry associated to this node
    auto re = Dune::Functions::flatVectorView(nodeToRangeEntry(node, tp, y));

    for (std::size_t i = 0; i < localBasis.size(); ++i) {
      auto&& multiIndex = localView_.index(node.localIndex(i));

      // Get coefficient associated to i-th shape function
      auto c = Dune::Functions::flatVectorView(coefficients[multiIndex]);

      // Get value of i-th shape function
      auto v = Dune::Functions::flatVectorView(shapeFunctionValues[i]);

      std::size_t dimC = c.size();
      std::size_t dimV = v.size();
      assert(dimC*dimV == std::size_t(re.size()));
      for(std::size_t j = 0; j < dimC; ++j) {
        auto&& c_j = c[j];
        for(std::size_t k = 0; k < dimV; ++k)
          re[j*dimV + k] += c_j*v[k];
      }
    }
  });

  return y;
}


template <class GB, class RT, class TP>
typename DOFVectorConstView<GB,RT,TP>::DerivativeRange DOFVectorConstView<GB,RT,TP>::
GradientLocalFunction::operator()(LocalDomain const& x) const
{
  assert( bound_ );
  Range dy;
  for (std::size_t j = 0; j < dy.size(); ++j)
    dy[j] = 0;

  auto&& coefficients = *globalFunction_->dofVector_;
  auto&& nodeToRangeEntry = globalFunction_->nodeToRangeEntry_;
  forEachLeafNode_(*subTree_, [&,this](auto const& node, auto const& tp)
  {
    // TODO: may DOFVectorView::Range to FieldVector type if necessary
    using LocalDerivativeTraits = Dune::Functions::DefaultDerivativeTraits<Dune::FieldVector<double,1>(Domain)>;
    using GradientBlock = typename LocalDerivativeTraits::Range;

    auto&& fe = node.finiteElement();
    auto&& localBasis = fe.localBasis();

    LocalBasisCache<std::decay_t<decltype(localBasis)>> localBasisCache(localBasis);
    auto const& referenceGradients = localBasisCache.evaluateJacobian(localView_.element(),x);

    // The transposed inverse Jacobian of the map from the reference element to the element
    auto&& jacobian = geometry_.value().jacobianInverseTransposed(x);

    // Compute the shape function gradients on the real element
    std::vector<GradientBlock> gradients(referenceGradients.size());
    for (std::size_t i = 0; i < gradients.size(); ++i)
      multiplies_ABt(referenceGradients[i], jacobian, gradients[i]); // D[phi] * J^(-1) -> grad

    // Get range entry associated to this node
    auto re = Dune::Functions::flatVectorView(nodeToRangeEntry(node, tp, dy));

    for (std::size_t i = 0; i < localBasis.size(); ++i) {
      auto&& multiIndex = localView_.index(node.localIndex(i));

      // Get coefficient associated to i-th shape function
      auto c = Dune::Functions::flatVectorView(coefficients[multiIndex]);

      // Get value of i-th transformed reference gradient
      auto grad = Dune::Functions::flatVectorView(gradients[i]);

      std::size_t dimC = c.size();
      std::size_t dimV = grad.size();
      assert(dimC*dimV == std::size_t(re.size()));
      for(std::size_t j = 0; j < dimC; ++j) {
        auto&& c_j = c[j];
        for(std::size_t k = 0; k < dimV; ++k)
          re[j*dimV + k] += c_j*grad[k];
      }
    }
  });

  return dy;
}

} // end namespace AMDiS