Simple example code demonstrating calling netgen from libmesh

netgen_test_01.cc

#include <iostream>
#include <fstream>
#include <map>
#include <cstdlib> // std::getenv()

#include "libmesh/libmesh.h"
#include "libmesh/mesh.h"
#include "libmesh/point.h"
#include "libmesh/cell_tet4.h"
// #include "libmesh/exodusII_io.h"

// NOTE: This example only works if you have built and installed
// netgen, (see http://sourceforge.net/projects/netgen-mesher) and you
// are using a Makefile that sets the appropriate include and library
// paths.
namespace nglib
{
// The cmake build of netgen does not install nglib.h???  Is this a new thing?
// There do not appear to be any other working headers installed...
#include <nglib.h>

// The nginterface.h header is invalid, it uses a naked string without prepending std::
// #include <nginterface.h>
}

// This header already wraps all of its contents in namespace netgen,
// but it can't be included directly because it uses stuff like
// NG_ELEMENT_TYPE which has not been defined.
// #include <nginterface_v2.hpp>

using namespace libMesh;

// This code reads cube.surf, which consists of a set of points and
// triangles which define the surface of a cube.  This code is copied
// from the ng_vol.cpp file which comes with the netgen distribution.
int main (int argc, char** argv)
{
  // Initialize the library.  This is necessary because the library
  // may depend on a number of other libraries (i.e. MPI  and Petsc)
  // that require initialization before use.
  LibMeshInit init(argc, argv);
  {
    // Must be called before any other nglib functions.
    nglib::Ng_Init();
    nglib::Ng_Mesh * netgen_mesh = nglib::Ng_NewMesh();

    // Read points and triangles from input surface file and store
    // them in mesh.
    {
      std::ifstream in("cube.surf");

      int np;
      in >> np;

      double point[3];
      std::cout << "Reading " << np  << " points...";
      for (int i = 1; i <= np; i++)
        {
          in >> point[0] >> point[1] >> point[2];
          nglib::Ng_AddPoint (netgen_mesh, point);
        }
      std::cout << "done" << std::endl;

      int nfaces;
      in >> nfaces;

      int trig[3];
      std::cout << "Reading " << nfaces  << " faces...";
      for (int i = 1; i <= nfaces; i++)
        {
          in >> trig[0] >> trig[1] >> trig[2];
          nglib::Ng_AddSurfaceElement (netgen_mesh, nglib::NG_TRIG, trig);
        }
      std::cout << "done" << std::endl;
    }

    // generate volume mesh
    nglib::Ng_Meshing_Parameters mp;
    // maxh is the "maximum global mesh size allowed".  It sounds like
    // this would allow one to generate a finer grid, but there is a
    // thread (with no responses) about how this does not work...
    // https://sourceforge.net/p/netgen-mesher/discussion/905307/thread/cd48bbc3/
    // Changing minh has no effect either.  Perhaps these parameters
    // are not allowed to insert additional points in the Mesh?
    mp.maxh = .5;
    mp.minh = .5;
    mp.fineness = 1;
    mp.second_order = 0;

    // We can also try to call an API that restricts the size of the
    // mesh globally...  This also has no effect on the volume mesh,
    // no matter how small I pick h to be.
    // nglib::Ng_RestrictMeshSizeGlobal(netgen_mesh, /*h=*/.001);

    // There's also an interface to restrict the mesh size near a point.  This seemed
    // to work, but only if I picked a number less than .5, like .1.  I also could not
    // get this to trigger refinement near the boundary, so maybe we are not allowed
    // to change the surface mesh?
    // double restriction_point[3] = {.5, .5, .5}; // triggers refinement
    // // double restriction_point[3] = {1, 1, 1}; // no refinement triggered
    // // double restriction_point[3] = {.5, .5, 1}; // no refinement triggered
    // // double restriction_point[3] = {.49, .49, .99}; // no refinement triggered
    // nglib::Ng_RestrictMeshSizePoint(netgen_mesh, restriction_point, /*h=*/.01);

    // It is also possible to restrict the maximum mesh size in a box
    // with the following interface.  This still does not have any
    // effect on the boundary mesh, however.
    // Notes:
    // .) bottom=.1, top=.9, h=.1 results in 1067 elements
    // .) bottom=.1, top=.9, h=.05 results in 9818 elements
    // .) bottom=.1, top=.9, h=.025 results in 122551 elements
    // double
    //   bottom = .1, top = .9;
    // double
    //   pmin[3] = {bottom, bottom, bottom},
    //   pmax[3] = {top, top, top};
    // nglib::Ng_RestrictMeshSizeBox(netgen_mesh, pmin, pmax, /*h=*/.025);

    std::cout << "start meshing" << std::endl;
    nglib::Ng_GenerateVolumeMesh (netgen_mesh, &mp);
    std::cout << "meshing done" << std::endl;

    // Build libmesh mesh to transfer points into
    Mesh mesh(init.comm(), 3);

    // Print lists of points and elements in volume mesh.
    {
      int np = nglib::Ng_GetNP(netgen_mesh);
      std::cout << "Points: " << np << std::endl;

      double point[3];
      for (int i = 1; i <= np; i++)
        {
          nglib::Ng_GetPoint (netgen_mesh, i, point);

          // Debugging: print the points as they are extracted from the netgen mesh
          // std::cout << i
          //           << ": " << point[0]
          //           << " " << point[1]
          //           << " " << point[2]
          //           << std::endl;

          // Add points to libmesh mesh, use 1-based numbering like netgen.
          mesh.add_point(Point(point[0], point[1], point[2]), i);
        }

      int ne = nglib::Ng_GetNE(netgen_mesh);
      std::cout << "Elements: " << ne << std::endl;

      int tet[4];
      for (int i = 1; i <= ne; i++)
        {
          nglib::Ng_GetVolumeElement (netgen_mesh, i, tet);

          // Debugging: print connectivity
          // std::cout << i
          //           << ": " << tet[0]
          //           << " " << tet[1]
          //           << " " << tet[2]
          //           << " " << tet[3]
          //           << std::endl;

          // Add tet to libmesh mesh, use 1-based numbering like
          // netgen.  Note that we have to swap the last two entries
          // of the connectivity array, because Netgen uses a
          // clockwise numbering for some reason.
          Elem* elem = mesh.add_elem(new Tet4);
          elem->set_node(0) = mesh.node_ptr(tet[0]);
          elem->set_node(1) = mesh.node_ptr(tet[1]);
          elem->set_node(2) = mesh.node_ptr(tet[3]);
          elem->set_node(3) = mesh.node_ptr(tet[2]);
        }
    }

    // Write mesh to Netgen-style file.
    nglib::Ng_SaveMesh(netgen_mesh, "cube.vol");

    // We're done with the Netgen mesh
    nglib::Ng_DeleteMesh(netgen_mesh);

    // We're now done calling Netgen functions
    nglib::Ng_Exit();

    // This is necessary, otherwise I get very strange results for
    // some reason (points with different coordinates than they are
    // supposed to have!
    mesh.prepare_for_use();

    // Write an Exodus file with the libmesh mesh
    mesh.write("netgen_mesh.e");
  }

  return 0;
}