1 #include <exploragram/hexdom/cavity.h> 2 #include <exploragram/hexdom/hex.h> 3 #include <exploragram/hexdom/PGP.h> 4 #include <exploragram/hexdom/basic.h> 5 #include <exploragram/hexdom/extra_connectivity.h> 6 #include <geogram/numerics/matrix_util.h> 7 #include <geogram/basic/permutation.h> 8 #include <algorithm> 9 #include <geogram/mesh/triangle_intersection.h> 10 #include <geogram/mesh/mesh_tetrahedralize.h> 11 #include <geogram/delaunay/delaunay.h> 12 13 #include <geogram/points/nn_search.h> 14 #include <geogram/points/colocate.h> 15 #include <queue> 16 17 #include <exploragram/hexdom/mesh_inspector.h> 18 #include <exploragram/hexdom/intersect_tools.h> 19 #include <exploragram/hexdom/polygon.h> 20 21 namespace GEO { 22 23 24 25 26 27 merge_hex_boundary_and_quadtri(Mesh * hex,Mesh * quad)28 void merge_hex_boundary_and_quadtri(Mesh* hex, Mesh* quad) { 29 if (hex->cells.nb() == 0) return; 30 double eps = 1e-3*get_cell_average_edge_size(hex); 31 32 // add hex surface to quadt 33 index_t off_v = quad->vertices.create_vertices(hex->vertices.nb()); 34 FOR(v, hex->vertices.nb()) X(quad)[off_v + v] = X(hex)[v]; 35 36 index_t off_f = quad->facets.create_facets(hex->facets.nb(), 4); 37 FOR(f, hex->facets.nb()) FOR(fc, 4) 38 quad->facets.set_vertex(off_f + f, (3 - fc), off_v + hex->facets.vertex(f, fc)); 39 40 41 // merge vertices from hex faces to quatri vertices 42 43 { 44 vector<index_t> to_kill(quad->vertices.nb(), 0); 45 vector<index_t> old2new(quad->vertices.nb()); 46 FOR(v, quad->vertices.nb()) old2new[v] = v; 47 48 NearestNeighborSearch_var NN = NearestNeighborSearch::create(3); 49 NN->set_points(off_v, quad->vertices.point_ptr(0)); 50 51 for (index_t v = off_v; v < quad->vertices.nb(); v++) { 52 index_t nearest = NN->get_nearest_neighbor(X(quad)[v].data()); 53 if ((X(quad)[v] - X(quad)[nearest]).length2() == 0) { 54 old2new[v] = nearest; 55 } 56 } 57 58 FOR(f, quad->facets.nb()) FOR(fv, quad->facets.nb_vertices(f)) 59 quad->facets.set_vertex(f, fv, old2new[quad->facets.vertex(f, fv)]); 60 FOR(v, quad->vertices.nb()) if (old2new[v] != v) to_kill[v] = NOT_AN_ID; 61 quad->vertices.delete_elements(to_kill); 62 } 63 64 65 66 67 // remove useless quads 68 vector<vec3> f_bary(quad->facets.nb(), vec3(0, 0, 0)); 69 70 FOR(f, quad->facets.nb()) FOR(fc, quad->facets.nb_vertices(f)) 71 f_bary[f] = f_bary[f] + (1. / quad->facets.nb_vertices(f))*quad->vertices.point(quad->facets.vertex(f, fc)); 72 73 vector<index_t> bary_old2new(quad->facets.nb()); 74 Geom::colocate((double*)(f_bary.data()), 3, quad->facets.nb(), bary_old2new, eps); 75 76 vector<index_t> to_kill(quad->facets.nb(), 0); 77 FOR(f, quad->facets.nb()) if (bary_old2new[f] != f) { 78 if (quad->facets.nb_vertices(f) != 4) continue; 79 if (quad->facets.nb_vertices(bary_old2new[f]) != 4) continue; 80 81 bool have_same_vertices = true; 82 FOR(fv1, quad->facets.nb_vertices(f)) { 83 bool isin = false; 84 FOR(fv2, quad->facets.nb_vertices(bary_old2new[f])) { 85 isin = isin || quad->facets.vertex(f, fv1) == quad->facets.vertex(bary_old2new[f], fv2); 86 } 87 have_same_vertices = have_same_vertices && isin; 88 } 89 if (!have_same_vertices) continue; 90 to_kill[f] = true; 91 to_kill[bary_old2new[f]] = true; 92 } 93 quad->facets.delete_elements(to_kill); 94 create_non_manifold_facet_adjacence(quad); 95 96 97 } 98 99 100 } 101