diff --git a/include/geom/cell_c0polyhedron.h b/include/geom/cell_c0polyhedron.h index e4d762232a..41ba0d8d3e 100644 --- a/include/geom/cell_c0polyhedron.h +++ b/include/geom/cell_c0polyhedron.h @@ -170,6 +170,14 @@ class C0Polyhedron : public Polyhedron Point side_vertex_average_normal(const unsigned int s) const override final; + /** + * \returns the local side-indices of subelement sides of the + * polyhedron + * + * Each subelement here is a tetrahedron + */ + virtual std::array subelement_sides_to_poly_sides(unsigned int i) const; + /** * Create a triangulation (tetrahedralization) based on the current * sides' triangulations. diff --git a/src/geom/cell_c0polyhedron.C b/src/geom/cell_c0polyhedron.C index 58be64ae42..08e7d09f41 100644 --- a/src/geom/cell_c0polyhedron.C +++ b/src/geom/cell_c0polyhedron.C @@ -270,6 +270,39 @@ C0Polyhedron::side_vertex_average_normal(const unsigned int s) const } +std::array +C0Polyhedron::subelement_sides_to_poly_sides(unsigned int tri_i) const +{ + std::array sides = {invalid_int, invalid_int, invalid_int, invalid_int}; + const auto & tet_nodes = _triangulation[tri_i]; + for (const auto poly_side : make_range(n_sides())) + { + const auto sd_nodes = nodes_on_side(poly_side); + bool zero_in = std::find(sd_nodes.begin(), sd_nodes.end(), tet_nodes[0]) != sd_nodes.end(); + bool one_in = std::find(sd_nodes.begin(), sd_nodes.end(), tet_nodes[1]) != sd_nodes.end(); + bool two_in = std::find(sd_nodes.begin(), sd_nodes.end(), tet_nodes[2]) != sd_nodes.end(); + bool three_in = std::find(sd_nodes.begin(), sd_nodes.end(), tet_nodes[3]) != sd_nodes.end(); + + // Take advantage of the fact that the poly side can only contain one tet side + // Note: we could optimize by replacing the booleans with the searches above + if (zero_in) + { + if (one_in) + { + if (two_in) + sides[0] = poly_side; + else if (three_in) + sides[1] = poly_side; + } + else if (two_in && three_in) + sides[3] = poly_side; + } + else if (three_in && two_in && one_in) + sides[2] = poly_side; + } + return sides; +} + void C0Polyhedron::retriangulate() { @@ -373,7 +406,7 @@ void C0Polyhedron::retriangulate() // sometimes with a similarly-sorted vector of surrounding elements // to go with it. auto surroundings_of = - [&nodes_to_elem_map, & surface] + [& nodes_to_elem_map, & surface] (const Node & node, std::vector * surrounding_elems) { @@ -439,7 +472,11 @@ void C0Polyhedron::retriangulate() v02 = static_cast(*surrounding_nodes[n+1]) - node, v03 = static_cast(*surrounding_nodes[n+2]) - node; - total_solid_angle += solid_angle(v01, v02, v03); + // Solid angles are all negative because the surface triangulation triangles were pointing + // outwards and the 4th node is inward because the polyhedron is convex + // Some are zero when all the nodes selected are from the same (planar) polygon + // Make it positive to store smallest angle first + total_solid_angle += std::abs(solid_angle(v01, v02, v03)); } return std::make_pair(n_surrounding, total_solid_angle); @@ -456,8 +493,11 @@ void C0Polyhedron::retriangulate() // in parallel. typedef std::multimap, Node*> node_map_type; node_map_type nodes_by_geometry; + // Keep track of index in map std::map node_index; + // The map is first sorted by the valence + // then by the angle. Max angle, if they are negative (positive right now), first for (auto node : surface.node_ptr_range()) node_index[node] = nodes_by_geometry.emplace(geometry_at(*node), node); @@ -466,6 +506,9 @@ void C0Polyhedron::retriangulate() // each vertex, and an inner loop to remove multiple tetrahedra in // cases where the vertex has more than 3 neighboring triangles. + // One of the checks does not work if we skip adding a 0 volume tet + bool has_skipped_adding_tets_and_retriangulating = false; + // We'll be done when there are only three "unremoved" nodes left, // so they don't actually enclose any volume. for (auto i : make_range(nodes_by_geometry.size()-3)) @@ -500,19 +543,56 @@ void C0Polyhedron::retriangulate() // predecessor and successor nodes in order afterward. auto find_valid_nodes_around = [n_surrounding, & surrounding_nodes] - (unsigned int j) + (unsigned int j, unsigned int skip=1e6) { unsigned int jnext = (j+1)%n_surrounding; - while (!surrounding_nodes[jnext]) + while (!surrounding_nodes[jnext] || jnext == skip) jnext = (jnext+1)%n_surrounding; unsigned int jprev = (j+n_surrounding-1)%n_surrounding; - while (!surrounding_nodes[jprev]) + while (!surrounding_nodes[jprev] || jprev == skip) jprev = (jprev+n_surrounding-1)%n_surrounding; return std::make_pair(jprev, jnext); }; + + // Vectors from the center node to each of its surrounding + // nodes are helpful for calculating prospective tet + // quality. + std::vector v0s(n_surrounding); + for (auto j : make_range(n_surrounding)) + v0s[j] = *(Point *)surrounding_nodes[j] - *node; + + // Find the tet quality we'd potentially get from each + // possible choice of tet + auto local_tet_quality_of = + [& surrounding_nodes, & v0s, & find_valid_nodes_around] + (unsigned int j, unsigned int skip=1e6) + { + auto [jminus, jplus] = find_valid_nodes_around(j, skip); + + // Anything proportional to the ratio of volume to + // total-edge-length-cubed should peak for perfect tets + // but hit 0 for pancakes and slivers. + + const Real total_len = + v0s[j].norm() + v0s[jminus].norm() + v0s[jplus].norm() + + (*(Point *)surrounding_nodes[jplus] - + *(Point *)surrounding_nodes[j]).norm() + + (*(Point *)surrounding_nodes[j] - + *(Point *)surrounding_nodes[jminus]).norm() + + (*(Point *)surrounding_nodes[jminus] - + *(Point *)surrounding_nodes[jplus]).norm(); + + // Orientation here is tricky. Think of the triple + // product as (v1 cross v2) dot v3, with right hand rule. + const Real six_vol = + triple_product(v0s[jminus], v0s[jplus], v0s[j]); + + return six_vol / (total_len * total_len * total_len); + }; + // We may have too many surrounding nodes to handle with // just one tet. In that case we'll keep a cache of the // element qualities that we'd get by making a tet with the @@ -525,8 +605,9 @@ void C0Polyhedron::retriangulate() std::vector local_tet_quality(n_surrounding, 1); // From our center node with N surrounding nodes we can make N-2 - // tetrahedra. The first N-3 each replace two surface tets with - // two new surface tets. + // tetrahedra. + // For the first N-3 tets we build, we may have to replace two surface + // triangles with two new surface triangles. // // My first idea was to greedily pick nodes with the smallest // local (solid) angles to get the best quality. This works in @@ -539,15 +620,39 @@ void C0Polyhedron::retriangulate() // My third idea is to try to fix the lowest quality tets first, // by picking cases where they have higher quality neighbors, // and creating those neighbors so as to change them. + // + // A fourth idea was to make sure to fix as many neighbor tets as + // possible. auto find_new_tet_nodes = - [& local_tet_quality, & find_valid_nodes_around] + [& local_tet_quality, & local_tet_quality_of, & find_valid_nodes_around, & surrounding_nodes] () { unsigned int jbest = 0; auto [jminus, jplus] = find_valid_nodes_around(jbest); Real qneighbest = std::min(local_tet_quality[jminus], local_tet_quality[jplus]); + // Count number of bad neighbors + auto num_bad_neigh_best = (local_tet_quality[jminus] <= 0) + (local_tet_quality[jplus] <= 0); + // Skip our heuristics if we are considering a bad tet + if (local_tet_quality[jbest] == -1e6) + { + num_bad_neigh_best = 0; + qneighbest = 1e6; + } + + // Count the number of valid surrounding nodes + unsigned int n_valid_surrounding = 0; + for (const auto surr_node : surrounding_nodes) + if (surr_node) + n_valid_surrounding++; + + bool bad_future_neighbors = false; + if ((n_valid_surrounding > 3) && local_tet_quality_of(jminus, /*without*/jbest) <= 0 && + local_tet_quality_of(jplus, /*without*/jbest) <= 0) + bad_future_neighbors = true; + + // Expected qualities post removing jbest from the surrounding nodes for (auto j : make_range(std::size_t(1), local_tet_quality.size())) { @@ -565,17 +670,35 @@ void C0Polyhedron::retriangulate() qneighj > 0) continue; + // Avoid chosing a tet that when constructed would leave both neighbors in bad shape + // TODO: simply pre-compute that for every vertex and exclude those nodes + // Note: It does not matter if we are down to 3 surrounding nodes, as we are building the final tet + if ((n_valid_surrounding > 3) && local_tet_quality_of(jminus, /*without*/j) <= 0 && + local_tet_quality_of(jplus, /*without*/j) <= 0) + continue; + + const auto num_bad_neigh = (local_tet_quality[jminus] <= 0) + + (local_tet_quality[jplus] <= 0); // We want to try for the best possible fix. - if ((local_tet_quality[j] - qneighj) > - (local_tet_quality[jbest] - qneighj)) + // - The more bad (volume <= 0) neighbors we fix the better + // - If same number of bad neighbors, pick best quality one + // - If the current pick for best tet will make the two neighbor tets bad, + // don't keep it and always try another + if (((num_bad_neigh > num_bad_neigh_best) || ((num_bad_neigh == num_bad_neigh_best) + && ((local_tet_quality[j] - qneighj) > + (local_tet_quality[jbest] - qneighj)))) + || bad_future_neighbors) { jbest = j; qneighbest = qneighj; + num_bad_neigh_best = num_bad_neigh; + bad_future_neighbors = false; } } + // If there are only 3 neighors and the tet is bad, we will skip it downstream libmesh_error_msg_if - (local_tet_quality[jbest] <= 0, + (n_valid_surrounding > 3 && local_tet_quality[jbest] <= 0, "Cannot build non-singular non-inverted tet"); std::tie(jminus, jplus) = find_valid_nodes_around(jbest); @@ -583,53 +706,22 @@ void C0Polyhedron::retriangulate() return std::make_tuple(jbest, jminus, jplus); }; - if (n_surrounding > 3) - { - // We'll be searching local_tet_quality even after tet - // extraction disconnects us from some nodes; when we do we - // don't want to get one. - constexpr Real far_node = -1e6; - - // Vectors from the center node to each of its surrounding - // nodes are helpful for calculating prospective tet - // quality. - std::vector v0s(n_surrounding); - for (auto j : make_range(n_surrounding)) - v0s[j] = *(Point *)surrounding_nodes[j] - *node; - - // Find the tet quality we'd potentially get from each - // possible choice of tet - auto local_tet_quality_of = - [& surrounding_nodes, & v0s, & find_valid_nodes_around] - (unsigned int j) - { - auto [jminus, jplus] = find_valid_nodes_around(j); - - // Anything proportional to the ratio of volume to - // total-edge-length-cubed should peak for perfect tets - // but hit 0 for pancakes and slivers. - - const Real total_len = - v0s[j].norm() + v0s[jminus].norm() + v0s[jplus].norm() + - (*(Point *)surrounding_nodes[jplus] - - *(Point *)surrounding_nodes[j]).norm() + - (*(Point *)surrounding_nodes[j] - - *(Point *)surrounding_nodes[jminus]).norm() + - (*(Point *)surrounding_nodes[jminus] - - *(Point *)surrounding_nodes[jplus]).norm(); - - // Orientation here is tricky. Think of the triple - // product as (v1 cross v2) dot v3, with right hand rule. - const Real six_vol = - triple_product(v0s[jminus], v0s[jplus], v0s[j]); + // Get the quality of the tets for every node around the node: + // we have chosen 1 node. If we chose another node for the tet, there + // are only two other nodes that are also neighbors, hence 1 tet per node + // Note: we always need to know the quality, because as we delete vertices + // we could end up with a a node with 3 surrounding nodes, but all coplanar! + for (auto j : make_range(n_surrounding)) + local_tet_quality[j] = local_tet_quality_of(j); - return six_vol / (total_len * total_len * total_len); - }; + // We'll be searching local_tet_quality even after tet + // extraction disconnects us from some nodes; when we do we + // don't want to re-use it so we mark the quality as -1e6 + constexpr Real far_node = -1e6; - for (auto j : make_range(n_surrounding)) - local_tet_quality[j] = local_tet_quality_of(j); - - // If we have N surrounding nodes, we can make N tets and + if (n_surrounding > 3) + { + // If we have N surrounding nodes, we can make N-3 tets and // that'll bring us back to the 3-surrounding-node case to // finish. for (auto t : make_range(n_surrounding-3)) @@ -697,6 +789,7 @@ void C0Polyhedron::retriangulate() // unchanged until we're out of this inner loop; let's // recalculate it here, and then we'll be done with it. Node * & nbestref = surrounding_nodes[jbest]; + libmesh_assert(nbestref); nodes_by_geometry.erase(node_index[nbestref]); node_index[nbestref] = nodes_by_geometry.emplace(geometry_at(*nbestref), nbestref); @@ -724,6 +817,18 @@ void C0Polyhedron::retriangulate() // counterclockwise" we get from the lambda. auto [j2, j1, j3] = find_new_tet_nodes(); + // If the final tet has zero volume, skip, we have covered all we need for this node + // This occurs when considering a polygonal prism. + // The opposite node on the other extruded face was deleted first + // and depending on the triangulation, this leaves the node only attached + // to 3 other coplanar nodes on one of the side + if (local_tet_quality[j2] == 0 || local_tet_quality[j2] == far_node) + { + nodes_by_geometry.erase(geometry_it); + has_skipped_adding_tets_and_retriangulating = true; + continue; + } + // Turn these last four nodes into a tet Node * n1 = surrounding_nodes[j1], * n2 = surrounding_nodes[j2], @@ -767,6 +872,26 @@ void C0Polyhedron::retriangulate() surface.delete_elem(oldtri2); surface.delete_elem(oldtri3); + // We've used up our center node, so it's not something we can + // eliminate again. + nodes_by_geometry.erase(geometry_it); + + // Recompute the valence and angles of the nodes we used + // The idea is that if one node is being isolated on one side, + // we need to treat it asap + if (nodes_by_geometry.size() > 4) + { + Node * & node_j2 = surrounding_nodes[j2]; + nodes_by_geometry.erase(node_index[node_j2]); + node_index[node_j2] = nodes_by_geometry.emplace(geometry_at(*node_j2), node_j2); + Node * & node_j1 = surrounding_nodes[j1]; + nodes_by_geometry.erase(node_index[node_j1]); + node_index[node_j1] = nodes_by_geometry.emplace(geometry_at(*node_j1), node_j1); + Node * & node_j3 = surrounding_nodes[j3]; + nodes_by_geometry.erase(node_index[node_j3]); + node_index[node_j3] = nodes_by_geometry.emplace(geometry_at(*node_j3), node_j3); + } + // We should have used up all our surrounding nodes now, and we // shouldn't have messed up our surface in the process, and our // center node should no longer be part of the surface. @@ -786,14 +911,11 @@ void C0Polyhedron::retriangulate() libmesh_assert_not_equal_to (elem->node_ptr(p), node); #endif - - // We've used up our center node, so it's not something we can - // eliminate again. - nodes_by_geometry.erase(geometry_it); } // At this point our surface should just have two triangles left. - libmesh_assert_equal_to(surface.n_elem(), 2); + if (!has_skipped_adding_tets_and_retriangulating) + libmesh_assert_equal_to(surface.n_elem(), 2); } diff --git a/src/geom/cell_polyhedron.C b/src/geom/cell_polyhedron.C index 23e4a4c0a4..35106f9fd5 100644 --- a/src/geom/cell_polyhedron.C +++ b/src/geom/cell_polyhedron.C @@ -115,6 +115,7 @@ Polyhedron::Polyhedron (const std::vector> & sides, { const Polygon & side = *sides[s]; const Point x_i = side.point(0); + // opposite of normal const Point n_i = (side.point(1) - side.point(0)).cross (side.point(0) - side.point(side.n_sides()-1)).unit(); @@ -123,7 +124,7 @@ Polyhedron::Polyhedron (const std::vector> & sides, inward_normal = (n_i * (center - x_i) > TOLERANCE); } - // We're betting a lot on "our polyhedra are all convex", so let's + // We're betting a lot on "our polygon sides are all convex", so let's // check that if we have time. #ifdef DEBUG for (unsigned int s : index_range(sides)) diff --git a/src/geom/face_c0polygon.C b/src/geom/face_c0polygon.C index 0437e51074..955af2133e 100644 --- a/src/geom/face_c0polygon.C +++ b/src/geom/face_c0polygon.C @@ -396,11 +396,16 @@ void C0Polygon::retriangulate() { Real min_cos_angle = 1; int best_vertex = -1; + std::cout << "Remaining nodes : "; + for (auto n : index_range(remaining_nodes)) + std::cout << remaining_nodes[n] << " "; + std::cout << std::endl; + const auto max_n = remaining_nodes.size(); for (auto n : index_range(remaining_nodes)) { const Point & pn = this->point(remaining_nodes[n]); - const Point & pnext = this->point(remaining_nodes[(n+1)%ns]); - const Point & pprev = this->point(remaining_nodes[(n+ns-1)%ns]); + const Point & pnext = this->point(remaining_nodes[(n+1)%max_n]); + const Point & pprev = this->point(remaining_nodes[(n+max_n-1)%max_n]); const Point vprev = (pn - pprev).unit(); const Point vnext = (pnext - pn).unit(); @@ -418,9 +423,9 @@ void C0Polygon::retriangulate() libmesh_assert(best_vertex >= 0); - this->_triangulation.push_back({remaining_nodes[(best_vertex+ns-1)%ns], + this->_triangulation.push_back({remaining_nodes[(best_vertex+max_n-1)%max_n], remaining_nodes[best_vertex], - remaining_nodes[(best_vertex+1)%ns]}); + remaining_nodes[(best_vertex+1)%max_n]}); remaining_nodes.erase(remaining_nodes.begin()+best_vertex); } } 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); }