geometry.hpp

#pragma once
 
#include "mfem.hpp"
 
namespace serac {
 
template <int d>
struct Dimension {
  constexpr operator int() { return d; }
};
 
constexpr int num_quadrature_points(mfem::Geometry::Type g, int Q)
{
  if (g == mfem::Geometry::SEGMENT) {
    return Q;
  }
  if (g == mfem::Geometry::TRIANGLE) {
    return (Q * (Q + 1)) / 2;
  }
  if (g == mfem::Geometry::SQUARE) {
    return Q * Q;
  }
  if (g == mfem::Geometry::TETRAHEDRON) {
    return (Q * (Q + 1) * (Q + 2)) / 6;
  }
  if (g == mfem::Geometry::CUBE) {
    return Q * Q * Q;
  }
  return -1;
}
 
constexpr int dimension_of(mfem::Geometry::Type g)
{
  if (g == mfem::Geometry::SEGMENT) {
    return 1;
  }
 
  if (g == mfem::Geometry::TRIANGLE || g == mfem::Geometry::SQUARE) {
    return 2;
  }
 
  if (g == mfem::Geometry::TETRAHEDRON || g == mfem::Geometry::CUBE) {
    return 3;
  }
 
  return -1;
}
 
inline std::array<uint32_t, mfem::Geometry::NUM_GEOMETRIES> geometry_counts(const mfem::Mesh& mesh)
{
  std::array<uint32_t, mfem::Geometry::NUM_GEOMETRIES> counts{};
  for (int i = 0; i < mesh.GetNE(); i++) {
    counts[uint64_t(mesh.GetElementGeometry(i))]++;
  }
  return counts;
}
 
inline std::array<uint32_t, mfem::Geometry::NUM_GEOMETRIES> boundary_geometry_counts(const mfem::Mesh& mesh)
{
  std::array<uint32_t, mfem::Geometry::NUM_GEOMETRIES> counts{};
  for (int f = 0; f < mesh.GetNumFaces(); f++) {
    // skip interior faces
    if (mesh.GetFaceInformation(f).IsInterior()) continue;
 
    counts[uint64_t(mesh.GetFaceGeometry(f))]++;
  }
  return counts;
}
 
}  // namespace serac