Skip to content

Instantly share code, notes, and snippets.

@Mathias-Fuchs Mathias-Fuchs/main.c

Last active Jun 3, 2020
Embed
What would you like to do?
Meshing with anisotropic metric in MMG2D without a msh/sol file
#include <stdio.h>
#include <stdlib.h>
#include <mmg/mmg2d/libmmg2d.h>
#include <math.h>
int main() {
MMG5_pMesh mesh2 = NULL;
MMG5_pSol met2 = NULL;
int ier = 0;
MMG2D_Init_mesh(MMG5_ARG_start, MMG5_ARG_ppMesh, &mesh2, MMG5_ARG_ppMet, &met2, MMG5_ARG_end);
MMG2D_Set_iparameter(mesh2, met2, MMG2D_IPARAM_verbose, 10);
// 4 vertices, no tris, no quads, no edges so far
MMG2D_Set_meshSize(mesh2, 4, 0, 0, 0);
// a square [0, 2 * M_PI] x [-M_PI / 2, M_PI / 2] for spherical coordinates in usual lon/lat way
double v[8] = {0, -M_PI/2, 2 * M_PI, -M_PI/2, 2 * M_PI, M_PI/2, 0, M_PI/2};
MMG2D_Set_vertices(mesh2, v, NULL);
MMG2D_Set_dparameter(mesh2, met2, MMG2D_DPARAM_hsiz, 0.01);
// generate a regular fine mesh of the square in meshing mode
ier = MMG2D_mmg2dmesh(mesh2, met2);
if (ier) {
fprintf(stderr, "error %i during meshing.\n", ier);
exit(1);
}
// save the "computational geometry" mesh
MMG2D_saveMshMesh(mesh2, NULL, "cg.msh");
// remesh with anisotropic metric
int np, nt, nquad, na;
MMG2D_Get_meshSize(mesh2, &np, &nt, &nquad, &na);
double* verts = (double*) malloc(2 * np * sizeof(double));
MMG2D_Get_vertices(mesh2, verts, NULL, NULL, NULL);
MMG2D_Set_solSize(mesh2, met2, MMG5_Vertex, np, MMG5_Tensor);
for (int i = 0; i < np; i++) {
// latitude
double y = verts[2 * i + 1];
// metric on the sphere, see standard textbooks on elementary differential geometry
// Gaussian fundamental quantities
double E = cos(y) * cos(y);
double F = 0.0;
double G = 1.0;
// the following factor was found by trial and error
// why does it not work with a factor of one?
// how to compute the factor in real life?
double factor = 1000.0;
// we add a small amount to the 1-1-entry to make the metric non-singular everywhere
MMG2D_Set_tensorSol(met2, factor * E + 0.1, factor * F, factor * G, i+1);
}
// disable hsiz because it is incompatible with an anisotropic metric
MMG2D_Set_dparameter(mesh2, met2, MMG2D_DPARAM_hsiz, -1);
// set gradation to a not too restrictive value
MMG2D_Set_dparameter(mesh2, met2, MMG2D_DPARAM_hgrad, 1.5);
// do it, remesh!
ier = MMG2D_mmg2dlib(mesh2, met2);
if (ier != 0) {
fprintf(stdout, "error %i \n", ier);
exit(1);
}
MMG2D_saveMshMesh(mesh2, NULL, "out2.msh");
// map to 3d and save as obj file
MMG2D_Get_meshSize(mesh2, &np, &nt, &nquad, &na);
verts = realloc(verts, 2 * np * sizeof(double));
int* tris = (int*) malloc(3 * nt * sizeof(int));
MMG2D_Get_vertices(mesh2, verts, NULL, NULL, NULL);
MMG2D_Get_triangles(mesh2, tris, NULL, NULL);
FILE *fp = fopen("sphere.obj", "w+");
for (int i = 0; i < np; i++) {
double x = verts[2 * i];
double y = verts[2 * i + 1];
double newx = cos(x) * cos(y);
double newy = sin(x) * cos(y);
double newz = sin(y);
fprintf(fp, "v %f %f %f\n", newx, newy, newz);
}
for (int i = 0; i < nt; i++)
fprintf(fp, "f %i %i %i\n", tris[3 * i], tris[3 * i +1], tris[3 * i + 2]);
fclose(fp);
free(verts);
free(tris);
MMG2D_Free_all(MMG5_ARG_start, MMG5_ARG_ppMesh, &mesh2, MMG5_ARG_ppMet, &met2, MMG5_ARG_end);
return 0;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.