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