v0.8.19
CutMeshInterface.cpp
Go to the documentation of this file.
1 /** \file CutMeshInterface.cpp
2  * \brief Cut mesh by surface
3  *
4  * MoFEM is free software: you can redistribute it and/or modify it under
5  * the terms of the GNU Lesser General Public License as published by the
6  * Free Software Foundation, either version 3 of the License, or (at your
7  * option) any later version.
8  *
9  * MoFEM is distributed in the hope that it will be useful, but WITHOUT
10  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
11  * FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
12  * License for more details.
13  *
14  * You should have received a copy of the GNU Lesser General Public
15  * License along with MoFEM. If not, see <http://www.gnu.org/licenses/>
16  */
17 
18 namespace MoFEM {
19 
20 struct SaveData {
21  moab::Interface &moab;
22  SaveData(moab::Interface &moab) : moab(moab) {}
23  MoFEMErrorCode operator()(const std::string name, const Range &ents) {
25  EntityHandle meshset;
26  CHKERR moab.create_meshset(MESHSET_SET, meshset);
27  CHKERR moab.add_entities(meshset, ents);
28  CHKERR moab.write_file(name.c_str(), "VTK", "", &meshset, 1);
29  CHKERR moab.delete_entities(&meshset, 1);
31  }
32 };
33 
36  UnknownInterface **iface) const {
38  *iface = NULL;
39  if (uuid == IDD_MOFEMCutMesh) {
40  *iface = const_cast<CutMeshInterface *>(this);
42  }
43  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "unknown interface");
45 }
46 
48  : cOre(const_cast<Core &>(core)) {
49  lineSearchSteps = 10;
50  nbMaxMergingCycles = 200;
52 }
53 
56  treeSurfPtr.reset();
58 }
59 
62  sUrface = surface;
64 }
65 
66 MoFEMErrorCode CutMeshInterface::copySurface(const Range &surface, Tag th,
67  double *shift, double *origin,
68  double *transform,
69  const std::string save_mesh) {
70  CoreInterface &m_field = cOre;
71  moab::Interface &moab = m_field.get_moab();
73  std::map<EntityHandle,EntityHandle> verts_map;
74  for (Range::const_iterator tit = surface.begin(); tit != surface.end();
75  tit++) {
76  int num_nodes;
77  const EntityHandle *conn;
78  CHKERR moab.get_connectivity(*tit, conn, num_nodes, true);
79  MatrixDouble coords(num_nodes, 3);
80  if (th) {
81  CHKERR moab.tag_get_data(th, conn, num_nodes, &coords(0, 0));
82  } else {
83  CHKERR moab.get_coords(conn, num_nodes, &coords(0, 0));
84  }
85  EntityHandle new_verts[num_nodes];
86  for (int nn = 0; nn != num_nodes; nn++) {
87  if(verts_map.find(conn[nn])!=verts_map.end()) {
88  new_verts[nn] = verts_map[conn[nn]];
89  } else {
90  if (transform) {
91  ublas::matrix_row<MatrixDouble> mr(coords, nn);
92  if (origin) {
93  VectorAdaptor vec_origin(
94  3, ublas::shallow_array_adaptor<double>(3, origin));
95  mr = mr - vec_origin;
96  }
97  MatrixAdaptor mat_transform = MatrixAdaptor(
98  3, 3, ublas::shallow_array_adaptor<double>(9, transform));
99  mr = prod(mat_transform, mr);
100  if (origin) {
101  VectorAdaptor vec_origin(
102  3, ublas::shallow_array_adaptor<double>(3, origin));
103  mr = mr + vec_origin;
104  }
105  }
106  if (shift) {
107  ublas::matrix_row<MatrixDouble> mr(coords, nn);
108  VectorAdaptor vec_shift(
109  3, ublas::shallow_array_adaptor<double>(3, shift));
110  mr = mr + vec_shift;
111  }
112  CHKERR moab.create_vertex(&coords(nn, 0), new_verts[nn]);
113  verts_map[conn[nn]] = new_verts[nn];
114  }
115  }
116  EntityHandle ele;
117  CHKERR moab.create_element(MBTRI, new_verts, num_nodes, ele);
118  sUrface.insert(ele);
119  }
120  if (!save_mesh.empty())
121  CHKERR SaveData(m_field.get_moab())(save_mesh, sUrface);
123 }
124 
127  vOlume = volume;
129 }
130 
133  sUrface.merge(surface);
135 }
136 
139  vOlume.merge(volume);
141 }
142 
144  CoreInterface &m_field = cOre;
145  moab::Interface &moab = m_field.get_moab();
147  treeSurfPtr = boost::shared_ptr<OrientedBoxTreeTool>(
148  new OrientedBoxTreeTool(&moab, "ROOTSETSURF", true));
151 }
152 
154  MoFEMErrorCode operator()(Core &core, const BitRefLevel &bit) const {
157  ->updateAllMeshsetsByEntitiesChildren(bit);
158  CHKERRG(ierr);
159 
161  }
162 };
163 
165  const BitRefLevel &bit_level1, const BitRefLevel &bit_level2, Tag th,
166  const double tol_cut, const double tol_cut_close, const double tol_trim,
167  const double tol_trim_close, Range *fixed_edges, Range *corner_nodes,
168  const bool update_meshsets,const bool debug) {
169  CoreInterface &m_field = cOre;
170  moab::Interface &moab = m_field.get_moab();
172 
173  // cut mesh
174  CHKERR findEdgesToCut(fixed_edges, corner_nodes, tol_cut, QUIET, debug);
175  CHKERR projectZeroDistanceEnts(fixed_edges, corner_nodes, tol_cut_close,
176  QUIET, debug);
177  CHKERR cutEdgesInMiddle(bit_level1, debug);
178  if (fixed_edges) {
179  CHKERR cOre.getInterface<BitRefManager>()->updateRange(*fixed_edges,
180  *fixed_edges);
181  }
182  if (corner_nodes) {
183  CHKERR cOre.getInterface<BitRefManager>()->updateRange(*corner_nodes,
184  *corner_nodes);
185  }
186  if (update_meshsets) {
187  CHKERR UpdateMeshsets()(cOre, bit_level1);
188  }
190 
191  auto get_min_quality = [&m_field](const BitRefLevel bit, Tag th) {
192  Range tets_level; // test at level
193  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
194  bit, BitRefLevel().set(), MBTET, tets_level);
195  double min_q = 1;
196  CHKERR m_field.getInterface<Tools>()->minTetsQuality(tets_level, min_q, th);
197  return min_q;
198  };
199 
200  PetscPrintf(PETSC_COMM_WORLD, "Min quality cut %6.4g\n",
201  get_min_quality(bit_level1, th));
202 
203  if(debug) {
205  }
206 
207  // trim mesh
208  CHKERR findEdgesToTrim(fixed_edges, corner_nodes, th, tol_trim, debug);
209  CHKERR trimEdgesInTheMiddle(bit_level2, th, tol_trim_close, debug);
210  if (fixed_edges) {
211  CHKERR cOre.getInterface<BitRefManager>()->updateRange(*fixed_edges,
212  *fixed_edges);
213  }
214  if (corner_nodes) {
215  CHKERR cOre.getInterface<BitRefManager>()->updateRange(*corner_nodes,
216  *corner_nodes);
217  }
218  if (update_meshsets) {
219  CHKERR UpdateMeshsets()(cOre, bit_level2);
220  }
222 
223  PetscPrintf(PETSC_COMM_WORLD, "Min quality trim %3.2g\n",
224  get_min_quality(bit_level2, th));
225 
226  if(debug) {
228  Range bit2_edges;
229  CHKERR cOre.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
230  bit_level2, BitRefLevel().set(), MBEDGE, bit2_edges);
231  CHKERR SaveData(moab)("trim_fixed_edges.vtk",
232  intersect(*fixed_edges, bit2_edges));
233  }
234 
236 }
237 
239  const int fraction_level, const BitRefLevel &bit_level1,
240  const BitRefLevel &bit_level2, const BitRefLevel &bit_level3, Tag th,
241  const double tol_cut, const double tol_cut_close, const double tol_trim,
242  const double tol_trim_close, Range &fixed_edges, Range &corner_nodes,
243  const bool update_meshsets, const bool debug) {
244  CoreInterface &m_field = cOre;
246  if(debug) {
247  CHKERR cOre.getInterface<BitRefManager>()->writeEntitiesNotInDatabase(
248  "ents_not_in_database.vtk", "VTK", "");
249  }
250  CHKERR cutAndTrim(bit_level1, bit_level2, th, tol_cut, tol_cut_close,
251  tol_trim, tol_trim_close, &fixed_edges, &corner_nodes,
252  update_meshsets, debug);
253  if(debug) {
254  CHKERR cOre.getInterface<BitRefManager>()->writeEntitiesNotInDatabase(
255  "cut_trim_ents_not_in_database.vtk", "VTK", "");
256  }
257 
258  CHKERR mergeBadEdges(fraction_level, bit_level2, bit_level1, bit_level3,
259  getNewTrimSurfaces(), fixed_edges, corner_nodes, th,
260  update_meshsets, debug);
261  // CHKERR removePathologicalFrontTris(bit_level3,
262  // const_cast<Range &>(getMergedSurfaces()));
263 
264  if(debug) {
265  CHKERR cOre.getInterface<BitRefManager>()->writeBitLevelByType(
266  bit_level3, BitRefLevel().set(), MBTET, "out_tets_merged.vtk", "VTK",
267  "");
268  CHKERR cOre.getInterface<BitRefManager>()->writeEntitiesNotInDatabase(
269  "cut_trim_merge_ents_not_in_database.vtk", "VTK", "");
270  }
271 
272  auto get_min_quality = [&m_field, debug](const BitRefLevel bit, Tag th) {
273  Range tets_level; // test at level
274  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
275  bit, BitRefLevel().set(), MBTET, tets_level);
276  double min_q = 1;
277  CHKERR m_field.getInterface<Tools>()->minTetsQuality(tets_level, min_q, th);
278  if (min_q < 0 && debug) {
279  CHKERR m_field.getInterface<Tools>()->writeTetsWithQuality(
280  "negative_tets.vtk", "VTK", "", tets_level, th);
281  }
282  return min_q;
283  };
284 
285  PetscPrintf(PETSC_COMM_WORLD, "Min quality node merge %6.4g\n",
286  get_min_quality(bit_level1, th));
287 
288  CHKERR
289  cOre.getInterface<BitRefManager>()->updateRange(fixed_edges, fixed_edges);
290  CHKERR cOre.getInterface<BitRefManager>()->updateRange(corner_nodes,
291  corner_nodes);
292 
294 }
295 
296 static double get_ave_edge_length(const EntityHandle ent,
297  const Range &vol_edges,
298  moab::Interface &moab) {
299  Range adj_edges;
300  if (moab.type_from_handle(ent) == MBVERTEX) {
301  CHKERR moab.get_adjacencies(&ent, 1, 1, false, adj_edges);
302  } else {
303  Range nodes;
304  CHKERR moab.get_connectivity(&ent, 1, nodes);
305  CHKERR moab.get_adjacencies(&ent, 1, 1, false, adj_edges,
306  moab::Interface::UNION);
307  }
308  adj_edges = intersect(adj_edges, vol_edges);
309  double ave_l = 0;
310  for (auto e : adj_edges) {
311  int num_nodes;
312  const EntityHandle *conn;
313  CHKERR moab.get_connectivity(e, conn, num_nodes, true);
314  VectorDouble6 coords(6);
315  CHKERR moab.get_coords(conn, num_nodes, &coords[0]);
317  &coords[0], &coords[1], &coords[2]);
319  &coords[3], &coords[4], &coords[5]);
321  t_n0(i) -= t_n1(i);
322  ave_l += sqrt(t_n0(i) * t_n0(i));
323  }
324  return ave_l / adj_edges.size();
325 };
326 
328  Range *corner_nodes,
329  const double low_tol, int verb,
330  const bool debug) {
331  CoreInterface &m_field = cOre;
332  moab::Interface &moab = m_field.get_moab();
334 
335  edgesToCut.clear();
336  cutEdges.clear();
337 
338  zeroDistanceVerts.clear();
339  zeroDistanceEnts.clear();
340  verticesOnCutEdges.clear();
341 
342  double ray_length;
343  double ray_point[3], unit_ray_dir[3];
344  VectorAdaptor vec_unit_ray_dir(
345  3, ublas::shallow_array_adaptor<double>(3, unit_ray_dir));
346  VectorAdaptor vec_ray_point(
347  3, ublas::shallow_array_adaptor<double>(3, ray_point));
348 
349  Tag th_dist;
350  rval = moab.tag_get_handle("DIST",th_dist);
351  if(rval == MB_SUCCESS) {
352  CHKERR moab.tag_delete(th_dist);
353  } else {
354  rval = MB_SUCCESS;
355  }
356  Tag th_dist_normal;
357  rval = moab.tag_get_handle("DIST_NORMAL", th_dist_normal);
358  if(rval == MB_SUCCESS) {
359  CHKERR moab.tag_delete(th_dist_normal);
360  } else {
361  rval = MB_SUCCESS;
362  }
363 
364  double def_val[] = {0, 0, 0};
365  CHKERR moab.tag_get_handle("DIST", 1, MB_TYPE_DOUBLE, th_dist,
366  MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
367  CHKERR moab.tag_get_handle("DIST_NORMAL", 3, MB_TYPE_DOUBLE, th_dist_normal,
368  MB_TAG_CREAT | MB_TAG_SPARSE, def_val);
369 
370  Range vol_vertices;
371  CHKERR moab.get_connectivity(vOlume, vol_vertices, true);
372  for (auto v : vol_vertices) {
373  VectorDouble3 point_in(3);
374  CHKERR moab.get_coords(&v, 1, &point_in[0]);
375  VectorDouble3 point_out(3);
376  EntityHandle facets_out;
377  CHKERR treeSurfPtr->closest_to_location(&point_in[0], rootSetSurf,
378  &point_out[0], facets_out);
379  VectorDouble3 n(3);
380  Util::normal(&moab, facets_out, n[0], n[1], n[2]);
381  VectorDouble3 delta = point_out - point_in;
382  double dist = norm_2(delta);
383  VectorDouble3 dist_normal = inner_prod(delta, n) * n;
384  CHKERR moab.tag_set_data(th_dist, &v, 1, &dist);
385  CHKERR moab.tag_set_data(th_dist_normal, &v, 1, &dist_normal[0]);
386  }
387 
388  auto get_normal_dist = [](const double *normal) {
389  FTensor::Tensor1<double, 3> t_n(normal[0], normal[1], normal[2]);
391  return sqrt(t_n(i) * t_n(i));
392  };
393 
394  auto get_edge_crossed = [&moab, get_normal_dist,
395  th_dist_normal](EntityHandle v0, EntityHandle v1) {
396  VectorDouble3 dist_normal0(3);
397  CHKERR moab.tag_get_data(th_dist_normal, &v0, 1, &dist_normal0[0]);
398  VectorDouble3 dist_normal1(3);
399  CHKERR moab.tag_get_data(th_dist_normal, &v1, 1, &dist_normal1[0]);
400  return (inner_prod(dist_normal0, dist_normal1) < 0);
401  };
402 
403  auto get_normal_dist_from_conn = [&moab, get_normal_dist,
404  th_dist_normal](EntityHandle v) {
405  double dist_normal[3];
406  CHKERR moab.tag_get_data(th_dist_normal, &v, 1, dist_normal);
407  return get_normal_dist(dist_normal);
408  };
409 
410  auto project_node = [this, &moab, th_dist_normal](const EntityHandle v) {
412  VectorDouble3 dist_normal(3);
413  rval = moab.tag_get_data(th_dist_normal, &v, 1, &dist_normal[0]);
414  VectorDouble3 s0(3);
415  CHKERR moab.get_coords(&v, 1, &s0[0]);
416  double dist = norm_2(dist_normal);
417  verticesOnCutEdges[v].dIst = dist;
418  verticesOnCutEdges[v].lEngth = dist;
419  verticesOnCutEdges[v].unitRayDir =
420  dist > 0 ? dist_normal / dist : dist_normal;
421  verticesOnCutEdges[v].rayPoint = s0;
423  };
424 
425  auto not_project_node = [this, &moab](const EntityHandle v) {
427  VectorDouble3 s0(3);
428  CHKERR moab.get_coords(&v, 1, &s0[0]);
429  verticesOnCutEdges[v].dIst = 0;
430  verticesOnCutEdges[v].lEngth = 0;
431  verticesOnCutEdges[v].unitRayDir = s0;
432  verticesOnCutEdges[v].rayPoint = s0;
434  };
435 
436  auto check_if_is_on_fixed_edge = [fixed_edges](const EntityHandle e) {
437  if(fixed_edges) {
438  if(fixed_edges->find(e)!=fixed_edges->end()) {
439  return true;
440  } else {
441  return false;
442  }
443  } else {
444  return false;
445  }
446  };
447 
448  auto check_if_is_on_cornet_node = [corner_nodes](const EntityHandle v) {
449  if (corner_nodes) {
450  if (corner_nodes->find(v) != corner_nodes->end()) {
451  return true;
452  } else {
453  return false;
454  }
455  } else {
456  return false;
457  }
458  };
459 
460  Range vol_edges;
461  CHKERR moab.get_adjacencies(vOlume, 1, true, vol_edges,
462  moab::Interface::UNION);
463 
464  aveLength = 0;
465  maxLength = 0;
466  int nb_ave_length = 0;
467  for (auto e : vol_edges) {
468  int num_nodes;
469  const EntityHandle *conn;
470  CHKERR moab.get_connectivity(e, conn, num_nodes, true);
471  double dist[num_nodes];
472  CHKERR moab.tag_get_data(th_dist, conn, num_nodes, dist);
473  CHKERR getRayForEdge(e, vec_ray_point, vec_unit_ray_dir, ray_length);
474  const double tol = ray_length * low_tol;
475  if (get_edge_crossed(conn[0], conn[1])) {
476  std::vector<double> distances_out;
477  std::vector<EntityHandle> facets_out;
478  CHKERR treeSurfPtr->ray_intersect_triangles(distances_out, facets_out,
479  rootSetSurf, tol, ray_point,
480  unit_ray_dir, &ray_length);
481  if (!distances_out.empty()) {
482  const auto dist_ptr =
483  std::min_element(distances_out.begin(), distances_out.end());
484  const double dist = *dist_ptr;
485  if (dist <= ray_length) {
486  aveLength += ray_length;
487  maxLength = fmax(maxLength, ray_length);
488  nb_ave_length++;
489  edgesToCut[e].dIst = dist;
490  edgesToCut[e].lEngth = ray_length;
491  edgesToCut[e].unitRayDir = vec_unit_ray_dir;
492  edgesToCut[e].rayPoint = vec_ray_point;
493  cutEdges.insert(e);
494  }
495  }
496  }
497 
498  if (fabs(dist[0]) < tol && fabs(dist[1]) < tol) {
499  aveLength += ray_length;
500  maxLength = fmax(maxLength, ray_length);
501  if (check_if_is_on_fixed_edge(e)) {
502  CHKERR not_project_node(conn[0]);
503  CHKERR not_project_node(conn[1]);
504  } else {
505  CHKERR project_node(conn[0]);
506  CHKERR project_node(conn[1]);
507  }
508  zeroDistanceEnts.insert(e);
509  }
510  }
511  aveLength /= nb_ave_length;
512 
513  Range cut_edges_verts;
514  CHKERR moab.get_connectivity(unite(cutEdges, zeroDistanceEnts),
515  cut_edges_verts, true);
516  vol_vertices = subtract(vol_vertices, cut_edges_verts);
517 
518  for (auto v : vol_vertices) {
519  double dist;
520  CHKERR moab.tag_get_data(th_dist, &v, 1, &dist);
521  const double tol = get_ave_edge_length(v, vol_edges, moab) * low_tol;
522  if(fabs(dist) < tol) {
523 
524  if (check_if_is_on_cornet_node(v)) {
525  CHKERR not_project_node(v);
526  } else {
527  CHKERR project_node(v);
528  }
529 
530  zeroDistanceVerts.insert(v);
531  }
532  }
533 
534  cutVolumes.clear();
535  // take all volumes adjacent to cut edges
536  CHKERR moab.get_adjacencies(cutEdges, 3, false, cutVolumes,
537  moab::Interface::UNION);
538  CHKERR moab.get_adjacencies(zeroDistanceVerts, 3, false, cutVolumes,
539  moab::Interface::UNION);
540  {
541  Range verts;
542  CHKERR moab.get_connectivity(unite(cutEdges, zeroDistanceEnts), verts,
543  true);
544  CHKERR moab.get_adjacencies(verts, 3, false, cutVolumes,
545  moab::Interface::UNION);
546  }
547  cutVolumes = intersect(cutVolumes,vOlume);
548 
549  // get edges on the cut volumes
550  Range edges;
551  CHKERR moab.get_adjacencies(cutVolumes, 1, false, edges,
552  moab::Interface::UNION);
553  edges = subtract(edges, cutEdges);
554 
555  // add to cut set edges which are cutted by extension of cutting surface
556  for (auto e : edges) {
557  int num_nodes;
558  const EntityHandle *conn;
559  CHKERR moab.get_connectivity(e, conn, num_nodes, true);
560  const double tol = get_ave_edge_length(e, vol_edges, moab) * low_tol;
561  double dist_normal[2];
562  dist_normal[0] = get_normal_dist_from_conn(conn[0]);
563  dist_normal[1] = get_normal_dist_from_conn(conn[1]);
564  if (get_edge_crossed(conn[0], conn[1])) {
565  CHKERR getRayForEdge(e, vec_ray_point, vec_unit_ray_dir, ray_length);
566  double s =
567  fabs(dist_normal[0]) / (fabs(dist_normal[0]) + fabs(dist_normal[1]));
568  edgesToCut[e].dIst = s * ray_length;
569  edgesToCut[e].lEngth = ray_length;
570  edgesToCut[e].unitRayDir = vec_unit_ray_dir;
571  edgesToCut[e].rayPoint = vec_ray_point;
572  cutEdges.insert(e);
573  } else if (fabs(dist_normal[0]) < tol && fabs(dist_normal[1]) < tol) {
574  if (check_if_is_on_fixed_edge(e)) {
575  CHKERR not_project_node(conn[0]);
576  CHKERR not_project_node(conn[1]);
577  } else {
578  CHKERR project_node(conn[0]);
579  CHKERR project_node(conn[1]);
580  }
581  zeroDistanceEnts.insert(e);
582  }
583  }
584 
585  CHKERR moab.get_adjacencies(cutVolumes, 1, false, edges,
586  moab::Interface::UNION);
587  Range add_verts;
588  CHKERR moab.get_connectivity(edges, add_verts, true);
589  add_verts = subtract(add_verts, zeroDistanceVerts);
590  CHKERR moab.get_connectivity(unite(cutEdges, zeroDistanceEnts),
591  cut_edges_verts, true);
592  add_verts = subtract(add_verts, cut_edges_verts);
593 
594  for (auto v : add_verts) {
595  double dist_normal = get_normal_dist_from_conn(v);
596  const double tol = get_ave_edge_length(v, vol_edges, moab) * low_tol;
597  if (fabs(dist_normal) < tol) {
598 
599  if (check_if_is_on_cornet_node(v)) {
600  CHKERR not_project_node(v);
601  } else {
602  CHKERR project_node(v);
603  }
604 
605  zeroDistanceVerts.insert(v);
606  }
607  }
608 
609  for (auto f : zeroDistanceEnts) {
610  int num_nodes;
611  const EntityHandle *conn;
612  CHKERR moab.get_connectivity(f, conn, num_nodes, true);
613  Range adj_edges;
614  CHKERR moab.get_adjacencies(conn, num_nodes, 1, false, adj_edges,
615  moab::Interface::UNION);
616  for (auto e : adj_edges) {
617  cutEdges.erase(e);
618  edgesToCut.erase(e);
619  }
620  }
621 
622  for (auto v : zeroDistanceVerts) {
623  Range adj_edges;
624  CHKERR moab.get_adjacencies(&v, 1, 1, false, adj_edges,
625  moab::Interface::UNION);
626  for (auto e : adj_edges) {
627  cutEdges.erase(e);
628  edgesToCut.erase(e);
629  }
630  }
631 
632  if(debug) {
633  CHKERR SaveData(m_field.get_moab())("cut_edges.vtk", cutEdges);
634  CHKERR SaveData(m_field.get_moab())(
635  "zero_dist_ents.vtk", unite(zeroDistanceVerts, zeroDistanceVerts));
636  }
637 
639 }
640 
642  Range *corner_nodes,
643  const double low_tol,
644  const int verb,
645  const bool debug) {
646  CoreInterface &m_field = cOre;
647  moab::Interface &moab = m_field.get_moab();
649 
650  // Get entities on body skin
651  Skinner skin(&moab);
652  Range tets_skin;
653  rval = skin.find_skin(0, vOlume, false, tets_skin);
654  Range tets_skin_edges;
655  CHKERR moab.get_adjacencies(tets_skin, 1, false, tets_skin_edges,
656  moab::Interface::UNION);
657  Range tets_skin_verts;
658  CHKERR moab.get_connectivity(tets_skin, tets_skin_verts, true);
659 
660  // Get entities in volume
661  Range vol_faces, vol_edges, vol_nodes;
662  CHKERR moab.get_adjacencies(vOlume, 2, false, vol_faces,
663  moab::Interface::UNION);
664  CHKERR moab.get_adjacencies(vOlume, 1, false, vol_edges,
665  moab::Interface::UNION);
666  CHKERR moab.get_connectivity(vOlume, vol_nodes, true);
667 
668  // Get nodes on cut edges
669  Range cut_edge_verts;
670  CHKERR moab.get_connectivity(cutEdges, cut_edge_verts, true);
671 
672  Range fixed_edges_nodes;
673  if(fixed_edges) {
674  CHKERR moab.get_connectivity(*fixed_edges, fixed_edges_nodes, true);
675  }
676 
677  // Get faces and edges
678  Range cut_edges_faces;
679  CHKERR moab.get_adjacencies(cut_edge_verts, 2, true, cut_edges_faces,
680  moab::Interface::UNION);
681  cut_edges_faces = intersect(cut_edges_faces, vol_faces);
682  Range cut_edges_faces_verts;
683  CHKERR moab.get_connectivity(cut_edges_faces, cut_edges_faces_verts, true);
684  cut_edges_faces_verts = subtract(cut_edges_faces_verts, cut_edge_verts);
685  Range to_remove_cut_edges_faces;
686  CHKERR moab.get_adjacencies(cut_edges_faces_verts, 2, true,
687  to_remove_cut_edges_faces,
688  moab::Interface::UNION);
689  cut_edges_faces = subtract(cut_edges_faces, to_remove_cut_edges_faces);
690  cut_edges_faces.merge(cutEdges);
691 
692  Tag th_dist_normal;
693  CHKERR moab.tag_get_handle("DIST_NORMAL", th_dist_normal);
694 
695  auto get_quality_change =
696  [this, &m_field,
697  &moab](const Range &adj_tets,
698  map<EntityHandle, TreeData> vertices_on_cut_edges) {
699  vertices_on_cut_edges.insert(verticesOnCutEdges.begin(),
700  verticesOnCutEdges.end());
701  double q0 = 1;
702  CHKERR m_field.getInterface<Tools>()->minTetsQuality(adj_tets, q0);
703  double q = 1;
704  for (auto t : adj_tets) {
705  int num_nodes;
706  const EntityHandle *conn;
707  CHKERR m_field.get_moab().get_connectivity(t, conn, num_nodes, true);
708  VectorDouble12 coords(12);
709  CHKERR moab.get_coords(conn, num_nodes, &coords[0]);
710  // cerr << coords << endl;
711  for (int n = 0; n != 4; ++n) {
712  bool ray_found = false;
713  auto mit = vertices_on_cut_edges.find(conn[n]);
714  if (mit != vertices_on_cut_edges.end()) {
715  ray_found = true;
716  }
717  if (ray_found) {
718  auto n_coords = getVectorAdaptor(&coords[3 * n], 3);
719  double dist = mit->second.dIst;
720  noalias(n_coords) =
721  mit->second.rayPoint + dist * mit->second.unitRayDir;
722  }
723  }
724  q = std::min(q, Tools::volumeLengthQuality(&coords[0]));
725  }
726  if (std::isnormal(q))
727  return q / q0;
728  else
729  return -1.;
730  };
731 
732  auto get_conn = [&moab](const EntityHandle &e,
733  const EntityHandle *&conn, int &num_nodes) {
735  EntityType type = moab.type_from_handle(e);
736  if (type == MBVERTEX) {
737  conn = &e;
738  num_nodes = 1;
739  } else {
740  CHKERR moab.get_connectivity(e, conn, num_nodes, true);
741  }
743  };
744 
745  auto get_normal_dist = [](const double *normal) {
746  FTensor::Tensor1<double, 3> t_n(normal[0], normal[1], normal[2]);
748  return sqrt(t_n(i) * t_n(i));
749  };
750 
751  auto get_normal_dist_from_conn = [&moab, get_normal_dist,
752  th_dist_normal](EntityHandle v) {
753  double dist_normal[3];
754  CHKERR moab.tag_get_data(th_dist_normal, &v, 1, dist_normal);
755  return get_normal_dist(dist_normal);
756  };
757 
758  auto project_node = [&moab, th_dist_normal](
759  const EntityHandle v,
760  map<EntityHandle, TreeData> &vertices_on_cut_edges) {
762  VectorDouble3 dist_normal(3);
763  rval = moab.tag_get_data(th_dist_normal, &v, 1, &dist_normal[0]);
764  VectorDouble3 s0(3);
765  CHKERR moab.get_coords(&v, 1, &s0[0]);
766  double dist = norm_2(dist_normal);
767  vertices_on_cut_edges[v].dIst = dist;
768  vertices_on_cut_edges[v].lEngth = dist;
769  vertices_on_cut_edges[v].unitRayDir =
770  dist > 0 ? dist_normal / dist : dist_normal;
771  vertices_on_cut_edges[v].rayPoint = s0;
773  };
774 
775  for (int d = 2; d >= 0; --d) {
776 
777  Range ents;
778  if (d > 0)
779  ents = cut_edges_faces.subset_by_dimension(d);
780  else
781  ents = cut_edge_verts;
782 
783  // make list of entities
784  multimap<double, EntityHandle> ents_to_check;
785  for (auto f : ents) {
786  int num_nodes;
787  const EntityHandle *conn;
788  CHKERR get_conn(f, conn, num_nodes);
789  VectorDouble3 dist(3);
790  for (int n = 0; n != num_nodes; ++n) {
791  dist[n] = get_normal_dist_from_conn(conn[n]);
792  }
793  double max_dist = 0;
794  for (int n = 0; n != num_nodes; ++n) {
795  max_dist = std::max(max_dist, fabs(dist[n]));
796  }
797  if (max_dist < low_tol * get_ave_edge_length(f, vol_edges, moab)) {
798  ents_to_check.insert(std::pair<double, EntityHandle>(max_dist, f));
799  }
800  }
801 
802  double ray_point[3], unit_ray_dir[3];
803  VectorAdaptor vec_unit_ray_dir(
804  3, ublas::shallow_array_adaptor<double>(3, unit_ray_dir));
805  VectorAdaptor vec_ray_point(
806  3, ublas::shallow_array_adaptor<double>(3, ray_point));
807 
808  for (auto m : ents_to_check) {
809 
810  EntityHandle f = m.second;
811 
812  int num_nodes;
813  const EntityHandle *conn;
814  CHKERR get_conn(f, conn, num_nodes);
815  VectorDouble9 coords(9);
816  CHKERR moab.get_coords(conn, num_nodes, &coords[0]);
817 
818  Range adj_tets;
819  CHKERR moab.get_adjacencies(conn, num_nodes, 3, false, adj_tets,
820  moab::Interface::UNION);
821  adj_tets = intersect(adj_tets, vOlume);
822 
823  map<EntityHandle, TreeData> vertices_on_cut_edges;
824  for (int n = 0; n != num_nodes; ++n) {
825  const EntityHandle node = conn[n];
826  CHKERR project_node(node, vertices_on_cut_edges);
827  }
828  if (static_cast<int>(vertices_on_cut_edges.size()) != num_nodes) {
829  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
830  "Data inconsistency");
831  }
832 
833  double q = get_quality_change(adj_tets, vertices_on_cut_edges);
834  if (q > 0.75) {
835  bool check_q_again = false;
836  for (auto &m : vertices_on_cut_edges) {
837  EntityHandle node = m.first;
838  if (tets_skin_verts.find(node) != tets_skin_verts.end()) {
839 
840  check_q_again = true;
841 
842  // check if node is at the corner
843  bool zero_disp_node = false;
844  if (corner_nodes) {
845  if (corner_nodes->find(node) != corner_nodes->end()) {
846  zero_disp_node = true;
847  }
848  }
849 
850  // check node is on the fixed edge
851  Range adj_edges;
852  CHKERR moab.get_adjacencies(&node, 1, 1, false, adj_edges);
853  adj_edges = intersect(adj_edges, tets_skin_edges);
854  if (fixed_edges) {
855  Range e;
856  // check if node is on fixed edge
857  e = intersect(adj_edges, *fixed_edges);
858  if (!e.empty()) {
859  adj_edges.swap(e);
860  }
861  // check if split edge is fixed edge
862  e = intersect(adj_edges, cutEdges);
863  if (!e.empty()) {
864  adj_edges.swap(e);
865  } else {
866  zero_disp_node = true;
867  }
868  }
869 
870  VectorDouble3 s0(3);
871  CHKERR moab.get_coords(&node, 1, &s0[0]);
872 
873  if (zero_disp_node) {
874  VectorDouble3 z(3);
875  z.clear();
876  m.second.dIst = 0;
877  m.second.lEngth = 0;
878  m.second.unitRayDir = z;
879  m.second.rayPoint = s0;
880  } else {
881  if (adj_edges.empty()) {
882  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
883  "Data inconsistency");
884  }
885  for (auto e : adj_edges) {
886  if (edgesToCut.find(e) != edgesToCut.end()) {
887  auto d = edgesToCut.at(e);
888  VectorDouble3 new_pos = d.rayPoint + d.dIst * d.unitRayDir;
889  VectorDouble3 ray = new_pos - s0;
890  double dist0 = norm_2(ray);
891  m.second.dIst = dist0;
892  m.second.lEngth = dist0;
893  m.second.unitRayDir = dist0 > 0 ? ray / dist0 : ray;
894  m.second.rayPoint = s0;
895  break;
896  } else {
897  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
898  "Data inconsistency");
899  }
900  }
901  }
902  }
903  }
904 
905  if (check_q_again) {
906  q = get_quality_change(adj_tets, vertices_on_cut_edges);
907  }
908  if (q > 0.75) {
909  verticesOnCutEdges.insert(vertices_on_cut_edges.begin(),
910  vertices_on_cut_edges.end());
911  EntityHandle type = moab.type_from_handle(f);
912  if (type == MBVERTEX) {
913  zeroDistanceVerts.insert(f);
914  } else {
915  zeroDistanceEnts.insert(f);
916  }
917  }
918  }
919  }
920 
921  }
922 
923  for (auto f : unite(zeroDistanceEnts, zeroDistanceVerts)) {
924  int num_nodes;
925  const EntityHandle *conn;
926  CHKERR get_conn(f, conn, num_nodes);
927  Range adj_edges;
928  CHKERR moab.get_adjacencies(conn, num_nodes, 1, false, adj_edges,
929  moab::Interface::UNION);
930  for (auto e : adj_edges) {
931  cutEdges.erase(e);
932  edgesToCut.erase(e);
933  }
934  }
935 
936  if (debug) {
937  SaveData(m_field.get_moab())("projected_zero_distance_ents.vtk",
939  }
940 
942 }
943 
945  const bool debug) {
946  CoreInterface &m_field = cOre;
947  moab::Interface &moab = m_field.get_moab();
948  MeshRefinement *refiner;
949  const RefEntity_multiIndex *ref_ents_ptr;
951  if(cutEdges.size() != edgesToCut.size()) {
952  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY, "Data inconsistency");
953  }
954  CHKERR m_field.getInterface(refiner);
955  CHKERR m_field.get_ref_ents(&ref_ents_ptr);
957  CHKERR refiner->refine_TET(vOlume, bit, false, QUIET,
958  debug ? &cutEdges : NULL);
959  cutNewVolumes.clear();
960  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
961  bit, BitRefLevel().set(), MBTET, cutNewVolumes);
962  cutNewSurfaces.clear();
963  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
964  bit, BitRefLevel().set(), MBTRI, cutNewSurfaces);
965  // Find new vertices on cut edges
966  cutNewVertices.clear();
967  CHKERR moab.get_connectivity(zeroDistanceEnts, cutNewVertices, true);
969  for (map<EntityHandle, TreeData>::iterator mit = edgesToCut.begin();
970  mit != edgesToCut.end(); ++mit) {
971  RefEntity_multiIndex::index<
972  Composite_ParentEnt_And_EntType_mi_tag>::type::iterator vit =
973  ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().find(
974  boost::make_tuple(MBVERTEX, mit->first));
975  if (vit ==
976  ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().end()) {
977  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
978  "No vertex on cut edges, that make no sense");
979  }
980  const boost::shared_ptr<RefEntity> &ref_ent = *vit;
981  if ((ref_ent->getBitRefLevel() & bit).any()) {
982  EntityHandle vert = ref_ent->getRefEnt();
983  cutNewVertices.insert(vert);
984  verticesOnCutEdges[vert] = mit->second;
985  } else {
986  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
987  "Vertex has wrong bit ref level");
988  }
989  }
990  // Add zero distance entities faces
991  Range tets_skin;
992  Skinner skin(&moab);
993  CHKERR skin.find_skin(0, cutNewVolumes, false, tets_skin);
994  cutNewSurfaces.merge(zeroDistanceEnts.subset_by_type(MBTRI));
995  // At that point cutNewSurfaces has all newly created faces, now take all
996  // nodes on those faces and subtract nodes on catted edges. Faces adjacent to
997  // nodes which left are not part of surface.
998  Range diff_verts;
999  CHKERR moab.get_connectivity(unite(cutNewSurfaces, zeroDistanceEnts),
1000  diff_verts, true);
1001  diff_verts = subtract(diff_verts, cutNewVertices);
1002  Range subtract_faces;
1003  CHKERR moab.get_adjacencies(diff_verts, 2, false, subtract_faces,
1004  moab::Interface::UNION);
1005  cutNewSurfaces = subtract(cutNewSurfaces, unite(subtract_faces, tets_skin));
1006  cutNewVertices.clear();
1007  CHKERR moab.get_connectivity(cutNewSurfaces, cutNewVertices, true);
1008 
1010 }
1011 
1014 
1015  CoreInterface &m_field = cOre;
1016  moab::Interface &moab = m_field.get_moab();
1018 
1019  // Range out_side_vertices;
1020  for (map<EntityHandle, TreeData>::iterator mit = verticesOnCutEdges.begin();
1021  mit != verticesOnCutEdges.end(); mit++) {
1022  double dist = mit->second.dIst;
1023  VectorDouble3 new_coors =
1024  mit->second.rayPoint + dist * mit->second.unitRayDir;
1025  if (th) {
1026  CHKERR moab.tag_set_data(th, &mit->first, 1, &new_coors[0]);
1027  } else {
1028  CHKERR moab.set_coords(&mit->first, 1, &new_coors[0]);
1029  }
1030  }
1031 
1033 }
1034 
1036  CoreInterface &m_field = cOre;
1037  moab::Interface &moab = m_field.get_moab();
1039  // Range out_side_vertices;
1040  for (auto &v : verticesOnTrimEdges) {
1041  double dist = v.second.dIst;
1042  VectorDouble3 new_coors = v.second.rayPoint + dist * v.second.unitRayDir;
1043  if (th) {
1044  CHKERR moab.tag_set_data(th, &v.first, 1, &new_coors[0]);
1045  } else {
1046  CHKERR moab.set_coords(&v.first, 1, &new_coors[0]);
1047  }
1048  }
1050 }
1051 
1053  Range *corner_nodes, Tag th,
1054  const double tol,
1055  const bool debug) {
1056  CoreInterface &m_field = cOre;
1057  moab::Interface &moab = m_field.get_moab();
1059 
1060  // takes body skin
1061  Skinner skin(&moab);
1062  Range tets_skin;
1063  CHKERR skin.find_skin(0, cutNewVolumes, false, tets_skin);
1064  // vertives on the skin
1065  Range tets_skin_verts;
1066  CHKERR moab.get_connectivity(tets_skin, tets_skin_verts, true);
1067  // edges on the skin
1068  Range tets_skin_edges;
1069  CHKERR moab.get_adjacencies(tets_skin, 1, false, tets_skin_edges,
1070  moab::Interface::UNION);
1071  // get edges on new surface
1072  Range edges;
1073  CHKERR moab.get_adjacencies(cutNewSurfaces, 1, false, edges,
1074  moab::Interface::UNION);
1075 
1076  Range cut_surface_edges_on_fixed_edges;
1077  if (fixed_edges) {
1078  cut_surface_edges_on_fixed_edges = intersect(edges, *fixed_edges);
1079  }
1080  Range cut_surface_edges_on_fixed_edges_verts;
1081  CHKERR moab.get_connectivity(cut_surface_edges_on_fixed_edges,
1082  cut_surface_edges_on_fixed_edges_verts, true);
1083 
1084  Range fixed_edges_verts_on_corners;
1085  if (fixed_edges) {
1086  CHKERR moab.get_connectivity(*fixed_edges, fixed_edges_verts_on_corners,
1087  true);
1088  }
1089  fixed_edges_verts_on_corners = subtract(
1090  fixed_edges_verts_on_corners, cut_surface_edges_on_fixed_edges_verts);
1091  if (corner_nodes) {
1092  fixed_edges_verts_on_corners.merge(*corner_nodes);
1093  }
1094 
1095  auto get_point_coords = [&th, &moab](EntityHandle v) {
1096  VectorDouble3 point(3);
1097  if (th) {
1098  CHKERR moab.tag_get_data(th, &v, 1, &point[0]);
1099  } else {
1100  CHKERR moab.get_coords(&v, 1, &point[0]);
1101  }
1102  return point;
1103  };
1104 
1105  // clear data ranges
1106  trimEdges.clear();
1107  edgesToTrim.clear();
1108  verticesOnTrimEdges.clear();
1109  trimNewVertices.clear();
1110 
1111  // iterate over entities on new cut surface
1112  std::multimap<double, std::pair<EntityHandle, EntityHandle>> verts_map;
1113  for (auto e : edges) {
1114  // Get edge connectivity and coords
1115  int num_nodes;
1116  const EntityHandle *conn;
1117  CHKERR moab.get_connectivity(e, conn, num_nodes, true);
1118  double coords[3 * num_nodes];
1119  if (th) {
1120  CHKERR moab.tag_get_data(th, conn, num_nodes, coords);
1121  } else {
1122  CHKERR moab.get_coords(conn, num_nodes, coords);
1123  }
1124  // Put edges coords into boost vectors
1125  auto get_s_adaptor = [&coords](const int n) {
1126  return VectorAdaptor(3,
1127  ublas::shallow_array_adaptor<double>(3, &coords[n]));
1128  };
1129  VectorAdaptor s0 = get_s_adaptor(0);
1130  VectorAdaptor s1 = get_s_adaptor(3);
1131  // get edge length
1132  double length = norm_2(s1 - s0);
1133 
1134  // Find point on surface closet to surface
1135  auto get_closets_delta = [this, &moab](const VectorAdaptor &s) {
1136  VectorDouble3 p(3);
1137  EntityHandle facets_out;
1138  // find closet point on the surface from first node
1139  CHKERR treeSurfPtr->closest_to_location(&s[0], rootSetSurf, &p[0],
1140  facets_out);
1141  VectorDouble3 n(3);
1142  Util::normal(&moab, facets_out, n[0], n[1], n[2]);
1143  VectorDouble3 w = p - s;
1144  VectorDouble3 normal = inner_prod(w, n) * n;
1145  w -= normal;
1146  return w;
1147  };
1148 
1149  // Calculate deltas, i.e. vectors from edges to closet point on surface
1150  VectorDouble3 delta0(3), delta1(3);
1151  noalias(delta0) = get_closets_delta(s0);
1152  noalias(delta1) = get_closets_delta(s1);
1153 
1154  // Calculate distances
1155  double dist0 = norm_2(delta0);
1156  double dist1 = norm_2(delta1);
1157  double min_dist = fmin(dist0, dist1);
1158  double max_dist = fmax(dist0, dist1);
1159 
1160  // add edge to trim
1161  double dist;
1162  VectorDouble3 ray;
1163  VectorDouble3 trimmed_end;
1164  VectorDouble3 itersection_point;
1165 
1166  if (min_dist < 1e-6 * aveLength && max_dist >= 1e-6 * aveLength) {
1167  if (max_dist == dist0) {
1168  // move mid node in reference to node 0
1169  trimmed_end = s0;
1170  ray = s1 - trimmed_end;
1171  } else {
1172  // move node in reference to node 1
1173  trimmed_end = s1;
1174  ray = s0 - trimmed_end;
1175  }
1176 
1177  // Solve nonlinera problem of finding point on surface front
1178  auto closest_point_projection =
1179  [this, &moab](VectorDouble3 ray, VectorDouble3 trimmed_end,
1180  const int max_it, const double tol) {
1181  VectorDouble3 n(3), w(3), normal(3);
1182  double length = norm_2(ray);
1183  ray /= length;
1184  for (int ii = 0; ii != max_it; ii++) {
1185  EntityHandle facets_out;
1186  VectorDouble3 point_out(3);
1187  treeSurfPtr->closest_to_location(&trimmed_end[0], rootSetSurf,
1188  &point_out[0], facets_out);
1189  Util::normal(&moab, facets_out, n[0], n[1], n[2]);
1190  noalias(w) = point_out - trimmed_end;
1191  noalias(normal) = inner_prod(w, n) * n;
1192  double s = inner_prod(ray, w - normal);
1193  trimmed_end += s * ray;
1194  if ((s / length) < tol)
1195  break;
1196  }
1197  return trimmed_end;
1198  };
1199 
1200  itersection_point = closest_point_projection(
1201  ray, trimmed_end, nbMaxTrimSearchIterations, 1e-12);
1202 
1203  ray = itersection_point - trimmed_end;
1204  dist = norm_2(ray);
1205 
1206  if ((1 - dist / length) > 0) {
1207 
1208  // check if edges should be trimmed, i.e. if edge is trimmed at very
1209  // end simply move closed node rather than trim
1210  edgesToTrim[e].dIst = dist;
1211  edgesToTrim[e].lEngth = dist;
1212  edgesToTrim[e].unitRayDir = ray / dist;
1213  edgesToTrim[e].rayPoint = trimmed_end;
1214  trimEdges.insert(e);
1215 
1216  auto add_vert = [&verts_map, e](EntityHandle v, double dist) {
1217  verts_map.insert(
1218  std::pair<double, std::pair<EntityHandle, EntityHandle>>(
1219  dist, std::pair<EntityHandle, EntityHandle>(v, e)));
1220  };
1221 
1222  double dist0_to_intersection =
1223  norm_2(itersection_point - s0) / aveLength;
1224  double dist1_to_intersection =
1225  norm_2(itersection_point - s1) / aveLength;
1226  if (dist0_to_intersection < dist1_to_intersection) {
1227  add_vert(conn[0], dist0_to_intersection);
1228  } else {
1229  add_vert(conn[1], dist1_to_intersection);
1230  }
1231  }
1232  }
1233  }
1234 
1235  // Iterate over all vertives close to surface foront and check if those can be
1236  // moved
1237  for (auto m : verts_map) {
1238 
1239  if (m.first < tol) {
1240 
1241  EntityHandle v = m.second.first;
1242  if (verticesOnTrimEdges.find(v) != verticesOnTrimEdges.end()) {
1243  continue;
1244  }
1245 
1246  VectorDouble3 ray_point = get_point_coords(v);
1247  Range adj_edges;
1248  CHKERR moab.get_adjacencies(&v, 1, 1, false, adj_edges);
1249  adj_edges = intersect(adj_edges, edges);
1250  Range w = intersect(adj_edges, tets_skin_edges);
1251  if (!w.empty()) {
1252  adj_edges.swap(w);
1253  }
1254  if (fixed_edges) {
1255  Range r = intersect(adj_edges, *fixed_edges);
1256  if (!r.empty()) {
1257  adj_edges.swap(r);
1258  }
1259  }
1260  if (adj_edges.empty()) {
1261  SETERRQ(PETSC_COMM_WORLD, MOFEM_DATA_INCONSISTENCY, "Imposible case");
1262  }
1263 
1264  EntityHandle e = m.second.second;
1265  if (adj_edges.find(e) == adj_edges.end()) {
1266  continue;
1267  }
1268 
1269  if (fixed_edges_verts_on_corners.find(v) !=
1270  fixed_edges_verts_on_corners.end()) {
1271 
1272  if (edgesToTrim.find(e) != edgesToTrim.end()) {
1273  auto &m = edgesToTrim.at(e);
1274  verticesOnTrimEdges[v] = m;
1275  verticesOnTrimEdges[v].rayPoint = get_point_coords(v);
1276  verticesOnTrimEdges[v].dIst = 0;
1277  trimNewVertices.insert(v);
1278  }
1279 
1280  } else {
1281 
1282  VectorDouble3 new_pos(3);
1283  new_pos.clear();
1284  if (edgesToTrim.find(e) != edgesToTrim.end()) {
1285 
1286  auto &r = edgesToTrim.at(e);
1287  noalias(new_pos) = r.rayPoint + r.dIst * r.unitRayDir;
1288  VectorDouble3 unit_ray_dir = ray_point - new_pos;
1289  double dist = norm_2(unit_ray_dir);
1290  unit_ray_dir /= dist;
1291 
1292  auto get_quality_change = [this, &m_field, &moab, &new_pos, v,
1293  th](const Range &adj_tets) {
1294  double q0 = 1;
1295  CHKERR m_field.getInterface<Tools>()->minTetsQuality(adj_tets, q0);
1296  double q = 1;
1297  for (auto t : adj_tets) {
1298  int num_nodes;
1299  const EntityHandle *conn;
1300  CHKERR m_field.get_moab().get_connectivity(t, conn, num_nodes,
1301  true);
1302  VectorDouble12 coords(12);
1303  if (th) {
1304  CHKERR moab.tag_get_data(th, conn, num_nodes, &coords[0]);
1305  } else {
1306  CHKERR moab.get_coords(conn, num_nodes, &coords[0]);
1307  }
1308  for (int n = 0; n != 4; ++n) {
1309  auto n_coords = getVectorAdaptor(&coords[3 * n], 3);
1310  if (conn[n] == v) {
1311  noalias(n_coords) = new_pos;
1312  } else {
1313  auto m = verticesOnTrimEdges.find(conn[n]);
1314  if (m != verticesOnTrimEdges.end()) {
1315  auto r = m->second;
1316  noalias(n_coords) = r.rayPoint + r.dIst * r.unitRayDir;
1317  }
1318  }
1319  }
1320  q = std::min(q, Tools::volumeLengthQuality(&coords[0]));
1321  }
1322  return q / q0;
1323  };
1324 
1325  Range adj_tets;
1326  CHKERR moab.get_adjacencies(&v, 1, 3, false, adj_tets);
1327  adj_tets = intersect(adj_tets, cutNewVolumes);
1328  double q = get_quality_change(adj_tets);
1329  if (q > 0.75) {
1330  VectorDouble3 unit_ray_dir = new_pos - ray_point;
1331  double dist = norm_2(unit_ray_dir);
1332  unit_ray_dir /= dist;
1333  verticesOnTrimEdges[v].dIst = dist;
1334  verticesOnTrimEdges[v].lEngth = dist;
1335  verticesOnTrimEdges[v].unitRayDir = unit_ray_dir;
1336  verticesOnTrimEdges[v].rayPoint = ray_point;
1337  trimNewVertices.insert(v);
1338  }
1339  }
1340  }
1341  }
1342  }
1343 
1344  for (auto m : verticesOnTrimEdges) {
1345  EntityHandle v = m.first;
1346  Range adj_edges;
1347  CHKERR moab.get_adjacencies(&v, 1, 1, false, adj_edges);
1348  adj_edges = intersect(adj_edges, edges);
1349  for (auto e : adj_edges) {
1350  edgesToTrim.erase(e);
1351  trimEdges.erase(e);
1352  }
1353  }
1354 
1355  if (debug) {
1356  CHKERR SaveData(moab)("fixed_edges_verts_on_corners.vtk",
1357  fixed_edges_verts_on_corners);
1358  }
1359 
1361 }
1362 
1364  Tag th, const double tol,
1365  const bool debug) {
1366  CoreInterface &m_field = cOre;
1367  moab::Interface &moab = m_field.get_moab();
1368  MeshRefinement *refiner;
1369  const RefEntity_multiIndex *ref_ents_ptr;
1371 
1372  CHKERR m_field.getInterface(refiner);
1373  CHKERR m_field.get_ref_ents(&ref_ents_ptr);
1375  CHKERR refiner->refine_TET(cutNewVolumes, bit, false, QUIET,
1376  debug ? &trimEdges : NULL);
1377 
1378  trimNewVolumes.clear();
1379  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1380  bit, BitRefLevel().set(), MBTET, trimNewVolumes);
1381 
1382  for (map<EntityHandle, TreeData>::iterator mit = edgesToTrim.begin();
1383  mit != edgesToTrim.end(); mit++) {
1384  auto vit = ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().find(
1385  boost::make_tuple(MBVERTEX, mit->first));
1386  if (vit ==
1387  ref_ents_ptr->get<Composite_ParentEnt_And_EntType_mi_tag>().end()) {
1388  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1389  "No vertex on trim edges, that make no sense");
1390  }
1391  const boost::shared_ptr<RefEntity> &ref_ent = *vit;
1392  if ((ref_ent->getBitRefLevel() & bit).any()) {
1393  EntityHandle vert = ref_ent->getRefEnt();
1394  trimNewVertices.insert(vert);
1395  verticesOnTrimEdges[vert] = mit->second;
1396  }
1397  }
1398 
1399  // Get faces which are trimmed
1400  trimNewSurfaces.clear();
1401  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1402  bit, BitRefLevel().set(), MBTRI, trimNewSurfaces);
1403 
1404  Range trim_new_surfaces_nodes;
1405  CHKERR moab.get_connectivity(trimNewSurfaces, trim_new_surfaces_nodes, true);
1406  trim_new_surfaces_nodes = subtract(trim_new_surfaces_nodes, trimNewVertices);
1407  trim_new_surfaces_nodes = subtract(trim_new_surfaces_nodes, cutNewVertices);
1408  Range faces_not_on_surface;
1409  CHKERR moab.get_adjacencies(trim_new_surfaces_nodes, 2, false,
1410  faces_not_on_surface, moab::Interface::UNION);
1411  trimNewSurfaces = subtract(trimNewSurfaces, faces_not_on_surface);
1412 
1413  // Get surfaces which are not trimmed and add them to surface
1414  Range all_surfaces_on_bit_level;
1415  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1416  bit, BitRefLevel().set(), MBTRI, all_surfaces_on_bit_level);
1417  all_surfaces_on_bit_level =
1418  intersect(all_surfaces_on_bit_level, cutNewSurfaces);
1419  trimNewSurfaces = unite(trimNewSurfaces, all_surfaces_on_bit_level);
1420 
1421  Range trim_surface_edges;
1422  CHKERR moab.get_adjacencies(trimNewSurfaces, 1, false, trim_surface_edges,
1423  moab::Interface::UNION);
1424 
1425  // check of nodes are outside surface and if it are remove adjacent faces to
1426  // those nodes.
1427  Range check_verts;
1428  CHKERR moab.get_connectivity(trimNewSurfaces, check_verts, true);
1429  check_verts = subtract(check_verts, trimNewVertices);
1430  for (auto v : check_verts) {
1431 
1432  VectorDouble3 s(3);
1433  if (th) {
1434  CHKERR moab.tag_get_data(th, &v, 1, &s[0]);
1435  } else {
1436  CHKERR moab.get_coords(&v, 1, &s[0]);
1437  }
1438 
1439  VectorDouble3 p(3);
1440  EntityHandle facets_out;
1441  CHKERR treeSurfPtr->closest_to_location(&s[0], rootSetSurf, &p[0],
1442  facets_out);
1443  VectorDouble3 n(3);
1444  Util::normal(&moab, facets_out, n[0], n[1], n[2]);
1445  VectorDouble3 delta = s - p;
1446  VectorDouble3 normal = inner_prod(delta, n) * n;
1447  if (norm_2(delta - normal) > tol * aveLength) {
1448  Range adj;
1449  CHKERR moab.get_adjacencies(&v, 1, 2, false, adj);
1450  trimNewSurfaces = subtract(trimNewSurfaces, adj);
1451  }
1452  }
1453 
1455 }
1456 
1458  VectorAdaptor& ray_point,
1459  VectorAdaptor& unit_ray_dir,
1460  double &ray_length) const {
1461  const CoreInterface &m_field = cOre;
1462  const moab::Interface &moab = m_field.get_moab();
1464  int num_nodes;
1465  const EntityHandle *conn;
1466  CHKERR moab.get_connectivity(ent, conn, num_nodes, true);
1467  double coords[6];
1468  CHKERR moab.get_coords(conn, num_nodes, coords);
1469  VectorAdaptor s0(3, ublas::shallow_array_adaptor<double>(3, &coords[0]));
1470  VectorAdaptor s1(3, ublas::shallow_array_adaptor<double>(3, &coords[3]));
1471  noalias(ray_point) = s0;
1472  noalias(unit_ray_dir) = s1 - s0;
1473  ray_length = norm_2(unit_ray_dir);
1474  unit_ray_dir /= ray_length;
1476 }
1477 
1478 // int CutMeshInterface::segmentPlane(
1479 // VectorAdaptor s0,
1480 // VectorAdaptor s1,
1481 // VectorAdaptor x0,
1482 // VectorAdaptor n,
1483 // double &s
1484 // ) const {
1485 // VectorDouble3 u = s1 - s0;
1486 // VectorDouble3 w = s0 - x0;
1487 // double nu = inner_prod(n,u);
1488 // double nw = -inner_prod(n,w);
1489 // const double tol = 1e-4;
1490 // if (fabs(nu) < tol) { // segment is parallel to plane
1491 // if (nw == 0) // segment lies in plane
1492 // return 2;
1493 // else
1494 // return 0; // no intersection
1495 // }
1496 // // they are not parallel
1497 // // compute intersect param
1498 // s = nw / nu;
1499 // if (s < 0 || s > 1)
1500 // return 0; // no intersection
1501 // return 1;
1502 // }
1503 
1506  Range &ents) {
1507  CoreInterface &m_field = cOre;
1508  moab::Interface &moab = m_field.get_moab();
1509  PrismInterface *interface;
1511  CHKERR m_field.getInterface(interface);
1512  // Remove tris on surface front
1513  {
1514  Range front_tris;
1515  EntityHandle meshset;
1516  CHKERR moab.create_meshset(MESHSET_SET, meshset);
1517  CHKERR moab.add_entities(meshset, ents);
1519  meshset, split_bit, true, front_tris);
1520  CHKERR moab.delete_entities(&meshset, 1);
1521  ents = subtract(ents, front_tris);
1522  }
1523  // Remove entities on skin
1524  Range tets;
1525  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1526  split_bit,BitRefLevel().set(),MBTET,tets
1527  );
1528  // Remove entities on skin
1529  Skinner skin(&moab);
1530  Range tets_skin;
1531  rval = skin.find_skin(0, tets, false, tets_skin);
1532  ents = subtract(ents, tets_skin);
1533 
1535 }
1536 
1538  const BitRefLevel bit,
1539  const Range &ents, Tag th) {
1540  CoreInterface &m_field = cOre;
1541  moab::Interface &moab = m_field.get_moab();
1542  PrismInterface *interface;
1544  CHKERR m_field.getInterface(interface);
1545  EntityHandle meshset_volume;
1546  CHKERR moab.create_meshset(MESHSET_SET, meshset_volume);
1547  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1548  split_bit, BitRefLevel().set(), MBTET, meshset_volume);
1549  EntityHandle meshset_surface;
1550  CHKERR moab.create_meshset(MESHSET_SET, meshset_surface);
1551  CHKERR moab.add_entities(meshset_surface, ents);
1552  CHKERR interface->getSides(meshset_surface, split_bit, true);
1553  CHKERR interface->splitSides(meshset_volume, bit, meshset_surface, true,
1554  true);
1555  CHKERR moab.delete_entities(&meshset_volume, 1);
1556  CHKERR moab.delete_entities(&meshset_surface, 1);
1557  if (th) {
1558  Range prisms;
1559  ierr = m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
1560  bit, BitRefLevel().set(), MBPRISM, prisms);
1561  for (Range::iterator pit = prisms.begin(); pit != prisms.end(); pit++) {
1562  int num_nodes;
1563  const EntityHandle *conn;
1564  CHKERR moab.get_connectivity(*pit, conn, num_nodes, true);
1565  MatrixDouble data(3, 3);
1566  CHKERR moab.tag_get_data(th, conn, 3, &data(0, 0));
1567  // cerr << data << endl;
1568  CHKERR moab.tag_set_data(th, &conn[3], 3, &data(0, 0));
1569  }
1570  }
1572 }
1573 
1575  double lEngth;
1576  double qUality;
1578  bool skip;
1579  LengthMapData(const double l, double q, const EntityHandle e)
1580  : lEngth(l), qUality(q), eDge(e), skip(false) {}
1581 };
1582 
1583 typedef multi_index_container<
1584  LengthMapData,
1585  indexed_by<
1586  ordered_non_unique<
1587  member<LengthMapData, double, &LengthMapData::lEngth>
1588  >,
1589  hashed_unique<
1590  member<LengthMapData, EntityHandle, &LengthMapData::eDge>
1591  >
1592  >
1594 
1595 
1596 
1598  const int fraction_level, const Range &tets, const Range &surface,
1599  const Range &fixed_edges, const Range &corner_nodes, Range &edges_to_merge,
1600  Range &out_tets, Range &new_surf, Tag th, const bool update_meshsets,
1601  const BitRefLevel *bit_ptr,const bool debug) {
1602  CoreInterface &m_field = cOre;
1603  moab::Interface &moab = m_field.get_moab();
1605 
1606  /**
1607  * \brief Merge nodes
1608  */
1609  struct MergeNodes {
1610  CoreInterface &mField;
1611  const bool onlyIfImproveQuality;
1612  const int lineSearch;
1613  Tag tH;
1614  bool updateMehsets;
1615 
1616  MergeNodes(CoreInterface &m_field,
1617  const bool only_if_improve_quality, const int line_search,
1618  Tag th, bool update_mehsets)
1619  : mField(m_field), onlyIfImproveQuality(only_if_improve_quality),
1620  lineSearch(line_search), tH(th), updateMehsets(update_mehsets) {
1621  mField.getInterface(nodeMergerPtr);
1622  }
1623  NodeMergerInterface *nodeMergerPtr;
1624  MoFEMErrorCode operator()(EntityHandle father, EntityHandle mother,
1625  Range &proc_tets, Range &new_surf,
1626  Range &edges_to_merge, Range &not_merged_edges,
1627  bool add_child = true) const {
1628  moab::Interface &moab = mField.get_moab();
1630  const EntityHandle conn[] = {father, mother};
1631  Range vert_tets;
1632  CHKERR moab.get_adjacencies(conn, 2, 3, false, vert_tets,
1633  moab::Interface::UNION);
1634  vert_tets = intersect(vert_tets, proc_tets);
1635  Range out_tets;
1636  CHKERR nodeMergerPtr->mergeNodes(father, mother, out_tets, &vert_tets,
1637  onlyIfImproveQuality, 0, lineSearch, tH);
1638  out_tets.merge(subtract(proc_tets, vert_tets));
1639  proc_tets.swap(out_tets);
1640 
1641  if (add_child && nodeMergerPtr->getSucessMerge()) {
1642 
1643  NodeMergerInterface::ParentChildMap &parent_child_map =
1644  nodeMergerPtr->getParentChildMap();
1645 
1646  Range child_ents;
1647  NodeMergerInterface::ParentChildMap::iterator it;
1648  for (it = parent_child_map.begin(); it != parent_child_map.end();
1649  it++) {
1650  child_ents.insert(it->pArent);
1651  }
1652 
1653  Range new_surf_child_ents = intersect(new_surf, child_ents);
1654  new_surf = subtract(new_surf, new_surf_child_ents);
1655  Range child_surf_ents;
1656  CHKERR updateRangeByChilds(parent_child_map, new_surf_child_ents,
1657  child_surf_ents);
1658  new_surf.merge(child_surf_ents);
1659 
1660  Range edges_to_merge_child_ents = intersect(edges_to_merge, child_ents);
1661  edges_to_merge = subtract(edges_to_merge, edges_to_merge_child_ents);
1662  Range merged_child_edge_ents;
1663  CHKERR updateRangeByChilds(parent_child_map, edges_to_merge_child_ents,
1664  merged_child_edge_ents);
1665 
1666  Range not_merged_edges_child_ents =
1667  intersect(not_merged_edges, child_ents);
1668  not_merged_edges =
1669  subtract(not_merged_edges, not_merged_edges_child_ents);
1670  Range not_merged_child_edge_ents;
1671  CHKERR updateRangeByChilds(parent_child_map, not_merged_edges_child_ents,
1672  not_merged_child_edge_ents);
1673 
1674  merged_child_edge_ents =
1675  subtract(merged_child_edge_ents, not_merged_child_edge_ents);
1676  edges_to_merge.merge(merged_child_edge_ents);
1677  not_merged_edges.merge(not_merged_child_edge_ents);
1678 
1679  if (updateMehsets) {
1681  (*mField.getInterface<MeshsetsManager>()), cubit_it)) {
1682  EntityHandle cubit_meshset = cubit_it->meshset;
1683  Range parent_ents;
1684  CHKERR moab.get_entities_by_handle(cubit_meshset, parent_ents, true);
1685  Range child_ents;
1686  CHKERR updateRangeByChilds(parent_child_map, parent_ents, child_ents);
1687  CHKERR moab.add_entities(cubit_meshset, child_ents);
1688  }
1689  }
1690  }
1692  }
1693 
1694  private:
1695  MoFEMErrorCode updateRangeByChilds(
1696  const NodeMergerInterface::ParentChildMap &parent_child_map,
1697  const Range &parents, Range &childs) const {
1699  NodeMergerInterface::ParentChildMap::nth_index<0>::type::iterator it;
1700  for (Range::const_iterator eit = parents.begin(); eit != parents.end();
1701  eit++) {
1702  it = parent_child_map.get<0>().find(*eit);
1703  if (it == parent_child_map.get<0>().end())
1704  continue;
1705  childs.insert(it->cHild);
1706  }
1708  }
1709  };
1710 
1711  /**
1712  * \brief Calculate edge length
1713  */
1714  struct LengthMap {
1715  Tag tH;
1716  CoreInterface &mField;
1717  moab::Interface &moab;
1718  const double maxLength;
1719  LengthMap(CoreInterface &m_field, Tag th, double max_length)
1720  : tH(th), mField(m_field), moab(m_field.get_moab()),
1721  maxLength(max_length) {
1722  gammaL = 1.;
1723  gammaQ = 3.;
1724  acceptedThrasholdMergeQuality = 0.5;
1725  }
1726 
1727  double
1728  gammaL; ///< Controls importance of length when ranking edges for merge
1729  double
1730  gammaQ; ///< Controls importance of quality when ranking edges for merge
1731  double acceptedThrasholdMergeQuality; ///< Do not merge quality if quality
1732  ///< above accepted thrashold
1733 
1734  MoFEMErrorCode operator()(const Range &tets, const Range &edges,
1735  LengthMapData_multi_index &length_map,
1736  double &ave) const {
1737  int num_nodes;
1738  const EntityHandle *conn;
1739  double coords[6];
1741  VectorAdaptor s0(3, ublas::shallow_array_adaptor<double>(3, &coords[0]));
1742  VectorAdaptor s1(3, ublas::shallow_array_adaptor<double>(3, &coords[3]));
1743  VectorDouble3 delta(3);
1744  for (auto edge : edges) {
1745  CHKERR moab.get_connectivity(edge, conn, num_nodes, true);
1746  Range adj_tet;
1747  CHKERR moab.get_adjacencies(conn, num_nodes, 3, false, adj_tet);
1748  adj_tet = intersect(adj_tet, tets);
1749  if (tH) {
1750  CHKERR moab.tag_get_data(tH, conn, num_nodes, coords);
1751  } else {
1752  CHKERR moab.get_coords(conn, num_nodes, coords);
1753  }
1754  double q = 1;
1755  auto abs_min = [](double a, double b) {
1756  return std::min(fabs(a), fabs(b));
1757  };
1758  CHKERR mField.getInterface<Tools>()->minTetsQuality(adj_tet, q, tH,
1759  abs_min);
1760  if (q != q)
1761  SETERRQ(PETSC_COMM_SELF, MOFEM_DATA_INCONSISTENCY,
1762  "Quality not a number");
1763  if (fabs(q) > acceptedThrasholdMergeQuality)
1764  continue;
1765  noalias(delta) = (s0 - s1) / maxLength;
1766  double dot = inner_prod(delta, delta);
1767  double val = pow(q, gammaQ) * pow(dot, 0.5 * gammaL);
1768  length_map.insert(LengthMapData(val, q, edge));
1769  }
1770  ave = 0;
1771  for (LengthMapData_multi_index::nth_index<0>::type::iterator mit =
1772  length_map.get<0>().begin();
1773  mit != length_map.get<0>().end(); mit++) {
1774  ave += mit->qUality;
1775  }
1776  ave /= length_map.size();
1778  }
1779  };
1780 
1781  /**
1782  * \brief Topological relations
1783  */
1784  struct Toplogy {
1785 
1786  CoreInterface &mField;
1787  Tag tH;
1788  const double tOL;
1789  Toplogy(CoreInterface &m_field, Tag th, const double tol)
1790  : mField(m_field), tH(th), tOL(tol) {}
1791 
1792  enum TYPE {
1793  FREE = 0,
1794  SKIN = 1 << 0,
1795  SURFACE = 1 << 1,
1796  SURFACE_SKIN = 1 << 2,
1797  FRONT_ENDS = 1 << 3,
1798  FIX_EDGES = 1 << 4,
1799  FIX_CORNERS = 1 << 5
1800  };
1801 
1802  typedef map<int, Range> SetsMap;
1803 
1804  MoFEMErrorCode classifyVerts(const Range &surface, const Range &tets,
1805  const Range &fixed_edges,
1806  const Range &corner_nodes,
1807  SetsMap &sets_map) const {
1808  moab::Interface &moab(mField.get_moab());
1809  Skinner skin(&moab);
1811 
1812  sets_map[FIX_CORNERS].merge(corner_nodes);
1813  Range fixed_verts;
1814  CHKERR moab.get_connectivity(fixed_edges, fixed_verts, true);
1815  sets_map[FIX_EDGES].swap(fixed_verts);
1816 
1817  Range tets_skin;
1818  CHKERR skin.find_skin(0, tets, false, tets_skin);
1819  Range tets_skin_edges;
1820  CHKERR moab.get_adjacencies(tets_skin, 1, false, tets_skin_edges,
1821  moab::Interface::UNION);
1822 
1823  // surface skin
1824  Range surface_skin;
1825  CHKERR skin.find_skin(0, surface, false, surface_skin);
1826  Range front_in_the_body;
1827  front_in_the_body = subtract(surface_skin, tets_skin_edges);
1828  Range front_ends;
1829  CHKERR skin.find_skin(0, front_in_the_body, false, front_ends);
1830  sets_map[FRONT_ENDS].swap(front_ends);
1831 
1832  Range surface_skin_verts;
1833  CHKERR moab.get_connectivity(surface_skin, surface_skin_verts, true);
1834  sets_map[SURFACE_SKIN].swap(surface_skin_verts);
1835 
1836  // surface
1837  Range surface_verts;
1838  CHKERR moab.get_connectivity(surface, surface_verts, true);
1839  sets_map[SURFACE].swap(surface_verts);
1840 
1841  // skin
1842  Range tets_skin_verts;
1843  CHKERR moab.get_connectivity(tets_skin, tets_skin_verts, true);
1844  sets_map[SKIN].swap(tets_skin_verts);
1845 
1846  Range tets_verts;
1847  CHKERR moab.get_connectivity(tets, tets_verts, true);
1848  sets_map[FREE].swap(tets_verts);
1849 
1851  }
1852 
1853  MoFEMErrorCode getProcTets(const Range &tets, const Range &edges_to_merge,
1854  Range &proc_tets) const {
1855  moab::Interface &moab(mField.get_moab());
1857  Range edges_to_merge_verts;
1858  CHKERR moab.get_connectivity(edges_to_merge, edges_to_merge_verts, true);
1859  Range edges_to_merge_verts_tets;
1860  CHKERR moab.get_adjacencies(edges_to_merge_verts, 3, false,
1861  edges_to_merge_verts_tets,
1862  moab::Interface::UNION);
1863  edges_to_merge_verts_tets = intersect(edges_to_merge_verts_tets, tets);
1864  proc_tets.swap(edges_to_merge_verts_tets);
1866  }
1867 
1868  MoFEMErrorCode edgesToMerge(const Range &surface, const Range &tets,
1869  Range &edges_to_merge) const {
1870  moab::Interface &moab(mField.get_moab());
1872 
1873  Range surface_verts;
1874  CHKERR moab.get_connectivity(surface, surface_verts, true);
1875  Range surface_verts_edges;
1876  CHKERR moab.get_adjacencies(surface_verts, 1, false, surface_verts_edges,
1877  moab::Interface::UNION);
1878  edges_to_merge.merge(surface_verts_edges);
1879  Range tets_edges;
1880  CHKERR moab.get_adjacencies(tets, 1, false, tets_edges,
1881  moab::Interface::UNION);
1882  edges_to_merge = intersect(edges_to_merge, tets_edges);
1884  }
1885 
1886  MoFEMErrorCode removeBadEdges(const Range &surface, const Range &tets,
1887  const Range &fixed_edges,
1888  const Range &corner_nodes,
1889  Range &edges_to_merge,
1890  Range &not_merged_edges) {
1891  moab::Interface &moab(mField.get_moab());
1893 
1894  // find skin
1895  Skinner skin(&moab);
1896  Range tets_skin;
1897  CHKERR skin.find_skin(0, tets, false, tets_skin);
1898  Range surface_skin;
1899  CHKERR skin.find_skin(0, surface, false, surface_skin);
1900 
1901  // end nodes
1902  Range tets_skin_edges;
1903  CHKERR moab.get_adjacencies(tets_skin, 1, false, tets_skin_edges,
1904  moab::Interface::UNION);
1905  Range surface_front;
1906  surface_front = subtract(surface_skin, tets_skin_edges);
1907  Range surface_front_nodes;
1908  CHKERR moab.get_connectivity(surface_front, surface_front_nodes, true);
1909  Range ends_nodes;
1910  CHKERR skin.find_skin(0, surface_front, false, ends_nodes);
1911 
1912  // remove bad merges
1913 
1914  // get surface and body skin verts
1915  Range surface_edges;
1916  CHKERR moab.get_adjacencies(surface, 1, false, surface_edges,
1917  moab::Interface::UNION);
1918  // get nodes on the surface
1919  Range surface_edges_verts;
1920  CHKERR moab.get_connectivity(surface_edges, surface_edges_verts, true);
1921  // get vertices on the body skin
1922  Range tets_skin_edges_verts;
1923  CHKERR moab.get_connectivity(tets_skin_edges, tets_skin_edges_verts, true);
1924 
1925  Range edges_to_remove;
1926 
1927  // remove edges self connected to body skin
1928  {
1929  Range ents_nodes_and_edges;
1930  ents_nodes_and_edges.merge(tets_skin_edges_verts);
1931  ents_nodes_and_edges.merge(tets_skin_edges);
1932  CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
1933  0, false);
1934  }
1935  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
1936  not_merged_edges.merge(edges_to_remove);
1937 
1938  // remove edges self connected to surface
1939  {
1940  Range ents_nodes_and_edges;
1941  ents_nodes_and_edges.merge(surface_edges_verts);
1942  ents_nodes_and_edges.merge(surface_edges);
1943  ents_nodes_and_edges.merge(tets_skin_edges_verts);
1944  ents_nodes_and_edges.merge(tets_skin_edges);
1945  CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
1946  0, false);
1947  }
1948  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
1949  not_merged_edges.merge(edges_to_remove);
1950 
1951  // remove edges adjacent corner_nodes execpt those on fixed edges
1952  Range fixed_edges_nodes;
1953  CHKERR moab.get_connectivity(fixed_edges, fixed_edges_nodes, true);
1954  {
1955  Range ents_nodes_and_edges;
1956  ents_nodes_and_edges.merge(fixed_edges_nodes);
1957  ents_nodes_and_edges.merge(ends_nodes);
1958  ents_nodes_and_edges.merge(corner_nodes);
1959  ents_nodes_and_edges.merge(fixed_edges);
1960  CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
1961  0, false);
1962  }
1963  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
1964  not_merged_edges.merge(edges_to_remove);
1965 
1966  // remove edges self connected to surface
1967  CHKERR removeSelfConectingEdges(surface_edges, edges_to_remove, 0, false);
1968  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
1969  not_merged_edges.merge(edges_to_remove);
1970 
1971  // remove edges self contented on surface skin
1972  {
1973  Range ents_nodes_and_edges;
1974  ents_nodes_and_edges.merge(surface_skin);
1975  ents_nodes_and_edges.merge(fixed_edges_nodes);
1976  CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
1977  0, false);
1978  }
1979  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
1980  not_merged_edges.merge(edges_to_remove);
1981 
1982  // remove crack front nodes connected to the surface
1983  {
1984  Range ents_nodes_and_edges;
1985  ents_nodes_and_edges.merge(surface_front_nodes);
1986  ents_nodes_and_edges.merge(surface_front);
1987  ents_nodes_and_edges.merge(tets_skin_edges_verts);
1988  ents_nodes_and_edges.merge(tets_skin_edges);
1989  CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
1990  0, false);
1991  }
1992  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
1993  not_merged_edges.merge(edges_to_remove);
1994 
1995  // remove edges connecting crack front and fixed edges, except those
1996  {
1997  Range ents_nodes_and_edges;
1998  ents_nodes_and_edges.merge(surface_skin.subset_by_type(MBEDGE));
1999  ents_nodes_and_edges.merge(fixed_edges.subset_by_type(MBEDGE));
2000  CHKERR removeSelfConectingEdges(ents_nodes_and_edges, edges_to_remove,
2001  tOL, false);
2002  }
2003  edges_to_merge = subtract(edges_to_merge, edges_to_remove);
2004  not_merged_edges.merge(edges_to_remove);
2005 
2006 
2008  }
2009 
2010  private:
2011  MoFEMErrorCode removeSelfConectingEdges(const Range &ents,
2012  Range &edges_to_remove,
2013  const bool length,
2014  bool debug) const {
2015  moab::Interface &moab(mField.get_moab());
2017  // get nodes
2018  Range ents_nodes = ents.subset_by_type(MBVERTEX);
2019  if (ents_nodes.empty()) {
2020  CHKERR moab.get_connectivity(ents, ents_nodes, true);
2021  }
2022  // edges adj. to nodes
2023  Range ents_nodes_edges;
2024  CHKERR moab.get_adjacencies(ents_nodes, 1, false, ents_nodes_edges,
2025  moab::Interface::UNION);
2026  // nodes of adj. edges
2027  Range ents_nodes_edges_nodes;
2028  CHKERR moab.get_connectivity(ents_nodes_edges, ents_nodes_edges_nodes,
2029  true);
2030  // hanging nodes
2031  ents_nodes_edges_nodes = subtract(ents_nodes_edges_nodes, ents_nodes);
2032  Range ents_nodes_edges_nodes_edges;
2033  CHKERR moab.get_adjacencies(ents_nodes_edges_nodes, 1, false,
2034  ents_nodes_edges_nodes_edges,
2035  moab::Interface::UNION);
2036  // remove edges adj. to hanging edges
2037  ents_nodes_edges =
2038  subtract(ents_nodes_edges, ents_nodes_edges_nodes_edges);
2039  ents_nodes_edges =
2040  subtract(ents_nodes_edges, ents.subset_by_type(MBEDGE));
2041  if(length>0) {
2042  Range::iterator eit = ents_nodes_edges.begin();
2043  for (; eit != ents_nodes_edges.end();) {
2044 
2045  int num_nodes;
2046  const EntityHandle *conn;
2047  rval = moab.get_connectivity(*eit, conn, num_nodes, true);
2048  double coords[6];
2049  if(tH) {
2050  CHKERR moab.tag_get_data(tH, conn, num_nodes, coords);
2051  } else {
2052  CHKERR moab.get_coords(conn, num_nodes, coords);
2053  }
2054 
2055  auto get_edge_length = [coords]() {
2057  &coords[0], &coords[1], &coords[2]);
2060  t_delta(i) = t_coords(i);
2061  ++t_coords;
2062  t_delta(i) -= t_coords(i);
2063  return sqrt(t_delta(i) * t_delta(i));
2064  };
2065 
2066  if (get_edge_length() < tOL) {
2067  eit = ents_nodes_edges.erase(eit);
2068  } else {
2069  eit++;
2070  }
2071  }
2072  }
2073  edges_to_remove.swap(ents_nodes_edges);
2074  if (debug) {
2075  CHKERR SaveData(moab)("edges_to_remove.vtk", edges_to_remove);
2076  }
2078  }
2079  };
2080 
2081  Range not_merged_edges;
2082  const double tol = 0.05;
2083  CHKERR Toplogy(m_field, th, tol * aveLength)
2084  .edgesToMerge(surface, tets, edges_to_merge);
2085  CHKERR Toplogy(m_field, th, tol * aveLength)
2086  .removeBadEdges(surface, tets, fixed_edges, corner_nodes, edges_to_merge,
2087  not_merged_edges);
2088  Toplogy::SetsMap sets_map;
2089  CHKERR Toplogy(m_field, th, tol * aveLength)
2090  .classifyVerts(surface, tets, fixed_edges, corner_nodes, sets_map);
2091  Range proc_tets;
2092  CHKERR Toplogy(m_field, th, tol * aveLength)
2093  .getProcTets(tets, edges_to_merge, proc_tets);
2094  out_tets = subtract(tets, proc_tets);
2095  if (bit_ptr) {
2096  for (int dd = 2; dd >= 0; dd--) {
2097  CHKERR moab.get_adjacencies(out_tets.subset_by_dimension(3), dd, false,
2098  out_tets, moab::Interface::UNION);
2099  }
2100  CHKERR m_field.getInterface<BitRefManager>()->addBitRefLevel(out_tets,
2101  *bit_ptr);
2102  }
2103 
2104  int nb_nodes_merged = 0;
2105  LengthMapData_multi_index length_map;
2106  new_surf = surface;
2107 
2108  double ave0 = 0, ave = 0, min = 0, min_p = 0, min_pp;
2109  for (int pp = 0; pp != nbMaxMergingCycles; ++pp) {
2110 
2111  int nb_nodes_merged_p = nb_nodes_merged;
2112  length_map.clear();
2113  min_pp = min_p;
2114  min_p = min;
2115  CHKERR LengthMap(m_field, th, aveLength)(proc_tets, edges_to_merge,
2116  length_map, ave);
2117  min = length_map.get<0>().begin()->qUality;
2118  if(pp == 0) {
2119  ave0 = ave;
2120  }
2121 
2122  int nn = 0;
2123  Range collapsed_edges;
2124  for (LengthMapData_multi_index::nth_index<0>::type::iterator
2125  mit = length_map.get<0>().begin();
2126  mit != length_map.get<0>().end(); mit++, nn++) {
2127  if(mit->skip) continue;
2128  int num_nodes;
2129  const EntityHandle *conn;
2130  CHKERR moab.get_connectivity(mit->eDge, conn, num_nodes, true);
2131  int conn_type[2] = {0, 0};
2132  for (int nn = 0; nn != 2; nn++) {
2133  conn_type[nn] = 0;
2134  for (Toplogy::SetsMap::reverse_iterator sit = sets_map.rbegin();
2135  sit != sets_map.rend(); sit++) {
2136  if (sit->second.find(conn[nn]) != sit->second.end()) {
2137  conn_type[nn] |= sit->first;
2138  }
2139  }
2140  }
2141  int type_father, type_mother;
2142  EntityHandle father, mother;
2143  if (conn_type[0] > conn_type[1]) {
2144  father = conn[0];
2145  mother = conn[1];
2146  type_father = conn_type[0];
2147  type_mother = conn_type[1];
2148  } else {
2149  father = conn[1];
2150  mother = conn[0];
2151  type_father = conn_type[1];
2152  type_mother = conn_type[0];
2153  }
2154  int line_search = 0;
2155  if (type_father == type_mother) {
2156  line_search = lineSearchSteps;
2157  }
2158 
2159  CHKERR MergeNodes(m_field, true, line_search, th,
2160  update_meshsets)(father, mother, proc_tets, new_surf,
2161  edges_to_merge, not_merged_edges);
2162 
2163  if (m_field.getInterface<NodeMergerInterface>()->getSucessMerge()) {
2164  Range adj_mother_tets;
2165  CHKERR moab.get_adjacencies(&mother, 1, 3, false, adj_mother_tets);
2166  Range adj_mother_tets_nodes;
2167  CHKERR moab.get_connectivity(adj_mother_tets, adj_mother_tets_nodes,
2168  true);
2169  Range adj_edges;
2170  CHKERR moab.get_adjacencies(adj_mother_tets_nodes, 1, false, adj_edges,
2171  moab::Interface::UNION);
2172  CHKERR moab.get_adjacencies(&father, 1, 1, false, adj_edges,
2173  moab::Interface::UNION);
2174  for (Range::iterator ait = adj_edges.begin(); ait != adj_edges.end();
2175  ait++) {
2176  LengthMapData_multi_index::nth_index<1>::type::iterator miit =
2177  length_map.get<1>().find(*ait);
2178  if (miit != length_map.get<1>().end()) {
2179  (const_cast<LengthMapData &>(*miit)).skip = true;
2180  }
2181  }
2182  nb_nodes_merged++;
2183  collapsed_edges.insert(mit->eDge);
2184  }
2185 
2186  if (nn > static_cast<int>(length_map.size() / fraction_level))
2187  break;
2188  if (mit->qUality > ave)
2189  break;
2190  }
2191 
2192  Range adj_faces, adj_edges;
2193  CHKERR moab.get_adjacencies(proc_tets, 2, false, adj_faces,
2194  moab::Interface::UNION);
2195  new_surf = intersect(new_surf, adj_faces);
2196 
2197  PetscPrintf(m_field.get_comm(),
2198  "(%d) Number of nodes merged %d ave q %3.4e min q %3.4e\n", pp,
2199  nb_nodes_merged, ave, min);
2200 
2201  if (debug) {
2202  // current surface and merged edges in step
2203  std::string name;
2204  name = "node_merge_step_" + boost::lexical_cast<std::string>(pp) + ".vtk";
2205  CHKERR SaveData(moab)(name, unite(new_surf, collapsed_edges));
2206  name = "edges_to_merge_step_" + boost::lexical_cast<std::string>(pp) +
2207  ".vtk";
2208  CHKERR SaveData(moab)(name, unite(new_surf, edges_to_merge));
2209  }
2210 
2211  if (nb_nodes_merged == nb_nodes_merged_p)
2212  break;
2213  if (min > 1e-2 && min == min_pp)
2214  break;
2215  if (min > ave0)
2216  break;
2217 
2218  CHKERR moab.get_adjacencies(proc_tets, 1, false, adj_edges,
2219  moab::Interface::UNION);
2220  edges_to_merge = intersect(edges_to_merge, adj_edges);
2221  CHKERR Toplogy(m_field, th, tol * aveLength)
2222  .removeBadEdges(new_surf, proc_tets, fixed_edges, corner_nodes,
2223  edges_to_merge, not_merged_edges);
2224  }
2225 
2226  if (bit_ptr) {
2227  CHKERR m_field.getInterface<BitRefManager>()->setBitRefLevel(proc_tets,
2228  *bit_ptr);
2229  }
2230  out_tets.merge(proc_tets);
2231 
2233 }
2234 
2236  const int fraction_level, const BitRefLevel trim_bit,
2237  const BitRefLevel cut_bit, const BitRefLevel bit, const Range &surface,
2238  const Range &fixed_edges, const Range &corner_nodes, Tag th,
2239  const bool update_meshsets, const bool debug) {
2240  CoreInterface &m_field = cOre;
2242  Range tets_level;
2243  ierr = m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2244  trim_bit, BitRefLevel().set(), MBTET, tets_level);
2245  CHKERRG(ierr);
2246 
2247  Range edges_to_merge;
2248  ierr = m_field.getInterface<BitRefManager>()->getEntitiesByParentType(
2249  trim_bit, trim_bit | cut_bit, MBEDGE, edges_to_merge);
2250  CHKERRG(ierr);
2251  edges_to_merge = edges_to_merge.subset_by_type(MBEDGE);
2252 
2253  // get all entities not in database
2254  Range all_ents_not_in_database_before;
2255  ierr = cOre.getInterface<BitRefManager>()->getAllEntitiesNotInDatabase(
2256  all_ents_not_in_database_before);
2257  CHKERRG(ierr);
2258 
2259  Range out_new_tets, out_new_surf;
2260  ierr = mergeBadEdges(fraction_level, tets_level, surface, fixed_edges,
2261  corner_nodes, edges_to_merge, out_new_tets, out_new_surf,
2262  th, update_meshsets, &bit, debug);
2263  CHKERRG(ierr);
2264 
2265  // get all entities not in database after merge
2266  Range all_ents_not_in_database_after;
2267  ierr = cOre.getInterface<BitRefManager>()->getAllEntitiesNotInDatabase(
2268  all_ents_not_in_database_after);
2269  CHKERRG(ierr);
2270  // delete hanging entities
2271  all_ents_not_in_database_after =
2272  subtract(all_ents_not_in_database_after, all_ents_not_in_database_before);
2273  Range meshsets;
2274  CHKERR m_field.get_moab().get_entities_by_type(0, MBENTITYSET, meshsets,
2275  false);
2276  for (Range::iterator mit = meshsets.begin(); mit != meshsets.end(); mit++) {
2277  CHKERR m_field.get_moab().remove_entities(*mit,
2278  all_ents_not_in_database_after);
2279  }
2280  m_field.get_moab().delete_entities(all_ents_not_in_database_after);
2281 
2282  mergedVolumes.swap(out_new_tets);
2283  mergedSurfaces.swap(out_new_surf);
2285 }
2286 
2287 #ifdef WITH_TETGEN
2288 
2290  vector<string> &switches, const BitRefLevel &mesh_bit,
2291  const BitRefLevel &bit, const Range &surface, const Range &fixed_edges,
2292  const Range &corner_nodes, Tag th, const bool debug) {
2293  CoreInterface &m_field = cOre;
2294  moab::Interface &moab = m_field.get_moab();
2295  TetGenInterface *tetgen_iface;
2297  CHKERR m_field.getInterface(tetgen_iface);
2298 
2299  tetGenData.clear();
2300  moabTetGenMap.clear();
2301  tetGenMoabMap.clear();
2302 
2303  if (tetGenData.size() < 1) {
2304  tetGenData.push_back(new tetgenio);
2305  }
2306  tetgenio &in = tetGenData.back();
2307 
2308  struct BitEnts {
2309 
2310  CoreInterface &mField;
2311  const BitRefLevel &bIt;
2312  BitEnts(CoreInterface &m_field, const BitRefLevel &bit)
2313  : mField(m_field), bIt(bit) {}
2314 
2315  Range mTets;
2316  Range mTris;
2317  Range mEdges;
2318  Range mNodes;
2319 
2320  MoFEMErrorCode getLevelEnts() {
2322  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2323  bIt, BitRefLevel().set(), MBTET, mTets);
2324  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2325  bIt, BitRefLevel().set(), MBTRI, mTris);
2326  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2327  bIt, BitRefLevel().set(), MBEDGE, mEdges);
2328  CHKERR mField.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2329  bIt, BitRefLevel().set(), MBVERTEX, mNodes);
2331  }
2332 
2333  Range mSkin;
2334  Range mSkinNodes;
2335  Range mSkinEdges;
2336 
2337  MoFEMErrorCode getSkin() {
2338  moab::Interface &moab = mField.get_moab();
2340  Skinner skin(&moab);
2341  CHKERR skin.find_skin(0, mTets, false, mSkin);
2342  CHKERR mField.get_moab().get_connectivity(mSkin, mSkinNodes, true);
2343  CHKERR mField.get_moab().get_adjacencies(mSkin, 1, false, mSkinEdges,
2344  moab::Interface::UNION);
2346  }
2347  };
2348 
2349  struct SurfaceEnts {
2350 
2351  CoreInterface &mField;
2352  SurfaceEnts(CoreInterface &m_field) : mField(m_field) {}
2353 
2354  Range sNodes;
2355  Range sEdges;
2356  Range sVols;
2357  Range vNodes;
2358 
2359  MoFEMErrorCode getVolume(const BitEnts &bit_ents, const Range &tris) {
2360  moab::Interface &moab = mField.get_moab();
2362  CHKERR moab.get_connectivity(tris, sNodes, true);
2363  CHKERR moab.get_adjacencies(tris, 1, false, sEdges,
2364  moab::Interface::UNION);
2365  CHKERR moab.get_adjacencies(sNodes, 3, false, sVols,
2366  moab::Interface::UNION);
2367  sVols = intersect(sVols, bit_ents.mTets);
2368  CHKERR moab.get_connectivity(sVols, vNodes, true);
2370  }
2371 
2372  Range sSkin;
2373  Range sSkinNodes;
2374  Range vSkin;
2375  Range vSkinNodes;
2376  Range vSkinWithoutBodySkin;
2377  Range vSkinNodesWithoutBodySkin;
2378  Range vSkinOnBodySkin;
2379  Range vSkinOnBodySkinNodes;
2380 
2381  MoFEMErrorCode getSkin(const BitEnts &bit_ents, const Range &tris,
2382  const int levels = 3) {
2383  moab::Interface &moab = mField.get_moab();
2385  Skinner skin(&moab);
2386  rval = skin.find_skin(0, sVols, false, vSkin);
2387  for (int ll = 0; ll != levels; ll++) {
2388  CHKERR moab.get_adjacencies(vSkin, 3, false, sVols,
2389  moab::Interface::UNION);
2390  sVols = intersect(sVols, bit_ents.mTets);
2391  vSkin.clear();
2392  CHKERR skin.find_skin(0, sVols, false, vSkin);
2393  }
2394  vSkinWithoutBodySkin = subtract(vSkin, bit_ents.mSkin);
2395  vSkinOnBodySkin = intersect(vSkin, bit_ents.mSkin);
2396  CHKERR moab.get_connectivity(vSkinOnBodySkin, vSkinOnBodySkinNodes, true);
2397  CHKERR moab.get_connectivity(sVols, vNodes, true);
2398  CHKERR moab.get_connectivity(vSkin, vSkinNodes, true);
2399  vSkinNodesWithoutBodySkin = subtract(vSkinNodes, bit_ents.mSkinNodes);
2400  CHKERR skin.find_skin(0, tris, false, sSkin);
2401  CHKERR moab.get_connectivity(sSkin, sSkinNodes, true);
2402  tVols = sVols;
2404  }
2405 
2406  Range tVols;
2407 
2408  MoFEMErrorCode getTetsForRemesh(const BitEnts &bit_ents, Tag th = NULL) {
2409  moab::Interface &moab = mField.get_moab();
2411 
2412  Range tets_with_four_nodes_on_skin;
2413  rval = moab.get_adjacencies(vSkinOnBodySkinNodes, 3, false,
2414  tets_with_four_nodes_on_skin,
2415  moab::Interface::UNION);
2416  Range tets_nodes;
2417  CHKERR moab.get_connectivity(tets_with_four_nodes_on_skin, tets_nodes,
2418  true);
2419  tets_nodes = subtract(tets_nodes, vSkinOnBodySkinNodes);
2420  Range other_tets;
2421  CHKERR moab.get_adjacencies(tets_nodes, 3, false, other_tets,
2422  moab::Interface::UNION);
2423  tets_with_four_nodes_on_skin =
2424  subtract(tets_with_four_nodes_on_skin, other_tets);
2425  Range to_remove;
2426  for (Range::iterator tit = tets_with_four_nodes_on_skin.begin();
2427  tit != tets_with_four_nodes_on_skin.end(); tit++) {
2428  int num_nodes;
2429  const EntityHandle *conn;
2430  CHKERR moab.get_connectivity(*tit, conn, num_nodes, true);
2431  double coords[12];
2432  if (th) {
2433  CHKERR moab.tag_get_data(th, conn, num_nodes, coords);
2434  } else {
2435  CHKERR moab.get_coords(conn, num_nodes, coords);
2436  }
2437  double quality = Tools::volumeLengthQuality(coords);
2438  if (quality < 1e-2) {
2439  to_remove.insert(*tit);
2440  }
2441  }
2442 
2443  sVols = subtract(sVols, to_remove);
2444 
2445  Skinner skin(&moab);
2446  vSkin.clear();
2447  CHKERR skin.find_skin(0, sVols, false, vSkin);
2448  Range m_skin;
2449  CHKERR
2450  skin.find_skin(0, subtract(bit_ents.mSkin, to_remove), false, m_skin);
2451  vSkinWithoutBodySkin = subtract(vSkin, m_skin);
2452  vSkinOnBodySkin = intersect(vSkin, m_skin);
2453  vNodes.clear();
2454  vSkinNodes.clear();
2455  vSkinOnBodySkinNodes.clear();
2456  CHKERR moab.get_connectivity(sVols, vNodes, true);
2457  CHKERR moab.get_connectivity(vSkinOnBodySkin, vSkinOnBodySkinNodes, true);
2458  CHKERR moab.get_connectivity(vSkin, vSkinNodes, true);
2460  }
2461  };
2462 
2463  BitEnts bit_ents(m_field, mesh_bit);
2464  CHKERR bit_ents.getLevelEnts();
2465  CHKERR bit_ents.getSkin();
2466  SurfaceEnts surf_ents(m_field);
2467  CHKERR surf_ents.getVolume(bit_ents, surface);
2468  CHKERR surf_ents.getSkin(bit_ents, surface);
2469  CHKERR surf_ents.getTetsForRemesh(bit_ents);
2470 
2471  map<int, Range> types_ents;
2472  types_ents[TetGenInterface::RIDGEVERTEX].merge(
2473  surf_ents.vSkinNodesWithoutBodySkin);
2474  // FREESEGVERTEX
2475  types_ents[TetGenInterface::FREESEGVERTEX].merge(surf_ents.sSkinNodes);
2476  types_ents[TetGenInterface::FREESEGVERTEX] =
2477  subtract(types_ents[TetGenInterface::FREESEGVERTEX],
2478  types_ents[TetGenInterface::RIDGEVERTEX]);
2479  // FREEFACETVERTEX
2480  types_ents[TetGenInterface::FREEFACETVERTEX].merge(surf_ents.sNodes);
2481  types_ents[TetGenInterface::FREEFACETVERTEX] =
2482  subtract(types_ents[TetGenInterface::FREEFACETVERTEX],
2483  types_ents[TetGenInterface::RIDGEVERTEX]);
2484  types_ents[TetGenInterface::FREEFACETVERTEX] =
2485  subtract(types_ents[TetGenInterface::FREEFACETVERTEX],
2486  types_ents[TetGenInterface::FREESEGVERTEX]);
2487  // FREEVOLVERTEX
2488  types_ents[TetGenInterface::FREEVOLVERTEX].merge(surf_ents.vNodes);
2489  types_ents[TetGenInterface::FREEVOLVERTEX] =
2490  subtract(types_ents[TetGenInterface::FREEVOLVERTEX],
2491  types_ents[TetGenInterface::RIDGEVERTEX]);
2492  types_ents[TetGenInterface::FREEVOLVERTEX] =
2493  subtract(types_ents[TetGenInterface::FREEVOLVERTEX],
2494  types_ents[TetGenInterface::FREESEGVERTEX]);
2495  types_ents[TetGenInterface::FREEVOLVERTEX] =
2496  subtract(types_ents[TetGenInterface::FREEVOLVERTEX],
2497  types_ents[TetGenInterface::FREEFACETVERTEX]);
2498 
2499  Tag th_marker;
2500  // Clean markers
2501  rval = m_field.get_moab().tag_get_handle("TETGEN_MARKER", th_marker);
2502  if(rval == MB_SUCCESS) {
2503  CHKERR m_field.get_moab().tag_delete(th_marker);
2504  rval = MB_SUCCESS;
2505  }
2506 
2507  int def_marker = 0;
2508  CHKERR m_field.get_moab().tag_get_handle(
2509  "TETGEN_MARKER", 1, MB_TYPE_INTEGER, th_marker,
2510  MB_TAG_CREAT | MB_TAG_SPARSE, &def_marker);
2511 
2512  // Mark surface with id = 1
2513  vector<int> markers(surface.size(), 1);
2514  CHKERR m_field.get_moab().tag_set_data(th_marker, surface, &*markers.begin());
2515  // Mark all side sets
2516  int shift = 1;
2517  map<int, int> id_shift_map; // each meshset has set unique bit
2519  (*cOre.getInterface<MeshsetsManager>()), SIDESET, it)) {
2520  int ms_id = it->getMeshsetId();
2521  id_shift_map[ms_id] = 1 << shift; // shift bit
2522  ++shift;
2523  Range sideset_faces;
2524  CHKERR m_field.getInterface<MeshsetsManager>()->getEntitiesByDimension(
2525  ms_id, SIDESET, 2, sideset_faces, true);
2526  markers.resize(sideset_faces.size());
2527  CHKERR m_field.get_moab().tag_get_data(th_marker, sideset_faces,
2528  &*markers.begin());
2529  for (unsigned int ii = 0; ii < markers.size(); ii++) {
2530  markers[ii] |= id_shift_map[ms_id]; // add bit to marker
2531  }
2532  CHKERR m_field.get_moab().tag_set_data(th_marker, sideset_faces,
2533  &*markers.begin());
2534  }
2535  Range nodes_to_remove; // none
2536  markers.resize(nodes_to_remove.size());
2537  fill(markers.begin(), markers.end(), -1);
2538  CHKERR m_field.get_moab().tag_set_data(th_marker, nodes_to_remove,
2539  &*markers.begin());
2540 
2541  // nodes
2542  if (tetGenData.size() == 1) {
2543 
2544  Range ents_to_tetgen = surf_ents.sVols;
2545  CHKERR m_field.get_moab().get_connectivity(surf_ents.sVols, ents_to_tetgen,
2546  true);
2547  for (int dd = 2; dd >= 1; dd--) {
2548  CHKERR m_field.get_moab().get_adjacencies(
2549  surf_ents.sVols, dd, false, ents_to_tetgen, moab::Interface::UNION);
2550  }
2551 
2552  // Load mesh to TetGen data structures
2553  CHKERR tetgen_iface->inData(ents_to_tetgen, in, moabTetGenMap,
2554  tetGenMoabMap, th);
2555  CHKERR tetgen_iface->setGeomData(in, moabTetGenMap, tetGenMoabMap,
2556  types_ents);
2557  std::vector<pair<Range, int> > markers;
2558  for (Range::iterator tit = surface.begin(); tit != surface.end(); tit++) {
2559  Range facet;
2560  facet.insert(*tit);
2561  markers.push_back(pair<Range, int>(facet, 2));
2562  }
2563  for (Range::iterator tit = surf_ents.vSkinWithoutBodySkin.begin();
2564  tit != surf_ents.vSkinWithoutBodySkin.end(); tit++) {
2565  Range facet;
2566  facet.insert(*tit);
2567  markers.push_back(pair<Range, int>(facet, 1));
2568  }
2569  Range other_facets;
2570  other_facets = subtract(surf_ents.vSkin, surf_ents.vSkinWithoutBodySkin);
2571  for (Range::iterator tit = other_facets.begin(); tit != other_facets.end();
2572  tit++) {
2573  Range facet;
2574  facet.insert(*tit);
2575  markers.push_back(pair<Range, int>(facet, 0));
2576  }
2577  CHKERR tetgen_iface->setFaceData(markers, in, moabTetGenMap, tetGenMoabMap);
2578  }
2579 
2580  if (debug) {
2581  if (m_field.get_comm_rank() == 0) {
2582  char tetgen_in_file_name[] = "in";
2583  in.save_nodes(tetgen_in_file_name);
2584  in.save_elements(tetgen_in_file_name);
2585  in.save_faces(tetgen_in_file_name);
2586  in.save_edges(tetgen_in_file_name);
2587  in.save_poly(tetgen_in_file_name);
2588  }
2589  }
2590 
2591  // generate new mesh
2592  {
2593  vector<string>::iterator sw = switches.begin();
2594  for (int ii = 0; sw != switches.end(); sw++, ii++) {
2595  tetgenio &_in_ = tetGenData.back();
2596  tetGenData.push_back(new tetgenio);
2597  tetgenio &_out_ = tetGenData.back();
2598  char *s = const_cast<char *>(sw->c_str());
2599  CHKERR tetgen_iface->tetRahedralize(s, _in_, _out_);
2600  }
2601  }
2602  tetgenio &out = tetGenData.back();
2603  // save elems
2604  if (debug) {
2605  char tetgen_out_file_name[] = "out";
2606  out.save_nodes(tetgen_out_file_name);
2607  out.save_elements(tetgen_out_file_name);
2608  out.save_faces(tetgen_out_file_name);
2609  out.save_edges(tetgen_out_file_name);
2610  out.save_poly(tetgen_out_file_name);
2611  }
2612 
2613  CHKERR tetgen_iface->outData(in, out, moabTetGenMap, tetGenMoabMap, bit,
2614  false, false);
2615 
2616  Range rest_of_ents = subtract(bit_ents.mTets, surf_ents.tVols);
2617  for (int dd = 2; dd >= 0; dd--) {
2618  CHKERR moab.get_adjacencies(rest_of_ents.subset_by_dimension(3), dd, false,
2619  rest_of_ents, moab::Interface::UNION);
2620  }
2621  CHKERR m_field.getInterface<BitRefManager>()->addBitRefLevel(rest_of_ents,
2622  bit);
2623 
2624  Range tetgen_faces;
2625  map<int, Range> face_markers_map;
2626  CHKERR tetgen_iface->getTriangleMarkers(tetGenMoabMap, out, &tetgen_faces,
2627  &face_markers_map);
2628 
2629  tetgenSurfaces = face_markers_map[1];
2630  for (map<int, Range>::iterator mit = face_markers_map.begin();
2631  mit != face_markers_map.end(); mit++) {
2633  (*cOre.getInterface<MeshsetsManager>()), SIDESET, it)) {
2634  int msId = it->getMeshsetId();
2635  if (id_shift_map[msId] & mit->first) {
2636  EntityHandle meshset = it->getMeshset();
2637  CHKERR m_field.get_moab().add_entities(
2638  meshset, mit->second.subset_by_type(MBTRI));
2639  }
2640  }
2641  }
2642 
2644 }
2645 
2646 #endif // WITH_TETGEN
2647 
2649  CoreInterface &m_field = cOre;
2650  moab::Interface &moab = m_field.get_moab();
2652  Range nodes;
2653  if(bit.none()) {
2654  CHKERR moab.get_entities_by_type(0, MBVERTEX, nodes);
2655  } else {
2656  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2657  bit,BitRefLevel().set(),MBVERTEX,nodes
2658  );
2659  }
2660  std::vector<double> coords(3 * nodes.size());
2661  CHKERR moab.get_coords(nodes, &coords[0]);
2662  CHKERR moab.tag_set_data(th, nodes, &coords[0]);
2664 }
2665 
2667  const BitRefLevel mask) {
2668  CoreInterface &m_field = cOre;
2669  moab::Interface &moab = m_field.get_moab();
2671  Range nodes;
2672  if(bit.none()) {
2673  CHKERR moab.get_entities_by_type(0, MBVERTEX, nodes);
2674  } else {
2675  CHKERR m_field.getInterface<BitRefManager>()->getEntitiesByTypeAndRefLevel(
2676  bit, mask, MBVERTEX, nodes);
2677  }
2678  std::vector<double> coords(3 * nodes.size());
2679  CHKERR moab.tag_get_data(th, nodes, &coords[0]);
2680  CHKERR moab.set_coords(nodes, &coords[0]);
2682 }
2683 
2685  CoreInterface &m_field = cOre;
2686  moab::Interface &moab = m_field.get_moab();
2688  CHKERR SaveData(moab)("out_vol.vtk", vOlume);
2689  CHKERR SaveData(moab)("out_surface.vtk", sUrface);
2690  CHKERR SaveData(moab)("out_cut_edges.vtk", cutEdges);
2691  CHKERR SaveData(moab)("out_cut_volumes.vtk", cutVolumes);
2692  CHKERR SaveData(moab)("out_cut_new_volumes.vtk", cutNewVolumes);
2693  CHKERR SaveData(moab)("out_cut_new_surfaces.vtk", cutNewSurfaces);
2694  CHKERR SaveData(moab)("out_cut_zero_distance_ents.vtk", zeroDistanceEnts);
2696 }
2697 
2699  moab::Interface &moab = cOre.getInterface<CoreInterface>()->get_moab();
2701  CHKERR SaveData(moab)("out_trim_new_volumes.vtk", trimNewVolumes);
2702  CHKERR SaveData(moab)("out_trim_new_surfaces.vtk", trimNewSurfaces);
2703  CHKERR SaveData(moab)("out_trim_edges.vtk", trimEdges);
2705 }
2706 
2707 } // namespace MoFEM
MoFEMErrorCode rebuildMeshWithTetGen(vector< string > &switches, const BitRefLevel &mesh_bit, const BitRefLevel &bit, const Range &surface, const Range &fixed_edges, const Range &corner_nodes, Tag th=NULL, const bool debug=false)
MoFEMErrorCode buildTree()
build tree
map< EntityHandle, TreeData > edgesToTrim
MoFEMErrorCode tetRahedralize(char switches[], tetgenio &in, tetgenio &out)
run tetgen
virtual MoFEMErrorCode get_ref_ents(const RefEntity_multiIndex **refined_ents_ptr) const =0
Get ref entities multi-index from database.
MoFEM interface unique ID.
Interface to cut meshes.
MoFEMErrorCode findEdgesToCut(Range *fixed_edges, Range *corner_nodes, const double low_tol=0, int verb=0, const bool debug=false)
find edges to cut
MoFEMErrorCode setTagData(Tag th, const BitRefLevel bit=BitRefLevel())
set coords to tag
MoFEMErrorCode getTriangleMarkers(tetGenMoab_Map &tetgen_moab_map, tetgenio &out, Range *ents=NULL, idxRange_Map *ents_map=NULL, bool only_non_zero=true)
get markers to faces
LengthMapData(const double l, double q, const EntityHandle e)
MoFEMErrorCode findIfTringleHasThreeNodesOnInternalSurfaceSkin(const EntityHandle sideset, const BitRefLevel mesh_bit_level, const bool recursive, Range &faces_with_three_nodes_on_front, int verb=QUIET)
Find if triangle has three nodes on internal surface skin.
static double volumeLengthQuality(const double *coords)
Calculate tetrahedron volume length quality.
Definition: Tools.cpp:32
VectorBoundedArray< double, 6 > VectorDouble6
Definition: Common.hpp:221
PetscErrorCode MoFEMErrorCode
MoFEM/PETSc error code.
Definition: Common.hpp:60
virtual moab::Interface & get_moab()=0
const Tensor2_symmetric_Expr< const ddTensor0< T, Dim, i, j >, typename promote< T, double >::V, Dim, i, j > dd(const Tensor0< T *> &a, const Index< i, Dim > index1, const Index< j, Dim > index2, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: ddTensor0.hpp:33
SaveData(moab::Interface &moab)
double maxLength
Maximal edge length.
MoFEMErrorCode copySurface(const Range &surface, Tag th=NULL, double *shift=NULL, double *origin=NULL, double *transform=NULL, const std::string save_mesh="")
copy surface entities
use TetGen to generate mesh
MoFEMErrorCode getRayForEdge(const EntityHandle ent, VectorAdaptor &ray_point, VectorAdaptor &unit_ray_dir, double &ray_length) const
#define MoFEMFunctionBeginHot
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:499
MoFEMErrorCode setSurface(const Range &surface)
set surface entities
VectorBoundedArray< double, 9 > VectorDouble9
Definition: Common.hpp:222
map< EntityHandle, TreeData > verticesOnTrimEdges
static double get_ave_edge_length(const EntityHandle ent, const Range &vol_edges, moab::Interface &moab)
const Tensor1_Expr< const dTensor0< T, Dim, i >, typename promote< T, double >::V, Dim, i > d(const Tensor0< T *> &a, const Index< i, Dim > index, const Tensor1< int, Dim > &d_ijk, const Tensor1< double, Dim > &d_xyz)
Definition: dTensor0.hpp:27
base class for all interface classes
boost::shared_ptr< OrientedBoxTreeTool > treeSurfPtr
#define MoFEMFunctionReturn(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:475
map< EntityHandle, unsigned long > moabTetGenMap
MoFEMErrorCode mergeBadEdges(const int fraction_level, const Range &tets, const Range &surface, const Range &fixed_edges, const Range &corner_nodes, Range &merged_nodes, Range &out_tets, Range &new_surf, Tag th, const bool update_meshsets=false, const BitRefLevel *bit_ptr=NULL, const bool debug=false)
Merge edges.
MoFEMErrorCode splitSides(const EntityHandle meshset, const BitRefLevel &bit, const int msId, const CubitBCType cubit_bc_type, const bool add_interface_entities, const bool recursive=false, int verb=QUIET)
split nodes and other entities of tetrahedral in children sets and add prism elements ...
#define CHKERRG(n)
Check error code of MoFEM/MOAB/PETSc function.
Definition: definitions.h:542
double aveLength
Average edge length.
Auxiliary tools.
Definition: Tools.hpp:30
Interface for managing meshsets containing materials and boundary conditions.
Core (interface) class.
Definition: Core.hpp:50
#define MoFEMFunctionReturnHot(a)
Last executable line of each PETSc function used for error handling. Replaces return() ...
Definition: definitions.h:506
MoFEMErrorCode cutTrimAndMerge(const int fraction_level, const BitRefLevel &bit_level1, const BitRefLevel &bit_level2, const BitRefLevel &bit_level3, Tag th, const double tol_cut, const double tol_cut_close, const double tol_trim, const double tol_trim_close, Range &fixed_edges, Range &corner_nodes, const bool update_meshsets=false, const bool debug=false)
MoFEMErrorCode removePathologicalFrontTris(const BitRefLevel split_bit, Range &ents)
Remove pathological elements on surface internal front.
implementation of Data Operators for Forces and Sources
Definition: Common.hpp:21
MoFEMErrorCode setFaceData(std::vector< std::pair< Range, int > > &markers, tetgenio &in, moabTetGen_Map &moab_tetgen_map, tetGenMoab_Map &tetgen_moab_map)
set markers to faces
VectorBoundedArray< double, 12 > VectorDouble12
Definition: Common.hpp:223
VectorShallowArrayAdaptor< double > VectorAdaptor
Definition: Common.hpp:234
MoFEMErrorCode cutAndTrim(const BitRefLevel &bit_level1, const BitRefLevel &bit_level2, Tag th, const double tol_cut, const double tol_cut_close, const double tol_trim, const double tol_trim_close, Range *fixed_edges=NULL, Range *corner_nodes=NULL, const bool update_meshsets=false, const bool debug=true)
MatrixShallowArrayAdaptor< double > MatrixAdaptor
Definition: Common.hpp:266
CutMeshInterface(const MoFEM::Core &core)
std::bitset< BITREFLEVEL_SIZE > BitRefLevel
Bit structure attached to each entity identifying to what mesh entity is attached.
Definition: Common.hpp:147
static MoFEMErrorCodeGeneric< PetscErrorCode > ierr
Definition: Common.hpp:80
moab::Interface & moab
const Range & projectZeroDistanceEnts() const
MoFEMErrorCode getInterface(const MOFEMuuid &uuid, IFACE *&iface) const
Get interface by uuid and return reference to pointer of interface.
virtual MoFEMErrorCode refine_TET(const EntityHandle meshset, const BitRefLevel &bit, const bool respect_interface=false, int verb=0, Range *ref_edges=NULL)
refine TET in the meshset
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 from two bit levels
Definition: NodeMerger.hpp:29
virtual int get_comm_rank() const =0
const Range & getNewTrimSurfaces() const
double tol
static MoFEMErrorCodeGeneric< moab::ErrorCode > rval
Definition: Common.hpp:78
MoFEMErrorCode operator()(Core &core, const BitRefLevel &bit) const
ublas::matrix< double, ublas::row_major, DoubleAllocator > MatrixDouble
Definition: Common.hpp:212
static const bool debug
InterfaceThis interface is used by user to:
Definition: Interface.hpp:42
const Range & getVolume() const
boost::ptr_vector< tetgenio > tetGenData
MoFEMErrorCode setCoords(Tag th, const BitRefLevel bit=BitRefLevel(), const BitRefLevel mask=BitRefLevel().set())
set coords from tag
Managing BitRefLevels.
MoFEMErrorCode mergeSurface(const Range &surface)
merge surface entities
MoFEMErrorCode moveMidNodesOnCutEdges(Tag th=NULL)
projecting of mid edge nodes on new mesh on surface
multi_index_container< LengthMapData, indexed_by< ordered_non_unique< member< LengthMapData, double, &LengthMapData::lEngth > >, hashed_unique< member< LengthMapData, EntityHandle, &LengthMapData::eDge > > >> LengthMapData_multi_index
MoFEMErrorCode mergeVolumes(const Range &volume)
merge volume entities
MoFEMErrorCode splitSides(const BitRefLevel split_bit, const BitRefLevel bit, const Range &ents, Tag th=NULL)
split sides
MoFEMErrorCode query_interface(const MOFEMuuid &uuid, UnknownInterface **iface) const
Make interface on given faces and create flat prism in that space.
ParentChildMap & getParentChildMap()
Get map of parent cand child.
Definition: NodeMerger.hpp:156
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:66
Mesh refinement interface.
MoFEMErrorCode inData(Range &ents, tetgenio &in, moabTetGen_Map &moab_tetgen_map, tetGenMoab_Map &tetgen_moab_map, Tag th=NULL)
create TetGen data structure form range of moab entities
MoFEMErrorCode getSides(const int msId, const CubitBCType cubit_bc_type, const BitRefLevel mesh_bit_level, const bool recursive, int verb=QUIET)
create two children meshsets in the meshset containing tetrahedral on two sides of faces ...
MoFEMErrorCode cutEdgesInMiddle(const BitRefLevel bit, const bool debug=false)
cut edges
#define _IT_CUBITMESHSETS_BY_SET_TYPE_FOR_LOOP_(MESHSET_MANAGER, CUBITBCTYPE, IT)
Iterator that loops over a specific Cubit MeshSet having a particular BC meshset in a moFEM field...
#define CHKERR
Inline error check.
Definition: definitions.h:594
static const MOFEMuuid IDD_MOFEMCutMesh
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 outData(tetgenio &in, tetgenio &out, moabTetGen_Map &moab_tetgen_map, tetGenMoab_Map &tetgen_moab_map, Range *ents=NULL, bool id_in_tags=false, bool error_if_created=false)
get entities for TetGen data structure
virtual MoFEMErrorCode add_verices_in_the_middel_of_edges(const EntityHandle meshset, const BitRefLevel &bit, const bool recursive=false, int verb=-1, EntityHandle start_v=0)
make vertices in the middle of edges in meshset and add them to refinement levels defined by bit ...
multi_index_container< boost::shared_ptr< RefEntity >, indexed_by< ordered_unique< tag< Ent_mi_tag >, member< RefEntity::BasicEntity, EntityHandle, &RefEntity::ent > >, ordered_non_unique< tag< Ent_Ent_mi_tag >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > >, ordered_non_unique< tag< EntType_mi_tag >, const_mem_fun< RefEntity::BasicEntity, EntityType, &RefEntity::getEntType > >, ordered_non_unique< tag< ParentEntType_mi_tag >, const_mem_fun< RefEntity, EntityType, &RefEntity::getParentEntType > >, ordered_non_unique< tag< Composite_EntType_and_ParentEntType_mi_tag >, composite_key< RefEntity, const_mem_fun< RefEntity::BasicEntity, EntityType, &RefEntity::getEntType >, const_mem_fun< RefEntity, EntityType, &RefEntity::getParentEntType > > >, ordered_non_unique< tag< Composite_ParentEnt_And_EntType_mi_tag >, composite_key< RefEntity, const_mem_fun< RefEntity::BasicEntity, EntityType, &RefEntity::getEntType >, const_mem_fun< RefEntity, EntityHandle, &RefEntity::getParentEnt > > > > > RefEntity_multiIndex
#define _IT_CUBITMESHSETS_FOR_LOOP_(MESHSET_MANAGER, IT)
Iterator that loops over all the Cubit MeshSets in a moFEM field.
auto getVectorAdaptor
Get Vector adaptor.
Definition: Common.hpp:256
MoFEMErrorCode operator()(const std::string name, const Range &ents)
map< EntityHandle, TreeData > verticesOnCutEdges
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...
Definition: definitions.h:405
virtual MPI_Comm & get_comm() const =0
MoFEMErrorCode saveTrimEdges()
multi_index_container< ParentChild, indexed_by< hashed_unique< member< ParentChild, EntityHandle,&ParentChild::pArent > >, hashed_non_unique< member< ParentChild, EntityHandle,&ParentChild::cHild > > > > ParentChildMap
Definition: NodeMerger.hpp:150
map< unsigned long, EntityHandle > tetGenMoabMap
MoFEMErrorCode trimEdgesInTheMiddle(const BitRefLevel bit, Tag th=NULL, const double tol=1e-4, const bool debug=false)
trim edges
map< EntityHandle, TreeData > edgesToCut
MoFEMErrorCode setGeomData(tetgenio &in, moabTetGen_Map &moab_tetgen_map, tetGenMoab_Map &tetgen_moab_map, std::map< int, Range > &type_ents)
set point tags and type
VectorBoundedArray< double, 3 > VectorDouble3
Definition: Common.hpp:220
MoFEMErrorCode setVolume(const Range &volume)
set volume entities
bool getSucessMerge()
Return true if successful merge.
Definition: NodeMerger.hpp:44
MoFEMErrorCode moveMidNodesOnTrimmedEdges(Tag th=NULL)
move trimmed edges mid nodes