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