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