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