quadrature.hpp

Go to the documentation of this file.

// Copyright (c) 2019-2024, Lawrence Livermore National Security, LLC and
// other Serac Project Developers. See the top-level LICENSE file for
// details.
//
// SPDX-License-Identifier: (BSD-3-Clause)
 
#pragma once
 
#include "tensor.hpp"
#include "polynomials.hpp"
#include "finite_element.hpp"
 
namespace serac {
 
template <int n, int dim>
struct QuadratureRule {
  tensor<double, n> weights;
 
  tensor<double, n, dim> points;
 
  SERAC_HOST_DEVICE constexpr std::size_t size() const { return n; }
};
 
template <mfem::Geometry::Type g, int Q>
SERAC_HOST_DEVICE constexpr auto GaussQuadratureRule()
{
  auto x = GaussLegendreNodes<Q, mfem::Geometry::SEGMENT>();
  auto w = GaussLegendreWeights<Q, mfem::Geometry::SEGMENT>();
 
  if constexpr (g == mfem::Geometry::SEGMENT) {
    return QuadratureRule<Q, 1>{w, make_tensor<Q, 1>([&x](int i, int /*j*/) { return x[i]; })};
  }
 
  if constexpr (g == mfem::Geometry::SQUARE) {
    QuadratureRule<Q * Q, 2> rule{};
    int                      count = 0;
    for (int j = 0; j < Q; j++) {
      for (int i = 0; i < Q; i++) {
        rule.points[count]    = {x[i], x[j]};
        rule.weights[count++] = w[i] * w[j];
      }
    }
    return rule;
  }
 
  if constexpr (g == mfem::Geometry::CUBE) {
    QuadratureRule<Q * Q * Q, 3> rule{};
    int                          count = 0;
    for (int k = 0; k < Q; k++) {
      for (int j = 0; j < Q; j++) {
        for (int i = 0; i < Q; i++) {
          rule.points[count]    = {x[i], x[j], x[k]};
          rule.weights[count++] = w[i] * w[j] * w[k];
        }
      }
    }
    return rule;
  }
}
 
}  // namespace serac