Hello,
the following program:
- defines a mesh of the unit square with two triangles
- remeshes the line x=0.3 with the level set function
- then remeshes the line x=0.7 with the level set function
- tries to set the edges from the first remeshing as required
but … these required edges are blatantly ignored.
It is a minimal example boiled down to 2d of the strategy discussed here BUT please don’t just dismiss it as “it’s not mmg’s job to do meshing, only remeshing”. If there is an API function called “MMG2D_Set_requiredEdge” and mmg does not respect it then it should at least print a warning, regardless of what one is ultimately trying to do with mmg.
#include <mmg/mmg2d/libmmg2d.h>
#include <stdlib.h> // for malloc
#include <stdio.h> // for stdout
#include <math.h> // for fabs
int main()
{
//
// step 1 :: mesh the line x=0.3
//
MMG5_pMesh mmg2dMesh = NULL;
MMG5_pSol mmg2dSol = NULL;
MMG2D_Init_mesh(MMG5_ARG_start, MMG5_ARG_ppMesh, &mmg2dMesh, MMG5_ARG_ppMet, &mmg2dSol, MMG5_ARG_end);
MMG2D_Set_meshSize(mmg2dMesh, 4, 2, 0, 0);
MMG2D_Set_vertices(mmg2dMesh, (double[]) { 0, 0, 1, 0, 1, 1, 0, 1 }, NULL);
MMG2D_Set_triangles(mmg2dMesh, (int[]) { 1, 2, 4, 2, 3, 4 }, NULL);
MMG2D_Set_solSize(mmg2dMesh, mmg2dSol, MMG5_Vertex, 4, MMG5_Scalar);
for (int i = 1; i < 5; i++) {
double x, y;
MMG2D_Get_vertex(mmg2dMesh, &x, &y, NULL, NULL, NULL);
MMG2D_Set_scalarSol(mmg2dSol, x - 0.3, i);
}
MMG2D_Set_iparameter(mmg2dMesh, mmg2dSol, MMG2D_IPARAM_iso, 1);
MMG2D_Set_iparameter(mmg2dMesh, mmg2dSol, MMG2D_IPARAM_verbose, 10);
MMG2D_Set_iparameter(mmg2dMesh, mmg2dSol, MMG2D_IPARAM_angle, 1);
MMG2D_Set_dparameter(mmg2dMesh, mmg2dSol, MMG2D_DPARAM_ls, 0);
MMG2D_Set_dparameter(mmg2dMesh, mmg2dSol, MMG2D_DPARAM_hmin, 0.025);
MMG2D_Set_dparameter(mmg2dMesh, mmg2dSol, MMG2D_DPARAM_hmax, 0.05);
MMG2D_mmg2dls(mmg2dMesh, mmg2dSol, NULL);
MMG2D_saveMesh(mmg2dMesh, "first.mesh");
//
// step 2: remesh also the line x 0.7 and require also the edges from the last step
//
int np, ne, na;
MMG2D_Get_meshSize(mmg2dMesh, &np, &ne, NULL, &na);
int* verts = malloc(2 * np * sizeof(double));
int* elements = malloc(3 * ne * sizeof(int));
int* edges = malloc(2 * na * sizeof(int));
MMG2D_Get_vertices(mmg2dMesh, verts, NULL, NULL, NULL);
MMG2D_Get_triangles(mmg2dMesh, elements, NULL, NULL);
MMG2D_Get_edges(mmg2dMesh, edges, NULL, NULL, NULL);
MMG2D_Free_all(MMG5_ARG_start, MMG5_ARG_ppMesh, &mmg2dMesh, MMG5_ARG_ppMet, &mmg2dSol, MMG5_ARG_end);
MMG2D_Init_mesh(MMG5_ARG_start, MMG5_ARG_ppMesh, &mmg2dMesh, MMG5_ARG_ppMet, &mmg2dSol, MMG5_ARG_end);
MMG2D_Set_meshSize(mmg2dMesh, np, ne, 0, na);
MMG2D_Set_vertices(mmg2dMesh, verts, NULL);
MMG2D_Set_triangles(mmg2dMesh, elements, NULL);
MMG2D_Set_edges(mmg2dMesh, edges, NULL);
for (int i = 0; i < na; i++) {
MMG2D_Set_requiredEdge(mmg2dMesh, i + 1);
double fx, fy, tx, ty;
MMG2D_GetByIdx_vertex(mmg2dMesh, &fx, &fy, NULL, NULL, NULL, edges[2 * i]);
MMG2D_GetByIdx_vertex(mmg2dMesh, &tx, &ty, NULL, NULL, NULL, edges[2 * i + 1]);
// just to make sure check that the edges at x=0.3 are really set to required
if (fabs(fx - 0.3) < 0.0001 && fabs(tx - 0.3) < 0.0001)
fprintf(stdout, "Edge from %.3f %.3f to %.3f %.3f is set to required.\n", fx, fy, tx, ty);
}
MMG2D_Set_solSize(mmg2dMesh, mmg2dSol, MMG5_Vertex, np, MMG5_Scalar);
for (int i = 1; i < np + 1; i++) {
double x, y;
MMG2D_Get_vertex(mmg2dMesh, &x, &y, NULL, NULL, NULL);
MMG2D_Set_scalarSol(mmg2dSol, x - 0.7, i);
}
MMG2D_Set_iparameter(mmg2dMesh, mmg2dSol, MMG2D_IPARAM_iso, 1);
MMG2D_Set_iparameter(mmg2dMesh, mmg2dSol, MMG2D_IPARAM_verbose, 10);
MMG2D_Set_iparameter(mmg2dMesh, mmg2dSol, MMG2D_IPARAM_angle, 1);
MMG2D_Set_dparameter(mmg2dMesh, mmg2dSol, MMG2D_DPARAM_ls, 0);
MMG2D_Set_dparameter(mmg2dMesh, mmg2dSol, MMG2D_DPARAM_hmin, 0.025);
MMG2D_Set_dparameter(mmg2dMesh, mmg2dSol, MMG2D_DPARAM_hmax, 0.05);
MMG2D_mmg2dls(mmg2dMesh, mmg2dSol, NULL);
MMG2D_saveMesh(mmg2dMesh, "second.mesh");
free(verts);
free(elements);
free(edges);
MMG2D_Free_all(MMG5_ARG_start, MMG5_ARG_ppMesh, &mmg2dMesh, MMG5_ARG_ppMet, &mmg2dSol, MMG5_ARG_end);
return 0;
}
Please someone tell me what’s going wrong here.
Mathias
).