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