7#define CutMeshFunctionBegin \
9 MOFEM_LOG_CHANNEL("WORLD"); \
10 MOFEM_LOG_FUNCTION(); \
11 MOFEM_LOG_TAG("WORLD", "CutMesh");
24 : cOre(const_cast<
Core &>(core)) {
49 double *shift,
double *origin,
51 const std::string save_mesh) {
53 moab::Interface &moab = m_field.
get_moab();
56 std::map<EntityHandle, EntityHandle> verts_map;
57 for (Range::const_iterator tit =
surface.begin(); tit !=
surface.end();
61 CHKERR moab.get_connectivity(*tit, conn, num_nodes,
true);
64 CHKERR moab.tag_get_data(
th, conn, num_nodes, &coords(0, 0));
66 CHKERR moab.get_coords(conn, num_nodes, &coords(0, 0));
69 for (
int nn = 0; nn != num_nodes; nn++) {
70 if (verts_map.find(conn[nn]) != verts_map.end()) {
71 new_verts[nn] = verts_map[conn[nn]];
74 ublas::matrix_row<MatrixDouble> mr(coords, nn);
77 3, ublas::shallow_array_adaptor<double>(3, origin));
81 3, 3, ublas::shallow_array_adaptor<double>(9, transform));
82 mr = prod(mat_transform, mr);
85 3, ublas::shallow_array_adaptor<double>(3, origin));
90 ublas::matrix_row<MatrixDouble> mr(coords, nn);
92 3, ublas::shallow_array_adaptor<double>(3, shift));
95 CHKERR moab.create_vertex(&coords(nn, 0), new_verts[nn]);
96 verts_map[conn[nn]] = new_verts[nn];
100 CHKERR moab.create_element(MBTRI, new_verts, num_nodes, ele);
103 if (!save_mesh.empty())
133 const double rel_tol,
134 const double abs_tol,
138 moab::Interface &moab = m_field.
get_moab();
153 const Range fixed_edges,
154 const double rel_tol,
155 const double abs_tol,
158 moab::Interface &moab = m_field.
get_moab();
162 map<EntityHandle, double> map_verts_length;
164 for (
auto f : surface_edges) {
167 CHKERR moab.get_connectivity(f, conn_skin, num_nodes,
true);
170 CHKERR moab.tag_get_data(
th, conn_skin, num_nodes, &coords_skin[0]);
172 CHKERR moab.get_coords(conn_skin, num_nodes, &coords_skin[0]);
174 &coords_skin[0], &coords_skin[1], &coords_skin[2]);
176 &coords_skin[3], &coords_skin[4], &coords_skin[5]);
178 t_skin_delta(
i) = t_n1(
i) - t_n0(
i);
179 const double skin_edge_length = sqrt(t_skin_delta(
i) * t_skin_delta(
i));
180 for (
int nn = 0; nn != 2; ++nn) {
181 auto it = map_verts_length.find(conn_skin[nn]);
182 if (it != map_verts_length.end())
183 it->second = std::min(it->second, skin_edge_length);
185 map_verts_length[conn_skin[nn]] = skin_edge_length;
189 for (
auto m : map_verts_length) {
193 CHKERR moab.tag_get_data(
th, &
m.first, 1, &t_n(0));
195 CHKERR moab.get_coords(&
m.first, 1, &t_n(0));
197 double min_dist = rel_tol *
m.second;
200 &t_n(0), 1, fixed_edges, &min_dist, &t_min_coords(0));
202 if (min_dist < rel_tol *
m.second || min_dist < abs_tol) {
204 MOFEM_LOG(
"WORLD", Sev::noisy) <<
"Snap " << min_dist;
206 CHKERR moab.tag_set_data(
th, &
m.first, 1, &t_min_coords(0));
208 CHKERR moab.set_coords(&
m.first, 1, &t_min_coords(0));
217 moab::Interface &moab = m_field.
get_moab();
219 treeSurfPtr = boost::shared_ptr<OrientedBoxTreeTool>(
220 new OrientedBoxTreeTool(&moab,
"ROOTSETSURF",
true));
227 const double tol_geometry,
const double tol_cut_close,
229 const bool update_meshsets,
const bool debug) {
231 moab::Interface &moab = m_field.
get_moab();
251 ->updateAllMeshsetsByEntitiesChildren(cut_bit);
257 CHKERR SaveData(moab)(
"out_cut_new_fixed_edges.vtk", *fixed_edges);
264 const double tol_trim_close,
267 const bool update_meshsets,
270 moab::Interface &moab = m_field.
get_moab();
289 ->updateAllMeshsetsByEntitiesChildren(trim_bit);
301 trim_bit,
BitRefLevel().set(), MBEDGE, bit2_edges);
303 intersect(*fixed_edges, bit2_edges));
310 int &first_bit, Tag
th,
const double tol_geometry,
311 const double tol_cut_close,
const double tol_trim_close,
Range *fixed_edges,
312 Range *corner_nodes,
const bool update_meshsets,
const bool debug) {
316 std::vector<BitRefLevel> bit_levels;
318 auto get_back_bit_levels = [&]() {
320 bit_levels.back().set(first_bit + bit_levels.size() - 1);
321 return bit_levels.back();
327 tol_geometry, tol_cut_close, fixed_edges, corner_nodes,
328 update_meshsets,
debug);
339 MOFEM_LOG_C(
"WORLD", Sev::inform,
"Min quality cut %6.4g",
340 get_min_quality(cut_bit,
th));
343 Range starting_volume_nodes;
345 starting_volume_nodes,
true);
346 std::vector<double> staring_volume_coords(3 * starting_volume_nodes.size());
348 &*staring_volume_coords.begin());
351 std::vector<double> staring_volume_th_coords(3 *
352 starting_volume_nodes.size());
354 &*staring_volume_th_coords.begin());
356 &*staring_volume_th_coords.begin());
365 update_meshsets,
debug);
367 MOFEM_LOG_C(
"WORLD", Sev::inform,
"Min quality trim %3.2g",
368 get_min_quality(trim_bit,
th));
370 first_bit += bit_levels.size() - 1;
374 &*staring_volume_coords.begin());
380 int &first_bit,
const int fraction_level, Tag
th,
const double tol_geometry,
381 const double tol_cut_close,
const double tol_trim_close,
Range &fixed_edges,
382 Range &corner_nodes,
const bool update_meshsets,
const bool debug) {
386 std::vector<BitRefLevel> bit_levels;
388 auto get_back_bit_levels = [&]() {
390 bit_levels.back().set(first_bit + bit_levels.size() - 1);
391 return bit_levels.back();
396 "ents_not_in_database.vtk",
"VTK",
"");
400 &fixed_edges, &corner_nodes, update_meshsets,
debug);
403 "cut_trim_ents_not_in_database.vtk",
"VTK",
"");
411 update_meshsets,
debug);
419 if (min_q < 0 &&
debug) {
421 "negative_tets.vtk",
"VTK",
"", tets_level,
th);
426 MOFEM_LOG_C(
"WORLD", Sev::inform,
"Min quality node merge %6.4g",
427 get_min_quality(bit_level3,
th));
436 first_bit += bit_levels.size() - 1;
440 bit_level3,
BitRefLevel().set(), MBTET,
"out_tets_merged.vtk",
"VTK",
443 "cut_trim_merge_ents_not_in_database.vtk",
"VTK",
"");
452 moab::Interface &moab = m_field.
get_moab();
457 Range tets_skin_edges;
458 ErrorCode tmp_result;
459 tmp_result = moab.get_adjacencies(tets_skin, 1,
false, tets_skin_edges,
460 moab::Interface::UNION);
462 if (MB_SUCCESS != tmp_result)
464 "Duplicated edges: most likely the source of error is comming from "
465 "adding the vertices of the cracking "
466 "volume to a BLOCKSET rather than NODESET (corresponding to the "
467 "input parameter-vertex_block_set)");
471 fRont = subtract(surface_skin, tets_skin_edges);
480 moab::Interface &moab = m_field.
get_moab();
484 auto create_tag = [&](
const std::string name,
const int dim) {
486 rval = moab.tag_get_handle(name.c_str(),
th);
487 if (
rval == MB_SUCCESS)
489 std::vector<double> def_val(
dim, 0);
490 CHKERR moab.tag_get_handle(name.c_str(),
dim, MB_TYPE_DOUBLE,
th,
491 MB_TAG_CREAT | MB_TAG_SPARSE, &*def_val.begin());
496 auto set_vol = [&](
const Range &vol_verts, std::vector<double> &coords,
497 std::vector<double> &dist_surface_vec,
498 std::vector<double> &dist_surface_normal_vec,
499 std::vector<double> &dist_surface_signed_dist_vec) {
502 coords.resize(3 * vol_verts.size());
503 dist_surface_vec.resize(3 * vol_verts.size());
504 dist_surface_normal_vec.resize(3 * vol_verts.size());
505 dist_surface_signed_dist_vec.resize(vol_verts.size());
507 CHKERR moab.get_coords(vol_verts, &*coords.begin());
510 for (
auto v : vol_verts) {
512 const int index = vol_verts.index(
v);
517 &point_out[0], facets_out);
521 CHKERR tools_interface->getTriNormal(facets_out, &*n_first.begin());
525 auto check_triangle_orientation = [&](
auto n) {
528 CHKERR moab.get_connectivity(facets_out, conn, num_nodes,
true);
530 CHKERR moab.get_coords(conn, 3, &coords(0, 0));
533 for (
auto ii : {0, 1, 2})
534 for (
auto jj : {0, 1, 2})
535 center(jj) += coords(ii, jj) / 3;
537 std::vector<EntityHandle> triangles_out;
538 std::vector<double> distance_out;
539 auto a_max = norm_2(
n);
542 const double ray_length = 2 * sqrt(a_max);
546 std::numeric_limits<float>::epsilon(), &pt[0], &ray[0],
549 if (triangles_out.size() > 1) {
551 for (
auto t : triangles_out) {
552 CHKERR tools_interface->getTriNormal(
t, &*nt.begin());
553 double at = norm_2(nt);
564 auto n = check_triangle_orientation(n_first);
566 const double dot = inner_prod(
delta,
n);
569 if (std::abs(dot) < std::numeric_limits<double>::epsilon()) {
570 if (std::rand() % 2 == 0)
571 delta +=
n * std::numeric_limits<double>::epsilon();
573 delta -=
n * std::numeric_limits<double>::epsilon();
577 noalias(dist_vec) =
delta;
579 auto dist_normal_vec =
582 dist_surface_signed_dist_vec[index] = dot;
583 noalias(dist_normal_vec) = dot *
n;
589 std::vector<double> coords;
590 std::vector<double> dist_surface_vec;
591 std::vector<double> dist_surface_normal_vec;
592 std::vector<double> dist_surface_signed_dist_vec;
596 CHKERR set_vol(vol_verts, coords, dist_surface_vec, dist_surface_normal_vec,
597 dist_surface_signed_dist_vec);
599 CHKERR moab.tag_set_data(create_tag(
"DIST_SURFACE_VECTOR", 3), vol_verts,
600 &*dist_surface_vec.begin());
601 CHKERR moab.tag_set_data(create_tag(
"DIST_SURFACE_NORMAL_VECTOR", 3),
602 vol_verts, &*dist_surface_normal_vec.begin());
603 CHKERR moab.tag_set_data(create_tag(
"DIST_SURFACE_NORMAL_SIGNED", 1),
604 vol_verts, &*dist_surface_signed_dist_vec.begin());
613 moab::Interface &moab = m_field.
get_moab();
616 auto create_tag = [&](
const std::string name,
const int dim) {
618 rval = moab.tag_get_handle(name.c_str(),
th);
619 if (
rval == MB_SUCCESS)
621 std::vector<double> def_val(
dim, 0);
622 CHKERR moab.tag_get_handle(name.c_str(),
dim, MB_TYPE_DOUBLE,
th,
623 MB_TAG_CREAT | MB_TAG_SPARSE, &*def_val.begin());
628 CHKERR moab.get_connectivity(vol, vol_vertices,
true);
630 std::vector<double> min_distances_from_front(vol_vertices.size(), -1);
631 std::vector<double> points_on_edges(3 * vol_vertices.size(), 0);
632 std::vector<EntityHandle> closest_edges(vol_vertices.size(), 0);
634 std::vector<double> coords(3 * vol_vertices.size());
636 CHKERR moab.tag_get_data(
th, vol_vertices, &*coords.begin());
638 CHKERR moab.get_coords(vol_vertices, &*coords.begin());
641 &*coords.begin(), vol_vertices.size(),
fRont,
642 &*min_distances_from_front.begin(), &*points_on_edges.begin(),
643 &*closest_edges.begin());
645 if (!points_on_edges.empty()) {
646 for (
int i = 0;
i != min_distances_from_front.size(); ++
i) {
648 CHKERR moab.get_adjacencies(&closest_edges[0], 1, 2,
false, faces);
651 point_out -= point_in;
655 auto th_dist_front_vec = create_tag(
"DIST_FRONT_VECTOR", 3);
656 CHKERR moab.tag_set_data(th_dist_front_vec, vol_vertices,
657 &*points_on_edges.begin());
663 Tag
th,
Range &vol_edges,
const bool remove_adj_prims_edges,
int verb,
664 const bool debug,
const std::string edges_file_name) {
666 moab::Interface &moab = m_field.
get_moab();
669 auto get_tag_dim = [&](
auto th) {
674 auto dim = get_tag_dim(
th);
676 auto get_tag_data = [&](
auto th,
auto conn) {
678 CHKERR moab.tag_get_by_ptr(
th, &conn, 1, &ptr);
680 const_cast<double *
>(
static_cast<const double *
>(ptr)), 3);
683 auto get_edge_ray = [&](
auto conn) {
685 CHKERR moab.get_coords(conn, 2, &coords[0]);
693 CHKERR moab.get_adjacencies(
vOlume, 1,
true, edges, moab::Interface::UNION);
695 auto remove_prisms_edges = [&](
Range &edges) {
698 CHKERR moab.get_adjacencies(edges, 3,
false, prisms,
699 moab::Interface::UNION);
700 prisms = prisms.subset_by_type(MBPRISM);
702 CHKERR moab.get_connectivity(prisms, prisms_verts,
true);
704 CHKERR moab.get_adjacencies(prisms_verts, 1,
false, prism_edges,
705 moab::Interface::UNION);
706 edges = subtract(edges, prism_edges);
709 if (remove_adj_prims_edges)
710 CHKERR remove_prisms_edges(edges);
712 std::vector<EntityHandle> cut_edges;
713 cut_edges.reserve(edges.size());
715 auto get_cut_edges_vec = [&](
auto th,
auto &cut_edges,
auto e,
auto &&conn) {
718 auto ray = get_edge_ray(conn);
719 const double length = norm_2(ray);
721 const auto dist0 = get_tag_data(
th, conn[0]);
722 const auto dist1 = get_tag_data(
th, conn[1]);
723 const double max_dist = std::max(norm_2(dist0), norm_2(dist1));
724 if (max_dist < length) {
725 cut_edges.push_back(e);
731 auto get_cut_edges_signed_dist = [&](
auto th,
auto &cut_edges,
auto e,
734 auto get_tag_signed_dist = [&](
auto conn) {
736 CHKERR moab.tag_get_data(
th, &conn, 1, &dist);
739 const auto dist0 = get_tag_signed_dist(conn[0]);
740 const auto dist1 = get_tag_signed_dist(conn[1]);
741 const double opposite = dist0 * dist1;
743 cut_edges.push_back(e);
747 auto get_conn = [&](
auto e) {
750 CHKERR moab.get_connectivity(e, conn, num_nodes,
true);
756 CHKERR get_cut_edges_vec(
th, cut_edges, e, get_conn(e));
759 CHKERR get_cut_edges_signed_dist(
th, cut_edges, e, get_conn(e));
761 CHKERR moab.get_adjacencies(&*cut_edges.begin(), cut_edges.size(), 3,
false,
762 vol_edges, moab::Interface::UNION);
763 vol_edges = intersect(vol_edges,
vOlume);
765 auto convert_to_ranger = [](
auto &
v) {
767 r.insert_list(
v.begin(),
v.end());
771 if (
debug && !edges_file_name.empty())
773 convert_to_ranger(cut_edges));
781 moab::Interface &moab = m_field.
get_moab();
785 Tag th_dist_front_vec;
786 CHKERR moab.tag_get_handle(
"DIST_FRONT_VECTOR", th_dist_front_vec);
788 debug,
"cutFrontEdges.vtk");
792 Tag th_dist_surface_vec;
793 CHKERR moab.tag_get_handle(
"DIST_SURFACE_VECTOR", th_dist_surface_vec);
796 debug,
"cutSurfaceEdges.vtk");
810 const int surf_levels,
811 const int front_levels,
812 Range *fixed_edges,
int verb,
815 moab::Interface &moab = m_field.
get_moab();
822 auto add_bit = [&](
const int bit) {
827 for (
auto d : {2, 1, 0}) {
829 moab::Interface::UNION);
835 CHKERR add_bit(init_bit_level);
837 auto update_range = [&](
Range *r_ptr) {
842 r_ptr->merge(childs);
850 CHKERR moab.get_connectivity(tets, verts,
true);
852 CHKERR moab.get_adjacencies(verts, 1,
true, ref_edges,
853 moab::Interface::UNION);
858 CHKERR update_range(fixed_edges);
861 ->updateAllMeshsetsByEntitiesChildren(
bit);
872 *fixed_edges = intersect(*fixed_edges, bit_edges);
877 for (
int ll = init_bit_level; ll != init_bit_level + surf_levels; ++ll) {
883 for (
int ll = init_bit_level + surf_levels;
884 ll != init_bit_level + surf_levels + front_levels; ++ll) {
900 moab::Interface &moab = m_field.
get_moab();
907 CHKERR moab.tag_get_handle(
"DIST_SURFACE_NORMAL_SIGNED", th_signed_dist);
909 auto get_tag_edge_dist = [&](
auto th,
auto conn) {
910 std::array<double, 2> r;
911 CHKERR moab.tag_get_data(
th, conn, 2, r.data());
915 auto get_conn = [&](
const auto e) {
918 CHKERR moab.get_connectivity(e, conn, num_nodes,
true);
922 auto get_adj = [&moab](
const Range r,
const int dim) {
925 CHKERR moab.get_adjacencies(r,
dim,
false,
a, moab::Interface::UNION);
927 CHKERR moab.get_connectivity(r,
a,
true);
931 auto vol_edges = get_adj(vol, 1);
935 int nb_ave_length = 0;
936 for (
auto e : vol_edges) {
938 auto conn = get_conn(e);
939 auto dist = get_tag_edge_dist(th_signed_dist, conn);
940 const auto dist_max = std::max(dist[0], dist[1]);
941 const auto dist_min = std::min(dist[0], dist[1]);
942 const auto dot = dist[0] * dist[1];
945 CHKERR moab.get_coords(conn, 2, &coords[0]);
949 const double ray_length = norm_2(ray);
954 (dot < 0 && dist_max > std::numeric_limits<double>::epsilon()) ||
955 (std::abs(dist_min) < std::numeric_limits<double>::epsilon() &&
956 dist_max > std::numeric_limits<double>::epsilon())
961 const double s = dist[0] / (dist[0] - dist[1]);
962 const double dist = s * ray_length;
964 auto add_edge = [&](
auto dist) {
992 Range *fixed_edges,
Range *corner_nodes,
const double geometry_tol,
993 const double close_tol,
const int verb,
const bool debug) {
995 moab::Interface &moab = m_field.
get_moab();
1006 auto get_adj = [&moab](
const Range r,
const int dim) {
1009 CHKERR moab.get_adjacencies(r,
dim,
false,
a, moab::Interface::UNION);
1011 CHKERR moab.get_connectivity(r,
a,
true);
1015 auto get_skin = [&skin](
const auto r) {
1017 CHKERR skin.find_skin(0, r,
false, s);
1022 auto get_range = [](std::vector<EntityHandle>
v) {
1024 r.insert_list(
v.begin(),
v.end());
1028 auto get_coords = [&](
const auto v) {
1030 CHKERR moab.get_coords(&
v, 1, &coords[0]);
1035 CHKERR moab.tag_get_handle(
"DIST_SURFACE_NORMAL_VECTOR", th_normal_dist);
1039 CHKERR moab.tag_get_data(th_normal_dist, &
v, 1, &*dist_vec.begin());
1044 const double geometry_tol) {
1045 auto get_tuple = [](
const EntityHandle e,
const double dist,
1046 const double l) {
return std::make_tuple(e, dist,
l); };
1048 std::tuple<EntityHandle, double, double> min_tuple{0, 0, 0};
1050 for (
auto e : edges) {
1054 auto get_dist = [&](
auto eit) {
1057 CHKERR moab.get_connectivity(eit->first, conn, num_nodes,
true);
1059 return eit->second.dIst;
1060 else if (conn[1] ==
v)
1061 return (eit->second.lEngth - eit->second.dIst);
1066 const auto d = get_dist(eit);
1067 if (std::get<0>(min_tuple) == 0) {
1068 min_tuple = get_tuple(e, d, eit->second.lEngth);
1069 }
else if (std::get<1>(min_tuple) > d) {
1070 min_tuple = get_tuple(e, d, eit->second.lEngth);
1075 const auto geom_dist_vec = get_normal_dits(
v);
1076 const auto geom_tol = norm_2(geom_dist_vec);
1077 const auto l = std::get<2>(min_tuple);
1079 if (geom_tol <
l * geometry_tol) {
1080 return std::make_pair(get_coords(
v),
l);
1083 const auto &d =
edgesToCut.at(std::get<0>(min_tuple));
1084 return std::make_pair(
VectorDouble3(d.rayPoint + d.dIst * d.unitRayDir),
1089 auto get_in_range = [](
auto v,
auto &r) {
return (r.find(
v) != r.end()); };
1091 auto project_nodes = [&](
auto nodes_to_check) {
1095 auto get_fix_e = [](
auto fixed_edges) {
1097 return *fixed_edges;
1102 const auto fix_e = get_fix_e(fixed_edges);
1104 const auto cut_fix = intersect(
cutEdges, fix_e);
1105 const auto cut_skin = intersect(
cutEdges, skin_e);
1109 corner_n = intersect(*corner_nodes, nodes_to_check);
1110 auto fixe_n = intersect(get_adj(fix_e, 0), nodes_to_check);
1111 auto skin_n = intersect(get_adj(skin_e, 0), nodes_to_check);
1113 std::vector<std::pair<EntityHandle, TreeData>> vertices_on_cut_edges;
1114 vertices_on_cut_edges.reserve(nodes_to_check.size());
1116 auto add_zero_vertex = [&](
const EntityHandle v,
auto &&new_pos) {
1117 auto coords = get_coords(
v);
1118 auto ray = new_pos - coords;
1119 auto l = norm_2(ray);
1122 if (
l > std::numeric_limits<double>::epsilon()) {
1129 return std::make_pair(
v, data);
1132 auto intersect_v = [&](
const auto v,
const auto r) {
1133 return intersect(r, get_ent_adj(
v, 1));
1136 for (
auto v : nodes_to_check) {
1138 const auto e = intersect_v(
v,
cutEdges);
1141 if (get_in_range(
v, corner_n)) {
1142 auto p = get_prj_point(
v, e, 0);
1143 if (norm_2(get_coords(
v) -
p.first) < close_tol *
p.second) {
1144 vertices_on_cut_edges.push_back(add_zero_vertex(
v, get_coords(
v)));
1148 }
else if (get_in_range(
v, fixe_n)) {
1149 const auto b = intersect_v(
v, cut_fix);
1151 auto p = get_prj_point(
v, b, geometry_tol);
1152 if (norm_2(get_coords(
v) -
p.first) < close_tol *
p.second) {
1153 vertices_on_cut_edges.push_back(add_zero_vertex(
v,
p.first));
1157 }
else if (get_in_range(
v, skin_n)) {
1158 const auto b = intersect_v(
v, cut_skin);
1160 auto p = get_prj_point(
v, b, geometry_tol);
1161 if (norm_2(get_coords(
v) -
p.first) < close_tol *
p.second) {
1162 vertices_on_cut_edges.push_back(add_zero_vertex(
v,
p.first));
1170 auto p = get_prj_point(
v, e, geometry_tol);
1171 if (norm_2(get_coords(
v) -
p.first) < close_tol *
p.second) {
1172 if (get_in_range(
v, fixe_n) || get_in_range(
v, skin_n))
1173 vertices_on_cut_edges.push_back(add_zero_vertex(
v, get_coords(
v)));
1175 vertices_on_cut_edges.push_back(add_zero_vertex(
v,
p.first));
1181 auto get_distances = [&](
auto &data) {
1182 std::map<EntityHandle, double> dist_map;
1183 if (!data.empty()) {
1185 CHKERR moab.tag_get_handle(
"DIST_SURFACE_NORMAL_SIGNED",
th);
1187 std::vector<EntityHandle> verts;
1188 verts.reserve(data.size());
1190 verts.emplace_back(d.first);
1191 std::vector<double> distances(verts.size());
1192 CHKERR moab.tag_get_data(
th, &*verts.begin(), verts.size(),
1193 &*distances.begin());
1194 for (
size_t i = 0;
i != verts.size(); ++
i)
1195 dist_map[verts[
i]] = distances[
i];
1200 auto dist_map = get_distances(vertices_on_cut_edges);
1201 if (!dist_map.empty()) {
1202 auto cmp = [&dist_map](
const auto &
a,
const auto &b) {
1203 return dist_map[
a.first] < dist_map[b.first];
1205 std::sort(vertices_on_cut_edges.begin(), vertices_on_cut_edges.end(),
1209 return vertices_on_cut_edges;
1212 auto get_min_quality =
1213 [&](
const Range &adj_tets,
1214 const map<EntityHandle, TreeData> &vertices_on_cut_edges,
1215 const std::pair<EntityHandle, TreeData> *test_data_ptr) {
1217 for (
auto t : adj_tets) {
1222 CHKERR moab.get_coords(conn, num_nodes, &coords[0]);
1224 auto set_new_coord = [](
auto d) {
1225 return d.rayPoint + d.dIst * d.unitRayDir;
1228 for (
auto n : {0, 1, 2, 3}) {
1230 if (test_data_ptr) {
1231 if (test_data_ptr->first == conn[
n]) {
1232 noalias(n_coords) = set_new_coord(test_data_ptr->second);
1235 auto mit = vertices_on_cut_edges.find(conn[
n]);
1236 if (mit != vertices_on_cut_edges.end()) {
1237 noalias(n_coords) = set_new_coord(mit->second);
1245 auto get_zero_distance_verts = [&](
const auto &&vertices_on_cut_edges) {
1246 std::vector<EntityHandle> zero_dist_vec;
1247 zero_dist_vec.reserve(vertices_on_cut_edges.size());
1248 for (
auto t : vertices_on_cut_edges) {
1250 auto adj_tet = intersect(
vOlume, get_ent_adj(
t.first, 3));
1255 zero_dist_vec.push_back(
t.first);
1258 return zero_dist_vec;
1263 auto get_zero_distant_ents = [&](
auto bridge_ents,
const int dim,
1264 const int bridge_dim) {
1266 intersect(get_adj(bridge_ents,
dim), get_adj(vol_cut_ents,
dim));
1267 auto ents_to_remove =
1268 get_adj(subtract(get_adj(ents, bridge_dim), bridge_ents),
dim);
1269 return subtract(ents, ents_to_remove);
1275 get_range(get_zero_distance_verts(project_nodes(get_adj(
cutEdges, 0)))));
1285 auto adj_edges = get_ent_adj(f, 1);
1286 for (
auto e : adj_edges) {
1304 moab::Interface &moab = m_field.
get_moab();
1312 auto refine_mesh = [&]() {
1331 bit,
bit, MBTRI, cut_surf);
1342 boost::make_tuple(MBVERTEX,
m.first));
1346 "No vertex on cut edges, that make no sense");
1348 const boost::shared_ptr<RefEntity> &ref_ent = *vit;
1349 if ((ref_ent->getBitRefLevel() &
bit).any()) {
1351 cut_verts.insert(vert);
1356 "Vertex has wrong bit ref level %s",
1357 boost::lexical_cast<std::string>(ref_ent->getBitRefLevel()).c_str());
1363 Skinner skin(&moab);
1364 CHKERR skin.find_skin(0, cut_vols,
false, tets_skin);
1372 CHKERR moab.get_connectivity(cut_surf, diff_verts,
true);
1373 diff_verts = subtract(diff_verts, cut_verts);
1375 Range subtract_faces;
1376 CHKERR moab.get_adjacencies(diff_verts, 2,
false, subtract_faces,
1377 moab::Interface::UNION);
1378 cut_surf = subtract(cut_surf, unite(subtract_faces, tets_skin));
1381 CHKERR moab.get_connectivity(cut_surf, cut_verts,
true);
1384 auto check_for_non_minfold = [&]() {
1387 CHKERR moab.get_adjacencies(cut_surf, 1,
false, surf_edges,
1388 moab::Interface::UNION);
1389 for (
auto e : surf_edges) {
1392 CHKERR moab.get_adjacencies(&e, 1, 2,
false, faces);
1393 faces = intersect(faces, cut_surf);
1394 if (faces.size() > 2) {
1396 bool resolved =
false;
1400 CHKERR moab.get_connectivity(faces, nodes,
true);
1401 for (
auto n : nodes) {
1403 CHKERR moab.get_adjacencies(&
n, 1, 2,
false, adj_faces);
1404 adj_faces = intersect(adj_faces, cut_surf);
1405 if (adj_faces.size() == 1) {
1406 cut_surf.erase(adj_faces[0]);
1413 CHKERR moab.get_adjacencies(faces, 1,
false, adj_edges,
1414 moab::Interface::UNION);
1415 adj_edges = intersect(adj_edges, surf_edges);
1417 for (
auto other_e : adj_edges) {
1419 CHKERR moab.get_adjacencies(&other_e, 1, 2,
false, other_faces);
1420 other_faces = intersect(other_faces, cut_surf);
1421 if (other_faces.size() > 2) {
1422 other_faces = intersect(other_faces, faces);
1423 cut_surf = subtract(cut_surf, other_faces);
1428 if (!resolved && !
debug)
1430 "Non-mainfold surface");
1433 CHKERR moab.get_connectivity(cut_surf, cut_verts,
true);
1439 CHKERR check_for_non_minfold();
1449 moab::Interface &moab = m_field.
get_moab();
1454 double dist =
m.second.dIst;
1455 VectorDouble3 new_coors =
m.second.rayPoint + dist *
m.second.unitRayDir;
1457 CHKERR moab.tag_set_data(
th, &
m.first, 1, &new_coors[0]);
1459 CHKERR moab.set_coords(&
m.first, 1, &new_coors[0]);
1467 moab::Interface &moab = m_field.
get_moab();
1470 double dist =
v.second.dIst;
1471 VectorDouble3 new_coors =
v.second.rayPoint + dist *
v.second.unitRayDir;
1473 CHKERR moab.tag_set_data(
th, &
v.first, 1, &new_coors[0]);
1475 CHKERR moab.set_coords(&
v.first, 1, &new_coors[0]);
1482 const double tol_trim_close,
1485 moab::Interface &moab = m_field.
get_moab();
1489 Skinner skin(&moab);
1495 Range tets_skin_verts;
1496 CHKERR moab.get_connectivity(tets_skin, tets_skin_verts,
true);
1498 Range tets_skin_edges;
1499 CHKERR moab.get_adjacencies(tets_skin, 1,
false, tets_skin_edges,
1500 moab::Interface::UNION);
1503 Range cut_surface_edges;
1505 moab::Interface::UNION);
1506 Range cut_surface_verts;
1511 corners = *corner_nodes;
1515 fix_edges = *fixed_edges;
1517 Range fixed_edges_verts;
1519 CHKERR moab.get_connectivity(*fixed_edges, fixed_edges_verts,
true);
1525 surface_skin =
fRont;
1530 CHKERR moab.tag_get_data(
th, &
v, 1, &point[0]);
1532 CHKERR moab.get_coords(&
v, 1, &point[0]);
1541 : d(d),
v(
v), e(e) {}
1544 typedef multi_index_container<
1547 ordered_non_unique<member<VertMap, const double, &VertMap::d>>,
1548 ordered_non_unique<member<VertMap, const EntityHandle, &VertMap::v>>,
1549 ordered_non_unique<member<VertMap, const EntityHandle, &VertMap::e>>
1554 VertMapMultiIndex verts_map;
1557 const double dist) {
1558 verts_map.insert(VertMap(dist,
v, e));
1567 auto cut_this_edge = [&](
const EntityHandle e,
const double length,
auto &ray,
1574 for (
int ii = 0; ii != 3; ++ii)
1576 for (
int ii = 0; ii != 3; ++ii)
1586 auto get_tag_data = [&](
auto th,
auto conn) {
1588 CHKERR moab.tag_get_data(
th, &conn, 1, &
t(0));
1592 double max_edge_length = 0;
1595 if (!
fRont.empty()) {
1597 treeSurfPtr = boost::shared_ptr<OrientedBoxTreeTool>(
1598 new OrientedBoxTreeTool(&moab,
"ROOTSETSURF",
true));
1601 for (
auto s : surface_skin) {
1603 auto project_on_surface = [&](
auto &
t) {
1613 t_normal(
i) /= sqrt(t_normal(
i) * t_normal(
i));
1614 const double dot = t_normal(
i) * (t_p(
i) -
t(
i));
1615 t_p(
i) =
t(
i) + dot * t_normal(
i);
1621 CHKERR moab.get_connectivity(s, conn, num_nodes,
true);
1624 CHKERR moab.get_coords(conn, num_nodes, &coords_front[0]);
1631 auto t_p_f0 = project_on_surface(t_f0);
1632 auto t_p_f1 = project_on_surface(t_f1);
1634 CHKERR moab.set_coords(&conn[0], 1, &t_p_f0(0));
1635 CHKERR moab.set_coords(&conn[1], 1, &t_p_f1(0));
1643 Tag th_dist_front_vec;
1644 CHKERR moab.tag_get_handle(
"DIST_FRONT_VECTOR", th_dist_front_vec);
1647 for (
auto e : cut_surface_edges) {
1649 auto get_conn_edge = [&moab](
auto e) {
1653 CHKERR moab.get_connectivity(e, conn_edge, num_nodes,
true);
1657 auto get_coords_edge = [&moab](
auto conn_edge) {
1658 std::array<double, 6> coords_edge;
1659 CHKERR moab.get_coords(conn_edge, 2, coords_edge.data());
1663 auto get_ftensor_coords = [](
const double *ptr) {
1667 auto conn_edge = get_conn_edge(e);
1668 auto coords_edge = get_coords_edge(conn_edge);
1670 auto t_dist0 = get_tag_data(th_dist_front_vec, conn_edge[0]);
1671 auto t_dist1 = get_tag_data(th_dist_front_vec, conn_edge[1]);
1672 auto t_e0 = get_ftensor_coords(&coords_edge[0]);
1673 auto t_e1 = get_ftensor_coords(&coords_edge[3]);
1676 t_edge_delta(
i) = t_e1(
i) - t_e0(
i);
1677 const double edge_length2 = t_edge_delta(
i) * t_edge_delta(
i);
1678 const double edge_length = sqrt(edge_length2);
1679 if (edge_length == 0)
1682 auto add_edge = [&](
auto t_cut) {
1685 t_ray(
i) = t_cut * t_edge_delta(
i);
1686 add_vert(conn_edge[0], e, fabs(t_cut));
1687 add_vert(conn_edge[1], e, fabs(t_cut - 1));
1688 cut_this_edge(e, edge_length, t_ray, t_e0);
1691 t_edge_point(
i) = t_e0(
i) + t_cut * t_edge_delta(
i);
1692 t_ray(
i) = t_edge_point(
i) - t_e1(
i);
1693 add_vert(conn_edge[0], e, fabs(t_cut));
1694 add_vert(conn_edge[1], e, fabs(t_cut - 1));
1695 cut_this_edge(e, edge_length, t_ray, t_e1);
1697 max_edge_length = std::max(max_edge_length, edge_length);
1700 auto trim_cross_edges = [&]() {
1701 if (std::min(t_dist0(
i) * t_dist0(
i), t_dist1(
i) * t_dist1(
i)) <
1704 auto make_pair = [&](
const double d,
const double te) {
1705 return std::make_pair(d, te);
1708 auto min_pair = make_pair(-1, -1);
1710 for (
auto f : surface_skin) {
1712 auto conn_front = get_conn_edge(f);
1713 auto coords_front = get_coords_edge(conn_front);
1714 auto t_f0 = get_ftensor_coords(&coords_front[0]);
1715 auto t_f1 = get_ftensor_coords(&coords_front[3]);
1721 &t_e0(0), &t_e1(0), &t_f0(0), &t_f1(0), &te, &tf);
1725 auto check = [](
auto v) {
1726 return (
v > -std::numeric_limits<double>::epsilon() &&
1727 (
v - 1) < std::numeric_limits<double>::epsilon());
1730 if (check(te) && check(tf)) {
1733 t_delta(
i) = (t_e0(
i) + te * t_edge_delta(
i)) -
1734 (t_f0(0) + tf * (t_f1(
i) - t_f0(
i)));
1736 t_delta(
i) /= edge_length;
1738 const double dot = t_delta(
i) * t_delta(
i);
1740 if (min_pair.first < 0 || min_pair.first > dot) {
1742 if (dot < tol_trim_close * tol_trim_close)
1743 min_pair = make_pair(dot, te);
1749 if (min_pair.first > -std::numeric_limits<double>::epsilon()) {
1750 add_edge(min_pair.second);
1758 if (!trim_cross_edges()) {
1760 const double dot = t_dist0(
i) * t_dist1(
i);
1761 const double dot_direction0 = t_dist0(
i) * t_edge_delta(
i);
1762 const double dot_direction1 = t_dist1(
i) * t_edge_delta(
i);
1764 if (dot < std::numeric_limits<double>::epsilon() &&
1765 dot_direction0 > -std::numeric_limits<double>::epsilon() &&
1766 dot_direction1 < std::numeric_limits<double>::epsilon()) {
1768 const double s0 = t_dist0(
i) * t_edge_delta(
i) / edge_length;
1769 const double s1 = t_dist1(
i) * t_edge_delta(
i) / edge_length;
1770 const double t_cut = s0 / (s0 - s1);
1778 CHKERR SaveData(moab)(
"edges_potentially_to_trim.vtk", cut_surface_edges);
1789 for (
auto t : adj_tets) {
1795 CHKERR moab.tag_get_data(
th, conn, num_nodes, &coords[0]);
1797 CHKERR moab.get_coords(conn, num_nodes, &coords[0]);
1799 for (
int n = 0;
n != 4; ++
n) {
1802 noalias(n_coords) = new_pos;
1807 m->second.rayPoint +
m->second.dIst *
m->second.unitRayDir;
1812 if (!std::isnormal(q)) {
1816 min_q = std::min(min_q, q);
1821 Range checked_verts;
1829 CHKERR moab.get_adjacencies(&
v, 1, 3,
false, adj_tets);
1831 double q = get_quality_change(adj_tets,
v, new_pos);
1834 double dist = norm_2(ray_dir);
1841 checked_verts.insert(
v);
1853 checked_verts.insert(
v);
1861 auto hi = verts_map.get<0>().upper_bound(tol_trim_close);
1862 verts_map.get<0>().erase(hi, verts_map.end());
1864 auto remove_verts = [&](
Range nodes) {
1865 for (
auto n : nodes) {
1866 auto r = verts_map.get<1>().equal_range(
n);
1867 verts_map.get<1>().erase(r.first, r.second);
1871 auto trim_verts = [&](
const Range verts,
const bool move) {
1872 VertMapMultiIndex verts_map_tmp;
1873 for (
auto p = corners.pair_begin();
p != corners.pair_end(); ++
p) {
1874 auto lo = verts_map.get<1>().lower_bound(
p->first);
1875 auto hi = verts_map.get<1>().upper_bound(
p->second);
1876 verts_map_tmp.insert(lo, hi);
1879 for (
auto &
m : verts_map_tmp.get<0>())
1880 add_trim_vert(
m.v,
m.e);
1882 for (
auto &
m : verts_map_tmp.get<0>())
1883 add_no_move_trim(
m.v,
m.e);
1887 auto trim_edges = [&](
const Range ents,
const bool move) {
1888 VertMapMultiIndex verts_map_tmp;
1889 for (
auto p = ents.pair_begin();
p != ents.pair_end(); ++
p) {
1890 auto lo = verts_map.get<2>().lower_bound(
p->first);
1891 auto hi = verts_map.get<2>().upper_bound(
p->second);
1892 verts_map_tmp.insert(lo, hi);
1893 verts_map.get<2>().erase(lo, hi);
1896 for (
auto &
m : verts_map_tmp.get<0>())
1897 add_trim_vert(
m.v,
m.e);
1899 for (
auto &
m : verts_map_tmp.get<0>())
1900 add_no_move_trim(
m.v,
m.e);
1904 auto intersect_self = [&](
Range &
a,
const Range b) {
a = intersect(
a, b); };
1906 map<std::string, Range> range_maps;
1911 intersect_self(range_maps[
"surface_skin"],
trimEdges);
1913 range_maps[
"fixed_edges_on_surface_skin"] =
1914 intersect(range_maps[
"surface_skin"], fix_edges);
1917 CHKERR moab.get_adjacencies(fixed_edges_verts, 1,
false,
1918 range_maps[
"fixed_edges_verts_edges"],
1919 moab::Interface::UNION);
1920 intersect_self(range_maps[
"fixed_edges_verts_edges"],
trimEdges);
1922 CHKERR moab.get_connectivity(
1923 range_maps[
"fixed_edges_verts_edges"],
1924 range_maps[
"fixed_edges_verts_edges_verts_on_the_skin"],
false);
1925 intersect_self(range_maps[
"fixed_edges_verts_edges_verts_on_the_skin"],
1929 trim_verts(corners,
false);
1930 remove_verts(corners);
1933 trim_edges(range_maps[
"fixed_edges_on_surface_skin"],
true);
1934 remove_verts(range_maps[
"fixed_edges_on_surface_skin_verts"]);
1937 trim_verts(range_maps[
"fixed_edges_verts_edges_verts_on_the_skin"],
false);
1938 remove_verts(range_maps[
"fixed_edges_verts_edges_verts_on_the_skin"]);
1941 trim_edges(range_maps[
"surface_skin"],
true);
1942 trim_verts(tets_skin_verts,
false);
1943 remove_verts(tets_skin_verts);
1945 for (
auto &
m : verts_map.get<0>())
1946 add_trim_vert(
m.v,
m.e);
1951 CHKERR moab.get_adjacencies(&
v, 1, 1,
false, adj);
1952 for (
auto e : adj) {
1972 moab::Interface &moab = m_field.
get_moab();
1985 for (map<EntityHandle, TreeData>::iterator mit =
edgesToTrim.begin();
1988 boost::make_tuple(MBVERTEX, mit->first));
1992 "No vertex on trim edges, that make no sense");
1994 const boost::shared_ptr<RefEntity> &ref_ent = *vit;
1995 if ((ref_ent->getBitRefLevel() &
bit).any()) {
2007 Range trim_new_surfaces_nodes;
2009 trim_new_surfaces_nodes = subtract(trim_new_surfaces_nodes,
trimNewVertices);
2010 trim_new_surfaces_nodes = subtract(trim_new_surfaces_nodes,
cutNewVertices);
2011 Range faces_not_on_surface;
2012 CHKERR moab.get_adjacencies(trim_new_surfaces_nodes, 2,
false,
2013 faces_not_on_surface, moab::Interface::UNION);
2017 Range all_surfaces_on_bit_level;
2020 all_surfaces_on_bit_level =
2029 Range *corner_nodes,
2030 const double tol, Tag
th,
2033 moab::Interface &moab = m_field.
get_moab();
2034 Skinner skin(&moab);
2036 auto get_adj = [&moab](
const Range r,
const int dim) {
2039 CHKERR moab.get_adjacencies(r,
dim,
false,
a, moab::Interface::UNION);
2041 CHKERR moab.get_connectivity(r,
a,
true);
2045 auto get_skin = [&skin](
const auto r) {
2047 CHKERR skin.find_skin(0, r,
false, s);
2054 auto trim_tets_skin_edges = get_adj(trim_tets_skin, 1);
2057 auto contarain_edges =
2060 contarain_edges.merge(
2061 intersect(fixed_edges->subset_by_type(MBEDGE), trim_surface_edges));
2066 barrier_vertices.merge(get_adj(
2067 intersect(get_adj(intersect(*corner_nodes, barrier_vertices), 2),
2072 auto get_nodes_with_one_node_on_fixed_edge_other_not = [&]() {
2073 const auto fixed_edges_on_surface =
2074 intersect(*fixed_edges, trim_surface_edges);
2075 const auto skin_fixed_edges_on_surface = get_skin(fixed_edges_on_surface);
2076 const auto barrier_nodes = subtract(skin_fixed_edges_on_surface,
2078 return barrier_nodes;
2082 barrier_vertices.merge(get_nodes_with_one_node_on_fixed_edge_other_not());
2084 auto add_close_surface_barrier = [&]() {
2096 test_edges = subtract(test_edges, *fixed_edges);
2097 auto trim_surface_nodes = get_adj(test_edges, 0);
2098 trim_surface_nodes = subtract(trim_surface_nodes, barrier_vertices);
2100 trim_surface_nodes =
2101 subtract(trim_surface_nodes, get_adj(*fixed_edges, 0));
2104 trim_skin.merge(contarain_edges);
2106 trim_skin.merge(get_skin(
sUrface));
2108 trim_skin.merge(
fRont);
2110 trim_skin.merge(intersect(*fixed_edges, trim_surface_edges));
2112 if (
debug && !trim_skin.empty())
2117 CHKERR moab.tag_get_data(
th, trim_surface_nodes, &*coords.begin());
2119 CHKERR moab.get_coords(trim_surface_nodes, &*coords.begin());
2121 if (!trim_skin.empty()) {
2123 std::vector<double> min_distances(trim_surface_nodes.size(), -1);
2125 &*coords.begin(), trim_surface_nodes.size(), trim_skin,
2126 &*min_distances.begin());
2128 auto coords_ptr = &*coords.begin();
2129 auto min_dist = &*min_distances.begin();
2131 std::vector<EntityHandle> add_nodes;
2132 add_nodes.reserve(trim_surface_nodes.size());
2134 for (
auto n : trim_surface_nodes) {
2136 if ((*min_dist) > std::numeric_limits<double>::epsilon()) {
2141 &point_out[0], facets_out);
2148 normal /= norm_2(normal);
2151 double dist = norm_2(
delta);
2152 if (dist <
tol * (*min_dist))
2153 add_nodes.emplace_back(
n);
2160 barrier_vertices.insert_list(add_nodes.begin(), add_nodes.end());
2166 auto remove_faces_on_skin = [&]() {
2169 if (!skin_faces.empty()) {
2170 barrier_vertices.merge(get_adj(skin_faces, 0));
2171 for (
auto f : skin_faces)
2177 auto get_trim_free_edges = [&]() {
2179 Range trim_surf_skin;
2181 trim_surf_skin = subtract(trim_surf_skin, trim_tets_skin_edges);
2183 Range trim_surf_skin_verts;
2184 CHKERR moab.get_connectivity(trim_surf_skin, trim_surf_skin_verts,
true);
2185 for (
auto e : barrier_vertices)
2186 trim_surf_skin_verts.erase(e);
2189 CHKERR moab.get_adjacencies(trim_surf_skin_verts, 1,
false, free_edges,
2190 moab::Interface::UNION);
2192 subtract(intersect(free_edges, trim_surf_skin), contarain_edges);
2197 CHKERR remove_faces_on_skin();
2198 CHKERR add_close_surface_barrier();
2200 if (
debug && !barrier_vertices.empty())
2206 while (!(out_edges = get_trim_free_edges()).empty()) {
2210 "trimNewSurfaces_" + boost::lexical_cast<std::string>(nn) +
".vtk",
2213 if (
debug && !out_edges.empty())
2215 "trimNewSurfacesOutsideVerts_" +
2216 boost::lexical_cast<std::string>(nn) +
".vtk",
2219 Range outside_faces;
2220 CHKERR moab.get_adjacencies(out_edges, 2,
false, outside_faces,
2221 moab::Interface::UNION);
2228 "trimNewSurfaces_" + boost::lexical_cast<std::string>(nn) +
".vtk",
2238 moab::Interface &moab = m_field.
get_moab();
2246 CHKERR moab.create_meshset(MESHSET_SET, meshset);
2247 CHKERR moab.add_entities(meshset, ents);
2249 meshset, split_bit,
true, front_tris);
2250 CHKERR moab.delete_entities(&meshset, 1);
2251 ents = subtract(ents, front_tris);
2258 Skinner skin(&moab);
2260 CHKERR skin.find_skin(0, tets,
false, tets_skin);
2261 ents = subtract(ents, tets_skin);
2270 moab::Interface &moab = m_field.
get_moab();
2275 CHKERR moab.create_meshset(MESHSET_SET, meshset_volume);
2277 split_bit,
BitRefLevel().set(), MBTET, meshset_volume);
2279 CHKERR moab.create_meshset(MESHSET_SET, meshset_surface);
2280 CHKERR moab.add_entities(meshset_surface, ents);
2284 CHKERR moab.delete_entities(&meshset_volume, 1);
2285 CHKERR moab.delete_entities(&meshset_surface, 1);
2290 for (Range::iterator pit = prisms.begin(); pit != prisms.end(); pit++) {
2293 CHKERR moab.get_connectivity(*pit, conn, num_nodes,
true);
2295 CHKERR moab.tag_get_data(
th, conn, 3, &data(0, 0));
2296 CHKERR moab.tag_set_data(
th, &conn[3], 3, &data(0, 0));
2304 const Range &fixed_edges,
const Range &corner_nodes,
Range &edges_to_merge,
2305 Range &out_tets,
Range &new_surf, Tag
th,
const bool update_meshsets,
2308 moab::Interface &moab = m_field.
get_moab();
2316 const bool onlyIfImproveQuality;
2320 MergeNodes(
CoreInterface &m_field,
const bool only_if_improve_quality,
2321 Tag
th,
bool update_mehsets)
2322 : mField(m_field), onlyIfImproveQuality(only_if_improve_quality),
2323 tH(
th), updateMehsets(update_mehsets) {
2329 bool add_child =
true) {
2330 moab::Interface &moab = mField.
get_moab();
2334 CHKERR moab.get_adjacencies(conn, 2, 3,
false, vert_tets,
2335 moab::Interface::UNION);
2336 vert_tets = intersect(vert_tets, proc_tets);
2339 onlyIfImproveQuality, 0, line_search,
2344 Range::iterator lo, hi = proc_tets.begin();
2345 for (
auto pt = vert_tets.pair_begin(); pt != vert_tets.pair_end();
2347 lo = proc_tets.lower_bound(hi, proc_tets.end(), pt->first);
2348 if (lo != proc_tets.end()) {
2349 hi = proc_tets.upper_bound(lo, proc_tets.end(), pt->second);
2350 proc_tets.erase(lo, hi);
2354 proc_tets.merge(out_tets);
2358 struct ChangeChild {
2360 ChangeChild(
const EntityHandle child) : child(child) {}
2366 std::vector<
decltype(parentsChildMap.get<0>().begin())> it_vec;
2367 it_vec.reserve(parentsChildMap.size());
2369 for (
auto &
p : parent_child_map) {
2372 for (
auto it = parentsChildMap.get<0>().equal_range(
p.pArent);
2373 it.first != it.second; ++it.first)
2374 it_vec.emplace_back(it.first);
2376 for (
auto it = parentsChildMap.get<1>().equal_range(
p.pArent);
2377 it.first != it.second; ++it.first)
2378 it_vec.emplace_back(parentsChildMap.project<0>(it.first));
2380 for (
auto &it : it_vec)
2381 parentsChildMap.modify(it, ChangeChild(
p.cHild));
2384 parentsChildMap.insert(parent_child_map.begin(),
2385 parent_child_map.end());
2391 Range ¬_merged_edges,
2393 moab::Interface &moab = mField.
get_moab();
2397 std::vector<EntityHandle> parents_ents_vec(parentsChildMap.size());
2398 for (
auto &it : parentsChildMap)
2399 parents_ents_vec.emplace_back(it.pArent);
2401 parent_ents.insert_list(parents_ents_vec.begin(),
2402 parents_ents_vec.end());
2404 Range surf_parent_ents = intersect(new_surf, parent_ents);
2405 new_surf = subtract(new_surf, surf_parent_ents);
2406 Range child_surf_ents;
2407 CHKERR updateRangeByChilds(parentsChildMap, surf_parent_ents,
2409 new_surf.merge(child_surf_ents);
2411 Range edges_to_merge_parent_ents =
2412 intersect(edges_to_merge, parent_ents);
2413 edges_to_merge = subtract(edges_to_merge, edges_to_merge_parent_ents);
2414 Range merged_child_edge_ents;
2415 CHKERR updateRangeByChilds(parentsChildMap, edges_to_merge_parent_ents,
2416 merged_child_edge_ents);
2418 Range not_merged_edges_child_ents =
2419 intersect(not_merged_edges, parent_ents);
2421 subtract(not_merged_edges, not_merged_edges_child_ents);
2422 Range not_merged_child_edge_ents;
2423 CHKERR updateRangeByChilds(parentsChildMap, not_merged_edges_child_ents,
2424 not_merged_child_edge_ents);
2426 merged_child_edge_ents =
2427 subtract(merged_child_edge_ents, not_merged_child_edge_ents);
2428 edges_to_merge.merge(merged_child_edge_ents);
2429 not_merged_edges.merge(not_merged_child_edge_ents);
2431 if (updateMehsets) {
2436 CHKERR moab.get_entities_by_handle(cubit_meshset, meshset_ents,
2439 CHKERR updateRangeByChilds(parentsChildMap, meshset_ents,
2441 CHKERR moab.add_entities(cubit_meshset, child_ents);
2451 std::vector<EntityHandle> childsVec;
2458 childsVec.reserve(parents.size());
2459 for (
auto pit = parents.pair_begin(); pit != parents.pair_end(); pit++) {
2460 auto it = parent_child_map.lower_bound(pit->first);
2461 if (it != parent_child_map.end()) {
2462 for (
auto hi_it = parent_child_map.upper_bound(pit->second);
2464 childsVec.emplace_back(it->cHild);
2467 childs.insert_list(childsVec.begin(), childsVec.end());
2478 moab::Interface &moab;
2479 const double maxLength;
2481 : tH(
th), mField(m_field), moab(m_field.
get_moab()),
2482 maxLength(max_length) {
2485 acceptedThrasholdMergeQuality = 0.5;
2492 double acceptedThrasholdMergeQuality;
2497 double &ave)
const {
2500 std::array<double, 12> coords;
2502 VectorAdaptor s0(3, ublas::shallow_array_adaptor<double>(3, &coords[0]));
2503 VectorAdaptor s1(3, ublas::shallow_array_adaptor<double>(3, &coords[3]));
2506 struct NodeQuality {
2509 NodeQuality(
const EntityHandle node) : node(node), quality(1) {}
2512 typedef multi_index_container<
2513 NodeQuality, indexed_by<
2517 hashed_non_unique<tag<Ent_mi_tag>,
2519 &NodeQuality::node>>
2522 NodeQuality_sequence;
2524 NodeQuality_sequence node_quality_sequence;
2527 CHKERR moab.get_connectivity(tets, edges_nodes,
false);
2529 CHKERR moab.get_adjacencies(edges, 3,
false, edges_tets,
2530 moab::Interface::UNION);
2531 edges_tets = intersect(edges_tets, tets);
2533 for (
auto node : edges_nodes)
2534 node_quality_sequence.emplace_back(node);
2536 for (
auto tet : edges_tets) {
2538 CHKERR moab.get_connectivity(tet, conn, num_nodes,
true);
2540 CHKERR moab.tag_get_data(tH, conn, num_nodes, coords.data());
2542 CHKERR moab.get_coords(conn, num_nodes, coords.data());
2544 const double q = Tools::volumeLengthQuality(coords.data());
2546 for (
auto n : {0, 1, 2, 3}) {
2547 auto it = node_quality_sequence.get<1>().find(conn[
n]);
2548 if (it != node_quality_sequence.get<1>().end())
2549 const_cast<double &
>(it->quality) = std::min(q, it->quality);
2553 for (
auto edge : edges) {
2555 CHKERR moab.get_connectivity(edge, conn, num_nodes,
true);
2558 CHKERR moab.tag_get_data(tH, conn, num_nodes, coords.data());
2560 CHKERR moab.get_coords(conn, num_nodes, coords.data());
2563 for (
auto n : {0, 1}) {
2564 auto it = node_quality_sequence.get<1>().find(conn[
n]);
2565 if (it != node_quality_sequence.get<1>().end())
2566 q = std::min(q, it->quality);
2569 if (q < acceptedThrasholdMergeQuality) {
2570 noalias(
delta) = (s0 - s1) / maxLength;
2572 double val = pow(q, gammaQ) * pow(dot, 0.5 * gammaL);
2578 for (LengthMapData_multi_index::nth_index<0>::type::iterator mit =
2579 length_map.get<0>().begin();
2580 mit != length_map.get<0>().end(); mit++) {
2581 ave += mit->qUality;
2583 ave /= length_map.size();
2601 SURFACE_SKIN = 1 << 2,
2602 FRONT_ENDS = 1 << 3,
2604 FIX_CORNERS = 1 << 5
2607 typedef map<int, Range> SetsMap;
2610 const Range &fixed_edges,
2611 const Range &corner_nodes,
2612 const Range &constrain_surface,
2613 SetsMap &sets_map)
const {
2614 moab::Interface &moab(mField.
get_moab());
2615 Skinner skin(&moab);
2618 sets_map[FIX_CORNERS].merge(corner_nodes);
2620 CHKERR moab.get_connectivity(fixed_edges, fixed_verts,
true);
2621 sets_map[FIX_EDGES].swap(fixed_verts);
2624 CHKERR skin.find_skin(0, tets,
false, tets_skin);
2625 Range tets_skin_edges;
2626 CHKERR moab.get_adjacencies(tets_skin, 1,
false, tets_skin_edges,
2627 moab::Interface::UNION);
2630 Range constrain_surface_verts;
2631 CHKERR moab.get_connectivity(constrain_surface, constrain_surface_verts,
2633 Range constrain_surface_edges;
2634 CHKERR moab.get_adjacencies(constrain_surface, 1,
false,
2635 constrain_surface_edges,
2636 moab::Interface::UNION);
2641 Range front_in_the_body;
2642 front_in_the_body = subtract(surface_skin, tets_skin_edges);
2643 Range front_in_the_body_verts;
2644 CHKERR moab.get_connectivity(front_in_the_body, front_in_the_body_verts,
2647 CHKERR skin.find_skin(0, front_in_the_body,
false, front_ends);
2649 intersect(front_in_the_body_verts, constrain_surface_verts));
2650 sets_map[FRONT_ENDS].swap(front_ends);
2652 Range surface_skin_verts;
2653 CHKERR moab.get_connectivity(surface_skin, surface_skin_verts,
true);
2654 sets_map[SURFACE_SKIN].swap(surface_skin_verts);
2657 Range surface_verts;
2659 sets_map[SURFACE_SKIN].merge(
2660 intersect(constrain_surface_verts, surface_verts));
2661 sets_map[SURFACE].swap(surface_verts);
2664 Range tets_skin_verts;
2665 CHKERR moab.get_connectivity(tets_skin, tets_skin_verts,
true);
2666 sets_map[SKIN].swap(tets_skin_verts);
2667 sets_map[SKIN].merge(constrain_surface_verts);
2670 CHKERR moab.get_connectivity(tets, tets_verts,
true);
2671 sets_map[FREE].swap(tets_verts);
2677 Range &proc_tets)
const {
2678 moab::Interface &moab(mField.
get_moab());
2680 Range edges_to_merge_verts;
2681 CHKERR moab.get_connectivity(edges_to_merge, edges_to_merge_verts,
true);
2682 Range edges_to_merge_verts_tets;
2683 CHKERR moab.get_adjacencies(edges_to_merge_verts, 3,
false,
2684 edges_to_merge_verts_tets,
2685 moab::Interface::UNION);
2686 edges_to_merge_verts_tets = intersect(edges_to_merge_verts_tets, tets);
2687 proc_tets.swap(edges_to_merge_verts_tets);
2692 const Range &fixed_edges,
2693 const Range &corner_nodes,
2694 const Range &constrain_surface,
2695 Range &edges_to_merge,
2696 Range ¬_merged_edges) {
2697 moab::Interface &moab(mField.
get_moab());
2701 Skinner skin(&moab);
2703 CHKERR skin.find_skin(0, tets,
false, tets_skin);
2708 Range constrain_surface_verts;
2709 CHKERR moab.get_connectivity(constrain_surface, constrain_surface_verts,
2711 Range constrain_surface_edges;
2712 CHKERR moab.get_adjacencies(constrain_surface, 1,
false,
2713 constrain_surface_edges,
2714 moab::Interface::UNION);
2717 Range tets_skin_edges;
2718 CHKERR moab.get_adjacencies(tets_skin, 1,
false, tets_skin_edges,
2719 moab::Interface::UNION);
2721 Range surface_front;
2722 surface_front = subtract(surface_skin, tets_skin_edges);
2723 Range surface_front_nodes;
2724 CHKERR moab.get_connectivity(surface_front, surface_front_nodes,
true);
2727 CHKERR skin.find_skin(0, surface_front,
false, ends_nodes);
2728 ends_nodes.merge(intersect(surface_front_nodes, constrain_surface_verts));
2731 surface_skin.merge(constrain_surface);
2732 tets_skin_edges.merge(constrain_surface_edges);
2735 Range surface_edges;
2736 CHKERR moab.get_adjacencies(
surface, 1,
false, surface_edges,
2737 moab::Interface::UNION);
2739 Range surface_edges_verts;
2740 CHKERR moab.get_connectivity(surface_edges, surface_edges_verts,
true);
2742 Range tets_skin_edges_verts;
2743 CHKERR moab.get_connectivity(tets_skin_edges, tets_skin_edges_verts,
2746 Range edges_to_remove;
2750 Range ents_nodes_and_edges;
2751 ents_nodes_and_edges.merge(tets_skin_edges_verts);
2752 ents_nodes_and_edges.merge(tets_skin_edges);
2753 CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2756 edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2757 not_merged_edges.merge(edges_to_remove);
2761 Range ents_nodes_and_edges;
2762 ents_nodes_and_edges.merge(surface_edges_verts);
2763 ents_nodes_and_edges.merge(surface_edges);
2764 ents_nodes_and_edges.merge(tets_skin_edges_verts);
2765 ents_nodes_and_edges.merge(tets_skin_edges);
2766 CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2769 edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2770 not_merged_edges.merge(edges_to_remove);
2773 Range fixed_edges_nodes;
2774 CHKERR moab.get_connectivity(fixed_edges, fixed_edges_nodes,
true);
2776 Range ents_nodes_and_edges;
2777 ents_nodes_and_edges.merge(fixed_edges_nodes);
2778 ents_nodes_and_edges.merge(ends_nodes);
2779 ents_nodes_and_edges.merge(corner_nodes);
2780 ents_nodes_and_edges.merge(fixed_edges);
2781 CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2784 edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2785 not_merged_edges.merge(edges_to_remove);
2788 CHKERR removeSelfConectingEdges(surface_edges, edges_to_remove,
false);
2789 edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2790 not_merged_edges.merge(edges_to_remove);
2794 Range ents_nodes_and_edges;
2795 ents_nodes_and_edges.merge(surface_skin);
2796 ents_nodes_and_edges.merge(fixed_edges_nodes);
2797 CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2800 edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2801 not_merged_edges.merge(edges_to_remove);
2805 Range ents_nodes_and_edges;
2806 ents_nodes_and_edges.merge(surface_skin.subset_by_type(MBEDGE));
2807 ents_nodes_and_edges.merge(fixed_edges.subset_by_type(MBEDGE));
2808 CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2811 edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2812 not_merged_edges.merge(edges_to_remove);
2816 Range ents_nodes_and_edges;
2817 ents_nodes_and_edges.merge(surface_front_nodes);
2818 ents_nodes_and_edges.merge(surface_front);
2819 ents_nodes_and_edges.merge(tets_skin_edges_verts);
2820 ents_nodes_and_edges.merge(tets_skin_edges);
2821 CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2824 edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2825 not_merged_edges.merge(edges_to_remove);
2832 Range &edges_to_remove,
2834 moab::Interface &moab(mField.
get_moab());
2837 Range ents_nodes = ents.subset_by_type(MBVERTEX);
2838 if (ents_nodes.empty()) {
2839 CHKERR moab.get_connectivity(ents, ents_nodes,
true);
2842 Range ents_nodes_edges;
2843 CHKERR moab.get_adjacencies(ents_nodes, 1,
false, ents_nodes_edges,
2844 moab::Interface::UNION);
2846 Range ents_nodes_edges_nodes;
2847 CHKERR moab.get_connectivity(ents_nodes_edges, ents_nodes_edges_nodes,
2850 ents_nodes_edges_nodes = subtract(ents_nodes_edges_nodes, ents_nodes);
2851 Range ents_nodes_edges_nodes_edges;
2852 CHKERR moab.get_adjacencies(ents_nodes_edges_nodes, 1,
false,
2853 ents_nodes_edges_nodes_edges,
2854 moab::Interface::UNION);
2857 subtract(ents_nodes_edges, ents_nodes_edges_nodes_edges);
2859 subtract(ents_nodes_edges, ents.subset_by_type(MBEDGE));
2861 edges_to_remove.swap(ents_nodes_edges);
2869 Range not_merged_edges;
2871 .removeBadEdges(
surface, tets, fixed_edges, corner_nodes,
2873 Toplogy::SetsMap sets_map;
2878 for (Toplogy::SetsMap::reverse_iterator sit = sets_map.rbegin();
2879 sit != sets_map.rend(); sit++) {
2880 std::string name =
"classification_verts_" +
2881 boost::lexical_cast<std::string>(sit->first) +
".vtk";
2882 if (!sit->second.empty())
2887 CHKERR Toplogy(m_field,
th).getProcTets(tets, edges_to_merge, proc_tets);
2888 out_tets = subtract(tets, proc_tets);
2891 Range all_out_ents = out_tets;
2892 for (
int dd = 2; dd >= 0; dd--) {
2893 CHKERR moab.get_adjacencies(out_tets, dd,
false, all_out_ents,
2894 moab::Interface::UNION);
2900 int nb_nodes_merged = 0;
2904 auto save_merge_step = [&](
const int pp,
const Range collapsed_edges) {
2907 CHKERR moab.get_adjacencies(proc_tets, 2,
false, adj_faces,
2908 moab::Interface::UNION);
2910 name =
"node_merge_step_" + boost::lexical_cast<std::string>(pp) +
".vtk";
2912 name, unite(intersect(new_surf, adj_faces), collapsed_edges));
2914 "edges_to_merge_step_" + boost::lexical_cast<std::string>(pp) +
".vtk";
2916 name, unite(intersect(new_surf, adj_faces), edges_to_merge));
2923 double ave0 = 0, ave = 0, min = 0, min_p = 0, min_pp;
2926 int nb_nodes_merged_p = nb_nodes_merged;
2933 if(!length_map.empty())
2934 min = length_map.get<2>().begin()->qUality;
2940 Range collapsed_edges;
2941 MergeNodes merge_nodes(m_field,
true,
th, update_meshsets);
2943 for (
auto mit = length_map.get<0>().begin();
2944 mit != length_map.get<0>().end(); mit++, nn++) {
2948 auto get_father_and_mother =
2953 CHKERR moab.get_connectivity(mit->eDge, conn, num_nodes,
true);
2954 std::array<int, 2> conn_type = {0, 0};
2955 for (
int nn = 0; nn != 2; nn++) {
2956 for (Toplogy::SetsMap::reverse_iterator sit = sets_map.rbegin();
2957 sit != sets_map.rend(); sit++) {
2958 if (sit->second.find(conn[nn]) != sit->second.end()) {
2959 conn_type[nn] |= sit->first;
2963 int type_father, type_mother;
2964 if (conn_type[0] > conn_type[1]) {
2967 type_father = conn_type[0];
2968 type_mother = conn_type[1];
2972 type_father = conn_type[1];
2973 type_mother = conn_type[0];
2975 if (type_father == type_mother) {
2981 int line_search = 0;
2983 CHKERR get_father_and_mother(father, mother, line_search);
2984 CHKERR merge_nodes.mergeNodes(line_search, father, mother, proc_tets);
2986 const EntityHandle father_and_mother[] = {father, mother};
2988 CHKERR moab.get_adjacencies(father_and_mother, 1, 3,
false, adj_tets);
2989 Range adj_tets_nodes;
2990 CHKERR moab.get_connectivity(adj_tets, adj_tets_nodes,
true);
2992 CHKERR moab.get_adjacencies(adj_tets_nodes, 1,
false, adj_edges,
2993 moab::Interface::UNION);
2994 for (
auto ait : adj_edges) {
2995 auto miit = length_map.get<1>().find(ait);
2996 if (miit != length_map.get<1>().end())
3000 collapsed_edges.insert(mit->eDge);
3003 if (nn >
static_cast<int>(length_map.size() / fraction_level))
3005 if (mit->qUality > ave)
3010 CHKERR merge_nodes.updateRangeByChilds(new_surf, edges_to_merge,
3011 not_merged_edges,
true);
3014 "(%d) Number of nodes merged %d ave q %3.4e min q %3.4e", pp,
3015 nb_nodes_merged, ave, min);
3018 CHKERR save_merge_step(pp + 1, collapsed_edges);
3020 if (nb_nodes_merged == nb_nodes_merged_p)
3022 if (min > 1e-2 && min == min_pp)
3028 CHKERR moab.get_adjacencies(proc_tets, 1,
false, adj_edges,
3029 moab::Interface::UNION);
3030 edges_to_merge = intersect(edges_to_merge, adj_edges);
3032 .removeBadEdges(new_surf, proc_tets, fixed_edges, corner_nodes,
3036 auto reconstruct_refined_ents = [&]() {
3047 CHKERR reconstruct_refined_ents();
3053 out_tets.merge(proc_tets);
3055 CHKERR moab.get_adjacencies(out_tets, 2,
false, adj_faces,
3056 moab::Interface::UNION);
3057 new_surf = intersect(new_surf, adj_faces);
3063 const int fraction_level,
const BitRefLevel cut_bit,
3065 const Range &fixed_edges,
const Range &corner_nodes, Tag
th,
3066 const bool update_meshsets,
const bool debug) {
3071 trim_bit,
BitRefLevel().set(), MBTET, tets_level);
3073 Range edges_to_merge;
3075 trim_bit, cut_bit | trim_bit, edges_to_merge);
3078 Range all_ents_not_in_database_before;
3080 all_ents_not_in_database_before);
3082 edges_to_merge = edges_to_merge.subset_by_type(MBEDGE);
3086 Range out_new_tets, out_new_surf;
3088 corner_nodes, edges_to_merge, out_new_tets, out_new_surf,
3092 Range all_ents_not_in_database_after;
3094 all_ents_not_in_database_after);
3097 all_ents_not_in_database_after =
3098 subtract(all_ents_not_in_database_after, all_ents_not_in_database_before);
3100 CHKERR m_field.
get_moab().get_entities_by_type(0, MBENTITYSET, meshsets,
3102 for (
auto m : meshsets)
3104 all_ents_not_in_database_after);
3106 m_field.
get_moab().delete_entities(all_ents_not_in_database_after);
3116 moab::Interface &moab = m_field.
get_moab();
3120 CHKERR moab.get_entities_by_type(0, MBVERTEX, nodes);
3124 std::vector<double> coords(3 * nodes.size());
3125 CHKERR moab.get_coords(nodes, &coords[0]);
3126 CHKERR moab.tag_set_data(
th, nodes, &coords[0]);
3133 moab::Interface &moab = m_field.
get_moab();
3137 CHKERR moab.get_entities_by_type(0, MBVERTEX, nodes);
3140 bit, mask, MBVERTEX, nodes);
3141 std::vector<double> coords(3 * nodes.size());
3142 CHKERR moab.tag_get_data(
th, nodes, &coords[0]);
3143 CHKERR moab.set_coords(nodes, &coords[0]);
3149 moab::Interface &moab = m_field.
get_moab();
#define CutMeshFunctionBegin
#define MOFEM_LOG_C(channel, severity, format,...)
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
@ MOFEM_DATA_INCONSISTENCY
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
FTensor::Index< 'n', SPACE_DIM > n
FTensor::Index< 'm', SPACE_DIM > m
virtual const RefEntity_multiIndex * get_ref_ents() const =0
Get the ref ents object.
MoFEMErrorCode addBitRefLevel(const Range &ents, const BitRefLevel &bit, int verb=QUIET) const
add bit ref level to ref entity
MoFEMErrorCode updateRangeByChildren(const Range &parent, Range &child, MoFEMTypes bh=MF_ZERO)
Update range by childrens.
#define MOFEM_LOG(channel, severity)
Log.
#define _IT_CUBITMESHSETS_FOR_LOOP_(MESHSET_MANAGER, IT)
Iterator that loops over all the Cubit MeshSets in a moFEM field.
FTensor::Index< 'i', SPACE_DIM > i
const double v
phase velocity of light in medium (cm/ns)
FTensor::Index< 'l', 3 > l
FTensor::Index< 'j', 3 > j
FTensor::Index< 'k', 3 > k
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
VectorBoundedArray< double, 3 > VectorDouble3
VectorBoundedArray< double, 6 > VectorDouble6
VectorShallowArrayAdaptor< double > VectorAdaptor
VectorBoundedArray< double, 12 > VectorDouble12
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
MatrixShallowArrayAdaptor< double > MatrixAdaptor
Matrix adaptor.
implementation of Data Operators for Forces and Sources
auto getVectorAdaptor(T1 ptr, const size_t n)
Get Vector adaptor.
MoFEMErrorCode reconstructMultiIndex(const MI &mi, MO &&mo=Modify_change_nothing())
Template used to reconstruct multi-index.
constexpr double t
plate stiffness
static constexpr double delta
virtual moab::Interface & get_moab()=0
MoFEMErrorCode findEdgesToTrim(Range *fixed_edges, Range *corner_nodes, Tag th=NULL, const double tol=1e-4, const bool debug=false)
Find edges to trimEdges.
MoFEMErrorCode copySurface(const Range surface, Tag th=NULL, double *shift=NULL, double *origin=NULL, double *transform=NULL, const std::string save_mesh="")
copy surface entities
MoFEMErrorCode setFront(const Range surface)
MoFEMErrorCode findLevelSetVolumes(Tag th, Range &vol_edges, const bool remove_adj_prims_edges, int verb=QUIET, const bool debug=false, const std::string edges_file_name=string())
Find level set volumes.
MoFEMErrorCode moveMidNodesOnTrimmedEdges(Tag th=NULL)
move trimmed edges mid nodes
MoFEMErrorCode trimOnly(const BitRefLevel trim_bit, Tag th, const double tol_trim_close, Range *fixed_edges=NULL, Range *corner_nodes=NULL, const bool update_meshsets=false, const bool debug=false)
Trim mesh only.
MoFEMErrorCode trimSurface(Range *fixed_edge, Range *corner_nodes, const double tol=1e-4, Tag th=NULL, const bool debug=false)
Trim surface from faces beyond front.
double aveLength
Average edge length.
MoFEMErrorCode trimEdgesInTheMiddle(const BitRefLevel bit, const bool debug=false)
trim edges
MoFEMErrorCode setCoords(Tag th, const BitRefLevel bit=BitRefLevel(), const BitRefLevel mask=BitRefLevel().set())
set coords from tag
const Range & projectZeroDistanceEnts() const
MoFEMErrorCode findEdgesToCut(Range vol, int verb=QUIET, const bool debug=false)
find edges to cut
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
MoFEMErrorCode cutEdgesInMiddle(const BitRefLevel bit, Range &cut_vols, Range &cut_surf, Range &cut_verts, const bool debug=false)
cut edges
MoFEMErrorCode setSurface(const Range surface)
set surface entities
map< EntityHandle, TreeData > edgesToTrim
map< EntityHandle, TreeData > verticesOnTrimEdges
MoFEMErrorCode buildTree()
build tree
double projectEntitiesQualityTrashold
double maxLength
Maximal edge length.
MoFEMErrorCode splitSides(const BitRefLevel split_bit, const BitRefLevel bit, const Range &ents, Tag th=NULL)
split sides
MoFEMErrorCode mergeBadEdges(const int fraction_level, const Range &tets, const Range &surface, const Range &fixed_edges, const Range &corner_nodes, Range &merged_nodes, Range &out_tets, Range &new_surf, Tag th, const bool update_meshsets=false, const BitRefLevel *bit_ptr=NULL, const bool debug=false)
Merge edges.
MoFEMErrorCode mergeVolumes(const Range volume)
merge volume entities
MoFEMErrorCode moveMidNodesOnCutEdges(Tag th=NULL)
projecting of mid edge nodes on new mesh on surface
map< EntityHandle, TreeData > verticesOnCutEdges
multi_index_container< LengthMapData, indexed_by< ordered_non_unique< member< LengthMapData, double, &LengthMapData::lEngth > >, hashed_unique< member< LengthMapData, EntityHandle, &LengthMapData::eDge > >, ordered_non_unique< member< LengthMapData, double, &LengthMapData::qUality > > > > LengthMapData_multi_index
MoFEMErrorCode saveCutEdges(const std::string prefix="")
MoFEMErrorCode saveTrimEdges()
MoFEMErrorCode snapSurfaceSkinToEdges(const Range fixed_edges, const double rel_tol, const double abs_tol, Tag th=nullptr, const bool debug=false)
MoFEMErrorCode clearMap()
MoFEMErrorCode removePathologicalFrontTris(const BitRefLevel split_bit, Range &ents)
Remove pathological elements on surface internal front.
MoFEMErrorCode setConstrainSurface(const Range surf)
Set the constrain surface object.
MoFEMErrorCode cutAndTrim(int &first_bit, Tag th, const double tol_geometry, const double tol_cut_close, const double tol_trim_close, Range *fixed_edges=NULL, Range *corner_nodes=NULL, const bool update_meshsets=false, const bool debug=false)
Cut and trim.
const Range & getNewTrimSurfaces() const
MoFEMErrorCode cutOnly(Range vol, const BitRefLevel cut_bit, Tag th, const double tol_geometry, const double tol_cut_close, Range *fixed_edges=NULL, Range *corner_nodes=NULL, const bool update_meshsets=false, const bool debug=false)
Cut mesh only.
MoFEMErrorCode setVolume(const Range volume)
set volume entities
MoFEMErrorCode createFrontLevelSets(Range vol, Tag th=nullptr, int verb=QUIET, const bool debug=false)
Calculate distance from mesh nodes to surface front.
MoFEMErrorCode setTagData(Tag th, const BitRefLevel bit=BitRefLevel())
set coords to tag
CutMeshInterface(const MoFEM::Core &core)
MoFEMErrorCode cutTrimAndMerge(int &first_bit, const int fraction_level, Tag th, const double tol_geometry, const double tol_cut_close, const double tol_trim_close, Range &fixed_edges, Range &corner_nodes, const bool update_meshsets=false, const bool debug=false)
Cut, trim and merge.
map< EntityHandle, TreeData > edgesToCut
MoFEMErrorCode makeFront(const bool debug=false)
Create front from the surface.
boost::shared_ptr< OrientedBoxTreeTool > treeSurfPtr
MoFEMErrorCode mergeSurface(const Range surface)
merge surface entities
MoFEMErrorCode refineMesh(const int init_bit_level, const int surf_levels, const int front_levels, Range *fixed_edges=nullptr, int verb=QUIET, const bool debug=false)
Refine and set level sets.
MoFEMErrorCode snapSurfaceToEdges(const Range surface_edges, const Range fixed_edges, const double rel_tol, const double abs_tol, Tag th=nullptr, const bool debug=false)
MoFEMErrorCode createSurfaceLevelSets(int verb=QUIET, const bool debug=false)
Calculate distance from mesh nodes to cut surface.
Mesh refinement interface.
MoFEMErrorCode refineTets(const EntityHandle meshset, const BitRefLevel &bit, int verb=QUIET, const bool debug=false)
refine TET in the meshset
MoFEMErrorCode addVerticesInTheMiddleOfEdges(const EntityHandle meshset, const BitRefLevel &bit, const bool recursive=false, int verb=QUIET, EntityHandle start_v=0)
make vertices in the middle of edges in meshset and add them to refinement levels defined by bit
Interface for managing meshsets containing materials and boundary conditions.
Merge node by collapsing edge between them.
multi_index_container< ParentChild, indexed_by< ordered_unique< member< ParentChild, EntityHandle, &ParentChild::pArent > >, ordered_non_unique< member< ParentChild, EntityHandle, &ParentChild::cHild > > > > ParentChildMap
bool getSuccessMerge()
Return true if successful merge.
MoFEMErrorCode mergeNodes(EntityHandle father, EntityHandle mother, Range &out_tets, Range *tets_ptr=NULL, const bool only_if_improve_quality=false, const double move=0, const int line_search=0, Tag th=NULL, const int verb=0)
merge nodes which sharing edge
ParentChildMap & getParentChildMap()
Get map of parent cand child.
Create interface from given surface and insert flat prisms in-between.
MoFEMErrorCode findFacesWithThreeNodesOnInternalSurfaceSkin(const EntityHandle sideset, const BitRefLevel mesh_bit_level, const bool recursive, Range &faces_with_three_nodes_on_front, int verb=QUIET)
Find triangles which have three nodes on internal surface skin.
MoFEMErrorCode getSides(const int msId, const CubitBCType cubit_bc_type, const BitRefLevel mesh_bit_level, const bool recursive, int verb=QUIET)
Store tetrahedra from each side of the interface separately in two child meshsets of the parent meshs...
MoFEMErrorCode splitSides(const EntityHandle meshset, const BitRefLevel &bit, const int msId, const CubitBCType cubit_bc_type, const bool add_interface_entities, const bool recursive=false, int verb=QUIET)
Split nodes and other entities of tetrahedra on both sides of the interface and insert flat prisms in...
base class for all interface classes
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface refernce to pointer of interface.