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