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