Add unit test for Netgen interface

We need a lot more here but this is a start

tests/Makefile.am

@@ -74,6 +74,7 @@ unit_tests_sources = \
   mesh/mesh_generation_test.C \
   mesh/mesh_input.C \
   mesh/mesh_stitch.C \
+  mesh/mesh_tet_test.C \
   mesh/mesh_triangulation.C \
   mesh/mixed_dim_mesh_test.C \
   mesh/mixed_order_test.C \

tests/mesh/mesh_tet_test.C

@@ -0,0 +1,305 @@
+#include <libmesh/boundary_info.h>
+#include <libmesh/elem.h>
+#include <libmesh/mesh.h>
+#include <libmesh/mesh_generation.h>
+#include <libmesh/mesh_netgen_interface.h>
+#include <libmesh/mesh_tetgen_interface.h>
+#include <libmesh/parallel_implementation.h> // max()
+
+#include "test_comm.h"
+#include "libmesh_cppunit.h"
+
+#include <algorithm>
+#include <cmath>
+#include <regex>
+
+
+using namespace libMesh;
+
+class MeshTetTest : public CppUnit::TestCase
+{
+  /**
+   * The goal of this test is to verify proper operation of the
+   * interfaces to tetrahedralization libraries
+   */
+public:
+  LIBMESH_CPPUNIT_TEST_SUITE( MeshTetTest );
+
+#ifdef LIBMESH_HAVE_NETGEN
+  // The most basic test to start with
+  CPPUNIT_TEST( testNetGen );
+
+  // We'll get to more advanced features later
+  /*
+  CPPUNIT_TEST( testNetGenInterp );
+  CPPUNIT_TEST( testNetGenInterp2 );
+  CPPUNIT_TEST( testNetGenRefined );
+  CPPUNIT_TEST( testNetGenNonRefined );
+  CPPUNIT_TEST( testNetGenExtraRefined );
+  */
+#endif
+
+  // Still need to work out the basics here - non-convex domains,
+  // precise control of refinement, etc.
+#ifdef LIBMESH_HAVE_TETGEN
+  /*
+  CPPUNIT_TEST( testTetGen );
+  CPPUNIT_TEST( testTetGenInterp );
+  CPPUNIT_TEST( testTetGenInterp2 );
+  */
+#endif
+
+  CPPUNIT_TEST_SUITE_END();
+
+public:
+  void setUp() {}
+
+  void tearDown() {}
+
+  void commonSettings (MeshTetInterface & tetinterface)
+  {
+    // Don't try to insert points unless we're requested to later
+    tetinterface.desired_volume() = 0;
+    tetinterface.smooth_after_generating() = false;
+  }
+
+  void testExceptionBase(MeshBase & mesh,
+                         MeshTetInterface & tetinterface,
+                         const char * re)
+  {
+#ifdef LIBMESH_ENABLE_EXCEPTIONS
+    // We can't just CPPUNIT_ASSERT_THROW, because we want to make
+    // sure we were thrown from the right place with the right error
+    // message!
+    bool threw_desired_exception = false;
+    try {
+      this->testTetInterfaceBase(mesh, tetinterface);
+    }
+    catch (libMesh::LogicError & e) {
+      std::regex msg_regex(re);
+      CPPUNIT_ASSERT(std::regex_search(e.what(), msg_regex));
+      threw_desired_exception = true;
+    }
+    catch (CppUnit::Exception & e) {
+      throw e;
+    }
+    catch (...) {
+      CPPUNIT_ASSERT_MESSAGE("Unexpected exception type thrown", false);
+    }
+    CPPUNIT_ASSERT(threw_desired_exception);
+#endif
+  }
+
+
+  void testTetInterfaceBase(MeshBase & mesh,
+                            MeshTetInterface & triangulator)
+  {
+    commonSettings(triangulator);
+
+    triangulator.triangulate();
+
+    CPPUNIT_ASSERT_EQUAL(mesh.n_elem(), static_cast<dof_id_type>(4));
+    CPPUNIT_ASSERT_EQUAL(mesh.n_nodes(), static_cast<dof_id_type>(6));
+    for (const auto & elem : mesh.element_ptr_range())
+      {
+        CPPUNIT_ASSERT_EQUAL(elem->type(), TET4);
+
+        // Make sure we're not getting any inverted elements
+        CPPUNIT_ASSERT(!elem->is_flipped());
+      }
+  }
+
+
+  void testTetInterface(MeshBase & mesh,
+                        MeshTetInterface & triangulator)
+  {
+    // An asymmetric octahedron, so we hopefully have an unambiguous
+    // choice of shortest diagonal for a Delaunay algorithm to pick.
+    mesh.add_point(Point(0,0,-0.1), 0);
+    mesh.add_point(Point(1,0,0), 1);
+    mesh.add_point(Point(0,1,0), 2);
+    mesh.add_point(Point(-1,0,0), 3);
+    mesh.add_point(Point(0,-1,0), 4);
+    mesh.add_point(Point(0,0,0.1), 5);
+
+    auto add_tri = [&mesh](std::array<dof_id_type,3> nodes)
+    {
+      auto elem = mesh.add_elem(Elem::build(TRI3));
+      elem->set_node(0) = mesh.node_ptr(nodes[0]);
+      elem->set_node(1) = mesh.node_ptr(nodes[1]);
+      elem->set_node(2) = mesh.node_ptr(nodes[2]);
+    };
+
+    add_tri({0,1,2});
+    add_tri({0,2,3});
+    add_tri({0,3,4});
+    add_tri({0,4,1});
+    add_tri({5,1,4});
+    add_tri({5,4,3});
+    add_tri({5,3,2});
+    add_tri({5,2,1});
+
+    mesh.prepare_for_use();
+
+    this->testTetInterfaceBase(mesh, triangulator);
+  }
+
+
+#ifdef LIBMESH_HAVE_TRIANGLE
+  void testTetGen()
+  {
+    LOG_UNIT_TEST;
+
+    Mesh mesh(*TestCommWorld);
+    TetGenMeshInterface tet_tet(mesh);
+    testTetInterface(mesh, tet_tet);
+  }
+
+
+  /*
+  void testTetGenInterp()
+  {
+    LOG_UNIT_TEST;
+
+    Mesh mesh(*TestCommWorld);
+    TetGenMeshInterface tet_tet(mesh);
+    testTetInterfaceInterp(mesh, tet_tet, 1, 6);
+  }
+
+
+  void testTetGenInterp2()
+  {
+    LOG_UNIT_TEST;
+
+    Mesh mesh(*TestCommWorld);
+    TetGenMeshInterface tet_tet(mesh);
+    testTetInterfaceInterp(mesh, tet_tet, 2, 10);
+  }
+  */
+
+#endif // LIBMESH_HAVE_TETGEN
+
+
+#ifdef LIBMESH_HAVE_NETGEN
+  void testNetGen()
+  {
+    LOG_UNIT_TEST;
+
+    Mesh mesh(*TestCommWorld);
+    NetGenMeshInterface net_tet(mesh);
+    testTetInterface(mesh, net_tet);
+  }
+
+
+  /*
+  void testNetGenInterp()
+  {
+    LOG_UNIT_TEST;
+
+    Mesh mesh(*TestCommWorld);
+    NetGenMeshInterface net_tet(mesh);
+    testTetInterfaceInterp(mesh, net_tet, 1, 6);
+  }
+
+
+  void testNetGenInterp2()
+  {
+    LOG_UNIT_TEST;
+
+    Mesh mesh(*TestCommWorld);
+    NetGenMeshInterface net_tet(mesh);
+    testTetInterfaceInterp(mesh, net_tet, 2, 10);
+  }
+
+
+  void testNetGenRefinementBase
+    (UnstructuredMesh & mesh,
+     const std::vector<MeshTetInterface::Hole*> * holes,
+     Real expected_total_area,
+     dof_id_type n_original_elem,
+     Real desired_area = 0.1,
+     FunctionBase<Real> * area_func = nullptr)
+  {
+    NetGenMeshInterface triangulator(mesh);
+
+    commonSettings(triangulator);
+
+    if (holes)
+      triangulator.attach_hole_list(holes);
+
+    // Try to insert points!
+    triangulator.desired_area() = desired_area;
+    triangulator.set_desired_area_function(area_func);
+
+    triangulator.triangulate();
+
+    // If refinement should have increased our element count, check it
+    if (desired_area || area_func)
+      CPPUNIT_ASSERT_GREATER(n_original_elem, mesh.n_elem()); // n_elem+++
+    else
+      CPPUNIT_ASSERT_EQUAL(mesh.n_elem(), n_original_elem);
+
+    Real area = 0;
+    for (const auto & elem : mesh.active_local_element_ptr_range())
+      {
+        CPPUNIT_ASSERT_EQUAL(elem->level(), 0u);
+        CPPUNIT_ASSERT_EQUAL(elem->type(), TRI3);
+
+        const Real my_area = elem->volume();
+
+        // my_area <= desired_area, wow this macro ordering hurts
+        if (desired_area != 0)
+          CPPUNIT_ASSERT_LESSEQUAL(desired_area, my_area);
+
+        if (area_func != nullptr)
+          for (auto v : make_range(elem->n_vertices()))
+            {
+              const Real local_desired_area =
+                (*area_func)(elem->point(v));
+              CPPUNIT_ASSERT_LESSEQUAL(local_desired_area, my_area);
+            }
+
+        area += my_area;
+      }
+
+    mesh.comm().sum(area);
+
+    LIBMESH_ASSERT_FP_EQUAL(area, expected_total_area, TOLERANCE*TOLERANCE);
+  }
+
+  void testNetGenRefined()
+  {
+    LOG_UNIT_TEST;
+
+    Mesh mesh(*TestCommWorld);
+    testTriangulatorTrapMesh(mesh);
+    testPoly2TriRefinementBase(mesh, nullptr, 1.5, 15);
+  }
+
+  void testNetGenNonRefined()
+  {
+    LOG_UNIT_TEST;
+
+    Mesh mesh(*TestCommWorld);
+    testTriangulatorTrapMesh(mesh);
+    // Make sure we see 0 as "don't refine", not "infinitely refine"
+    testPoly2TriRefinementBase(mesh, nullptr, 1.5, 2, 0);
+  }
+
+
+  void testNetGenExtraRefined()
+  {
+    LOG_UNIT_TEST;
+
+    Mesh mesh(*TestCommWorld);
+    testTriangulatorTrapMesh(mesh);
+    testPoly2TriRefinementBase(mesh, nullptr, 1.5, 150, 0.01);
+  }
+*/
+
+#endif // LIBMESH_HAVE_NETGEN
+
+};
+
+
+CPPUNIT_TEST_SUITE_REGISTRATION( MeshTetTest );