contains_point.C

Go to the documentation of this file.

#include <libmesh/libmesh.h>
#include <libmesh/elem.h>

#include "test_comm.h"
#include "libmesh_cppunit.h"


using namespace libMesh;

class ContainsPointTest : public CppUnit::TestCase
{
public:
  LIBMESH_CPPUNIT_TEST_SUITE( ContainsPointTest );

#if LIBMESH_DIM > 2
  CPPUNIT_TEST( testContainsPointNodeElem );
  CPPUNIT_TEST( testContainsPointTri3 );
  CPPUNIT_TEST( testContainsPointTet4 );
#endif

  CPPUNIT_TEST_SUITE_END();

public:
  void setUp() {}

  void tearDown() {}

  // NodeElem test
  void testContainsPointNodeElem()
  {
    Node node  (1., 1., 1., /*id=*/0);
    std::unique_ptr<Elem> elem = Elem::build(NODEELEM);
    elem->set_node(0, &node);

    Real epsilon = 1.e-4;
    CPPUNIT_ASSERT(elem->contains_point(Point(1.+epsilon/2, 1.-epsilon/2, 1+epsilon/2), /*tol=*/epsilon));
    CPPUNIT_ASSERT(!elem->contains_point(Point(1.+epsilon/2, 1.-epsilon/2, 1+epsilon/2), /*tol=*/epsilon/2));
  }

  // TRI3 test
  void testContainsPointTri3()
  {
    LOG_UNIT_TEST;

    Point a(3,1,2), b(1,2,3), c(2,3,1);
    containsPointTri3Helper(a, b, c, a);

    // Ben's 1st "small triangle" test case.
    containsPointTri3Helper(a/5000, b/5000, c/5000, c/5000);

    // Ben's 2nd "small triangle" test case.
    containsPointTri3Helper(Point(0.000808805, 0.0047284,  0.),
                            Point(0.0011453,   0.00472831, 0.),
                            Point(0.000982249, 0.00508037, 0.),
                            Point(0.001,       0.005,      0.));
  }



  // TET4 test
  void testContainsPointTet4()
  {
    LOG_UNIT_TEST;

    // Construct unit Tet.
    {
      Node zero  (0., 0., 0., 0);
      Node one   (1., 0., 0., 1);
      Node two   (0., 1., 0., 2);
      Node three (0., 0., 1., 3);

      std::unique_ptr<Elem> elem = Elem::build(TET4);
      elem->set_node(0, &zero);
      elem->set_node(1, &one);
      elem->set_node(2, &two);
      elem->set_node(3, &three);

      // The centroid (equal to vertex average for Tet4) must be inside the element.
      CPPUNIT_ASSERT (elem->contains_point(elem->vertex_average()));

      // The vertices must be contained in the element.
      CPPUNIT_ASSERT (elem->contains_point(zero));
      CPPUNIT_ASSERT (elem->contains_point(one));
      CPPUNIT_ASSERT (elem->contains_point(two));
      CPPUNIT_ASSERT (elem->contains_point(three));

      // Make sure that outside points are not contained.
      CPPUNIT_ASSERT (!elem->contains_point(Point(.34, .34, .34)));
      CPPUNIT_ASSERT (!elem->contains_point(Point(.33, .33, -.1)));
      CPPUNIT_ASSERT (!elem->contains_point(Point(0., -.1, .5)));
    }


    // Construct a nearly degenerate (sliver) tet.  A unit tet with
    // nodes "one" and "two" moved to within a distance of epsilon
    // from the origin.
    {
      Real epsilon = 1.e-4;

      Node zero  (0., 0., 0., 0);
      Node one   (epsilon, 0., 0., 1);
      Node two   (0., epsilon, 0., 2);
      Node three (0., 0., 1., 3);

      std::unique_ptr<Elem> elem = Elem::build(TET4);
      elem->set_node(0, &zero);
      elem->set_node(1, &one);
      elem->set_node(2, &two);
      elem->set_node(3, &three);

      // The centroid (equal to vertex average for Tet4) must be inside the element.
      CPPUNIT_ASSERT (elem->contains_point(elem->vertex_average()));

      // The vertices must be contained in the element.
      CPPUNIT_ASSERT (elem->contains_point(zero));
      CPPUNIT_ASSERT (elem->contains_point(one));
      CPPUNIT_ASSERT (elem->contains_point(two));
      CPPUNIT_ASSERT (elem->contains_point(three));

      // Verify that a point which is on a mid-edge is contained in the element.
      CPPUNIT_ASSERT (elem->contains_point(Point(epsilon/2, 0, 0.5)));

      // Make sure that "just barely" outside points are outside.
      CPPUNIT_ASSERT (!elem->contains_point(Point(epsilon, epsilon, epsilon/2)));
      CPPUNIT_ASSERT (!elem->contains_point(Point(epsilon/10, epsilon/10, 1.0)));
      CPPUNIT_ASSERT (!elem->contains_point(Point(epsilon/2, -epsilon/10, 0.5)));
    }
  }

protected:
  // The first 3 arguments are the points of the triangle, the last argument is a
  // custom point that should be inside the Tri.
  void containsPointTri3Helper(Point a_in, Point b_in, Point c_in, Point p)
  {
    // vertices
    Node a(a_in, 0);
    Node b(b_in, 1);
    Node c(c_in, 2);

    // Build the test Elem
    std::unique_ptr<Elem> elem = Elem::build(TRI3);

    elem->set_node(0, &a);
    elem->set_node(1, &b);
    elem->set_node(2, &c);

    // helper vectors to span the tri and point out of plane
    Point va(a-c);
    Point vb(b-c);
    Point oop(va.cross(vb));
    Point oop_unit = oop.unit();

    // The centroid (equal to vertex average for Tri3) must be inside the element.
    CPPUNIT_ASSERT (elem->contains_point(elem->vertex_average()));

    // triangle should contain all its vertices
    CPPUNIT_ASSERT (elem->contains_point(a));
    CPPUNIT_ASSERT (elem->contains_point(b));
    CPPUNIT_ASSERT (elem->contains_point(c));

    // out of plane from the centroid, should *not* be in the element
    CPPUNIT_ASSERT (!elem->contains_point(elem->vertex_average() + std::sqrt(TOLERANCE) * oop_unit));

    // These points should be outside the triangle
    CPPUNIT_ASSERT (!elem->contains_point(a + va * TOLERANCE * 10));
    CPPUNIT_ASSERT (!elem->contains_point(b + vb * TOLERANCE * 10));
    CPPUNIT_ASSERT (!elem->contains_point(c - (va + vb) * TOLERANCE * 10));

    // Test the custom point
    CPPUNIT_ASSERT (elem->contains_point(p));
  }
};


CPPUNIT_TEST_SUITE_REGISTRATION( ContainsPointTest );