Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 :
4 : // This library is free software; you can redistribute it and/or
5 : // modify it under the terms of the GNU Lesser General Public
6 : // License as published by the Free Software Foundation; either
7 : // version 2.1 of the License, or (at your option) any later version.
8 :
9 : // This library is distributed in the hope that it will be useful,
10 : // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 : // Lesser General Public License for more details.
13 :
14 : // You should have received a copy of the GNU Lesser General Public
15 : // License along with this library; if not, write to the Free Software
16 : // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 :
18 :
19 :
20 : // C++ includes
21 : #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
22 : #include <cmath> // for std::acos()
23 : #include <algorithm>
24 : #include <limits>
25 : #include <map>
26 : #include <array>
27 :
28 : // Local includes
29 : #include "libmesh/boundary_info.h"
30 : #include "libmesh/function_base.h"
31 : #include "libmesh/cell_tet4.h"
32 : #include "libmesh/cell_tet10.h"
33 : #include "libmesh/cell_c0polyhedron.h"
34 : #include "libmesh/cell_polyhedron.h"
35 : #include "libmesh/elem_range.h"
36 : #include "libmesh/face_c0polygon.h"
37 : #include "libmesh/face_polygon.h"
38 : #include "libmesh/face_tri3.h"
39 : #include "libmesh/face_tri6.h"
40 : #include "libmesh/libmesh_logging.h"
41 : #include "libmesh/mesh_communication.h"
42 : #include "libmesh/mesh_modification.h"
43 : #include "libmesh/mesh_tools.h"
44 : #include "libmesh/node_range.h"
45 : #include "libmesh/parallel.h"
46 : #include "libmesh/parallel_ghost_sync.h"
47 : #include "libmesh/remote_elem.h"
48 : #include "libmesh/surface.h"
49 : #include "libmesh/enum_to_string.h"
50 : #include "libmesh/unstructured_mesh.h"
51 : #include "libmesh/elem_side_builder.h"
52 : #include "libmesh/tensor_value.h"
53 :
54 : namespace
55 : {
56 : using namespace libMesh;
57 :
58 2576 : bool split_first_diagonal(const Elem * elem,
59 : unsigned int diag_1_node_1,
60 : unsigned int diag_1_node_2,
61 : unsigned int diag_2_node_1,
62 : unsigned int diag_2_node_2)
63 : {
64 432 : return ((elem->node_id(diag_1_node_1) > elem->node_id(diag_2_node_1) &&
65 2936 : elem->node_id(diag_1_node_1) > elem->node_id(diag_2_node_2)) ||
66 2292 : (elem->node_id(diag_1_node_2) > elem->node_id(diag_2_node_1) &&
67 2652 : elem->node_id(diag_1_node_2) > elem->node_id(diag_2_node_2)));
68 : }
69 :
70 :
71 : // Return the local index of the vertex on \p elem with the highest
72 : // node id.
73 44518 : unsigned int highest_vertex_on(const Elem * elem)
74 : {
75 1792 : unsigned int highest_n = 0;
76 3584 : dof_id_type highest_n_id = elem->node_id(0);
77 356144 : for (auto n : make_range(1u, elem->n_vertices()))
78 : {
79 25088 : const dof_id_type n_id = elem->node_id(n);
80 311626 : if (n_id > highest_n_id)
81 : {
82 6928 : highest_n = n;
83 6928 : highest_n_id = n_id;
84 : }
85 : }
86 :
87 44518 : return highest_n;
88 : }
89 :
90 :
91 : static const std::array<std::array<unsigned int, 3>, 8> opposing_nodes =
92 : {{ {2,5,7},{3,4,6},{0,5,7},{1,4,6},{1,3,6},{0,2,7},{1,3,4},{0,2,5} }};
93 :
94 :
95 : // Find the highest id on these side nodes of this element
96 : std::pair<unsigned int, unsigned int>
97 201321 : split_diagonal(const Elem * elem,
98 : const std::vector<unsigned int> & nodes_on_side)
99 : {
100 6552 : libmesh_assert_equal_to(elem->type(), HEX8);
101 :
102 201321 : unsigned int highest_n = nodes_on_side.front();
103 13104 : dof_id_type highest_n_id = elem->node_id(nodes_on_side.front());
104 1006605 : for (auto n : nodes_on_side)
105 : {
106 26208 : const dof_id_type n_id = elem->node_id(n);
107 805284 : if (n_id > highest_n_id)
108 : {
109 11056 : highest_n = n;
110 11056 : highest_n_id = n_id;
111 : }
112 : }
113 :
114 409682 : for (auto n : nodes_on_side)
115 : {
116 1114909 : for (auto n2 : opposing_nodes[highest_n])
117 906548 : if (n2 == n)
118 6552 : return std::make_pair(highest_n, n2);
119 : }
120 :
121 0 : libmesh_error();
122 :
123 : return std::make_pair(libMesh::invalid_uint, libMesh::invalid_uint);
124 : }
125 :
126 :
127 : // Reconstruct a C++20 feature in C++14
128 : template <typename T>
129 : struct reversion_wrapper { T& iterable; };
130 :
131 : template <typename T>
132 4928 : auto begin (reversion_wrapper<T> w) {return std::rbegin(w.iterable);}
133 :
134 : template <typename T>
135 4928 : auto end (reversion_wrapper<T> w) {return std::rend(w.iterable);}
136 :
137 : template <typename T>
138 4928 : reversion_wrapper<T> reverse(T&& iterable) {return {iterable};}
139 :
140 : }
141 :
142 :
143 : namespace libMesh
144 : {
145 :
146 :
147 : // ------------------------------------------------------------
148 : // MeshTools::Modification functions for mesh modification
149 426 : void MeshTools::Modification::distort (MeshBase & mesh,
150 : const Real factor,
151 : const bool perturb_boundary)
152 : {
153 12 : libmesh_assert (mesh.n_nodes());
154 12 : libmesh_assert (mesh.n_elem());
155 12 : libmesh_assert ((factor >= 0.) && (factor <= 1.));
156 :
157 36 : LOG_SCOPE("distort()", "MeshTools::Modification");
158 :
159 : // If we are not perturbing boundary nodes, make a
160 : // quickly-searchable list of node ids we can check against.
161 24 : std::unordered_set<dof_id_type> boundary_node_ids;
162 426 : if (!perturb_boundary)
163 840 : boundary_node_ids = MeshTools::find_boundary_nodes (mesh);
164 :
165 : // Now calculate the minimum distance to
166 : // neighboring nodes for each node.
167 : // hmin holds these distances.
168 438 : std::vector<float> hmin (mesh.max_node_id(),
169 438 : std::numeric_limits<float>::max());
170 :
171 13674 : for (const auto & elem : mesh.active_element_ptr_range())
172 43425 : for (auto & n : elem->node_ref_range())
173 38700 : hmin[n.id()] = std::min(hmin[n.id()],
174 48831 : static_cast<float>(elem->hmin()));
175 :
176 : // Now actually move the nodes
177 : {
178 12 : const unsigned int seed = 123456;
179 :
180 : // seed the random number generator.
181 : // We'll loop from 1 to n_nodes on every processor, even those
182 : // that don't have a particular node, so that the pseudorandom
183 : // numbers will be the same everywhere.
184 426 : std::srand(seed);
185 :
186 : // If the node is on the boundary or
187 : // the node is not used by any element (hmin[n]<1.e20)
188 : // then we should not move it.
189 : // [Note: Testing for (in)equality might be wrong
190 : // (different types, namely float and double)]
191 10437 : for (auto n : make_range(mesh.max_node_id()))
192 11432 : if ((perturb_boundary || !boundary_node_ids.count(n)) && hmin[n] < 1.e20)
193 : {
194 : // the direction, random but unit normalized
195 1207 : Point dir (static_cast<Real>(std::rand())/static_cast<Real>(RAND_MAX),
196 2414 : (mesh.mesh_dimension() > 1) ? static_cast<Real>(std::rand())/static_cast<Real>(RAND_MAX) : 0.,
197 4828 : ((mesh.mesh_dimension() == 3) ? static_cast<Real>(std::rand())/static_cast<Real>(RAND_MAX) : 0.));
198 :
199 1207 : dir(0) = (dir(0)-.5)*2.;
200 : #if LIBMESH_DIM > 1
201 1207 : if (mesh.mesh_dimension() > 1)
202 1207 : dir(1) = (dir(1)-.5)*2.;
203 : #endif
204 : #if LIBMESH_DIM > 2
205 1207 : if (mesh.mesh_dimension() == 3)
206 852 : dir(2) = (dir(2)-.5)*2.;
207 : #endif
208 :
209 1207 : dir = dir.unit();
210 :
211 1207 : Node * node = mesh.query_node_ptr(n);
212 1207 : if (!node)
213 0 : continue;
214 :
215 1207 : (*node)(0) += dir(0)*factor*hmin[n];
216 : #if LIBMESH_DIM > 1
217 1207 : if (mesh.mesh_dimension() > 1)
218 1241 : (*node)(1) += dir(1)*factor*hmin[n];
219 : #endif
220 : #if LIBMESH_DIM > 2
221 1207 : if (mesh.mesh_dimension() == 3)
222 876 : (*node)(2) += dir(2)*factor*hmin[n];
223 : #endif
224 : }
225 : }
226 :
227 : // We haven't changed any topology, but just changing geometry could
228 : // have invalidated a point locator.
229 426 : mesh.clear_point_locator();
230 426 : }
231 :
232 :
233 :
234 298576 : void MeshTools::Modification::permute_elements(MeshBase & mesh)
235 : {
236 25176 : LOG_SCOPE("permute_elements()", "MeshTools::Modification");
237 :
238 : // We don't yet support doing permute() on a parent element, which
239 : // would require us to consistently permute all its children and
240 : // give them different local child numbers.
241 298576 : unsigned int n_levels = MeshTools::n_levels(mesh);
242 298576 : if (n_levels > 1)
243 0 : libmesh_error();
244 :
245 8392 : const unsigned int seed = 123456;
246 :
247 : // seed the random number generator.
248 : // We'll loop from 1 to max_elem_id on every processor, even those
249 : // that don't have a particular element, so that the pseudorandom
250 : // numbers will be the same everywhere.
251 298576 : std::srand(seed);
252 :
253 :
254 5260210 : for (auto e_id : make_range(mesh.max_elem_id()))
255 : {
256 4961634 : int my_rand = std::rand();
257 :
258 4961634 : Elem * elem = mesh.query_elem_ptr(e_id);
259 :
260 4961634 : if (!elem)
261 2882276 : continue;
262 :
263 2079358 : const unsigned int max_permutation = elem->n_permutations();
264 2079358 : if (!max_permutation)
265 8436 : continue;
266 :
267 2069678 : const unsigned int perm = my_rand % max_permutation;
268 :
269 2069678 : elem->permute(perm);
270 : }
271 298576 : }
272 :
273 :
274 2236 : void MeshTools::Modification::orient_elements(MeshBase & mesh)
275 : {
276 228 : LOG_SCOPE("orient_elements()", "MeshTools::Modification");
277 :
278 : // We don't yet support doing orient() on a parent element, which
279 : // would require us to consistently orient all its children and
280 : // give them different local child numbers.
281 2236 : unsigned int n_levels = MeshTools::n_levels(mesh);
282 2236 : if (n_levels > 1)
283 0 : libmesh_not_implemented_msg("orient_elements() does not support refined meshes");
284 :
285 76 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
286 92940 : for (auto elem : mesh.element_ptr_range())
287 45412 : elem->orient(&boundary_info);
288 2236 : }
289 :
290 :
291 :
292 300259 : void MeshTools::Modification::redistribute (MeshBase & mesh,
293 : const FunctionBase<Real> & mapfunc)
294 : {
295 8458 : libmesh_assert (mesh.n_nodes());
296 8458 : libmesh_assert (mesh.n_elem());
297 :
298 25374 : LOG_SCOPE("redistribute()", "MeshTools::Modification");
299 :
300 : Threads::parallel_for
301 300259 : (mesh.node_stored_range(),
302 600518 : [&mapfunc](const NodeRange & range)
303 : {
304 300259 : DenseVector<Real> output_vec(LIBMESH_DIM);
305 308717 : std::unique_ptr<FunctionBase<Real>> myfunc = mapfunc.clone();
306 5474265 : for (Node * node : range)
307 : {
308 5174006 : (*myfunc)(*node, output_vec);
309 :
310 5174006 : (*node)(0) = output_vec(0);
311 : #if LIBMESH_DIM > 1
312 5174006 : (*node)(1) = output_vec(1);
313 : #endif
314 : #if LIBMESH_DIM > 2
315 5174006 : (*node)(2) = output_vec(2);
316 : #endif
317 : }
318 583602 : });
319 :
320 : // If we just moved a mesh in or out out of the X axis or XY plane
321 : // then we might have changed its spatial_dimension()
322 8458 : mesh.unset_has_cached_elem_data();
323 :
324 : // We haven't changed any topology, but just changing geometry could
325 : // have invalidated a point locator.
326 300259 : mesh.clear_point_locator();
327 300259 : }
328 :
329 :
330 :
331 217 : void MeshTools::Modification::translate (MeshBase & mesh,
332 : const Real xt,
333 : const Real yt,
334 : const Real zt)
335 : {
336 8 : const Point p(xt, yt, zt);
337 :
338 : Threads::parallel_for_each
339 217 : (mesh.node_stored_range(),
340 1694 : [p](Node * node)
341 : {
342 1292 : *node += p;
343 1292 : });
344 :
345 : // If we just moved a mesh in or out out of the X axis or XY plane
346 : // then we might have changed its spatial_dimension()
347 8 : mesh.unset_has_cached_elem_data();
348 :
349 : // We haven't changed any topology, but just changing geometry could
350 : // have invalidated a point locator.
351 217 : mesh.clear_point_locator();
352 217 : }
353 :
354 :
355 : // void MeshTools::Modification::rotate2D (MeshBase & mesh,
356 : // const Real alpha)
357 : // {
358 : // libmesh_assert_not_equal_to (mesh.mesh_dimension(), 1);
359 :
360 : // const Real pi = std::acos(-1);
361 : // const Real a = alpha/180.*pi;
362 : // for (unsigned int n=0; n<mesh.n_nodes(); n++)
363 : // {
364 : // const Point p = mesh.node_ref(n);
365 : // const Real x = p(0);
366 : // const Real y = p(1);
367 : // const Real z = p(2);
368 : // mesh.node_ref(n) = Point(std::cos(a)*x - std::sin(a)*y,
369 : // std::sin(a)*x + std::cos(a)*y,
370 : // z);
371 : // }
372 :
373 : // }
374 :
375 :
376 :
377 : RealTensorValue
378 253633 : MeshTools::Modification::rotate (MeshBase & mesh,
379 : const Real phi,
380 : const Real theta,
381 : const Real psi)
382 : {
383 : // We won't change any topology, but just changing geometry could
384 : // invalidate a point locator.
385 253633 : mesh.clear_point_locator();
386 :
387 : #if LIBMESH_DIM == 3
388 253633 : const auto R = RealTensorValue::intrinsic_rotation_matrix(phi, theta, psi);
389 :
390 : Threads::parallel_for_each
391 253633 : (mesh.node_stored_range(),
392 723709 : [R](Node * node)
393 : {
394 242164 : Point & pt = *node;
395 5930479 : pt = R * pt;
396 5930479 : });
397 :
398 : // If we just moved a mesh in or out out of the X axis or XY plane
399 : // then we might have changed its spatial_dimension()
400 7126 : mesh.unset_has_cached_elem_data();
401 :
402 253633 : return R;
403 :
404 : #else
405 : libmesh_ignore(mesh, phi, theta, psi);
406 : libmesh_error_msg("MeshTools::Modification::rotate() requires libMesh to be compiled with LIBMESH_DIM==3");
407 : // We'll never get here
408 : return RealTensorValue();
409 : #endif
410 239381 : }
411 :
412 :
413 71 : void MeshTools::Modification::scale (MeshBase & mesh,
414 : const Real xs,
415 : const Real ys,
416 : const Real zs)
417 : {
418 2 : const Real x_scale = xs;
419 2 : Real y_scale = ys;
420 2 : Real z_scale = zs;
421 :
422 71 : if (ys == 0.)
423 : {
424 0 : libmesh_assert_equal_to (zs, 0.);
425 :
426 0 : y_scale = z_scale = x_scale;
427 : }
428 :
429 : // Scale the x coordinate in all dimensions
430 : Threads::parallel_for
431 71 : (mesh.node_stored_range(),
432 233 : [x_scale, y_scale, z_scale](const NodeRange & range)
433 : {
434 : if (LIBMESH_DIM < 2)
435 : for (auto & node : range)
436 : {
437 : (*node)(0) *= x_scale;
438 : }
439 : else if (LIBMESH_DIM < 3)
440 : for (auto & node : range)
441 : {
442 : (*node)(0) *= x_scale;
443 : (*node)(1) *= y_scale;
444 : }
445 : else
446 426 : for (auto & node : range)
447 : {
448 365 : (*node)(0) *= x_scale;
449 365 : (*node)(1) *= y_scale;
450 355 : (*node)(2) *= z_scale;
451 : }
452 69 : });
453 :
454 : // If we just collapsed a manifold onto the X axis or XY plane
455 : // then we might have changed its spatial_dimension()
456 2 : mesh.unset_has_cached_elem_data();
457 :
458 : // We haven't changed any topology, but just changing geometry could
459 : // have invalidated a point locator.
460 71 : mesh.clear_point_locator();
461 71 : }
462 :
463 :
464 :
465 3550 : void MeshTools::Modification::all_tri (MeshBase & mesh)
466 : {
467 300 : LOG_SCOPE("all_tri()", "MeshTools::Modification");
468 :
469 3550 : if (!mesh.is_replicated() && !mesh.is_prepared())
470 0 : mesh.prepare_for_use();
471 :
472 : // The number of elements in the original mesh before any additions
473 : // or deletions.
474 3550 : const dof_id_type n_orig_elem = mesh.n_elem();
475 3550 : const dof_id_type max_orig_id = mesh.max_elem_id();
476 :
477 : // We store pointers to the newly created elements in a vector
478 : // until they are ready to be added to the mesh. This is because
479 : // adding new elements on the fly can cause reallocation and invalidation
480 : // of existing mesh element_iterators.
481 300 : std::vector<std::unique_ptr<Elem>> new_elements;
482 :
483 3550 : unsigned int max_subelems = 1; // in 1D nothing needs to change
484 3550 : if (mesh.mesh_dimension() == 2) // in 2D quads can split into 2 tris
485 2840 : max_subelems = 2;
486 3550 : if (mesh.mesh_dimension() == 3) // in 3D hexes can split into 6 tets
487 710 : max_subelems = 6;
488 :
489 : // 2D polygons and 3D polyhedra can be split into an arbitrary
490 : // number of triangles/tetrahedra depending on their topology, so we
491 : // have to scan the mesh to find the largest split we will need.
492 174918 : for (const Elem * elem : mesh.element_ptr_range())
493 : {
494 87581 : if (const Polygon * poly = dynamic_cast<const Polygon *>(elem))
495 919 : max_subelems = std::max(max_subelems, poly->n_subtriangles());
496 86800 : else if (const Polyhedron * polyhedron = dynamic_cast<const Polyhedron *>(elem))
497 211 : max_subelems = std::max(max_subelems, polyhedron->n_subelements());
498 3350 : }
499 3550 : mesh.comm().max(max_subelems);
500 :
501 3550 : new_elements.reserve (max_subelems*n_orig_elem);
502 :
503 : // If the original mesh has *side* boundary data, we carry that over
504 : // to the new mesh with triangular elements. We currently only
505 : // support bringing over side-based BCs to the all-tri mesh, but
506 : // that could probably be extended to node and edge-based BCs as
507 : // well.
508 3550 : const bool mesh_has_boundary_data = (mesh.get_boundary_info().n_boundary_conds() > 0);
509 :
510 : // Temporary vectors to store the new boundary element pointers, side numbers, and boundary ids
511 200 : std::vector<Elem *> new_bndry_elements;
512 200 : std::vector<unsigned short int> new_bndry_sides;
513 200 : std::vector<boundary_id_type> new_bndry_ids;
514 :
515 : // We may need to add new points if we run into a 1.5th order
516 : // element; if we do that on a DistributedMesh in a ghost element then
517 : // we will need to fix their ids / unique_ids
518 3550 : bool added_new_ghost_point = false;
519 :
520 : // Iterate over the elements, splitting:
521 : // QUADs into pairs of conforming triangles
522 : // PYRAMIDs into pairs of conforming tets,
523 : // PRISMs into triplets of conforming tets, and
524 : // HEXs into quintets or sextets of conforming tets.
525 : // We split on the shortest diagonal to give us better
526 : // triangle quality in 2D, and we split based on node ids
527 : // to guarantee consistency in 3D.
528 : // C0POLYGONs into their sub-triangles
529 : // C0POLYHEDRA into their sub-elements (currently only tets)
530 :
531 : // FIXME: This algorithm does not work on refined grids!
532 : {
533 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
534 3550 : unique_id_type max_unique_id = mesh.parallel_max_unique_id();
535 : #endif
536 :
537 : // For avoiding extraneous allocation when building side elements
538 3550 : std::unique_ptr<const Elem> elem_side, subside_elem;
539 :
540 174918 : for (auto & elem : mesh.element_ptr_range())
541 : {
542 87581 : const ElemType etype = elem->type();
543 :
544 : // all_tri currently only works on coarse meshes
545 91203 : if (elem->parent())
546 0 : libmesh_not_implemented_msg("Cannot convert a refined element into simplices\n");
547 :
548 : // The new elements we will split the original into. Reserving
549 : // for the maximum number of sub-elements created for each element
550 90269 : std::vector<std::unique_ptr<Elem>> subelem(max_subelems);
551 :
552 16237 : auto set_nodes = [&elem, &subelem]
553 188368 : (const std::initializer_list<std::initializer_list<int>> & node_ids) {
554 870 : int i=0;
555 54499 : for (auto row : node_ids)
556 : {
557 1756 : int j=0;
558 36522 : Elem * sub = subelem[i++].get();
559 1756 : libmesh_assert(sub);
560 186628 : for (auto node_id : row)
561 157478 : sub->set_node(j++, elem->node_ptr(node_id));
562 : }
563 101936 : };
564 :
565 87581 : switch (etype)
566 : {
567 12451 : case QUADSHELL4:
568 : case QUAD4:
569 : {
570 12451 : const ElemType tri_type = (etype == QUAD4) ? TRI3 : TRISHELL3;
571 23794 : subelem[0] = Elem::build(tri_type);
572 23794 : subelem[1] = Elem::build(tri_type);
573 :
574 : // Check for possible edge swap
575 13005 : if ((elem->point(0) - elem->point(2)).norm() <
576 12451 : (elem->point(1) - elem->point(3)).norm())
577 4102 : set_nodes({{0,1,2},{0,2,3}});
578 : else
579 8349 : set_nodes({{0,1,3},{1,2,3}});
580 :
581 554 : break;
582 : }
583 :
584 142 : case QUAD8:
585 : {
586 146 : if (elem->processor_id() != mesh.processor_id())
587 118 : added_new_ghost_point = true;
588 :
589 276 : subelem[0] = Elem::build(TRI6);
590 276 : subelem[1] = Elem::build(TRI6);
591 :
592 : // Add a new node at the center (vertex average) of the element.
593 150 : Node * new_node = mesh.add_point((mesh.point(elem->node_id(0)) +
594 150 : mesh.point(elem->node_id(1)) +
595 150 : mesh.point(elem->node_id(2)) +
596 146 : mesh.point(elem->node_id(3)))/4,
597 : DofObject::invalid_id,
598 142 : elem->processor_id());
599 :
600 : // Check for possible edge swap
601 146 : if ((elem->point(0) - elem->point(2)).norm() <
602 142 : (elem->point(1) - elem->point(3)).norm())
603 : {
604 0 : set_nodes({{0,1,2,4,5},{0,2,3,3,6,7}});
605 0 : subelem[0]->set_node(5, new_node);
606 0 : subelem[1]->set_node(3, new_node);
607 : }
608 : else
609 : {
610 142 : set_nodes({{3,0,1,7,4},{1,2,3,5,6}});
611 146 : subelem[0]->set_node(5, new_node);
612 146 : subelem[1]->set_node(5, new_node);
613 : }
614 :
615 4 : break;
616 : }
617 :
618 3964 : case QUAD9:
619 : {
620 7384 : subelem[0] = Elem::build(TRI6);
621 7384 : subelem[1] = Elem::build(TRI6);
622 :
623 : // Check for possible edge swap
624 4236 : if ((elem->point(0) - elem->point(2)).norm() <
625 3964 : (elem->point(1) - elem->point(3)).norm())
626 1648 : set_nodes({{0,1,2,4,5,8},{0,2,3,8,6,7}});
627 : else
628 2316 : set_nodes({{0,1,3,4,8,7},{1,2,3,5,6,8}});
629 :
630 272 : break;
631 : }
632 :
633 1792 : case HEX8:
634 : {
635 1792 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
636 :
637 : // Hexes all split into six tetrahedra
638 85452 : subelem[0] = Elem::build(TET4);
639 85452 : subelem[1] = Elem::build(TET4);
640 85452 : subelem[2] = Elem::build(TET4);
641 85452 : subelem[3] = Elem::build(TET4);
642 85452 : subelem[4] = Elem::build(TET4);
643 85452 : subelem[5] = Elem::build(TET4);
644 :
645 : // On faces, we choose the node with the highest
646 : // global id, and we split on the diagonal which
647 : // includes that node. This ensures that (even in
648 : // parallel, even on distributed meshes) the same
649 : // diagonal split will be chosen for elements on either
650 : // side of the same quad face.
651 44518 : const unsigned int highest_n = highest_vertex_on(elem);
652 :
653 : // opposing_node[n] is the local node number of the node
654 : // on the farthest corner of a hex8 from local node n
655 : static const std::array<unsigned int, 8> opposing_node =
656 : {6, 7, 4, 5, 2, 3, 0, 1};
657 :
658 : static const std::vector<std::vector<unsigned int>> sides_opposing_highest =
659 45153 : {{2,3,5},{3,4,5},{1,4,5},{1,2,5},{0,2,3},{0,3,4},{0,1,4},{0,1,2}};
660 : static const std::vector<std::vector<unsigned int>> nodes_neighboring_highest =
661 45153 : {{1,3,4},{0,2,5},{1,3,6},{0,2,7},{0,5,7},{1,4,6},{2,5,7},{3,4,6}};
662 :
663 : // Start by looking in three directions away from the
664 : // highest-id node. In each direction there will be two
665 : // different possibilities for the split depending on
666 : // how the opposing face nodes are numbered.
667 : //
668 : // This is tricky enough that I'm not going to worry
669 : // about manually keeping tets oriented; we'll just call
670 : // orient() on each as we go.
671 :
672 1792 : unsigned int next_subelem = 0;
673 179864 : for (auto side : sides_opposing_highest[highest_n])
674 : {
675 : const std::vector<unsigned int> nodes_on_side =
676 138930 : elem->nodes_on_side(side);
677 :
678 133554 : auto [dn, dn2] = split_diagonal(elem, nodes_on_side);
679 :
680 5376 : unsigned int split_on_neighbor = false;
681 341933 : for (auto n : nodes_neighboring_highest[highest_n])
682 302503 : if (dn == n || dn2 == n)
683 : {
684 4928 : split_on_neighbor = true;
685 4928 : break;
686 : }
687 :
688 : // Add one or two elements for each opposing side,
689 : // depending on whether the diagonal split there
690 : // connects to the neighboring diagonal split or
691 : // not.
692 133554 : if (split_on_neighbor)
693 : {
694 109356 : subelem[next_subelem]->set_node(0, elem->node_ptr(highest_n));
695 109356 : subelem[next_subelem]->set_node(1, elem->node_ptr(dn));
696 109356 : subelem[next_subelem]->set_node(2, elem->node_ptr(dn2));
697 150520 : for (auto n : nodes_on_side)
698 150520 : if (n != dn && n != dn2)
699 : {
700 109356 : subelem[next_subelem]->set_node(3, elem->node_ptr(n));
701 4928 : break;
702 : }
703 99500 : subelem[next_subelem]->orient(&boundary_info);
704 99500 : ++next_subelem;
705 :
706 109356 : subelem[next_subelem]->set_node(0, elem->node_ptr(highest_n));
707 109356 : subelem[next_subelem]->set_node(1, elem->node_ptr(dn));
708 109356 : subelem[next_subelem]->set_node(2, elem->node_ptr(dn2));
709 147980 : for (auto n : reverse(nodes_on_side))
710 147980 : if (n != dn && n != dn2)
711 : {
712 109356 : subelem[next_subelem]->set_node(3, elem->node_ptr(n));
713 4928 : break;
714 : }
715 99500 : subelem[next_subelem]->orient(&boundary_info);
716 99500 : ++next_subelem;
717 : }
718 : else
719 : {
720 34950 : subelem[next_subelem]->set_node(0, elem->node_ptr(highest_n));
721 34950 : subelem[next_subelem]->set_node(1, elem->node_ptr(dn));
722 34950 : subelem[next_subelem]->set_node(2, elem->node_ptr(dn2));
723 101267 : for (auto n : nodes_on_side)
724 337087 : for (auto n2 : nodes_neighboring_highest[highest_n])
725 268406 : if (n == n2)
726 : {
727 34950 : subelem[next_subelem]->set_node(3, elem->node_ptr(n));
728 33606 : goto break_both_loops;
729 : }
730 :
731 34054 : break_both_loops:
732 34054 : subelem[next_subelem]->orient(&boundary_info);
733 34054 : ++next_subelem;
734 : }
735 : }
736 :
737 : // At this point we've created between 3 and 6 tets.
738 : // What's left to do depends on how many.
739 :
740 : // If we just chopped off three vertices into three
741 : // tets, then the best way to split this hex would be
742 : // the symmetric five-split. Chop off the opposing
743 : // vertex too, and then the remaining interior is our
744 : // final tet.
745 44518 : if (next_subelem == 3)
746 : {
747 2470 : subelem[next_subelem]->set_node(0, elem->node_ptr(opposing_nodes[highest_n][0]));
748 2470 : subelem[next_subelem]->set_node(1, elem->node_ptr(opposing_nodes[highest_n][1]));
749 2470 : subelem[next_subelem]->set_node(2, elem->node_ptr(opposing_nodes[highest_n][2]));
750 2470 : subelem[next_subelem]->set_node(3, elem->node_ptr(opposing_node[highest_n]));
751 2470 : subelem[next_subelem]->orient(&boundary_info);
752 0 : ++next_subelem;
753 :
754 2470 : subelem[next_subelem]->set_node(0, elem->node_ptr(opposing_nodes[highest_n][0]));
755 2470 : subelem[next_subelem]->set_node(1, elem->node_ptr(opposing_nodes[highest_n][1]));
756 2470 : subelem[next_subelem]->set_node(2, elem->node_ptr(opposing_nodes[highest_n][2]));
757 2470 : subelem[next_subelem]->set_node(3, elem->node_ptr(highest_n));
758 2470 : subelem[next_subelem]->orient(&boundary_info);
759 0 : ++next_subelem;
760 :
761 : // We don't need the 6th tet after all
762 0 : subelem[next_subelem].reset();
763 0 : ++next_subelem;
764 : }
765 :
766 : // If we just chopped off one (or two) vertices into
767 : // tets, then the remaining gap is best (or only) filled
768 : // by pairing another tet with each.
769 44518 : if (next_subelem == 4 ||
770 : next_subelem == 5)
771 : {
772 90748 : for (auto side : sides_opposing_highest[highest_n])
773 : {
774 : const std::vector<unsigned int> nodes_on_side =
775 68943 : elem->nodes_on_side(side);
776 :
777 67767 : auto [dn, dn2] = split_diagonal(elem, nodes_on_side);
778 :
779 1176 : unsigned int split_on_neighbor = false;
780 191339 : for (auto n : nodes_neighboring_highest[highest_n])
781 163519 : if (dn == n || dn2 == n)
782 : {
783 728 : split_on_neighbor = true;
784 728 : break;
785 : }
786 :
787 : // The two !split_on_neighbor sides are where we
788 : // need the two remaining tets
789 67767 : if (!split_on_neighbor)
790 : {
791 27540 : subelem[next_subelem]->set_node(0, elem->node_ptr(highest_n));
792 27540 : subelem[next_subelem]->set_node(1, elem->node_ptr(dn));
793 27540 : subelem[next_subelem]->set_node(2, elem->node_ptr(dn2));
794 27540 : subelem[next_subelem]->set_node(3, elem->node_ptr(opposing_node[highest_n]));
795 26644 : subelem[next_subelem]->orient(&boundary_info);
796 26644 : ++next_subelem;
797 : }
798 : }
799 : }
800 :
801 : // Whether we got there by creating six tets from the
802 : // first for loop or by patching up the split afterward,
803 : // we should have considered six tets (possibly
804 : // including one deleted one...) at this point.
805 1792 : libmesh_assert(next_subelem == 6);
806 :
807 1792 : break;
808 : }
809 :
810 142 : case PRISM6:
811 : {
812 : // Prisms all split into three tetrahedra
813 276 : subelem[0] = Elem::build(TET4);
814 276 : subelem[1] = Elem::build(TET4);
815 276 : subelem[2] = Elem::build(TET4);
816 :
817 : // Triangular faces are not split.
818 :
819 : // On quad faces, we choose the node with the highest
820 : // global id, and we split on the diagonal which
821 : // includes that node. This ensures that (even in
822 : // parallel, even on distributed meshes) the same
823 : // diagonal split will be chosen for elements on either
824 : // side of the same quad face. It also ensures that we
825 : // always have a mix of "clockwise" and
826 : // "counterclockwise" split faces (two of one and one
827 : // of the other on each prism; this is useful since the
828 : // alternative all-clockwise or all-counterclockwise
829 : // face splittings can't be turned into tets without
830 : // adding more nodes
831 :
832 : // Split on 0-4 diagonal
833 142 : if (split_first_diagonal(elem, 0,4, 1,3))
834 : {
835 : // Split on 0-5 diagonal
836 142 : if (split_first_diagonal(elem, 0,5, 2,3))
837 : {
838 : // Split on 1-5 diagonal
839 142 : if (split_first_diagonal(elem, 1,5, 2,4))
840 71 : set_nodes({{0,4,5,3},{0,4,1,5},{0,1,2,5}});
841 : else // Split on 2-4 diagonal
842 : {
843 2 : libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
844 71 : set_nodes({{0,4,5,3},{0,4,2,5},{0,1,2,4}});
845 : }
846 : }
847 : else // Split on 2-3 diagonal
848 : {
849 0 : libmesh_assert (split_first_diagonal(elem, 2,3, 0,5));
850 :
851 : // 0-4 and 2-3 split implies 2-4 split
852 0 : libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
853 :
854 0 : set_nodes({{0,4,2,3},{3,4,2,5},{0,1,2,4}});
855 : }
856 : }
857 : else // Split on 1-3 diagonal
858 : {
859 0 : libmesh_assert (split_first_diagonal(elem, 1,3, 0,4));
860 :
861 : // Split on 0-5 diagonal
862 0 : if (split_first_diagonal(elem, 0,5, 2,3))
863 : {
864 : // 1-3 and 0-5 split implies 1-5 split
865 0 : libmesh_assert (split_first_diagonal(elem, 1,5, 2,4));
866 :
867 0 : set_nodes({{1,3,4,5},{1,0,3,5},{0,1,2,5}});
868 : }
869 : else // Split on 2-3 diagonal
870 : {
871 0 : libmesh_assert (split_first_diagonal(elem, 2,3, 0,5));
872 :
873 : // Split on 1-5 diagonal
874 0 : if (split_first_diagonal(elem, 1,5, 2,4))
875 0 : set_nodes({{0,1,2,3},{3,1,2,5},{1,3,4,5}});
876 : else // Split on 2-4 diagonal
877 : {
878 0 : libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
879 0 : set_nodes({{0,1,2,3},{2,3,4,5},{3,1,2,4}});
880 : }
881 : }
882 : }
883 :
884 4 : break;
885 : }
886 :
887 426 : case PRISM20:
888 : case PRISM21:
889 : libmesh_experimental(); // We should upgrade this to TET14...
890 : libmesh_fallthrough();
891 : case PRISM18:
892 : {
893 828 : subelem[0] = Elem::build(TET10);
894 828 : subelem[1] = Elem::build(TET10);
895 828 : subelem[2] = Elem::build(TET10);
896 :
897 : // Split on 0-4 diagonal
898 426 : if (split_first_diagonal(elem, 0,4, 1,3))
899 : {
900 : // Split on 0-5 diagonal
901 426 : if (split_first_diagonal(elem, 0,5, 2,3))
902 : {
903 : // Split on 1-5 diagonal
904 426 : if (split_first_diagonal(elem, 1,5, 2,4))
905 213 : set_nodes({{0,4,5,3,15,13,17,9,12,14},
906 : {0,4,1,5,15,10,6,17,13,16},
907 : {0,1,2,5,6,7,8,17,16,11}});
908 : else // Split on 2-4 diagonal
909 : {
910 6 : libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
911 :
912 213 : set_nodes({{0,4,5,3,15,13,17,9,12,14},
913 : {0,4,2,5,15,16,8,17,13,11},
914 : {0,1,2,4,6,7,8,15,10,16}});
915 : }
916 : }
917 : else // Split on 2-3 diagonal
918 : {
919 0 : libmesh_assert (split_first_diagonal(elem, 2,3, 0,5));
920 :
921 : // 0-4 and 2-3 split implies 2-4 split
922 0 : libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
923 :
924 0 : set_nodes({{0,4,2,3,15,16,8,9,12,17},
925 : {3,4,2,5,12,16,17,14,13,11},
926 : {0,1,2,4,6,7,8,15,10,16}});
927 : }
928 : }
929 : else // Split on 1-3 diagonal
930 : {
931 0 : libmesh_assert (split_first_diagonal(elem, 1,3, 0,4));
932 :
933 : // Split on 0-5 diagonal
934 0 : if (split_first_diagonal(elem, 0,5, 2,3))
935 : {
936 : // 1-3 and 0-5 split implies 1-5 split
937 0 : libmesh_assert (split_first_diagonal(elem, 1,5, 2,4));
938 :
939 0 : set_nodes({{1,3,4,5,15,12,10,16,14,13},
940 : {1,0,3,5,6,9,15,16,17,14},
941 : {0,1,2,5,6,7,8,17,16,11}});
942 : }
943 : else // Split on 2-3 diagonal
944 : {
945 0 : libmesh_assert (split_first_diagonal(elem, 2,3, 0,5));
946 :
947 : // Split on 1-5 diagonal
948 0 : if (split_first_diagonal(elem, 1,5, 2,4))
949 0 : set_nodes({{0,1,2,3,6,7,8,9,15,17},
950 : {3,1,2,5,15,7,17,14,16,11},
951 : {1,3,4,5,15,12,10,16,14,13}});
952 : else // Split on 2-4 diagonal
953 : {
954 0 : libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
955 :
956 0 : set_nodes({{0,1,2,3,6,7,8,9,15,17},
957 : {2,3,4,5,17,12,16,11,14,13},
958 : {3,1,2,4,15,7,17,12,10,16}});
959 : }
960 : }
961 : }
962 :
963 12 : break;
964 : }
965 :
966 426 : case PYRAMID5:
967 : {
968 : // Pyramids all split into two tetrahedra
969 828 : subelem[0] = Elem::build(TET4);
970 828 : subelem[1] = Elem::build(TET4);
971 :
972 : // Choose how to split the quad face in a way that will
973 : // be consistent from possibly-different-type elements
974 : // splitting from the other side
975 : //
976 : // Split on 0-2 diagonal
977 426 : if (split_first_diagonal(elem, 0,2, 1,3))
978 213 : set_nodes({{0,1,2,4},{0,2,3,4}});
979 : // Split on 1-3 diagonal
980 : else
981 : {
982 6 : libmesh_assert (split_first_diagonal(elem, 1,3, 0,2));
983 213 : set_nodes({{0,1,3,4},{1,2,3,4}});
984 : }
985 :
986 12 : break;
987 : }
988 :
989 426 : case PYRAMID14:
990 : {
991 : // Pyramids all split into two tetrahedra
992 828 : subelem[0] = Elem::build(TET10);
993 828 : subelem[1] = Elem::build(TET10);
994 :
995 : // Choose how to split the quad face in a way that will
996 : // be consistent from possibly-different-type elements
997 : // splitting from the other side
998 : //
999 : // Split on 0-2 diagonal
1000 426 : if (split_first_diagonal(elem, 0,2, 1,3))
1001 213 : set_nodes({{0,1,2,4,5,6,13,9,10,11},
1002 : {0,2,3,4,13,7,8,9,11,12}});
1003 : // Split on 1-3 diagonal
1004 : else
1005 : {
1006 6 : libmesh_assert (split_first_diagonal(elem, 1,3, 0,2));
1007 213 : set_nodes({{0,1,3,4,5,13,8,9,10,12},
1008 : {1,2,3,4,6,7,13,10,11,12}});
1009 : }
1010 :
1011 12 : break;
1012 : }
1013 :
1014 781 : case C0POLYGON:
1015 : {
1016 : // Split a C0Polygon into the triangles defined by its
1017 : // current triangulation. This relies on the user having
1018 : // a valid triangulation (the constructor sets a default
1019 : // one, and the user can refresh it via retriangulate()
1020 : // after moving nodes).
1021 781 : const C0Polygon * polygon = cast_ptr<const C0Polygon *>(elem);
1022 22 : const unsigned int n_subtri = polygon->n_subtriangles();
1023 2698 : for (unsigned int t = 0; t != n_subtri; ++t)
1024 : {
1025 1917 : const std::array<int, 3> tri = polygon->subtriangle(t);
1026 1917 : if (tri[0] < 0 || tri[1] < 0 || tri[2] < 0)
1027 0 : libmesh_not_implemented_msg
1028 : ("Cannot convert a C0Polygon whose triangulation\n"
1029 : "introduces special (non-vertex) points");
1030 1917 : subelem[t] = Elem::build(TRI3);
1031 2025 : subelem[t]->set_node(0, elem->node_ptr(tri[0]));
1032 2025 : subelem[t]->set_node(1, elem->node_ptr(tri[1]));
1033 2025 : subelem[t]->set_node(2, elem->node_ptr(tri[2]));
1034 : }
1035 :
1036 22 : break;
1037 : }
1038 :
1039 142 : case C0POLYHEDRON:
1040 : {
1041 : // Split a C0Polyhedron into the tetrahedra defined by its
1042 : // current tetrahedralization. If the polyhedron required
1043 : // a mid-element node, the user is expected to have added
1044 : // that node to the mesh during construction; we just
1045 : // reference it via the polyhedron's node pointers.
1046 : const C0Polyhedron * polyhedron =
1047 142 : cast_ptr<const C0Polyhedron *>(elem);
1048 4 : const unsigned int n_sub = polyhedron->n_subelements();
1049 1917 : for (unsigned int t = 0; t != n_sub; ++t)
1050 : {
1051 1775 : const std::array<int, 4> tet = polyhedron->subelement(t);
1052 1775 : if (tet[0] < 0 || tet[1] < 0 || tet[2] < 0 || tet[3] < 0)
1053 0 : libmesh_not_implemented_msg
1054 : ("Cannot convert a C0Polyhedron whose triangulation\n"
1055 : "introduces special (non-vertex) points");
1056 1775 : subelem[t] = Elem::build(TET4);
1057 1875 : subelem[t]->set_node(0, elem->node_ptr(tet[0]));
1058 1875 : subelem[t]->set_node(1, elem->node_ptr(tet[1]));
1059 1875 : subelem[t]->set_node(2, elem->node_ptr(tet[2]));
1060 1875 : subelem[t]->set_node(3, elem->node_ptr(tet[3]));
1061 : }
1062 : // There is a concern that two neighbor polyhedra might have
1063 : // a triangulation of a side that does not match. But the
1064 : // default triangulation is based on the side's triangulation
1065 : // and the side element is supposed to be shared (that's why
1066 : // shared pointers to polygons are used to build the polyhedra).
1067 : // So the default one should work.
1068 :
1069 4 : break;
1070 : }
1071 :
1072 : // No need to split elements that are already simplicial:
1073 24163 : case EDGE2:
1074 : case EDGE3:
1075 : case EDGE4:
1076 : case TRI3:
1077 : case TRI6:
1078 : case TRI7:
1079 : case TET4:
1080 : case TET10:
1081 : case TET14:
1082 : case INFEDGE2:
1083 : // No way to split infinite quad/prism elements, so
1084 : // hopefully no need to
1085 : case INFQUAD4:
1086 : case INFQUAD6:
1087 : case INFPRISM6:
1088 : case INFPRISM12:
1089 1868 : continue;
1090 : // If we're left with an unimplemented element we're
1091 : // probably out of luck. TODO: implement hex20, hex27,
1092 : // pyramid13,...
1093 0 : default:
1094 0 : libmesh_not_implemented_msg
1095 : ("Error, encountered unimplemented element "
1096 : << Utility::enum_to_string<ElemType>(etype)
1097 934 : << " in MeshTools::Modification::all_tri()...");
1098 22295 : } // end switch (etype)
1099 :
1100 : // Be sure the correct data is set for all subelems.
1101 63418 : const unsigned int nei = elem->n_extra_integers();
1102 377272 : for (unsigned int i=0; i != max_subelems; ++i)
1103 326650 : if (subelem[i]) {
1104 304852 : subelem[i]->processor_id() = elem->processor_id();
1105 304852 : subelem[i]->subdomain_id() = elem->subdomain_id();
1106 :
1107 : // Copy any extra element data. Since the subelements
1108 : // haven't been added to the mesh yet any allocation has
1109 : // to be done manually.
1110 304852 : subelem[i]->add_extra_integers(nei);
1111 311936 : for (unsigned int ei=0; ei != nei; ++ei)
1112 7588 : subelem[ei]->set_extra_integer(ei, elem->get_extra_integer(ei));
1113 :
1114 :
1115 : // Copy any mapping data.
1116 317464 : subelem[i]->set_mapping_type(elem->mapping_type());
1117 25224 : subelem[i]->set_mapping_data(elem->mapping_data());
1118 : }
1119 :
1120 : // On a mesh with boundary data, we need to move that data to
1121 : // the new elements.
1122 :
1123 : // On a mesh which is distributed, we need to move
1124 : // remote_elem links to the new elements.
1125 63418 : bool mesh_is_serial = mesh.is_serial();
1126 :
1127 63418 : if (mesh_has_boundary_data || !mesh_is_serial)
1128 : {
1129 : // Container to key boundary IDs handed back by the BoundaryInfo object.
1130 500 : std::vector<boundary_id_type> bc_ids;
1131 :
1132 123224 : for (auto sn : elem->side_index_range())
1133 : {
1134 102649 : mesh.get_boundary_info().boundary_ids(elem, sn, bc_ids);
1135 :
1136 102649 : if (bc_ids.empty() && elem->neighbor_ptr(sn) != remote_elem)
1137 81617 : continue;
1138 :
1139 : // Make a sorted list of node ids for elem->side(sn)
1140 21032 : elem->build_side_ptr(elem_side, sn);
1141 21418 : std::vector<dof_id_type> elem_side_nodes(elem_side->n_nodes());
1142 23254 : for (unsigned int esn=0,
1143 772 : n_esn = cast_int<unsigned int>(elem_side_nodes.size());
1144 97424 : esn != n_esn; ++esn)
1145 77696 : elem_side_nodes[esn] = elem_side->node_id(esn);
1146 21032 : std::sort(elem_side_nodes.begin(), elem_side_nodes.end());
1147 :
1148 119750 : for (unsigned int i=0; i != max_subelems; ++i)
1149 100286 : if (subelem[i])
1150 : {
1151 406589 : for (auto subside : subelem[i]->side_index_range())
1152 : {
1153 325232 : subelem[i]->build_side_ptr(subside_elem, subside);
1154 :
1155 : // Make a list of *vertex* node ids for this subside, see if they are all present
1156 : // in elem->side(sn). Note 1: we can't just compare elem->key(sn) to
1157 : // subelem[i]->key(subside) in the Prism cases, since the new side is
1158 : // a different type. Note 2: we only use vertex nodes since, in the future,
1159 : // a Hex20 or Prism15's QUAD8 face may be split into two Tri6 faces, and the
1160 : // original face will not contain the mid-edge node.
1161 325232 : std::vector<dof_id_type> subside_nodes(subside_elem->n_vertices());
1162 338704 : for (unsigned int ssn=0,
1163 8512 : n_ssn = cast_int<unsigned int>(subside_nodes.size());
1164 1219476 : ssn != n_ssn; ++ssn)
1165 909492 : subside_nodes[ssn] = subside_elem->node_id(ssn);
1166 320976 : std::sort(subside_nodes.begin(), subside_nodes.end());
1167 :
1168 : // std::includes returns true if every element of the second sorted range is
1169 : // contained in the first sorted range.
1170 320976 : if (std::includes(elem_side_nodes.begin(), elem_side_nodes.end(),
1171 : subside_nodes.begin(), subside_nodes.end()))
1172 : {
1173 44250 : for (const auto & b_id : bc_ids)
1174 12853 : if (b_id != BoundaryInfo::invalid_id)
1175 : {
1176 12853 : new_bndry_ids.push_back(b_id);
1177 13307 : new_bndry_elements.push_back(subelem[i].get());
1178 12853 : new_bndry_sides.push_back(subside);
1179 : }
1180 :
1181 : // If the original element had a RemoteElem neighbor on side 'sn',
1182 : // then the subelem has one on side 'subside'.
1183 31875 : if (elem->neighbor_ptr(sn) == remote_elem)
1184 72 : subelem[i]->set_neighbor(subside, const_cast<RemoteElem*>(remote_elem));
1185 : }
1186 : }
1187 : } // end for loop over subelem
1188 : } // end for loop over sides
1189 :
1190 : // Remove the original element from the BoundaryInfo structure.
1191 20325 : mesh.get_boundary_info().remove(elem);
1192 :
1193 : } // end if (mesh_has_boundary_data)
1194 :
1195 : // Determine new IDs for the split elements which will be
1196 : // the same on all processors, therefore keeping the Mesh
1197 : // in sync. Note: we offset the new IDs by max_orig_id to
1198 : // avoid overwriting any of the original IDs.
1199 377272 : for (unsigned int i=0; i != max_subelems; ++i)
1200 326650 : if (subelem[i])
1201 : {
1202 : // Determine new IDs for the split elements which will be
1203 : // the same on all processors, therefore keeping the Mesh
1204 : // in sync. Note: we offset the new IDs by the max of the
1205 : // pre-existing ids to avoid conflicting with originals.
1206 304852 : subelem[i]->set_id( max_orig_id + max_subelems*elem->id() + i );
1207 :
1208 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
1209 304852 : subelem[i]->set_unique_id(max_unique_id + max_subelems*elem->unique_id() + i);
1210 : #endif
1211 :
1212 : // Prepare to add the newly-created simplices
1213 12612 : new_elements.push_back(std::move(subelem[i]));
1214 : }
1215 :
1216 : // Delete the original element
1217 63418 : mesh.delete_elem(elem);
1218 83687 : } // End for loop over elements
1219 3350 : } // end scope
1220 :
1221 :
1222 : // Now, iterate over the new elements vector, and add them each to
1223 : // the Mesh.
1224 308402 : for (auto & elem : new_elements)
1225 330076 : mesh.add_elem(std::move(elem));
1226 :
1227 3550 : if (mesh_has_boundary_data)
1228 : {
1229 : // If the old mesh had boundary data, the new mesh better have
1230 : // some. However, we can't assert that the size of
1231 : // new_bndry_elements vector is > 0, since we may not have split
1232 : // any elements actually on the boundary. We also can't assert
1233 : // that the original number of boundary sides is equal to the
1234 : // sum of the boundary sides currently in the mesh and the
1235 : // newly-added boundary sides, since in 3D, we may have split a
1236 : // boundary QUAD into two boundary TRIs. Therefore, we won't be
1237 : // too picky about the actual number of BCs, and just assert that
1238 : // there are some, somewhere.
1239 : #ifndef NDEBUG
1240 50 : bool nbe_nonempty = new_bndry_elements.size();
1241 50 : mesh.comm().max(nbe_nonempty);
1242 50 : libmesh_assert(nbe_nonempty ||
1243 : mesh.get_boundary_info().n_boundary_conds()>0);
1244 : #endif
1245 :
1246 : // We should also be sure that the lengths of the new boundary data vectors
1247 : // are all the same.
1248 50 : libmesh_assert_equal_to (new_bndry_elements.size(), new_bndry_sides.size());
1249 50 : libmesh_assert_equal_to (new_bndry_sides.size(), new_bndry_ids.size());
1250 :
1251 : // Add the new boundary info to the mesh
1252 14628 : for (auto s : index_range(new_bndry_elements))
1253 13761 : mesh.get_boundary_info().add_side(new_bndry_elements[s],
1254 908 : new_bndry_sides[s],
1255 908 : new_bndry_ids[s]);
1256 : }
1257 :
1258 : // In a DistributedMesh any newly added ghost node ids may be
1259 : // inconsistent, and unique_ids of newly added ghost nodes remain
1260 : // unset.
1261 : // make_nodes_parallel_consistent() will fix all this.
1262 3550 : if (!mesh.is_serial())
1263 : {
1264 1287 : mesh.comm().max(added_new_ghost_point);
1265 :
1266 1287 : if (added_new_ghost_point)
1267 0 : MeshCommunication().make_nodes_parallel_consistent (mesh);
1268 : }
1269 :
1270 : // Prepare the newly created mesh for use.
1271 3550 : mesh.prepare_for_use();
1272 :
1273 : // Let the new_elements and new_bndry_elements vectors go out of scope.
1274 3818 : }
1275 :
1276 :
1277 :
1278 2130 : void MeshTools::Modification::all_rbb (MeshBase & mesh)
1279 : {
1280 180 : LOG_SCOPE("all_rbb()", "MeshTools::Modification");
1281 :
1282 : // By default, use 1.0 as the weight on every RATIONAL_BERNSTEIN
1283 : // mapped node
1284 2130 : const Real default_weight = 1.0;
1285 :
1286 : const auto weight_index =
1287 4200 : (mesh.add_node_datum<Real>("rational_weight", true,
1288 : &default_weight));
1289 :
1290 60 : mesh.set_default_mapping_type(RATIONAL_BERNSTEIN_MAP);
1291 2130 : mesh.set_default_mapping_data(weight_index);
1292 :
1293 : // Out of loop to reduce heap allocations
1294 2130 : std::unique_ptr<Elem> edge_ptr, face_ptr;
1295 :
1296 57426 : for (auto & elem : mesh.element_ptr_range())
1297 : {
1298 28397 : if (elem->level())
1299 0 : libmesh_not_implemented_msg
1300 : ("all_rbb() currently only supports flat meshes with no refinement levels");
1301 :
1302 : #ifdef LIBMESH_ENABLE_INFINITE_ELEMENTS
1303 4460 : if (elem->infinite())
1304 0 : libmesh_not_implemented_msg
1305 : ("all_rbb() currently only supports finite geometric elements");
1306 : #endif
1307 :
1308 28397 : elem->set_mapping_type(RATIONAL_BERNSTEIN_MAP);
1309 1784 : elem->set_mapping_data(weight_index);
1310 :
1311 : // Nothing to do unless we have curves to correct
1312 28397 : if (elem->default_order() == FIRST)
1313 2920 : continue;
1314 :
1315 : // Modify the center node of an "edge" - possibly an actual edge
1316 : // element's node, possibly a center node between points on a
1317 : // face's or cell's edge - for RBB interpolation. This should fit
1318 : // a circular curve exactly in cases where the original nodes are
1319 : // equispaced and the outer nodes' weights are equal, and should
1320 : // be a good fit otherwise.
1321 : //
1322 : // We want to use this to interpolate "internal" conceptual
1323 : // edges of a Hex27 too, so we'll handle the cases where w0 and
1324 : // w1 aren't 1, as well as the cases where the Nodes n0 and n1
1325 : // are already control points which don't match their
1326 : // corresponding physical points.
1327 129773 : auto make_edge_rbb = [default_weight, weight_index]
1328 : (const Node & n0, const Node & n1, Node & n_center,
1329 167144 : const Point & p0, const Point & p1)
1330 : {
1331 : // Skip edges we've already modified; the center node for
1332 : // these is no longer at the curve point we wish to
1333 : // interpolate, it should already be at the control point that
1334 : // accomplishes the interpolation.
1335 118136 : const Real old_weight = n_center.get_extra_datum<Real>(weight_index);
1336 118136 : if (old_weight != default_weight)
1337 1588 : return;
1338 :
1339 5332 : Point & p2 = n_center;
1340 :
1341 97742 : const Real w0 = n0.get_extra_datum<Real>(weight_index);
1342 97742 : const Real w1 = n1.get_extra_datum<Real>(weight_index);
1343 :
1344 5332 : const Point e02 = p2-p0,
1345 5332 : e21 = p1-p2;
1346 5332 : const Real chord_02_len_sq = e02.norm_sq(),
1347 5332 : chord_21_len_sq = e21.norm_sq();
1348 :
1349 : // First find the cosine of phi, the angle between our two
1350 : // subchords (turning from the direction of one to the
1351 : // direction of the other; this is the supplementary angle to
1352 : // the angle at the midpoint). This is the same as half of
1353 : // the angle of our circular arc, which nicely enough is also
1354 : // the angle we take cos and sec of in NURBS formulae
1355 97742 : const Real cos_phi = (e02*e21)/std::sqrt(chord_02_len_sq*chord_21_len_sq);
1356 :
1357 : // There's a way to do really large arcs using negative
1358 : // weights, but we're going to get lousy approximation quality
1359 : // from isoparametric elements if we go too low, as well as
1360 : // bad numerics here, so let's just disallow it.
1361 97742 : if (cos_phi < 0.5)
1362 0 : libmesh_not_implemented_msg
1363 : ("all_rbb() is not recommended for extremely sharp curves on one edge");
1364 :
1365 97742 : const Real w_center = cos_phi*std::sqrt(w0*w1);
1366 :
1367 92410 : n_center.set_extra_datum<Real>(weight_index, w_center);
1368 :
1369 : // Now let's get the control point location. This comes from
1370 : // a lot of back-and-forth with Gemini, but fortunately I'm
1371 : // rewriting it after I've already added unit tests that
1372 : // should scream if it's badly wrong.
1373 97742 : const Real w_mid = w0/4 + w1/4 + w_center/2;
1374 97742 : n_center *= 2*w_mid;
1375 97742 : n_center -= (w0 * p0 + w1 * p1)/2;
1376 5332 : n_center /= w_center;
1377 25477 : };
1378 :
1379 128314 : auto make_face_rbb = [weight_index] (Elem & face)
1380 : {
1381 : // Prisms and pyramids may need to skip some faces while
1382 : // adjusting others
1383 20659 : if (face.type() == TRI6)
1384 0 : return;
1385 :
1386 20659 : if (face.type() != QUAD9)
1387 0 : libmesh_not_implemented_msg
1388 : ("all_rbb() currently only supports mid-face nodes on Quad9 faces");
1389 :
1390 : // We only use [4,8) but matching indices is nice and stack is
1391 : // cheap.
1392 : Real w[9];
1393 :
1394 103295 : for (unsigned int i : make_range(4u, 8u))
1395 86996 : w[i] = face.node_ref(i).get_extra_datum<Real>(weight_index);
1396 :
1397 : // We can't currently handle arbitrary vertex weights
1398 : #ifndef NDEBUG
1399 5450 : for (unsigned int i : make_range(4u))
1400 4360 : libmesh_assert_equal_to
1401 : (face.node_ref(i).get_extra_datum<Real>(weight_index), 1);
1402 : #endif
1403 :
1404 : // For the mid-face point, if we want to exactly match
1405 : // any cylinders and cones and spheres, we're actually already
1406 : // entirely constrained by the other points.
1407 : //
1408 : // This formula gives the minimum-energy Steiner surface based
1409 : // on the outer 8 points.
1410 : //
1411 : // That's an isogeometric representation of a cylinder aligned
1412 : // to either axis, or of a sphere where the quad edges are on
1413 : // latitude/longitude lines, or of a cone where two edges are
1414 : // segments of cone generating lines and the other two are
1415 : // arcs perpendicular to the axis.
1416 : //
1417 : // It's not perfectly isogeometric for the spheres we generate
1418 : // (where the quad edges are all great circles), but it should
1419 : // still converge asymptotically faster than non-rational
1420 : // quadratic Lagrange.
1421 2180 : const Point xi_avg = (face.point(7) + face.point(5))/2;
1422 1090 : const Point eta_avg = (face.point(4) + face.point(6))/2;
1423 1090 : const Point vertex_avg = (face.point(0) + face.point(1) +
1424 2180 : face.point(2) + face.point(3))/4;
1425 :
1426 20659 : const Real w_xi = (w[7] + w[5])/2;
1427 20659 : const Real w_eta = (w[4] + w[6])/2;
1428 20659 : const Real w_mid = w_xi * w_eta;
1429 :
1430 1090 : Node & midnode = face.node_ref(8);
1431 20659 : midnode.set_extra_datum<Real>(weight_index, w_mid);
1432 20659 : midnode = ((1+w_mid)/(w_xi+w_eta) * (w_xi*xi_avg + w_eta*eta_avg) - vertex_avg)/w_mid;
1433 25477 : };
1434 :
1435 : // If we're on a Hex27, our formula for the mid-volume node
1436 : // relies on the locations of the mid-face points. We could
1437 : // re-calculate those later but let's just save them now.
1438 178339 : Point midfacepts[6];
1439 25477 : if (elem->type() == HEX27)
1440 16366 : for (auto i : make_range(6))
1441 14652 : midfacepts[i] = elem->point(20+i);
1442 :
1443 : // Check each edge for a curve, and adjust it if needed.
1444 137599 : for (auto e : elem->edge_index_range())
1445 : {
1446 110454 : elem->build_edge_ptr(edge_ptr, e);
1447 :
1448 : // We should add EDGE4 once we have QUAD16/TRI10/HEX64 to
1449 : // use it
1450 110454 : if (edge_ptr->type() != EDGE3)
1451 0 : libmesh_not_implemented_msg
1452 : ("all_rbb() currently only supports meshes with 2- and/or 3-node edges");
1453 :
1454 123550 : make_edge_rbb(edge_ptr->node_ref(0), edge_ptr->node_ref(1),
1455 : edge_ptr->node_ref(2),
1456 6548 : edge_ptr->node_ref(0), edge_ptr->node_ref(1));
1457 :
1458 : }
1459 :
1460 : // If we're in 3D, we may have face nodes that also need to be
1461 : // adjusted to replace an interpolated curve with a spline
1462 : // curve. We know what to do with a quad face, but we'll have
1463 : // to scream and die if we see a Tri7 face node.
1464 25681 : bool check_face_points = (elem->dim() > 2) &&
1465 5035 : (elem->n_nodes() > elem->n_edges() + elem->n_vertices());
1466 :
1467 1668 : if (check_face_points)
1468 16470 : for (auto f : elem->side_index_range())
1469 : {
1470 : // Prisms and pyramids may need to skip some faces while
1471 : // adjusting others
1472 14028 : if (elem->side_type(f) == TRI6)
1473 0 : continue;
1474 :
1475 14028 : elem->build_side_ptr(face_ptr, f);
1476 :
1477 14028 : make_face_rbb(*face_ptr);
1478 : }
1479 :
1480 : bool check_interior_points =
1481 25477 : elem->n_nodes() > elem->n_edges() + elem->n_vertices() + elem->n_faces();
1482 :
1483 25477 : if (check_interior_points)
1484 : {
1485 9637 : if (elem->type() == EDGE3)
1486 : {
1487 788 : make_edge_rbb(elem->node_ref(0), elem->node_ref(1),
1488 : elem->node_ref(2),
1489 668 : elem->node_ref(0), elem->node_ref(1));
1490 : }
1491 8969 : else if (elem->dim() == 2)
1492 : {
1493 6631 : make_face_rbb(*elem);
1494 : }
1495 2338 : else if (elem->type() == HEX27)
1496 : {
1497 : // We still have the midnode left to go. We want
1498 : // something here that will preserve the tensor product
1499 : // structure for 2.5D extrusions of IGA faces, but also
1500 : // be at least near to the minimum-energy control point
1501 : // and weight for general cases. We'll treat opposing
1502 : // mid-face nodes as the endpoints of a (more general
1503 : // than our edges, since they might have non-1 weights)
1504 : // Edge3, and see what we'd need on the midnode to
1505 : // interpolate the center point with them. If we've got
1506 : // something isogeometric like an extrusion then our
1507 : // results should agree; for a quick-but-good output in
1508 : // general we'll take an average.
1509 2338 : const int opposite_sides[3][2] = {{0,5}, {1,3}, {2,4}};
1510 :
1511 2338 : Node & midnode = elem->node_ref(26);
1512 2338 : const Point original_midpoint = midnode;
1513 :
1514 : // Averaging in projective space
1515 104 : Point sum_weighted_point = 0;
1516 104 : Real sum_weight = 0;
1517 :
1518 9352 : for (int i : make_range(3))
1519 : {
1520 7014 : Node & n0 = elem->node_ref(20+opposite_sides[i][0]);
1521 7014 : Node & n1 = elem->node_ref(20+opposite_sides[i][1]);
1522 :
1523 7014 : make_edge_rbb(n0, n1, midnode,
1524 7014 : midfacepts[opposite_sides[i][0]],
1525 7014 : midfacepts[opposite_sides[i][1]]);
1526 :
1527 : const Real midweight =
1528 7014 : midnode.get_extra_datum<Real>(weight_index);
1529 7014 : sum_weight += midweight;
1530 312 : sum_weighted_point += midweight * midnode;
1531 :
1532 : // Reset for next run
1533 312 : midnode = original_midpoint;
1534 6702 : midnode.set_extra_datum<Real>(weight_index,
1535 : default_weight);
1536 :
1537 : }
1538 :
1539 2338 : const Real midweight = sum_weight/3;
1540 2338 : midnode.set_extra_datum<Real>(weight_index,
1541 : midweight);
1542 :
1543 104 : midnode = sum_weighted_point / 3 / midweight;
1544 : }
1545 : else
1546 0 : libmesh_not_implemented_msg
1547 : ("all_rbb() doesn't yet support " << elem->type());
1548 : }
1549 2010 : }
1550 2130 : }
1551 :
1552 :
1553 :
1554 0 : void MeshTools::Modification::smooth (MeshBase & mesh,
1555 : const unsigned int n_iterations,
1556 : const Real power)
1557 : {
1558 : /**
1559 : * This implementation assumes every element "side" has only 2 nodes.
1560 : */
1561 0 : libmesh_assert_equal_to (mesh.mesh_dimension(), 2);
1562 :
1563 : /*
1564 : * Create a quickly-searchable list of boundary nodes.
1565 : */
1566 : std::unordered_set<dof_id_type> boundary_node_ids =
1567 0 : MeshTools::find_boundary_nodes (mesh);
1568 :
1569 : // For avoiding extraneous element side allocation
1570 0 : ElemSideBuilder side_builder;
1571 :
1572 0 : for (unsigned int iter=0; iter<n_iterations; iter++)
1573 : {
1574 : /*
1575 : * loop over the mesh refinement level
1576 : */
1577 0 : unsigned int n_levels = MeshTools::n_levels(mesh);
1578 0 : for (unsigned int refinement_level=0; refinement_level != n_levels;
1579 : refinement_level++)
1580 : {
1581 : // initialize the storage (have to do it on every level to get empty vectors
1582 0 : std::vector<Point> new_positions;
1583 0 : std::vector<Real> weight;
1584 0 : new_positions.resize(mesh.n_nodes());
1585 0 : weight.resize(mesh.n_nodes());
1586 :
1587 : {
1588 : // Loop over the elements to calculate new node positions
1589 0 : for (const auto & elem : as_range(mesh.level_elements_begin(refinement_level),
1590 0 : mesh.level_elements_end(refinement_level)))
1591 : {
1592 : /*
1593 : * We relax all nodes on level 0 first
1594 : * If the element is refined (level > 0), we interpolate the
1595 : * parents nodes with help of the embedding matrix
1596 : */
1597 0 : if (refinement_level == 0)
1598 : {
1599 0 : for (auto s : elem->side_index_range())
1600 : {
1601 : /*
1602 : * Only operate on sides which are on the
1603 : * boundary or for which the current element's
1604 : * id is greater than its neighbor's.
1605 : * Sides get only built once.
1606 : */
1607 0 : if ((elem->neighbor_ptr(s) != nullptr) &&
1608 0 : (elem->id() > elem->neighbor_ptr(s)->id()))
1609 : {
1610 0 : const Elem & side = side_builder(*elem, s);
1611 0 : const Node & node0 = side.node_ref(0);
1612 0 : const Node & node1 = side.node_ref(1);
1613 :
1614 0 : Real node_weight = 1.;
1615 : // calculate the weight of the nodes
1616 0 : if (power > 0)
1617 : {
1618 0 : Point diff = node0-node1;
1619 0 : node_weight = std::pow(diff.norm(), power);
1620 : }
1621 :
1622 0 : const dof_id_type id0 = node0.id(), id1 = node1.id();
1623 0 : new_positions[id0].add_scaled( node1, node_weight );
1624 0 : new_positions[id1].add_scaled( node0, node_weight );
1625 0 : weight[id0] += node_weight;
1626 0 : weight[id1] += node_weight;
1627 : }
1628 : } // element neighbor loop
1629 : }
1630 : #ifdef LIBMESH_ENABLE_AMR
1631 : else // refinement_level > 0
1632 : {
1633 : /*
1634 : * Find the positions of the hanging nodes of refined elements.
1635 : * We do this by calculating their position based on the parent
1636 : * (one level less refined) element, and the embedding matrix
1637 : */
1638 :
1639 0 : const Elem * parent = elem->parent();
1640 :
1641 : /*
1642 : * find out which child I am
1643 : */
1644 0 : unsigned int c = parent->which_child_am_i(elem);
1645 : /*
1646 : *loop over the childs (that is, the current elements) nodes
1647 : */
1648 0 : for (auto nc : elem->node_index_range())
1649 : {
1650 : /*
1651 : * the new position of the node
1652 : */
1653 0 : Point point;
1654 0 : for (auto n : parent->node_index_range())
1655 : {
1656 : /*
1657 : * The value from the embedding matrix
1658 : */
1659 0 : const Real em_val = parent->embedding_matrix(c,nc,n);
1660 :
1661 0 : if (em_val != 0.)
1662 0 : point.add_scaled (parent->point(n), em_val);
1663 : }
1664 :
1665 0 : const dof_id_type id = elem->node_ptr(nc)->id();
1666 0 : new_positions[id] = point;
1667 0 : weight[id] = 1.;
1668 : }
1669 : } // if element refinement_level
1670 : #endif // #ifdef LIBMESH_ENABLE_AMR
1671 :
1672 0 : } // element loop
1673 :
1674 : /*
1675 : * finally reposition the vertex nodes
1676 : */
1677 0 : for (auto nid : make_range(mesh.n_nodes()))
1678 0 : if (!boundary_node_ids.count(nid) && weight[nid] > 0.)
1679 0 : mesh.node_ref(nid) = new_positions[nid]/weight[nid];
1680 : }
1681 :
1682 : // Now handle the additional second_order nodes by calculating
1683 : // their position based on the vertex positions
1684 : // we do a second loop over the level elements
1685 0 : for (auto & elem : as_range(mesh.level_elements_begin(refinement_level),
1686 0 : mesh.level_elements_end(refinement_level)))
1687 : {
1688 0 : const unsigned int son_begin = elem->n_vertices();
1689 0 : const unsigned int son_end = elem->n_nodes();
1690 0 : for (unsigned int n=son_begin; n<son_end; n++)
1691 : {
1692 : const unsigned int n_adjacent_vertices =
1693 0 : elem->n_second_order_adjacent_vertices(n);
1694 :
1695 0 : Point point;
1696 0 : for (unsigned int v=0; v<n_adjacent_vertices; v++)
1697 0 : point.add(elem->point( elem->second_order_adjacent_vertex(n,v) ));
1698 :
1699 0 : const dof_id_type id = elem->node_ptr(n)->id();
1700 0 : mesh.node_ref(id) = point/n_adjacent_vertices;
1701 : }
1702 0 : }
1703 : } // refinement_level loop
1704 : } // end iteration
1705 :
1706 : // We haven't changed any topology, but just changing geometry could
1707 : // have invalidated a point locator.
1708 0 : mesh.clear_point_locator();
1709 0 : }
1710 :
1711 :
1712 :
1713 3209 : void MeshTools::Modification::interpolate_surface (MeshBase & mesh,
1714 : const Surface & surface,
1715 : std::set<std::size_t> ids,
1716 : bool use_boundary_nodes)
1717 : {
1718 3209 : const bool is_serial = mesh.is_serial();
1719 192 : const processor_id_type mesh_pid = mesh.processor_id();
1720 :
1721 : // We might have to move ghost nodes on a distributed mesh if their
1722 : // owners don't see a requisite element or boundary they're on.
1723 192 : std::unordered_set<dof_id_type> moved_ghost_nodes;
1724 :
1725 180143 : auto move_node = [& moved_ghost_nodes, & surface, is_serial, mesh_pid]
1726 95311 : (Node & node) {
1727 198687 : node = surface.closest_point(node);
1728 :
1729 198687 : if (!is_serial && node.processor_id() != mesh_pid)
1730 45991 : moved_ghost_nodes.insert(node.id());
1731 12385 : };
1732 :
1733 96 : const bool no_ids = ids.empty();
1734 96 : const BoundaryInfo & boundary_info = mesh.get_boundary_info();
1735 :
1736 436360 : for (const auto & elem : mesh.active_element_ptr_range())
1737 : {
1738 225335 : if (elem->mapping_type() != LAGRANGE_MAP)
1739 0 : libmesh_not_implemented();
1740 :
1741 225335 : if (use_boundary_nodes)
1742 : {
1743 1500603 : for (auto s : elem->side_index_range())
1744 : {
1745 1275268 : if (no_ids)
1746 : {
1747 : // If we're not using boundary ids, we're
1748 : // interpolating all external and no internal
1749 : // boundaries
1750 1333676 : if (elem->neighbor_ptr(s))
1751 1165408 : continue;
1752 : }
1753 : else
1754 : {
1755 0 : if (std::none_of(ids.begin(), ids.end(),
1756 0 : [&boundary_info,elem,s](std::size_t bcid)
1757 0 : {return boundary_info.has_boundary_id(elem, s, bcid);}))
1758 0 : continue;
1759 : }
1760 :
1761 255211 : for (auto n : elem->nodes_on_side(s))
1762 207959 : move_node(elem->node_ref(n));
1763 : }
1764 : }
1765 : else
1766 : {
1767 0 : if (no_ids || ids.count(elem->subdomain_id()))
1768 0 : for (Node & node : elem->node_ref_range())
1769 0 : move_node(node);
1770 : }
1771 3017 : }
1772 :
1773 3209 : if (!is_serial)
1774 : {
1775 24 : std::map<processor_id_type, std::vector<dof_id_type>> moved_nodes_map;
1776 21284 : for (auto id : moved_ghost_nodes)
1777 : {
1778 19216 : const Node & node = mesh.node_ref(id);
1779 19216 : moved_nodes_map[node.processor_id()].push_back(node.id());
1780 : }
1781 :
1782 : auto action_functor =
1783 5645 : [& mesh, & surface]
1784 : (processor_id_type /* pid */,
1785 38800 : const std::vector<dof_id_type> & my_moved_nodes)
1786 : {
1787 24893 : for (auto id : my_moved_nodes)
1788 : {
1789 19216 : Node & node = mesh.node_ref(id);
1790 19216 : node = surface.closest_point(node);
1791 : }
1792 2076 : };
1793 :
1794 : // First get new node positions to their owners
1795 : Parallel::push_parallel_vector_data
1796 2068 : (mesh.comm(), moved_nodes_map, action_functor);
1797 :
1798 : // Then get node positions to anyone else with them ghosted
1799 2068 : SyncNodalPositions sync_object(mesh);
1800 : Parallel::sync_dofobject_data_by_id
1801 4112 : (mesh.comm(), mesh.nodes_begin(), mesh.nodes_end(),
1802 : sync_object);
1803 : }
1804 :
1805 : // We haven't changed any topology, but just changing geometry could
1806 : // have invalidated a point locator.
1807 3209 : mesh.clear_point_locator();
1808 3209 : }
1809 :
1810 :
1811 :
1812 : #ifdef LIBMESH_ENABLE_AMR
1813 2352 : void MeshTools::Modification::flatten(MeshBase & mesh)
1814 : {
1815 70 : libmesh_assert(mesh.is_prepared() || mesh.is_replicated());
1816 :
1817 : // Algorithm:
1818 : // .) For each active element in the mesh: construct a
1819 : // copy which is the same in every way *except* it is
1820 : // a level 0 element. Store the pointers to these in
1821 : // a separate vector. Save any boundary information as well.
1822 : // Delete the active element from the mesh.
1823 : // .) Loop over all (remaining) elements in the mesh, delete them.
1824 : // .) Add the level-0 copies back to the mesh
1825 :
1826 : // Temporary storage for new element pointers
1827 210 : std::vector<std::unique_ptr<Elem>> new_elements;
1828 :
1829 : // BoundaryInfo Storage for element ids, sides, and BC ids
1830 140 : std::vector<Elem *> saved_boundary_elements;
1831 140 : std::vector<boundary_id_type> saved_bc_ids;
1832 140 : std::vector<unsigned short int> saved_bc_sides;
1833 :
1834 : // Container to catch boundary ids passed back by BoundaryInfo
1835 140 : std::vector<boundary_id_type> bc_ids;
1836 :
1837 : // Reserve a reasonable amt. of space for each
1838 2352 : new_elements.reserve(mesh.n_active_elem());
1839 2352 : saved_boundary_elements.reserve(mesh.get_boundary_info().n_boundary_conds());
1840 2352 : saved_bc_ids.reserve(mesh.get_boundary_info().n_boundary_conds());
1841 2352 : saved_bc_sides.reserve(mesh.get_boundary_info().n_boundary_conds());
1842 :
1843 408622 : for (auto & elem : mesh.active_element_ptr_range())
1844 : {
1845 : // Make a new element of the same type
1846 221226 : auto copy = Elem::build(elem->type());
1847 :
1848 : // Set node pointers (they still point to nodes in the original mesh)
1849 1707972 : for (auto n : elem->node_index_range())
1850 1555544 : copy->set_node(n, elem->node_ptr(n));
1851 :
1852 : // Copy over ids
1853 211610 : copy->processor_id() = elem->processor_id();
1854 211610 : copy->subdomain_id() = elem->subdomain_id();
1855 :
1856 : // Retain the original element's ID(s) as well, otherwise
1857 : // the Mesh may try to create them for you...
1858 19232 : copy->set_id( elem->id() );
1859 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
1860 19232 : copy->set_unique_id(elem->unique_id());
1861 : #endif
1862 :
1863 : // This element could have boundary info or DistributedMesh
1864 : // remote_elem links as well. We need to save the (elem,
1865 : // side, bc_id) triples and those links
1866 1366452 : for (auto s : elem->side_index_range())
1867 : {
1868 1208122 : if (elem->neighbor_ptr(s) == remote_elem)
1869 744 : copy->set_neighbor(s, const_cast<RemoteElem *>(remote_elem));
1870 :
1871 1154842 : mesh.get_boundary_info().boundary_ids(elem, s, bc_ids);
1872 1155176 : for (const auto & bc_id : bc_ids)
1873 334 : if (bc_id != BoundaryInfo::invalid_id)
1874 : {
1875 334 : saved_boundary_elements.push_back(copy.get());
1876 334 : saved_bc_ids.push_back(bc_id);
1877 334 : saved_bc_sides.push_back(s);
1878 : }
1879 : }
1880 :
1881 : // Copy any extra element data. Since the copy hasn't been
1882 : // added to the mesh yet any allocation has to be done manually.
1883 211610 : const unsigned int nei = elem->n_extra_integers();
1884 211610 : copy->add_extra_integers(nei);
1885 211610 : for (unsigned int i=0; i != nei; ++i)
1886 0 : copy->set_extra_integer(i, elem->get_extra_integer(i));
1887 :
1888 : // Copy any mapping data.
1889 211610 : copy->set_mapping_type(elem->mapping_type());
1890 19232 : copy->set_mapping_data(elem->mapping_data());
1891 :
1892 : // We're done with this element
1893 211610 : mesh.delete_elem(elem);
1894 :
1895 : // But save the copy
1896 9616 : new_elements.push_back(std::move(copy));
1897 194590 : }
1898 :
1899 : // Make sure we saved the same number of boundary conditions
1900 : // in each vector.
1901 70 : libmesh_assert_equal_to (saved_boundary_elements.size(), saved_bc_ids.size());
1902 70 : libmesh_assert_equal_to (saved_bc_ids.size(), saved_bc_sides.size());
1903 :
1904 : // Loop again, delete any remaining elements
1905 92992 : for (auto & elem : mesh.element_ptr_range())
1906 48135 : mesh.delete_elem(elem);
1907 :
1908 : // Add the copied (now level-0) elements back to the mesh
1909 213962 : for (auto & new_elem : new_elements)
1910 : {
1911 : // Save the original ID, because the act of adding the Elem can
1912 : // change new_elem's id!
1913 9616 : dof_id_type orig_id = new_elem->id();
1914 :
1915 230842 : Elem * added_elem = mesh.add_elem(std::move(new_elem));
1916 :
1917 : // If the Elem, as it was re-added to the mesh, now has a
1918 : // different ID (this is unlikely, so it's just an assert)
1919 : // the boundary information will no longer be correct.
1920 9616 : libmesh_assert_equal_to (orig_id, added_elem->id());
1921 :
1922 : // Avoid compiler warnings in opt mode.
1923 9616 : libmesh_ignore(added_elem, orig_id);
1924 : }
1925 :
1926 : // Finally, also add back the saved boundary information
1927 2686 : for (auto e : index_range(saved_boundary_elements))
1928 358 : mesh.get_boundary_info().add_side(saved_boundary_elements[e],
1929 24 : saved_bc_sides[e],
1930 24 : saved_bc_ids[e]);
1931 :
1932 : // Trim unused and renumber nodes and elements
1933 2352 : mesh.prepare_for_use();
1934 2352 : }
1935 : #endif // #ifdef LIBMESH_ENABLE_AMR
1936 :
1937 :
1938 :
1939 3408 : void MeshTools::Modification::change_boundary_id (MeshBase & mesh,
1940 : const boundary_id_type old_id,
1941 : const boundary_id_type new_id)
1942 : {
1943 : // This is just a shim around the member implementation, now
1944 3408 : mesh.get_boundary_info().renumber_id(old_id, new_id);
1945 3408 : }
1946 :
1947 :
1948 :
1949 0 : void MeshTools::Modification::change_subdomain_id (MeshBase & mesh,
1950 : const subdomain_id_type old_id,
1951 : const subdomain_id_type new_id)
1952 : {
1953 0 : if (old_id == new_id)
1954 : {
1955 : // If the IDs are the same, this is a no-op.
1956 0 : return;
1957 : }
1958 :
1959 : Threads::parallel_for_each
1960 0 : (mesh.element_stored_range(),
1961 0 : [old_id, new_id](Elem * elem)
1962 : {
1963 0 : if (elem->subdomain_id() == old_id)
1964 0 : elem->subdomain_id() = new_id;
1965 0 : });
1966 :
1967 : // We just invalidated mesh.get_subdomain_ids(), but it might not be
1968 : // efficient to fix that here.
1969 0 : mesh.unset_has_cached_elem_data();
1970 : }
1971 :
1972 :
1973 : } // namespace libMesh
|
