libMesh::MeshTools::Generation Namespace Reference

Namespaces

 Private
 

Functions

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)
 
void build_point (UnstructuredMesh &mesh, const ElemType type=INVALID_ELEM, const bool gauss_lobatto_grid=false)
 
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)
 
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)
 
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)
 
void build_extrusion (UnstructuredMesh &mesh, const MeshBase &cross_section, const unsigned int nz, RealVectorValue extrusion_vector)
 
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=NULL)
 

Detailed Description

Tools for Mesh generation.

Author
Benjamin S. Kirk
Date
2004

Function Documentation

void libMesh::MeshTools::Generation::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 $ nx \times ny \times nz $ (elements) cube. Defaults to a unit cube (or line in 1D, square in 2D), but the dimensions can be specified through the optional arguments.

Boundary ids are set to be equal to the side indexing on a master hex

Definition at line 150 of file mesh_generation.C.

References libMesh::MeshBase::add_elem(), libMesh::MeshBase::add_point(), libMesh::UnstructuredMesh::all_second_order(), libMesh::MeshBase::boundary_info, libMesh::Elem::build_side(), libMesh::MeshBase::clear(), libMesh::MeshBase::delete_elem(), libMeshEnums::EDGE2, libMeshEnums::EDGE3, libMeshEnums::EDGE4, libMesh::MeshBase::elements_begin(), libMesh::MeshBase::elements_end(), libMesh::err, libMesh::Elem::get_node(), libMeshEnums::HEX20, libMeshEnums::HEX27, libMeshEnums::HEX8, libMesh::MeshTools::Generation::Private::idx(), libMeshEnums::INVALID_ELEM, libMesh::BoundaryInfo::invalid_id, libMesh::libmesh_assert(), libMesh::MeshBase::mesh_dimension(), libMesh::MeshBase::n_elem(), libMesh::MeshBase::n_nodes(), libMesh::Elem::n_sides(), libMesh::MeshBase::node(), libMesh::MeshBase::node_ptr(), libMeshEnums::NODEELEM, libMesh::pi, libMesh::MeshBase::prepare_for_use(), libMeshEnums::PRISM15, libMeshEnums::PRISM18, libMeshEnums::PRISM6, libMeshEnums::PYRAMID14, libMeshEnums::PYRAMID5, libMeshEnums::QUAD4, libMeshEnums::QUAD8, libMeshEnums::QUAD9, libMesh::Real, libMesh::MeshBase::reserve_elem(), libMesh::MeshBase::reserve_nodes(), libMesh::MeshBase::set_mesh_dimension(), libMesh::Elem::set_node(), side, libMesh::START_LOG(), libMesh::STOP_LOG(), libMeshEnums::TET10, libMeshEnums::TET4, libMeshEnums::TRI3, and libMeshEnums::TRI6.

Referenced by build_line(), build_point(), and build_square().

159 {
160  START_LOG("build_cube()", "MeshTools::Generation");
161 
162  // Declare that we are using the indexing utility routine
163  // in the "Private" part of our current namespace. If this doesn't
164  // work in GCC 2.95.3 we can either remove it or stop supporting
165  // 2.95.3 altogether.
166  // Changing this to import the whole namespace... just importing idx
167  // causes an internal compiler error for Intel Compiler 11.0 on Linux
168  // in debug mode.
169  using namespace MeshTools::Generation::Private;
170 
171  // Clear the mesh and start from scratch
172  mesh.clear();
173 
174  if (nz != 0)
175  mesh.set_mesh_dimension(3);
176  else if (ny != 0)
177  mesh.set_mesh_dimension(2);
178  else if (nx != 0)
179  mesh.set_mesh_dimension(1);
180  else
181  mesh.set_mesh_dimension(0);
182 
183  switch (mesh.mesh_dimension())
184  {
185  //---------------------------------------------------------------------
186  // Build a 0D point
187  case 0:
188  {
189  libmesh_assert_equal_to (nx, 0);
190  libmesh_assert_equal_to (ny, 0);
191  libmesh_assert_equal_to (nz, 0);
192 
193  libmesh_assert (type == INVALID_ELEM || type == NODEELEM);
194 
195  // Build one nodal element for the mesh
196  mesh.add_point (Point(0, 0, 0), 0);
197  Elem* elem = mesh.add_elem (new NodeElem);
198  elem->set_node(0) = mesh.node_ptr(0);
199 
200  break;
201  }
202 
203 
204 
205  //---------------------------------------------------------------------
206  // Build a 1D line
207  case 1:
208  {
209  libmesh_assert_not_equal_to (nx, 0);
210  libmesh_assert_equal_to (ny, 0);
211  libmesh_assert_equal_to (nz, 0);
212  libmesh_assert_less (xmin, xmax);
213 
214  // Reserve elements
215  switch (type)
216  {
217  case INVALID_ELEM:
218  case EDGE2:
219  case EDGE3:
220  case EDGE4:
221  {
222  mesh.reserve_elem (nx);
223  break;
224  }
225 
226  default:
227  {
228  libMesh::err << "ERROR: Unrecognized 1D element type." << std::endl;
229  libmesh_error();
230  }
231  }
232 
233  // Reserve nodes
234  switch (type)
235  {
236  case INVALID_ELEM:
237  case EDGE2:
238  {
239  mesh.reserve_nodes(nx+1);
240  break;
241  }
242 
243  case EDGE3:
244  {
245  mesh.reserve_nodes(2*nx+1);
246  break;
247  }
248 
249  case EDGE4:
250  {
251  mesh.reserve_nodes(3*nx+1);
252  break;
253  }
254 
255  default:
256  {
257  libMesh::err << "ERROR: Unrecognized 1D element type." << std::endl;
258  libmesh_error();
259  }
260  }
261 
262 
263  // Build the nodes, depends on whether we're using linears,
264  // quadratics or cubics and whether using uniform grid or Gauss-Lobatto
265  unsigned int node_id = 0;
266  switch(type)
267  {
268  case INVALID_ELEM:
269  case EDGE2:
270  {
271  for (unsigned int i=0; i<=nx; i++)
272  {
273  if (gauss_lobatto_grid)
274  mesh.add_point (Point(0.5*(std::cos(libMesh::pi*static_cast<Real>(nx-i)/static_cast<Real>(nx))+1.0),
275  0,
276  0), node_id++);
277  else
278  mesh.add_point (Point(static_cast<Real>(i)/static_cast<Real>(nx),
279  0,
280  0), node_id++);
281  }
282  break;
283  }
284 
285  case EDGE3:
286  {
287  for (unsigned int i=0; i<=2*nx; i++)
288  {
289  if (gauss_lobatto_grid)
290  {
291  // The x location of the point.
292  Real x=0.;
293 
294  // Shortcut quantities (do not depend on i)
295  const Real c = std::cos( libMesh::pi*i / static_cast<Real>(2*nx) );
296 
297  // If i is even, compute a normal Gauss-Lobatto point
298  if (i%2 == 0)
299  x = 0.5*(1.0 - c);
300 
301  // Otherwise, it is the average of the previous and next points
302  else
303  {
304  Real cmin = std::cos( libMesh::pi*(i-1) / static_cast<Real>(2*nx) );
305  Real cmax = std::cos( libMesh::pi*(i+1) / static_cast<Real>(2*nx) );
306 
307  Real gl_xmin = 0.5*(1.0 - cmin);
308  Real gl_xmax = 0.5*(1.0 - cmax);
309  x = 0.5*(gl_xmin + gl_xmax);
310  }
311 
312  mesh.add_point (Point(x,0.,0.), node_id++);
313  }
314  else
315  mesh.add_point (Point(static_cast<Real>(i)/static_cast<Real>(2*nx),
316  0,
317  0), node_id++);
318  }
319  break;
320  }
321 
322  case EDGE4:
323  {
324  for (unsigned int i=0; i<=3*nx; i++)
325  {
326  if (gauss_lobatto_grid)
327  {
328  // The x location of the point
329  Real x=0.;
330 
331  // Shortcut quantities
332  const Real c = std::cos( libMesh::pi*i / static_cast<Real>(3*nx) );
333 
334  // If i is multiple of 3, compute a normal Gauss-Lobatto point
335  if (i%3 == 0)
336  x = 0.5*(1.0 - c);
337 
338  // Otherwise, distribute points evenly within the element
339  else
340  {
341  if(i%3 == 1)
342  {
343  Real cmin = std::cos( libMesh::pi*(i-1) / static_cast<Real>(3*nx) );
344  Real cmax = std::cos( libMesh::pi*(i+2) / static_cast<Real>(3*nx) );
345 
346  Real gl_xmin = 0.5*(1.0 - cmin);
347  Real gl_xmax = 0.5*(1.0 - cmax);
348 
349  x = (2.*gl_xmin + gl_xmax)/3.;
350  }
351  else
352  if(i%3 == 2)
353  {
354  Real cmin = std::cos( libMesh::pi*(i-2) / static_cast<Real>(3*nx) );
355  Real cmax = std::cos( libMesh::pi*(i+1) / static_cast<Real>(3*nx) );
356 
357  Real gl_xmin = 0.5*(1.0 - cmin);
358  Real gl_xmax = 0.5*(1.0 - cmax);
359 
360  x = (gl_xmin + 2.*gl_xmax)/3.;
361  }
362 
363  }
364 
365  mesh.add_point (Point(x,0.,0.), node_id++);
366  }
367  else
368  mesh.add_point (Point(static_cast<Real>(i)/static_cast<Real>(3*nx),
369  0,
370  0), node_id++);
371  }
372 
373 
374 
375  break;
376  }
377 
378  default:
379  {
380  libMesh::err << "ERROR: Unrecognized 1D element type." << std::endl;
381  libmesh_error();
382  }
383 
384  }
385 
386  // Build the elements of the mesh
387  switch(type)
388  {
389  case INVALID_ELEM:
390  case EDGE2:
391  {
392  for (unsigned int i=0; i<nx; i++)
393  {
394  Elem* elem = mesh.add_elem (new Edge2);
395  elem->set_node(0) = mesh.node_ptr(i);
396  elem->set_node(1) = mesh.node_ptr(i+1);
397 
398  if (i == 0)
399  mesh.boundary_info->add_side(elem, 0, 0);
400 
401  if (i == (nx-1))
402  mesh.boundary_info->add_side(elem, 1, 1);
403  }
404  break;
405  }
406 
407  case EDGE3:
408  {
409  for (unsigned int i=0; i<nx; i++)
410  {
411  Elem* elem = mesh.add_elem (new Edge3);
412  elem->set_node(0) = mesh.node_ptr(2*i);
413  elem->set_node(2) = mesh.node_ptr(2*i+1);
414  elem->set_node(1) = mesh.node_ptr(2*i+2);
415 
416  if (i == 0)
417  mesh.boundary_info->add_side(elem, 0, 0);
418 
419  if (i == (nx-1))
420  mesh.boundary_info->add_side(elem, 1, 1);
421  }
422  break;
423  }
424 
425  case EDGE4:
426  {
427  for (unsigned int i=0; i<nx; i++)
428  {
429  Elem* elem = mesh.add_elem (new Edge4);
430  elem->set_node(0) = mesh.node_ptr(3*i);
431  elem->set_node(2) = mesh.node_ptr(3*i+1);
432  elem->set_node(3) = mesh.node_ptr(3*i+2);
433  elem->set_node(1) = mesh.node_ptr(3*i+3);
434 
435  if (i == 0)
436  mesh.boundary_info->add_side(elem, 0, 0);
437 
438  if (i == (nx-1))
439  mesh.boundary_info->add_side(elem, 1, 1);
440  }
441  break;
442  }
443 
444  default:
445  {
446  libMesh::err << "ERROR: Unrecognized 1D element type." << std::endl;
447  libmesh_error();
448  }
449  }
450 
451  // Scale the nodal positions
452  for (unsigned int p=0; p<mesh.n_nodes(); p++)
453  mesh.node(p)(0) = (mesh.node(p)(0))*(xmax-xmin) + xmin;
454 
455  // Add sideset names to boundary info
456  mesh.boundary_info->sideset_name(0) = "left";
457  mesh.boundary_info->sideset_name(1) = "right";
458 
459  // Add nodeset names to boundary info
460  mesh.boundary_info->nodeset_name(0) = "left";
461  mesh.boundary_info->nodeset_name(1) = "right";
462 
463  break;
464  }
465 
466 
467 
468 
469 
470 
471 
472 
473 
474 
475  //---------------------------------------------------------------------
476  // Build a 2D quadrilateral
477  case 2:
478  {
479  libmesh_assert_not_equal_to (nx, 0);
480  libmesh_assert_not_equal_to (ny, 0);
481  libmesh_assert_equal_to (nz, 0);
482  libmesh_assert_less (xmin, xmax);
483  libmesh_assert_less (ymin, ymax);
484 
485  // Reserve elements. The TRI3 and TRI6 meshes
486  // have twice as many elements...
487  switch (type)
488  {
489  case INVALID_ELEM:
490  case QUAD4:
491  case QUAD8:
492  case QUAD9:
493  {
494  mesh.reserve_elem (nx*ny);
495  break;
496  }
497 
498  case TRI3:
499  case TRI6:
500  {
501  mesh.reserve_elem (2*nx*ny);
502  break;
503  }
504 
505  default:
506  {
507  libMesh::err << "ERROR: Unrecognized 2D element type." << std::endl;
508  libmesh_error();
509  }
510  }
511 
512 
513 
514  // Reserve nodes. The quadratic element types
515  // need to reserve more nodes than the linear types.
516  switch (type)
517  {
518  case INVALID_ELEM:
519  case QUAD4:
520  case TRI3:
521  {
522  mesh.reserve_nodes( (nx+1)*(ny+1) );
523  break;
524  }
525 
526  case QUAD8:
527  case QUAD9:
528  case TRI6:
529  {
530  mesh.reserve_nodes( (2*nx+1)*(2*ny+1) );
531  break;
532  }
533 
534 
535  default:
536  {
537  libMesh::err << "ERROR: Unrecognized 2D element type." << std::endl;
538  libmesh_error();
539  }
540  }
541 
542 
543 
544  // Build the nodes. Depends on whether you are using a linear
545  // or quadratic element, and whether you are using a uniform
546  // grid or the Gauss-Lobatto grid points.
547  unsigned int node_id = 0;
548  switch (type)
549  {
550  case INVALID_ELEM:
551  case QUAD4:
552  case TRI3:
553  {
554  for (unsigned int j=0; j<=ny; j++)
555  for (unsigned int i=0; i<=nx; i++)
556  {
557  if (gauss_lobatto_grid)
558  {
559  mesh.add_point (Point(0.5*(1.0 - std::cos(libMesh::pi*static_cast<Real>(i)/static_cast<Real>(nx))),
560  0.5*(1.0 - std::cos(libMesh::pi*static_cast<Real>(j)/static_cast<Real>(ny))),
561  0.), node_id++);
562  }
563 
564  else
565  mesh.add_point (Point(static_cast<Real>(i)/static_cast<Real>(nx),
566  static_cast<Real>(j)/static_cast<Real>(ny),
567  0.), node_id++);
568  }
569 
570  break;
571  }
572 
573  case QUAD8:
574  case QUAD9:
575  case TRI6:
576  {
577  for (unsigned int j=0; j<=(2*ny); j++)
578  for (unsigned int i=0; i<=(2*nx); i++)
579  {
580  if (gauss_lobatto_grid)
581  {
582  // The x,y locations of the point.
583  Real x=0., y=0.;
584 
585  // Shortcut quantities (do not depend on i,j)
586  const Real a = std::cos( libMesh::pi / static_cast<Real>(2*nx) );
587  const Real b = std::cos( libMesh::pi / static_cast<Real>(2*ny) );
588 
589  // Shortcut quantities (depend on i,j)
590  const Real c = std::cos( libMesh::pi*i / static_cast<Real>(2*nx) );
591  const Real d = std::cos( libMesh::pi*j / static_cast<Real>(2*ny) );
592 
593  // If i is even, compute a normal Gauss-Lobatto point
594  if (i%2 == 0)
595  x = 0.5*(1.0 - c);
596 
597  // Otherwise, it is the average of the previous and next points
598  else
599  x = 0.5*(1.0 - a*c);
600 
601  // If j is even, compute a normal Gauss-Lobatto point
602  if (j%2 == 0)
603  y = 0.5*(1.0 - d);
604 
605  // Otherwise, it is the average of the previous and next points
606  else
607  y = 0.5*(1.0 - b*d);
608 
609 
610  mesh.add_point (Point(x,y,0.), node_id++);
611  }
612 
613 
614  else
615  mesh.add_point (Point(static_cast<Real>(i)/static_cast<Real>(2*nx),
616  static_cast<Real>(j)/static_cast<Real>(2*ny),
617  0), node_id++);
618  }
619 
620  break;
621  }
622 
623 
624  default:
625  {
626  libMesh::err << "ERROR: Unrecognized 2D element type." << std::endl;
627  libmesh_error();
628  }
629  }
630 
631 
632 
633 
634 
635 
636  // Build the elements. Each one is a bit different.
637  switch (type)
638  {
639 
640  case INVALID_ELEM:
641  case QUAD4:
642  {
643  for (unsigned int j=0; j<ny; j++)
644  for (unsigned int i=0; i<nx; i++)
645  {
646  Elem* elem = mesh.add_elem(new Quad4);
647 
648  elem->set_node(0) = mesh.node_ptr(idx(type,nx,i,j) );
649  elem->set_node(1) = mesh.node_ptr(idx(type,nx,i+1,j) );
650  elem->set_node(2) = mesh.node_ptr(idx(type,nx,i+1,j+1));
651  elem->set_node(3) = mesh.node_ptr(idx(type,nx,i,j+1) );
652 
653  if (j == 0)
654  mesh.boundary_info->add_side(elem, 0, 0);
655 
656  if (j == (ny-1))
657  mesh.boundary_info->add_side(elem, 2, 2);
658 
659  if (i == 0)
660  mesh.boundary_info->add_side(elem, 3, 3);
661 
662  if (i == (nx-1))
663  mesh.boundary_info->add_side(elem, 1, 1);
664  }
665  break;
666  }
667 
668 
669  case TRI3:
670  {
671  for (unsigned int j=0; j<ny; j++)
672  for (unsigned int i=0; i<nx; i++)
673  {
674  Elem* elem = NULL;
675 
676  // Add first Tri3
677  elem = mesh.add_elem(new Tri3);
678 
679  elem->set_node(0) = mesh.node_ptr(idx(type,nx,i,j) );
680  elem->set_node(1) = mesh.node_ptr(idx(type,nx,i+1,j) );
681  elem->set_node(2) = mesh.node_ptr(idx(type,nx,i+1,j+1));
682 
683  if (j == 0)
684  mesh.boundary_info->add_side(elem, 0, 0);
685 
686  if (i == (nx-1))
687  mesh.boundary_info->add_side(elem, 1, 1);
688 
689  // Add second Tri3
690  elem = mesh.add_elem(new Tri3);
691 
692  elem->set_node(0) = mesh.node_ptr(idx(type,nx,i,j) );
693  elem->set_node(1) = mesh.node_ptr(idx(type,nx,i+1,j+1));
694  elem->set_node(2) = mesh.node_ptr(idx(type,nx,i,j+1) );
695 
696  if (j == (ny-1))
697  mesh.boundary_info->add_side(elem, 1, 2);
698 
699  if (i == 0)
700  mesh.boundary_info->add_side(elem, 2, 3);
701  }
702  break;
703  }
704 
705 
706 
707  case QUAD8:
708  case QUAD9:
709  {
710  for (unsigned int j=0; j<(2*ny); j += 2)
711  for (unsigned int i=0; i<(2*nx); i += 2)
712  {
713  Elem* elem = (type == QUAD8) ?
714  mesh.add_elem(new Quad8) :
715  mesh.add_elem(new Quad9);
716 
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+2,j) );
720  elem->set_node(2) = mesh.node_ptr(idx(type,nx,i+2,j+2));
721  elem->set_node(3) = mesh.node_ptr(idx(type,nx,i,j+2) );
722  elem->set_node(4) = mesh.node_ptr(idx(type,nx,i+1,j) );
723  elem->set_node(5) = mesh.node_ptr(idx(type,nx,i+2,j+1));
724  elem->set_node(6) = mesh.node_ptr(idx(type,nx,i+1,j+2));
725  elem->set_node(7) = mesh.node_ptr(idx(type,nx,i,j+1) );
726  if (type == QUAD9)
727  elem->set_node(8) = mesh.node_ptr(idx(type,nx,i+1,j+1));
728 
729 
730  if (j == 0)
731  mesh.boundary_info->add_side(elem, 0, 0);
732 
733  if (j == 2*(ny-1))
734  mesh.boundary_info->add_side(elem, 2, 2);
735 
736  if (i == 0)
737  mesh.boundary_info->add_side(elem, 3, 3);
738 
739  if (i == 2*(nx-1))
740  mesh.boundary_info->add_side(elem, 1, 1);
741  }
742  break;
743  }
744 
745 
746  case TRI6:
747  {
748  for (unsigned int j=0; j<(2*ny); j += 2)
749  for (unsigned int i=0; i<(2*nx); i += 2)
750  {
751  Elem* elem = NULL;
752 
753  // Add first Tri6
754  elem = mesh.add_elem(new Tri6);
755 
756  elem->set_node(0) = mesh.node_ptr(idx(type,nx,i,j) );
757  elem->set_node(1) = mesh.node_ptr(idx(type,nx,i+2,j) );
758  elem->set_node(2) = mesh.node_ptr(idx(type,nx,i+2,j+2));
759  elem->set_node(3) = mesh.node_ptr(idx(type,nx,i+1,j) );
760  elem->set_node(4) = mesh.node_ptr(idx(type,nx,i+2,j+1));
761  elem->set_node(5) = mesh.node_ptr(idx(type,nx,i+1,j+1));
762 
763  if (j == 0)
764  mesh.boundary_info->add_side(elem, 0, 0);
765 
766  if (i == 2*(nx-1))
767  mesh.boundary_info->add_side(elem, 1, 1);
768 
769  // Add second Tri6
770  elem = mesh.add_elem(new Tri6);
771 
772  elem->set_node(0) = mesh.node_ptr(idx(type,nx,i,j) );
773  elem->set_node(1) = mesh.node_ptr(idx(type,nx,i+2,j+2));
774  elem->set_node(2) = mesh.node_ptr(idx(type,nx,i,j+2) );
775  elem->set_node(3) = mesh.node_ptr(idx(type,nx,i+1,j+1));
776  elem->set_node(4) = mesh.node_ptr(idx(type,nx,i+1,j+2));
777  elem->set_node(5) = mesh.node_ptr(idx(type,nx,i,j+1) );
778 
779  if (j == 2*(ny-1))
780  mesh.boundary_info->add_side(elem, 1, 2);
781 
782  if (i == 0)
783  mesh.boundary_info->add_side(elem, 2, 3);
784 
785  }
786  break;
787  };
788 
789 
790  default:
791  {
792  libMesh::err << "ERROR: Unrecognized 2D element type." << std::endl;
793  libmesh_error();
794  }
795  }
796 
797 
798 
799 
800  // Scale the nodal positions
801  for (unsigned int p=0; p<mesh.n_nodes(); p++)
802  {
803  mesh.node(p)(0) = (mesh.node(p)(0))*(xmax-xmin) + xmin;
804  mesh.node(p)(1) = (mesh.node(p)(1))*(ymax-ymin) + ymin;
805  }
806 
807  // Add sideset names to boundary info
808  mesh.boundary_info->sideset_name(0) = "bottom";
809  mesh.boundary_info->sideset_name(1) = "right";
810  mesh.boundary_info->sideset_name(2) = "top";
811  mesh.boundary_info->sideset_name(3) = "left";
812 
813  // Add nodeset names to boundary info
814  mesh.boundary_info->nodeset_name(0) = "bottom";
815  mesh.boundary_info->nodeset_name(1) = "right";
816  mesh.boundary_info->nodeset_name(2) = "top";
817  mesh.boundary_info->nodeset_name(3) = "left";
818 
819  break;
820  }
821 
822 
823 
824 
825 
826 
827 
828 
829 
830 
831 
832  //---------------------------------------------------------------------
833  // Build a 3D mesh using hexes, tets, prisms, or pyramids.
834  case 3:
835  {
836  libmesh_assert_not_equal_to (nx, 0);
837  libmesh_assert_not_equal_to (ny, 0);
838  libmesh_assert_not_equal_to (nz, 0);
839  libmesh_assert_less (xmin, xmax);
840  libmesh_assert_less (ymin, ymax);
841  libmesh_assert_less (zmin, zmax);
842 
843 
844  // Reserve elements. Meshes with prismatic elements require
845  // twice as many elements.
846  switch (type)
847  {
848  case INVALID_ELEM:
849  case HEX8:
850  case HEX20:
851  case HEX27:
852  case TET4: // TET4's are created from an initial HEX27 discretization
853  case TET10: // TET10's are created from an initial HEX27 discretization
854  case PYRAMID5: // PYRAMID5's are created from an initial HEX27 discretization
855  case PYRAMID14:
856  {
857  mesh.reserve_elem(nx*ny*nz);
858  break;
859  }
860 
861  case PRISM6:
862  case PRISM15:
863  case PRISM18:
864  {
865  mesh.reserve_elem(2*nx*ny*nz);
866  break;
867  }
868 
869  default:
870  {
871  libMesh::err << "ERROR: Unrecognized 3D element type." << std::endl;
872  libmesh_error();
873  }
874  }
875 
876 
877 
878 
879 
880  // Reserve nodes. Quadratic elements need twice as many nodes as linear elements.
881  switch (type)
882  {
883  case INVALID_ELEM:
884  case HEX8:
885  case PRISM6:
886  {
887  mesh.reserve_nodes( (nx+1)*(ny+1)*(nz+1) );
888  break;
889  }
890 
891  case HEX20:
892  case HEX27:
893  case TET4: // TET4's are created from an initial HEX27 discretization
894  case TET10: // TET10's are created from an initial HEX27 discretization
895  case PYRAMID5: // PYRAMID5's are created from an initial HEX27 discretization
896  case PYRAMID14:
897  case PRISM15:
898  case PRISM18:
899  {
900  // FYI: The resulting TET4 mesh will have exactly
901  // 5*(nx*ny*nz) + 2*(nx*ny + nx*nz + ny*nz) + (nx+ny+nz) + 1
902  // nodes once the additional mid-edge nodes for the HEX27 discretization
903  // have been deleted.
904  mesh.reserve_nodes( (2*nx+1)*(2*ny+1)*(2*nz+1) );
905  break;
906  }
907 
908  default:
909  {
910  libMesh::err << "ERROR: Unrecognized 3D element type." << std::endl;
911  libmesh_error();
912  }
913  }
914 
915 
916 
917 
918  // Build the nodes.
919  unsigned int node_id = 0;
920  switch (type)
921  {
922  case INVALID_ELEM:
923  case HEX8:
924  case PRISM6:
925  {
926  for (unsigned int k=0; k<=nz; k++)
927  for (unsigned int j=0; j<=ny; j++)
928  for (unsigned int i=0; i<=nx; i++)
929  {
930  if (gauss_lobatto_grid)
931  {
932  mesh.add_point (Point(0.5*(1.0 - std::cos(libMesh::pi*static_cast<Real>(i)/static_cast<Real>(nx))),
933  0.5*(1.0 - std::cos(libMesh::pi*static_cast<Real>(j)/static_cast<Real>(ny))),
934  0.5*(1.0 - std::cos(libMesh::pi*static_cast<Real>(k)/static_cast<Real>(nz)))), node_id++);
935  }
936 
937  else
938  mesh.add_point(Point(static_cast<Real>(i)/static_cast<Real>(nx),
939  static_cast<Real>(j)/static_cast<Real>(ny),
940  static_cast<Real>(k)/static_cast<Real>(nz)), node_id++);
941  }
942  break;
943  }
944 
945  case HEX20:
946  case HEX27:
947  case TET4: // TET4's are created from an initial HEX27 discretization
948  case TET10: // TET10's are created from an initial HEX27 discretization
949  case PYRAMID5: // PYRAMID5's are created from an initial HEX27 discretization
950  case PYRAMID14:
951  case PRISM15:
952  case PRISM18:
953  {
954  for (unsigned int k=0; k<=(2*nz); k++)
955  for (unsigned int j=0; j<=(2*ny); j++)
956  for (unsigned int i=0; i<=(2*nx); i++)
957  {
958  if (gauss_lobatto_grid)
959  {
960  // The x,y locations of the point.
961  Real x=0., y=0., z=0.;
962 
963  // Shortcut quantities (do not depend on i,j)
964  const Real a = std::cos( libMesh::pi / static_cast<Real>(2*nx) );
965  const Real b = std::cos( libMesh::pi / static_cast<Real>(2*ny) );
966 
967  // Shortcut quantities (depend on i,j)
968  const Real c = std::cos( libMesh::pi*i / static_cast<Real>(2*nx) );
969  const Real d = std::cos( libMesh::pi*j / static_cast<Real>(2*ny) );
970 
971  // Additional shortcut quantities (for 3D)
972  const Real e = std::cos( libMesh::pi / static_cast<Real>(2*nz) );
973  const Real f = std::cos( libMesh::pi*k / static_cast<Real>(2*nz) );
974 
975  // If i is even, compute a normal Gauss-Lobatto point
976  if (i%2 == 0)
977  x = 0.5*(1.0 - c);
978 
979  // Otherwise, it is the average of the previous and next points
980  else
981  x = 0.5*(1.0 - a*c);
982 
983  // If j is even, compute a normal Gauss-Lobatto point
984  if (j%2 == 0)
985  y = 0.5*(1.0 - d);
986 
987  // Otherwise, it is the average of the previous and next points
988  else
989  y = 0.5*(1.0 - b*d);
990 
991  // If k is even, compute a normal Gauss-Lobatto point
992  if (k%2 == 0)
993  z = 0.5*(1.0 - f);
994 
995  // Otherwise, it is the average of the previous and next points
996  else
997  z = 0.5*(1.0 - e*f);
998 
999 
1000  mesh.add_point (Point(x,y,z), node_id++);
1001  }
1002 
1003  else
1004  mesh.add_point(Point(static_cast<Real>(i)/static_cast<Real>(2*nx),
1005  static_cast<Real>(j)/static_cast<Real>(2*ny),
1006  static_cast<Real>(k)/static_cast<Real>(2*nz)), node_id++);
1007  }
1008  break;
1009  }
1010 
1011 
1012  default:
1013  {
1014  libMesh::err << "ERROR: Unrecognized 3D element type." << std::endl;
1015  libmesh_error();
1016  }
1017  }
1018 
1019 
1020 
1021 
1022  // Build the elements.
1023  switch (type)
1024  {
1025  case INVALID_ELEM:
1026  case HEX8:
1027  {
1028  for (unsigned int k=0; k<nz; k++)
1029  for (unsigned int j=0; j<ny; j++)
1030  for (unsigned int i=0; i<nx; i++)
1031  {
1032  Elem* elem = mesh.add_elem(new Hex8);
1033 
1034  elem->set_node(0) = mesh.node_ptr(idx(type,nx,ny,i,j,k) );
1035  elem->set_node(1) = mesh.node_ptr(idx(type,nx,ny,i+1,j,k) );
1036  elem->set_node(2) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) );
1037  elem->set_node(3) = mesh.node_ptr(idx(type,nx,ny,i,j+1,k) );
1038  elem->set_node(4) = mesh.node_ptr(idx(type,nx,ny,i,j,k+1) );
1039  elem->set_node(5) = mesh.node_ptr(idx(type,nx,ny,i+1,j,k+1) );
1040  elem->set_node(6) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1));
1041  elem->set_node(7) = mesh.node_ptr(idx(type,nx,ny,i,j+1,k+1) );
1042 
1043  if (k == 0)
1044  mesh.boundary_info->add_side(elem, 0, 0);
1045 
1046  if (k == (nz-1))
1047  mesh.boundary_info->add_side(elem, 5, 5);
1048 
1049  if (j == 0)
1050  mesh.boundary_info->add_side(elem, 1, 1);
1051 
1052  if (j == (ny-1))
1053  mesh.boundary_info->add_side(elem, 3, 3);
1054 
1055  if (i == 0)
1056  mesh.boundary_info->add_side(elem, 4, 4);
1057 
1058  if (i == (nx-1))
1059  mesh.boundary_info->add_side(elem, 2, 2);
1060  }
1061  break;
1062  }
1063 
1064 
1065 
1066 
1067  case PRISM6:
1068  {
1069  for (unsigned int k=0; k<nz; k++)
1070  for (unsigned int j=0; j<ny; j++)
1071  for (unsigned int i=0; i<nx; i++)
1072  {
1073  // First Prism
1074  Elem* elem = NULL;
1075  elem = mesh.add_elem(new Prism6);
1076 
1077  elem->set_node(0) = mesh.node_ptr(idx(type,nx,ny,i,j,k) );
1078  elem->set_node(1) = mesh.node_ptr(idx(type,nx,ny,i+1,j,k) );
1079  elem->set_node(2) = mesh.node_ptr(idx(type,nx,ny,i,j+1,k) );
1080  elem->set_node(3) = mesh.node_ptr(idx(type,nx,ny,i,j,k+1) );
1081  elem->set_node(4) = mesh.node_ptr(idx(type,nx,ny,i+1,j,k+1) );
1082  elem->set_node(5) = mesh.node_ptr(idx(type,nx,ny,i,j+1,k+1) );
1083 
1084  // Add sides for first prism to boundary info object
1085  if (i==0)
1086  mesh.boundary_info->add_side(elem, 3, 4);
1087 
1088  if (j==0)
1089  mesh.boundary_info->add_side(elem, 1, 1);
1090 
1091  if (k==0)
1092  mesh.boundary_info->add_side(elem, 0, 0);
1093 
1094  if (k == (nz-1))
1095  mesh.boundary_info->add_side(elem, 4, 5);
1096 
1097  // Second Prism
1098  elem = mesh.add_elem(new Prism6);
1099 
1100  elem->set_node(0) = mesh.node_ptr(idx(type,nx,ny,i+1,j,k) );
1101  elem->set_node(1) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) );
1102  elem->set_node(2) = mesh.node_ptr(idx(type,nx,ny,i,j+1,k) );
1103  elem->set_node(3) = mesh.node_ptr(idx(type,nx,ny,i+1,j,k+1) );
1104  elem->set_node(4) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1));
1105  elem->set_node(5) = mesh.node_ptr(idx(type,nx,ny,i,j+1,k+1) );
1106 
1107  // Add sides for second prism to boundary info object
1108  if (i == (nx-1))
1109  mesh.boundary_info->add_side(elem, 1, 2);
1110 
1111  if (j == (ny-1))
1112  mesh.boundary_info->add_side(elem, 2, 3);
1113 
1114  if (k==0)
1115  mesh.boundary_info->add_side(elem, 0, 0);
1116 
1117  if (k == (nz-1))
1118  mesh.boundary_info->add_side(elem, 4, 5);
1119  }
1120  break;
1121  }
1122 
1123 
1124 
1125 
1126 
1127 
1128  case HEX20:
1129  case HEX27:
1130  case TET4: // TET4's are created from an initial HEX27 discretization
1131  case TET10: // TET10's are created from an initial HEX27 discretization
1132  case PYRAMID5: // PYRAMID5's are created from an initial HEX27 discretization
1133  case PYRAMID14:
1134  {
1135  for (unsigned int k=0; k<(2*nz); k += 2)
1136  for (unsigned int j=0; j<(2*ny); j += 2)
1137  for (unsigned int i=0; i<(2*nx); i += 2)
1138  {
1139  Elem* elem = (type == HEX20) ?
1140  mesh.add_elem(new Hex20) :
1141  mesh.add_elem(new Hex27);
1142 
1143  elem->set_node(0) = mesh.node_ptr(idx(type,nx,ny,i, j, k) );
1144  elem->set_node(1) = mesh.node_ptr(idx(type,nx,ny,i+2,j, k) );
1145  elem->set_node(2) = mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k) );
1146  elem->set_node(3) = mesh.node_ptr(idx(type,nx,ny,i, j+2,k) );
1147  elem->set_node(4) = mesh.node_ptr(idx(type,nx,ny,i, j, k+2));
1148  elem->set_node(5) = mesh.node_ptr(idx(type,nx,ny,i+2,j, k+2));
1149  elem->set_node(6) = mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k+2));
1150  elem->set_node(7) = mesh.node_ptr(idx(type,nx,ny,i, j+2,k+2));
1151  elem->set_node(8) = mesh.node_ptr(idx(type,nx,ny,i+1,j, k) );
1152  elem->set_node(9) = mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k) );
1153  elem->set_node(10) = mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k) );
1154  elem->set_node(11) = mesh.node_ptr(idx(type,nx,ny,i, j+1,k) );
1155  elem->set_node(12) = mesh.node_ptr(idx(type,nx,ny,i, j, k+1));
1156  elem->set_node(13) = mesh.node_ptr(idx(type,nx,ny,i+2,j, k+1));
1157  elem->set_node(14) = mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k+1));
1158  elem->set_node(15) = mesh.node_ptr(idx(type,nx,ny,i, j+2,k+1));
1159  elem->set_node(16) = mesh.node_ptr(idx(type,nx,ny,i+1,j, k+2));
1160  elem->set_node(17) = mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k+2));
1161  elem->set_node(18) = mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k+2));
1162  elem->set_node(19) = mesh.node_ptr(idx(type,nx,ny,i, j+1,k+2));
1163  if ((type == HEX27) || (type == TET4) || (type == TET10) || (type == PYRAMID5) || (type == PYRAMID14))
1164  {
1165  elem->set_node(20) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) );
1166  elem->set_node(21) = mesh.node_ptr(idx(type,nx,ny,i+1,j, k+1));
1167  elem->set_node(22) = mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k+1));
1168  elem->set_node(23) = mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k+1));
1169  elem->set_node(24) = mesh.node_ptr(idx(type,nx,ny,i, j+1,k+1));
1170  elem->set_node(25) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+2));
1171  elem->set_node(26) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1));
1172  }
1173 
1174 
1175  if (k == 0)
1176  mesh.boundary_info->add_side(elem, 0, 0);
1177 
1178  if (k == 2*(nz-1))
1179  mesh.boundary_info->add_side(elem, 5, 5);
1180 
1181  if (j == 0)
1182  mesh.boundary_info->add_side(elem, 1, 1);
1183 
1184  if (j == 2*(ny-1))
1185  mesh.boundary_info->add_side(elem, 3, 3);
1186 
1187  if (i == 0)
1188  mesh.boundary_info->add_side(elem, 4, 4);
1189 
1190  if (i == 2*(nx-1))
1191  mesh.boundary_info->add_side(elem, 2, 2);
1192  }
1193  break;
1194  }
1195 
1196 
1197 
1198 
1199  case PRISM15:
1200  case PRISM18:
1201  {
1202  for (unsigned int k=0; k<(2*nz); k += 2)
1203  for (unsigned int j=0; j<(2*ny); j += 2)
1204  for (unsigned int i=0; i<(2*nx); i += 2)
1205  {
1206  // First Prism
1207  Elem* elem = NULL;
1208  elem = ((type == PRISM15) ?
1209  mesh.add_elem(new Prism15) :
1210  mesh.add_elem(new Prism18));
1211 
1212  elem->set_node(0) = mesh.node_ptr(idx(type,nx,ny,i, j, k) );
1213  elem->set_node(1) = mesh.node_ptr(idx(type,nx,ny,i+2,j, k) );
1214  elem->set_node(2) = mesh.node_ptr(idx(type,nx,ny,i, j+2,k) );
1215  elem->set_node(3) = mesh.node_ptr(idx(type,nx,ny,i, j, k+2));
1216  elem->set_node(4) = mesh.node_ptr(idx(type,nx,ny,i+2,j, k+2));
1217  elem->set_node(5) = mesh.node_ptr(idx(type,nx,ny,i, j+2,k+2));
1218  elem->set_node(6) = mesh.node_ptr(idx(type,nx,ny,i+1,j, k) );
1219  elem->set_node(7) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) );
1220  elem->set_node(8) = mesh.node_ptr(idx(type,nx,ny,i, j+1,k) );
1221  elem->set_node(9) = mesh.node_ptr(idx(type,nx,ny,i, j, k+1));
1222  elem->set_node(10) = mesh.node_ptr(idx(type,nx,ny,i+2,j, k+1));
1223  elem->set_node(11) = mesh.node_ptr(idx(type,nx,ny,i, j+2,k+1));
1224  elem->set_node(12) = mesh.node_ptr(idx(type,nx,ny,i+1,j, k+2));
1225  elem->set_node(13) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+2));
1226  elem->set_node(14) = mesh.node_ptr(idx(type,nx,ny,i, j+1,k+2));
1227  if (type == PRISM18)
1228  {
1229  elem->set_node(15) = mesh.node_ptr(idx(type,nx,ny,i+1,j, k+1));
1230  elem->set_node(16) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1));
1231  elem->set_node(17) = mesh.node_ptr(idx(type,nx,ny,i, j+1,k+1));
1232  }
1233 
1234  // Add sides for first prism to boundary info object
1235  if (i==0)
1236  mesh.boundary_info->add_side(elem, 3, 4);
1237 
1238  if (j==0)
1239  mesh.boundary_info->add_side(elem, 1, 1);
1240 
1241  if (k==0)
1242  mesh.boundary_info->add_side(elem, 0, 0);
1243 
1244  if (k == 2*(nz-1))
1245  mesh.boundary_info->add_side(elem, 4, 5);
1246 
1247 
1248  // Second Prism
1249  elem = ((type == PRISM15) ?
1250  mesh.add_elem(new Prism15) :
1251  mesh.add_elem(new Prism18));
1252 
1253  elem->set_node(0) = mesh.node_ptr(idx(type,nx,ny,i+2,j,k) );
1254  elem->set_node(1) = mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k) );
1255  elem->set_node(2) = mesh.node_ptr(idx(type,nx,ny,i,j+2,k) );
1256  elem->set_node(3) = mesh.node_ptr(idx(type,nx,ny,i+2,j,k+2) );
1257  elem->set_node(4) = mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k+2) );
1258  elem->set_node(5) = mesh.node_ptr(idx(type,nx,ny,i,j+2,k+2) );
1259  elem->set_node(6) = mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k) );
1260  elem->set_node(7) = mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k) );
1261  elem->set_node(8) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) );
1262  elem->set_node(9) = mesh.node_ptr(idx(type,nx,ny,i+2,j,k+1) );
1263  elem->set_node(10) = mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k+1));
1264  elem->set_node(11) = mesh.node_ptr(idx(type,nx,ny,i,j+2,k+1) );
1265  elem->set_node(12) = mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k+2));
1266  elem->set_node(13) = mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k+2));
1267  elem->set_node(14) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+2));
1268  if (type == PRISM18)
1269  {
1270  elem->set_node(15) = mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k+1));
1271  elem->set_node(16) = mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k+1));
1272  elem->set_node(17) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1));
1273  }
1274 
1275  // Add sides for second prism to boundary info object
1276  if (i == 2*(nx-1))
1277  mesh.boundary_info->add_side(elem, 1, 2);
1278 
1279  if (j == 2*(ny-1))
1280  mesh.boundary_info->add_side(elem, 2, 3);
1281 
1282  if (k==0)
1283  mesh.boundary_info->add_side(elem, 0, 0);
1284 
1285  if (k == 2*(nz-1))
1286  mesh.boundary_info->add_side(elem, 4, 5);
1287 
1288  }
1289  break;
1290  }
1291 
1292 
1293 
1294 
1295 
1296  default:
1297  {
1298  libMesh::err << "ERROR: Unrecognized 3D element type." << std::endl;
1299  libmesh_error();
1300  }
1301  }
1302 
1303 
1304 
1305 
1306  //.......................................
1307  // Scale the nodal positions
1308  for (unsigned int p=0; p<mesh.n_nodes(); p++)
1309  {
1310  mesh.node(p)(0) = (mesh.node(p)(0))*(xmax-xmin) + xmin;
1311  mesh.node(p)(1) = (mesh.node(p)(1))*(ymax-ymin) + ymin;
1312  mesh.node(p)(2) = (mesh.node(p)(2))*(zmax-zmin) + zmin;
1313  }
1314 
1315 
1316 
1317 
1318  // Additional work for tets and pyramids: we take the existing
1319  // HEX27 discretization and split each element into 24
1320  // sub-tets or 6 sub-pyramids.
1321  //
1322  // 24 isn't the minimum-possible number of tets, but it
1323  // obviates any concerns about the edge orientations between
1324  // the various elements.
1325  if ((type == TET4) ||
1326  (type == TET10) ||
1327  (type == PYRAMID5) ||
1328  (type == PYRAMID14))
1329  {
1330  // Temporary storage for new elements. (24 tets per hex, 6 pyramids)
1331  std::vector<Elem*> new_elements;
1332 
1333  if ((type == TET4) || (type == TET10))
1334  new_elements.reserve(24*mesh.n_elem());
1335  else
1336  new_elements.reserve(6*mesh.n_elem());
1337 
1338  // Create tetrahedra or pyramids
1339  {
1340  MeshBase::element_iterator el = mesh.elements_begin();
1341  const MeshBase::element_iterator end_el = mesh.elements_end();
1342 
1343  for ( ; el != end_el; ++el)
1344  {
1345  // Get a pointer to the HEX27 element.
1346  Elem* base_hex = *el;
1347 
1348  // Get a pointer to the node located at the HEX27 centroid
1349  Node* apex_node = base_hex->get_node(26);
1350 
1351  for (unsigned int s=0; s<base_hex->n_sides(); ++s)
1352  {
1353  // Get the boundary ID for this side
1354  boundary_id_type b_id = mesh.boundary_info->boundary_id(*el, s);
1355 
1356  // Need to build the full-ordered side!
1357  AutoPtr<Elem> side = base_hex->build_side(s);
1358 
1359  if ((type == TET4) || (type == TET10))
1360  {
1361  // Build 4 sub-tets per side
1362  for (unsigned int sub_tet=0; sub_tet<4; ++sub_tet)
1363  {
1364  new_elements.push_back( new Tet4 );
1365  Elem* sub_elem = new_elements.back();
1366  sub_elem->set_node(0) = side->get_node(sub_tet);
1367  sub_elem->set_node(1) = side->get_node(8); // centroid of the face
1368  sub_elem->set_node(2) = side->get_node(sub_tet==3 ? 0 : sub_tet+1 ); // wrap-around
1369  sub_elem->set_node(3) = apex_node; // apex node always used!
1370 
1371  // If the original hex was a boundary hex, add the new sub_tet's side
1372  // 0 with the same b_id. Note: the tets are all aligned so that their
1373  // side 0 is on the boundary.
1374  if (b_id != BoundaryInfo::invalid_id)
1375  mesh.boundary_info->add_side(sub_elem, 0, b_id);
1376  }
1377  } // end if ((type == TET4) || (type == TET10))
1378 
1379  else // type==PYRAMID5 || type==PYRAMID14
1380  {
1381  // Build 1 sub-pyramid per side.
1382  new_elements.push_back(new Pyramid5);
1383  Elem* sub_elem = new_elements.back();
1384 
1385  // Set the base. Note that since the apex is *inside* the base_hex,
1386  // and the pyramid uses a counter-clockwise base numbering, we need to
1387  // reverse the [1] and [3] node indices.
1388  sub_elem->set_node(0) = side->get_node(0);
1389  sub_elem->set_node(1) = side->get_node(3);
1390  sub_elem->set_node(2) = side->get_node(2);
1391  sub_elem->set_node(3) = side->get_node(1);
1392 
1393  // Set the apex
1394  sub_elem->set_node(4) = apex_node;
1395 
1396  // If the original hex was a boundary hex, add the new sub_pyr's side
1397  // 4 (the square base) with the same b_id.
1398  if (b_id != BoundaryInfo::invalid_id)
1399  mesh.boundary_info->add_side(sub_elem, 4, b_id);
1400  } // end else type==PYRAMID5 || type==PYRAMID14
1401  }
1402  }
1403  }
1404 
1405 
1406  // Delete the original HEX27 elements from the mesh, and the boundary info structure.
1407  {
1408  MeshBase::element_iterator el = mesh.elements_begin();
1409  const MeshBase::element_iterator end_el = mesh.elements_end();
1410 
1411  for ( ; el != end_el; ++el)
1412  {
1413  mesh.boundary_info->remove(*el); // Safe even if *el has no boundary info.
1414  mesh.delete_elem(*el);
1415  }
1416  }
1417 
1418  // Add the new elements
1419  for (unsigned int i=0; i<new_elements.size(); ++i)
1420  mesh.add_elem(new_elements[i]);
1421 
1422  } // end if (type == TET4,TET10,PYRAMID5,PYRAMID14
1423 
1424 
1425  // Use all_second_order to convert the TET4's to TET10's or PYRAMID5's to PYRAMID14's
1426  if ((type == TET10) || (type == PYRAMID14))
1427  {
1428  mesh.all_second_order();
1429  }
1430 
1431  // Add sideset names to boundary info (Z axis out of the screen)
1432  mesh.boundary_info->sideset_name(0) = "back";
1433  mesh.boundary_info->sideset_name(1) = "bottom";
1434  mesh.boundary_info->sideset_name(2) = "right";
1435  mesh.boundary_info->sideset_name(3) = "top";
1436  mesh.boundary_info->sideset_name(4) = "left";
1437  mesh.boundary_info->sideset_name(5) = "front";
1438 
1439  // Add nodeset names to boundary info
1440  mesh.boundary_info->nodeset_name(0) = "back";
1441  mesh.boundary_info->nodeset_name(1) = "bottom";
1442  mesh.boundary_info->nodeset_name(2) = "right";
1443  mesh.boundary_info->nodeset_name(3) = "top";
1444  mesh.boundary_info->nodeset_name(4) = "left";
1445  mesh.boundary_info->nodeset_name(5) = "front";
1446 
1447  break;
1448  } // end case dim==3
1449 
1450  default:
1451  {
1452  libmesh_error();
1453  }
1454  }
1455 
1456  STOP_LOG("build_cube()", "MeshTools::Generation");
1457 
1458 
1459 
1460  // Done building the mesh. Now prepare it for use.
1461  mesh.prepare_for_use (/*skip_renumber =*/ false);
1462 }
void libMesh::MeshTools::Generation::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 = NULL 
)

Meshes a rectangular (2D) region (with or without holes) with a Delaunay triangulation. This function internally calls the triangle library written by J.R. Shewchuk.

Definition at line 2212 of file mesh_generation.C.

References libMesh::MeshBase::add_point(), libMesh::TriangleInterface::attach_hole_list(), bc_id, libMesh::MeshBase::boundary_info, libMesh::Elem::build_side(), libMesh::MeshBase::clear(), libMesh::TriangleInterface::desired_area(), libMesh::TriangleInterface::elem_type(), libMesh::MeshBase::elements_begin(), libMesh::MeshBase::elements_end(), libMesh::DofObject::id(), libMesh::MeshBase::n_nodes(), libMesh::Elem::n_sides(), libMesh::Elem::neighbor(), libMesh::TriangleInterface::PSLG, libMesh::Real, libMesh::MeshBase::set_mesh_dimension(), side, libMesh::TOLERANCE, libMesh::TriangleInterface::triangulate(), and libMesh::TriangleInterface::triangulation_type().

2219 {
2220  // Check for reasonable size
2221  libmesh_assert_greater_equal (nx, 1); // need at least 1 element in x-direction
2222  libmesh_assert_greater_equal (ny, 1); // need at least 1 element in y-direction
2223  libmesh_assert_less (xmin, xmax);
2224  libmesh_assert_less (ymin, ymax);
2225 
2226  // Clear out any data which may have been in the Mesh
2227  mesh.clear();
2228 
2229  // Make sure the new Mesh will be 2D
2230  mesh.set_mesh_dimension(2);
2231 
2232  // The x and y spacing between boundary points
2233  const Real delta_x = (xmax-xmin) / static_cast<Real>(nx);
2234  const Real delta_y = (ymax-ymin) / static_cast<Real>(ny);
2235 
2236  // Bottom
2237  for (unsigned int p=0; p<=nx; ++p)
2238  mesh.add_point(Point(xmin + p*delta_x, ymin));
2239 
2240  // Right side
2241  for (unsigned int p=1; p<ny; ++p)
2242  mesh.add_point(Point(xmax, ymin + p*delta_y));
2243 
2244  // Top
2245  for (unsigned int p=0; p<=nx; ++p)
2246  mesh.add_point(Point(xmax - p*delta_x, ymax));
2247 
2248  // Left side
2249  for (unsigned int p=1; p<ny; ++p)
2250  mesh.add_point(Point(xmin, ymax - p*delta_y));
2251 
2252  // Be sure we added as many points as we thought we did
2253  libmesh_assert_equal_to (mesh.n_nodes(), 2*(nx+ny));
2254 
2255  // Construct the Triangle Interface object
2256  TriangleInterface t(mesh);
2257 
2258  // Set custom variables for the triangulation
2259  t.desired_area() = 0.5 * (xmax-xmin)*(ymax-ymin) / static_cast<Real>(nx*ny);
2260  t.triangulation_type() = TriangleInterface::PSLG;
2261  t.elem_type() = type;
2262 
2263  if (holes != NULL)
2264  t.attach_hole_list(holes);
2265 
2266  // Triangulate!
2267  t.triangulate();
2268 
2269  // The mesh is now generated, but we still need to mark the boundaries
2270  // to be consistent with the other build_square routines. Note that all
2271  // hole boundary elements get the same ID, 4.
2272  MeshBase::element_iterator el = mesh.elements_begin();
2273  const MeshBase::element_iterator end_el = mesh.elements_end();
2274  for ( ; el != end_el; ++el)
2275  {
2276  const Elem* elem = *el;
2277 
2278  for (unsigned int s=0; s<elem->n_sides(); s++)
2279  if (elem->neighbor(s) == NULL)
2280  {
2281  AutoPtr<Elem> side (elem->build_side(s));
2282 
2283  // Check the location of the side's midpoint. Since
2284  // the square has straight sides, the midpoint is not
2285  // on the corner and thus it is uniquely on one of the
2286  // sides.
2287  Point side_midpoint= 0.5f*( (*side->get_node(0)) + (*side->get_node(1)) );
2288 
2289  // The boundary ids are set following the same convention as Quad4 sides
2290  // bottom = 0
2291  // right = 1
2292  // top = 2
2293  // left = 3
2294  // hole = 4
2296 
2297  // bottom
2298  if (std::fabs(side_midpoint(1) - ymin) < TOLERANCE)
2299  bc_id=0;
2300 
2301  // right
2302  else if (std::fabs(side_midpoint(0) - xmax) < TOLERANCE)
2303  bc_id=1;
2304 
2305  // top
2306  else if (std::fabs(side_midpoint(1) - ymax) < TOLERANCE)
2307  bc_id=2;
2308 
2309  // left
2310  else if (std::fabs(side_midpoint(0) - xmin) < TOLERANCE)
2311  bc_id=3;
2312 
2313  // If the point is not on any of the external boundaries, it
2314  // is on one of the holes....
2315 
2316  // Finally, add this element's information to the boundary info object.
2317  mesh.boundary_info->add_side(elem->id(), s, bc_id);
2318  }
2319  }
2320 
2321 } // end build_delaunay_square
void libMesh::MeshTools::Generation::build_extrusion ( UnstructuredMesh &  mesh,
const MeshBase &  cross_section,
const unsigned int  nz,
RealVectorValue  extrusion_vector 
)

Meshes the tensor product of a 1D and a 1D-or-2D domain.

Definition at line 1992 of file mesh_generation.C.

References libMesh::MeshBase::add_elem(), libMesh::MeshBase::add_point(), libMesh::MeshBase::boundary_info, libMesh::MeshBase::delete_remote_elements(), libMeshEnums::EDGE2, libMeshEnums::EDGE3, libMesh::MeshBase::elements_begin(), libMesh::MeshBase::elements_end(), end, libMesh::Elem::get_node(), libMesh::DofObject::id(), libMesh::MeshBase::is_serial(), libMesh::libmesh_assert(), libMesh::MeshBase::n_elem(), libMesh::MeshBase::n_nodes(), libMesh::Elem::n_sides(), libMesh::MeshBase::node_ptr(), libMesh::MeshBase::nodes_begin(), libMesh::MeshBase::nodes_end(), libMesh::Elem::parent(), libMesh::MeshBase::prepare_for_use(), libMesh::DofObject::processor_id(), libMeshEnums::QUAD4, libMeshEnums::QUAD9, libMesh::MeshBase::reserve_elem(), libMesh::MeshBase::reserve_nodes(), libMeshEnums::SECOND, libMesh::DofObject::set_id(), libMesh::Elem::set_node(), libMesh::START_LOG(), libMesh::STOP_LOG(), libMesh::Elem::subdomain_id(), libMeshEnums::TRI3, libMeshEnums::TRI6, and libMesh::Elem::type().

1996 {
1997  if (!cross_section.n_elem())
1998  return;
1999 
2000  START_LOG("build_extrusion()", "MeshTools::Generation");
2001 
2002  dof_id_type orig_elem = cross_section.n_elem();
2003  dof_id_type orig_nodes = cross_section.n_nodes();
2004 
2005  unsigned int order = 1;
2006 
2007  // If cross_section is distributed, so is its extrusion
2008  if (!cross_section.is_serial())
2009  mesh.delete_remote_elements();
2010 
2011  // We know a priori how many elements we'll need
2012  mesh.reserve_elem(nz*orig_elem);
2013 
2014  // For straightforward meshes we need one or two additional layers per
2015  // element.
2016  if ((*cross_section.elements_begin())->default_order() == SECOND)
2017  order = 2;
2018 
2019  mesh.reserve_nodes((order*nz+1)*orig_nodes);
2020 
2021  MeshBase::const_node_iterator nd = cross_section.nodes_begin();
2022  const MeshBase::const_node_iterator nend = cross_section.nodes_end();
2023  for (; nd!=nend; ++nd)
2024  {
2025  const Node* node = *nd;
2026 
2027  for (unsigned int k=0; k != order*nz+1; ++k)
2028  {
2029  Node *new_node =
2030  mesh.add_point(*node +
2031  (extrusion_vector * k / nz / order),
2032  node->id() + (k * orig_nodes),
2033  node->processor_id());
2034 
2035  const std::vector<boundary_id_type> ids_to_copy =
2036  cross_section.boundary_info->boundary_ids(node);
2037 
2038  mesh.boundary_info->add_node(new_node, ids_to_copy);
2039  }
2040  }
2041 
2042  const std::set<boundary_id_type> &side_ids =
2043  cross_section.boundary_info->get_side_boundary_ids();
2044  const boundary_id_type next_side_id = side_ids.empty() ?
2045  0 : *side_ids.rbegin() + 1;
2046 
2047  MeshBase::const_element_iterator el = cross_section.elements_begin();
2048  const MeshBase::const_element_iterator end = cross_section.elements_end();
2049  for (; el!=end; ++el)
2050  {
2051  const Elem* elem = *el;
2052  const ElemType etype = elem->type();
2053 
2054  // build_extrusion currently only works on coarse meshes
2055  libmesh_assert (!elem->parent());
2056 
2057  // We need a map from low-D to high-D sides for boundary id
2058  // setting
2059  std::vector<unsigned char> sidemap(4);
2060 
2061  for (unsigned int k=0; k != nz; ++k)
2062  {
2063  Elem *new_elem;
2064  switch (etype)
2065  {
2066  case EDGE2:
2067  {
2068  new_elem = new Quad4;
2069  new_elem->set_node(0) = mesh.node_ptr(elem->get_node(0)->id() + (k * orig_nodes));
2070  new_elem->set_node(1) = mesh.node_ptr(elem->get_node(1)->id() + (k * orig_nodes));
2071  new_elem->set_node(2) = mesh.node_ptr(elem->get_node(1)->id() + ((k+1) * orig_nodes));
2072  new_elem->set_node(3) = mesh.node_ptr(elem->get_node(0)->id() + ((k+1) * orig_nodes));
2073  break;
2074  }
2075  case EDGE3:
2076  {
2077  new_elem = new Quad9;
2078  new_elem->set_node(0) = mesh.node_ptr(elem->get_node(0)->id() + (2*k * orig_nodes));
2079  new_elem->set_node(1) = mesh.node_ptr(elem->get_node(1)->id() + (2*k * orig_nodes));
2080  new_elem->set_node(2) = mesh.node_ptr(elem->get_node(1)->id() + ((2*k+2) * orig_nodes));
2081  new_elem->set_node(3) = mesh.node_ptr(elem->get_node(0)->id() + ((2*k+2) * orig_nodes));
2082  new_elem->set_node(4) = mesh.node_ptr(elem->get_node(2)->id() + (2*k * orig_nodes));
2083  new_elem->set_node(5) = mesh.node_ptr(elem->get_node(1)->id() + ((2*k+1) * orig_nodes));
2084  new_elem->set_node(6) = mesh.node_ptr(elem->get_node(2)->id() + ((2*k+2) * orig_nodes));
2085  new_elem->set_node(7) = mesh.node_ptr(elem->get_node(0)->id() + ((2*k+1) * orig_nodes));
2086  new_elem->set_node(8) = mesh.node_ptr(elem->get_node(2)->id() + ((2*k+1) * orig_nodes));
2087  break;
2088  }
2089  case TRI3:
2090  {
2091  new_elem = new Prism6;
2092  new_elem->set_node(0) = mesh.node_ptr(elem->get_node(0)->id() + (k * orig_nodes));
2093  new_elem->set_node(1) = mesh.node_ptr(elem->get_node(1)->id() + (k * orig_nodes));
2094  new_elem->set_node(2) = mesh.node_ptr(elem->get_node(2)->id() + (k * orig_nodes));
2095  new_elem->set_node(3) = mesh.node_ptr(elem->get_node(0)->id() + ((k+1) * orig_nodes));
2096  new_elem->set_node(4) = mesh.node_ptr(elem->get_node(1)->id() + ((k+1) * orig_nodes));
2097  new_elem->set_node(5) = mesh.node_ptr(elem->get_node(2)->id() + ((k+1) * orig_nodes));
2098  break;
2099  }
2100  case TRI6:
2101  {
2102  new_elem = new Prism18;
2103  new_elem->set_node(0) = mesh.node_ptr(elem->get_node(0)->id() + (2*k * orig_nodes));
2104  new_elem->set_node(1) = mesh.node_ptr(elem->get_node(1)->id() + (2*k * orig_nodes));
2105  new_elem->set_node(2) = mesh.node_ptr(elem->get_node(2)->id() + (2*k * orig_nodes));
2106  new_elem->set_node(3) = mesh.node_ptr(elem->get_node(0)->id() + ((2*k+2) * orig_nodes));
2107  new_elem->set_node(4) = mesh.node_ptr(elem->get_node(1)->id() + ((2*k+2) * orig_nodes));
2108  new_elem->set_node(5) = mesh.node_ptr(elem->get_node(2)->id() + ((2*k+2) * orig_nodes));
2109  new_elem->set_node(6) = mesh.node_ptr(elem->get_node(3)->id() + (2*k * orig_nodes));
2110  new_elem->set_node(7) = mesh.node_ptr(elem->get_node(4)->id() + (2*k * orig_nodes));
2111  new_elem->set_node(8) = mesh.node_ptr(elem->get_node(5)->id() + (2*k * orig_nodes));
2112  new_elem->set_node(9) = mesh.node_ptr(elem->get_node(0)->id() + ((2*k+1) * orig_nodes));
2113  new_elem->set_node(10) = mesh.node_ptr(elem->get_node(1)->id() + ((2*k+1) * orig_nodes));
2114  new_elem->set_node(11) = mesh.node_ptr(elem->get_node(2)->id() + ((2*k+1) * orig_nodes));
2115  new_elem->set_node(12) = mesh.node_ptr(elem->get_node(3)->id() + ((2*k+2) * orig_nodes));
2116  new_elem->set_node(13) = mesh.node_ptr(elem->get_node(4)->id() + ((2*k+2) * orig_nodes));
2117  new_elem->set_node(14) = mesh.node_ptr(elem->get_node(5)->id() + ((2*k+2) * orig_nodes));
2118  new_elem->set_node(15) = mesh.node_ptr(elem->get_node(3)->id() + ((2*k+1) * orig_nodes));
2119  new_elem->set_node(16) = mesh.node_ptr(elem->get_node(4)->id() + ((2*k+1) * orig_nodes));
2120  new_elem->set_node(17) = mesh.node_ptr(elem->get_node(5)->id() + ((2*k+1) * orig_nodes));
2121  break;
2122  }
2123  case QUAD4:
2124  {
2125  new_elem = new Hex8;
2126  new_elem->set_node(0) = mesh.node_ptr(elem->get_node(0)->id() + (k * orig_nodes));
2127  new_elem->set_node(1) = mesh.node_ptr(elem->get_node(1)->id() + (k * orig_nodes));
2128  new_elem->set_node(2) = mesh.node_ptr(elem->get_node(2)->id() + (k * orig_nodes));
2129  new_elem->set_node(3) = mesh.node_ptr(elem->get_node(3)->id() + (k * orig_nodes));
2130  new_elem->set_node(4) = mesh.node_ptr(elem->get_node(0)->id() + ((k+1) * orig_nodes));
2131  new_elem->set_node(5) = mesh.node_ptr(elem->get_node(1)->id() + ((k+1) * orig_nodes));
2132  new_elem->set_node(6) = mesh.node_ptr(elem->get_node(2)->id() + ((k+1) * orig_nodes));
2133  new_elem->set_node(7) = mesh.node_ptr(elem->get_node(3)->id() + ((k+1) * orig_nodes));
2134  break;
2135  }
2136  case QUAD9:
2137  {
2138  new_elem = new Hex27;
2139  new_elem->set_node(0) = mesh.node_ptr(elem->get_node(0)->id() + (2*k * orig_nodes));
2140  new_elem->set_node(1) = mesh.node_ptr(elem->get_node(1)->id() + (2*k * orig_nodes));
2141  new_elem->set_node(2) = mesh.node_ptr(elem->get_node(2)->id() + (2*k * orig_nodes));
2142  new_elem->set_node(3) = mesh.node_ptr(elem->get_node(3)->id() + (2*k * orig_nodes));
2143  new_elem->set_node(4) = mesh.node_ptr(elem->get_node(0)->id() + ((2*k+2) * orig_nodes));
2144  new_elem->set_node(5) = mesh.node_ptr(elem->get_node(1)->id() + ((2*k+2) * orig_nodes));
2145  new_elem->set_node(6) = mesh.node_ptr(elem->get_node(2)->id() + ((2*k+2) * orig_nodes));
2146  new_elem->set_node(7) = mesh.node_ptr(elem->get_node(3)->id() + ((2*k+2) * orig_nodes));
2147  new_elem->set_node(8) = mesh.node_ptr(elem->get_node(4)->id() + (2*k * orig_nodes));
2148  new_elem->set_node(9) = mesh.node_ptr(elem->get_node(5)->id() + (2*k * orig_nodes));
2149  new_elem->set_node(10) = mesh.node_ptr(elem->get_node(6)->id() + (2*k * orig_nodes));
2150  new_elem->set_node(11) = mesh.node_ptr(elem->get_node(7)->id() + (2*k * orig_nodes));
2151  new_elem->set_node(12) = mesh.node_ptr(elem->get_node(0)->id() + ((2*k+1) * orig_nodes));
2152  new_elem->set_node(13) = mesh.node_ptr(elem->get_node(1)->id() + ((2*k+1) * orig_nodes));
2153  new_elem->set_node(14) = mesh.node_ptr(elem->get_node(2)->id() + ((2*k+1) * orig_nodes));
2154  new_elem->set_node(15) = mesh.node_ptr(elem->get_node(3)->id() + ((2*k+1) * orig_nodes));
2155  new_elem->set_node(16) = mesh.node_ptr(elem->get_node(4)->id() + ((2*k+2) * orig_nodes));
2156  new_elem->set_node(17) = mesh.node_ptr(elem->get_node(5)->id() + ((2*k+2) * orig_nodes));
2157  new_elem->set_node(18) = mesh.node_ptr(elem->get_node(6)->id() + ((2*k+2) * orig_nodes));
2158  new_elem->set_node(19) = mesh.node_ptr(elem->get_node(7)->id() + ((2*k+2) * orig_nodes));
2159  new_elem->set_node(20) = mesh.node_ptr(elem->get_node(8)->id() + (2*k * orig_nodes));
2160  new_elem->set_node(21) = mesh.node_ptr(elem->get_node(4)->id() + ((2*k+1) * orig_nodes));
2161  new_elem->set_node(22) = mesh.node_ptr(elem->get_node(5)->id() + ((2*k+1) * orig_nodes));
2162  new_elem->set_node(23) = mesh.node_ptr(elem->get_node(6)->id() + ((2*k+1) * orig_nodes));
2163  new_elem->set_node(24) = mesh.node_ptr(elem->get_node(7)->id() + ((2*k+1) * orig_nodes));
2164  new_elem->set_node(25) = mesh.node_ptr(elem->get_node(8)->id() + ((2*k+2) * orig_nodes));
2165  new_elem->set_node(26) = mesh.node_ptr(elem->get_node(8)->id() + ((2*k+1) * orig_nodes));
2166  break;
2167  }
2168  default:
2169  {
2170  libmesh_not_implemented();
2171  break;
2172  }
2173  }
2174 
2175  new_elem->set_id(elem->id() + (k * orig_elem));
2176  new_elem->processor_id() = elem->processor_id();
2177 
2178  // maintain the subdomain_id
2179  new_elem->subdomain_id() = elem->subdomain_id();
2180 
2181  new_elem = mesh.add_elem(new_elem);
2182 
2183  // Copy any old boundary ids on all sides
2184  for (unsigned int s = 0; s != elem->n_sides(); ++s)
2185  {
2186  const std::vector<boundary_id_type> ids_to_copy =
2187  cross_section.boundary_info->boundary_ids(elem, s);
2188 
2189  mesh.boundary_info->add_side(new_elem, s+1, ids_to_copy);
2190  }
2191 
2192  // Give new boundary ids to bottom and top
2193  if (k == 0)
2194  mesh.boundary_info->add_side(new_elem, 0, next_side_id);
2195  if (k == nz-1)
2196  mesh.boundary_info->add_side(new_elem, elem->n_sides()+1, next_side_id+1);
2197  }
2198  }
2199 
2200  STOP_LOG("build_extrusion()", "MeshTools::Generation");
2201 
2202  // Done building the mesh. Now prepare it for use.
2203  mesh.prepare_for_use(/*skip_renumber =*/ false);
2204 }
void libMesh::MeshTools::Generation::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

Boundary ids are set to be equal to the side indexing on a master edge

Definition at line 1484 of file mesh_generation.C.

References build_cube().

1489 {
1490  // This method only makes sense in 1D!
1491  // But we now just turn a non-1D mesh into a 1D mesh
1492  //libmesh_assert_equal_to (mesh.mesh_dimension(), 1);
1493 
1494  build_cube(mesh,
1495  nx, 0, 0,
1496  xmin, xmax,
1497  0., 0.,
1498  0., 0.,
1499  type,
1500  gauss_lobatto_grid);
1501 }
void libMesh::MeshTools::Generation::build_point ( UnstructuredMesh &  mesh,
const ElemType  type = INVALID_ELEM,
const bool  gauss_lobatto_grid = false 
)

A specialized build_cube() for 0D meshes. The resulting mesh is a single NodeElem suitable for ODE tests

Definition at line 1466 of file mesh_generation.C.

References build_cube().

1469 {
1470  // This method only makes sense in 0D!
1471  // But we now just turn a non-0D mesh into a 0D mesh
1472  //libmesh_assert_equal_to (mesh.mesh_dimension(), 1);
1473 
1474  build_cube(mesh,
1475  0, 0, 0,
1476  0., 0.,
1477  0., 0.,
1478  0., 0.,
1479  type,
1480  gauss_lobatto_grid);
1481 }
void libMesh::MeshTools::Generation::build_sphere ( UnstructuredMesh &  mesh,
const Real  rad = 1,
const unsigned int  nr = 2,
const ElemType  type = INVALID_ELEM,
const unsigned int  n_smooth = 2,
const bool  flat = true 
)

Meshes a spherical or mapped-spherical domain.

Definition at line 1536 of file mesh_generation.C.

References libMesh::out.

1542 {
1543  libMesh::out << "Building a circle/sphere only works with AMR." << std::endl;
1544  libmesh_error();
1545 }
void libMesh::MeshTools::Generation::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.

Boundary ids are set to be equal to the side indexing on a master quad

Definition at line 1505 of file mesh_generation.C.

References build_cube().

1512 {
1513  // This method only makes sense in 2D!
1514  // But we now just turn a non-2D mesh into a 2D mesh
1515  //libmesh_assert_equal_to (mesh.mesh_dimension(), 2);
1516 
1517  // Call the build_cube() member to actually do the work for us.
1518  build_cube (mesh,
1519  nx, ny, 0,
1520  xmin, xmax,
1521  ymin, ymax,
1522  0., 0.,
1523  type,
1524  gauss_lobatto_grid);
1525 }

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:33 UTC

Hosted By:
SourceForge.net Logo