v0.15.0
Loading...
Searching...
No Matches
NodeMerger.cpp
Go to the documentation of this file.
1/** \file NodeMerger.cpp
2 * \brief Interface for merging nodes
3 */
4
5
6
7namespace MoFEM {
8
10NodeMergerInterface::query_interface(boost::typeindex::type_index type_index,
11 UnknownInterface **iface) const {
12 *iface = const_cast<NodeMergerInterface *>(this);
13 return 0;
14}
15
16static auto min_non_abs(const double a, const double b) {
17 return std::min(a, b);
18};
19
21 : cOre(const_cast<MoFEM::Core &>(core)), minQualityFunction(min_non_abs),
22 successMerge(false), errorIfNoCommonEdge(false) {}
23
25 Interface &m_field = cOre;
27 PetscOptionsBegin(m_field.get_comm(), "node_merge",
28 "Node merge options", "none");
29 PetscOptionsEnd();
30
32}
33
35 EntityHandle father, EntityHandle mother, Range &out_tets, Range *tets_ptr,
36 const bool only_if_improve_quality, const double move,
37 const int line_search, Tag th, int verb) {
38 Interface &m_field = cOre;
39 FTensor::Index<'i', 3> i;
41
42 // Get edges adjacent to father and mother, i.e. mother is merged to father.
43 Range father_edges;
44 CHKERR m_field.get_moab().get_adjacencies(&father, 1, 1, false, father_edges);
45 Range mother_edges;
46 CHKERR m_field.get_moab().get_adjacencies(&mother, 1, 1, false, mother_edges);
47
48 // Get tets adjacent to mother and father
49 Range father_tets;
50 CHKERR m_field.get_moab().get_adjacencies(&father, 1, 3, false, father_tets);
51 Range mother_tets;
52 CHKERR m_field.get_moab().get_adjacencies(&mother, 1, 3, false, mother_tets);
53 if (tets_ptr != NULL) {
54 mother_tets = intersect(mother_tets, *tets_ptr);
55 father_tets = intersect(father_tets, *tets_ptr);
56 }
57
58 // Find common edge
59 Range common_edge;
60 common_edge = intersect(father_edges, mother_edges);
61 if (tets_ptr != NULL) {
62 Range tets = unite(father_tets, mother_tets);
63 Range tets_edges;
64 CHKERR m_field.get_moab().get_adjacencies(tets, 1, false, tets_edges,
65 moab::Interface::UNION);
66 common_edge = intersect(common_edge, tets_edges);
67 father_edges = intersect(father_edges, tets_edges);
68 mother_edges = intersect(mother_edges, tets_edges);
69 }
70
71 // No common edge, merge no possible
72 if (errorIfNoCommonEdge && common_edge.empty()) {
73 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
74 "no common edge between nodes");
75 } else if (common_edge.empty()) {
76 Range seed_tets;
77 if (tets_ptr != NULL) {
78 seed_tets.merge(*tets_ptr);
79 }
80 out_tets = seed_tets;
81 successMerge = false;
83 }
84
85 // Common edge tets, that tests will be squashed
86 Range edge_tets;
87 CHKERR m_field.get_moab().get_adjacencies(common_edge, 3, true, edge_tets);
88 // Intersect with ptr_tets (usually associated with some bit level)
89 if (tets_ptr != NULL) {
90 edge_tets = intersect(edge_tets, *tets_ptr);
91 }
92 // Mother tets, has only one mother vertex and no father vertex.
93 mother_tets = subtract(mother_tets, edge_tets);
94 father_tets = subtract(father_tets, edge_tets);
95
96 auto get_coords = [&m_field](Tag th, const EntityHandle *conn,
97 const int num_nodes) {
98 VectorDouble coords(3 * num_nodes);
99 if (th == NULL) {
100 CHKERR m_field.get_moab().get_coords(conn, num_nodes, &coords[0]);
101 } else {
102 CHKERR m_field.get_moab().tag_get_data(th, conn, num_nodes, &coords[0]);
103 }
104 return coords;
105 };
106
107 auto get_tensor = [](VectorDouble &coords, const int shift) {
109 &coords[shift], &coords[shift + 1], &coords[shift + 2]);
110 };
111
112 // move father coord is move > 0
114 if (move > 0) {
115 EntityHandle conn[] = {father, mother};
116 VectorDouble coords = get_coords(th, conn, 2);
117 auto t_n0 = get_tensor(coords, 0);
118 auto t_n1 = get_tensor(coords, 3);
119 t_move(i) = t_n0(i) + move * (t_n1(i) - t_n0(i));
120 }
121
122 if (line_search > 0) {
123 Range check_tests = unite(father_tets, mother_tets);
124 CHKERR lineSearch(check_tests, father, mother, line_search, t_move, th);
125 }
126
127 if (only_if_improve_quality) {
128 Range check_tests = mother_tets;
129 // no need to check father tets since no change of quality for them
130 if (move > 0 || line_search) {
131 check_tests.merge(father_tets);
132 }
133
134 double min_quality0 = 1;
135 CHKERR minQuality(edge_tets, 0, 0, NULL, min_quality0, th,
137 CHKERR minQuality(check_tests, 0, 0, NULL, min_quality0, th,
139 double min_quality = 1;
140 CHKERR minQuality(check_tests, father, mother,
141 ((move > 0) || line_search) ? &t_move(0) : NULL,
142 min_quality, th, minQualityFunction);
143 if (min_quality < min_quality0) {
144 if (tets_ptr != NULL) {
145 out_tets = *tets_ptr;
146 }
147 successMerge = false;
149 }
150 }
151
152 // Move node
153 if (move > 0 || line_search) {
154 if (th == NULL) {
155 CHKERR m_field.get_moab().set_coords(&father, 1, &t_move(0));
156 } else {
157 CHKERR m_field.get_moab().tag_set_data(th, &father, 1, &t_move(0));
158 }
159 }
160
161 auto get_conn = [&m_field](const EntityHandle ent,
162 int *ret_num_nodes = nullptr) {
163 int num_nodes;
164 const EntityHandle *conn;
165 CHKERR m_field.get_moab().get_connectivity(ent, conn, num_nodes, true);
166 if (ret_num_nodes)
167 *ret_num_nodes = num_nodes;
168 return conn;
169 };
170
171 auto create_tet = [this, &m_field](const EntityHandle *new_conn,
172 const EntityHandle parent) {
173 EntityHandle tet;
174 Range tets;
175 CHKERR m_field.get_moab().get_adjacencies(new_conn, 4, 3, false, tets);
176 bool tet_found = false;
177 for (auto it_tet : tets) {
178 const EntityHandle *tet_conn;
179 int num_nodes;
180 CHKERR m_field.get_moab().get_connectivity(it_tet, tet_conn, num_nodes,
181 true);
182 const EntityHandle *p = std::find(tet_conn, &tet_conn[4], new_conn[0]);
183 if (p != tet_conn + 4) {
184 int s = std::distance(tet_conn, p);
185 int n = 0;
186 for (; n != 4; ++n) {
187 const int idx[] = {0, 1, 2, 3, 0, 1, 2, 3};
188 if (tet_conn[idx[s + n]] != new_conn[n])
189 break;
190 }
191 if (n == 4 && !tet_found) {
192 tet = it_tet;
193 tet_found = true;
194 } else if (n == 4) {
195 THROW_MESSAGE("More that one tet with the same connectivity");
196 }
197 }
198 }
199 if (!tet_found) {
200 // Create tet with new connectivity
201 CHKERR m_field.get_moab().create_element(MBTET, new_conn, 4, tet);
202 CHKERR m_field.get_moab().tag_set_data(cOre.get_th_RefParentHandle(),
203 &tet, 1, &parent);
204 parentChildMap.insert(ParentChild(parent, tet));
205 }
206 return tet;
207 };
208
209 // clear map
210 parentChildMap.clear();
211
212 Range created_tets;
213 for (auto m_tet : mother_tets) {
214 const EntityHandle *conn = get_conn(m_tet);
215 EntityHandle new_conn[4];
216 // Replace mother vertices by father vertices
217 int nb_mother_verts = 0;
218 for (int nn = 0; nn != 4; ++nn) {
219 if (conn[nn] == father) {
220 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
221 "Tet has father vertex, impossible but here it is");
222 }
223 if (conn[nn] == mother) {
224 new_conn[nn] = father;
225 nb_mother_verts++;
226 } else {
227 new_conn[nn] = conn[nn];
228 }
229 }
230 if (nb_mother_verts != 1) {
231 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
232 "Tet should have only one vertex but have %d", nb_mother_verts);
233 }
234
235 VectorDouble new_coords = get_coords(th, new_conn, 4);
236
237 // add tet to range
238 created_tets.insert(create_tet(new_conn, m_tet));
239 }
240
241 // Loop over mother adjacent entities to use them as parents
242 auto get_adj_ents = [&](const Range &ents, const bool create) {
243 Range adj;
244 for (int dd = 1; dd <= 2; dd++)
245 CHKERR m_field.get_moab().get_adjacencies(ents, dd, create, adj,
246 moab::Interface::UNION);
247 return adj;
248 };
249 auto adj_crated_ents = get_adj_ents(created_tets, true);
250 adj_crated_ents.erase(common_edge[0]);
251
252 FaceMapIdx face_map;
253 for (auto ent : adj_crated_ents) {
254 int num_nodes;
255 const EntityHandle *conn = get_conn(ent, &num_nodes);
256 EntityHandle small_conn[num_nodes];
257 int ii = 0;
258 bool father_node = false;
259 for (int nn = 0; nn != num_nodes; nn++) {
260 if (conn[nn] == father)
261 father_node = true;
262 else
263 small_conn[ii++] = conn[nn];
264 }
265 if (father_node) {
266 if (ii > 1)
267 std::sort(&small_conn[0], &small_conn[ii]);
268 if (ii == 2)
269 face_map.insert(FaceMap(ent, small_conn[0], small_conn[1]));
270 else
271 face_map.insert(FaceMap(ent, small_conn[0], 0));
272 }
273 }
274
275 auto adj_mother_ents = get_adj_ents(mother_tets, false);
276 adj_mother_ents.erase(common_edge[0]);
277 for (auto ent : adj_mother_ents) {
278 int num_nodes;
279 const EntityHandle *conn = get_conn(ent, &num_nodes);
280 EntityHandle small_conn[num_nodes];
281 int nb_new_node = 0;
282 int nn = 0;
283 int ii = 0;
284 for (; nn != num_nodes; ++nn) {
285 if (conn[nn] == mother) {
286 nb_new_node++;
287 } else {
288 small_conn[ii++] = conn[nn];
289 }
290 }
291 if (nb_new_node > 0) {
292 if (ii > 1)
293 std::sort(&small_conn[0], &small_conn[ii]);
294
295 EntityHandle n0 = small_conn[0], n1 = 0;
296 if (ii == 2)
297 n1 = small_conn[1];
298
299 FaceMapIdx::iterator fit = face_map.find(boost::make_tuple(n0, n1));
300 if (fit != face_map.end()) {
301 const EntityHandle child = fit->e;
302 const EntityHandle parent = ent;
303 if (m_field.get_moab().dimension_from_handle(parent) !=
304 m_field.get_moab().dimension_from_handle(child))
305 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
306 "Huston we have a problem!");
307
308 // Set parent child relation
309 CHKERR m_field.get_moab().tag_set_data(cOre.get_th_RefParentHandle(),
310 &child, 1, &parent);
311 // Create map
312 parentChildMap.insert(ParentChild(parent, child));
313 }
314 }
315 }
316
317 // Seed tets to given bit level
318 Range seed_tets;
319 if (tets_ptr != NULL) {
320 seed_tets.merge(*tets_ptr);
321 seed_tets = subtract(seed_tets, edge_tets);
322 seed_tets = subtract(seed_tets, mother_tets);
323 }
324 seed_tets.merge(created_tets);
325 out_tets.swap(seed_tets);
326
327 successMerge = true;
328
329 if (verb > VERY_VERBOSE)
330 std::cerr << "Nodes merged" << endl;
331
333}
334
337 EntityHandle mother, double *coords_move,
338 double &min_quality, Tag th,
339 boost::function<double(double, double)> f) {
340 Interface &m_field = cOre;
341 double coords[12];
342 FTensor::Index<'i', 3> i;
344 for (Range::iterator tit = check_tests.begin(); tit != check_tests.end();
345 tit++) {
346 const EntityHandle *conn;
347 int num_nodes;
348 CHKERR m_field.get_moab().get_connectivity(*tit, conn, num_nodes, true);
349 if (mother > 0) {
350 EntityHandle new_conn[4];
351 // Replace mother vertices by father vertices
352 int nb_mother_verts = 0;
353 int father_nn = 0;
354 for (int nn = 0; nn < 4; nn++) {
355 if (conn[nn] == father) {
356 father_nn = nn;
357 }
358 if (conn[nn] == mother) {
359 new_conn[nn] = father;
360 father_nn = nn;
361 nb_mother_verts++;
362 } else {
363 new_conn[nn] = conn[nn];
364 }
365 }
366 if (nb_mother_verts > 1) {
367 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
368 "Tet should have no more than one mother vertex but have %d",
369 nb_mother_verts);
370 }
371 if (th == NULL) {
372 CHKERR m_field.get_moab().get_coords(new_conn, num_nodes, coords);
373 } else {
374 CHKERR m_field.get_moab().tag_get_data(th, new_conn, num_nodes, coords);
375 }
376 if (coords_move) {
377 int shift = 3 * father_nn;
378 for (int nn = 0; nn != 3; nn++) {
379 coords[shift + nn] = coords_move[nn];
380 }
381 }
382 } else {
383 if (th == NULL) {
384 CHKERR m_field.get_moab().get_coords(conn, num_nodes, coords);
385 } else {
386 CHKERR m_field.get_moab().tag_get_data(th, conn, num_nodes, coords);
387 }
388 }
389 double quality = Tools::volumeLengthQuality(coords);
390 if (!std::isnormal(quality)) {
391 min_quality = 0;
393 }
394 min_quality = f(min_quality, quality);
395 }
397};
398
401 EntityHandle mother, int line_search,
403 Interface &m_field = cOre;
405
406 EntityHandle conn[] = {father, mother};
407
408 double coords[6];
409 if (th == NULL) {
410 CHKERR m_field.get_moab().get_coords(conn, 2, coords);
411 } else {
412 CHKERR m_field.get_moab().tag_get_data(th, conn, 2, coords);
413 }
414
415 FTensor::Index<'i', 3> i;
416 FTensor::Tensor1<double, 3> t_coords(coords[0], coords[1], coords[2]);
418 for (int nn = 0; nn != 3; nn++) {
419 t_delta(nn) = coords[3 + nn] - t_coords(nn);
420 }
421
422 t_move(i) = t_coords(i);
423 double min_quality_i = 1;
424 CHKERR minQuality(check_tests, father, mother, &t_move(0), min_quality_i, th,
426
427 t_move(i) = t_coords(i) + t_delta(i);
428 double min_quality_k = 1;
429 CHKERR minQuality(check_tests, father, mother, &t_move(0), min_quality_k, th,
431
432 double alpha_i = 0;
433 double alpha_k = 1;
434
435 for (int ii = 0; ii != line_search; ii++) {
436 double min_quality = 1;
437 double alpha = (alpha_i + alpha_k) * 0.5;
438 t_move(i) = t_coords(i) + alpha * t_delta(i);
439 CHKERR minQuality(check_tests, father, mother, &t_move(0), min_quality, th,
441 if (min_quality_i >= min_quality_k) {
442 min_quality_k = min_quality;
443 alpha_k = alpha;
444 } else {
445 min_quality_i = min_quality;
446 alpha_i = alpha;
447 }
448 }
449
450 if (min_quality_i > min_quality_k) {
451 t_move(i) = t_coords(i) + alpha_i * t_delta(i);
452 } else {
453 t_move(i) = t_coords(i) + alpha_k * t_delta(i);
454 }
455
457}
458
460 EntityHandle father, EntityHandle mother, BitRefLevel bit, Range *tets_ptr,
461 const bool only_if_improve_quality, const double move, Tag th) {
462 Interface &m_field = cOre;
464 Range out_tets;
465 CHKERR mergeNodes(father, mother, out_tets, tets_ptr, only_if_improve_quality,
466 move, 0, th);
467 CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevel(out_tets, bit);
469}
470
473 BitRefLevel tets_from_bit_ref_level, const bool only_if_improve_quality,
474 const double move, Tag th) {
475 Interface &m_field = cOre;
477 Range level_tets;
478 CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
479 tets_from_bit_ref_level, BitRefLevel().set(), MBTET, level_tets);
480 CHKERR mergeNodes(father, mother, bit, &level_tets, only_if_improve_quality,
481 move, th);
483}
484
485} // namespace MoFEM
constexpr double a
@ VERY_VERBOSE
#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
Definition definitions.h:31
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
#define CHKERR
Inline error check.
#define THROW_MESSAGE(msg)
Throw MoFEM exception.
auto bit
set bit
FTensor::Index< 'i', SPACE_DIM > i
const double n
refractive index of diffusive medium
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition Types.hpp:40
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
static auto min_non_abs(const double a, const double b)
Managing BitRefLevels.
virtual moab::Interface & get_moab()=0
virtual MPI_Comm & get_comm() const =0
Core (interface) class.
Definition Core.hpp:82
Tag get_th_RefParentHandle() const
Definition Core.hpp:197
Deprecated interface functions.
Merge node by collapsing edge between them.
ParentChildMap parentChildMap
MoFEMErrorCode minQuality(Range &check_tests, EntityHandle father, EntityHandle mother, double *coords_move, double &min_quality, Tag th=NULL, boost::function< double(double, double)> f=[](double a, double b) -> double { return std::min(a, b);})
Calualte quality if nodes merged.
bool successMerge
True if marge is success.
boost::function< double(const double a, const double b)> minQualityFunction
NodeMergerInterface(const MoFEM::Core &core)
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
multi_index_container< FaceMap, indexed_by< hashed_unique< composite_key< FaceMap, member< FaceMap, EntityHandle, &FaceMap::n0 >, member< FaceMap, EntityHandle, &FaceMap::n1 > > > > > FaceMapIdx
bool errorIfNoCommonEdge
Send error if no common edge.
MoFEMErrorCode getSubInterfaceOptions()
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
MoFEMErrorCode lineSearch(Range &check_tests, EntityHandle father, EntityHandle mother, int line_search, FTensor::Tensor1< double, 3 > &t_move, Tag th=NULL)
Use bisection method to find point of edge collapse.
static double volumeLengthQuality(const double *coords)
Calculate tetrahedron volume length quality.
Definition Tools.cpp:15
base class for all interface classes
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.