Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2025 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/face_tri3.h"
34 : #include "libmesh/face_tri6.h"
35 : #include "libmesh/libmesh_logging.h"
36 : #include "libmesh/mesh_communication.h"
37 : #include "libmesh/mesh_modification.h"
38 : #include "libmesh/mesh_tools.h"
39 : #include "libmesh/parallel.h"
40 : #include "libmesh/remote_elem.h"
41 : #include "libmesh/enum_to_string.h"
42 : #include "libmesh/unstructured_mesh.h"
43 : #include "libmesh/elem_side_builder.h"
44 : #include "libmesh/tensor_value.h"
45 :
46 : namespace
47 : {
48 : using namespace libMesh;
49 :
50 1712 : bool split_first_diagonal(const Elem * elem,
51 : unsigned int diag_1_node_1,
52 : unsigned int diag_1_node_2,
53 : unsigned int diag_2_node_1,
54 : unsigned int diag_2_node_2)
55 : {
56 372 : return ((elem->node_id(diag_1_node_1) > elem->node_id(diag_2_node_1) &&
57 2044 : elem->node_id(diag_1_node_1) > elem->node_id(diag_2_node_2)) ||
58 1712 : (elem->node_id(diag_1_node_2) > elem->node_id(diag_2_node_1) &&
59 1768 : elem->node_id(diag_1_node_2) > elem->node_id(diag_2_node_2)));
60 : }
61 :
62 :
63 : // Return the local index of the vertex on \p elem with the highest
64 : // node id.
65 44535 : unsigned int highest_vertex_on(const Elem * elem)
66 : {
67 1792 : unsigned int highest_n = 0;
68 3584 : dof_id_type highest_n_id = elem->node_id(0);
69 356280 : for (auto n : make_range(1u, elem->n_vertices()))
70 : {
71 25088 : const dof_id_type n_id = elem->node_id(n);
72 311745 : if (n_id > highest_n_id)
73 : {
74 6928 : highest_n = n;
75 6928 : highest_n_id = n_id;
76 : }
77 : }
78 :
79 44535 : return highest_n;
80 : }
81 :
82 :
83 : static const std::array<std::array<unsigned int, 3>, 8> opposing_nodes =
84 : {{ {2,5,7},{3,4,6},{0,5,7},{1,4,6},{1,3,6},{0,2,7},{1,3,4},{0,2,5} }};
85 :
86 :
87 : // Find the highest id on these side nodes of this element
88 : std::pair<unsigned int, unsigned int>
89 201543 : split_diagonal(const Elem * elem,
90 : const std::vector<unsigned int> & nodes_on_side)
91 : {
92 6552 : libmesh_assert_equal_to(elem->type(), HEX8);
93 :
94 201543 : unsigned int highest_n = nodes_on_side.front();
95 13104 : dof_id_type highest_n_id = elem->node_id(nodes_on_side.front());
96 1007715 : for (auto n : nodes_on_side)
97 : {
98 26208 : const dof_id_type n_id = elem->node_id(n);
99 806172 : if (n_id > highest_n_id)
100 : {
101 11056 : highest_n = n;
102 11056 : highest_n_id = n_id;
103 : }
104 : }
105 :
106 409839 : for (auto n : nodes_on_side)
107 : {
108 1114271 : for (auto n2 : opposing_nodes[highest_n])
109 905975 : if (n2 == n)
110 6552 : return std::make_pair(highest_n, n2);
111 : }
112 :
113 0 : libmesh_error();
114 :
115 : return std::make_pair(libMesh::invalid_uint, libMesh::invalid_uint);
116 : }
117 :
118 :
119 : // Reconstruct a C++20 feature in C++14
120 : template <typename T>
121 : struct reversion_wrapper { T& iterable; };
122 :
123 : template <typename T>
124 4928 : auto begin (reversion_wrapper<T> w) {return std::rbegin(w.iterable);}
125 :
126 : template <typename T>
127 4928 : auto end (reversion_wrapper<T> w) {return std::rend(w.iterable);}
128 :
129 : template <typename T>
130 4928 : reversion_wrapper<T> reverse(T&& iterable) {return {iterable};}
131 :
132 : }
133 :
134 :
135 : namespace libMesh
136 : {
137 :
138 :
139 : // ------------------------------------------------------------
140 : // MeshTools::Modification functions for mesh modification
141 426 : void MeshTools::Modification::distort (MeshBase & mesh,
142 : const Real factor,
143 : const bool perturb_boundary)
144 : {
145 12 : libmesh_assert (mesh.n_nodes());
146 12 : libmesh_assert (mesh.n_elem());
147 12 : libmesh_assert ((factor >= 0.) && (factor <= 1.));
148 :
149 24 : LOG_SCOPE("distort()", "MeshTools::Modification");
150 :
151 : // If we are not perturbing boundary nodes, make a
152 : // quickly-searchable list of node ids we can check against.
153 24 : std::unordered_set<dof_id_type> boundary_node_ids;
154 426 : if (!perturb_boundary)
155 840 : boundary_node_ids = MeshTools::find_boundary_nodes (mesh);
156 :
157 : // Now calculate the minimum distance to
158 : // neighboring nodes for each node.
159 : // hmin holds these distances.
160 438 : std::vector<float> hmin (mesh.max_node_id(),
161 438 : std::numeric_limits<float>::max());
162 :
163 13674 : for (const auto & elem : mesh.active_element_ptr_range())
164 43425 : for (auto & n : elem->node_ref_range())
165 38700 : hmin[n.id()] = std::min(hmin[n.id()],
166 48831 : static_cast<float>(elem->hmin()));
167 :
168 : // Now actually move the nodes
169 : {
170 12 : const unsigned int seed = 123456;
171 :
172 : // seed the random number generator.
173 : // We'll loop from 1 to n_nodes on every processor, even those
174 : // that don't have a particular node, so that the pseudorandom
175 : // numbers will be the same everywhere.
176 426 : std::srand(seed);
177 :
178 : // If the node is on the boundary or
179 : // the node is not used by any element (hmin[n]<1.e20)
180 : // then we should not move it.
181 : // [Note: Testing for (in)equality might be wrong
182 : // (different types, namely float and double)]
183 10437 : for (auto n : make_range(mesh.max_node_id()))
184 11432 : if ((perturb_boundary || !boundary_node_ids.count(n)) && hmin[n] < 1.e20)
185 : {
186 : // the direction, random but unit normalized
187 1207 : Point dir (static_cast<Real>(std::rand())/static_cast<Real>(RAND_MAX),
188 2414 : (mesh.mesh_dimension() > 1) ? static_cast<Real>(std::rand())/static_cast<Real>(RAND_MAX) : 0.,
189 4828 : ((mesh.mesh_dimension() == 3) ? static_cast<Real>(std::rand())/static_cast<Real>(RAND_MAX) : 0.));
190 :
191 1207 : dir(0) = (dir(0)-.5)*2.;
192 : #if LIBMESH_DIM > 1
193 1207 : if (mesh.mesh_dimension() > 1)
194 1207 : dir(1) = (dir(1)-.5)*2.;
195 : #endif
196 : #if LIBMESH_DIM > 2
197 1207 : if (mesh.mesh_dimension() == 3)
198 852 : dir(2) = (dir(2)-.5)*2.;
199 : #endif
200 :
201 1207 : dir = dir.unit();
202 :
203 1207 : Node * node = mesh.query_node_ptr(n);
204 1207 : if (!node)
205 0 : continue;
206 :
207 1207 : (*node)(0) += dir(0)*factor*hmin[n];
208 : #if LIBMESH_DIM > 1
209 1207 : if (mesh.mesh_dimension() > 1)
210 1241 : (*node)(1) += dir(1)*factor*hmin[n];
211 : #endif
212 : #if LIBMESH_DIM > 2
213 1207 : if (mesh.mesh_dimension() == 3)
214 876 : (*node)(2) += dir(2)*factor*hmin[n];
215 : #endif
216 : }
217 : }
218 426 : }
219 :
220 :
221 :
222 186659 : void MeshTools::Modification::permute_elements(MeshBase & mesh)
223 : {
224 10516 : LOG_SCOPE("permute_elements()", "MeshTools::Modification");
225 :
226 : // We don't yet support doing permute() on a parent element, which
227 : // would require us to consistently permute all its children and
228 : // give them different local child numbers.
229 186659 : unsigned int n_levels = MeshTools::n_levels(mesh);
230 186659 : if (n_levels > 1)
231 0 : libmesh_error();
232 :
233 5258 : const unsigned int seed = 123456;
234 :
235 : // seed the random number generator.
236 : // We'll loop from 1 to max_elem_id on every processor, even those
237 : // that don't have a particular element, so that the pseudorandom
238 : // numbers will be the same everywhere.
239 186659 : std::srand(seed);
240 :
241 :
242 4852992 : for (auto e_id : make_range(mesh.max_elem_id()))
243 : {
244 4666333 : int my_rand = std::rand();
245 :
246 4666333 : Elem * elem = mesh.query_elem_ptr(e_id);
247 :
248 4666333 : if (!elem)
249 2768542 : continue;
250 :
251 1897791 : const unsigned int max_permutation = elem->n_permutations();
252 1897791 : if (!max_permutation)
253 4627 : continue;
254 :
255 1892506 : const unsigned int perm = my_rand % max_permutation;
256 :
257 1892506 : elem->permute(perm);
258 : }
259 186659 : }
260 :
261 :
262 2236 : void MeshTools::Modification::orient_elements(MeshBase & mesh)
263 : {
264 152 : LOG_SCOPE("orient_elements()", "MeshTools::Modification");
265 :
266 : // We don't yet support doing orient() on a parent element, which
267 : // would require us to consistently orient all its children and
268 : // give them different local child numbers.
269 2236 : unsigned int n_levels = MeshTools::n_levels(mesh);
270 2236 : if (n_levels > 1)
271 0 : libmesh_not_implemented_msg("orient_elements() does not support refined meshes");
272 :
273 76 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
274 92940 : for (auto elem : mesh.element_ptr_range())
275 45412 : elem->orient(&boundary_info);
276 2236 : }
277 :
278 :
279 :
280 184174 : void MeshTools::Modification::redistribute (MeshBase & mesh,
281 : const FunctionBase<Real> & mapfunc)
282 : {
283 5188 : libmesh_assert (mesh.n_nodes());
284 5188 : libmesh_assert (mesh.n_elem());
285 :
286 10376 : LOG_SCOPE("redistribute()", "MeshTools::Modification");
287 :
288 184174 : DenseVector<Real> output_vec(LIBMESH_DIM);
289 :
290 : // FIXME - we should thread this later.
291 189362 : std::unique_ptr<FunctionBase<Real>> myfunc = mapfunc.clone();
292 :
293 4390348 : for (auto & node : mesh.node_ptr_range())
294 : {
295 4027188 : (*myfunc)(*node, output_vec);
296 :
297 4027188 : (*node)(0) = output_vec(0);
298 : #if LIBMESH_DIM > 1
299 4027188 : (*node)(1) = output_vec(1);
300 : #endif
301 : #if LIBMESH_DIM > 2
302 4027188 : (*node)(2) = output_vec(2);
303 : #endif
304 173798 : }
305 357972 : }
306 :
307 :
308 :
309 217 : void MeshTools::Modification::translate (MeshBase & mesh,
310 : const Real xt,
311 : const Real yt,
312 : const Real zt)
313 : {
314 8 : const Point p(xt, yt, zt);
315 :
316 72638 : for (auto & node : mesh.node_ptr_range())
317 37599 : *node += p;
318 217 : }
319 :
320 :
321 : // void MeshTools::Modification::rotate2D (MeshBase & mesh,
322 : // const Real alpha)
323 : // {
324 : // libmesh_assert_not_equal_to (mesh.mesh_dimension(), 1);
325 :
326 : // const Real pi = std::acos(-1);
327 : // const Real a = alpha/180.*pi;
328 : // for (unsigned int n=0; n<mesh.n_nodes(); n++)
329 : // {
330 : // const Point p = mesh.node_ref(n);
331 : // const Real x = p(0);
332 : // const Real y = p(1);
333 : // const Real z = p(2);
334 : // mesh.node_ref(n) = Point(std::cos(a)*x - std::sin(a)*y,
335 : // std::sin(a)*x + std::cos(a)*y,
336 : // z);
337 : // }
338 :
339 : // }
340 :
341 :
342 :
343 : RealTensorValue
344 162519 : MeshTools::Modification::rotate (MeshBase & mesh,
345 : const Real phi,
346 : const Real theta,
347 : const Real psi)
348 : {
349 : #if LIBMESH_DIM == 3
350 162519 : const auto R = RealTensorValue::intrinsic_rotation_matrix(phi, theta, psi);
351 :
352 162519 : if (theta)
353 88253 : mesh.set_spatial_dimension(3);
354 :
355 10235306 : for (auto & node : mesh.node_ptr_range())
356 : {
357 5267113 : Point & pt = *node;
358 5267113 : pt = R * pt;
359 153363 : }
360 :
361 162519 : return R;
362 :
363 : #else
364 : libmesh_ignore(mesh, phi, theta, psi);
365 : libmesh_error_msg("MeshTools::Modification::rotate() requires libMesh to be compiled with LIBMESH_DIM==3");
366 : // We'll never get here
367 : return RealTensorValue();
368 : #endif
369 : }
370 :
371 :
372 71 : void MeshTools::Modification::scale (MeshBase & mesh,
373 : const Real xs,
374 : const Real ys,
375 : const Real zs)
376 : {
377 2 : const Real x_scale = xs;
378 2 : Real y_scale = ys;
379 2 : Real z_scale = zs;
380 :
381 71 : if (ys == 0.)
382 : {
383 0 : libmesh_assert_equal_to (zs, 0.);
384 :
385 0 : y_scale = z_scale = x_scale;
386 : }
387 :
388 : // Scale the x coordinate in all dimensions
389 495 : for (auto & node : mesh.node_ptr_range())
390 422 : (*node)(0) *= x_scale;
391 :
392 : // Only scale the y coordinate in 2 and 3D
393 : if (LIBMESH_DIM < 2)
394 : return;
395 :
396 495 : for (auto & node : mesh.node_ptr_range())
397 422 : (*node)(1) *= y_scale;
398 :
399 : // Only scale the z coordinate in 3D
400 : if (LIBMESH_DIM < 3)
401 : return;
402 :
403 495 : for (auto & node : mesh.node_ptr_range())
404 422 : (*node)(2) *= z_scale;
405 : }
406 :
407 :
408 :
409 2272 : void MeshTools::Modification::all_tri (MeshBase & mesh)
410 : {
411 2272 : if (!mesh.is_replicated() && !mesh.is_prepared())
412 0 : mesh.prepare_for_use();
413 :
414 : // The number of elements in the original mesh before any additions
415 : // or deletions.
416 2272 : const dof_id_type n_orig_elem = mesh.n_elem();
417 2272 : const dof_id_type max_orig_id = mesh.max_elem_id();
418 :
419 : // We store pointers to the newly created elements in a vector
420 : // until they are ready to be added to the mesh. This is because
421 : // adding new elements on the fly can cause reallocation and invalidation
422 : // of existing mesh element_iterators.
423 192 : std::vector<std::unique_ptr<Elem>> new_elements;
424 :
425 64 : unsigned int max_subelems = 1; // in 1D nothing needs to change
426 2272 : if (mesh.mesh_dimension() == 2) // in 2D quads can split into 2 tris
427 52 : max_subelems = 2;
428 2272 : if (mesh.mesh_dimension() == 3) // in 3D hexes can split into 6 tets
429 12 : max_subelems = 6;
430 :
431 2272 : new_elements.reserve (max_subelems*n_orig_elem);
432 :
433 : // If the original mesh has *side* boundary data, we carry that over
434 : // to the new mesh with triangular elements. We currently only
435 : // support bringing over side-based BCs to the all-tri mesh, but
436 : // that could probably be extended to node and edge-based BCs as
437 : // well.
438 2272 : const bool mesh_has_boundary_data = (mesh.get_boundary_info().n_boundary_conds() > 0);
439 :
440 : // Temporary vectors to store the new boundary element pointers, side numbers, and boundary ids
441 128 : std::vector<Elem *> new_bndry_elements;
442 128 : std::vector<unsigned short int> new_bndry_sides;
443 128 : std::vector<boundary_id_type> new_bndry_ids;
444 :
445 : // We may need to add new points if we run into a 1.5th order
446 : // element; if we do that on a DistributedMesh in a ghost element then
447 : // we will need to fix their ids / unique_ids
448 2272 : bool added_new_ghost_point = false;
449 :
450 : // Iterate over the elements, splitting:
451 : // QUADs into pairs of conforming triangles
452 : // PYRAMIDs into pairs of conforming tets,
453 : // PRISMs into triplets of conforming tets, and
454 : // HEXs into quintets or sextets of conforming tets.
455 : // We split on the shortest diagonal to give us better
456 : // triangle quality in 2D, and we split based on node ids
457 : // to guarantee consistency in 3D.
458 :
459 : // FIXME: This algorithm does not work on refined grids!
460 : {
461 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
462 2272 : unique_id_type max_unique_id = mesh.parallel_max_unique_id();
463 : #endif
464 :
465 : // For avoiding extraneous allocation when building side elements
466 2272 : std::unique_ptr<const Elem> elem_side, subside_elem;
467 :
468 151078 : for (auto & elem : mesh.element_ptr_range())
469 : {
470 76325 : const ElemType etype = elem->type();
471 :
472 : // all_tri currently only works on coarse meshes
473 79351 : if (elem->parent())
474 0 : libmesh_not_implemented_msg("Cannot convert a refined element into simplices\n");
475 :
476 : // The new elements we will split the quad into.
477 : // In 3D we may need as many as 6 tets per hex
478 73299 : std::array<std::unique_ptr<Elem>, 6> subelem {};
479 :
480 76325 : switch (etype)
481 : {
482 9323 : case QUAD4:
483 : {
484 17966 : subelem[0] = Elem::build(TRI3);
485 17966 : subelem[1] = Elem::build(TRI3);
486 :
487 : // Check for possible edge swap
488 9663 : if ((elem->point(0) - elem->point(2)).norm() <
489 9323 : (elem->point(1) - elem->point(3)).norm())
490 : {
491 3118 : subelem[0]->set_node(0, elem->node_ptr(0));
492 3118 : subelem[0]->set_node(1, elem->node_ptr(1));
493 3118 : subelem[0]->set_node(2, elem->node_ptr(2));
494 :
495 3118 : subelem[1]->set_node(0, elem->node_ptr(0));
496 3118 : subelem[1]->set_node(1, elem->node_ptr(2));
497 3118 : subelem[1]->set_node(2, elem->node_ptr(3));
498 : }
499 :
500 : else
501 : {
502 6545 : subelem[0]->set_node(0, elem->node_ptr(0));
503 6545 : subelem[0]->set_node(1, elem->node_ptr(1));
504 6545 : subelem[0]->set_node(2, elem->node_ptr(3));
505 :
506 6545 : subelem[1]->set_node(0, elem->node_ptr(1));
507 6545 : subelem[1]->set_node(1, elem->node_ptr(2));
508 6545 : subelem[1]->set_node(2, elem->node_ptr(3));
509 : }
510 :
511 :
512 340 : break;
513 : }
514 :
515 142 : case QUAD8:
516 : {
517 146 : if (elem->processor_id() != mesh.processor_id())
518 118 : added_new_ghost_point = true;
519 :
520 276 : subelem[0] = Elem::build(TRI6);
521 276 : subelem[1] = Elem::build(TRI6);
522 :
523 : // Add a new node at the center (vertex average) of the element.
524 150 : Node * new_node = mesh.add_point((mesh.point(elem->node_id(0)) +
525 150 : mesh.point(elem->node_id(1)) +
526 150 : mesh.point(elem->node_id(2)) +
527 146 : mesh.point(elem->node_id(3)))/4,
528 : DofObject::invalid_id,
529 142 : elem->processor_id());
530 :
531 : // Check for possible edge swap
532 146 : if ((elem->point(0) - elem->point(2)).norm() <
533 142 : (elem->point(1) - elem->point(3)).norm())
534 : {
535 0 : subelem[0]->set_node(0, elem->node_ptr(0));
536 0 : subelem[0]->set_node(1, elem->node_ptr(1));
537 0 : subelem[0]->set_node(2, elem->node_ptr(2));
538 0 : subelem[0]->set_node(3, elem->node_ptr(4));
539 0 : subelem[0]->set_node(4, elem->node_ptr(5));
540 0 : subelem[0]->set_node(5, new_node);
541 :
542 0 : subelem[1]->set_node(0, elem->node_ptr(0));
543 0 : subelem[1]->set_node(1, elem->node_ptr(2));
544 0 : subelem[1]->set_node(2, elem->node_ptr(3));
545 0 : subelem[1]->set_node(3, new_node);
546 0 : subelem[1]->set_node(4, elem->node_ptr(6));
547 0 : subelem[1]->set_node(5, elem->node_ptr(7));
548 :
549 : }
550 :
551 : else
552 : {
553 146 : subelem[0]->set_node(0, elem->node_ptr(3));
554 146 : subelem[0]->set_node(1, elem->node_ptr(0));
555 146 : subelem[0]->set_node(2, elem->node_ptr(1));
556 146 : subelem[0]->set_node(3, elem->node_ptr(7));
557 146 : subelem[0]->set_node(4, elem->node_ptr(4));
558 142 : subelem[0]->set_node(5, new_node);
559 :
560 146 : subelem[1]->set_node(0, elem->node_ptr(1));
561 146 : subelem[1]->set_node(1, elem->node_ptr(2));
562 146 : subelem[1]->set_node(2, elem->node_ptr(3));
563 146 : subelem[1]->set_node(3, elem->node_ptr(5));
564 146 : subelem[1]->set_node(4, elem->node_ptr(6));
565 142 : subelem[1]->set_node(5, new_node);
566 : }
567 :
568 4 : break;
569 : }
570 :
571 422 : case QUAD9:
572 : {
573 804 : subelem[0] = Elem::build(TRI6);
574 804 : subelem[1] = Elem::build(TRI6);
575 :
576 : // Check for possible edge swap
577 442 : if ((elem->point(0) - elem->point(2)).norm() <
578 422 : (elem->point(1) - elem->point(3)).norm())
579 : {
580 0 : subelem[0]->set_node(0, elem->node_ptr(0));
581 0 : subelem[0]->set_node(1, elem->node_ptr(1));
582 0 : subelem[0]->set_node(2, elem->node_ptr(2));
583 0 : subelem[0]->set_node(3, elem->node_ptr(4));
584 0 : subelem[0]->set_node(4, elem->node_ptr(5));
585 0 : subelem[0]->set_node(5, elem->node_ptr(8));
586 :
587 0 : subelem[1]->set_node(0, elem->node_ptr(0));
588 0 : subelem[1]->set_node(1, elem->node_ptr(2));
589 0 : subelem[1]->set_node(2, elem->node_ptr(3));
590 0 : subelem[1]->set_node(3, elem->node_ptr(8));
591 0 : subelem[1]->set_node(4, elem->node_ptr(6));
592 0 : subelem[1]->set_node(5, elem->node_ptr(7));
593 : }
594 :
595 : else
596 : {
597 442 : subelem[0]->set_node(0, elem->node_ptr(0));
598 442 : subelem[0]->set_node(1, elem->node_ptr(1));
599 442 : subelem[0]->set_node(2, elem->node_ptr(3));
600 442 : subelem[0]->set_node(3, elem->node_ptr(4));
601 442 : subelem[0]->set_node(4, elem->node_ptr(8));
602 442 : subelem[0]->set_node(5, elem->node_ptr(7));
603 :
604 442 : subelem[1]->set_node(0, elem->node_ptr(1));
605 442 : subelem[1]->set_node(1, elem->node_ptr(2));
606 442 : subelem[1]->set_node(2, elem->node_ptr(3));
607 442 : subelem[1]->set_node(3, elem->node_ptr(5));
608 442 : subelem[1]->set_node(4, elem->node_ptr(6));
609 442 : subelem[1]->set_node(5, elem->node_ptr(8));
610 : }
611 :
612 20 : break;
613 : }
614 :
615 1792 : case HEX8:
616 : {
617 1792 : BoundaryInfo & boundary_info = mesh.get_boundary_info();
618 :
619 : // Hexes all split into six tetrahedra
620 85486 : subelem[0] = Elem::build(TET4);
621 85486 : subelem[1] = Elem::build(TET4);
622 85486 : subelem[2] = Elem::build(TET4);
623 85486 : subelem[3] = Elem::build(TET4);
624 85486 : subelem[4] = Elem::build(TET4);
625 85486 : subelem[5] = Elem::build(TET4);
626 :
627 : // On faces, we choose the node with the highest
628 : // global id, and we split on the diagonal which
629 : // includes that node. This ensures that (even in
630 : // parallel, even on distributed meshes) the same
631 : // diagonal split will be chosen for elements on either
632 : // side of the same quad face.
633 44535 : const unsigned int highest_n = highest_vertex_on(elem);
634 :
635 : // opposing_node[n] is the local node number of the node
636 : // on the farthest corner of a hex8 from local node n
637 : static const std::array<unsigned int, 8> opposing_node =
638 : {6, 7, 4, 5, 2, 3, 0, 1};
639 :
640 : static const std::vector<std::vector<unsigned int>> sides_opposing_highest =
641 45170 : {{2,3,5},{3,4,5},{1,4,5},{1,2,5},{0,2,3},{0,3,4},{0,1,4},{0,1,2}};
642 : static const std::vector<std::vector<unsigned int>> nodes_neighboring_highest =
643 45170 : {{1,3,4},{0,2,5},{1,3,6},{0,2,7},{0,5,7},{1,4,6},{2,5,7},{3,4,6}};
644 :
645 : // Start by looking in three directions away from the
646 : // highest-id node. In each direction there will be two
647 : // different possibilities for the split depending on
648 : // how the opposing face nodes are numbered.
649 : //
650 : // This is tricky enough that I'm not going to worry
651 : // about manually keeping tets oriented; we'll just call
652 : // orient() on each as we go.
653 :
654 1792 : unsigned int next_subelem = 0;
655 179932 : for (auto side : sides_opposing_highest[highest_n])
656 : {
657 : const std::vector<unsigned int> nodes_on_side =
658 138981 : elem->nodes_on_side(side);
659 :
660 133605 : auto [dn, dn2] = split_diagonal(elem, nodes_on_side);
661 :
662 5376 : unsigned int split_on_neighbor = false;
663 342347 : for (auto n : nodes_neighboring_highest[highest_n])
664 302750 : if (dn == n || dn2 == n)
665 : {
666 4928 : split_on_neighbor = true;
667 4928 : break;
668 : }
669 :
670 : // Add one or two elements for each opposing side,
671 : // depending on whether the diagonal split there
672 : // connects to the neighboring diagonal split or
673 : // not.
674 133605 : if (split_on_neighbor)
675 : {
676 109240 : subelem[next_subelem]->set_node(0, elem->node_ptr(highest_n));
677 104312 : subelem[next_subelem]->set_node(1, elem->node_ptr(dn));
678 104312 : subelem[next_subelem]->set_node(2, elem->node_ptr(dn2));
679 150357 : for (auto n : nodes_on_side)
680 150357 : if (n != dn && n != dn2)
681 : {
682 104312 : subelem[next_subelem]->set_node(3, elem->node_ptr(n));
683 4928 : break;
684 : }
685 94456 : subelem[next_subelem]->orient(&boundary_info);
686 99384 : ++next_subelem;
687 :
688 109240 : subelem[next_subelem]->set_node(0, elem->node_ptr(highest_n));
689 104312 : subelem[next_subelem]->set_node(1, elem->node_ptr(dn));
690 104312 : subelem[next_subelem]->set_node(2, elem->node_ptr(dn2));
691 147795 : for (auto n : reverse(nodes_on_side))
692 147795 : if (n != dn && n != dn2)
693 : {
694 104312 : subelem[next_subelem]->set_node(3, elem->node_ptr(n));
695 4928 : break;
696 : }
697 94456 : subelem[next_subelem]->orient(&boundary_info);
698 99384 : ++next_subelem;
699 : }
700 : else
701 : {
702 35117 : subelem[next_subelem]->set_node(0, elem->node_ptr(highest_n));
703 34669 : subelem[next_subelem]->set_node(1, elem->node_ptr(dn));
704 34669 : subelem[next_subelem]->set_node(2, elem->node_ptr(dn2));
705 101845 : for (auto n : nodes_on_side)
706 339087 : for (auto n2 : nodes_neighboring_highest[highest_n])
707 269995 : if (n == n2)
708 : {
709 34669 : subelem[next_subelem]->set_node(3, elem->node_ptr(n));
710 33773 : goto break_both_loops;
711 : }
712 :
713 34221 : break_both_loops:
714 33773 : subelem[next_subelem]->orient(&boundary_info);
715 34221 : ++next_subelem;
716 : }
717 : }
718 :
719 : // At this point we've created between 3 and 6 tets.
720 : // What's left to do depends on how many.
721 :
722 : // If we just chopped off three vertices into three
723 : // tets, then the best way to split this hex would be
724 : // the symmetric five-split. Chop off the opposing
725 : // vertex too, and then the remaining interior is our
726 : // final tet.
727 44535 : if (next_subelem == 3)
728 : {
729 2525 : subelem[next_subelem]->set_node(0, elem->node_ptr(opposing_nodes[highest_n][0]));
730 2525 : subelem[next_subelem]->set_node(1, elem->node_ptr(opposing_nodes[highest_n][1]));
731 2525 : subelem[next_subelem]->set_node(2, elem->node_ptr(opposing_nodes[highest_n][2]));
732 2525 : subelem[next_subelem]->set_node(3, elem->node_ptr(opposing_node[highest_n]));
733 2525 : subelem[next_subelem]->orient(&boundary_info);
734 0 : ++next_subelem;
735 :
736 2525 : subelem[next_subelem]->set_node(0, elem->node_ptr(opposing_nodes[highest_n][0]));
737 2525 : subelem[next_subelem]->set_node(1, elem->node_ptr(opposing_nodes[highest_n][1]));
738 2525 : subelem[next_subelem]->set_node(2, elem->node_ptr(opposing_nodes[highest_n][2]));
739 2525 : subelem[next_subelem]->set_node(3, elem->node_ptr(highest_n));
740 2525 : subelem[next_subelem]->orient(&boundary_info);
741 0 : ++next_subelem;
742 :
743 : // We don't need the 6th tet after all
744 0 : subelem[next_subelem].reset();
745 0 : ++next_subelem;
746 : }
747 :
748 : // If we just chopped off one (or two) vertices into
749 : // tets, then the remaining gap is best (or only) filled
750 : // by pairing another tet with each.
751 44535 : if (next_subelem == 4 ||
752 : next_subelem == 5)
753 : {
754 90976 : for (auto side : sides_opposing_highest[highest_n])
755 : {
756 : const std::vector<unsigned int> nodes_on_side =
757 69114 : elem->nodes_on_side(side);
758 :
759 67938 : auto [dn, dn2] = split_diagonal(elem, nodes_on_side);
760 :
761 1176 : unsigned int split_on_neighbor = false;
762 191663 : for (auto n : nodes_neighboring_highest[highest_n])
763 163841 : if (dn == n || dn2 == n)
764 : {
765 728 : split_on_neighbor = true;
766 728 : break;
767 : }
768 :
769 : // The two !split_on_neighbor sides are where we
770 : // need the two remaining tets
771 67938 : if (!split_on_neighbor)
772 : {
773 27542 : subelem[next_subelem]->set_node(0, elem->node_ptr(highest_n));
774 27094 : subelem[next_subelem]->set_node(1, elem->node_ptr(dn));
775 27094 : subelem[next_subelem]->set_node(2, elem->node_ptr(dn2));
776 27094 : subelem[next_subelem]->set_node(3, elem->node_ptr(opposing_node[highest_n]));
777 26198 : subelem[next_subelem]->orient(&boundary_info);
778 26646 : ++next_subelem;
779 : }
780 : }
781 : }
782 :
783 : // Whether we got there by creating six tets from the
784 : // first for loop or by patching up the split afterward,
785 : // we should have considered six tets (possibly
786 : // including one deleted one...) at this point.
787 1792 : libmesh_assert(next_subelem == 6);
788 :
789 1792 : break;
790 : }
791 :
792 142 : case PRISM6:
793 : {
794 : // Prisms all split into three tetrahedra
795 276 : subelem[0] = Elem::build(TET4);
796 276 : subelem[1] = Elem::build(TET4);
797 276 : subelem[2] = Elem::build(TET4);
798 :
799 : // Triangular faces are not split.
800 :
801 : // On quad faces, we choose the node with the highest
802 : // global id, and we split on the diagonal which
803 : // includes that node. This ensures that (even in
804 : // parallel, even on distributed meshes) the same
805 : // diagonal split will be chosen for elements on either
806 : // side of the same quad face. It also ensures that we
807 : // always have a mix of "clockwise" and
808 : // "counterclockwise" split faces (two of one and one
809 : // of the other on each prism; this is useful since the
810 : // alternative all-clockwise or all-counterclockwise
811 : // face splittings can't be turned into tets without
812 : // adding more nodes
813 :
814 : // Split on 0-4 diagonal
815 142 : if (split_first_diagonal(elem, 0,4, 1,3))
816 : {
817 : // Split on 0-5 diagonal
818 142 : if (split_first_diagonal(elem, 0,5, 2,3))
819 : {
820 : // Split on 1-5 diagonal
821 142 : if (split_first_diagonal(elem, 1,5, 2,4))
822 : {
823 73 : subelem[0]->set_node(0, elem->node_ptr(0));
824 73 : subelem[0]->set_node(1, elem->node_ptr(4));
825 73 : subelem[0]->set_node(2, elem->node_ptr(5));
826 73 : subelem[0]->set_node(3, elem->node_ptr(3));
827 :
828 73 : subelem[1]->set_node(0, elem->node_ptr(0));
829 73 : subelem[1]->set_node(1, elem->node_ptr(4));
830 73 : subelem[1]->set_node(2, elem->node_ptr(1));
831 73 : subelem[1]->set_node(3, elem->node_ptr(5));
832 :
833 73 : subelem[2]->set_node(0, elem->node_ptr(0));
834 73 : subelem[2]->set_node(1, elem->node_ptr(1));
835 73 : subelem[2]->set_node(2, elem->node_ptr(2));
836 73 : subelem[2]->set_node(3, elem->node_ptr(5));
837 : }
838 : else // Split on 2-4 diagonal
839 : {
840 2 : libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
841 :
842 73 : subelem[0]->set_node(0, elem->node_ptr(0));
843 73 : subelem[0]->set_node(1, elem->node_ptr(4));
844 73 : subelem[0]->set_node(2, elem->node_ptr(5));
845 73 : subelem[0]->set_node(3, elem->node_ptr(3));
846 :
847 73 : subelem[1]->set_node(0, elem->node_ptr(0));
848 73 : subelem[1]->set_node(1, elem->node_ptr(4));
849 73 : subelem[1]->set_node(2, elem->node_ptr(2));
850 73 : subelem[1]->set_node(3, elem->node_ptr(5));
851 :
852 73 : subelem[2]->set_node(0, elem->node_ptr(0));
853 73 : subelem[2]->set_node(1, elem->node_ptr(1));
854 73 : subelem[2]->set_node(2, elem->node_ptr(2));
855 73 : subelem[2]->set_node(3, elem->node_ptr(4));
856 : }
857 : }
858 : else // Split on 2-3 diagonal
859 : {
860 0 : libmesh_assert (split_first_diagonal(elem, 2,3, 0,5));
861 :
862 : // 0-4 and 2-3 split implies 2-4 split
863 0 : libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
864 :
865 0 : subelem[0]->set_node(0, elem->node_ptr(0));
866 0 : subelem[0]->set_node(1, elem->node_ptr(4));
867 0 : subelem[0]->set_node(2, elem->node_ptr(2));
868 0 : subelem[0]->set_node(3, elem->node_ptr(3));
869 :
870 0 : subelem[1]->set_node(0, elem->node_ptr(3));
871 0 : subelem[1]->set_node(1, elem->node_ptr(4));
872 0 : subelem[1]->set_node(2, elem->node_ptr(2));
873 0 : subelem[1]->set_node(3, elem->node_ptr(5));
874 :
875 0 : subelem[2]->set_node(0, elem->node_ptr(0));
876 0 : subelem[2]->set_node(1, elem->node_ptr(1));
877 0 : subelem[2]->set_node(2, elem->node_ptr(2));
878 0 : subelem[2]->set_node(3, elem->node_ptr(4));
879 : }
880 : }
881 : else // Split on 1-3 diagonal
882 : {
883 0 : libmesh_assert (split_first_diagonal(elem, 1,3, 0,4));
884 :
885 : // Split on 0-5 diagonal
886 0 : if (split_first_diagonal(elem, 0,5, 2,3))
887 : {
888 : // 1-3 and 0-5 split implies 1-5 split
889 0 : libmesh_assert (split_first_diagonal(elem, 1,5, 2,4));
890 :
891 0 : subelem[0]->set_node(0, elem->node_ptr(1));
892 0 : subelem[0]->set_node(1, elem->node_ptr(3));
893 0 : subelem[0]->set_node(2, elem->node_ptr(4));
894 0 : subelem[0]->set_node(3, elem->node_ptr(5));
895 :
896 0 : subelem[1]->set_node(0, elem->node_ptr(1));
897 0 : subelem[1]->set_node(1, elem->node_ptr(0));
898 0 : subelem[1]->set_node(2, elem->node_ptr(3));
899 0 : subelem[1]->set_node(3, elem->node_ptr(5));
900 :
901 0 : subelem[2]->set_node(0, elem->node_ptr(0));
902 0 : subelem[2]->set_node(1, elem->node_ptr(1));
903 0 : subelem[2]->set_node(2, elem->node_ptr(2));
904 0 : subelem[2]->set_node(3, elem->node_ptr(5));
905 : }
906 : else // Split on 2-3 diagonal
907 : {
908 0 : libmesh_assert (split_first_diagonal(elem, 2,3, 0,5));
909 :
910 : // Split on 1-5 diagonal
911 0 : if (split_first_diagonal(elem, 1,5, 2,4))
912 : {
913 0 : subelem[0]->set_node(0, elem->node_ptr(0));
914 0 : subelem[0]->set_node(1, elem->node_ptr(1));
915 0 : subelem[0]->set_node(2, elem->node_ptr(2));
916 0 : subelem[0]->set_node(3, elem->node_ptr(3));
917 :
918 0 : subelem[1]->set_node(0, elem->node_ptr(3));
919 0 : subelem[1]->set_node(1, elem->node_ptr(1));
920 0 : subelem[1]->set_node(2, elem->node_ptr(2));
921 0 : subelem[1]->set_node(3, elem->node_ptr(5));
922 :
923 0 : subelem[2]->set_node(0, elem->node_ptr(1));
924 0 : subelem[2]->set_node(1, elem->node_ptr(3));
925 0 : subelem[2]->set_node(2, elem->node_ptr(4));
926 0 : subelem[2]->set_node(3, elem->node_ptr(5));
927 : }
928 : else // Split on 2-4 diagonal
929 : {
930 0 : libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
931 :
932 0 : subelem[0]->set_node(0, elem->node_ptr(0));
933 0 : subelem[0]->set_node(1, elem->node_ptr(1));
934 0 : subelem[0]->set_node(2, elem->node_ptr(2));
935 0 : subelem[0]->set_node(3, elem->node_ptr(3));
936 :
937 0 : subelem[1]->set_node(0, elem->node_ptr(2));
938 0 : subelem[1]->set_node(1, elem->node_ptr(3));
939 0 : subelem[1]->set_node(2, elem->node_ptr(4));
940 0 : subelem[1]->set_node(3, elem->node_ptr(5));
941 :
942 0 : subelem[2]->set_node(0, elem->node_ptr(3));
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(4));
946 : }
947 : }
948 : }
949 :
950 4 : break;
951 : }
952 :
953 426 : case PRISM20:
954 : case PRISM21:
955 : libmesh_experimental(); // We should upgrade this to TET14...
956 : libmesh_fallthrough();
957 : case PRISM18:
958 : {
959 828 : subelem[0] = Elem::build(TET10);
960 828 : subelem[1] = Elem::build(TET10);
961 828 : subelem[2] = Elem::build(TET10);
962 :
963 : // Split on 0-4 diagonal
964 426 : if (split_first_diagonal(elem, 0,4, 1,3))
965 : {
966 : // Split on 0-5 diagonal
967 426 : if (split_first_diagonal(elem, 0,5, 2,3))
968 : {
969 : // Split on 1-5 diagonal
970 426 : if (split_first_diagonal(elem, 1,5, 2,4))
971 : {
972 219 : subelem[0]->set_node(0, elem->node_ptr(0));
973 219 : subelem[0]->set_node(1, elem->node_ptr(4));
974 219 : subelem[0]->set_node(2, elem->node_ptr(5));
975 219 : subelem[0]->set_node(3, elem->node_ptr(3));
976 :
977 219 : subelem[0]->set_node(4, elem->node_ptr(15));
978 219 : subelem[0]->set_node(5, elem->node_ptr(13));
979 219 : subelem[0]->set_node(6, elem->node_ptr(17));
980 219 : subelem[0]->set_node(7, elem->node_ptr(9));
981 219 : subelem[0]->set_node(8, elem->node_ptr(12));
982 219 : subelem[0]->set_node(9, elem->node_ptr(14));
983 :
984 219 : subelem[1]->set_node(0, elem->node_ptr(0));
985 219 : subelem[1]->set_node(1, elem->node_ptr(4));
986 219 : subelem[1]->set_node(2, elem->node_ptr(1));
987 219 : subelem[1]->set_node(3, elem->node_ptr(5));
988 :
989 219 : subelem[1]->set_node(4, elem->node_ptr(15));
990 219 : subelem[1]->set_node(5, elem->node_ptr(10));
991 219 : subelem[1]->set_node(6, elem->node_ptr(6));
992 219 : subelem[1]->set_node(7, elem->node_ptr(17));
993 219 : subelem[1]->set_node(8, elem->node_ptr(13));
994 219 : subelem[1]->set_node(9, elem->node_ptr(16));
995 :
996 219 : subelem[2]->set_node(0, elem->node_ptr(0));
997 219 : subelem[2]->set_node(1, elem->node_ptr(1));
998 219 : subelem[2]->set_node(2, elem->node_ptr(2));
999 219 : subelem[2]->set_node(3, elem->node_ptr(5));
1000 :
1001 219 : subelem[2]->set_node(4, elem->node_ptr(6));
1002 219 : subelem[2]->set_node(5, elem->node_ptr(7));
1003 219 : subelem[2]->set_node(6, elem->node_ptr(8));
1004 219 : subelem[2]->set_node(7, elem->node_ptr(17));
1005 219 : subelem[2]->set_node(8, elem->node_ptr(16));
1006 219 : subelem[2]->set_node(9, elem->node_ptr(11));
1007 : }
1008 : else // Split on 2-4 diagonal
1009 : {
1010 6 : libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
1011 :
1012 219 : subelem[0]->set_node(0, elem->node_ptr(0));
1013 219 : subelem[0]->set_node(1, elem->node_ptr(4));
1014 219 : subelem[0]->set_node(2, elem->node_ptr(5));
1015 219 : subelem[0]->set_node(3, elem->node_ptr(3));
1016 :
1017 219 : subelem[0]->set_node(4, elem->node_ptr(15));
1018 219 : subelem[0]->set_node(5, elem->node_ptr(13));
1019 219 : subelem[0]->set_node(6, elem->node_ptr(17));
1020 219 : subelem[0]->set_node(7, elem->node_ptr(9));
1021 219 : subelem[0]->set_node(8, elem->node_ptr(12));
1022 219 : subelem[0]->set_node(9, elem->node_ptr(14));
1023 :
1024 219 : subelem[1]->set_node(0, elem->node_ptr(0));
1025 219 : subelem[1]->set_node(1, elem->node_ptr(4));
1026 219 : subelem[1]->set_node(2, elem->node_ptr(2));
1027 219 : subelem[1]->set_node(3, elem->node_ptr(5));
1028 :
1029 219 : subelem[1]->set_node(4, elem->node_ptr(15));
1030 219 : subelem[1]->set_node(5, elem->node_ptr(16));
1031 219 : subelem[1]->set_node(6, elem->node_ptr(8));
1032 219 : subelem[1]->set_node(7, elem->node_ptr(17));
1033 219 : subelem[1]->set_node(8, elem->node_ptr(13));
1034 219 : subelem[1]->set_node(9, elem->node_ptr(11));
1035 :
1036 219 : subelem[2]->set_node(0, elem->node_ptr(0));
1037 219 : subelem[2]->set_node(1, elem->node_ptr(1));
1038 219 : subelem[2]->set_node(2, elem->node_ptr(2));
1039 219 : subelem[2]->set_node(3, elem->node_ptr(4));
1040 :
1041 219 : subelem[2]->set_node(4, elem->node_ptr(6));
1042 219 : subelem[2]->set_node(5, elem->node_ptr(7));
1043 219 : subelem[2]->set_node(6, elem->node_ptr(8));
1044 219 : subelem[2]->set_node(7, elem->node_ptr(15));
1045 219 : subelem[2]->set_node(8, elem->node_ptr(10));
1046 219 : subelem[2]->set_node(9, elem->node_ptr(16));
1047 : }
1048 : }
1049 : else // Split on 2-3 diagonal
1050 : {
1051 0 : libmesh_assert (split_first_diagonal(elem, 2,3, 0,5));
1052 :
1053 : // 0-4 and 2-3 split implies 2-4 split
1054 0 : libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
1055 :
1056 0 : subelem[0]->set_node(0, elem->node_ptr(0));
1057 0 : subelem[0]->set_node(1, elem->node_ptr(4));
1058 0 : subelem[0]->set_node(2, elem->node_ptr(2));
1059 0 : subelem[0]->set_node(3, elem->node_ptr(3));
1060 :
1061 0 : subelem[0]->set_node(4, elem->node_ptr(15));
1062 0 : subelem[0]->set_node(5, elem->node_ptr(16));
1063 0 : subelem[0]->set_node(6, elem->node_ptr(8));
1064 0 : subelem[0]->set_node(7, elem->node_ptr(9));
1065 0 : subelem[0]->set_node(8, elem->node_ptr(12));
1066 0 : subelem[0]->set_node(9, elem->node_ptr(17));
1067 :
1068 0 : subelem[1]->set_node(0, elem->node_ptr(3));
1069 0 : subelem[1]->set_node(1, elem->node_ptr(4));
1070 0 : subelem[1]->set_node(2, elem->node_ptr(2));
1071 0 : subelem[1]->set_node(3, elem->node_ptr(5));
1072 :
1073 0 : subelem[1]->set_node(4, elem->node_ptr(12));
1074 0 : subelem[1]->set_node(5, elem->node_ptr(16));
1075 0 : subelem[1]->set_node(6, elem->node_ptr(17));
1076 0 : subelem[1]->set_node(7, elem->node_ptr(14));
1077 0 : subelem[1]->set_node(8, elem->node_ptr(13));
1078 0 : subelem[1]->set_node(9, elem->node_ptr(11));
1079 :
1080 0 : subelem[2]->set_node(0, elem->node_ptr(0));
1081 0 : subelem[2]->set_node(1, elem->node_ptr(1));
1082 0 : subelem[2]->set_node(2, elem->node_ptr(2));
1083 0 : subelem[2]->set_node(3, elem->node_ptr(4));
1084 :
1085 0 : subelem[2]->set_node(4, elem->node_ptr(6));
1086 0 : subelem[2]->set_node(5, elem->node_ptr(7));
1087 0 : subelem[2]->set_node(6, elem->node_ptr(8));
1088 0 : subelem[2]->set_node(7, elem->node_ptr(15));
1089 0 : subelem[2]->set_node(8, elem->node_ptr(10));
1090 0 : subelem[2]->set_node(9, elem->node_ptr(16));
1091 : }
1092 : }
1093 : else // Split on 1-3 diagonal
1094 : {
1095 0 : libmesh_assert (split_first_diagonal(elem, 1,3, 0,4));
1096 :
1097 : // Split on 0-5 diagonal
1098 0 : if (split_first_diagonal(elem, 0,5, 2,3))
1099 : {
1100 : // 1-3 and 0-5 split implies 1-5 split
1101 0 : libmesh_assert (split_first_diagonal(elem, 1,5, 2,4));
1102 :
1103 0 : subelem[0]->set_node(0, elem->node_ptr(1));
1104 0 : subelem[0]->set_node(1, elem->node_ptr(3));
1105 0 : subelem[0]->set_node(2, elem->node_ptr(4));
1106 0 : subelem[0]->set_node(3, elem->node_ptr(5));
1107 :
1108 0 : subelem[0]->set_node(4, elem->node_ptr(15));
1109 0 : subelem[0]->set_node(5, elem->node_ptr(12));
1110 0 : subelem[0]->set_node(6, elem->node_ptr(10));
1111 0 : subelem[0]->set_node(7, elem->node_ptr(16));
1112 0 : subelem[0]->set_node(8, elem->node_ptr(14));
1113 0 : subelem[0]->set_node(9, elem->node_ptr(13));
1114 :
1115 0 : subelem[1]->set_node(0, elem->node_ptr(1));
1116 0 : subelem[1]->set_node(1, elem->node_ptr(0));
1117 0 : subelem[1]->set_node(2, elem->node_ptr(3));
1118 0 : subelem[1]->set_node(3, elem->node_ptr(5));
1119 :
1120 0 : subelem[1]->set_node(4, elem->node_ptr(6));
1121 0 : subelem[1]->set_node(5, elem->node_ptr(9));
1122 0 : subelem[1]->set_node(6, elem->node_ptr(15));
1123 0 : subelem[1]->set_node(7, elem->node_ptr(16));
1124 0 : subelem[1]->set_node(8, elem->node_ptr(17));
1125 0 : subelem[1]->set_node(9, elem->node_ptr(14));
1126 :
1127 0 : subelem[2]->set_node(0, elem->node_ptr(0));
1128 0 : subelem[2]->set_node(1, elem->node_ptr(1));
1129 0 : subelem[2]->set_node(2, elem->node_ptr(2));
1130 0 : subelem[2]->set_node(3, elem->node_ptr(5));
1131 :
1132 0 : subelem[2]->set_node(4, elem->node_ptr(6));
1133 0 : subelem[2]->set_node(5, elem->node_ptr(7));
1134 0 : subelem[2]->set_node(6, elem->node_ptr(8));
1135 0 : subelem[2]->set_node(7, elem->node_ptr(17));
1136 0 : subelem[2]->set_node(8, elem->node_ptr(16));
1137 0 : subelem[2]->set_node(9, elem->node_ptr(11));
1138 : }
1139 : else // Split on 2-3 diagonal
1140 : {
1141 0 : libmesh_assert (split_first_diagonal(elem, 2,3, 0,5));
1142 :
1143 : // Split on 1-5 diagonal
1144 0 : if (split_first_diagonal(elem, 1,5, 2,4))
1145 : {
1146 0 : subelem[0]->set_node(0, elem->node_ptr(0));
1147 0 : subelem[0]->set_node(1, elem->node_ptr(1));
1148 0 : subelem[0]->set_node(2, elem->node_ptr(2));
1149 0 : subelem[0]->set_node(3, elem->node_ptr(3));
1150 :
1151 0 : subelem[0]->set_node(4, elem->node_ptr(6));
1152 0 : subelem[0]->set_node(5, elem->node_ptr(7));
1153 0 : subelem[0]->set_node(6, elem->node_ptr(8));
1154 0 : subelem[0]->set_node(7, elem->node_ptr(9));
1155 0 : subelem[0]->set_node(8, elem->node_ptr(15));
1156 0 : subelem[0]->set_node(9, elem->node_ptr(17));
1157 :
1158 0 : subelem[1]->set_node(0, elem->node_ptr(3));
1159 0 : subelem[1]->set_node(1, elem->node_ptr(1));
1160 0 : subelem[1]->set_node(2, elem->node_ptr(2));
1161 0 : subelem[1]->set_node(3, elem->node_ptr(5));
1162 :
1163 0 : subelem[1]->set_node(4, elem->node_ptr(15));
1164 0 : subelem[1]->set_node(5, elem->node_ptr(7));
1165 0 : subelem[1]->set_node(6, elem->node_ptr(17));
1166 0 : subelem[1]->set_node(7, elem->node_ptr(14));
1167 0 : subelem[1]->set_node(8, elem->node_ptr(16));
1168 0 : subelem[1]->set_node(9, elem->node_ptr(11));
1169 :
1170 0 : subelem[2]->set_node(0, elem->node_ptr(1));
1171 0 : subelem[2]->set_node(1, elem->node_ptr(3));
1172 0 : subelem[2]->set_node(2, elem->node_ptr(4));
1173 0 : subelem[2]->set_node(3, elem->node_ptr(5));
1174 :
1175 0 : subelem[2]->set_node(4, elem->node_ptr(15));
1176 0 : subelem[2]->set_node(5, elem->node_ptr(12));
1177 0 : subelem[2]->set_node(6, elem->node_ptr(10));
1178 0 : subelem[2]->set_node(7, elem->node_ptr(16));
1179 0 : subelem[2]->set_node(8, elem->node_ptr(14));
1180 0 : subelem[2]->set_node(9, elem->node_ptr(13));
1181 : }
1182 : else // Split on 2-4 diagonal
1183 : {
1184 0 : libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
1185 :
1186 0 : subelem[0]->set_node(0, elem->node_ptr(0));
1187 0 : subelem[0]->set_node(1, elem->node_ptr(1));
1188 0 : subelem[0]->set_node(2, elem->node_ptr(2));
1189 0 : subelem[0]->set_node(3, elem->node_ptr(3));
1190 :
1191 0 : subelem[0]->set_node(4, elem->node_ptr(6));
1192 0 : subelem[0]->set_node(5, elem->node_ptr(7));
1193 0 : subelem[0]->set_node(6, elem->node_ptr(8));
1194 0 : subelem[0]->set_node(7, elem->node_ptr(9));
1195 0 : subelem[0]->set_node(8, elem->node_ptr(15));
1196 0 : subelem[0]->set_node(9, elem->node_ptr(17));
1197 :
1198 0 : subelem[1]->set_node(0, elem->node_ptr(2));
1199 0 : subelem[1]->set_node(1, elem->node_ptr(3));
1200 0 : subelem[1]->set_node(2, elem->node_ptr(4));
1201 0 : subelem[1]->set_node(3, elem->node_ptr(5));
1202 :
1203 0 : subelem[1]->set_node(4, elem->node_ptr(17));
1204 0 : subelem[1]->set_node(5, elem->node_ptr(12));
1205 0 : subelem[1]->set_node(6, elem->node_ptr(16));
1206 0 : subelem[1]->set_node(7, elem->node_ptr(11));
1207 0 : subelem[1]->set_node(8, elem->node_ptr(14));
1208 0 : subelem[1]->set_node(9, elem->node_ptr(13));
1209 :
1210 0 : subelem[2]->set_node(0, elem->node_ptr(3));
1211 0 : subelem[2]->set_node(1, elem->node_ptr(1));
1212 0 : subelem[2]->set_node(2, elem->node_ptr(2));
1213 0 : subelem[2]->set_node(3, elem->node_ptr(4));
1214 :
1215 0 : subelem[2]->set_node(4, elem->node_ptr(15));
1216 0 : subelem[2]->set_node(5, elem->node_ptr(7));
1217 0 : subelem[2]->set_node(6, elem->node_ptr(17));
1218 0 : subelem[2]->set_node(7, elem->node_ptr(12));
1219 0 : subelem[2]->set_node(8, elem->node_ptr(10));
1220 0 : subelem[2]->set_node(9, elem->node_ptr(16));
1221 : }
1222 : }
1223 : }
1224 :
1225 12 : break;
1226 : }
1227 :
1228 : // No need to split elements that are already simplicial:
1229 20481 : case EDGE2:
1230 : case EDGE3:
1231 : case EDGE4:
1232 : case TRI3:
1233 : case TRI6:
1234 : case TRI7:
1235 : case TET4:
1236 : case TET10:
1237 : case TET14:
1238 : case INFEDGE2:
1239 : // No way to split infinite quad/prism elements, so
1240 : // hopefully no need to
1241 : case INFQUAD4:
1242 : case INFQUAD6:
1243 : case INFPRISM6:
1244 : case INFPRISM12:
1245 854 : continue;
1246 : // If we're left with an unimplemented hex we're probably
1247 : // out of luck. TODO: implement hexes
1248 0 : default:
1249 : {
1250 0 : libMesh::err << "Error, encountered unimplemented element "
1251 0 : << Utility::enum_to_string<ElemType>(etype)
1252 0 : << " in MeshTools::Modification::all_tri()..."
1253 0 : << std::endl;
1254 0 : libmesh_not_implemented();
1255 854 : }
1256 19627 : } // end switch (etype)
1257 :
1258 : // Be sure the correct data is set for all subelems.
1259 54990 : const unsigned int nei = elem->n_extra_integers();
1260 345382 : for (unsigned int i=0; i != max_subelems; ++i)
1261 301968 : if (subelem[i]) {
1262 286163 : subelem[i]->processor_id() = elem->processor_id();
1263 286163 : subelem[i]->subdomain_id() = elem->subdomain_id();
1264 :
1265 : // Copy any extra element data. Since the subelements
1266 : // haven't been added to the mesh yet any allocation has
1267 : // to be done manually.
1268 286163 : subelem[i]->add_extra_integers(nei);
1269 286163 : for (unsigned int ei=0; ei != nei; ++ei)
1270 0 : subelem[ei]->set_extra_integer(ei, elem->get_extra_integer(ei));
1271 :
1272 :
1273 : // Copy any mapping data.
1274 286163 : subelem[i]->set_mapping_type(elem->mapping_type());
1275 23056 : subelem[i]->set_mapping_data(elem->mapping_data());
1276 : }
1277 :
1278 : // On a mesh with boundary data, we need to move that data to
1279 : // the new elements.
1280 :
1281 : // On a mesh which is distributed, we need to move
1282 : // remote_elem links to the new elements.
1283 54990 : bool mesh_is_serial = mesh.is_serial();
1284 :
1285 54990 : if (mesh_has_boundary_data || !mesh_is_serial)
1286 : {
1287 : // Container to key boundary IDs handed back by the BoundaryInfo object.
1288 392 : std::vector<boundary_id_type> bc_ids;
1289 :
1290 87516 : for (auto sn : elem->side_index_range())
1291 : {
1292 73806 : mesh.get_boundary_info().boundary_ids(elem, sn, bc_ids);
1293 :
1294 73806 : if (bc_ids.empty() && elem->neighbor_ptr(sn) != remote_elem)
1295 59665 : continue;
1296 :
1297 : // Make a sorted list of node ids for elem->side(sn)
1298 14141 : elem->build_side_ptr(elem_side, sn);
1299 14425 : std::vector<dof_id_type> elem_side_nodes(elem_side->n_nodes());
1300 15713 : for (unsigned int esn=0,
1301 568 : n_esn = cast_int<unsigned int>(elem_side_nodes.size());
1302 68855 : esn != n_esn; ++esn)
1303 55642 : elem_side_nodes[esn] = elem_side->node_id(esn);
1304 14141 : std::sort(elem_side_nodes.begin(), elem_side_nodes.end());
1305 :
1306 79967 : for (unsigned int i=0; i != max_subelems; ++i)
1307 66650 : if (subelem[i])
1308 : {
1309 283200 : for (auto subside : subelem[i]->side_index_range())
1310 : {
1311 224658 : subelem[i]->build_side_ptr(subside_elem, subside);
1312 :
1313 : // Make a list of *vertex* node ids for this subside, see if they are all present
1314 : // in elem->side(sn). Note 1: we can't just compare elem->key(sn) to
1315 : // subelem[i]->key(subside) in the Prism cases, since the new side is
1316 : // a different type. Note 2: we only use vertex nodes since, in the future,
1317 : // a Hex20 or Prism15's QUAD8 face may be split into two Tri6 faces, and the
1318 : // original face will not contain the mid-edge node.
1319 226746 : std::vector<dof_id_type> subside_nodes(subside_elem->n_vertices());
1320 232458 : for (unsigned int ssn=0,
1321 4176 : n_ssn = cast_int<unsigned int>(subside_nodes.size());
1322 870102 : ssn != n_ssn; ++ssn)
1323 650388 : subside_nodes[ssn] = subside_elem->node_id(ssn);
1324 224658 : std::sort(subside_nodes.begin(), subside_nodes.end());
1325 :
1326 : // std::includes returns true if every element of the second sorted range is
1327 : // contained in the first sorted range.
1328 224658 : if (std::includes(elem_side_nodes.begin(), elem_side_nodes.end(),
1329 : subside_nodes.begin(), subside_nodes.end()))
1330 : {
1331 29493 : for (const auto & b_id : bc_ids)
1332 7102 : if (b_id != BoundaryInfo::invalid_id)
1333 : {
1334 7102 : new_bndry_ids.push_back(b_id);
1335 7102 : new_bndry_elements.push_back(subelem[i].get());
1336 7102 : new_bndry_sides.push_back(subside);
1337 : }
1338 :
1339 : // If the original element had a RemoteElem neighbor on side 'sn',
1340 : // then the subelem has one on side 'subside'.
1341 22707 : if (elem->neighbor_ptr(sn) == remote_elem)
1342 48 : subelem[i]->set_neighbor(subside, const_cast<RemoteElem*>(remote_elem));
1343 : }
1344 : }
1345 : } // end for loop over subelem
1346 : } // end for loop over sides
1347 :
1348 : // Remove the original element from the BoundaryInfo structure.
1349 13514 : mesh.get_boundary_info().remove(elem);
1350 :
1351 : } // end if (mesh_has_boundary_data)
1352 :
1353 : // Determine new IDs for the split elements which will be
1354 : // the same on all processors, therefore keeping the Mesh
1355 : // in sync. Note: we offset the new IDs by max_orig_id to
1356 : // avoid overwriting any of the original IDs.
1357 345382 : for (unsigned int i=0; i != max_subelems; ++i)
1358 301968 : if (subelem[i])
1359 : {
1360 : // Determine new IDs for the split elements which will be
1361 : // the same on all processors, therefore keeping the Mesh
1362 : // in sync. Note: we offset the new IDs by the max of the
1363 : // pre-existing ids to avoid conflicting with originals.
1364 286163 : subelem[i]->set_id( max_orig_id + 6*elem->id() + i );
1365 :
1366 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
1367 286163 : subelem[i]->set_unique_id(max_unique_id + max_subelems*elem->unique_id() + i);
1368 : #endif
1369 :
1370 : // Prepare to add the newly-created simplices
1371 11528 : new_elements.push_back(std::move(subelem[i]));
1372 : }
1373 :
1374 : // Delete the original element
1375 54990 : mesh.delete_elem(elem);
1376 2144 : } // End for loop over elements
1377 2144 : } // end scope
1378 :
1379 :
1380 : // Now, iterate over the new elements vector, and add them each to
1381 : // the Mesh.
1382 288435 : for (auto & elem : new_elements)
1383 309219 : mesh.add_elem(std::move(elem));
1384 :
1385 2272 : if (mesh_has_boundary_data)
1386 : {
1387 : // If the old mesh had boundary data, the new mesh better have
1388 : // some. However, we can't assert that the size of
1389 : // new_bndry_elements vector is > 0, since we may not have split
1390 : // any elements actually on the boundary. We also can't assert
1391 : // that the original number of boundary sides is equal to the
1392 : // sum of the boundary sides currently in the mesh and the
1393 : // newly-added boundary sides, since in 3D, we may have split a
1394 : // boundary QUAD into two boundary TRIs. Therefore, we won't be
1395 : // too picky about the actual number of BCs, and just assert that
1396 : // there are some, somewhere.
1397 : #ifndef NDEBUG
1398 36 : bool nbe_nonempty = new_bndry_elements.size();
1399 36 : mesh.comm().max(nbe_nonempty);
1400 36 : libmesh_assert(nbe_nonempty ||
1401 : mesh.get_boundary_info().n_boundary_conds()>0);
1402 : #endif
1403 :
1404 : // We should also be sure that the lengths of the new boundary data vectors
1405 : // are all the same.
1406 36 : libmesh_assert_equal_to (new_bndry_elements.size(), new_bndry_sides.size());
1407 36 : libmesh_assert_equal_to (new_bndry_sides.size(), new_bndry_ids.size());
1408 :
1409 : // Add the new boundary info to the mesh
1410 8380 : for (auto s : index_range(new_bndry_elements))
1411 7686 : mesh.get_boundary_info().add_side(new_bndry_elements[s],
1412 584 : new_bndry_sides[s],
1413 584 : new_bndry_ids[s]);
1414 : }
1415 :
1416 : // In a DistributedMesh any newly added ghost node ids may be
1417 : // inconsistent, and unique_ids of newly added ghost nodes remain
1418 : // unset.
1419 : // make_nodes_parallel_consistent() will fix all this.
1420 2272 : if (!mesh.is_serial())
1421 : {
1422 903 : mesh.comm().max(added_new_ghost_point);
1423 :
1424 903 : if (added_new_ghost_point)
1425 0 : MeshCommunication().make_nodes_parallel_consistent (mesh);
1426 : }
1427 :
1428 :
1429 :
1430 : // Prepare the newly created mesh for use.
1431 2272 : mesh.prepare_for_use();
1432 :
1433 : // Let the new_elements and new_bndry_elements vectors go out of scope.
1434 2540 : }
1435 :
1436 :
1437 0 : void MeshTools::Modification::smooth (MeshBase & mesh,
1438 : const unsigned int n_iterations,
1439 : const Real power)
1440 : {
1441 : /**
1442 : * This implementation assumes every element "side" has only 2 nodes.
1443 : */
1444 0 : libmesh_assert_equal_to (mesh.mesh_dimension(), 2);
1445 :
1446 : /*
1447 : * Create a quickly-searchable list of boundary nodes.
1448 : */
1449 : std::unordered_set<dof_id_type> boundary_node_ids =
1450 0 : MeshTools::find_boundary_nodes (mesh);
1451 :
1452 : // For avoiding extraneous element side allocation
1453 0 : ElemSideBuilder side_builder;
1454 :
1455 0 : for (unsigned int iter=0; iter<n_iterations; iter++)
1456 : {
1457 : /*
1458 : * loop over the mesh refinement level
1459 : */
1460 0 : unsigned int n_levels = MeshTools::n_levels(mesh);
1461 0 : for (unsigned int refinement_level=0; refinement_level != n_levels;
1462 : refinement_level++)
1463 : {
1464 : // initialize the storage (have to do it on every level to get empty vectors
1465 0 : std::vector<Point> new_positions;
1466 0 : std::vector<Real> weight;
1467 0 : new_positions.resize(mesh.n_nodes());
1468 0 : weight.resize(mesh.n_nodes());
1469 :
1470 : {
1471 : // Loop over the elements to calculate new node positions
1472 0 : for (const auto & elem : as_range(mesh.level_elements_begin(refinement_level),
1473 0 : mesh.level_elements_end(refinement_level)))
1474 : {
1475 : /*
1476 : * We relax all nodes on level 0 first
1477 : * If the element is refined (level > 0), we interpolate the
1478 : * parents nodes with help of the embedding matrix
1479 : */
1480 0 : if (refinement_level == 0)
1481 : {
1482 0 : for (auto s : elem->side_index_range())
1483 : {
1484 : /*
1485 : * Only operate on sides which are on the
1486 : * boundary or for which the current element's
1487 : * id is greater than its neighbor's.
1488 : * Sides get only built once.
1489 : */
1490 0 : if ((elem->neighbor_ptr(s) != nullptr) &&
1491 0 : (elem->id() > elem->neighbor_ptr(s)->id()))
1492 : {
1493 0 : const Elem & side = side_builder(*elem, s);
1494 0 : const Node & node0 = side.node_ref(0);
1495 0 : const Node & node1 = side.node_ref(1);
1496 :
1497 0 : Real node_weight = 1.;
1498 : // calculate the weight of the nodes
1499 0 : if (power > 0)
1500 : {
1501 0 : Point diff = node0-node1;
1502 0 : node_weight = std::pow(diff.norm(), power);
1503 : }
1504 :
1505 0 : const dof_id_type id0 = node0.id(), id1 = node1.id();
1506 0 : new_positions[id0].add_scaled( node1, node_weight );
1507 0 : new_positions[id1].add_scaled( node0, node_weight );
1508 0 : weight[id0] += node_weight;
1509 0 : weight[id1] += node_weight;
1510 : }
1511 : } // element neighbor loop
1512 : }
1513 : #ifdef LIBMESH_ENABLE_AMR
1514 : else // refinement_level > 0
1515 : {
1516 : /*
1517 : * Find the positions of the hanging nodes of refined elements.
1518 : * We do this by calculating their position based on the parent
1519 : * (one level less refined) element, and the embedding matrix
1520 : */
1521 :
1522 0 : const Elem * parent = elem->parent();
1523 :
1524 : /*
1525 : * find out which child I am
1526 : */
1527 0 : unsigned int c = parent->which_child_am_i(elem);
1528 : /*
1529 : *loop over the childs (that is, the current elements) nodes
1530 : */
1531 0 : for (auto nc : elem->node_index_range())
1532 : {
1533 : /*
1534 : * the new position of the node
1535 : */
1536 0 : Point point;
1537 0 : for (auto n : parent->node_index_range())
1538 : {
1539 : /*
1540 : * The value from the embedding matrix
1541 : */
1542 0 : const Real em_val = parent->embedding_matrix(c,nc,n);
1543 :
1544 0 : if (em_val != 0.)
1545 0 : point.add_scaled (parent->point(n), em_val);
1546 : }
1547 :
1548 0 : const dof_id_type id = elem->node_ptr(nc)->id();
1549 0 : new_positions[id] = point;
1550 0 : weight[id] = 1.;
1551 : }
1552 : } // if element refinement_level
1553 : #endif // #ifdef LIBMESH_ENABLE_AMR
1554 :
1555 0 : } // element loop
1556 :
1557 : /*
1558 : * finally reposition the vertex nodes
1559 : */
1560 0 : for (auto nid : make_range(mesh.n_nodes()))
1561 0 : if (!boundary_node_ids.count(nid) && weight[nid] > 0.)
1562 0 : mesh.node_ref(nid) = new_positions[nid]/weight[nid];
1563 : }
1564 :
1565 : // Now handle the additional second_order nodes by calculating
1566 : // their position based on the vertex positions
1567 : // we do a second loop over the level elements
1568 0 : for (auto & elem : as_range(mesh.level_elements_begin(refinement_level),
1569 0 : mesh.level_elements_end(refinement_level)))
1570 : {
1571 0 : const unsigned int son_begin = elem->n_vertices();
1572 0 : const unsigned int son_end = elem->n_nodes();
1573 0 : for (unsigned int n=son_begin; n<son_end; n++)
1574 : {
1575 : const unsigned int n_adjacent_vertices =
1576 0 : elem->n_second_order_adjacent_vertices(n);
1577 :
1578 0 : Point point;
1579 0 : for (unsigned int v=0; v<n_adjacent_vertices; v++)
1580 0 : point.add(elem->point( elem->second_order_adjacent_vertex(n,v) ));
1581 :
1582 0 : const dof_id_type id = elem->node_ptr(n)->id();
1583 0 : mesh.node_ref(id) = point/n_adjacent_vertices;
1584 : }
1585 0 : }
1586 : } // refinement_level loop
1587 : } // end iteration
1588 0 : }
1589 :
1590 :
1591 :
1592 : #ifdef LIBMESH_ENABLE_AMR
1593 1145 : void MeshTools::Modification::flatten(MeshBase & mesh)
1594 : {
1595 36 : libmesh_assert(mesh.is_prepared() || mesh.is_replicated());
1596 :
1597 : // Algorithm:
1598 : // .) For each active element in the mesh: construct a
1599 : // copy which is the same in every way *except* it is
1600 : // a level 0 element. Store the pointers to these in
1601 : // a separate vector. Save any boundary information as well.
1602 : // Delete the active element from the mesh.
1603 : // .) Loop over all (remaining) elements in the mesh, delete them.
1604 : // .) Add the level-0 copies back to the mesh
1605 :
1606 : // Temporary storage for new element pointers
1607 108 : std::vector<std::unique_ptr<Elem>> new_elements;
1608 :
1609 : // BoundaryInfo Storage for element ids, sides, and BC ids
1610 72 : std::vector<Elem *> saved_boundary_elements;
1611 72 : std::vector<boundary_id_type> saved_bc_ids;
1612 72 : std::vector<unsigned short int> saved_bc_sides;
1613 :
1614 : // Container to catch boundary ids passed back by BoundaryInfo
1615 72 : std::vector<boundary_id_type> bc_ids;
1616 :
1617 : // Reserve a reasonable amt. of space for each
1618 1145 : new_elements.reserve(mesh.n_active_elem());
1619 1145 : saved_boundary_elements.reserve(mesh.get_boundary_info().n_boundary_conds());
1620 1145 : saved_bc_ids.reserve(mesh.get_boundary_info().n_boundary_conds());
1621 1145 : saved_bc_sides.reserve(mesh.get_boundary_info().n_boundary_conds());
1622 :
1623 329648 : for (auto & elem : mesh.active_element_ptr_range())
1624 : {
1625 : // Make a new element of the same type
1626 177481 : auto copy = Elem::build(elem->type());
1627 :
1628 : // Set node pointers (they still point to nodes in the original mesh)
1629 1374787 : for (auto n : elem->node_index_range())
1630 1247032 : copy->set_node(n, elem->node_ptr(n));
1631 :
1632 : // Copy over ids
1633 170589 : copy->processor_id() = elem->processor_id();
1634 170589 : copy->subdomain_id() = elem->subdomain_id();
1635 :
1636 : // Retain the original element's ID(s) as well, otherwise
1637 : // the Mesh may try to create them for you...
1638 13784 : copy->set_id( elem->id() );
1639 : #ifdef LIBMESH_ENABLE_UNIQUE_ID
1640 13784 : copy->set_unique_id(elem->unique_id());
1641 : #endif
1642 :
1643 : // This element could have boundary info or DistributedMesh
1644 : // remote_elem links as well. We need to save the (elem,
1645 : // side, bc_id) triples and those links
1646 1098669 : for (auto s : elem->side_index_range())
1647 : {
1648 966376 : if (elem->neighbor_ptr(s) == remote_elem)
1649 744 : copy->set_neighbor(s, const_cast<RemoteElem *>(remote_elem));
1650 :
1651 928080 : mesh.get_boundary_info().boundary_ids(elem, s, bc_ids);
1652 928414 : for (const auto & bc_id : bc_ids)
1653 334 : if (bc_id != BoundaryInfo::invalid_id)
1654 : {
1655 334 : saved_boundary_elements.push_back(copy.get());
1656 334 : saved_bc_ids.push_back(bc_id);
1657 334 : saved_bc_sides.push_back(s);
1658 : }
1659 : }
1660 :
1661 : // Copy any extra element data. Since the copy hasn't been
1662 : // added to the mesh yet any allocation has to be done manually.
1663 170589 : const unsigned int nei = elem->n_extra_integers();
1664 170589 : copy->add_extra_integers(nei);
1665 170589 : for (unsigned int i=0; i != nei; ++i)
1666 0 : copy->set_extra_integer(i, elem->get_extra_integer(i));
1667 :
1668 : // Copy any mapping data.
1669 170589 : copy->set_mapping_type(elem->mapping_type());
1670 13784 : copy->set_mapping_data(elem->mapping_data());
1671 :
1672 : // We're done with this element
1673 170589 : mesh.delete_elem(elem);
1674 :
1675 : // But save the copy
1676 6892 : new_elements.push_back(std::move(copy));
1677 157878 : }
1678 :
1679 : // Make sure we saved the same number of boundary conditions
1680 : // in each vector.
1681 36 : libmesh_assert_equal_to (saved_boundary_elements.size(), saved_bc_ids.size());
1682 36 : libmesh_assert_equal_to (saved_bc_ids.size(), saved_bc_sides.size());
1683 :
1684 : // Loop again, delete any remaining elements
1685 71004 : for (auto & elem : mesh.element_ptr_range())
1686 36722 : mesh.delete_elem(elem);
1687 :
1688 : // Add the copied (now level-0) elements back to the mesh
1689 171734 : for (auto & new_elem : new_elements)
1690 : {
1691 : // Save the original ID, because the act of adding the Elem can
1692 : // change new_elem's id!
1693 6892 : dof_id_type orig_id = new_elem->id();
1694 :
1695 184373 : Elem * added_elem = mesh.add_elem(std::move(new_elem));
1696 :
1697 : // If the Elem, as it was re-added to the mesh, now has a
1698 : // different ID (this is unlikely, so it's just an assert)
1699 : // the boundary information will no longer be correct.
1700 6892 : libmesh_assert_equal_to (orig_id, added_elem->id());
1701 :
1702 : // Avoid compiler warnings in opt mode.
1703 6892 : libmesh_ignore(added_elem, orig_id);
1704 : }
1705 :
1706 : // Finally, also add back the saved boundary information
1707 1479 : for (auto e : index_range(saved_boundary_elements))
1708 358 : mesh.get_boundary_info().add_side(saved_boundary_elements[e],
1709 24 : saved_bc_sides[e],
1710 24 : saved_bc_ids[e]);
1711 :
1712 : // Trim unused and renumber nodes and elements
1713 1145 : mesh.prepare_for_use();
1714 1145 : }
1715 : #endif // #ifdef LIBMESH_ENABLE_AMR
1716 :
1717 :
1718 :
1719 3408 : void MeshTools::Modification::change_boundary_id (MeshBase & mesh,
1720 : const boundary_id_type old_id,
1721 : const boundary_id_type new_id)
1722 : {
1723 : // This is just a shim around the member implementation, now
1724 3408 : mesh.get_boundary_info().renumber_id(old_id, new_id);
1725 3408 : }
1726 :
1727 :
1728 :
1729 0 : void MeshTools::Modification::change_subdomain_id (MeshBase & mesh,
1730 : const subdomain_id_type old_id,
1731 : const subdomain_id_type new_id)
1732 : {
1733 0 : if (old_id == new_id)
1734 : {
1735 : // If the IDs are the same, this is a no-op.
1736 0 : return;
1737 : }
1738 :
1739 0 : for (auto & elem : mesh.element_ptr_range())
1740 : {
1741 0 : if (elem->subdomain_id() == old_id)
1742 0 : elem->subdomain_id() = new_id;
1743 0 : }
1744 : }
1745 :
1746 :
1747 : } // namespace libMesh
|
