detect_slit.cc

// This utility code can be used to detect slits in a Mesh.  A "slit"
// is a zero-width gap between two elements which should actually be
// neighbors.  This can be caused by having duplicate nodes in a mesh
// generated by Cubit.
//
// The following journal file should work in any recent version of Cubit
// to generate a mesh with a bunch of duplicated nodes along an interface.
// This is caused by the user forgetting to call "merge all" after setting
// up his geometry.

// reset
// brick x 1 y 1 z 1
// brick x 1 y 1 z 1
// move body 2 x 1
// imprint all # does this do anything? just prints "group imprint finished"
//
// # Leaving out this line creates the "slit" in the mesh, otherwise prints
// # "Consolidated 1 pair of surfaces".  If you leave this out, you
// # should have 2662 nodes rather than 2541!
// # merge all
//
// mesh surf 4
// mesh volume all
// block 1 volume 1
// block 2 volume 2
// export mesh "~/Desktop/test_slit.e" dimension 3 block all overwrite

#include "libmesh/libmesh.h"
#include "libmesh/mesh.h"
#include "libmesh/elem.h"
#include "libmesh/mesh_generation.h"
#include "libmesh/mesh_modification.h"
#ifdef LIBMESH_HAVE_NANOFLANN
#  include "libmesh/nanoflann.hpp"
#endif
#include "libmesh/meshfree_interpolation.h"
#include "libmesh/perf_log.h"
#include "libmesh/getpot.h"

// Bring in the libmesh namespace
using namespace libMesh;


// Nanoflann uses "duck typing" to allow users to define their own adaptors...
template <unsigned int Dim>
class NanoflannMeshAdaptor
{
private:
  // Constant reference to the Mesh we are adapting for use in Nanoflann
  const MeshBase & _mesh;

public:
  NanoflannMeshAdaptor (const MeshBase & mesh) :
    _mesh(mesh)
  {}

  /**
   * libMesh \p Point coordinate type
   */
  typedef Real coord_t;

  /**
   * Must return the number of data points
   */
  inline size_t kdtree_get_point_count() const { return _mesh.n_nodes(); }

  /**
   * Returns the distance between the vector "p1[0:size-1]"
   * and the data point with index "idx_p2" stored in _mesh
   */
  inline coord_t kdtree_distance(const coord_t * p1, const size_t idx_p2, size_t size) const
  {
    libmesh_assert_equal_to (size, Dim);
    libmesh_assert_less (idx_p2, _pts.size());

    // Construct a libmesh Point object from the input coord_t.  This
    // assumes LIBMESH_DIM==3.
    Point point1(p1[0],
                 size > 1 ? p1[1] : 0.,
                 size > 2 ? p1[2] : 0.);

    // Get the referred-to point from the Mesh
    const Point & point2 = _mesh.point(idx_p2);

    // Compute Euclidean distance, squared
    return (point1 - point2).norm_sq();
  }

  /**
   * Returns the dim'th component of the idx'th point in the class:
   * Since this is inlined and the "dim" argument is typically an immediate value, the
   *  "if's" are actually solved at compile time.
   */
  inline coord_t kdtree_get_pt(const size_t idx, int dim) const
  {
    libmesh_assert_less (dim, (int) Dim);
    libmesh_assert_less (idx, _mesh.n_nodes());
    libmesh_assert_less (dim, 3);

    return _mesh.point(idx)(dim);
  }

  /**
   * Optional bounding-box computation: return false to default to a standard bbox computation loop.
   * Return true if the BBOX was already computed by the class and returned in "bb" so it can be
   * avoided to redo it again. Look at bb.size() to find out the expected dimensionality
   * (e.g. 2 or 3 for point clouds)
   */
  template <class BBOX>
  bool kdtree_get_bbox(BBOX & /* bb */) const { return false; }
};



int main (int argc, char ** argv)
{
  LibMeshInit init(argc, argv);

  PerfLog perf_main("Main Program");

  // Create a GetPot object to parse the command line
  GetPot command_line (argc, argv);

  std::string input_filename = "";
  if ( command_line.search(2, "-i", "--input-file") )
    input_filename = command_line.next(input_filename);

  // The user must provide an input filename
  if (input_filename == "")
    libmesh_error_msg("Please run with -i or --input-filename to set the input mesh filename.");

  Mesh mesh(init.comm());
  mesh.read(input_filename);

  // Log how long it takes to do the duplicate node detection using Nanoflann
  {
    libMesh::out << "Starting Nanoflann tree construction and search algorithm." << std::endl;

    perf_main.push("Nanoflann");

    // Loop over nodes to try and detect duplicates.  We use nanoflann
    // for this, inspired by
    // contrib/nanoflann/examples/pointcloud_adaptor_example.cpp

    // Declare a type templated on NanoflannMeshAdaptor
    typedef nanoflann::L2_Simple_Adaptor<Real, NanoflannMeshAdaptor<3> > adatper_t;

    // Declare a KDTree type based on NanoflannMeshAdaptor
    typedef nanoflann::KDTreeSingleIndexAdaptor<adatper_t, NanoflannMeshAdaptor<3>, 3> kd_tree_t;

    // Build adaptor and tree objects
    NanoflannMeshAdaptor<3> mesh_adaptor(mesh);
    kd_tree_t kd_tree(3, mesh_adaptor, nanoflann::KDTreeSingleIndexAdaptorParams(/*max leaf=*/10));

    // Construct the tree
    kd_tree.buildIndex();

    // Count the number of duplicate nodes found
    unsigned dupe_count = 0;

    // For every node in the mesh, search the KDtree and find any
    // zero-distance neighbors with a different index from the current
    // node being searched... this indicates a duplicate node.
    for (const auto & node : mesh.node_ptr_range())
      {
        Real query_pt[3] = {(*node)(0), (*node)(1), (*node)(2)};

        // The number of results we want to get.  We'll look for the
        // two closest neighbors, since even 1 duplicate indicates
        // failure.
        const size_t num_results = 2;

        // To catch values returned by the search
        std::vector<size_t> ret_index(num_results);
        std::vector<Real> out_dist_sqr(num_results);
        nanoflann::KNNResultSet<Real> result_set(num_results);

        // Initialize the result_set
        result_set.init(&ret_index[0], &out_dist_sqr[0]);

        // Do the search
        kd_tree.findNeighbors(result_set, &query_pt[0], nanoflann::SearchParams(10));

        // There should be one zero-distance node found, but not more.
        // Don't print out both halves of the duplicate pair, only
        // when ret_index[r] < node->id().
        for (unsigned r=0; r<result_set.size(); ++r)
          if (ret_index[r] < node->id() && out_dist_sqr[r] < TOLERANCE)
            {
              // libMesh::out << "Node " << ret_index[r] << " is a duplicate of node " << node->id() << std::endl;
              dupe_count++;
            }
      }

    perf_main.pop("Nanoflann");

    // Only print the warning if there were duplicates found, and if so, make the output message stand out.
    if (dupe_count)
      {
        libMesh::out << std::endl;
        libMesh::out << "********************************************************************************" << std::endl;
        libMesh::out << "Warning! Nanoflann found " << dupe_count << " duplicate nodes." << std::endl;
        libMesh::out << "********************************************************************************" << std::endl;
        libMesh::out << std::endl;
      }
  }

  // Log how long the naive N^2 algorithm takes.  Since it can take
  // a long time, also give the option to skip it.
  bool test_n2 = false;
  if (command_line.search(1, "--test-n2"))
    test_n2 = true;

  if (test_n2)
    {
      libMesh::out << "Starting naive N^2 algorithm." << std::endl;

      perf_main.push("N^2 algorithm");

      // Count the number of duplicate nodes found
      unsigned dupe_count = 0;

      // For every node in the mesh, search the KDtree and find any
      // zero-distance neighbors with a different index from the current
      // node being searched... this indicates a duplicate node.
      for (const auto & node : mesh.node_ptr_range())
        for (const auto & inner_node : mesh.node_ptr_range())
          {
            // Don't even bother to make the comparison unless node->id() < inner_node->id().
            // The duplicate node relation is symmetric, and every node is a duplicate to itself.
            if (node->id() < inner_node->id() && ((*node) - (*inner_node)).norm_sq() < TOLERANCE)
              {
                // libMesh::out << "Node " << node->id() << " is a duplicate of node " << inner_node->id() << std::endl;
                dupe_count++;

                // We only need to find at most 1 duplicate, so we can break.
                break;
              }
          }

      perf_main.pop("N^2 algorithm");

      if (dupe_count)
        {
          libMesh::out << std::endl;
          libMesh::out << "********************************************************************************" << std::endl;
          libMesh::out << "Naive N^2 algorihtm found " << dupe_count << " duplicate nodes." << std::endl;
          libMesh::out << "********************************************************************************" << std::endl;
          libMesh::out << std::endl;
        }

      // Compute the timing ratio
      libMesh::out << "Ratio of N^2 algorithm to Nanoflann = "
                   << perf_main.get_perf_data("N^2 algorithm").tot_time / perf_main.get_perf_data("Nanoflann").tot_time
                   << "X"
                   << std::endl;
    }


  return 0;
}