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