FirstOrderPartialTest.hpp

#pragma once

#include <type_traits>

#include <amdis/GridFunctionOperator.hpp>
#include <amdis/common/StaticSize.hpp>
#include <amdis/utility/LocalToGlobalAdapter.hpp>

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


  class FirstOrderPartialTest
  {
    int comp_;

  public:
    FirstOrderPartialTest(tag::partialtest tag)
      : comp_(tag.comp)
    {}

    template <class CG, class Node, class Quad, class LocalFct, class Vec>
    void assemble(CG const& contextGeo, Node const& node, Quad const& quad,
                  LocalFct const& localFct, Vec& elementVector) const
    {
      static_assert(static_size_v<typename LocalFct::Range> == 1, "Expression must be of scalar type." );
      static_assert(Node::isLeaf, "Operator can be applied to Leaf-Nodes only.");

      LocalToGlobalBasisAdapter localBasis(node, contextGeo.elementGeometry());

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

        // The multiplicative factor in the integral transformation formula
        auto dx = contextGeo.integrationElement(qp.position()) * qp.weight();
        dx *= localFct(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];
        }
      }
    }
  };

  template <class LC>
  struct GridFunctionOperatorRegistry<tag::partialtest, LC>
  {
    static constexpr int degree = 1;
    using type = FirstOrderPartialTest;
  };

} // end namespace AMDiS