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