v0.15.0
Loading...
Searching...
No Matches
MeshRefinement.cpp
Go to the documentation of this file.
1/** \file MeshRefinement.cpp
2 * \brief Mesh refinement interface
3 */
4
5
6
7#include <EntityRefine.hpp>
8
9namespace MoFEM {
10
11typedef multi_index_container<
12 boost::shared_ptr<RefElement>,
13 // ptrWrapperRefElement,
14 indexed_by<
15 ordered_unique<tag<Ent_mi_tag>,
16 const_mem_fun<RefElement::interface_type_RefEntity,
17 EntityHandle, &RefElement::getEnt>>,
18 ordered_non_unique<
19 tag<Ent_Ent_mi_tag>,
20 const_mem_fun<RefElement::interface_type_RefEntity, EntityHandle,
21 &RefElement::getParentEnt>>,
22 ordered_non_unique<
23 tag<Composite_ParentEnt_And_BitsOfRefinedEdges_mi_tag>,
24 composite_key<
25 RefElement,
26 const_mem_fun<RefElement::interface_type_RefEntity,
27 EntityHandle, &RefElement::getParentEnt>,
28 const_mem_fun<RefElement, int,
29 &RefElement::getBitRefEdgesUlong>>>>>
30 RefElement_multiIndex_parents_view;
31
32MoFEMErrorCode
33MeshRefinement::query_interface(boost::typeindex::type_index type_index,
34 UnknownInterface **iface) const {
35 *iface = const_cast<MeshRefinement *>(this);
36 return 0;
37}
38
39MeshRefinement::MeshRefinement(const Core &core)
40 : cOre(const_cast<Core &>(core)) {}
41
42MoFEMErrorCode MeshRefinement::addVerticesInTheMiddleOfEdges(
43 const EntityHandle meshset, const BitRefLevel &bit, const bool recursive,
44 int verb, EntityHandle start_v) {
45 Interface &m_field = cOre;
46 moab::Interface &moab = m_field.get_moab();
47 MoFEMFunctionBegin;
48 Range edges;
49 CHKERR moab.get_entities_by_type(meshset, MBEDGE, edges, recursive);
50 if (edges.empty()) {
51 Range tets;
52 CHKERR moab.get_entities_by_type(meshset, MBTET, tets, recursive);
53 CHKERR moab.get_adjacencies(tets, 1, true, edges, moab::Interface::UNION);
54 if (tets.empty()) {
55 Range prisms;
56 CHKERR moab.get_entities_by_type(meshset, MBPRISM, prisms, recursive);
57 for (Range::iterator pit = prisms.begin(); pit != prisms.end(); pit++) {
58 const EntityHandle *conn;
59 int num_nodes;
60 CHKERR moab.get_connectivity(*pit, conn, num_nodes, true);
61 assert(num_nodes == 6);
62 //
63 Range edge;
64 CHKERR moab.get_adjacencies(&conn[0], 2, 1, true, edge);
65 assert(edge.size() == 1);
66 edges.insert(edge[0]);
67 edge.clear();
68 CHKERR moab.get_adjacencies(&conn[1], 2, 1, true, edge);
69 assert(edge.size() == 1);
70 edges.insert(edge[0]);
71 EntityHandle conn_edge2[] = {conn[2], conn[0]};
72 edge.clear();
73 CHKERR moab.get_adjacencies(conn_edge2, 2, 1, true, edge);
74 assert(edge.size() == 1);
75 edges.insert(edge[0]);
76 //
77 edge.clear();
78 CHKERR moab.get_adjacencies(&conn[3], 2, 1, true, edge);
79 assert(edge.size() == 1);
80 edges.insert(edge[0]);
81 edge.clear();
82 CHKERR moab.get_adjacencies(&conn[4], 2, 1, true, edge);
83 assert(edge.size() == 1);
84 edges.insert(edge[0]);
85 EntityHandle conn_edge8[] = {conn[5], conn[3]};
86 edge.clear();
87 CHKERR moab.get_adjacencies(conn_edge8, 2, 1, true, edge);
88 assert(edge.size() == 1);
89 edges.insert(edge[0]);
90 }
91 }
92 }
93 CHKERR addVerticesInTheMiddleOfEdges(edges, bit, verb, start_v);
94 MoFEMFunctionReturn(0);
95}
96MoFEMErrorCode MeshRefinement::addVerticesInTheMiddleOfEdges(
97 const Range &ents, const BitRefLevel &bit, int verb, EntityHandle start_v) {
98 Interface &m_field = cOre;
99 moab::Interface &moab = m_field.get_moab();
100 auto refined_ents_ptr = m_field.get_ref_ents();
101 MoFEMFunctionBegin;
102 auto miit =
103 refined_ents_ptr->get<Composite_EntType_and_ParentEntType_mi_tag>()
104 .lower_bound(boost::make_tuple(MBVERTEX, MBEDGE));
105 auto hi_miit =
106 refined_ents_ptr->get<Composite_EntType_and_ParentEntType_mi_tag>()
107 .upper_bound(boost::make_tuple(MBVERTEX, MBEDGE));
108 RefEntity_multiIndex_view_by_ordered_parent_entity ref_parent_ents_view;
109 ref_parent_ents_view.insert(miit, hi_miit);
110
111 Range edges = ents.subset_by_type(MBEDGE);
112 if (verb >= VERBOSE) {
113 MOFEM_TAG_AND_LOG("WORLD", Sev::verbose, "MeshRefinemnt")
114 << "ref level " << bit << " nb. edges to refine " << edges.size();
115 }
116
117 std::array<std::vector<double>, 3> vert_coords;
118 for (auto &vc : vert_coords)
119 vc.reserve(edges.size());
120
121 std::vector<EntityHandle> parent_edge;
122 parent_edge.reserve(edges.size());
123
124 std::array<double, 9> coords;
125 FTensor::Tensor1<FTensor::PackPtr<double *, 0>, 3> t_0 = {
126 &coords[0], &coords[1], &coords[2]};
127 FTensor::Tensor1<FTensor::PackPtr<double *, 0>, 3> t_1 = {
128 &coords[3], &coords[4], &coords[5]};
129 FTensor::Tensor1<FTensor::PackPtr<double *, 0>, 3> t_2 = {
130 &coords[6], &coords[7], &coords[8]};
131
132 FTensor::Index<'i', 3> i;
133
134 Range add_bit;
135 for (auto p_eit = edges.pair_begin(); p_eit != edges.pair_end(); ++p_eit) {
136 auto miit_view = ref_parent_ents_view.lower_bound(p_eit->first);
137
138 Range edge_having_parent_vertex;
139 if (miit_view != ref_parent_ents_view.end()) {
140 for (auto hi_miit_view = ref_parent_ents_view.upper_bound(p_eit->second);
141 miit_view != hi_miit_view; ++miit_view) {
142 edge_having_parent_vertex.insert(edge_having_parent_vertex.end(),
143 miit_view->get()->getParentEnt());
144 add_bit.insert(add_bit.end(), miit_view->get()->getEnt());
145 }
146 }
147
148 Range add_vertex_edges =
149 subtract(Range(p_eit->first, p_eit->second), edge_having_parent_vertex);
150
151 for (auto e : add_vertex_edges)
152 parent_edge.emplace_back(e);
153
154 for (auto e : add_vertex_edges) {
155
156 const EntityHandle *conn;
157 int num_nodes;
158 CHKERR moab.get_connectivity(e, conn, num_nodes, false);
159 if(num_nodes == 2) {
160 CHKERR moab.get_coords(conn, num_nodes, coords.data());
161 t_0(i) += t_1(i);
162 t_0(i) *= 0.5;
163 for (auto j : {0, 1, 2})
164 vert_coords[j].emplace_back(t_0(j));
165
166 } else if (num_nodes == 3) {
167 for (auto j : {0, 1, 2})
168 vert_coords[j].emplace_back(t_2(j));
169 } else {
170 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
171 "edge should have 2 or 3 nodes but have %d", num_nodes);
172 }
173 }
174 }
175
176 CHKERR m_field.getInterface<BitRefManager>()->addBitRefLevel(add_bit, bit);
177
178 if (!vert_coords[0].empty()) {
179 ReadUtilIface *read_util;
180 CHKERR moab.query_interface(read_util);
181 int num_nodes = vert_coords[0].size();
182 vector<double *> arrays_coord;
183 CHKERR read_util->get_node_coords(3, num_nodes, 0, start_v, arrays_coord);
184 Range verts(start_v, start_v + num_nodes - 1);
185 for (auto dd : {0, 1, 2}) {
186 std::copy(vert_coords[dd].begin(), vert_coords[dd].end(),
187 arrays_coord[dd]);
188 }
189 CHKERR moab.tag_set_data(cOre.get_th_RefParentHandle(), verts,
190 &*parent_edge.begin());
191 CHKERR m_field.getInterface<BitRefManager>()->setEntitiesBitRefLevel(
192 verts, bit, verb);
193 }
194 MoFEMFunctionReturn(0);
195}
196
197MoFEMErrorCode MeshRefinement::refineTets(const EntityHandle meshset,
198 const BitRefLevel &bit, int verb,
199 const bool debug) {
200 Interface &m_field = cOre;
201 moab::Interface &moab = m_field.get_moab();
202 MoFEMFunctionBegin;
203 Range tets;
204 CHKERR moab.get_entities_by_type(meshset, MBTET, tets, false);
205 CHKERR refineTets(tets, bit, verb, debug);
206 MoFEMFunctionReturn(0);
207}
208
209MoFEMErrorCode MeshRefinement::SetParent::operator()(
210 const EntityHandle ent, const EntityHandle parent,
211 const RefEntity_multiIndex *ref_ents_ptr, MoFEM::Core &cOre) {
212 MoFEM::Interface &m_field = cOre;
213 MoFEMFunctionBegin;
214 auto it = ref_ents_ptr->find(ent);
215 if (it != ref_ents_ptr->end()) {
216 if (it->get()->getParentEnt() != parent && ent != parent) {
217 parentsToChange[ent] = parent;
218 }
219 } else {
220 if (ent != parent) {
221 CHKERR m_field.get_moab().tag_set_data(cOre.get_th_RefParentHandle(),
222 &ent, 1, &parent);
223 }
224 }
225 MoFEMFunctionReturn(0);
226}
227
228MoFEMErrorCode MeshRefinement::SetParent::operator()(
229 const RefEntity_multiIndex *ref_ents_ptr) {
230 MoFEMFunctionBegin;
231 for (auto mit = parentsToChange.begin(); mit != parentsToChange.end();
232 ++mit) {
233 auto it = ref_ents_ptr->find(mit->first);
234 bool success = const_cast<RefEntity_multiIndex *>(ref_ents_ptr)
235 ->modify(it, RefEntity_change_parent(mit->second));
236 if (!success) {
237 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
238 "impossible to set parent");
239 }
240 }
241 MoFEMFunctionReturn(0);
242}
243
244MoFEMErrorCode MeshRefinement::refineTets(const Range &_tets,
245 const BitRefLevel &bit, int verb,
246 const bool debug) {
247
248 MoFEMFunctionBegin;
249
250 auto set_edge_bits = [](moab::Interface &moab,
251 RefEntity_multiIndex_view_by_ordered_parent_entity
252 &ref_parent_ents_view,
253 EntityHandle tet, BitRefEdges &parent_edges_bit,
254 EntityHandle *edge_new_nodes, int *split_edges) {
255 MoFEMFunctionBeginHot;
256
257 for (int ee = 0; ee != 6; ++ee) {
258 EntityHandle edge;
259 CHKERR moab.side_element(tet, 1, ee, edge);
260 auto eit = ref_parent_ents_view.find(edge);
261 if (eit != ref_parent_ents_view.end()) {
262 edge_new_nodes[ee] = (*eit)->getEnt();
263 split_edges[parent_edges_bit.count()] = ee;
264 parent_edges_bit.set(ee, 1);
265 }
266 }
267
268 MoFEMFunctionReturnHot(0);
269 };
270
271 CHKERR refineTets(_tets, bit, set_edge_bits, verb, debug);
272
273 MoFEMFunctionReturn(0);
274}
275
276MoFEMErrorCode MeshRefinement::refineTetsHangingNodes(const Range &_tets,
277 const BitRefLevel &bit,
278 int verb,
279 const bool debug) {
280
281 MoFEMFunctionBegin;
282
283 auto set_edge_bits = [](moab::Interface &moab,
284 RefEntity_multiIndex_view_by_ordered_parent_entity
285 &ref_parent_ents_view,
286 EntityHandle tet, BitRefEdges &parent_edges_bit,
287 EntityHandle *edge_new_nodes, int *split_edges) {
288 MoFEMFunctionBeginHot;
289
290 for (int ee = 0; ee != 6; ++ee) {
291 EntityHandle edge;
292 CHKERR moab.side_element(tet, 1, ee, edge);
293 auto eit = ref_parent_ents_view.find(edge);
294 if (eit != ref_parent_ents_view.end()) {
295 edge_new_nodes[ee] = (*eit)->getEnt();
296 split_edges[parent_edges_bit.count()] = ee;
297 parent_edges_bit.set(ee, 1);
298 }
299 }
300
301 if(parent_edges_bit.count() < 6)
302 parent_edges_bit.reset();
303
304 MoFEMFunctionReturnHot(0);
305 };
306
307 CHKERR refineTets(_tets, bit, set_edge_bits, verb, debug);
308
309 MoFEMFunctionReturn(0);
310}
311
312MoFEMErrorCode
313MeshRefinement::refineTetsHangingNodes(const EntityHandle meshset,
314 const BitRefLevel &bit, int verb,
315 const bool debug) {
316 Interface &m_field = cOre;
317 moab::Interface &moab = m_field.get_moab();
318 MoFEMFunctionBegin;
319 Range tets;
320 CHKERR moab.get_entities_by_type(meshset, MBTET, tets, false);
321 CHKERR refineTetsHangingNodes(tets, bit, verb, debug);
322 MoFEMFunctionReturn(0);
323}
324
325MoFEMErrorCode MeshRefinement::refineTets(const Range &_tets,
326 const BitRefLevel &bit,
327 SetEdgeBitsFun set_edge_bits,
328 int verb, const bool debug) {
329
330 Interface &m_field = cOre;
331 moab::Interface &moab = m_field.get_moab();
332 auto refined_ents_ptr = m_field.get_ref_ents();
333 auto refined_finite_elements_ptr = m_field.get_ref_finite_elements();
334 ReadUtilIface *read_util;
335 MoFEMFunctionBegin;
336
337 CHKERR m_field.get_moab().query_interface(read_util);
338
339 Range ents_to_set_bit;
340
341 // FIXME: refinement is based on entity handlers, should work on global ids of
342 // nodes, this will allow parallelise algorithm in the future
343
344 auto get_parent_ents_view = [&](const auto type_child,
345 const auto type_parent) {
346 auto &ref_ents =
347 refined_ents_ptr->get<Composite_EntType_and_ParentEntType_mi_tag>();
348 auto range =
349 ref_ents.equal_range(boost::make_tuple(type_child, type_parent));
350
351 RefEntity_multiIndex_view_by_ordered_parent_entity ref_parent_ents_view;
352
353 using I = decltype(range.first);
354
355 boost::function<bool(I it)> tester = [&](I it) {
356 return ((*it)->getBitRefLevel() & bit).any();
357 };
358
359 boost::function<MoFEMErrorCode(I f, I s)> inserter = [&](I f, I s) {
360 ref_parent_ents_view.insert(f, s);
361 return 0;
362 };
363
364 rangeInserter(range.first, range.second, tester, inserter);
365
366 return ref_parent_ents_view;
367 };
368
369 auto ref_parent_ents_view = get_parent_ents_view(MBVERTEX, MBEDGE);
370
371 auto tets = _tets.subset_by_type(MBTET);
372
373 auto get_ele_parent_view = [&]() {
374 auto &ref_ents = refined_finite_elements_ptr->get<Ent_mi_tag>();
375 RefElement_multiIndex_parents_view ref_ele_parent_view;
376 for (auto p = tets.pair_begin(); p != tets.pair_end(); ++p) {
377 auto tmi = ref_ents.lower_bound(p->first);
378 auto tme = ref_ents.upper_bound(p->second);
379 ref_ele_parent_view.insert(tmi, tme);
380 }
381 return ref_ele_parent_view;
382 };
383
384 auto ref_ele_parent_view = get_ele_parent_view();
385
386 std::vector<EntityHandle> parent_tets_refinded;
387 std::vector<EntityHandle> vertices_of_created_tets;
388 std::vector<BitRefEdges> parent_edges_bit_vec;
389 std::vector<int> nb_new_tets_vec;
390 std::vector<int> sub_type_vec;
391
392 parent_tets_refinded.reserve(tets.size());
393 vertices_of_created_tets.reserve(4 * tets.size());
394 parent_edges_bit_vec.reserve(tets.size());
395 nb_new_tets_vec.reserve(tets.size());
396 sub_type_vec.reserve(tets.size());
397
398 // make loop over all tets which going to be refined
399 for (auto tit : tets) {
400
401 // get tet connectivity
402 const EntityHandle *conn;
403 int num_nodes;
404 CHKERR moab.get_connectivity(tit, conn, num_nodes, true);
405
406 // get edges
407 BitRefEdges parent_edges_bit(0);
408 EntityHandle edge_new_nodes[] = {0, 0, 0, 0, 0, 0};
409 int split_edges[] = {-1, -1, -1, -1, -1, -1};
410
411 CHKERR set_edge_bits(moab, ref_parent_ents_view, tit, parent_edges_bit,
412 edge_new_nodes, split_edges);
413
414 if (parent_edges_bit.count()) {
415
416 // swap nodes forward
417 EntityHandle _conn_[4];
418 std::copy(&conn[0], &conn[4], &_conn_[0]);
419
420 // build connectivity for ref tets
421 EntityHandle new_tets_conns[8 * 4];
422 std::fill(&new_tets_conns[0], &new_tets_conns[8 * 4], no_handle);
423
424 int sub_type = -1, nb_new_tets = 0;
425 switch (parent_edges_bit.count()) {
426 case 1:
427 sub_type = 0;
428 tet_type_1(_conn_, split_edges[0], edge_new_nodes[split_edges[0]],
429 new_tets_conns);
430 nb_new_tets = 2;
431 break;
432 case 2:
433 sub_type =
434 tet_type_2(_conn_, split_edges, edge_new_nodes, new_tets_conns);
435 if (sub_type & (4 | 8 | 16)) {
436 nb_new_tets = 3;
437 break;
438 } else if (sub_type == 1) {
439 nb_new_tets = 4;
440 break;
441 };
442 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "Impossible case");
443 break;
444 case 3:
445 sub_type =
446 tet_type_3(_conn_, split_edges, edge_new_nodes, new_tets_conns);
447 if (sub_type <= 4) {
448 nb_new_tets = 4;
449 break;
450 } else if (sub_type <= 7) {
451 nb_new_tets = 5;
452 break;
453 }
454 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "Impossible case");
455 case 4:
456 sub_type =
457 tet_type_4(_conn_, split_edges, edge_new_nodes, new_tets_conns);
458 if (sub_type == 0) {
459 nb_new_tets = 5;
460 break;
461 } else if (sub_type <= 7) {
462 nb_new_tets = 6;
463 break;
464 }
465 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "Impossible case");
466 case 5:
467 sub_type = tet_type_5(moab, _conn_, edge_new_nodes, new_tets_conns);
468 nb_new_tets = 7;
469 break;
470 case 6:
471 sub_type = 0;
472 tet_type_6(moab, _conn_, edge_new_nodes, new_tets_conns);
473 nb_new_tets = 8;
474 break;
475 default:
476 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "Impossible case");
477 }
478
479 // find that tets
480 auto &ref_ele_by_parent_and_ref_edges =
481 ref_ele_parent_view
482 .get<Composite_ParentEnt_And_BitsOfRefinedEdges_mi_tag>();
483 auto it_by_ref_edges = ref_ele_by_parent_and_ref_edges.equal_range(
484 boost::make_tuple(tit, parent_edges_bit.to_ulong()));
485 // check if tet with this refinement scheme already exits
486 std::vector<EntityHandle> ents_to_set_bit_vec;
487 if (std::distance(it_by_ref_edges.first, it_by_ref_edges.second) ==
488 static_cast<size_t>(nb_new_tets)) {
489 ents_to_set_bit_vec.reserve(nb_new_tets);
490 for (int tt = 0; it_by_ref_edges.first != it_by_ref_edges.second;
491 it_by_ref_edges.first++, tt++) {
492 auto tet = (*it_by_ref_edges.first)->getEnt();
493 ents_to_set_bit_vec.emplace_back(tet);
494 // set ref tets entities
495 if (debug) {
496 // add this tet if exist to this ref level
497 auto ref_tet_it = refined_ents_ptr->find(tet);
498 if (ref_tet_it == refined_ents_ptr->end()) {
499 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
500 "Tetrahedron should be in database");
501 }
502 }
503 }
504 ents_to_set_bit.insert_list(ents_to_set_bit_vec.begin(),
505 ents_to_set_bit_vec.end());
506
507 } else {
508 // if this element was not refined or was refined with different
509 // patterns of split edges create new elements
510
511 parent_tets_refinded.emplace_back(tit);
512 for (int tt = 0; tt != nb_new_tets; ++tt) {
513 for (auto nn : {0, 1, 2, 3})
514 vertices_of_created_tets.emplace_back(new_tets_conns[4 * tt + nn]);
515 }
516 parent_edges_bit_vec.emplace_back(parent_edges_bit);
517 nb_new_tets_vec.emplace_back(nb_new_tets);
518 sub_type_vec.emplace_back(sub_type);
519 }
520 } else {
521
522 if (debug) {
523 RefEntity_multiIndex::iterator tit_miit;
524 tit_miit = refined_ents_ptr->find(tit);
525 if (tit_miit == refined_ents_ptr->end())
526 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
527 "can not find this tet");
528 }
529
530 ents_to_set_bit.insert(tit);
531 }
532 }
533
534 // Create tets
535 EntityHandle start_e = 0;
536 EntityHandle *conn_moab;
537 const int nb_new_tets = vertices_of_created_tets.size() / 4;
538 read_util->get_element_connect(nb_new_tets, 4, MBTET, 0, start_e, conn_moab);
539 std::copy(vertices_of_created_tets.begin(), vertices_of_created_tets.end(),
540 conn_moab);
541 CHKERR read_util->update_adjacencies(start_e, nb_new_tets, 4, conn_moab);
542 ents_to_set_bit.insert(start_e, start_e + nb_new_tets - 1);
543
544 // Create adj entities
545 Range ents;
546 for (auto d : {1, 2}) {
547 CHKERR moab.get_adjacencies(ents_to_set_bit, d, true, ents,
548 moab::Interface::UNION);
549 }
550
551 SetParent set_parent;
552
553 // Set parents and adjacencies
554 for (int idx = 0; idx != parent_tets_refinded.size(); ++idx) {
555
556 const EntityHandle tit = parent_tets_refinded[idx];
557 const BitRefEdges &parent_edges_bit = parent_edges_bit_vec[idx];
558 const int nb_new_tets = nb_new_tets_vec[idx];
559
560 std::array<EntityHandle, 8> ref_tets;
561 for (int tt = 0; tt != nb_new_tets; ++tt, ++start_e)
562 ref_tets[tt] = start_e;
563
564 if (nb_new_tets) {
565
566 CHKERR moab.tag_clear_data(cOre.get_th_RefBitEdge(), &ref_tets[0],
567 nb_new_tets, &parent_edges_bit);
568 CHKERR moab.tag_clear_data(cOre.get_th_RefParentHandle(), &ref_tets[0],
569 nb_new_tets, &tit);
570
571 // hash map of nodes (RefEntity) by edges (EntityHandle)
572 std::map<EntityHandle /*edge*/, EntityHandle /*node*/>
573 map_ref_nodes_by_edges;
574
575 Range tit_edges;
576 CHKERR moab.get_adjacencies(&tit, 1, 1, false, tit_edges);
577 for (auto edge : tit_edges) {
578 auto eit = ref_parent_ents_view.find(edge);
579 if (eit != ref_parent_ents_view.end()) {
580 map_ref_nodes_by_edges[(*eit)->getParentEnt()] = eit->get()->getEnt();
581 }
582 }
583
584 // find parents for new edges and faces
585 // get tet edges and faces
586 Range tit_faces;
587 CHKERR moab.get_adjacencies(&tit, 1, 2, false, tit_faces);
588 Range edges_nodes[6], faces_nodes[4];
589 // for edges - add ref nodes
590 // edges_nodes[ee] - contains all nodes on edge ee including mid nodes if
591 // exist
592 Range::iterator eit = tit_edges.begin();
593 for (int ee = 0; eit != tit_edges.end(); eit++, ee++) {
594 CHKERR moab.get_connectivity(&*eit, 1, edges_nodes[ee], true);
595 auto map_miit = map_ref_nodes_by_edges.find(*eit);
596 if (map_miit != map_ref_nodes_by_edges.end()) {
597 edges_nodes[ee].insert(map_miit->second);
598 }
599 }
600 // for faces - add ref nodes
601 // faces_nodes[ff] - contains all nodes on face ff including mid nodes if
602 // exist
603 Range::iterator fit = tit_faces.begin();
604 for (int ff = 0; fit != tit_faces.end(); fit++, ff++) {
605 CHKERR moab.get_connectivity(&*fit, 1, faces_nodes[ff], true);
606 // Get edges on face and loop over those edges to add mid-nodes to range
607 Range fit_edges;
608 CHKERR moab.get_adjacencies(&*fit, 1, 1, false, fit_edges);
609 for (Range::iterator eit2 = fit_edges.begin(); eit2 != fit_edges.end();
610 eit2++) {
611 auto map_miit = map_ref_nodes_by_edges.find(*eit2);
612 if (map_miit != map_ref_nodes_by_edges.end()) {
613 faces_nodes[ff].insert(map_miit->second);
614 }
615 }
616 }
617 // add ref nodes to tet
618 // tet_nodes contains all nodes on tet including mid edge nodes
619 Range tet_nodes;
620 CHKERR moab.get_connectivity(&tit, 1, tet_nodes, true);
621 for (std::map<EntityHandle, EntityHandle>::iterator map_miit =
622 map_ref_nodes_by_edges.begin();
623 map_miit != map_ref_nodes_by_edges.end(); map_miit++) {
624 tet_nodes.insert(map_miit->second);
625 }
626 Range ref_edges;
627 // Get all all edges of refined tets
628 CHKERR moab.get_adjacencies(ref_tets.data(), nb_new_tets, 1, false,
629 ref_edges, moab::Interface::UNION);
630 // Check for all ref edge and set parents
631 for (Range::iterator reit = ref_edges.begin(); reit != ref_edges.end();
632 reit++) {
633 Range ref_edges_nodes;
634 CHKERR moab.get_connectivity(&*reit, 1, ref_edges_nodes, true);
635
636 // Check if ref edge is an coarse edge (loop over coarse tet edges)
637 int ee = 0;
638 for (; ee < 6; ee++) {
639 // Two nodes are common (node[0],node[1],ref_node (if exist))
640 // this tests if given edge is contained by edge of refined
641 // tetrahedral
642 if (intersect(edges_nodes[ee], ref_edges_nodes).size() == 2) {
643 EntityHandle edge = tit_edges[ee];
644 CHKERR set_parent(*reit, edge, refined_ents_ptr, cOre);
645 break;
646 }
647 }
648 if (ee < 6)
649 continue; // this refined edge is contained by edge of tetrahedral
650 // check if ref edge is in coarse face
651 int ff = 0;
652 for (; ff < 4; ff++) {
653 // two nodes are common (node[0],node[1],ref_node (if exist))
654 // this tests if given face is contained by face of tetrahedral
655 if (intersect(faces_nodes[ff], ref_edges_nodes).size() == 2) {
656 EntityHandle face = tit_faces[ff];
657 CHKERR set_parent(*reit, face, refined_ents_ptr, cOre);
658 break;
659 }
660 }
661 if (ff < 4)
662 continue; // this refined edge is contained by face of tetrahedral
663
664 // check if ref edge is in coarse tetrahedral (i.e. that is internal
665 // edge of refined tetrahedral)
666 if (intersect(tet_nodes, ref_edges_nodes).size() == 2) {
667 CHKERR set_parent(*reit, tit, refined_ents_ptr, cOre);
668 continue;
669 }
670
671 // Refined edge is not child of any edge, face or tetrahedral, this is
672 // Impossible edge
673 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
674 "impossible refined edge");
675 }
676
677 Range ref_faces;
678 CHKERR moab.get_adjacencies(ref_tets.data(), nb_new_tets, 2, false,
679 ref_faces, moab::Interface::UNION);
680 Tag th_interface_side;
681 const int def_side[] = {0};
682 CHKERR moab.tag_get_handle("INTERFACE_SIDE", 1, MB_TYPE_INTEGER,
683 th_interface_side,
684 MB_TAG_CREAT | MB_TAG_SPARSE, def_side);
685 // Check for all ref faces
686 for (auto rfit : ref_faces) {
687 Range ref_faces_nodes;
688 CHKERR moab.get_connectivity(&rfit, 1, ref_faces_nodes, true);
689 // Check if ref face is in coarse face
690 int ff = 0;
691 for (; ff < 4; ff++) {
692 // Check if refined triangle is contained by face of tetrahedral
693 if (intersect(faces_nodes[ff], ref_faces_nodes).size() == 3) {
694 EntityHandle face = tit_faces[ff];
695 CHKERR set_parent(rfit, face, refined_ents_ptr, cOre);
696 int side = 0;
697 // Set face side if it is on interface
698 CHKERR moab.tag_get_data(th_interface_side, &face, 1, &side);
699 CHKERR moab.tag_set_data(th_interface_side, &rfit, 1, &side);
700 break;
701 }
702 }
703 if (ff < 4)
704 continue; // this face is contained by one of tetrahedrons
705 // check if ref face is in coarse tetrahedral
706 // this is ref face which is contained by tetrahedral volume
707 if (intersect(tet_nodes, ref_faces_nodes).size() == 3) {
708 CHKERR set_parent(rfit, tit, refined_ents_ptr, cOre);
709 continue;
710 }
711 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
712 "impossible refined face");
713 }
714 }
715 }
716
717 CHKERR set_parent(refined_ents_ptr);
718 CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevel(
719 ents_to_set_bit, bit, true, verb, &ents);
720 MoFEMFunctionReturn(0);
721}
722MoFEMErrorCode MeshRefinement::refinePrisms(const EntityHandle meshset,
723 const BitRefLevel &bit, int verb) {
724
725 Interface &m_field = cOre;
726 moab::Interface &moab = m_field.get_moab();
727 auto refined_ents_ptr = m_field.get_ref_ents();
728 auto refined_finite_elements_ptr = m_field.get_ref_finite_elements();
729
730 // FIXME: refinement is based on entity handlers, should work on global ids of
731 // nodes, this will allow parallelise algorithm in the future
732
733 MoFEMFunctionBegin;
734
735 RefElement_multiIndex_parents_view ref_ele_parent_view;
736 ref_ele_parent_view.insert(
737 refined_finite_elements_ptr->get<Ent_mi_tag>().lower_bound(
738 get_id_for_min_type<MBPRISM>()),
739 refined_finite_elements_ptr->get<Ent_mi_tag>().upper_bound(
740 get_id_for_max_type<MBPRISM>()));
741 auto &ref_ele_by_parent_and_ref_edges =
742 ref_ele_parent_view
743 .get<Composite_ParentEnt_And_BitsOfRefinedEdges_mi_tag>();
744 // find all vertices which parent is edge
745 auto &ref_ents_by_comp =
746 refined_ents_ptr->get<Composite_EntType_and_ParentEntType_mi_tag>();
747 RefEntity_multiIndex_view_by_hashed_parent_entity ref_parent_ents_view;
748 ref_parent_ents_view.insert(
749 ref_ents_by_comp.lower_bound(boost::make_tuple(MBVERTEX, MBEDGE)),
750 ref_ents_by_comp.upper_bound(boost::make_tuple(MBVERTEX, MBEDGE)));
751 Range prisms;
752 CHKERR moab.get_entities_by_type(meshset, MBPRISM, prisms, false);
753 Range::iterator pit = prisms.begin();
754 for (; pit != prisms.end(); pit++) {
755 auto miit_prism = refined_ents_ptr->get<Ent_mi_tag>().find(*pit);
756 if (miit_prism == refined_ents_ptr->end()) {
757 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
758 "this prism is not in ref database");
759 }
760 if (verb >= NOISY) {
761 MOFEM_TAG_AND_LOG("WORLD", Sev::noisy, "MeshRefinment")
762 << "ref prism " << **miit_prism;
763 }
764 // prism connectivity
765 int num_nodes;
766 const EntityHandle *conn;
767 CHKERR moab.get_connectivity(*pit, conn, num_nodes, true);
768 assert(num_nodes == 6);
769 // edges connectivity
770 EntityHandle edges[6];
771 for (int ee = 0; ee < 3; ee++) {
772 CHKERR moab.side_element(*pit, 1, ee, edges[ee]);
773 }
774 for (int ee = 6; ee < 9; ee++) {
775 CHKERR moab.side_element(*pit, 1, ee, edges[ee - 3]);
776 }
777 // detect split edges
778 BitRefEdges split_edges(0);
779 EntityHandle edge_nodes[6];
780 std::fill(&edge_nodes[0], &edge_nodes[6], no_handle);
781 for (int ee = 0; ee < 6; ee++) {
782 RefEntity_multiIndex_view_by_hashed_parent_entity::iterator miit_view =
783 ref_parent_ents_view.find(edges[ee]);
784 if (miit_view != ref_parent_ents_view.end()) {
785 if (((*miit_view)->getBitRefLevel() & bit).any()) {
786 edge_nodes[ee] = (*miit_view)->getEnt();
787 split_edges.set(ee);
788 }
789 }
790 }
791 if (split_edges.count() == 0) {
792 *(const_cast<RefEntity *>(miit_prism->get())->getBitRefLevelPtr()) |= bit;
793 if (verb >= VERY_NOISY)
794 MOFEM_TAG_AND_LOG("WORLD", Sev::noisy, "MeshRefinment")
795 << "no refinement";
796 continue;
797 }
798 // check consistency
799 if (verb >= NOISY) {
800 std::ostringstream ss;
801 ss << "prism split edges " << split_edges << " count "
802 << split_edges.count() << std::endl;
803 PetscPrintf(m_field.get_comm(), "%s", ss.str().c_str());
804 }
805 // prism ref
806 EntityHandle new_prism_conn[4 * 6];
807 std::fill(&new_prism_conn[0], &new_prism_conn[4 * 6], no_handle);
808 int nb_new_prisms = 0;
809 switch (split_edges.count()) {
810 case 0:
811 break;
812 case 2:
813 CHKERR prism_type_1(conn, split_edges, edge_nodes, new_prism_conn);
814 nb_new_prisms = 2;
815 break;
816 case 4:
817 CHKERR prism_type_2(conn, split_edges, edge_nodes, new_prism_conn);
818 nb_new_prisms = 3;
819 break;
820 case 6:
821 CHKERR prism_type_3(conn, split_edges, edge_nodes, new_prism_conn);
822 nb_new_prisms = 4;
823 break;
824 default:
825 std::ostringstream ss;
826 ss << split_edges << " : [ " << conn[0] << " " << conn[1] << " "
827 << conn[2] << " " << conn[3] << " " << conn[4] << " " << conn[5]
828 << " ]";
829 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY, "%s",
830 ss.str().c_str());
831 }
832 // find that prism
833 std::bitset<4> ref_prism_bit(0);
834 auto it_by_ref_edges = ref_ele_by_parent_and_ref_edges.lower_bound(
835 boost::make_tuple(*pit, split_edges.to_ulong()));
836 auto hi_it_by_ref_edges = ref_ele_by_parent_and_ref_edges.upper_bound(
837 boost::make_tuple(*pit, split_edges.to_ulong()));
838 auto it_by_ref_edges2 = it_by_ref_edges;
839 for (int pp = 0; it_by_ref_edges2 != hi_it_by_ref_edges;
840 it_by_ref_edges2++, pp++) {
841 // Add this tet to this ref
842 *(const_cast<RefElement *>(it_by_ref_edges2->get())
843 ->getBitRefLevelPtr()) |= bit;
844 ref_prism_bit.set(pp, 1);
845 if (verb > 2) {
846 std::ostringstream ss;
847 ss << "is refined " << *it_by_ref_edges2 << std::endl;
848 PetscPrintf(m_field.get_comm(), "%s", ss.str().c_str());
849 }
850 }
851 if (it_by_ref_edges != hi_it_by_ref_edges) {
852 if (ref_prism_bit.count() != (unsigned int)nb_new_prisms)
853 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
854 "data inconsistency");
855 } else {
856 EntityHandle ref_prisms[4];
857 // create prism
858 for (int pp = 0; pp < nb_new_prisms; pp++) {
859 if (verb > 3) {
860 std::ostringstream ss;
861 ss << "ref prism " << ref_prism_bit << std::endl;
862 PetscPrintf(m_field.get_comm(), "%s", ss.str().c_str());
863 }
864 if (!ref_prism_bit.test(pp)) {
865 CHKERR moab.create_element(MBPRISM, &new_prism_conn[6 * pp], 6,
866 ref_prisms[pp]);
867 CHKERR moab.tag_set_data(cOre.get_th_RefParentHandle(),
868 &ref_prisms[pp], 1, &*pit);
869 CHKERR moab.tag_set_data(cOre.get_th_RefBitLevel(), &ref_prisms[pp],
870 1, &bit);
871 CHKERR moab.tag_set_data(cOre.get_th_RefBitEdge(), &ref_prisms[pp], 1,
872 &split_edges);
873 auto p_ent =
874 const_cast<RefEntity_multiIndex *>(refined_ents_ptr)
875 ->insert(boost::shared_ptr<RefEntity>(new RefEntity(
876 m_field.get_basic_entity_data_ptr(), ref_prisms[pp])));
877 std::pair<RefElement_multiIndex::iterator, bool> p_fe;
878 try {
879 p_fe =
880 const_cast<RefElement_multiIndex *>(refined_finite_elements_ptr)
881 ->insert(boost::shared_ptr<RefElement>(
882 new RefElement_PRISM(*p_ent.first)));
883 } catch (MoFEMException const &e) {
884 SETERRQ(PETSC_COMM_SELF, e.errorCode, "%s", e.errorMessage);
885 }
886 ref_prism_bit.set(pp);
887 CHKERR cOre.addPrismToDatabase(ref_prisms[pp]);
888 if (verb > 2) {
889 std::ostringstream ss;
890 ss << "add prism: " << **p_fe.first << std::endl;
891 if (verb > 7) {
892 for (int nn = 0; nn < 6; nn++) {
893 ss << new_prism_conn[nn] << " ";
894 }
895 ss << std::endl;
896 }
897 PetscPrintf(m_field.get_comm(), "%s", ss.str().c_str());
898 }
899 }
900 }
901 }
902 }
903 MoFEMFunctionReturn(0);
904}
905MoFEMErrorCode MeshRefinement::refineMeshset(const EntityHandle meshset,
906 const BitRefLevel &bit,
907 const bool recursive, int verb) {
908 Interface &m_field = cOre;
909 auto refined_ents_ptr = m_field.get_ref_ents();
910 MoFEMFunctionBegin;
911 auto miit = refined_ents_ptr->find(meshset);
912 if (miit == refined_ents_ptr->end()) {
913 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
914 "this meshset is not in ref database");
915 }
916 CHKERR m_field.getInterface<BitRefManager>()->updateMeshsetByEntitiesChildren(
917 meshset, bit, meshset, MBMAXTYPE, recursive, verb);
918 *(const_cast<RefEntity *>(miit->get())->getBitRefLevelPtr()) |= bit;
919 MoFEMFunctionReturn(0);
920}
921
922MoFEMErrorCode MeshRefinement::refineTris(const EntityHandle meshset,
923 const BitRefLevel &bit, int verb,
924 const bool debug) {
925 Interface &m_field = cOre;
926 moab::Interface &moab = m_field.get_moab();
927 MoFEMFunctionBegin;
928 Range tris;
929 CHKERR moab.get_entities_by_type(meshset, MBTRI, tris, false);
930 CHKERR refineTris(tris, bit, verb, debug);
931 MoFEMFunctionReturn(0);
932}
933
934MoFEMErrorCode MeshRefinement::refineTris(const Range &_tris,
935 const BitRefLevel &bit, int verb,
936 const bool debug) {
937 MoFEMFunctionBegin;
938
939 auto set_edge_bits = [](moab::Interface &moab,
940 RefEntity_multiIndex_view_by_ordered_parent_entity
941 &ref_parent_ents_view,
942 EntityHandle tet, BitRefEdges &parent_edges_bit,
943 EntityHandle *edge_new_nodes, int *split_edges) {
944 MoFEMFunctionBeginHot;
945
946 for (int ee = 0; ee != 3; ++ee) {
947 EntityHandle edge;
948 CHKERR moab.side_element(tet, 1, ee, edge);
949 auto eit = ref_parent_ents_view.find(edge);
950 if (eit != ref_parent_ents_view.end()) {
951 edge_new_nodes[ee] = (*eit)->getEnt();
952 split_edges[parent_edges_bit.count()] = ee;
953 parent_edges_bit.set(ee, 1);
954 }
955 }
956
957 MoFEMFunctionReturnHot(0);
958 };
959
960 CHKERR refineTris(_tris, bit, set_edge_bits, verb, debug);
961
962 MoFEMFunctionReturn(0);
963}
964
965MoFEMErrorCode MeshRefinement::refineTrisHangingNodes(const Range &_tris,
966 const BitRefLevel &bit,
967 int verb,
968 const bool debug) {
969
970 MoFEMFunctionBegin;
971
972 auto set_edge_bits = [](moab::Interface &moab,
973 RefEntity_multiIndex_view_by_ordered_parent_entity
974 &ref_parent_ents_view,
975 EntityHandle tet, BitRefEdges &parent_edges_bit,
976 EntityHandle *edge_new_nodes, int *split_edges) {
977 MoFEMFunctionBeginHot;
978
979 for (int ee = 0; ee != 3; ++ee) {
980 EntityHandle edge;
981 CHKERR moab.side_element(tet, 1, ee, edge);
982 auto eit = ref_parent_ents_view.find(edge);
983 if (eit != ref_parent_ents_view.end()) {
984 edge_new_nodes[ee] = (*eit)->getEnt();
985 split_edges[parent_edges_bit.count()] = ee;
986 parent_edges_bit.set(ee, 1);
987 }
988 }
989
990 if (parent_edges_bit.count() < 3) {
991 parent_edges_bit.reset();
992 }
993
994 MoFEMFunctionReturnHot(0);
995 };
996
997 CHKERR refineTris(_tris, bit, set_edge_bits, verb, debug);
998
999 MoFEMFunctionReturn(0);
1000};
1001
1002MoFEMErrorCode
1003MeshRefinement::refineTrisHangingNodes(const EntityHandle meshset,
1004 const BitRefLevel &bit, int verb,
1005 const bool debug) {
1006 Interface &m_field = cOre;
1007 moab::Interface &moab = m_field.get_moab();
1008 MoFEMFunctionBegin;
1009 Range tris;
1010 CHKERR moab.get_entities_by_type(meshset, MBTRI, tris, false);
1011 CHKERR refineTrisHangingNodes(tris, bit, verb, debug);
1012 MoFEMFunctionReturn(0);
1013}
1014
1015MoFEMErrorCode MeshRefinement::refineTris(const Range &_tris,
1016 const BitRefLevel &bit,
1017 SetEdgeBitsFun set_edge_bits,
1018 int verb, const bool debug) {
1019 Interface &m_field = cOre;
1020 moab::Interface &moab = m_field.get_moab();
1021 auto refined_ents_ptr = m_field.get_ref_ents();
1022 auto refined_finite_elements_ptr = m_field.get_ref_finite_elements();
1023 ReadUtilIface *read_util;
1024
1025 MoFEMFunctionBegin;
1026
1027 auto get_parent_ents_view = [&](const auto type_child,
1028 const auto type_parent) {
1029 auto &ref_ents =
1030 refined_ents_ptr->get<Composite_EntType_and_ParentEntType_mi_tag>();
1031 auto range =
1032 ref_ents.equal_range(boost::make_tuple(type_child, type_parent));
1033
1034 RefEntity_multiIndex_view_by_ordered_parent_entity ref_parent_ents_view;
1035
1036 using I = decltype(range.first);
1037
1038 boost::function<bool(I it)> tester = [&](I it) {
1039 return ((*it)->getBitRefLevel() & bit).any();
1040 };
1041
1042 boost::function<MoFEMErrorCode(I f, I s)> inserter = [&](I f, I s) {
1043 ref_parent_ents_view.insert(f, s);
1044 return 0;
1045 };
1046
1047 rangeInserter(range.first, range.second, tester, inserter);
1048
1049 return ref_parent_ents_view;
1050 };
1051
1052 auto ref_parent_ents_view = get_parent_ents_view(MBVERTEX, MBEDGE);
1053
1054 auto tris = _tris.subset_by_type(MBTRI);
1055
1056 auto get_ele_parent_view = [&]() {
1057 auto &ref_ents = refined_finite_elements_ptr->get<Ent_mi_tag>();
1058 RefElement_multiIndex_parents_view ref_ele_parent_view;
1059 for (auto p = tris.pair_begin(); p != tris.pair_end(); ++p) {
1060 auto tmi = ref_ents.lower_bound(p->first);
1061 auto tme = ref_ents.upper_bound(p->second);
1062 ref_ele_parent_view.insert(tmi, tme);
1063 }
1064 return ref_ele_parent_view;
1065 };
1066
1067 auto ref_ele_parent_view = get_ele_parent_view();
1068
1069 Range ents_to_set_bit;
1070
1071 std::vector<EntityHandle> parent_tris_refined; // entities refined
1072 std::vector<EntityHandle> vertices_of_created_tris;
1073 std::vector<BitRefEdges> parent_edges_bit_vec;
1074 std::vector<int> nb_new_tris_vec;
1075
1076 parent_tris_refined.reserve(tris.size());
1077 vertices_of_created_tris.reserve(3 * tris.size());
1078 parent_edges_bit_vec.reserve(tris.size());
1079 nb_new_tris_vec.reserve(tris.size());
1080
1081 // make loop over all tets which going to be refined
1082 for (auto tit : tris) {
1083
1084 // get tet connectivity
1085 const EntityHandle *conn;
1086 int num_nodes;
1087 CHKERR moab.get_connectivity(tit, conn, num_nodes, true);
1088
1089 // get edges
1090 BitRefEdges parent_edges_bit(0);
1091 EntityHandle edge_new_nodes[] = {0, 0, 0};
1092 int split_edges[] = {-1, -1, -1};
1093
1094 CHKERR set_edge_bits(moab, ref_parent_ents_view, tit, parent_edges_bit,
1095 edge_new_nodes, split_edges);
1096
1097 if (parent_edges_bit.count()) {
1098
1099 // swap nodes forward
1100 EntityHandle _conn_[3];
1101 std::copy(&conn[0], &conn[3], &_conn_[0]);
1102
1103 // build connectivity for ref tets
1104 EntityHandle new_tets_conns[4 * 3];
1105 std::fill(&new_tets_conns[0], &new_tets_conns[4 * 3], no_handle);
1106
1107 int nb_new_tris = 0;
1108 switch (parent_edges_bit.count()) {
1109 case 1:
1110 CHKERR tri_type_1(conn, split_edges[0], edge_new_nodes[split_edges[0]],
1111 new_tets_conns);
1112 nb_new_tris += 2;
1113 break;
1114 case 2:
1115 CHKERR tri_type_2(conn, split_edges, edge_new_nodes, new_tets_conns);
1116 nb_new_tris += 3;
1117 break;
1118 case 3:
1119 CHKERR tri_type_3(conn, edge_new_nodes, new_tets_conns);
1120 nb_new_tris += 4;
1121 break;
1122 default:
1123 SETERRQ(PETSC_COMM_SELF, MOFEM_IMPOSSIBLE_CASE, "Impossible case");
1124 }
1125
1126 // find that tets
1127 auto &ref_ele_by_parent_and_ref_edges =
1128 ref_ele_parent_view
1129 .get<Composite_ParentEnt_And_BitsOfRefinedEdges_mi_tag>();
1130 auto it_by_ref_edges = ref_ele_by_parent_and_ref_edges.equal_range(
1131 boost::make_tuple(tit, parent_edges_bit.to_ulong()));
1132 // check if tet with this refinement scheme already exits
1133 std::vector<EntityHandle> ents_to_set_bit_vec;
1134 if (std::distance(it_by_ref_edges.first, it_by_ref_edges.second) ==
1135 static_cast<size_t>(nb_new_tris)) {
1136 ents_to_set_bit_vec.reserve(nb_new_tris);
1137 for (int tt = 0; it_by_ref_edges.first != it_by_ref_edges.second;
1138 it_by_ref_edges.first++, tt++) {
1139 auto tet = (*it_by_ref_edges.first)->getEnt();
1140 ents_to_set_bit_vec.emplace_back(tet);
1141 // set ref tets entities
1142 if (debug) {
1143 // add this tet if exist to this ref level
1144 auto ref_tet_it = refined_ents_ptr->find(tet);
1145 if (ref_tet_it == refined_ents_ptr->end()) {
1146 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1147 "Tetrahedron should be in database");
1148 }
1149 }
1150 }
1151 ents_to_set_bit.insert_list(ents_to_set_bit_vec.begin(),
1152 ents_to_set_bit_vec.end());
1153
1154 } else {
1155 // if this element was not refined or was refined with different
1156 // patterns of split edges create new elements
1157
1158 parent_tris_refined.emplace_back(tit);
1159 for (int tt = 0; tt != nb_new_tris; ++tt) {
1160 for (auto nn : {0, 1, 2})
1161 vertices_of_created_tris.emplace_back(new_tets_conns[3 * tt + nn]);
1162 }
1163 parent_edges_bit_vec.emplace_back(parent_edges_bit);
1164 nb_new_tris_vec.emplace_back(nb_new_tris);
1165 }
1166 } else {
1167
1168 if (debug) {
1169 RefEntity_multiIndex::iterator tit_miit;
1170 tit_miit = refined_ents_ptr->find(tit);
1171 if (tit_miit == refined_ents_ptr->end())
1172 SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1173 "can not find this tet");
1174 }
1175
1176 ents_to_set_bit.insert(tit);
1177 }
1178 }
1179
1180 const int nb_new_tris = vertices_of_created_tris.size() / 3;
1181 EntityHandle start_e = 0;
1182
1183 if (nb_new_tris) {
1184
1185 // Create tets
1186 EntityHandle *conn_moab;
1187 CHKERR m_field.get_moab().query_interface(read_util);
1188 CHKERR read_util->get_element_connect(nb_new_tris, 3, MBTRI, 0, start_e,
1189 conn_moab);
1190 std::copy(vertices_of_created_tris.begin(), vertices_of_created_tris.end(),
1191 conn_moab);
1192 CHKERR read_util->update_adjacencies(start_e, nb_new_tris, 3, conn_moab);
1193 ents_to_set_bit.insert(start_e, start_e + nb_new_tris - 1);
1194
1195 if (debug) {
1196 auto meshset_ptr = get_temp_meshset_ptr(moab);
1197 CHKERR moab.add_entities(*meshset_ptr, ents_to_set_bit);
1198 CHKERR moab.write_file("out_refinded_debug.vtk", "VTK", "",
1199 meshset_ptr->get_ptr(), 1);
1200 }
1201 }
1202
1203 // Create adj entities
1204 Range ents;
1205 for (auto d : {1}) {
1206 CHKERR moab.get_adjacencies(ents_to_set_bit, d, true, ents,
1207 moab::Interface::UNION);
1208 }
1209
1210 SetParent set_parent;
1211
1212 // Set parrents and adjacencies. Iterate over all parent triangles.
1213 for (int idx = 0; idx != parent_tris_refined.size(); ++idx) {
1214
1215 const auto tit = parent_tris_refined[idx]; // parent triangle
1216 const auto &parent_edges_bit =
1217 parent_edges_bit_vec[idx]; // marks what edges are refined
1218 const auto nb_new_tris =
1219 nb_new_tris_vec[idx]; // number of triangles created by refinement
1220
1221 // Iterate over refined triangles
1222 std::array<EntityHandle, 4> ref_tris;
1223 for (int tt = 0; tt != nb_new_tris; ++tt, ++start_e)
1224 ref_tris[tt] = start_e;
1225
1226 if (nb_new_tris) {
1227
1228 CHKERR moab.tag_clear_data(cOre.get_th_RefBitEdge(), &ref_tris[0],
1229 nb_new_tris, &parent_edges_bit);
1230 CHKERR moab.tag_clear_data(cOre.get_th_RefParentHandle(), &ref_tris[0],
1231 nb_new_tris, &tit);
1232
1233 // Hash map of nodes (RefEntity) by edges (EntityHandle)
1234 std::map<EntityHandle /*edge*/, EntityHandle /*node*/>
1235 map_ref_nodes_by_edges;
1236
1237 // find parents for new edges and faces
1238 // get tet edges and faces
1239 Range tit_edges; // parent edges
1240 CHKERR moab.get_adjacencies(&tit, 1, 1, false, tit_edges);
1241 for (auto edge : tit_edges) {
1242 auto eit = ref_parent_ents_view.find(edge);
1243 if (eit != ref_parent_ents_view.end()) {
1244 map_ref_nodes_by_edges[(*eit)->getParentEnt()] = eit->get()->getEnt();
1245 }
1246 }
1247
1248 Range edges_nodes[3];
1249 // for edges - add ref nodes
1250 // edges_nodes[ee] - contains all nodes on edge ee including mid nodes if
1251 // exist
1252 Range::iterator eit = tit_edges.begin();
1253 for (int ee = 0; eit != tit_edges.end(); eit++, ee++) {
1254 CHKERR moab.get_connectivity(&*eit, 1, edges_nodes[ee], true);
1255 auto map_miit = map_ref_nodes_by_edges.find(*eit);
1256 if (map_miit != map_ref_nodes_by_edges.end()) {
1257 edges_nodes[ee].insert(map_miit->second);
1258 }
1259 }
1260
1261 // add ref nodes to tri
1262 // tet_nodes contains all nodes on tet including mid edge nodes
1263 Range tri_nodes;
1264 CHKERR moab.get_connectivity(&tit, 1, tri_nodes, true);
1265 for (auto map_miit = map_ref_nodes_by_edges.begin();
1266 map_miit != map_ref_nodes_by_edges.end(); map_miit++) {
1267 tri_nodes.insert(map_miit->second);
1268 }
1269
1270 Range ref_edges; // edges on all new tets
1271 // Get all all edges of refined tets
1272 CHKERR moab.get_adjacencies(ref_tris.data(), nb_new_tris, 1, false,
1273 ref_edges, moab::Interface::UNION);
1274
1275 // Check for all ref edge and set parents
1276 for (auto reit = ref_edges.begin(); reit != ref_edges.end(); ++reit) {
1277 Range ref_edges_nodes;
1278 CHKERR moab.get_connectivity(&*reit, 1, ref_edges_nodes, true);
1279 // Check if ref edge is an coarse edge (loop over coarse tet edges)
1280 int ee = 0;
1281 for (; ee < 3; ee++) {
1282 // Two nodes are common (node[0],node[1],ref_node (if exist))
1283 // this tests if given edge is contained by edge of refined
1284 // triangle
1285 if (intersect(edges_nodes[ee], ref_edges_nodes).size() == 2) {
1286 EntityHandle edge = tit_edges[ee];
1287 CHKERR set_parent(*reit, edge, refined_ents_ptr, cOre);
1288 break;
1289 }
1290 }
1291 if (ee < 3)
1292 continue; // this refined edge is contained by edge of triangle
1293
1294 // check if ref edge is in coarse tetrahedral (i.e. that is internal
1295 // edge of refined tetrahedral)
1296 if (intersect(tri_nodes, ref_edges_nodes).size() == 2) {
1297 CHKERR set_parent(*reit, tit, refined_ents_ptr, cOre);
1298 continue;
1299 }
1300
1301 // Refined edge is not child of any edge, face or tetrahedral, this is
1302 // Impossible edge
1303 SETERRQ(m_field.get_comm(), MOFEM_DATA_INCONSISTENCY,
1304 "impossible refined edge");
1305 }
1306 }
1307 }
1308
1309 CHKERR set_parent(refined_ents_ptr);
1310 CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevel(
1311 ents_to_set_bit, bit, true, verb, &ents);
1312
1313 MoFEMFunctionReturn(0);
1314}
1315
1316} // namespace MoFEM
MOFEM_TAG_AND_LOG
#define MOFEM_TAG_AND_LOG(channel, severity, tag)
Tag and log in channel.
Definition LogManager.hpp:376
NOISY
@ NOISY
Definition definitions.h:224
VERBOSE
@ VERBOSE
Definition definitions.h:222
VERY_NOISY
@ VERY_NOISY
Definition definitions.h:225
MoFEMFunctionReturnHot
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:463
MoFEMFunctionBegin
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:359
MOFEM_IMPOSSIBLE_CASE
@ MOFEM_IMPOSSIBLE_CASE
Definition definitions.h:35
MOFEM_DATA_INCONSISTENCY
@ MOFEM_DATA_INCONSISTENCY
Definition definitions.h:31
MoFEMFunctionReturn
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition definitions.h:432
CHKERR
#define CHKERR
Inline error check.
Definition definitions.h:551
MoFEMFunctionBeginHot
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition definitions.h:456
MoFEM::RefEntity_multiIndex_view_by_hashed_parent_entity
multi_index_container< boost::shared_ptr< RefEntity >, indexed_by< hashed_non_unique< const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > >, hashed_unique< tag< Composite_EntType_and_ParentEntType_mi_tag >, composite_key< boost::shared_ptr< RefEntity >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getEnt >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > > > > > RefEntity_multiIndex_view_by_hashed_parent_entity
multi-index view of RefEntity by parent entity
Definition RefEntsMultiIndices.hpp:777
MoFEM::CoreInterface::get_ref_ents
virtual const RefEntity_multiIndex * get_ref_ents() const =0
Get the ref ents object.
MoFEM::CoreInterface::get_ref_finite_elements
virtual const RefElement_multiIndex * get_ref_finite_elements() const =0
Get the ref finite elements object.
bit
auto bit
set bit
Definition hanging_node_approx.cpp:75
i
FTensor::Index< 'i', SPACE_DIM > i
Definition hcurl_divergence_operator_2d.cpp:27
j
FTensor::Index< 'j', 3 > j
Definition matrix_function.cpp:19
MoFEM::Exceptions::MoFEMErrorCode
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition Exceptions.hpp:56
MoFEM::Types::BitRefEdges
std::bitset< BITREFEDGES_SIZE > BitRefEdges
Definition Types.hpp:34
MoFEM::Types::BitRefLevel
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition Types.hpp:40
MoFEM
implementation of Data Operators for Forces and Sources
Definition Common.hpp:10
MoFEM::tet_type_3
int tet_type_3(const EntityHandle *conn, const int *split_edges, const EntityHandle *edge_new_nodes, EntityHandle *new_tets_conn)
Definition EntityRefine.cpp:565
MoFEM::RefEntity_multiIndex
multi_index_container< boost::shared_ptr< RefEntity >, indexed_by< ordered_unique< tag< Ent_mi_tag >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getEnt > >, ordered_non_unique< tag< Ent_Ent_mi_tag >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > >, ordered_non_unique< tag< Composite_EntType_and_ParentEntType_mi_tag >, composite_key< RefEntity, const_mem_fun< RefEntity, EntityType, &RefEntity::getEntType >, const_mem_fun< RefEntity, EntityType, &RefEntity::getParentEntType > > >, ordered_non_unique< tag< Composite_ParentEnt_And_EntType_mi_tag >, composite_key< RefEntity, const_mem_fun< RefEntity, EntityType, &RefEntity::getEntType >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > > > > > RefEntity_multiIndex
Definition RefEntsMultiIndices.hpp:760
MoFEM::tet_type_6
void tet_type_6(moab::Interface &moab, const EntityHandle *conn, const EntityHandle *edge_new_nodes, EntityHandle *new_tets_conn)
Definition EntityRefine.cpp:31
MoFEM::tri_type_2
MoFEMErrorCode tri_type_2(const EntityHandle *conn, const int *split_edges, const EntityHandle *edge_new_nodes, EntityHandle *new_tris_conn)
Definition EntityRefine.cpp:1177
MoFEM::debug
static const bool debug
Definition ConstrainMatrixCtx.cpp:9
MoFEM::prism_type_2
MoFEMErrorCode prism_type_2(const EntityHandle *conn, const BitRefEdges split_edges, const EntityHandle *edge_new_nodes, EntityHandle *new_prism_conn)
Definition EntityRefine.cpp:1266
MoFEM::tri_type_1
MoFEMErrorCode tri_type_1(const EntityHandle *conn, const int split_edge, const EntityHandle edge_new_node, EntityHandle *new_tris_conn)
Definition EntityRefine.cpp:1137
MoFEM::RefEntity
RefEntityTmp< 0 > RefEntity
Definition RefEntsMultiIndices.hpp:566
MoFEM::prism_type_3
MoFEMErrorCode prism_type_3(const EntityHandle *conn, const BitRefEdges split_edges, const EntityHandle *edge_new_nodes, EntityHandle *new_prism_conn)
Definition EntityRefine.cpp:1339
MoFEM::tri_type_3
MoFEMErrorCode tri_type_3(const EntityHandle *conn, const EntityHandle *edge_new_nodes, EntityHandle *new_tris_conn)
Definition EntityRefine.cpp:1111
MoFEM::prism_type_1
MoFEMErrorCode prism_type_1(const EntityHandle *conn, const BitRefEdges split_edges, const EntityHandle *edge_new_nodes, EntityHandle *new_prism_conn)
Definition EntityRefine.cpp:1230
MoFEM::RefElement_multiIndex
multi_index_container< boost::shared_ptr< RefElement >, indexed_by< ordered_unique< tag< Ent_mi_tag >, const_mem_fun< RefElement::interface_type_RefEntity, EntityHandle, &RefElement::getEnt > > > > RefElement_multiIndex
Definition RefElementMultiIndices.hpp:181
MoFEM::tet_type_2
int tet_type_2(const EntityHandle *conn, const int *split_edges, const EntityHandle *edge_new_nodes, EntityHandle *new_tets_conn)
Definition EntityRefine.cpp:988
MoFEM::RefElement_multiIndex_parents_view
multi_index_container< boost::shared_ptr< RefElement >, indexed_by< ordered_unique< tag< Ent_mi_tag >, const_mem_fun< RefElement::interface_type_RefEntity, EntityHandle, &RefElement::getEnt > >, ordered_non_unique< tag< Ent_Ent_mi_tag >, const_mem_fun< RefElement::interface_type_RefEntity, EntityHandle, &RefElement::getParentEnt > >, ordered_non_unique< tag< Composite_ParentEnt_And_BitsOfRefinedEdges_mi_tag >, composite_key< RefElement, const_mem_fun< RefElement::interface_type_RefEntity, EntityHandle, &RefElement::getParentEnt >, const_mem_fun< RefElement, int, &RefElement::getBitRefEdgesUlong > > > > > RefElement_multiIndex_parents_view
Definition MeshRefinement.cpp:30
MoFEM::tet_type_5
int tet_type_5(moab::Interface &moab, const EntityHandle *conn, const EntityHandle *edge_new_nodes, EntityHandle *new_tets_conn)
Definition EntityRefine.cpp:133
MoFEM::no_handle
const EntityHandle no_handle
No entity handle is indicated by zero handle, i.e. root meshset.
Definition Common.hpp:12
MoFEM::tet_type_4
int tet_type_4(const EntityHandle *conn, const int *split_edges, const EntityHandle *edge_new_nodes, EntityHandle *new_tets_conn)
Definition EntityRefine.cpp:239
MoFEM::get_temp_meshset_ptr
auto get_temp_meshset_ptr(moab::Interface &moab)
Create smart pointer to temporary meshset.
Definition Templates.hpp:1980
MoFEM::tet_type_1
void tet_type_1(const EntityHandle *conn, const int split_edge, const EntityHandle edge_new_node, EntityHandle *new_tets_conn)
Definition EntityRefine.cpp:1096
MoFEM::rangeInserter
auto rangeInserter(const I f, const I s, boost::function< bool(I it)> tester, boost::function< MoFEMErrorCode(I f, I s)> inserter)
Insert ranges.
Definition Templates.hpp:2039
MoFEM::RefEntity_multiIndex_view_by_ordered_parent_entity
multi_index_container< boost::shared_ptr< RefEntity >, indexed_by< ordered_non_unique< const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > >, hashed_unique< tag< Composite_EntType_and_ParentEntType_mi_tag >, composite_key< boost::shared_ptr< RefEntity >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getEnt >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > > > > > RefEntity_multiIndex_view_by_ordered_parent_entity
Definition RefEntsMultiIndices.hpp:792
I
constexpr IntegrationType I
Definition operators_tests.cpp:31
MoFEM::BitRefManager
Managing BitRefLevels.
Definition BitRefManager.hpp:21
MoFEM::Composite_EntType_and_ParentEntType_mi_tag
Definition TagMultiIndices.hpp:77
MoFEM::Composite_ParentEnt_And_BitsOfRefinedEdges_mi_tag
Definition TagMultiIndices.hpp:79
MoFEM::CoreInterface::get_moab
virtual moab::Interface & get_moab()=0
MoFEM::CoreInterface::get_comm
virtual MPI_Comm & get_comm() const =0
MoFEM::CoreInterface::get_basic_entity_data_ptr
virtual boost::shared_ptr< BasicEntityData > & get_basic_entity_data_ptr()=0
Get pointer to basic entity data.
MoFEM::CoreTmp< 0 >
Core (interface) class.
Definition Core.hpp:82
MoFEM::CoreTmp< 0 >::addPrismToDatabase
MoFEMErrorCode addPrismToDatabase(const EntityHandle prism, int verb=DEFAULT_VERBOSITY)
add prim element
Definition Core.cpp:544
MoFEM::CoreTmp< 0 >::get_th_RefBitLevel
Tag get_th_RefBitLevel() const
Definition Core.hpp:198
MoFEM::CoreTmp< 0 >::get_th_RefParentHandle
Tag get_th_RefParentHandle() const
Definition Core.hpp:197
MoFEM::CoreTmp< 0 >::get_th_RefBitEdge
Tag get_th_RefBitEdge() const
Definition Core.hpp:199
MoFEM::DeprecatedCoreInterface
Deprecated interface functions.
Definition DeprecatedCoreInterface.hpp:16
MoFEM::Ent_mi_tag
Definition TagMultiIndices.hpp:21
MoFEM::MeshRefinement::SetParent::operator()
MoFEMErrorCode operator()(const EntityHandle ent, const EntityHandle parent, const RefEntity_multiIndex *ref_ents_ptr, MoFEM::Core &cOre)
Definition MeshRefinement.cpp:209
MoFEM::MeshRefinement::SetParent::parentsToChange
map< EntityHandle, EntityHandle > parentsToChange
Definition MeshRefinement.hpp:182
MoFEM::MeshRefinement
Mesh refinement interface.
Definition MeshRefinement.hpp:26
MoFEM::MeshRefinement::refineTetsHangingNodes
MoFEMErrorCode refineTetsHangingNodes(const Range &tets, const BitRefLevel &bit, int verb=QUIET, const bool debug=false)
refine TET in the meshset
Definition MeshRefinement.cpp:276
MoFEM::MeshRefinement::refineTets
MoFEMErrorCode refineTets(const EntityHandle meshset, const BitRefLevel &bit, int verb=QUIET, const bool debug=false)
refine TET in the meshset
Definition MeshRefinement.cpp:197
MoFEM::MeshRefinement::query_interface
MoFEMErrorCode query_interface(boost::typeindex::type_index type_index, UnknownInterface **iface) const
Definition MeshRefinement.cpp:33
MoFEM::MeshRefinement::addVerticesInTheMiddleOfEdges
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
Definition MeshRefinement.cpp:42
MoFEM::MeshRefinement::refinePrisms
MoFEMErrorCode refinePrisms(const EntityHandle meshset, const BitRefLevel &bit, int verb=QUIET)
refine PRISM in the meshset
Definition MeshRefinement.cpp:722
MoFEM::MeshRefinement::refineTrisHangingNodes
MoFEMErrorCode refineTrisHangingNodes(const EntityHandle meshset, const BitRefLevel &bit, int verb=QUIET, const bool debug=false)
refine TRI in the meshset
Definition MeshRefinement.cpp:1003
MoFEM::MeshRefinement::MeshRefinement
MeshRefinement(const MoFEM::Core &core)
Definition MeshRefinement.cpp:39
MoFEM::MeshRefinement::SetEdgeBitsFun
boost::function< MoFEMErrorCode(moab::Interface &moab, RefEntity_multiIndex_view_by_ordered_parent_entity &ref_parent_ents_view, EntityHandle tet, BitRefEdges &parent_edges_bit, EntityHandle *edge_new_nodes, int *split_edges)> SetEdgeBitsFun
Functions setting edges for refinemnt on enetity level.
Definition MeshRefinement.hpp:203
MoFEM::MeshRefinement::refineTris
MoFEMErrorCode refineTris(const EntityHandle meshset, const BitRefLevel &bit, int verb=QUIET, const bool debug=false)
refine triangles in the meshset
Definition MeshRefinement.cpp:922
MoFEM::MeshRefinement::refineMeshset
MoFEMErrorCode refineMeshset(const EntityHandle meshset, const BitRefLevel &bit, const bool recursive=false, int verb=QUIET)
refine meshset, i.e. add child of refined entities to meshset
Definition MeshRefinement.cpp:905
MoFEM::RefElement_PRISM
keeps data about abstract PRISM finite element
Definition RefElementMultiIndices.hpp:67
MoFEM::RefElement
keeps data about abstract refined finite element
Definition RefElementMultiIndices.hpp:16
MoFEM::RefElement::interface_type_RefEntity
interface_RefEntity< RefEntity > interface_type_RefEntity
Definition RefElementMultiIndices.hpp:18
MoFEM::RefElement::getBitRefEdgesUlong
virtual int getBitRefEdgesUlong() const
Definition RefElementMultiIndices.hpp:25
MoFEM::RefEntityTmp< 0 >
Struct keeps handle to refined handle.
Definition RefEntsMultiIndices.hpp:141
MoFEM::RefEntity_change_parent
change parent
Definition RefEntsMultiIndices.hpp:812
MoFEM::UnknownInterface
base class for all interface classes
Definition UnknownInterface.hpp:34
MoFEM::UnknownInterface::getInterface
MoFEMErrorCode getInterface(IFACE *&iface) const
Get interface reference to pointer of interface.
Definition UnknownInterface.hpp:93
MoFEM::interface_RefEntity< RefEntity >::getParentEnt
EntityHandle getParentEnt() const
Definition RefEntsMultiIndices.hpp:615
MoFEM::interface_RefEntity< RefEntity >::getEnt
EntityHandle getEnt() const
Definition RefEntsMultiIndices.hpp:603
Generated by Doxygen 1.9.8 and hosted at University of Glasgow