Mesh generation misses surface components depending on node ordering

[roystgnr]

I can take even a simple domain (one cube inside another: 16 points connected by 24 triangles) and get Netgen to either give me a correct tetrahedral mesh or an incorrect one (omitting the hole defined by the internal cube) depending on the numbering (implicit in the ordering of the calls to Ng_AddPoint, then kept consistent for the data handed to Ng_AddSurfaceElement) of the boundary vertices.

This was triggered when wrapped by another library wrapped by another code ( https://github.com/idaholab/moose/pull/28298 ), but I've distilled it into an independent (g++ test_netgen.C -o test_netgen.x -I$MY_NETGEN_SRCDIR -L $MY_NETGEN_BUILDDIR -Wl,-rpath -Wl,$MY_NETGEN_BUILDDIR -lnglib on Linux) executable with the following test_netgen.C:

#include <array>
#include <cassert>
#include <cstddef>
#include <iostream>
#include <numeric>
#include <set>
#include <vector>

namespace nglib {
#include "netgen/nglib/nglib.h"
}

using namespace nglib;

  std::size_t testPermutedHole (const std::vector<std::size_t> & permutation)
  { 
    auto ngmesh = Ng_NewMesh();

    Ng_Meshing_Parameters params;
    params.elementsperedge = 1;
    params.elementspercurve = 1;

    std::vector<std::array<double, 3>> point_coords =
      {{-4, 4, -4},
      {-4, 4, 4},
      {4, 4, 4},
      {4, -4, 4},
      {4, 4, -4},
      {4, -4, -4},
      {-4, -4, 4},
      {-4, -4, -4},

      {3, -3, -3},
      {-3, -3, 3},
      {-3, -3, -3},
      {-3, 3, -3},
      {-3, 3, 3},
      {3, 3, 3},
      {3, -3, 3},
      {3, 3, -3}};

    auto n_input_points = point_coords.size();

    assert (permutation.size() == n_input_points);

    std::vector<std::size_t> inverse_permutation (n_input_points);

    for (std::size_t i=0; i != n_input_points; ++i)
      {
        inverse_permutation[permutation[i]] = i;
        Ng_AddPoint(ngmesh, point_coords[permutation[i]].data());
      }

    std::vector<std::array<int, 3>> tri_nodes =
      {{1, 2, 3},
      {4, 3, 2},
      {5, 1, 3},
      {5, 3, 4},
      {1, 5, 6},
      {6, 5, 4},
      {6, 4, 7},
      {7, 4, 2},
      {8, 1, 6},
      {8, 6, 7},
      {1, 8, 2},
      {8, 7, 2},

      {11, 10, 9},
      {11, 9, 12},
      {10, 11, 13},
      {13, 12, 14},
      {12, 13, 11},
      {9, 10, 15},
      {13, 14, 15},
      {13, 15, 10},
      {15, 14, 16},
      {16, 14, 12},
      {15, 16, 9},
      {9, 16, 12}};

    for (auto & nodes : tri_nodes)
      {
        std::array<int, 3> permuted_nodes = nodes;
        for (auto & node_i : permuted_nodes)  
          node_i = inverse_permutation[node_i-1]+1;
        Ng_AddSurfaceElement(ngmesh, NG_TRIG, permuted_nodes.data());
      }

    Ng_GenerateVolumeMesh(ngmesh, &params);   
    const std::size_t n_elem = Ng_GetNE(ngmesh);
    const std::size_t n_points = Ng_GetNP(ngmesh);
    assert(n_points >= n_input_points);

    std::set<int> nodes_seen;
    for (std::size_t i = 0; i != n_elem; ++i) 
      {
        // Avoid segfault even if ngtype isn't a tet4
        int ngnodes[11];
        Ng_Volume_Element_Type ngtype =
          Ng_GetVolumeElement(ngmesh, i+1, ngnodes);
        assert(ngtype == NG_TET);
        for (int n = 0; n != 4; ++n)
          if (ngnodes[n] <= int(n_input_points) && ngnodes[n] > 0)
            nodes_seen.insert(ngnodes[n]);
      }

    Ng_DeleteMesh(ngmesh);

    return nodes_seen.size();
  }

  int main(void)
  {
    const std::size_t n_input_points = 16;

    std::vector<std::size_t> permutation (n_input_points);
    std::iota(permutation.begin(), permutation.end(), 0);

    int fails=0, successes = 0;

    auto run_test = [&fails, &successes, &permutation]()
    {
      int n_original_nodes = testPermutedHole(permutation);
      fails += (n_original_nodes != n_input_points);
      successes += (n_original_nodes == n_input_points);
      std::cout << "Found " << n_original_nodes << "/" << n_input_points << " original nodes\n";
      std::cout << "Fails = " << fails << ", successes = " << successes << std::endl;
    };

    // Heap's Algorithm, non-recursive version
    run_test();

    std::vector<std::size_t> c(permutation.size());

    for (std::size_t i = 1; i != n_input_points ; ++i)
      {
        if (c[i] < i)
          {
            if (i%2)
              std::swap(permutation[c[i]], permutation[i]);
            else
              std::swap(permutation[0], permutation[i]);
            run_test();

            ++c[i];
            i = 0; // Yeah, this runs in factorial time
          }
        else
          c[i] = 0;
      }

    return 0;
  }

With the nodes in their original numbering in the code (hand-copied from my real code's input in the simplest failure I could easily reproduce), Netgen meshes the outer cube with 12 tets, none of which connect to any nodes from the inner cube. With most permutations (75%? But Heap's algorithm isn't random sampling and I'm not going to wait for 16! cases to finish) I instead get more like what I expect: many more tets (36 whenever I check), with all 16 nodes of the inner and outer cube accounted for among them.

Are there any requirements for the ordering of input nodes (and/or surface elements? my original failure case pseudorandomizes those too, and I was going to try permutations there if permuting nodes hadn't been sufficient) for Netgen to generate a tetrahedralization respecting each surface?

Or is there an Ng_Meshing_Parameters setting I'm missing that might cause some surface components to be removed? I've tried turning off optsurfmeshenable and optvolmeshenable to no effect, and nothing else looked like it might be responsible for any form of deliberate decimation.

[roystgnr]

Seriously? A malware bot comment within seconds of issue creation time? My apologies for however I attracted its attention; I hope this doesn't spread to any repos I work in...

[ChrLackner]

Interesting catch! Thanks for the easy testcase.

Should be fixed with https://github.com/NGSolve/netgen/commit/4964691fd97f735300030ed7d8e4764d3b478f5a

Note: it does not necessarily manage to fill without additional nodes. this is a bit randomly depending on the node order... But now it does not ignore the inner cube any more.

best Christopher

[roystgnr]

4964691 is fixing all my unit tests, thanks! It'll take me a little longer to run application tests through all the configurations I can, but those runs are already well past the cases that were breaking last week, so I'm not expecting to reopen any more bug reports on this.

Thanks for the easy testcase.

I think that super-quick fix is gratitude enough, thank you.

But if I'm wrong, could I trade the rest for a little tech support? Some of our users could really benefit from being able to generate a tetrahedralization with zero additional nodes inserted, using only a prespecified set of interior nodes. I failed to find parameters that would reliably prevent any creation of new nodes ... and it looks in the code like you're using an advancing front method that creates more nodes from the get-go, not just starting from a coarse Delaunay mesh and only refining/smoothing afterward. Good for getting a nice boundary layer for most use cases, not so good for their particular use case. Is there any way to force Netgen to output a mesh where tet vertices are limited to input nodes only?

[ChrLackner]

I think in general you can always create counterexamples where this is not possible. I agree in this case it would be, but in netgen there is kind of a tradeoff between speed, quality and elements number and netgen does not "try too hard" to make it work directly in the advancing front, but instead relatively quickly adds inner nodes in meshing, and often the optimizer can then optimize them away. Reason is mainly checking for these vis a vis rules is quite expensive. Change of this behaviour is currently not possible without hacking source code. First place to play around is the rules/tetrules.rls file where putting free tet rules in higher quality classes makes them later to be used and less prefered. This would probably lead meshes with less nodes, but become less stable. If you know which cases you are considering there can be some tuning with meshing parameters, or rewarding of swaps that remove elements, but they all come with downsides in general purpose cases I think...

one simple thing you can do is do more volume optimizations, and see if that helps.

[roystgnr]

I'll give that a try; if it doesn't help, then at least knowing that I wasn't misinterpreting the code design helps. Thanks again!

I'm afraid I may be back soon. I've got many users reporting what looks to me like evidence for UB in the static constructors in profiler.cpp and/or paje_trace.cpp ever since we pushed out our first Netgen-enabled update. I don't want to waste your time with it until I can confirm (maybe some of our other new or newly-integrated code is writing OOB earlier for Netgen to trip over later), though. I looked for obvious culprits (a Static Initialization Order Fiasco, since it works for some builds but not others?) but didn't see anything, and whatever it is seems hard to replicate.

[roystgnr]

Oh, I just realized I never reported back here. There was no bug in the code, just me failing to notice some aggressive optimization in the build system. You default USE_NATIVE_ARCH to ON in CMakeLists.txt, and I didn't think to change that in our build scripts, and yet -march=native does not play well with binary distribution when a few users have CPU architectures older than the build box architecture. We weren't crashing in Netgen constructors because of initialization order, we were crashing there because that was where we first reached CPU instructions that those users' CPUs didn't support.