libMesh::MeshTools::Generation Namespace Reference
Namespaces | |
| namespace | 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.
- 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 
(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 141 of file mesh_generation.C.
References libMesh::MeshBase::add_elem(), libMesh::MeshBase::add_point(), libMesh::UnstructuredMesh::all_second_order(), libMesh::Elem::build_side(), libMesh::MeshBase::clear(), libMesh::MeshBase::delete_elem(), libMesh::EDGE2, libMesh::EDGE3, libMesh::EDGE4, libMesh::MeshBase::elements_begin(), libMesh::MeshBase::elements_end(), libMesh::MeshBase::get_boundary_info(), libMesh::Elem::get_node(), libMesh::HEX20, libMesh::HEX27, libMesh::HEX8, libMesh::MeshTools::Generation::Private::idx(), libMesh::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(), libMesh::NODEELEM, libMesh::pi, libMesh::MeshBase::prepare_for_use(), libMesh::PRISM15, libMesh::PRISM18, libMesh::PRISM6, libMesh::PYRAMID13, libMesh::PYRAMID14, libMesh::PYRAMID5, libMesh::QUAD4, libMesh::QUAD8, libMesh::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::TET10, libMesh::TET4, libMesh::TRI3, libMesh::TRI6, and libMesh::x.
Referenced by build_line(), build_point(), and build_square().
{
START_LOG("build_cube()", "MeshTools::Generation");
// Declare that we are using the indexing utility routine
// in the "Private" part of our current namespace. If this doesn't
// work in GCC 2.95.3 we can either remove it or stop supporting
// 2.95.3 altogether.
// Changing this to import the whole namespace... just importing idx
// causes an internal compiler error for Intel Compiler 11.0 on Linux
// in debug mode.
using namespace MeshTools::Generation::Private;
// Clear the mesh and start from scratch
mesh.clear();
BoundaryInfo& boundary_info = mesh.get_boundary_info();
if (nz != 0)
mesh.set_mesh_dimension(3);
else if (ny != 0)
mesh.set_mesh_dimension(2);
else if (nx != 0)
mesh.set_mesh_dimension(1);
else
mesh.set_mesh_dimension(0);
switch (mesh.mesh_dimension())
{
//---------------------------------------------------------------------
// Build a 0D point
case 0:
{
libmesh_assert_equal_to (nx, 0);
libmesh_assert_equal_to (ny, 0);
libmesh_assert_equal_to (nz, 0);
libmesh_assert (type == INVALID_ELEM || type == NODEELEM);
// Build one nodal element for the mesh
mesh.add_point (Point(0, 0, 0), 0);
Elem* elem = mesh.add_elem (new NodeElem);
elem->set_node(0) = mesh.node_ptr(0);
break;
}
//---------------------------------------------------------------------
// Build a 1D line
case 1:
{
libmesh_assert_not_equal_to (nx, 0);
libmesh_assert_equal_to (ny, 0);
libmesh_assert_equal_to (nz, 0);
libmesh_assert_less (xmin, xmax);
// Reserve elements
switch (type)
{
case INVALID_ELEM:
case EDGE2:
case EDGE3:
case EDGE4:
{
mesh.reserve_elem (nx);
break;
}
default:
libmesh_error_msg("ERROR: Unrecognized 1D element type.");
}
// Reserve nodes
switch (type)
{
case INVALID_ELEM:
case EDGE2:
{
mesh.reserve_nodes(nx+1);
break;
}
case EDGE3:
{
mesh.reserve_nodes(2*nx+1);
break;
}
case EDGE4:
{
mesh.reserve_nodes(3*nx+1);
break;
}
default:
libmesh_error_msg("ERROR: Unrecognized 1D element type.");
}
// Build the nodes, depends on whether we're using linears,
// quadratics or cubics and whether using uniform grid or Gauss-Lobatto
unsigned int node_id = 0;
switch(type)
{
case INVALID_ELEM:
case EDGE2:
{
for (unsigned int i=0; i<=nx; i++)
{
if (gauss_lobatto_grid)
mesh.add_point (Point(0.5*(std::cos(libMesh::pi*static_cast<Real>(nx-i)/static_cast<Real>(nx))+1.0),
0,
0), node_id++);
else
mesh.add_point (Point(static_cast<Real>(i)/static_cast<Real>(nx),
0,
0), node_id++);
}
break;
}
case EDGE3:
{
for (unsigned int i=0; i<=2*nx; i++)
{
if (gauss_lobatto_grid)
{
// The x location of the point.
Real x=0.;
// Shortcut quantities (do not depend on i)
const Real c = std::cos( libMesh::pi*i / static_cast<Real>(2*nx) );
// If i is even, compute a normal Gauss-Lobatto point
if (i%2 == 0)
x = 0.5*(1.0 - c);
// Otherwise, it is the average of the previous and next points
else
{
Real cmin = std::cos( libMesh::pi*(i-1) / static_cast<Real>(2*nx) );
Real cmax = std::cos( libMesh::pi*(i+1) / static_cast<Real>(2*nx) );
Real gl_xmin = 0.5*(1.0 - cmin);
Real gl_xmax = 0.5*(1.0 - cmax);
x = 0.5*(gl_xmin + gl_xmax);
}
mesh.add_point (Point(x,0.,0.), node_id++);
}
else
mesh.add_point (Point(static_cast<Real>(i)/static_cast<Real>(2*nx),
0,
0), node_id++);
}
break;
}
case EDGE4:
{
for (unsigned int i=0; i<=3*nx; i++)
{
if (gauss_lobatto_grid)
{
// The x location of the point
Real x=0.;
// Shortcut quantities
const Real c = std::cos( libMesh::pi*i / static_cast<Real>(3*nx) );
// If i is multiple of 3, compute a normal Gauss-Lobatto point
if (i%3 == 0)
x = 0.5*(1.0 - c);
// Otherwise, distribute points evenly within the element
else
{
if(i%3 == 1)
{
Real cmin = std::cos( libMesh::pi*(i-1) / static_cast<Real>(3*nx) );
Real cmax = std::cos( libMesh::pi*(i+2) / static_cast<Real>(3*nx) );
Real gl_xmin = 0.5*(1.0 - cmin);
Real gl_xmax = 0.5*(1.0 - cmax);
x = (2.*gl_xmin + gl_xmax)/3.;
}
else
if(i%3 == 2)
{
Real cmin = std::cos( libMesh::pi*(i-2) / static_cast<Real>(3*nx) );
Real cmax = std::cos( libMesh::pi*(i+1) / static_cast<Real>(3*nx) );
Real gl_xmin = 0.5*(1.0 - cmin);
Real gl_xmax = 0.5*(1.0 - cmax);
x = (gl_xmin + 2.*gl_xmax)/3.;
}
}
mesh.add_point (Point(x,0.,0.), node_id++);
}
else
mesh.add_point (Point(static_cast<Real>(i)/static_cast<Real>(3*nx),
0,
0), node_id++);
}
break;
}
default:
libmesh_error_msg("ERROR: Unrecognized 1D element type.");
}
// Build the elements of the mesh
switch(type)
{
case INVALID_ELEM:
case EDGE2:
{
for (unsigned int i=0; i<nx; i++)
{
Elem* elem = mesh.add_elem (new Edge2);
elem->set_node(0) = mesh.node_ptr(i);
elem->set_node(1) = mesh.node_ptr(i+1);
if (i == 0)
boundary_info.add_side(elem, 0, 0);
if (i == (nx-1))
boundary_info.add_side(elem, 1, 1);
}
break;
}
case EDGE3:
{
for (unsigned int i=0; i<nx; i++)
{
Elem* elem = mesh.add_elem (new Edge3);
elem->set_node(0) = mesh.node_ptr(2*i);
elem->set_node(2) = mesh.node_ptr(2*i+1);
elem->set_node(1) = mesh.node_ptr(2*i+2);
if (i == 0)
boundary_info.add_side(elem, 0, 0);
if (i == (nx-1))
boundary_info.add_side(elem, 1, 1);
}
break;
}
case EDGE4:
{
for (unsigned int i=0; i<nx; i++)
{
Elem* elem = mesh.add_elem (new Edge4);
elem->set_node(0) = mesh.node_ptr(3*i);
elem->set_node(2) = mesh.node_ptr(3*i+1);
elem->set_node(3) = mesh.node_ptr(3*i+2);
elem->set_node(1) = mesh.node_ptr(3*i+3);
if (i == 0)
boundary_info.add_side(elem, 0, 0);
if (i == (nx-1))
boundary_info.add_side(elem, 1, 1);
}
break;
}
default:
libmesh_error_msg("ERROR: Unrecognized 1D element type.");
}
// Scale the nodal positions
for (unsigned int p=0; p<mesh.n_nodes(); p++)
mesh.node(p)(0) = (mesh.node(p)(0))*(xmax-xmin) + xmin;
// Add sideset names to boundary info
boundary_info.sideset_name(0) = "left";
boundary_info.sideset_name(1) = "right";
// Add nodeset names to boundary info
boundary_info.nodeset_name(0) = "left";
boundary_info.nodeset_name(1) = "right";
break;
}
//---------------------------------------------------------------------
// Build a 2D quadrilateral
case 2:
{
libmesh_assert_not_equal_to (nx, 0);
libmesh_assert_not_equal_to (ny, 0);
libmesh_assert_equal_to (nz, 0);
libmesh_assert_less (xmin, xmax);
libmesh_assert_less (ymin, ymax);
// Reserve elements. The TRI3 and TRI6 meshes
// have twice as many elements...
switch (type)
{
case INVALID_ELEM:
case QUAD4:
case QUAD8:
case QUAD9:
{
mesh.reserve_elem (nx*ny);
break;
}
case TRI3:
case TRI6:
{
mesh.reserve_elem (2*nx*ny);
break;
}
default:
libmesh_error_msg("ERROR: Unrecognized 2D element type.");
}
// Reserve nodes. The quadratic element types
// need to reserve more nodes than the linear types.
switch (type)
{
case INVALID_ELEM:
case QUAD4:
case TRI3:
{
mesh.reserve_nodes( (nx+1)*(ny+1) );
break;
}
case QUAD8:
case QUAD9:
case TRI6:
{
mesh.reserve_nodes( (2*nx+1)*(2*ny+1) );
break;
}
default:
libmesh_error_msg("ERROR: Unrecognized 2D element type.");
}
// Build the nodes. Depends on whether you are using a linear
// or quadratic element, and whether you are using a uniform
// grid or the Gauss-Lobatto grid points.
unsigned int node_id = 0;
switch (type)
{
case INVALID_ELEM:
case QUAD4:
case TRI3:
{
for (unsigned int j=0; j<=ny; j++)
for (unsigned int i=0; i<=nx; i++)
{
if (gauss_lobatto_grid)
{
mesh.add_point (Point(0.5*(1.0 - std::cos(libMesh::pi*static_cast<Real>(i)/static_cast<Real>(nx))),
0.5*(1.0 - std::cos(libMesh::pi*static_cast<Real>(j)/static_cast<Real>(ny))),
0.), node_id++);
}
else
mesh.add_point (Point(static_cast<Real>(i)/static_cast<Real>(nx),
static_cast<Real>(j)/static_cast<Real>(ny),
0.), node_id++);
}
break;
}
case QUAD8:
case QUAD9:
case TRI6:
{
for (unsigned int j=0; j<=(2*ny); j++)
for (unsigned int i=0; i<=(2*nx); i++)
{
if (gauss_lobatto_grid)
{
// The x,y locations of the point.
Real x=0., y=0.;
// Shortcut quantities (do not depend on i,j)
const Real a = std::cos( libMesh::pi / static_cast<Real>(2*nx) );
const Real b = std::cos( libMesh::pi / static_cast<Real>(2*ny) );
// Shortcut quantities (depend on i,j)
const Real c = std::cos( libMesh::pi*i / static_cast<Real>(2*nx) );
const Real d = std::cos( libMesh::pi*j / static_cast<Real>(2*ny) );
// If i is even, compute a normal Gauss-Lobatto point
if (i%2 == 0)
x = 0.5*(1.0 - c);
// Otherwise, it is the average of the previous and next points
else
x = 0.5*(1.0 - a*c);
// If j is even, compute a normal Gauss-Lobatto point
if (j%2 == 0)
y = 0.5*(1.0 - d);
// Otherwise, it is the average of the previous and next points
else
y = 0.5*(1.0 - b*d);
mesh.add_point (Point(x,y,0.), node_id++);
}
else
mesh.add_point (Point(static_cast<Real>(i)/static_cast<Real>(2*nx),
static_cast<Real>(j)/static_cast<Real>(2*ny),
0), node_id++);
}
break;
}
default:
libmesh_error_msg("ERROR: Unrecognized 2D element type.");
}
// Build the elements. Each one is a bit different.
switch (type)
{
case INVALID_ELEM:
case QUAD4:
{
for (unsigned int j=0; j<ny; j++)
for (unsigned int i=0; i<nx; i++)
{
Elem* elem = mesh.add_elem(new Quad4);
elem->set_node(0) = mesh.node_ptr(idx(type,nx,i,j) );
elem->set_node(1) = mesh.node_ptr(idx(type,nx,i+1,j) );
elem->set_node(2) = mesh.node_ptr(idx(type,nx,i+1,j+1));
elem->set_node(3) = mesh.node_ptr(idx(type,nx,i,j+1) );
if (j == 0)
boundary_info.add_side(elem, 0, 0);
if (j == (ny-1))
boundary_info.add_side(elem, 2, 2);
if (i == 0)
boundary_info.add_side(elem, 3, 3);
if (i == (nx-1))
boundary_info.add_side(elem, 1, 1);
}
break;
}
case TRI3:
{
for (unsigned int j=0; j<ny; j++)
for (unsigned int i=0; i<nx; i++)
{
Elem* elem = NULL;
// Add first Tri3
elem = mesh.add_elem(new Tri3);
elem->set_node(0) = mesh.node_ptr(idx(type,nx,i,j) );
elem->set_node(1) = mesh.node_ptr(idx(type,nx,i+1,j) );
elem->set_node(2) = mesh.node_ptr(idx(type,nx,i+1,j+1));
if (j == 0)
boundary_info.add_side(elem, 0, 0);
if (i == (nx-1))
boundary_info.add_side(elem, 1, 1);
// Add second Tri3
elem = mesh.add_elem(new Tri3);
elem->set_node(0) = mesh.node_ptr(idx(type,nx,i,j) );
elem->set_node(1) = mesh.node_ptr(idx(type,nx,i+1,j+1));
elem->set_node(2) = mesh.node_ptr(idx(type,nx,i,j+1) );
if (j == (ny-1))
boundary_info.add_side(elem, 1, 2);
if (i == 0)
boundary_info.add_side(elem, 2, 3);
}
break;
}
case QUAD8:
case QUAD9:
{
for (unsigned int j=0; j<(2*ny); j += 2)
for (unsigned int i=0; i<(2*nx); i += 2)
{
Elem* elem = (type == QUAD8) ?
mesh.add_elem(new Quad8) :
mesh.add_elem(new Quad9);
elem->set_node(0) = mesh.node_ptr(idx(type,nx,i,j) );
elem->set_node(1) = mesh.node_ptr(idx(type,nx,i+2,j) );
elem->set_node(2) = mesh.node_ptr(idx(type,nx,i+2,j+2));
elem->set_node(3) = mesh.node_ptr(idx(type,nx,i,j+2) );
elem->set_node(4) = mesh.node_ptr(idx(type,nx,i+1,j) );
elem->set_node(5) = mesh.node_ptr(idx(type,nx,i+2,j+1));
elem->set_node(6) = mesh.node_ptr(idx(type,nx,i+1,j+2));
elem->set_node(7) = mesh.node_ptr(idx(type,nx,i,j+1) );
if (type == QUAD9)
elem->set_node(8) = mesh.node_ptr(idx(type,nx,i+1,j+1));
if (j == 0)
boundary_info.add_side(elem, 0, 0);
if (j == 2*(ny-1))
boundary_info.add_side(elem, 2, 2);
if (i == 0)
boundary_info.add_side(elem, 3, 3);
if (i == 2*(nx-1))
boundary_info.add_side(elem, 1, 1);
}
break;
}
case TRI6:
{
for (unsigned int j=0; j<(2*ny); j += 2)
for (unsigned int i=0; i<(2*nx); i += 2)
{
Elem* elem = NULL;
// Add first Tri6
elem = mesh.add_elem(new Tri6);
elem->set_node(0) = mesh.node_ptr(idx(type,nx,i,j) );
elem->set_node(1) = mesh.node_ptr(idx(type,nx,i+2,j) );
elem->set_node(2) = mesh.node_ptr(idx(type,nx,i+2,j+2));
elem->set_node(3) = mesh.node_ptr(idx(type,nx,i+1,j) );
elem->set_node(4) = mesh.node_ptr(idx(type,nx,i+2,j+1));
elem->set_node(5) = mesh.node_ptr(idx(type,nx,i+1,j+1));
if (j == 0)
boundary_info.add_side(elem, 0, 0);
if (i == 2*(nx-1))
boundary_info.add_side(elem, 1, 1);
// Add second Tri6
elem = mesh.add_elem(new Tri6);
elem->set_node(0) = mesh.node_ptr(idx(type,nx,i,j) );
elem->set_node(1) = mesh.node_ptr(idx(type,nx,i+2,j+2));
elem->set_node(2) = mesh.node_ptr(idx(type,nx,i,j+2) );
elem->set_node(3) = mesh.node_ptr(idx(type,nx,i+1,j+1));
elem->set_node(4) = mesh.node_ptr(idx(type,nx,i+1,j+2));
elem->set_node(5) = mesh.node_ptr(idx(type,nx,i,j+1) );
if (j == 2*(ny-1))
boundary_info.add_side(elem, 1, 2);
if (i == 0)
boundary_info.add_side(elem, 2, 3);
}
break;
};
default:
libmesh_error_msg("ERROR: Unrecognized 2D element type.");
}
// Scale the nodal positions
for (unsigned int p=0; p<mesh.n_nodes(); p++)
{
mesh.node(p)(0) = (mesh.node(p)(0))*(xmax-xmin) + xmin;
mesh.node(p)(1) = (mesh.node(p)(1))*(ymax-ymin) + ymin;
}
// Add sideset names to boundary info
boundary_info.sideset_name(0) = "bottom";
boundary_info.sideset_name(1) = "right";
boundary_info.sideset_name(2) = "top";
boundary_info.sideset_name(3) = "left";
// Add nodeset names to boundary info
boundary_info.nodeset_name(0) = "bottom";
boundary_info.nodeset_name(1) = "right";
boundary_info.nodeset_name(2) = "top";
boundary_info.nodeset_name(3) = "left";
break;
}
//---------------------------------------------------------------------
// Build a 3D mesh using hexes, tets, prisms, or pyramids.
case 3:
{
libmesh_assert_not_equal_to (nx, 0);
libmesh_assert_not_equal_to (ny, 0);
libmesh_assert_not_equal_to (nz, 0);
libmesh_assert_less (xmin, xmax);
libmesh_assert_less (ymin, ymax);
libmesh_assert_less (zmin, zmax);
// Reserve elements. Meshes with prismatic elements require
// twice as many elements.
switch (type)
{
case INVALID_ELEM:
case HEX8:
case HEX20:
case HEX27:
case TET4: // TET4's are created from an initial HEX27 discretization
case TET10: // TET10's are created from an initial HEX27 discretization
case PYRAMID5: // PYRAMIDs are created from an initial HEX27 discretization
case PYRAMID13:
case PYRAMID14:
{
mesh.reserve_elem(nx*ny*nz);
break;
}
case PRISM6:
case PRISM15:
case PRISM18:
{
mesh.reserve_elem(2*nx*ny*nz);
break;
}
default:
libmesh_error_msg("ERROR: Unrecognized 3D element type.");
}
// Reserve nodes. Quadratic elements need twice as many nodes as linear elements.
switch (type)
{
case INVALID_ELEM:
case HEX8:
case PRISM6:
{
mesh.reserve_nodes( (nx+1)*(ny+1)*(nz+1) );
break;
}
case HEX20:
case HEX27:
case TET4: // TET4's are created from an initial HEX27 discretization
case TET10: // TET10's are created from an initial HEX27 discretization
case PYRAMID5: // PYRAMIDs are created from an initial HEX27 discretization
case PYRAMID13:
case PYRAMID14:
case PRISM15:
case PRISM18:
{
// FYI: The resulting TET4 mesh will have exactly
// 5*(nx*ny*nz) + 2*(nx*ny + nx*nz + ny*nz) + (nx+ny+nz) + 1
// nodes once the additional mid-edge nodes for the HEX27 discretization
// have been deleted.
mesh.reserve_nodes( (2*nx+1)*(2*ny+1)*(2*nz+1) );
break;
}
default:
libmesh_error_msg("ERROR: Unrecognized 3D element type.");
}
// Build the nodes.
unsigned int node_id = 0;
switch (type)
{
case INVALID_ELEM:
case HEX8:
case PRISM6:
{
for (unsigned int k=0; k<=nz; k++)
for (unsigned int j=0; j<=ny; j++)
for (unsigned int i=0; i<=nx; i++)
{
if (gauss_lobatto_grid)
{
mesh.add_point (Point(0.5*(1.0 - std::cos(libMesh::pi*static_cast<Real>(i)/static_cast<Real>(nx))),
0.5*(1.0 - std::cos(libMesh::pi*static_cast<Real>(j)/static_cast<Real>(ny))),
0.5*(1.0 - std::cos(libMesh::pi*static_cast<Real>(k)/static_cast<Real>(nz)))), node_id++);
}
else
mesh.add_point(Point(static_cast<Real>(i)/static_cast<Real>(nx),
static_cast<Real>(j)/static_cast<Real>(ny),
static_cast<Real>(k)/static_cast<Real>(nz)), node_id++);
}
break;
}
case HEX20:
case HEX27:
case TET4: // TET4's are created from an initial HEX27 discretization
case TET10: // TET10's are created from an initial HEX27 discretization
case PYRAMID5: // PYRAMIDs are created from an initial HEX27 discretization
case PYRAMID13:
case PYRAMID14:
case PRISM15:
case PRISM18:
{
for (unsigned int k=0; k<=(2*nz); k++)
for (unsigned int j=0; j<=(2*ny); j++)
for (unsigned int i=0; i<=(2*nx); i++)
{
if (gauss_lobatto_grid)
{
// The x,y locations of the point.
Real x=0., y=0., z=0.;
// Shortcut quantities (do not depend on i,j)
const Real a = std::cos( libMesh::pi / static_cast<Real>(2*nx) );
const Real b = std::cos( libMesh::pi / static_cast<Real>(2*ny) );
// Shortcut quantities (depend on i,j)
const Real c = std::cos( libMesh::pi*i / static_cast<Real>(2*nx) );
const Real d = std::cos( libMesh::pi*j / static_cast<Real>(2*ny) );
// Additional shortcut quantities (for 3D)
const Real e = std::cos( libMesh::pi / static_cast<Real>(2*nz) );
const Real f = std::cos( libMesh::pi*k / static_cast<Real>(2*nz) );
// If i is even, compute a normal Gauss-Lobatto point
if (i%2 == 0)
x = 0.5*(1.0 - c);
// Otherwise, it is the average of the previous and next points
else
x = 0.5*(1.0 - a*c);
// If j is even, compute a normal Gauss-Lobatto point
if (j%2 == 0)
y = 0.5*(1.0 - d);
// Otherwise, it is the average of the previous and next points
else
y = 0.5*(1.0 - b*d);
// If k is even, compute a normal Gauss-Lobatto point
if (k%2 == 0)
z = 0.5*(1.0 - f);
// Otherwise, it is the average of the previous and next points
else
z = 0.5*(1.0 - e*f);
mesh.add_point (Point(x,y,z), node_id++);
}
else
mesh.add_point(Point(static_cast<Real>(i)/static_cast<Real>(2*nx),
static_cast<Real>(j)/static_cast<Real>(2*ny),
static_cast<Real>(k)/static_cast<Real>(2*nz)), node_id++);
}
break;
}
default:
libmesh_error_msg("ERROR: Unrecognized 3D element type.");
}
// Build the elements.
switch (type)
{
case INVALID_ELEM:
case HEX8:
{
for (unsigned int k=0; k<nz; k++)
for (unsigned int j=0; j<ny; j++)
for (unsigned int i=0; i<nx; i++)
{
Elem* elem = mesh.add_elem(new Hex8);
elem->set_node(0) = mesh.node_ptr(idx(type,nx,ny,i,j,k) );
elem->set_node(1) = mesh.node_ptr(idx(type,nx,ny,i+1,j,k) );
elem->set_node(2) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) );
elem->set_node(3) = mesh.node_ptr(idx(type,nx,ny,i,j+1,k) );
elem->set_node(4) = mesh.node_ptr(idx(type,nx,ny,i,j,k+1) );
elem->set_node(5) = mesh.node_ptr(idx(type,nx,ny,i+1,j,k+1) );
elem->set_node(6) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1));
elem->set_node(7) = mesh.node_ptr(idx(type,nx,ny,i,j+1,k+1) );
if (k == 0)
boundary_info.add_side(elem, 0, 0);
if (k == (nz-1))
boundary_info.add_side(elem, 5, 5);
if (j == 0)
boundary_info.add_side(elem, 1, 1);
if (j == (ny-1))
boundary_info.add_side(elem, 3, 3);
if (i == 0)
boundary_info.add_side(elem, 4, 4);
if (i == (nx-1))
boundary_info.add_side(elem, 2, 2);
}
break;
}
case PRISM6:
{
for (unsigned int k=0; k<nz; k++)
for (unsigned int j=0; j<ny; j++)
for (unsigned int i=0; i<nx; i++)
{
// First Prism
Elem* elem = NULL;
elem = mesh.add_elem(new Prism6);
elem->set_node(0) = mesh.node_ptr(idx(type,nx,ny,i,j,k) );
elem->set_node(1) = mesh.node_ptr(idx(type,nx,ny,i+1,j,k) );
elem->set_node(2) = mesh.node_ptr(idx(type,nx,ny,i,j+1,k) );
elem->set_node(3) = mesh.node_ptr(idx(type,nx,ny,i,j,k+1) );
elem->set_node(4) = mesh.node_ptr(idx(type,nx,ny,i+1,j,k+1) );
elem->set_node(5) = mesh.node_ptr(idx(type,nx,ny,i,j+1,k+1) );
// Add sides for first prism to boundary info object
if (i==0)
boundary_info.add_side(elem, 3, 4);
if (j==0)
boundary_info.add_side(elem, 1, 1);
if (k==0)
boundary_info.add_side(elem, 0, 0);
if (k == (nz-1))
boundary_info.add_side(elem, 4, 5);
// Second Prism
elem = mesh.add_elem(new Prism6);
elem->set_node(0) = mesh.node_ptr(idx(type,nx,ny,i+1,j,k) );
elem->set_node(1) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) );
elem->set_node(2) = mesh.node_ptr(idx(type,nx,ny,i,j+1,k) );
elem->set_node(3) = mesh.node_ptr(idx(type,nx,ny,i+1,j,k+1) );
elem->set_node(4) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1));
elem->set_node(5) = mesh.node_ptr(idx(type,nx,ny,i,j+1,k+1) );
// Add sides for second prism to boundary info object
if (i == (nx-1))
boundary_info.add_side(elem, 1, 2);
if (j == (ny-1))
boundary_info.add_side(elem, 2, 3);
if (k==0)
boundary_info.add_side(elem, 0, 0);
if (k == (nz-1))
boundary_info.add_side(elem, 4, 5);
}
break;
}
case HEX20:
case HEX27:
case TET4: // TET4's are created from an initial HEX27 discretization
case TET10: // TET10's are created from an initial HEX27 discretization
case PYRAMID5: // PYRAMIDs are created from an initial HEX27 discretization
case PYRAMID13:
case PYRAMID14:
{
for (unsigned int k=0; k<(2*nz); k += 2)
for (unsigned int j=0; j<(2*ny); j += 2)
for (unsigned int i=0; i<(2*nx); i += 2)
{
Elem* elem = (type == HEX20) ?
mesh.add_elem(new Hex20) :
mesh.add_elem(new Hex27);
elem->set_node(0) = mesh.node_ptr(idx(type,nx,ny,i, j, k) );
elem->set_node(1) = mesh.node_ptr(idx(type,nx,ny,i+2,j, k) );
elem->set_node(2) = mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k) );
elem->set_node(3) = mesh.node_ptr(idx(type,nx,ny,i, j+2,k) );
elem->set_node(4) = mesh.node_ptr(idx(type,nx,ny,i, j, k+2));
elem->set_node(5) = mesh.node_ptr(idx(type,nx,ny,i+2,j, k+2));
elem->set_node(6) = mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k+2));
elem->set_node(7) = mesh.node_ptr(idx(type,nx,ny,i, j+2,k+2));
elem->set_node(8) = mesh.node_ptr(idx(type,nx,ny,i+1,j, k) );
elem->set_node(9) = mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k) );
elem->set_node(10) = mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k) );
elem->set_node(11) = mesh.node_ptr(idx(type,nx,ny,i, j+1,k) );
elem->set_node(12) = mesh.node_ptr(idx(type,nx,ny,i, j, k+1));
elem->set_node(13) = mesh.node_ptr(idx(type,nx,ny,i+2,j, k+1));
elem->set_node(14) = mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k+1));
elem->set_node(15) = mesh.node_ptr(idx(type,nx,ny,i, j+2,k+1));
elem->set_node(16) = mesh.node_ptr(idx(type,nx,ny,i+1,j, k+2));
elem->set_node(17) = mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k+2));
elem->set_node(18) = mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k+2));
elem->set_node(19) = mesh.node_ptr(idx(type,nx,ny,i, j+1,k+2));
if ((type == HEX27) || (type == TET4) || (type == TET10) ||
(type == PYRAMID5) || (type == PYRAMID13) || (type == PYRAMID14))
{
elem->set_node(20) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) );
elem->set_node(21) = mesh.node_ptr(idx(type,nx,ny,i+1,j, k+1));
elem->set_node(22) = mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k+1));
elem->set_node(23) = mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k+1));
elem->set_node(24) = mesh.node_ptr(idx(type,nx,ny,i, j+1,k+1));
elem->set_node(25) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+2));
elem->set_node(26) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1));
}
if (k == 0)
boundary_info.add_side(elem, 0, 0);
if (k == 2*(nz-1))
boundary_info.add_side(elem, 5, 5);
if (j == 0)
boundary_info.add_side(elem, 1, 1);
if (j == 2*(ny-1))
boundary_info.add_side(elem, 3, 3);
if (i == 0)
boundary_info.add_side(elem, 4, 4);
if (i == 2*(nx-1))
boundary_info.add_side(elem, 2, 2);
}
break;
}
case PRISM15:
case PRISM18:
{
for (unsigned int k=0; k<(2*nz); k += 2)
for (unsigned int j=0; j<(2*ny); j += 2)
for (unsigned int i=0; i<(2*nx); i += 2)
{
// First Prism
Elem* elem = NULL;
elem = ((type == PRISM15) ?
mesh.add_elem(new Prism15) :
mesh.add_elem(new Prism18));
elem->set_node(0) = mesh.node_ptr(idx(type,nx,ny,i, j, k) );
elem->set_node(1) = mesh.node_ptr(idx(type,nx,ny,i+2,j, k) );
elem->set_node(2) = mesh.node_ptr(idx(type,nx,ny,i, j+2,k) );
elem->set_node(3) = mesh.node_ptr(idx(type,nx,ny,i, j, k+2));
elem->set_node(4) = mesh.node_ptr(idx(type,nx,ny,i+2,j, k+2));
elem->set_node(5) = mesh.node_ptr(idx(type,nx,ny,i, j+2,k+2));
elem->set_node(6) = mesh.node_ptr(idx(type,nx,ny,i+1,j, k) );
elem->set_node(7) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) );
elem->set_node(8) = mesh.node_ptr(idx(type,nx,ny,i, j+1,k) );
elem->set_node(9) = mesh.node_ptr(idx(type,nx,ny,i, j, k+1));
elem->set_node(10) = mesh.node_ptr(idx(type,nx,ny,i+2,j, k+1));
elem->set_node(11) = mesh.node_ptr(idx(type,nx,ny,i, j+2,k+1));
elem->set_node(12) = mesh.node_ptr(idx(type,nx,ny,i+1,j, k+2));
elem->set_node(13) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+2));
elem->set_node(14) = mesh.node_ptr(idx(type,nx,ny,i, j+1,k+2));
if (type == PRISM18)
{
elem->set_node(15) = mesh.node_ptr(idx(type,nx,ny,i+1,j, k+1));
elem->set_node(16) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1));
elem->set_node(17) = mesh.node_ptr(idx(type,nx,ny,i, j+1,k+1));
}
// Add sides for first prism to boundary info object
if (i==0)
boundary_info.add_side(elem, 3, 4);
if (j==0)
boundary_info.add_side(elem, 1, 1);
if (k==0)
boundary_info.add_side(elem, 0, 0);
if (k == 2*(nz-1))
boundary_info.add_side(elem, 4, 5);
// Second Prism
elem = ((type == PRISM15) ?
mesh.add_elem(new Prism15) :
mesh.add_elem(new Prism18));
elem->set_node(0) = mesh.node_ptr(idx(type,nx,ny,i+2,j,k) );
elem->set_node(1) = mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k) );
elem->set_node(2) = mesh.node_ptr(idx(type,nx,ny,i,j+2,k) );
elem->set_node(3) = mesh.node_ptr(idx(type,nx,ny,i+2,j,k+2) );
elem->set_node(4) = mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k+2) );
elem->set_node(5) = mesh.node_ptr(idx(type,nx,ny,i,j+2,k+2) );
elem->set_node(6) = mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k) );
elem->set_node(7) = mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k) );
elem->set_node(8) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k) );
elem->set_node(9) = mesh.node_ptr(idx(type,nx,ny,i+2,j,k+1) );
elem->set_node(10) = mesh.node_ptr(idx(type,nx,ny,i+2,j+2,k+1));
elem->set_node(11) = mesh.node_ptr(idx(type,nx,ny,i,j+2,k+1) );
elem->set_node(12) = mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k+2));
elem->set_node(13) = mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k+2));
elem->set_node(14) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+2));
if (type == PRISM18)
{
elem->set_node(15) = mesh.node_ptr(idx(type,nx,ny,i+2,j+1,k+1));
elem->set_node(16) = mesh.node_ptr(idx(type,nx,ny,i+1,j+2,k+1));
elem->set_node(17) = mesh.node_ptr(idx(type,nx,ny,i+1,j+1,k+1));
}
// Add sides for second prism to boundary info object
if (i == 2*(nx-1))
boundary_info.add_side(elem, 1, 2);
if (j == 2*(ny-1))
boundary_info.add_side(elem, 2, 3);
if (k==0)
boundary_info.add_side(elem, 0, 0);
if (k == 2*(nz-1))
boundary_info.add_side(elem, 4, 5);
}
break;
}
default:
libmesh_error_msg("ERROR: Unrecognized 3D element type.");
}
//.......................................
// Scale the nodal positions
for (unsigned int p=0; p<mesh.n_nodes(); p++)
{
mesh.node(p)(0) = (mesh.node(p)(0))*(xmax-xmin) + xmin;
mesh.node(p)(1) = (mesh.node(p)(1))*(ymax-ymin) + ymin;
mesh.node(p)(2) = (mesh.node(p)(2))*(zmax-zmin) + zmin;
}
// Additional work for tets and pyramids: we take the existing
// HEX27 discretization and split each element into 24
// sub-tets or 6 sub-pyramids.
//
// 24 isn't the minimum-possible number of tets, but it
// obviates any concerns about the edge orientations between
// the various elements.
if ((type == TET4) ||
(type == TET10) ||
(type == PYRAMID5) ||
(type == PYRAMID13) ||
(type == PYRAMID14))
{
// Temporary storage for new elements. (24 tets per hex, 6 pyramids)
std::vector<Elem*> new_elements;
if ((type == TET4) || (type == TET10))
new_elements.reserve(24*mesh.n_elem());
else
new_elements.reserve(6*mesh.n_elem());
// Create tetrahedra or pyramids
{
MeshBase::element_iterator el = mesh.elements_begin();
const MeshBase::element_iterator end_el = mesh.elements_end();
for ( ; el != end_el; ++el)
{
// Get a pointer to the HEX27 element.
Elem* base_hex = *el;
// Get a pointer to the node located at the HEX27 centroid
Node* apex_node = base_hex->get_node(26);
for (unsigned short s=0; s<base_hex->n_sides(); ++s)
{
// Get the boundary ID for this side
boundary_id_type b_id = boundary_info.boundary_id(*el, s);
// Need to build the full-ordered side!
UniquePtr<Elem> side = base_hex->build_side(s);
if ((type == TET4) || (type == TET10))
{
// Build 4 sub-tets per side
for (unsigned int sub_tet=0; sub_tet<4; ++sub_tet)
{
new_elements.push_back( new Tet4 );
Elem* sub_elem = new_elements.back();
sub_elem->set_node(0) = side->get_node(sub_tet);
sub_elem->set_node(1) = side->get_node(8); // centroid of the face
sub_elem->set_node(2) = side->get_node(sub_tet==3 ? 0 : sub_tet+1 ); // wrap-around
sub_elem->set_node(3) = apex_node; // apex node always used!
// If the original hex was a boundary hex, add the new sub_tet's side
// 0 with the same b_id. Note: the tets are all aligned so that their
// side 0 is on the boundary.
if (b_id != BoundaryInfo::invalid_id)
boundary_info.add_side(sub_elem, 0, b_id);
}
} // end if ((type == TET4) || (type == TET10))
else // type==PYRAMID5 || type==PYRAMID13 || type==PYRAMID14
{
// Build 1 sub-pyramid per side.
new_elements.push_back(new Pyramid5);
Elem* sub_elem = new_elements.back();
// Set the base. Note that since the apex is *inside* the base_hex,
// and the pyramid uses a counter-clockwise base numbering, we need to
// reverse the [1] and [3] node indices.
sub_elem->set_node(0) = side->get_node(0);
sub_elem->set_node(1) = side->get_node(3);
sub_elem->set_node(2) = side->get_node(2);
sub_elem->set_node(3) = side->get_node(1);
// Set the apex
sub_elem->set_node(4) = apex_node;
// If the original hex was a boundary hex, add the new sub_pyr's side
// 4 (the square base) with the same b_id.
if (b_id != BoundaryInfo::invalid_id)
boundary_info.add_side(sub_elem, 4, b_id);
} // end else type==PYRAMID5 || type==PYRAMID13 || type==PYRAMID14
}
}
}
// Delete the original HEX27 elements from the mesh, and the boundary info structure.
{
MeshBase::element_iterator el = mesh.elements_begin();
const MeshBase::element_iterator end_el = mesh.elements_end();
for ( ; el != end_el; ++el)
{
boundary_info.remove(*el); // Safe even if *el has no boundary info.
mesh.delete_elem(*el);
}
}
// Add the new elements
for (unsigned int i=0; i<new_elements.size(); ++i)
mesh.add_elem(new_elements[i]);
} // end if (type == TET4,TET10,PYRAMID5,PYRAMID13,PYRAMID14
// Use all_second_order to convert the TET4's to TET10's or PYRAMID5's to PYRAMID14's
if ((type == TET10) || (type == PYRAMID14))
mesh.all_second_order();
else if (type == PYRAMID13)
mesh.all_second_order(/*full_ordered=*/false);
// Add sideset names to boundary info (Z axis out of the screen)
boundary_info.sideset_name(0) = "back";
boundary_info.sideset_name(1) = "bottom";
boundary_info.sideset_name(2) = "right";
boundary_info.sideset_name(3) = "top";
boundary_info.sideset_name(4) = "left";
boundary_info.sideset_name(5) = "front";
// Add nodeset names to boundary info
boundary_info.nodeset_name(0) = "back";
boundary_info.nodeset_name(1) = "bottom";
boundary_info.nodeset_name(2) = "right";
boundary_info.nodeset_name(3) = "top";
boundary_info.nodeset_name(4) = "left";
boundary_info.nodeset_name(5) = "front";
break;
} // end case dim==3
default:
libmesh_error_msg("Unknown dimension " << mesh.mesh_dimension());
}
STOP_LOG("build_cube()", "MeshTools::Generation");
// Done building the mesh. Now prepare it for use.
mesh.prepare_for_use (/*skip_renumber =*/ false);
}
| 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 2201 of file mesh_generation.C.
References libMesh::MeshBase::add_point(), libMesh::BoundaryInfo::add_side(), libMesh::TriangleInterface::attach_hole_list(), bc_id, libMesh::Elem::build_side(), libMesh::MeshBase::clear(), libMesh::TriangleInterface::desired_area(), libMesh::TriangleInterface::elem_type(), libMesh::MeshBase::elements_begin(), libMesh::MeshBase::elements_end(), libMesh::MeshBase::get_boundary_info(), 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().
{
// Check for reasonable size
libmesh_assert_greater_equal (nx, 1); // need at least 1 element in x-direction
libmesh_assert_greater_equal (ny, 1); // need at least 1 element in y-direction
libmesh_assert_less (xmin, xmax);
libmesh_assert_less (ymin, ymax);
// Clear out any data which may have been in the Mesh
mesh.clear();
BoundaryInfo& boundary_info = mesh.get_boundary_info();
// Make sure the new Mesh will be 2D
mesh.set_mesh_dimension(2);
// The x and y spacing between boundary points
const Real delta_x = (xmax-xmin) / static_cast<Real>(nx);
const Real delta_y = (ymax-ymin) / static_cast<Real>(ny);
// Bottom
for (unsigned int p=0; p<=nx; ++p)
mesh.add_point(Point(xmin + p*delta_x, ymin));
// Right side
for (unsigned int p=1; p<ny; ++p)
mesh.add_point(Point(xmax, ymin + p*delta_y));
// Top
for (unsigned int p=0; p<=nx; ++p)
mesh.add_point(Point(xmax - p*delta_x, ymax));
// Left side
for (unsigned int p=1; p<ny; ++p)
mesh.add_point(Point(xmin, ymax - p*delta_y));
// Be sure we added as many points as we thought we did
libmesh_assert_equal_to (mesh.n_nodes(), 2*(nx+ny));
// Construct the Triangle Interface object
TriangleInterface t(mesh);
// Set custom variables for the triangulation
t.desired_area() = 0.5 * (xmax-xmin)*(ymax-ymin) / static_cast<Real>(nx*ny);
t.triangulation_type() = TriangleInterface::PSLG;
t.elem_type() = type;
if (holes != NULL)
t.attach_hole_list(holes);
// Triangulate!
t.triangulate();
// The mesh is now generated, but we still need to mark the boundaries
// to be consistent with the other build_square routines. Note that all
// hole boundary elements get the same ID, 4.
MeshBase::element_iterator el = mesh.elements_begin();
const MeshBase::element_iterator end_el = mesh.elements_end();
for ( ; el != end_el; ++el)
{
const Elem* elem = *el;
for (unsigned int s=0; s<elem->n_sides(); s++)
if (elem->neighbor(s) == NULL)
{
UniquePtr<Elem> side (elem->build_side(s));
// Check the location of the side's midpoint. Since
// the square has straight sides, the midpoint is not
// on the corner and thus it is uniquely on one of the
// sides.
Point side_midpoint= 0.5f*( (*side->get_node(0)) + (*side->get_node(1)) );
// The boundary ids are set following the same convention as Quad4 sides
// bottom = 0
// right = 1
// top = 2
// left = 3
// hole = 4
boundary_id_type bc_id=4;
// bottom
if (std::fabs(side_midpoint(1) - ymin) < TOLERANCE)
bc_id=0;
// right
else if (std::fabs(side_midpoint(0) - xmax) < TOLERANCE)
bc_id=1;
// top
else if (std::fabs(side_midpoint(1) - ymax) < TOLERANCE)
bc_id=2;
// left
else if (std::fabs(side_midpoint(0) - xmin) < TOLERANCE)
bc_id=3;
// If the point is not on any of the external boundaries, it
// is on one of the holes....
// Finally, add this element's information to the boundary info object.
boundary_info.add_side(elem->id(), s, bc_id);
}
}
} // 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 1954 of file mesh_generation.C.
References libMesh::MeshBase::add_elem(), libMesh::BoundaryInfo::add_node(), libMesh::MeshBase::add_point(), libMesh::BoundaryInfo::add_side(), libMesh::BoundaryInfo::boundary_ids(), libMesh::MeshBase::delete_remote_elements(), libMesh::Elem::dim(), libMesh::EDGE2, libMesh::EDGE3, libMesh::MeshBase::elements_begin(), libMesh::MeshBase::elements_end(), end, libMesh::MeshBase::get_boundary_info(), libMesh::Elem::get_node(), libMesh::BoundaryInfo::get_side_boundary_ids(), 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(), libMesh::QUAD4, libMesh::QUAD9, libMesh::MeshBase::reserve_elem(), libMesh::MeshBase::reserve_nodes(), libMesh::SECOND, libMesh::DofObject::set_id(), libMesh::Elem::set_node(), libMesh::START_LOG(), libMesh::Elem::subdomain_id(), libMesh::TRI3, libMesh::TRI6, and libMesh::Elem::type().
{
if (!cross_section.n_elem())
return;
START_LOG("build_extrusion()", "MeshTools::Generation");
dof_id_type orig_elem = cross_section.n_elem();
dof_id_type orig_nodes = cross_section.n_nodes();
unsigned int order = 1;
BoundaryInfo& boundary_info = mesh.get_boundary_info();
const BoundaryInfo& cross_section_boundary_info = cross_section.get_boundary_info();
// If cross_section is distributed, so is its extrusion
if (!cross_section.is_serial())
mesh.delete_remote_elements();
// We know a priori how many elements we'll need
mesh.reserve_elem(nz*orig_elem);
// For straightforward meshes we need one or two additional layers per
// element.
if ((*cross_section.elements_begin())->default_order() == SECOND)
order = 2;
mesh.reserve_nodes((order*nz+1)*orig_nodes);
MeshBase::const_node_iterator nd = cross_section.nodes_begin();
const MeshBase::const_node_iterator nend = cross_section.nodes_end();
for (; nd!=nend; ++nd)
{
const Node* node = *nd;
for (unsigned int k=0; k != order*nz+1; ++k)
{
Node *new_node =
mesh.add_point(*node +
(extrusion_vector * k / nz / order),
node->id() + (k * orig_nodes),
node->processor_id());
const std::vector<boundary_id_type> ids_to_copy =
cross_section_boundary_info.boundary_ids(node);
boundary_info.add_node(new_node, ids_to_copy);
}
}
const std::set<boundary_id_type> &side_ids =
cross_section_boundary_info.get_side_boundary_ids();
const boundary_id_type next_side_id = side_ids.empty() ?
0 : cast_int<boundary_id_type>(*side_ids.rbegin() + 1);
MeshBase::const_element_iterator el = cross_section.elements_begin();
const MeshBase::const_element_iterator end = cross_section.elements_end();
for (; el!=end; ++el)
{
const Elem* elem = *el;
const ElemType etype = elem->type();
// build_extrusion currently only works on coarse meshes
libmesh_assert (!elem->parent());
for (unsigned int k=0; k != nz; ++k)
{
Elem *new_elem;
switch (etype)
{
case EDGE2:
{
new_elem = new Quad4;
new_elem->set_node(0) = mesh.node_ptr(elem->get_node(0)->id() + (k * orig_nodes));
new_elem->set_node(1) = mesh.node_ptr(elem->get_node(1)->id() + (k * orig_nodes));
new_elem->set_node(2) = mesh.node_ptr(elem->get_node(1)->id() + ((k+1) * orig_nodes));
new_elem->set_node(3) = mesh.node_ptr(elem->get_node(0)->id() + ((k+1) * orig_nodes));
break;
}
case EDGE3:
{
new_elem = new Quad9;
new_elem->set_node(0) = mesh.node_ptr(elem->get_node(0)->id() + (2*k * orig_nodes));
new_elem->set_node(1) = mesh.node_ptr(elem->get_node(1)->id() + (2*k * orig_nodes));
new_elem->set_node(2) = mesh.node_ptr(elem->get_node(1)->id() + ((2*k+2) * orig_nodes));
new_elem->set_node(3) = mesh.node_ptr(elem->get_node(0)->id() + ((2*k+2) * orig_nodes));
new_elem->set_node(4) = mesh.node_ptr(elem->get_node(2)->id() + (2*k * orig_nodes));
new_elem->set_node(5) = mesh.node_ptr(elem->get_node(1)->id() + ((2*k+1) * orig_nodes));
new_elem->set_node(6) = mesh.node_ptr(elem->get_node(2)->id() + ((2*k+2) * orig_nodes));
new_elem->set_node(7) = mesh.node_ptr(elem->get_node(0)->id() + ((2*k+1) * orig_nodes));
new_elem->set_node(8) = mesh.node_ptr(elem->get_node(2)->id() + ((2*k+1) * orig_nodes));
break;
}
case TRI3:
{
new_elem = new Prism6;
new_elem->set_node(0) = mesh.node_ptr(elem->get_node(0)->id() + (k * orig_nodes));
new_elem->set_node(1) = mesh.node_ptr(elem->get_node(1)->id() + (k * orig_nodes));
new_elem->set_node(2) = mesh.node_ptr(elem->get_node(2)->id() + (k * orig_nodes));
new_elem->set_node(3) = mesh.node_ptr(elem->get_node(0)->id() + ((k+1) * orig_nodes));
new_elem->set_node(4) = mesh.node_ptr(elem->get_node(1)->id() + ((k+1) * orig_nodes));
new_elem->set_node(5) = mesh.node_ptr(elem->get_node(2)->id() + ((k+1) * orig_nodes));
break;
}
case TRI6:
{
new_elem = new Prism18;
new_elem->set_node(0) = mesh.node_ptr(elem->get_node(0)->id() + (2*k * orig_nodes));
new_elem->set_node(1) = mesh.node_ptr(elem->get_node(1)->id() + (2*k * orig_nodes));
new_elem->set_node(2) = mesh.node_ptr(elem->get_node(2)->id() + (2*k * orig_nodes));
new_elem->set_node(3) = mesh.node_ptr(elem->get_node(0)->id() + ((2*k+2) * orig_nodes));
new_elem->set_node(4) = mesh.node_ptr(elem->get_node(1)->id() + ((2*k+2) * orig_nodes));
new_elem->set_node(5) = mesh.node_ptr(elem->get_node(2)->id() + ((2*k+2) * orig_nodes));
new_elem->set_node(6) = mesh.node_ptr(elem->get_node(3)->id() + (2*k * orig_nodes));
new_elem->set_node(7) = mesh.node_ptr(elem->get_node(4)->id() + (2*k * orig_nodes));
new_elem->set_node(8) = mesh.node_ptr(elem->get_node(5)->id() + (2*k * orig_nodes));
new_elem->set_node(9) = mesh.node_ptr(elem->get_node(0)->id() + ((2*k+1) * orig_nodes));
new_elem->set_node(10) = mesh.node_ptr(elem->get_node(1)->id() + ((2*k+1) * orig_nodes));
new_elem->set_node(11) = mesh.node_ptr(elem->get_node(2)->id() + ((2*k+1) * orig_nodes));
new_elem->set_node(12) = mesh.node_ptr(elem->get_node(3)->id() + ((2*k+2) * orig_nodes));
new_elem->set_node(13) = mesh.node_ptr(elem->get_node(4)->id() + ((2*k+2) * orig_nodes));
new_elem->set_node(14) = mesh.node_ptr(elem->get_node(5)->id() + ((2*k+2) * orig_nodes));
new_elem->set_node(15) = mesh.node_ptr(elem->get_node(3)->id() + ((2*k+1) * orig_nodes));
new_elem->set_node(16) = mesh.node_ptr(elem->get_node(4)->id() + ((2*k+1) * orig_nodes));
new_elem->set_node(17) = mesh.node_ptr(elem->get_node(5)->id() + ((2*k+1) * orig_nodes));
break;
}
case QUAD4:
{
new_elem = new Hex8;
new_elem->set_node(0) = mesh.node_ptr(elem->get_node(0)->id() + (k * orig_nodes));
new_elem->set_node(1) = mesh.node_ptr(elem->get_node(1)->id() + (k * orig_nodes));
new_elem->set_node(2) = mesh.node_ptr(elem->get_node(2)->id() + (k * orig_nodes));
new_elem->set_node(3) = mesh.node_ptr(elem->get_node(3)->id() + (k * orig_nodes));
new_elem->set_node(4) = mesh.node_ptr(elem->get_node(0)->id() + ((k+1) * orig_nodes));
new_elem->set_node(5) = mesh.node_ptr(elem->get_node(1)->id() + ((k+1) * orig_nodes));
new_elem->set_node(6) = mesh.node_ptr(elem->get_node(2)->id() + ((k+1) * orig_nodes));
new_elem->set_node(7) = mesh.node_ptr(elem->get_node(3)->id() + ((k+1) * orig_nodes));
break;
}
case QUAD9:
{
new_elem = new Hex27;
new_elem->set_node(0) = mesh.node_ptr(elem->get_node(0)->id() + (2*k * orig_nodes));
new_elem->set_node(1) = mesh.node_ptr(elem->get_node(1)->id() + (2*k * orig_nodes));
new_elem->set_node(2) = mesh.node_ptr(elem->get_node(2)->id() + (2*k * orig_nodes));
new_elem->set_node(3) = mesh.node_ptr(elem->get_node(3)->id() + (2*k * orig_nodes));
new_elem->set_node(4) = mesh.node_ptr(elem->get_node(0)->id() + ((2*k+2) * orig_nodes));
new_elem->set_node(5) = mesh.node_ptr(elem->get_node(1)->id() + ((2*k+2) * orig_nodes));
new_elem->set_node(6) = mesh.node_ptr(elem->get_node(2)->id() + ((2*k+2) * orig_nodes));
new_elem->set_node(7) = mesh.node_ptr(elem->get_node(3)->id() + ((2*k+2) * orig_nodes));
new_elem->set_node(8) = mesh.node_ptr(elem->get_node(4)->id() + (2*k * orig_nodes));
new_elem->set_node(9) = mesh.node_ptr(elem->get_node(5)->id() + (2*k * orig_nodes));
new_elem->set_node(10) = mesh.node_ptr(elem->get_node(6)->id() + (2*k * orig_nodes));
new_elem->set_node(11) = mesh.node_ptr(elem->get_node(7)->id() + (2*k * orig_nodes));
new_elem->set_node(12) = mesh.node_ptr(elem->get_node(0)->id() + ((2*k+1) * orig_nodes));
new_elem->set_node(13) = mesh.node_ptr(elem->get_node(1)->id() + ((2*k+1) * orig_nodes));
new_elem->set_node(14) = mesh.node_ptr(elem->get_node(2)->id() + ((2*k+1) * orig_nodes));
new_elem->set_node(15) = mesh.node_ptr(elem->get_node(3)->id() + ((2*k+1) * orig_nodes));
new_elem->set_node(16) = mesh.node_ptr(elem->get_node(4)->id() + ((2*k+2) * orig_nodes));
new_elem->set_node(17) = mesh.node_ptr(elem->get_node(5)->id() + ((2*k+2) * orig_nodes));
new_elem->set_node(18) = mesh.node_ptr(elem->get_node(6)->id() + ((2*k+2) * orig_nodes));
new_elem->set_node(19) = mesh.node_ptr(elem->get_node(7)->id() + ((2*k+2) * orig_nodes));
new_elem->set_node(20) = mesh.node_ptr(elem->get_node(8)->id() + (2*k * orig_nodes));
new_elem->set_node(21) = mesh.node_ptr(elem->get_node(4)->id() + ((2*k+1) * orig_nodes));
new_elem->set_node(22) = mesh.node_ptr(elem->get_node(5)->id() + ((2*k+1) * orig_nodes));
new_elem->set_node(23) = mesh.node_ptr(elem->get_node(6)->id() + ((2*k+1) * orig_nodes));
new_elem->set_node(24) = mesh.node_ptr(elem->get_node(7)->id() + ((2*k+1) * orig_nodes));
new_elem->set_node(25) = mesh.node_ptr(elem->get_node(8)->id() + ((2*k+2) * orig_nodes));
new_elem->set_node(26) = mesh.node_ptr(elem->get_node(8)->id() + ((2*k+1) * orig_nodes));
break;
}
default:
{
libmesh_not_implemented();
break;
}
}
new_elem->set_id(elem->id() + (k * orig_elem));
new_elem->processor_id() = elem->processor_id();
// maintain the subdomain_id
new_elem->subdomain_id() = elem->subdomain_id();
new_elem = mesh.add_elem(new_elem);
// Copy any old boundary ids on all sides
for (unsigned short s = 0; s != elem->n_sides(); ++s)
{
const std::vector<boundary_id_type> ids_to_copy =
cross_section_boundary_info.boundary_ids(elem, s);
if (new_elem->dim() == 3)
{
// For 2D->3D extrusion, we give the boundary IDs
// for side s on the old element to side s+1 on the
// new element. This is just a happy coincidence as
// far as I can tell...
boundary_info.add_side
(new_elem, cast_int<unsigned short>(s+1),
ids_to_copy);
}
else
{
// For 1D->2D extrusion, the boundary IDs map as:
// Old elem -> New elem
// 0 -> 3
// 1 -> 1
libmesh_assert_less(s, 2);
const unsigned short sidemap[2] = {3, 1};
boundary_info.add_side(new_elem, sidemap[s], ids_to_copy);
}
}
// Give new boundary ids to bottom and top
if (k == 0)
boundary_info.add_side(new_elem, 0, next_side_id);
if (k == nz-1)
{
// For 2D->3D extrusion, the "top" ID is 1+the original
// element's number of sides. For 1D->2D extrusion, the
// "top" ID is side 2.
const unsigned short top_id = new_elem->dim() == 3 ?
cast_int<unsigned short>(elem->n_sides()+1) : 2;
boundary_info.add_side
(new_elem, top_id,
cast_int<boundary_id_type>(next_side_id+1));
}
}
}
STOP_LOG("build_extrusion()", "MeshTools::Generation");
// Done building the mesh. Now prepare it for use.
mesh.prepare_for_use(/*skip_renumber =*/ false);
}
| 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 1447 of file mesh_generation.C.
References build_cube().
{
// This method only makes sense in 1D!
// But we now just turn a non-1D mesh into a 1D mesh
//libmesh_assert_equal_to (mesh.mesh_dimension(), 1);
build_cube(mesh,
nx, 0, 0,
xmin, xmax,
0., 0.,
0., 0.,
type,
gauss_lobatto_grid);
}
| 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 1429 of file mesh_generation.C.
References build_cube().
{
// This method only makes sense in 0D!
// But we now just turn a non-0D mesh into a 0D mesh
//libmesh_assert_equal_to (mesh.mesh_dimension(), 1);
build_cube(mesh,
0, 0, 0,
0., 0.,
0., 0.,
0., 0.,
type,
gauss_lobatto_grid);
}
| 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 1499 of file mesh_generation.C.
{
libmesh_error_msg("Building a circle/sphere only works with AMR.");
}
| 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 1468 of file mesh_generation.C.
References build_cube().
{
// This method only makes sense in 2D!
// But we now just turn a non-2D mesh into a 2D mesh
//libmesh_assert_equal_to (mesh.mesh_dimension(), 2);
// Call the build_cube() member to actually do the work for us.
build_cube (mesh,
nx, ny, 0,
xmin, xmax,
ymin, ymax,
0., 0.,
type,
gauss_lobatto_grid);
}
generated by
