Test demonstrating meshing a 2D surface with Netgen from within libmesh

netgen_test_02.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/face_tri3.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
{
#include <nglib.h>
}

using namespace libMesh;

// This code adds 2D points and line segments to a Netgen Mesh object.
int main (int argc, char** argv)
{
  // Initialize libmesh.  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();

    // Let netgen allocate the mesh below.
    nglib::Ng_Mesh * netgen_mesh = NULL;

    // There does not appear to be a function that generates a 2D mesh given
    // a set of boundary points and segments. There are:
    //
    // Ng_LoadGeometry_2D
    // Ng_GenerateMesh_2D
    //
    // so one approach is to create a simple file in the right format
    // and then call Ng_LoadGeometry_2D(), as this spits back a Ng_Geometry_2D
    // pointer.  The following input file is from:
    // http://www.asc.tuwien.ac.at/~lehrenfeld/xfem/exercise/ex1_poisson/poisson.pdf
    //
    // splinecurves2dv2
    // 5
    // points
    // 1 -1 -1
    // 2 1 -1
    // 3 1 1
    // 4 -1 1
    // segments
    // # domain-left, domain-right, num-points, point1, point2
    // # domain 0 == outside
    // # domain 1 == inside
    // 1 0 2 1 2
    // 1 0 2 2 3
    // 1 0 2 3 4
    // 1 0 2 4 1
    // materials
    // 1 inner -maxh=0.2
    std::string filename = "netgen_io.tmp";
    {
      std::ofstream netgen_io(filename.c_str());
      netgen_io << "splinecurves2dv2\n";
      netgen_io << "5\n";
      netgen_io << "points\n";
      netgen_io << "1 -1 -1\n";
      netgen_io << "2 1 -1\n";
      netgen_io << "3 1 1\n";
      netgen_io << "4 -1 1\n";
      netgen_io << "segments\n";
      netgen_io << "1 0 2 1 2\n";
      netgen_io << "1 0 2 2 3\n";
      netgen_io << "1 0 2 3 4\n";
      netgen_io << "1 0 2 4 1\n";
      netgen_io << "materials\n";
      netgen_io << "1 inner -maxh=0.2\n";
    }

    nglib::Ng_Geometry_2D * geom = nglib::Ng_LoadGeometry_2D (filename.c_str());

    // Remove the temporary file now that we are done with it.
    std::remove(filename.c_str());

    // Call netgen to generate the 2D mesh.  According to the docs,
    // the function is supposed to allocate the Netgen mesh object,
    // but there is a bug in upstream master -- the first thing this
    // function tries to do is call mesh->SetDimension(2) on a NULL
    // pointer...
    //
    // The Mesh class is defined in libsrc/meshing/meshclass.hpp
    // MeshFromSpline2D is in libsrc/geom2d//genmesh2d.cpp
    // This is now fixed by my commit: 3867621 in
    // https://github.com/jwpeterson/netgen
    nglib::Ng_Meshing_Parameters mp;
    nglib::Ng_Result res = nglib::Ng_GenerateMesh_2D(geom, &netgen_mesh, &mp);

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

   // Print lists of points and elements in volume mesh.  Note: we are using
   // the 2D version of the C interfaces in what follows below.
   {
     int np = nglib::Ng_GetNP_2D(netgen_mesh);
     std::cout << "Points: " << np << std::endl;

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

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

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

     int tri[3];
     for (int i = 1; i <= ne; i++)
       {
         nglib::Ng_GetElement_2D (netgen_mesh, i, tri);

         // Add tri to libmesh mesh, use 1-based numbering like
         // netgen.  I think the triangles are numbered the same
         // way as in libmesh.
         Elem* elem = mesh.add_elem(new Tri3);
         elem->set_node(0) = mesh.node_ptr(tri[0]);
         elem->set_node(1) = mesh.node_ptr(tri[1]);
         elem->set_node(2) = mesh.node_ptr(tri[2]);
       }
   }

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

    // 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_test_02.e");
  }

  return 0;
}