17 return std::min(
a, b);
22 successMerge(false), errorIfNoCommonEdge(false) {}
27 ierr = PetscOptionsBegin(m_field.
get_comm(),
"node_merge",
28 "Node merge options",
"none");
31 ierr = PetscOptionsEnd();
39 const bool only_if_improve_quality,
const double move,
40 const int line_search, Tag
th,
int verb) {
47 CHKERR m_field.
get_moab().get_adjacencies(&father, 1, 1,
false, father_edges);
49 CHKERR m_field.
get_moab().get_adjacencies(&mother, 1, 1,
false, mother_edges);
53 CHKERR m_field.
get_moab().get_adjacencies(&father, 1, 3,
false, father_tets);
55 CHKERR m_field.
get_moab().get_adjacencies(&mother, 1, 3,
false, mother_tets);
56 if (tets_ptr != NULL) {
57 mother_tets = intersect(mother_tets, *tets_ptr);
58 father_tets = intersect(father_tets, *tets_ptr);
63 common_edge = intersect(father_edges, mother_edges);
64 if (tets_ptr != NULL) {
65 Range tets = unite(father_tets, mother_tets);
67 CHKERR m_field.
get_moab().get_adjacencies(tets, 1,
false, tets_edges,
68 moab::Interface::UNION);
69 common_edge = intersect(common_edge, tets_edges);
70 father_edges = intersect(father_edges, tets_edges);
71 mother_edges = intersect(mother_edges, tets_edges);
77 "no common edge between nodes");
78 }
else if (common_edge.empty()) {
80 if (tets_ptr != NULL) {
81 seed_tets.merge(*tets_ptr);
90 CHKERR m_field.
get_moab().get_adjacencies(common_edge, 3,
true, edge_tets);
92 if (tets_ptr != NULL) {
93 edge_tets = intersect(edge_tets, *tets_ptr);
96 mother_tets = subtract(mother_tets, edge_tets);
97 father_tets = subtract(father_tets, edge_tets);
100 const int num_nodes) {
103 CHKERR m_field.
get_moab().get_coords(conn, num_nodes, &coords[0]);
110 auto get_tensor = [](
VectorDouble &coords,
const int shift) {
112 &coords[shift], &coords[shift + 1], &coords[shift + 2]);
120 auto t_n0 = get_tensor(coords, 0);
121 auto t_n1 = get_tensor(coords, 3);
122 t_move(
i) = t_n0(
i) + move * (t_n1(
i) - t_n0(
i));
125 if (line_search > 0) {
130 if (only_if_improve_quality) {
133 if (move > 0 || line_search) {
137 double min_quality0 = 1;
142 double min_quality = 1;
144 ((move > 0) || line_search) ? &t_move(0) : NULL,
146 if (min_quality < min_quality0) {
147 if (tets_ptr != NULL) {
148 out_tets = *tets_ptr;
156 if (move > 0 || line_search) {
165 int *ret_num_nodes =
nullptr) {
168 CHKERR m_field.
get_moab().get_connectivity(ent, conn, num_nodes,
true);
170 *ret_num_nodes = num_nodes;
174 auto create_tet = [
this, &m_field](
const EntityHandle *new_conn,
178 CHKERR m_field.
get_moab().get_adjacencies(new_conn, 4, 3,
false, tets);
179 bool tet_found =
false;
180 for (
auto it_tet : tets) {
183 CHKERR m_field.
get_moab().get_connectivity(it_tet, tet_conn, num_nodes,
185 const EntityHandle *p = std::find(tet_conn, &tet_conn[4], new_conn[0]);
186 if (p != tet_conn + 4) {
187 int s = std::distance(tet_conn, p);
189 for (;
n != 4; ++
n) {
190 const int idx[] = {0, 1, 2, 3, 0, 1, 2, 3};
191 if (tet_conn[idx[s +
n]] != new_conn[
n])
194 if (
n == 4 && !tet_found) {
198 THROW_MESSAGE(
"More that one tet with the same connectivity");
204 CHKERR m_field.
get_moab().create_element(MBTET, new_conn, 4, tet);
216 for (
auto m_tet : mother_tets) {
220 int nb_mother_verts = 0;
221 for (
int nn = 0; nn != 4; ++nn) {
222 if (conn[nn] == father) {
224 "Tet has father vertex, impossible but here it is");
226 if (conn[nn] == mother) {
227 new_conn[nn] = father;
230 new_conn[nn] = conn[nn];
233 if (nb_mother_verts != 1) {
235 "Tet should have only one vertex but have %d", nb_mother_verts);
241 created_tets.insert(create_tet(new_conn, m_tet));
247 for (
int dd = 1;
dd <= 2;
dd++)
249 moab::Interface::UNION);
252 auto adj_crated_ents =
get_adj_ents(created_tets,
true);
253 adj_crated_ents.erase(common_edge[0]);
256 for (
auto ent : adj_crated_ents) {
261 bool father_node =
false;
262 for (
int nn = 0; nn != num_nodes; nn++) {
263 if (conn[nn] == father)
266 small_conn[ii++] = conn[nn];
270 std::sort(&small_conn[0], &small_conn[ii]);
272 face_map.insert(
FaceMap(ent, small_conn[0], small_conn[1]));
274 face_map.insert(
FaceMap(ent, small_conn[0], 0));
278 auto adj_mother_ents =
get_adj_ents(mother_tets,
false);
279 adj_mother_ents.erase(common_edge[0]);
280 for (
auto ent : adj_mother_ents) {
287 for (; nn != num_nodes; ++nn) {
288 if (conn[nn] == mother) {
291 small_conn[ii++] = conn[nn];
294 if (nb_new_node > 0) {
296 std::sort(&small_conn[0], &small_conn[ii]);
302 FaceMapIdx::iterator fit = face_map.find(boost::make_tuple(n0, n1));
303 if (fit != face_map.end()) {
306 if (m_field.
get_moab().dimension_from_handle(parent) !=
307 m_field.
get_moab().dimension_from_handle(child))
309 "Huston we have a problem!");
322 if (tets_ptr != NULL) {
323 seed_tets.merge(*tets_ptr);
324 seed_tets = subtract(seed_tets, edge_tets);
325 seed_tets = subtract(seed_tets, mother_tets);
327 seed_tets.merge(created_tets);
328 out_tets.swap(seed_tets);
333 std::cerr <<
"Nodes merged" << endl;
341 double &min_quality, Tag
th,
342 boost::function<
double(
double,
double)>
f) {
351 CHKERR m_field.
get_moab().get_connectivity(*tit, conn, num_nodes,
true);
355 int nb_mother_verts = 0;
357 for (
int nn = 0; nn < 4; nn++) {
358 if (conn[nn] == father) {
361 if (conn[nn] == mother) {
362 new_conn[nn] = father;
366 new_conn[nn] = conn[nn];
369 if (nb_mother_verts > 1) {
371 "Tet should have no more than one mother vertex but have %d",
375 CHKERR m_field.
get_moab().get_coords(new_conn, num_nodes, coords);
380 int shift = 3 * father_nn;
381 for (
int nn = 0; nn != 3; nn++) {
382 coords[shift + nn] = coords_move[nn];
393 if (!std::isnormal(quality)) {
397 min_quality =
f(min_quality, quality);
421 for (
int nn = 0; nn != 3; nn++) {
422 t_delta(nn) = coords[3 + nn] - t_coords(nn);
425 t_move(
i) = t_coords(
i);
426 double min_quality_i = 1;
430 t_move(
i) = t_coords(
i) + t_delta(
i);
431 double min_quality_k = 1;
438 for (
int ii = 0; ii != line_search; ii++) {
439 double min_quality = 1;
440 double alpha = (alpha_i + alpha_k) * 0.5;
441 t_move(
i) = t_coords(
i) + alpha * t_delta(
i);
444 if (min_quality_i >= min_quality_k) {
445 min_quality_k = min_quality;
448 min_quality_i = min_quality;
453 if (min_quality_i > min_quality_k) {
454 t_move(
i) = t_coords(
i) + alpha_i * t_delta(
i);
456 t_move(
i) = t_coords(
i) + alpha_k * t_delta(
i);
464 const bool only_if_improve_quality,
const double move, Tag
th) {
468 CHKERR mergeNodes(father, mother, out_tets, tets_ptr, only_if_improve_quality,
476 BitRefLevel tets_from_bit_ref_level,
const bool only_if_improve_quality,
477 const double move, Tag
th) {
482 tets_from_bit_ref_level,
BitRefLevel().set(), MBTET, level_tets);