From 71cc19e94fb1eac1c25025fc24c2ac15ae56c354 Mon Sep 17 00:00:00 2001 From: Sebastian Grimberg Date: Wed, 22 May 2024 14:45:30 -0700 Subject: [PATCH 01/14] Preliminary cleanup for mesh cracking (remove unused subdomain boundary element addition) --- palace/main.cpp | 2 +- palace/utils/geodata.cpp | 219 +++++++++++++++------------------------ palace/utils/geodata.hpp | 5 +- 3 files changed, 85 insertions(+), 141 deletions(-) diff --git a/palace/main.cpp b/palace/main.cpp index ac3cfebf2..f76af96da 100644 --- a/palace/main.cpp +++ b/palace/main.cpp @@ -308,7 +308,7 @@ int main(int argc, char *argv[]) std::vector> mesh; { std::vector> mfem_mesh; - mfem_mesh.push_back(mesh::ReadMesh(world_comm, iodata, false, true, true, false)); + mfem_mesh.push_back(mesh::ReadMesh(world_comm, iodata, false, true, true)); iodata.NondimensionalizeInputs(*mfem_mesh[0]); mesh::RefineMesh(iodata, mfem_mesh); for (auto &m : mfem_mesh) diff --git a/palace/utils/geodata.cpp b/palace/utils/geodata.cpp index 487ef4c80..a4bbfbfaf 100644 --- a/palace/utils/geodata.cpp +++ b/palace/utils/geodata.cpp @@ -9,9 +9,10 @@ #include #include #include -#include #include #include +#include +#include #include #include "utils/communication.hpp" #include "utils/diagnostic.hpp" @@ -46,11 +47,10 @@ void ReorderMesh(mfem::Mesh &, bool = true); std::unique_ptr GetMeshPartitioning(const mfem::Mesh &, int, const std::string & = "", bool = true); -// Clean the provided serial mesh by removing unnecessary domain and elements, adding -// boundary elements for material interfaces and exterior boundaries, and adding boundary -// elements for subdomain interfaces. +// Clean the provided serial mesh by removing unnecessary domain and elements, and adding +// boundary elements for material interfaces and exterior boundaries. std::map> CheckMesh(std::unique_ptr &, const IoData &, - bool, bool, bool, const int *); + bool, bool); // Given a serial mesh on the root processor and element partitioning, create a parallel // mesh over the given communicator. The serial mesh is destroyed when no longer needed. @@ -66,8 +66,9 @@ void RebalanceConformalMesh(std::unique_ptr &); namespace mesh { -std::unique_ptr ReadMesh(MPI_Comm comm, const IoData &iodata, bool reorder, - bool clean_elem, bool add_bdr, bool add_subdomain) +std::unique_ptr ReadMesh(MPI_Comm comm, const IoData &iodata, + bool reorder_elem, bool clean_elem, + bool add_bdr_elem) { // If possible on root, read the serial mesh (converting format if necessary), and do all // necessary serial preprocessing. When finished, distribute the mesh to all processes. @@ -103,7 +104,7 @@ std::unique_ptr ReadMesh(MPI_Comm comm, const IoData &iodata, boo // Optionally reorder elements (and vertices) based on spatial location after loading // the serial mesh. - if (reorder) + if (reorder_elem) { ReorderMesh(*smesh); } @@ -111,20 +112,13 @@ std::unique_ptr ReadMesh(MPI_Comm comm, const IoData &iodata, boo // Clean up unused domain elements from the mesh, add new boundary elements for material // interfaces if not present. Can only clean up conforming meshes, assumes that any // nonconformal mesh was generated by adaptation and thus does not need checking. - if ((clean_elem || add_bdr) && smesh->Conforming()) + if ((clean_elem || add_bdr_elem) && smesh->Conforming()) { - static_cast(CheckMesh(smesh, iodata, clean_elem, add_bdr, false, nullptr)); + static_cast(CheckMesh(smesh, iodata, clean_elem, add_bdr_elem)); } // Generate the mesh partitioning. partitioning = GetMeshPartitioning(*smesh, Mpi::Size(comm), iodata.model.partition); - - // Optionally add subdomain interface boundary elements. - if (add_subdomain && smesh->Conforming()) - { - static_cast( - CheckMesh(smesh, iodata, false, false, add_subdomain, partitioning.get())); - } } std::unique_ptr pmesh; @@ -739,7 +733,7 @@ int CollectPointCloudOnRoot(const mfem::ParMesh &mesh, const mfem::Array &m // Linear mesh, work with element vertices directly. PalacePragmaOmp(parallel) { - std::set vertex_indices; + std::unordered_set vertex_indices; mfem::Array v; if (bdr) { @@ -1769,8 +1763,7 @@ std::unique_ptr GetMeshPartitioning(const mfem::Mesh &mesh, int size, std::map> CheckMesh(std::unique_ptr &orig_mesh, const IoData &iodata, bool clean_elem, - bool add_bdr, bool add_subdomain, - const int *partitioning) + bool add_bdr_elem) { // - Check that all external boundaries of the mesh have a corresponding boundary // condition. @@ -1790,7 +1783,7 @@ std::map> CheckMesh(std::unique_ptr &orig_me orig_mesh->bdr_attributes.Size() ? orig_mesh->bdr_attributes.Max() : 0, iodata.boundaries.attributes, true); { - std::set bdr_warn_list; + std::unordered_set bdr_warn_list; for (int be = 0; be < orig_mesh->GetNBE(); be++) { int attr = orig_mesh->GetBdrAttribute(be); @@ -1822,30 +1815,25 @@ std::map> CheckMesh(std::unique_ptr &orig_me // Mapping from new interface boundary attribute tags to vector of neighboring domain // attributes (when adding new boundary elements). - if (!clean_elem && !add_bdr && !add_subdomain) + if (!clean_elem && !add_bdr_elem) { return {}; } // Count deleted or added domain and boundary elements. int new_ne = orig_mesh->GetNE(); - int new_nbdr = orig_mesh->GetNBE(); - mfem::Array elem_delete, bdr_delete; - mfem::Array orig_bdr_faces, add_bdr_faces; - elem_delete.SetSize(orig_mesh->GetNE(), false); - bdr_delete.SetSize(orig_mesh->GetNBE(), false); - orig_bdr_faces.SetSize(orig_mesh->GetNumFaces(), -1); + int new_nbe = orig_mesh->GetNBE(); + std::vector elem_delete(orig_mesh->GetNE(), false), + bdr_elem_delete(orig_mesh->GetNBE(), false); + std::unordered_map face_to_be, new_face_be; + face_to_be.reserve(orig_mesh->GetNBE()); for (int be = 0; be < orig_mesh->GetNBE(); be++) { int f, o; orig_mesh->GetBdrElementFace(be, &f, &o); - MFEM_VERIFY(orig_bdr_faces[f] < 0, + MFEM_VERIFY(face_to_be.find(f) == face_to_be.end(), "Mesh should not define boundary elements multiple times!"); - orig_bdr_faces[f] = be; - } - if (add_bdr || add_subdomain) - { - add_bdr_faces.SetSize(orig_mesh->GetNumFaces(), -1); + face_to_be[f] = be; } if (clean_elem) @@ -1864,21 +1852,18 @@ std::map> CheckMesh(std::unique_ptr &orig_me // Make sure to remove any boundary elements which are no longer attached to elements in // the domain. - for (int f = 0; f < orig_mesh->GetNumFaces(); f++) + for (int be = 0; be < orig_mesh->GetNBE(); be++) { - const int &be = orig_bdr_faces[f]; - if (be >= 0) + int f, o, e1, e2; + orig_mesh->GetBdrElementFace(be, &f, &o); + orig_mesh->GetFaceElements(f, &e1, &e2); + bool no_e1 = (e1 < 0 || elem_delete[e1]); + bool no_e2 = (e2 < 0 || elem_delete[e2]); + if (no_e1 && no_e2) { - int e1, e2; - orig_mesh->GetFaceElements(f, &e1, &e2); - bool no_e1 = (e1 < 0 || elem_delete[e1]); - bool no_e2 = (e2 < 0 || elem_delete[e2]); - if (no_e1 && no_e2) - { - // Mpi::Print("Deleting an unattached boundary element!\n"); - bdr_delete[be] = true; - new_nbdr--; - } + // Mpi::Print("Deleting an unattached boundary element!\n"); + bdr_elem_delete[be] = true; + new_nbe--; } } if (new_ne < orig_mesh->GetNE()) @@ -1886,110 +1871,71 @@ std::map> CheckMesh(std::unique_ptr &orig_me Mpi::Print("Removed {:d} unmarked domain elements from the mesh\n", orig_mesh->GetNE() - new_ne); } - if (new_nbdr < orig_mesh->GetNBE()) + if (new_nbe < orig_mesh->GetNBE()) { Mpi::Print("Removed {:d} unattached boundary elements from the mesh\n", - orig_mesh->GetNBE() - new_nbdr); + orig_mesh->GetNBE() - new_nbe); } } - int new_ne_step1 = new_ne; - int new_nbdr_step1 = new_nbdr; + int new_ne_step_1 = new_ne; + int new_nbe_step_1 = new_nbe; - if (add_bdr) + if (add_bdr_elem) { // Add new boundary elements at material interfaces or on the exterior boundary of the // simulation domain, if there is not already a boundary element present. MFEM_VERIFY(!orig_mesh->Nonconforming(), "Adding material interface boundary elements " "is not supported for nonconforming meshes!"); - int add_bdr_ext = 0, add_bdr_int = 0; + int new_nbe_ext = 0, new_nbe_int = 0; for (int f = 0; f < orig_mesh->GetNumFaces(); f++) { - const int &be = orig_bdr_faces[f]; - if (be < 0 && add_bdr_faces[f] < 0) + if (face_to_be.find(f) != face_to_be.end()) { - int e1, e2; - orig_mesh->GetFaceElements(f, &e1, &e2); - bool no_e1 = (e1 < 0 || elem_delete[e1]); - bool no_e2 = (e2 < 0 || elem_delete[e2]); - if ((no_e1 || no_e2) && !(no_e1 && no_e2)) - { - // Mpi::Print("Adding exterior boundary element!\n"); - add_bdr_faces[f] = 1; - add_bdr_ext++; - } - else if (orig_mesh->GetAttribute(e1) != orig_mesh->GetAttribute(e2)) - { - // Add new boundary element at material interface between two domains. - // Mpi::Print("Adding material interface boundary element!\n"); - add_bdr_faces[f] = 1; - add_bdr_int++; - } + continue; + } + int e1, e2; + orig_mesh->GetFaceElements(f, &e1, &e2); + bool no_e1 = (e1 < 0 || elem_delete[e1]); + bool no_e2 = (e2 < 0 || elem_delete[e2]); + if ((no_e1 || no_e2) && !(no_e1 && no_e2)) + { + // Mpi::Print("Adding exterior boundary element!\n"); + new_face_be[f] = 1; + new_nbe_ext++; + } + else if (orig_mesh->GetAttribute(e1) != orig_mesh->GetAttribute(e2)) + { + // Add new boundary element at material interface between two domains. + // Mpi::Print("Adding material interface boundary element!\n"); + new_face_be[f] = 1; + new_nbe_int++; } } - new_nbdr += (add_bdr_ext + add_bdr_int); - if (add_bdr_ext > 0) + new_nbe += (new_nbe_ext + new_nbe_int); + if (new_nbe_ext > 0) { Mpi::Print("Added {:d} boundary elements for exterior boundaries to the mesh\n", - add_bdr_ext); + new_nbe_ext); } - if (add_bdr_int > 0) + if (new_nbe_int > 0) { Mpi::Print("Added {:d} boundary elements for material interfaces to the mesh\n", - add_bdr_int); - } - } - int new_ne_step2 = new_ne; - int new_nbdr_step2 = new_nbdr; - - if (add_subdomain) - { - // Add new boundary elements at interfaces between elements belonging to different - // subdomains. This uses similar code to mfem::Mesh::PrintWithPartitioning. - MFEM_VERIFY(partitioning, "Cannot add subdomain interface boundary elements without " - "supplied mesh partitioning!"); - MFEM_VERIFY(!orig_mesh->Nonconforming(), "Adding subdomain interface boundary elements " - "is not supported for nonconforming meshes!"); - for (int f = 0; f < orig_mesh->GetNumFaces(); f++) - { - const int &be = orig_bdr_faces[f]; - if (be < 0 && add_bdr_faces[f] < 0) - { - int e1, e2; - orig_mesh->GetFaceElements(f, &e1, &e2); - bool no_e1 = (e1 < 0 || elem_delete[e1]); - bool no_e2 = (e2 < 0 || elem_delete[e2]); - if (!no_e1 && !no_e2 && partitioning[e1] != partitioning[e2]) - { - // Internal face is connected to two elements belonging to different subdomains - // (this works for conforming meshes). - add_bdr_faces[f] = 2; - new_nbdr += 2; - } - } - // else - // { - // // This face is attached to a boundary element. We could define a new boundary - // // element with opposite orientation to ensure both subdomains in the distributed - // // ParMesh have the boundary element. - // } - } - if (new_nbdr > new_nbdr_step2) - { - Mpi::Print("Added {:d} boundary elements for subdomain interfaces to the mesh\n", - new_nbdr - new_nbdr_step2); + new_nbe_int); } } + int new_ne_step_2 = new_ne; + int new_nbe_step_2 = new_nbe; // Create the new mesh. - if (new_ne == new_ne_step1 && new_ne_step1 == new_ne_step2 && - new_ne_step2 == orig_mesh->GetNE() && new_nbdr == new_nbdr_step1 && - new_nbdr_step1 == new_nbdr_step2 && new_nbdr_step2 == orig_mesh->GetNBE()) + if (new_ne == new_ne_step_1 && new_ne_step_1 == new_ne_step_2 && + new_ne_step_2 == orig_mesh->GetNE() && new_nbe == new_nbe_step_1 && + new_nbe_step_1 == new_nbe_step_2 && new_nbe_step_2 == orig_mesh->GetNBE()) { return {}; } auto new_mesh = std::make_unique(orig_mesh->Dimension(), orig_mesh->GetNV(), new_ne, - new_nbdr, orig_mesh->SpaceDimension()); + new_nbe, orig_mesh->SpaceDimension()); std::map> new_attr_map; // Copy vertices and non-deleted domain and boundary elements. @@ -2007,7 +1953,7 @@ std::map> CheckMesh(std::unique_ptr &orig_me } for (int be = 0; be < orig_mesh->GetNBE(); be++) { - if (!bdr_delete[be]) + if (!bdr_elem_delete[be]) { mfem::Element *el = orig_mesh->GetBdrElement(be)->Duplicate(new_mesh.get()); new_mesh->AddBdrElement(el); @@ -2015,16 +1961,8 @@ std::map> CheckMesh(std::unique_ptr &orig_me } // Add new boundary elements. - if (add_bdr || add_subdomain) + if (add_bdr_elem) { - auto FlipVertices = [](mfem::Element *el) - { - mfem::Array v; - el->GetVertices(v); - std::reverse(v.begin(), v.end()); - el->SetVertices(v.HostRead()); - }; - // Some (1-based) boundary attributes may be empty since they were removed from the // original mesh, but to keep attributes the same as config file we don't compress the // list. @@ -2032,7 +1970,7 @@ std::map> CheckMesh(std::unique_ptr &orig_me orig_mesh->bdr_attributes.Size() ? orig_mesh->bdr_attributes.Max() : 0; for (int f = 0; f < orig_mesh->GetNumFaces(); f++) { - if (add_bdr_faces[f] > 0) + if (new_face_be[f] > 0) { // Assign new unique attribute based on attached elements. Save so that the // attributes of e1 and e2 can be easily referenced using the new attribute. Since @@ -2070,18 +2008,23 @@ std::map> CheckMesh(std::unique_ptr &orig_me mfem::Element *el = orig_mesh->GetFace(f)->Duplicate(new_mesh.get()); el->SetAttribute(new_attr); new_mesh->AddBdrElement(el); - if (add_bdr_faces[f] > 1) + if (new_face_be[f] > 1) { // Flip order of vertices to reverse normal direction of second added element. el = orig_mesh->GetFace(f)->Duplicate(new_mesh.get()); - FlipVertices(el); el->SetAttribute(new_attr); + { + mfem::Array v; + el->GetVertices(v); + std::reverse(v.begin(), v.end()); + el->SetVertices(v.HostRead()); + }; new_mesh->AddBdrElement(el); if constexpr (false) { - Mpi::Print( - "Adding two boundary elements with attr {:d} from elements {:d} and {:d}\n", - new_attr, a, b); + Mpi::Print("Adding two boundary elements with attribute {:d} from elements " + "{:d} and {:d}\n", + new_attr, a, b); } } } diff --git a/palace/utils/geodata.hpp b/palace/utils/geodata.hpp index 9c9ef3c38..862c372bd 100644 --- a/palace/utils/geodata.hpp +++ b/palace/utils/geodata.hpp @@ -24,8 +24,9 @@ namespace mesh // Read and partition a serial mesh from file, returning a pointer to the new parallel mesh // object, which should be destroyed by the user. -std::unique_ptr ReadMesh(MPI_Comm comm, const IoData &iodata, bool reorder, - bool clean_elem, bool add_bdr, bool add_subdomain); +std::unique_ptr ReadMesh(MPI_Comm comm, const IoData &iodata, + bool reorder_elem, bool clean_elem, + bool add_bdr_elem); // Refine the provided mesh according to the data in the input file. If levels of refinement // are requested, the refined meshes are stored in order of increased refinement. Ownership From 651af79ec4b14c6b0c7498f8536b4f6b483e0b34 Mon Sep 17 00:00:00 2001 From: Sebastian Grimberg Date: Fri, 24 May 2024 13:51:45 -0700 Subject: [PATCH 02/14] Add configuration file options for mesh reordering, mesh cleaning, and interface boundary element additions (rather than hard-coding) --- docs/src/config/model.md | 16 +++++++---- palace/main.cpp | 2 +- palace/models/postoperator.cpp | 2 +- palace/utils/configfile.cpp | 21 ++++++++++---- palace/utils/configfile.hpp | 22 +++++++++++---- palace/utils/geodata.cpp | 47 ++++++++++++++++---------------- palace/utils/geodata.hpp | 4 +-- palace/utils/iodata.cpp | 40 +++++++++++++-------------- palace/utils/iodata.hpp | 2 +- scripts/schema/config/model.json | 7 +++-- 10 files changed, 94 insertions(+), 69 deletions(-) diff --git a/docs/src/config/model.md b/docs/src/config/model.md index 1c3c3078d..f31daddd3 100644 --- a/docs/src/config/model.md +++ b/docs/src/config/model.md @@ -22,11 +22,14 @@ with `"Mesh" [None]` : Input mesh file path, an absolute path is recommended. -`"L0" [1.0e-6]` : Mesh vertex coordinate length unit, m. +`"L0" [1.0e-6]` : Unit, relative to m, for mesh vertex coordinates. For example, a value +of `1.0e-6` implies the mesh coordinates are in μm. `"Lc" [0.0]` : Characteristic length scale used for nondimensionalization, specified in -mesh length units. A value less than or equal to zero uses an internally calculated length -scale based on the bounding box of the computational domain. +mesh length units. This keyword should typically not be specified by the user. A value less +than or equal to zero uses an internally calculated length scale based on the bounding box +of the computational domain. A value of 1.0 will disable nondimensionalization of lengths +in the model and all computations will take place in the same units as the mesh. `"Refinement"` : Top-level object for configuring mesh refinement. @@ -110,9 +113,12 @@ mesh length units. ### Advanced model options - - `"Partition" [""]` - - `"ReorientTetMesh" [false]` - `"RemoveCurvature" [false]` + - `"ReorderElements" [false]` + - `"CleanUnusedElements" [true]` + - `"AddInterfaceBoundaryElements" [true]` + - `"ReorientTetMesh" [false]` + - `"Partitioning" [""]` - `"MaxNCLevels" [1]` - `"MaximumImbalance" [1.1]` - `"SaveAdaptIterations" [true]` diff --git a/palace/main.cpp b/palace/main.cpp index f76af96da..3e0c671d1 100644 --- a/palace/main.cpp +++ b/palace/main.cpp @@ -308,7 +308,7 @@ int main(int argc, char *argv[]) std::vector> mesh; { std::vector> mfem_mesh; - mfem_mesh.push_back(mesh::ReadMesh(world_comm, iodata, false, true, true)); + mfem_mesh.push_back(mesh::ReadMesh(world_comm, iodata)); iodata.NondimensionalizeInputs(*mfem_mesh[0]); mesh::RefineMesh(iodata, mfem_mesh); for (auto &m : mfem_mesh) diff --git a/palace/models/postoperator.cpp b/palace/models/postoperator.cpp index a1a906d79..c5c5e5538 100644 --- a/palace/models/postoperator.cpp +++ b/palace/models/postoperator.cpp @@ -160,7 +160,7 @@ void PostOperator::InitializeDataCollection(const IoData &iodata) const bool use_ho = true; const int refine_ho = HasE() ? E->ParFESpace()->GetMaxElementOrder() : B->ParFESpace()->GetMaxElementOrder(); - mesh_Lc0 = iodata.GetLengthScale(); + mesh_Lc0 = iodata.GetMeshLengthScale(); // Output mesh coordinate units same as input. paraview.SetCycle(-1); diff --git a/palace/utils/configfile.cpp b/palace/utils/configfile.cpp index 6b1fbc32c..7db0cff3d 100644 --- a/palace/utils/configfile.cpp +++ b/palace/utils/configfile.cpp @@ -462,18 +462,24 @@ void ModelData::SetUp(json &config) mesh = model->at("Mesh"); // Required L0 = model->value("L0", L0); Lc = model->value("Lc", Lc); - partition = model->value("Partition", partition); - reorient_tet = model->value("ReorientTetMesh", reorient_tet); remove_curvature = model->value("RemoveCurvature", remove_curvature); + 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); + reorient_tet_mesh = model->value("ReorientTetMesh", reorient_tet_mesh); + partitioning = model->value("Partitioning", partitioning); refinement.SetUp(*model); // Cleanup model->erase("Mesh"); model->erase("L0"); model->erase("Lc"); - model->erase("Partition"); - model->erase("ReorientTetMesh"); model->erase("RemoveCurvature"); + model->erase("ReorderElements"); + model->erase("CleanUnusedElements"); + model->erase("AddInterfaceBoundaryElements"); + model->erase("ReorientTetMesh"); + model->erase("Partitioning"); model->erase("Refinement"); MFEM_VERIFY(model->empty(), "Found an unsupported configuration file keyword under \"Model\"!\n" @@ -485,9 +491,12 @@ void ModelData::SetUp(json &config) std::cout << "Mesh: " << mesh << '\n'; std::cout << "L0: " << L0 << '\n'; std::cout << "Lc: " << Lc << '\n'; - std::cout << "Partition: " << partition << '\n'; - std::cout << "ReorientTetMesh: " << reorient_tet << '\n'; std::cout << "RemoveCurvature: " << remove_curvature << '\n'; + std::cout << "ReorderElements: " << reorder_elements << '\n'; + std::cout << "CleanUnusedElements: " << clean_unused_elements << '\n'; + std::cout << "AddInterfaceBoundaryElements: " << add_bdr_elements << '\n'; + std::cout << "ReorientTetMesh: " << reorient_tet_mesh << '\n'; + std::cout << "Partitioning: " << partitioning << '\n'; } } diff --git a/palace/utils/configfile.hpp b/palace/utils/configfile.hpp index 0014dc8ef..19e529996 100644 --- a/palace/utils/configfile.hpp +++ b/palace/utils/configfile.hpp @@ -199,14 +199,26 @@ struct ModelData double L0 = 1.0e-6; double Lc = -1.0; - // Partitioning file (if specified, does not compute a new partitioning). - std::string partition = ""; + // Remove high-order curvature information from the mesh. + bool remove_curvature = false; + + // Reorder elements based on spatial location after loading the serial mesh, which can + // potentially increase memory coherency. + bool reorder_elements = false; + + // Remove elements (along with any associated unattached boundary elements) from the mesh + // which do not have any material properties specified. + bool clean_unused_elements = true; + + // Add new boundary elements for faces are on the computational domain boundary or which + // have attached elements on either side with different domain attributes. + bool add_bdr_elements = true; // Call MFEM's ReorientTetMesh as a check of mesh orientation after partitioning. - bool reorient_tet = false; + bool reorient_tet_mesh = false; - // Remove high-order curvature information from the mesh. - bool remove_curvature = false; + // Partitioning file (if specified, does not compute a new partitioning). + std::string partitioning = ""; // Object controlling mesh refinement. RefinementData refinement = {}; diff --git a/palace/utils/geodata.cpp b/palace/utils/geodata.cpp index a4bbfbfaf..809a8feb8 100644 --- a/palace/utils/geodata.cpp +++ b/palace/utils/geodata.cpp @@ -66,9 +66,7 @@ void RebalanceConformalMesh(std::unique_ptr &); namespace mesh { -std::unique_ptr ReadMesh(MPI_Comm comm, const IoData &iodata, - bool reorder_elem, bool clean_elem, - bool add_bdr_elem) +std::unique_ptr ReadMesh(MPI_Comm comm, const IoData &iodata) { // If possible on root, read the serial mesh (converting format if necessary), and do all // necessary serial preprocessing. When finished, distribute the mesh to all processes. @@ -104,7 +102,7 @@ std::unique_ptr ReadMesh(MPI_Comm comm, const IoData &iodata, // Optionally reorder elements (and vertices) based on spatial location after loading // the serial mesh. - if (reorder_elem) + if (iodata.model.reorder_elements) { ReorderMesh(*smesh); } @@ -112,13 +110,14 @@ std::unique_ptr ReadMesh(MPI_Comm comm, const IoData &iodata, // Clean up unused domain elements from the mesh, add new boundary elements for material // interfaces if not present. Can only clean up conforming meshes, assumes that any // nonconformal mesh was generated by adaptation and thus does not need checking. - if ((clean_elem || add_bdr_elem) && smesh->Conforming()) + if (smesh->Conforming()) { - static_cast(CheckMesh(smesh, iodata, clean_elem, add_bdr_elem)); + static_cast(CheckMesh(smesh, iodata, iodata.model.clean_unused_elements, + iodata.model.add_bdr_elements)); } // Generate the mesh partitioning. - partitioning = GetMeshPartitioning(*smesh, Mpi::Size(comm), iodata.model.partition); + partitioning = GetMeshPartitioning(*smesh, Mpi::Size(comm), iodata.model.partitioning); } std::unique_ptr pmesh; @@ -206,7 +205,7 @@ void RefineMesh(const IoData &iodata, std::vector // reorient only the coarse mesh so that the refinements are still true refinements of // the original mesh (required for geometric multigrid). Otherwise, it happens after // refinement. - if (iodata.model.reorient_tet && mesh.capacity() > 1) + if (iodata.model.reorient_tet_mesh && mesh.capacity() > 1) { PalacePragmaDiagnosticPush PalacePragmaDiagnosticDisableDeprecated @@ -354,7 +353,7 @@ void RefineMesh(const IoData &iodata, std::vector // not required after MFEM's PR #1046, but in case you want to be absolutely sure, we // reorient only the mesh after refinement if there is a single mesh (doesn't work with // h-refinement geometric multigrid). - if (iodata.model.reorient_tet && mesh.size() == 1) + if (iodata.model.reorient_tet_mesh && mesh.capacity() == 1) { PalacePragmaDiagnosticPush PalacePragmaDiagnosticDisableDeprecated @@ -1524,7 +1523,7 @@ double RebalanceMesh(std::unique_ptr &mesh, const IoData &iodata) // fo << std::fixed; fo << std::scientific; fo.precision(MSH_FLT_PRECISION); - mesh::DimensionalizeMesh(smesh, iodata.GetLengthScale()); + mesh::DimensionalizeMesh(smesh, iodata.GetMeshLengthScale()); smesh.Mesh::Print(fo); // Do not need to nondimensionalize the temporary mesh } Mpi::Barrier(comm); @@ -1588,7 +1587,7 @@ template void AttrToMarker(int, const std::vector &, mfem::Array &, bo namespace { -std::unique_ptr LoadMesh(const std::string &path, bool remove_curvature) +std::unique_ptr LoadMesh(const std::string &mesh_file, bool remove_curvature) { // Read the (serial) mesh from the given mesh file. Handle preparation for refinement and // orientations here to avoid possible reorientations and reordering later on. MFEM @@ -1597,9 +1596,9 @@ std::unique_ptr LoadMesh(const std::string &path, bool remove_curvat // or error out if not supported. constexpr bool generate_edges = true, refine = true, fix_orientation = true; std::unique_ptr mesh; - std::filesystem::path mfile(path); - if (mfile.extension() == ".mphtxt" || mfile.extension() == ".mphbin" || - mfile.extension() == ".nas" || mfile.extension() == ".bdf") + std::filesystem::path mesh_path(mesh_file); + if (mesh_path.extension() == ".mphtxt" || mesh_path.extension() == ".mphbin" || + mesh_path.extension() == ".nas" || mesh_path.extension() == ".bdf") { // Put translated mesh in temporary string buffer. std::stringstream fi(std::stringstream::in | std::stringstream::out); @@ -1623,15 +1622,15 @@ std::unique_ptr LoadMesh(const std::string &path, bool remove_curvat fo.precision(MSH_FLT_PRECISION); #endif - if (mfile.extension() == ".mphtxt" || mfile.extension() == ".mphbin") + if (mesh_path.extension() == ".mphtxt" || mesh_path.extension() == ".mphbin") { - mesh::ConvertMeshComsol(path, fi, remove_curvature); - // mesh::ConvertMeshComsol(path, fo, remove_curvature); + mesh::ConvertMeshComsol(mesh_file, fi, remove_curvature); + // mesh::ConvertMeshComsol(mesh_file, fo, remove_curvature); } else { - mesh::ConvertMeshNastran(path, fi, remove_curvature); - // mesh::ConvertMeshNastran(path, fo, remove_curvature); + mesh::ConvertMeshNastran(mesh_file, fi, remove_curvature); + // mesh::ConvertMeshNastran(mesh_file, fo, remove_curvature); } #if 0 @@ -1648,10 +1647,10 @@ std::unique_ptr LoadMesh(const std::string &path, bool remove_curvat else { // Otherwise, just rely on MFEM load the mesh. - std::ifstream fi(path); + std::ifstream fi(mesh_file); if (!fi.good()) { - MFEM_ABORT("Unable to open mesh file \"" << path << "\"!"); + MFEM_ABORT("Unable to open mesh file \"" << mesh_file << "\"!"); } mesh = std::make_unique(fi, generate_edges, refine, fix_orientation); } @@ -1708,11 +1707,11 @@ void ReorderMesh(mfem::Mesh &mesh, bool print) } std::unique_ptr GetMeshPartitioning(const mfem::Mesh &mesh, int size, - const std::string &partition, bool print) + const std::string &part_file, bool print) { MFEM_VERIFY(size <= mesh.GetNE(), "Mesh partitioning must have parts <= mesh elements (" << size << " vs. " << mesh.GetNE() << ")!"); - if (partition.length() == 0) + if (part_file.length() == 0) { const int part_method = 1; std::unique_ptr partitioning( @@ -1734,7 +1733,7 @@ std::unique_ptr GetMeshPartitioning(const mfem::Mesh &mesh, int size, // // int ne, np; - std::ifstream part_ifs(partition); + std::ifstream part_ifs(part_file); part_ifs.ignore(std::numeric_limits::max(), ' '); part_ifs >> ne; if (ne != mesh.GetNE()) diff --git a/palace/utils/geodata.hpp b/palace/utils/geodata.hpp index 862c372bd..285ebd52f 100644 --- a/palace/utils/geodata.hpp +++ b/palace/utils/geodata.hpp @@ -24,9 +24,7 @@ namespace mesh // Read and partition a serial mesh from file, returning a pointer to the new parallel mesh // object, which should be destroyed by the user. -std::unique_ptr ReadMesh(MPI_Comm comm, const IoData &iodata, - bool reorder_elem, bool clean_elem, - bool add_bdr_elem); +std::unique_ptr ReadMesh(MPI_Comm comm, const IoData &iodata); // Refine the provided mesh according to the data in the input file. If levels of refinement // are requested, the refined meshes are stored in order of increased refinement. Ownership diff --git a/palace/utils/iodata.cpp b/palace/utils/iodata.cpp index 4409abf69..f2c5b6818 100644 --- a/palace/utils/iodata.cpp +++ b/palace/utils/iodata.cpp @@ -460,55 +460,53 @@ void IoData::NondimensionalizeInputs(mfem::ParMesh &mesh) MFEM_VERIFY(!init, "NondimensionalizeInputs should only be called once!"); init = true; - // Calculate the reference length and time. - if (model.Lc > 0.0) - { - // User specified Lc in mesh length units. - Lc = model.Lc * model.L0; // [m] - } - else + // Calculate the reference length and time. A user specified model.Lc is in mesh length + // units. + if (model.Lc <= 0.0) { mfem::Vector bbmin, bbmax; mesh::GetAxisAlignedBoundingBox(mesh, bbmin, bbmax); bbmax -= bbmin; - bbmax *= model.L0; // [m] - Lc = *std::max_element(bbmax.begin(), bbmax.end()); + model.Lc = *std::max_element(bbmax.begin(), bbmax.end()); } + Lc = model.Lc * model.L0; // [m] tc = 1.0e9 * Lc / electromagnetics::c0_; // [ns] // Mesh refinement parameters. - auto DivideLength = [this](double val) { return val / (Lc / model.L0); }; + auto DivideLengthScale = [Lc0 = GetMeshLengthScale()](double val) { return val / Lc0; }; for (auto &box : model.refinement.GetBoxes()) { - std::transform(box.bbmin.begin(), box.bbmin.end(), box.bbmin.begin(), DivideLength); - std::transform(box.bbmax.begin(), box.bbmax.end(), box.bbmax.begin(), DivideLength); + std::transform(box.bbmin.begin(), box.bbmin.end(), box.bbmin.begin(), + DivideLengthScale); + std::transform(box.bbmax.begin(), box.bbmax.end(), box.bbmax.begin(), + DivideLengthScale); } for (auto &sphere : model.refinement.GetSpheres()) { - sphere.r /= Lc / model.L0; + sphere.r /= GetMeshLengthScale(); std::transform(sphere.center.begin(), sphere.center.end(), sphere.center.begin(), - DivideLength); + DivideLengthScale); } // Materials: conductivity and London penetration depth. for (auto &data : domains.materials) { data.sigma /= 1.0 / (electromagnetics::Z0_ * Lc); - data.lambda_L /= Lc / model.L0; + data.lambda_L /= GetMeshLengthScale(); } // Probe location coordinates. for (auto &[idx, data] : domains.postpro.probe) { std::transform(data.center.begin(), data.center.end(), data.center.begin(), - DivideLength); + DivideLengthScale); } // Finite conductivity boundaries. for (auto &data : boundaries.conductivity) { data.sigma /= 1.0 / (electromagnetics::Z0_ * Lc); - data.h /= Lc / model.L0; + data.h /= GetMeshLengthScale(); } // Impedance boundaries and lumped ports. @@ -531,20 +529,20 @@ void IoData::NondimensionalizeInputs(mfem::ParMesh &mesh) // Wave port offset distance. for (auto &[idx, data] : boundaries.waveport) { - data.d_offset /= Lc / model.L0; + data.d_offset /= GetMeshLengthScale(); } // Center coordinates for surface flux. for (auto &[idx, data] : boundaries.postpro.flux) { std::transform(data.center.begin(), data.center.end(), data.center.begin(), - DivideLength); + DivideLengthScale); } // Dielectric interface thickness. for (auto &[idx, data] : boundaries.postpro.dielectric) { - data.t /= Lc / model.L0; + data.t /= GetMeshLengthScale(); } // For eigenmode simulations: @@ -562,7 +560,7 @@ void IoData::NondimensionalizeInputs(mfem::ParMesh &mesh) solver.transient.delta_t /= tc; // Scale mesh vertices for correct nondimensionalization. - mesh::NondimensionalizeMesh(mesh, GetLengthScale()); + mesh::NondimensionalizeMesh(mesh, GetMeshLengthScale()); // Print some information. Mpi::Print(mesh.GetComm(), diff --git a/palace/utils/iodata.hpp b/palace/utils/iodata.hpp index b57d3082f..39e8a5fa6 100644 --- a/palace/utils/iodata.hpp +++ b/palace/utils/iodata.hpp @@ -46,7 +46,7 @@ class IoData void NondimensionalizeInputs(mfem::ParMesh &mesh); // Return the mesh scaling factor in units model.L0 x [m] for mesh IO. - double GetLengthScale() const { return Lc / model.L0; } + double GetMeshLengthScale() const { return Lc / model.L0; } // Redimensionalize values for output. Outputs which depend on the fields assume a // characteristic reference magnetic field strength Hc such that Pc = 1 W, where Pc is the diff --git a/scripts/schema/config/model.json b/scripts/schema/config/model.json index d6e263f6e..148ab890f 100644 --- a/scripts/schema/config/model.json +++ b/scripts/schema/config/model.json @@ -9,9 +9,12 @@ "Mesh": { "type": "string" }, "L0": { "type": "number", "exclusiveMinimum": 0.0 }, "Lc": { "type": "number", "exclusiveMinimum": 0.0 }, - "Partition": { "type": "string" }, - "ReorientTetMesh": { "type": "boolean" }, "RemoveCurvature": { "type": "boolean" }, + "ReorderElements": { "type": "boolean" }, + "CleanUnusedElements": { "type": "boolean" }, + "AddInterfaceBoundaryElements": { "type": "boolean" }, + "ReorientTetMesh": { "type": "boolean" }, + "Partitioning": { "type": "string" }, "Refinement": { "type": "object", From f3e307f3ce36d8b0edcf83909f018582c18bd34b Mon Sep 17 00:00:00 2001 From: Sebastian Grimberg Date: Fri, 24 May 2024 13:51:55 -0700 Subject: [PATCH 03/14] Minor docs cleanup --- examples/rings/mesh/mesh.jl | 2 +- examples/rings/rings.json | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/examples/rings/mesh/mesh.jl b/examples/rings/mesh/mesh.jl index b7aca9ec4..71c16337f 100644 --- a/examples/rings/mesh/mesh.jl +++ b/examples/rings/mesh/mesh.jl @@ -54,7 +54,7 @@ function generate_ring_mesh(; end gmsh.model.add("rings") - # Geometry parameters (in um) + # Geometry parameters (in μm) farfield_radius = 10.0 * outer_radius # Mesh parameters diff --git a/examples/rings/rings.json b/examples/rings/rings.json index 0c6134cc2..c0d96c232 100644 --- a/examples/rings/rings.json +++ b/examples/rings/rings.json @@ -8,7 +8,7 @@ "Model": { "Mesh": "mesh/rings.msh", - "L0": 1.0e-6 // um + "L0": 1.0e-6 // μm }, "Domains": { From cab5f61fce1134350071a443a21e19674703b69f Mon Sep 17 00:00:00 2001 From: Sebastian Grimberg Date: Wed, 29 May 2024 20:46:02 -0700 Subject: [PATCH 04/14] Add function interpolation between meshes, which enables preservation of curvature when splitting mesh into simplices --- docs/src/config/model.md | 1 + palace/fem/interpolator.cpp | 182 +++++++++++++++++++++++++++++-- palace/fem/interpolator.hpp | 9 ++ palace/utils/configfile.cpp | 3 + palace/utils/configfile.hpp | 3 + palace/utils/geodata.cpp | 165 ++++++++++++++++++---------- scripts/schema/config/model.json | 1 + 7 files changed, 298 insertions(+), 66 deletions(-) diff --git a/docs/src/config/model.md b/docs/src/config/model.md index f31daddd3..b727c9250 100644 --- a/docs/src/config/model.md +++ b/docs/src/config/model.md @@ -114,6 +114,7 @@ mesh length units. ### Advanced model options - `"RemoveCurvature" [false]` + - `"MakeSimplex" [false]` - `"ReorderElements" [false]` - `"CleanUnusedElements" [true]` - `"AddInterfaceBoundaryElements" [true]` diff --git a/palace/fem/interpolator.cpp b/palace/fem/interpolator.cpp index a24abc25f..29664f0e0 100644 --- a/palace/fem/interpolator.cpp +++ b/palace/fem/interpolator.cpp @@ -24,13 +24,14 @@ InterpolationOperator::InterpolationOperator(const IoData &iodata, mfem::ParMesh { return; } - const double bb_t = 0.1; // MFEM defaults - const double newton_tol = 1.0e-12; - const int npts = static_cast(iodata.domains.postpro.probe.size()); + constexpr double bb_t = 0.1; // MFEM defaults + constexpr double newton_tol = 1.0e-12; + const int dim = mesh.SpaceDimension(); MFEM_VERIFY( - mesh.Dimension() == mesh.SpaceDimension(), + mesh.Dimension() == dim, "Probe postprocessing functionality requires mesh dimension == space dimension!"); - mfem::Vector xyz(npts * mesh.SpaceDimension()); + const int npts = static_cast(iodata.domains.postpro.probe.size()); + mfem::Vector xyz(npts * dim); op_idx.resize(npts); int i = 0; for (const auto &[idx, data] : iodata.domains.postpro.probe) @@ -38,7 +39,7 @@ InterpolationOperator::InterpolationOperator(const IoData &iodata, mfem::ParMesh // Use default ordering byNODES. xyz(i) = data.center[0]; xyz(npts + i) = data.center[1]; - if (mesh.SpaceDimension() == 3) + if (dim == 3) { xyz(2 * npts + i) = data.center[2]; } @@ -72,15 +73,15 @@ std::vector InterpolationOperator::ProbeField(const mfem::ParGridFunctio // Interpolated vector values are returned from GSLIB interpolator byNODES, which we // transform to byVDIM for output. const int npts = op.GetCode().Size(); - const int dim = U.VectorDim(); - std::vector vals(npts * dim); - mfem::Vector v(npts * dim); + const int vdim = U.VectorDim(); + std::vector vals(npts * vdim); + mfem::Vector v(npts * vdim); op.Interpolate(U, v); - for (int d = 0; d < dim; d++) + for (int d = 0; d < vdim; d++) { for (int i = 0; i < npts; i++) { - vals[i * dim + d] = v(d * npts + i); + vals[i * vdim + d] = v(d * npts + i); } } return vals; @@ -107,4 +108,163 @@ std::vector> InterpolationOperator::ProbeField(const GridFu } } +namespace fem +{ + +void InterpolateFunction(const mfem::GridFunction &U, mfem::GridFunction &V) +{ +#if defined(MFEM_USE_GSLIB) + // Generate list of points where the grid function will be evaluated. + auto &dest_mesh = *V.FESpace()->GetMesh(); + MFEM_VERIFY(dest_mesh.GetNodes(), "Destination mesh has no nodal FE space!"); + const int dim = dest_mesh.SpaceDimension(); + mfem::Vector xyz; + mfem::Ordering::Type ordering; + const auto *dest_fec_h1 = + dynamic_cast(V.FESpace()->FEColl()); + const auto *dest_fec_l2 = + dynamic_cast(V.FESpace()->FEColl()); + const auto *dest_nodes_h1 = dynamic_cast( + dest_mesh.GetNodes()->FESpace()->FEColl()); + const auto *dest_nodes_l2 = dynamic_cast( + dest_mesh.GetNodes()->FESpace()->FEColl()); + int dest_fespace_order = V.FESpace()->GetMaxElementOrder(); + int dest_nodes_order = dest_mesh.GetNodes()->FESpace()->GetMaxElementOrder(); + dest_mesh.GetNodes()->HostRead(); + if (((dest_fec_h1 && dest_nodes_h1) || (dest_fec_l2 && dest_nodes_l2)) && + dest_fespace_order == dest_nodes_order) + { + xyz.MakeRef(*dest_mesh.GetNodes(), 0, dest_mesh.GetNodes()->Size()); + ordering = dest_mesh.GetNodes()->FESpace()->GetOrdering(); + } + else + { + int npts = 0, offset = 0; + for (int i = 0; i < dest_mesh.GetNE(); i++) + { + npts += V.FESpace()->GetFE(i)->GetNodes().GetNPoints(); + } + xyz.SetSize(npts * dim); + mfem::DenseMatrix pointmat; + for (int i = 0; i < dest_mesh.GetNE(); i++) + { + const mfem::FiniteElement &fe = *V.FESpace()->GetFE(i); + mfem::ElementTransformation &T = *V.FESpace()->GetElementTransformation(i); + T.Transform(fe.GetNodes(), pointmat); + for (int j = 0; j < pointmat.Width(); j++) + { + // Fill with ordering byNODES. + xyz(offset + j) = pointmat(0, j); + xyz(npts + offset + j) = pointmat(1, j); + if (dim == 3) + { + xyz(2 * npts + offset + j) = pointmat(2, j); + } + } + offset += pointmat.Width(); + } + ordering = mfem::Ordering::byNODES; + } + const int npts = xyz.Size() / dim; + + // Set up the interpolator. + auto &src_mesh = *U.FESpace()->GetMesh(); + MFEM_VERIFY(src_mesh.GetNodes(), "Source mesh has no nodal FE space!"); + auto *src_pmesh = dynamic_cast(&src_mesh); + MPI_Comm comm = (src_pmesh) ? src_pmesh->GetComm() : MPI_COMM_SELF; + mfem::FindPointsGSLIB op(comm); + constexpr double bb_t = 0.1; // MFEM defaults + constexpr double newton_tol = 1.0e-12; + op.Setup(src_mesh, bb_t, newton_tol, npts); + + // Perform the interpolation and fill the target GridFunction (see MFEM's field-interp + // miniapp). + const int vdim = V.VectorDim(); + mfem::Vector vals(npts * vdim); + op.SetDefaultInterpolationValue(0.0); + op.SetL2AvgType(mfem::FindPointsGSLIB::NONE); + op.Interpolate(xyz, U, vals, ordering); + if (dest_fec_h1 || dest_fec_l2) + { + if ((dest_fec_h1 && dest_fespace_order == dest_nodes_order) || dest_fec_l2) + { + // H1 and L2 (copy interpolated values to vdofs). + MFEM_ASSERT(V.Size() == vals.Size(), + "Unexpected size mismatch for interpolated values and grid function!"); + V = vals; + } + else + { + // H1 with order != mesh order needs to handle duplicated points. + mfem::Vector elem_vals; + mfem::Array vdofs; + int offset = 0; + for (int i = 0; i < dest_mesh.GetNE(); i++) + { + const mfem::FiniteElement &fe = *V.FESpace()->GetFE(i); + const int elem_npts = fe.GetNodes().GetNPoints(); + elem_vals.SetSize(elem_npts * vdim); + + for (int j = 0; j < elem_npts; j++) + { + for (int d = 0; d < vdim; d++) + { + // Arrange values byNODES to align with GetElementVDofs. + int idx = (U.FESpace()->GetOrdering() == mfem::Ordering::byNODES) + ? d * npts + offset + j + : (offset + j) * dim + d; + elem_vals(d * elem_npts + j) = vals(idx); + } + } + offset += elem_npts; + const auto *dof_trans = V.FESpace()->GetElementVDofs(i, vdofs); + if (dof_trans) + { + dof_trans->TransformPrimal(elem_vals); + } + V.SetSubVector(vdofs, elem_vals); + } + } + } + else + { + // H(div) or H(curl) use ProjectFromNodes. + mfem::Vector elem_vals, v; + mfem::Array vdofs; + int offset = 0; + for (int i = 0; i < dest_mesh.GetNE(); i++) + { + const mfem::FiniteElement &fe = *V.FESpace()->GetFE(i); + mfem::ElementTransformation &T = *V.FESpace()->GetElementTransformation(i); + const int elem_npts = fe.GetNodes().GetNPoints(); + elem_vals.SetSize(elem_npts * vdim); + for (int j = 0; j < elem_npts; j++) + { + for (int d = 0; d < vdim; d++) + { + // Arrange values byVDIM. + int idx = (U.FESpace()->GetOrdering() == mfem::Ordering::byNODES) + ? d * npts + offset + j + : (offset + j) * dim + d; + elem_vals(j * vdim + d) = vals(idx); + } + } + offset += elem_npts; + const auto *dof_trans = V.FESpace()->GetElementVDofs(i, vdofs); + v.SetSize(vdofs.Size()); + fe.ProjectFromNodes(elem_vals, T, v); + if (dof_trans) + { + dof_trans->TransformPrimal(v); + } + V.SetSubVector(vdofs, v); + } + } +#else + MFEM_ABORT("InterpolateFunction requires MFEM_USE_GSLIB!"); +#endif +} + +} // namespace fem + } // namespace palace diff --git a/palace/fem/interpolator.hpp b/palace/fem/interpolator.hpp index 460caa2e9..757eb4108 100644 --- a/palace/fem/interpolator.hpp +++ b/palace/fem/interpolator.hpp @@ -35,6 +35,15 @@ class InterpolationOperator std::vector> ProbeField(const GridFunction &U); }; +namespace fem +{ + +// Interpolate a function on a serial or parallel mesh to a different mesh, using GSLIB. +// Similar to MFEM's field-interp miniapp. +void InterpolateFunction(const mfem::GridFunction &U, mfem::GridFunction &V); + +} // namespace fem + } // namespace palace #endif // PALACE_FEM_INTERPOLATOR_HPP diff --git a/palace/utils/configfile.cpp b/palace/utils/configfile.cpp index 7db0cff3d..53c4104f7 100644 --- a/palace/utils/configfile.cpp +++ b/palace/utils/configfile.cpp @@ -463,6 +463,7 @@ void ModelData::SetUp(json &config) L0 = model->value("L0", L0); Lc = model->value("Lc", Lc); remove_curvature = model->value("RemoveCurvature", remove_curvature); + make_simplex = model->value("MakeSimplex", make_simplex); 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); @@ -475,6 +476,7 @@ void ModelData::SetUp(json &config) model->erase("L0"); model->erase("Lc"); model->erase("RemoveCurvature"); + model->erase("MakeSimplex"); model->erase("ReorderElements"); model->erase("CleanUnusedElements"); model->erase("AddInterfaceBoundaryElements"); @@ -492,6 +494,7 @@ void ModelData::SetUp(json &config) std::cout << "L0: " << L0 << '\n'; std::cout << "Lc: " << Lc << '\n'; std::cout << "RemoveCurvature: " << remove_curvature << '\n'; + std::cout << "MakeSimplex: " << make_simplex << '\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 19e529996..40c502eb6 100644 --- a/palace/utils/configfile.hpp +++ b/palace/utils/configfile.hpp @@ -202,6 +202,9 @@ struct ModelData // Remove high-order curvature information from the mesh. bool remove_curvature = false; + // Convert mesh to simplex elements. + bool make_simplex = 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 809a8feb8..ec936e356 100644 --- a/palace/utils/geodata.cpp +++ b/palace/utils/geodata.cpp @@ -14,6 +14,7 @@ #include #include #include +#include "fem/interpolator.hpp" #include "utils/communication.hpp" #include "utils/diagnostic.hpp" #include "utils/filesystem.hpp" @@ -43,15 +44,20 @@ std::unique_ptr LoadMesh(const std::string &, bool); // cache usage. void ReorderMesh(mfem::Mesh &, bool = true); -// Generate element-based mesh partitioning, using either a provided file or METIS. -std::unique_ptr GetMeshPartitioning(const mfem::Mesh &, int, - const std::string & = "", 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); // Clean the provided serial mesh by removing unnecessary domain and elements, and adding // boundary elements for material interfaces and exterior boundaries. std::map> CheckMesh(std::unique_ptr &, const IoData &, bool, bool); +// Generate element-based mesh partitioning, using either a provided file or METIS. +std::unique_ptr GetMeshPartitioning(const mfem::Mesh &, int, + const std::string & = "", bool = true); + // Given a serial mesh on the root processor and element partitioning, create a parallel // mesh over the given communicator. The serial mesh is destroyed when no longer needed. std::unique_ptr DistributeMesh(MPI_Comm, std::unique_ptr &, @@ -93,13 +99,23 @@ std::unique_ptr ReadMesh(MPI_Comm comm, const IoData &iodata) { // Check the the AMR specification and the mesh elements are compatible. const auto element_types = CheckElements(*smesh); - MFEM_VERIFY(!use_amr || !element_types.has_hexahedra || refinement.nonconformal, + MFEM_VERIFY(!use_amr || iodata.model.make_simplex || !element_types.has_hexahedra || + refinement.nonconformal, "If there are tensor elements, AMR must be nonconformal!"); - MFEM_VERIFY(!use_amr || !element_types.has_pyramids || refinement.nonconformal, + MFEM_VERIFY(!use_amr || iodata.model.make_simplex || !element_types.has_pyramids || + refinement.nonconformal, "If there are pyramid elements, AMR must be nonconformal!"); - MFEM_VERIFY(!use_amr || !element_types.has_prisms || refinement.nonconformal, + MFEM_VERIFY(!use_amr || iodata.model.make_simplex || !element_types.has_prisms || + refinement.nonconformal, "If there are wedge 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) + { + MakeSimplicial(*smesh); + } + // Optionally reorder elements (and vertices) based on spatial location after loading // the serial mesh. if (iodata.model.reorder_elements) @@ -1563,7 +1579,6 @@ double RebalanceMesh(std::unique_ptr &mesh, const IoData &iodata) } if (ratio > tol) { - mesh->ExchangeFaceNbrData(); if (mesh->Nonconforming()) { mesh->Rebalance(); @@ -1574,7 +1589,6 @@ double RebalanceMesh(std::unique_ptr &mesh, const IoData &iodata) // processor and then redistributed. RebalanceConformalMesh(mesh); } - mesh->ExchangeFaceNbrData(); } return ratio; } @@ -1706,58 +1720,46 @@ void ReorderMesh(mfem::Mesh &mesh, bool print) mesh.ReorderElements(ordering); } -std::unique_ptr GetMeshPartitioning(const mfem::Mesh &mesh, int size, - const std::string &part_file, bool print) +void MakeSimplicial(mfem::Mesh &orig_mesh, bool preserve_curvature) { - MFEM_VERIFY(size <= mesh.GetNE(), "Mesh partitioning must have parts <= mesh elements (" - << size << " vs. " << mesh.GetNE() << ")!"); - if (part_file.length() == 0) + // Back up the original high-order nodes. + mfem::GridFunction *nodes = nullptr; + int own_nodes = 1; + if (preserve_curvature && orig_mesh.GetNodes()) { - const int part_method = 1; - std::unique_ptr partitioning( - const_cast(mesh).GeneratePartitioning(size, part_method)); - if (print) - { - Mpi::Print("Finished partitioning mesh into {:d} subdomain{}\n", size, - (size > 1) ? "s" : ""); - } - return partitioning; + orig_mesh.SwapNodes(nodes, own_nodes); } - // User can optionally specify a mesh partitioning file as generated from the MFEM - // mesh-explorer miniapp, for example. It has the format: - // - // number_of_elements - // number_of_processors - // - // ... - // - // - int ne, np; - std::ifstream part_ifs(part_file); - part_ifs.ignore(std::numeric_limits::max(), ' '); - part_ifs >> ne; - if (ne != mesh.GetNE()) - { - MFEM_ABORT("Invalid partitioning file (number of elements)!"); - } - part_ifs.ignore(std::numeric_limits::max(), ' '); - part_ifs >> np; - if (np != size) + + // MFEM's function removes curvature information from the new mesh. So, if needed, we + // interpolate it onto the new mesh with GSLIB. + mfem::Mesh new_mesh = mfem::Mesh::MakeSimplicial(orig_mesh); + if (preserve_curvature && nodes) { - MFEM_ABORT("Invalid partitioning file (number of processors)!"); + // Prepare to interpolate the grid function for high-order nodes from the old mesh + // on the new one. + orig_mesh.EnsureNodes(); + new_mesh.EnsureNodes(); + const mfem::FiniteElementSpace *fespace = nodes->FESpace(); + mfem::Ordering::Type ordering = fespace->GetOrdering(); + int order = fespace->GetMaxElementOrder(); + int sdim = orig_mesh.SpaceDimension(); + bool discont = + (dynamic_cast(fespace->FEColl()) != nullptr); + mfem::FiniteElementSpace new_fespace(&new_mesh, fespace->FEColl(), sdim, ordering); + mfem::GridFunction new_nodes(&new_fespace); + fem::InterpolateFunction(*nodes, new_nodes); + + // Finally, copy the nodal grid function to the new mesh. + new_mesh.SetCurvature(order, discont, sdim, ordering); + MFEM_VERIFY(new_mesh.GetNodes()->Size() == new_nodes.Size(), + "Unexpected size mismatch for nodes!"); + new_mesh.SetNodes(new_nodes); } - auto partitioning = std::make_unique(mesh.GetNE()); - int i = 0; - while (i < mesh.GetNE()) + if (own_nodes) { - part_ifs >> partitioning[i++]; + delete nodes; } - if (print) - { - Mpi::Print("Read mesh partitioning into {:d} subdomain{} from disk\n", size, - (size > 1) ? "s" : ""); - } - return partitioning; + orig_mesh = std::move(new_mesh); } std::map> CheckMesh(std::unique_ptr &orig_mesh, @@ -2043,12 +2045,11 @@ std::map> CheckMesh(std::unique_ptr &orig_me { const mfem::GridFunction *nodes = orig_mesh->GetNodes(); const mfem::FiniteElementSpace *fespace = nodes->FESpace(); - mfem::Ordering::Type ordering = fespace->GetOrdering(); int order = fespace->GetMaxElementOrder(); int sdim = orig_mesh->SpaceDimension(); bool discont = - dynamic_cast(fespace->FEColl()) != nullptr; + (dynamic_cast(fespace->FEColl()) != nullptr); new_mesh->SetCurvature(order, discont, sdim, ordering); mfem::GridFunction *new_nodes = new_mesh->GetNodes(); const mfem::FiniteElementSpace *new_fespace = new_nodes->FESpace(); @@ -2076,6 +2077,60 @@ std::map> CheckMesh(std::unique_ptr &orig_me return new_attr_map; } +std::unique_ptr GetMeshPartitioning(const mfem::Mesh &mesh, int size, + const std::string &part_file, bool print) +{ + MFEM_VERIFY(size <= mesh.GetNE(), "Mesh partitioning must have parts <= mesh elements (" + << size << " vs. " << mesh.GetNE() << ")!"); + if (part_file.length() == 0) + { + const int part_method = 1; + std::unique_ptr partitioning( + const_cast(mesh).GeneratePartitioning(size, part_method)); + if (print) + { + Mpi::Print("Finished partitioning mesh into {:d} subdomain{}\n", size, + (size > 1) ? "s" : ""); + } + return partitioning; + } + // User can optionally specify a mesh partitioning file as generated from the MFEM + // mesh-explorer miniapp, for example. It has the format: + // + // number_of_elements + // number_of_processors + // + // ... + // + // + int ne, np; + std::ifstream part_ifs(part_file); + part_ifs.ignore(std::numeric_limits::max(), ' '); + part_ifs >> ne; + if (ne != mesh.GetNE()) + { + MFEM_ABORT("Invalid partitioning file (number of elements)!"); + } + part_ifs.ignore(std::numeric_limits::max(), ' '); + part_ifs >> np; + if (np != size) + { + MFEM_ABORT("Invalid partitioning file (number of processors)!"); + } + auto partitioning = std::make_unique(mesh.GetNE()); + int i = 0; + while (i < mesh.GetNE()) + { + part_ifs >> partitioning[i++]; + } + if (print) + { + Mpi::Print("Read mesh partitioning into {:d} subdomain{} from disk\n", size, + (size > 1) ? "s" : ""); + } + return partitioning; +} + std::unique_ptr DistributeMesh(MPI_Comm comm, std::unique_ptr &smesh, const int *partitioning, diff --git a/scripts/schema/config/model.json b/scripts/schema/config/model.json index 148ab890f..4de2435a3 100644 --- a/scripts/schema/config/model.json +++ b/scripts/schema/config/model.json @@ -10,6 +10,7 @@ "L0": { "type": "number", "exclusiveMinimum": 0.0 }, "Lc": { "type": "number", "exclusiveMinimum": 0.0 }, "RemoveCurvature": { "type": "boolean" }, + "MakeSimplex": { "type": "boolean" }, "ReorderElements": { "type": "boolean" }, "CleanUnusedElements": { "type": "boolean" }, "AddInterfaceBoundaryElements": { "type": "boolean" }, From f4bdc6137412a016772a1f1b6925f9c9035352d5 Mon Sep 17 00:00:00 2001 From: Sebastian Grimberg Date: Thu, 30 May 2024 09:17:11 -0700 Subject: [PATCH 05/14] Fix interpolation for non-H1 spaces --- palace/fem/interpolator.cpp | 37 ++++++++++++++++++------------------- 1 file changed, 18 insertions(+), 19 deletions(-) diff --git a/palace/fem/interpolator.cpp b/palace/fem/interpolator.cpp index 29664f0e0..712f8c501 100644 --- a/palace/fem/interpolator.cpp +++ b/palace/fem/interpolator.cpp @@ -114,7 +114,10 @@ namespace fem void InterpolateFunction(const mfem::GridFunction &U, mfem::GridFunction &V) { #if defined(MFEM_USE_GSLIB) - // Generate list of points where the grid function will be evaluated. + // Generate list of points where the grid function will be evaluated. If the grid function + // to interpolate is an H1 space of the same order as the mesh nodes, we can use the + // mesh node points directly. Otherwise, for a different basis order or type, we generate + // the interpolation points in the physical space manually. auto &dest_mesh = *V.FESpace()->GetMesh(); MFEM_VERIFY(dest_mesh.GetNodes(), "Destination mesh has no nodal FE space!"); const int dim = dest_mesh.SpaceDimension(); @@ -122,17 +125,12 @@ void InterpolateFunction(const mfem::GridFunction &U, mfem::GridFunction &V) mfem::Ordering::Type ordering; const auto *dest_fec_h1 = dynamic_cast(V.FESpace()->FEColl()); - const auto *dest_fec_l2 = - dynamic_cast(V.FESpace()->FEColl()); const auto *dest_nodes_h1 = dynamic_cast( dest_mesh.GetNodes()->FESpace()->FEColl()); - const auto *dest_nodes_l2 = dynamic_cast( - dest_mesh.GetNodes()->FESpace()->FEColl()); int dest_fespace_order = V.FESpace()->GetMaxElementOrder(); int dest_nodes_order = dest_mesh.GetNodes()->FESpace()->GetMaxElementOrder(); dest_mesh.GetNodes()->HostRead(); - if (((dest_fec_h1 && dest_nodes_h1) || (dest_fec_l2 && dest_nodes_l2)) && - dest_fespace_order == dest_nodes_order) + if (dest_fec_h1 && dest_nodes_h1 && dest_fespace_order == dest_nodes_order) { xyz.MakeRef(*dest_mesh.GetNodes(), 0, dest_mesh.GetNodes()->Size()); ordering = dest_mesh.GetNodes()->FESpace()->GetOrdering(); @@ -184,18 +182,13 @@ void InterpolateFunction(const mfem::GridFunction &U, mfem::GridFunction &V) op.SetDefaultInterpolationValue(0.0); op.SetL2AvgType(mfem::FindPointsGSLIB::NONE); op.Interpolate(xyz, U, vals, ordering); + const auto *dest_fec_l2 = + dynamic_cast(V.FESpace()->FEColl()); if (dest_fec_h1 || dest_fec_l2) { - if ((dest_fec_h1 && dest_fespace_order == dest_nodes_order) || dest_fec_l2) - { - // H1 and L2 (copy interpolated values to vdofs). - MFEM_ASSERT(V.Size() == vals.Size(), - "Unexpected size mismatch for interpolated values and grid function!"); - V = vals; - } - else + if (dest_fec_h1 && dest_fespace_order != dest_nodes_order) { - // H1 with order != mesh order needs to handle duplicated points. + // H1 with order != mesh order needs to handle duplicated interpolation points. mfem::Vector elem_vals; mfem::Array vdofs; int offset = 0; @@ -204,12 +197,11 @@ void InterpolateFunction(const mfem::GridFunction &U, mfem::GridFunction &V) const mfem::FiniteElement &fe = *V.FESpace()->GetFE(i); const int elem_npts = fe.GetNodes().GetNPoints(); elem_vals.SetSize(elem_npts * vdim); - for (int j = 0; j < elem_npts; j++) { for (int d = 0; d < vdim; d++) { - // Arrange values byNODES to align with GetElementVDofs. + // Arrange element values byNODES to align with GetElementVDofs. int idx = (U.FESpace()->GetOrdering() == mfem::Ordering::byNODES) ? d * npts + offset + j : (offset + j) * dim + d; @@ -225,6 +217,13 @@ void InterpolateFunction(const mfem::GridFunction &U, mfem::GridFunction &V) V.SetSubVector(vdofs, elem_vals); } } + else + { + // Otherwise, H1 and L2 copy interpolated values to vdofs. + MFEM_ASSERT(V.Size() == vals.Size(), + "Unexpected size mismatch for interpolated values and grid function!"); + V = vals; + } } else { @@ -242,7 +241,7 @@ void InterpolateFunction(const mfem::GridFunction &U, mfem::GridFunction &V) { for (int d = 0; d < vdim; d++) { - // Arrange values byVDIM. + // Arrange element values byVDIM for ProjectFromNodes. int idx = (U.FESpace()->GetOrdering() == mfem::Ordering::byNODES) ? d * npts + offset + j : (offset + j) * dim + d; From 22af1ff56fd698d6fcadc6cf0f1ec971b34062fe Mon Sep 17 00:00:00 2001 From: Sebastian Grimberg Date: Thu, 30 May 2024 09:17:27 -0700 Subject: [PATCH 06/14] Add checks for pyramid elements when converting to simplices --- palace/utils/geodata.cpp | 35 +++++++++++++++++++++-------------- 1 file changed, 21 insertions(+), 14 deletions(-) diff --git a/palace/utils/geodata.cpp b/palace/utils/geodata.cpp index ec936e356..1178c30af 100644 --- a/palace/utils/geodata.cpp +++ b/palace/utils/geodata.cpp @@ -102,12 +102,12 @@ std::unique_ptr ReadMesh(MPI_Comm comm, const IoData &iodata) MFEM_VERIFY(!use_amr || iodata.model.make_simplex || !element_types.has_hexahedra || refinement.nonconformal, "If there are tensor elements, AMR must be nonconformal!"); - MFEM_VERIFY(!use_amr || iodata.model.make_simplex || !element_types.has_pyramids || - refinement.nonconformal, - "If there are pyramid elements, AMR must be nonconformal!"); MFEM_VERIFY(!use_amr || iodata.model.make_simplex || !element_types.has_prisms || refinement.nonconformal, "If there are wedge elements, AMR must be nonconformal!"); + MFEM_VERIFY(!use_amr || iodata.model.make_simplex || !element_types.has_pyramids || + refinement.nonconformal, + "If there are pyramid elements, AMR must be nonconformal!"); // Optionally convert mesh elements to simplices, for example in order to enable // conformal mesh refinement. @@ -1722,19 +1722,26 @@ void ReorderMesh(mfem::Mesh &mesh, bool print) void MakeSimplicial(mfem::Mesh &orig_mesh, bool preserve_curvature) { - // Back up the original high-order nodes. - mfem::GridFunction *nodes = nullptr; - int own_nodes = 1; - if (preserve_curvature && orig_mesh.GetNodes()) + // 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) { - orig_mesh.SwapNodes(nodes, own_nodes); + 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's function removes curvature information from the new mesh. So, if needed, we // interpolate it onto the new mesh with GSLIB. - mfem::Mesh new_mesh = mfem::Mesh::MakeSimplicial(orig_mesh); - if (preserve_curvature && nodes) + if (preserve_curvature && orig_mesh.GetNodes()) { + // Back up the original high-order nodes. + mfem::GridFunction *nodes = nullptr; + int own_nodes = 1; + orig_mesh.SwapNodes(nodes, own_nodes); + // Prepare to interpolate the grid function for high-order nodes from the old mesh // on the new one. orig_mesh.EnsureNodes(); @@ -1748,6 +1755,10 @@ void MakeSimplicial(mfem::Mesh &orig_mesh, bool preserve_curvature) mfem::FiniteElementSpace new_fespace(&new_mesh, fespace->FEColl(), sdim, ordering); mfem::GridFunction new_nodes(&new_fespace); fem::InterpolateFunction(*nodes, new_nodes); + if (own_nodes) + { + delete nodes; + } // Finally, copy the nodal grid function to the new mesh. new_mesh.SetCurvature(order, discont, sdim, ordering); @@ -1755,10 +1766,6 @@ void MakeSimplicial(mfem::Mesh &orig_mesh, bool preserve_curvature) "Unexpected size mismatch for nodes!"); new_mesh.SetNodes(new_nodes); } - if (own_nodes) - { - delete nodes; - } orig_mesh = std::move(new_mesh); } From 0bbd579182f752c969912de7c4df656e77c46288 Mon Sep 17 00:00:00 2001 From: Sebastian Grimberg Date: Thu, 30 May 2024 11:40:49 -0700 Subject: [PATCH 07/14] Add tet-to-hex mesh conversion option, also supporting high-order meshes --- docs/src/config/model.md | 1 + palace/utils/configfile.cpp | 3 + palace/utils/configfile.hpp | 3 + palace/utils/geodata.cpp | 231 +++++++++++++++++++++++++------ scripts/schema/config/model.json | 1 + 5 files changed, 200 insertions(+), 39 deletions(-) diff --git a/docs/src/config/model.md b/docs/src/config/model.md index b727c9250..41f05de9f 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 53c4104f7..e6dd79ed9 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 40c502eb6..27947e8ac 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 1178c30af..b212fc192 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 @@ -1688,53 +1688,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. @@ -1769,6 +1890,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" }, From f0c9bcf1bbe48d079c71f94d6ea300af3290e28d Mon Sep 17 00:00:00 2001 From: Sebastian Grimberg Date: Tue, 4 Jun 2024 16:34:06 -0700 Subject: [PATCH 08/14] Add serial refinement capability for debugging --- docs/src/config/model.md | 1 + palace/utils/configfile.cpp | 5 ++++- palace/utils/configfile.hpp | 3 +++ palace/utils/geodata.cpp | 16 ++++++++++++++++ scripts/schema/config/model.json | 1 + 5 files changed, 25 insertions(+), 1 deletion(-) diff --git a/docs/src/config/model.md b/docs/src/config/model.md index 41f05de9f..3a7da2d88 100644 --- a/docs/src/config/model.md +++ b/docs/src/config/model.md @@ -125,3 +125,4 @@ mesh length units. - `"MaximumImbalance" [1.1]` - `"SaveAdaptIterations" [true]` - `"SaveAdaptMesh" [false]` + - `"SerialUniformLevels" [0]` diff --git a/palace/utils/configfile.cpp b/palace/utils/configfile.cpp index e6dd79ed9..daddd69d7 100644 --- a/palace/utils/configfile.cpp +++ b/palace/utils/configfile.cpp @@ -337,7 +337,8 @@ void RefinementData::SetUp(json &model) // Options for a priori refinement. uniform_ref_levels = refinement->value("UniformLevels", uniform_ref_levels); - MFEM_VERIFY(uniform_ref_levels >= 0, + ser_uniform_ref_levels = refinement->value("SerialUniformLevels", ser_uniform_ref_levels); + MFEM_VERIFY(uniform_ref_levels >= 0 && ser_uniform_ref_levels >= 0, "Number of uniform mesh refinement levels must be non-negative!"); auto boxes = refinement->find("Boxes"); if (boxes != refinement->end()) @@ -430,6 +431,7 @@ void RefinementData::SetUp(json &model) refinement->erase("SaveAdaptIterations"); refinement->erase("SaveAdaptMesh"); refinement->erase("UniformLevels"); + refinement->erase("SerialUniformLevels"); refinement->erase("Boxes"); refinement->erase("Spheres"); MFEM_VERIFY(refinement->empty(), @@ -449,6 +451,7 @@ void RefinementData::SetUp(json &model) std::cout << "SaveAdaptIterations: " << save_adapt_iterations << '\n'; std::cout << "SaveAdaptMesh: " << save_adapt_mesh << '\n'; std::cout << "UniformLevels: " << uniform_ref_levels << '\n'; + std::cout << "SerialUniformLevels: " << ser_uniform_ref_levels << '\n'; } } diff --git a/palace/utils/configfile.hpp b/palace/utils/configfile.hpp index 27947e8ac..765f4e4af 100644 --- a/palace/utils/configfile.hpp +++ b/palace/utils/configfile.hpp @@ -171,6 +171,9 @@ struct RefinementData // Parallel uniform mesh refinement levels. int uniform_ref_levels = 0; + // Serial uniform mesh refinement levels. + int ser_uniform_ref_levels = 0; + private: // Refinement data for mesh regions. std::vector box_list = {}; diff --git a/palace/utils/geodata.cpp b/palace/utils/geodata.cpp index b212fc192..631456117 100644 --- a/palace/utils/geodata.cpp +++ b/palace/utils/geodata.cpp @@ -123,6 +123,22 @@ std::unique_ptr ReadMesh(MPI_Comm comm, const IoData &iodata) ReorderMeshElements(*smesh); } + // Refine the serial mesh (not typically used, prefer parallel uniform refinement + // instead). + { + int ne = smesh->GetNE(); + for (int l = 0; l < iodata.model.refinement.ser_uniform_ref_levels; l++) + { + smesh->UniformRefinement(); + } + if (iodata.model.refinement.ser_uniform_ref_levels > 0) + { + Mpi::Print("Serial uniform mesh refinement levels added {:d} elements (initial = " + "{:d}, final = {:d})\n", + smesh->GetNE() - ne, ne, smesh->GetNE()); + } + } + // Clean up unused domain elements from the mesh, add new boundary elements for material // interfaces if not present. Can only clean up conforming meshes, assumes that any // nonconformal mesh was generated by adaptation and thus does not need checking. diff --git a/scripts/schema/config/model.json b/scripts/schema/config/model.json index 556464dcd..cab779794 100644 --- a/scripts/schema/config/model.json +++ b/scripts/schema/config/model.json @@ -34,6 +34,7 @@ "SaveAdaptIterations": {"type": "boolean"}, "SaveAdaptMesh": {"type": "boolean"}, "UniformLevels": { "type": "integer", "minimum": 0 }, + "SerialUniformLevels": { "type": "integer", "minimum": 0 }, "Boxes": { "type": "array", From 4c188840cbb9027d76cde55c54ac73fbf568146e Mon Sep 17 00:00:00 2001 From: Sebastian Grimberg Date: Tue, 11 Jun 2024 11:27:16 -0700 Subject: [PATCH 09/14] Improve mesh partitioning for nonconforming meshes/meshes not using mesh partitioner class --- palace/utils/geodata.cpp | 230 ++++++++++++++++----------------------- 1 file changed, 95 insertions(+), 135 deletions(-) diff --git a/palace/utils/geodata.cpp b/palace/utils/geodata.cpp index 631456117..e876e1f0a 100644 --- a/palace/utils/geodata.cpp +++ b/palace/utils/geodata.cpp @@ -83,9 +83,17 @@ std::unique_ptr ReadMesh(MPI_Comm comm, const IoData &iodata) const auto &refinement = iodata.model.refinement; const bool use_amr = refinement.max_it > 0; const bool use_mesh_partitioner = !use_amr || !refinement.nonconformal; + MPI_Comm node_comm = nullptr; + if (!use_mesh_partitioner) + { + MPI_Comm_split_type(comm, MPI_COMM_TYPE_SHARED, Mpi::Rank(comm), MPI_INFO_NULL, + &node_comm); + } + + // Only one process per node reads the serial mesh. { BlockTimer bt(Timer::IO); - if (Mpi::Root(comm) || !use_mesh_partitioner) + if ((use_mesh_partitioner && Mpi::Root(comm)) || Mpi::Root(node_comm)) { smesh = LoadMesh(iodata.model.mesh, iodata.model.remove_curvature); MFEM_VERIFY(!(smesh->Nonconforming() && use_mesh_partitioner), @@ -94,8 +102,9 @@ std::unique_ptr ReadMesh(MPI_Comm comm, const IoData &iodata) Mpi::Barrier(comm); } + // Do some mesh preprocessing, and generate the partitioning. std::unique_ptr partitioning; - if (Mpi::Root(comm) || !use_mesh_partitioner) + if (smesh) { // Check the the AMR specification and the mesh elements are compatible. const auto element_types = CheckElements(*smesh); @@ -152,6 +161,7 @@ std::unique_ptr ReadMesh(MPI_Comm comm, const IoData &iodata) partitioning = GetMeshPartitioning(*smesh, Mpi::Size(comm), iodata.model.partitioning); } + // Distribute the mesh. std::unique_ptr pmesh; if (use_mesh_partitioner) { @@ -159,10 +169,47 @@ std::unique_ptr ReadMesh(MPI_Comm comm, const IoData &iodata) } else { + // Send the preprocessed serial mesh and partitioning as a byte string. + constexpr bool generate_edges = false, refine = true, fix_orientation = false; + std::string so; + int slen = 0; + if (smesh) + { + std::ostringstream fo(std::stringstream::out); + // fo << std::fixed; + fo << std::scientific; + fo.precision(MSH_FLT_PRECISION); + smesh->Print(fo); + smesh.reset(); // Root process needs to rebuild the mesh to ensure consistency with + // the saved serial mesh (refinement marking, for example) + so = fo.str(); + // so = zlib::CompressString(fo.str()); + slen = static_cast(so.size()); + MFEM_VERIFY(so.size() == (std::size_t)slen, "Overflow in stringbuffer size!"); + } + Mpi::Broadcast(1, &slen, 0, node_comm); + if (so.empty()) + { + so.resize(slen); + } + Mpi::Broadcast(slen, so.data(), 0, node_comm); + if (!smesh) + { + std::istringstream fi(so); + // std::istringstream fi(zlib::DecompressString(so)); + smesh = std::make_unique(fi, generate_edges, refine, fix_orientation); + } + so.clear(); if (refinement.nonconformal && use_amr) { smesh->EnsureNCMesh(true); } + if (!partitioning) + { + partitioning = std::make_unique(smesh->GetNE()); + } + Mpi::Broadcast(smesh->GetNE(), partitioning.get(), 0, node_comm); + MPI_Comm_free(&node_comm); pmesh = std::make_unique(comm, *smesh, partitioning.get()); smesh.reset(); } @@ -1635,23 +1682,6 @@ std::unique_ptr LoadMesh(const std::string &mesh_file, bool remove_c // fi << std::fixed; fi << std::scientific; fi.precision(MSH_FLT_PRECISION); - -#if 0 - // Put translated mesh in temporary storage (directory is created and destroyed in - // calling function). - std::string tmp = iodata.problem.output; - if (tmp.back() != '/') - { - tmp += '/'; - } - tmp += "tmp/serial.msh"; - std::ofstream fo(tmp); - // mfem::ofgzstream fo(tmp, true); // Use zlib compression if available - // fo << std::fixed; - fo << std::scientific; - fo.precision(MSH_FLT_PRECISION); -#endif - if (mesh_path.extension() == ".mphtxt" || mesh_path.extension() == ".mphbin") { mesh::ConvertMeshComsol(mesh_file, fi, remove_curvature); @@ -1662,16 +1692,6 @@ std::unique_ptr LoadMesh(const std::string &mesh_file, bool remove_c mesh::ConvertMeshNastran(mesh_file, fi, remove_curvature); // mesh::ConvertMeshNastran(mesh_file, fo, remove_curvature); } - -#if 0 - std::ifstream fi(tmp); - // mfem::ifgzstream fi(tmp); - if (!fi.good()) - { - MFEM_ABORT("Unable to open translated mesh file \"" << tmp << "\"!"); - } -#endif - mesh = std::make_unique(fi, generate_edges, refine, fix_orientation); } else @@ -2313,120 +2333,60 @@ std::unique_ptr DistributeMesh(MPI_Comm comm, const std::string &output_dir) { // Take a serial mesh and partitioning on the root process and construct the global - // parallel mesh. For now, prefer the MPI-based version. When constructing the ParMesh, we - // pass arguments to ensure no topological changes (this isn't required since the serial - // mesh was marked for refinement). - constexpr bool generate_edges = true, refine = true, fix_orientation = true; - if constexpr (false) + // parallel mesh. For now, prefer the MPI-based version to the file IO one. When + // constructing the ParMesh, we mark for refinement since refinement flags are not copied + // from the serial mesh. Beware that mfem::ParMesh constructor argument order is not the + // same as mfem::Mesh! Each processor's component gets sent as a byte string. + constexpr bool generate_edges = false, refine = true, fix_orientation = false; + std::unique_ptr pmesh; + if (Mpi::Root(comm)) { - // Write each processor's component to file. - std::string tmp = output_dir; - if (tmp.back() != '/') - { - tmp += '/'; - } - tmp += "tmp/"; - int width = 1 + static_cast(std::log10(Mpi::Size(comm) - 1)); - if (Mpi::Root(comm)) - { - if (!std::filesystem::exists(tmp)) - { - std::filesystem::create_directories(tmp); - } - mfem::MeshPartitioner partitioner(*smesh, Mpi::Size(comm), - const_cast(partitioning)); - for (int i = 0; i < Mpi::Size(comm); i++) + mfem::MeshPartitioner partitioner(*smesh, Mpi::Size(comm), + const_cast(partitioning)); + std::vector send_requests(Mpi::Size(comm) - 1, MPI_REQUEST_NULL); + std::vector so; + so.reserve(Mpi::Size(comm)); + for (int i = 0; i < Mpi::Size(comm); i++) + { + mfem::MeshPart part; + partitioner.ExtractPart(i, part); + std::ostringstream fo(std::stringstream::out); + // fo << std::fixed; + fo << std::scientific; + fo.precision(MSH_FLT_PRECISION); + part.Print(fo); + so.push_back(fo.str()); + // so.push_back((i > 0) ? zlib::CompressString(fo.str()) : fo.str()); + if (i > 0) { - mfem::MeshPart part; - partitioner.ExtractPart(i, part); - std::string pfile = mfem::MakeParFilename(tmp + "part.", i, ".mesh", width); - std::ofstream fo(pfile); - // mfem::ofgzstream fo(pfile, true); // Use zlib compression if available - // fo << std::fixed; - fo << std::scientific; - fo.precision(MSH_FLT_PRECISION); - part.Print(fo); + int slen = static_cast(so[i].length()); + MFEM_VERIFY(so[i].length() == (std::size_t)slen, + "Overflow error distributing parallel mesh!"); + MPI_Isend(so[i].data(), slen, MPI_CHAR, i, i, comm, &send_requests[i - 1]); } - smesh.reset(); - } - - // Each process loads its own partitioned mesh file and constructs the parallel mesh. - std::string pfile = - mfem::MakeParFilename(tmp + "part.", Mpi::Rank(comm), ".mesh", width); - int exists = 0; - while (!exists) // Wait for root to finish writing all files - { - exists = std::filesystem::exists(pfile); - Mpi::GlobalMax(1, &exists, comm); - } - std::ifstream fi(pfile); - // mfem::ifgzstream fi(pfile); - if (!fi.good()) - { - MFEM_ABORT("Unable to open partitioned mesh file \"" << pfile << "\"!"); - } - auto pmesh = - std::make_unique(comm, fi, generate_edges, refine, fix_orientation); - Mpi::Barrier(comm); // Wait for all processes to read the mesh from file - if (Mpi::Root(comm)) - { - std::filesystem::remove_all(tmp); // Remove the temporary directory } - return pmesh; + smesh.reset(); + std::istringstream fi(so[0]); // This is never compressed + pmesh = + std::make_unique(comm, fi, refine, generate_edges, fix_orientation); + MPI_Waitall(static_cast(send_requests.size()), send_requests.data(), + MPI_STATUSES_IGNORE); } else { - // Send each processor's component as a byte string. - std::unique_ptr pmesh; - if (Mpi::Root(comm)) - { - mfem::MeshPartitioner partitioner(*smesh, Mpi::Size(comm), - const_cast(partitioning)); - std::vector send_requests(Mpi::Size(comm) - 1, MPI_REQUEST_NULL); - std::vector so; - so.reserve(Mpi::Size(comm)); - for (int i = 0; i < Mpi::Size(comm); i++) - { - mfem::MeshPart part; - partitioner.ExtractPart(i, part); - std::ostringstream fo(std::stringstream::out); - // fo << std::fixed; - fo << std::scientific; - fo.precision(MSH_FLT_PRECISION); - part.Print(fo); - so.push_back(fo.str()); - // so.push_back((i > 0) ? zlib::CompressString(fo.str()) : fo.str()); - if (i > 0) - { - int slen = static_cast(so[i].length()); - MFEM_VERIFY(so[i].length() == (std::size_t)slen, - "Overflow error distributing parallel mesh!"); - MPI_Isend(so[i].data(), slen, MPI_CHAR, i, i, comm, &send_requests[i - 1]); - } - } - smesh.reset(); - std::istringstream fi(so[0]); // This is never compressed - pmesh = std::make_unique(comm, fi, generate_edges, refine, - fix_orientation); - MPI_Waitall(static_cast(send_requests.size()), send_requests.data(), - MPI_STATUSES_IGNORE); - } - else - { - MPI_Status status; - int rlen; - std::string si; - MPI_Probe(0, Mpi::Rank(comm), comm, &status); - MPI_Get_count(&status, MPI_CHAR, &rlen); - si.resize(rlen); - MPI_Recv(si.data(), rlen, MPI_CHAR, 0, Mpi::Rank(comm), comm, MPI_STATUS_IGNORE); - std::istringstream fi(si); - // std::istringstream fi(zlib::DecompressString(si)); - pmesh = std::make_unique(comm, fi, generate_edges, refine, - fix_orientation); - } - return pmesh; + MPI_Status status; + int rlen; + std::string si; + MPI_Probe(0, Mpi::Rank(comm), comm, &status); + MPI_Get_count(&status, MPI_CHAR, &rlen); + si.resize(rlen); + MPI_Recv(si.data(), rlen, MPI_CHAR, 0, Mpi::Rank(comm), comm, MPI_STATUS_IGNORE); + std::istringstream fi(si); + // std::istringstream fi(zlib::DecompressString(si)); + pmesh = + std::make_unique(comm, fi, refine, generate_edges, fix_orientation); } + return pmesh; } void RebalanceConformalMesh(std::unique_ptr &pmesh) From 67fd635092ee64cdbd262066b83b66a01b873855 Mon Sep 17 00:00:00 2001 From: Sebastian Grimberg Date: Wed, 19 Jun 2024 10:19:02 -0700 Subject: [PATCH 10/14] Improve high-order node transfer for mesh splitting --- palace/fem/interpolator.cpp | 102 +++++++++++++------- palace/fem/interpolator.hpp | 6 ++ palace/utils/geodata.cpp | 185 +++++++++++++++++++++++++++++++----- 3 files changed, 234 insertions(+), 59 deletions(-) diff --git a/palace/fem/interpolator.cpp b/palace/fem/interpolator.cpp index 712f8c501..50b017882 100644 --- a/palace/fem/interpolator.cpp +++ b/palace/fem/interpolator.cpp @@ -11,6 +11,14 @@ namespace palace { +namespace +{ + +constexpr auto GSLIB_BB_TOL = 0.01; // MFEM defaults, slightly reduced bounding box +constexpr auto GSLIB_NEWTON_TOL = 1.0e-12; + +} // namespace + #if defined(MFEM_USE_GSLIB) InterpolationOperator::InterpolationOperator(const IoData &iodata, mfem::ParMesh &mesh) : op(mesh.GetComm()) @@ -24,8 +32,6 @@ InterpolationOperator::InterpolationOperator(const IoData &iodata, mfem::ParMesh { return; } - constexpr double bb_t = 0.1; // MFEM defaults - constexpr double newton_tol = 1.0e-12; const int dim = mesh.SpaceDimension(); MFEM_VERIFY( mesh.Dimension() == dim, @@ -36,16 +42,14 @@ InterpolationOperator::InterpolationOperator(const IoData &iodata, mfem::ParMesh int i = 0; for (const auto &[idx, data] : iodata.domains.postpro.probe) { - // Use default ordering byNODES. - xyz(i) = data.center[0]; - xyz(npts + i) = data.center[1]; - if (dim == 3) + for (int d = 0; d < dim; d++) { - xyz(2 * npts + i) = data.center[2]; + // Use default ordering byNODES. + xyz(d * npts + i) = data.center[d]; } op_idx[i++] = idx; } - op.Setup(mesh, bb_t, newton_tol, npts); + op.Setup(mesh, GSLIB_BB_TOL, GSLIB_NEWTON_TOL, npts); op.FindPoints(xyz, mfem::Ordering::byNODES); op.SetDefaultInterpolationValue(0.0); i = 0; @@ -70,18 +74,26 @@ InterpolationOperator::InterpolationOperator(const IoData &iodata, mfem::ParMesh std::vector InterpolationOperator::ProbeField(const mfem::ParGridFunction &U) { #if defined(MFEM_USE_GSLIB) - // Interpolated vector values are returned from GSLIB interpolator byNODES, which we - // transform to byVDIM for output. + // Interpolated vector values are returned from GSLIB interpolator with the same ordering + // as the source grid function, which we transform to byVDIM for output. const int npts = op.GetCode().Size(); const int vdim = U.VectorDim(); std::vector vals(npts * vdim); - mfem::Vector v(npts * vdim); - op.Interpolate(U, v); - for (int d = 0; d < vdim; d++) + if (U.FESpace()->GetOrdering() == mfem::Ordering::byVDIM) + { + mfem::Vector v(vals.data(), npts * vdim); + op.Interpolate(U, v); + } + else { - for (int i = 0; i < npts; i++) + mfem::Vector v(npts * vdim); + op.Interpolate(U, v); + for (int d = 0; d < vdim; d++) { - vals[i * vdim + d] = v(d * npts + i); + for (int i = 0; i < npts; i++) + { + vals[i * vdim + d] = v(d * npts + i); + } } } return vals; @@ -147,16 +159,14 @@ void InterpolateFunction(const mfem::GridFunction &U, mfem::GridFunction &V) for (int i = 0; i < dest_mesh.GetNE(); i++) { const mfem::FiniteElement &fe = *V.FESpace()->GetFE(i); - mfem::ElementTransformation &T = *V.FESpace()->GetElementTransformation(i); + mfem::ElementTransformation &T = *dest_mesh.GetElementTransformation(i); T.Transform(fe.GetNodes(), pointmat); for (int j = 0; j < pointmat.Width(); j++) { - // Fill with ordering byNODES. - xyz(offset + j) = pointmat(0, j); - xyz(npts + offset + j) = pointmat(1, j); - if (dim == 3) + for (int d = 0; d < dim; d++) { - xyz(2 * npts + offset + j) = pointmat(2, j); + // Use default ordering byNODES. + xyz(d * npts + offset + j) = pointmat(d, j); } } offset += pointmat.Width(); @@ -171,13 +181,11 @@ void InterpolateFunction(const mfem::GridFunction &U, mfem::GridFunction &V) auto *src_pmesh = dynamic_cast(&src_mesh); MPI_Comm comm = (src_pmesh) ? src_pmesh->GetComm() : MPI_COMM_SELF; mfem::FindPointsGSLIB op(comm); - constexpr double bb_t = 0.1; // MFEM defaults - constexpr double newton_tol = 1.0e-12; - op.Setup(src_mesh, bb_t, newton_tol, npts); + op.Setup(src_mesh, GSLIB_BB_TOL, GSLIB_NEWTON_TOL, npts); // Perform the interpolation and fill the target GridFunction (see MFEM's field-interp // miniapp). - const int vdim = V.VectorDim(); + const int vdim = U.VectorDim(); mfem::Vector vals(npts * vdim); op.SetDefaultInterpolationValue(0.0); op.SetL2AvgType(mfem::FindPointsGSLIB::NONE); @@ -197,24 +205,24 @@ void InterpolateFunction(const mfem::GridFunction &U, mfem::GridFunction &V) const mfem::FiniteElement &fe = *V.FESpace()->GetFE(i); const int elem_npts = fe.GetNodes().GetNPoints(); elem_vals.SetSize(elem_npts * vdim); - for (int j = 0; j < elem_npts; j++) + for (int d = 0; d < vdim; d++) { - for (int d = 0; d < vdim; d++) + for (int j = 0; j < elem_npts; j++) { // Arrange element values byNODES to align with GetElementVDofs. int idx = (U.FESpace()->GetOrdering() == mfem::Ordering::byNODES) ? d * npts + offset + j - : (offset + j) * dim + d; + : (offset + j) * vdim + d; elem_vals(d * elem_npts + j) = vals(idx); } } - offset += elem_npts; const auto *dof_trans = V.FESpace()->GetElementVDofs(i, vdofs); if (dof_trans) { dof_trans->TransformPrimal(elem_vals); } V.SetSubVector(vdofs, elem_vals); + offset += elem_npts; } } else @@ -234,21 +242,20 @@ void InterpolateFunction(const mfem::GridFunction &U, mfem::GridFunction &V) for (int i = 0; i < dest_mesh.GetNE(); i++) { const mfem::FiniteElement &fe = *V.FESpace()->GetFE(i); - mfem::ElementTransformation &T = *V.FESpace()->GetElementTransformation(i); + mfem::ElementTransformation &T = *dest_mesh.GetElementTransformation(i); const int elem_npts = fe.GetNodes().GetNPoints(); elem_vals.SetSize(elem_npts * vdim); - for (int j = 0; j < elem_npts; j++) + for (int d = 0; d < vdim; d++) { - for (int d = 0; d < vdim; d++) + for (int j = 0; j < elem_npts; j++) { // Arrange element values byVDIM for ProjectFromNodes. int idx = (U.FESpace()->GetOrdering() == mfem::Ordering::byNODES) ? d * npts + offset + j - : (offset + j) * dim + d; + : (offset + j) * vdim + d; elem_vals(j * vdim + d) = vals(idx); } } - offset += elem_npts; const auto *dof_trans = V.FESpace()->GetElementVDofs(i, vdofs); v.SetSize(vdofs.Size()); fe.ProjectFromNodes(elem_vals, T, v); @@ -257,6 +264,7 @@ void InterpolateFunction(const mfem::GridFunction &U, mfem::GridFunction &V) dof_trans->TransformPrimal(v); } V.SetSubVector(vdofs, v); + offset += elem_npts; } } #else @@ -264,6 +272,32 @@ void InterpolateFunction(const mfem::GridFunction &U, mfem::GridFunction &V) #endif } +void InterpolateFunction(const mfem::Vector &xyz, const mfem::GridFunction &U, + mfem::Vector &vals, mfem::Ordering::Type ordering) +{ +#if defined(MFEM_USE_GSLIB) + // Set up the interpolator. + auto &src_mesh = *U.FESpace()->GetMesh(); + MFEM_VERIFY(src_mesh.GetNodes(), "Source mesh has no nodal FE space!"); + const int dim = src_mesh.SpaceDimension(); + const int npts = xyz.Size() / dim; + auto *src_pmesh = dynamic_cast(&src_mesh); + MPI_Comm comm = (src_pmesh) ? src_pmesh->GetComm() : MPI_COMM_SELF; + mfem::FindPointsGSLIB op(comm); + op.Setup(src_mesh, GSLIB_BB_TOL, GSLIB_NEWTON_TOL, npts); + + // Perform the interpolation, with the ordering of the returned values matching the + // ordering of the source grid function. + const int vdim = U.VectorDim(); + MFEM_VERIFY(vals.Size() == npts * vdim, "Incorrect size for interpolated values vector!"); + op.SetDefaultInterpolationValue(0.0); + op.SetL2AvgType(mfem::FindPointsGSLIB::NONE); + op.Interpolate(xyz, U, vals, ordering); +#else + MFEM_ABORT("InterpolateFunction requires MFEM_USE_GSLIB!"); +#endif +} + } // namespace fem } // namespace palace diff --git a/palace/fem/interpolator.hpp b/palace/fem/interpolator.hpp index 757eb4108..468512f20 100644 --- a/palace/fem/interpolator.hpp +++ b/palace/fem/interpolator.hpp @@ -42,6 +42,12 @@ namespace fem // Similar to MFEM's field-interp miniapp. void InterpolateFunction(const mfem::GridFunction &U, mfem::GridFunction &V); +// Interpolate a function at a specific list of points, specified using the provided +// ordering. The output vector values are always arranged byVDIM. +void InterpolateFunction(const mfem::Vector &xyz, const mfem::GridFunction &U, + mfem::Vector &V, + mfem::Ordering::Type ordering = mfem::Ordering::byNODES); + } // namespace fem } // namespace palace diff --git a/palace/utils/geodata.cpp b/palace/utils/geodata.cpp index e876e1f0a..729228ea6 100644 --- a/palace/utils/geodata.cpp +++ b/palace/utils/geodata.cpp @@ -43,7 +43,7 @@ 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); +void SplitMeshElements(std::unique_ptr &, bool, bool, bool = true); // Optionally reorder mesh elements based on MFEM's internal reordeing tools for improved // cache usage. @@ -122,7 +122,7 @@ std::unique_ptr ReadMesh(MPI_Comm comm, const IoData &iodata) // conformal mesh refinement, or hexes. if (iodata.model.make_simplex || iodata.model.make_hex) { - SplitMeshElements(*smesh, iodata.model.make_simplex, iodata.model.make_hex); + SplitMeshElements(smesh, iodata.model.make_simplex, iodata.model.make_hex); } // Optionally reorder elements (and vertices) based on spatial location after loading @@ -1848,18 +1848,31 @@ mfem::Mesh MeshTetToHex(const mfem::Mesh &orig_mesh) return hex_mesh; } -void SplitMeshElements(mfem::Mesh &orig_mesh, bool make_simplex, bool make_hex, - bool preserve_curvature) +void SplitMeshElements(std::unique_ptr &orig_mesh, bool make_simplex, + bool make_hex, bool preserve_curvature) { if (!make_simplex && !make_hex) { return; } - MFEM_VERIFY(!orig_mesh.Nonconforming(), - "Mesh element splitting is not supported for nonconforming meshes!"); - mfem::Mesh *mesh = &orig_mesh; + mfem::Mesh *mesh = orig_mesh.get(); mfem::Mesh new_mesh; + // In order to track the mapping from parent mesh element to split mesh elements for + // high-order node interpolation, we use a unique attribute per parent mesh element. No + // need to call mfem::Mesh::SetAttributes at the end, since we won't use the attribute + // list data structure. + mfem::Array orig_mesh_attr; + auto BackUpAttributes = [&orig_mesh_attr](mfem::Mesh &mesh) + { + orig_mesh_attr.SetSize(mesh.GetNE()); + for (int e = 0; e < mesh.GetNE(); e++) + { + orig_mesh_attr[e] = mesh.GetAttribute(e); + mesh.SetAttribute(e, 1 + e); + } + }; + // Convert all element types to simplices. if (make_simplex) { @@ -1867,10 +1880,19 @@ void SplitMeshElements(mfem::Mesh &orig_mesh, bool make_simplex, bool make_hex, if (element_types.has_hexahedra || element_types.has_prisms || element_types.has_pyramids) { + MFEM_VERIFY(!mesh->Nonconforming(), + "Mesh element splitting is not supported for nonconforming meshes!"); MFEM_VERIFY( !element_types.has_pyramids, "Splitting mesh elements to simplices does not support pyramid elements yet!"); + if (preserve_curvature && mesh->GetNodes() && !orig_mesh_attr.Size()) + { + BackUpAttributes(*mesh); + } + int ne = mesh->GetNE(); new_mesh = mfem::Mesh::MakeSimplicial(*mesh); + Mpi::Print("Added {:d} elements to the mesh during conversion to simplices\n", + new_mesh.GetNE() - ne); mesh = &new_mesh; } } @@ -1882,40 +1904,153 @@ void SplitMeshElements(mfem::Mesh &orig_mesh, bool make_simplex, bool make_hex, if (element_types.has_simplices || element_types.has_prisms || element_types.has_pyramids) { + MFEM_VERIFY(!mesh->Nonconforming(), + "Mesh element splitting is not supported for nonconforming meshes!"); MFEM_VERIFY(!element_types.has_prisms && !element_types.has_pyramids, "Splitting mesh elements to hexahedra only supports simplex elements " "(tetrahedra) for now!"); + if (preserve_curvature && mesh->GetNodes() && !orig_mesh_attr.Size()) + { + BackUpAttributes(*mesh); + } + int ne = mesh->GetNE(); new_mesh = MeshTetToHex(*mesh); + Mpi::Print("Added {:d} elements to the mesh during conversion to hexahedra\n", + new_mesh.GetNE() - ne); 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()) + // Return if no modifications were made. + if (mesh == orig_mesh.get()) { - // Back up the original high-order nodes. - mfem::GridFunction *nodes = nullptr; - int own_nodes = 1; - orig_mesh.SwapNodes(nodes, own_nodes); + return; + } - // Prepare to interpolate the grid function for high-order nodes from the old mesh - // on the new one. - orig_mesh.EnsureNodes(); + // The previous splitting functions remove curvature information from the new mesh. So, if + // needed, we interpolate it onto the new mesh with GSLIB. The topology of the mesh must + // be finalized for this to work (assumes FinalizeTopology has been called on the new + // mesh). + if (preserve_curvature && orig_mesh->GetNodes()) + { + // Prepare to interpolate the grid function for high-order nodes from the old mesh to + // the new one. This first sets up the new mesh as a linear mesh and constructs a + // separate high-order grid function for storing the interpolated nodes. new_mesh.EnsureNodes(); + const mfem::GridFunction *nodes = orig_mesh->GetNodes(); const mfem::FiniteElementSpace *fespace = nodes->FESpace(); mfem::Ordering::Type ordering = fespace->GetOrdering(); int order = fespace->GetMaxElementOrder(); - int sdim = orig_mesh.SpaceDimension(); + int sdim = orig_mesh->SpaceDimension(); bool discont = (dynamic_cast(fespace->FEColl()) != nullptr); mfem::FiniteElementSpace new_fespace(&new_mesh, fespace->FEColl(), sdim, ordering); mfem::GridFunction new_nodes(&new_fespace); - fem::InterpolateFunction(*nodes, new_nodes); - if (own_nodes) + + mfem::Array vdofs; + mfem::Vector vals, elem_vals, xyz; + mfem::DenseMatrix pointmat; + int start = 0; + for (int e = 0; e < orig_mesh->GetNE(); e++) { - delete nodes; + // Get the high-order nodes restricted to this parent element, always returned + // byNODES. + fespace->GetElementVDofs(e, vdofs); + nodes->GetSubVector(vdofs, vals); + + // Find all child elements of this parent in the split mesh and restore their correct + // original attribute. + int attr = new_mesh.GetAttribute(start); + int end = start; + while (end < new_mesh.GetNE() && new_mesh.GetAttribute(end) == attr) + { + new_mesh.SetAttribute(end++, orig_mesh_attr[e]); + } + + // Special case: Parent element is unsplit, no interpolation is needed. + if (end == start + 1) + { + new_fespace.GetElementVDofs(start, vdofs); + new_nodes.SetSubVector(vdofs, vals); + start = end; + continue; + } + + // Create a list of points at which to interpolate the high order node information. + int npts = 0, offset = 0; + for (int i = start; i < end; i++) + { + npts += new_fespace.GetFE(i)->GetNodes().GetNPoints(); + } + xyz.SetSize(npts * sdim); + for (int i = start; i < end; i++) + { + const mfem::FiniteElement &fe = *new_fespace.GetFE(i); + mfem::ElementTransformation &T = *new_mesh.GetElementTransformation(i); + T.Transform(fe.GetNodes(), pointmat); + for (int j = 0; j < pointmat.Width(); j++) + { + for (int d = 0; d < sdim; d++) + { + // Use default ordering byNODES. + xyz(d * npts + offset + j) = pointmat(d, j); + } + } + offset += pointmat.Width(); + } + + // Create a single element linear mesh for interpolation. + const mfem::Element *el = orig_mesh->GetElement(e); + mfem::Mesh parent_mesh(orig_mesh->Dimension(), el->GetNVertices(), 1, 0, + orig_mesh->SpaceDimension()); + parent_mesh.AddElement(el->Duplicate(&parent_mesh)); + for (int i = 0; i < el->GetNVertices(); i++) + { + int v = parent_mesh.AddVertex(orig_mesh->GetVertex(el->GetVertices()[i])); + parent_mesh.GetElement(0)->GetVertices()[i] = v; + } + constexpr bool generate_bdr = false; + parent_mesh.FinalizeTopology(generate_bdr); + + // Create a grid function on the single element parent mesh with the high-order nodal + // grid function to interpolate. + parent_mesh.EnsureNodes(); + mfem::FiniteElementSpace parent_fespace(&parent_mesh, fespace->FEColl(), sdim, + mfem::Ordering::byNODES); + mfem::GridFunction parent_nodes(&parent_fespace); + MFEM_ASSERT(parent_nodes.Size() == vals.Size(), + "Unexpected size mismatch for high-order node interpolation!"); + parent_nodes = vals; + + // Interpolate the high-order nodes grid function and copy into the new mesh nodes + // grid function. + vals.SetSize(npts * sdim); + fem::InterpolateFunction(xyz, parent_nodes, vals, mfem::Ordering::byNODES); + offset = 0; + for (int i = start; i < end; i++) + { + const int elem_npts = new_fespace.GetFE(i)->GetNodes().GetNPoints(); + elem_vals.SetSize(elem_npts * sdim); + for (int d = 0; d < sdim; d++) + { + for (int j = 0; j < elem_npts; j++) + { + // Arrange element values byNODES to align with GetElementVDofs. + elem_vals(d * elem_npts + j) = vals(d * npts + offset + j); + } + } + new_fespace.GetElementVDofs(i, vdofs); + MFEM_ASSERT(vdofs.Size() == elem_vals.Size(), + "Unexpected size mismatch for high-order node interpolation!"); + new_nodes.SetSubVector(vdofs, elem_vals); + offset += elem_npts; + } + + // Prepare for next parent element. + start = end; } + MFEM_VERIFY(start == new_mesh.GetNE(), + "Premature abort for high-order curvature in mesh splitting!"); // Finally, copy the nodal grid function to the new mesh. new_mesh.SetCurvature(order, discont, sdim, ordering); @@ -1923,7 +2058,7 @@ void SplitMeshElements(mfem::Mesh &orig_mesh, bool make_simplex, bool make_hex, "Unexpected size mismatch for nodes!"); new_mesh.SetNodes(new_nodes); } - orig_mesh = std::move(new_mesh); + orig_mesh = std::make_unique(std::move(new_mesh)); // Call move constructor } void ReorderMeshElements(mfem::Mesh &mesh, bool print) @@ -2252,7 +2387,7 @@ std::map> CheckMesh(std::unique_ptr &orig_me // The element loop works because we know the mapping from old_mesh to new_mesh element // indices from the insertion order. - mfem::Array vdofs, new_vdofs; + mfem::Array vdofs; mfem::Vector loc_vec; int te = 0; for (int e = 0; e < orig_mesh->GetNE(); e++) @@ -2262,8 +2397,8 @@ std::map> CheckMesh(std::unique_ptr &orig_me // No need for DofTransformation here since spaces are H1 or L2. fespace->GetElementVDofs(e, vdofs); nodes->GetSubVector(vdofs, loc_vec); - new_fespace->GetElementVDofs(te, new_vdofs); - new_nodes->SetSubVector(new_vdofs, loc_vec); + new_fespace->GetElementVDofs(te, vdofs); + new_nodes->SetSubVector(vdofs, loc_vec); te++; } } From f3192000218bf4771daaee66cc49e1d5fed90f50 Mon Sep 17 00:00:00 2001 From: Sebastian Grimberg Date: Wed, 19 Jun 2024 18:14:31 -0700 Subject: [PATCH 11/14] Fix bug for certain MPI distributions --- palace/utils/geodata.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/palace/utils/geodata.cpp b/palace/utils/geodata.cpp index 729228ea6..e31e1bd60 100644 --- a/palace/utils/geodata.cpp +++ b/palace/utils/geodata.cpp @@ -83,7 +83,7 @@ std::unique_ptr ReadMesh(MPI_Comm comm, const IoData &iodata) const auto &refinement = iodata.model.refinement; const bool use_amr = refinement.max_it > 0; const bool use_mesh_partitioner = !use_amr || !refinement.nonconformal; - MPI_Comm node_comm = nullptr; + MPI_Comm node_comm; if (!use_mesh_partitioner) { MPI_Comm_split_type(comm, MPI_COMM_TYPE_SHARED, Mpi::Rank(comm), MPI_INFO_NULL, @@ -93,7 +93,8 @@ std::unique_ptr ReadMesh(MPI_Comm comm, const IoData &iodata) // Only one process per node reads the serial mesh. { BlockTimer bt(Timer::IO); - if ((use_mesh_partitioner && Mpi::Root(comm)) || Mpi::Root(node_comm)) + if ((use_mesh_partitioner && Mpi::Root(comm)) || + (!use_mesh_partitioner && Mpi::Root(node_comm))) { smesh = LoadMesh(iodata.model.mesh, iodata.model.remove_curvature); MFEM_VERIFY(!(smesh->Nonconforming() && use_mesh_partitioner), From 54d6631cc40158b44d8950c2b9cce323ff9e32dd Mon Sep 17 00:00:00 2001 From: Sebastian Grimberg Date: Fri, 21 Jun 2024 11:41:17 -0700 Subject: [PATCH 12/14] MakeHex -> MakeHexahedral --- docs/src/config/model.md | 2 +- palace/utils/configfile.cpp | 6 +++--- scripts/schema/config/model.json | 2 +- 3 files changed, 5 insertions(+), 5 deletions(-) diff --git a/docs/src/config/model.md b/docs/src/config/model.md index 3a7da2d88..9f644bf8e 100644 --- a/docs/src/config/model.md +++ b/docs/src/config/model.md @@ -115,7 +115,7 @@ mesh length units. - `"RemoveCurvature" [false]` - `"MakeSimplex" [false]` - - `"MakeHex" [false]` + - `"MakeHexahedral" [false]` - `"ReorderElements" [false]` - `"CleanUnusedElements" [true]` - `"AddInterfaceBoundaryElements" [true]` diff --git a/palace/utils/configfile.cpp b/palace/utils/configfile.cpp index daddd69d7..62ab01fd3 100644 --- a/palace/utils/configfile.cpp +++ b/palace/utils/configfile.cpp @@ -467,7 +467,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); + make_hex = model->value("MakeHexahedral", 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); @@ -481,7 +481,7 @@ void ModelData::SetUp(json &config) model->erase("Lc"); model->erase("RemoveCurvature"); model->erase("MakeSimplex"); - model->erase("MakeHex"); + model->erase("MakeHexahedral"); model->erase("ReorderElements"); model->erase("CleanUnusedElements"); model->erase("AddInterfaceBoundaryElements"); @@ -500,7 +500,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 << "MakeHexahedral: " << 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/scripts/schema/config/model.json b/scripts/schema/config/model.json index cab779794..9894835ae 100644 --- a/scripts/schema/config/model.json +++ b/scripts/schema/config/model.json @@ -11,7 +11,7 @@ "Lc": { "type": "number", "exclusiveMinimum": 0.0 }, "RemoveCurvature": { "type": "boolean" }, "MakeSimplex": { "type": "boolean" }, - "MakeHex": { "type": "boolean" }, + "MakeHexahedral": { "type": "boolean" }, "ReorderElements": { "type": "boolean" }, "CleanUnusedElements": { "type": "boolean" }, "AddInterfaceBoundaryElements": { "type": "boolean" }, From 5053453caa4ec364951d616300a14fbea28c13cd Mon Sep 17 00:00:00 2001 From: Sebastian Grimberg Date: Fri, 21 Jun 2024 11:43:19 -0700 Subject: [PATCH 13/14] Update changelog --- CHANGELOG.md | 3 +++ 1 file changed, 3 insertions(+) diff --git a/CHANGELOG.md b/CHANGELOG.md index 88275bef2..600a50b9e 100644 --- a/CHANGELOG.md +++ b/CHANGELOG.md @@ -22,6 +22,9 @@ The format of this changelog is based on `"WavePortPEC"`). - Fixed a bug in divergence-free projection for problems without essential or mixed boundary conditions. + - Added `"MakeSimplex"` and `"MakeHexahedral"` mesh options to convert an input mesh to + all tetrahedra or all hexahedra. Also adds `"SerialUniformLevels"` option to + `config["Model"]["Refinement"]` for testing or debugging. ## [0.13.0] - 2024-05-20 From 3e03a64504d2660ff5fd72bf0b09527f6d700caf Mon Sep 17 00:00:00 2001 From: Sebastian Grimberg Date: Wed, 26 Jun 2024 10:27:56 -0700 Subject: [PATCH 14/14] Address PR feedback --- palace/utils/geodata.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/palace/utils/geodata.cpp b/palace/utils/geodata.cpp index e31e1bd60..e69f2e049 100644 --- a/palace/utils/geodata.cpp +++ b/palace/utils/geodata.cpp @@ -194,13 +194,12 @@ std::unique_ptr ReadMesh(MPI_Comm comm, const IoData &iodata) so.resize(slen); } Mpi::Broadcast(slen, so.data(), 0, node_comm); - if (!smesh) { std::istringstream fi(so); // std::istringstream fi(zlib::DecompressString(so)); smesh = std::make_unique(fi, generate_edges, refine, fix_orientation); + so.clear(); } - so.clear(); if (refinement.nonconformal && use_amr) { smesh->EnsureNCMesh(true); @@ -1960,7 +1959,8 @@ void SplitMeshElements(std::unique_ptr &orig_mesh, bool make_simplex nodes->GetSubVector(vdofs, vals); // Find all child elements of this parent in the split mesh and restore their correct - // original attribute. + // original attribute. This works because child elements are added in contiguous + // batches in the same order as the parents. int attr = new_mesh.GetAttribute(start); int end = start; while (end < new_mesh.GetNE() && new_mesh.GetAttribute(end) == attr) @@ -1989,9 +1989,9 @@ void SplitMeshElements(std::unique_ptr &orig_mesh, bool make_simplex const mfem::FiniteElement &fe = *new_fespace.GetFE(i); mfem::ElementTransformation &T = *new_mesh.GetElementTransformation(i); T.Transform(fe.GetNodes(), pointmat); - for (int j = 0; j < pointmat.Width(); j++) + for (int d = 0; d < sdim; d++) { - for (int d = 0; d < sdim; d++) + for (int j = 0; j < pointmat.Width(); j++) { // Use default ordering byNODES. xyz(d * npts + offset + j) = pointmat(d, j);