Oversized commit. Adds num_vertices. Implements MOAB ID<-->index mapping using intrinsic IDs

Author: pshriwise

include/xdg/libmesh/mesh_manager.h

@@ -80,6 +80,10 @@ class LibMeshManager : public MeshManager {
     }
   }
 
+  int num_vertices() const override {
+    return mesh()->n_nodes();
+  }
+
   int num_volume_elements(MeshID volume) const override {
     return get_volume_elements(volume).size();
   }

include/xdg/mesh_manager_interface.h

@@ -67,6 +67,8 @@ class MeshManager {
                const Position& u) const;
 
   // Mesh
+  virtual int num_vertices() const = 0;
+
   virtual int num_volume_elements(MeshID volume) const = 0;
 
   virtual int num_volume_elements() const;

include/xdg/moab/direct_access.h

@@ -47,35 +47,13 @@ class MBDirectAccess {
 
   //! \brief Determine the index of an element in the managed data
   inline size_t element_index(EntityHandle element) {
-    size_t element_index = 0;
-    auto fe = element_data_.first_elements.begin();
-    while(true) {
-      if (element - fe->first < fe->second) { break; }
-      element_index += fe->second;
-      fe++;
-      if (fe == element_data_.first_elements.end()) {
-        throw std::runtime_error("Element not found in MBDirectAccess::element_index");
-      }
-    }
-    element_index += element - fe->first;
-    return element_index;
+    return element_data_.entity_range.index(element);
   }
 
   //! \brief Determine the element handle of an element based on its index in
   //! the managed data
   inline EntityHandle element_handle(size_t index) {
-    size_t running_index = 0;
-    auto fe = element_data_.first_elements.begin();
-    while(true) {
-      if (index < running_index + fe->second) {
-        return fe->first + (index - running_index);
-      }
-      running_index += fe->second;
-      fe++;
-      if (fe == element_data_.first_elements.end()) {
-        throw std::runtime_error("Index out of range in MBDirectAccess::element_handle");
-      }
-    }
+    return element_data_.entity_range[index];
   }
 
   //! \brief Get the coordinates of a triangle as XDG Vertices
@@ -102,33 +80,11 @@ class MBDirectAccess {
   }
 
   int vertex_index(EntityHandle vertex) {
-    int vertex_index = 0;
-    auto fe = vertex_data_.first_vertices.begin();
-    while(true) {
-      if (vertex - fe->first < fe->second) { break; }
-      vertex_index += fe->second;
-      fe++;
-      if (fe == vertex_data_.first_vertices.end()) {
-        throw std::runtime_error("Vertex not found in MBDirectAccess::vertex_index");
-      }
-    }
-    vertex_index += vertex - fe->first;
-    return vertex_index;
+    return vertex_data_.vertex_range.index(vertex);
   }
 
   EntityHandle vertex_handle(int vertex_index) {
-    int running_index = 0;
-    auto fe = vertex_data_.first_vertices.begin();
-    while(true) {
-      if (vertex_index < running_index + fe->second) {
-        return fe->first + (vertex_index - running_index);
-      }
-      running_index += fe->second;
-      fe++;
-      if (fe == vertex_data_.first_vertices.end()) {
-        throw std::runtime_error("Index out of range in MBDirectAccess::vertex_handle");
-      }
-    }
+    return vertex_data_.vertex_range[vertex_index];
   }
 
   //! \brief Get the adjacent element
@@ -233,33 +189,33 @@ class MBDirectAccess {
     int element_stride {-1}; //!< Number of vertices used by each element
     std::vector<std::pair<EntityHandle, size_t>> first_elements; //!< Pairs of first element and length pairs for contiguous blocks of memory
     std::vector<const EntityHandle*> vconn; //!< Storage array(s) for the connectivity array
+    moab::Range entity_range; //!< Range of entities managed
 
     void setup(Interface * mbi) {
       ErrorCode rval;
 
       // setup face connectivity data
-      Range faces;
-      rval = mbi->get_entities_by_type(0, entity_type, faces, true);
-      MB_CHK_SET_ERR_CONT(rval, "Failed to get all elements of dimension 2 (faces)");
-      num_entities = faces.size();
+      rval = mbi->get_entities_by_type(0, entity_type, entity_range, true);
+      MB_CHK_SET_ERR_CONT(rval, "Failed to get all elements for the given entity type");
+      num_entities = entity_range.size();
 
       // only supporting triangle elements for now
-      if (!faces.all_of_type(entity_type)) { throw std::runtime_error("Not all 2D elements are triangles"); }
+      if (!entity_range.all_of_type(entity_type)) { throw std::runtime_error("Not all 2D elements are triangles"); }
 
-      moab::Range::iterator faces_it = faces.begin();
-      while(faces_it != faces.end()) {
+      moab::Range::iterator entity_it = entity_range.begin();
+      while(entity_it != entity_range.end()) {
         // set connectivity pointer, element stride and the number of elements
         EntityHandle* conntmp;
         int n_elements;
-        rval = mbi->connect_iterate(faces_it, faces.end(), conntmp, element_stride, n_elements);
+        rval = mbi->connect_iterate(entity_it, entity_range.end(), conntmp, element_stride, n_elements);
         MB_CHK_SET_ERR_CONT(rval, "Failed to get direct access to triangle elements");
 
         // set const pointers for the connectivity array and add first element/length pair to the set of first elements
         vconn.push_back(conntmp);
-        first_elements.push_back({*faces_it, n_elements});
+        first_elements.push_back({*entity_it, n_elements});
 
         // move iterator forward by the number of triangles in this contiguous memory block
-        faces_it += n_elements;
+        entity_it += n_elements;
       }
     }
 
@@ -295,19 +251,18 @@ class MBDirectAccess {
     void setup(Interface* mbi) {
       ErrorCode rval;
       // setup vertices
-      Range verts;
-      rval = mbi->get_entities_by_dimension(0, 0, verts, true);
+      rval = mbi->get_entities_by_dimension(0, 0, vertex_range, true);
       MB_CHK_SET_ERR_CONT(rval, "Failed to get all elements of dimension 0 (vertices)");
-      num_vertices = verts.size();
+      num_vertices = vertex_range.size();
 
-      moab::Range::iterator verts_it = verts.begin();
-      while (verts_it != verts.end()) {
+      moab::Range::iterator verts_it = vertex_range.begin();
+      while (verts_it != vertex_range.end()) {
         // set vertex coordinate pointers
         double* xtmp;
         double* ytmp;
         double* ztmp;
         int n_vertices;
-        rval = mbi->coords_iterate(verts_it, verts.end(), xtmp, ytmp, ztmp, n_vertices);
+        rval = mbi->coords_iterate(verts_it, vertex_range.end(), xtmp, ytmp, ztmp, n_vertices);
         MB_CHK_SET_ERR_CONT(rval, "Failed to get direct access to vertex elements");
 
         // add the vertex coordinate arrays to their corresponding vector of array pointers
@@ -351,6 +306,7 @@ class MBDirectAccess {
     std::vector<const double*> ty; //!< Storage array(s) for vertex y coordinates
     std::vector<const double*> tz; //!< Storage array(s) for vertex z coordinates
     std::vector<std::pair<EntityHandle, size_t>> first_vertices; //!< Pairs of first vertex and length pairs for contiguous blocks of memory
+    moab::Range vertex_range; //!< Range of entities managed
   };
 
   ConnectivityData face_data_;

include/xdg/moab/mesh_manager.h

@@ -54,6 +54,10 @@ class MOABMeshManager : public MeshManager {
   void add_surface_to_volume(MeshID volume, MeshID surface, Sense sense, bool overwrite=false) override;
 
   // Mesh
+  int num_vertices() const override {
+    return this->mb_direct()->n_vertices();
+  }
+
   int num_volume_elements(MeshID volume) const override;
 
   int num_volume_elements() const override;
@@ -83,7 +87,7 @@ class MOABMeshManager : public MeshManager {
   inline MeshID element_id(size_t element_idx) const override {
     // MOAB element IDs start at one and may not be contiguous,
     // so we need to map from index to ID
-    return tag_data<MeshID>(global_id_tag_, mb_direct()->element_handle(element_idx));
+    return moab_interface()->id_from_handle(mb_direct()->element_handle(element_idx));
   }
 
   inline int element_index(MeshID element) const override {
@@ -96,7 +100,7 @@ class MOABMeshManager : public MeshManager {
 
   inline MeshID vertex_id(size_t vertex_idx) const override {
     moab::EntityHandle vertex_handle = mb_direct()->vertex_handle(vertex_idx);
-    return tag_data<MeshID>(global_id_tag_, vertex_handle);
+    return moab_interface()->id_from_handle(vertex_handle);
   }
 
   inline int vertex_index(MeshID vertex) const override {

tests/mesh_mock.h

@@ -56,6 +56,10 @@ class MeshMock : public MeshManager {
     return -1;
   }
 
+  virtual int num_vertices() const override {
+    return vertices_.size();
+  }
+
   virtual int num_volume_elements(MeshID volume) const override {
     if (!volumetric_elements_) return 0;
     return 12;

tests/test_files

@@ -388,7 +388,8 @@ TEST_CASE("LibMesh Element ID and Index Mapping")
       REQUIRE(element_index == static_cast<int>(i));
     }
 
-    size_t num_vertices = static_cast<LibMeshManager*>(mesh_manager.get())->mesh()->n_nodes();
+    size_t num_vertices = mesh_manager->num_vertices();
+    REQUIRE(num_vertices == 2067);
     for (size_t i = 0; i < num_vertices; i++) {
       MeshID vertex_id = mesh_manager->vertex_id(i);
       int vertex_index = mesh_manager->vertex_index(vertex_id);
@@ -398,12 +399,12 @@ TEST_CASE("LibMesh Element ID and Index Mapping")
 
   // test mapping for non-contiguous IDs via manual modification
   {
-    // now create a new mesh manager
-    std::unique_ptr<LibMeshManager> mesh_manager2  {std::make_unique<LibMeshManager>()};
-    mesh_manager2->load_file("jezebel.exo");
+    // now create a new mesh manager (create explicitly so we can modify the mesh before init)
+    std::unique_ptr<LibMeshManager> mesh_manager  {std::make_unique<LibMeshManager>()};
+    mesh_manager->load_file("jezebel.exo");
 
     // tweak some of the element IDs to create gaps
-    auto* mesh = mesh_manager2->mesh();
+    auto* mesh = mesh_manager->mesh();
     int next_id = 0;
     std::vector<MeshID> modified_element_ids;
     for (auto* elem : mesh->active_element_ptr_range()) {
@@ -437,38 +438,38 @@ TEST_CASE("LibMesh Element ID and Index Mapping")
     mesh->allow_renumbering(false);
 
     // now initialize the mesh manager
-    mesh_manager2->init();
-    REQUIRE(mesh_manager2->num_volume_elements() == 10333);
+    mesh_manager->init();
+    REQUIRE(mesh_manager->num_volume_elements() == 10333);
 
       // check element ID to index mapping
     int expected_index = 0;
-    for (const auto* elem : mesh_manager2->mesh()->active_element_ptr_range()) {
+    for (const auto* elem : mesh_manager->mesh()->active_element_ptr_range()) {
       MeshID element_id = elem->id();
-      int element_index = mesh_manager2->element_index(element_id);
+      int element_index = mesh_manager->element_index(element_id);
       REQUIRE(element_index == expected_index);
       expected_index++;
     }
 
     // check node ID to index mapping
     expected_index = 0;
-    for (const auto* node : mesh_manager2->mesh()->node_ptr_range()) {
+    for (const auto* node : mesh_manager->mesh()->node_ptr_range()) {
       MeshID vertex_id = node->id();
-      int vertex_index = mesh_manager2->vertex_index(vertex_id);
+      int vertex_index = mesh_manager->vertex_index(vertex_id);
       REQUIRE(vertex_index == expected_index);
       expected_index++;
     }
 
     // check index to element ID mapping
-    size_t num_elements = mesh_manager2->num_volume_elements();
+    size_t num_elements = mesh_manager->num_volume_elements();
     for (size_t i = 0; i < num_elements; i++) {
-      MeshID element_id = mesh_manager2->element_id(i);
+      MeshID element_id = mesh_manager->element_id(i);
       REQUIRE(element_id == modified_element_ids[i]);
     }
 
     // check index to vertex ID mapping
-    size_t num_vertices = mesh_manager2->mesh()->n_nodes();
+    size_t num_vertices = mesh_manager->mesh()->n_nodes();
     for (size_t i = 0; i < num_vertices; i++) {
-      MeshID vertex_id = mesh_manager2->vertex_id(i);
+      MeshID vertex_id = mesh_manager->vertex_id(i);
       REQUIRE(vertex_id == modified_vertex_ids[i]);
     }
   }

tests/test_libmesh.cpp

@@ -282,20 +282,95 @@ TEST_CASE("TEST MOAB Find Element Method")
 
 TEST_CASE("MOAB Element ID and Index Mapping")
 {
-  std::shared_ptr<XDG> xdg = XDG::create(MeshLibrary::MOAB);
-  REQUIRE(xdg->mesh_manager()->mesh_library() == MeshLibrary::MOAB);
-  const auto& mesh_manager = xdg->mesh_manager();
-  mesh_manager->load_file("jezebel.h5m");
-  mesh_manager->init();
+  // test mapping for contiguous MOAB IDs using jezebel model
+  {
+    std::shared_ptr<XDG> xdg = XDG::create(MeshLibrary::MOAB);
+    REQUIRE(xdg->mesh_manager()->mesh_library() == MeshLibrary::MOAB);
+    const auto& mesh_manager = xdg->mesh_manager();
+    mesh_manager->load_file("jezebel.h5m");
+    mesh_manager->init();
+
+    size_t num_elements = mesh_manager->num_volume_elements();
+    REQUIRE(num_elements == 10333);
+
+    for (size_t idx = 0; idx < num_elements; ++idx) {
+      MeshID element_id = mesh_manager->element_id(idx);
+      // MOAB element IDs start at 1 and, for this model, are contiguous
+      REQUIRE(element_id == idx + 1);
+      int mapped_idx = mesh_manager->element_index(element_id);
+      REQUIRE(mapped_idx == static_cast<int>(idx));
+    }
 
-  size_t num_elements = mesh_manager->num_volume_elements();
-  REQUIRE(num_elements == 10333);
+    size_t num_vertices = mesh_manager->num_vertices();
+    REQUIRE(num_vertices == 2067);
+    for (size_t idx = 0; idx < num_vertices; ++idx) {
+      MeshID vertex_id = mesh_manager->vertex_id(idx);
+      // MOAB vertex IDs start at 1 and, for this model, are contiguous
+      REQUIRE(vertex_id == idx + 1);
+      int mapped_idx = mesh_manager->vertex_index(vertex_id);
+      REQUIRE(mapped_idx == static_cast<int>(idx));
+    }
+  }
+
+  // test mapping for non-contiguous MOAB IDs via modification
+  {
+    std::shared_ptr<MOABMeshManager> mesh_manager = std::make_shared<MOABMeshManager>();
+    mesh_manager->load_file("jezebel.h5m");
+
+    moab::Interface* mbi = mesh_manager->moab_interface();
+
+    // create gaps in the intrinsic ID space by deleteing some elements
+    moab::Range elem_range;
+    mbi->get_entities_by_type(0, moab::MBTET, elem_range);
+    int next_id = 0;
+    std::vector<MeshID> modified_ids;
+    for (const auto& elem : elem_range) {
+      if (next_id % 100 == 0) {
+        mbi->delete_entities(&elem, 1);
+        next_id++;
+        continue;
+      }
+      next_id ++;
+      modified_ids.push_back(mbi->id_from_handle(elem));
+    }
+
+    mesh_manager->init();
+
+    size_t num_elements = mesh_manager->num_volume_elements();
+    REQUIRE(num_elements == modified_ids.size());
+    for (size_t idx = 0; idx < num_elements; ++idx) {
+      MeshID expected_id = modified_ids[idx];
+      MeshID element_id = mesh_manager->element_id(idx);
+      REQUIRE(element_id == expected_id);
+      int mapped_idx = mesh_manager->element_index(element_id);
+      REQUIRE(mapped_idx == static_cast<int>(idx));
+    }
 
-  for (size_t idx = 0; idx < num_elements; ++idx) {
-    MeshID element_id = mesh_manager->element_id(idx);
-    // MOAB element IDs start at 1 and, for this model, are contiguous
-    REQUIRE(element_id == idx + 1);
-    int mapped_idx = mesh_manager->element_index(element_id);
-    REQUIRE(mapped_idx == static_cast<int>(idx));
+    //     // modify the vertex global IDs to be non-contiguous
+    // next_id = 0;
+    // std::vector<MeshID> modified_vertex_ids;
+    // moab::Range vertex_range;
+    // mbi->get_entities_by_type(0, moab::MBVERTEX, vertex_range);
+    // for (const auto& vertex : vertex_range) {
+    //   if (next_id % 50 == 0) {
+    //     mbi->delete_entities(&vertex, 1);
+    //     next_id++;
+    //     continue;
+    //   }
+    //   next_id++;
+    // }
+
+    // // reset the direct access manager to update vertex data
+    // mesh_manager->mb_direct()->update();
+
+    // size_t num_vertices = mesh_manager->num_vertices();
+    // REQUIRE(num_vertices == modified_vertex_ids.size());
+    // for (size_t idx = 0; idx < num_vertices; ++idx) {
+    //   MeshID expected_id = modified_vertex_ids[idx];
+    //   MeshID vertex_id = mesh_manager->vertex_id(idx);
+    //   REQUIRE(vertex_id == expected_id);
+    //   int mapped_idx = mesh_manager->vertex_index(vertex_id);
+    //   REQUIRE(mapped_idx == static_cast<int>(idx));
+    // }
   }
 }