mesh_hierarchy.h

#ifndef INCG_LF_REFINEMENT_HIER_H_
#define INCG_LF_REFINEMENT_HIER_H_
 
#include <iostream>
 
#include "hybrid2d_refinement_pattern.h"
 
namespace lf::refinement {
 
struct PointChildInfo {
  explicit PointChildInfo() = default;
  RefPat ref_pat{RefPat::rp_nil};
  glb_idx_t child_point_idx{idx_nil};
};
 
struct EdgeChildInfo {
  explicit EdgeChildInfo() = default;
  RefPat ref_pat_{RefPat::rp_nil};
  std::vector<glb_idx_t> child_edge_idx;
  std::vector<glb_idx_t> child_point_idx;
};
 
struct CellChildInfo {
  explicit CellChildInfo() = default;
  RefPat ref_pat_{RefPat::rp_nil};
  sub_idx_t anchor_{idx_nil};
  std::vector<glb_idx_t> child_cell_idx;
  std::vector<glb_idx_t> child_edge_idx;
  std::vector<glb_idx_t> child_point_idx;
};
 
struct ParentInfo {
  explicit ParentInfo() = default;
  // Data members
  const mesh::Entity *parent_ptr{
      nullptr}; 
  glb_idx_t parent_index{
      idx_nil}; 
  sub_idx_t child_number{idx_nil}; 
  std::unique_ptr<const lf::geometry::Geometry> rel_ref_geo_{nullptr};
};
 
class MeshHierarchy {
 public:
  MeshHierarchy(const std::shared_ptr<mesh::Mesh> &base_mesh,
                std::unique_ptr<mesh::MeshFactory> mesh_factory);
  MeshHierarchy(const MeshHierarchy &) = delete;
  MeshHierarchy &operator=(const MeshHierarchy &) = delete;
  MeshHierarchy(MeshHierarchy &&) = delete;
  MeshHierarchy &operator=(MeshHierarchy &&) = delete;
 
  [[nodiscard]] size_type NumLevels() const { return meshes_.size(); }
 
  [[nodiscard]] std::shared_ptr<const mesh::Mesh> getMesh(
      size_type level) const {
    LF_VERIFY_MSG(level < meshes_.size(),
                  "Level " << level << " outside scope");
    return meshes_.at(level);
  }
  [[nodiscard]] std::shared_ptr<mesh::Mesh> getMesh(size_type level) {
    LF_VERIFY_MSG(level < meshes_.size(),
                  "Level " << level << " outside scope");
    return meshes_.at(level);
  }
 
  [[nodiscard]] std::vector<std::shared_ptr<const mesh::Mesh>> getMeshes()
      const {
    return {meshes_.begin(), meshes_.end()};
  }
 
  [[nodiscard]] const std::vector<PointChildInfo> &PointChildInfos(
      size_type level) const {
    LF_VERIFY_MSG(level < NumLevels(), "Illegal level " << level);
    return point_child_infos_[level];
  }
  [[nodiscard]] const std::vector<EdgeChildInfo> &EdgeChildInfos(
      size_type level) const {
    LF_VERIFY_MSG(level < NumLevels(), "Illegal level " << level);
    return edge_child_infos_[level];
  }
  [[nodiscard]] const std::vector<CellChildInfo> &CellChildInfos(
      size_type level) const {
    LF_VERIFY_MSG(level < NumLevels(), "Illegal level " << level);
    return cell_child_infos_[level];
  }
  [[nodiscard]] const std::vector<ParentInfo> &ParentInfos(size_type level,
                                                           dim_t codim) const {
    LF_VERIFY_MSG(level < NumLevels(), "Illegal level " << level);
    LF_VERIFY_MSG(codim < 3, "Codim = " << codim << " illegal");
    return parent_infos_[level][codim];
  }
  [[nodiscard]] const std::vector<sub_idx_t> &RefinementEdges(
      size_type level) const {
    LF_VERIFY_MSG(level < NumLevels(), "Illegal level " << level);
    return refinement_edges_[level];
  }
 
  void RefineRegular(RefPat ref_pat = RefPat::rp_regular);
  template <typename Marker>
  void MarkEdges(Marker &&marker);
 
  void RefineMarked();
  void Coarsen();
  [[nodiscard]] const lf::geometry::Geometry *GeometryInParent(
      size_type level, const lf::mesh::Entity &e) const;
 
  [[nodiscard]] const lf::mesh::Entity *ParentEntity(
      size_type level, const lf::mesh::Entity &e) const;
 
  std::ostream &PrintInfo(std::ostream &o, unsigned ctrl = 0) const;
 
  virtual ~MeshHierarchy() = default;
 
 private:
  void PerformRefinement();
 
  void initGeometryInParent();
 
  std::vector<std::shared_ptr<mesh::Mesh>> meshes_;
  std::unique_ptr<mesh::MeshFactory> mesh_factory_;
  std::vector<std::vector<PointChildInfo>> point_child_infos_;
  std::vector<std::vector<EdgeChildInfo>> edge_child_infos_;
  std::vector<std::vector<CellChildInfo>> cell_child_infos_;
  std::vector<std::array<std::vector<ParentInfo>, 3>> parent_infos_;
  std::vector<std::vector<bool>> edge_marked_;
  std::vector<std::vector<sub_idx_t>> refinement_edges_;
 
  [[nodiscard]] static sub_idx_t LongestEdge(const lf::mesh::Entity &T);
 
 public:
  static std::shared_ptr<spdlog::logger> &Logger();
};
 
template <typename Marker>
void MeshHierarchy::MarkEdges(Marker &&marker) {
  // Retrieve the finest mesh in the hierarchy
  const mesh::Mesh &finest_mesh(*meshes_.back());
 
  LF_VERIFY_MSG(edge_marked_.back().size() == finest_mesh.NumEntities(1),
                "Length  mismatch for edge flag array");
 
  // Run through the edges = entities of co-dimension 1
  for (const mesh::Entity *edge : finest_mesh.Entities(1)) {
    const glb_idx_t edge_index = finest_mesh.Index(*edge);
    (edge_marked_.back())[edge_index] = marker(finest_mesh, *edge);
  }
}  // end MeshHierarchy::MarkEdges
 
std::shared_ptr<MeshHierarchy> GenerateMeshHierarchyByUniformRefinemnt(
    const std::shared_ptr<lf::mesh::Mesh> &mesh_p, lf::base::size_type ref_lev,
    RefPat ref_pat = RefPat::rp_regular);
 
template <typename POSPRED>
class EntityCenterPositionSelector {
 public:
  EntityCenterPositionSelector(const EntityCenterPositionSelector &) = default;
  EntityCenterPositionSelector(EntityCenterPositionSelector &&) noexcept =
      default;
  EntityCenterPositionSelector &operator=(
      const EntityCenterPositionSelector &) = default;
  EntityCenterPositionSelector &operator=(
      EntityCenterPositionSelector &&) noexcept = default;
  explicit EntityCenterPositionSelector(POSPRED pos_pred)
      : pos_pred_(std::move(pos_pred)) {}
  bool operator()(const lf::mesh::Entity &ent) const {
    const lf::base::RefEl ref_el_type = ent.RefEl();
    // Obtain shape of entity
    const lf::geometry::Geometry *geo_ptr = ent.Geometry();
    LF_ASSERT_MSG(geo_ptr != nullptr, "Missing geometry for " << ent);
    switch (ref_el_type) {
      case lf::base::RefEl::kPoint(): {
        const Eigen::MatrixXd pos(geo_ptr->Global(kpoint_center_));
        return pos_pred_(pos);
      }
      case lf::base::RefEl::kSegment(): {
        const Eigen::MatrixXd pos(geo_ptr->Global(kedge_center_));
        return pos_pred_(pos);
      }
      case lf::base::RefEl::kTria(): {
        const Eigen::MatrixXd pos(geo_ptr->Global(ktria_center_));
        return pos_pred_(pos);
      }
      case lf::base::RefEl::kQuad(): {
        const Eigen::MatrixXd pos(geo_ptr->Global(kquad_center_));
        return pos_pred_(pos);
      }
      default: {
        LF_ASSERT_MSG(false, "Illegal entity type");
        break;
      }
    }  // end switch
    return false;
  }
 
  virtual ~EntityCenterPositionSelector() = default;
 
 private:
  POSPRED pos_pred_;
 
  static const Eigen::MatrixXd kpoint_center_;
  static const Eigen::MatrixXd kedge_center_;
  static const Eigen::MatrixXd ktria_center_;
  static const Eigen::MatrixXd kquad_center_;
};
 
template <typename POSPRED>
const Eigen::MatrixXd EntityCenterPositionSelector<POSPRED>::kpoint_center_ =
    Eigen::MatrixXd::Zero(0, 1);
 
template <typename POSPRED>
const Eigen::MatrixXd EntityCenterPositionSelector<POSPRED>::kedge_center_ =
    (Eigen::MatrixXd(1, 1) << 0.5).finished();
 
template <typename POSPRED>
const Eigen::MatrixXd EntityCenterPositionSelector<POSPRED>::ktria_center_ =
    (Eigen::MatrixXd(2, 1) << 0.33, 0.33).finished();
 
template <typename POSPRED>
const Eigen::MatrixXd EntityCenterPositionSelector<POSPRED>::kquad_center_ =
    (Eigen::MatrixXd(2, 1) << 0.5, 0.5).finished();
 
}  // namespace lf::refinement
 
#endif