FirstOrderPartialTest.hpp

#pragma once

#include <type_traits>

#include <amdis/localoperators/FirstOrderTestPartialTrial.hpp>

namespace AMDiS
{
  namespace tag
  {
    struct partialtest
    {
      std::size_t comp;
    };
  }



  template <class LC, class GridFct>
  class GridFunctionOperator<tag::partialtest, LC, GridFct>
      : public GridFunctionOperatorBase<GridFunctionOperator<tag::partialtest, LC, GridFct>, LC, GridFct>
  {
    using Self = GridFunctionOperator;
    using Super = GridFunctionOperatorBase<Self, LC, GridFct>;

    static_assert( static_size_v<typename GridFct::Range> == 1, "Expression must be of scalar type." );

  public:
    GridFunctionOperator(tag::partialtest tag, GridFct const& expr)
      : Super(expr, 1)
      , comp_(tag.comp)
    {}

    template <class CG, class Node, class Vec>
    void getElementVector(CG const& contextGeo, Node const& node, Vec& elementVector)
    {
      static_assert(Node::isLeaf, "Operator can be applied to Leaf-Nodes only.");

      auto const& quad = this->getQuadratureRule(contextGeo.type(), node);
      LocalToGlobalBasisAdapter localBasis(node, contextGeo.geometry());

      for (auto const& qp : quad) {
        // Position of the current quadrature point in the reference element
        decltype(auto) local = contextGeo.local(qp.position());

        // The multiplicative factor in the integral transformation formula
        auto dx = contextGeo.integrationElement(qp.position()) * qp.weight();
        dx *= Super::coefficient(local);

        // Compute the shape function gradients on the real element
        auto const& partial = localBasis.partialsAt(local, comp_);

        for (std::size_t j = 0; j < localBasis.size(); ++j) {
          const auto local_j = node.localIndex(j);
          elementVector[local_j] += dx * partial[j];
        }
      }
    }

  private:
    std::size_t comp_;
  };

} // end namespace AMDiS