ZeroOrderTestTrialvec.hpp

#pragma once

#include <type_traits>

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

namespace AMDiS
{
  namespace tag
  {
    struct test_trialvec {};
  }


  class ZeroOrderTestTrialvec
  {
  public:
    ZeroOrderTestTrialvec(tag::test_trialvec) {}

    template <class CG, class RN, class CN, class Quad, class LocalFct, class Mat>
    void assemble(CG const& contextGeo, RN const& rowNode, CN const& colNode,
                  Quad const& quad, LocalFct const& localFct, Mat& elementMatrix) const
    {
      static_assert(static_size_v<typename LocalFct::Range> == CG::dow,
        "Expression must be of vector type." );
      static_assert(RN::isLeaf && CN::isPower,
        "RN must be Leaf-Node and CN must be a Power-Node.");

      assert(colNode.degree() == CG::dow);

      std::size_t rowSize = rowNode.size();
      std::size_t colSize = colNode.child(0).size();

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

        // The multiplicative factor in the integral transformation formula
        const auto factor = contextGeo.integrationElement(qp.position()) * qp.weight();
        const auto b = localFct(local);

        auto const& rowShapeValues = rowNode.localBasisValuesAt(local);
        auto const& colShapeValues = colNode.child(0).localBasisValuesAt(local);

        for (std::size_t i = 0; i < rowSize; ++i) {
          const auto local_i = rowNode.localIndex(i);

          for (std::size_t j = 0; j < colSize; ++j) {
            const auto value = b * (factor * rowShapeValues[i] * colShapeValues[j]);

            for (std::size_t k = 0; k < colNode.degree(); ++k) {
              const auto local_kj = colNode.child(k).localIndex(j);
              elementMatrix[local_i][local_kj] += at(value,k);
            }
          }
        }
      }
    }
  };

  template <class LC>
  struct GridFunctionOperatorRegistry<tag::test_trialvec, LC>
  {
    static constexpr int degree = 0;
    using type = ZeroOrderTestTrialvec;
  };

} // end namespace AMDiS