v0.8.20
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 
20 // FIXME: Some definiton of classes, like SaveData or LengthMapData, can
21 // generate conflicts. Those classes has to be packaged, such that will not
22 // generate conflicts someware else in the code. In princliple this type
23 // implementation is wrong and should be avoided.
24 
25 struct SaveData {
26  moab::Interface &moab;
27  SaveData(moab::Interface &moab) : moab(moab) {}
28  MoFEMErrorCode operator()(const std::string name, const Range &ents) {
30  EntityHandle meshset;
31  CHKERR moab.create_meshset(MESHSET_SET, meshset);
32  CHKERR moab.add_entities(meshset, ents);
33  CHKERR moab.write_file(name.c_str(), "VTK", "", &meshset, 1);
34  CHKERR moab.delete_entities(&meshset, 1);
36  }
37 };
38 
41  UnknownInterface **iface) const {
43  *iface = NULL;
44  if (uuid == IDD_MOFEMCutMesh) {
45  *iface = const_cast<CutMeshInterface *>(this);
47  }
48  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "unknown interface");
50 }
51 
53  : cOre(const_cast<Core &>(core)) {
54  lineSearchSteps = 10;
55  nbMaxMergingCycles = 200;
57 }
58 
61  treeSurfPtr.reset();
63 }
64 
67  fRont = front;
69 }
70 
73  sUrface = surface;
75 }
76 
77 MoFEMErrorCode CutMeshInterface::copySurface(const Range &surface, Tag th,
78  double *shift, double *origin,
79  double *transform,
80  const std::string save_mesh) {
81  CoreInterface &m_field = cOre;
82  moab::Interface &moab = m_field.get_moab();
84  std::map<EntityHandle, EntityHandle> verts_map;
85  for (Range::const_iterator tit = surface.begin(); tit != surface.end();
86  tit++) {
87  int num_nodes;
88  const EntityHandle *conn;
89  CHKERR moab.get_connectivity(*tit, conn, num_nodes, true);
90  MatrixDouble coords(num_nodes, 3);
91  if (th) {
92  CHKERR moab.tag_get_data(th, conn, num_nodes, &coords(0, 0));
93  } else {
94  CHKERR moab.get_coords(conn, num_nodes, &coords(0, 0));
95  }
96  EntityHandle new_verts[num_nodes];
97  for (int nn = 0; nn != num_nodes; nn++) {
98  if (verts_map.find(conn[nn]) != verts_map.end()) {
99  new_verts[nn] = verts_map[conn[nn]];
100  } else {
101  if (transform) {
102  ublas::matrix_row<MatrixDouble> mr(coords, nn);
103  if (origin) {
104  VectorAdaptor vec_origin(
105  3, ublas::shallow_array_adaptor<double>(3, origin));
106  mr = mr - vec_origin;
107  }
108  MatrixAdaptor mat_transform = MatrixAdaptor(
109  3, 3, ublas::shallow_array_adaptor<double>(9, transform));
110  mr = prod(mat_transform, mr);
111  if (origin) {
112  VectorAdaptor vec_origin(
113  3, ublas::shallow_array_adaptor<double>(3, origin));
114  mr = mr + vec_origin;
115  }
116  }
117  if (shift) {
118  ublas::matrix_row<MatrixDouble> mr(coords, nn);
119  VectorAdaptor vec_shift(
120  3, ublas::shallow_array_adaptor<double>(3, shift));
121  mr = mr + vec_shift;
122  }
123  CHKERR moab.create_vertex(&coords(nn, 0), new_verts[nn]);
124  verts_map[conn[nn]] = new_verts[nn];
125  }
126  }
127  EntityHandle ele;
128  CHKERR moab.create_element(MBTRI, new_verts, num_nodes, ele);
129  sUrface.insert(ele);
130  }
131  if (!save_mesh.empty())
132  CHKERR SaveData(m_field.get_moab())(save_mesh, sUrface);
134 }
135 
138  vOlume = volume;
140 }
141 
144  sUrface.merge(surface);
146 }
147 
150  vOlume.merge(volume);
152 }
153 
155  const Range &fixed_edges, const double rel_tol, const double abs_tol,
156  Tag th, const bool debug) {
157  CoreInterface &m_field = cOre;
158  moab::Interface &moab = m_field.get_moab();
160 
161  // Get cutting surface skin
162  Skinner skin(&moab);
163  Range surface_skin;
164  CHKERR skin.find_skin(0, sUrface, false, surface_skin);
165 
166  CHKERR snapSurfaceToEdges(surface_skin, fixed_edges, rel_tol, abs_tol, th,
167  debug);
168 
170 }
171 
173  const Range &fixed_edges,
174  const double rel_tol,
175  const double abs_tol,
176  Tag th, const bool debug) {
177  CoreInterface &m_field = cOre;
178  moab::Interface &moab = m_field.get_moab();
181 
182  auto get_point = [i](auto &t_w, auto &t_delta, auto t) {
184  t = std::max(0., std::min(1., t));
185  t_p(i) = t_w(i) + t * t_delta(i);
186  return t_p;
187  };
188 
189  auto get_distance = [i](auto &t_p, auto &t_n) {
190  FTensor::Tensor1<double, 3> t_dist_vector;
191  t_dist_vector(i) = t_p(i) - t_n(i);
192  return sqrt(t_dist_vector(i) * t_dist_vector(i));
193  };
194 
195  map<EntityHandle, double> map_verts_length;
196 
197  for (auto f : surface_edges) {
198  int num_nodes;
199  const EntityHandle *conn_skin;
200  CHKERR moab.get_connectivity(f, conn_skin, num_nodes, true);
201  VectorDouble6 coords_skin(6);
202  if (th)
203  CHKERR moab.tag_get_data(th, conn_skin, num_nodes, &coords_skin[0]);
204  else
205  CHKERR moab.get_coords(conn_skin, num_nodes, &coords_skin[0]);
207  &coords_skin[0], &coords_skin[1], &coords_skin[2]);
209  &coords_skin[3], &coords_skin[4], &coords_skin[5]);
210  FTensor::Tensor1<double, 3> t_skin_delta;
211  t_skin_delta(i) = t_n1(i) - t_n0(i);
212  const double skin_edge_length = sqrt(t_skin_delta(i) * t_skin_delta(i));
213  for (int nn = 0; nn != 2; ++nn) {
214  auto it = map_verts_length.find(conn_skin[nn]);
215  if (it != map_verts_length.end())
216  it->second = std::min(it->second, skin_edge_length);
217  else
218  map_verts_length[conn_skin[nn]] = skin_edge_length;
219  }
220  }
221 
222  for (auto m : map_verts_length) {
223 
225  if (th)
226  CHKERR moab.tag_get_data(th, &m.first, 1, &t_n(0));
227  else
228  CHKERR moab.get_coords(&m.first, 1, &t_n(0));
229 
230  double min_dist = rel_tol * m.second;
231  FTensor::Tensor1<double, 3> t_min_coords;
232  CHKERR cOre.getInterface<Tools>()->findMinDistanceFromTheEdges(
233  &t_n(0), fixed_edges, &min_dist, &t_min_coords(0));
234 
235  if (min_dist < rel_tol * m.second || min_dist < abs_tol) {
236  if (th)
237  CHKERR moab.tag_set_data(th, &m.first, 1, &t_min_coords(0));
238  else
239  CHKERR moab.set_coords(&m.first, 1, &t_min_coords(0));
240  }
241  }
242 
244 }
245 
247  CoreInterface &m_field = cOre;
248  moab::Interface &moab = m_field.get_moab();
250  treeSurfPtr = boost::shared_ptr<OrientedBoxTreeTool>(
251  new OrientedBoxTreeTool(&moab, "ROOTSETSURF", true));
254 }
255 
257  MoFEMErrorCode operator()(Core &core, const BitRefLevel &bit) const {
260  ->updateAllMeshsetsByEntitiesChildren(bit);
262  }
263 };
264 
266  int &first_bit, const int before_trim_levels, Tag th, const double tol_cut,
267  const double tol_cut_close, const double tol_trim,
268  const double tol_trim_close, Range *fixed_edges, Range *corner_nodes,
269  const bool update_meshsets, const bool debug) {
270  CoreInterface &m_field = cOre;
271  moab::Interface &moab = m_field.get_moab();
273 
274  std::vector<BitRefLevel> bit_levels;
275 
276  auto get_back_bit_levels = [&]() {
277  bit_levels.push_back(BitRefLevel());
278  bit_levels.back().set(first_bit + bit_levels.size() - 1);
279  return bit_levels.back();
280  };
281 
282  BitRefLevel bit_cut = get_back_bit_levels();
283 
284  // cut mesh
285  CHKERR findEdgesToCut(fixed_edges, corner_nodes, tol_cut, QUIET, debug);
286  CHKERR projectZeroDistanceEnts(fixed_edges, corner_nodes, tol_cut_close,
287  QUIET, debug);
288  CHKERR cutEdgesInMiddle(bit_cut, debug);
289  if (fixed_edges) {
290  CHKERR cOre.getInterface<BitRefManager>()->updateRange(*fixed_edges,
291  *fixed_edges);
292  }
293  if (corner_nodes) {
294  CHKERR cOre.getInterface<BitRefManager>()->updateRange(*corner_nodes,
295  *corner_nodes);
296  }
297  if (update_meshsets) {
298  CHKERR UpdateMeshsets()(cOre, bit_cut);
299  }
301 
302  auto get_min_quality = [&m_field](const BitRefLevel bit, Tag th) {
303  Range tets_level; // test at level
304  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
305  bit, BitRefLevel().set(), MBTET, tets_level);
306  double min_q = 1;
307  CHKERR m_field.getInterface<Tools>()->minTetsQuality(tets_level, min_q, th);
308  return min_q;
309  };
310 
311  PetscPrintf(PETSC_COMM_WORLD, "Min quality cut %6.4g\n",
312  get_min_quality(bit_cut, th));
313 
314  if (debug)
316 
317  Range starting_volume = cutNewVolumes;
318  Range starting_volume_nodes;
319  CHKERR m_field.get_moab().get_connectivity(starting_volume,
320  starting_volume_nodes, true);
321  std::vector<double> staring_volume_coords(3 * starting_volume_nodes.size());
322  CHKERR m_field.get_moab().get_coords(starting_volume_nodes,
323  &*staring_volume_coords.begin());
324 
325  if (th) {
326  std::vector<double> staring_volume_th_coords(3 *
327  starting_volume_nodes.size());
328  CHKERR m_field.get_moab().tag_get_data(th, starting_volume_nodes,
329  &*staring_volume_th_coords.begin());
330  CHKERR m_field.get_moab().set_coords(starting_volume_nodes,
331  &*staring_volume_th_coords.begin());
332  }
333 
334  if (th)
335  CHKERR setTagData(th);
336 
337  for (int ll = 0; ll != before_trim_levels; ++ll) {
338  CHKERR findEdgesToTrim(fixed_edges, corner_nodes, th, tol_trim,
339  tol_trim_close, debug);
340  BitRefLevel bit = get_back_bit_levels();
341  CHKERR refineBeforeTrim(bit, fixed_edges, true);
342  if (th)
343  CHKERR setTagData(th);
344  }
345 
346  BitRefLevel bit_trim = get_back_bit_levels();
347  // trim mesh
348  CHKERR findEdgesToTrim(fixed_edges, corner_nodes, th, tol_trim,
349  tol_trim_close, debug);
350  CHKERR trimEdgesInTheMiddle(bit_trim, th, tol_trim_close, debug);
351  if (fixed_edges) {
352  CHKERR cOre.getInterface<BitRefManager>()->updateRange(*fixed_edges,
353  *fixed_edges);
354  }
355  if (corner_nodes) {
356  CHKERR cOre.getInterface<BitRefManager>()->updateRange(*corner_nodes,
357  *corner_nodes);
358  }
359  if (update_meshsets) {
360  CHKERR UpdateMeshsets()(cOre, bit_trim);
361  }
363 
364  PetscPrintf(PETSC_COMM_WORLD, "Min quality trim %3.2g\n",
365  get_min_quality(bit_trim, th));
366 
367  if (debug) {
369  Range bit2_edges;
370  CHKERR cOre.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
371  bit_trim, BitRefLevel().set(), MBEDGE, bit2_edges);
372  CHKERR SaveData(moab)("trim_fixed_edges.vtk",
373  intersect(*fixed_edges, bit2_edges));
374  }
375 
376  first_bit += bit_levels.size() - 1;
377 
378  if (th)
379  CHKERR m_field.get_moab().set_coords(starting_volume_nodes,
380  &*staring_volume_coords.begin());
381 
383 }
384 
386  int &first_bit, const int fraction_level, const int before_trim_levels,
387  Tag th, const double tol_cut, const double tol_cut_close,
388  const double tol_trim, const double tol_trim_close, Range &fixed_edges,
389  Range &corner_nodes, const bool update_meshsets, const bool debug) {
390  CoreInterface &m_field = cOre;
392 
393  std::vector<BitRefLevel> bit_levels;
394 
395  auto get_back_bit_levels = [&]() {
396  bit_levels.push_back(BitRefLevel());
397  bit_levels.back().set(first_bit + bit_levels.size() - 1);
398  return bit_levels.back();
399  };
400 
401  if (debug) {
402  CHKERR cOre.getInterface<BitRefManager>()->writeEntitiesNotInDatabase(
403  "ents_not_in_database.vtk", "VTK", "");
404  }
405 
406  CHKERR cutAndTrim(first_bit, before_trim_levels, th, tol_cut, tol_cut_close,
407  tol_trim, tol_trim_close, &fixed_edges, &corner_nodes,
408  update_meshsets, debug);
409  if (debug) {
410  CHKERR cOre.getInterface<BitRefManager>()->writeEntitiesNotInDatabase(
411  "cut_trim_ents_not_in_database.vtk", "VTK", "");
412  }
413 
414  BitRefLevel bit_level1 = BitRefLevel().set(first_bit - 1);
415  for (int ll = 1; ll <= before_trim_levels; ++ll)
416  bit_level1.set(first_bit - 1 - ll);
417  BitRefLevel bit_level2 = get_back_bit_levels();
418  BitRefLevel bit_level3 = get_back_bit_levels();
419 
420  CHKERR mergeBadEdges(fraction_level, bit_level1, bit_level2, bit_level3,
421  getNewTrimSurfaces(), fixed_edges, corner_nodes, th,
422  update_meshsets, debug);
423 
424  if (debug) {
425  CHKERR cOre.getInterface<BitRefManager>()->writeBitLevelByType(
426  bit_level3, BitRefLevel().set(), MBTET, "out_tets_merged.vtk", "VTK",
427  "");
428  CHKERR cOre.getInterface<BitRefManager>()->writeEntitiesNotInDatabase(
429  "cut_trim_merge_ents_not_in_database.vtk", "VTK", "");
430  CHKERR SaveData(m_field.get_moab())("merged_surface.vtk", mergedSurfaces);
431  }
432 
433  auto get_min_quality = [&m_field, debug](const BitRefLevel bit, Tag th) {
434  Range tets_level; // test at level
435  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
436  bit, BitRefLevel().set(), MBTET, tets_level);
437  double min_q = 1;
438  CHKERR m_field.getInterface<Tools>()->minTetsQuality(tets_level, min_q, th);
439  if (min_q < 0 && debug) {
440  CHKERR m_field.getInterface<Tools>()->writeTetsWithQuality(
441  "negative_tets.vtk", "VTK", "", tets_level, th);
442  }
443  return min_q;
444  };
445 
446  PetscPrintf(PETSC_COMM_WORLD, "Min quality node merge %6.4g\n",
447  get_min_quality(bit_level3, th));
448 
449  CHKERR cOre.getInterface<BitRefManager>()->updateRange(fixed_edges,
450  fixed_edges);
451  CHKERR cOre.getInterface<BitRefManager>()->updateRange(corner_nodes,
452  corner_nodes);
453 
454  first_bit += bit_levels.size() - 1;
455 
457 }
458 
460  int &first_bit, const int fraction_level, const int ref_before_cut_levels,
461  const int befor_trim_levels, Tag th, const double tol_cut,
462  const double tol_cut_close, const double tol_trim,
463  const double tol_trim_close, Range &fixed_edges, Range &corner_nodes,
464  const bool update_meshsets, const bool debug) {
465  CoreInterface &m_field = cOre;
467 
468  std::vector<BitRefLevel> bit_levels;
469 
470  auto get_back_bit_levels = [&]() {
471  bit_levels.push_back(BitRefLevel());
472  bit_levels.back().set(first_bit + bit_levels.size() - 1);
473  return bit_levels.back();
474  };
475 
476  Range starting_volume = vOlume;
477  Range starting_volume_nodes;
478  CHKERR m_field.get_moab().get_connectivity(starting_volume,
479  starting_volume_nodes, true);
480  std::vector<double> staring_volume_coords(3 * starting_volume_nodes.size());
481  CHKERR m_field.get_moab().get_coords(starting_volume_nodes,
482  &*staring_volume_coords.begin());
483 
484  if (th) {
485  std::vector<double> staring_volume_th_coords(3 *
486  starting_volume_nodes.size());
487  CHKERR m_field.get_moab().tag_get_data(th, starting_volume_nodes,
488  &*staring_volume_th_coords.begin());
489  CHKERR m_field.get_moab().set_coords(starting_volume_nodes,
490  &*staring_volume_th_coords.begin());
491  }
492 
493  for (int ll = 0; ll != ref_before_cut_levels; ++ll) {
494  CHKERR findEdgesToCut(&fixed_edges, &corner_nodes, tol_trim, QUIET);
495  BitRefLevel bit = get_back_bit_levels();
496  CHKERR refineBeforeCut(bit, &fixed_edges, true);
497  if (th)
498  CHKERR setTagData(th);
499  }
500 
501  if (debug && ref_before_cut_levels)
502  CHKERR m_field.getInterface<BitRefManager>()->writeBitLevelByType(
503  bit_levels.back(), BitRefLevel().set(), MBTET,
504  "out_after_refine_before_cut.vtk", "VTK", "");
505 
506  first_bit += bit_levels.size();
507  CHKERR cutTrimAndMerge(first_bit, fraction_level, befor_trim_levels, th,
508  tol_cut, tol_cut_close, tol_trim, tol_trim_close,
509  fixed_edges, corner_nodes, update_meshsets, debug);
510 
511  if (th)
512  CHKERR m_field.get_moab().set_coords(starting_volume_nodes,
513  &*staring_volume_coords.begin());
514 
516 }
517 
518 
520  Range *corner_nodes,
521  const double low_tol, int verb,
522  const bool debug) {
523  CoreInterface &m_field = cOre;
524  moab::Interface &moab = m_field.get_moab();
526 
527  edgesToCut.clear();
528  cutEdges.clear();
529 
530  zeroDistanceVerts.clear();
531  zeroDistanceEnts.clear();
532  verticesOnCutEdges.clear();
533 
534  double ray_length;
535  double ray_point[3], unit_ray_dir[3];
536  VectorAdaptor vec_unit_ray_dir(
537  3, ublas::shallow_array_adaptor<double>(3, unit_ray_dir));
538  VectorAdaptor vec_ray_point(
539  3, ublas::shallow_array_adaptor<double>(3, ray_point));
540 
541  Tag th_dist;
542  rval = moab.tag_get_handle("DIST", th_dist);
543  if (rval == MB_SUCCESS)
544  CHKERR moab.tag_delete(th_dist);
545  else
546  rval = MB_SUCCESS;
547 
548  Tag th_dist_normal;
549  rval = moab.tag_get_handle("DIST_NORMAL", th_dist_normal);
550  if (rval == MB_SUCCESS)
551  CHKERR moab.tag_delete(th_dist_normal);
552  else
553  rval = MB_SUCCESS;
554 
555  double def_val[] = {0, 0, 0};
556  CHKERR moab.tag_get_handle("DIST", 1, MB_TYPE_DOUBLE, th_dist,
557  MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
558  CHKERR moab.tag_get_handle("DIST_NORMAL", 3, MB_TYPE_DOUBLE, th_dist_normal,
559  MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
560 
561  // Calculate not signed and not signed distances from nodes to surface
562  Range vol_vertices;
563  CHKERR moab.get_connectivity(vOlume, vol_vertices, true);
564  for (auto v : vol_vertices) {
565  VectorDouble3 point_in(3);
566  CHKERR moab.get_coords(&v, 1, &point_in[0]);
567  VectorDouble3 point_out(3);
568  EntityHandle facets_out;
569  CHKERR treeSurfPtr->closest_to_location(&point_in[0], rootSetSurf,
570  &point_out[0], facets_out);
571  VectorDouble3 n(3);
572  Util::normal(&moab, facets_out, n[0], n[1], n[2]);
573  VectorDouble3 delta = point_out - point_in;
574  double dist = norm_2(delta);
575  n /= norm_2(n);
576  VectorDouble3 dist_normal = inner_prod(delta, n) * n;
577  // not signed distance
578  CHKERR moab.tag_set_data(th_dist, &v, 1, &dist);
579  // signed distance
580  CHKERR moab.tag_set_data(th_dist_normal, &v, 1, &dist_normal[0]);
581  }
582 
583  auto get_normal_dist = [](const double *normal) {
584  FTensor::Tensor1<const double *, 3> t_n(normal, &normal[1], &normal[2]);
586  return sqrt(t_n(i) * t_n(i));
587  };
588 
589  auto get_edge_crossed = [&moab, get_normal_dist,
590  th_dist_normal](EntityHandle v0, EntityHandle v1) {
591  VectorDouble3 dist_normal0(3);
592  CHKERR moab.tag_get_data(th_dist_normal, &v0, 1, &dist_normal0[0]);
593  VectorDouble3 dist_normal1(3);
594  CHKERR moab.tag_get_data(th_dist_normal, &v1, 1, &dist_normal1[0]);
595  return (inner_prod(dist_normal0, dist_normal1) < 0);
596  };
597 
598  auto get_normal_dist_from_conn = [&moab, get_normal_dist,
599  th_dist_normal](EntityHandle v) {
600  double dist_normal[3];
601  CHKERR moab.tag_get_data(th_dist_normal, &v, 1, dist_normal);
602  return get_normal_dist(dist_normal);
603  };
604 
605  auto get_ray_for_edge = [&](const EntityHandle ent, VectorAdaptor &ray_point,
606  VectorAdaptor &unit_ray_dir, double &ray_length) {
608  int num_nodes;
609  const EntityHandle *conn;
610  CHKERR moab.get_connectivity(ent, conn, num_nodes, true);
611  double coords[6];
612  CHKERR moab.get_coords(conn, num_nodes, coords);
613  VectorAdaptor s0(3, ublas::shallow_array_adaptor<double>(3, &coords[0]));
614  VectorAdaptor s1(3, ublas::shallow_array_adaptor<double>(3, &coords[3]));
615  noalias(ray_point) = s0;
616  noalias(unit_ray_dir) = s1 - s0;
617  ray_length = norm_2(unit_ray_dir);
618  if (ray_length > 0)
619  unit_ray_dir /= ray_length;
621  };
622 
623  auto set_edge_cut_point_from_level_set = [&](const EntityHandle e,
624  const auto dist_normal) {
626  CHKERR get_ray_for_edge(e, vec_ray_point, vec_unit_ray_dir, ray_length);
627  const double s = dist_normal[0] / (dist_normal[1] + dist_normal[0]);
628  edgesToCut[e].dIst = s * ray_length;
629  edgesToCut[e].lEngth = ray_length;
630  edgesToCut[e].unitRayDir = vec_unit_ray_dir;
631  edgesToCut[e].rayPoint = vec_ray_point;
632  cutEdges.insert(e);
634  };
635 
636  auto not_project_node = [this, &moab](const EntityHandle v) {
638  VectorDouble3 s0(3);
639  CHKERR moab.get_coords(&v, 1, &s0[0]);
640  verticesOnCutEdges[v].dIst = 0;
641  verticesOnCutEdges[v].lEngth = 0;
642  verticesOnCutEdges[v].unitRayDir.resize(3, false);
643  verticesOnCutEdges[v].unitRayDir.clear();
644  verticesOnCutEdges[v].rayPoint = s0;
646  };
647 
648  auto get_ave_edge_length = [&](const EntityHandle ent,
649  const Range &vol_edges) {
650  Range adj_edges;
651  if (moab.type_from_handle(ent) == MBVERTEX) {
652  CHKERR moab.get_adjacencies(&ent, 1, 1, false, adj_edges);
653  } else {
654  adj_edges.insert(ent);
655  }
656  adj_edges = intersect(adj_edges, vol_edges);
657  double ave_l = 0;
658  for (auto e : adj_edges) {
659  int num_nodes;
660  const EntityHandle *conn;
661  CHKERR moab.get_connectivity(e, conn, num_nodes, true);
662  VectorDouble6 coords(6);
663  CHKERR moab.get_coords(conn, num_nodes, &coords[0]);
665  &coords[0], &coords[1], &coords[2]);
667  &coords[3], &coords[4], &coords[5]);
669  t_n0(i) -= t_n1(i);
670  ave_l += sqrt(t_n0(i) * t_n0(i));
671  }
672  return ave_l / adj_edges.size();
673  };
674 
675  Range vol_edges;
676  CHKERR moab.get_adjacencies(vOlume, 1, true, vol_edges,
677  moab::Interface::UNION);
678 
679  aveLength = 0;
680  maxLength = 0;
681  int nb_ave_length = 0;
682  for (auto e : vol_edges) {
683  int num_nodes;
684  const EntityHandle *conn;
685  CHKERR moab.get_connectivity(e, conn, num_nodes, true);
686  double dist[num_nodes];
687  CHKERR moab.tag_get_data(th_dist, conn, num_nodes, dist);
688  CHKERR get_ray_for_edge(e, vec_ray_point, vec_unit_ray_dir, ray_length);
689  const double tol = ray_length * low_tol;
690 
691  // Edges is on two sides of the surface
692  if (get_edge_crossed(conn[0], conn[1])) {
693  std::vector<double> distances_out;
694  std::vector<EntityHandle> facets_out;
695  // sen the ray in the direction of the edge
696  CHKERR treeSurfPtr->ray_intersect_triangles(distances_out, facets_out,
697  rootSetSurf, tol, ray_point,
698  unit_ray_dir, &ray_length);
699  // if ray punching through the surface
700  if (!distances_out.empty()) {
701  // minimal distance
702  const auto dist_ptr =
703  std::min_element(distances_out.begin(), distances_out.end());
704  const double dist = *dist_ptr;
705  // if distance is shorter than edge length
706  if (dist <= ray_length) {
707  aveLength += ray_length;
708  maxLength = std::max(maxLength, ray_length);
709  nb_ave_length++;
710  edgesToCut[e].dIst = dist;
711  edgesToCut[e].lEngth = ray_length;
712  edgesToCut[e].unitRayDir = vec_unit_ray_dir;
713  edgesToCut[e].rayPoint = vec_ray_point;
714  cutEdges.insert(e);
715  }
716  }
717  }
718 
719  // check if cutting point is close to the end of the edges
720  if (fabs(dist[0]) < tol && fabs(dist[1]) < tol) {
721  aveLength += ray_length;
722  maxLength = fmax(maxLength, ray_length);
723  nb_ave_length++;
724  for (int n : {0, 1})
725  CHKERR not_project_node(conn[n]);
726  zeroDistanceEnts.insert(e);
727  }
728  }
729  aveLength /= nb_ave_length;
730 
731  for (auto v : vol_vertices) {
732  double dist;
733  CHKERR moab.tag_get_data(th_dist, &v, 1, &dist);
734  const double tol = get_ave_edge_length(v, vol_edges) * low_tol;
735  if (fabs(dist) < tol) {
736  CHKERR not_project_node(v);
737  zeroDistanceVerts.insert(v);
738  }
739  }
740 
741  cutVolumes.clear();
742  Range cut_edges_verts;
743  CHKERR moab.get_connectivity(unite(cutEdges, zeroDistanceEnts),
744  cut_edges_verts, true);
745  CHKERR moab.get_adjacencies(unite(cut_edges_verts, zeroDistanceVerts), 3,
746  false, cutVolumes, moab::Interface::UNION);
747  cutVolumes = intersect(cutVolumes, vOlume);
748 
749  // get edges on the cut volumes
750  Range edges;
751  CHKERR moab.get_adjacencies(cutVolumes, 1, false, edges,
752  moab::Interface::UNION);
753  edges = subtract(edges, cutEdges);
754 
755  // add to cut set edges which are cutted by extension of cutting surface
756  for (auto e : edges) {
757  int num_nodes;
758  const EntityHandle *conn;
759  CHKERR moab.get_connectivity(e, conn, num_nodes, true);
760  const double tol = get_ave_edge_length(e, vol_edges) * low_tol;
761  double dist_normal[2];
762  dist_normal[0] = get_normal_dist_from_conn(conn[0]);
763  dist_normal[1] = get_normal_dist_from_conn(conn[1]);
764  if (get_edge_crossed(conn[0], conn[1]))
765  CHKERR set_edge_cut_point_from_level_set(e, dist_normal);
766  else if (fabs(dist_normal[0]) < tol && fabs(dist_normal[1]) < tol) {
767  for (int n : {0, 1})
768  CHKERR not_project_node(conn[n]);
769  zeroDistanceEnts.insert(e);
770  }
771  }
772 
773  for (auto f : zeroDistanceEnts) {
774  int num_nodes;
775  const EntityHandle *conn;
776  CHKERR moab.get_connectivity(f, conn, num_nodes, true);
777  Range adj_edges;
778  CHKERR moab.get_adjacencies(conn, num_nodes, 1, false, adj_edges,
779  moab::Interface::UNION);
780  for (auto e : adj_edges) {
781  cutEdges.erase(e);
782  edgesToCut.erase(e);
783  }
784  }
785 
786  for (auto v : zeroDistanceVerts) {
787  Range adj_edges;
788  CHKERR moab.get_adjacencies(&v, 1, 1, false, adj_edges,
789  moab::Interface::UNION);
790  for (auto e : adj_edges) {
791  cutEdges.erase(e);
792  edgesToCut.erase(e);
793  }
794  }
795 
796  if (debug)
797  CHKERR SaveData(m_field.get_moab())("cut_edges.vtk", cutEdges);
798 
799  if (debug)
800  CHKERR SaveData(m_field.get_moab())(
801  "cut_zero_dist.vtk", unite(zeroDistanceVerts, zeroDistanceEnts));
802 
804 }
805 
807  Range *fixed_edges,
808  const bool update_meshsets,
809  const bool debug) {
810  CoreInterface &m_field = cOre;
811  moab::Interface &moab = m_field.get_moab();
812  MeshRefinement *refiner;
813  BitRefManager *bit_ref_manager;
815 
816  if (cutEdges.size() != edgesToCut.size()) {
817  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
818  }
819  CHKERR m_field.getInterface(refiner);
820  CHKERR m_field.getInterface(bit_ref_manager);
821 
822  Range ref_edges_tets;
823  CHKERR moab.get_adjacencies(cutEdges, 3, false, ref_edges_tets,
824  moab::Interface::UNION);
825  CHKERR moab.get_adjacencies(zeroDistanceVerts, 3, false, ref_edges_tets,
826  moab::Interface::UNION);
827  Range ref_edges;
828  CHKERR moab.get_adjacencies(ref_edges_tets, 1, true, ref_edges,
829  moab::Interface::UNION);
830  CHKERR refiner->add_verices_in_the_middel_of_edges(ref_edges, bit);
831  CHKERR refiner->refine_TET(vOlume, bit, false, QUIET,
832  debug ? &cutEdges : NULL);
833 
834  auto update_range = [bit_ref_manager](Range *r_ptr) {
836  if (r_ptr) {
837  Range childs;
838  CHKERR bit_ref_manager->updateRange(*r_ptr, childs);
839  r_ptr->merge(childs);
840  }
842  };
843  CHKERR update_range(fixed_edges);
844  CHKERR update_range(&vOlume);
845 
846  if (update_meshsets)
847  CHKERR UpdateMeshsets()(cOre, bit);
848 
849  Range bit_tets;
850  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
851  bit, BitRefLevel().set(), MBTET, bit_tets);
852  vOlume = intersect(vOlume, bit_tets);
853 
854  Range bit_edges;
855  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
856  bit, BitRefLevel().set(), MBEDGE, bit_edges);
857  if (fixed_edges)
858  *fixed_edges = intersect(*fixed_edges, bit_edges);
859 
861 }
862 
864  Range *fixed_edges,
865  const bool update_meshsets,
866  const bool debug) {
867  CoreInterface &m_field = cOre;
868  moab::Interface &moab = m_field.get_moab();
869  MeshRefinement *refiner;
870  BitRefManager *bit_ref_manager;
872 
873  // refine before trim
874  Range vol_trim_edges;
875  CHKERR moab.get_adjacencies(trimEdges, 3, true, vol_trim_edges,
876  moab::Interface::UNION);
877  CHKERR moab.get_adjacencies(trimNewVertices, 3, true, vol_trim_edges,
878  moab::Interface::UNION);
879  vol_trim_edges = intersect(vol_trim_edges, cutNewVolumes);
880 
881  CHKERR m_field.getInterface(refiner);
882  CHKERR m_field.getInterface(bit_ref_manager);
883 
884  Range ref_edges;
885  CHKERR moab.get_adjacencies(vol_trim_edges, 1, true, ref_edges,
886  moab::Interface::UNION);
887  CHKERR refiner->add_verices_in_the_middel_of_edges(ref_edges, bit);
888  CHKERR refiner->refine_TET(cutNewVolumes, bit, false, QUIET,
889  debug ? &cutEdges : NULL);
890 
891  auto update_range = [bit_ref_manager](Range *r_ptr) {
893  if (r_ptr) {
894  Range childs;
895  CHKERR bit_ref_manager->updateRange(*r_ptr, childs);
896  r_ptr->merge(childs);
897  }
899  };
900  CHKERR update_range(fixed_edges);
901  CHKERR update_range(&cutNewVolumes);
902  CHKERR update_range(&cutNewSurfaces);
903 
904  if (update_meshsets)
905  CHKERR UpdateMeshsets()(cOre, bit);
906 
907  Range bit_tets;
908  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
909  bit, BitRefLevel().set(), MBTET, bit_tets);
910  cutNewVolumes = intersect(cutNewVolumes, bit_tets);
911 
912  Range bit_faces;
913  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
914  bit, BitRefLevel().set(), MBTRI, bit_faces);
915  cutNewSurfaces = intersect(cutNewSurfaces, bit_faces);
916  cutNewVertices.clear();
917  CHKERR moab.get_connectivity(cutNewSurfaces, cutNewVertices, true);
918 
919  Range bit_edges;
920  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
921  bit, BitRefLevel().set(), MBEDGE, bit_edges);
922  if (fixed_edges)
923  *fixed_edges = intersect(*fixed_edges, bit_edges);
924 
926 }
927 
929  Range *corner_nodes,
930  const double low_tol,
931  const int verb,
932  const bool debug) {
933  CoreInterface &m_field = cOre;
934  moab::Interface &moab = m_field.get_moab();
936 
937  // Get entities on body skin
938  Skinner skin(&moab);
939  Range tets_skin;
940  rval = skin.find_skin(0, vOlume, false, tets_skin);
941  Range tets_skin_edges;
942  CHKERR moab.get_adjacencies(tets_skin, 1, false, tets_skin_edges,
943  moab::Interface::UNION);
944  Range tets_skin_verts;
945  CHKERR moab.get_connectivity(tets_skin, tets_skin_verts, true);
946 
947  // Get entities in volume
948  Range vol_faces, vol_edges, vol_nodes;
949  CHKERR moab.get_adjacencies(vOlume, 2, false, vol_faces,
950  moab::Interface::UNION);
951  CHKERR moab.get_adjacencies(vOlume, 1, false, vol_edges,
952  moab::Interface::UNION);
953  CHKERR moab.get_connectivity(vOlume, vol_nodes, true);
954 
955  // Get nodes on cut edges
956  Range cut_edge_verts;
957  CHKERR moab.get_connectivity(cutEdges, cut_edge_verts, true);
958 
959  // Get faces and edges
960  Range cut_edges_faces;
961  CHKERR moab.get_adjacencies(cut_edge_verts, 2, true, cut_edges_faces,
962  moab::Interface::UNION);
963  cut_edges_faces = intersect(cut_edges_faces, vol_faces);
964  Range cut_edges_faces_verts;
965  CHKERR moab.get_connectivity(cut_edges_faces, cut_edges_faces_verts, true);
966  cut_edges_faces_verts = subtract(cut_edges_faces_verts, cut_edge_verts);
967  Range to_remove_cut_edges_faces;
968  CHKERR moab.get_adjacencies(cut_edges_faces_verts, 2, true,
969  to_remove_cut_edges_faces,
970  moab::Interface::UNION);
971  // Those are faces which have vertices adjacent to cut edges vertices without
972  // hanging nodes (nodes not adjacent to cutting edge)
973  cut_edges_faces = subtract(cut_edges_faces, to_remove_cut_edges_faces);
974  if (debug)
975  CHKERR SaveData(moab)("cut_edges_faces.vtk", cut_edges_faces);
976  cut_edges_faces.merge(cutEdges);
977 
978  Range fixed_edges_nodes;
979  if (fixed_edges)
980  CHKERR moab.get_connectivity(*fixed_edges, fixed_edges_nodes, true);
981 
982  Tag th_dist_normal;
983  CHKERR moab.tag_get_handle("DIST_NORMAL", th_dist_normal);
984 
985  // Create map of closes points to the surface
986  enum TYPE { FREE = 0, SKIN, FIXED, CORNER, NOT_MOVE };
987  map<EntityHandle, std::pair<std::pair<TYPE, EntityHandle>, TreeData>>
988  min_dist_map;
989  double ave_cut_edge_length = 0;
990  for (auto e : cutEdges) {
991 
992  TYPE edge_type = FREE;
993  if (tets_skin_edges.find(e) != tets_skin_edges.end())
994  edge_type = SKIN;
995  if (fixed_edges)
996  if (fixed_edges->find(e) != fixed_edges->end())
997  edge_type = FIXED;
998 
999  auto eit = edgesToCut.find(e);
1000  if (eit != edgesToCut.end()) {
1001 
1002  int num_nodes;
1003  const EntityHandle *conn;
1004  CHKERR moab.get_connectivity(e, conn, num_nodes, true);
1005  VectorDouble6 pos(6);
1006  CHKERR moab.get_coords(conn, num_nodes, &pos[0]);
1007  VectorDouble3 p[2];
1008  p[0] = getVectorAdaptor(&pos[0], 3);
1009  p[1] = getVectorAdaptor(&pos[3], 3);
1010  ave_cut_edge_length += norm_2(p[0] - p[1]);
1011 
1012  VectorDouble6 dist_normal(6);
1013  CHKERR moab.tag_get_data(th_dist_normal, conn, num_nodes,
1014  &dist_normal[0]);
1015  auto get_normal_dist = [](const double *normal) {
1016  FTensor::Tensor1<const double *, 3> t_n(normal, &normal[1], &normal[2]);
1018  return sqrt(t_n(i) * t_n(i));
1019  };
1020 
1021  auto &d = eit->second;
1022  VectorDouble3 new_pos = d.rayPoint + d.dIst * d.unitRayDir;
1023  for (int nn = 0; nn != 2; ++nn) {
1024 
1025  VectorDouble3 ray = new_pos - p[nn];
1026  const double dist = norm_2(ray);
1027  const double length = get_normal_dist(&dist_normal[3 * nn]);
1028 
1029  bool add_node = true;
1030  auto vit = min_dist_map.find(conn[nn]);
1031  if (vit != min_dist_map.end()) {
1032  if (vit->second.first.first > edge_type)
1033  add_node = false;
1034  else if (vit->second.first.first == edge_type) {
1035  if (vit->second.second.dIst < dist)
1036  add_node = false;
1037  }
1038  }
1039 
1040  if (add_node) {
1041  min_dist_map[conn[nn]].first.first = edge_type;
1042  min_dist_map[conn[nn]].first.second = e;
1043  auto &data = min_dist_map[conn[nn]].second;
1044  data.lEngth = length;
1045  data.rayPoint = p[nn];
1046  data.dIst = dist;
1047  if (dist > 0)
1048  data.unitRayDir = ray / dist;
1049  else {
1050  data.unitRayDir.resize(3);
1051  data.unitRayDir.clear();
1052  }
1053  }
1054  }
1055 
1056  } else
1057  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Edge not found");
1058  }
1059 
1060  ave_cut_edge_length /= cutEdges.size();
1061 
1062  auto get_quality_change =
1063  [this, &m_field,
1064  &moab](const Range &adj_tets,
1065  map<EntityHandle, TreeData> vertices_on_cut_edges) {
1066  vertices_on_cut_edges.insert(verticesOnCutEdges.begin(),
1067  verticesOnCutEdges.end());
1068  double q0 = 1;
1069  CHKERR m_field.getInterface<Tools>()->minTetsQuality(adj_tets, q0);
1070  double min_q = 1;
1071  for (auto t : adj_tets) {
1072  int num_nodes;
1073  const EntityHandle *conn;
1074  CHKERR m_field.get_moab().get_connectivity(t, conn, num_nodes, true);
1075  VectorDouble12 coords(12);
1076  CHKERR moab.get_coords(conn, num_nodes, &coords[0]);
1077  for (int n = 0; n != 4; ++n) {
1078  bool ray_found = false;
1079  auto mit = vertices_on_cut_edges.find(conn[n]);
1080  if (mit != vertices_on_cut_edges.end()) {
1081  ray_found = true;
1082  }
1083  if (ray_found) {
1084  auto n_coords = getVectorAdaptor(&coords[3 * n], 3);
1085  double dist = mit->second.dIst;
1086  noalias(n_coords) =
1087  mit->second.rayPoint + dist * mit->second.unitRayDir;
1088  }
1089  }
1090  const double q = Tools::volumeLengthQuality(&coords[0]);
1091  if (!std::isnormal(q)) {
1092  min_q = -2;
1093  break;
1094  }
1095  min_q = std::min(min_q, q);
1096  }
1097  return min_q / q0;
1098  };
1099 
1100  auto get_conn = [&moab](const EntityHandle &e, const EntityHandle *&conn,
1101  int &num_nodes) {
1103  EntityType type = moab.type_from_handle(e);
1104  if (type == MBVERTEX) {
1105  conn = &e;
1106  num_nodes = 1;
1107  } else {
1108  CHKERR moab.get_connectivity(e, conn, num_nodes, true);
1109  }
1111  };
1112 
1113  auto project_node = [&](const EntityHandle v,
1114  map<EntityHandle, TreeData> &vertices_on_cut_edges) {
1116  vertices_on_cut_edges[v].dIst = min_dist_map[v].second.dIst;
1117  vertices_on_cut_edges[v].lEngth = min_dist_map[v].second.lEngth;
1118  vertices_on_cut_edges[v].unitRayDir = min_dist_map[v].second.unitRayDir;
1119  vertices_on_cut_edges[v].rayPoint = min_dist_map[v].second.rayPoint;
1121  };
1122 
1123  auto remove_surface_tets = [&](Range &zero_dist_tets,
1124  Range zero_distance_ents,
1125  Range zero_distance_verts) {
1127  Range zero_dist_all_verts;
1128  CHKERR moab.get_connectivity(zero_distance_ents, zero_dist_all_verts, true);
1129  zero_dist_all_verts.merge(zero_distance_verts);
1130  CHKERR moab.get_adjacencies(zero_dist_all_verts, 3, false, zero_dist_tets,
1131  moab::Interface::UNION);
1132  Range zero_tets_verts;
1133  CHKERR moab.get_connectivity(zero_dist_tets, zero_tets_verts, true);
1134  zero_tets_verts = subtract(zero_tets_verts, zero_dist_all_verts);
1135  Range free_tets;
1136  CHKERR moab.get_adjacencies(zero_tets_verts, 3, false, free_tets,
1137  moab::Interface::UNION);
1138  zero_dist_tets = subtract(zero_dist_tets, free_tets);
1139 
1141  };
1142 
1143  for (int d = 2; d >= 0; --d) {
1144 
1145  Range ents;
1146  if (d > 0)
1147  ents = cut_edges_faces.subset_by_dimension(d);
1148  else
1149  ents = cut_edge_verts;
1150 
1151  // make list of entities
1152  multimap<double, EntityHandle> ents_to_check;
1153  for (auto f : ents) {
1154  int num_nodes;
1155  const EntityHandle *conn;
1156  CHKERR get_conn(f, conn, num_nodes);
1157  VectorDouble3 dist(3);
1158 
1159  bool move = true;
1160  for (int n = 0; n != num_nodes; ++n) {
1161  auto &d = min_dist_map[conn[n]];
1162  if (d.first.first == NOT_MOVE) {
1163  move = false;
1164  break;
1165  }
1166  dist[n] = d.second.lEngth;
1167  }
1168 
1169  if (move) {
1170  double max_dist = 0;
1171  for (int n = 0; n != num_nodes; ++n)
1172  max_dist = std::max(max_dist, fabs(dist[n]));
1173  if (max_dist < low_tol * ave_cut_edge_length)
1174  ents_to_check.insert(std::pair<double, EntityHandle>(max_dist, f));
1175  }
1176  }
1177 
1178  double ray_point[3], unit_ray_dir[3];
1179  VectorAdaptor vec_unit_ray_dir(
1180  3, ublas::shallow_array_adaptor<double>(3, unit_ray_dir));
1181  VectorAdaptor vec_ray_point(
1182  3, ublas::shallow_array_adaptor<double>(3, ray_point));
1183 
1184  for (auto m : ents_to_check) {
1185 
1186  EntityHandle f = m.second;
1187 
1188  int num_nodes;
1189  const EntityHandle *conn;
1190  CHKERR get_conn(f, conn, num_nodes);
1191  VectorDouble9 coords(9);
1192  CHKERR moab.get_coords(conn, num_nodes, &coords[0]);
1193 
1194  Range adj_tets;
1195  CHKERR moab.get_adjacencies(conn, num_nodes, 3, false, adj_tets,
1196  moab::Interface::UNION);
1197  adj_tets = intersect(adj_tets, vOlume);
1198 
1199  map<EntityHandle, TreeData> vertices_on_cut_edges;
1200  for (int n = 0; n != num_nodes; ++n)
1201  CHKERR project_node(conn[n], vertices_on_cut_edges);
1202 
1203  const double q = get_quality_change(adj_tets, vertices_on_cut_edges);
1204 
1206  EntityHandle type = moab.type_from_handle(f);
1207 
1208  Range zero_dist_tets;
1209  if (type == MBVERTEX) {
1210  Range zero_distance_verts_test = zeroDistanceVerts;
1211  zero_distance_verts_test.insert(f);
1212  CHKERR remove_surface_tets(zero_dist_tets, zeroDistanceEnts,
1213  zero_distance_verts_test);
1214  } else {
1215  Range zero_distance_ents_test = zeroDistanceEnts;
1216  zero_distance_ents_test.insert(f);
1217  CHKERR remove_surface_tets(zero_dist_tets, zero_distance_ents_test,
1219  }
1220 
1221  if (zero_dist_tets.empty()) {
1222  verticesOnCutEdges.insert(vertices_on_cut_edges.begin(),
1223  vertices_on_cut_edges.end());
1224 
1225  if (type == MBVERTEX) {
1226  zeroDistanceVerts.insert(f);
1227  } else {
1228  zeroDistanceEnts.insert(f);
1229  }
1230  }
1231  }
1232  }
1233  }
1234 
1235  for (auto &v : verticesOnCutEdges) {
1236 
1237  TYPE node_type;
1238 
1239  if (corner_nodes->find(v.first) != corner_nodes->end())
1240  node_type = CORNER;
1241  else if (fixed_edges_nodes.find(v.first) != fixed_edges_nodes.end())
1242  node_type = FIXED;
1243  else if (tets_skin_verts.find(v.first) != tets_skin_verts.end())
1244  node_type = SKIN;
1245  else
1246  node_type = FREE;
1247 
1248  if (node_type > min_dist_map[v.first].first.first)
1249  v.second.unitRayDir.clear();
1250  }
1251 
1252  for (auto f : unite(zeroDistanceEnts, zeroDistanceVerts)) {
1253  int num_nodes;
1254  const EntityHandle *conn;
1255  CHKERR get_conn(f, conn, num_nodes);
1256  Range adj_edges;
1257  CHKERR moab.get_adjacencies(conn, num_nodes, 1, false, adj_edges,
1258  moab::Interface::UNION);
1259  for (auto e : adj_edges) {
1260  cutEdges.erase(e);
1261  edgesToCut.erase(e);
1262  }
1263  }
1264 
1265  if (debug)
1266  SaveData(m_field.get_moab())("projected_zero_distance_ents.vtk",
1268 
1270 }
1271 
1273  const bool debug) {
1274  CoreInterface &m_field = cOre;
1275  moab::Interface &moab = m_field.get_moab();
1276  MeshRefinement *refiner;
1277  const RefEntity_multiIndex *ref_ents_ptr;
1279  if (cutEdges.size() != edgesToCut.size()) {
1280  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
1281  }
1282  CHKERR m_field.getInterface(refiner);
1283  CHKERR m_field.get_ref_ents(&ref_ents_ptr);
1285  CHKERR refiner->refine_TET(vOlume, bit, false, QUIET,
1286  debug ? &cutEdges : NULL);
1287  cutNewVolumes.clear();
1288  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1289  bit, BitRefLevel().set(), MBTET, cutNewVolumes);
1290  cutNewSurfaces.clear();
1291  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1292  bit, BitRefLevel().set(), MBTRI, cutNewSurfaces);
1293  // Find new vertices on cut edges
1294  cutNewVertices.clear();
1295  CHKERR moab.get_connectivity(zeroDistanceEnts, cutNewVertices, true);
1297  for (auto &m : edgesToCut) {
1298  RefEntity_multiIndex::index<
1299  Composite_ParentEnt_And_EntType_mi_tag>::type::iterator vit =
1300  ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().find(
1301  boost::make_tuple(MBVERTEX, m.first));
1302  if (vit ==
1303  ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().end()) {
1304  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1305  "No vertex on cut edges, that make no sense");
1306  }
1307  const boost::shared_ptr<RefEntity> &ref_ent = *vit;
1308  if ((ref_ent->getBitRefLevel() & bit).any()) {
1309  EntityHandle vert = ref_ent->getRefEnt();
1310  cutNewVertices.insert(vert);
1311  verticesOnCutEdges[vert] = m.second;
1312  } else {
1313  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1314  "Vertex has wrong bit ref level");
1315  }
1316  }
1317 
1318  // Add zero distance entities faces
1319  Range tets_skin;
1320  Skinner skin(&moab);
1321  CHKERR skin.find_skin(0, cutNewVolumes, false, tets_skin);
1322  cutNewSurfaces.merge(zeroDistanceEnts.subset_by_type(MBTRI));
1323 
1324  // At that point cutNewSurfaces has all newly created faces, now take all
1325  // nodes on those faces and subtract nodes on catted edges. Faces adjacent to
1326  // nodes which left are not part of surface.
1327  Range diff_verts;
1328  CHKERR moab.get_connectivity(unite(cutNewSurfaces, zeroDistanceEnts),
1329  diff_verts, true);
1330  diff_verts = subtract(diff_verts, cutNewVertices);
1331  Range subtract_faces;
1332  CHKERR moab.get_adjacencies(diff_verts, 2, false, subtract_faces,
1333  moab::Interface::UNION);
1334  cutNewSurfaces = subtract(cutNewSurfaces, unite(subtract_faces, tets_skin));
1335  cutNewVertices.clear();
1336  CHKERR moab.get_connectivity(cutNewSurfaces, cutNewVertices, true);
1337 
1339 }
1340 
1343 
1344  CoreInterface &m_field = cOre;
1345  moab::Interface &moab = m_field.get_moab();
1347 
1348  // Range out_side_vertices;
1349  for (auto m : verticesOnCutEdges) {
1350  double dist = m.second.dIst;
1351  VectorDouble3 new_coors = m.second.rayPoint + dist * m.second.unitRayDir;
1352  if (th)
1353  CHKERR moab.tag_set_data(th, &m.first, 1, &new_coors[0]);
1354  else
1355  CHKERR moab.set_coords(&m.first, 1, &new_coors[0]);
1356  }
1357 
1359 }
1360 
1362  CoreInterface &m_field = cOre;
1363  moab::Interface &moab = m_field.get_moab();
1365  // Range out_side_vertices;
1366  for (auto &v : verticesOnTrimEdges) {
1367  double dist = v.second.dIst;
1368  VectorDouble3 new_coors = v.second.rayPoint + dist * v.second.unitRayDir;
1369  if (th) {
1370  CHKERR moab.tag_set_data(th, &v.first, 1, &new_coors[0]);
1371  } else {
1372  CHKERR moab.set_coords(&v.first, 1, &new_coors[0]);
1373  }
1374  }
1376 }
1377 
1379  Range *corner_nodes, Tag th,
1380  const double tol,
1381  const double tol_close,
1382  const bool debug) {
1383  CoreInterface &m_field = cOre;
1384  moab::Interface &moab = m_field.get_moab();
1386 
1387  // takes body skin
1388  Skinner skin(&moab);
1389  Range tets_skin;
1390  CHKERR skin.find_skin(0, cutNewVolumes, false, tets_skin);
1391  // vertices on the skin
1392  Range tets_skin_verts;
1393  CHKERR moab.get_connectivity(tets_skin, tets_skin_verts, true);
1394  // edges on the skin
1395  Range tets_skin_edges;
1396  CHKERR moab.get_adjacencies(tets_skin, 1, false, tets_skin_edges,
1397  moab::Interface::UNION);
1398  // get edges on new surface
1399  Range edges;
1400  CHKERR moab.get_adjacencies(cutNewSurfaces, 1, false, edges,
1401  moab::Interface::UNION);
1402 
1403  Range cut_surface_edges_on_fixed_edges;
1404  if (fixed_edges)
1405  cut_surface_edges_on_fixed_edges = intersect(edges, *fixed_edges);
1406 
1407  Range cut_surface_edges_on_fixed_edges_verts;
1408  CHKERR moab.get_connectivity(cut_surface_edges_on_fixed_edges,
1409  cut_surface_edges_on_fixed_edges_verts, true);
1410 
1411  Range fixed_edges_verts;
1412  if (fixed_edges)
1413  CHKERR moab.get_connectivity(*fixed_edges, fixed_edges_verts, true);
1414 
1415  Range surface_skin;
1416  if (fRont.empty())
1417  CHKERR skin.find_skin(0, sUrface, false, surface_skin);
1418  else
1419  surface_skin = fRont;
1420 
1421  auto get_point_coords = [&th, &moab](EntityHandle v) {
1422  VectorDouble3 point(3);
1423  if (th)
1424  CHKERR moab.tag_get_data(th, &v, 1, &point[0]);
1425  else
1426  CHKERR moab.get_coords(&v, 1, &point[0]);
1427  return point;
1428  };
1429 
1430  struct VertMap {
1431  const double d;
1432  const EntityHandle v;
1433  const EntityHandle e;
1434  VertMap(const double d, const EntityHandle v, const EntityHandle e)
1435  : d(d), v(v), e(e) {}
1436  };
1437 
1438  typedef multi_index_container<
1439  VertMap,
1440  indexed_by<
1441  ordered_non_unique<member<VertMap, const double, &VertMap::d>>,
1442  ordered_non_unique<member<VertMap, const EntityHandle, &VertMap::v>>
1443 
1444  >>
1445  VertMapMultiIndex;
1446 
1447  VertMapMultiIndex verts_map;
1448 
1449  auto add_vert = [&](const EntityHandle v, const EntityHandle e,
1450  const double dist) {
1451  verts_map.insert(VertMap(dist, v, e));
1452  };
1453 
1454  // clear data ranges
1455  trimEdges.clear();
1456  edgesToTrim.clear();
1457  verticesOnTrimEdges.clear();
1458  trimNewVertices.clear();
1459 
1460  if (debug)
1461  CHKERR SaveData(moab)("edges_potentially_to_trim.vtk", edges);
1462 
1463  if (debug)
1464  CHKERR SaveData(moab)("surface_skin_to_trim.vtk", surface_skin);
1465 
1466  auto cut_this_edge = [&](const EntityHandle e, const double length, auto &ray,
1467  auto &ray_point) {
1468  trimEdges.insert(e);
1469  edgesToTrim[e].dIst = 1;
1470  edgesToTrim[e].lEngth = length;
1471  edgesToTrim[e].unitRayDir.resize(3, false);
1472  edgesToTrim[e].rayPoint.resize(3, false);
1473  for (int ii = 0; ii != 3; ++ii)
1474  edgesToTrim[e].unitRayDir[ii] = ray(ii);
1475  for (int ii = 0; ii != 3; ++ii)
1476  edgesToTrim[e].rayPoint[ii] = ray_point(ii);
1477  };
1478 
1480  int num_nodes;
1481 
1482  MatrixDouble edge_face_normal(3, surface_skin.size());
1483  FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 3> t_edge_face_normal =
1484  getFTensor1FromMat<3>(edge_face_normal);
1485  for (auto s : surface_skin) {
1486  Range adj_face;
1487  moab.get_adjacencies(&s, 1, 2, false, adj_face, moab::Interface::UNION);
1488  adj_face = intersect(adj_face, sUrface);
1489  if (adj_face.size() == 1)
1490  Util::normal(&moab, adj_face[0], t_edge_face_normal(0),
1491  t_edge_face_normal(1), t_edge_face_normal(2));
1492  else
1493  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1494  "Should be only one face in the range");
1495 
1496  t_edge_face_normal(i) /=
1497  sqrt(t_edge_face_normal(i) * t_edge_face_normal(i));
1498  ++t_edge_face_normal;
1499  }
1500 
1501  for (auto e : edges) {
1502 
1503  // Get edge connectivity and coords
1504  const EntityHandle *conn_edge;
1505  CHKERR moab.get_connectivity(e, conn_edge, num_nodes, true);
1506  double coords_edge[3 * num_nodes];
1507  CHKERR moab.get_coords(conn_edge, num_nodes, coords_edge);
1508 
1509  FTensor::Tensor1<double *, 3> t_e0(&coords_edge[0], &coords_edge[1],
1510  &coords_edge[2]);
1511  FTensor::Tensor1<double *, 3> t_e1(&coords_edge[3], &coords_edge[4],
1512  &coords_edge[5]);
1513 
1514  FTensor::Tensor1<double, 3> t_edge_delta;
1515  t_edge_delta(i) = t_e1(i) - t_e0(i);
1516  const double edge_length2 = t_edge_delta(i) * t_edge_delta(i);
1517  const double edge_length = sqrt(edge_length2);
1518  if (edge_length == 0)
1519  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Zero edge length");
1520 
1521  FTensor::Tensor1<FTensor::PackPtr<double *, 1>, 3> t_edge_face_normal =
1522  getFTensor1FromMat<3>(edge_face_normal);
1523 
1524  auto t_project = [&](auto &t_vec) {
1526  t_prj(i) =
1527  t_vec(i) - (t_edge_face_normal(i) * t_vec(i)) * t_edge_face_normal(i);
1528  return t_prj;
1529  };
1530 
1531  auto get_edge_coors = [&](const auto e) {
1532  const EntityHandle *conn;
1533  CHKERR moab.get_connectivity(e, conn, num_nodes, true);
1534  VectorDouble6 coords(6);
1535  if (th)
1536  CHKERR moab.tag_get_data(th, conn, num_nodes, &coords[0]);
1537  else
1538  CHKERR moab.get_coords(conn, num_nodes, &coords[0]);
1539  return coords;
1540  };
1541 
1542  // iterate entities on inner surface font to find that edge cross
1543  for (auto s : surface_skin) {
1544 
1545  auto coords_front = get_edge_coors(s);
1546 
1547  FTensor::Tensor1<double *, 3> t_f0(&coords_front[0], &coords_front[1],
1548  &coords_front[2]);
1549  FTensor::Tensor1<double *, 3> t_f1(&coords_front[3], &coords_front[4],
1550  &coords_front[5]);
1551 
1552  auto t_e0_prj = t_project(t_e0);
1553  auto t_e1_prj = t_project(t_e1);
1554  auto t_f0_prj = t_project(t_f0);
1555  auto t_f1_prj = t_project(t_f1);
1556 
1557  double t_edge, t_front;
1558  auto res = Tools::minDistanceFromSegments(&t_e0_prj(0), &t_e1_prj(0),
1559  &t_f0_prj(0), &t_f1_prj(0),
1560  &t_edge, &t_front);
1561 
1562  if (res != Tools::NO_SOLUTION) {
1563 
1564  const double overlap_tol = 1e-2;
1565  if (t_edge >= 0 && t_edge <= 1 && t_front >= -overlap_tol &&
1566  t_front <= 1 + overlap_tol) {
1567 
1568  FTensor::Tensor1<double, 3> t_front_delta_prj;
1569  t_front_delta_prj(i) = t_f1_prj(i) - t_f0_prj(i);
1570  FTensor::Tensor1<double, 3> t_edge_delta_prj;
1571  t_edge_delta_prj(i) = t_e1_prj(i) - t_e0_prj(i);
1572  FTensor::Tensor1<double, 3> t_edge_point_prj, t_front_point_prj;
1573  t_edge_point_prj(i) = t_e0_prj(i) + t_edge * t_edge_delta_prj(i);
1574  t_front_point_prj(i) = t_f0_prj(i) + t_front * t_front_delta_prj(i);
1575  FTensor::Tensor1<double, 3> t_ray_prj;
1576  t_ray_prj(i) = t_front_point_prj(i) - t_edge_point_prj(i);
1577  const double dist = sqrt(t_ray_prj(i) * t_ray_prj(i));
1578 
1579  if ((dist / edge_length) < 0.1) {
1580 
1581  auto check_to_add_edge = [&](const EntityHandle e,
1582  const double dist) {
1583  auto eit = edgesToTrim.find(e);
1584  if (eit != edgesToTrim.end())
1585  if (eit->second.dIst < dist)
1586  return false;
1587  return true;
1588  };
1589 
1591  if (t_edge < 0.5) {
1592  t_ray(i) = t_edge * t_edge_delta(i);
1593  const double ray_length = sqrt(t_ray(i) * t_ray(i));
1594  if (check_to_add_edge(e, ray_length)) {
1595  add_vert(conn_edge[0], e, fabs(t_edge));
1596  add_vert(conn_edge[1], e, fabs(t_edge - 1));
1597  cut_this_edge(e, edge_length, t_ray, t_e0);
1598  }
1599  } else {
1600  FTensor::Tensor1<double, 3> t_edge_point;
1601  t_edge_point(i) = t_e0(i) + t_edge * t_edge_delta(i);
1602  t_ray(i) = t_edge_point(i) - t_e1(i);
1603  const double ray_length = sqrt(t_ray(i) * t_ray(i));
1604  if (check_to_add_edge(e, ray_length)) {
1605  add_vert(conn_edge[0], e, fabs(t_edge));
1606  add_vert(conn_edge[1], e, fabs(t_edge - 1));
1607  cut_this_edge(e, edge_length, t_ray, t_e1);
1608  }
1609  }
1610  }
1611  }
1612  }
1613 
1614  ++t_edge_face_normal;
1615  }
1616  }
1617 
1618  if (debug)
1619  CHKERR SaveData(moab)("edges_selected_to_trim.vtk", trimEdges);
1620 
1621  Range corners;
1622  if (corner_nodes)
1623  corners = *corner_nodes;
1624  Range fix_edges;
1625  if (fixed_edges)
1626  fix_edges = *fixed_edges;
1627 
1628  // Iterate over all vertives close to surface front and check if those can be
1629  // moved
1630  auto lo_m = verts_map.get<0>().lower_bound(0);
1631  auto hi_m = verts_map.get<0>().lower_bound(tol);
1632  for (; lo_m != hi_m; ++lo_m) {
1633 
1634  const double dist = lo_m->d;
1635  EntityHandle v = lo_m->v;
1636  if (verticesOnTrimEdges.find(v) == verticesOnTrimEdges.end()) {
1637 
1638  VectorDouble3 ray_point = get_point_coords(v);
1639  EntityHandle e = lo_m->e;
1640 
1641  VectorDouble3 new_pos(3);
1642  new_pos.clear();
1643 
1644  auto get_quality_change = [&](const Range &adj_tets) {
1645  double q0 = 1;
1646  CHKERR m_field.getInterface<Tools>()->minTetsQuality(adj_tets, q0, th);
1647 
1648  double min_q = 1;
1649  for (auto t : adj_tets) {
1650  int num_nodes;
1651  const EntityHandle *conn;
1652  CHKERR m_field.get_moab().get_connectivity(t, conn, num_nodes, true);
1653  VectorDouble12 coords(12);
1654  if (th)
1655  CHKERR moab.tag_get_data(th, conn, num_nodes, &coords[0]);
1656  else
1657  CHKERR moab.get_coords(conn, num_nodes, &coords[0]);
1658 
1659  for (int n = 0; n != 4; ++n) {
1660  auto n_coords = getVectorAdaptor(&coords[3 * n], 3);
1661  if (conn[n] == v) {
1662  noalias(n_coords) = new_pos;
1663  } else {
1664  auto m = verticesOnTrimEdges.find(conn[n]);
1665  if (m != verticesOnTrimEdges.end())
1666  noalias(n_coords) =
1667  m->second.rayPoint + m->second.dIst * m->second.unitRayDir;
1668  }
1669  }
1670 
1671  const double q = Tools::volumeLengthQuality(&coords[0]);
1672  if (!std::isnormal(q)) {
1673  min_q = -2;
1674  break;
1675  }
1676  min_q = std::min(min_q, q);
1677  }
1678  return min_q / q0;
1679  };
1680 
1681  auto add_trim_vert = [&]() {
1682  auto &r = edgesToTrim.at(e);
1683  noalias(new_pos) = r.rayPoint + r.dIst * r.unitRayDir;
1684  VectorDouble3 unit_ray_dir = ray_point - new_pos;
1685 
1686  Range adj_tets;
1687  CHKERR moab.get_adjacencies(&v, 1, 3, false, adj_tets);
1688  adj_tets = intersect(adj_tets, cutNewVolumes);
1689 
1690  double q = get_quality_change(adj_tets);
1692  VectorDouble3 ray_dir = new_pos - ray_point;
1693  double dist = norm_2(ray_dir);
1694  verticesOnTrimEdges[v].dIst = 1;
1695  verticesOnTrimEdges[v].lEngth = dist;
1696  verticesOnTrimEdges[v].unitRayDir = ray_dir;
1697  verticesOnTrimEdges[v].rayPoint = ray_point;
1698  trimNewVertices.insert(v);
1699  }
1700  };
1701 
1702  if (trimNewVertices.find(v) == trimNewVertices.end()) {
1703  if (edgesToTrim.find(e) != edgesToTrim.end()) {
1704  if (dist < tol_close) {
1705 
1706  int rank_vertex = 0;
1707  if (corners.find(v) != corners.end())
1708  ++rank_vertex;
1709  if (fixed_edges_verts.find(v) != fixed_edges_verts.end())
1710  ++rank_vertex;
1711  if (tets_skin_verts.find(v) != tets_skin_verts.end())
1712  ++rank_vertex;
1713 
1714  int rank_edge = 0;
1715  if (fix_edges.find(e) != fix_edges.end())
1716  ++rank_edge;
1717  if (tets_skin_edges.find(e) != tets_skin_edges.end())
1718  ++rank_edge;
1719 
1720  if (rank_vertex > rank_edge) {
1721  auto &m = edgesToTrim.at(e);
1722  verticesOnTrimEdges[v] = m;
1723  verticesOnTrimEdges[v].rayPoint = get_point_coords(v);
1724  verticesOnTrimEdges[v].lEngth = 0;
1725  verticesOnTrimEdges[v].dIst = 0;
1726  trimNewVertices.insert(v);
1727  } else if (rank_vertex == rank_edge)
1728  add_trim_vert();
1729  else
1730  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1731  "Imposible case");
1732 
1733  } else if (corners.find(v) != corners.end()) {
1734  auto &m = edgesToTrim.at(e);
1735  verticesOnTrimEdges[v] = m;
1736  verticesOnTrimEdges[v].rayPoint = get_point_coords(v);
1737  verticesOnTrimEdges[v].lEngth = 0;
1738  verticesOnTrimEdges[v].dIst = 0;
1739  trimNewVertices.insert(v);
1740  } else if (fixed_edges_verts.find(v) != fixed_edges_verts.end()) {
1741  if (fix_edges.find(e) != fix_edges.end()) {
1742  add_trim_vert();
1743  }
1744  } else if (tets_skin_verts.find(v) != tets_skin_verts.end()) {
1745  if (tets_skin_edges.find(e) != tets_skin_edges.end()) {
1746  add_trim_vert();
1747  }
1748  } else {
1749  add_trim_vert();
1750  }
1751  }
1752  }
1753  }
1754  }
1755 
1756  for (auto m : verticesOnTrimEdges) {
1757  EntityHandle v = m.first;
1758  // Range adj_vertex_edges;
1759  // CHKERR moab.get_adjacencies(&v, 1, 1, false, adj_vertex_edges,
1760  // moab::Interface::UNION);
1761  // for(auto e : adj_vertex_edges) {
1762  // edgesToTrim.erase(e);
1763  // trimEdges.erase(e);
1764  // }
1765  auto lo = verts_map.get<1>().lower_bound(v);
1766  auto hi = verts_map.get<1>().upper_bound(v);
1767  for (; lo != hi; ++lo) {
1768  edgesToTrim.erase(lo->e);
1769  trimEdges.erase(lo->e);
1770  }
1771  }
1772 
1773  if (!trimNewVertices.empty()) {
1774  Range adj_trim_vertex_edges;
1775  CHKERR moab.get_adjacencies(trimNewVertices, 1, false,
1776  adj_trim_vertex_edges, moab::Interface::UNION);
1777  adj_trim_vertex_edges = intersect(adj_trim_vertex_edges, trimEdges);
1778  Range adj_trim_vertex_edges_verts;
1779  CHKERR moab.get_connectivity(adj_trim_vertex_edges,
1780  adj_trim_vertex_edges_verts, true);
1781  adj_trim_vertex_edges_verts =
1782  subtract(adj_trim_vertex_edges_verts, trimNewVertices);
1783  Range edges_connected_to_close_verts;
1784  CHKERR moab.get_adjacencies(adj_trim_vertex_edges_verts, 1, false,
1785  edges_connected_to_close_verts,
1786  moab::Interface::UNION);
1787  edges_connected_to_close_verts =
1788  subtract(adj_trim_vertex_edges, edges_connected_to_close_verts);
1789 
1790  if (debug)
1791  if (!edges_connected_to_close_verts.empty())
1792  CHKERR SaveData(moab)("trim_edges_connected_to_close_verts.vtk",
1793  edges_connected_to_close_verts);
1794 
1795  for (auto e : edges_connected_to_close_verts) {
1796  edgesToTrim.erase(e);
1797  trimEdges.erase(e);
1798  }
1799  }
1800 
1801  if (debug)
1802  if (!trimNewVertices.empty())
1803  CHKERR SaveData(moab)("trim_close_vertices.vtk", trimNewVertices);
1804 
1805  if (debug)
1806  if (!trimEdges.empty())
1807  CHKERR SaveData(moab)("trim_edges.vtk", trimEdges);
1808 
1810 } // namespace MoFEM
1811 
1813  Tag th, const double tol,
1814  const bool debug) {
1815  CoreInterface &m_field = cOre;
1816  moab::Interface &moab = m_field.get_moab();
1817  MeshRefinement *refiner;
1818  const RefEntity_multiIndex *ref_ents_ptr;
1820 
1821  CHKERR m_field.getInterface(refiner);
1822  CHKERR m_field.get_ref_ents(&ref_ents_ptr);
1824  CHKERR refiner->refine_TET(cutNewVolumes, bit, false, QUIET,
1825  debug ? &trimEdges : NULL);
1826 
1827  trimNewVolumes.clear();
1828  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1829  bit, BitRefLevel().set(), MBTET, trimNewVolumes);
1830 
1831  for (map<EntityHandle, TreeData>::iterator mit = edgesToTrim.begin();
1832  mit != edgesToTrim.end(); mit++) {
1833  auto vit = ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().find(
1834  boost::make_tuple(MBVERTEX, mit->first));
1835  if (vit ==
1836  ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().end()) {
1837  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1838  "No vertex on trim edges, that make no sense");
1839  }
1840  const boost::shared_ptr<RefEntity> &ref_ent = *vit;
1841  if ((ref_ent->getBitRefLevel() & bit).any()) {
1842  EntityHandle vert = ref_ent->getRefEnt();
1843  trimNewVertices.insert(vert);
1844  verticesOnTrimEdges[vert] = mit->second;
1845  }
1846  }
1847 
1848  // Get faces which are trimmed
1849  trimNewSurfaces.clear();
1850  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1851  bit, BitRefLevel().set(), MBTRI, trimNewSurfaces);
1852 
1853  Range trim_new_surfaces_nodes;
1854  CHKERR moab.get_connectivity(trimNewSurfaces, trim_new_surfaces_nodes, true);
1855  trim_new_surfaces_nodes = subtract(trim_new_surfaces_nodes, trimNewVertices);
1856  trim_new_surfaces_nodes = subtract(trim_new_surfaces_nodes, cutNewVertices);
1857  Range faces_not_on_surface;
1858  CHKERR moab.get_adjacencies(trim_new_surfaces_nodes, 2, false,
1859  faces_not_on_surface, moab::Interface::UNION);
1860  trimNewSurfaces = subtract(trimNewSurfaces, faces_not_on_surface);
1861 
1862  // Get surfaces which are not trimmed and add them to surface
1863  Range all_surfaces_on_bit_level;
1864  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1865  bit, BitRefLevel().set(), MBTRI, all_surfaces_on_bit_level);
1866  all_surfaces_on_bit_level =
1867  intersect(all_surfaces_on_bit_level, cutNewSurfaces);
1868  trimNewSurfaces = unite(trimNewSurfaces, all_surfaces_on_bit_level);
1869 
1870  Range trim_surface_edges;
1871  CHKERR moab.get_adjacencies(trimNewSurfaces, 1, false, trim_surface_edges,
1872  moab::Interface::UNION);
1873 
1874  // check of nodes are outside surface and if it are remove adjacent faces to
1875  // those nodes.
1876  Range check_verts;
1877  CHKERR moab.get_connectivity(trimNewSurfaces, check_verts, true);
1878  check_verts = subtract(check_verts, trimNewVertices);
1879  for (auto v : check_verts) {
1880 
1881  VectorDouble3 s(3);
1882  if (th) {
1883  CHKERR moab.tag_get_data(th, &v, 1, &s[0]);
1884  } else {
1885  CHKERR moab.get_coords(&v, 1, &s[0]);
1886  }
1887 
1888  VectorDouble3 p(3);
1889  EntityHandle facets_out;
1890  CHKERR treeSurfPtr->closest_to_location(&s[0], rootSetSurf, &p[0],
1891  facets_out);
1892  VectorDouble3 n(3);
1893  Util::normal(&moab, facets_out, n[0], n[1], n[2]);
1894  n /= norm_2(n);
1895  VectorDouble3 delta = s - p;
1896  VectorDouble3 normal = inner_prod(delta, n) * n;
1897  if (norm_2(delta - normal) > tol * aveLength) {
1898  Range adj;
1899  CHKERR moab.get_adjacencies(&v, 1, 2, false, adj);
1900  trimNewSurfaces = subtract(trimNewSurfaces, adj);
1901  }
1902  }
1903 
1905 }
1906 
1909  Range &ents) {
1910  CoreInterface &m_field = cOre;
1911  moab::Interface &moab = m_field.get_moab();
1912  PrismInterface *interface;
1914  CHKERR m_field.getInterface(interface);
1915  // Remove tris on surface front
1916  {
1917  Range front_tris;
1918  EntityHandle meshset;
1919  CHKERR moab.create_meshset(MESHSET_SET, meshset);
1920  CHKERR moab.add_entities(meshset, ents);
1922  meshset, split_bit, true, front_tris);
1923  CHKERR moab.delete_entities(&meshset, 1);
1924  ents = subtract(ents, front_tris);
1925  }
1926  // Remove entities on skin
1927  Range tets;
1928  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1929  split_bit, BitRefLevel().set(), MBTET, tets);
1930  // Remove entities on skin
1931  Skinner skin(&moab);
1932  Range tets_skin;
1933  rval = skin.find_skin(0, tets, false, tets_skin);
1934  ents = subtract(ents, tets_skin);
1935 
1937 }
1938 
1940  const BitRefLevel bit,
1941  const Range &ents, Tag th) {
1942  CoreInterface &m_field = cOre;
1943  moab::Interface &moab = m_field.get_moab();
1944  PrismInterface *interface;
1946  CHKERR m_field.getInterface(interface);
1947  EntityHandle meshset_volume;
1948  CHKERR moab.create_meshset(MESHSET_SET, meshset_volume);
1949  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1950  split_bit, BitRefLevel().set(), MBTET, meshset_volume);
1951  EntityHandle meshset_surface;
1952  CHKERR moab.create_meshset(MESHSET_SET, meshset_surface);
1953  CHKERR moab.add_entities(meshset_surface, ents);
1954  CHKERR interface->getSides(meshset_surface, split_bit, true);
1955  CHKERR interface->splitSides(meshset_volume, bit, meshset_surface, true,
1956  true);
1957  CHKERR moab.delete_entities(&meshset_volume, 1);
1958  CHKERR moab.delete_entities(&meshset_surface, 1);
1959  if (th) {
1960  Range prisms;
1961  ierr = m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1962  bit, BitRefLevel().set(), MBPRISM, prisms);
1963  for (Range::iterator pit = prisms.begin(); pit != prisms.end(); pit++) {
1964  int num_nodes;
1965  const EntityHandle *conn;
1966  CHKERR moab.get_connectivity(*pit, conn, num_nodes, true);
1967  MatrixDouble data(3, 3);
1968  CHKERR moab.tag_get_data(th, conn, 3, &data(0, 0));
1969  CHKERR moab.tag_set_data(th, &conn[3], 3, &data(0, 0));
1970  }
1971  }
1973 }
1974 
1976  double lEngth;
1977  double qUality;
1979  bool skip;
1980  LengthMapData(const double l, double q, const EntityHandle e)
1981  : lEngth(l), qUality(q), eDge(e), skip(false) {}
1982 };
1983 
1984 typedef multi_index_container<
1985  LengthMapData,
1986  indexed_by<
1987 
1988  ordered_non_unique<
1989  member<LengthMapData, double, &LengthMapData::lEngth>>,
1990 
1991  hashed_unique<
1992  member<LengthMapData, EntityHandle, &LengthMapData::eDge>>,
1993 
1994  ordered_non_unique<
1995  member<LengthMapData, double, &LengthMapData::qUality>>
1996 
1997  >>
1999 
2001  const int fraction_level, const Range &tets, const Range &surface,
2002  const Range &fixed_edges, const Range &corner_nodes, Range &edges_to_merge,
2003  Range &out_tets, Range &new_surf, Tag th, const bool update_meshsets,
2004  const BitRefLevel *bit_ptr, const bool debug) {
2005  CoreInterface &m_field = cOre;
2006  moab::Interface &moab = m_field.get_moab();
2008 
2009  /**
2010  * \brief Merge nodes
2011  */
2012  struct MergeNodes {
2013  CoreInterface &mField;
2014  const bool onlyIfImproveQuality;
2015  const int lineSearch;
2016  Tag tH;
2017  bool updateMehsets;
2018 
2019  MergeNodes(CoreInterface &m_field, const bool only_if_improve_quality,
2020  const int line_search, Tag th, bool update_mehsets)
2021  : mField(m_field), onlyIfImproveQuality(only_if_improve_quality),
2022  lineSearch(line_search), tH(th), updateMehsets(update_mehsets) {
2023  mField.getInterface(nodeMergerPtr);
2024  }
2025  NodeMergerInterface *nodeMergerPtr;
2026  MoFEMErrorCode operator()(EntityHandle father, EntityHandle mother,
2027  Range &proc_tets, Range &new_surf,
2028  Range &edges_to_merge, Range &not_merged_edges,
2029  bool add_child = true) const {
2030  moab::Interface &moab = mField.get_moab();
2032  const EntityHandle conn[] = {father, mother};
2033  Range vert_tets;
2034  CHKERR moab.get_adjacencies(conn, 2, 3, false, vert_tets,
2035  moab::Interface::UNION);
2036  vert_tets = intersect(vert_tets, proc_tets);
2037  Range out_tets;
2038  CHKERR nodeMergerPtr->mergeNodes(father, mother, out_tets, &vert_tets,
2039  onlyIfImproveQuality, 0, lineSearch, tH);
2040  out_tets.merge(subtract(proc_tets, vert_tets));
2041  proc_tets.swap(out_tets);
2042 
2043  if (add_child && nodeMergerPtr->getSucessMerge()) {
2044 
2045  auto &parent_child_map = nodeMergerPtr->getParentChildMap();
2046 
2047  Range child_ents;
2048  for (auto &it : parent_child_map)
2049  child_ents.insert(it.pArent);
2050 
2051  Range new_surf_child_ents = intersect(new_surf, child_ents);
2052  new_surf = subtract(new_surf, new_surf_child_ents);
2053  Range child_surf_ents;
2054  CHKERR updateRangeByChilds(parent_child_map, new_surf_child_ents,
2055  child_surf_ents);
2056  new_surf.merge(child_surf_ents);
2057 
2058  Range edges_to_merge_child_ents = intersect(edges_to_merge, child_ents);
2059  edges_to_merge = subtract(edges_to_merge, edges_to_merge_child_ents);
2060  Range merged_child_edge_ents;
2061  CHKERR updateRangeByChilds(parent_child_map, edges_to_merge_child_ents,
2062  merged_child_edge_ents);
2063 
2064  Range not_merged_edges_child_ents =
2065  intersect(not_merged_edges, child_ents);
2066  not_merged_edges =
2067  subtract(not_merged_edges, not_merged_edges_child_ents);
2068  Range not_merged_child_edge_ents;
2069  CHKERR updateRangeByChilds(parent_child_map,
2070  not_merged_edges_child_ents,
2071  not_merged_child_edge_ents);
2072 
2073  merged_child_edge_ents =
2074  subtract(merged_child_edge_ents, not_merged_child_edge_ents);
2075  edges_to_merge.merge(merged_child_edge_ents);
2076  not_merged_edges.merge(not_merged_child_edge_ents);
2077 
2078  if (updateMehsets) {
2080  (*mField.getInterface<MeshsetsManager>()), cubit_it)) {
2081  EntityHandle cubit_meshset = cubit_it->meshset;
2082  Range parent_ents;
2083  CHKERR moab.get_entities_by_handle(cubit_meshset, parent_ents,
2084  true);
2085  Range child_ents;
2086  CHKERR updateRangeByChilds(parent_child_map, parent_ents,
2087  child_ents);
2088  CHKERR moab.add_entities(cubit_meshset, child_ents);
2089  }
2090  }
2091  }
2093  }
2094 
2095  private:
2096  MoFEMErrorCode updateRangeByChilds(
2097  const NodeMergerInterface::ParentChildMap &parent_child_map,
2098  const Range &parents, Range &childs) const {
2100  NodeMergerInterface::ParentChildMap::nth_index<0>::type::iterator it;
2101  for (Range::const_iterator eit = parents.begin(); eit != parents.end();
2102  eit++) {
2103  it = parent_child_map.get<0>().find(*eit);
2104  if (it == parent_child_map.get<0>().end())
2105  continue;
2106  childs.insert(it->cHild);
2107  }
2109  }
2110  };
2111 
2112  /**
2113  * \brief Calculate edge length
2114  */
2115  struct LengthMap {
2116  Tag tH;
2117  CoreInterface &mField;
2118  moab::Interface &moab;
2119  const double maxLength;
2120  LengthMap(CoreInterface &m_field, Tag th, double max_length)
2121  : tH(th), mField(m_field), moab(m_field.get_moab()),
2122  maxLength(max_length) {
2123  gammaL = 1.;
2124  gammaQ = 3.;
2125  acceptedThrasholdMergeQuality = 0.5;
2126  }
2127 
2128  double
2129  gammaL; ///< Controls importance of length when ranking edges for merge
2130  double
2131  gammaQ; ///< Controls importance of quality when ranking edges for merge
2132  double acceptedThrasholdMergeQuality; ///< Do not merge quality if quality
2133  ///< above accepted thrashold
2134 
2135  MoFEMErrorCode operator()(const Range &tets, const Range &edges,
2136  LengthMapData_multi_index &length_map,
2137  double &ave) const {
2138  int num_nodes;
2139  const EntityHandle *conn;
2140  double coords[6];
2142  VectorAdaptor s0(3, ublas::shallow_array_adaptor<double>(3, &coords[0]));
2143  VectorAdaptor s1(3, ublas::shallow_array_adaptor<double>(3, &coords[3]));
2144  VectorDouble3 delta(3);
2145  for (auto edge : edges) {
2146  CHKERR moab.get_connectivity(edge, conn, num_nodes, true);
2147  Range adj_tet;
2148  CHKERR moab.get_adjacencies(conn, num_nodes, 3, false, adj_tet);
2149  adj_tet = intersect(adj_tet, tets);
2150  if (tH) {
2151  CHKERR moab.tag_get_data(tH, conn, num_nodes, coords);
2152  } else {
2153  CHKERR moab.get_coords(conn, num_nodes, coords);
2154  }
2155  double q = 1;
2156  CHKERR mField.getInterface<Tools>()->minTetsQuality(adj_tet, q, tH);
2157  if (q < acceptedThrasholdMergeQuality) {
2158  noalias(delta) = (s0 - s1) / maxLength;
2159  double dot = inner_prod(delta, delta);
2160  double val = pow(q, gammaQ) * pow(dot, 0.5 * gammaL);
2161  length_map.insert(LengthMapData(val, q, edge));
2162  }
2163  }
2164  ave = 0;
2165  for (LengthMapData_multi_index::nth_index<0>::type::iterator mit =
2166  length_map.get<0>().begin();
2167  mit != length_map.get<0>().end(); mit++) {
2168  ave += mit->qUality;
2169  }
2170  ave /= length_map.size();
2172  }
2173  };
2174 
2175  /**
2176  * \brief Topological relations
2177  */
2178  struct Toplogy {
2179 
2180  CoreInterface &mField;
2181  Tag tH;
2182  const double tOL;
2183  Toplogy(CoreInterface &m_field, Tag th, const double tol)
2184  : mField(m_field), tH(th), tOL(tol) {}
2185 
2186  enum TYPE {
2187  FREE = 0,
2188  SKIN = 1 << 0,
2189  SURFACE = 1 << 1,
2190  SURFACE_SKIN = 1 << 2,
2191  FRONT_ENDS = 1 << 3,
2192  FIX_EDGES = 1 << 4,
2193  FIX_CORNERS = 1 << 5
2194  };
2195 
2196  typedef map<int, Range> SetsMap;
2197 
2198  MoFEMErrorCode classifyVerts(const Range &surface, const Range &tets,
2199  const Range &fixed_edges,
2200  const Range &corner_nodes,
2201  SetsMap &sets_map) const {
2202  moab::Interface &moab(mField.get_moab());
2203  Skinner skin(&moab);
2205 
2206  sets_map[FIX_CORNERS].merge(corner_nodes);
2207  Range fixed_verts;
2208  CHKERR moab.get_connectivity(fixed_edges, fixed_verts, true);
2209  sets_map[FIX_EDGES].swap(fixed_verts);
2210 
2211  Range tets_skin;
2212  CHKERR skin.find_skin(0, tets, false, tets_skin);
2213  Range tets_skin_edges;
2214  CHKERR moab.get_adjacencies(tets_skin, 1, false, tets_skin_edges,
2215  moab::Interface::UNION);
2216 
2217  // surface skin
2218  Range surface_skin;
2219  CHKERR skin.find_skin(0, surface, false, surface_skin);
2220  Range front_in_the_body;
2221  front_in_the_body = subtract(surface_skin, tets_skin_edges);
2222  Range front_ends;
2223  CHKERR skin.find_skin(0, front_in_the_body, false, front_ends);
2224  sets_map[FRONT_ENDS].swap(front_ends);
2225 
2226  Range surface_skin_verts;
2227  CHKERR moab.get_connectivity(surface_skin, surface_skin_verts, true);
2228  sets_map[SURFACE_SKIN].swap(surface_skin_verts);
2229 
2230  // surface
2231  Range surface_verts;
2232  CHKERR moab.get_connectivity(surface, surface_verts, true);
2233  sets_map[SURFACE].swap(surface_verts);
2234 
2235  // skin
2236  Range tets_skin_verts;
2237  CHKERR moab.get_connectivity(tets_skin, tets_skin_verts, true);
2238  sets_map[SKIN].swap(tets_skin_verts);
2239 
2240  Range tets_verts;
2241  CHKERR moab.get_connectivity(tets, tets_verts, true);
2242  sets_map[FREE].swap(tets_verts);
2243 
2245  }
2246 
2247  MoFEMErrorCode getProcTets(const Range &tets, const Range &edges_to_merge,
2248  Range &proc_tets) const {
2249  moab::Interface &moab(mField.get_moab());
2251  Range edges_to_merge_verts;
2252  CHKERR moab.get_connectivity(edges_to_merge, edges_to_merge_verts, true);
2253  Range edges_to_merge_verts_tets;
2254  CHKERR moab.get_adjacencies(edges_to_merge_verts, 3, false,
2255  edges_to_merge_verts_tets,
2256  moab::Interface::UNION);
2257  edges_to_merge_verts_tets = intersect(edges_to_merge_verts_tets, tets);
2258  proc_tets.swap(edges_to_merge_verts_tets);
2260  }
2261 
2262  MoFEMErrorCode edgesToMerge(const Range &surface, const Range &tets,
2263  Range &edges_to_merge) const {
2264  moab::Interface &moab(mField.get_moab());
2266  Range surface_verts;
2267  CHKERR moab.get_connectivity(surface, surface_verts, true);
2268  Range surface_verts_edges;
2269  CHKERR moab.get_adjacencies(surface_verts, 1, false, surface_verts_edges,
2270  moab::Interface::UNION);
2271  edges_to_merge.merge(surface_verts_edges);
2272  Range tets_edges;
2273  CHKERR moab.get_adjacencies(tets, 1, false, tets_edges,
2274  moab::Interface::UNION);
2275  edges_to_merge = intersect(edges_to_merge, tets_edges);
2277  }
2278 
2279  MoFEMErrorCode removeBadEdges(const Range &surface, const Range &tets,
2280  const Range &fixed_edges,
2281  const Range &corner_nodes,
2282  Range &edges_to_merge,
2283  Range &not_merged_edges) {
2284  moab::Interface &moab(mField.get_moab());
2286 
2287  // find skin
2288  Skinner skin(&moab);
2289  Range tets_skin;
2290  CHKERR skin.find_skin(0, tets, false, tets_skin);
2291  Range surface_skin;
2292  CHKERR skin.find_skin(0, surface, false, surface_skin);
2293 
2294  // end nodes
2295  Range tets_skin_edges;
2296  CHKERR moab.get_adjacencies(tets_skin, 1, false, tets_skin_edges,
2297  moab::Interface::UNION);
2298  Range surface_front;
2299  surface_front = subtract(surface_skin, tets_skin_edges);
2300  Range surface_front_nodes;
2301  CHKERR moab.get_connectivity(surface_front, surface_front_nodes, true);
2302  Range ends_nodes;
2303  CHKERR skin.find_skin(0, surface_front, false, ends_nodes);
2304 
2305  // remove bad merges
2306 
2307  // get surface and body skin verts
2308  Range surface_edges;
2309  CHKERR moab.get_adjacencies(surface, 1, false, surface_edges,
2310  moab::Interface::UNION);
2311  // get nodes on the surface
2312  Range surface_edges_verts;
2313  CHKERR moab.get_connectivity(surface_edges, surface_edges_verts, true);
2314  // get vertices on the body skin
2315  Range tets_skin_edges_verts;
2316  CHKERR moab.get_connectivity(tets_skin_edges, tets_skin_edges_verts,
2317  true);
2318 
2319  Range edges_to_remove;
2320 
2321  // remove edges self connected to body skin
2322  {
2323  Range ents_nodes_and_edges;
2324  ents_nodes_and_edges.merge(tets_skin_edges_verts);
2325  ents_nodes_and_edges.merge(tets_skin_edges);
2326  CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2327  0, false);
2328  }
2329  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2330  not_merged_edges.merge(edges_to_remove);
2331 
2332  // remove edges self connected to surface
2333  {
2334  Range ents_nodes_and_edges;
2335  ents_nodes_and_edges.merge(surface_edges_verts);
2336  ents_nodes_and_edges.merge(surface_edges);
2337  ents_nodes_and_edges.merge(tets_skin_edges_verts);
2338  ents_nodes_and_edges.merge(tets_skin_edges);
2339  CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2340  0, false);
2341  }
2342  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2343  not_merged_edges.merge(edges_to_remove);
2344 
2345  // remove edges adjacent corner_nodes execpt those on fixed edges
2346  Range fixed_edges_nodes;
2347  CHKERR moab.get_connectivity(fixed_edges, fixed_edges_nodes, true);
2348  {
2349  Range ents_nodes_and_edges;
2350  ents_nodes_and_edges.merge(fixed_edges_nodes);
2351  ents_nodes_and_edges.merge(ends_nodes);
2352  ents_nodes_and_edges.merge(corner_nodes);
2353  ents_nodes_and_edges.merge(fixed_edges);
2354  CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2355  0, false);
2356  }
2357  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2358  not_merged_edges.merge(edges_to_remove);
2359 
2360  // remove edges self connected to surface
2361  CHKERR removeSelfConectingEdges(surface_edges, edges_to_remove, 0, false);
2362  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2363  not_merged_edges.merge(edges_to_remove);
2364 
2365  // remove edges self contented on surface skin
2366  {
2367  Range ents_nodes_and_edges;
2368  ents_nodes_and_edges.merge(surface_skin);
2369  ents_nodes_and_edges.merge(fixed_edges_nodes);
2370  CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2371  0, false);
2372  }
2373  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2374  not_merged_edges.merge(edges_to_remove);
2375 
2376  // remove crack front nodes connected to the surface, except those short
2377  {
2378  Range ents_nodes_and_edges;
2379  ents_nodes_and_edges.merge(surface_front_nodes);
2380  ents_nodes_and_edges.merge(surface_front);
2381  ents_nodes_and_edges.merge(tets_skin_edges_verts);
2382  ents_nodes_and_edges.merge(tets_skin_edges);
2383  CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2384  tOL, false);
2385  }
2386  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2387  not_merged_edges.merge(edges_to_remove);
2388 
2389  // remove edges connecting crack front and fixed edges, except those short
2390  {
2391  Range ents_nodes_and_edges;
2392  ents_nodes_and_edges.merge(surface_skin.subset_by_type(MBEDGE));
2393  ents_nodes_and_edges.merge(fixed_edges.subset_by_type(MBEDGE));
2394  CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2395  tOL, false);
2396  }
2397  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2398  not_merged_edges.merge(edges_to_remove);
2399 
2401  }
2402 
2403  private:
2404  MoFEMErrorCode removeSelfConectingEdges(const Range &ents,
2405  Range &edges_to_remove,
2406  const bool length,
2407  bool debug) const {
2408  moab::Interface &moab(mField.get_moab());
2410  // get nodes
2411  Range ents_nodes = ents.subset_by_type(MBVERTEX);
2412  if (ents_nodes.empty()) {
2413  CHKERR moab.get_connectivity(ents, ents_nodes, true);
2414  }
2415  // edges adj. to nodes
2416  Range ents_nodes_edges;
2417  CHKERR moab.get_adjacencies(ents_nodes, 1, false, ents_nodes_edges,
2418  moab::Interface::UNION);
2419  // nodes of adj. edges
2420  Range ents_nodes_edges_nodes;
2421  CHKERR moab.get_connectivity(ents_nodes_edges, ents_nodes_edges_nodes,
2422  true);
2423  // hanging nodes
2424  ents_nodes_edges_nodes = subtract(ents_nodes_edges_nodes, ents_nodes);
2425  Range ents_nodes_edges_nodes_edges;
2426  CHKERR moab.get_adjacencies(ents_nodes_edges_nodes, 1, false,
2427  ents_nodes_edges_nodes_edges,
2428  moab::Interface::UNION);
2429  // remove edges adj. to hanging edges
2430  ents_nodes_edges =
2431  subtract(ents_nodes_edges, ents_nodes_edges_nodes_edges);
2432  ents_nodes_edges =
2433  subtract(ents_nodes_edges, ents.subset_by_type(MBEDGE));
2434  if (length > 0) {
2435  Range::iterator eit = ents_nodes_edges.begin();
2436  for (; eit != ents_nodes_edges.end();) {
2437 
2438  int num_nodes;
2439  const EntityHandle *conn;
2440  rval = moab.get_connectivity(*eit, conn, num_nodes, true);
2441  double coords[6];
2442  if (tH)
2443  CHKERR moab.tag_get_data(tH, conn, num_nodes, coords);
2444  else
2445  CHKERR moab.get_coords(conn, num_nodes, coords);
2446 
2447  auto get_edge_length = [coords]() {
2449  &coords[0], &coords[1], &coords[2]);
2452  t_delta(i) = t_coords(i);
2453  ++t_coords;
2454  t_delta(i) -= t_coords(i);
2455  return sqrt(t_delta(i) * t_delta(i));
2456  };
2457 
2458  if (get_edge_length() < tOL) {
2459  eit = ents_nodes_edges.erase(eit);
2460  } else {
2461  eit++;
2462  }
2463  }
2464  }
2465  edges_to_remove.swap(ents_nodes_edges);
2466  if (debug) {
2467  CHKERR SaveData(moab)("edges_to_remove.vtk", edges_to_remove);
2468  }
2470  }
2471  };
2472 
2473  Range not_merged_edges;
2474  const double tol = 0.1;
2475  CHKERR Toplogy(m_field, th, tol * aveLength)
2476  .edgesToMerge(surface, tets, edges_to_merge);
2477  CHKERR Toplogy(m_field, th, tol * aveLength)
2478  .removeBadEdges(surface, tets, fixed_edges, corner_nodes, edges_to_merge,
2479  not_merged_edges);
2480  Toplogy::SetsMap sets_map;
2481  CHKERR Toplogy(m_field, th, tol * aveLength)
2482  .classifyVerts(surface, tets, fixed_edges, corner_nodes, sets_map);
2483  Range proc_tets;
2484  CHKERR Toplogy(m_field, th, tol * aveLength)
2485  .getProcTets(tets, edges_to_merge, proc_tets);
2486  out_tets = subtract(tets, proc_tets);
2487 
2488  if (bit_ptr) {
2489  Range all_out_ents = out_tets;
2490  for (int dd = 2; dd >= 0; dd--) {
2491  CHKERR moab.get_adjacencies(out_tets, dd, false, all_out_ents,
2492  moab::Interface::UNION);
2493  }
2494  CHKERR m_field.getInterface<BitRefManager>()->addBitRefLevel(all_out_ents,
2495  *bit_ptr);
2496  }
2497 
2498  int nb_nodes_merged = 0;
2499  LengthMapData_multi_index length_map;
2500  new_surf = surface;
2501 
2502  double ave0 = 0, ave = 0, min = 0, min_p = 0, min_pp;
2503  for (int pp = 0; pp != nbMaxMergingCycles; ++pp) {
2504 
2505  int nb_nodes_merged_p = nb_nodes_merged;
2506  length_map.clear();
2507  min_pp = min_p;
2508  min_p = min;
2509  CHKERR LengthMap(m_field, th, aveLength)(proc_tets, edges_to_merge,
2510  length_map, ave);
2511  min = length_map.get<2>().begin()->qUality;
2512  if (pp == 0) {
2513  ave0 = ave;
2514  }
2515 
2516  int nn = 0;
2517  Range collapsed_edges;
2518  for (auto mit = length_map.get<0>().begin();
2519  mit != length_map.get<0>().end(); mit++, nn++) {
2520 
2521  if (!mit->skip) {
2522 
2523  auto get_father_and_mother =
2524  [&](EntityHandle &father, EntityHandle &mother, int &line_search) {
2526  int num_nodes;
2527  const EntityHandle *conn;
2528  CHKERR moab.get_connectivity(mit->eDge, conn, num_nodes, true);
2529  int conn_type[2] = {0, 0};
2530  for (int nn = 0; nn != 2; nn++) {
2531  conn_type[nn] = 0;
2532  for (Toplogy::SetsMap::reverse_iterator sit = sets_map.rbegin();
2533  sit != sets_map.rend(); sit++) {
2534  if (sit->second.find(conn[nn]) != sit->second.end()) {
2535  conn_type[nn] |= sit->first;
2536  }
2537  }
2538  }
2539  int type_father, type_mother;
2540  if (conn_type[0] > conn_type[1]) {
2541  father = conn[0];
2542  mother = conn[1];
2543  type_father = conn_type[0];
2544  type_mother = conn_type[1];
2545  } else {
2546  father = conn[1];
2547  mother = conn[0];
2548  type_father = conn_type[1];
2549  type_mother = conn_type[0];
2550  }
2551  if (type_father == type_mother) {
2552  line_search = lineSearchSteps;
2553  }
2555  };
2556 
2557  int line_search = 0;
2558  EntityHandle father, mother;
2559  CHKERR get_father_and_mother(father, mother, line_search);
2560  CHKERR MergeNodes(m_field, true, line_search, th,
2561  update_meshsets)(father, mother, proc_tets, new_surf,
2562  edges_to_merge, not_merged_edges);
2563 
2564  if (m_field.getInterface<NodeMergerInterface>()->getSucessMerge()) {
2565  Range adj_mother_tets;
2566  CHKERR moab.get_adjacencies(&mother, 1, 3, false, adj_mother_tets);
2567  Range adj_mother_tets_nodes;
2568  CHKERR moab.get_connectivity(adj_mother_tets, adj_mother_tets_nodes,
2569  true);
2570  Range adj_edges;
2571  CHKERR moab.get_adjacencies(adj_mother_tets_nodes, 1, false,
2572  adj_edges, moab::Interface::UNION);
2573  CHKERR moab.get_adjacencies(&father, 1, 1, false, adj_edges,
2574  moab::Interface::UNION);
2575  Range proc_edges;
2576  CHKERR moab.get_adjacencies(proc_tets, 1, true, proc_edges);
2577  adj_edges = intersect(proc_edges, adj_edges);
2578  for (Range::iterator ait = adj_edges.begin(); ait != adj_edges.end();
2579  ++ait) {
2580  auto miit = length_map.get<1>().find(*ait);
2581  if (miit != length_map.get<1>().end())
2582  (const_cast<LengthMapData &>(*miit)).skip = true;
2583  }
2584  nb_nodes_merged++;
2585  collapsed_edges.insert(mit->eDge);
2586  }
2587 
2588  if (nn > static_cast<int>(length_map.size() / fraction_level))
2589  break;
2590  if (mit->qUality > ave)
2591  break;
2592  }
2593  }
2594 
2595  Range adj_faces, adj_edges;
2596  CHKERR moab.get_adjacencies(proc_tets, 2, false, adj_faces,
2597  moab::Interface::UNION);
2598  new_surf = intersect(new_surf, adj_faces);
2599 
2600  PetscPrintf(m_field.get_comm(),
2601  "(%d) Number of nodes merged %d ave q %3.4e min q %3.4e\n", pp,
2602  nb_nodes_merged, ave, min);
2603 
2604  if (debug) {
2605  // current surface and merged edges in step
2606  std::string name;
2607  name = "node_merge_step_" + boost::lexical_cast<std::string>(pp) + ".vtk";
2608  CHKERR SaveData(moab)(name, unite(new_surf, collapsed_edges));
2609  name = "edges_to_merge_step_" + boost::lexical_cast<std::string>(pp) +
2610  ".vtk";
2611  CHKERR SaveData(moab)(name, unite(new_surf, edges_to_merge));
2612  }
2613 
2614  if (nb_nodes_merged == nb_nodes_merged_p)
2615  break;
2616  if (min > 1e-2 && min == min_pp)
2617  break;
2618  if (min > ave0)
2619  break;
2620 
2621  CHKERR moab.get_adjacencies(proc_tets, 1, false, adj_edges,
2622  moab::Interface::UNION);
2623  edges_to_merge = intersect(edges_to_merge, adj_edges);
2624  CHKERR Toplogy(m_field, th, tol * aveLength)
2625  .removeBadEdges(new_surf, proc_tets, fixed_edges, corner_nodes,
2626  edges_to_merge, not_merged_edges);
2627  }
2628 
2629  if (bit_ptr) {
2630  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevel(proc_tets,
2631  *bit_ptr);
2632  }
2633  out_tets.merge(proc_tets);
2634 
2636 }
2637 
2639  const int fraction_level, const BitRefLevel cut_bit,
2640  const BitRefLevel trim_bit, const BitRefLevel bit, const Range &surface,
2641  const Range &fixed_edges, const Range &corner_nodes, Tag th,
2642  const bool update_meshsets, const bool debug) {
2643  CoreInterface &m_field = cOre;
2645  Range tets_level;
2646  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2647  trim_bit, BitRefLevel().set(), MBTET, tets_level);
2648 
2649  Range edges_to_merge;
2650  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByParentType(
2651  trim_bit, cut_bit | trim_bit, MBEDGE, edges_to_merge);
2652 
2653  // get all entities not in database
2654  Range all_ents_not_in_database_before;
2655  CHKERR cOre.getInterface<BitRefManager>()->getAllEntitiesNotInDatabase(
2656  all_ents_not_in_database_before);
2657 
2658  Range out_new_tets, out_new_surf;
2659  CHKERR mergeBadEdges(fraction_level, tets_level, surface, fixed_edges,
2660  corner_nodes, edges_to_merge, out_new_tets, out_new_surf,
2661  th, update_meshsets, &bit, debug);
2662 
2663  // get all entities not in database after merge
2664  Range all_ents_not_in_database_after;
2665  CHKERR cOre.getInterface<BitRefManager>()->getAllEntitiesNotInDatabase(
2666  all_ents_not_in_database_after);
2667 
2668  // delete hanging entities
2669  all_ents_not_in_database_after =
2670  subtract(all_ents_not_in_database_after, all_ents_not_in_database_before);
2671  Range meshsets;
2672  CHKERR m_field.get_moab().get_entities_by_type(0, MBENTITYSET, meshsets,
2673  false);
2674  for (auto m : meshsets)
2675  CHKERR m_field.get_moab().remove_entities(m,
2676  all_ents_not_in_database_after);
2677 
2678  m_field.get_moab().delete_entities(all_ents_not_in_database_after);
2679 
2680  mergedVolumes.swap(out_new_tets);
2681  mergedSurfaces.swap(out_new_surf);
2683 }
2684 
2685 #ifdef WITH_TETGEN
2686 
2688  vector<string> &switches, const BitRefLevel &mesh_bit,
2689  const BitRefLevel &bit, const Range &surface, const Range &fixed_edges,
2690  const Range &corner_nodes, Tag th, const bool debug) {
2691  CoreInterface &m_field = cOre;
2692  moab::Interface &moab = m_field.get_moab();
2693  TetGenInterface *tetgen_iface;
2695  CHKERR m_field.getInterface(tetgen_iface);
2696 
2697  tetGenData.clear();
2698  moabTetGenMap.clear();
2699  tetGenMoabMap.clear();
2700 
2701  if (tetGenData.size() < 1) {
2702  tetGenData.push_back(new tetgenio);
2703  }
2704  tetgenio &in = tetGenData.back();
2705 
2706  struct BitEnts {
2707 
2708  CoreInterface &mField;
2709  const BitRefLevel &bIt;
2710  BitEnts(CoreInterface &m_field, const BitRefLevel &bit)
2711  : mField(m_field), bIt(bit) {}
2712 
2713  Range mTets;
2714  Range mTris;
2715  Range mEdges;
2716  Range mNodes;
2717 
2718  MoFEMErrorCode getLevelEnts() {
2720  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2721  bIt, BitRefLevel().set(), MBTET, mTets);
2722  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2723  bIt, BitRefLevel().set(), MBTRI, mTris);
2724  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2725  bIt, BitRefLevel().set(), MBEDGE, mEdges);
2726  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2727  bIt, BitRefLevel().set(), MBVERTEX, mNodes);
2729  }
2730 
2731  Range mSkin;
2732  Range mSkinNodes;
2733  Range mSkinEdges;
2734 
2735  MoFEMErrorCode getSkin() {
2736  moab::Interface &moab = mField.get_moab();
2738  Skinner skin(&moab);
2739  CHKERR skin.find_skin(0, mTets, false, mSkin);
2740  CHKERR mField.get_moab().get_connectivity(mSkin, mSkinNodes, true);
2741  CHKERR mField.get_moab().get_adjacencies(mSkin, 1, false, mSkinEdges,
2742  moab::Interface::UNION);
2744  }
2745  };
2746 
2747  struct SurfaceEnts {
2748 
2749  CoreInterface &mField;
2750  SurfaceEnts(CoreInterface &m_field) : mField(m_field) {}
2751 
2752  Range sNodes;
2753  Range sEdges;
2754  Range sVols;
2755  Range vNodes;
2756 
2757  MoFEMErrorCode getVolume(const BitEnts &bit_ents, const Range &tris) {
2758  moab::Interface &moab = mField.get_moab();
2760  CHKERR moab.get_connectivity(tris, sNodes, true);
2761  CHKERR moab.get_adjacencies(tris, 1, false, sEdges,
2762  moab::Interface::UNION);
2763  CHKERR moab.get_adjacencies(sNodes, 3, false, sVols,
2764  moab::Interface::UNION);
2765  sVols = intersect(sVols, bit_ents.mTets);
2766  CHKERR moab.get_connectivity(sVols, vNodes, true);
2768  }
2769 
2770  Range sSkin;
2771  Range sSkinNodes;
2772  Range vSkin;
2773  Range vSkinNodes;
2774  Range vSkinWithoutBodySkin;
2775  Range vSkinNodesWithoutBodySkin;
2776  Range vSkinOnBodySkin;
2777  Range vSkinOnBodySkinNodes;
2778 
2779  MoFEMErrorCode getSkin(const BitEnts &bit_ents, const Range &tris,
2780  const int levels = 3) {
2781  moab::Interface &moab = mField.get_moab();
2783  Skinner skin(&moab);
2784  rval = skin.find_skin(0, sVols, false, vSkin);
2785  for (int ll = 0; ll != levels; ll++) {
2786  CHKERR moab.get_adjacencies(vSkin, 3, false, sVols,
2787  moab::Interface::UNION);
2788  sVols = intersect(sVols, bit_ents.mTets);
2789  vSkin.clear();
2790  CHKERR skin.find_skin(0, sVols, false, vSkin);
2791  }
2792  vSkinWithoutBodySkin = subtract(vSkin, bit_ents.mSkin);
2793  vSkinOnBodySkin = intersect(vSkin, bit_ents.mSkin);
2794  CHKERR moab.get_connectivity(vSkinOnBodySkin, vSkinOnBodySkinNodes, true);
2795  CHKERR moab.get_connectivity(sVols, vNodes, true);
2796  CHKERR moab.get_connectivity(vSkin, vSkinNodes, true);
2797  vSkinNodesWithoutBodySkin = subtract(vSkinNodes, bit_ents.mSkinNodes);
2798  CHKERR skin.find_skin(0, tris, false, sSkin);
2799  CHKERR moab.get_connectivity(sSkin, sSkinNodes, true);
2800  tVols = sVols;
2802  }
2803 
2804  Range tVols;
2805 
2806  MoFEMErrorCode getTetsForRemesh(const BitEnts &bit_ents, Tag th = NULL) {
2807  moab::Interface &moab = mField.get_moab();
2809 
2810  Range tets_with_four_nodes_on_skin;
2811  CHKERR moab.get_adjacencies(vSkinOnBodySkinNodes, 3, false,
2812  tets_with_four_nodes_on_skin,
2813  moab::Interface::UNION);
2814  Range tets_nodes;
2815  CHKERR moab.get_connectivity(tets_with_four_nodes_on_skin, tets_nodes,
2816  true);
2817  tets_nodes = subtract(tets_nodes, vSkinOnBodySkinNodes);
2818  Range other_tets;
2819  CHKERR moab.get_adjacencies(tets_nodes, 3, false, other_tets,
2820  moab::Interface::UNION);
2821  tets_with_four_nodes_on_skin =
2822  subtract(tets_with_four_nodes_on_skin, other_tets);
2823  Range to_remove;
2824  for (Range::iterator tit = tets_with_four_nodes_on_skin.begin();
2825  tit != tets_with_four_nodes_on_skin.end(); tit++) {
2826  int num_nodes;
2827  const EntityHandle *conn;
2828  CHKERR moab.get_connectivity(*tit, conn, num_nodes, true);
2829  double coords[12];
2830  if (th) {
2831  CHKERR moab.tag_get_data(th, conn, num_nodes, coords);
2832  } else {
2833  CHKERR moab.get_coords(conn, num_nodes, coords);
2834  }
2835  double quality = Tools::volumeLengthQuality(coords);
2836  if (quality < 1e-2) {
2837  to_remove.insert(*tit);
2838  }
2839  }
2840 
2841  sVols = subtract(sVols, to_remove);
2842 
2843  Skinner skin(&moab);
2844  vSkin.clear();
2845  CHKERR skin.find_skin(0, sVols, false, vSkin);
2846  Range m_skin;
2847  CHKERR
2848  skin.find_skin(0, subtract(bit_ents.mSkin, to_remove), false, m_skin);
2849  vSkinWithoutBodySkin = subtract(vSkin, m_skin);
2850  vSkinOnBodySkin = intersect(vSkin, m_skin);
2851  vNodes.clear();
2852  vSkinNodes.clear();
2853  vSkinOnBodySkinNodes.clear();
2854  CHKERR moab.get_connectivity(sVols, vNodes, true);
2855  CHKERR moab.get_connectivity(vSkinOnBodySkin, vSkinOnBodySkinNodes, true);
2856  CHKERR moab.get_connectivity(vSkin, vSkinNodes, true);
2858  }
2859  };
2860 
2861  BitEnts bit_ents(m_field, mesh_bit);
2862  CHKERR bit_ents.getLevelEnts();
2863  CHKERR bit_ents.getSkin();
2864  SurfaceEnts surf_ents(m_field);
2865  CHKERR surf_ents.getVolume(bit_ents, surface);
2866  CHKERR surf_ents.getSkin(bit_ents, surface);
2867  CHKERR surf_ents.getTetsForRemesh(bit_ents);
2868 
2869  map<int, Range> types_ents;
2870  types_ents[TetGenInterface::RIDGEVERTEX].merge(
2871  surf_ents.vSkinNodesWithoutBodySkin);
2872  // FREESEGVERTEX
2873  types_ents[TetGenInterface::FREESEGVERTEX].merge(surf_ents.sSkinNodes);
2874  types_ents[TetGenInterface::FREESEGVERTEX] =
2875  subtract(types_ents[TetGenInterface::FREESEGVERTEX],
2876  types_ents[TetGenInterface::RIDGEVERTEX]);
2877  // FREEFACETVERTEX
2878  types_ents[TetGenInterface::FREEFACETVERTEX].merge(surf_ents.sNodes);
2879  types_ents[TetGenInterface::FREEFACETVERTEX] =
2880  subtract(types_ents[TetGenInterface::FREEFACETVERTEX],
2881  types_ents[TetGenInterface::RIDGEVERTEX]);
2882  types_ents[TetGenInterface::FREEFACETVERTEX] =
2883  subtract(types_ents[TetGenInterface::FREEFACETVERTEX],
2884  types_ents[TetGenInterface::FREESEGVERTEX]);
2885  // FREEVOLVERTEX
2886  types_ents[TetGenInterface::FREEVOLVERTEX].merge(surf_ents.vNodes);
2887  types_ents[TetGenInterface::FREEVOLVERTEX] =
2888  subtract(types_ents[TetGenInterface::FREEVOLVERTEX],
2889  types_ents[TetGenInterface::RIDGEVERTEX]);
2890  types_ents[TetGenInterface::FREEVOLVERTEX] =
2891  subtract(types_ents[TetGenInterface::FREEVOLVERTEX],
2892  types_ents[TetGenInterface::FREESEGVERTEX]);
2893  types_ents[TetGenInterface::FREEVOLVERTEX] =
2894  subtract(types_ents[TetGenInterface::FREEVOLVERTEX],
2895  types_ents[TetGenInterface::FREEFACETVERTEX]);
2896 
2897  Tag th_marker;
2898  // Clean markers
2899  rval = m_field.get_moab().tag_get_handle("TETGEN_MARKER", th_marker);
2900  if (rval == MB_SUCCESS) {
2901  CHKERR m_field.get_moab().tag_delete(th_marker);
2902  rval = MB_SUCCESS;
2903  }
2904 
2905  int def_marker = 0;
2906  CHKERR m_field.get_moab().tag_get_handle(
2907  "TETGEN_MARKER", 1, MB_TYPE_INTEGER, th_marker,
2908  MB_TAG_CREAT | MB_TAG_SPARSE, &def_marker);
2909 
2910  // Mark surface with id = 1
2911  vector<int> markers(surface.size(), 1);
2912  CHKERR m_field.get_moab().tag_set_data(th_marker, surface, &*markers.begin());
2913  // Mark all side sets
2914  int shift = 1;
2915  map<int, int> id_shift_map; // each meshset has set unique bit
2917  (*cOre.getInterface<MeshsetsManager>()), SIDESET, it)) {
2918  int ms_id = it->getMeshsetId();
2919  id_shift_map[ms_id] = 1 << shift; // shift bit
2920  ++shift;
2921  Range sideset_faces;
2922  CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
2923  ms_id, SIDESET, 2, sideset_faces, true);
2924  markers.resize(sideset_faces.size());
2925  CHKERR m_field.get_moab().tag_get_data(th_marker, sideset_faces,
2926  &*markers.begin());
2927  for (unsigned int ii = 0; ii < markers.size(); ii++) {
2928  markers[ii] |= id_shift_map[ms_id]; // add bit to marker
2929  }
2930  CHKERR m_field.get_moab().tag_set_data(th_marker, sideset_faces,
2931  &*markers.begin());
2932  }
2933  Range nodes_to_remove; // none
2934  markers.resize(nodes_to_remove.size());
2935  fill(markers.begin(), markers.end(), -1);
2936  CHKERR m_field.get_moab().tag_set_data(th_marker, nodes_to_remove,
2937  &*markers.begin());
2938 
2939  // nodes
2940  if (tetGenData.size() == 1) {
2941 
2942  Range ents_to_tetgen = surf_ents.sVols;
2943  CHKERR m_field.get_moab().get_connectivity(surf_ents.sVols, ents_to_tetgen,
2944  true);
2945  for (int dd = 2; dd >= 1; dd--) {
2946  CHKERR m_field.get_moab().get_adjacencies(
2947  surf_ents.sVols, dd, false, ents_to_tetgen, moab::Interface::UNION);
2948  }
2949 
2950  // Load mesh to TetGen data structures
2951  CHKERR tetgen_iface->inData(ents_to_tetgen, in, moabTetGenMap,
2952  tetGenMoabMap, th);
2953  CHKERR tetgen_iface->setGeomData(in, moabTetGenMap, tetGenMoabMap,
2954  types_ents);
2955  std::vector<pair<Range, int>> markers;
2956  for (Range::iterator tit = surface.begin(); tit != surface.end(); tit++) {
2957  Range facet;
2958  facet.insert(*tit);
2959  markers.push_back(pair<Range, int>(facet, 2));
2960  }
2961  for (Range::iterator tit = surf_ents.vSkinWithoutBodySkin.begin();
2962  tit != surf_ents.vSkinWithoutBodySkin.end(); tit++) {
2963  Range facet;
2964  facet.insert(*tit);
2965  markers.push_back(pair<Range, int>(facet, 1));
2966  }
2967  Range other_facets;
2968  other_facets = subtract(surf_ents.vSkin, surf_ents.vSkinWithoutBodySkin);
2969  for (Range::iterator tit = other_facets.begin(); tit != other_facets.end();
2970  tit++) {
2971  Range facet;
2972  facet.insert(*tit);
2973  markers.push_back(pair<Range, int>(facet, 0));
2974  }
2975  CHKERR tetgen_iface->setFaceData(markers, in, moabTetGenMap, tetGenMoabMap);
2976  }
2977 
2978  if (debug) {
2979  if (m_field.get_comm_rank() == 0) {
2980  char tetgen_in_file_name[] = "in";
2981  in.save_nodes(tetgen_in_file_name);
2982  in.save_elements(tetgen_in_file_name);
2983  in.save_faces(tetgen_in_file_name);
2984  in.save_edges(tetgen_in_file_name);
2985  in.save_poly(tetgen_in_file_name);
2986  }
2987  }
2988 
2989  // generate new mesh
2990  {
2991  vector<string>::iterator sw = switches.begin();
2992  for (int ii = 0; sw != switches.end(); sw++, ii++) {
2993  tetgenio &_in_ = tetGenData.back();
2994  tetGenData.push_back(new tetgenio);
2995  tetgenio &_out_ = tetGenData.back();
2996  char *s = const_cast<char *>(sw->c_str());
2997  CHKERR tetgen_iface->tetRahedralize(s, _in_, _out_);
2998  }
2999  }
3000  tetgenio &out = tetGenData.back();
3001  // save elems
3002  if (debug) {
3003  char tetgen_out_file_name[] = "out";
3004  out.save_nodes(tetgen_out_file_name);
3005  out.save_elements(tetgen_out_file_name);
3006  out.save_faces(tetgen_out_file_name);
3007  out.save_edges(tetgen_out_file_name);
3008  out.save_poly(tetgen_out_file_name);
3009  }
3010 
3011  CHKERR tetgen_iface->outData(in, out, moabTetGenMap, tetGenMoabMap, bit,
3012  false, false, true, th);
3013 
3014  Range rest_of_ents = subtract(bit_ents.mTets, surf_ents.tVols);
3015  for (int dd = 2; dd >= 0; dd--) {
3016  CHKERR moab.get_adjacencies(rest_of_ents.subset_by_dimension(3), dd, false,
3017  rest_of_ents, moab::Interface::UNION);
3018  }
3019  CHKERR m_field.getInterface<BitRefManager>()->addBitRefLevel(rest_of_ents,
3020  bit);
3021 
3022  Range tetgen_faces;
3023  map<int, Range> face_markers_map;
3024  CHKERR tetgen_iface->getTriangleMarkers(tetGenMoabMap, out, &tetgen_faces,
3025  &face_markers_map);
3026 
3027  tetgenSurfaces = face_markers_map[1];
3028  for (map<int, Range>::iterator mit = face_markers_map.begin();
3029  mit != face_markers_map.end(); mit++) {
3031  (*cOre.getInterface<MeshsetsManager>()), SIDESET, it)) {
3032  int msId = it->getMeshsetId();
3033  if (id_shift_map[msId] & mit->first) {
3034  EntityHandle meshset = it->getMeshset();
3035  CHKERR m_field.get_moab().add_entities(
3036  meshset, mit->second.subset_by_type(MBTRI));
3037  }
3038  }
3039  }
3040 
3041  auto get_min_quality = [&m_field, debug](const BitRefLevel bit, Tag th) {
3042  Range tets_level; // test at level
3043  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
3044  bit, BitRefLevel().set(), MBTET, tets_level);
3045  double min_q = 1;
3046  CHKERR m_field.getInterface<Tools>()->minTetsQuality(tets_level, min_q, th);
3047  return min_q;
3048  };
3049 
3050  PetscPrintf(PETSC_COMM_WORLD, "Min quality tetgen %6.4g\n",
3051  get_min_quality(bit, th));
3052 
3054 }
3055 
3056 #endif // WITH_TETGEN
3057 
3059  CoreInterface &m_field = cOre;
3060  moab::Interface &moab = m_field.get_moab();
3062  Range nodes;
3063  if (bit.none()) {
3064  CHKERR moab.get_entities_by_type(0, MBVERTEX, nodes);
3065  } else {
3066  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
3067  bit, BitRefLevel().set(), MBVERTEX, nodes);
3068  }
3069  std::vector<double> coords(3 * nodes.size());
3070  CHKERR moab.get_coords(nodes, &coords[0]);
3071  CHKERR moab.tag_set_data(th, nodes, &coords[0]);
3073 }
3074 
3076  const BitRefLevel mask) {
3077  CoreInterface &m_field = cOre;
3078  moab::Interface &moab = m_field.get_moab();
3080  Range nodes;
3081  if (bit.none()) {
3082  CHKERR moab.get_entities_by_type(0, MBVERTEX, nodes);
3083  } else {
3084  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
3085  bit, mask, MBVERTEX, nodes);
3086  }
3087  std::vector<double> coords(3 * nodes.size());
3088  CHKERR moab.tag_get_data(th, nodes, &coords[0]);
3089  CHKERR moab.set_coords(nodes, &coords[0]);
3091 }
3092 
3094  CoreInterface &m_field = cOre;
3095  moab::Interface &moab = m_field.get_moab();
3097  CHKERR SaveData(moab)(prefix + "out_vol.vtk", vOlume);
3098  CHKERR SaveData(moab)(prefix + "out_surface.vtk", sUrface);
3099  CHKERR SaveData(moab)(prefix + "out_cut_edges.vtk", cutEdges);
3100  CHKERR SaveData(moab)(prefix + "out_cut_volumes.vtk", cutVolumes);
3101  CHKERR SaveData(moab)(prefix + "out_cut_new_volumes.vtk", cutNewVolumes);
3102  CHKERR SaveData(moab)(prefix + "out_cut_new_surfaces.vtk", cutNewSurfaces);
3103  CHKERR SaveData(moab)(prefix + "out_cut_zero_distance_ents.vtk",
3106 }
3107 
3109  moab::Interface &moab = cOre.getInterface<CoreInterface>()->get_moab();
3111  CHKERR SaveData(moab)("out_trim_surface.vtk", sUrface);
3112  CHKERR SaveData(moab)("out_trim_new_volumes.vtk", trimNewVolumes);
3113  CHKERR SaveData(moab)("out_trim_new_surfaces.vtk", trimNewSurfaces);
3114  CHKERR SaveData(moab)("out_trim_edges.vtk", trimEdges);
3116 }
3117 
3118 } // namespace MoFEM
MoFEMErrorCode rebuildMeshWithTetGen(vector< string > &switches, const BitRefLevel &mesh_bit, const BitRefLevel &bit, const Range &surface, const Range &fixed_edges, const Range &corner_nodes, Tag th=NULL, const bool debug=false)
MoFEMErrorCode buildTree()
build tree
virtual MoFEMErrorCode refine_TET(const EntityHandle meshset, const BitRefLevel &bit, const bool respect_interface=false, int verb=QUIET, Range *ref_edges=NULL)
refine TET in the meshset
map< EntityHandle, TreeData > edgesToTrim
MoFEMErrorCode tetRahedralize(char switches[], tetgenio &in, tetgenio &out)
run tetgen
virtual MoFEMErrorCode get_ref_ents(const RefEntity_multiIndex **refined_ents_ptr) const =0
Get ref entities multi-index from database.
MoFEMErrorCode cutAndTrim(int &first_bit, const int before_trim_levels, Tag th, const double tol_cut, const double tol_cut_close, const double tol_trim, const double tol_trim_close, Range *fixed_edges=NULL, Range *corner_nodes=NULL, const bool update_meshsets=false, const bool debug=false)
MoFEM interface unique ID.
Interface to cut meshes.
MoFEMErrorCode outData(tetgenio &in, tetgenio &out, moabTetGen_Map &moab_tetgen_map, tetGenMoab_Map &tetgen_moab_map, Range *ents=NULL, bool id_in_tags=false, bool error_if_created=false, bool assume_first_nodes_the_same=false, Tag th=nullptr)
get entities for TetGen data structure
MoFEMErrorCode setTagData(Tag th, const BitRefLevel bit=BitRefLevel())
set coords to tag
MoFEMErrorCode getTriangleMarkers(tetGenMoab_Map &tetgen_moab_map, tetgenio &out, Range *ents=NULL, idxRange_Map *ents_map=NULL, bool only_non_zero=true)
get markers to faces
LengthMapData(const double l, double q, const EntityHandle e)
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:32
virtual moab::Interface & get_moab()=0
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
SaveData(moab::Interface &moab)
double maxLength
Maximal edge length.
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
use TetGen to generate mesh
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)
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:499
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:284
map< EntityHandle, TreeData > verticesOnTrimEdges
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Types.hpp:77
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
base class for all interface classes
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:475
MoFEMErrorCode setFront(const Range &surface)
map< EntityHandle, unsigned long > moabTetGenMap
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 ...
double aveLength
Average edge length.
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:506
MoFEMErrorCode removePathologicalFrontTris(const BitRefLevel split_bit, Range &ents)
Remove pathological elements on surface internal front.
MoFEMErrorCode findEdgesToTrim(Range *fixed_edges, Range *corner_nodes, Tag th=NULL, const double tol=1e-4, const double tol_close=0, const bool debug=false)
Find edges to trimEdges.
MoFEMErrorCode refineBeforeCut(const BitRefLevel &bit, Range *fixed_edges, const bool update_meshsets, const bool debug=false)
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
MoFEMErrorCode setFaceData(std::vector< std::pair< Range, int > > &markers, tetgenio &in, moabTetGen_Map &moab_tetgen_map, tetGenMoab_Map &tetgen_moab_map)
set markers to faces
VectorBoundedArray< double, 6 > VectorDouble6
Definition: Types.hpp:86
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Exceptions.hpp:84
auto getVectorAdaptor
Get Vector adaptor.
Definition: Types.hpp:121
CutMeshInterface(const MoFEM::Core &core)
moab::Interface & moab
const Range & projectZeroDistanceEnts() const
VectorBoundedArray< double, 12 > VectorDouble12
Definition: Types.hpp:88
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
MoFEMErrorCode saveCutEdges(const std::string prefix="")
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
virtual int get_comm_rank() const =0
const Range & getNewTrimSurfaces() const
MoFEMErrorCode findEdgesToCut(Range *fixed_edges, Range *corner_nodes, const double low_tol, int verb=QUIET, const bool debug=false)
find edges to cut
double tol
MatrixShallowArrayAdaptor< double > MatrixAdaptor
Matrix adaptor.
Definition: Types.hpp:141
MoFEMErrorCode operator()(Core &core, const BitRefLevel &bit) const
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
const Range & getVolume() const
boost::ptr_vector< tetgenio > tetGenData
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
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:99
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 inData(Range &ents, tetgenio &in, moabTetGen_Map &moab_tetgen_map, tetGenMoab_Map &tetgen_moab_map, Tag th=NULL)
create TetGen data structure form range of moab entities
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 cutEdgesInMiddle(const BitRefLevel bit, const bool debug=false)
cut edges
MoFEMErrorCode updateRange(const Range &parent, Range &child)
Update range by prents.
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field...
#define CHKERR
Inline error check.
Definition: definitions.h:594
static const MOFEMuuid IDD_MOFEMCutMesh
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
#define _IT_CUBITMESHSETS_FOR_LOOP_(MESHSET_MANAGER, IT)
Iterator that loops over all the Cubit MeshSets in a moFEM field.
MoFEMErrorCode operator()(const std::string name, const Range &ents)
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:405
virtual MPI_Comm & get_comm() const =0
MoFEMErrorCode saveTrimEdges()
VectorBoundedArray< double, 9 > VectorDouble9
Definition: Types.hpp:87
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
map< unsigned long, EntityHandle > tetGenMoabMap
MoFEMErrorCode trimEdgesInTheMiddle(const BitRefLevel bit, Tag th=NULL, const double tol=1e-4, const bool debug=false)
trim edges
MoFEMErrorCode refineBeforeTrim(const BitRefLevel &bit, Range *fixed_edges, const bool update_meshsets, const bool debug=false)
MoFEMErrorCode cutTrimAndMerge(int &first_bit, const int fraction_level, const int before_trim_levels, Tag th, const double tol_cut, const double tol_cut_close, const double tol_trim, const double tol_trim_close, Range &fixed_edges, Range &corner_nodes, const bool update_meshsets=false, const bool debug=false)
map< EntityHandle, TreeData > edgesToCut
MoFEMErrorCode setGeomData(tetgenio &in, moabTetGen_Map &moab_tetgen_map, tetGenMoab_Map &tetgen_moab_map, std::map< int, Range > &type_ents)
set point tags and type
multi_index_container< ParentChild, indexed_by< hashed_unique< member< ParentChild, EntityHandle, &ParentChild::pArent > >, hashed_non_unique< member< ParentChild, EntityHandle, &ParentChild::cHild > > > > ParentChildMap
Definition: NodeMerger.hpp:128
MoFEMErrorCode setVolume(const Range &volume)
set volume entities
virtual MoFEMErrorCode add_verices_in_the_middel_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 ...
bool getSucessMerge()
Return true if successful merge.
Definition: NodeMerger.hpp:44
MoFEMErrorCode moveMidNodesOnTrimmedEdges(Tag th=NULL)
move trimmed edges mid nodes
MoFEMErrorCode refCutTrimAndMerge(int &first_bit, const int fraction_level, const int ref_before_cut_levels, const int befor_trim_levels, Tag th, const double tol_cut, const double tol_cut_close, const double tol_trim, const double tol_trim_close, Range &fixed_edges, Range &corner_nodes, const bool update_meshsets=false, const bool debug=false)