v0.9.0
CutMeshInterface.cpp
Go to the documentation of this file.
1 /** \file CutMeshInterface.cpp
2  * \brief Cut mesh by surface
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 namespace MoFEM {
19 
22  UnknownInterface **iface) const {
24  *iface = NULL;
25  if (uuid == IDD_MOFEMCutMesh) {
26  *iface = const_cast<CutMeshInterface *>(this);
28  }
29  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "unknown interface");
31 }
32 
34  : cOre(const_cast<Core &>(core)) {
35  lineSearchSteps = 10;
36  nbMaxMergingCycles = 200;
38  trimFixedEdges = true;
39 }
40 
43  treeSurfPtr.reset();
45 }
46 
49  fRont = front;
51 }
52 
55  sUrface = surface;
57 }
58 
59 MoFEMErrorCode CutMeshInterface::copySurface(const Range &surface, Tag th,
60  double *shift, double *origin,
61  double *transform,
62  const std::string save_mesh) {
63  CoreInterface &m_field = cOre;
64  moab::Interface &moab = m_field.get_moab();
66  sUrface.clear();
67  std::map<EntityHandle, EntityHandle> verts_map;
68  for (Range::const_iterator tit = surface.begin(); tit != surface.end();
69  tit++) {
70  int num_nodes;
71  const EntityHandle *conn;
72  CHKERR moab.get_connectivity(*tit, conn, num_nodes, true);
73  MatrixDouble coords(num_nodes, 3);
74  if (th) {
75  CHKERR moab.tag_get_data(th, conn, num_nodes, &coords(0, 0));
76  } else {
77  CHKERR moab.get_coords(conn, num_nodes, &coords(0, 0));
78  }
79  EntityHandle new_verts[num_nodes];
80  for (int nn = 0; nn != num_nodes; nn++) {
81  if (verts_map.find(conn[nn]) != verts_map.end()) {
82  new_verts[nn] = verts_map[conn[nn]];
83  } else {
84  if (transform) {
85  ublas::matrix_row<MatrixDouble> mr(coords, nn);
86  if (origin) {
87  VectorAdaptor vec_origin(
88  3, ublas::shallow_array_adaptor<double>(3, origin));
89  mr = mr - vec_origin;
90  }
91  MatrixAdaptor mat_transform = MatrixAdaptor(
92  3, 3, ublas::shallow_array_adaptor<double>(9, transform));
93  mr = prod(mat_transform, mr);
94  if (origin) {
95  VectorAdaptor vec_origin(
96  3, ublas::shallow_array_adaptor<double>(3, origin));
97  mr = mr + vec_origin;
98  }
99  }
100  if (shift) {
101  ublas::matrix_row<MatrixDouble> mr(coords, nn);
102  VectorAdaptor vec_shift(
103  3, ublas::shallow_array_adaptor<double>(3, shift));
104  mr = mr + vec_shift;
105  }
106  CHKERR moab.create_vertex(&coords(nn, 0), new_verts[nn]);
107  verts_map[conn[nn]] = new_verts[nn];
108  }
109  }
110  EntityHandle ele;
111  CHKERR moab.create_element(MBTRI, new_verts, num_nodes, ele);
112  sUrface.insert(ele);
113  }
114  if (!save_mesh.empty())
115  CHKERR SaveData(m_field.get_moab())(save_mesh, sUrface);
117 }
118 
121  vOlume = volume;
123 }
124 
127  sUrface.merge(surface);
129 }
130 
133  vOlume.merge(volume);
135 }
136 
138  const Range &fixed_edges, const double rel_tol, const double abs_tol,
139  Tag th, const bool debug) {
140  CoreInterface &m_field = cOre;
141  moab::Interface &moab = m_field.get_moab();
143 
144  // Get cutting surface skin
145  Skinner skin(&moab);
146  Range surface_skin;
147  CHKERR skin.find_skin(0, sUrface, false, surface_skin);
148 
149  CHKERR snapSurfaceToEdges(surface_skin, fixed_edges, rel_tol, abs_tol, th,
150  debug);
151 
153 }
154 
156  const Range &fixed_edges,
157  const double rel_tol,
158  const double abs_tol,
159  Tag th, const bool debug) {
160  CoreInterface &m_field = cOre;
161  moab::Interface &moab = m_field.get_moab();
164 
165  map<EntityHandle, double> map_verts_length;
166 
167  for (auto f : surface_edges) {
168  int num_nodes;
169  const EntityHandle *conn_skin;
170  CHKERR moab.get_connectivity(f, conn_skin, num_nodes, true);
171  VectorDouble6 coords_skin(6);
172  if (th)
173  CHKERR moab.tag_get_data(th, conn_skin, num_nodes, &coords_skin[0]);
174  else
175  CHKERR moab.get_coords(conn_skin, num_nodes, &coords_skin[0]);
177  &coords_skin[0], &coords_skin[1], &coords_skin[2]);
179  &coords_skin[3], &coords_skin[4], &coords_skin[5]);
180  FTensor::Tensor1<double, 3> t_skin_delta;
181  t_skin_delta(i) = t_n1(i) - t_n0(i);
182  const double skin_edge_length = sqrt(t_skin_delta(i) * t_skin_delta(i));
183  for (int nn = 0; nn != 2; ++nn) {
184  auto it = map_verts_length.find(conn_skin[nn]);
185  if (it != map_verts_length.end())
186  it->second = std::min(it->second, skin_edge_length);
187  else
188  map_verts_length[conn_skin[nn]] = skin_edge_length;
189  }
190  }
191 
192  for (auto m : map_verts_length) {
193 
195  if (th)
196  CHKERR moab.tag_get_data(th, &m.first, 1, &t_n(0));
197  else
198  CHKERR moab.get_coords(&m.first, 1, &t_n(0));
199 
200  double min_dist = rel_tol * m.second;
201  FTensor::Tensor1<double, 3> t_min_coords;
202  CHKERR cOre.getInterface<Tools>()->findMinDistanceFromTheEdges(
203  &t_n(0), 1, fixed_edges, &min_dist, &t_min_coords(0));
204 
205  if (min_dist < rel_tol * m.second || min_dist < abs_tol) {
206  if(debug)
207  cerr << "Snap " << min_dist << endl;
208  if (th)
209  CHKERR moab.tag_set_data(th, &m.first, 1, &t_min_coords(0));
210  else
211  CHKERR moab.set_coords(&m.first, 1, &t_min_coords(0));
212  }
213  }
214 
216 }
217 
219  CoreInterface &m_field = cOre;
220  moab::Interface &moab = m_field.get_moab();
222  treeSurfPtr = boost::shared_ptr<OrientedBoxTreeTool>(
223  new OrientedBoxTreeTool(&moab, "ROOTSETSURF", true));
226 }
227 
229 CutMeshInterface::cutOnly(Range vol, const BitRefLevel cut_bit, Tag th,
230  const double tol_cut, const double tol_cut_close,
231  Range *fixed_edges, Range *corner_nodes,
232  const bool update_meshsets, const bool debug) {
233  CoreInterface &m_field = cOre;
234  moab::Interface &moab = m_field.get_moab();
236 
237  // cut mesh
238  CHKERR findEdgesToCut(vol, fixed_edges, corner_nodes, tol_cut, QUIET, debug);
239  CHKERR projectZeroDistanceEnts(fixed_edges, corner_nodes, tol_cut_close,
240  QUIET, debug);
243  if (fixed_edges)
244  CHKERR cOre.getInterface<BitRefManager>()->updateRange(*fixed_edges,
245  *fixed_edges);
246  if (corner_nodes)
247  CHKERR cOre.getInterface<BitRefManager>()->updateRange(*corner_nodes,
248  *corner_nodes);
249  if (update_meshsets)
251  ->updateAllMeshsetsByEntitiesChildren(cut_bit);
253 
254  if (debug) {
256  if (fixed_edges)
257  CHKERR SaveData(moab)("out_cut_new_fixed_edges.vtk", *fixed_edges);
258  }
259 
261 }
262 
264  const double tol_trim_close,
265  Range *fixed_edges,
266  Range *corner_nodes,
267  const bool update_meshsets,
268  const bool debug) {
269  CoreInterface &m_field = cOre;
270  moab::Interface &moab = m_field.get_moab();
272 
273  // trim mesh
274  CHKERR findEdgesToTrim(fixed_edges, corner_nodes, th, tol_trim_close, debug);
275  CHKERR trimEdgesInTheMiddle(trim_bit, debug);
276  if (fixed_edges)
277  CHKERR cOre.getInterface<BitRefManager>()->updateRange(*fixed_edges,
278  *fixed_edges);
279 
280  if (corner_nodes)
281  CHKERR cOre.getInterface<BitRefManager>()->updateRange(*corner_nodes,
282  *corner_nodes);
283 
284  if (update_meshsets)
286  ->updateAllMeshsetsByEntitiesChildren(trim_bit);
287 
288  // move nodes
290 
291  // remove faces
292  CHKERR trimSurface(fixed_edges, corner_nodes, debug);
293 
294  if (debug) {
296  Range bit2_edges;
297  CHKERR cOre.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
298  trim_bit, BitRefLevel().set(), MBEDGE, bit2_edges);
299  CHKERR SaveData(moab)("trim_fixed_edges.vtk",
300  intersect(*fixed_edges, bit2_edges));
301  }
302 
304 }
305 
307  int &first_bit, Tag th, const double tol_cut, const double tol_cut_close,
308  const double tol_trim_close, Range *fixed_edges, Range *corner_nodes,
309  const bool update_meshsets, const bool debug) {
310  CoreInterface &m_field = cOre;
312 
313  std::vector<BitRefLevel> bit_levels;
314 
315  auto get_back_bit_levels = [&]() {
316  bit_levels.push_back(BitRefLevel());
317  bit_levels.back().set(first_bit + bit_levels.size() - 1);
318  return bit_levels.back();
319  };
320 
321  BitRefLevel cut_bit = get_back_bit_levels();
322 
323  CHKERR cutOnly(unite(cutSurfaceVolumes, cutFrontVolumes), cut_bit, th,
324  tol_cut, tol_cut_close, fixed_edges, corner_nodes,
325  update_meshsets, debug);
326 
327  auto get_min_quality = [&m_field](const BitRefLevel bit, Tag th) {
328  Range tets_level; // test at level
329  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
330  bit, BitRefLevel().set(), MBTET, tets_level);
331  double min_q = 1;
332  CHKERR m_field.getInterface<Tools>()->minTetsQuality(tets_level, min_q, th);
333  return min_q;
334  };
335 
336  PetscPrintf(PETSC_COMM_WORLD, "Min quality cut %6.4g\n",
337  get_min_quality(cut_bit, th));
338 
339  Range starting_volume = cutNewVolumes;
340  Range starting_volume_nodes;
341  CHKERR m_field.get_moab().get_connectivity(starting_volume,
342  starting_volume_nodes, true);
343  std::vector<double> staring_volume_coords(3 * starting_volume_nodes.size());
344  CHKERR m_field.get_moab().get_coords(starting_volume_nodes,
345  &*staring_volume_coords.begin());
346 
347  if (th) {
348  std::vector<double> staring_volume_th_coords(3 *
349  starting_volume_nodes.size());
350  CHKERR m_field.get_moab().tag_get_data(th, starting_volume_nodes,
351  &*staring_volume_th_coords.begin());
352  CHKERR m_field.get_moab().set_coords(starting_volume_nodes,
353  &*staring_volume_th_coords.begin());
354  }
355 
356  if (th)
357  CHKERR setTagData(th);
358 
359  BitRefLevel trim_bit = get_back_bit_levels();
360 
361  CHKERR trimOnly(trim_bit, th, tol_trim_close, fixed_edges, corner_nodes,
362  update_meshsets, debug);
363 
364  PetscPrintf(PETSC_COMM_WORLD, "Min quality trim %3.2g\n",
365  get_min_quality(trim_bit, th));
366 
367  first_bit += bit_levels.size() - 1;
368 
369  if (th)
370  CHKERR m_field.get_moab().set_coords(starting_volume_nodes,
371  &*staring_volume_coords.begin());
372 
374 }
375 
377  int &first_bit, const int fraction_level, Tag th, const double tol_cut,
378  const double tol_cut_close, const double tol_trim_close, Range &fixed_edges,
379  Range &corner_nodes, const bool update_meshsets, const bool debug) {
380  CoreInterface &m_field = cOre;
382 
383  std::vector<BitRefLevel> bit_levels;
384 
385  auto get_back_bit_levels = [&]() {
386  bit_levels.push_back(BitRefLevel());
387  bit_levels.back().set(first_bit + bit_levels.size() - 1);
388  return bit_levels.back();
389  };
390 
391  if (debug) {
392  CHKERR cOre.getInterface<BitRefManager>()->writeEntitiesNotInDatabase(
393  "ents_not_in_database.vtk", "VTK", "");
394  }
395 
396  CHKERR cutAndTrim(first_bit, th, tol_cut, tol_cut_close, tol_trim_close,
397  &fixed_edges, &corner_nodes, update_meshsets, debug);
398  if (debug)
399  CHKERR cOre.getInterface<BitRefManager>()->writeEntitiesNotInDatabase(
400  "cut_trim_ents_not_in_database.vtk", "VTK", "");
401 
402  BitRefLevel bit_level1 = BitRefLevel().set(first_bit - 1);
403  BitRefLevel bit_level2 = get_back_bit_levels();
404  BitRefLevel bit_level3 = get_back_bit_levels();
405 
406  CHKERR mergeBadEdges(fraction_level, bit_level1, bit_level2, bit_level3,
407  getNewTrimSurfaces(), fixed_edges, corner_nodes, th,
408  update_meshsets, debug);
409 
410  auto get_min_quality = [&m_field, debug](const BitRefLevel bit, Tag th) {
411  Range tets_level; // test at level
412  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
413  bit, BitRefLevel().set(), MBTET, tets_level);
414  double min_q = 1;
415  CHKERR m_field.getInterface<Tools>()->minTetsQuality(tets_level, min_q, th);
416  if (min_q < 0 && debug) {
417  CHKERR m_field.getInterface<Tools>()->writeTetsWithQuality(
418  "negative_tets.vtk", "VTK", "", tets_level, th);
419  }
420  return min_q;
421  };
422 
423  PetscPrintf(PETSC_COMM_WORLD, "Min quality node merge %6.4g\n",
424  get_min_quality(bit_level3, th));
425 
426  CHKERR cOre.getInterface<BitRefManager>()->updateRange(fixed_edges,
427  fixed_edges);
428  CHKERR cOre.getInterface<BitRefManager>()->updateRange(corner_nodes,
429  corner_nodes);
430 
431  first_bit += bit_levels.size() - 1;
432 
433  if (debug) {
434  CHKERR cOre.getInterface<BitRefManager>()->writeBitLevelByType(
435  bit_level3, BitRefLevel().set(), MBTET, "out_tets_merged.vtk", "VTK",
436  "");
437  CHKERR cOre.getInterface<BitRefManager>()->writeEntitiesNotInDatabase(
438  "cut_trim_merge_ents_not_in_database.vtk", "VTK", "");
439  CHKERR SaveData(m_field.get_moab())("merged_surface.vtk", mergedSurfaces);
440  }
441 
443 }
444 
446  CoreInterface &m_field = cOre;
447  moab::Interface &moab = m_field.get_moab();
449  Skinner skin(&moab);
450  Range tets_skin;
451  CHKERR skin.find_skin(0, vOlume, false, tets_skin);
452  Range tets_skin_edges;
453  CHKERR moab.get_adjacencies(tets_skin, 1, false, tets_skin_edges,
454  moab::Interface::UNION);
455  Range surface_skin;
456  CHKERR skin.find_skin(0, sUrface, false, surface_skin);
457  fRont = subtract(surface_skin, tets_skin_edges);
458  if (debug)
459  CHKERR SaveData(moab)("front.vtk", fRont);
461 }
462 
464  const bool debug) {
465  CoreInterface &m_field = cOre;
466  moab::Interface &moab = m_field.get_moab();
468  auto tools_interface = m_field.getInterface<Tools>();
469 
470  auto create_tag = [&](const std::string name, const int dim) {
471  Tag th;
472  rval = moab.tag_get_handle(name.c_str(), th);
473  if (rval == MB_SUCCESS)
474  return th;
475  std::vector<double> def_val(dim, 0);
476  CHKERR moab.tag_get_handle(name.c_str(), dim, MB_TYPE_DOUBLE, th,
477  MB_TAG_CREAT | MB_TAG_SPARSE, &*def_val.begin());
478 
479  return th;
480  };
481 
482  auto set_vol = [&](const Range &vol_verts, std::vector<double> &coords,
483  std::vector<double> &dist_surface_vec,
484  std::vector<double> &dist_surface_normal_vec) {
486 
487  coords.resize(3 * vol_verts.size());
488  dist_surface_vec.resize(3 * vol_verts.size());
489  dist_surface_normal_vec.resize(3 * vol_verts.size());
490  CHKERR moab.get_coords(vol_verts, &*coords.begin());
491  std::srand(0);
492 
493  for (auto v : vol_verts) {
494 
495  const int index = vol_verts.index(v);
496  auto point_in = getVectorAdaptor(&coords[3 * index], 3);
497  VectorDouble3 point_out(3);
498  EntityHandle facets_out;
499  CHKERR treeSurfPtr->closest_to_location(&point_in[0], rootSetSurf,
500  &point_out[0], facets_out);
501 
502  VectorDouble3 n(3);
503  CHKERR tools_interface->getTriNormal(facets_out, &*n.begin());
504  n /= norm_2(n);
505 
506  VectorDouble3 delta = point_out - point_in;
507  if (norm_2(delta) < std::numeric_limits<double>::epsilon()) {
508  if (std::rand() % 2 == 0)
509  delta += n * std::numeric_limits<double>::epsilon();
510  else
511  delta -= n * std::numeric_limits<double>::epsilon();
512  }
513 
514  auto dist_vec = getVectorAdaptor(&dist_surface_vec[3 * index], 3);
515  noalias(dist_vec) = delta;
516 
517  auto dist_normal_vec =
518  getVectorAdaptor(&dist_surface_normal_vec[3 * index], 3);
519  noalias(dist_normal_vec) = inner_prod(delta, n) * n;
520  }
521 
523  };
524 
525  std::vector<double> coords;
526  std::vector<double> dist_surface_vec;
527  std::vector<double> dist_surface_normal_vec;
528  Range vol_verts;
529  CHKERR moab.get_connectivity(vOlume, vol_verts, true);
530 
531  CHKERR set_vol(vol_verts, coords, dist_surface_vec, dist_surface_normal_vec);
532 
533  auto th_dist_surface_vec = create_tag("DIST_SURFACE_VECTOR", 3);
534  auto th_dist_surface_normal_vec = create_tag("DIST_SURFACE_NORMAL_VECTOR", 3);
535  CHKERR moab.tag_set_data(th_dist_surface_vec, vol_verts,
536  &*dist_surface_vec.begin());
537  CHKERR moab.tag_set_data(th_dist_surface_normal_vec, vol_verts,
538  &*dist_surface_normal_vec.begin());
539 
541 }
542 
544  int verb,
545  const bool debug) {
546  CoreInterface &m_field = cOre;
547  moab::Interface &moab = m_field.get_moab();
549 
550  auto create_tag = [&](const std::string name, const int dim) {
551  Tag th;
552  rval = moab.tag_get_handle(name.c_str(), th);
553  if (rval == MB_SUCCESS)
554  return th;
555  std::vector<double> def_val(dim, 0);
556  CHKERR moab.tag_get_handle(name.c_str(), dim, MB_TYPE_DOUBLE, th,
557  MB_TAG_CREAT | MB_TAG_SPARSE, &*def_val.begin());
558  return th;
559  };
560 
561  Range vol_vertices;
562  CHKERR moab.get_connectivity(vol, vol_vertices, true);
563 
564  std::vector<double> min_distances_from_front(vol_vertices.size(), -1);
565  std::vector<double> points_on_edges(3 * vol_vertices.size(), 0);
566  std::vector<EntityHandle> closest_edges(vol_vertices.size(), 0);
567 
568  std::vector<double> coords(3 * vol_vertices.size());
569  if (th)
570  CHKERR moab.tag_get_data(th, vol_vertices, &*coords.begin());
571  else
572  CHKERR moab.get_coords(vol_vertices, &*coords.begin());
573 
574  CHKERR cOre.getInterface<Tools>()->findMinDistanceFromTheEdges(
575  &*coords.begin(), vol_vertices.size(), fRont,
576  &*min_distances_from_front.begin(), &*points_on_edges.begin(),
577  &*closest_edges.begin());
578 
579  if (!points_on_edges.empty()) {
580  for (int i = 0; i != min_distances_from_front.size(); ++i) {
581  Range faces;
582  CHKERR moab.get_adjacencies(&closest_edges[0], 1, 2, false, faces);
583  auto point_in = getVectorAdaptor(&coords[3 * i], 3);
584  auto point_out = getVectorAdaptor(&points_on_edges[3 * i], 3);
585  point_out -= point_in;
586  }
587  }
588 
589  auto th_dist_front_vec = create_tag("DIST_FRONT_VECTOR", 3);
590  CHKERR moab.tag_set_data(th_dist_front_vec, vol_vertices,
591  &*points_on_edges.begin());
592 
594 }
595 
597  Tag th, Range &vol_edges, const bool remove_adj_prims_edges, int verb,
598  const bool debug, const std::string edges_file_name) {
599  CoreInterface &m_field = cOre;
600  moab::Interface &moab = m_field.get_moab();
602 
603  auto get_tag_data = [&](auto th, auto conn) {
604  const void *ptr;
605  CHKERR moab.tag_get_by_ptr(th, &conn, 1, &ptr);
606  return getVectorAdaptor(
607  const_cast<double *>(static_cast<const double *>(ptr)), 3);
608  };
609 
610  auto get_edge_ray = [&](auto conn) {
611  VectorDouble6 coords(6);
612  CHKERR moab.get_coords(conn, 2, &coords[0]);
613  VectorAdaptor n0 = getVectorAdaptor(&coords[0], 3);
614  VectorAdaptor n1 = getVectorAdaptor(&coords[3], 3);
615  VectorDouble3 ray = n1 - n0;
616  return ray;
617  };
618 
619  Range cut_edges;
620 
621  Range edges;
622  CHKERR moab.get_adjacencies(vOlume, 1, true, edges, moab::Interface::UNION);
623 
624  auto remove_prisms_edges = [&](Range &edges) {
626  Range prisms;
627  CHKERR moab.get_adjacencies(edges, 3, false, prisms,
628  moab::Interface::UNION);
629  prisms = prisms.subset_by_type(MBPRISM);
630  Range prisms_verts;
631  CHKERR moab.get_connectivity(prisms, prisms_verts, true);
632  Range prism_edges;
633  CHKERR moab.get_adjacencies(prisms_verts, 1, false, prism_edges,
634  moab::Interface::UNION);
635  edges = subtract(edges, prism_edges);
637  };
638  if (remove_adj_prims_edges)
639  CHKERR remove_prisms_edges(edges);
640 
641  for (auto e : edges) {
642 
643  int num_nodes;
644  const EntityHandle *conn;
645  CHKERR moab.get_connectivity(e, conn, num_nodes, true);
646  auto ray = get_edge_ray(conn);
647  const double length = norm_2(ray);
648  ray /= length;
649 
650  auto signed_norm = [&](const auto &v) { return inner_prod(ray, v); };
651 
652  auto get_cut_edges = [&](auto th, Range &cut_edges) {
654  const auto dist0 = get_tag_data(th, conn[0]);
655  const auto dist1 = get_tag_data(th, conn[1]);
656  const double min_dist = std::min(norm_2(dist0), norm_2(dist1));
657  if (min_dist < length) {
658  auto opposite = inner_prod(dist0, dist1);
659  if (opposite <= 0) {
660  const double sign_dist0 = signed_norm(dist0);
661  const double sign_dist1 = signed_norm(dist1);
662  if (sign_dist0 > -std::numeric_limits<double>::epsilon() &&
663  sign_dist1 < std::numeric_limits<double>::epsilon())
664  cut_edges.insert(e);
665  }
666  }
668  };
669 
670  CHKERR get_cut_edges(th, cut_edges);
671  }
672 
673  CHKERR moab.get_adjacencies(cut_edges, 3, false, vol_edges,
674  moab::Interface::UNION);
675 
676  vol_edges = intersect(vol_edges, vOlume);
677 
678  if (debug && !edges_file_name.empty())
679  CHKERR SaveData(m_field.get_moab())(edges_file_name, cut_edges);
680 
682 }
683 
685  CoreInterface &m_field = cOre;
686  moab::Interface &moab = m_field.get_moab();
688 
689  CHKERR createFrontLevelSets(vOlume, nullptr, verb, debug);
690  Tag th_dist_front_vec;
691  CHKERR moab.tag_get_handle("DIST_FRONT_VECTOR", th_dist_front_vec);
692  CHKERR createLevelSets(th_dist_front_vec, cutFrontVolumes, true, verb, debug,
693  "cutFrontEdges.vtk");
694 
696 
697  Tag th_dist_surface_vec;
698  CHKERR moab.tag_get_handle("DIST_SURFACE_VECTOR", th_dist_surface_vec);
699  cutSurfaceVolumes.clear();
700  CHKERR createLevelSets(th_dist_surface_vec, cutSurfaceVolumes, true, verb,
701  debug, "cutSurfaceEdges.vtk");
702 
703  if (debug)
704  CHKERR SaveData(m_field.get_moab())("level_sets.vtk", vOlume);
705  if (debug)
706  CHKERR SaveData(m_field.get_moab())("cutSurfaceVolumes.vtk",
708  if (debug)
709  CHKERR SaveData(m_field.get_moab())("cutFrontVolumes.vtk", cutFrontVolumes);
710 
712 }
713 
715  const int surf_levels,
716  const int front_levels,
717  Range *fixed_edges, int verb,
718  const bool debug) {
719  CoreInterface &m_field = cOre;
720  moab::Interface &moab = m_field.get_moab();
721  MeshRefinement *refiner;
722  BitRefManager *bit_ref_manager;
724  CHKERR m_field.getInterface(refiner);
725  CHKERR m_field.getInterface(bit_ref_manager);
726 
727  auto add_bit = [&](const int bit) {
729  CHKERR bit_ref_manager->addBitRefLevel(vOlume, BitRefLevel().set(bit),
730  verb);
731  Range adj_ents;
732  for (auto d : {2, 1, 0})
733  CHKERR moab.get_adjacencies(vOlume, d, true, adj_ents,
734  moab::Interface::UNION);
735  CHKERR bit_ref_manager->addBitRefLevel(vOlume, BitRefLevel().set(bit),
736  verb);
738  };
739  CHKERR add_bit(init_bit_level);
740 
741  auto update_range = [&](Range *r_ptr) {
743  if (r_ptr) {
744  Range childs;
745  CHKERR bit_ref_manager->updateRange(*r_ptr, childs);
746  r_ptr->merge(childs);
747  }
749  };
750 
751  auto refine = [&](const BitRefLevel bit, const Range tets) {
753  Range verts;
754  CHKERR moab.get_connectivity(tets, verts, true);
755  Range ref_edges;
756  CHKERR moab.get_adjacencies(verts, 1, true, ref_edges,
757  moab::Interface::UNION);
758 
759  CHKERR refiner->add_vertices_in_the_middle_of_edges(ref_edges, bit);
760  CHKERR refiner->refine_TET(vOlume, bit, false, verb);
761 
762  CHKERR update_range(fixed_edges);
763  CHKERR update_range(&vOlume);
765  ->updateAllMeshsetsByEntitiesChildren(bit);
766 
767  Range bit_tets;
768  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
769  bit, BitRefLevel().set(), MBTET, bit_tets);
770  vOlume = intersect(vOlume, bit_tets);
771 
772  Range bit_edges;
773  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
774  bit, BitRefLevel().set(), MBEDGE, bit_edges);
775  if (fixed_edges)
776  *fixed_edges = intersect(*fixed_edges, bit_edges);
777 
779  };
780 
781  for (int ll = init_bit_level; ll != init_bit_level + surf_levels; ++ll) {
783  CHKERR refine(BitRefLevel().set(ll + 1),
785  }
786 
787  for (int ll = init_bit_level + surf_levels;
788  ll != init_bit_level + surf_levels + front_levels; ++ll) {
790  CHKERR refine(BitRefLevel().set(ll + 1), cutFrontVolumes);
791  }
792 
794 
795  if (debug)
796  CHKERR SaveData(m_field.get_moab())("refinedVolume.vtk", vOlume);
797 
799 }
800 
801 MoFEMErrorCode CutMeshInterface::findEdgesToCut(Range vol, Range *fixed_edges,
802  Range *corner_nodes,
803  const double geometry_tol,
804  int verb, const bool debug) {
805  CoreInterface &m_field = cOre;
806  moab::Interface &moab = m_field.get_moab();
808 
809  edgesToCut.clear();
810  cutEdges.clear();
811 
812  zeroDistanceVerts.clear();
813  zeroDistanceEnts.clear();
814  verticesOnCutEdges.clear();
815 
816  Tag th_test_dist;
817  if (debug) {
818  double def_val[] = {0, 0, 0};
819  CHKERR moab.tag_get_handle("TETS_DIST", 1, MB_TYPE_DOUBLE, th_test_dist,
820  MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
821  }
822 
823  Tag th_dist_normal;
824  CHKERR moab.tag_get_handle("DIST_SURFACE_NORMAL_VECTOR", th_dist_normal);
825 
826  auto get_tag_data = [&](auto th, auto conn) {
827  const void *ptr;
828  CHKERR moab.tag_get_by_ptr(th, &conn, 1, &ptr);
829  return getVectorAdaptor(
830  const_cast<double *>(static_cast<const double *>(ptr)), 3);
831  };
832 
833  Range vol_edges;
834  CHKERR moab.get_adjacencies(vol, 1, true, vol_edges, moab::Interface::UNION);
835 
836  aveLength = 0;
837  maxLength = 0;
838  int nb_ave_length = 0;
839  for (auto e : vol_edges) {
840 
841  int num_nodes;
842  const EntityHandle *conn;
843  CHKERR moab.get_connectivity(e, conn, num_nodes, true);
844 
845  VectorDouble6 coords(6);
846  CHKERR moab.get_coords(conn, 2, &coords[0]);
847  VectorAdaptor n0 = getVectorAdaptor(&coords[0], 3);
848  VectorAdaptor n1 = getVectorAdaptor(&coords[3], 3);
849  VectorDouble3 ray = n1 - n0;
850  const double ray_length = norm_2(ray);
851  ray /= ray_length;
852 
853  auto dist_vec0 = get_tag_data(th_dist_normal, conn[0]);
854  auto dist_vec1 = get_tag_data(th_dist_normal, conn[1]);
855 
856  const double s0 = norm_2(dist_vec0);
857  const double s1 = norm_2(dist_vec1);
858 
859  if (inner_prod(dist_vec0, dist_vec1) < 0) {
860 
861  // Edges is on two sides of the surface
862  const double s = s0 / (s0 + s1);
863  const double dist = s * ray_length;
864 
865  VectorDouble3 p = n0 + dist * ray;
866  VectorDouble3 w = n0 + dist_vec0;
867  VectorDouble3 v = n1 + dist_vec1;
868  double t;
869  auto res = Tools::minDistancePointFromOnSegment(&w[0], &v[0], &p[0], &t);
870  t = std::max(0., std::min(t, 1.));
871  double d = 0;
872  if (res == Tools::SOLUTION_EXIST) {
873  VectorDouble3 o = w + t * (v - w);
874  d = norm_2(o - p) / ray_length;
875  }
876 
877  if (debug)
878  CHKERR moab.tag_set_data(th_test_dist, &e, 1, &d);
879 
880  if (d < 0.25) {
881 
882  edgesToCut[e].dIst = dist;
883  edgesToCut[e].lEngth = ray_length;
884  edgesToCut[e].unitRayDir = ray;
885  edgesToCut[e].rayPoint = n0;
886  cutEdges.insert(e);
887 
888  aveLength += norm_2(ray);
889  maxLength = fmax(maxLength, norm_2(ray));
890  ++nb_ave_length;
891  }
892  }
893  }
894  aveLength /= nb_ave_length;
895 
896  auto not_project_node = [this, &moab](const EntityHandle v) {
898  VectorDouble3 s0(3);
899  CHKERR moab.get_coords(&v, 1, &s0[0]);
900  verticesOnCutEdges[v].dIst = 0;
901  verticesOnCutEdges[v].lEngth = 0;
902  verticesOnCutEdges[v].unitRayDir.resize(3, false);
903  verticesOnCutEdges[v].unitRayDir.clear();
904  verticesOnCutEdges[v].rayPoint = s0;
906  };
907 
908  auto project_node = [this, &moab](const EntityHandle v,
909  VectorDouble3 dist_normal) {
911  VectorDouble3 s0(3);
912  CHKERR moab.get_coords(&v, 1, &s0[0]);
913  verticesOnCutEdges[v].dIst = 1;
914  verticesOnCutEdges[v].lEngth = norm_2(dist_normal);
915  verticesOnCutEdges[v].unitRayDir = dist_normal;
916  verticesOnCutEdges[v].rayPoint = s0;
918  };
919 
920  auto get_ave_edge_length = [&](const EntityHandle ent,
921  const Range &vol_edges) {
922 
923  Range adj_edges;
924  if (moab.type_from_handle(ent) == MBVERTEX)
925  CHKERR moab.get_adjacencies(&ent, 1, 1, false, adj_edges);
926  else
927  adj_edges.insert(ent);
928  adj_edges = intersect(adj_edges, vol_edges);
929 
930  double ave_l = 0;
931  for (auto e : adj_edges) {
932  int num_nodes;
933  const EntityHandle *conn;
934  CHKERR moab.get_connectivity(e, conn, num_nodes, true);
935  VectorDouble6 coords(6);
936  CHKERR moab.get_coords(conn, num_nodes, &coords[0]);
938  &coords[0], &coords[1], &coords[2]);
940  &coords[3], &coords[4], &coords[5]);
942  t_n0(i) -= t_n1(i);
943  const double l = sqrt(t_n0(i) * t_n0(i));
944  ave_l += l;
945  }
946  return ave_l / adj_edges.size();
947  };
948 
949  auto project_vertices_close_to_geometry_features = [&]() {
951 
952  Range vol_vertices;
953  CHKERR moab.get_connectivity(vol, vol_vertices, true);
954 
955  Range fixed_edges_verts;
956  if (fixed_edges)
957  CHKERR moab.get_connectivity(*fixed_edges, fixed_edges_verts, true);
958  if (corner_nodes)
959  fixed_edges_verts.merge(*corner_nodes);
960 
961  Range fix_vertices = intersect(fixed_edges_verts, vol_vertices);
962 
963  for (auto v : fix_vertices) {
964 
965  VectorDouble3 dist_normal(3);
966  CHKERR moab.tag_get_data(th_dist_normal, &v, 1, &*dist_normal.begin());
967  const double dist = norm_2(dist_normal);
968 
969  const double tol = get_ave_edge_length(v, vol_edges) * geometry_tol;
970  if (dist < tol) {
971  CHKERR not_project_node(v);
972  zeroDistanceVerts.insert(v);
973  }
974  }
975 
976  Skinner skin(&moab);
977  Range tets_skin;
978  CHKERR skin.find_skin(0, vOlume, false, tets_skin);
979  Range tets_skin_verts;
980  CHKERR moab.get_connectivity(tets_skin, tets_skin_verts, true);
981 
982  for (auto v : subtract(tets_skin_verts, fix_vertices)) {
983 
984  VectorDouble3 dist_normal(3);
985  CHKERR moab.tag_get_data(th_dist_normal, &v, 1, &*dist_normal.begin());
986  const double dist = norm_2(dist_normal);
987 
988  const double tol =
989  get_ave_edge_length(v, vol_edges) * pow(geometry_tol, 2);
990  if (dist < tol) {
991  CHKERR not_project_node(v);
992  zeroDistanceVerts.insert(v);
993  }
994  }
995 
996  for (auto v : subtract(vol_vertices, tets_skin_verts)) {
997 
998  VectorDouble3 dist_normal(3);
999  CHKERR moab.tag_get_data(th_dist_normal, &v, 1, &*dist_normal.begin());
1000  const double dist = norm_2(dist_normal);
1001 
1002  const double tol =
1003  get_ave_edge_length(v, vol_edges) * pow(geometry_tol, 3);
1004  if (dist < tol) {
1005  CHKERR project_node(v, dist_normal);
1006  zeroDistanceVerts.insert(v);
1007  }
1008  }
1009 
1011  };
1012 
1013  CHKERR project_vertices_close_to_geometry_features();
1014 
1015  for (auto v : zeroDistanceVerts) {
1016  Range adj_edges;
1017  CHKERR moab.get_adjacencies(&v, 1, 1, false, adj_edges);
1018  for (auto e : adj_edges) {
1019  cutEdges.erase(e);
1020  edgesToCut.erase(e);
1021  }
1022  }
1023 
1024  if (debug)
1025  CHKERR SaveData(m_field.get_moab())("cut_edges.vtk", cutEdges);
1026 
1027  if (debug)
1028  CHKERR SaveData(m_field.get_moab())("cut_edges_zero_distance_verts.vtk",
1030 
1032 }
1033 
1035  Range *corner_nodes,
1036  const double close_tol,
1037  const int verb,
1038  const bool debug) {
1039  CoreInterface &m_field = cOre;
1040  moab::Interface &moab = m_field.get_moab();
1042 
1043  // Get entities on body skin
1044  Skinner skin(&moab);
1045  Range tets_skin;
1046  CHKERR skin.find_skin(0, vOlume, false, tets_skin);
1047  Range tets_skin_edges;
1048  CHKERR moab.get_adjacencies(tets_skin, 1, false, tets_skin_edges,
1049  moab::Interface::UNION);
1050  Range tets_skin_verts;
1051  CHKERR moab.get_connectivity(tets_skin, tets_skin_verts, true);
1052 
1053  // Get entities in volume
1054  Range vol_faces, vol_edges, vol_nodes;
1055  CHKERR moab.get_adjacencies(vOlume, 2, false, vol_faces,
1056  moab::Interface::UNION);
1057  CHKERR moab.get_adjacencies(vOlume, 1, false, vol_edges,
1058  moab::Interface::UNION);
1059  CHKERR moab.get_connectivity(vOlume, vol_nodes, true);
1060 
1061  // Get nodes on cut edges
1062  Range cut_edge_verts;
1063  CHKERR moab.get_connectivity(cutEdges, cut_edge_verts, true);
1064 
1065  // Get faces and edges
1066  Range cut_edges_faces;
1067  CHKERR moab.get_adjacencies(cut_edge_verts, 2, true, cut_edges_faces,
1068  moab::Interface::UNION);
1069  cut_edges_faces = intersect(cut_edges_faces, vol_faces);
1070  Range cut_edges_faces_verts;
1071  CHKERR moab.get_connectivity(cut_edges_faces, cut_edges_faces_verts, true);
1072  cut_edges_faces_verts = subtract(cut_edges_faces_verts, cut_edge_verts);
1073  Range to_remove_cut_edges_faces;
1074  CHKERR moab.get_adjacencies(cut_edges_faces_verts, 2, true,
1075  to_remove_cut_edges_faces,
1076  moab::Interface::UNION);
1077  // Those are faces which have vertices adjacent to cut edges vertices without
1078  // hanging nodes (nodes not adjacent to cutting edge)
1079  cut_edges_faces = subtract(cut_edges_faces, to_remove_cut_edges_faces);
1080  if (debug)
1081  CHKERR SaveData(moab)("cut_edges_faces.vtk", cut_edges_faces);
1082  cut_edges_faces.merge(cutEdges);
1083 
1084  Range fixed_edges_nodes;
1085  if (fixed_edges)
1086  CHKERR moab.get_connectivity(*fixed_edges, fixed_edges_nodes, true);
1087 
1088  Tag th_dist_normal;
1089  CHKERR moab.tag_get_handle("DIST_SURFACE_NORMAL_VECTOR", th_dist_normal);
1090 
1091  // Create map of closes points to the surface
1092  enum TYPE { FREE = 0, SKIN, FIXED, CORNER, NOT_MOVE };
1093  map<EntityHandle, std::pair<std::pair<TYPE, EntityHandle>, TreeData>>
1094  min_dist_map;
1095  double ave_cut_edge_length = 0;
1096  for (auto e : cutEdges) {
1097 
1098  auto eit = edgesToCut.find(e);
1099  if (eit != edgesToCut.end()) {
1100 
1101  TYPE edge_type = FREE;
1102  if (tets_skin_edges.find(e) != tets_skin_edges.end())
1103  edge_type = SKIN;
1104  if (fixed_edges)
1105  if (fixed_edges->find(e) != fixed_edges->end())
1106  edge_type = FIXED;
1107 
1108  int num_nodes;
1109  const EntityHandle *conn;
1110  CHKERR moab.get_connectivity(e, conn, num_nodes, true);
1111  VectorDouble6 pos(6);
1112  CHKERR moab.get_coords(conn, num_nodes, &pos[0]);
1113  VectorDouble3 p[2];
1114  p[0] = getVectorAdaptor(&pos[0], 3);
1115  p[1] = getVectorAdaptor(&pos[3], 3);
1116  ave_cut_edge_length += norm_2(p[0] - p[1]);
1117 
1118  auto &d = eit->second;
1119  VectorDouble3 new_pos = d.rayPoint + d.dIst * d.unitRayDir;
1120  for (int nn = 0; nn != 2; ++nn) {
1121 
1122  VectorDouble3 ray = new_pos - p[nn];
1123  const double dist = norm_2(ray);
1124  const double length = dist;
1125 
1126  bool add_node = true;
1127  auto vit = min_dist_map.find(conn[nn]);
1128  if (vit != min_dist_map.end()) {
1129  if (vit->second.second.dIst < dist)
1130  add_node = false;
1131  }
1132 
1133  if (add_node) {
1134  min_dist_map[conn[nn]].first.first = edge_type;
1135  min_dist_map[conn[nn]].first.second = e;
1136  auto &data = min_dist_map[conn[nn]].second;
1137  data.lEngth = length;
1138  data.rayPoint = p[nn];
1139  data.dIst = dist;
1140  if (dist > 0)
1141  data.unitRayDir = ray / dist;
1142  else {
1143  data.unitRayDir.resize(3);
1144  data.unitRayDir.clear();
1145  }
1146  }
1147  }
1148 
1149  } else
1150  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Edge not found");
1151  }
1152 
1153  ave_cut_edge_length /= cutEdges.size();
1154 
1155  auto get_min_quality =
1156  [&](const Range &adj_tets,
1157  const map<EntityHandle, TreeData> &vertices_on_cut_edges) {
1158  double min_q = 1;
1159  for (auto t : adj_tets) {
1160  int num_nodes;
1161  const EntityHandle *conn;
1162  CHKERR m_field.get_moab().get_connectivity(t, conn, num_nodes, true);
1163  VectorDouble12 coords(12);
1164  CHKERR moab.get_coords(conn, num_nodes, &coords[0]);
1165  for (auto n : {0, 1, 2, 3}) {
1166  auto mit = vertices_on_cut_edges.find(conn[n]);
1167  if (mit != vertices_on_cut_edges.end()) {
1168  auto n_coords = getVectorAdaptor(&coords[3 * n], 3);
1169  const double dist = mit->second.dIst;
1170  noalias(n_coords) =
1171  mit->second.rayPoint + dist * mit->second.unitRayDir;
1172  }
1173  }
1174  min_q = std::min(min_q, Tools::volumeLengthQuality(&coords[0]));
1175  }
1176  return min_q;
1177  };
1178 
1179  auto get_quality_change =
1180  [&](const Range &adj_tets,
1181  map<EntityHandle, TreeData> vertices_on_cut_edges) {
1182  double q0 = get_min_quality(adj_tets, verticesOnCutEdges);
1183  vertices_on_cut_edges.insert(verticesOnCutEdges.begin(),
1184  verticesOnCutEdges.end());
1185  double q = get_min_quality(adj_tets, vertices_on_cut_edges);
1186  return q / q0;
1187  };
1188 
1189  auto get_conn = [&moab](const EntityHandle &e, const EntityHandle *&conn,
1190  int &num_nodes) {
1192  EntityType type = moab.type_from_handle(e);
1193  if (type == MBVERTEX) {
1194  conn = &e;
1195  num_nodes = 1;
1196  } else {
1197  CHKERR moab.get_connectivity(e, conn, num_nodes, true);
1198  }
1200  };
1201 
1202  auto project_node = [&](const EntityHandle v,
1203  map<EntityHandle, TreeData> &vertices_on_cut_edges) {
1205 
1206  vertices_on_cut_edges[v].dIst = min_dist_map[v].second.dIst;
1207  vertices_on_cut_edges[v].lEngth = min_dist_map[v].second.lEngth;
1208  vertices_on_cut_edges[v].unitRayDir = min_dist_map[v].second.unitRayDir;
1209  vertices_on_cut_edges[v].rayPoint = min_dist_map[v].second.rayPoint;
1210 
1212  };
1213 
1214  auto remove_surface_tets = [&](Range &zero_dist_tets,
1215  Range zero_distance_ents,
1216  Range zero_distance_verts) {
1218  Range zero_dist_all_verts;
1219  CHKERR moab.get_connectivity(zero_distance_ents, zero_dist_all_verts, true);
1220  zero_dist_all_verts.merge(zero_distance_verts);
1221  CHKERR moab.get_adjacencies(zero_dist_all_verts, 3, false, zero_dist_tets,
1222  moab::Interface::UNION);
1223  zero_dist_tets = intersect(zero_dist_tets, vOlume);
1224  Range zero_tets_verts;
1225  CHKERR moab.get_connectivity(zero_dist_tets, zero_tets_verts, true);
1226  zero_tets_verts = subtract(zero_tets_verts, zero_dist_all_verts);
1227  Range free_tets;
1228  CHKERR moab.get_adjacencies(zero_tets_verts, 3, false, free_tets,
1229  moab::Interface::UNION);
1230  zero_dist_tets = subtract(zero_dist_tets, free_tets);
1231 
1233  };
1234 
1235  for (int d = 2; d >= 0; --d) {
1236 
1237  Range ents;
1238  if (d > 0)
1239  ents = cut_edges_faces.subset_by_dimension(d);
1240  else
1241  ents = cut_edge_verts;
1242 
1243  // make list of entities
1244  multimap<double, EntityHandle> ents_to_check;
1245  for (auto f : ents) {
1246  int num_nodes;
1247  const EntityHandle *conn;
1248  CHKERR get_conn(f, conn, num_nodes);
1249  VectorDouble3 dist(3);
1250 
1251  for (int n = 0; n != num_nodes; ++n) {
1252  auto &d = min_dist_map[conn[n]];
1253  dist[n] = d.second.lEngth;
1254  }
1255 
1256  double max_dist = 0;
1257  for (int n = 0; n != num_nodes; ++n)
1258  max_dist = std::max(max_dist, fabs(dist[n]));
1259  if (max_dist < close_tol * ave_cut_edge_length)
1260  ents_to_check.insert(std::pair<double, EntityHandle>(max_dist, f));
1261  }
1262 
1263  if (debug) {
1264  Range ents;
1265  for (auto m : ents_to_check)
1266  ents.insert(m.second);
1267  CHKERR SaveData(moab)("ents_to_check_to_project.vtk", ents);
1268  }
1269 
1270  double ray_point[3], unit_ray_dir[3];
1271  VectorAdaptor vec_unit_ray_dir(
1272  3, ublas::shallow_array_adaptor<double>(3, unit_ray_dir));
1273  VectorAdaptor vec_ray_point(
1274  3, ublas::shallow_array_adaptor<double>(3, ray_point));
1275 
1276  for (auto m : ents_to_check) {
1277 
1278  EntityHandle f = m.second;
1279 
1280  int num_nodes;
1281  const EntityHandle *conn;
1282  CHKERR get_conn(f, conn, num_nodes);
1283  VectorDouble9 coords(9);
1284  CHKERR moab.get_coords(conn, num_nodes, &coords[0]);
1285 
1286  Range adj_tets;
1287  CHKERR moab.get_adjacencies(conn, num_nodes, 3, false, adj_tets,
1288  moab::Interface::UNION);
1289  adj_tets = intersect(adj_tets, vOlume);
1290 
1291  map<EntityHandle, TreeData> vertices_on_cut_edges;
1292  for (int n = 0; n != num_nodes; ++n)
1293  CHKERR project_node(conn[n], vertices_on_cut_edges);
1294 
1295  const double q = get_quality_change(adj_tets, vertices_on_cut_edges);
1296  if (std::isnormal(q) && q > projectEntitiesQualityTrashold) {
1297  EntityHandle type = moab.type_from_handle(f);
1298 
1299  Range zero_dist_tets;
1300  if (type == MBVERTEX) {
1301  Range zero_distance_verts_test = zeroDistanceVerts;
1302  zero_distance_verts_test.insert(f);
1303  CHKERR remove_surface_tets(zero_dist_tets, zeroDistanceEnts,
1304  zero_distance_verts_test);
1305  } else {
1306  Range zero_distance_ents_test = zeroDistanceEnts;
1307  zero_distance_ents_test.insert(f);
1308  CHKERR remove_surface_tets(zero_dist_tets, zero_distance_ents_test,
1310  }
1311 
1312  if (zero_dist_tets.empty()) {
1313  verticesOnCutEdges.insert(vertices_on_cut_edges.begin(),
1314  vertices_on_cut_edges.end());
1315  if (type == MBVERTEX) {
1316  zeroDistanceVerts.insert(f);
1317  } else {
1318  zeroDistanceEnts.insert(f);
1319  }
1320  }
1321  }
1322  }
1323  }
1324 
1325  for (auto &v : verticesOnCutEdges) {
1326 
1327  TYPE node_type;
1328 
1329  if (corner_nodes->find(v.first) != corner_nodes->end())
1330  node_type = CORNER;
1331  else if (fixed_edges_nodes.find(v.first) != fixed_edges_nodes.end())
1332  node_type = FIXED;
1333  else if (tets_skin_verts.find(v.first) != tets_skin_verts.end())
1334  node_type = SKIN;
1335  else
1336  node_type = FREE;
1337 
1338  if (node_type > min_dist_map[v.first].first.first)
1339  v.second.unitRayDir.clear();
1340  }
1341 
1342  for (auto f : unite(zeroDistanceEnts, zeroDistanceVerts)) {
1343  int num_nodes;
1344  const EntityHandle *conn;
1345  CHKERR get_conn(f, conn, num_nodes);
1346  Range adj_edges;
1347  CHKERR moab.get_adjacencies(conn, num_nodes, 1, false, adj_edges,
1348  moab::Interface::UNION);
1349  for (auto e : adj_edges) {
1350  cutEdges.erase(e);
1351  edgesToCut.erase(e);
1352  }
1353  }
1354 
1355  if (debug)
1356  SaveData(m_field.get_moab())("projected_zero_distance_ents.vtk",
1358 
1360 }
1361 
1363  Range &cut_vols,
1364  Range &cut_surf,
1365  Range &cut_verts,
1366  const bool debug) {
1367  CoreInterface &m_field = cOre;
1368  moab::Interface &moab = m_field.get_moab();
1369  MeshRefinement *refiner;
1370  const RefEntity_multiIndex *ref_ents_ptr;
1372 
1373  if (cutEdges.size() != edgesToCut.size())
1374  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
1375 
1376  CHKERR m_field.getInterface(refiner);
1377  CHKERR m_field.get_ref_ents(&ref_ents_ptr);
1379  CHKERR refiner->refine_TET(vOlume, bit, false, QUIET,
1380  debug ? &cutEdges : NULL);
1381 
1382  if (debug)
1383  if (cutEdges.size() != edgesToCut.size())
1384  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
1385 
1386  cut_vols.clear();
1387  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1388  bit, BitRefLevel().set(), MBTET, cut_vols);
1389  cut_surf.clear();
1390  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1391  bit, BitRefLevel().set(), MBTRI, cut_surf);
1392 
1393  // Find new vertices on cut edges
1394  cut_verts.clear();
1395  CHKERR moab.get_connectivity(zeroDistanceEnts, cut_verts, true);
1396  cut_verts.merge(zeroDistanceVerts);
1397  for (auto &m : edgesToCut) {
1398  auto vit = ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().find(
1399  boost::make_tuple(MBVERTEX, m.first));
1400  if (vit ==
1401  ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().end()) {
1402  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1403  "No vertex on cut edges, that make no sense");
1404  }
1405  const boost::shared_ptr<RefEntity> &ref_ent = *vit;
1406  if ((ref_ent->getBitRefLevel() & bit).any()) {
1407  EntityHandle vert = ref_ent->getRefEnt();
1408  cut_verts.insert(vert);
1409  verticesOnCutEdges[vert] = m.second;
1410  } else {
1411  SETERRQ1(
1412  PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1413  "Vertex has wrong bit ref level %s",
1414  boost::lexical_cast<std::string>(ref_ent->getBitRefLevel()).c_str());
1415  }
1416  }
1417 
1418  // Add zero distance entities faces
1419  Range tets_skin;
1420  Skinner skin(&moab);
1421  CHKERR skin.find_skin(0, cut_vols, false, tets_skin);
1422  cut_surf.merge(zeroDistanceEnts.subset_by_type(MBTRI));
1423 
1424  // At that point cut_surf has all newly created faces, now take all
1425  // nodes on those faces and subtract nodes on cut edges. Faces adjacent to
1426  // nodes which left are not part of surface.
1427  Range diff_verts;
1428  CHKERR moab.get_connectivity(unite(cut_surf, zeroDistanceEnts), diff_verts,
1429  true);
1430  diff_verts = subtract(diff_verts, cut_verts);
1431  Range subtract_faces;
1432  CHKERR moab.get_adjacencies(diff_verts, 2, false, subtract_faces,
1433  moab::Interface::UNION);
1434  cut_surf = subtract(cut_surf, unite(subtract_faces, tets_skin));
1435  cut_verts.clear();
1436  CHKERR moab.get_connectivity(cut_surf, cut_verts, true);
1437 
1438  // Check non-mainfolds
1439  auto check_for_non_minfold = [&]() {
1441  Range surf_edges;
1442  CHKERR moab.get_adjacencies(cut_surf, 1, false, surf_edges,
1443  moab::Interface::UNION);
1444  for (auto e : surf_edges) {
1445 
1446  Range faces;
1447  CHKERR moab.get_adjacencies(&e, 1, 2, false, faces);
1448  faces = intersect(faces, cut_surf);
1449  if (faces.size() > 2) {
1450 
1451  bool resolved = false;
1452 
1453  // Check for haning node
1454  Range nodes;
1455  CHKERR moab.get_connectivity(faces, nodes, true);
1456  for (auto n : nodes) {
1457  Range adj_faces;
1458  CHKERR moab.get_adjacencies(&n, 1, 2, false, adj_faces);
1459  adj_faces = intersect(adj_faces, cut_surf);
1460  if (adj_faces.size() == 1) {
1461  cut_surf.erase(adj_faces[0]);
1462  resolved = true;
1463  }
1464  }
1465 
1466  // Check for two edges minfold
1467  Range adj_edges;
1468  CHKERR moab.get_adjacencies(faces, 1, false, adj_edges,
1469  moab::Interface::UNION);
1470  adj_edges = intersect(adj_edges, surf_edges);
1471  adj_edges.erase(e);
1472  for (auto other_e : adj_edges) {
1473  Range other_faces;
1474  CHKERR moab.get_adjacencies(&other_e, 1, 2, false, other_faces);
1475  other_faces = intersect(other_faces, cut_surf);
1476  if (other_faces.size() > 2) {
1477  other_faces = intersect(other_faces, faces);
1478  cut_surf = subtract(cut_surf, other_faces);
1479  resolved = true;
1480  }
1481  }
1482 
1483  if (!resolved && !debug)
1484  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY,
1485  "Non-mainfold surface");
1486 
1487  cut_verts.clear();
1488  CHKERR moab.get_connectivity(cut_surf, cut_verts, true);
1489  }
1490  }
1492  };
1493 
1494  CHKERR check_for_non_minfold();
1495 
1497 }
1498 
1501 
1502  CoreInterface &m_field = cOre;
1503  moab::Interface &moab = m_field.get_moab();
1505 
1506  // Range out_side_vertices;
1507  for (auto m : verticesOnCutEdges) {
1508  double dist = m.second.dIst;
1509  VectorDouble3 new_coors = m.second.rayPoint + dist * m.second.unitRayDir;
1510  if (th)
1511  CHKERR moab.tag_set_data(th, &m.first, 1, &new_coors[0]);
1512  else
1513  CHKERR moab.set_coords(&m.first, 1, &new_coors[0]);
1514  }
1515 
1517 }
1518 
1520  CoreInterface &m_field = cOre;
1521  moab::Interface &moab = m_field.get_moab();
1523  for (auto &v : verticesOnTrimEdges) {
1524  double dist = v.second.dIst;
1525  VectorDouble3 new_coors = v.second.rayPoint + dist * v.second.unitRayDir;
1526  if (th)
1527  CHKERR moab.tag_set_data(th, &v.first, 1, &new_coors[0]);
1528  else
1529  CHKERR moab.set_coords(&v.first, 1, &new_coors[0]);
1530  }
1532 }
1533 
1535  Range *corner_nodes, Tag th,
1536  const double tol,
1537  const bool debug) {
1538  CoreInterface &m_field = cOre;
1539  moab::Interface &moab = m_field.get_moab();
1541 
1542  // takes body skin
1543  Skinner skin(&moab);
1544  Range tets_skin;
1545  CHKERR skin.find_skin(0, cutNewVolumes, false, tets_skin);
1546 
1547  // vertices on the skin
1548  Range tets_skin_verts;
1549  CHKERR moab.get_connectivity(tets_skin, tets_skin_verts, true);
1550  // edges on the skin
1551  Range tets_skin_edges;
1552  CHKERR moab.get_adjacencies(tets_skin, 1, false, tets_skin_edges,
1553  moab::Interface::UNION);
1554  // get edges on new surface
1555  Range cut_surface_edges;
1556  CHKERR moab.get_adjacencies(cutNewSurfaces, 1, false, cut_surface_edges,
1557  moab::Interface::UNION);
1558  Range cut_surface_verts;
1559  CHKERR moab.get_connectivity(cutNewSurfaces, cut_surface_verts, true);
1560 
1561  Range corners;
1562  if (corner_nodes)
1563  corners = *corner_nodes;
1564 
1565  Range fix_edges;
1566  if (fixed_edges)
1567  fix_edges = *fixed_edges;
1568 
1569  Range fixed_edges_verts;
1570  if (fixed_edges)
1571  CHKERR moab.get_connectivity(*fixed_edges, fixed_edges_verts, true);
1572 
1573  Range surface_skin;
1574  if (fRont.empty())
1575  CHKERR skin.find_skin(0, sUrface, false, surface_skin);
1576  else
1577  surface_skin = fRont;
1578 
1579  auto get_point_coords = [&](EntityHandle v) {
1580  VectorDouble3 point(3);
1581  if (th)
1582  CHKERR moab.tag_get_data(th, &v, 1, &point[0]);
1583  else
1584  CHKERR moab.get_coords(&v, 1, &point[0]);
1585  return point;
1586  };
1587 
1588  struct VertMap {
1589  const double d;
1590  const EntityHandle v;
1591  const EntityHandle e;
1592  VertMap(const double d, const EntityHandle v, const EntityHandle e)
1593  : d(d), v(v), e(e) {}
1594  };
1595 
1596  typedef multi_index_container<
1597  VertMap,
1598  indexed_by<
1599  ordered_non_unique<member<VertMap, const double, &VertMap::d>>,
1600  ordered_non_unique<member<VertMap, const EntityHandle, &VertMap::v>>,
1601  ordered_non_unique<member<VertMap, const EntityHandle, &VertMap::e>>
1602 
1603  >>
1604  VertMapMultiIndex;
1605 
1606  VertMapMultiIndex verts_map;
1607 
1608  auto add_vert = [&](const EntityHandle v, const EntityHandle e,
1609  const double dist) {
1610  verts_map.insert(VertMap(dist, v, e));
1611  };
1612 
1613  // clear data ranges
1614  trimEdges.clear();
1615  edgesToTrim.clear();
1616  verticesOnTrimEdges.clear();
1617  trimNewVertices.clear();
1618 
1619  auto cut_this_edge = [&](const EntityHandle e, const double length, auto &ray,
1620  auto &ray_point) {
1621  trimEdges.insert(e);
1622  edgesToTrim[e].dIst = 1;
1623  edgesToTrim[e].lEngth = length;
1624  edgesToTrim[e].unitRayDir.resize(3, false);
1625  edgesToTrim[e].rayPoint.resize(3, false);
1626  for (int ii = 0; ii != 3; ++ii)
1627  edgesToTrim[e].unitRayDir[ii] = ray(ii);
1628  for (int ii = 0; ii != 3; ++ii)
1629  edgesToTrim[e].rayPoint[ii] = ray_point(ii);
1630  };
1631 
1633  int num_nodes;
1634 
1635  auto get_tag_data = [&](auto th, auto conn) {
1637  CHKERR moab.tag_get_data(th, &conn, 1, &t(0));
1638  return t;
1639  };
1640 
1641  double max_edge_length = 0;
1642 
1643  if (!fRont.empty()) {
1644  // Calculate distances
1645  treeSurfPtr = boost::shared_ptr<OrientedBoxTreeTool>(
1646  new OrientedBoxTreeTool(&moab, "ROOTSETSURF", true));
1648 
1649  for (auto s : surface_skin) {
1650 
1651  auto project_on_surface = [&](auto &t) {
1653 
1654  EntityHandle facets_out;
1655  FTensor::Tensor1<double, 3> t_point_on_cutting_surface;
1656  CHKERR treeSurfPtr->closest_to_location(&t(0), rootSetSurf, &t_p(0),
1657  facets_out);
1658 
1659  FTensor::Tensor1<double,3> t_normal;
1660  CHKERR m_field.getInterface<Tools>()->getTriNormal(facets_out,
1661  &t_normal(0));
1662  t_normal(i) /= sqrt(t_normal(i) * t_normal(i));
1663  const double dot = t_normal(i) * (t_p(i) - t(i));
1664  t_p(i) = t(i) + dot * t_normal(i);
1665 
1666  return t_p;
1667  };
1668 
1669  const EntityHandle *conn;
1670  CHKERR moab.get_connectivity(s, conn, num_nodes, true);
1671 
1672  VectorDouble6 coords_front(6);
1673  CHKERR moab.get_coords(conn, num_nodes, &coords_front[0]);
1674 
1675  FTensor::Tensor1<double *, 3> t_f0(&coords_front[0], &coords_front[1],
1676  &coords_front[2]);
1677  FTensor::Tensor1<double *, 3> t_f1(&coords_front[3], &coords_front[4],
1678  &coords_front[5]);
1679 
1680  auto t_p_f0 = project_on_surface(t_f0);
1681  auto t_p_f1 = project_on_surface(t_f1);
1682 
1683  CHKERR moab.set_coords(&conn[0], 1, &t_p_f0(0));
1684  CHKERR moab.set_coords(&conn[1], 1, &t_p_f1(0));
1685  }
1686  }
1687 
1688  if (debug)
1689  CHKERR SaveData(moab)("surface_skin_to_trim.vtk", surface_skin);
1690 
1692  Tag th_dist_front_vec;
1693  CHKERR moab.tag_get_handle("DIST_FRONT_VECTOR", th_dist_front_vec);
1694 
1695  if (debug)
1696  CHKERR SaveData(moab)("edges_potentially_to_trim.vtk", cut_surface_edges);
1697 
1698  // iterate over edges on cut surface
1699  for (auto e : cut_surface_edges) {
1700 
1701  // Get edge connectivity and coords
1702  const EntityHandle *conn_edge;
1703  CHKERR moab.get_connectivity(e, conn_edge, num_nodes, true);
1704 
1705  auto t_dist0 = get_tag_data(th_dist_front_vec, conn_edge[0]);
1706  auto t_dist1 = get_tag_data(th_dist_front_vec, conn_edge[1]);
1707 
1708  double coords_edge[3 * num_nodes];
1709  CHKERR moab.get_coords(conn_edge, num_nodes, coords_edge);
1710 
1711  FTensor::Tensor1<double *, 3> t_e0(&coords_edge[0], &coords_edge[1],
1712  &coords_edge[2]);
1713  FTensor::Tensor1<double *, 3> t_e1(&coords_edge[3], &coords_edge[4],
1714  &coords_edge[5]);
1715 
1716  FTensor::Tensor1<double, 3> t_edge_delta;
1717  t_edge_delta(i) = t_e1(i) - t_e0(i);
1718  const double edge_length2 = t_edge_delta(i) * t_edge_delta(i);
1719  const double edge_length = sqrt(edge_length2);
1720  if (edge_length == 0)
1721  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Zero edge length");
1722 
1723  max_edge_length = std::max(max_edge_length, edge_length);
1724  const double s0 = t_dist0(i) * t_edge_delta(i) / edge_length;
1725  const double s1 = t_dist1(i) * t_edge_delta(i) / edge_length;
1726  const double dot = t_dist0(i) * t_dist1(i);
1727 
1728  // iterate entities surface front to find cross to trim
1729  for (auto s : surface_skin) {
1730 
1731  auto get_edge_coors = [&](const auto e) {
1732  const EntityHandle *conn;
1733  CHKERR moab.get_connectivity(e, conn, num_nodes, true);
1734  VectorDouble6 coords(6);
1735  CHKERR moab.get_coords(conn, num_nodes, &coords[0]);
1736  return coords;
1737  };
1738 
1739  auto coords_front = get_edge_coors(s);
1740 
1741  FTensor::Tensor1<double *, 3> t_f0(&coords_front[0], &coords_front[1],
1742  &coords_front[2]);
1743  FTensor::Tensor1<double *, 3> t_f1(&coords_front[3], &coords_front[4],
1744  &coords_front[5]);
1745 
1746  // find point of minilam distance between front and cut surface edge
1747  double t_edge = -1, t_front = -1;
1748  auto res = Tools::minDistanceFromSegments(&t_e0(0), &t_e1(0), &t_f0(0),
1749  &t_f1(0), &t_edge, &t_front);
1750 
1751  if (res != Tools::NO_SOLUTION) {
1752  // check if edges crossing each other in the middle (it not imply that
1753  // have common point)
1754  const double overlap_tol = 1e-2;
1755  if (
1756 
1757  (t_edge > -std::numeric_limits<float>::epsilon() &&
1758  t_edge < 1 + std::numeric_limits<float>::epsilon())
1759 
1760  &&
1761 
1762  (t_front >= -overlap_tol && t_front <= 1 + overlap_tol)
1763 
1764  ) {
1765 
1766  t_front = std::max(0., std::min(t_front, 1.));
1767 
1768  FTensor::Tensor1<double, 3> t_front_delta;
1769  t_front_delta(i) = t_f1(i) - t_f0(i);
1770  FTensor::Tensor1<double, 3> t_edge_point, t_front_point;
1771  t_edge_point(i) = t_e0(i) + t_edge * t_edge_delta(i);
1772  t_front_point(i) = t_f0(i) + t_front * t_front_delta(i);
1773 
1775  t_ray(i) = t_front_point(i) - t_edge_point(i);
1776  const double dist = sqrt(t_ray(i) * t_ray(i));
1777 
1778  // that imply that edges have common point
1779  if ((dist < 0.1 * edge_length) ||
1780  (s0 >= 0 && s1 < -std::numeric_limits<double>::epsilon() &&
1781  dot <= 0)) {
1782 
1783  auto check_to_add_edge = [&](const EntityHandle e,
1784  const double dist) {
1785  auto eit = edgesToTrim.find(e);
1786  if (eit != edgesToTrim.end())
1787  if (eit->second.dIst < dist)
1788  return false;
1789  return true;
1790  };
1791 
1792  if (t_edge < 0.5) {
1793  t_ray(i) = t_edge * t_edge_delta(i);
1794  if (check_to_add_edge(e, dist)) {
1795  add_vert(conn_edge[0], e, fabs(t_edge));
1796  add_vert(conn_edge[1], e, fabs(t_edge - 1));
1797  cut_this_edge(e, edge_length, t_ray, t_e0);
1798  }
1799  } else {
1800  FTensor::Tensor1<double, 3> t_edge_point;
1801  t_edge_point(i) = t_e0(i) + t_edge * t_edge_delta(i);
1802  t_ray(i) = t_edge_point(i) - t_e1(i);
1803  if (check_to_add_edge(e, dist)) {
1804  add_vert(conn_edge[0], e, fabs(t_edge));
1805  add_vert(conn_edge[1], e, fabs(t_edge - 1));
1806  cut_this_edge(e, edge_length, t_ray, t_e1);
1807  }
1808  }
1809  }
1810  }
1811  }
1812  }
1813  }
1814 
1815  if (debug)
1816  CHKERR SaveData(moab)("edges_selected_to_trim.vtk", trimEdges);
1817 
1818  auto get_quality_change = [&](const Range &adj_tets, const EntityHandle &v,
1819  const VectorDouble3 &new_pos) {
1820  double q0 = 1;
1821  CHKERR m_field.getInterface<Tools>()->minTetsQuality(adj_tets, q0, th);
1822 
1823  double min_q = 1;
1824  for (auto t : adj_tets) {
1825  int num_nodes;
1826  const EntityHandle *conn;
1827  CHKERR m_field.get_moab().get_connectivity(t, conn, num_nodes, true);
1828  VectorDouble12 coords(12);
1829  if (th)
1830  CHKERR moab.tag_get_data(th, conn, num_nodes, &coords[0]);
1831  else
1832  CHKERR moab.get_coords(conn, num_nodes, &coords[0]);
1833 
1834  for (int n = 0; n != 4; ++n) {
1835  auto n_coords = getVectorAdaptor(&coords[3 * n], 3);
1836  if (conn[n] == v) {
1837  noalias(n_coords) = new_pos;
1838  } else {
1839  auto m = verticesOnTrimEdges.find(conn[n]);
1840  if (m != verticesOnTrimEdges.end())
1841  noalias(n_coords) =
1842  m->second.rayPoint + m->second.dIst * m->second.unitRayDir;
1843  }
1844  }
1845 
1846  const double q = Tools::volumeLengthQuality(&coords[0]);
1847  if (!std::isnormal(q)) {
1848  min_q = -2;
1849  break;
1850  }
1851  min_q = std::min(min_q, q);
1852  }
1853  return min_q / q0;
1854  };
1855 
1856  Range checked_verts;
1857  auto add_trim_vert = [&](const EntityHandle v, const EntityHandle e) {
1858  if (!(trimNewVertices.find(v) != trimNewVertices.end())) {
1859  auto &r = edgesToTrim.at(e);
1860  VectorDouble3 ray_point = get_point_coords(v);
1861  VectorDouble3 new_pos = r.rayPoint + r.dIst * r.unitRayDir;
1862  VectorDouble3 unit_ray_dir = ray_point - new_pos;
1863  Range adj_tets;
1864  CHKERR moab.get_adjacencies(&v, 1, 3, false, adj_tets);
1865  adj_tets = intersect(adj_tets, cutNewVolumes);
1866  double q = get_quality_change(adj_tets, v, new_pos);
1868  VectorDouble3 ray_dir = new_pos - ray_point;
1869  double dist = norm_2(ray_dir);
1870  verticesOnTrimEdges[v].dIst = 1;
1871  verticesOnTrimEdges[v].lEngth = dist;
1872  verticesOnTrimEdges[v].unitRayDir = ray_dir;
1873  verticesOnTrimEdges[v].rayPoint = ray_point;
1874  trimNewVertices.insert(v);
1875  }
1876  checked_verts.insert(v);
1877  }
1878  };
1879 
1880  auto add_no_move_trim = [&](const EntityHandle v, const EntityHandle e) {
1881  if (!(trimNewVertices.find(v) != trimNewVertices.end())) {
1882  auto &m = edgesToTrim.at(e);
1883  verticesOnTrimEdges[v] = m;
1884  verticesOnTrimEdges[v].rayPoint = get_point_coords(v);
1885  verticesOnTrimEdges[v].lEngth = 0;
1886  verticesOnTrimEdges[v].dIst = 0;
1887  trimNewVertices.insert(v);
1888  checked_verts.insert(v);
1889  }
1890  };
1891 
1892  // Iterate over all vertives close to surface front and check if those can
1893  // be moved
1894 
1895  // filer nodes which distance is in given tolerance
1896  auto hi = verts_map.get<0>().upper_bound(tol);
1897  verts_map.get<0>().erase(hi, verts_map.end());
1898 
1899  auto remove_verts = [&](Range nodes) {
1900  for (auto n : nodes) {
1901  auto r = verts_map.get<1>().equal_range(n);
1902  verts_map.get<1>().erase(r.first, r.second);
1903  }
1904  };
1905 
1906  auto trim_verts = [&](const Range verts, const bool move) {
1907  VertMapMultiIndex verts_map_tmp;
1908  for (auto p = corners.pair_begin(); p != corners.pair_end(); ++p) {
1909  auto lo = verts_map.get<1>().lower_bound(p->first);
1910  auto hi = verts_map.get<1>().upper_bound(p->second);
1911  verts_map_tmp.insert(lo, hi);
1912  }
1913  if (move) {
1914  for (auto &m : verts_map_tmp.get<0>())
1915  add_trim_vert(m.v, m.e);
1916  } else {
1917  for (auto &m : verts_map_tmp.get<0>())
1918  add_no_move_trim(m.v, m.e);
1919  }
1920  };
1921 
1922  auto trim_edges = [&](const Range ents, const bool move) {
1923  VertMapMultiIndex verts_map_tmp;
1924  for (auto p = ents.pair_begin(); p != ents.pair_end(); ++p) {
1925  auto lo = verts_map.get<2>().lower_bound(p->first);
1926  auto hi = verts_map.get<2>().upper_bound(p->second);
1927  verts_map_tmp.insert(lo, hi);
1928  verts_map.get<2>().erase(lo, hi);
1929  }
1930  if (move) {
1931  for (auto &m : verts_map_tmp.get<0>())
1932  add_trim_vert(m.v, m.e);
1933  } else {
1934  for (auto &m : verts_map_tmp.get<0>())
1935  add_no_move_trim(m.v, m.e);
1936  }
1937  };
1938 
1939  auto intersect_self = [&](Range &a, const Range b) { a = intersect(a, b); };
1940 
1941  map<std::string, Range> range_maps;
1942  CHKERR skin.find_skin(0, cutNewSurfaces, false, range_maps["surface_skin"]);
1943  intersect_self(range_maps["surface_skin"], trimEdges);
1944  range_maps["fixed_edges_on_surface_skin"] =
1945  intersect(range_maps["surface_skin"], fix_edges);
1946  CHKERR moab.get_adjacencies(range_maps["fixed_edges_verts"], 1, false,
1947  range_maps["fixed_edges_verts_edges"],
1948  moab::Interface::UNION);
1949  intersect_self(range_maps["fixed_edges_verts_edges"], trimEdges);
1950  CHKERR moab.get_connectivity(
1951  range_maps["fixed_edges_verts_edges"],
1952  range_maps["fixed_edges_verts_edges_verts_on_the_skin"], false);
1953  intersect_self(range_maps["fixed_edges_verts_edges_verts_on_the_skin"],
1954  tets_skin_verts);
1955 
1956  // do not move nodes at the corners
1957  trim_verts(corners, false);
1958  remove_verts(corners);
1959  trim_edges(range_maps["fixed_edges_on_surface_skin"], true);
1960  remove_verts(range_maps["fixed_edges_on_surface_skin_verts"]);
1961  trim_verts(range_maps["fixed_edges_verts_edges_verts_on_the_skin"], false);
1962  remove_verts(range_maps["fixed_edges_verts_edges_verts_on_the_skin"]);
1963  trim_edges(range_maps["surface_skin"], true);
1964  trim_verts(tets_skin_verts, false);
1965  remove_verts(tets_skin_verts);
1966 
1967  for (auto &m : verts_map.get<0>())
1968  add_trim_vert(m.v, m.e);
1969 
1970  for (auto v : subtract(cut_surface_verts, checked_verts)) {
1971 
1972  if (!(trimNewVertices.find(v) != trimNewVertices.end())) {
1973 
1974  auto get_tag_data = [&](auto th, auto conn) {
1976  CHKERR moab.tag_get_data(th, &conn, 1, &t(0));
1977  return t;
1978  };
1979 
1980  auto t_dist = get_tag_data(th_dist_front_vec, v);
1981  const double d = sqrt(t_dist(i) * t_dist(i));
1982  if (d < tol * max_edge_length) {
1983 
1984  Range adj;
1985  CHKERR moab.get_adjacencies(&v, 1, 1, false, adj);
1986  adj = intersect(adj, cut_surface_edges);
1987  double min_length = max_edge_length;
1988  for (auto e : adj) {
1989 
1990  // Get edge connectivity and coords
1991  const EntityHandle *conn_edge;
1992  CHKERR moab.get_connectivity(e, conn_edge, num_nodes, true);
1993  double coords_edge[3 * num_nodes];
1994  CHKERR moab.get_coords(conn_edge, num_nodes, coords_edge);
1995  FTensor::Tensor1<double *, 3> t_e0(&coords_edge[0], &coords_edge[1],
1996  &coords_edge[2]);
1997  FTensor::Tensor1<double *, 3> t_e1(&coords_edge[3], &coords_edge[4],
1998  &coords_edge[5]);
1999  FTensor::Tensor1<double, 3> t_edge_delta;
2000  t_edge_delta(i) = t_e1(i) - t_e0(i);
2001 
2002  const double length = sqrt(t_edge_delta(i) * t_edge_delta(i));
2003  min_length = std::min(min_length, length);
2004  }
2005 
2006  if (d < (tol * tol * min_length)) {
2007  verticesOnTrimEdges[v].dIst = 0;
2008  verticesOnTrimEdges[v].lEngth = 0;
2009  verticesOnTrimEdges[v].unitRayDir.resize(3);
2010  verticesOnTrimEdges[v].unitRayDir.clear();
2011  verticesOnTrimEdges[v].rayPoint = get_point_coords(v);
2012  trimNewVertices.insert(v);
2013  }
2014  }
2015  }
2016  }
2017 
2018  for (auto m : verticesOnTrimEdges) {
2019  EntityHandle v = m.first;
2020  Range adj;
2021  CHKERR moab.get_adjacencies(&v, 1, 1, false, adj);
2022  for (auto e : adj) {
2023  edgesToTrim.erase(e);
2024  trimEdges.erase(e);
2025  }
2026  }
2027 
2028  if (debug)
2029  if (!trimNewVertices.empty())
2030  CHKERR SaveData(moab)("trim_close_vertices.vtk", trimNewVertices);
2031 
2032  if (debug)
2033  if (!trimEdges.empty())
2034  CHKERR SaveData(moab)("trim_edges.vtk", trimEdges);
2035 
2037 }
2038 
2040  const bool debug) {
2041  CoreInterface &m_field = cOre;
2042  moab::Interface &moab = m_field.get_moab();
2043  MeshRefinement *refiner;
2044  const RefEntity_multiIndex *ref_ents_ptr;
2046 
2047  CHKERR m_field.getInterface(refiner);
2048  CHKERR m_field.get_ref_ents(&ref_ents_ptr);
2050  CHKERR refiner->refine_TET(cutNewVolumes, bit, false, QUIET,
2051  debug ? &trimEdges : NULL);
2052 
2053  trimNewVolumes.clear();
2054  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2055  bit, BitRefLevel().set(), MBTET, trimNewVolumes);
2056 
2057  for (map<EntityHandle, TreeData>::iterator mit = edgesToTrim.begin();
2058  mit != edgesToTrim.end(); mit++) {
2059  auto vit = ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().find(
2060  boost::make_tuple(MBVERTEX, mit->first));
2061  if (vit ==
2062  ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().end())
2063  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
2064  "No vertex on trim edges, that make no sense");
2065 
2066  const boost::shared_ptr<RefEntity> &ref_ent = *vit;
2067  if ((ref_ent->getBitRefLevel() & bit).any()) {
2068  EntityHandle vert = ref_ent->getRefEnt();
2069  trimNewVertices.insert(vert);
2070  verticesOnTrimEdges[vert] = mit->second;
2071  }
2072  }
2073 
2074  // Get faces which are trimmed
2075  trimNewSurfaces.clear();
2076  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2077  bit, bit, MBTRI, trimNewSurfaces);
2078 
2079  Range trim_new_surfaces_nodes;
2080  CHKERR moab.get_connectivity(trimNewSurfaces, trim_new_surfaces_nodes, true);
2081  trim_new_surfaces_nodes = subtract(trim_new_surfaces_nodes, trimNewVertices);
2082  trim_new_surfaces_nodes = subtract(trim_new_surfaces_nodes, cutNewVertices);
2083  Range faces_not_on_surface;
2084  CHKERR moab.get_adjacencies(trim_new_surfaces_nodes, 2, false,
2085  faces_not_on_surface, moab::Interface::UNION);
2086  trimNewSurfaces = subtract(trimNewSurfaces, faces_not_on_surface);
2087 
2088  // Get surfaces which are not trimmed and add them to surface
2089  Range all_surfaces_on_bit_level;
2090  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2091  bit, BitRefLevel().set(), MBTRI, all_surfaces_on_bit_level);
2092  all_surfaces_on_bit_level =
2093  intersect(all_surfaces_on_bit_level, cutNewSurfaces);
2094 
2095  trimNewSurfaces = unite(trimNewSurfaces, all_surfaces_on_bit_level);
2096 
2098 }
2099 
2101  Range *corner_nodes,
2102  const bool debug) {
2103  CoreInterface &m_field = cOre;
2104  moab::Interface &moab = m_field.get_moab();
2106 
2107  Skinner skin(&moab);
2108  Range trim_tets_skin;
2109  CHKERR skin.find_skin(0, trimNewVolumes, false, trim_tets_skin);
2110  Range trim_tets_skin_edges;
2111  CHKERR moab.get_adjacencies(trim_tets_skin, 1, false, trim_tets_skin_edges,
2112  moab::Interface::UNION);
2113 
2114  Range barrier_vertices(trimNewVertices);
2115 
2116  if (fixed_edges && trimFixedEdges) {
2117 
2118  // get all vertices on fixed edges and surface
2119  Range trim_surface_edges;
2120  CHKERR moab.get_adjacencies(trimNewSurfaces, 1, false, trim_surface_edges,
2121  moab::Interface::UNION);
2122  Range fixed_edges_on_trim_surface;
2123  fixed_edges_on_trim_surface = intersect(trim_surface_edges, *fixed_edges);
2124  Range fixed_edges_on_trim_surface_verts;
2125  CHKERR moab.get_connectivity(fixed_edges_on_trim_surface,
2126  fixed_edges_on_trim_surface_verts, false);
2127 
2128  // get faces adjacent to barrier_vertices
2129  Range barrier_vertices_faces;
2130  CHKERR moab.get_adjacencies(barrier_vertices, 2, false,
2131  barrier_vertices_faces, moab::Interface::UNION);
2132  barrier_vertices_faces = intersect(barrier_vertices_faces, trimNewSurfaces);
2133 
2134  // get vertices on fixed edges
2135  Range fixed_edges_vertices;
2136  CHKERR moab.get_connectivity(*fixed_edges, fixed_edges_vertices, false);
2137  fixed_edges_vertices = intersect(barrier_vertices, fixed_edges_vertices);
2138  fixed_edges_vertices =
2139  subtract(fixed_edges_vertices, fixed_edges_on_trim_surface_verts);
2140  if (corner_nodes)
2141  fixed_edges_vertices.merge(intersect(barrier_vertices, *corner_nodes));
2142 
2143  // get faces adjacent to vertices on fixed edges
2144  Range fixed_edges_faces;
2145  CHKERR moab.get_adjacencies(fixed_edges_vertices, 2, false,
2146  fixed_edges_faces, moab::Interface::UNION);
2147  fixed_edges_faces = intersect(fixed_edges_faces, barrier_vertices_faces);
2148 
2149  if (debug && !fixed_edges_faces.empty())
2150  CHKERR SaveData(m_field.get_moab())("fixed_edges_faces.vtk",
2151  fixed_edges_faces);
2152 
2153  // get nodes on faces
2154  Range fixed_edges_faces_vertices;
2155  CHKERR moab.get_connectivity(fixed_edges_faces, fixed_edges_faces_vertices,
2156  false);
2157  barrier_vertices.merge(fixed_edges_faces_vertices);
2158  }
2159 
2160  auto remove_faces_on_skin = [&]() {
2162  Range skin_faces = intersect(trimNewSurfaces, trim_tets_skin);
2163  if (!skin_faces.empty()) {
2164  Range add_to_barrier;
2165  CHKERR moab.get_connectivity(skin_faces, add_to_barrier, false);
2166  barrier_vertices.merge(add_to_barrier);
2167  for (auto f : skin_faces)
2168  trimNewSurfaces.erase(f);
2169  }
2171  };
2172 
2173  auto get_trim_free_edges = [&]() {
2174  // get current surface skin
2175  Range trim_surf_skin;
2176  CHKERR skin.find_skin(0, trimNewSurfaces, false, trim_surf_skin);
2177  trim_surf_skin = subtract(trim_surf_skin, trim_tets_skin_edges);
2178 
2179  Range trim_surf_skin_verts;
2180  CHKERR moab.get_connectivity(trim_surf_skin, trim_surf_skin_verts, true);
2181  for (auto e : barrier_vertices)
2182  trim_surf_skin_verts.erase(e);
2183 
2184  Range free_edges;
2185  CHKERR moab.get_adjacencies(trim_surf_skin_verts, 1, false, free_edges,
2186  moab::Interface::UNION);
2187  free_edges = intersect(free_edges, trim_surf_skin);
2188 
2189  return free_edges;
2190  };
2191 
2192  CHKERR remove_faces_on_skin();
2193 
2194  if (debug && !barrier_vertices.empty())
2195  CHKERR SaveData(m_field.get_moab())("barrier_vertices.vtk",
2196  barrier_vertices);
2197 
2198  int nn = 0;
2199  Range out_edges;
2200  while (!(out_edges = get_trim_free_edges()).empty()) {
2201 
2202  if (debug && !trimNewSurfaces.empty())
2203  CHKERR SaveData(m_field.get_moab())(
2204  "trimNewSurfaces_" + boost::lexical_cast<std::string>(nn) + ".vtk",
2205  trimNewSurfaces);
2206 
2207  if (debug && !out_edges.empty())
2208  CHKERR SaveData(m_field.get_moab())(
2209  "trimNewSurfacesOutsideVerts_" +
2210  boost::lexical_cast<std::string>(nn) + ".vtk",
2211  out_edges);
2212 
2213  Range outside_faces;
2214  CHKERR moab.get_adjacencies(out_edges, 2, false, outside_faces,
2215  moab::Interface::UNION);
2216  trimNewSurfaces = subtract(trimNewSurfaces, outside_faces);
2217  ++nn;
2218  }
2219 
2220  if (debug && !trimNewSurfaces.empty())
2221  CHKERR SaveData(m_field.get_moab())(
2222  "trimNewSurfaces_" + boost::lexical_cast<std::string>(nn) + ".vtk",
2223  trimNewSurfaces);
2224 
2226 }
2227 
2230  Range &ents) {
2231  CoreInterface &m_field = cOre;
2232  moab::Interface &moab = m_field.get_moab();
2233  PrismInterface *interface;
2235  CHKERR m_field.getInterface(interface);
2236  // Remove tris on surface front
2237  {
2238  Range front_tris;
2239  EntityHandle meshset;
2240  CHKERR moab.create_meshset(MESHSET_SET, meshset);
2241  CHKERR moab.add_entities(meshset, ents);
2243  meshset, split_bit, true, front_tris);
2244  CHKERR moab.delete_entities(&meshset, 1);
2245  ents = subtract(ents, front_tris);
2246  }
2247  // Remove entities on skin
2248  Range tets;
2249  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2250  split_bit, BitRefLevel().set(), MBTET, tets);
2251  // Remove entities on skin
2252  Skinner skin(&moab);
2253  Range tets_skin;
2254  CHKERR skin.find_skin(0, tets, false, tets_skin);
2255  ents = subtract(ents, tets_skin);
2256 
2258 }
2259 
2261  const BitRefLevel bit,
2262  const Range &ents, Tag th) {
2263  CoreInterface &m_field = cOre;
2264  moab::Interface &moab = m_field.get_moab();
2265  PrismInterface *interface;
2267  CHKERR m_field.getInterface(interface);
2268  EntityHandle meshset_volume;
2269  CHKERR moab.create_meshset(MESHSET_SET, meshset_volume);
2270  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2271  split_bit, BitRefLevel().set(), MBTET, meshset_volume);
2272  EntityHandle meshset_surface;
2273  CHKERR moab.create_meshset(MESHSET_SET, meshset_surface);
2274  CHKERR moab.add_entities(meshset_surface, ents);
2275  CHKERR interface->getSides(meshset_surface, split_bit, true);
2276  CHKERR interface->splitSides(meshset_volume, bit, meshset_surface, true,
2277  true);
2278  CHKERR moab.delete_entities(&meshset_volume, 1);
2279  CHKERR moab.delete_entities(&meshset_surface, 1);
2280  if (th) {
2281  Range prisms;
2282  ierr = m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2283  bit, BitRefLevel().set(), MBPRISM, prisms);
2284  for (Range::iterator pit = prisms.begin(); pit != prisms.end(); pit++) {
2285  int num_nodes;
2286  const EntityHandle *conn;
2287  CHKERR moab.get_connectivity(*pit, conn, num_nodes, true);
2288  MatrixDouble data(3, 3);
2289  CHKERR moab.tag_get_data(th, conn, 3, &data(0, 0));
2290  CHKERR moab.tag_set_data(th, &conn[3], 3, &data(0, 0));
2291  }
2292  }
2294 }
2295 
2297  const int fraction_level, const Range &tets, const Range &surface,
2298  const Range &fixed_edges, const Range &corner_nodes, Range &edges_to_merge,
2299  Range &out_tets, Range &new_surf, Tag th, const bool update_meshsets,
2300  const BitRefLevel *bit_ptr, const bool debug) {
2301  CoreInterface &m_field = cOre;
2302  moab::Interface &moab = m_field.get_moab();
2304 
2305  /**
2306  * \brief Merge nodes
2307  */
2308  struct MergeNodes {
2309  CoreInterface &mField;
2310  const bool onlyIfImproveQuality;
2311  Tag tH;
2312  bool updateMehsets;
2313 
2314  MergeNodes(CoreInterface &m_field, const bool only_if_improve_quality,
2315  Tag th, bool update_mehsets)
2316  : mField(m_field), onlyIfImproveQuality(only_if_improve_quality),
2317  tH(th), updateMehsets(update_mehsets) {
2318  mField.getInterface(nodeMergerPtr);
2319  }
2320  NodeMergerInterface *nodeMergerPtr;
2321  MoFEMErrorCode mergeNodes(int line_search, EntityHandle father,
2322  EntityHandle mother, Range &proc_tets,
2323  bool add_child = true) {
2324  moab::Interface &moab = mField.get_moab();
2326  const EntityHandle conn[] = {father, mother};
2327  Range vert_tets;
2328  CHKERR moab.get_adjacencies(conn, 2, 3, false, vert_tets,
2329  moab::Interface::UNION);
2330  vert_tets = intersect(vert_tets, proc_tets);
2331  Range out_tets;
2332  CHKERR nodeMergerPtr->mergeNodes(father, mother, out_tets, &vert_tets,
2333  onlyIfImproveQuality, 0, line_search,
2334  tH);
2335 
2336  if (add_child && nodeMergerPtr->getSuccessMerge()) {
2337 
2338  Range::iterator lo, hi = proc_tets.begin();
2339  for (auto pt = vert_tets.pair_begin(); pt != vert_tets.pair_end();
2340  ++pt) {
2341  lo = proc_tets.lower_bound(hi, proc_tets.end(), pt->first);
2342  if (lo != proc_tets.end()) {
2343  hi = proc_tets.upper_bound(lo, proc_tets.end(), pt->second);
2344  proc_tets.erase(lo, hi);
2345  } else
2346  break;
2347  }
2348  proc_tets.merge(out_tets);
2349 
2350  auto &parent_child_map = nodeMergerPtr->getParentChildMap();
2351 
2352  struct ChangeChild {
2353  EntityHandle child;
2354  ChangeChild(const EntityHandle child) : child(child) {}
2355  void operator()(NodeMergerInterface::ParentChild &p) {
2356  p.cHild = child;
2357  }
2358  };
2359 
2360  std::vector<decltype(parentsChildMap.get<0>().begin())> it_vec;
2361  it_vec.reserve(parentsChildMap.size());
2362 
2363  for (auto &p : parent_child_map) {
2364 
2365  it_vec.clear();
2366  for (auto it = parentsChildMap.get<0>().equal_range(p.pArent);
2367  it.first != it.second; ++it.first)
2368  it_vec.emplace_back(it.first);
2369 
2370  for (auto it = parentsChildMap.get<1>().equal_range(p.pArent);
2371  it.first != it.second; ++it.first)
2372  it_vec.emplace_back(parentsChildMap.project<0>(it.first));
2373 
2374  for (auto &it : it_vec)
2375  parentsChildMap.modify(it, ChangeChild(p.cHild));
2376  }
2377 
2378  parentsChildMap.insert(parent_child_map.begin(),
2379  parent_child_map.end());
2380  }
2382  }
2383 
2384  MoFEMErrorCode updateRangeByChilds(Range &new_surf, Range &edges_to_merge,
2385  Range &not_merged_edges,
2386  bool add_child) {
2387  moab::Interface &moab = mField.get_moab();
2389  if (add_child) {
2390 
2391  std::vector<EntityHandle> parents_ents_vec(parentsChildMap.size());
2392  for (auto &it : parentsChildMap)
2393  parents_ents_vec.emplace_back(it.pArent);
2394  Range parent_ents;
2395  parent_ents.insert_list(parents_ents_vec.begin(),
2396  parents_ents_vec.end());
2397 
2398  Range surf_parent_ents = intersect(new_surf, parent_ents);
2399  new_surf = subtract(new_surf, surf_parent_ents);
2400  Range child_surf_ents;
2401  CHKERR updateRangeByChilds(parentsChildMap, surf_parent_ents,
2402  child_surf_ents);
2403  new_surf.merge(child_surf_ents);
2404 
2405  Range edges_to_merge_parent_ents =
2406  intersect(edges_to_merge, parent_ents);
2407  edges_to_merge = subtract(edges_to_merge, edges_to_merge_parent_ents);
2408  Range merged_child_edge_ents;
2409  CHKERR updateRangeByChilds(parentsChildMap, edges_to_merge_parent_ents,
2410  merged_child_edge_ents);
2411 
2412  Range not_merged_edges_child_ents =
2413  intersect(not_merged_edges, parent_ents);
2414  not_merged_edges =
2415  subtract(not_merged_edges, not_merged_edges_child_ents);
2416  Range not_merged_child_edge_ents;
2417  CHKERR updateRangeByChilds(parentsChildMap, not_merged_edges_child_ents,
2418  not_merged_child_edge_ents);
2419 
2420  merged_child_edge_ents =
2421  subtract(merged_child_edge_ents, not_merged_child_edge_ents);
2422  edges_to_merge.merge(merged_child_edge_ents);
2423  not_merged_edges.merge(not_merged_child_edge_ents);
2424 
2425  if (updateMehsets) {
2427  (*mField.getInterface<MeshsetsManager>()), cubit_it)) {
2428  EntityHandle cubit_meshset = cubit_it->meshset;
2429  Range meshset_ents;
2430  CHKERR moab.get_entities_by_handle(cubit_meshset, meshset_ents,
2431  true);
2432  Range child_ents;
2433  CHKERR updateRangeByChilds(parentsChildMap, meshset_ents,
2434  child_ents);
2435  CHKERR moab.add_entities(cubit_meshset, child_ents);
2436  }
2437  }
2438  }
2439 
2441  };
2442 
2443  private:
2444  NodeMergerInterface::ParentChildMap parentsChildMap;
2445  std::vector<EntityHandle> childsVec;
2446 
2447  inline MoFEMErrorCode updateRangeByChilds(
2448  const NodeMergerInterface::ParentChildMap &parent_child_map,
2449  const Range &parents, Range &childs) {
2451  childsVec.clear();
2452  childsVec.reserve(parents.size());
2453  for (auto pit = parents.pair_begin(); pit != parents.pair_end(); pit++) {
2454  auto it = parent_child_map.lower_bound(pit->first);
2455  if (it != parent_child_map.end()) {
2456  for (auto hi_it = parent_child_map.upper_bound(pit->second);
2457  it != hi_it; ++it)
2458  childsVec.emplace_back(it->cHild);
2459  }
2460  }
2461  childs.insert_list(childsVec.begin(), childsVec.end());
2463  }
2464  };
2465 
2466  /**
2467  * \brief Calculate edge length
2468  */
2469  struct LengthMap {
2470  Tag tH;
2471  CoreInterface &mField;
2472  moab::Interface &moab;
2473  const double maxLength;
2474  LengthMap(CoreInterface &m_field, Tag th, double max_length)
2475  : tH(th), mField(m_field), moab(m_field.get_moab()),
2476  maxLength(max_length) {
2477  gammaL = 1.;
2478  gammaQ = 3.;
2479  acceptedThrasholdMergeQuality = 0.5;
2480  }
2481 
2482  double
2483  gammaL; ///< Controls importance of length when ranking edges for merge
2484  double
2485  gammaQ; ///< Controls importance of quality when ranking edges for merge
2486  double acceptedThrasholdMergeQuality; ///< Do not merge quality if quality
2487  ///< above accepted thrashold
2488 
2489  MoFEMErrorCode operator()(const Range &tets, const Range &edges,
2490  LengthMapData_multi_index &length_map,
2491  double &ave) const {
2492  int num_nodes;
2493  const EntityHandle *conn;
2494  std::array<double, 12> coords;
2496  VectorAdaptor s0(3, ublas::shallow_array_adaptor<double>(3, &coords[0]));
2497  VectorAdaptor s1(3, ublas::shallow_array_adaptor<double>(3, &coords[3]));
2498  VectorDouble3 delta(3);
2499 
2500  struct NodeQuality {
2501  EntityHandle node;
2502  double quality;
2503  NodeQuality(const EntityHandle node) : node(node), quality(1) {}
2504  };
2505 
2506  typedef multi_index_container<
2507  NodeQuality, indexed_by<
2508 
2509  sequenced<>,
2510 
2511  hashed_non_unique<tag<Ent_mi_tag>,
2512  member<NodeQuality, EntityHandle,
2513  &NodeQuality::node>>
2514 
2515  >>
2516  NodeQuality_sequence;
2517 
2518  NodeQuality_sequence node_quality_sequence;
2519 
2520  Range edges_nodes;
2521  CHKERR moab.get_connectivity(tets, edges_nodes, false);
2522  Range edges_tets;
2523  CHKERR moab.get_adjacencies(edges, 3, false, edges_tets,
2524  moab::Interface::UNION);
2525  edges_tets = intersect(edges_tets, tets);
2526 
2527  for (auto node : edges_nodes)
2528  node_quality_sequence.emplace_back(node);
2529 
2530  for (auto tet : edges_tets) {
2531 
2532  CHKERR moab.get_connectivity(tet, conn, num_nodes, true);
2533  if (tH)
2534  CHKERR moab.tag_get_data(tH, conn, num_nodes, coords.data());
2535  else
2536  CHKERR moab.get_coords(conn, num_nodes, coords.data());
2537 
2538  const double q = Tools::volumeLengthQuality(coords.data());
2539 
2540  for (auto n : {0, 1, 2, 3}) {
2541  auto it = node_quality_sequence.get<1>().find(conn[n]);
2542  if (it != node_quality_sequence.get<1>().end())
2543  const_cast<double &>(it->quality) = std::min(q, it->quality);
2544  }
2545  }
2546 
2547  for (auto edge : edges) {
2548 
2549  CHKERR moab.get_connectivity(edge, conn, num_nodes, true);
2550 
2551  if (tH)
2552  CHKERR moab.tag_get_data(tH, conn, num_nodes, coords.data());
2553  else
2554  CHKERR moab.get_coords(conn, num_nodes, coords.data());
2555 
2556  double q = 1;
2557  for (auto n : {0, 1}) {
2558  auto it = node_quality_sequence.get<1>().find(conn[n]);
2559  if (it != node_quality_sequence.get<1>().end())
2560  q = std::min(q, it->quality);
2561  }
2562 
2563  if (q < acceptedThrasholdMergeQuality) {
2564  noalias(delta) = (s0 - s1) / maxLength;
2565  double dot = inner_prod(delta, delta);
2566  double val = pow(q, gammaQ) * pow(dot, 0.5 * gammaL);
2567  length_map.insert(LengthMapData(val, q, edge));
2568  }
2569  }
2570 
2571  ave = 0;
2572  for (LengthMapData_multi_index::nth_index<0>::type::iterator mit =
2573  length_map.get<0>().begin();
2574  mit != length_map.get<0>().end(); mit++) {
2575  ave += mit->qUality;
2576  }
2577  ave /= length_map.size();
2579  }
2580  };
2581 
2582  /**
2583  * \brief Topological relations
2584  */
2585  struct Toplogy {
2586 
2587  CoreInterface &mField;
2588  Tag tH;
2589  const double tOL;
2590  Toplogy(CoreInterface &m_field, Tag th, const double tol)
2591  : mField(m_field), tH(th), tOL(tol) {}
2592 
2593  enum TYPE {
2594  FREE = 0,
2595  SKIN = 1 << 0,
2596  SURFACE = 1 << 1,
2597  SURFACE_SKIN = 1 << 2,
2598  FRONT_ENDS = 1 << 3,
2599  FIX_EDGES = 1 << 4,
2600  FIX_CORNERS = 1 << 5
2601  };
2602 
2603  typedef map<int, Range> SetsMap;
2604 
2605  MoFEMErrorCode classifyVerts(const Range &surface, const Range &tets,
2606  const Range &fixed_edges,
2607  const Range &corner_nodes,
2608  SetsMap &sets_map) const {
2609  moab::Interface &moab(mField.get_moab());
2610  Skinner skin(&moab);
2612 
2613  sets_map[FIX_CORNERS].merge(corner_nodes);
2614  Range fixed_verts;
2615  CHKERR moab.get_connectivity(fixed_edges, fixed_verts, true);
2616  sets_map[FIX_EDGES].swap(fixed_verts);
2617 
2618  Range tets_skin;
2619  CHKERR skin.find_skin(0, tets, false, tets_skin);
2620  Range tets_skin_edges;
2621  CHKERR moab.get_adjacencies(tets_skin, 1, false, tets_skin_edges,
2622  moab::Interface::UNION);
2623 
2624  // surface skin
2625  Range surface_skin;
2626  CHKERR skin.find_skin(0, surface, false, surface_skin);
2627  Range front_in_the_body;
2628  front_in_the_body = subtract(surface_skin, tets_skin_edges);
2629  Range front_ends;
2630  CHKERR skin.find_skin(0, front_in_the_body, false, front_ends);
2631  sets_map[FRONT_ENDS].swap(front_ends);
2632 
2633  Range surface_skin_verts;
2634  CHKERR moab.get_connectivity(surface_skin, surface_skin_verts, true);
2635  sets_map[SURFACE_SKIN].swap(surface_skin_verts);
2636 
2637  // surface
2638  Range surface_verts;
2639  CHKERR moab.get_connectivity(surface, surface_verts, true);
2640  sets_map[SURFACE].swap(surface_verts);
2641 
2642  // skin
2643  Range tets_skin_verts;
2644  CHKERR moab.get_connectivity(tets_skin, tets_skin_verts, true);
2645  sets_map[SKIN].swap(tets_skin_verts);
2646 
2647  Range tets_verts;
2648  CHKERR moab.get_connectivity(tets, tets_verts, true);
2649  sets_map[FREE].swap(tets_verts);
2650 
2652  }
2653 
2654  MoFEMErrorCode getProcTets(const Range &tets, const Range &edges_to_merge,
2655  Range &proc_tets) const {
2656  moab::Interface &moab(mField.get_moab());
2658  Range edges_to_merge_verts;
2659  CHKERR moab.get_connectivity(edges_to_merge, edges_to_merge_verts, true);
2660  Range edges_to_merge_verts_tets;
2661  CHKERR moab.get_adjacencies(edges_to_merge_verts, 3, false,
2662  edges_to_merge_verts_tets,
2663  moab::Interface::UNION);
2664  edges_to_merge_verts_tets = intersect(edges_to_merge_verts_tets, tets);
2665  proc_tets.swap(edges_to_merge_verts_tets);
2667  }
2668 
2669  MoFEMErrorCode removeBadEdges(const Range &surface, const Range &tets,
2670  const Range &fixed_edges,
2671  const Range &corner_nodes,
2672  Range &edges_to_merge,
2673  Range &not_merged_edges) {
2674  moab::Interface &moab(mField.get_moab());
2676 
2677  // find skin
2678  Skinner skin(&moab);
2679  Range tets_skin;
2680  CHKERR skin.find_skin(0, tets, false, tets_skin);
2681  Range surface_skin;
2682  CHKERR skin.find_skin(0, surface, false, surface_skin);
2683 
2684  // end nodes
2685  Range tets_skin_edges;
2686  CHKERR moab.get_adjacencies(tets_skin, 1, false, tets_skin_edges,
2687  moab::Interface::UNION);
2688  Range surface_front;
2689  surface_front = subtract(surface_skin, tets_skin_edges);
2690  Range surface_front_nodes;
2691  CHKERR moab.get_connectivity(surface_front, surface_front_nodes, true);
2692  Range ends_nodes;
2693  CHKERR skin.find_skin(0, surface_front, false, ends_nodes);
2694 
2695  // remove bad merges
2696 
2697  // get surface and body skin verts
2698  Range surface_edges;
2699  CHKERR moab.get_adjacencies(surface, 1, false, surface_edges,
2700  moab::Interface::UNION);
2701  // get nodes on the surface
2702  Range surface_edges_verts;
2703  CHKERR moab.get_connectivity(surface_edges, surface_edges_verts, true);
2704  // get vertices on the body skin
2705  Range tets_skin_edges_verts;
2706  CHKERR moab.get_connectivity(tets_skin_edges, tets_skin_edges_verts,
2707  true);
2708 
2709  Range edges_to_remove;
2710 
2711  // remove edges self connected to body skin
2712  {
2713  Range ents_nodes_and_edges;
2714  ents_nodes_and_edges.merge(tets_skin_edges_verts);
2715  ents_nodes_and_edges.merge(tets_skin_edges);
2716  CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2717  0, false);
2718  }
2719  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2720  not_merged_edges.merge(edges_to_remove);
2721 
2722  // remove edges self connected to surface
2723  {
2724  Range ents_nodes_and_edges;
2725  ents_nodes_and_edges.merge(surface_edges_verts);
2726  ents_nodes_and_edges.merge(surface_edges);
2727  ents_nodes_and_edges.merge(tets_skin_edges_verts);
2728  ents_nodes_and_edges.merge(tets_skin_edges);
2729  CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2730  0, false);
2731  }
2732  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2733  not_merged_edges.merge(edges_to_remove);
2734 
2735  // remove edges adjacent corner_nodes execpt those on fixed edges
2736  Range fixed_edges_nodes;
2737  CHKERR moab.get_connectivity(fixed_edges, fixed_edges_nodes, true);
2738  {
2739  Range ents_nodes_and_edges;
2740  ents_nodes_and_edges.merge(fixed_edges_nodes);
2741  ents_nodes_and_edges.merge(ends_nodes);
2742  ents_nodes_and_edges.merge(corner_nodes);
2743  ents_nodes_and_edges.merge(fixed_edges);
2744  CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2745  0, false);
2746  }
2747  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2748  not_merged_edges.merge(edges_to_remove);
2749 
2750  // remove edges self connected to surface
2751  CHKERR removeSelfConectingEdges(surface_edges, edges_to_remove, 0, false);
2752  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2753  not_merged_edges.merge(edges_to_remove);
2754 
2755  // remove edges self contented on surface skin
2756  {
2757  Range ents_nodes_and_edges;
2758  ents_nodes_and_edges.merge(surface_skin);
2759  ents_nodes_and_edges.merge(fixed_edges_nodes);
2760  CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2761  0, false);
2762  }
2763  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2764  not_merged_edges.merge(edges_to_remove);
2765 
2766  // remove edges connecting crack front and fixed edges, except those short
2767  {
2768  Range ents_nodes_and_edges;
2769  ents_nodes_and_edges.merge(surface_skin.subset_by_type(MBEDGE));
2770  ents_nodes_and_edges.merge(fixed_edges.subset_by_type(MBEDGE));
2771  CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2772  0, false);
2773  }
2774  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2775  not_merged_edges.merge(edges_to_remove);
2776 
2777  // remove crack front nodes connected to the surface, except those short
2778  {
2779  Range ents_nodes_and_edges;
2780  ents_nodes_and_edges.merge(surface_front_nodes);
2781  ents_nodes_and_edges.merge(surface_front);
2782  ents_nodes_and_edges.merge(tets_skin_edges_verts);
2783  ents_nodes_and_edges.merge(tets_skin_edges);
2784  CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2785  tOL, false);
2786  }
2787  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2788  not_merged_edges.merge(edges_to_remove);
2789 
2791  }
2792 
2793  private:
2794  MoFEMErrorCode removeSelfConectingEdges(const Range &ents,
2795  Range &edges_to_remove,
2796  const bool length,
2797  bool debug) const {
2798  moab::Interface &moab(mField.get_moab());
2800  // get nodes
2801  Range ents_nodes = ents.subset_by_type(MBVERTEX);
2802  if (ents_nodes.empty()) {
2803  CHKERR moab.get_connectivity(ents, ents_nodes, true);
2804  }
2805  // edges adj. to nodes
2806  Range ents_nodes_edges;
2807  CHKERR moab.get_adjacencies(ents_nodes, 1, false, ents_nodes_edges,
2808  moab::Interface::UNION);
2809  // nodes of adj. edges
2810  Range ents_nodes_edges_nodes;
2811  CHKERR moab.get_connectivity(ents_nodes_edges, ents_nodes_edges_nodes,
2812  true);
2813  // hanging nodes
2814  ents_nodes_edges_nodes = subtract(ents_nodes_edges_nodes, ents_nodes);
2815  Range ents_nodes_edges_nodes_edges;
2816  CHKERR moab.get_adjacencies(ents_nodes_edges_nodes, 1, false,
2817  ents_nodes_edges_nodes_edges,
2818  moab::Interface::UNION);
2819  // remove edges adj. to hanging edges
2820  ents_nodes_edges =
2821  subtract(ents_nodes_edges, ents_nodes_edges_nodes_edges);
2822  ents_nodes_edges =
2823  subtract(ents_nodes_edges, ents.subset_by_type(MBEDGE));
2824  if (length > 0) {
2825  Range::iterator eit = ents_nodes_edges.begin();
2826  for (; eit != ents_nodes_edges.end();) {
2827 
2828  int num_nodes;
2829  const EntityHandle *conn;
2830  CHKERR moab.get_connectivity(*eit, conn, num_nodes, true);
2831  double coords[6];
2832  if (tH)
2833  CHKERR moab.tag_get_data(tH, conn, num_nodes, coords);
2834  else
2835  CHKERR moab.get_coords(conn, num_nodes, coords);
2836 
2837  auto get_edge_length = [coords]() {
2839  &coords[0], &coords[1], &coords[2]);
2842  t_delta(i) = t_coords(i);
2843  ++t_coords;
2844  t_delta(i) -= t_coords(i);
2845  return sqrt(t_delta(i) * t_delta(i));
2846  };
2847 
2848  if (get_edge_length() < tOL) {
2849  eit = ents_nodes_edges.erase(eit);
2850  } else {
2851  eit++;
2852  }
2853  }
2854  }
2855  edges_to_remove.swap(ents_nodes_edges);
2856  if (debug) {
2857  CHKERR SaveData(moab)("edges_to_remove.vtk", edges_to_remove);
2858  }
2860  }
2861  };
2862 
2863  Range not_merged_edges;
2864  const double tol = 1e-3;
2865  CHKERR Toplogy(m_field, th, tol * aveLength)
2866  .removeBadEdges(surface, tets, fixed_edges, corner_nodes, edges_to_merge,
2867  not_merged_edges);
2868  Toplogy::SetsMap sets_map;
2869  CHKERR Toplogy(m_field, th, tol * aveLength)
2870  .classifyVerts(surface, tets, fixed_edges, corner_nodes, sets_map);
2871  if (debug) {
2872  for (Toplogy::SetsMap::reverse_iterator sit = sets_map.rbegin();
2873  sit != sets_map.rend(); sit++) {
2874  std::string name = "classification_verts_" +
2875  boost::lexical_cast<std::string>(sit->first) + ".vtk";
2876  if (!sit->second.empty())
2877  CHKERR SaveData(moab)(name, sit->second);
2878  }
2879  }
2880  Range proc_tets;
2881  CHKERR Toplogy(m_field, th, tol * aveLength)
2882  .getProcTets(tets, edges_to_merge, proc_tets);
2883  out_tets = subtract(tets, proc_tets);
2884 
2885  if (bit_ptr) {
2886  Range all_out_ents = out_tets;
2887  for (int dd = 2; dd >= 0; dd--) {
2888  CHKERR moab.get_adjacencies(out_tets, dd, false, all_out_ents,
2889  moab::Interface::UNION);
2890  }
2891  CHKERR m_field.getInterface<BitRefManager>()->addBitRefLevel(all_out_ents,
2892  *bit_ptr);
2893  }
2894 
2895  int nb_nodes_merged = 0;
2896  LengthMapData_multi_index length_map;
2897  new_surf = surface;
2898 
2899  auto save_merge_step = [&](const int pp, const Range collapsed_edges) {
2901  Range adj_faces;
2902  CHKERR moab.get_adjacencies(proc_tets, 2, false, adj_faces,
2903  moab::Interface::UNION);
2904  std::string name;
2905  name = "node_merge_step_" + boost::lexical_cast<std::string>(pp) + ".vtk";
2906  CHKERR SaveData(moab)(
2907  name, unite(intersect(new_surf, adj_faces), collapsed_edges));
2908  name =
2909  "edges_to_merge_step_" + boost::lexical_cast<std::string>(pp) + ".vtk";
2910  CHKERR SaveData(moab)(
2911  name, unite(intersect(new_surf, adj_faces), edges_to_merge));
2913  };
2914 
2915  if (debug)
2916  CHKERR save_merge_step(0, Range());
2917 
2918  double ave0 = 0, ave = 0, min = 0, min_p = 0, min_pp;
2919  for (int pp = 0; pp != nbMaxMergingCycles; ++pp) {
2920 
2921  int nb_nodes_merged_p = nb_nodes_merged;
2922  length_map.clear();
2923  min_pp = min_p;
2924  min_p = min;
2925  CHKERR LengthMap(m_field, th, aveLength)(proc_tets, edges_to_merge,
2926  length_map, ave);
2927  min = length_map.get<2>().begin()->qUality;
2928  if (pp == 0) {
2929  ave0 = ave;
2930  }
2931 
2932  int nn = 0;
2933  Range collapsed_edges;
2934  MergeNodes merge_nodes(m_field, true, th, update_meshsets);
2935 
2936  for (auto mit = length_map.get<0>().begin();
2937  mit != length_map.get<0>().end(); mit++, nn++) {
2938 
2939  if (!mit->skip) {
2940 
2941  auto get_father_and_mother =
2942  [&](EntityHandle &father, EntityHandle &mother, int &line_search) {
2944  int num_nodes;
2945  const EntityHandle *conn;
2946  CHKERR moab.get_connectivity(mit->eDge, conn, num_nodes, true);
2947  std::array<int, 2> conn_type = {0, 0};
2948  for (int nn = 0; nn != 2; nn++) {
2949  for (Toplogy::SetsMap::reverse_iterator sit = sets_map.rbegin();
2950  sit != sets_map.rend(); sit++) {
2951  if (sit->second.find(conn[nn]) != sit->second.end()) {
2952  conn_type[nn] |= sit->first;
2953  }
2954  }
2955  }
2956  int type_father, type_mother;
2957  if (conn_type[0] > conn_type[1]) {
2958  father = conn[0];
2959  mother = conn[1];
2960  type_father = conn_type[0];
2961  type_mother = conn_type[1];
2962  } else {
2963  father = conn[1];
2964  mother = conn[0];
2965  type_father = conn_type[1];
2966  type_mother = conn_type[0];
2967  }
2968  if (type_father == type_mother) {
2969  line_search = lineSearchSteps;
2970  }
2972  };
2973 
2974  int line_search = 0;
2975  EntityHandle father, mother;
2976  CHKERR get_father_and_mother(father, mother, line_search);
2977  CHKERR merge_nodes.mergeNodes(line_search, father, mother, proc_tets);
2978  if (m_field.getInterface<NodeMergerInterface>()->getSuccessMerge()) {
2979  const EntityHandle father_and_mother[] = {father, mother};
2980  Range adj_tets;
2981  CHKERR moab.get_adjacencies(father_and_mother, 1, 3, false, adj_tets);
2982  Range adj_tets_nodes;
2983  CHKERR moab.get_connectivity(adj_tets, adj_tets_nodes, true);
2984  Range adj_edges;
2985  CHKERR moab.get_adjacencies(adj_tets_nodes, 1, false, adj_edges,
2986  moab::Interface::UNION);
2987  for (auto ait : adj_edges) {
2988  auto miit = length_map.get<1>().find(ait);
2989  if (miit != length_map.get<1>().end())
2990  (const_cast<LengthMapData &>(*miit)).skip = true;
2991  }
2992  nb_nodes_merged++;
2993  collapsed_edges.insert(mit->eDge);
2994  }
2995 
2996  if (nn > static_cast<int>(length_map.size() / fraction_level))
2997  break;
2998  if (mit->qUality > ave)
2999  break;
3000  }
3001  }
3002 
3003  CHKERR merge_nodes.updateRangeByChilds(new_surf, edges_to_merge,
3004  not_merged_edges, true);
3005 
3006  PetscPrintf(m_field.get_comm(),
3007  "(%d) Number of nodes merged %d ave q %3.4e min q %3.4e\n", pp,
3008  nb_nodes_merged, ave, min);
3009 
3010  if (debug)
3011  CHKERR save_merge_step(pp + 1, collapsed_edges);
3012 
3013  if (nb_nodes_merged == nb_nodes_merged_p)
3014  break;
3015  if (min > 1e-2 && min == min_pp)
3016  break;
3017  if (min > ave0)
3018  break;
3019 
3020  Range adj_edges;
3021  CHKERR moab.get_adjacencies(proc_tets, 1, false, adj_edges,
3022  moab::Interface::UNION);
3023  edges_to_merge = intersect(edges_to_merge, adj_edges);
3024  CHKERR Toplogy(m_field, th, tol * aveLength)
3025  .removeBadEdges(new_surf, proc_tets, fixed_edges, corner_nodes,
3026  edges_to_merge, not_merged_edges);
3027  }
3028 
3029  if (bit_ptr)
3030  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevel(proc_tets,
3031  *bit_ptr);
3032 
3033  out_tets.merge(proc_tets);
3034  Range adj_faces;
3035  CHKERR moab.get_adjacencies(out_tets, 2, false, adj_faces,
3036  moab::Interface::UNION);
3037  new_surf = intersect(new_surf, adj_faces);
3038 
3040 }
3041 
3043  const int fraction_level, const BitRefLevel cut_bit,
3044  const BitRefLevel trim_bit, const BitRefLevel bit, const Range &surface,
3045  const Range &fixed_edges, const Range &corner_nodes, Tag th,
3046  const bool update_meshsets, const bool debug) {
3047  CoreInterface &m_field = cOre;
3049  Range tets_level;
3050  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
3051  trim_bit, BitRefLevel().set(), MBTET, tets_level);
3052 
3053  Range edges_to_merge;
3054  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByRefLevel(
3055  trim_bit, cut_bit | trim_bit, edges_to_merge);
3056 
3057  // get all entities not in database
3058  Range all_ents_not_in_database_before;
3059  CHKERR cOre.getInterface<BitRefManager>()->getAllEntitiesNotInDatabase(
3060  all_ents_not_in_database_before);
3061 
3062  edges_to_merge = edges_to_merge.subset_by_type(MBEDGE);
3063  if (debug)
3064  CHKERR SaveData(m_field.get_moab())("edges_to_merge.vtk", edges_to_merge);
3065 
3066  Range out_new_tets, out_new_surf;
3067  CHKERR mergeBadEdges(fraction_level, tets_level, surface, fixed_edges,
3068  corner_nodes, edges_to_merge, out_new_tets, out_new_surf,
3069  th, update_meshsets, &bit, debug);
3070 
3071  // get all entities not in database after merge
3072  Range all_ents_not_in_database_after;
3073  CHKERR cOre.getInterface<BitRefManager>()->getAllEntitiesNotInDatabase(
3074  all_ents_not_in_database_after);
3075 
3076  // delete hanging entities
3077  all_ents_not_in_database_after =
3078  subtract(all_ents_not_in_database_after, all_ents_not_in_database_before);
3079  Range meshsets;
3080  CHKERR m_field.get_moab().get_entities_by_type(0, MBENTITYSET, meshsets,
3081  false);
3082  for (auto m : meshsets)
3083  CHKERR m_field.get_moab().remove_entities(m,
3084  all_ents_not_in_database_after);
3085 
3086  m_field.get_moab().delete_entities(all_ents_not_in_database_after);
3087 
3088  mergedVolumes.swap(out_new_tets);
3089  mergedSurfaces.swap(out_new_surf);
3091 }
3092 
3094  CoreInterface &m_field = cOre;
3095  moab::Interface &moab = m_field.get_moab();
3097  Range nodes;
3098  if (bit.none())
3099  CHKERR moab.get_entities_by_type(0, MBVERTEX, nodes);
3100  else
3101  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
3102  bit, BitRefLevel().set(), MBVERTEX, nodes);
3103  std::vector<double> coords(3 * nodes.size());
3104  CHKERR moab.get_coords(nodes, &coords[0]);
3105  CHKERR moab.tag_set_data(th, nodes, &coords[0]);
3107 }
3108 
3110  const BitRefLevel mask) {
3111  CoreInterface &m_field = cOre;
3112  moab::Interface &moab = m_field.get_moab();
3114  Range nodes;
3115  if (bit.none())
3116  CHKERR moab.get_entities_by_type(0, MBVERTEX, nodes);
3117  else
3118  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
3119  bit, mask, MBVERTEX, nodes);
3120  std::vector<double> coords(3 * nodes.size());
3121  CHKERR moab.tag_get_data(th, nodes, &coords[0]);
3122  CHKERR moab.set_coords(nodes, &coords[0]);
3124 }
3125 
3127  CoreInterface &m_field = cOre;
3128  moab::Interface &moab = m_field.get_moab();
3130  CHKERR SaveData(moab)(prefix + "out_vol.vtk", vOlume);
3131  CHKERR SaveData(moab)(prefix + "out_surface.vtk", sUrface);
3132  CHKERR SaveData(moab)(prefix + "out_cut_edges.vtk", cutEdges);
3133  CHKERR SaveData(moab)(prefix + "out_cut_new_volumes.vtk", cutNewVolumes);
3134  CHKERR SaveData(moab)(prefix + "out_cut_new_surfaces.vtk", cutNewSurfaces);
3135  CHKERR SaveData(moab)(prefix + "out_cut_zero_distance_ents.vtk",
3138 }
3139 
3141  moab::Interface &moab = cOre.getInterface<CoreInterface>()->get_moab();
3143  CHKERR SaveData(moab)("out_trim_surface.vtk", sUrface);
3144  CHKERR SaveData(moab)("out_trim_new_volumes.vtk", trimNewVolumes);
3145  CHKERR SaveData(moab)("out_trim_new_surfaces.vtk", trimNewSurfaces);
3146  CHKERR SaveData(moab)("out_trim_edges.vtk", trimEdges);
3148 }
3149 
3150 } // namespace MoFEM
MoFEMErrorCode buildTree()
build tree
MoFEMErrorCode trimSurface(Range *fixed_edge, Range *corner_nodes, const bool debug=false)
Trim surface from faces beyond front.
map< EntityHandle, TreeData > edgesToTrim
virtual MoFEMErrorCode get_ref_ents(const RefEntity_multiIndex **refined_ents_ptr) const =0
Get ref entities multi-index from database.
MoFEM interface unique ID.
MoFEMErrorCode setTagData(Tag th, const BitRefLevel bit=BitRefLevel())
set coords to tag
MoFEMErrorCode findIfTringleHasThreeNodesOnInternalSurfaceSkin(const EntityHandle sideset, const BitRefLevel mesh_bit_level, const bool recursive, Range &faces_with_three_nodes_on_front, int verb=QUIET)
Find if triangle has three nodes on internal surface skin.
static double volumeLengthQuality(const double *coords)
Calculate tetrahedron volume length quality.
Definition: Tools.cpp:34
MoFEMErrorCode trimEdgesInTheMiddle(const BitRefLevel bit, const bool debug=false)
trim edges
virtual moab::Interface & get_moab()=0
double maxLength
Maximal edge length.
virtual MoFEMErrorCode refine_TET(const EntityHandle meshset, const BitRefLevel &bit, const bool respect_interface=false, int verb=QUIET, Range *ref_edges=NULL, const bool debug=false)
refine TET in the meshset
virtual MoFEMErrorCode add_vertices_in_the_middle_of_edges(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 copySurface(const Range &surface, Tag th=NULL, double *shift=NULL, double *origin=NULL, double *transform=NULL, const std::string save_mesh="")
copy surface entities
MoFEMErrorCode snapSurfaceToEdges(const Range &surface_edges, const Range &fixed_edges, const double rel_tol, const double abs_tol, Tag th=nullptr, const bool debug=false)
MoFEMErrorCode addBitRefLevel(const Range &ents, const BitRefLevel bit, int verb=QUIET) const
add bit ref level to ref entity
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:501
MoFEMErrorCode setSurface(const Range &surface)
set surface entities
static SEGMENT_MIN_DISTANCE minDistanceFromSegments(const double *w_ptr, const double *v_ptr, const double *k_ptr, const double *l_ptr, double *const tvw_ptr=nullptr, double *const tlk_ptr=nullptr)
Find points on two segments in closest distance.
Definition: Tools.cpp:345
map< EntityHandle, TreeData > verticesOnTrimEdges
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:74
MoFEMErrorCode makeFront(const bool debug=false)
Create front from the surface.
base class for all interface classes
static constexpr double delta
boost::shared_ptr< OrientedBoxTreeTool > treeSurfPtr
MoFEMErrorCode snapSurfaceSkinToEdges(const Range &fixed_edges, const double rel_tol, const double abs_tol, Tag th=nullptr, const bool debug=false)
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:477
MoFEMErrorCode setFront(const Range &surface)
MoFEMErrorCode mergeBadEdges(const int fraction_level, const Range &tets, const Range &surface, const Range &fixed_edges, const Range &corner_nodes, Range &merged_nodes, Range &out_tets, Range &new_surf, Tag th, const bool update_meshsets=false, const BitRefLevel *bit_ptr=NULL, const bool debug=false)
Merge edges.
MoFEMErrorCode splitSides(const EntityHandle meshset, const BitRefLevel &bit, const int msId, const CubitBCType cubit_bc_type, const bool add_interface_entities, const bool recursive=false, int verb=QUIET)
split nodes and other entities of tetrahedral in children sets and add prism elements
MoFEMErrorCode createSurfaceLevelSets(int verb=QUIET, const bool debug=false)
Calculate distance from mesh nodes to cut surface.
MoFEMErrorCode cutTrimAndMerge(int &first_bit, const int fraction_level, Tag th, const double tol_cut, const double tol_cut_close, const double tol_trim_close, Range &fixed_edges, Range &corner_nodes, const bool update_meshsets=false, const bool debug=false)
Cut, trim and merge.
double aveLength
Average edge length.
bool getSuccessMerge()
Return true if successful merge.
Definition: NodeMerger.hpp:44
Auxiliary tools.
Definition: Tools.hpp:30
Interface for managing meshsets containing materials and boundary conditions.
Core (interface) class.
Definition: Core.hpp:50
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return()
Definition: definitions.h:508
MoFEMErrorCode removePathologicalFrontTris(const BitRefLevel split_bit, Range &ents)
Remove pathological elements on surface internal front.
MoFEMErrorCode cutEdgesInMiddle(const BitRefLevel bit, Range &cut_vols, Range &cut_surf, Range &cut_verts, const bool debug=false)
cut edges
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
multi_index_container< ParentChild, indexed_by< ordered_unique< member< ParentChild, EntityHandle, &ParentChild::pArent > >, ordered_non_unique< member< ParentChild, EntityHandle, &ParentChild::cHild > > > > ParentChildMap
Definition: NodeMerger.hpp:128
VectorBoundedArray< double, 6 > VectorDouble6
Definition: Types.hpp:88
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:84
CutMeshInterface(const MoFEM::Core &core)
const Range & projectZeroDistanceEnts() const
VectorBoundedArray< double, 12 > VectorDouble12
Definition: Types.hpp:90
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
MoFEMErrorCode trimOnly(const BitRefLevel trim_bit, Tag th, const double tol_cut_close, Range *fixed_edges=NULL, Range *corner_nodes=NULL, const bool update_meshsets=false, const bool debug=false)
Trim mesh only.
MoFEMErrorCode saveCutEdges(const std::string prefix="")
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
Tensor2_Expr< transform_Tensor2< A, B, T, Dim0, Dim1, i, j >, T, Dim0, Dim1, i, j > transform(const Tensor2_Expr< A, T, Dim0, Dim1, i, j > &a, B function)
Merge node by collapsing edge between them.
Definition: NodeMerger.hpp:30
const Range & getNewTrimSurfaces() const
double tol
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
MatrixShallowArrayAdaptor< double > MatrixAdaptor
Matrix adaptor.
Definition: Types.hpp:118
MoFEMErrorCode cutAndTrim(int &first_bit, Tag th, const double tol_cut, const double tol_cut_close, const double tol_trim_close, Range *fixed_edges=NULL, Range *corner_nodes=NULL, const bool update_meshsets=false, const bool debug=false)
Cut and trim.
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Exceptions.hpp:66
static const bool debug
InterfaceThis interface is used by user to:
Definition: Interface.hpp:42
multi_index_container< LengthMapData, indexed_by< ordered_non_unique< member< LengthMapData, double, &LengthMapData::lEngth > >, hashed_unique< member< LengthMapData, EntityHandle, &LengthMapData::eDge > >, ordered_non_unique< member< LengthMapData, double, &LengthMapData::qUality > > > > LengthMapData_multi_index
MoFEMErrorCode setCoords(Tag th, const BitRefLevel bit=BitRefLevel(), const BitRefLevel mask=BitRefLevel().set())
set coords from tag
Managing BitRefLevels.
MoFEMErrorCode mergeSurface(const Range &surface)
merge surface entities
MoFEMErrorCode moveMidNodesOnCutEdges(Tag th=NULL)
projecting of mid edge nodes on new mesh on surface
MoFEMErrorCode mergeVolumes(const Range &volume)
merge volume entities
MoFEMErrorCode splitSides(const BitRefLevel split_bit, const BitRefLevel bit, const Range &ents, Tag th=NULL)
split sides
auto getVectorAdaptor(T1 ptr, const size_t n)
Get Vector adaptor.
Definition: Templates.hpp:43
MoFEMErrorCode query_interface(const MOFEMuuid &uuid, UnknownInterface **iface) const
Make interface on given faces and create flat prism in that space.
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Types.hpp:101
ParentChildMap & getParentChildMap()
Get map of parent cand child.
Definition: NodeMerger.hpp:134
MoFEMErrorCode mergeNodes(EntityHandle father, EntityHandle mother, Range &out_tets, Range *tets_ptr=NULL, const bool only_if_improve_quality=false, const double move=0, const int line_search=0, Tag th=NULL, const int verb=0)
merge nodes which sharing edge
Definition: NodeMerger.cpp:89
Mesh refinement interface.
MoFEMErrorCode getSides(const int msId, const CubitBCType cubit_bc_type, const BitRefLevel mesh_bit_level, const bool recursive, int verb=QUIET)
create two children meshsets in the meshset containing tetrahedral on two sides of faces
MoFEMErrorCode updateRange(const Range &parent, Range &child)
Update range by prents.
#define CHKERR
Inline error check.
Definition: definitions.h:596
static SEGMENT_MIN_DISTANCE minDistancePointFromOnSegment(const double *w_ptr, const double *v_ptr, const double *p_ptr, double *const t_ptr=nullptr)
Find closet point on the segment from the point.
Definition: Tools.cpp:322
static const MOFEMuuid IDD_MOFEMCutMesh
MoFEMErrorCode findEdgesToTrim(Range *fixed_edges, Range *corner_nodes, Tag th=NULL, const double tol=1e-4, const bool debug=false)
Find edges to trimEdges.
MoFEMErrorCode refineMesh(const int init_bit_level, const int surf_levels, const int front_levels, Range *fixed_edges=nullptr, int verb=QUIET, const bool debug=false)
Refine and set level sets.
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Types.hpp:85
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Types.hpp:51
multi_index_container< boost::shared_ptr< RefEntity >, indexed_by< ordered_unique< tag< Ent_mi_tag >, member< RefEntity::BasicEntity, EntityHandle, &RefEntity::ent > >, ordered_non_unique< tag< Ent_Ent_mi_tag >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > >, ordered_non_unique< tag< EntType_mi_tag >, const_mem_fun< RefEntity::BasicEntity, EntityType, &RefEntity::getEntType > >, ordered_non_unique< tag< ParentEntType_mi_tag >, const_mem_fun< RefEntity, EntityType, &RefEntity::getParentEntType > >, ordered_non_unique< tag< Composite_EntType_and_ParentEntType_mi_tag >, composite_key< RefEntity, const_mem_fun< RefEntity::BasicEntity, 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::BasicEntity, EntityType, &RefEntity::getEntType >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > > > > > RefEntity_multiIndex
MoFEMErrorCode createLevelSets(Tag th, Range &vol_edges, const bool remove_adj_prims_edges, int verb=QUIET, const bool debug=false, const std::string edges_file_name=string())
Create a level sets, i.e. distances from surface and surface front.
#define _IT_CUBITMESHSETS_FOR_LOOP_(MESHSET_MANAGER, IT)
Iterator that loops over all the Cubit MeshSets in a moFEM field.
MoFEMErrorCode findEdgesToCut(Range vol, Range *fixed_edges, Range *corner_nodes, const double geometry_tol, int verb=QUIET, const bool debug=false)
find edges to cut
map< EntityHandle, TreeData > verticesOnCutEdges
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:407
virtual MPI_Comm & get_comm() const =0
MoFEMErrorCode saveTrimEdges()
VectorBoundedArray< double, 9 > VectorDouble9
Definition: Types.hpp:89
MoFEMErrorCode createFrontLevelSets(Range vol, Tag th=nullptr, int verb=QUIET, const bool debug=false)
Calculate distance from mesh nodes to surface front.
map< EntityHandle, TreeData > edgesToCut
MoFEMErrorCode setVolume(const Range &volume)
set volume entities
MoFEMErrorCode cutOnly(Range vol, const BitRefLevel cut_bit, Tag th, const double tol_cut, const double tol_cut_close, Range *fixed_edges=NULL, Range *corner_nodes=NULL, const bool update_meshsets=false, const bool debug=false)
Cut mesh onlu.
MoFEMErrorCode moveMidNodesOnTrimmedEdges(Tag th=NULL)
move trimmed edges mid nodes