quadraturerulehelper.hh

Go to the documentation of this file.

// SPDX-FileCopyrightText: 2021-2025 The Ikarus Developers ikarus@ibb.uni-stuttgart.de
// SPDX-License-Identifier: LGPL-3.0-or-later
 
#pragma once
 
#include <dune/geometry/quadraturerules.hh>
namespace Ikarus {
 
template <class BaseQuadrature, class Quadrature>
auto tensorProductQuadrature(const BaseQuadrature& baseQuad, const Quadrature& onedQuad) {
  constexpr int baseQuadDim       = BaseQuadrature::d;
  auto rule                       = Dune::QuadratureRule<double, baseQuadDim + 1>();
  const unsigned int baseQuadSize = baseQuad.size();
  for (unsigned int bqi = 0; bqi < baseQuadSize; ++bqi) {
    const typename Dune::QuadraturePoint<double, baseQuadDim>::Vector& basePoint = baseQuad[bqi].position();
    const double& baseWeight                                                     = baseQuad[bqi].weight();
 
    typename Dune::QuadraturePoint<double, baseQuadDim + 1>::Vector point;
    for (unsigned int i = 0; i < baseQuadDim; ++i)
      point[i] = basePoint[i];
 
    const unsigned int onedQuadSize = onedQuad.size();
    for (unsigned int oqi = 0; oqi < onedQuadSize; ++oqi) {
      point[baseQuadDim] = onedQuad[oqi].position()[0];
      rule.emplace_back(Dune::QuadraturePoint(point, baseWeight * onedQuad[oqi].weight()));
    }
  }
  return rule;
}
 
template <int dim, typename Element>
auto defaultQuadratureRule(const Element& element, int order) {
  const auto defaultRule = Dune::QuadratureRules<double, dim>::rule(element.type(), order);
  if constexpr (requires { element.impl().getQuadratureRule(order); })
    if (element.impl().isTrimmed())
      return element.impl().getQuadratureRule(order);
    else
      return defaultRule;
  else
    return defaultRule;
}
 
int numberOfGaussPointsToOrder(int nGP) { return 2 * nGP - 1; }
} // namespace Ikarus