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