diff --git a/docs/src/config/model.md b/docs/src/config/model.md index 7b3c59ccd..7a674e804 100644 --- a/docs/src/config/model.md +++ b/docs/src/config/model.md @@ -115,6 +115,7 @@ mesh length units. - `"RemoveCurvature" [false]` - `"MakeSimplex" [false]` + - `"MakeHex" [false]` - `"ReorderElements" [false]` - `"CleanUnusedElements" [true]` - `"AddInterfaceBoundaryElements" [true]` diff --git a/palace/utils/configfile.cpp b/palace/utils/configfile.cpp index 1e489d3a9..72e8486a5 100644 --- a/palace/utils/configfile.cpp +++ b/palace/utils/configfile.cpp @@ -464,6 +464,7 @@ void ModelData::SetUp(json &config) Lc = model->value("Lc", Lc); remove_curvature = model->value("RemoveCurvature", remove_curvature); make_simplex = model->value("MakeSimplex", make_simplex); + make_hex = model->value("MakeHex", make_hex); reorder_elements = model->value("ReorderElements", reorder_elements); clean_unused_elements = model->value("CleanUnusedElements", clean_unused_elements); add_bdr_elements = model->value("AddInterfaceBoundaryElements", add_bdr_elements); @@ -477,6 +478,7 @@ void ModelData::SetUp(json &config) model->erase("Lc"); model->erase("RemoveCurvature"); model->erase("MakeSimplex"); + model->erase("MakeHex"); model->erase("ReorderElements"); model->erase("CleanUnusedElements"); model->erase("AddInterfaceBoundaryElements"); @@ -495,6 +497,7 @@ void ModelData::SetUp(json &config) std::cout << "Lc: " << Lc << '\n'; std::cout << "RemoveCurvature: " << remove_curvature << '\n'; std::cout << "MakeSimplex: " << make_simplex << '\n'; + std::cout << "MakeHex: " << make_hex << '\n'; std::cout << "ReorderElements: " << reorder_elements << '\n'; std::cout << "CleanUnusedElements: " << clean_unused_elements << '\n'; std::cout << "AddInterfaceBoundaryElements: " << add_bdr_elements << '\n'; diff --git a/palace/utils/configfile.hpp b/palace/utils/configfile.hpp index 194967933..3e65af1cc 100644 --- a/palace/utils/configfile.hpp +++ b/palace/utils/configfile.hpp @@ -205,6 +205,9 @@ struct ModelData // Convert mesh to simplex elements. bool make_simplex = false; + // Convert mesh to hexahedral elements (using tet-to-hex algorithm). + bool make_hex = false; + // Reorder elements based on spatial location after loading the serial mesh, which can // potentially increase memory coherency. bool reorder_elements = false; diff --git a/palace/utils/geodata.cpp b/palace/utils/geodata.cpp index e59fd37b2..baecd8a3c 100644 --- a/palace/utils/geodata.cpp +++ b/palace/utils/geodata.cpp @@ -40,14 +40,14 @@ constexpr auto MSH_FLT_PRECISION = std::numeric_limits::max_digits10; // Load the serial mesh from disk. std::unique_ptr LoadMesh(const std::string &, bool); +// Create a new mesh by splitting all elements of the mesh into simplices or hexes +// (using tet-to-hex). Optionally preserves curvature of the original mesh by interpolating +// the high-order nodes with GSLIB. +void SplitMeshElements(mfem::Mesh &, bool, bool, bool = true); + // Optionally reorder mesh elements based on MFEM's internal reordeing tools for improved // cache usage. -void ReorderMesh(mfem::Mesh &, bool = true); - -// Create a new mesh by splitting all elements of the mesh into simplices. Optionally -// preserves curvature of the original mesh by interpolating the high-order nodes with -// GSLIB. -void MakeSimplicial(mfem::Mesh &, bool = true); +void ReorderMeshElements(mfem::Mesh &, bool = true); // Clean the provided serial mesh by removing unnecessary domain and elements, and adding // boundary elements for material interfaces and exterior boundaries. @@ -110,17 +110,17 @@ std::unique_ptr ReadMesh(MPI_Comm comm, const IoData &iodata) "If there are pyramid elements, AMR must be nonconformal!"); // Optionally convert mesh elements to simplices, for example in order to enable - // conformal mesh refinement. - if (iodata.model.make_simplex) + // conformal mesh refinement, or hexes. + if (iodata.model.make_simplex || iodata.model.make_hex) { - MakeSimplicial(*smesh); + SplitMeshElements(*smesh, iodata.model.make_simplex, iodata.model.make_hex); } // Optionally reorder elements (and vertices) based on spatial location after loading // the serial mesh. if (iodata.model.reorder_elements) { - ReorderMesh(*smesh); + ReorderMeshElements(*smesh); } // Clean up unused domain elements from the mesh, add new boundary elements for material @@ -1578,53 +1578,174 @@ std::unique_ptr LoadMesh(const std::string &mesh_file, bool remove_c return mesh; } -void ReorderMesh(mfem::Mesh &mesh, bool print) +mfem::Mesh MeshTetToHex(const mfem::Mesh &orig_mesh) { - mfem::Array ordering; + // Courtesy of https://gist.github.com/pazner/e9376f77055c0918d7c43e034e9e5888, only + // supports tetrahedral elements for now. Eventually should be expanded to support prism + // and pyramid elements but this mixed mesh support requires a bit more work. + MFEM_VERIFY(orig_mesh.Dimension() == 3, "Tet-to-hex conversion only supports 3D meshes!"); + { + mfem::Array geoms; + orig_mesh.GetGeometries(3, geoms); + MFEM_VERIFY(geoms.Size() == 1 && geoms[0] == mfem::Geometry::TETRAHEDRON, + "Tet-to-hex conversion only works for pure tetrahedral meshes!"); + } - if constexpr (false) + // Add new vertices in every edge, face, and volume. Each tet is subdivided into 4 hexes, + // and each triangular face subdivided into 3 quads. + const int nv_tet = orig_mesh.GetNV(); + const int nedge_tet = orig_mesh.GetNEdges(); + const int nface_tet = orig_mesh.GetNFaces(); + const int ne_tet = orig_mesh.GetNE(); + const int nbe_tet = orig_mesh.GetNBE(); + const int nv = nv_tet + nedge_tet + nface_tet + ne_tet; + const int ne = 4 * ne_tet; + const int nbe = 3 * nbe_tet; + mfem::Mesh hex_mesh(orig_mesh.Dimension(), nv, ne, nbe, orig_mesh.SpaceDimension()); + + // Add original vertices. + for (int v = 0; v < nv_tet; v++) { - // Gecko reordering. - mfem::Array tentative; - int outer = 3, inner = 3, window = 4, period = 2; - double best_cost = mfem::infinity(); - for (int i = 0; i < outer; i++) + hex_mesh.AddVertex(orig_mesh.GetVertex(v)); + } + + // Add midpoints of edges, faces, and elements. + auto AddCentroid = [&orig_mesh, &hex_mesh](const int *verts, int nv) + { + double coord[3] = {0.0, 0.0, 0.0}; + for (int i = 0; i < nv; i++) { - int seed = i + 1; - double cost = - mesh.GetGeckoElementOrdering(tentative, inner, window, period, seed, true); - if (cost < best_cost) + for (int d = 0; d < orig_mesh.SpaceDimension(); d++) { - ordering = tentative; - best_cost = cost; + coord[d] += orig_mesh.GetVertex(verts[i])[d] / nv; } } - if (print) + hex_mesh.AddVertex(coord); + }; + { + mfem::Array verts; + for (int e = 0; e < nedge_tet; ++e) { - Mpi::Print("Final cost: {:e}\n", best_cost); + orig_mesh.GetEdgeVertices(e, verts); + AddCentroid(verts.GetData(), verts.Size()); } } + for (int f = 0; f < nface_tet; ++f) + { + AddCentroid(orig_mesh.GetFace(f)->GetVertices(), orig_mesh.GetFace(f)->GetNVertices()); + } + for (int e = 0; e < ne_tet; ++e) + { + AddCentroid(orig_mesh.GetElement(e)->GetVertices(), + orig_mesh.GetElement(e)->GetNVertices()); + } - // (Faster) Hilbert reordering. - mesh.GetHilbertElementOrdering(ordering); - mesh.ReorderElements(ordering); + // Connectivity of tetrahedron vertices to the edges. + constexpr int tet_vertex_edge_map[4 * 3] = {0, 1, 2, 3, 0, 4, 1, 3, 5, 5, 4, 2}; + constexpr int tet_vertex_face_map[4 * 3] = {3, 2, 1, 3, 0, 2, 3, 1, 0, 0, 1, 2}; + constexpr int tri_vertex_edge_map[3 * 2] = {0, 2, 1, 0, 2, 1}; + + // Add four hexahedra for each tetrahedron. + { + mfem::Array edges, faces, orients; + for (int e = 0; e < ne_tet; ++e) + { + const int *verts = orig_mesh.GetElement(e)->GetVertices(); + orig_mesh.GetElementEdges(e, edges, orients); + orig_mesh.GetElementFaces(e, faces, orients); + + // One hex for each vertex of the tet. + for (int i = 0; i < 4; ++i) + { + int hex_v[8]; + hex_v[0] = verts[i]; + hex_v[1] = nv_tet + edges[tet_vertex_edge_map[3 * i + 0]]; + hex_v[2] = nv_tet + nedge_tet + faces[tet_vertex_face_map[3 * i + 0]]; + hex_v[3] = nv_tet + edges[tet_vertex_edge_map[3 * i + 1]]; + hex_v[4] = nv_tet + edges[tet_vertex_edge_map[3 * i + 2]]; + hex_v[5] = nv_tet + nedge_tet + faces[tet_vertex_face_map[3 * i + 1]]; + hex_v[6] = nv_tet + nedge_tet + nface_tet + e; + hex_v[7] = nv_tet + nedge_tet + faces[tet_vertex_face_map[3 * i + 2]]; + hex_mesh.AddHex(hex_v, orig_mesh.GetAttribute(e)); + } + } + } + + // Add the boundary elements. + { + mfem::Array edges, orients; + for (int be = 0; be < nbe_tet; ++be) + { + int f, o; + const int *verts = orig_mesh.GetBdrElement(be)->GetVertices(); + orig_mesh.GetBdrElementEdges(be, edges, orients); + orig_mesh.GetBdrElementFace(be, &f, &o); + + // One quad for each vertex of the tri. + for (int i = 0; i < 3; ++i) + { + int quad_v[4]; + quad_v[0] = verts[i]; + quad_v[1] = nv_tet + edges[tri_vertex_edge_map[2 * i + 0]]; + quad_v[2] = nv_tet + nedge_tet + f; + quad_v[3] = nv_tet + edges[tri_vertex_edge_map[2 * i + 1]]; + hex_mesh.AddBdrQuad(quad_v, orig_mesh.GetBdrAttribute(be)); + } + } + } + + // Finalize the hex mesh. No need to generate boundary elements or mark for refinement, + // but we fix orientations for the new elements as needed. + constexpr bool generate_bdr = false, refine = false, fix_orientation = true; + hex_mesh.FinalizeTopology(generate_bdr); + hex_mesh.Finalize(refine, fix_orientation); + return hex_mesh; } -void MakeSimplicial(mfem::Mesh &orig_mesh, bool preserve_curvature) +void SplitMeshElements(mfem::Mesh &orig_mesh, bool make_simplex, bool make_hex, + bool preserve_curvature) { - // Convert all element types to simplices. - const auto element_types = mesh::CheckElements(orig_mesh); - if (element_types.has_simplices && !element_types.has_hexahedra && - !element_types.has_prisms && !element_types.has_pyramids) + if (!make_simplex && !make_hex) { return; } - MFEM_VERIFY(!element_types.has_pyramids, - "mfem::Mesh::MakeSimplicial does not support pyramid elements yet!"); - mfem::Mesh new_mesh = mfem::Mesh::MakeSimplicial(orig_mesh); + MFEM_VERIFY(!orig_mesh.Nonconforming(), + "Mesh element splitting is not supported for nonconforming meshes!"); + mfem::Mesh *mesh = &orig_mesh; + mfem::Mesh new_mesh; + + // Convert all element types to simplices. + if (make_simplex) + { + const auto element_types = mesh::CheckElements(*mesh); + if (element_types.has_hexahedra || element_types.has_prisms || + element_types.has_pyramids) + { + MFEM_VERIFY( + !element_types.has_pyramids, + "Splitting mesh elements to simplices does not support pyramid elements yet!"); + new_mesh = mfem::Mesh::MakeSimplicial(*mesh); + mesh = &new_mesh; + } + } - // MFEM's function removes curvature information from the new mesh. So, if needed, we - // interpolate it onto the new mesh with GSLIB. + // Convert all element types to hexahedra (currently only tet-to-hex). + if (make_hex) + { + const auto element_types = mesh::CheckElements(*mesh); + if (element_types.has_simplices || element_types.has_prisms || + element_types.has_pyramids) + { + MFEM_VERIFY(!element_types.has_prisms && !element_types.has_pyramids, + "Splitting mesh elements to hexahedra only supports simplex elements " + "(tetrahedra) for now!"); + new_mesh = MeshTetToHex(*mesh); + mesh = &new_mesh; + } + } + + // The previous splitting functions remove curvature information from the new mesh. So, if + // needed, we interpolate it onto the new mesh with GSLIB. if (preserve_curvature && orig_mesh.GetNodes()) { // Back up the original high-order nodes. @@ -1659,6 +1780,38 @@ void MakeSimplicial(mfem::Mesh &orig_mesh, bool preserve_curvature) orig_mesh = std::move(new_mesh); } +void ReorderMeshElements(mfem::Mesh &mesh, bool print) +{ + mfem::Array ordering; + + if constexpr (false) + { + // Gecko reordering. + mfem::Array tentative; + int outer = 3, inner = 3, window = 4, period = 2; + double best_cost = mfem::infinity(); + for (int i = 0; i < outer; i++) + { + int seed = i + 1; + double cost = + mesh.GetGeckoElementOrdering(tentative, inner, window, period, seed, true); + if (cost < best_cost) + { + ordering = tentative; + best_cost = cost; + } + } + if (print) + { + Mpi::Print("Final cost: {:e}\n", best_cost); + } + } + + // (Faster) Hilbert reordering. + mesh.GetHilbertElementOrdering(ordering); + mesh.ReorderElements(ordering); +} + std::map> CheckMesh(std::unique_ptr &orig_mesh, const IoData &iodata, bool clean_elem, bool add_bdr_elem) diff --git a/scripts/schema/config/model.json b/scripts/schema/config/model.json index 4de2435a3..556464dcd 100644 --- a/scripts/schema/config/model.json +++ b/scripts/schema/config/model.json @@ -11,6 +11,7 @@ "Lc": { "type": "number", "exclusiveMinimum": 0.0 }, "RemoveCurvature": { "type": "boolean" }, "MakeSimplex": { "type": "boolean" }, + "MakeHex": { "type": "boolean" }, "ReorderElements": { "type": "boolean" }, "CleanUnusedElements": { "type": "boolean" }, "AddInterfaceBoundaryElements": { "type": "boolean" },