tet_contains_point_test.cc

#include "libmesh/libmesh.h"
#include "libmesh/mesh.h"
#include "libmesh/elem.h"
#include "libmesh/mesh_generation.h"
#include "libmesh/mesh_modification.h"

using namespace libMesh;

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

    // Create the Mesh
    Mesh mesh(init.comm());

    // Note that there is an extra factor of 24 in each direction for Tets.
    const unsigned n_elem = 10;

    MeshTools::Generation::build_cube (mesh,
                                       n_elem, n_elem, n_elem,
                                       0., 1,
                                       0., 1,
                                       0., 1,
                                       TET4);

    // Introduce some irregularity in the Mesh.
    MeshTools::Modification::distort(mesh, 0.3, /*perturb_boundary=*/false);

    // Build a PointLocator object.  This is the most common way that contains_point() calls are triggered.
    UniquePtr<PointLocatorBase> pl_ptr = PointLocatorBase::build(TREE, mesh);
    PointLocatorBase & pl = *pl_ptr;

    // How many random points do we want to find?  Note: each point
    // will trigger multiple contains_point() calls.
    const unsigned n_points = 1e6;
    unsigned points_found = 0;

    PerfLog perf_log ("Locate points");
    perf_log.push("locate points");

    for (unsigned i=0; i<n_points; ++i)
      {
        Point p(random()/static_cast<Real>(RAND_MAX),
                random()/static_cast<Real>(RAND_MAX),
                random()/static_cast<Real>(RAND_MAX));

        const Elem * elem = pl(p);

        if (!elem)
          libmesh_error_msg("Point not found!");

        // std::cout << "Point found in elem " << elem->id() << std::endl;
        points_found++;
      }

    perf_log.pop("locate points");

    libMesh::out << "points_found = " << points_found << "." << std::endl;

    // Inspect the mesh with a clip plane to see how much distortion there is.
    // mesh.write("tet_contains_point_test.e");

  } // LibMeshInit


  return 0;
}


// Results 1e6 nodes:
// Original contains_point() routine from Elem base class (sorted):
// 46.5478,46.6196,46.9555,47.1207,47.4358; median=46.9555
//
// Optimized contains_point() routine:
// 15.5872,15.7150,15.7725,15.7873,15.7888; median=15.7725