18 template <
class T1,
class T2>
21 det = +
t(0, 0) *
t(1, 1) *
t(2, 2) +
t(1, 0) *
t(2, 1) *
t(0, 2) +
22 t(2, 0) *
t(0, 1) *
t(1, 2) -
t(0, 0) *
t(2, 1) *
t(1, 2) -
23 t(2, 0) *
t(1, 1) *
t(0, 2) -
t(1, 0) *
t(0, 1) *
t(2, 2);
39 std::map<EntityHandle, unsigned long> &moab_tetgen_map,
40 std::map<unsigned long, EntityHandle> &tetgen_moab_map,
50 "TETGEN_MARKER", 1, MB_TYPE_INTEGER, th_marker,
51 MB_TAG_CREAT | MB_TAG_SPARSE, &def_marker);
56 Range points = ents.subset_by_dimension(0);
57 in.numberofpoints = points.size();
58 if (points.size() > 0) {
59 in.pointlist =
new double[in.numberofpoints * 3];
60 in.pointmarkerlist =
new int[in.numberofpoints];
69 for (
int ii = 0; it != points.end(); it++, ii++) {
71 tetgen_moab_map[iii] = *it;
72 moab_tetgen_map[*it] = iii;
76 Range tets = ents.subset_by_type(MBTET);
77 in.numberoftetrahedra = tets.size();
78 if (in.numberoftetrahedra > 0) {
79 in.tetrahedronlist =
new int[4 * ents.subset_by_type(MBTET).size()];
81 for (
int ii = 0; it != tets.end(); it++, ii++) {
84 CHKERR m_field.
get_moab().get_connectivity(*it, conn, num_nodes,
true);
89 CHKERR m_field.tag_get_data(
th, conn, num_nodes, coords);
103 for (
int nn = 0; nn != 4; nn++) {
104 jac(
i,
j) += t_coords(
i) * t_diff_n(
j);
112 "Negative volume det = %6.4e", det);
118 for (
int nn = 0; nn != 4; nn++) {
119 if (moab_tetgen_map.find(conn[nn]) == moab_tetgen_map.end()) {
121 "data inconsistency between TetGen and MoAB");
123 in.tetrahedronlist[4 * ii + nn] =
129 Range tris = ents.subset_by_type(MBTRI);
130 in.numberoftrifaces = tris.size();
131 if (in.numberoftrifaces) {
132 in.trifacelist =
new int[3 * in.numberoftrifaces];
133 in.trifacemarkerlist =
new int[in.numberoftrifaces];
136 in.trifacemarkerlist);
138 for (
int ii = 0; it != tris.end(); it++, ii++) {
141 CHKERR m_field.
get_moab().get_connectivity(*it, conn, num_nodes,
true);
142 int order[] = {0, 1, 2};
144 CHKERR m_field.
get_moab().get_adjacencies(&*it, 1, 3,
true, adj_tets);
145 adj_tets = intersect(adj_tets, tets);
146 if (adj_tets.size() == 1) {
150 CHKERR m_field.
get_moab().side_number(adj_tets[0], *it, side_number,
159 for (
int nn = 0; nn < 3; nn++) {
160 if (moab_tetgen_map.find(conn[
order[nn]]) == moab_tetgen_map.end()) {
162 "data inconsistency between TetGen and MoAB");
164 in.trifacelist[3 * ii + nn] =
170 Range edges = ents.subset_by_type(MBEDGE);
171 in.numberofedges = edges.size();
172 if (in.numberofedges > 0) {
173 in.edgelist =
new int[2 * in.numberofedges];
174 in.edgemarkerlist =
new int[in.numberofedges];
176 CHKERR m_field.
get_moab().tag_get_data(th_marker, edges, in.edgemarkerlist);
178 for (
int ii = 0; it != edges.end(); it++, ii++) {
181 CHKERR m_field.
get_moab().get_connectivity(*it, conn, num_nodes,
true);
184 for (
int nn = 0; nn < 2; nn++) {
185 if (moab_tetgen_map.find(conn[nn]) == moab_tetgen_map.end()) {
187 "data inconsistency between TetGen and MoAB");
189 in.edgelist[2 * ii + nn] = moab_tetgen_map[conn[nn]] >>
MB_TYPE_WIDTH;
200 std::map<int, Range> &type_ents) {
206 in.pointparamlist =
new tetgenio::pointparam[in.numberofpoints];
208 std::map<int, Range>::iterator mit = type_ents.begin();
209 for (; mit != type_ents.end(); mit++) {
210 if (mit->first < 0 && mit->first > 3) {
212 "Wrong TetGen point type");
214 Range::iterator it = mit->second.begin();
215 for (; it != mit->second.end(); it++) {
216 moabTetGen_Map::iterator miit = moab_tetgen_map.find(*it);
217 if (miit == moab_tetgen_map.end()) {
219 "Data inconsistency between TetGen and MoAB");
223 in.pointparamlist[id].uv[0] = 0;
224 in.pointparamlist[id].uv[1] = 0;
225 in.pointparamlist[id].type = mit->first;
226 in.pointparamlist[id].tag = m_field.
get_moab().id_from_handle(*it) + 1;
247 std::map<EntityHandle, unsigned long> &moab_tetgen_map,
248 std::map<unsigned long, EntityHandle> &tetgen_moab_map,
249 Range *ents,
bool id_in_tags,
bool error_if_created,
250 bool assume_first_nodes_the_same, Tag
th) {
258 "TETGEN_MARKER", 1, MB_TYPE_INTEGER, th_marker,
259 MB_TAG_CREAT | MB_TAG_SPARSE, &def_marker);
264 for (; ii < out.numberofpoints; ii++) {
265 if (ii < in.numberofpoints) {
266 bool mem_the_same = memcmp(&in.pointlist[3 * ii], &out.pointlist[3 * ii],
267 3 *
sizeof(
double)) == 0;
268 if (assume_first_nodes_the_same || mem_the_same) {
270 if (tetgen_moab_map.find(iii) != tetgen_moab_map.end()) {
274 1, &out.pointlist[3 * ii]);
277 &out.pointlist[3 * ii]);
280 th_marker, &tetgen_moab_map[iii], 1, &out.pointmarkerlist[ii]);
284 "data inconsistency between TetGen and MoAB");
289 if (out.pointparamlist[ii].tag > 0) {
292 MBVERTEX, in.pointparamlist[ii].tag - 1, node);
293 if (moab_tetgen_map.find(node) != moab_tetgen_map.end()) {
298 if (error_if_created) {
300 "node should not be created");
305 ReadUtilIface *iface;
309 vector<double *> arrays;
311 CHKERR iface->get_node_coords(3, num_nodes, 0, startv, arrays);
312 Range verts(startv, startv + num_nodes - 1);
313 int ii = in.numberofpoints;
314 for (Range::iterator vit = verts.begin(); vit != verts.end(); vit++, ii++) {
315 for (
int nn = 0; nn != 3; nn++)
316 arrays[nn][ii - in.numberofpoints] = out.pointlist[3 * ii + nn];
323 th_marker, verts, &out.pointmarkerlist[in.numberofpoints]);
326 th, verts, &out.pointlist[3 * in.numberofpoints]);
329 std::vector<int> tetgen_ii;
332 std::vector<EntityHandle> conn_seq_tet;
333 conn_seq_tet.reserve(4 * out.numberoftetrahedra);
334 tetgen_ii.reserve(out.numberoftetrahedra);
335 conn_seq_tet.clear();
338 for (; ii < out.numberoftetrahedra; ii++) {
340 if (ii < in.numberoftetrahedra) {
341 if (memcmp(&in.tetrahedronlist[4 * ii], &out.tetrahedronlist[4 * ii],
342 4 *
sizeof(
int)) == 0) {
343 if (tetgen_moab_map.find(iii) != tetgen_moab_map.end()) {
345 ents->insert(tetgen_moab_map[iii]);
349 "data inconsistency between TetGen and MoAB");
354 for (
int nn = 0; nn < 4; nn++) {
355 int nnn = out.tetrahedronlist[4 * ii + nn];
356 if (tetgen_moab_map.find(MBVERTEX | (nnn <<
MB_TYPE_WIDTH)) ==
357 tetgen_moab_map.end()) {
359 "data inconsistency between TetGen and MoAB");
361 conn[nn] = tetgen_moab_map.at(MBVERTEX | (nnn <<
MB_TYPE_WIDTH));
377 CHKERR m_field.
get_moab().get_adjacencies(conn, 4, 3,
false, tets);
378 bool tet_found =
false;
379 for (
auto tet : tets) {
382 CHKERR m_field.
get_moab().get_connectivity(tet, tet_conn, num_nodes,
384 const EntityHandle *p = std::find(tet_conn, &tet_conn[4], conn[0]);
385 if (p != &tet_conn[4]) {
386 int s = std::distance(tet_conn, p);
388 for (;
n != 4; ++
n) {
389 const int idx[] = {0, 1, 2, 3, 0, 1, 2, 3};
390 if (tet_conn[idx[s +
n]] != conn[
n])
393 if (
n == 4 && !tet_found) {
394 moab_tetgen_map[tet] = iii;
395 tetgen_moab_map[iii] = tet;
401 "More that one tet with the same connectivity");
408 for (
int nn = 0; nn != 4; nn++) {
409 conn_seq_tet.push_back(conn[nn]);
411 tetgen_ii.push_back(ii);
416 if (!conn_seq_tet.empty()) {
419 int num_el = tetgen_ii.size();
420 CHKERR iface->get_element_connect(num_el, 4, MBTET, 0, starte, conn);
421 std::copy(conn_seq_tet.begin(), conn_seq_tet.end(), conn);
422 CHKERR iface->update_adjacencies(starte, num_el, 4, conn);
423 new_tets =
Range(starte, starte + num_el - 1);
424 std::vector<int>::iterator ii_it = tetgen_ii.begin();
426 for (Range::iterator tit = new_tets.begin(); tit != new_tets.end();
427 tit++, ii_it++, ii++) {
429 moab_tetgen_map[*tit] = iii;
430 tetgen_moab_map[iii] = *tit;
433 ents->merge(new_tets);
440 std::map<EntityHandle, unsigned long> &moab_tetgen_map,
441 std::map<unsigned long, EntityHandle> &tetgen_moab_map,
443 bool error_if_created,
444 bool assume_first_nodes_the_same, Tag
th) {
448 CHKERR outData(in, out, moab_tetgen_map, tetgen_moab_map, &ents, id_in_tags,
449 error_if_created, assume_first_nodes_the_same,
th);
451 ents.subset_by_type(MBTET),
bit);
455 std::vector<std::pair<Range, int>> &markers, tetgenio &in,
456 std::map<EntityHandle, unsigned long> &moab_tetgen_map,
457 std::map<unsigned long, EntityHandle> &tetgen_moab_map) {
461 in.numberoffacets = markers.size();
462 in.facetlist =
new tetgenio::facet[in.numberoffacets];
463 in.facetmarkerlist =
new int[in.numberoffacets];
464 std::vector<std::pair<Range, int>>::iterator mit = markers.begin();
465 for (
int ii = 0; mit != markers.end(); mit++, ii++) {
466 in.facetmarkerlist[ii] = mit->second;
467 Range &faces = mit->first;
468 tetgenio::facet *
f = &(in.facetlist[ii]);
469 f->numberofpolygons = faces.size();
470 f->polygonlist =
new tetgenio::polygon[
f->numberofpolygons];
472 for (
int dd = 3;
dd >= 0;
dd--) {
473 Range dd_faces = faces.subset_by_dimension(
dd);
474 if (dd_faces.empty())
476 Range::iterator it = dd_faces.begin();
477 for (; it != dd_faces.end(); it++, jj++) {
480 tetgenio::polygon *p = &(
f->polygonlist[jj]);
483 p->numberofvertices = 1;
488 m_field.
get_moab().get_connectivity(*it, conn, num_nodes,
true);
490 p->numberofvertices = num_nodes;
493 p->vertexlist =
new int[p->numberofvertices];
494 for (
int nn = 0; nn < p->numberofvertices; nn++) {
495 if (moab_tetgen_map.find(conn[nn]) == moab_tetgen_map.end()) {
497 "data inconsistency between TetGen and MoAB");
499 p->vertexlist[nn] = moab_tetgen_map[conn[nn]] >>
MB_TYPE_WIDTH;
504 f->numberofholes = 0;
510 std::map<EntityHandle, unsigned long> &tetgen_moab_map, tetgenio &out,
518 "TETGEN_MARKER", 1, MB_TYPE_INTEGER, th_marker,
519 MB_TAG_CREAT | MB_TAG_SPARSE, &def_marker);
522 for (; ii < out.numberoftrifaces; ii++) {
524 if (out.trifacemarkerlist[ii] == 0) {
529 for (
int nn = 0; nn < 3; nn++) {
530 int iii = MBVERTEX | (out.trifacelist[3 * ii + nn] <<
MB_TYPE_WIDTH);
531 if (tetgen_moab_map.find(iii) == tetgen_moab_map.end()) {
533 "data inconsistency between TetGen and MoAB");
535 conn[nn] = tetgen_moab_map[iii];
539 rval = m_field.
get_moab().get_adjacencies(conn, 3, 2,
false, face);
541 face = face.subset_by_type(MBTRI);
542 if (face.size() != 1) {
544 "data inconsistency between TetGen and MoAB, %u", face.size());
548 if (ents_map != NULL)
549 (*ents_map)[out.trifacemarkerlist[ii]].merge(face);
550 rval = m_field.
get_moab().tag_set_data(th_marker, &*face.begin(), 1,
551 &out.trifacemarkerlist[ii]);
557 std::vector<std::pair<EntityHandle, int>> ®ions, tetgenio &in, Tag
th) {
561 in.numberofregions = regions.size();
562 in.regionlist =
new double[5 * in.numberofregions];
564 std::vector<std::pair<EntityHandle, int>>::iterator it = regions.begin();
565 for (
int ii = 0; it != regions.end(); it++, ii++) {
570 rval = m_field.
get_moab().tag_get_data(
th, &it->first, 1, coords);
573 rval = m_field.
get_moab().get_coords(&it->first, 1, coords);
581 m_field.
get_moab().get_connectivity(it->first, conn, num_nodes,
true);
585 rval = m_field.
get_moab().tag_get_data(
th, conn, num_nodes, _coords);
588 rval = m_field.
get_moab().get_coords(conn, num_nodes, _coords);
591 coords[0] = (_coords[0] + _coords[3] + _coords[6] + _coords[9]) / 4.;
592 coords[1] = (_coords[1] + _coords[4] + _coords[7] + _coords[10]) / 4.;
593 coords[2] = (_coords[2] + _coords[5] + _coords[8] + _coords[11]) / 4.;
598 for (
int nn = 0; nn < 3; nn++) {
599 in.regionlist[kk++] = coords[nn];
601 in.regionlist[kk++] = it->second;
602 in.regionlist[kk++] = it->second;
607 std::map<EntityHandle, unsigned long> &tetgen_moab_map, tetgenio &out,
611 int nbattributes = out.numberoftetrahedronattributes;
612 if (nbattributes == 0) {
614 "tetgen has no regions attributes");
617 rval = m_field.
get_moab().tag_get_handle(
"TETGEN_REGION", th_region);
618 if (
rval == MB_SUCCESS) {
622 double def_marker = 0;
624 "TETGEN_REGION", nbattributes, MB_TYPE_DOUBLE, th_region,
625 MB_TAG_CREAT | MB_TAG_SPARSE, &def_marker);
627 for (; ii < out.numberoftetrahedra; ii++) {
629 for (; jj < nbattributes; jj++) {
630 double id = out.tetrahedronattributelist[ii * nbattributes + jj];
632 if (tetgen_moab_map.find(iii) == tetgen_moab_map.end()) {
634 "data inconsistency between TetGen and MoAB");
640 if (ents_map != NULL)
641 (*ents_map)[id].insert(ent);
650 a.parse_commandline(switches);
651 tetrahedralize(&
a, &in, &out);
656 if (!in.load_poly(file_name)) {
658 "can not read TetGen poly file");
666 double *pa = &coords[0];
667 double *pb = &coords[3];
668 double *pc = &coords[6];
669 double *pd = &coords[9];
670 double adx = pa[0] - pd[0];
671 double bdx = pb[0] - pd[0];
672 double cdx = pc[0] - pd[0];
673 double ady = pa[1] - pd[1];
674 double bdy = pb[1] - pd[1];
675 double cdy = pc[1] - pd[1];
676 double adz = pa[2] - pd[2];
677 double bdz = pb[2] - pd[2];
678 double cdz = pc[2] - pd[2];
679 double v = adx * (bdy * cdz - bdz * cdy) + bdx * (cdy * adz - cdz * ady) +
680 cdx * (ady * bdz - adz * bdy);
681 double l = sqrt(pow(pa[0] - pb[0], 2) + pow(pa[1] - pb[1], 2) +
682 pow(pa[2] - pb[2], 2)) +
683 sqrt(pow(pa[0] - pc[0], 2) + pow(pa[1] - pc[1], 2) +
684 pow(pa[2] - pc[2], 2)) +
685 sqrt(pow(pa[0] - pd[0], 2) + pow(pa[1] - pd[1], 2) +
686 pow(pa[2] - pd[2], 2)) +
687 sqrt(pow(pb[0] - pc[0], 2) + pow(pb[1] - pc[1], 2) +
688 pow(pb[2] - pc[2], 2)) +
689 sqrt(pow(pb[0] - pd[0], 2) + pow(pb[1] - pd[1], 2) +
690 pow(pb[2] - pd[2], 2)) +
691 sqrt(pow(pc[0] - pd[0], 2) + pow(pc[1] - pd[1], 2) +
692 pow(pc[2] - pd[2], 2));
694 *result = fabs(
v / pow(
l, 3)) <
eps ? true :
false;
698 std::vector<Range> &sorted,
699 const double eps, Tag
th) {
707 Range noplanar_to_anyother;
708 std::vector<Range>::iterator vit = sorted.begin();
716 CHKERR skin.find_skin(0, *vit,
false, skin_edges);
719 Range skin_edges_nodes;
720 CHKERR m_field.
get_moab().get_connectivity(skin_edges, skin_edges_nodes,
724 Range skin_edges_tris;
726 skin_edges, 2,
false, skin_edges_tris, moab::Interface::UNION);
728 Range inner_tris = intersect(skin_edges_tris, *vit);
729 Range outer_tris = intersect(skin_edges_tris, tris);
734 Range::iterator tit = outer_tris.begin();
735 for (; tit != outer_tris.end(); tit++) {
737 CHKERR m_field.
get_moab().get_connectivity(&*tit, 1, tit_conn,
true);
738 tit_conn = subtract(tit_conn, skin_edges_nodes);
739 if (tit_conn.empty()) {
740 coplanar_tris.insert(*tit);
741 noplanar_to_anyother.erase(*tit);
746 &*tit, 1, 1,
false, tit_edges, moab::Interface::UNION);
747 tit_edges = intersect(tit_edges, skin_edges);
748 if (tit_edges.size() != 1) {
750 "data inconsistency");
754 tit_edges, 2,
false, inner_tri, moab::Interface::UNION);
755 inner_tri = intersect(inner_tri, inner_tris);
756 if (inner_tri.size() != 1) {
758 "data inconsistency");
764 *inner_tri.begin(), inner_tri_conn, num_nodes,
true);
779 coplanar_tris.insert(*tit);
780 noplanar_to_anyother.erase(*tit);
783 noplanar_to_anyother.insert(*tit);
788 vit->merge(coplanar_tris);
789 tris = subtract(tris, *vit);
792 vit = sorted.begin();
797 }
while (vit != sorted.end());
799 if (noplanar_to_anyother.empty()) {
803 seed.insert(noplanar_to_anyother[0]);
804 tris.erase(noplanar_to_anyother[0]);
805 sorted.push_back(seed);
813 Range &tris, std::vector<std::vector<Range>> &sorted,
const double eps) {
819 seed.insert(tris[0]);
821 std::vector<Range> vec_seed;
822 vec_seed.push_back(seed);
823 sorted.push_back(vec_seed);
826 std::vector<Range> &vec = sorted.back();
832 seed.insert(tris[0]);
834 std::vector<Range> vec_seed;
835 vec_seed.push_back(seed);
836 sorted.push_back(vec_seed);
844 Range *not_reducable_nodes,
845 const double eps, Tag
th) {
851 "no ents to build polygon");
859 CHKERR skin.find_skin(0, ents,
false, skin_edges);
861 std::vector<EntityHandle> polygon_nodes;
864 seen_edges.insert(seed);
865 skin_edges.erase(seed);
868 CHKERR m_field.
get_moab().get_connectivity(seed, conn, num_nodes,
true);
869 polygon_nodes.push_back(conn[0]);
875 CHKERR m_field.
get_moab().get_adjacencies(&last_node, 1, 1,
false,
877 adj_edges = intersect(adj_edges, skin_edges);
878 if (adj_edges.size() == 0) {
881 if (adj_edges.size() != 1) {
883 "should be only one edge");
885 seen_edges.insert(adj_edges[0]);
886 skin_edges.erase(adj_edges[0]);
887 CHKERR m_field.
get_moab().get_connectivity(adj_edges[0], conn, num_nodes,
889 EntityHandle add_node = (last_node == conn[0]) ? conn[1] : conn[0];
890 polygon_nodes.push_back(add_node);
896 std::vector<EntityHandle>::iterator pit = polygon_nodes.begin();
898 for (; pit != polygon_nodes.end();) {
899 if (not_reducable_nodes != NULL) {
900 if (not_reducable_nodes->find(*pit) != not_reducable_nodes->end()) {
906 if (pit == polygon_nodes.begin()) {
907 mm = polygon_nodes.back();
913 if (polygon_nodes.back() == mc) {
914 mp = polygon_nodes[0];
928 cblas_daxpy(3, -1, &coords[3 * 1], 1, &coords[3 * 0], 1);
930 cblas_daxpy(3, -1, &coords[3 * 1], 1, &coords[3 * 2], 1);
934 double l0 = cblas_dnrm2(3, &coords[3 * 0], 1);
935 cblas_dgemv(CblasRowMajor, CblasNoTrans, 3, 3, 1. / l0, spin, 3,
936 &coords[3 * 2], 1, 0., &coords[3 * 1], 1);
937 double dot = cblas_dnrm2(3, &coords[3 * 1], 1);
940 polygon_nodes.erase(pit);
941 pit = polygon_nodes.begin();
949 Range existing_polygon;
951 &polygon_nodes[0], polygon_nodes.size(), 2,
true, existing_polygon);
952 if (existing_polygon.empty()) {
954 CHKERR m_field.
get_moab().create_element(MBPOLYGON, &polygon_nodes[0],
955 polygon_nodes.size(), polygon);
956 polygons.insert(polygon);
958 polygons.merge(existing_polygon);