Thanks a lot for the detailed answer, @Algiane.

I’m only thinking of the case where the surface is represented in a CAD program – for instance as a NURBS surface or whatever – , and the level function is its signed distance function. This is the approach you suggested here , if I understand correctly. In that case, if two vertices v, w are at distance || v-w || apart from each other, the level function f can only vary by the amount |f(v) - f(w)| \le ||v - w || by the triangle inequality (it is Lipschitz-continuous with Lipschitz constant 1). So, if |f(v)| > || v-w || or |f(w)| > || v - w || then f can’t be zero on the interior of the edge. So, in that case, mmg would not have to bother with querying the level function on the interior of the edge, or subdividing that edge.

In any case, I see there are good reasons why the level function should only be given at the vertices, I’m totally fine with that.

On the practical side, here is a simple example C program. Let’s take the unit ball to serve as an example of a solid given in a CAD program. It more-or-less works acceptably - iterative refinement where in each iteration, the distance to the surface is evaluated at the vertices, and then mmg is asked to refine it according to the distance function. The only problem is that there are small areas where the number of vertices explodes, disregarding the hmin setting. Unfortunately this explosion of vertices makes the result very hard to use for FEM - some tetrahedra have a very bad quality. See the screenshot below. Is that related to the oscillation problem you mention? It would be extremely cool if you had an idea how to improve the program, or point me to the place in the mmg source code where those excessive divisions occur.

If this level set approach is unusual, what is the recommended approach to generate a tet mesh of the interior of a solid given in CAD? Let’s assume I want to stick to mmg, not use other meshers. The following already works quite ok.

```
#include <mmg/mmg3d/libmmg3d.h>
#include <math.h>
int main() {
MMG5_pMesh mmg3dMesh=NULL;
MMG5_pSol mmg3dSol=NULL;
MMG3D_Init_mesh(MMG5_ARG_start, MMG5_ARG_ppMesh, &mmg3dMesh, MMG5_ARG_ppMet, &mmg3dSol, MMG5_ARG_end);
MMG3D_Set_meshSize(mmg3dMesh, 8, 6, 0, 0, 0, 0);
// the corners of the unit cube
MMG3D_Set_vertices(mmg3dMesh, (double[]) {0,0,0,1,0,0,1,1,0,0,1,0,0,0,1,1,0,1,1,1,1,0,1,1}, NULL);
// a tet mesh of the unit cube
MMG3D_Set_tetrahedra(mmg3dMesh, (int[]) {3,6,7,8,2,5,6,8,1,2,4,5,2,4,5,8,2,3,6,8,2,3,4,8}, NULL);
// refine this tet mesh to have max edge length 0.5 to split all edges
MMG3D_Set_dparameter(mmg3dMesh, mmg3dSol, MMG3D_DPARAM_hmax, 0.5);
MMG3D_mmg3dlib(mmg3dMesh, mmg3dSol);
// let's now run an iterative procedure ten times to refine the mesh
for (int j=0;j<10;j++) {
int np;
MMG3D_Get_meshSize(mmg3dMesh, &np, NULL, NULL, NULL, NULL, NULL);
MMG3D_Set_solSize(mmg3dMesh, mmg3dSol, MMG5_Vertex, np, MMG5_Scalar);
for (int i = 1; i < np+1; i++) {
double x,y,z;
MMG3D_Get_vertex(mmg3dMesh, &x, &y, &z, NULL, NULL, NULL);
// signed distance to the unit sphere inside the cube
MMG3D_Set_scalarSol(mmg3dSol, sqrt(pow(x-0.5,2)+pow(y-0.5,2)+pow(z-0.5,2))-0.5, i);
}
MMG3D_Set_dparameter(mmg3dMesh, mmg3dSol, MMG3D_DPARAM_ls, 0);
MMG3D_Set_dparameter(mmg3dMesh, mmg3dSol, MMG3D_DPARAM_hmin, 0.05);
MMG3D_Set_dparameter(mmg3dMesh, mmg3dSol, MMG3D_DPARAM_hmax, 0.1);
MMG3D_mmg3dls(mmg3dMesh, mmg3dSol, NULL);
}
MMG3D_saveMesh(mmg3dMesh, "x.meshb");
MMG3D_Free_all(MMG5_ARG_start, MMG5_ARG_ppMesh, &mmg3dMesh, MMG5_ARG_ppMet, &mmg3dSol, MMG5_ARG_end);
return 0;
}
```

(In reality, the problem is similar but worse - there are different materials in the model, each material being a solid. I have to run this iterative procedure for each solid so that in the end each tetrahedron is inside one and only one solid. Vertices on the surface of one solid should not be destroyed by level-setting the distances to the other solids, so they have to be set “required”).

Don’t get me wrong, mmg is really doing an amazing job here. Maybe the “local explosion problem” can be solved with some Gaussian or Laplacian smoothing?!

Thanks!