libMesh
mesh_generation.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2026 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 
20 // libmesh includes
21 #include "libmesh/mesh_generation.h"
22 #include "libmesh/unstructured_mesh.h"
23 #include "libmesh/mesh_refinement.h"
24 #include "libmesh/edge_edge2.h"
25 #include "libmesh/edge_edge3.h"
26 #include "libmesh/edge_edge4.h"
27 #include "libmesh/face_tri3.h"
28 #include "libmesh/face_tri6.h"
29 #include "libmesh/face_tri7.h"
30 #include "libmesh/face_quad4.h"
31 #include "libmesh/face_quad8.h"
32 #include "libmesh/face_quad9.h"
33 #include "libmesh/face_c0polygon.h"
34 #include "libmesh/cell_hex8.h"
35 #include "libmesh/cell_hex20.h"
36 #include "libmesh/cell_hex27.h"
37 #include "libmesh/cell_prism6.h"
38 #include "libmesh/cell_prism15.h"
39 #include "libmesh/cell_prism18.h"
40 #include "libmesh/cell_prism20.h"
41 #include "libmesh/cell_prism21.h"
42 #include "libmesh/cell_tet4.h"
43 #include "libmesh/cell_pyramid5.h"
44 #include "libmesh/libmesh_logging.h"
45 #include "libmesh/boundary_info.h"
46 #include "libmesh/remote_elem.h"
47 #include "libmesh/sphere.h"
48 #include "libmesh/mesh_modification.h"
49 #include "libmesh/mesh_smoother_laplace.h"
50 #include "libmesh/node_elem.h"
51 #include "libmesh/vector_value.h"
52 #include "libmesh/function_base.h"
53 #include "libmesh/enum_order.h"
54 #include "libmesh/int_range.h"
55 #include "libmesh/parallel.h"
56 #include "libmesh/parallel_ghost_sync.h"
57 #include "libmesh/enum_to_string.h"
58 
59 // C++ includes
60 #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
61 #include <cmath> // for std::sqrt
62 #include <unordered_set>
63 
64 
65 namespace libMesh
66 {
67 
68 namespace MeshTools {
69 namespace Generation {
70 namespace Private {
78 inline
79 unsigned int idx(const ElemType type,
80  const unsigned int nx,
81  const unsigned int i,
82  const unsigned int j)
83 {
84  switch(type)
85  {
86  case INVALID_ELEM:
87  case QUAD4:
88  case QUADSHELL4:
89  case TRI3:
90  case TRISHELL3:
91  {
92  return i + j*(nx+1);
93  }
94 
95  case QUAD8:
96  case QUADSHELL8:
97  case QUAD9:
98  case QUADSHELL9:
99  case TRI6:
100  case TRI7:
101  {
102  return i + j*(2*nx+1);
103  }
104 
105  default:
106  libmesh_error_msg("ERROR: Unrecognized 2D element type == " << Utility::enum_to_string(type));
107  }
108 
109  return libMesh::invalid_uint;
110 }
111 
112 
113 
114 // Same as the function above, but for 3D elements
115 inline
116 unsigned int idx(const ElemType type,
117  const unsigned int nx,
118  const unsigned int ny,
119  const unsigned int i,
120  const unsigned int j,
121  const unsigned int k)
122 {
123  switch(type)
124  {
125  case INVALID_ELEM:
126  case HEX8:
127  case PRISM6:
128  {
129  return i + (nx+1)*(j + k*(ny+1));
130  }
131 
132  case HEX20:
133  case HEX27:
134  case TET4: // TET4's are created from an initial HEX27 discretization
135  case TET10: // TET10's are created from an initial HEX27 discretization
136  case TET14: // TET14's are created from an initial HEX27 discretization
137  case PYRAMID5: // PYRAMID5's are created from an initial HEX27 discretization
138  case PYRAMID13:
139  case PYRAMID14:
140  case PYRAMID18:
141  case PRISM15:
142  case PRISM18:
143  case PRISM20:
144  case PRISM21:
145  {
146  return i + (2*nx+1)*(j + k*(2*ny+1));
147  }
148 
149  default:
150  libmesh_error_msg("ERROR: Unrecognized element type == " << Utility::enum_to_string(type));
151  }
152 
153  return libMesh::invalid_uint;
154 }
155 
156 
162 class GaussLobattoRedistributionFunction : public FunctionBase<Real>
163 {
164 public:
168  GaussLobattoRedistributionFunction(unsigned int nx,
169  Real xmin,
170  Real xmax,
171  unsigned int ny=0,
172  Real ymin=0,
173  Real ymax=0,
174  unsigned int nz=0,
175  Real zmin=0,
176  Real zmax=0) :
177  FunctionBase<Real>(nullptr)
178  {
179  _nelem.resize(3);
180  _nelem[0] = nx;
181  _nelem[1] = ny;
182  _nelem[2] = nz;
183 
184  _mins.resize(3);
185  _mins[0] = xmin;
186  _mins[1] = ymin;
187  _mins[2] = zmin;
188 
189  _widths.resize(3);
190  _widths[0] = xmax - xmin;
191  _widths[1] = ymax - ymin;
192  _widths[2] = zmax - zmin;
193 
194  // Precompute the cosine values.
195  _cosines.resize(3);
196  for (unsigned dir=0; dir<3; ++dir)
197  if (_nelem[dir] != 0)
198  {
199  _cosines[dir].resize(_nelem[dir]+1);
200  for (auto i : index_range(_cosines[dir]))
201  _cosines[dir][i] = std::cos(libMesh::pi * Real(i) / _nelem[dir]);
202  }
203  }
204 
208  GaussLobattoRedistributionFunction (GaussLobattoRedistributionFunction &&) = default;
209  GaussLobattoRedistributionFunction (const GaussLobattoRedistributionFunction &) = default;
210  GaussLobattoRedistributionFunction & operator= (const GaussLobattoRedistributionFunction &) = default;
211  GaussLobattoRedistributionFunction & operator= (GaussLobattoRedistributionFunction &&) = default;
212  virtual ~GaussLobattoRedistributionFunction () = default;
213 
218  virtual std::unique_ptr<FunctionBase<Real>> clone () const override
219  {
220  return std::make_unique<GaussLobattoRedistributionFunction>(*this);
221  }
222 
228  virtual void operator() (const Point & p,
229  const Real /*time*/,
230  DenseVector<Real> & output) override
231  {
232  output.resize(3);
233 
234  for (unsigned dir=0; dir<3; ++dir)
235  if (_nelem[dir] != 0)
236  {
237  // Figure out the index of the current point.
238  Real float_index = (p(dir) - _mins[dir]) * _nelem[dir] / _widths[dir];
239 
240  // std::modf separates the fractional and integer parts of the index.
241  Real integer_part_f = 0;
242  const Real fractional_part = std::modf(float_index, &integer_part_f);
243 
244  const int integer_part = int(integer_part_f);
245 
246  // Vertex node?
247  if (std::abs(fractional_part) < TOLERANCE || std::abs(fractional_part - 1.0) < TOLERANCE)
248  {
249  int index = int(round(float_index));
250 
251  // Move node to the Gauss-Lobatto position.
252  output(dir) = _mins[dir] + _widths[dir] * 0.5 * (1.0 - _cosines[dir][index]);
253  }
254 
255  // Mid-edge (quadratic) node?
256  else if (std::abs(fractional_part - 0.5) < TOLERANCE)
257  {
258  // Move node to the Gauss-Lobatto position, which is the average of
259  // the node to the left and the node to the right.
260  output(dir) = _mins[dir] + _widths[dir] * 0.5 *
261  (1.0 - 0.5*(_cosines[dir][integer_part] + _cosines[dir][integer_part+1]));
262  }
263 
264  // 1D only: Left interior (cubic) node?
265  else if (std::abs(fractional_part - 1./3.) < TOLERANCE)
266  {
267  // Move node to the Gauss-Lobatto position, which is
268  // 2/3*left_vertex + 1/3*right_vertex.
269  output(dir) = _mins[dir] + _widths[dir] * 0.5 *
270  (1.0 - 2./3.*_cosines[dir][integer_part] - 1./3.*_cosines[dir][integer_part+1]);
271  }
272 
273  // 1D only: Right interior (cubic) node?
274  else if (std::abs(fractional_part - 2./3.) < TOLERANCE)
275  {
276  // Move node to the Gauss-Lobatto position, which is
277  // 1/3*left_vertex + 2/3*right_vertex.
278  output(dir) = _mins[dir] + _widths[dir] * 0.5 *
279  (1.0 - 1./3.*_cosines[dir][integer_part] - 2./3.*_cosines[dir][integer_part+1]);
280  }
281 
282  else
283  libmesh_error_msg("Cannot redistribute node: " << p);
284  }
285  }
286 
291  virtual Real operator() (const Point & /*p*/,
292  const Real /*time*/) override
293  {
294  libmesh_not_implemented();
295  }
296 
297 protected:
298  // Stored data
299  std::vector<Real> _mins;
300  std::vector<unsigned int> _nelem;
301  std::vector<Real> _widths;
302 
303  // Precomputed values
304  std::vector<std::vector<Real>> _cosines;
305 };
306 
307 
308 } // namespace Private
309 } // namespace Generation
310 } // namespace MeshTools
311 
312 // ------------------------------------------------------------
313 // MeshTools::Generation function for mesh generation
314 void MeshTools::Generation::build_cube(UnstructuredMesh & mesh,
315  const unsigned int nx,
316  const unsigned int ny,
317  const unsigned int nz,
318  const Real xmin, const Real xmax,
319  const Real ymin, const Real ymax,
320  const Real zmin, const Real zmax,
321  const ElemType type,
322  const bool gauss_lobatto_grid)
323 {
324  LOG_SCOPE("build_cube()", "MeshTools::Generation");
325 
326  // Declare that we are using the indexing utility routine
327  // in the "Private" part of our current namespace. If this doesn't
328  // work in GCC 2.95.3 we can either remove it or stop supporting
329  // 2.95.3 altogether.
330  // Changing this to import the whole namespace... just importing idx
331  // causes an internal compiler error for Intel Compiler 11.0 on Linux
332  // in debug mode.
333  using namespace MeshTools::Generation::Private;
334 
335  // Clear the mesh and start from scratch
336  mesh.clear();
337 
338  BoundaryInfo & boundary_info = mesh.get_boundary_info();
339 
340  if (nz != 0)
341  {
342  mesh.set_mesh_dimension(3);
343  mesh.set_spatial_dimension(3);
344  }
345  else if (ny != 0)
346  {
347  mesh.set_mesh_dimension(2);
348  mesh.set_spatial_dimension(2);
349  }
350  else if (nx != 0)
351  {
352  mesh.set_mesh_dimension(1);
353  mesh.set_spatial_dimension(1);
354  }
355  else
356  {
357  // Will we get here?
358  mesh.set_mesh_dimension(0);
359  mesh.set_spatial_dimension(0);
360  }
361 
362  switch (mesh.mesh_dimension())
363  {
364  //---------------------------------------------------------------------
365  // Build a 0D point
366  case 0:
367  {
368  libmesh_assert_equal_to (nx, 0);
369  libmesh_assert_equal_to (ny, 0);
370  libmesh_assert_equal_to (nz, 0);
371 
372  libmesh_assert (type == INVALID_ELEM || type == NODEELEM);
373 
374  // Build one nodal element for the mesh
375  mesh.add_point (Point(0, 0, 0), 0);
376  Elem * elem = mesh.add_elem(Elem::build(NODEELEM));
377  elem->set_node(0, mesh.node_ptr(0));
378 
379  break;
380  }
381 
382 
383 
384  //---------------------------------------------------------------------
385  // Build a 1D line
386  case 1:
387  {
388  libmesh_assert_not_equal_to (nx, 0);
389  libmesh_assert_equal_to (ny, 0);
390  libmesh_assert_equal_to (nz, 0);
391  libmesh_assert_less (xmin, xmax);
392 
393  // Reserve elements
394  switch (type)
395  {
396  case INVALID_ELEM:
397  case EDGE2:
398  case EDGE3:
399  case EDGE4:
400  {
401  mesh.reserve_elem (nx);
402  break;
403  }
404 
405  default:
406  libmesh_error_msg("ERROR: Unrecognized 1D element type == " << Utility::enum_to_string(type));
407  }
408 
409  // Reserve nodes
410  switch (type)
411  {
412  case INVALID_ELEM:
413  case EDGE2:
414  {
415  mesh.reserve_nodes(nx+1);
416  break;
417  }
418 
419  case EDGE3:
420  {
421  mesh.reserve_nodes(2*nx+1);
422  break;
423  }
424 
425  case EDGE4:
426  {
427  mesh.reserve_nodes(3*nx+1);
428  break;
429  }
430 
431  default:
432  libmesh_error_msg("ERROR: Unrecognized 1D element type == " << Utility::enum_to_string(type));
433  }
434 
435 
436  // Build the nodes, depends on whether we're using linears,
437  // quadratics or cubics and whether using uniform grid or Gauss-Lobatto
438  unsigned int node_id = 0;
439  switch(type)
440  {
441  case INVALID_ELEM:
442  case EDGE2:
443  {
444  for (unsigned int i=0; i<=nx; i++)
445  {
446  const Node * const node = mesh.add_point (Point(static_cast<Real>(i)/nx, 0, 0), node_id++);
447  if (i == 0)
448  boundary_info.add_node(node, 0);
449  if (i == nx)
450  boundary_info.add_node(node, 1);
451  }
452 
453  break;
454  }
455 
456  case EDGE3:
457  {
458  for (unsigned int i=0; i<=2*nx; i++)
459  {
460  const Node * const node = mesh.add_point (Point(static_cast<Real>(i)/(2*nx), 0, 0), node_id++);
461  if (i == 0)
462  boundary_info.add_node(node, 0);
463  if (i == 2*nx)
464  boundary_info.add_node(node, 1);
465  }
466  break;
467  }
468 
469  case EDGE4:
470  {
471  for (unsigned int i=0; i<=3*nx; i++)
472  {
473  const Node * const node = mesh.add_point (Point(static_cast<Real>(i)/(3*nx), 0, 0), node_id++);
474  if (i == 0)
475  boundary_info.add_node(node, 0);
476  if (i == 3*nx)
477  boundary_info.add_node(node, 1);
478  }
479 
480  break;
481  }
482 
483  default:
484  libmesh_error_msg("ERROR: Unrecognized 1D element type == " << Utility::enum_to_string(type));
485 
486  }
487 
488  // Build the elements of the mesh
489  switch(type)
490  {
491  case INVALID_ELEM:
492  case EDGE2:
493  {
494  for (unsigned int i=0; i<nx; i++)
495  {
496  Elem * elem = mesh.add_elem(Elem::build_with_id(EDGE2, i));
497  elem->set_node(0, mesh.node_ptr(i));
498  elem->set_node(1, mesh.node_ptr(i+1));
499 
500  if (i == 0)
501  boundary_info.add_side(elem, 0, 0);
502 
503  if (i == (nx-1))
504  boundary_info.add_side(elem, 1, 1);
505 
506  }
507  break;
508  }
509 
510  case EDGE3:
511  {
512  for (unsigned int i=0; i<nx; i++)
513  {
514  Elem * elem = mesh.add_elem(Elem::build_with_id(EDGE3, i));
515  elem->set_node(0, mesh.node_ptr(2*i));
516  elem->set_node(2, mesh.node_ptr(2*i+1));
517  elem->set_node(1, mesh.node_ptr(2*i+2));
518 
519  if (i == 0)
520  boundary_info.add_side(elem, 0, 0);
521 
522  if (i == (nx-1))
523  boundary_info.add_side(elem, 1, 1);
524  }
525  break;
526  }
527 
528  case EDGE4:
529  {
530  for (unsigned int i=0; i<nx; i++)
531  {
532  Elem * elem = mesh.add_elem(Elem::build_with_id(EDGE4, i));
533  elem->set_node(0, mesh.node_ptr(3*i));
534  elem->set_node(2, mesh.node_ptr(3*i+1));
535  elem->set_node(3, mesh.node_ptr(3*i+2));
536  elem->set_node(1, mesh.node_ptr(3*i+3));
537 
538  if (i == 0)
539  boundary_info.add_side(elem, 0, 0);
540 
541  if (i == (nx-1))
542  boundary_info.add_side(elem, 1, 1);
543  }
544  break;
545  }
546 
547  default:
548  libmesh_error_msg("ERROR: Unrecognized 1D element type == " << Utility::enum_to_string(type));
549  }
550 
551  // Move the nodes to their final locations.
552  if (gauss_lobatto_grid)
553  {
554  GaussLobattoRedistributionFunction func(nx, xmin, xmax);
555  MeshTools::Modification::redistribute(mesh, func);
556  }
557  else // !gauss_lobatto_grid
558  {
559  for (Node * node : mesh.node_ptr_range())
560  (*node)(0) = (*node)(0)*(xmax-xmin) + xmin;
561  }
562 
563  // Add sideset names to boundary info
564  boundary_info.sideset_name(0) = "left";
565  boundary_info.sideset_name(1) = "right";
566 
567  // Add nodeset names to boundary info
568  boundary_info.nodeset_name(0) = "left";
569  boundary_info.nodeset_name(1) = "right";
570 
571  break;
572  }
573 
574 
575 
576 
577 
578 
579 
580 
581 
582 
583  //---------------------------------------------------------------------
584  // Build a 2D quadrilateral
585  case 2:
586  {
587  libmesh_assert_not_equal_to (nx, 0);
588  libmesh_assert_not_equal_to (ny, 0);
589  libmesh_assert_equal_to (nz, 0);
590  libmesh_assert_less (xmin, xmax);
591  libmesh_assert_less (ymin, ymax);
592 
593  // Reserve elements. The TRI3 and TRI6 meshes
594  // have twice as many elements...
595  switch (type)
596  {
597  case INVALID_ELEM:
598  case QUAD4:
599  case QUADSHELL4:
600  case QUAD8:
601  case QUADSHELL8:
602  case QUAD9:
603  case QUADSHELL9:
604  {
605  mesh.reserve_elem (nx*ny);
606  break;
607  }
608 
609  case TRI3:
610  case TRISHELL3:
611  case TRI6:
612  case TRI7:
613  {
614  mesh.reserve_elem (2*nx*ny);
615  break;
616  }
617 
618  case C0POLYGON:
619  {
620  mesh.reserve_elem ((nx + 1) * (ny + 1));
621  break;
622  }
623 
624  default:
625  libmesh_error_msg("ERROR: Unrecognized 2D element type == " << Utility::enum_to_string(type));
626  }
627 
628 
629 
630  // Reserve nodes. The quadratic element types
631  // need to reserve more nodes than the linear types.
632  switch (type)
633  {
634  case INVALID_ELEM:
635  case QUAD4:
636  case QUADSHELL4:
637  case TRI3:
638  case TRISHELL3:
639  {
640  mesh.reserve_nodes( (nx+1)*(ny+1) );
641  break;
642  }
643 
644  case QUAD8:
645  case QUADSHELL8:
646  case QUAD9:
647  case QUADSHELL9:
648  case TRI6:
649  {
650  mesh.reserve_nodes( (2*nx+1)*(2*ny+1) );
651  break;
652  }
653 
654  case TRI7:
655  {
656  mesh.reserve_nodes( (2*nx+1)*(2*ny+1) + 2*nx*ny );
657  break;
658  }
659  case C0POLYGON:
660  {
661  mesh.reserve_nodes (4 + 3*nx*ny + 2*nx + 2*ny);
662  break;
663  }
664 
665  default:
666  libmesh_error_msg("ERROR: Unrecognized 2D element type == " << Utility::enum_to_string(type));
667  }
668 
669 
670 
671  // Build the nodes. Depends on whether you are using a linear
672  // or quadratic element, and whether you are using a uniform
673  // grid or the Gauss-Lobatto grid points.
674  unsigned int node_id = 0;
675  switch (type)
676  {
677  case INVALID_ELEM:
678  case QUAD4:
679  case QUADSHELL4:
680  case TRI3:
681  case TRISHELL3:
682  {
683  for (unsigned int j=0; j<=ny; j++)
684  for (unsigned int i=0; i<=nx; i++)
685  {
686  const Node * const node =
687  mesh.add_point(Point(static_cast<Real>(i) / static_cast<Real>(nx),
688  static_cast<Real>(j) / static_cast<Real>(ny),
689  0.),
690  node_id++);
691  if (j == 0)
692  boundary_info.add_node(node, 0);
693  if (j == ny)
694  boundary_info.add_node(node, 2);
695  if (i == 0)
696  boundary_info.add_node(node, 3);
697  if (i == nx)
698  boundary_info.add_node(node, 1);
699  }
700 
701  break;
702  }
703 
704  case QUAD8:
705  case QUADSHELL8:
706  case QUAD9:
707  case QUADSHELL9:
708  case TRI6:
709  case TRI7:
710  {
711  for (unsigned int j=0; j<=(2*ny); j++)
712  for (unsigned int i=0; i<=(2*nx); i++)
713  {
714  const Node * const node =
715  mesh.add_point(Point(static_cast<Real>(i) / static_cast<Real>(2 * nx),
716  static_cast<Real>(j) / static_cast<Real>(2 * ny),
717  0),
718  node_id++);
719  if (j == 0)
720  boundary_info.add_node(node, 0);
721  if (j == 2*ny)
722  boundary_info.add_node(node, 2);
723  if (i == 0)
724  boundary_info.add_node(node, 3);
725  if (i == 2*nx)
726  boundary_info.add_node(node, 1);
727  }
728 
729  // We'll add any interior Tri7 nodes last, to keep from
730  // messing with our idx function
731  if (type == TRI7)
732  for (unsigned int j=0; j<(3*ny); j += 3)
733  for (unsigned int i=0; i<(3*nx); i += 3)
734  {
735  // The bottom-right triangle's center node
736  mesh.add_point(Point(static_cast<Real>(i+2) / static_cast<Real>(3 * nx),
737  static_cast<Real>(j+1) / static_cast<Real>(3 * ny),
738  0),
739  node_id++);
740  // The top-left triangle's center node
741  mesh.add_point(Point(static_cast<Real>(i+1) / static_cast<Real>(3 * nx),
742  static_cast<Real>(j+2) / static_cast<Real>(3 * ny),
743  0),
744  node_id++);
745  }
746 
747  break;
748  }
749 
750  case C0POLYGON:
751  {
752  // we create the nodes at the same time as the elements
753  break;
754  }
755 
756  default:
757  libmesh_error_msg("ERROR: Unrecognized 2D element type == " << Utility::enum_to_string(type));
758  }
759 
760 
761 
762 
763 
764 
765  // Build the elements. Each one is a bit different.
766  unsigned int elem_id = 0;
767  switch (type)
768  {
769 
770  case INVALID_ELEM:
771  case QUAD4:
772  case QUADSHELL4:
773  {
774  for (unsigned int j=0; j<ny; j++)
775  for (unsigned int i=0; i<nx; i++)
776  {
777  Elem * elem = mesh.add_elem(Elem::build_with_id(type == INVALID_ELEM ? QUAD4 : type, elem_id++));
778  elem->set_node(0, mesh.node_ptr(idx(type,nx,i,j) ));
779  elem->set_node(1, mesh.node_ptr(idx(type,nx,i+1,j) ));
780  elem->set_node(2, mesh.node_ptr(idx(type,nx,i+1,j+1)));
781  elem->set_node(3, mesh.node_ptr(idx(type,nx,i,j+1) ));
782 
783  if (j == 0)
784  boundary_info.add_side(elem, 0, 0);
785 
786  if (j == (ny-1))
787  boundary_info.add_side(elem, 2, 2);
788 
789  if (i == 0)
790  boundary_info.add_side(elem, 3, 3);
791 
792  if (i == (nx-1))
793  boundary_info.add_side(elem, 1, 1);
794  }
795  break;
796  }
797 
798 
799  case TRI3:
800  case TRISHELL3:
801  {
802  for (unsigned int j=0; j<ny; j++)
803  for (unsigned int i=0; i<nx; i++)
804  {
805  // Add first Tri3
806  Elem * elem = mesh.add_elem(Elem::build_with_id(type, elem_id++));
807  elem->set_node(0, mesh.node_ptr(idx(type,nx,i,j) ));
808  elem->set_node(1, mesh.node_ptr(idx(type,nx,i+1,j) ));
809  elem->set_node(2, mesh.node_ptr(idx(type,nx,i+1,j+1)));
810 
811  if (j == 0)
812  boundary_info.add_side(elem, 0, 0);
813 
814  if (i == (nx-1))
815  boundary_info.add_side(elem, 1, 1);
816 
817  // Add second Tri3
818  elem = mesh.add_elem(Elem::build_with_id(type, elem_id++));
819  elem->set_node(0, mesh.node_ptr(idx(type,nx,i,j) ));
820  elem->set_node(1, mesh.node_ptr(idx(type,nx,i+1,j+1)));
821  elem->set_node(2, mesh.node_ptr(idx(type,nx,i,j+1) ));
822 
823  if (j == (ny-1))
824  boundary_info.add_side(elem, 1, 2);
825 
826  if (i == 0)
827  boundary_info.add_side(elem, 2, 3);
828  }
829  break;
830  }
831 
832 
833 
834  case QUAD8:
835  case QUADSHELL8:
836  case QUAD9:
837  case QUADSHELL9:
838  {
839  for (unsigned int j=0; j<(2*ny); j += 2)
840  for (unsigned int i=0; i<(2*nx); i += 2)
841  {
842  Elem * elem = mesh.add_elem(Elem::build_with_id(type, elem_id++));
843  elem->set_node(0, mesh.node_ptr(idx(type,nx,i,j) ));
844  elem->set_node(1, mesh.node_ptr(idx(type,nx,i+2,j) ));
845  elem->set_node(2, mesh.node_ptr(idx(type,nx,i+2,j+2)));
846  elem->set_node(3, mesh.node_ptr(idx(type,nx,i,j+2) ));
847  elem->set_node(4, mesh.node_ptr(idx(type,nx,i+1,j) ));
848  elem->set_node(5, mesh.node_ptr(idx(type,nx,i+2,j+1)));
849  elem->set_node(6, mesh.node_ptr(idx(type,nx,i+1,j+2)));
850  elem->set_node(7, mesh.node_ptr(idx(type,nx,i,j+1) ));
851 
852  if (type == QUAD9 || type == QUADSHELL9)
853  elem->set_node(8, mesh.node_ptr(idx(type,nx,i+1,j+1)));
854 
855  if (j == 0)
856  boundary_info.add_side(elem, 0, 0);
857 
858  if (j == 2*(ny-1))
859  boundary_info.add_side(elem, 2, 2);
860 
861  if (i == 0)
862  boundary_info.add_side(elem, 3, 3);
863 
864  if (i == 2*(nx-1))
865  boundary_info.add_side(elem, 1, 1);
866  }
867  break;
868  }
869 
870 
871  case TRI6:
872  case TRI7:
873  {
874  for (unsigned int j=0; j<(2*ny); j += 2)
875  for (unsigned int i=0; i<(2*nx); i += 2)
876  {
877  // Add first Tri in the bottom-right of its quad
878  Elem * elem = mesh.add_elem(Elem::build_with_id(type, elem_id++));
879  elem->set_node(0, mesh.node_ptr(idx(type,nx,i,j) ));
880  elem->set_node(1, mesh.node_ptr(idx(type,nx,i+2,j) ));
881  elem->set_node(2, mesh.node_ptr(idx(type,nx,i+2,j+2)));
882  elem->set_node(3, mesh.node_ptr(idx(type,nx,i+1,j) ));
883  elem->set_node(4, mesh.node_ptr(idx(type,nx,i+2,j+1)));
884  elem->set_node(5, mesh.node_ptr(idx(type,nx,i+1,j+1)));
885 
886  if (type == TRI7)
887  elem->set_node(6, mesh.node_ptr(elem->id()+(2*nx+1)*(2*ny+1)));
888 
889  if (j == 0)
890  boundary_info.add_side(elem, 0, 0);
891 
892  if (i == 2*(nx-1))
893  boundary_info.add_side(elem, 1, 1);
894 
895  // Add second Tri in the top left of its quad
896  elem = mesh.add_elem(Elem::build_with_id(type, elem_id++));
897  elem->set_node(0, mesh.node_ptr(idx(type,nx,i,j) ));
898  elem->set_node(1, mesh.node_ptr(idx(type,nx,i+2,j+2)));
899  elem->set_node(2, mesh.node_ptr(idx(type,nx,i,j+2) ));
900  elem->set_node(3, mesh.node_ptr(idx(type,nx,i+1,j+1)));
901  elem->set_node(4, mesh.node_ptr(idx(type,nx,i+1,j+2)));
902  elem->set_node(5, mesh.node_ptr(idx(type,nx,i,j+1) ));
903 
904  if (type == TRI7)
905  elem->set_node(6, mesh.node_ptr(elem->id()+(2*nx+1)*(2*ny+1)));
906 
907  if (j == 2*(ny-1))
908  boundary_info.add_side(elem, 1, 2);
909 
910  if (i == 0)
911  boundary_info.add_side(elem, 2, 3);
912  }
913  break;
914  };
915 
916  case C0POLYGON:
917  {
918  // Build a 2D paving using hexagons (center), quads (part of y-boundary)
919  // and triangles (x-boundaries).
920  // Vector to re-use previously created nodes
921  std::vector<Node *> node_list;
922 
923  // Start with a layer of triangles on the boundary
924  const auto dx_tri = (xmax - xmin) / (nx);
925  const auto dy_tri = (ymax - ymin) / (ny + 1);
926  std::unique_ptr<Elem> new_elem;
927  for (const auto i : make_range(nx + 1))
928  {
929  // Make new nodes for bottom layer of triangles
930  Node *node0, *node1, *node2;
931  if (i == 0)
932  {
933  node0 = mesh.add_point(Point(0., 0, 0.));
934  node1 = mesh.add_point(Point(0., dy_tri / 2., 0.));
935  node2 = mesh.add_point(Point(dx_tri / 2., 0., 0.));
936  node_list.push_back(node0);
937  node_list.push_back(node1);
938  node_list.push_back(node2);
939  }
940  else if (i < nx)
941  {
942  node0 = node_list.back();
943  node1 = mesh.add_point(Point((i)*dx_tri, dy_tri / 2., 0.));
944  node2 = mesh.add_point(Point((i + 1. / 2.) * dx_tri, 0., 0.));
945  node_list.push_back(node1);
946  node_list.push_back(node2);
947  }
948  else
949  {
950  node0 = node_list.back();
951  node1 = mesh.add_point(Point((i)*dx_tri, dy_tri / 2., 0.));
952  node2 = mesh.add_point(Point((i)*dx_tri, 0., 0.));
953  node_list.push_back(node1);
954  node_list.push_back(node2);
955  }
956 
957  new_elem = std::make_unique<C0Polygon>(3);
958  // Switch to Tri3 when exodus default output supports element type mixes
959  new_elem->set_node(0, node0);
960  new_elem->set_node(1, node1);
961  new_elem->set_node(2, node2);
962  auto * elem = mesh.add_elem(std::move(new_elem));
963 
964  // Set boundaries
965  if (i == 0)
966  boundary_info.add_side(elem, 0, 3); // left
967  else if (i == nx)
968  boundary_info.add_side(elem, 1, 1); // right
969  boundary_info.add_side(elem, 2, 0); // bottom
970  }
971  // Start with the second node to build hexagons
972  unsigned int running_index = 1;
973 
974  // Build layers of hexagons
975  const auto hex_side =
976  (ymax - ymin - (ny == 1 ? dy_tri : (1. + (ny - 1) / 2.) * dy_tri)) / ny;
977  for (const auto j : make_range(ny))
978  {
979  for (const auto i : make_range(nx + (j % 2)))
980  {
981  if ((j % 2 == 0) || ((i > 0) && (i < nx)))
982  {
983  Node *n0, *n1, *n2, *n3, *n4, *n5;
984  n0 = node_list[running_index++];
985  n1 = node_list[running_index++];
986  n2 = node_list[running_index];
987 
988  if (i == 0)
989  {
990  n3 = mesh.add_point(Point(*n0) + RealVectorValue(0, hex_side, 0));
991  node_list.push_back(n3);
992  }
993  else
994  n3 = node_list.back();
995 
996  n4 = mesh.add_point(Point(*n1) + RealVectorValue(0, hex_side + dy_tri, 0));
997  n5 = mesh.add_point(Point(*n2) + RealVectorValue(0, hex_side, 0));
998  node_list.push_back(n4);
999  node_list.push_back(n5);
1000 
1001  new_elem = std::make_unique<libMesh::C0Polygon>(6);
1002  new_elem->set_node(0, n0);
1003  new_elem->set_node(1, n1);
1004  new_elem->set_node(2, n2);
1005  new_elem->set_node(3, n5);
1006  new_elem->set_node(4, n4);
1007  new_elem->set_node(5, n3);
1008  auto * elem = mesh.add_elem(std::move(new_elem));
1009 
1010  // Set boundaries
1011  if (i == 0)
1012  boundary_info.add_side(elem, 5, 3); // left
1013  else if (i == nx)
1014  boundary_info.add_side(elem, 2, 1); // right
1015  }
1016  // The hexagons are offset, so we build on a quad on each external side to fill
1017  else if (i == 0 || i == nx)
1018  {
1019  Node *n0, *n1, *n2, *n3;
1020  n0 = node_list[running_index++];
1021  n1 = node_list[running_index];
1022 
1023  if (i == 0)
1024  {
1025  n2 = mesh.add_point(Point(*n0) + RealVectorValue(0, hex_side + dy_tri, 0));
1026  node_list.push_back(n2);
1027  n3 = mesh.add_point(Point(*n1) + RealVectorValue(0, hex_side, 0));
1028  }
1029  else
1030  {
1031  n2 = node_list.back();
1032  n3 = mesh.add_point(Point(*n1) + RealVectorValue(0, hex_side + dy_tri, 0));
1033  }
1034  node_list.push_back(n3);
1035 
1036  new_elem = std::make_unique<C0Polygon>(4);
1037  // Switch to Quad4 when exodus default output supports element type mixes
1038  new_elem->set_node(0, n0);
1039  new_elem->set_node(1, n1);
1040  new_elem->set_node(3, n2);
1041  new_elem->set_node(2, n3);
1042  auto * elem = mesh.add_elem(std::move(new_elem));
1043 
1044  // Set boundaries
1045  if (i == 0)
1046  boundary_info.add_side(elem, 3, 3); // left
1047  else if (i == nx)
1048  boundary_info.add_side(elem, 1, 1); // right
1049  }
1050  else
1051  libmesh_assert(false);
1052  }
1053  // Increment once to switch to next 'row' of nodes
1054  running_index++;
1055 
1056  // Skip lower right corner node
1057  if (j == 0)
1058  running_index++;
1059  }
1060 
1061  // Build a final layer of triangles
1062  const bool ny_odd = (ny % 2 == 1);
1063  for (const auto i : make_range(nx + ny_odd))
1064  {
1065  // Use existing nodes, except at the corners
1066  Node *node0, *node1, *node2;
1067  if (i == 0 && ny_odd)
1068  {
1069  node0 = mesh.add_point(Point(0., ymax, 0.));
1070  node1 = node_list[running_index++];
1071  node2 = node_list[running_index];
1072  }
1073  else if (i < nx)
1074  {
1075  node0 = node_list[running_index++];
1076  node1 = node_list[running_index++];
1077  node2 = node_list[running_index];
1078  }
1079  // This case only reached if ny is odd and we are using a triangle in top right corner
1080  else
1081  {
1082  node0 = node_list[running_index++];
1083  node1 = node_list[running_index];
1084  node2 = mesh.add_point(Point(xmax, ymax, 0.));
1085  }
1086 
1087  new_elem = std::make_unique<C0Polygon>(3);
1088  // Switch to Tri3 when exodus default output supports element type mixes
1089  new_elem->set_node(0, node0);
1090  new_elem->set_node(1, node1);
1091  new_elem->set_node(2, node2);
1092  auto * elem = mesh.add_elem(std::move(new_elem));
1093 
1094  // Set boundaries
1095  if (i == 0)
1096  boundary_info.add_side(elem, 0, 3); // left
1097  else if (i == nx)
1098  boundary_info.add_side(elem, 1, 1); // right
1099  boundary_info.add_side(elem, 2, 2); // top
1100 
1101  }
1102  break;
1103  }
1104 
1105 
1106  default:
1107  libmesh_error_msg("ERROR: Unrecognized 2D element type == " << Utility::enum_to_string(type));
1108  }
1109 
1110 
1111 
1112 
1113  // Scale the nodal positions
1114  if (gauss_lobatto_grid)
1115  {
1116  GaussLobattoRedistributionFunction func(nx, xmin, xmax,
1117  ny, ymin, ymax);
1118  MeshTools::Modification::redistribute(mesh, func);
1119  }
1120  else // !gauss_lobatto_grid
1121  {
1122  for (Node * node : mesh.node_ptr_range())
1123  {
1124  (*node)(0) = ((*node)(0))*(xmax-xmin) + xmin;
1125  (*node)(1) = ((*node)(1))*(ymax-ymin) + ymin;
1126  }
1127  }
1128 
1129  // Add sideset names to boundary info
1130  boundary_info.sideset_name(0) = "bottom";
1131  boundary_info.sideset_name(1) = "right";
1132  boundary_info.sideset_name(2) = "top";
1133  boundary_info.sideset_name(3) = "left";
1134 
1135  // Add nodeset names to boundary info
1136  boundary_info.nodeset_name(0) = "bottom";
1137  boundary_info.nodeset_name(1) = "right";
1138  boundary_info.nodeset_name(2) = "top";
1139  boundary_info.nodeset_name(3) = "left";
1140 
1141  break;
1142  }
1143 
1144 
1145 
1146 
1147 
1148 
1149 
1150 
1151 
1152 
1153 
1154  //---------------------------------------------------------------------
1155  // Build a 3D mesh using hexes, tets, prisms, or pyramids.
1156  case 3:
1157  {
1158  libmesh_assert_not_equal_to (nx, 0);
1159  libmesh_assert_not_equal_to (ny, 0);
1160  libmesh_assert_not_equal_to (nz, 0);
1161  libmesh_assert_less (xmin, xmax);
1162  libmesh_assert_less (ymin, ymax);
1163  libmesh_assert_less (zmin, zmax);
1164 
1165 
1166  // Reserve elements. Meshes with prismatic elements require
1167  // twice as many elements.
1168  switch (type)
1169  {
1170  case INVALID_ELEM:
1171  case HEX8:
1172  case HEX20:
1173  case HEX27:
1174  case TET4: // TET4's are created from an initial HEX27 discretization
1175  case TET10: // TET10's are created from an initial HEX27 discretization
1176  case TET14: // TET14's are created from an initial HEX27 discretization
1177  case PYRAMID5: // PYRAMIDs are created from an initial HEX27 discretization
1178  case PYRAMID13:
1179  case PYRAMID14:
1180  case PYRAMID18:
1181  {
1182  mesh.reserve_elem(nx*ny*nz);
1183  break;
1184  }
1185 
1186  case PRISM6:
1187  case PRISM15:
1188  case PRISM18:
1189  case PRISM20:
1190  case PRISM21:
1191  {
1192  mesh.reserve_elem(2*nx*ny*nz);
1193  break;
1194  }
1195 
1196  default:
1197  libmesh_error_msg("ERROR: Unrecognized 3D element type == " << Utility::enum_to_string(type));
1198  }
1199 
1200 
1201 
1202 
1203 
1204  // Reserve nodes. Quadratic elements need twice as many nodes as linear elements.
1205  switch (type)
1206  {
1207  case INVALID_ELEM:
1208  case HEX8:
1209  case PRISM6:
1210  {
1211  mesh.reserve_nodes( (nx+1)*(ny+1)*(nz+1) );
1212  break;
1213  }
1214 
1215  case HEX20:
1216  case HEX27:
1217  case TET4: // TET4's are created from an initial HEX27 discretization
1218  case TET10: // TET10's are created from an initial HEX27 discretization
1219  case PYRAMID5: // PYRAMIDs are created from an initial HEX27 discretization
1220  case PYRAMID13:
1221  case PYRAMID14:
1222  case PYRAMID18:
1223  case PRISM15:
1224  case PRISM18:
1225  {
1226  // FYI: The resulting TET4 mesh will have exactly
1227  // 5*(nx*ny*nz) + 2*(nx*ny + nx*nz + ny*nz) + (nx+ny+nz) + 1
1228  // nodes once the additional mid-edge nodes for the HEX27 discretization
1229  // have been deleted.
1230  mesh.reserve_nodes( (2*nx+1)*(2*ny+1)*(2*nz+1) );
1231  break;
1232  }
1233 
1234  case TET14:
1235  {
1236  mesh.reserve_nodes( (2*nx+1)*(2*ny+1)*(2*nz+1) +
1237  24*nx*ny*nz +
1238  4*(nx*ny + ny*nz + nx*nz) );
1239  break;
1240  }
1241 
1242  case PRISM20:
1243  {
1244  mesh.reserve_nodes( (2*nx+1)*(2*ny+1)*(2*nz+1) +
1245  2*nx*ny*(nz+1) );
1246  break;
1247  }
1248 
1249  case PRISM21:
1250  {
1251  mesh.reserve_nodes( (2*nx+1)*(2*ny+1)*(2*nz+1) +
1252  2*nx*ny*(2*nz+1) );
1253  break;
1254  }
1255 
1256  default:
1257  libmesh_error_msg("ERROR: Unrecognized 3D element type == " << Utility::enum_to_string(type));
1258  }
1259 
1260 
1261 
1262 
1263  // Build the nodes.
1264  unsigned int node_id = 0;
1265  switch (type)
1266  {
1267  case INVALID_ELEM:
1268  case HEX8:
1269  case PRISM6:
1270  {
1271  for (unsigned int k=0; k<=nz; k++)
1272  for (unsigned int j=0; j<=ny; j++)
1273  for (unsigned int i=0; i<=nx; i++)
1274  {
1275  const Node * const node =
1276  mesh.add_point(Point(static_cast<Real>(i) / static_cast<Real>(nx),
1277  static_cast<Real>(j) / static_cast<Real>(ny),
1278  static_cast<Real>(k) / static_cast<Real>(nz)),
1279  node_id++);
1280  if (k == 0)
1281  boundary_info.add_node(node, 0);
1282  if (k == nz)
1283  boundary_info.add_node(node, 5);
1284  if (j == 0)
1285  boundary_info.add_node(node, 1);
1286  if (j == ny)
1287  boundary_info.add_node(node, 3);
1288  if (i == 0)
1289  boundary_info.add_node(node, 4);
1290  if (i == nx)
1291  boundary_info.add_node(node, 2);
1292  }
1293 
1294  break;
1295  }
1296 
1297  case HEX20:
1298  case HEX27:
1299  case TET4: // TET4's are created from an initial HEX27 discretization
1300  case TET10: // TET10's are created from an initial HEX27 discretization
1301  case TET14: // TET14's are created from an initial HEX27 discretization
1302  case PYRAMID5: // PYRAMIDs are created from an initial HEX27 discretization
1303  case PYRAMID13:
1304  case PYRAMID14:
1305  case PYRAMID18:
1306  case PRISM15:
1307  case PRISM18:
1308  case PRISM20:
1309  case PRISM21:
1310  {
1311  for (unsigned int k=0; k<=(2*nz); k++)
1312  for (unsigned int j=0; j<=(2*ny); j++)
1313  for (unsigned int i=0; i<=(2*nx); i++)
1314  {
1315  const Node * const node =
1316  mesh.add_point(Point(static_cast<Real>(i) / static_cast<Real>(2 * nx),
1317  static_cast<Real>(j) / static_cast<Real>(2 * ny),
1318  static_cast<Real>(k) / static_cast<Real>(2 * nz)),
1319  node_id++);
1320  if (k == 0)
1321  boundary_info.add_node(node, 0);
1322  if (k == 2*nz)
1323  boundary_info.add_node(node, 5);
1324  if (j == 0)
1325  boundary_info.add_node(node, 1);
1326  if (j == 2*ny)
1327  boundary_info.add_node(node, 3);
1328  if (i == 0)
1329  boundary_info.add_node(node, 4);
1330  if (i == 2*nx)
1331  boundary_info.add_node(node, 2);
1332  }
1333 
1334  if (type == PRISM20 ||
1335  type == PRISM21)
1336  {
1337  const unsigned int kmax = (type == PRISM20) ? nz : 2*nz;
1338  for (unsigned int k=0; k<=kmax; k++)
1339  for (unsigned int j=0; j<ny; j++)
1340  for (unsigned int i=0; i<nx; i++)
1341  {
1342  const Node * const node1 =
1343  mesh.add_point(Point((static_cast<Real>(i)+1/Real(3)) / static_cast<Real>(nx),
1344  (static_cast<Real>(j)+1/Real(3)) / static_cast<Real>(ny),
1345  static_cast<Real>(k) / static_cast<Real>(kmax)),
1346  node_id++);
1347  if (k == 0)
1348  boundary_info.add_node(node1, 0);
1349  if (k == kmax)
1350  boundary_info.add_node(node1, 5);
1351 
1352  const Node * const node2 =
1353  mesh.add_point(Point((static_cast<Real>(i)+2/Real(3)) / static_cast<Real>(nx),
1354  (static_cast<Real>(j)+2/Real(3)) / static_cast<Real>(ny),
1355  static_cast<Real>(k) / static_cast<Real>(kmax)),
1356  node_id++);
1357  if (k == 0)
1358  boundary_info.add_node(node2, 0);
1359  if (k == kmax)
1360  boundary_info.add_node(node2, 5);
1361  }
1362  }
1363 
1364  break;
1365  }
1366 
1367 
1368  default:
1369  libmesh_error_msg("ERROR: Unrecognized 3D element type == " << Utility::enum_to_string(type));
1370  }
1371 
1372 
1373 
1374 
1375  // Build the elements.
1376  unsigned int elem_id = 0;
1377  switch (type)
1378  {
1379  case INVALID_ELEM:
1380  case HEX8:
1381  {
1382  for (unsigned int k=0; k<nz; k++)
1383  for (unsigned int j=0; j<ny; j++)
1384  for (unsigned int i=0; i<nx; i++)
1385  {
1386  Elem * elem = mesh.add_elem(Elem::build_with_id(HEX8, elem_id++));
1387  elem->set_node(0, mesh.node_ptr(idx(type,nx,ny,i,j,k) ));
1388  elem->set_node(1, mesh.node_ptr(idx(type,nx,ny,i+1,j,k) ));
1389  elem->set_node(2, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) ));
1390  elem->set_node(3, mesh.node_ptr(idx(type,nx,ny,i,j+1,k) ));
1391  elem->set_node(4, mesh.node_ptr(idx(type,nx,ny,i,j,k+1) ));
1392  elem->set_node(5, mesh.node_ptr(idx(type,nx,ny,i+1,j,k+1) ));
1393  elem->set_node(6, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1)));
1394  elem->set_node(7, mesh.node_ptr(idx(type,nx,ny,i,j+1,k+1) ));
1395 
1396  if (k == 0)
1397  boundary_info.add_side(elem, 0, 0);
1398 
1399  if (k == (nz-1))
1400  boundary_info.add_side(elem, 5, 5);
1401 
1402  if (j == 0)
1403  boundary_info.add_side(elem, 1, 1);
1404 
1405  if (j == (ny-1))
1406  boundary_info.add_side(elem, 3, 3);
1407 
1408  if (i == 0)
1409  boundary_info.add_side(elem, 4, 4);
1410 
1411  if (i == (nx-1))
1412  boundary_info.add_side(elem, 2, 2);
1413  }
1414  break;
1415  }
1416 
1417 
1418 
1419 
1420  case PRISM6:
1421  {
1422  for (unsigned int k=0; k<nz; k++)
1423  for (unsigned int j=0; j<ny; j++)
1424  for (unsigned int i=0; i<nx; i++)
1425  {
1426  // First Prism
1427  Elem * elem = mesh.add_elem(Elem::build_with_id(PRISM6, elem_id++));
1428  elem->set_node(0, mesh.node_ptr(idx(type,nx,ny,i,j,k) ));
1429  elem->set_node(1, mesh.node_ptr(idx(type,nx,ny,i+1,j,k) ));
1430  elem->set_node(2, mesh.node_ptr(idx(type,nx,ny,i,j+1,k) ));
1431  elem->set_node(3, mesh.node_ptr(idx(type,nx,ny,i,j,k+1) ));
1432  elem->set_node(4, mesh.node_ptr(idx(type,nx,ny,i+1,j,k+1) ));
1433  elem->set_node(5, mesh.node_ptr(idx(type,nx,ny,i,j+1,k+1) ));
1434 
1435  // Add sides for first prism to boundary info object
1436  if (i==0)
1437  boundary_info.add_side(elem, 3, 4);
1438 
1439  if (j==0)
1440  boundary_info.add_side(elem, 1, 1);
1441 
1442  if (k==0)
1443  boundary_info.add_side(elem, 0, 0);
1444 
1445  if (k == (nz-1))
1446  boundary_info.add_side(elem, 4, 5);
1447 
1448  // Second Prism
1449  elem = mesh.add_elem(Elem::build_with_id(PRISM6, elem_id++));
1450  elem->set_node(0, mesh.node_ptr(idx(type,nx,ny,i+1,j,k) ));
1451  elem->set_node(1, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) ));
1452  elem->set_node(2, mesh.node_ptr(idx(type,nx,ny,i,j+1,k) ));
1453  elem->set_node(3, mesh.node_ptr(idx(type,nx,ny,i+1,j,k+1) ));
1454  elem->set_node(4, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1)));
1455  elem->set_node(5, mesh.node_ptr(idx(type,nx,ny,i,j+1,k+1) ));
1456 
1457  // Add sides for second prism to boundary info object
1458  if (i == (nx-1))
1459  boundary_info.add_side(elem, 1, 2);
1460 
1461  if (j == (ny-1))
1462  boundary_info.add_side(elem, 2, 3);
1463 
1464  if (k==0)
1465  boundary_info.add_side(elem, 0, 0);
1466 
1467  if (k == (nz-1))
1468  boundary_info.add_side(elem, 4, 5);
1469  }
1470  break;
1471  }
1472 
1473 
1474 
1475 
1476 
1477 
1478  case HEX20:
1479  case HEX27:
1480  case TET4: // TET4's are created from an initial HEX27 discretization
1481  case TET10: // TET10's are created from an initial HEX27 discretization
1482  case TET14: // TET14's are created from an initial HEX27 discretization
1483  case PYRAMID5: // PYRAMIDs are created from an initial HEX27 discretization
1484  case PYRAMID13:
1485  case PYRAMID14:
1486  case PYRAMID18:
1487  {
1488  for (unsigned int k=0; k<(2*nz); k += 2)
1489  for (unsigned int j=0; j<(2*ny); j += 2)
1490  for (unsigned int i=0; i<(2*nx); i += 2)
1491  {
1492  ElemType build_type = (type == HEX20) ? HEX20 : HEX27;
1493  Elem * elem = mesh.add_elem(Elem::build_with_id(build_type, elem_id++));
1494 
1495  elem->set_node(0, mesh.node_ptr(idx(type,nx,ny,i, j, k) ));
1496  elem->set_node(1, mesh.node_ptr(idx(type,nx,ny,i+2,j, k) ));
1497  elem->set_node(2, mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k) ));
1498  elem->set_node(3, mesh.node_ptr(idx(type,nx,ny,i, j+2,k) ));
1499  elem->set_node(4, mesh.node_ptr(idx(type,nx,ny,i, j, k+2)));
1500  elem->set_node(5, mesh.node_ptr(idx(type,nx,ny,i+2,j, k+2)));
1501  elem->set_node(6, mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k+2)));
1502  elem->set_node(7, mesh.node_ptr(idx(type,nx,ny,i, j+2,k+2)));
1503  elem->set_node(8, mesh.node_ptr(idx(type,nx,ny,i+1,j, k) ));
1504  elem->set_node(9, mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k) ));
1505  elem->set_node(10, mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k) ));
1506  elem->set_node(11, mesh.node_ptr(idx(type,nx,ny,i, j+1,k) ));
1507  elem->set_node(12, mesh.node_ptr(idx(type,nx,ny,i, j, k+1)));
1508  elem->set_node(13, mesh.node_ptr(idx(type,nx,ny,i+2,j, k+1)));
1509  elem->set_node(14, mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k+1)));
1510  elem->set_node(15, mesh.node_ptr(idx(type,nx,ny,i, j+2,k+1)));
1511  elem->set_node(16, mesh.node_ptr(idx(type,nx,ny,i+1,j, k+2)));
1512  elem->set_node(17, mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k+2)));
1513  elem->set_node(18, mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k+2)));
1514  elem->set_node(19, mesh.node_ptr(idx(type,nx,ny,i, j+1,k+2)));
1515 
1516  if ((type == HEX27) || (type == TET4) || (type == TET10) || (type == TET14) ||
1517  (type == PYRAMID5) || (type == PYRAMID13) || (type == PYRAMID14) ||
1518  (type == PYRAMID18))
1519  {
1520  elem->set_node(20, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) ));
1521  elem->set_node(21, mesh.node_ptr(idx(type,nx,ny,i+1,j, k+1)));
1522  elem->set_node(22, mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k+1)));
1523  elem->set_node(23, mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k+1)));
1524  elem->set_node(24, mesh.node_ptr(idx(type,nx,ny,i, j+1,k+1)));
1525  elem->set_node(25, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+2)));
1526  elem->set_node(26, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1)));
1527  }
1528 
1529  if (k == 0)
1530  boundary_info.add_side(elem, 0, 0);
1531 
1532  if (k == 2*(nz-1))
1533  boundary_info.add_side(elem, 5, 5);
1534 
1535  if (j == 0)
1536  boundary_info.add_side(elem, 1, 1);
1537 
1538  if (j == 2*(ny-1))
1539  boundary_info.add_side(elem, 3, 3);
1540 
1541  if (i == 0)
1542  boundary_info.add_side(elem, 4, 4);
1543 
1544  if (i == 2*(nx-1))
1545  boundary_info.add_side(elem, 2, 2);
1546  }
1547  break;
1548  }
1549 
1550 
1551 
1552 
1553  case PRISM15:
1554  case PRISM18:
1555  case PRISM20:
1556  case PRISM21:
1557  {
1558  for (unsigned int k=0; k<(2*nz); k += 2)
1559  for (unsigned int j=0; j<(2*ny); j += 2)
1560  for (unsigned int i=0; i<(2*nx); i += 2)
1561  {
1562  // First Prism
1563  Elem * elem = mesh.add_elem(Elem::build_with_id(type, elem_id++));
1564  elem->set_node(0, mesh.node_ptr(idx(type,nx,ny,i, j, k) ));
1565  elem->set_node(1, mesh.node_ptr(idx(type,nx,ny,i+2,j, k) ));
1566  elem->set_node(2, mesh.node_ptr(idx(type,nx,ny,i, j+2,k) ));
1567  elem->set_node(3, mesh.node_ptr(idx(type,nx,ny,i, j, k+2)));
1568  elem->set_node(4, mesh.node_ptr(idx(type,nx,ny,i+2,j, k+2)));
1569  elem->set_node(5, mesh.node_ptr(idx(type,nx,ny,i, j+2,k+2)));
1570  elem->set_node(6, mesh.node_ptr(idx(type,nx,ny,i+1,j, k) ));
1571  elem->set_node(7, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) ));
1572  elem->set_node(8, mesh.node_ptr(idx(type,nx,ny,i, j+1,k) ));
1573  elem->set_node(9, mesh.node_ptr(idx(type,nx,ny,i, j, k+1)));
1574  elem->set_node(10, mesh.node_ptr(idx(type,nx,ny,i+2,j, k+1)));
1575  elem->set_node(11, mesh.node_ptr(idx(type,nx,ny,i, j+2,k+1)));
1576  elem->set_node(12, mesh.node_ptr(idx(type,nx,ny,i+1,j, k+2)));
1577  elem->set_node(13, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+2)));
1578  elem->set_node(14, mesh.node_ptr(idx(type,nx,ny,i, j+1,k+2)));
1579 
1580  if (type == PRISM18 ||
1581  type == PRISM20 ||
1582  type == PRISM21)
1583  {
1584  elem->set_node(15, mesh.node_ptr(idx(type,nx,ny,i+1,j, k+1)));
1585  elem->set_node(16, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1)));
1586  elem->set_node(17, mesh.node_ptr(idx(type,nx,ny,i, j+1,k+1)));
1587  }
1588 
1589  if (type == PRISM20)
1590  {
1591  const dof_id_type base_idx = (2*nx+1)*(2*ny+1)*(2*nz+1);
1592  elem->set_node(18, mesh.node_ptr(base_idx+((k/2)*(nx*ny)+j/2*nx+i/2)*2));
1593  elem->set_node(19, mesh.node_ptr(base_idx+(((k/2)+1)*(nx*ny)+j/2*nx+i/2)*2));
1594  }
1595 
1596  if (type == PRISM21)
1597  {
1598  const dof_id_type base_idx = (2*nx+1)*(2*ny+1)*(2*nz+1);
1599  elem->set_node(18, mesh.node_ptr(base_idx+(k*(nx*ny)+j/2*nx+i/2)*2));
1600  elem->set_node(19, mesh.node_ptr(base_idx+((k+2)*(nx*ny)+j/2*nx+i/2)*2));
1601  elem->set_node(20, mesh.node_ptr(base_idx+((k+1)*(nx*ny)+j/2*nx+i/2)*2));
1602  }
1603 
1604  // Add sides for first prism to boundary info object
1605  if (i==0)
1606  boundary_info.add_side(elem, 3, 4);
1607 
1608  if (j==0)
1609  boundary_info.add_side(elem, 1, 1);
1610 
1611  if (k==0)
1612  boundary_info.add_side(elem, 0, 0);
1613 
1614  if (k == 2*(nz-1))
1615  boundary_info.add_side(elem, 4, 5);
1616 
1617 
1618  // Second Prism
1619  elem = mesh.add_elem(Elem::build_with_id(type, elem_id++));
1620  elem->set_node(0, mesh.node_ptr(idx(type,nx,ny,i+2,j,k) ));
1621  elem->set_node(1, mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k) ));
1622  elem->set_node(2, mesh.node_ptr(idx(type,nx,ny,i,j+2,k) ));
1623  elem->set_node(3, mesh.node_ptr(idx(type,nx,ny,i+2,j,k+2) ));
1624  elem->set_node(4, mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k+2) ));
1625  elem->set_node(5, mesh.node_ptr(idx(type,nx,ny,i,j+2,k+2) ));
1626  elem->set_node(6, mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k) ));
1627  elem->set_node(7, mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k) ));
1628  elem->set_node(8, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) ));
1629  elem->set_node(9, mesh.node_ptr(idx(type,nx,ny,i+2,j,k+1) ));
1630  elem->set_node(10, mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k+1)));
1631  elem->set_node(11, mesh.node_ptr(idx(type,nx,ny,i,j+2,k+1) ));
1632  elem->set_node(12, mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k+2)));
1633  elem->set_node(13, mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k+2)));
1634  elem->set_node(14, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+2)));
1635 
1636  if (type == PRISM18 ||
1637  type == PRISM20 ||
1638  type == PRISM21)
1639  {
1640  elem->set_node(15, mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k+1)));
1641  elem->set_node(16, mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k+1)));
1642  elem->set_node(17, mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1)));
1643  }
1644 
1645  if (type == PRISM20)
1646  {
1647  const dof_id_type base_idx = (2*nx+1)*(2*ny+1)*(2*nz+1);
1648  elem->set_node(18, mesh.node_ptr(base_idx+((k/2)*(nx*ny)+j/2*nx+i/2)*2+1));
1649  elem->set_node(19, mesh.node_ptr(base_idx+(((k/2)+1)*(nx*ny)+j/2*nx+i/2)*2+1));
1650  }
1651 
1652  if (type == PRISM21)
1653  {
1654  const dof_id_type base_idx = (2*nx+1)*(2*ny+1)*(2*nz+1);
1655  elem->set_node(18, mesh.node_ptr(base_idx+(k*(nx*ny)+j/2*nx+i/2)*2+1));
1656  elem->set_node(19, mesh.node_ptr(base_idx+((k+2)*(nx*ny)+j/2*nx+i/2)*2+1));
1657  elem->set_node(20, mesh.node_ptr(base_idx+((k+1)*(nx*ny)+j/2*nx+i/2)*2+1));
1658  }
1659 
1660  // Add sides for second prism to boundary info object
1661  if (i == 2*(nx-1))
1662  boundary_info.add_side(elem, 1, 2);
1663 
1664  if (j == 2*(ny-1))
1665  boundary_info.add_side(elem, 2, 3);
1666 
1667  if (k==0)
1668  boundary_info.add_side(elem, 0, 0);
1669 
1670  if (k == 2*(nz-1))
1671  boundary_info.add_side(elem, 4, 5);
1672 
1673  }
1674  break;
1675  }
1676 
1677 
1678 
1679 
1680 
1681  default:
1682  libmesh_error_msg("ERROR: Unrecognized 3D element type == " << Utility::enum_to_string(type));
1683  }
1684 
1685 
1686 
1687 
1688  //.......................................
1689  // Scale the nodal positions
1690  if (gauss_lobatto_grid)
1691  {
1692  GaussLobattoRedistributionFunction func(nx, xmin, xmax,
1693  ny, ymin, ymax,
1694  nz, zmin, zmax);
1695  MeshTools::Modification::redistribute(mesh, func);
1696  }
1697  else // !gauss_lobatto_grid
1698  {
1699  for (unsigned int p=0; p<mesh.n_nodes(); p++)
1700  {
1701  mesh.node_ref(p)(0) = (mesh.node_ref(p)(0))*(xmax-xmin) + xmin;
1702  mesh.node_ref(p)(1) = (mesh.node_ref(p)(1))*(ymax-ymin) + ymin;
1703  mesh.node_ref(p)(2) = (mesh.node_ref(p)(2))*(zmax-zmin) + zmin;
1704  }
1705  }
1706 
1707 
1708 
1709  // Additional work for tets and pyramids: we take the existing
1710  // HEX27 discretization and split each element into 24
1711  // sub-tets or 6 sub-pyramids.
1712  //
1713  // 24 isn't the minimum-possible number of tets, but it
1714  // obviates any concerns about the edge orientations between
1715  // the various elements.
1716  if ((type == TET4) ||
1717  (type == TET10) ||
1718  (type == TET14) ||
1719  (type == PYRAMID5) ||
1720  (type == PYRAMID13) ||
1721  (type == PYRAMID14) ||
1722  (type == PYRAMID18))
1723  {
1724  // Temporary storage for new elements. (24 tets per hex, 6 pyramids)
1725  std::vector<std::unique_ptr<Elem>> new_elements;
1726 
1727  // For avoiding extraneous construction of element sides
1728  std::unique_ptr<Elem> side;
1729 
1730  if ((type == TET4) || (type == TET10) || (type == TET14))
1731  new_elements.reserve(24*mesh.n_elem());
1732  else
1733  new_elements.reserve(6*mesh.n_elem());
1734 
1735  // Create tetrahedra or pyramids
1736  for (auto & base_hex : mesh.element_ptr_range())
1737  {
1738  // Get a pointer to the node located at the HEX27 center
1739  Node * apex_node = base_hex->node_ptr(26);
1740 
1741  // Container to catch ids handed back from BoundaryInfo
1742  std::vector<boundary_id_type> ids;
1743 
1744  for (auto s : base_hex->side_index_range())
1745  {
1746  // Get the boundary ID(s) for this side
1747  boundary_info.boundary_ids(base_hex, s, ids);
1748 
1749  // We're creating this Mesh, so there should be 0 or 1 boundary IDs.
1750  libmesh_assert(ids.size() <= 1);
1751 
1752  // A convenient name for the side's ID.
1753  boundary_id_type b_id = ids.empty() ? BoundaryInfo::invalid_id : ids[0];
1754 
1755  // Need to build the full-ordered side!
1756  base_hex->build_side_ptr(side, s);
1757 
1758  if ((type == TET4) || (type == TET10) || (type == TET14))
1759  {
1760  // Build 4 sub-tets per side
1761  for (unsigned int sub_tet=0; sub_tet<4; ++sub_tet)
1762  {
1763  new_elements.push_back( Elem::build(TET4) );
1764  auto & sub_elem = new_elements.back();
1765  sub_elem->set_node(0, side->node_ptr(sub_tet));
1766  sub_elem->set_node(1, side->node_ptr(8)); // center of the face
1767  sub_elem->set_node(2, side->node_ptr(sub_tet==3 ? 0 : sub_tet+1 )); // wrap-around
1768  sub_elem->set_node(3, apex_node); // apex node always used!
1769 
1770  // If the original hex was a boundary hex, add the new sub_tet's side
1771  // 0 with the same b_id. Note: the tets are all aligned so that their
1772  // side 0 is on the boundary.
1773  if (b_id != BoundaryInfo::invalid_id)
1774  boundary_info.add_side(sub_elem.get(), 0, b_id);
1775  }
1776  } // end if ((type == TET4) || (type == TET10) || (type == TET14))
1777 
1778  else // type==PYRAMID*
1779  {
1780  // Build 1 sub-pyramid per side.
1781  new_elements.push_back( Elem::build(PYRAMID5) );
1782  auto & sub_elem = new_elements.back();
1783 
1784  // Set the base. Note that since the apex is *inside* the base_hex,
1785  // and the pyramid uses a counter-clockwise base numbering, we need to
1786  // reverse the [1] and [3] node indices.
1787  sub_elem->set_node(0, side->node_ptr(0));
1788  sub_elem->set_node(1, side->node_ptr(3));
1789  sub_elem->set_node(2, side->node_ptr(2));
1790  sub_elem->set_node(3, side->node_ptr(1));
1791 
1792  // Set the apex
1793  sub_elem->set_node(4, apex_node);
1794 
1795  // If the original hex was a boundary hex, add the new sub_pyr's side
1796  // 4 (the square base) with the same b_id.
1797  if (b_id != BoundaryInfo::invalid_id)
1798  boundary_info.add_side(sub_elem.get(), 4, b_id);
1799  } // end else type==PYRAMID*
1800  }
1801  }
1802 
1803 
1804  // Delete the original HEX27 elements from the mesh, and the boundary info structure.
1805  for (auto & elem : mesh.element_ptr_range())
1806  {
1807  boundary_info.remove(elem); // Safe even if elem has no boundary info.
1808  mesh.delete_elem(elem);
1809  }
1810 
1811  // Add the new elements
1812  for (auto i : index_range(new_elements))
1813  {
1814  new_elements[i]->set_id(i);
1815  mesh.add_elem( std::move(new_elements[i]) );
1816  }
1817 
1818  } // end if (type == TET*,PYRAMID*)
1819 
1820 
1821  // Use all_second_order to convert the TET4's to TET10's or PYRAMID5's to PYRAMID14's
1822  if ((type == TET10) || (type == PYRAMID14))
1823  mesh.all_second_order();
1824 
1825  else if (type == PYRAMID13)
1826  mesh.all_second_order(/*full_ordered=*/false);
1827 
1828  else if ((type == TET14) || (type == PYRAMID18))
1829  mesh.all_complete_order();
1830 
1831 
1832  // Add sideset names to boundary info (Z axis out of the screen)
1833  boundary_info.sideset_name(0) = "back";
1834  boundary_info.sideset_name(1) = "bottom";
1835  boundary_info.sideset_name(2) = "right";
1836  boundary_info.sideset_name(3) = "top";
1837  boundary_info.sideset_name(4) = "left";
1838  boundary_info.sideset_name(5) = "front";
1839 
1840  // Add nodeset names to boundary info
1841  boundary_info.nodeset_name(0) = "back";
1842  boundary_info.nodeset_name(1) = "bottom";
1843  boundary_info.nodeset_name(2) = "right";
1844  boundary_info.nodeset_name(3) = "top";
1845  boundary_info.nodeset_name(4) = "left";
1846  boundary_info.nodeset_name(5) = "front";
1847 
1848  break;
1849  } // end case dim==3
1850 
1851  default:
1852  libmesh_error_msg("Unknown dimension " << mesh.mesh_dimension());
1853  }
1854 
1855  // Done building the mesh. Now prepare it for use.
1856  mesh.prepare_for_use ();
1857 }
1858 
1859 
1860 
1861 void MeshTools::Generation::build_point (UnstructuredMesh & mesh,
1862  const ElemType type,
1863  const bool gauss_lobatto_grid)
1864 {
1865  // This method only makes sense in 0D!
1866  // But we now just turn a non-0D mesh into a 0D mesh
1867  //libmesh_assert_equal_to (mesh.mesh_dimension(), 1);
1868 
1869  build_cube(mesh,
1870  0, 0, 0,
1871  0., 0.,
1872  0., 0.,
1873  0., 0.,
1874  type,
1875  gauss_lobatto_grid);
1876 }
1877 
1878 
1879 void MeshTools::Generation::build_line (UnstructuredMesh & mesh,
1880  const unsigned int nx,
1881  const Real xmin, const Real xmax,
1882  const ElemType type,
1883  const bool gauss_lobatto_grid)
1884 {
1885  // This method only makes sense in 1D!
1886  // But we now just turn a non-1D mesh into a 1D mesh
1887  //libmesh_assert_equal_to (mesh.mesh_dimension(), 1);
1888 
1889  build_cube(mesh,
1890  nx, 0, 0,
1891  xmin, xmax,
1892  0., 0.,
1893  0., 0.,
1894  type,
1895  gauss_lobatto_grid);
1896 }
1897 
1898 
1899 
1900 void MeshTools::Generation::build_square (UnstructuredMesh & mesh,
1901  const unsigned int nx,
1902  const unsigned int ny,
1903  const Real xmin, const Real xmax,
1904  const Real ymin, const Real ymax,
1905  const ElemType type,
1906  const bool gauss_lobatto_grid)
1907 {
1908  // This method only makes sense in 2D!
1909  // But we now just turn a non-2D mesh into a 2D mesh
1910  //libmesh_assert_equal_to (mesh.mesh_dimension(), 2);
1911 
1912  // Call the build_cube() member to actually do the work for us.
1913  build_cube (mesh,
1914  nx, ny, 0,
1915  xmin, xmax,
1916  ymin, ymax,
1917  0., 0.,
1918  type,
1919  gauss_lobatto_grid);
1920 }
1921 
1922 
1923 
1924 
1925 
1926 
1927 
1928 
1929 
1930 #ifndef LIBMESH_ENABLE_AMR
1931 void MeshTools::Generation::build_sphere (UnstructuredMesh &,
1932  const Real,
1933  const unsigned int,
1934  const ElemType,
1935  const unsigned int,
1936  const bool)
1937 {
1938  libmesh_error_msg("Building a circle/sphere only works with AMR.");
1939 }
1940 
1941 #else
1942 
1943 void MeshTools::Generation::build_sphere (UnstructuredMesh & mesh,
1944  const Real rad,
1945  const unsigned int nr,
1946  const ElemType type,
1947  const unsigned int n_smooth,
1948  const bool flat)
1949 {
1950  libmesh_assert_greater (rad, 0.);
1951  //libmesh_assert_greater (nr, 0); // must refine at least once otherwise will end up with a square/cube
1952 
1953  LOG_SCOPE("build_sphere()", "MeshTools::Generation");
1954 
1955  // Clear the mesh and start from scratch, but save the original
1956  // mesh_dimension, since the original intent of this function was to
1957  // allow the geometric entity (line, circle, ball, sphere)
1958  // constructed to be determined by the mesh's dimension.
1959  unsigned char orig_mesh_dimension =
1960  cast_int<unsigned char>(mesh.mesh_dimension());
1961  mesh.clear();
1962  mesh.set_mesh_dimension(orig_mesh_dimension);
1963 
1964  // If mesh.mesh_dimension()==1, it *could* be because the user
1965  // constructed a Mesh without specifying a dimension (since this is
1966  // allowed now) and hence it got the default dimension of 1. In
1967  // this case, we will try to infer the dimension they *really*
1968  // wanted from the requested ElemType, and if they don't match, go
1969  // with the ElemType.
1970  if (mesh.mesh_dimension() == 1)
1971  {
1972  switch (type)
1973  {
1974  case HEX8:
1975  case HEX27:
1976  case TET4:
1977  case TET10:
1978  case TET14:
1979  mesh.set_mesh_dimension(3);
1980  break;
1981  case TRI3:
1982  case TRI6:
1983  case TRI7:
1984  case QUAD4:
1985  case QUADSHELL4:
1986  case QUAD8:
1987  case QUADSHELL8:
1988  case QUAD9:
1989  case QUADSHELL9:
1990  mesh.set_mesh_dimension(2);
1991  break;
1992  case EDGE2:
1993  case EDGE3:
1994  case EDGE4:
1995  mesh.set_mesh_dimension(1);
1996  break;
1997  case INVALID_ELEM:
1998  // Just keep the existing dimension
1999  break;
2000  default:
2001  libmesh_error_msg("build_sphere(): Please specify a mesh dimension or a valid ElemType (EDGE{2,3,4}, TRI{3,6,7}, QUAD{4,8,9}, HEX{8,27}, TET{4,10,14})");
2002  }
2003  }
2004 
2005  BoundaryInfo & boundary_info = mesh.get_boundary_info();
2006 
2007  // Building while distributed is a little more complicated
2008  const bool is_replicated = mesh.is_replicated();
2009 
2010  // Sphere is centered at origin by default
2011  const Point cent;
2012 
2013  const Sphere sphere (cent, rad);
2014 
2015  switch (mesh.mesh_dimension())
2016  {
2017  //-----------------------------------------------------------------
2018  // Build a line in one dimension
2019  case 1:
2020  {
2021  build_line (mesh, 3, -rad, rad, type);
2022 
2023  break;
2024  }
2025 
2026 
2027 
2028 
2029  //-----------------------------------------------------------------
2030  // Build a circle or hollow sphere in two dimensions
2031  case 2:
2032  {
2033  // For DistributedMesh, if we don't specify node IDs the Mesh
2034  // will try to pick an appropriate (unique) one for us. But
2035  // since we are adding these nodes on all processors, we want
2036  // to be sure they have consistent IDs across all processors.
2037  unsigned node_id = 0;
2038 
2039  if (flat)
2040  {
2041  const Real sqrt_2 = std::sqrt(2.);
2042  const Real rad_2 = .25*rad;
2043  const Real rad_sqrt_2 = rad/sqrt_2;
2044 
2045  // (Temporary) convenient storage for node pointers
2046  std::vector<Node *> nodes(8);
2047 
2048  // Point 0
2049  nodes[0] = mesh.add_point (Point(-rad_2,-rad_2, 0.), node_id++);
2050 
2051  // Point 1
2052  nodes[1] = mesh.add_point (Point( rad_2,-rad_2, 0.), node_id++);
2053 
2054  // Point 2
2055  nodes[2] = mesh.add_point (Point( rad_2, rad_2, 0.), node_id++);
2056 
2057  // Point 3
2058  nodes[3] = mesh.add_point (Point(-rad_2, rad_2, 0.), node_id++);
2059 
2060  // Point 4
2061  nodes[4] = mesh.add_point (Point(-rad_sqrt_2,-rad_sqrt_2, 0.), node_id++);
2062 
2063  // Point 5
2064  nodes[5] = mesh.add_point (Point( rad_sqrt_2,-rad_sqrt_2, 0.), node_id++);
2065 
2066  // Point 6
2067  nodes[6] = mesh.add_point (Point( rad_sqrt_2, rad_sqrt_2, 0.), node_id++);
2068 
2069  // Point 7
2070  nodes[7] = mesh.add_point (Point(-rad_sqrt_2, rad_sqrt_2, 0.), node_id++);
2071 
2072  // Build the elements & set node pointers
2073 
2074  // Element 0
2075  {
2076  Elem * elem0 = mesh.add_elem (Elem::build(QUAD4));
2077  elem0->set_node(0, nodes[0]);
2078  elem0->set_node(1, nodes[1]);
2079  elem0->set_node(2, nodes[2]);
2080  elem0->set_node(3, nodes[3]);
2081  }
2082 
2083  // Element 1
2084  {
2085  Elem * elem1 = mesh.add_elem (Elem::build(QUAD4));
2086  elem1->set_node(0, nodes[4]);
2087  elem1->set_node(1, nodes[0]);
2088  elem1->set_node(2, nodes[3]);
2089  elem1->set_node(3, nodes[7]);
2090  }
2091 
2092  // Element 2
2093  {
2094  Elem * elem2 = mesh.add_elem (Elem::build(QUAD4));
2095  elem2->set_node(0, nodes[4]);
2096  elem2->set_node(1, nodes[5]);
2097  elem2->set_node(2, nodes[1]);
2098  elem2->set_node(3, nodes[0]);
2099  }
2100 
2101  // Element 3
2102  {
2103  Elem * elem3 = mesh.add_elem (Elem::build(QUAD4));
2104  elem3->set_node(0, nodes[1]);
2105  elem3->set_node(1, nodes[5]);
2106  elem3->set_node(2, nodes[6]);
2107  elem3->set_node(3, nodes[2]);
2108  }
2109 
2110  // Element 4
2111  {
2112  Elem * elem4 = mesh.add_elem (Elem::build(QUAD4));
2113  elem4->set_node(0, nodes[3]);
2114  elem4->set_node(1, nodes[2]);
2115  elem4->set_node(2, nodes[6]);
2116  elem4->set_node(3, nodes[7]);
2117  }
2118 
2119  }
2120  else
2121  {
2122  // Create the 12 vertices of a regular unit icosahedron
2123  Real t = 0.5 * (1 + std::sqrt(5.0));
2124  Real s = rad / std::sqrt(1 + t*t);
2125  t *= s;
2126 
2127  mesh.add_point (Point(-s, t, 0), node_id++);
2128  mesh.add_point (Point( s, t, 0), node_id++);
2129  mesh.add_point (Point(-s, -t, 0), node_id++);
2130  mesh.add_point (Point( s, -t, 0), node_id++);
2131 
2132  mesh.add_point (Point( 0, -s, t), node_id++);
2133  mesh.add_point (Point( 0, s, t), node_id++);
2134  mesh.add_point (Point( 0, -s, -t), node_id++);
2135  mesh.add_point (Point( 0, s, -t), node_id++);
2136 
2137  mesh.add_point (Point( t, 0, -s), node_id++);
2138  mesh.add_point (Point( t, 0, s), node_id++);
2139  mesh.add_point (Point(-t, 0, -s), node_id++);
2140  mesh.add_point (Point(-t, 0, s), node_id++);
2141 
2142  // Create the 20 triangles of the icosahedron
2143  static const unsigned int idx1 [6] = {11, 5, 1, 7, 10, 11};
2144  static const unsigned int idx2 [6] = {9, 4, 2, 6, 8, 9};
2145  static const unsigned int idx3 [6] = {1, 5, 11, 10, 7, 1};
2146 
2147  for (unsigned int i = 0; i < 5; ++i)
2148  {
2149  // 5 elems around point 0
2150  Elem * new_elem = mesh.add_elem(Elem::build(TRI3));
2151  new_elem->set_node(0, mesh.node_ptr(0));
2152  new_elem->set_node(1, mesh.node_ptr(idx1[i]));
2153  new_elem->set_node(2, mesh.node_ptr(idx1[i+1]));
2154 
2155  // 5 adjacent elems
2156  new_elem = mesh.add_elem(Elem::build(TRI3));
2157  new_elem->set_node(0, mesh.node_ptr(idx3[i]));
2158  new_elem->set_node(1, mesh.node_ptr(idx3[i+1]));
2159  new_elem->set_node(2, mesh.node_ptr(idx2[i]));
2160 
2161  // 5 elems around point 3
2162  new_elem = mesh.add_elem(Elem::build(TRI3));
2163  new_elem->set_node(0, mesh.node_ptr(3));
2164  new_elem->set_node(1, mesh.node_ptr(idx2[i]));
2165  new_elem->set_node(2, mesh.node_ptr(idx2[i+1]));
2166 
2167  // 5 adjacent elems
2168  new_elem = mesh.add_elem(Elem::build(TRI3));
2169  new_elem->set_node(0, mesh.node_ptr(idx2[i+1]));
2170  new_elem->set_node(1, mesh.node_ptr(idx2[i]));
2171  new_elem->set_node(2, mesh.node_ptr(idx3[i+1]));
2172  }
2173  }
2174 
2175  break;
2176  } // end case 2
2177 
2178 
2179 
2180 
2181 
2182  //-----------------------------------------------------------------
2183  // Build a sphere in three dimensions
2184  case 3:
2185  {
2186  // (Currently) supported types
2187  if (!((type == HEX8) || (type == HEX27) || (type == TET4) ||
2188  (type == TET10) || (type == TET14)))
2189  {
2190  libmesh_error_msg("Error: Only HEX8/27 and TET4/10/14 are currently supported in 3D.");
2191  }
2192 
2193 
2194  // 3D analog of 2D initial grid:
2195  const Real
2196  r_small = 0.25*rad, // 0.25 *radius
2197  r_med = (0.125*std::sqrt(2.)+0.5)*rad; // .67677*radius
2198 
2199  // (Temporary) convenient storage for node pointers
2200  std::vector<Node *> nodes(16);
2201 
2202  // For DistributedMesh, if we don't specify node IDs the Mesh
2203  // will try to pick an appropriate (unique) one for us. But
2204  // since we are adding these nodes on all processors, we want
2205  // to be sure they have consistent IDs across all processors.
2206  unsigned node_id = 0;
2207 
2208  // Points 0-7 are the initial HEX8
2209  nodes[0] = mesh.add_point (Point(-r_small,-r_small, -r_small), node_id++);
2210  nodes[1] = mesh.add_point (Point( r_small,-r_small, -r_small), node_id++);
2211  nodes[2] = mesh.add_point (Point( r_small, r_small, -r_small), node_id++);
2212  nodes[3] = mesh.add_point (Point(-r_small, r_small, -r_small), node_id++);
2213  nodes[4] = mesh.add_point (Point(-r_small,-r_small, r_small), node_id++);
2214  nodes[5] = mesh.add_point (Point( r_small,-r_small, r_small), node_id++);
2215  nodes[6] = mesh.add_point (Point( r_small, r_small, r_small), node_id++);
2216  nodes[7] = mesh.add_point (Point(-r_small, r_small, r_small), node_id++);
2217 
2218  // Points 8-15 are for the outer hexes, we number them in the same way
2219  nodes[8] = mesh.add_point (Point(-r_med,-r_med, -r_med), node_id++);
2220  nodes[9] = mesh.add_point (Point( r_med,-r_med, -r_med), node_id++);
2221  nodes[10] = mesh.add_point (Point( r_med, r_med, -r_med), node_id++);
2222  nodes[11] = mesh.add_point (Point(-r_med, r_med, -r_med), node_id++);
2223  nodes[12] = mesh.add_point (Point(-r_med,-r_med, r_med), node_id++);
2224  nodes[13] = mesh.add_point (Point( r_med,-r_med, r_med), node_id++);
2225  nodes[14] = mesh.add_point (Point( r_med, r_med, r_med), node_id++);
2226  nodes[15] = mesh.add_point (Point(-r_med, r_med, r_med), node_id++);
2227 
2228  // Now create the elements and add them to the mesh
2229  // Element 0 - center element
2230  {
2231  Elem * elem0 = mesh.add_elem(Elem::build(HEX8));
2232  elem0->set_node(0, nodes[0]);
2233  elem0->set_node(1, nodes[1]);
2234  elem0->set_node(2, nodes[2]);
2235  elem0->set_node(3, nodes[3]);
2236  elem0->set_node(4, nodes[4]);
2237  elem0->set_node(5, nodes[5]);
2238  elem0->set_node(6, nodes[6]);
2239  elem0->set_node(7, nodes[7]);
2240  }
2241 
2242  // Element 1 - "bottom"
2243  {
2244  Elem * elem1 = mesh.add_elem(Elem::build(HEX8));
2245  elem1->set_node(0, nodes[8]);
2246  elem1->set_node(1, nodes[9]);
2247  elem1->set_node(2, nodes[10]);
2248  elem1->set_node(3, nodes[11]);
2249  elem1->set_node(4, nodes[0]);
2250  elem1->set_node(5, nodes[1]);
2251  elem1->set_node(6, nodes[2]);
2252  elem1->set_node(7, nodes[3]);
2253  }
2254 
2255  // Element 2 - "front"
2256  {
2257  Elem * elem2 = mesh.add_elem(Elem::build(HEX8));
2258  elem2->set_node(0, nodes[8]);
2259  elem2->set_node(1, nodes[9]);
2260  elem2->set_node(2, nodes[1]);
2261  elem2->set_node(3, nodes[0]);
2262  elem2->set_node(4, nodes[12]);
2263  elem2->set_node(5, nodes[13]);
2264  elem2->set_node(6, nodes[5]);
2265  elem2->set_node(7, nodes[4]);
2266  }
2267 
2268  // Element 3 - "right"
2269  {
2270  Elem * elem3 = mesh.add_elem(Elem::build(HEX8));
2271  elem3->set_node(0, nodes[1]);
2272  elem3->set_node(1, nodes[9]);
2273  elem3->set_node(2, nodes[10]);
2274  elem3->set_node(3, nodes[2]);
2275  elem3->set_node(4, nodes[5]);
2276  elem3->set_node(5, nodes[13]);
2277  elem3->set_node(6, nodes[14]);
2278  elem3->set_node(7, nodes[6]);
2279  }
2280 
2281  // Element 4 - "back"
2282  {
2283  Elem * elem4 = mesh.add_elem(Elem::build(HEX8));
2284  elem4->set_node(0, nodes[3]);
2285  elem4->set_node(1, nodes[2]);
2286  elem4->set_node(2, nodes[10]);
2287  elem4->set_node(3, nodes[11]);
2288  elem4->set_node(4, nodes[7]);
2289  elem4->set_node(5, nodes[6]);
2290  elem4->set_node(6, nodes[14]);
2291  elem4->set_node(7, nodes[15]);
2292  }
2293 
2294  // Element 5 - "left"
2295  {
2296  Elem * elem5 = mesh.add_elem(Elem::build(HEX8));
2297  elem5->set_node(0, nodes[8]);
2298  elem5->set_node(1, nodes[0]);
2299  elem5->set_node(2, nodes[3]);
2300  elem5->set_node(3, nodes[11]);
2301  elem5->set_node(4, nodes[12]);
2302  elem5->set_node(5, nodes[4]);
2303  elem5->set_node(6, nodes[7]);
2304  elem5->set_node(7, nodes[15]);
2305  }
2306 
2307  // Element 6 - "top"
2308  {
2309  Elem * elem6 = mesh.add_elem(Elem::build(HEX8));
2310  elem6->set_node(0, nodes[4]);
2311  elem6->set_node(1, nodes[5]);
2312  elem6->set_node(2, nodes[6]);
2313  elem6->set_node(3, nodes[7]);
2314  elem6->set_node(4, nodes[12]);
2315  elem6->set_node(5, nodes[13]);
2316  elem6->set_node(6, nodes[14]);
2317  elem6->set_node(7, nodes[15]);
2318  }
2319 
2320  break;
2321  } // end case 3
2322 
2323  default:
2324  libmesh_error_msg("Unknown dimension " << mesh.mesh_dimension());
2325 
2326 
2327 
2328  } // end switch (dim)
2329 
2330  // Now we have the beginnings of a sphere.
2331  // Add some more elements by doing uniform refinements and
2332  // popping nodes to the boundary.
2333  MeshRefinement mesh_refinement (mesh);
2334 
2335  // For avoiding extraneous element side construction
2336  std::unique_ptr<Elem> side;
2337 
2338  // Loop over the elements, refine, pop nodes to boundary.
2339  for (unsigned int r=0; r<nr; r++)
2340  {
2341  // A DistributedMesh needs a little prep before refinement, and
2342  // may need us to keep track of ghost node movement.
2343  std::unordered_set<dof_id_type> moved_ghost_nodes;
2344  if (!is_replicated)
2345  mesh.prepare_for_use();
2346 
2347  mesh_refinement.uniformly_refine(1);
2348 
2349  for (const auto & elem : mesh.active_element_ptr_range())
2350  for (auto s : elem->side_index_range())
2351  if (elem->neighbor_ptr(s) == nullptr || (mesh.mesh_dimension() == 2 && !flat))
2352  {
2353  elem->build_side_ptr(side, s);
2354 
2355  // Pop each point to the sphere boundary. Keep track of
2356  // any points we don't own, so we can push their "moved"
2357  // status to their owners.
2358  for (auto n : side->node_index_range())
2359  {
2360  Node & side_node = side->node_ref(n);
2361  side_node =
2362  sphere.closest_point(side->point(n));
2363 
2364  if (!is_replicated &&
2365  side_node.processor_id() != mesh.processor_id())
2366  moved_ghost_nodes.insert(side_node.id());
2367  }
2368  }
2369 
2370  if (!is_replicated)
2371  {
2372  std::map<processor_id_type, std::vector<dof_id_type>> moved_nodes_map;
2373  for (auto id : moved_ghost_nodes)
2374  {
2375  const Node & node = mesh.node_ref(id);
2376  moved_nodes_map[node.processor_id()].push_back(node.id());
2377  }
2378 
2379  auto action_functor =
2380  [& mesh, & sphere]
2381  (processor_id_type /* pid */,
2382  const std::vector<dof_id_type> & my_moved_nodes)
2383  {
2384  for (auto id : my_moved_nodes)
2385  {
2386  Node & node = mesh.node_ref(id);
2387  node = sphere.closest_point(node);
2388  }
2389  };
2390 
2391  // First get new node positions to their owners
2392  Parallel::push_parallel_vector_data
2393  (mesh.comm(), moved_nodes_map, action_functor);
2394 
2395  // Then get node positions to anyone else with them ghosted
2396  SyncNodalPositions sync_object(mesh);
2397  Parallel::sync_dofobject_data_by_id
2398  (mesh.comm(), mesh.nodes_begin(), mesh.nodes_end(),
2399  sync_object);
2400  }
2401  }
2402 
2403  // A DistributedMesh needs a little prep before flattening
2404  if (!is_replicated)
2405  mesh.prepare_for_use();
2406 
2407  // The mesh now contains a refinement hierarchy due to the refinements
2408  // used to generate the grid. In order to call other support functions
2409  // like all_tri() and all_second_order, you need a "flat" mesh file (with no
2410  // refinement trees) so
2411  MeshTools::Modification::flatten(mesh);
2412 
2413  // Convert all the tensor product elements to simplices if requested
2414  if ((type == TRI7) || (type == TRI6) || (type == TRI3) ||
2415  (type == TET4) || (type == TET10) || (type == TET14))
2416  {
2417  // A DistributedMesh needs a little prep before all_tri()
2418  if (is_replicated)
2419  mesh.prepare_for_use();
2420 
2421  MeshTools::Modification::all_tri(mesh);
2422  }
2423 
2424  // Convert to second-order elements if the user requested it.
2425  if (Elem::build(type)->default_order() != FIRST)
2426  {
2427  if (type == TET14)
2428  mesh.all_complete_order();
2429  else
2430  {
2431  // type is second-order, determine if it is the
2432  // "full-ordered" second-order element, or the "serendipity"
2433  // second order element. Note also that all_second_order
2434  // can't be called once the mesh has been refined.
2435  bool full_ordered = !((type==QUAD8) || (type==HEX20));
2436  mesh.all_second_order(full_ordered);
2437  }
2438 
2439  // And pop to the boundary again...
2440  for (const auto & elem : mesh.active_element_ptr_range())
2441  for (auto s : elem->side_index_range())
2442  if (elem->neighbor_ptr(s) == nullptr)
2443  {
2444  elem->build_side_ptr(side, s);
2445 
2446  // Pop each point to the sphere boundary
2447  for (auto n : side->node_index_range())
2448  side->point(n) =
2449  sphere.closest_point(side->point(n));
2450  }
2451  }
2452 
2453 
2454  // The meshes could probably use some smoothing.
2455  if (mesh.mesh_dimension() > 1)
2456  {
2457  LaplaceMeshSmoother smoother(mesh, n_smooth);
2458  smoother.smooth();
2459  }
2460 
2461  // We'll give the whole sphere surface a boundary id of 0
2462  for (const auto & elem : mesh.active_element_ptr_range())
2463  for (auto s : elem->side_index_range())
2464  if (!elem->neighbor_ptr(s))
2465  boundary_info.add_side(elem, s, 0);
2466 
2467  // Done building the mesh. Now prepare it for use.
2468  mesh.prepare_for_use();
2469 }
2470 
2471 #endif // #ifndef LIBMESH_ENABLE_AMR
2472 
2473 
2474 // Meshes the tensor product of a 1D and a 1D-or-2D domain.
2475 void MeshTools::Generation::build_extrusion (UnstructuredMesh & mesh,
2476  const MeshBase & cross_section,
2477  const unsigned int nz,
2478  RealVectorValue extrusion_vector,
2479  QueryElemSubdomainIDBase * elem_subdomain)
2480 {
2481  LOG_SCOPE("build_extrusion()", "MeshTools::Generation");
2482 
2483  if (!cross_section.n_elem())
2484  return;
2485 
2486  dof_id_type orig_elem = cross_section.n_elem();
2487  dof_id_type orig_nodes = cross_section.n_nodes();
2488 
2489 #ifdef LIBMESH_ENABLE_UNIQUE_ID
2490  unique_id_type orig_unique_ids = cross_section.parallel_max_unique_id();
2491 #endif
2492 
2493  unsigned int order = 1;
2494 
2495  BoundaryInfo & boundary_info = mesh.get_boundary_info();
2496  const BoundaryInfo & cross_section_boundary_info = cross_section.get_boundary_info();
2497 
2498  // Copy name maps from old to new boundary. We won't copy the whole
2499  // BoundaryInfo because that copies bc ids too, and we need to set
2500  // those more carefully.
2501  boundary_info.set_sideset_name_map() = cross_section_boundary_info.get_sideset_name_map();
2502  boundary_info.set_nodeset_name_map() = cross_section_boundary_info.get_nodeset_name_map();
2503  boundary_info.set_edgeset_name_map() = cross_section_boundary_info.get_edgeset_name_map();
2504 
2505  // If cross_section is distributed, so is its extrusion
2506  if (!cross_section.is_serial())
2507  mesh.delete_remote_elements();
2508 
2509  // We know a priori how many elements we'll need
2510  mesh.reserve_elem(nz*orig_elem);
2511 
2512  // For straightforward meshes we need one or two additional layers per
2513  // element.
2514  if (cross_section.elements_begin() != cross_section.elements_end() &&
2515  (*cross_section.elements_begin())->default_order() == SECOND)
2516  order = 2;
2517  mesh.comm().max(order);
2518 
2519  mesh.reserve_nodes((order*nz+1)*orig_nodes);
2520 
2521  // Container to catch the boundary IDs handed back by the BoundaryInfo object
2522  std::vector<boundary_id_type> ids_to_copy;
2523 
2524  for (const auto & node : cross_section.node_ptr_range())
2525  {
2526  for (unsigned int k=0; k != order*nz+1; ++k)
2527  {
2528  const dof_id_type new_node_id = node->id() + k * orig_nodes;
2529  Node * my_node = mesh.query_node_ptr(new_node_id);
2530  if (!my_node)
2531  {
2532  std::unique_ptr<Node> new_node = Node::build
2533  (*node + (extrusion_vector * k / nz / order),
2534  new_node_id);
2535  new_node->processor_id() = node->processor_id();
2536 
2537 #ifdef LIBMESH_ENABLE_UNIQUE_ID
2538  // Let's give the base of the extruded mesh the same
2539  // unique_ids as the source mesh, in case anyone finds that
2540  // a useful map to preserve.
2541  const unique_id_type uid = (k == 0) ?
2542  node->unique_id() :
2543  orig_unique_ids + (k-1)*(orig_nodes + orig_elem) + node->id();
2544 
2545  new_node->set_unique_id(uid);
2546 #endif
2547 
2548  cross_section_boundary_info.boundary_ids(node, ids_to_copy);
2549  boundary_info.add_node(new_node.get(), ids_to_copy);
2550 
2551  mesh.add_node(std::move(new_node));
2552  }
2553  }
2554  }
2555 
2556  const std::set<boundary_id_type> & side_ids =
2557  cross_section_boundary_info.get_side_boundary_ids();
2558 
2559  boundary_id_type next_side_id = side_ids.empty() ?
2560  0 : cast_int<boundary_id_type>(*side_ids.rbegin() + 1);
2561 
2562  // side_ids may not include ids from remote elements, in which case
2563  // some processors may have underestimated the next_side_id; let's
2564  // fix that.
2565  cross_section.comm().max(next_side_id);
2566 
2567  for (const auto & elem : cross_section.element_ptr_range())
2568  {
2569  const ElemType etype = elem->type();
2570 
2571  // build_extrusion currently only works on coarse meshes
2572  libmesh_assert (!elem->parent());
2573 
2574  for (unsigned int k=0; k != nz; ++k)
2575  {
2576  std::unique_ptr<Elem> new_elem;
2577  switch (etype)
2578  {
2579  case EDGE2:
2580  {
2581  new_elem = Elem::build(QUAD4);
2582  new_elem->set_node(0, mesh.node_ptr(elem->node_ptr(0)->id() + (k * orig_nodes)));
2583  new_elem->set_node(1, mesh.node_ptr(elem->node_ptr(1)->id() + (k * orig_nodes)));
2584  new_elem->set_node(2, mesh.node_ptr(elem->node_ptr(1)->id() + ((k+1) * orig_nodes)));
2585  new_elem->set_node(3, mesh.node_ptr(elem->node_ptr(0)->id() + ((k+1) * orig_nodes)));
2586 
2587  if (elem->neighbor_ptr(0) == remote_elem)
2588  new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
2589  if (elem->neighbor_ptr(1) == remote_elem)
2590  new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
2591 
2592  break;
2593  }
2594  case EDGE3:
2595  {
2596  new_elem = Elem::build(QUAD9);
2597  new_elem->set_node(0, mesh.node_ptr(elem->node_ptr(0)->id() + (2*k * orig_nodes)));
2598  new_elem->set_node(1, mesh.node_ptr(elem->node_ptr(1)->id() + (2*k * orig_nodes)));
2599  new_elem->set_node(2, mesh.node_ptr(elem->node_ptr(1)->id() + ((2*k+2) * orig_nodes)));
2600  new_elem->set_node(3, mesh.node_ptr(elem->node_ptr(0)->id() + ((2*k+2) * orig_nodes)));
2601  new_elem->set_node(4, mesh.node_ptr(elem->node_ptr(2)->id() + (2*k * orig_nodes)));
2602  new_elem->set_node(5, mesh.node_ptr(elem->node_ptr(1)->id() + ((2*k+1) * orig_nodes)));
2603  new_elem->set_node(6, mesh.node_ptr(elem->node_ptr(2)->id() + ((2*k+2) * orig_nodes)));
2604  new_elem->set_node(7, mesh.node_ptr(elem->node_ptr(0)->id() + ((2*k+1) * orig_nodes)));
2605  new_elem->set_node(8, mesh.node_ptr(elem->node_ptr(2)->id() + ((2*k+1) * orig_nodes)));
2606 
2607  if (elem->neighbor_ptr(0) == remote_elem)
2608  new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
2609  if (elem->neighbor_ptr(1) == remote_elem)
2610  new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
2611 
2612  break;
2613  }
2614  case TRI3:
2615  {
2616  new_elem = Elem::build(PRISM6);
2617  new_elem->set_node(0, mesh.node_ptr(elem->node_ptr(0)->id() + (k * orig_nodes)));
2618  new_elem->set_node(1, mesh.node_ptr(elem->node_ptr(1)->id() + (k * orig_nodes)));
2619  new_elem->set_node(2, mesh.node_ptr(elem->node_ptr(2)->id() + (k * orig_nodes)));
2620  new_elem->set_node(3, mesh.node_ptr(elem->node_ptr(0)->id() + ((k+1) * orig_nodes)));
2621  new_elem->set_node(4, mesh.node_ptr(elem->node_ptr(1)->id() + ((k+1) * orig_nodes)));
2622  new_elem->set_node(5, mesh.node_ptr(elem->node_ptr(2)->id() + ((k+1) * orig_nodes)));
2623 
2624  if (elem->neighbor_ptr(0) == remote_elem)
2625  new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
2626  if (elem->neighbor_ptr(1) == remote_elem)
2627  new_elem->set_neighbor(2, const_cast<RemoteElem *>(remote_elem));
2628  if (elem->neighbor_ptr(2) == remote_elem)
2629  new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
2630 
2631  break;
2632  }
2633  case TRI6:
2634  {
2635  new_elem = Elem::build(PRISM18);
2636  new_elem->set_node(0, mesh.node_ptr(elem->node_ptr(0)->id() + (2*k * orig_nodes)));
2637  new_elem->set_node(1, mesh.node_ptr(elem->node_ptr(1)->id() + (2*k * orig_nodes)));
2638  new_elem->set_node(2, mesh.node_ptr(elem->node_ptr(2)->id() + (2*k * orig_nodes)));
2639  new_elem->set_node(3, mesh.node_ptr(elem->node_ptr(0)->id() + ((2*k+2) * orig_nodes)));
2640  new_elem->set_node(4, mesh.node_ptr(elem->node_ptr(1)->id() + ((2*k+2) * orig_nodes)));
2641  new_elem->set_node(5, mesh.node_ptr(elem->node_ptr(2)->id() + ((2*k+2) * orig_nodes)));
2642  new_elem->set_node(6, mesh.node_ptr(elem->node_ptr(3)->id() + (2*k * orig_nodes)));
2643  new_elem->set_node(7, mesh.node_ptr(elem->node_ptr(4)->id() + (2*k * orig_nodes)));
2644  new_elem->set_node(8, mesh.node_ptr(elem->node_ptr(5)->id() + (2*k * orig_nodes)));
2645  new_elem->set_node(9, mesh.node_ptr(elem->node_ptr(0)->id() + ((2*k+1) * orig_nodes)));
2646  new_elem->set_node(10, mesh.node_ptr(elem->node_ptr(1)->id() + ((2*k+1) * orig_nodes)));
2647  new_elem->set_node(11, mesh.node_ptr(elem->node_ptr(2)->id() + ((2*k+1) * orig_nodes)));
2648  new_elem->set_node(12, mesh.node_ptr(elem->node_ptr(3)->id() + ((2*k+2) * orig_nodes)));
2649  new_elem->set_node(13, mesh.node_ptr(elem->node_ptr(4)->id() + ((2*k+2) * orig_nodes)));
2650  new_elem->set_node(14, mesh.node_ptr(elem->node_ptr(5)->id() + ((2*k+2) * orig_nodes)));
2651  new_elem->set_node(15, mesh.node_ptr(elem->node_ptr(3)->id() + ((2*k+1) * orig_nodes)));
2652  new_elem->set_node(16, mesh.node_ptr(elem->node_ptr(4)->id() + ((2*k+1) * orig_nodes)));
2653  new_elem->set_node(17, mesh.node_ptr(elem->node_ptr(5)->id() + ((2*k+1) * orig_nodes)));
2654 
2655  if (elem->neighbor_ptr(0) == remote_elem)
2656  new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
2657  if (elem->neighbor_ptr(1) == remote_elem)
2658  new_elem->set_neighbor(2, const_cast<RemoteElem *>(remote_elem));
2659  if (elem->neighbor_ptr(2) == remote_elem)
2660  new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
2661 
2662  break;
2663  }
2664  case TRI7:
2665  {
2666  new_elem = Elem::build(PRISM21);
2667  new_elem->set_node(0, mesh.node_ptr(elem->node_ptr(0)->id() + (2*k * orig_nodes)));
2668  new_elem->set_node(1, mesh.node_ptr(elem->node_ptr(1)->id() + (2*k * orig_nodes)));
2669  new_elem->set_node(2, mesh.node_ptr(elem->node_ptr(2)->id() + (2*k * orig_nodes)));
2670  new_elem->set_node(3, mesh.node_ptr(elem->node_ptr(0)->id() + ((2*k+2) * orig_nodes)));
2671  new_elem->set_node(4, mesh.node_ptr(elem->node_ptr(1)->id() + ((2*k+2) * orig_nodes)));
2672  new_elem->set_node(5, mesh.node_ptr(elem->node_ptr(2)->id() + ((2*k+2) * orig_nodes)));
2673  new_elem->set_node(6, mesh.node_ptr(elem->node_ptr(3)->id() + (2*k * orig_nodes)));
2674  new_elem->set_node(7, mesh.node_ptr(elem->node_ptr(4)->id() + (2*k * orig_nodes)));
2675  new_elem->set_node(8, mesh.node_ptr(elem->node_ptr(5)->id() + (2*k * orig_nodes)));
2676  new_elem->set_node(9, mesh.node_ptr(elem->node_ptr(0)->id() + ((2*k+1) * orig_nodes)));
2677  new_elem->set_node(10, mesh.node_ptr(elem->node_ptr(1)->id() + ((2*k+1) * orig_nodes)));
2678  new_elem->set_node(11, mesh.node_ptr(elem->node_ptr(2)->id() + ((2*k+1) * orig_nodes)));
2679  new_elem->set_node(12, mesh.node_ptr(elem->node_ptr(3)->id() + ((2*k+2) * orig_nodes)));
2680  new_elem->set_node(13, mesh.node_ptr(elem->node_ptr(4)->id() + ((2*k+2) * orig_nodes)));
2681  new_elem->set_node(14, mesh.node_ptr(elem->node_ptr(5)->id() + ((2*k+2) * orig_nodes)));
2682  new_elem->set_node(15, mesh.node_ptr(elem->node_ptr(3)->id() + ((2*k+1) * orig_nodes)));
2683  new_elem->set_node(16, mesh.node_ptr(elem->node_ptr(4)->id() + ((2*k+1) * orig_nodes)));
2684  new_elem->set_node(17, mesh.node_ptr(elem->node_ptr(5)->id() + ((2*k+1) * orig_nodes)));
2685 
2686  new_elem->set_node(18, mesh.node_ptr(elem->node_ptr(6)->id() + (2*k * orig_nodes)));
2687  new_elem->set_node(19, mesh.node_ptr(elem->node_ptr(6)->id() + ((2*k+2) * orig_nodes)));
2688  new_elem->set_node(20, mesh.node_ptr(elem->node_ptr(6)->id() + ((2*k+1) * orig_nodes)));
2689 
2690  if (elem->neighbor_ptr(0) == remote_elem)
2691  new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
2692  if (elem->neighbor_ptr(1) == remote_elem)
2693  new_elem->set_neighbor(2, const_cast<RemoteElem *>(remote_elem));
2694  if (elem->neighbor_ptr(2) == remote_elem)
2695  new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
2696 
2697  break;
2698  }
2699  case QUAD4:
2700  {
2701  new_elem = Elem::build(HEX8);
2702  new_elem->set_node(0, mesh.node_ptr(elem->node_ptr(0)->id() + (k * orig_nodes)));
2703  new_elem->set_node(1, mesh.node_ptr(elem->node_ptr(1)->id() + (k * orig_nodes)));
2704  new_elem->set_node(2, mesh.node_ptr(elem->node_ptr(2)->id() + (k * orig_nodes)));
2705  new_elem->set_node(3, mesh.node_ptr(elem->node_ptr(3)->id() + (k * orig_nodes)));
2706  new_elem->set_node(4, mesh.node_ptr(elem->node_ptr(0)->id() + ((k+1) * orig_nodes)));
2707  new_elem->set_node(5, mesh.node_ptr(elem->node_ptr(1)->id() + ((k+1) * orig_nodes)));
2708  new_elem->set_node(6, mesh.node_ptr(elem->node_ptr(2)->id() + ((k+1) * orig_nodes)));
2709  new_elem->set_node(7, mesh.node_ptr(elem->node_ptr(3)->id() + ((k+1) * orig_nodes)));
2710 
2711  if (elem->neighbor_ptr(0) == remote_elem)
2712  new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
2713  if (elem->neighbor_ptr(1) == remote_elem)
2714  new_elem->set_neighbor(2, const_cast<RemoteElem *>(remote_elem));
2715  if (elem->neighbor_ptr(2) == remote_elem)
2716  new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
2717  if (elem->neighbor_ptr(3) == remote_elem)
2718  new_elem->set_neighbor(4, const_cast<RemoteElem *>(remote_elem));
2719 
2720  break;
2721  }
2722  case QUAD9:
2723  {
2724  new_elem = Elem::build(HEX27);
2725  new_elem->set_node(0, mesh.node_ptr(elem->node_ptr(0)->id() + (2*k * orig_nodes)));
2726  new_elem->set_node(1, mesh.node_ptr(elem->node_ptr(1)->id() + (2*k * orig_nodes)));
2727  new_elem->set_node(2, mesh.node_ptr(elem->node_ptr(2)->id() + (2*k * orig_nodes)));
2728  new_elem->set_node(3, mesh.node_ptr(elem->node_ptr(3)->id() + (2*k * orig_nodes)));
2729  new_elem->set_node(4, mesh.node_ptr(elem->node_ptr(0)->id() + ((2*k+2) * orig_nodes)));
2730  new_elem->set_node(5, mesh.node_ptr(elem->node_ptr(1)->id() + ((2*k+2) * orig_nodes)));
2731  new_elem->set_node(6, mesh.node_ptr(elem->node_ptr(2)->id() + ((2*k+2) * orig_nodes)));
2732  new_elem->set_node(7, mesh.node_ptr(elem->node_ptr(3)->id() + ((2*k+2) * orig_nodes)));
2733  new_elem->set_node(8, mesh.node_ptr(elem->node_ptr(4)->id() + (2*k * orig_nodes)));
2734  new_elem->set_node(9, mesh.node_ptr(elem->node_ptr(5)->id() + (2*k * orig_nodes)));
2735  new_elem->set_node(10, mesh.node_ptr(elem->node_ptr(6)->id() + (2*k * orig_nodes)));
2736  new_elem->set_node(11, mesh.node_ptr(elem->node_ptr(7)->id() + (2*k * orig_nodes)));
2737  new_elem->set_node(12, mesh.node_ptr(elem->node_ptr(0)->id() + ((2*k+1) * orig_nodes)));
2738  new_elem->set_node(13, mesh.node_ptr(elem->node_ptr(1)->id() + ((2*k+1) * orig_nodes)));
2739  new_elem->set_node(14, mesh.node_ptr(elem->node_ptr(2)->id() + ((2*k+1) * orig_nodes)));
2740  new_elem->set_node(15, mesh.node_ptr(elem->node_ptr(3)->id() + ((2*k+1) * orig_nodes)));
2741  new_elem->set_node(16, mesh.node_ptr(elem->node_ptr(4)->id() + ((2*k+2) * orig_nodes)));
2742  new_elem->set_node(17, mesh.node_ptr(elem->node_ptr(5)->id() + ((2*k+2) * orig_nodes)));
2743  new_elem->set_node(18, mesh.node_ptr(elem->node_ptr(6)->id() + ((2*k+2) * orig_nodes)));
2744  new_elem->set_node(19, mesh.node_ptr(elem->node_ptr(7)->id() + ((2*k+2) * orig_nodes)));
2745  new_elem->set_node(20, mesh.node_ptr(elem->node_ptr(8)->id() + (2*k * orig_nodes)));
2746  new_elem->set_node(21, mesh.node_ptr(elem->node_ptr(4)->id() + ((2*k+1) * orig_nodes)));
2747  new_elem->set_node(22, mesh.node_ptr(elem->node_ptr(5)->id() + ((2*k+1) * orig_nodes)));
2748  new_elem->set_node(23, mesh.node_ptr(elem->node_ptr(6)->id() + ((2*k+1) * orig_nodes)));
2749  new_elem->set_node(24, mesh.node_ptr(elem->node_ptr(7)->id() + ((2*k+1) * orig_nodes)));
2750  new_elem->set_node(25, mesh.node_ptr(elem->node_ptr(8)->id() + ((2*k+2) * orig_nodes)));
2751  new_elem->set_node(26, mesh.node_ptr(elem->node_ptr(8)->id() + ((2*k+1) * orig_nodes)));
2752 
2753  if (elem->neighbor_ptr(0) == remote_elem)
2754  new_elem->set_neighbor(1, const_cast<RemoteElem *>(remote_elem));
2755  if (elem->neighbor_ptr(1) == remote_elem)
2756  new_elem->set_neighbor(2, const_cast<RemoteElem *>(remote_elem));
2757  if (elem->neighbor_ptr(2) == remote_elem)
2758  new_elem->set_neighbor(3, const_cast<RemoteElem *>(remote_elem));
2759  if (elem->neighbor_ptr(3) == remote_elem)
2760  new_elem->set_neighbor(4, const_cast<RemoteElem *>(remote_elem));
2761 
2762  break;
2763  }
2764  default:
2765  {
2766  libmesh_not_implemented();
2767  break;
2768  }
2769  }
2770 
2771  new_elem->set_id(elem->id() + (k * orig_elem));
2772  new_elem->processor_id() = elem->processor_id();
2773 
2774 #ifdef LIBMESH_ENABLE_UNIQUE_ID
2775  // Let's give the base of the extruded mesh the same
2776  // unique_ids as the source mesh, in case anyone finds that
2777  // a useful map to preserve.
2778  const unique_id_type uid = (k == 0) ?
2779  elem->unique_id() :
2780  orig_unique_ids + (k-1)*(orig_nodes + orig_elem) + orig_nodes + elem->id();
2781 
2782  new_elem->set_unique_id(uid);
2783 #endif
2784 
2785  if (!elem_subdomain)
2786  // maintain the subdomain_id
2787  new_elem->subdomain_id() = elem->subdomain_id();
2788  else
2789  // Allow the user to choose new subdomain_ids
2790  new_elem->subdomain_id() = elem_subdomain->get_subdomain_for_layer(elem, k);
2791 
2792  Elem * added_elem = mesh.add_elem(std::move(new_elem));
2793 
2794  // Copy any old boundary ids on all sides
2795  for (auto s : elem->side_index_range())
2796  {
2797  cross_section_boundary_info.boundary_ids(elem, s, ids_to_copy);
2798 
2799  if (added_elem->dim() == 3)
2800  {
2801  // For 2D->3D extrusion, we give the boundary IDs
2802  // for side s on the old element to side s+1 on the
2803  // new element. This is just a happy coincidence as
2804  // far as I can tell...
2805  boundary_info.add_side(added_elem,
2806  cast_int<unsigned short>(s+1),
2807  ids_to_copy);
2808  }
2809  else
2810  {
2811  // For 1D->2D extrusion, the boundary IDs map as:
2812  // Old elem -> New elem
2813  // 0 -> 3
2814  // 1 -> 1
2815  libmesh_assert_less(s, 2);
2816  const unsigned short sidemap[2] = {3, 1};
2817  boundary_info.add_side(added_elem, sidemap[s], ids_to_copy);
2818  }
2819  }
2820 
2821  // Give new boundary ids to bottom and top
2822  if (k == 0)
2823  boundary_info.add_side(added_elem, 0, next_side_id);
2824  if (k == nz-1)
2825  {
2826  // For 2D->3D extrusion, the "top" ID is 1+the original
2827  // element's number of sides. For 1D->2D extrusion, the
2828  // "top" ID is side 2.
2829  const unsigned short top_id = added_elem->dim() == 3 ?
2830  cast_int<unsigned short>(elem->n_sides()+1) : 2;
2831  boundary_info.add_side
2832  (added_elem, top_id,
2833  cast_int<boundary_id_type>(next_side_id+1));
2834  }
2835  }
2836  }
2837 
2838  // Done building the mesh. Now prepare it for use.
2839  mesh.prepare_for_use();
2840 }
2841 
2842 
2843 
2844 
2845 #if defined(LIBMESH_HAVE_TRIANGLE) && LIBMESH_DIM > 1
2846 
2847 // Triangulates a 2D rectangular region with or without holes
2848 void MeshTools::Generation::build_delaunay_square(UnstructuredMesh & mesh,
2849  const unsigned int nx, // num. of elements in x-dir
2850  const unsigned int ny, // num. of elements in y-dir
2851  const Real xmin, const Real xmax,
2852  const Real ymin, const Real ymax,
2853  const ElemType type,
2854  const std::vector<TriangleInterface::Hole*> * holes)
2855 {
2856  // Check for reasonable size
2857  libmesh_assert_greater_equal (nx, 1); // need at least 1 element in x-direction
2858  libmesh_assert_greater_equal (ny, 1); // need at least 1 element in y-direction
2859  libmesh_assert_less (xmin, xmax);
2860  libmesh_assert_less (ymin, ymax);
2861 
2862  // Clear out any data which may have been in the Mesh
2863  mesh.clear();
2864 
2865  BoundaryInfo & boundary_info = mesh.get_boundary_info();
2866 
2867  // Make sure the new Mesh will be 2D
2868  mesh.set_mesh_dimension(2);
2869 
2870  // The x and y spacing between boundary points
2871  const Real delta_x = (xmax-xmin) / static_cast<Real>(nx);
2872  const Real delta_y = (ymax-ymin) / static_cast<Real>(ny);
2873 
2874  // Bottom
2875  for (unsigned int p=0; p<=nx; ++p)
2876  mesh.add_point(Point(xmin + p*delta_x, ymin));
2877 
2878  // Right side
2879  for (unsigned int p=1; p<ny; ++p)
2880  mesh.add_point(Point(xmax, ymin + p*delta_y));
2881 
2882  // Top
2883  for (unsigned int p=0; p<=nx; ++p)
2884  mesh.add_point(Point(xmax - p*delta_x, ymax));
2885 
2886  // Left side
2887  for (unsigned int p=1; p<ny; ++p)
2888  mesh.add_point(Point(xmin, ymax - p*delta_y));
2889 
2890  // Be sure we added as many points as we thought we did
2891  libmesh_assert_equal_to (mesh.n_nodes(), 2*(nx+ny));
2892 
2893  // Construct the Triangle Interface object
2894  TriangleInterface t(mesh);
2895 
2896  // Set custom variables for the triangulation
2897  t.desired_area() = 0.5 * (xmax-xmin)*(ymax-ymin) / static_cast<Real>(nx*ny);
2898  t.triangulation_type() = TriangleInterface::PSLG;
2899  t.elem_type() = type;
2900 
2901  if (holes != nullptr)
2902  t.attach_hole_list(holes);
2903 
2904  // Triangulate!
2905  t.triangulate();
2906 
2907  // For avoiding extraneous side element construction
2908  std::unique_ptr<const Elem> side;
2909 
2910  // The mesh is now generated, but we still need to mark the boundaries
2911  // to be consistent with the other build_square routines. Note that all
2912  // hole boundary elements get the same ID, 4.
2913  for (auto & elem : mesh.element_ptr_range())
2914  for (auto s : elem->side_index_range())
2915  if (elem->neighbor_ptr(s) == nullptr)
2916  {
2917  elem->build_side_ptr(side, s);
2918 
2919  // Check the location of the side's midpoint. Since
2920  // the square has straight sides, the midpoint is not
2921  // on the corner and thus it is uniquely on one of the
2922  // sides.
2923  Point side_midpoint= 0.5f*( side->point(0) + side->point(1) );
2924 
2925  // The boundary ids are set following the same convention as Quad4 sides
2926  // bottom = 0
2927  // right = 1
2928  // top = 2
2929  // left = 3
2930  // hole = 4
2931  boundary_id_type bc_id=4;
2932 
2933  // bottom
2934  if (std::fabs(side_midpoint(1) - ymin) < TOLERANCE)
2935  bc_id=0;
2936 
2937  // right
2938  else if (std::fabs(side_midpoint(0) - xmax) < TOLERANCE)
2939  bc_id=1;
2940 
2941  // top
2942  else if (std::fabs(side_midpoint(1) - ymax) < TOLERANCE)
2943  bc_id=2;
2944 
2945  // left
2946  else if (std::fabs(side_midpoint(0) - xmin) < TOLERANCE)
2947  bc_id=3;
2948 
2949  // If the point is not on any of the external boundaries, it
2950  // is on one of the holes....
2951 
2952  // Finally, add this element's information to the boundary info object.
2953  boundary_info.add_side(elem->id(), s, bc_id);
2954  }
2955 
2956 } // end build_delaunay_square
2957 
2958 #endif // LIBMESH_HAVE_TRIANGLE && LIBMESH_DIM > 1
2959 
2960 
2961 void MeshTools::Generation::surface_octahedron
2962  (UnstructuredMesh & mesh,
2963  Real xmin, Real xmax,
2964  Real ymin, Real ymax,
2965  Real zmin, Real zmax,
2966  bool flip_tris)
2967 {
2968  const Real xavg = (xmin + xmax)/2;
2969  const Real yavg = (ymin + ymax)/2;
2970  const Real zavg = (zmin + zmax)/2;
2971  mesh.add_point(Point(xavg,yavg,zmin), 0);
2972  mesh.add_point(Point(xmax,yavg,zavg), 1);
2973  mesh.add_point(Point(xavg,ymax,zavg), 2);
2974  mesh.add_point(Point(xmin,yavg,zavg), 3);
2975  mesh.add_point(Point(xavg,ymin,zavg), 4);
2976  mesh.add_point(Point(xavg,yavg,zmax), 5);
2977 
2978  auto add_tri = [&mesh, flip_tris](std::array<dof_id_type,3> nodes)
2979  {
2980  auto elem = mesh.add_elem(Elem::build(TRI3));
2981  elem->set_node(0, mesh.node_ptr(nodes[0]));
2982  elem->set_node(1, mesh.node_ptr(nodes[1]));
2983  elem->set_node(2, mesh.node_ptr(nodes[2]));
2984  if (flip_tris)
2985  elem->flip(&mesh.get_boundary_info());
2986  };
2987 
2988  add_tri({0,2,1});
2989  add_tri({0,3,2});
2990  add_tri({0,4,3});
2991  add_tri({0,1,4});
2992  add_tri({5,4,1});
2993  add_tri({5,3,4});
2994  add_tri({5,2,3});
2995  add_tri({5,1,2});
2996 
2997  mesh.prepare_for_use();
2998 }
2999 
3000 
3001 } // namespace libMesh
libMesh::ElemType
ElemType
Defines an enum for geometric element types.
Definition: enum_elem_type.h:33
libMesh::MeshTools::Generation::build_extrusion
void build_extrusion(UnstructuredMesh &mesh, const MeshBase &cross_section, const unsigned int nz, RealVectorValue extrusion_vector, QueryElemSubdomainIDBase *elem_subdomain=nullptr)
Meshes the tensor product of a 1D and a 1D-or-2D domain.
Definition: mesh_generation.C:2475
libMesh::MeshTools::Generation::QueryElemSubdomainIDBase
Class for receiving the callback during extrusion generation and providing user-defined subdomains ba...
Definition: mesh_generation.h:163
libMesh::MeshBase::reserve_nodes
virtual void reserve_nodes(const dof_id_type nn)=0
Reserves space for a known number of nodes.
libMesh::TriangleInterface::triangulate
virtual void triangulate() override
Internally, this calls Triangle&#39;s triangulate routine.
Definition: mesh_triangle_interface.C:58
libMesh::BoundaryInfo::get_side_boundary_ids
const std::set< boundary_id_type > & get_side_boundary_ids() const
Definition: boundary_info.h:807
libMesh::Elem::set_node
virtual Node *& set_node(const unsigned int i)
Definition: elem.h:2564
libMesh::Node
A Node is like a Point, but with more information.
Definition: node.h:52
libMesh::MeshTools::Generation::Private::GaussLobattoRedistributionFunction
This object is passed to MeshTools::Modification::redistribute() to redistribute the points on a unif...
Definition: mesh_generation.C:162
libMesh::MeshTools::Generation::build_point
void build_point(UnstructuredMesh &mesh, const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
A specialized build_cube() for 0D meshes.
Definition: mesh_generation.C:1861
libMesh::MeshBase::parallel_max_unique_id
virtual unique_id_type parallel_max_unique_id() const =0
libMesh::BoundaryInfo::nodeset_name
std::string & nodeset_name(boundary_id_type id)
Definition: boundary_info.C:3261
libMesh::invalid_uint
const unsigned int invalid_uint
A number which is used quite often to represent an invalid or uninitialized value for an unsigned int...
Definition: libmesh.h:303
libMesh::RealVectorValue
VectorValue< Real > RealVectorValue
Useful typedefs to allow transparent switching between Real and Complex data types.
Definition: hp_coarsentest.h:46
libMesh::TOLERANCE
static constexpr Real TOLERANCE
Definition: libmesh_common.h:151
libMesh::DenseVector::resize
void resize(const unsigned int n)
Resize the vector.
Definition: dense_vector.h:396
libMesh::MeshBase::prepare_for_use
void prepare_for_use(const bool skip_renumber_nodes_and_elements, const bool skip_find_neighbors)
Prepare a newly created (or read) mesh for use.
Definition: mesh_base.C:824
libMesh::BoundaryInfo::remove
void remove(const Node *node)
Removes the boundary conditions associated with node node, if any exist.
Definition: boundary_info.C:1786
libMesh::BoundaryInfo::get_sideset_name_map
const std::map< boundary_id_type, std::string > & get_sideset_name_map() const
Definition: boundary_info.h:889
libMesh::Elem
This is the base class from which all geometric element types are derived.
Definition: elem.h:94
mesh
MeshBase & mesh
Definition: mesh_communication.C:1458
libMesh::ParallelObject::comm
const Parallel::Communicator & comm() const
Definition: parallel_object.h:97
libMesh::MeshTools::Generation::build_square
void build_square(UnstructuredMesh &mesh, const unsigned int nx, const unsigned int ny, const Real xmin=0., const Real xmax=1., const Real ymin=0., const Real ymax=1., const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
A specialized build_cube() for 2D meshes.
Definition: mesh_generation.C:1900
libMesh::BoundaryInfo::boundary_ids
void boundary_ids(const Node *node, std::vector< boundary_id_type > &vec_to_fill) const
Fills a user-provided std::vector with the boundary ids associated with Node node.
Definition: boundary_info.C:1348
libMesh
The libMesh namespace provides an interface to certain functionality in the library.
libMesh::TriangulatorInterface::elem_type
ElemType & elem_type()
Sets and/or gets the desired element type.
Definition: triangulator_interface.h:164
libMesh::MeshBase::get_boundary_info
const BoundaryInfo & get_boundary_info() const
The information about boundary ids on the mesh.
Definition: mesh_base.h:170
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
This is the MeshBase class.
Definition: mesh_base.h:80
libMesh::TriangulatorInterface::triangulation_type
TriangulationType & triangulation_type()
Sets and/or gets the desired triangulation type.
Definition: triangulator_interface.h:205
libMesh::MeshTools::Generation::Private::GaussLobattoRedistributionFunction::clone
virtual std::unique_ptr< FunctionBase< Real > > clone() const override
We must provide a way to clone ourselves to satisfy the pure virtual interface.
Definition: mesh_generation.C:218
processor_id_type
uint8_t processor_id_type
libMesh::MeshTools::Modification::all_tri
void all_tri(MeshBase &mesh)
Subdivides any non-simplex elements in a Mesh to produce simplex (triangular in 2D, tetrahedral in 3D) elements.
Definition: mesh_modification.C:434
libMesh::BoundaryInfo::set_sideset_name_map
std::map< boundary_id_type, std::string > & set_sideset_name_map()
Definition: boundary_info.h:887
libMesh::MeshBase::is_serial
virtual bool is_serial() const
Definition: mesh_base.h:347
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:954
libMesh::BoundaryInfo::get_nodeset_name_map
const std::map< boundary_id_type, std::string > & get_nodeset_name_map() const
Definition: boundary_info.h:897
libMesh::boundary_id_type
int8_t boundary_id_type
Definition: id_types.h:51
libMesh::BoundaryInfo::invalid_id
static const boundary_id_type invalid_id
Number used for internal use.
Definition: boundary_info.h:912
libMesh::DofObject::id
dof_id_type id() const
Definition: dof_object.h:819
libMesh::TriangulatorInterface::attach_hole_list
void attach_hole_list(const std::vector< Hole *> *holes)
Attaches a vector of Hole* pointers which will be meshed around.
Definition: triangulator_interface.h:255
libMesh::UnstructuredMesh
The UnstructuredMesh class is derived from the MeshBase class.
Definition: unstructured_mesh.h:48
libMesh::MeshBase::add_elem
virtual Elem * add_elem(Elem *e)=0
Add elem e to the end of the element array.
libMesh::Elem::build
static std::unique_ptr< Elem > build(const ElemType type, Elem *p=nullptr)
Definition: elem.C:442
libMesh::MeshBase::query_node_ptr
virtual const Node * query_node_ptr(const dof_id_type i) const =0
libMesh::BoundaryInfo::get_edgeset_name_map
const std::map< boundary_id_type, std::string > & get_edgeset_name_map() const
Definition: boundary_info.h:905
libMesh::BoundaryInfo
The BoundaryInfo class contains information relevant to boundary conditions including storing faces...
Definition: boundary_info.h:57
libMesh::libmesh_assert
libmesh_assert(ctx)
libMesh::MeshBase::set_mesh_dimension
void set_mesh_dimension(unsigned char d)
Resets the logical dimension of the mesh.
Definition: mesh_base.h:413
libMesh::Parallel::sync_dofobject_data_by_id
void sync_dofobject_data_by_id(const Communicator &comm, const Iterator &range_begin, const Iterator &range_end, SyncFunctor &sync)
Request data about a range of ghost dofobjects uniquely identified by their id.
Definition: parallel_ghost_sync.h:357
libMesh::MeshTools::Generation::Private::GaussLobattoRedistributionFunction::GaussLobattoRedistributionFunction
GaussLobattoRedistributionFunction(unsigned int nx, Real xmin, Real xmax, unsigned int ny=0, Real ymin=0, Real ymax=0, unsigned int nz=0, Real zmin=0, Real zmax=0)
Constructor.
Definition: mesh_generation.C:168
libMesh::MeshTools::Generation::Private::GaussLobattoRedistributionFunction::operator=
GaussLobattoRedistributionFunction & operator=(const GaussLobattoRedistributionFunction &)=default
libMesh::MeshBase::clear
virtual void clear()
Deletes all the element and node data that is currently stored.
Definition: mesh_base.C:1029
libMesh::MeshTools::Generation::surface_octahedron
void surface_octahedron(UnstructuredMesh &mesh, Real xmin, Real xmax, Real ymin, Real ymax, Real zmin, Real zmax, bool flip_tris=false)
Meshes the surface of an octahedron with 8 Tri3 elements, with counter-clockwise (libMesh default) no...
Definition: mesh_generation.C:2962
libMesh::BoundaryInfo::sideset_name
std::string & sideset_name(boundary_id_type id)
Definition: boundary_info.C:3246
top_id
const boundary_id_type top_id
Definition: periodic_bc_test.C:156
libMesh::Utility::enum_to_string
std::string enum_to_string(const T e)
libMesh::TriangulatorInterface::desired_area
Real & desired_area()
Sets and/or gets the desired triangle area.
Definition: triangulator_interface.h:177
libMesh::MeshBase::add_node
virtual Node * add_node(Node *n)=0
Add Node n to the end of the vertex array.
libMesh::TriangulatorInterface::PSLG
Triangulate the interior of a Planar Straight Line Graph, which is defined implicitly by the order of...
Definition: triangulator_interface.h:130
libMesh::Node::build
static std::unique_ptr< Node > build(const Node &n)
Definition: node.h:315
libMesh::MeshTools::Generation::build_sphere
void build_sphere(UnstructuredMesh &mesh, const Real radius=1, const unsigned int n_refinements=2, const ElemType type=INVALID_ELEM, const unsigned int n_smooth=2, const bool flat=true)
Meshes a spherical or mapped-spherical domain.
Definition: mesh_generation.C:1931
libMesh::Real
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
Definition: libmesh_common.h:144
TIMPI::Communicator::max
void max(const T &r, T &o, Request &req) const
libMesh::Elem::dim
virtual unsigned short dim() const =0
libMesh::MeshBase::all_complete_order
virtual void all_complete_order()
Calls the range-based version of this function with a range consisting of all elements in the mesh...
Definition: mesh_base.C:1801
libMesh::BoundaryInfo::set_nodeset_name_map
std::map< boundary_id_type, std::string > & set_nodeset_name_map()
Definition: boundary_info.h:895
libMesh::Elem::build_with_id
static std::unique_ptr< Elem > build_with_id(const ElemType type, dof_id_type id)
Calls the build() method above with a nullptr parent, and additionally sets the newly-created Elem&#39;s ...
Definition: elem.C:556
libMesh::MeshBase::is_replicated
virtual bool is_replicated() const
Definition: mesh_base.h:369
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:1210
libMesh::make_range
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
Definition: int_range.h:176
libMesh::MeshBase::all_second_order
void all_second_order(const bool full_ordered=true)
Calls the range-based version of this function with a range consisting of all elements in the mesh...
Definition: mesh_base.C:1796
libMesh::MeshTools::Generation::build_delaunay_square
void build_delaunay_square(UnstructuredMesh &mesh, const unsigned int nx, const unsigned int ny, const Real xmin, const Real xmax, const Real ymin, const Real ymax, const ElemType type, const std::vector< TriangleInterface::Hole *> *holes=nullptr)
Meshes a rectangular (2D) region (with or without holes) with a Delaunay triangulation.
Definition: mesh_generation.C:2848
libMesh::MeshBase::mesh_dimension
unsigned int mesh_dimension() const
Definition: mesh_base.C:430
libMesh::MeshTools::Generation::build_line
void build_line(UnstructuredMesh &mesh, const unsigned int nx, const Real xmin=0., const Real xmax=1., const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
A specialized build_cube() for 1D meshes.
Definition: mesh_generation.C:1879
libMesh::MeshTools::Modification::redistribute
void redistribute(MeshBase &mesh, const FunctionBase< Real > &mapfunc)
Deterministically perturb the nodal locations.
Definition: mesh_modification.C:289
libMesh::MeshBase::node_ref
virtual const Node & node_ref(const dof_id_type i) const
Definition: mesh_base.h:735
libMesh::DenseVector
Defines a dense vector for use in Finite Element-type computations.
Definition: dof_map.h:74
libMesh::MeshBase::delete_remote_elements
virtual void delete_remote_elements()
When supported, deletes all nonlocal elements of the mesh except for "ghosts" which touch a local ele...
Definition: mesh_base.h:389
libMesh::MeshTools::Generation::QueryElemSubdomainIDBase::get_subdomain_for_layer
virtual subdomain_id_type get_subdomain_for_layer(const Elem *old_elem, unsigned int layer)=0
libMesh::MeshBase::n_elem
virtual dof_id_type n_elem() const =0
libMesh::MeshBase::node_ptr
virtual const Node * node_ptr(const dof_id_type i) const =0
libMesh::FunctionBase
Base class for functors that can be evaluated at a point and (optionally) time.
Definition: dirichlet_boundaries.h:44
libMesh::ParallelObject::processor_id
processor_id_type processor_id() const
Definition: parallel_object.h:114
libMesh::TriangleInterface
A C++ interface between LibMesh and the Triangle library written by J.R.
Definition: mesh_triangle_interface.h:41
libMesh::MeshTools::Generation::Private::GaussLobattoRedistributionFunction::operator()
virtual void operator()(const Point &p, const Real, DenseVector< Real > &output) override
This is the actual function that MeshTools::Modification::redistribute() calls.
Definition: mesh_generation.C:228
libMesh::Point
A Point defines a location in LIBMESH_DIM dimensional Real space.
Definition: point.h:39
libMesh::MeshBase::reserve_elem
virtual void reserve_elem(const dof_id_type ne)=0
Reserves space for a known number of elements.
libMesh::unique_id_type
uint8_t unique_id_type
Definition: id_types.h:86
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:1967
int
void ErrorVector unsigned int
Definition: adjoints_ex3.C:360
libMesh::index_range
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
Definition: int_range.h:153
libMesh::MeshTools::Generation::build_cube
void build_cube(UnstructuredMesh &mesh, const unsigned int nx=0, const unsigned int ny=0, const unsigned int nz=0, const Real xmin=0., const Real xmax=1., const Real ymin=0., const Real ymax=1., const Real zmin=0., const Real zmax=1., const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
Builds a (elements) cube.
Definition: mesh_generation.C:314
libMesh::MeshBase::n_nodes
virtual dof_id_type n_nodes() const =0
libMesh::MeshTools::Generation::Private::idx
unsigned int idx(const ElemType type, const unsigned int nx, const unsigned int i, const unsigned int j)
A useful inline function which replaces the macros used previously.
Definition: mesh_generation.C:79
libMesh::dof_id_type
uint8_t dof_id_type
Definition: id_types.h:67
libMesh::BoundaryInfo::set_edgeset_name_map
std::map< boundary_id_type, std::string > & set_edgeset_name_map()
Definition: boundary_info.h:903
libMesh::pi
const Real pi
.
Definition: libmesh.h:292
libMesh::remote_elem
const RemoteElem * remote_elem
Definition: remote_elem.C:57
Generated on Wed Jun 24 2026 08:04:14 for libMesh by   doxygen 1.8.14