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