ref_el.h

 
 
#ifndef __96e6ff0ee0034f4584fcdfc7e9c53f82
#define __96e6ff0ee0034f4584fcdfc7e9c53f82
 
#include <array>
#include <vector>
 
//#include <boost/range.hpp>
 
#include "base.h"
#include "lf_assert.h"
 
namespace lf::base {
 
enum class RefElType : unsigned char {
  kPoint = 1,
  kSegment = 2,
  kTria = 3,
  kQuad = 4,
};
 
namespace internal {
// Some utility methods that are needed to deduce compile time return types:
constexpr dim_t DimensionImpl(RefElType type) {
  switch (type) {
    case RefElType::kPoint:
      return 0;
    case RefElType::kSegment:
      return 1;
    case RefElType::kTria:
    case RefElType::kQuad:
      return 2;
    default:
      throw std::runtime_error(
          "RefEl::Dimension() not implemented for this RefEl type.");
  }
}
}  // namespace internal
 
class RefEl {
  // Node coordinates as dynamic eigen matrices:
  static const Eigen::MatrixXd ncoords_point_dynamic_;
  static const Eigen::MatrixXd ncoords_segment_dynamic_;
  static const Eigen::MatrixXd ncoords_tria_dynamic_;
  static const Eigen::MatrixXd ncoords_quad_dynamic_;
 
  // Node coordinates as fixed eigen matrices:
  static const std::vector<Eigen::Matrix<double, 0, 1>> ncoords_point_static_;
  static const std::vector<Eigen::Matrix<double, 1, 1>> ncoords_segment_static_;
  static const std::vector<Eigen::Vector2d> ncoords_tria_static_;
  static const std::vector<Eigen::Vector2d> ncoords_quad_static_;
 
  // Member variable
  RefElType type_;
 
  // subSubEntities, used by SubSubEntity2SubEntity
  static constexpr std::array<std::array<base::dim_t, 2>, 3>
      sub_sub_entity_index_tria_ = {{{0, 1}, {1, 2}, {2, 0}}};
  static constexpr std::array<std::array<base::dim_t, 2>, 4>
      sub_sub_entity_index_quad_ = {{{0, 1}, {1, 2}, {2, 3}, {3, 0}}};
 
 public:
  using dim_t = unsigned int;
  template <RefElType type>
  using NodeCoordVector =
      Eigen::Matrix<double, internal::DimensionImpl(type), 1>;
 
  static constexpr RefEl kPoint() { return RefEl(RefElType::kPoint); }
 
  static constexpr RefEl kSegment() { return RefEl(RefElType::kSegment); }
 
  static constexpr RefEl kTria() { return RefEl(RefElType::kTria); }
 
  static constexpr RefEl kQuad() { return RefEl(RefElType::kQuad); }
 
  explicit constexpr RefEl(RefElType type) noexcept : type_(type) {}
 
  constexpr RefEl(const RefEl&) = default;
 
  constexpr RefEl(RefEl&&) = default;
 
  // NOLINTNEXTLINE(modernize-use-equals-default,hicpp-use-equals-default,cert-oop54-cpp)
  constexpr RefEl& operator=(const RefEl& rhs) {
    type_ = rhs.type_;
    return *this;
  }
 
  constexpr RefEl& operator=(RefEl&& rhs) noexcept {
    type_ = rhs.type_;
    return *this;
  }
 
  [[nodiscard]] constexpr dim_t Dimension() const {
    return internal::DimensionImpl(type_);
  }
 
  [[nodiscard]] constexpr size_type NumNodes() const {
    switch (type_) {
      case RefElType::kPoint:
        return 1;
      case RefElType::kSegment:
        return 2;
      case RefElType::kTria:
        return 3;
      case RefElType::kQuad:
        return 4;
      default:
        throw std::runtime_error(
            "RefEl::NumNodes() not implemented for this RefEl type.");
    }
  }
 
  [[nodiscard]] const Eigen::MatrixXd& NodeCoords() const {
    switch (type_) {
      case RefElType::kPoint:
        return ncoords_point_dynamic_;
      case RefElType::kSegment:
        return ncoords_segment_dynamic_;
      case RefElType::kTria:
        return ncoords_tria_dynamic_;
      case RefElType::kQuad:
        return ncoords_quad_dynamic_;
      default:
        LF_VERIFY_MSG(
            false, "RefEl::NodeCoords() not implemented for this RefEl type.");
    }
  }
 
  template <RefElType type>
  static const std::vector<NodeCoordVector<type>>& NodeCoords() {
    // NOLINTNEXTLINE
    if constexpr (type == RefElType::kPoint) {
      return ncoords_point_static_;
    }
    // NOLINTNEXTLINE
    if constexpr (type == RefElType::kSegment) {
      return ncoords_segment_static_;
    }
    // NOLINTNEXTLINE
    if constexpr (type == RefElType::kTria) {
      return ncoords_tria_static_;
    }
    // NOLINTNEXTLINE
    if constexpr (type == RefElType::kQuad) {
      return ncoords_quad_static_;
    }
    LF_VERIFY_MSG(false,
                  "RefEl::NodeCoords<>() not implemented for this RefEl type.");
  }
 
  [[nodiscard]] constexpr size_type NumSubEntities(dim_t sub_codim) const {
    LF_ASSERT_MSG_CONSTEXPR(sub_codim >= 0, "sub_codim is negative");
    LF_ASSERT_MSG_CONSTEXPR(sub_codim <= Dimension(),
                            "sub_codim > Dimension()");
    if (sub_codim == 0) {
      return 1;
    }
    switch (type_) {
      case RefElType::kSegment:
        return 2;  // sub_codim=1
      case RefElType::kTria:
        return 3;  // sub_codim=1,2
      case RefElType::kQuad:
        return 4;  // sub_codim=1,2
      default:
        LF_ASSERT_MSG_CONSTEXPR(
            false,
            "RefEl::NumSubEntities() not implemented for this RefElType.");
    }
    return 0;  // prevent warnings from compilers
  }
 
  [[nodiscard]] constexpr RefEl SubType(dim_t sub_codim,
                                        dim_t sub_index) const {
    LF_ASSERT_MSG_CONSTEXPR(sub_codim >= 0, "sub_codim is negative");
    LF_ASSERT_MSG_CONSTEXPR(sub_codim <= Dimension(),
                            "sub_codim > Dimension()");
    LF_ASSERT_MSG_CONSTEXPR(sub_index >= 0, "sub_index is negative");
    LF_ASSERT_MSG_CONSTEXPR(sub_index < NumSubEntities(sub_codim),
                            "sub_index >= NumSubEntities");
 
    if (sub_codim == 0) {
      return *this;
    }
    if (sub_codim == Dimension()) {
      return kPoint();
    }
    if (Dimension() - sub_codim == 1) {
      return kSegment();
    } else {  // NOLINT(readability-else-after-return)
      LF_ASSERT_MSG_CONSTEXPR(false, "This code should never be reached.");
    }
 
    return kPoint();  // prevent warnings from compiler
  }
 
  [[nodiscard]] constexpr sub_idx_t SubSubEntity2SubEntity(
      dim_t sub_codim, sub_idx_t sub_index, dim_t sub_rel_codim,
      sub_idx_t sub_rel_index) const {
    LF_ASSERT_MSG_CONSTEXPR(sub_codim >= 0, "sub_codim negative");
    LF_ASSERT_MSG_CONSTEXPR(sub_codim <= Dimension(), "sub_codim > Dimension");
    LF_ASSERT_MSG_CONSTEXPR(sub_index >= 0, "sub_index negative");
    LF_ASSERT_MSG_CONSTEXPR(sub_index <= NumSubEntities(sub_codim),
                            "sub_index >= NumSubEntities");
    LF_ASSERT_MSG_CONSTEXPR(sub_rel_codim >= 0, "sub_rel_codim negative.");
    LF_ASSERT_MSG_CONSTEXPR(sub_rel_codim <= Dimension() - sub_codim,
                            "subSubCodim out of bounds.");
    LF_ASSERT_MSG_CONSTEXPR(sub_rel_index >= 0, "sub_rel_index negative.");
    LF_ASSERT_MSG_CONSTEXPR(
        sub_rel_index <
            SubType(sub_codim, sub_index).NumSubEntities(sub_rel_codim),
        "sub_sub_index out of bounds.");
 
    if (type_ == RefElType::kPoint) {
      return 0;
    }
    if (sub_codim == 0) {
      return sub_rel_index;
    }
    if (sub_codim == Dimension()) {
      return sub_index;
    }
 
    // from here on, it must be a segment
    switch (type_) {
      case RefElType::kTria:
        return sub_sub_entity_index_tria_[sub_index][sub_rel_index];
      case RefElType::kQuad:
        return sub_sub_entity_index_quad_[sub_index][sub_rel_index];
      default:
        LF_ASSERT_MSG_CONSTEXPR(false, "This code should never be reached.");
    }
 
    return 0;  // Prevent warnings from compiler...
  }
 
  [[nodiscard]] std::string ToString() const {
    switch (type_) {
      case RefElType::kPoint:
        return "NODE";
      case RefElType::kSegment:
        return "EDGE";
      case RefElType::kTria:
        return "TRIA";
      case RefElType::kQuad:
        return "QUAD";
      default:
        LF_VERIFY_MSG(false, "ToString() not implemented for this RefElType");
    }
  }
 
  // NOLINTNEXTLINE(google-explicit-constructor, hicpp-explicit-conversions)
  [[nodiscard]] constexpr operator RefElType() const { return type_; }
 
  // NOLINTNEXTLINE(google-explicit-constructor, hicpp-explicit-conversions)
  [[nodiscard]] constexpr unsigned int Id() const {
    return static_cast<unsigned int>(type_);
  }
 
  ~RefEl() = default;
 
};  // class RefEl
 
// Declare print function
void PrintInfo(std::ostream& o, const RefEl& ref_el, int output_ctrl = 0);
 
inline std::ostream& operator<<(std::ostream& stream, const RefEl& ref_el) {
  return stream << ref_el.ToString();
}
 
}  // namespace lf::base
 
#endif  // __96e6ff0ee0034f4584fcdfc7e9c53f82