diff --git a/src/mesh/mesh_generation.C b/src/mesh/mesh_generation.C index 2525efd46e..f95ca49f0f 100644 --- a/src/mesh/mesh_generation.C +++ b/src/mesh/mesh_generation.C @@ -30,6 +30,7 @@ #include "libmesh/face_quad4.h" #include "libmesh/face_quad8.h" #include "libmesh/face_quad9.h" +#include "libmesh/face_c0polygon.h" #include "libmesh/cell_hex8.h" #include "libmesh/cell_hex20.h" #include "libmesh/cell_hex27.h" @@ -614,6 +615,12 @@ void MeshTools::Generation::build_cube(UnstructuredMesh & mesh, break; } + case C0POLYGON: + { + mesh.reserve_elem ((nx + 1) * (ny + 1)); + break; + } + default: libmesh_error_msg("ERROR: Unrecognized 2D element type == " << Utility::enum_to_string(type)); } @@ -649,6 +656,11 @@ void MeshTools::Generation::build_cube(UnstructuredMesh & mesh, mesh.reserve_nodes( (2*nx+1)*(2*ny+1) + 2*nx*ny ); break; } + case C0POLYGON: + { + mesh.reserve_nodes (4 + 3*nx*ny + 2*nx + 2*ny); + break; + } default: libmesh_error_msg("ERROR: Unrecognized 2D element type == " << Utility::enum_to_string(type)); @@ -735,6 +747,11 @@ void MeshTools::Generation::build_cube(UnstructuredMesh & mesh, break; } + case C0POLYGON: + { + // we create the nodes at the same time as the elements + break; + } default: libmesh_error_msg("ERROR: Unrecognized 2D element type == " << Utility::enum_to_string(type)); @@ -896,6 +913,195 @@ void MeshTools::Generation::build_cube(UnstructuredMesh & mesh, break; }; + case C0POLYGON: + { + // Build a 2D paving using hexagons (center), quads (part of y-boundary) + // and triangles (x-boundaries). + // Vector to re-use previously created nodes + std::vector node_list; + + // Start with a layer of triangles on the boundary + const auto dx_tri = (xmax - xmin) / (nx); + const auto dy_tri = (ymax - ymin) / (ny + 1); + std::unique_ptr new_elem; + for (const auto i : make_range(nx + 1)) + { + // Make new nodes for bottom layer of triangles + Node *node0, *node1, *node2; + if (i == 0) + { + node0 = mesh.add_point(Point(0., 0, 0.)); + node1 = mesh.add_point(Point(0., dy_tri / 2., 0.)); + node2 = mesh.add_point(Point(dx_tri / 2., 0., 0.)); + node_list.push_back(node0); + node_list.push_back(node1); + node_list.push_back(node2); + } + else if (i < nx) + { + node0 = node_list.back(); + node1 = mesh.add_point(Point((i)*dx_tri, dy_tri / 2., 0.)); + node2 = mesh.add_point(Point((i + 1. / 2.) * dx_tri, 0., 0.)); + node_list.push_back(node1); + node_list.push_back(node2); + } + else + { + node0 = node_list.back(); + node1 = mesh.add_point(Point((i)*dx_tri, dy_tri / 2., 0.)); + node2 = mesh.add_point(Point((i)*dx_tri, 0., 0.)); + node_list.push_back(node1); + node_list.push_back(node2); + } + + new_elem = std::make_unique(3); + // Switch to Tri3 when exodus default output supports element type mixes + new_elem->set_node(0, node0); + new_elem->set_node(1, node1); + new_elem->set_node(2, node2); + auto * elem = mesh.add_elem(std::move(new_elem)); + + // Set boundaries + if (i == 0) + boundary_info.add_side(elem, 0, 3); // left + else if (i == nx) + boundary_info.add_side(elem, 1, 1); // right + boundary_info.add_side(elem, 2, 0); // bottom + } + // Start with the second node to build hexagons + unsigned int running_index = 1; + + // Build layers of hexagons + const auto hex_side = + (ymax - ymin - (ny == 1 ? dy_tri : (1. + (ny - 1) / 2.) * dy_tri)) / ny; + for (const auto j : make_range(ny)) + { + for (const auto i : make_range(nx + (j % 2))) + { + if ((j % 2 == 0) || ((i > 0) && (i < nx))) + { + Node *n0, *n1, *n2, *n3, *n4, *n5; + n0 = node_list[running_index++]; + n1 = node_list[running_index++]; + n2 = node_list[running_index]; + + if (i == 0) + { + n3 = mesh.add_point(Point(*n0) + RealVectorValue(0, hex_side, 0)); + node_list.push_back(n3); + } + else + n3 = node_list.back(); + + n4 = mesh.add_point(Point(*n1) + RealVectorValue(0, hex_side + dy_tri, 0)); + n5 = mesh.add_point(Point(*n2) + RealVectorValue(0, hex_side, 0)); + node_list.push_back(n4); + node_list.push_back(n5); + + new_elem = std::make_unique(6); + new_elem->set_node(0, n0); + new_elem->set_node(1, n1); + new_elem->set_node(2, n2); + new_elem->set_node(3, n5); + new_elem->set_node(4, n4); + new_elem->set_node(5, n3); + auto * elem = mesh.add_elem(std::move(new_elem)); + + // Set boundaries + if (i == 0) + boundary_info.add_side(elem, 5, 3); // left + else if (i == nx) + boundary_info.add_side(elem, 2, 1); // right + } + // The hexagons are offset, so we build on a quad on each external side to fill + else if (i == 0 || i == nx) + { + Node *n0, *n1, *n2, *n3; + n0 = node_list[running_index++]; + n1 = node_list[running_index]; + + if (i == 0) + { + n2 = mesh.add_point(Point(*n0) + RealVectorValue(0, hex_side + dy_tri, 0)); + node_list.push_back(n2); + n3 = mesh.add_point(Point(*n1) + RealVectorValue(0, hex_side, 0)); + } + else + { + n2 = node_list.back(); + n3 = mesh.add_point(Point(*n1) + RealVectorValue(0, hex_side + dy_tri, 0)); + } + node_list.push_back(n3); + + new_elem = std::make_unique(4); + // Switch to Quad4 when exodus default output supports element type mixes + new_elem->set_node(0, n0); + new_elem->set_node(1, n1); + new_elem->set_node(3, n2); + new_elem->set_node(2, n3); + auto * elem = mesh.add_elem(std::move(new_elem)); + + // Set boundaries + if (i == 0) + boundary_info.add_side(elem, 3, 3); // left + else if (i == nx) + boundary_info.add_side(elem, 1, 1); // right + } + else + libmesh_assert(false); + } + // Increment once to switch to next 'row' of nodes + running_index++; + + // Skip lower right corner node + if (j == 0) + running_index++; + } + + // Build a final layer of triangles + const bool ny_odd = (ny % 2 == 1); + for (const auto i : make_range(nx + ny_odd)) + { + // Use existing nodes, except at the corners + Node *node0, *node1, *node2; + if (i == 0 && ny_odd) + { + node0 = mesh.add_point(Point(0., ymax, 0.)); + node1 = node_list[running_index++]; + node2 = node_list[running_index]; + } + else if (i < nx) + { + node0 = node_list[running_index++]; + node1 = node_list[running_index++]; + node2 = node_list[running_index]; + } + // This case only reached if ny is odd and we are using a triangle in top right corner + else + { + node0 = node_list[running_index++]; + node1 = node_list[running_index]; + node2 = mesh.add_point(Point(xmax, ymax, 0.)); + } + + new_elem = std::make_unique(3); + // Switch to Tri3 when exodus default output supports element type mixes + new_elem->set_node(0, node0); + new_elem->set_node(1, node1); + new_elem->set_node(2, node2); + auto * elem = mesh.add_elem(std::move(new_elem)); + + // Set boundaries + if (i == 0) + boundary_info.add_side(elem, 0, 3); // left + else if (i == nx) + boundary_info.add_side(elem, 1, 1); // right + boundary_info.add_side(elem, 2, 2); // top + + } + break; + } + default: libmesh_error_msg("ERROR: Unrecognized 2D element type == " << Utility::enum_to_string(type)); diff --git a/tests/mesh/mesh_generation_test.C b/tests/mesh/mesh_generation_test.C index 987fedd71d..56c232f755 100644 --- a/tests/mesh/mesh_generation_test.C +++ b/tests/mesh/mesh_generation_test.C @@ -37,6 +37,8 @@ public: CPPUNIT_TEST( buildSquareQuad4 ); CPPUNIT_TEST( buildSquareQuad8 ); CPPUNIT_TEST( buildSquareQuad9 ); + CPPUNIT_TEST( buildSquareC0PolygonEven ); + CPPUNIT_TEST( buildSquareC0PolygonOdd ); # ifdef LIBMESH_ENABLE_AMR CPPUNIT_TEST( buildSphereTri3 ); CPPUNIT_TEST( buildSphereQuad4 ); @@ -107,7 +109,9 @@ public: void testBuildSquare(UnstructuredMesh & mesh, unsigned int n, ElemType type) { MeshTools::Generation::build_square (mesh, n, n, -2.0, 3.0, -4.0, 5.0, type); - if (Elem::type_to_n_sides_map[type] == 4) + if (type == C0POLYGON) + CPPUNIT_ASSERT_EQUAL(mesh.n_elem(), cast_int(n*n + 4 + 2 * (n - 1) + ((n - 1) / 2))); + else if (Elem::type_to_n_sides_map[type] == 4) CPPUNIT_ASSERT_EQUAL(mesh.n_elem(), cast_int(n*n)); else CPPUNIT_ASSERT_EQUAL(mesh.n_elem(), cast_int(n*n*2)); @@ -132,6 +136,10 @@ public: CPPUNIT_ASSERT_EQUAL(mesh.n_nodes(), cast_int((2*n+1)*(2*n+1) + 2*n*n)); break; + case C0POLYGON: + CPPUNIT_ASSERT_EQUAL(mesh.n_nodes(), + cast_int(4 + 2*n*n + (n - 1) + 2*n + 2 * (n%2))); + break; default: // Wait, what did we try to build? CPPUNIT_ASSERT(false); } @@ -146,8 +154,9 @@ public: // Do serial assertions *after* all parallel assertions, so we // stay in sync after failure on only some processor(s) - for (auto & elem : mesh.element_ptr_range()) - CPPUNIT_ASSERT(elem->has_affine_map()); + if (type != C0POLYGON) + for (auto & elem : mesh.element_ptr_range()) + CPPUNIT_ASSERT(elem->has_affine_map()); } void testBuildCube(UnstructuredMesh & mesh, unsigned int n, ElemType type) @@ -293,6 +302,8 @@ public: void buildSquareQuad4 () { LOG_UNIT_TEST; tester(&MeshGenerationTest::testBuildSquare, 4, QUAD4); } void buildSquareQuad8 () { LOG_UNIT_TEST; tester(&MeshGenerationTest::testBuildSquare, 4, QUAD8); } void buildSquareQuad9 () { LOG_UNIT_TEST; tester(&MeshGenerationTest::testBuildSquare, 4, QUAD9); } + void buildSquareC0PolygonOdd() { LOG_UNIT_TEST; tester(&MeshGenerationTest::testBuildSquare, 5, C0POLYGON); } + void buildSquareC0PolygonEven(){ LOG_UNIT_TEST; tester(&MeshGenerationTest::testBuildSquare, 6, C0POLYGON); } void buildSphereTri3 () { LOG_UNIT_TEST; testBuildSphere(2, TRI3); } void buildSphereQuad4 () { LOG_UNIT_TEST; testBuildSphere(2, QUAD4); }