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