Merge pull request #15 from aurora-multiphysics/AutoToleranceFinder

Merge in functionality to automatically set mesh-merging tolerance

.github/workflows/main.yml

@@ -14,7 +14,7 @@ jobs:
     steps:
       # First check out the repository to get the docker file
       - name: Checkout
-        uses: actions/checkout@v2
+        uses: actions/checkout@master
       # Now build in a container with all deps
       - name:  DockerBuildTestPush
         run: docker build --build-arg BUILD_GIT_SHA=$GITHUB_HEAD_REF -t ci-vacmesh-ubuntu docker/skinning-ubuntu/ 

.gitignore

@@ -14,7 +14,6 @@ build/
 /test/build/
 CMakeFiles/
 Meshes/
-test/testingMeshes/
 *.face
 *.node
 *.step

VacuumMeshing/include/BoundaryGeneration/BoundaryGenerator.hpp

@@ -4,7 +4,7 @@
  * to doing exactly that.
  */
 #pragma once
-#include "Utils/libmeshConversions.hpp"
+#include "Utils/utils.hpp"
 #include "igl/triangle/triangulate.h"
 #include "libmesh/boundary_info.h"
 #include "libmesh/elem.h"
@@ -26,7 +26,7 @@ class BoundaryGenerator {
    */
   BoundaryGenerator(libMesh::Mesh &mesh, libMesh::Mesh &surface_mesh,
                     libMesh::Mesh &boundary_mesh,
-                    const double &merge_tolerance = 1e-08);
+                    double mesh_merge_tolerance = 0);
 
   /** Adds a boundary to the \c skinnedMesh. This new mesh with the boundary
  added is stored in the \c boundaryMesh. The \c length and number of mesh \c
@@ -67,13 +67,17 @@ class BoundaryGenerator {
   void checkBoundary(const double &length) const;
 
 protected:
+  // Method used to get tolerances automatically
+  void setMergeToleranceAuto();
+
   // Libmesh meshes to store the meshes we need
   libMesh::Mesh &mesh_, &surface_mesh_, &boundary_mesh_;
 
   // Eigen data structures to store mesh data for our boundaryz
   // Eigen::MatrixXd boundary_verts;
   // Eigen::MatrixXi boundary_elems;
-  double merge_tolerance_;
+  double mesh_merge_tolerance_;
+  double boundary_face_merge_tolerance_ = 1e-06;
 
 private:
 };

VacuumMeshing/include/BoundaryGeneration/CoilBoundaryGenerator.hpp

@@ -8,7 +8,7 @@ class CoilBoundaryGenerator : public BoundaryGenerator {
                         libMesh::Mesh &boundary_mesh,
                         const int sideset_one_id = 1,
                         const int sideset_two_id = 2,
-                        const double &merge_tolerance = 1e-06);
+                        double merge_tolerance = 0);
 
   ~CoilBoundaryGenerator();
 

VacuumMeshing/include/SurfaceMeshing/SurfaceGenerator.hpp

@@ -18,6 +18,11 @@
 
 class SurfaceMeshGenerator {
 public:
+  /**
+   * Enum to define what type of neighbor matching the user wants to use.
+   * */
+  enum NEIGHBOURTYPE { FACE, EDGE, VERTEX };
+
   /**
    * Default constructor
    */
@@ -49,6 +54,22 @@ class SurfaceMeshGenerator {
   void getSurface(bool writeMesh = false,
                   std::string outputFilename = "surface_mesh.e");
 
+  /**
+   * Return const reference to the surface face map
+   */
+  const std::multimap<unsigned int, unsigned int> &getSurfaceMap() const {
+    return surface_face_map;
+  }
+
+  /**
+   * Method that organises data about a face of an element. This data inc
+   */
+  void getFaceInfo(
+      libMesh::Elem *elem, int &face_id, std::vector<int> &original_node_ids,
+      std::vector<int> &connectivity,
+      std::map<int, std::vector<libMesh::boundary_id_type>> &boundary_data);
+
+protected:
   /** Method for checking whether an element has sides which should be in the
    * skin. Looks at the sides (faces or edges, depends if 2D or 3D element) of
    * \p element, and checks whether they are null. If they are, this side of \p
@@ -74,7 +95,8 @@ class SurfaceMeshGenerator {
   /** Method for grouping a discontinuous mesh into its continuous chunks.
    */
   void groupElems(libMesh::Mesh mesh,
-                  std::vector<std::vector<libMesh::dof_id_type>> &groups);
+                  std::vector<std::vector<libMesh::dof_id_type>> &groups,
+                  NEIGHBOURTYPE neighbour_type = NEIGHBOURTYPE::FACE);
 
   /** Method for checking whether an element has sides which should be in the
    * skin.
@@ -83,20 +105,10 @@ class SurfaceMeshGenerator {
                         std::vector<std::vector<libMesh::dof_id_type>> &groups,
                         std::string componentFilename);
 
-  /**
-   * Method that organises data about a face of an element. This data inc
-   */
-  void getFaceInfo(
-      libMesh::Elem *elem, int &face_id, std::vector<int> &original_node_ids,
-      std::vector<int> &connectivity,
-      std::map<int, std::vector<libMesh::boundary_id_type>> &boundary_data);
-
   /**
    *
    */
   std::multimap<unsigned int, unsigned int> surface_face_map;
-
-protected:
   // Mesh references
   libMesh::Mesh &mesh, &surfaceMesh;
 

VacuumMeshing/include/Utils/divConq.hpp

@@ -0,0 +1,88 @@
+/** \class
+ * Once a vacuum region mesh has been generated, it needs to be recombined,
+ * with the mesh of the original part to create the full mesh. This will results
+ * in lots of duplicate nodes. We can store nodes in an rTree, and when we want
+ * to add nodes we query the tree to see if an identical node exists. To decide
+ * whether the node is a duplicate or not, a tolerance value is used. So if
+ * nodes are within the euclidian distance described by "tolerance", the rTree
+ * returns a "hit" to show that an identical node already exists.
+ *
+ * The functionality in this class is used to help identify a suitable
+ * tolerance. A divide and conquer algorithm is used to identify the closest
+ * pair of nodes in a mesh. Using this value, we can identify a suitable choice
+ * of tolerance for the rTree. To learn more about the divide and conquer check
+ * out this paper by Shamos and Bentley,
+ * http://euro.ecom.cmu.edu/people/faculty/mshamos/1976ShamosBentley.pdf
+ */
+
+#include "libmesh/libmesh.h"
+#include "libmesh/mesh.h"
+#include "libmesh/node.h"
+#include "stdlib.h"
+#include "string"
+
+namespace VacuumMesher {
+/** Class used to find the closest pair of nodes.*/
+
+enum AXIS { X_AXIS, Y_AXIS, Z_AXIS };
+
+class ClosestPairFinder {
+public:
+  // Default constructor
+  ClosestPairFinder();
+
+  // Constructor that should be used most of the time.
+  ClosestPairFinder(libMesh::Mesh &mesh_);
+
+  // Return the euclidian distance between the two closest points
+  double closestPairMagnitude(int dim = 3);
+
+protected:
+  // Type def to prevent headache
+  typedef std::pair<libMesh::Node *, libMesh::Node *> nodePair;
+
+  //
+  double closestPair2D(std::vector<libMesh::Node *> &X);
+
+  /** Method to find the closest distance between any pair of points in a 3D
+   * mesh. Uses divide and conquer to perform the operation in O(nlog(n)) time
+   * (almost). 
+   *
+   * The only argument this function takes is a vector of libMesh node
+   * pointers, \p X. This vector should be pre-sorted based on the x coordinate of
+   * the nodes.
+   */
+  double closestPair3D(std::vector<libMesh::Node *> &X);
+
+  /** Method to return list of pairs of nodes that need their distance checking
+   * in the 3D method.*/
+  void getPotentialPairs(std::vector<libMesh::Node *> &X, double delta,
+                         std::vector<nodePair> &pairs);
+
+  /** Method to find the euclidian distance between \p n1 and \p n2 */
+  double eucDist(libMesh::Node *n1, libMesh::Node *n2);
+
+  /** Method that sorts \p vec by its x, y or z coordainte. \p axis = 0 will
+   * sort by x, 1 by y and 2 by z*/
+  void sortByIthCoord(AXIS axis, std::vector<libMesh::Node *> &vec);
+
+  /** Method to intialise class member \p xPoints. xPoints is a vector
+   * of libMesh node pointers, sorted by their x coordinate value, from low to
+   * high
+   */
+  void initialiseArrays();
+
+  /** Vector of libmesh node pointers you can use for the closestPair methods.
+   * You might be wondering why the methods take an argument if they could just
+   * directly call this vector, as it is a class member. This is not done
+   * because the methods get called recursively, and need their own input vector
+   * in each recursive call. If they just all reference \p xPoints stuff would
+   * go wrong. This is here just out of convenience*/
+  std::vector<libMesh::Node *> xPoints;
+
+  /** shared_ptr to libmesh mesh, of which we want to find the closest pair of
+   * nodes (points).*/
+  std::shared_ptr<libMesh::Mesh> mesh;
+};
+
+} // namespace VacuumMesher

VacuumMeshing/include/Utils/utils.hpp

@@ -0,0 +1,5 @@
+#pragma once
+#include "Utils/divConq.hpp"
+#include "Utils/getElemInfo.hpp"
+#include "Utils/libmeshConversions.hpp"
+#include "Utils/removeDupeNodes.hpp"

VacuumMeshing/include/VacuumGeneration/VacuumGenerator.hpp

@@ -11,9 +11,10 @@
 
 class VacuumGenerator {
 public:
-  VacuumGenerator(libMesh::Mesh &mesh, libMesh::Mesh &surface_mesh,
-                  libMesh::Mesh &boundary_mesh, libMesh::Mesh &vacuum_mesh,
-                  std::multimap<unsigned int, unsigned int> &surface_face_map);
+  VacuumGenerator(
+      libMesh::Mesh &mesh, libMesh::Mesh &surface_mesh,
+      libMesh::Mesh &boundary_mesh, libMesh::Mesh &vacuum_mesh,
+      const std::multimap<unsigned int, unsigned int> &surface_face_map);
 
   ~VacuumGenerator();
 

VacuumMeshing/src/BoundaryGeneration/BoundaryGenerator.cpp

@@ -6,9 +6,9 @@
 BoundaryGenerator::BoundaryGenerator(libMesh::Mesh &mesh,
                                      libMesh::Mesh &surface_mesh,
                                      libMesh::Mesh &boundary_mesh,
-                                     const double &merge_tolerance)
+                                     double mesh_merge_tolerance)
     : mesh_(mesh), surface_mesh_(surface_mesh), boundary_mesh_(boundary_mesh),
-      merge_tolerance_(merge_tolerance) {}
+      mesh_merge_tolerance_(mesh_merge_tolerance) {}
 
 BoundaryGenerator::~BoundaryGenerator() {}
 
@@ -28,8 +28,14 @@ void BoundaryGenerator::addBoundary(double length, int subdivisions,
 
   // Turn IGL mesh into libmesh Mesh
   IGLToLibMesh(boundary_mesh_, boundary_verts, boundary_elems);
+
+  // If unspecified, get merge tolerance
+  if (mesh_merge_tolerance_ == 0) {
+    setMergeToleranceAuto();
+  }
+
   // Combine IGL mesh with boundary
-  combineMeshes(1e-06, boundary_mesh_, surface_mesh_);
+  combineMeshes(mesh_merge_tolerance_, boundary_mesh_, surface_mesh_);
 }
 
 void BoundaryGenerator::genBoundary(Eigen::MatrixXd &tri_vertices,
@@ -97,8 +103,8 @@ void BoundaryGenerator::genBoundary(Eigen::MatrixXd &tri_vertices,
   // Combine the 6 square meshes to create the cubic boundary, using the rTree
   //  to avoid having any duplicate nodes
   for (int i = 0; i < 6; i++) {
-    combineMeshes(merge_tolerance_, tri_vertices, tri_elements, new_faces[i],
-                  square_elems_2);
+    combineMeshes(boundary_face_merge_tolerance_, tri_vertices, tri_elements,
+                  new_faces[i], square_elems_2);
   }
 }
 
@@ -203,4 +209,13 @@ void BoundaryGenerator::checkBoundary(const double &length) const {
         "The cubic boundary will overlap with the mesh. Please provide a "
         "larger boundary length.");
   }
+}
+
+void BoundaryGenerator::setMergeToleranceAuto() {
+  VacuumMesher::ClosestPairFinder closestPairBoundary(boundary_mesh_);
+  VacuumMesher::ClosestPairFinder closestPairSkin(surface_mesh_);
+
+  mesh_merge_tolerance_ = std::min(closestPairSkin.closestPairMagnitude(),
+                                   closestPairBoundary.closestPairMagnitude());
+  mesh_merge_tolerance_ /= 10;
 }

VacuumMeshing/src/BoundaryGeneration/CoilBoundaryGenerator.cpp

@@ -9,7 +9,7 @@ CoilBoundaryGenerator::CoilBoundaryGenerator(libMesh::Mesh &mesh,
                                              libMesh::Mesh &boundary_mesh,
                                              const int sideset_one_id,
                                              const int sideset_two_id,
-                                             const double &merge_tolerance)
+                                             double merge_tolerance)
     : BoundaryGenerator(mesh, surface_mesh, boundary_mesh, merge_tolerance),
       coil_sideset_one_id(sideset_one_id), coil_sideset_two_id(sideset_two_id) {
 }
@@ -23,9 +23,14 @@ void CoilBoundaryGenerator::addBoundary(const double length,
   // generate the boundary and store it in boundary_mesh_
   generateCoilBoundary(length, subdivisions, tri_flags);
 
+  if (mesh_merge_tolerance_ == 0) {
+    setMergeToleranceAuto();
+    std::cout << "Mesh merge tolerance used: " << mesh_merge_tolerance_
+              << std::endl;
+  }
   // Combine the boundary mesh with the surface mesh to create a mesh ready for
   // tetrahedralisation
-  combineMeshes(merge_tolerance_, boundary_mesh_, surface_mesh_);
+  combineMeshes(mesh_merge_tolerance_, boundary_mesh_, surface_mesh_);
 }
 
 void CoilBoundaryGenerator::generateCoilBoundary(const double length,
@@ -179,7 +184,8 @@ void CoilBoundaryGenerator::generateCoilFaceBound(
   igl::triangle::triangulate(verts, elems, holes, tri_args, newVerts, newElems);
   IGLToLibMesh(output_mesh, newVerts, newElems);
   //
-  combineMeshes(merge_tolerance_, output_mesh, remaining_boundary);
+  combineMeshes(boundary_face_merge_tolerance_, output_mesh,
+                remaining_boundary);
 }
 
 void CoilBoundaryGenerator::genSidesetMesh(libMesh::Mesh &mesh,
@@ -365,8 +371,8 @@ void CoilBoundaryGenerator::genRemainingFiveBoundaryFaces(
   translateMesh(new_faces[4], {0, 0, length});
 
   for (int i = 0; i < 5; i++) {
-    combineMeshes(merge_tolerance_, tri_vertices, tri_elems, new_faces[i],
-                  square_elems);
+    combineMeshes(boundary_face_merge_tolerance_, tri_vertices, tri_elems,
+                  new_faces[i], square_elems);
   }
 }