v0.15.0
Loading...
Searching...
No Matches
Namespaces | Macros
CutMeshInterface.cpp File Reference

Cut mesh by surface. More...

Go to the source code of this file.

Namespaces

namespace  MoFEM
 implementation of Data Operators for Forces and Sources
 

Macros

#define CutMeshFunctionBegin
 

Detailed Description

Cut mesh by surface.

Definition in file CutMeshInterface.cpp.

Macro Definition Documentation

◆ CutMeshFunctionBegin

#define CutMeshFunctionBegin
Value:
MOFEM_LOG_CHANNEL("WORLD"); \
MOFEM_LOG_FUNCTION(); \
MOFEM_LOG_TAG("WORLD", "CutMesh");
#define MoFEMFunctionBegin
First executable line of each MoFEM function, used for error handling. Final line of MoFEM functions ...

Definition at line 7 of file CutMeshInterface.cpp.

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