mesh_modification.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 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 
27 // Local includes
28 #include "libmesh/boundary_info.h"
29 #include "libmesh/function_base.h"
30 #include "libmesh/cell_tet4.h"
31 #include "libmesh/cell_tet10.h"
32 #include "libmesh/face_tri3.h"
33 #include "libmesh/face_tri6.h"
37 #include "libmesh/mesh_tools.h"
38 #include "libmesh/parallel.h"
39 #include "libmesh/remote_elem.h"
40 #include "libmesh/string_to_enum.h"
42 
43 namespace
44 {
45 bool split_first_diagonal(const libMesh::Elem * elem,
46  unsigned int diag_1_node_1,
47  unsigned int diag_1_node_2,
48  unsigned int diag_2_node_1,
49  unsigned int diag_2_node_2)
50 {
51  return ((elem->node_id(diag_1_node_1) > elem->node_id(diag_2_node_1) &&
52  elem->node_id(diag_1_node_1) > elem->node_id(diag_2_node_2)) ||
53  (elem->node_id(diag_1_node_2) > elem->node_id(diag_2_node_1) &&
54  elem->node_id(diag_1_node_2) > elem->node_id(diag_2_node_2)));
55 }
56 
57 }
58 
59 namespace libMesh
60 {
61 
62 
63 // ------------------------------------------------------------
64 // MeshTools::Modification functions for mesh modification
66  const Real factor,
67  const bool perturb_boundary)
68 {
69  libmesh_assert (mesh.n_nodes());
70  libmesh_assert (mesh.n_elem());
71  libmesh_assert ((factor >= 0.) && (factor <= 1.));
72 
73  LOG_SCOPE("distort()", "MeshTools::Modification");
74 
75  // If we are not perturbing boundary nodes, make a
76  // quickly-searchable list of node ids we can check against.
77  std::unordered_set<dof_id_type> boundary_node_ids;
78  if (!perturb_boundary)
79  boundary_node_ids = MeshTools::find_boundary_nodes (mesh);
80 
81  // Now calculate the minimum distance to
82  // neighboring nodes for each node.
83  // hmin holds these distances.
84  std::vector<float> hmin (mesh.max_node_id(),
86 
87  for (const auto & elem : mesh.active_element_ptr_range())
88  for (auto & n : elem->node_ref_range())
89  hmin[n.id()] = std::min(hmin[n.id()],
90  static_cast<float>(elem->hmin()));
91 
92  // Now actually move the nodes
93  {
94  const unsigned int seed = 123456;
95 
96  // seed the random number generator.
97  // We'll loop from 1 to n_nodes on every processor, even those
98  // that don't have a particular node, so that the pseudorandom
99  // numbers will be the same everywhere.
100  std::srand(seed);
101 
102  // If the node is on the boundary or
103  // the node is not used by any element (hmin[n]<1.e20)
104  // then we should not move it.
105  // [Note: Testing for (in)equality might be wrong
106  // (different types, namely float and double)]
107  for (unsigned int n=0; n<mesh.max_node_id(); n++)
108  if ((perturb_boundary || !boundary_node_ids.count(n)) && hmin[n] < 1.e20)
109  {
110  // the direction, random but unit normalized
111  Point dir (static_cast<Real>(std::rand())/static_cast<Real>(RAND_MAX),
112  (mesh.mesh_dimension() > 1) ? static_cast<Real>(std::rand())/static_cast<Real>(RAND_MAX) : 0.,
113  ((mesh.mesh_dimension() == 3) ? static_cast<Real>(std::rand())/static_cast<Real>(RAND_MAX) : 0.));
114 
115  dir(0) = (dir(0)-.5)*2.;
116  if (mesh.mesh_dimension() > 1)
117  dir(1) = (dir(1)-.5)*2.;
118  if (mesh.mesh_dimension() == 3)
119  dir(2) = (dir(2)-.5)*2.;
120 
121  dir = dir.unit();
122 
123  Node * node = mesh.query_node_ptr(n);
124  if (!node)
125  continue;
126 
127  (*node)(0) += dir(0)*factor*hmin[n];
128  if (mesh.mesh_dimension() > 1)
129  (*node)(1) += dir(1)*factor*hmin[n];
130  if (mesh.mesh_dimension() == 3)
131  (*node)(2) += dir(2)*factor*hmin[n];
132  }
133  }
134 }
135 
136 
137 
139  const FunctionBase<Real> & mapfunc)
140 {
141  libmesh_assert (mesh.n_nodes());
142  libmesh_assert (mesh.n_elem());
143 
144  LOG_SCOPE("redistribute()", "MeshTools::Modification");
145 
146  DenseVector<Real> output_vec(LIBMESH_DIM);
147 
148  // FIXME - we should thread this later.
149  std::unique_ptr<FunctionBase<Real>> myfunc = mapfunc.clone();
150 
151  for (auto & node : mesh.node_ptr_range())
152  {
153  (*myfunc)(*node, output_vec);
154 
155  (*node)(0) = output_vec(0);
156 #if LIBMESH_DIM > 1
157  (*node)(1) = output_vec(1);
158 #endif
159 #if LIBMESH_DIM > 2
160  (*node)(2) = output_vec(2);
161 #endif
162  }
163 }
164 
165 
166 
168  const Real xt,
169  const Real yt,
170  const Real zt)
171 {
172  const Point p(xt, yt, zt);
173 
174  for (auto & node : mesh.node_ptr_range())
175  *node += p;
176 }
177 
178 
179 // void MeshTools::Modification::rotate2D (MeshBase & mesh,
180 // const Real alpha)
181 // {
182 // libmesh_assert_not_equal_to (mesh.mesh_dimension(), 1);
183 
184 // const Real pi = std::acos(-1);
185 // const Real a = alpha/180.*pi;
186 // for (unsigned int n=0; n<mesh.n_nodes(); n++)
187 // {
188 // const Point p = mesh.node_ref(n);
189 // const Real x = p(0);
190 // const Real y = p(1);
191 // const Real z = p(2);
192 // mesh.node_ref(n) = Point(std::cos(a)*x - std::sin(a)*y,
193 // std::sin(a)*x + std::cos(a)*y,
194 // z);
195 // }
196 
197 // }
198 
199 
200 
202  const Real phi,
203  const Real theta,
204  const Real psi)
205 {
206 #if LIBMESH_DIM == 3
207  const Real p = -phi/180.*libMesh::pi;
208  const Real t = -theta/180.*libMesh::pi;
209  const Real s = -psi/180.*libMesh::pi;
210  const Real sp = std::sin(p), cp = std::cos(p);
211  const Real st = std::sin(t), ct = std::cos(t);
212  const Real ss = std::sin(s), cs = std::cos(s);
213 
214  // We follow the convention described at http://mathworld.wolfram.com/EulerAngles.html
215  // (equations 6-14 give the entries of the composite transformation matrix).
216  // The rotations are performed sequentially about the z, x, and z axes, in that order.
217  // A positive angle yields a counter-clockwise rotation about the axis in question.
218  for (auto & node : mesh.node_ptr_range())
219  {
220  const Point pt = *node;
221  const Real x = pt(0);
222  const Real y = pt(1);
223  const Real z = pt(2);
224  *node = Point(( cp*cs-sp*ct*ss)*x + ( sp*cs+cp*ct*ss)*y + (st*ss)*z,
225  (-cp*ss-sp*ct*cs)*x + (-sp*ss+cp*ct*cs)*y + (st*cs)*z,
226  ( sp*st)*x + (-cp*st)*y + (ct)*z );
227  }
228 #else
229  libmesh_error_msg("MeshTools::Modification::rotate() requires libMesh to be compiled with LIBMESH_DIM==3");
230 #endif
231 }
232 
233 
235  const Real xs,
236  const Real ys,
237  const Real zs)
238 {
239  const Real x_scale = xs;
240  Real y_scale = ys;
241  Real z_scale = zs;
242 
243  if (ys == 0.)
244  {
245  libmesh_assert_equal_to (zs, 0.);
246 
247  y_scale = z_scale = x_scale;
248  }
249 
250  // Scale the x coordinate in all dimensions
251  for (auto & node : mesh.node_ptr_range())
252  (*node)(0) *= x_scale;
253 
254  // Only scale the y coordinate in 2 and 3D
255  if (LIBMESH_DIM < 2)
256  return;
257 
258  for (auto & node : mesh.node_ptr_range())
259  (*node)(1) *= y_scale;
260 
261  // Only scale the z coordinate in 3D
262  if (LIBMESH_DIM < 3)
263  return;
264 
265  for (auto & node : mesh.node_ptr_range())
266  (*node)(2) *= z_scale;
267 }
268 
269 
270 
272 {
273  // The number of elements in the original mesh before any additions
274  // or deletions.
275  const dof_id_type n_orig_elem = mesh.n_elem();
276  const dof_id_type max_orig_id = mesh.max_elem_id();
277 
278  // We store pointers to the newly created elements in a vector
279  // until they are ready to be added to the mesh. This is because
280  // adding new elements on the fly can cause reallocation and invalidation
281  // of existing mesh element_iterators.
282  std::vector<Elem *> new_elements;
283 
284  unsigned int max_subelems = 1; // in 1D nothing needs to change
285  if (mesh.mesh_dimension() == 2) // in 2D quads can split into 2 tris
286  max_subelems = 2;
287  if (mesh.mesh_dimension() == 3) // in 3D hexes can split into 6 tets
288  max_subelems = 6;
289 
290  new_elements.reserve (max_subelems*n_orig_elem);
291 
292  // If the original mesh has *side* boundary data, we carry that over
293  // to the new mesh with triangular elements. We currently only
294  // support bringing over side-based BCs to the all-tri mesh, but
295  // that could probably be extended to node and edge-based BCs as
296  // well.
297  const bool mesh_has_boundary_data = (mesh.get_boundary_info().n_boundary_conds() > 0);
298 
299  // Temporary vectors to store the new boundary element pointers, side numbers, and boundary ids
300  std::vector<Elem *> new_bndry_elements;
301  std::vector<unsigned short int> new_bndry_sides;
302  std::vector<boundary_id_type> new_bndry_ids;
303 
304  // We may need to add new points if we run into a 1.5th order
305  // element; if we do that on a DistributedMesh in a ghost element then
306  // we will need to fix their ids / unique_ids
307  bool added_new_ghost_point = false;
308 
309  // Iterate over the elements, splitting:
310  // QUADs into pairs of conforming triangles
311  // PYRAMIDs into pairs of conforming tets,
312  // PRISMs into triplets of conforming tets, and
313  // HEXs into quintets or sextets of conforming tets.
314  // We split on the shortest diagonal to give us better
315  // triangle quality in 2D, and we split based on node ids
316  // to guarantee consistency in 3D.
317 
318  // FIXME: This algorithm does not work on refined grids!
319  {
320 #ifdef LIBMESH_ENABLE_UNIQUE_ID
321  unique_id_type max_unique_id = mesh.parallel_max_unique_id();
322 #endif
323 
324  for (auto & elem : mesh.element_ptr_range())
325  {
326  const ElemType etype = elem->type();
327 
328  // all_tri currently only works on coarse meshes
329  libmesh_assert (!elem->parent());
330 
331  // The new elements we will split the quad into.
332  // In 3D we may need as many as 6 tets per hex
333  Elem * subelem[6];
334 
335  for (unsigned int i = 0; i != max_subelems; ++i)
336  subelem[i] = nullptr;
337 
338  switch (etype)
339  {
340  case QUAD4:
341  {
342  subelem[0] = new Tri3;
343  subelem[1] = new Tri3;
344 
345  // Check for possible edge swap
346  if ((elem->point(0) - elem->point(2)).norm() <
347  (elem->point(1) - elem->point(3)).norm())
348  {
349  subelem[0]->set_node(0) = elem->node_ptr(0);
350  subelem[0]->set_node(1) = elem->node_ptr(1);
351  subelem[0]->set_node(2) = elem->node_ptr(2);
352 
353  subelem[1]->set_node(0) = elem->node_ptr(0);
354  subelem[1]->set_node(1) = elem->node_ptr(2);
355  subelem[1]->set_node(2) = elem->node_ptr(3);
356  }
357 
358  else
359  {
360  subelem[0]->set_node(0) = elem->node_ptr(0);
361  subelem[0]->set_node(1) = elem->node_ptr(1);
362  subelem[0]->set_node(2) = elem->node_ptr(3);
363 
364  subelem[1]->set_node(0) = elem->node_ptr(1);
365  subelem[1]->set_node(1) = elem->node_ptr(2);
366  subelem[1]->set_node(2) = elem->node_ptr(3);
367  }
368 
369 
370  break;
371  }
372 
373  case QUAD8:
374  {
375  if (elem->processor_id() != mesh.processor_id())
376  added_new_ghost_point = true;
377 
378  subelem[0] = new Tri6;
379  subelem[1] = new Tri6;
380 
381  // Add a new node at the center (vertex average) of the element.
382  Node * new_node = mesh.add_point((mesh.point(elem->node_id(0)) +
383  mesh.point(elem->node_id(1)) +
384  mesh.point(elem->node_id(2)) +
385  mesh.point(elem->node_id(3)))/4,
387  elem->processor_id());
388 
389  // Check for possible edge swap
390  if ((elem->point(0) - elem->point(2)).norm() <
391  (elem->point(1) - elem->point(3)).norm())
392  {
393  subelem[0]->set_node(0) = elem->node_ptr(0);
394  subelem[0]->set_node(1) = elem->node_ptr(1);
395  subelem[0]->set_node(2) = elem->node_ptr(2);
396  subelem[0]->set_node(3) = elem->node_ptr(4);
397  subelem[0]->set_node(4) = elem->node_ptr(5);
398  subelem[0]->set_node(5) = new_node;
399 
400  subelem[1]->set_node(0) = elem->node_ptr(0);
401  subelem[1]->set_node(1) = elem->node_ptr(2);
402  subelem[1]->set_node(2) = elem->node_ptr(3);
403  subelem[1]->set_node(3) = new_node;
404  subelem[1]->set_node(4) = elem->node_ptr(6);
405  subelem[1]->set_node(5) = elem->node_ptr(7);
406 
407  }
408 
409  else
410  {
411  subelem[0]->set_node(0) = elem->node_ptr(3);
412  subelem[0]->set_node(1) = elem->node_ptr(0);
413  subelem[0]->set_node(2) = elem->node_ptr(1);
414  subelem[0]->set_node(3) = elem->node_ptr(7);
415  subelem[0]->set_node(4) = elem->node_ptr(4);
416  subelem[0]->set_node(5) = new_node;
417 
418  subelem[1]->set_node(0) = elem->node_ptr(1);
419  subelem[1]->set_node(1) = elem->node_ptr(2);
420  subelem[1]->set_node(2) = elem->node_ptr(3);
421  subelem[1]->set_node(3) = elem->node_ptr(5);
422  subelem[1]->set_node(4) = elem->node_ptr(6);
423  subelem[1]->set_node(5) = new_node;
424  }
425 
426  break;
427  }
428 
429  case QUAD9:
430  {
431  subelem[0] = new Tri6;
432  subelem[1] = new Tri6;
433 
434  // Check for possible edge swap
435  if ((elem->point(0) - elem->point(2)).norm() <
436  (elem->point(1) - elem->point(3)).norm())
437  {
438  subelem[0]->set_node(0) = elem->node_ptr(0);
439  subelem[0]->set_node(1) = elem->node_ptr(1);
440  subelem[0]->set_node(2) = elem->node_ptr(2);
441  subelem[0]->set_node(3) = elem->node_ptr(4);
442  subelem[0]->set_node(4) = elem->node_ptr(5);
443  subelem[0]->set_node(5) = elem->node_ptr(8);
444 
445  subelem[1]->set_node(0) = elem->node_ptr(0);
446  subelem[1]->set_node(1) = elem->node_ptr(2);
447  subelem[1]->set_node(2) = elem->node_ptr(3);
448  subelem[1]->set_node(3) = elem->node_ptr(8);
449  subelem[1]->set_node(4) = elem->node_ptr(6);
450  subelem[1]->set_node(5) = elem->node_ptr(7);
451  }
452 
453  else
454  {
455  subelem[0]->set_node(0) = elem->node_ptr(0);
456  subelem[0]->set_node(1) = elem->node_ptr(1);
457  subelem[0]->set_node(2) = elem->node_ptr(3);
458  subelem[0]->set_node(3) = elem->node_ptr(4);
459  subelem[0]->set_node(4) = elem->node_ptr(8);
460  subelem[0]->set_node(5) = elem->node_ptr(7);
461 
462  subelem[1]->set_node(0) = elem->node_ptr(1);
463  subelem[1]->set_node(1) = elem->node_ptr(2);
464  subelem[1]->set_node(2) = elem->node_ptr(3);
465  subelem[1]->set_node(3) = elem->node_ptr(5);
466  subelem[1]->set_node(4) = elem->node_ptr(6);
467  subelem[1]->set_node(5) = elem->node_ptr(8);
468  }
469 
470  break;
471  }
472 
473  case PRISM6:
474  {
475  // Prisms all split into three tetrahedra
476  subelem[0] = new Tet4;
477  subelem[1] = new Tet4;
478  subelem[2] = new Tet4;
479 
480  // Triangular faces are not split.
481 
482  // On quad faces, we choose the node with the highest
483  // global id, and we split on the diagonal which
484  // includes that node. This ensures that (even in
485  // parallel, even on distributed meshes) the same
486  // diagonal split will be chosen for elements on either
487  // side of the same quad face. It also ensures that we
488  // always have a mix of "clockwise" and
489  // "counterclockwise" split faces (two of one and one
490  // of the other on each prism; this is useful since the
491  // alternative all-clockwise or all-counterclockwise
492  // face splittings can't be turned into tets without
493  // adding more nodes
494 
495  // Split on 0-4 diagonal
496  if (split_first_diagonal(elem, 0,4, 1,3))
497  {
498  // Split on 0-5 diagonal
499  if (split_first_diagonal(elem, 0,5, 2,3))
500  {
501  // Split on 1-5 diagonal
502  if (split_first_diagonal(elem, 1,5, 2,4))
503  {
504  subelem[0]->set_node(0) = elem->node_ptr(0);
505  subelem[0]->set_node(1) = elem->node_ptr(4);
506  subelem[0]->set_node(2) = elem->node_ptr(5);
507  subelem[0]->set_node(3) = elem->node_ptr(3);
508 
509  subelem[1]->set_node(0) = elem->node_ptr(0);
510  subelem[1]->set_node(1) = elem->node_ptr(4);
511  subelem[1]->set_node(2) = elem->node_ptr(1);
512  subelem[1]->set_node(3) = elem->node_ptr(5);
513 
514  subelem[2]->set_node(0) = elem->node_ptr(0);
515  subelem[2]->set_node(1) = elem->node_ptr(1);
516  subelem[2]->set_node(2) = elem->node_ptr(2);
517  subelem[2]->set_node(3) = elem->node_ptr(5);
518  }
519  else // Split on 2-4 diagonal
520  {
521  libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
522 
523  subelem[0]->set_node(0) = elem->node_ptr(0);
524  subelem[0]->set_node(1) = elem->node_ptr(4);
525  subelem[0]->set_node(2) = elem->node_ptr(5);
526  subelem[0]->set_node(3) = elem->node_ptr(3);
527 
528  subelem[1]->set_node(0) = elem->node_ptr(0);
529  subelem[1]->set_node(1) = elem->node_ptr(4);
530  subelem[1]->set_node(2) = elem->node_ptr(2);
531  subelem[1]->set_node(3) = elem->node_ptr(5);
532 
533  subelem[2]->set_node(0) = elem->node_ptr(0);
534  subelem[2]->set_node(1) = elem->node_ptr(1);
535  subelem[2]->set_node(2) = elem->node_ptr(2);
536  subelem[2]->set_node(3) = elem->node_ptr(4);
537  }
538  }
539  else // Split on 2-3 diagonal
540  {
541  libmesh_assert (split_first_diagonal(elem, 2,3, 0,5));
542 
543  // 0-4 and 2-3 split implies 2-4 split
544  libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
545 
546  subelem[0]->set_node(0) = elem->node_ptr(0);
547  subelem[0]->set_node(1) = elem->node_ptr(4);
548  subelem[0]->set_node(2) = elem->node_ptr(2);
549  subelem[0]->set_node(3) = elem->node_ptr(3);
550 
551  subelem[1]->set_node(0) = elem->node_ptr(3);
552  subelem[1]->set_node(1) = elem->node_ptr(4);
553  subelem[1]->set_node(2) = elem->node_ptr(2);
554  subelem[1]->set_node(3) = elem->node_ptr(5);
555 
556  subelem[2]->set_node(0) = elem->node_ptr(0);
557  subelem[2]->set_node(1) = elem->node_ptr(1);
558  subelem[2]->set_node(2) = elem->node_ptr(2);
559  subelem[2]->set_node(3) = elem->node_ptr(4);
560  }
561  }
562  else // Split on 1-3 diagonal
563  {
564  libmesh_assert (split_first_diagonal(elem, 1,3, 0,4));
565 
566  // Split on 0-5 diagonal
567  if (split_first_diagonal(elem, 0,5, 2,3))
568  {
569  // 1-3 and 0-5 split implies 1-5 split
570  libmesh_assert (split_first_diagonal(elem, 1,5, 2,4));
571 
572  subelem[0]->set_node(0) = elem->node_ptr(1);
573  subelem[0]->set_node(1) = elem->node_ptr(3);
574  subelem[0]->set_node(2) = elem->node_ptr(4);
575  subelem[0]->set_node(3) = elem->node_ptr(5);
576 
577  subelem[1]->set_node(0) = elem->node_ptr(1);
578  subelem[1]->set_node(1) = elem->node_ptr(0);
579  subelem[1]->set_node(2) = elem->node_ptr(3);
580  subelem[1]->set_node(3) = elem->node_ptr(5);
581 
582  subelem[2]->set_node(0) = elem->node_ptr(0);
583  subelem[2]->set_node(1) = elem->node_ptr(1);
584  subelem[2]->set_node(2) = elem->node_ptr(2);
585  subelem[2]->set_node(3) = elem->node_ptr(5);
586  }
587  else // Split on 2-3 diagonal
588  {
589  libmesh_assert (split_first_diagonal(elem, 2,3, 0,5));
590 
591  // Split on 1-5 diagonal
592  if (split_first_diagonal(elem, 1,5, 2,4))
593  {
594  subelem[0]->set_node(0) = elem->node_ptr(0);
595  subelem[0]->set_node(1) = elem->node_ptr(1);
596  subelem[0]->set_node(2) = elem->node_ptr(2);
597  subelem[0]->set_node(3) = elem->node_ptr(3);
598 
599  subelem[1]->set_node(0) = elem->node_ptr(3);
600  subelem[1]->set_node(1) = elem->node_ptr(1);
601  subelem[1]->set_node(2) = elem->node_ptr(2);
602  subelem[1]->set_node(3) = elem->node_ptr(5);
603 
604  subelem[2]->set_node(0) = elem->node_ptr(1);
605  subelem[2]->set_node(1) = elem->node_ptr(3);
606  subelem[2]->set_node(2) = elem->node_ptr(4);
607  subelem[2]->set_node(3) = elem->node_ptr(5);
608  }
609  else // Split on 2-4 diagonal
610  {
611  libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
612 
613  subelem[0]->set_node(0) = elem->node_ptr(0);
614  subelem[0]->set_node(1) = elem->node_ptr(1);
615  subelem[0]->set_node(2) = elem->node_ptr(2);
616  subelem[0]->set_node(3) = elem->node_ptr(3);
617 
618  subelem[1]->set_node(0) = elem->node_ptr(2);
619  subelem[1]->set_node(1) = elem->node_ptr(3);
620  subelem[1]->set_node(2) = elem->node_ptr(4);
621  subelem[1]->set_node(3) = elem->node_ptr(5);
622 
623  subelem[2]->set_node(0) = elem->node_ptr(3);
624  subelem[2]->set_node(1) = elem->node_ptr(1);
625  subelem[2]->set_node(2) = elem->node_ptr(2);
626  subelem[2]->set_node(3) = elem->node_ptr(4);
627  }
628  }
629  }
630 
631  break;
632  }
633 
634  case PRISM18:
635  {
636  subelem[0] = new Tet10;
637  subelem[1] = new Tet10;
638  subelem[2] = new Tet10;
639 
640  // Split on 0-4 diagonal
641  if (split_first_diagonal(elem, 0,4, 1,3))
642  {
643  // Split on 0-5 diagonal
644  if (split_first_diagonal(elem, 0,5, 2,3))
645  {
646  // Split on 1-5 diagonal
647  if (split_first_diagonal(elem, 1,5, 2,4))
648  {
649  subelem[0]->set_node(0) = elem->node_ptr(0);
650  subelem[0]->set_node(1) = elem->node_ptr(4);
651  subelem[0]->set_node(2) = elem->node_ptr(5);
652  subelem[0]->set_node(3) = elem->node_ptr(3);
653 
654  subelem[0]->set_node(4) = elem->node_ptr(15);
655  subelem[0]->set_node(5) = elem->node_ptr(13);
656  subelem[0]->set_node(6) = elem->node_ptr(17);
657  subelem[0]->set_node(7) = elem->node_ptr(9);
658  subelem[0]->set_node(8) = elem->node_ptr(12);
659  subelem[0]->set_node(9) = elem->node_ptr(14);
660 
661  subelem[1]->set_node(0) = elem->node_ptr(0);
662  subelem[1]->set_node(1) = elem->node_ptr(4);
663  subelem[1]->set_node(2) = elem->node_ptr(1);
664  subelem[1]->set_node(3) = elem->node_ptr(5);
665 
666  subelem[1]->set_node(4) = elem->node_ptr(15);
667  subelem[1]->set_node(5) = elem->node_ptr(10);
668  subelem[1]->set_node(6) = elem->node_ptr(6);
669  subelem[1]->set_node(7) = elem->node_ptr(17);
670  subelem[1]->set_node(8) = elem->node_ptr(13);
671  subelem[1]->set_node(9) = elem->node_ptr(16);
672 
673  subelem[2]->set_node(0) = elem->node_ptr(0);
674  subelem[2]->set_node(1) = elem->node_ptr(1);
675  subelem[2]->set_node(2) = elem->node_ptr(2);
676  subelem[2]->set_node(3) = elem->node_ptr(5);
677 
678  subelem[2]->set_node(4) = elem->node_ptr(6);
679  subelem[2]->set_node(5) = elem->node_ptr(7);
680  subelem[2]->set_node(6) = elem->node_ptr(8);
681  subelem[2]->set_node(7) = elem->node_ptr(17);
682  subelem[2]->set_node(8) = elem->node_ptr(16);
683  subelem[2]->set_node(9) = elem->node_ptr(11);
684  }
685  else // Split on 2-4 diagonal
686  {
687  libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
688 
689  subelem[0]->set_node(0) = elem->node_ptr(0);
690  subelem[0]->set_node(1) = elem->node_ptr(4);
691  subelem[0]->set_node(2) = elem->node_ptr(5);
692  subelem[0]->set_node(3) = elem->node_ptr(3);
693 
694  subelem[0]->set_node(4) = elem->node_ptr(15);
695  subelem[0]->set_node(5) = elem->node_ptr(13);
696  subelem[0]->set_node(6) = elem->node_ptr(17);
697  subelem[0]->set_node(7) = elem->node_ptr(9);
698  subelem[0]->set_node(8) = elem->node_ptr(12);
699  subelem[0]->set_node(9) = elem->node_ptr(14);
700 
701  subelem[1]->set_node(0) = elem->node_ptr(0);
702  subelem[1]->set_node(1) = elem->node_ptr(4);
703  subelem[1]->set_node(2) = elem->node_ptr(2);
704  subelem[1]->set_node(3) = elem->node_ptr(5);
705 
706  subelem[1]->set_node(4) = elem->node_ptr(15);
707  subelem[1]->set_node(5) = elem->node_ptr(16);
708  subelem[1]->set_node(6) = elem->node_ptr(8);
709  subelem[1]->set_node(7) = elem->node_ptr(17);
710  subelem[1]->set_node(8) = elem->node_ptr(13);
711  subelem[1]->set_node(9) = elem->node_ptr(11);
712 
713  subelem[2]->set_node(0) = elem->node_ptr(0);
714  subelem[2]->set_node(1) = elem->node_ptr(1);
715  subelem[2]->set_node(2) = elem->node_ptr(2);
716  subelem[2]->set_node(3) = elem->node_ptr(4);
717 
718  subelem[2]->set_node(4) = elem->node_ptr(6);
719  subelem[2]->set_node(5) = elem->node_ptr(7);
720  subelem[2]->set_node(6) = elem->node_ptr(8);
721  subelem[2]->set_node(7) = elem->node_ptr(15);
722  subelem[2]->set_node(8) = elem->node_ptr(10);
723  subelem[2]->set_node(9) = elem->node_ptr(16);
724  }
725  }
726  else // Split on 2-3 diagonal
727  {
728  libmesh_assert (split_first_diagonal(elem, 2,3, 0,5));
729 
730  // 0-4 and 2-3 split implies 2-4 split
731  libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
732 
733  subelem[0]->set_node(0) = elem->node_ptr(0);
734  subelem[0]->set_node(1) = elem->node_ptr(4);
735  subelem[0]->set_node(2) = elem->node_ptr(2);
736  subelem[0]->set_node(3) = elem->node_ptr(3);
737 
738  subelem[0]->set_node(4) = elem->node_ptr(15);
739  subelem[0]->set_node(5) = elem->node_ptr(16);
740  subelem[0]->set_node(6) = elem->node_ptr(8);
741  subelem[0]->set_node(7) = elem->node_ptr(9);
742  subelem[0]->set_node(8) = elem->node_ptr(12);
743  subelem[0]->set_node(9) = elem->node_ptr(17);
744 
745  subelem[1]->set_node(0) = elem->node_ptr(3);
746  subelem[1]->set_node(1) = elem->node_ptr(4);
747  subelem[1]->set_node(2) = elem->node_ptr(2);
748  subelem[1]->set_node(3) = elem->node_ptr(5);
749 
750  subelem[1]->set_node(4) = elem->node_ptr(12);
751  subelem[1]->set_node(5) = elem->node_ptr(16);
752  subelem[1]->set_node(6) = elem->node_ptr(17);
753  subelem[1]->set_node(7) = elem->node_ptr(14);
754  subelem[1]->set_node(8) = elem->node_ptr(13);
755  subelem[1]->set_node(9) = elem->node_ptr(11);
756 
757  subelem[2]->set_node(0) = elem->node_ptr(0);
758  subelem[2]->set_node(1) = elem->node_ptr(1);
759  subelem[2]->set_node(2) = elem->node_ptr(2);
760  subelem[2]->set_node(3) = elem->node_ptr(4);
761 
762  subelem[2]->set_node(4) = elem->node_ptr(6);
763  subelem[2]->set_node(5) = elem->node_ptr(7);
764  subelem[2]->set_node(6) = elem->node_ptr(8);
765  subelem[2]->set_node(7) = elem->node_ptr(15);
766  subelem[2]->set_node(8) = elem->node_ptr(10);
767  subelem[2]->set_node(9) = elem->node_ptr(16);
768  }
769  }
770  else // Split on 1-3 diagonal
771  {
772  libmesh_assert (split_first_diagonal(elem, 1,3, 0,4));
773 
774  // Split on 0-5 diagonal
775  if (split_first_diagonal(elem, 0,5, 2,3))
776  {
777  // 1-3 and 0-5 split implies 1-5 split
778  libmesh_assert (split_first_diagonal(elem, 1,5, 2,4));
779 
780  subelem[0]->set_node(0) = elem->node_ptr(1);
781  subelem[0]->set_node(1) = elem->node_ptr(3);
782  subelem[0]->set_node(2) = elem->node_ptr(4);
783  subelem[0]->set_node(3) = elem->node_ptr(5);
784 
785  subelem[0]->set_node(4) = elem->node_ptr(15);
786  subelem[0]->set_node(5) = elem->node_ptr(12);
787  subelem[0]->set_node(6) = elem->node_ptr(10);
788  subelem[0]->set_node(7) = elem->node_ptr(16);
789  subelem[0]->set_node(8) = elem->node_ptr(14);
790  subelem[0]->set_node(9) = elem->node_ptr(13);
791 
792  subelem[1]->set_node(0) = elem->node_ptr(1);
793  subelem[1]->set_node(1) = elem->node_ptr(0);
794  subelem[1]->set_node(2) = elem->node_ptr(3);
795  subelem[1]->set_node(3) = elem->node_ptr(5);
796 
797  subelem[1]->set_node(4) = elem->node_ptr(6);
798  subelem[1]->set_node(5) = elem->node_ptr(9);
799  subelem[1]->set_node(6) = elem->node_ptr(15);
800  subelem[1]->set_node(7) = elem->node_ptr(16);
801  subelem[1]->set_node(8) = elem->node_ptr(17);
802  subelem[1]->set_node(9) = elem->node_ptr(14);
803 
804  subelem[2]->set_node(0) = elem->node_ptr(0);
805  subelem[2]->set_node(1) = elem->node_ptr(1);
806  subelem[2]->set_node(2) = elem->node_ptr(2);
807  subelem[2]->set_node(3) = elem->node_ptr(5);
808 
809  subelem[2]->set_node(4) = elem->node_ptr(6);
810  subelem[2]->set_node(5) = elem->node_ptr(7);
811  subelem[2]->set_node(6) = elem->node_ptr(8);
812  subelem[2]->set_node(7) = elem->node_ptr(17);
813  subelem[2]->set_node(8) = elem->node_ptr(16);
814  subelem[2]->set_node(9) = elem->node_ptr(11);
815  }
816  else // Split on 2-3 diagonal
817  {
818  libmesh_assert (split_first_diagonal(elem, 2,3, 0,5));
819 
820  // Split on 1-5 diagonal
821  if (split_first_diagonal(elem, 1,5, 2,4))
822  {
823  subelem[0]->set_node(0) = elem->node_ptr(0);
824  subelem[0]->set_node(1) = elem->node_ptr(1);
825  subelem[0]->set_node(2) = elem->node_ptr(2);
826  subelem[0]->set_node(3) = elem->node_ptr(3);
827 
828  subelem[0]->set_node(4) = elem->node_ptr(6);
829  subelem[0]->set_node(5) = elem->node_ptr(7);
830  subelem[0]->set_node(6) = elem->node_ptr(8);
831  subelem[0]->set_node(7) = elem->node_ptr(9);
832  subelem[0]->set_node(8) = elem->node_ptr(15);
833  subelem[0]->set_node(9) = elem->node_ptr(17);
834 
835  subelem[1]->set_node(0) = elem->node_ptr(3);
836  subelem[1]->set_node(1) = elem->node_ptr(1);
837  subelem[1]->set_node(2) = elem->node_ptr(2);
838  subelem[1]->set_node(3) = elem->node_ptr(5);
839 
840  subelem[1]->set_node(4) = elem->node_ptr(15);
841  subelem[1]->set_node(5) = elem->node_ptr(7);
842  subelem[1]->set_node(6) = elem->node_ptr(17);
843  subelem[1]->set_node(7) = elem->node_ptr(14);
844  subelem[1]->set_node(8) = elem->node_ptr(16);
845  subelem[1]->set_node(9) = elem->node_ptr(11);
846 
847  subelem[2]->set_node(0) = elem->node_ptr(1);
848  subelem[2]->set_node(1) = elem->node_ptr(3);
849  subelem[2]->set_node(2) = elem->node_ptr(4);
850  subelem[2]->set_node(3) = elem->node_ptr(5);
851 
852  subelem[2]->set_node(4) = elem->node_ptr(15);
853  subelem[2]->set_node(5) = elem->node_ptr(12);
854  subelem[2]->set_node(6) = elem->node_ptr(10);
855  subelem[2]->set_node(7) = elem->node_ptr(16);
856  subelem[2]->set_node(8) = elem->node_ptr(14);
857  subelem[2]->set_node(9) = elem->node_ptr(13);
858  }
859  else // Split on 2-4 diagonal
860  {
861  libmesh_assert (split_first_diagonal(elem, 2,4, 1,5));
862 
863  subelem[0]->set_node(0) = elem->node_ptr(0);
864  subelem[0]->set_node(1) = elem->node_ptr(1);
865  subelem[0]->set_node(2) = elem->node_ptr(2);
866  subelem[0]->set_node(3) = elem->node_ptr(3);
867 
868  subelem[0]->set_node(4) = elem->node_ptr(6);
869  subelem[0]->set_node(5) = elem->node_ptr(7);
870  subelem[0]->set_node(6) = elem->node_ptr(8);
871  subelem[0]->set_node(7) = elem->node_ptr(9);
872  subelem[0]->set_node(8) = elem->node_ptr(15);
873  subelem[0]->set_node(9) = elem->node_ptr(17);
874 
875  subelem[1]->set_node(0) = elem->node_ptr(2);
876  subelem[1]->set_node(1) = elem->node_ptr(3);
877  subelem[1]->set_node(2) = elem->node_ptr(4);
878  subelem[1]->set_node(3) = elem->node_ptr(5);
879 
880  subelem[1]->set_node(4) = elem->node_ptr(17);
881  subelem[1]->set_node(5) = elem->node_ptr(12);
882  subelem[1]->set_node(6) = elem->node_ptr(16);
883  subelem[1]->set_node(7) = elem->node_ptr(11);
884  subelem[1]->set_node(8) = elem->node_ptr(14);
885  subelem[1]->set_node(9) = elem->node_ptr(13);
886 
887  subelem[2]->set_node(0) = elem->node_ptr(3);
888  subelem[2]->set_node(1) = elem->node_ptr(1);
889  subelem[2]->set_node(2) = elem->node_ptr(2);
890  subelem[2]->set_node(3) = elem->node_ptr(4);
891 
892  subelem[2]->set_node(4) = elem->node_ptr(15);
893  subelem[2]->set_node(5) = elem->node_ptr(7);
894  subelem[2]->set_node(6) = elem->node_ptr(17);
895  subelem[2]->set_node(7) = elem->node_ptr(12);
896  subelem[2]->set_node(8) = elem->node_ptr(10);
897  subelem[2]->set_node(9) = elem->node_ptr(16);
898  }
899  }
900  }
901 
902  break;
903  }
904 
905  // No need to split elements that are already simplicial:
906  case EDGE2:
907  case EDGE3:
908  case EDGE4:
909  case TRI3:
910  case TRI6:
911  case TET4:
912  case TET10:
913  case INFEDGE2:
914  // No way to split infinite quad/prism elements, so
915  // hopefully no need to
916  case INFQUAD4:
917  case INFQUAD6:
918  case INFPRISM6:
919  case INFPRISM12:
920  continue;
921  // If we're left with an unimplemented hex we're probably
922  // out of luck. TODO: implement hexes
923  default:
924  {
925  libMesh::err << "Error, encountered unimplemented element "
926  << Utility::enum_to_string<ElemType>(etype)
927  << " in MeshTools::Modification::all_tri()..."
928  << std::endl;
929  libmesh_not_implemented();
930  }
931  } // end switch (etype)
932 
933 
934 
935  // Be sure the correct IDs are also set for all subelems.
936  for (unsigned int i=0; i != max_subelems; ++i)
937  if (subelem[i]) {
938  subelem[i]->processor_id() = elem->processor_id();
939  subelem[i]->subdomain_id() = elem->subdomain_id();
940  }
941 
942  // On a mesh with boundary data, we need to move that data to
943  // the new elements.
944 
945  // On a mesh which is distributed, we need to move
946  // remote_elem links to the new elements.
947  bool mesh_is_serial = mesh.is_serial();
948 
949  if (mesh_has_boundary_data || mesh_is_serial)
950  {
951  // Container to key boundary IDs handed back by the BoundaryInfo object.
952  std::vector<boundary_id_type> bc_ids;
953 
954  for (auto sn : elem->side_index_range())
955  {
956  mesh.get_boundary_info().boundary_ids(elem, sn, bc_ids);
957  for (const auto & b_id : bc_ids)
958  {
959  if (mesh_is_serial && b_id == BoundaryInfo::invalid_id)
960  continue;
961 
962  // Make a sorted list of node ids for elem->side(sn)
963  std::unique_ptr<Elem> elem_side = elem->build_side_ptr(sn);
964  std::vector<dof_id_type> elem_side_nodes(elem_side->n_nodes());
965  for (unsigned int esn=0,
966  n_esn = cast_int<unsigned int>(elem_side_nodes.size());
967  esn != n_esn; ++esn)
968  elem_side_nodes[esn] = elem_side->node_id(esn);
969  std::sort(elem_side_nodes.begin(), elem_side_nodes.end());
970 
971  for (unsigned int i=0; i != max_subelems; ++i)
972  if (subelem[i])
973  {
974  for (auto subside : subelem[i]->side_index_range())
975  {
976  std::unique_ptr<Elem> subside_elem = subelem[i]->build_side_ptr(subside);
977 
978  // Make a list of *vertex* node ids for this subside, see if they are all present
979  // in elem->side(sn). Note 1: we can't just compare elem->key(sn) to
980  // subelem[i]->key(subside) in the Prism cases, since the new side is
981  // a different type. Note 2: we only use vertex nodes since, in the future,
982  // a Hex20 or Prism15's QUAD8 face may be split into two Tri6 faces, and the
983  // original face will not contain the mid-edge node.
984  std::vector<dof_id_type> subside_nodes(subside_elem->n_vertices());
985  for (unsigned int ssn=0,
986  n_ssn = cast_int<unsigned int>(subside_nodes.size());
987  ssn != n_ssn; ++ssn)
988  subside_nodes[ssn] = subside_elem->node_id(ssn);
989  std::sort(subside_nodes.begin(), subside_nodes.end());
990 
991  // std::includes returns true if every element of the second sorted range is
992  // contained in the first sorted range.
993  if (std::includes(elem_side_nodes.begin(), elem_side_nodes.end(),
994  subside_nodes.begin(), subside_nodes.end()))
995  {
996  if (b_id != BoundaryInfo::invalid_id)
997  {
998  new_bndry_ids.push_back(b_id);
999  new_bndry_elements.push_back(subelem[i]);
1000  new_bndry_sides.push_back(subside);
1001  }
1002 
1003  // If the original element had a RemoteElem neighbor on side 'sn',
1004  // then the subelem has one on side 'subside'.
1005  if (elem->neighbor_ptr(sn) == remote_elem)
1006  subelem[i]->set_neighbor(subside, const_cast<RemoteElem*>(remote_elem));
1007  }
1008  }
1009  }
1010  } // end for loop over boundary IDs
1011  } // end for loop over sides
1012 
1013  // Remove the original element from the BoundaryInfo structure.
1014  mesh.get_boundary_info().remove(elem);
1015 
1016  } // end if (mesh_has_boundary_data)
1017 
1018  // Determine new IDs for the split elements which will be
1019  // the same on all processors, therefore keeping the Mesh
1020  // in sync. Note: we offset the new IDs by max_orig_id to
1021  // avoid overwriting any of the original IDs.
1022  for (unsigned int i=0; i != max_subelems; ++i)
1023  if (subelem[i])
1024  {
1025  // Determine new IDs for the split elements which will be
1026  // the same on all processors, therefore keeping the Mesh
1027  // in sync. Note: we offset the new IDs by the max of the
1028  // pre-existing ids to avoid conflicting with originals.
1029  subelem[i]->set_id( max_orig_id + 6*elem->id() + i );
1030 
1031 #ifdef LIBMESH_ENABLE_UNIQUE_ID
1032  subelem[i]->set_unique_id() = max_unique_id + 2*elem->unique_id() + i;
1033 #endif
1034 
1035  // Prepare to add the newly-created simplices
1036  new_elements.push_back(subelem[i]);
1037  }
1038 
1039  // Delete the original element
1040  mesh.delete_elem(elem);
1041  } // End for loop over elements
1042  } // end scope
1043 
1044 
1045  // Now, iterate over the new elements vector, and add them each to
1046  // the Mesh.
1047  {
1048  std::vector<Elem *>::iterator el = new_elements.begin();
1049  const std::vector<Elem *>::iterator end = new_elements.end();
1050  for (; el != end; ++el)
1051  mesh.add_elem(*el);
1052  }
1053 
1054  if (mesh_has_boundary_data)
1055  {
1056  // If the old mesh had boundary data, the new mesh better have
1057  // some. However, we can't assert that the size of
1058  // new_bndry_elements vector is > 0, since we may not have split
1059  // any elements actually on the boundary. We also can't assert
1060  // that the original number of boundary sides is equal to the
1061  // sum of the boundary sides currently in the mesh and the
1062  // newly-added boundary sides, since in 3D, we may have split a
1063  // boundary QUAD into two boundary TRIs. Therefore, we won't be
1064  // too picky about the actual number of BCs, and just assert that
1065  // there are some, somewhere.
1066 #ifndef NDEBUG
1067  bool nbe_nonempty = new_bndry_elements.size();
1068  mesh.comm().max(nbe_nonempty);
1069  libmesh_assert(nbe_nonempty ||
1071 #endif
1072 
1073  // We should also be sure that the lengths of the new boundary data vectors
1074  // are all the same.
1075  libmesh_assert_equal_to (new_bndry_elements.size(), new_bndry_sides.size());
1076  libmesh_assert_equal_to (new_bndry_sides.size(), new_bndry_ids.size());
1077 
1078  // Add the new boundary info to the mesh
1079  for (std::size_t s=0; s<new_bndry_elements.size(); ++s)
1080  mesh.get_boundary_info().add_side(new_bndry_elements[s],
1081  new_bndry_sides[s],
1082  new_bndry_ids[s]);
1083  }
1084 
1085  // In a DistributedMesh any newly added ghost node ids may be
1086  // inconsistent, and unique_ids of newly added ghost nodes remain
1087  // unset.
1088  // make_nodes_parallel_consistent() will fix all this.
1089  if (!mesh.is_serial())
1090  {
1091  mesh.comm().max(added_new_ghost_point);
1092 
1093  if (added_new_ghost_point)
1095  }
1096 
1097 
1098 
1099  // Prepare the newly created mesh for use.
1100  mesh.prepare_for_use(/*skip_renumber =*/ false);
1101 
1102  // Let the new_elements and new_bndry_elements vectors go out of scope.
1103 }
1104 
1105 
1107  const unsigned int n_iterations,
1108  const Real power)
1109 {
1113  libmesh_assert_equal_to (mesh.mesh_dimension(), 2);
1114 
1115  /*
1116  * Create a quickly-searchable list of boundary nodes.
1117  */
1118  std::unordered_set<dof_id_type> boundary_node_ids =
1120 
1121  for (unsigned int iter=0; iter<n_iterations; iter++)
1122  {
1123  /*
1124  * loop over the mesh refinement level
1125  */
1126  unsigned int n_levels = MeshTools::n_levels(mesh);
1127  for (unsigned int refinement_level=0; refinement_level != n_levels;
1128  refinement_level++)
1129  {
1130  // initialize the storage (have to do it on every level to get empty vectors
1131  std::vector<Point> new_positions;
1132  std::vector<Real> weight;
1133  new_positions.resize(mesh.n_nodes());
1134  weight.resize(mesh.n_nodes());
1135 
1136  {
1137  // Loop over the elements to calculate new node positions
1138  for (const auto & elem : as_range(mesh.level_elements_begin(refinement_level),
1139  mesh.level_elements_end(refinement_level)))
1140  {
1141  /*
1142  * We relax all nodes on level 0 first
1143  * If the element is refined (level > 0), we interpolate the
1144  * parents nodes with help of the embedding matrix
1145  */
1146  if (refinement_level == 0)
1147  {
1148  for (auto s : elem->side_index_range())
1149  {
1150  /*
1151  * Only operate on sides which are on the
1152  * boundary or for which the current element's
1153  * id is greater than its neighbor's.
1154  * Sides get only built once.
1155  */
1156  if ((elem->neighbor_ptr(s) != nullptr) &&
1157  (elem->id() > elem->neighbor_ptr(s)->id()))
1158  {
1159  std::unique_ptr<const Elem> side(elem->build_side_ptr(s));
1160 
1161  const Node & node0 = side->node_ref(0);
1162  const Node & node1 = side->node_ref(1);
1163 
1164  Real node_weight = 1.;
1165  // calculate the weight of the nodes
1166  if (power > 0)
1167  {
1168  Point diff = node0-node1;
1169  node_weight = std::pow(diff.norm(), power);
1170  }
1171 
1172  const dof_id_type id0 = node0.id(), id1 = node1.id();
1173  new_positions[id0].add_scaled( node1, node_weight );
1174  new_positions[id1].add_scaled( node0, node_weight );
1175  weight[id0] += node_weight;
1176  weight[id1] += node_weight;
1177  }
1178  } // element neighbor loop
1179  }
1180 #ifdef LIBMESH_ENABLE_AMR
1181  else // refinement_level > 0
1182  {
1183  /*
1184  * Find the positions of the hanging nodes of refined elements.
1185  * We do this by calculating their position based on the parent
1186  * (one level less refined) element, and the embedding matrix
1187  */
1188 
1189  const Elem * parent = elem->parent();
1190 
1191  /*
1192  * find out which child I am
1193  */
1194  unsigned int c = parent->which_child_am_i(elem);
1195  /*
1196  *loop over the childs (that is, the current elements) nodes
1197  */
1198  for (auto nc : elem->node_index_range())
1199  {
1200  /*
1201  * the new position of the node
1202  */
1203  Point point;
1204  for (auto n : parent->node_index_range())
1205  {
1206  /*
1207  * The value from the embedding matrix
1208  */
1209  const float em_val = parent->embedding_matrix(c,nc,n);
1210 
1211  if (em_val != 0.)
1212  point.add_scaled (parent->point(n), em_val);
1213  }
1214 
1215  const dof_id_type id = elem->node_ptr(nc)->id();
1216  new_positions[id] = point;
1217  weight[id] = 1.;
1218  }
1219  } // if element refinement_level
1220 #endif // #ifdef LIBMESH_ENABLE_AMR
1221 
1222  } // element loop
1223 
1224  /*
1225  * finally reposition the vertex nodes
1226  */
1227  for (unsigned int nid=0; nid<mesh.n_nodes(); ++nid)
1228  if (!boundary_node_ids.count(nid) && weight[nid] > 0.)
1229  mesh.node_ref(nid) = new_positions[nid]/weight[nid];
1230  }
1231 
1232  // Now handle the additional second_order nodes by calculating
1233  // their position based on the vertex positions
1234  // we do a second loop over the level elements
1235  for (auto & elem : as_range(mesh.level_elements_begin(refinement_level),
1236  mesh.level_elements_end(refinement_level)))
1237  {
1238  const unsigned int son_begin = elem->n_vertices();
1239  const unsigned int son_end = elem->n_nodes();
1240  for (unsigned int n=son_begin; n<son_end; n++)
1241  {
1242  const unsigned int n_adjacent_vertices =
1244 
1245  Point point;
1246  for (unsigned int v=0; v<n_adjacent_vertices; v++)
1247  point.add(elem->point( elem->second_order_adjacent_vertex(n,v) ));
1248 
1249  const dof_id_type id = elem->node_ptr(n)->id();
1250  mesh.node_ref(id) = point/n_adjacent_vertices;
1251  }
1252  }
1253  } // refinement_level loop
1254  } // end iteration
1255 }
1256 
1257 
1258 
1259 #ifdef LIBMESH_ENABLE_AMR
1261 {
1262  // Algorithm:
1263  // .) For each active element in the mesh: construct a
1264  // copy which is the same in every way *except* it is
1265  // a level 0 element. Store the pointers to these in
1266  // a separate vector. Save any boundary information as well.
1267  // Delete the active element from the mesh.
1268  // .) Loop over all (remaining) elements in the mesh, delete them.
1269  // .) Add the level-0 copies back to the mesh
1270 
1271  // Temporary storage for new element pointers
1272  std::vector<Elem *> new_elements;
1273 
1274  // BoundaryInfo Storage for element ids, sides, and BC ids
1275  std::vector<Elem *> saved_boundary_elements;
1276  std::vector<boundary_id_type> saved_bc_ids;
1277  std::vector<unsigned short int> saved_bc_sides;
1278 
1279  // Container to catch boundary ids passed back by BoundaryInfo
1280  std::vector<boundary_id_type> bc_ids;
1281 
1282  // Reserve a reasonable amt. of space for each
1283  new_elements.reserve(mesh.n_active_elem());
1284  saved_boundary_elements.reserve(mesh.get_boundary_info().n_boundary_conds());
1285  saved_bc_ids.reserve(mesh.get_boundary_info().n_boundary_conds());
1286  saved_bc_sides.reserve(mesh.get_boundary_info().n_boundary_conds());
1287 
1288  for (auto & elem : mesh.active_element_ptr_range())
1289  {
1290  // Make a new element of the same type
1291  Elem * copy = Elem::build(elem->type()).release();
1292 
1293  // Set node pointers (they still point to nodes in the original mesh)
1294  for (auto n : elem->node_index_range())
1295  copy->set_node(n) = elem->node_ptr(n);
1296 
1297  // Copy over ids
1298  copy->processor_id() = elem->processor_id();
1299  copy->subdomain_id() = elem->subdomain_id();
1300 
1301  // Retain the original element's ID(s) as well, otherwise
1302  // the Mesh may try to create them for you...
1303  copy->set_id( elem->id() );
1304 #ifdef LIBMESH_ENABLE_UNIQUE_ID
1305  copy->set_unique_id() = elem->unique_id();
1306 #endif
1307 
1308  // This element could have boundary info or DistributedMesh
1309  // remote_elem links as well. We need to save the (elem,
1310  // side, bc_id) triples and those links
1311  for (auto s : elem->side_index_range())
1312  {
1313  if (elem->neighbor_ptr(s) == remote_elem)
1314  copy->set_neighbor(s, const_cast<RemoteElem *>(remote_elem));
1315 
1316  mesh.get_boundary_info().boundary_ids(elem, s, bc_ids);
1317  for (const auto & bc_id : bc_ids)
1319  {
1320  saved_boundary_elements.push_back(copy);
1321  saved_bc_ids.push_back(bc_id);
1322  saved_bc_sides.push_back(s);
1323  }
1324  }
1325 
1326 
1327  // We're done with this element
1328  mesh.delete_elem(elem);
1329 
1330  // But save the copy
1331  new_elements.push_back(copy);
1332  }
1333 
1334  // Make sure we saved the same number of boundary conditions
1335  // in each vector.
1336  libmesh_assert_equal_to (saved_boundary_elements.size(), saved_bc_ids.size());
1337  libmesh_assert_equal_to (saved_bc_ids.size(), saved_bc_sides.size());
1338 
1339  // Loop again, delete any remaining elements
1340  for (auto & elem : mesh.element_ptr_range())
1341  mesh.delete_elem(elem);
1342 
1343  // Add the copied (now level-0) elements back to the mesh
1344  for (auto & new_elem : new_elements)
1345  {
1346  // Save the original ID, because the act of adding the Elem can
1347  // change new_elem's id!
1348  dof_id_type orig_id = new_elem->id();
1349 
1350  Elem * added_elem = mesh.add_elem(new_elem);
1351 
1352  // If the Elem, as it was re-added to the mesh, now has a
1353  // different ID (this is unlikely, so it's just an assert)
1354  // the boundary information will no longer be correct.
1355  libmesh_assert_equal_to (orig_id, added_elem->id());
1356 
1357  // Avoid compiler warnings in opt mode.
1358  libmesh_ignore(added_elem, orig_id);
1359  }
1360 
1361  // Finally, also add back the saved boundary information
1362  for (std::size_t e=0; e<saved_boundary_elements.size(); ++e)
1363  mesh.get_boundary_info().add_side(saved_boundary_elements[e],
1364  saved_bc_sides[e],
1365  saved_bc_ids[e]);
1366 
1367  // Trim unused and renumber nodes and elements
1368  mesh.prepare_for_use(/*skip_renumber =*/ false);
1369 }
1370 #endif // #ifdef LIBMESH_ENABLE_AMR
1371 
1372 
1373 
1375  const boundary_id_type old_id,
1376  const boundary_id_type new_id)
1377 {
1378  if (old_id == new_id)
1379  {
1380  // If the IDs are the same, this is a no-op.
1381  return;
1382  }
1383 
1384  // A reference to the Mesh's BoundaryInfo object, for convenience.
1386 
1387  {
1388  // Temporary vector to hold ids
1389  std::vector<boundary_id_type> bndry_ids;
1390 
1391  // build_node_list returns a vector of (node, bc) tuples.
1392  for (const auto & t : bi.build_node_list())
1393  if (std::get<1>(t) == old_id)
1394  {
1395  // Get the node in question
1396  const Node * node = mesh.node_ptr(std::get<0>(t));
1397 
1398  // Get all the current IDs for this node.
1399  bi.boundary_ids(node, bndry_ids);
1400 
1401  // Update the IDs accordingly
1402  std::replace(bndry_ids.begin(), bndry_ids.end(), old_id, new_id);
1403 
1404  // Remove all traces of that node from the BoundaryInfo object
1405  bi.remove(node);
1406 
1407  // Add it back with the updated IDs
1408  bi.add_node(node, bndry_ids);
1409  }
1410  }
1411 
1412  {
1413  // Temporary vector to hold ids
1414  std::vector<boundary_id_type> bndry_ids;
1415 
1416  // build_edge_list returns a vector of (elem, side, bc) tuples.
1417  for (const auto & t : bi.build_edge_list())
1418  if (std::get<2>(t) == old_id)
1419  {
1420  // Get the elem in question
1421  const Elem * elem = mesh.elem_ptr(std::get<0>(t));
1422 
1423  // The edge of the elem in question
1424  unsigned short int edge = std::get<1>(t);
1425 
1426  // Get all the current IDs for the edge in question.
1427  bi.edge_boundary_ids(elem, edge, bndry_ids);
1428 
1429  // Update the IDs accordingly
1430  std::replace(bndry_ids.begin(), bndry_ids.end(), old_id, new_id);
1431 
1432  // Remove all traces of that edge from the BoundaryInfo object
1433  bi.remove_edge(elem, edge);
1434 
1435  // Add it back with the updated IDs
1436  bi.add_edge(elem, edge, bndry_ids);
1437  }
1438  }
1439 
1440  {
1441  // Temporary vector to hold ids
1442  std::vector<boundary_id_type> bndry_ids;
1443 
1444  // build_shellface_list returns a vector of (elem, side, bc) tuples.
1445  for (const auto & t : bi.build_shellface_list())
1446  if (std::get<2>(t) == old_id)
1447  {
1448  // Get the elem in question
1449  const Elem * elem = mesh.elem_ptr(std::get<0>(t));
1450 
1451  // The shellface of the elem in question
1452  unsigned short int shellface = std::get<1>(t);
1453 
1454  // Get all the current IDs for the shellface in question.
1455  bi.shellface_boundary_ids(elem, shellface, bndry_ids);
1456 
1457  // Update the IDs accordingly
1458  std::replace(bndry_ids.begin(), bndry_ids.end(), old_id, new_id);
1459 
1460  // Remove all traces of that shellface from the BoundaryInfo object
1461  bi.remove_shellface(elem, shellface);
1462 
1463  // Add it back with the updated IDs
1464  bi.add_shellface(elem, shellface, bndry_ids);
1465  }
1466  }
1467 
1468  {
1469  // Temporary vector to hold ids
1470  std::vector<boundary_id_type> bndry_ids;
1471 
1472  // build_side_list returns a vector of (elem, side, bc) tuples.
1473  for (const auto & t : bi.build_side_list())
1474  if (std::get<2>(t) == old_id)
1475  {
1476  // Get the elem in question
1477  const Elem * elem = mesh.elem_ptr(std::get<0>(t));
1478 
1479  // The side of the elem in question
1480  unsigned short int side = std::get<1>(t);
1481 
1482  // Get all the current IDs for the side in question.
1483  bi.boundary_ids(elem, side, bndry_ids);
1484 
1485  // Update the IDs accordingly
1486  std::replace(bndry_ids.begin(), bndry_ids.end(), old_id, new_id);
1487 
1488  // Remove all traces of that side from the BoundaryInfo object
1489  bi.remove_side(elem, side);
1490 
1491  // Add it back with the updated IDs
1492  bi.add_side(elem, side, bndry_ids);
1493  }
1494  }
1495 
1496  // Remove any remaining references to the old_id so it does not show
1497  // up in lists of boundary ids, etc.
1498  bi.remove_id(old_id);
1499 }
1500 
1501 
1502 
1504  const subdomain_id_type old_id,
1505  const subdomain_id_type new_id)
1506 {
1507  if (old_id == new_id)
1508  {
1509  // If the IDs are the same, this is a no-op.
1510  return;
1511  }
1512 
1513  for (auto & elem : mesh.element_ptr_range())
1514  {
1515  if (elem->subdomain_id() == old_id)
1516  elem->subdomain_id() = new_id;
1517  }
1518 }
1519 
1520 
1521 } // namespace libMesh
std::size_t n_boundary_conds() const
A 2D triangular element with 3 nodes.
Definition: face_tri3.h:56
unique_id_type & set_unique_id()
Definition: dof_object.h:685
const Elem * parent() const
Definition: elem.h:2479
virtual dof_id_type n_active_elem() const =0
virtual Node *& set_node(const unsigned int i)
Definition: elem.h:2024
void remove_edge(const Elem *elem, const unsigned short int edge)
A geometric point in (x,y,z) space associated with a DOF.
Definition: node.h:52
void rotate(MeshBase &mesh, const Real phi, const Real theta=0., const Real psi=0.)
virtual unique_id_type parallel_max_unique_id() const =0
virtual element_iterator level_elements_begin(unsigned int level)=0
Real norm() const
Definition: type_vector.h:912
IntRange< unsigned short > side_index_range() const
Definition: elem.h:2166
void add_scaled(const TypeVector< T2 > &, const T)
Definition: type_vector.h:627
void remove_id(boundary_id_type id)
void scale(MeshBase &mesh, const Real xs, const Real ys=0., const Real zs=0.)
void remove_shellface(const Elem *elem, const unsigned short int shellface)
void distort(MeshBase &mesh, const Real factor, const bool perturb_boundary=false)
void remove(const Node *node)
unsigned short int side
Definition: xdr_io.C:50
The base class for all geometric element types.
Definition: elem.h:100
MeshBase & mesh
virtual float embedding_matrix(const unsigned int child_num, const unsigned int child_node_num, const unsigned int parent_node_num) const =0
unique_id_type unique_id() const
Definition: dof_object.h:672
const Parallel::Communicator & comm() const
void shellface_boundary_ids(const Elem *const elem, const unsigned short int shellface, std::vector< boundary_id_type > &vec_to_fill) const
IterBase * end
virtual SimpleRange< element_iterator > active_element_ptr_range()=0
const BoundaryInfo & get_boundary_info() const
Definition: mesh_base.h:131
long double max(long double a, double b)
virtual Node * add_point(const Point &p, const dof_id_type id=DofObject::invalid_id, const processor_id_type proc_id=DofObject::invalid_processor_id)=0
boundary_id_type bc_id
Definition: xdr_io.C:51
Base class for Mesh.
Definition: mesh_base.h:77
void add(const TypeVector< T2 > &)
Definition: type_vector.h:603
dof_id_type & set_id()
Definition: dof_object.h:664
void build_side_list(std::vector< dof_id_type > &element_id_list, std::vector< unsigned short int > &side_list, std::vector< boundary_id_type > &bc_id_list) const
A 2D triangular element with 6 nodes.
Definition: face_tri6.h:56
std::vector< boundary_id_type > boundary_ids(const Node *node) const
void change_subdomain_id(MeshBase &mesh, const subdomain_id_type old_id, const subdomain_id_type new_id)
virtual std::unique_ptr< Elem > build_side_ptr(const unsigned int i, bool proxy=true)=0
virtual element_iterator level_elements_end(unsigned int level)=0
void smooth(MeshBase &, unsigned int, Real)
virtual bool is_serial() const
Definition: mesh_base.h:154
TypeVector< T > unit() const
Definition: type_vector.C:37
void libmesh_ignore(const Args &...)
void build_node_list(std::vector< dof_id_type > &node_id_list, std::vector< boundary_id_type > &bc_id_list) const
void add_node(const Node *node, const boundary_id_type id)
int8_t boundary_id_type
Definition: id_types.h:51
static const boundary_id_type invalid_id
virtual void delete_elem(Elem *e)=0
virtual SimpleRange< element_iterator > element_ptr_range()=0
dof_id_type id() const
Definition: dof_object.h:655
virtual Real hmin() const
Definition: elem.C:356
virtual unsigned int n_nodes() const =0
dof_id_type weight(const MeshBase &mesh, const processor_id_type pid)
Definition: mesh_tools.C:233
virtual Elem * add_elem(Elem *e)=0
static std::unique_ptr< Elem > build(const ElemType type, Elem *p=nullptr)
Definition: elem.C:245
void change_boundary_id(MeshBase &mesh, const boundary_id_type old_id, const boundary_id_type new_id)
unsigned int n_levels(const MeshBase &mesh)
Definition: mesh_tools.C:653
virtual const Node * query_node_ptr(const dof_id_type i) const =0
virtual SimpleRange< node_iterator > node_ptr_range()=0
virtual dof_id_type max_elem_id() const =0
Used by the Mesh to keep track of boundary nodes and elements.
Definition: boundary_info.h:57
SimpleRange< I > as_range(const std::pair< I, I > &p)
Definition: simple_range.h:57
static const dof_id_type invalid_id
Definition: dof_object.h:347
void prepare_for_use(const bool skip_renumber_nodes_and_elements=false, const bool skip_find_neighbors=false)
Definition: mesh_base.C:152
OStreamProxy err(std::cerr)
virtual unsigned int n_second_order_adjacent_vertices(const unsigned int n) const
Definition: elem.C:2558
void set_neighbor(const unsigned int i, Elem *n)
Definition: elem.h:2083
void build_shellface_list(std::vector< dof_id_type > &element_id_list, std::vector< unsigned short int > &shellface_list, std::vector< boundary_id_type > &bc_id_list) const
void translate(MeshBase &mesh, const Real xt=0., const Real yt=0., const Real zt=0.)
double pow(double a, int b)
unsigned int which_child_am_i(const Elem *e) const
Definition: elem.h:2620
SimpleRange< NodeRefIter > node_ref_range()
Definition: elem.h:2130
virtual const Elem * elem_ptr(const dof_id_type i) const =0
const Elem * neighbor_ptr(unsigned int i) const
Definition: elem.h:2050
void remove_side(const Elem *elem, const unsigned short int side)
virtual unsigned int n_vertices() const =0
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
subdomain_id_type subdomain_id() const
Definition: elem.h:2034
A 3D tetrahedral element with 10 nodes.
Definition: cell_tet10.h:60
const Node * node_ptr(const unsigned int i) const
Definition: elem.h:1957
void build_edge_list(std::vector< dof_id_type > &element_id_list, std::vector< unsigned short int > &edge_list, std::vector< boundary_id_type > &bc_id_list) const
void make_nodes_parallel_consistent(MeshBase &)
void add_side(const dof_id_type elem, const unsigned short int side, const boundary_id_type id)
virtual unsigned short int second_order_adjacent_vertex(const unsigned int n, const unsigned int v) const
Definition: elem.C:2566
void add_shellface(const dof_id_type elem, const unsigned short int shellface, const boundary_id_type id)
unsigned int mesh_dimension() const
Definition: mesh_base.C:126
void redistribute(MeshBase &mesh, const FunctionBase< Real > &mapfunc)
virtual std::unique_ptr< FunctionBase< Output > > clone() const =0
A 3D tetrahedral element with 4 nodes.
Definition: cell_tet4.h:53
IntRange< unsigned short > node_index_range() const
Definition: elem.h:2148
virtual const Node & node_ref(const dof_id_type i) const
Definition: mesh_base.h:434
virtual const Point & point(const dof_id_type i) const =0
virtual dof_id_type max_node_id() const =0
virtual dof_id_type n_elem() const =0
virtual const Node * node_ptr(const dof_id_type i) const =0
processor_id_type processor_id() const
processor_id_type processor_id() const
Definition: dof_object.h:717
virtual ElemType type() const =0
long double min(long double a, double b)
A geometric point in (x,y,z) space.
Definition: point.h:38
dof_id_type node_id(const unsigned int i) const
Definition: elem.h:1914
const Point & point(const unsigned int i) const
Definition: elem.h:1892
uint8_t unique_id_type
Definition: id_types.h:79
virtual dof_id_type n_nodes() const =0
void find_boundary_nodes(const MeshBase &mesh, std::vector< bool > &on_boundary)
Definition: mesh_tools.C:303
uint8_t dof_id_type
Definition: id_types.h:64
void add_edge(const dof_id_type elem, const unsigned short int edge, const boundary_id_type id)
const Real pi
Definition: libmesh.h:233
std::vector< boundary_id_type > edge_boundary_ids(const Elem *const elem, const unsigned short int edge) const
const RemoteElem * remote_elem
Definition: remote_elem.C:57