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