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