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