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