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 diff --git a/docs/src/config/model.md b/docs/src/config/model.md index 1c3c3078d..9f644bf8e 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,10 +113,16 @@ mesh length units. ### Advanced model options - - `"Partition" [""]` - - `"ReorientTetMesh" [false]` - `"RemoveCurvature" [false]` + - `"MakeSimplex" [false]` + - `"MakeHexahedral" [false]` + - `"ReorderElements" [false]` + - `"CleanUnusedElements" [true]` + - `"AddInterfaceBoundaryElements" [true]` + - `"ReorientTetMesh" [false]` + - `"Partitioning" [""]` - `"MaxNCLevels" [1]` - `"MaximumImbalance" [1.1]` - `"SaveAdaptIterations" [true]` - `"SaveAdaptMesh" [false]` + - `"SerialUniformLevels" [0]` 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": { diff --git a/palace/fem/interpolator.cpp b/palace/fem/interpolator.cpp index a24abc25f..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,27 +32,24 @@ 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()); + 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) { - // Use default ordering byNODES. - xyz(i) = data.center[0]; - xyz(npts + i) = data.center[1]; - if (mesh.SpaceDimension() == 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; @@ -69,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 dim = U.VectorDim(); - std::vector vals(npts * dim); - mfem::Vector v(npts * dim); - op.Interpolate(U, v); - for (int d = 0; d < dim; d++) + const int vdim = U.VectorDim(); + std::vector vals(npts * vdim); + if (U.FESpace()->GetOrdering() == mfem::Ordering::byVDIM) { - for (int i = 0; i < npts; i++) + mfem::Vector v(vals.data(), npts * vdim); + op.Interpolate(U, v); + } + else + { + mfem::Vector v(npts * vdim); + op.Interpolate(U, v); + for (int d = 0; d < vdim; d++) { - vals[i * dim + d] = v(d * npts + i); + for (int i = 0; i < npts; i++) + { + vals[i * vdim + d] = v(d * npts + i); + } } } return vals; @@ -107,4 +120,184 @@ 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. 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(); + mfem::Vector xyz; + mfem::Ordering::Type ordering; + const auto *dest_fec_h1 = + dynamic_cast(V.FESpace()->FEColl()); + const auto *dest_nodes_h1 = 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_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 = *dest_mesh.GetElementTransformation(i); + T.Transform(fe.GetNodes(), pointmat); + for (int j = 0; j < pointmat.Width(); j++) + { + for (int d = 0; d < dim; d++) + { + // Use default ordering byNODES. + xyz(d * npts + offset + j) = pointmat(d, 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); + 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 = U.VectorDim(); + mfem::Vector vals(npts * vdim); + 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) + { + // H1 with order != mesh order needs to handle duplicated interpolation 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 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) * vdim + d; + elem_vals(d * elem_npts + j) = vals(idx); + } + } + 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 + { + // 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 + { + // 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 = *dest_mesh.GetElementTransformation(i); + const int elem_npts = fe.GetNodes().GetNPoints(); + elem_vals.SetSize(elem_npts * vdim); + 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) * vdim + d; + elem_vals(j * vdim + d) = vals(idx); + } + } + 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); + offset += elem_npts; + } + } +#else + MFEM_ABORT("InterpolateFunction requires MFEM_USE_GSLIB!"); +#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 460caa2e9..468512f20 100644 --- a/palace/fem/interpolator.hpp +++ b/palace/fem/interpolator.hpp @@ -35,6 +35,21 @@ 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); + +// 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 #endif // PALACE_FEM_INTERPOLATOR_HPP diff --git a/palace/main.cpp b/palace/main.cpp index ac3cfebf2..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, false)); + 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..62ab01fd3 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'; } } @@ -462,18 +465,28 @@ 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); + make_simplex = model->value("MakeSimplex", make_simplex); + 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); + 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("MakeSimplex"); + model->erase("MakeHexahedral"); + 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 +498,14 @@ 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 << "MakeSimplex: " << make_simplex << '\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'; + 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..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 = {}; @@ -199,14 +202,32 @@ 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; + + // 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; + + // 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 487ef4c80..e69f2e049 100644 --- a/palace/utils/geodata.cpp +++ b/palace/utils/geodata.cpp @@ -9,10 +9,12 @@ #include #include #include -#include #include #include +#include +#include #include +#include "fem/interpolator.hpp" #include "utils/communication.hpp" #include "utils/diagnostic.hpp" #include "utils/filesystem.hpp" @@ -38,20 +40,24 @@ 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(std::unique_ptr &, 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); +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. +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); -// 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. -std::map> CheckMesh(std::unique_ptr &, const IoData &, - bool, bool, bool, const int *); - // 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 &, @@ -66,8 +72,7 @@ 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) { // 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. @@ -78,9 +83,18 @@ std::unique_ptr ReadMesh(MPI_Comm comm, const IoData &iodata, boo 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; + 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)) || + (!use_mesh_partitioner && Mpi::Root(node_comm))) { smesh = LoadMesh(iodata.model.mesh, iodata.model.remove_curvature); MFEM_VERIFY(!(smesh->Nonconforming() && use_mesh_partitioner), @@ -89,44 +103,66 @@ std::unique_ptr ReadMesh(MPI_Comm comm, const IoData &iodata, boo 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); - 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, - "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!"); + 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, or hexes. + if (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 // the serial mesh. - if (reorder) + if (iodata.model.reorder_elements) { - ReorderMesh(*smesh); + 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. - if ((clean_elem || add_bdr) && smesh->Conforming()) + if (smesh->Conforming()) { - static_cast(CheckMesh(smesh, iodata, clean_elem, add_bdr, false, nullptr)); + 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); - - // Optionally add subdomain interface boundary elements. - if (add_subdomain && smesh->Conforming()) - { - static_cast( - CheckMesh(smesh, iodata, false, false, add_subdomain, partitioning.get())); - } + partitioning = GetMeshPartitioning(*smesh, Mpi::Size(comm), iodata.model.partitioning); } + // Distribute the mesh. std::unique_ptr pmesh; if (use_mesh_partitioner) { @@ -134,10 +170,46 @@ std::unique_ptr ReadMesh(MPI_Comm comm, const IoData &iodata, boo } 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); + { + 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(); } @@ -212,7 +284,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 @@ -360,7 +432,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 @@ -739,7 +811,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) { @@ -1530,7 +1602,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); @@ -1570,7 +1642,6 @@ double RebalanceMesh(std::unique_ptr &mesh, const IoData &iodata) } if (ratio > tol) { - mesh->ExchangeFaceNbrData(); if (mesh->Nonconforming()) { mesh->Rebalance(); @@ -1581,7 +1652,6 @@ double RebalanceMesh(std::unique_ptr &mesh, const IoData &iodata) // processor and then redistributed. RebalanceConformalMesh(mesh); } - mesh->ExchangeFaceNbrData(); } return ratio; } @@ -1594,7 +1664,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 @@ -1603,61 +1673,34 @@ 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); // 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() != '/') + if (mesh_path.extension() == ".mphtxt" || mesh_path.extension() == ".mphbin") { - 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 (mfile.extension() == ".mphtxt" || mfile.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 - 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 { // 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); } @@ -1681,7 +1724,345 @@ std::unique_ptr LoadMesh(const std::string &path, bool remove_curvat return mesh; } -void ReorderMesh(mfem::Mesh &mesh, bool print) +mfem::Mesh MeshTetToHex(const mfem::Mesh &orig_mesh) +{ + // 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!"); + } + + // 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++) + { + 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++) + { + for (int d = 0; d < orig_mesh.SpaceDimension(); d++) + { + coord[d] += orig_mesh.GetVertex(verts[i])[d] / nv; + } + } + hex_mesh.AddVertex(coord); + }; + { + mfem::Array verts; + for (int e = 0; e < nedge_tet; ++e) + { + 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()); + } + + // 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 SplitMeshElements(std::unique_ptr &orig_mesh, bool make_simplex, + bool make_hex, bool preserve_curvature) +{ + if (!make_simplex && !make_hex) + { + return; + } + 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) + { + const auto element_types = mesh::CheckElements(*mesh); + 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; + } + } + + // 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(!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; + } + } + + // Return if no modifications were made. + if (mesh == orig_mesh.get()) + { + return; + } + + // 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(); + bool discont = + (dynamic_cast(fespace->FEColl()) != nullptr); + mfem::FiniteElementSpace new_fespace(&new_mesh, fespace->FEColl(), sdim, ordering); + mfem::GridFunction new_nodes(&new_fespace); + + mfem::Array vdofs; + mfem::Vector vals, elem_vals, xyz; + mfem::DenseMatrix pointmat; + int start = 0; + for (int e = 0; e < orig_mesh->GetNE(); e++) + { + // 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. 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) + { + 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 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); + } + } + 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); + MFEM_VERIFY(new_mesh.GetNodes()->Size() == new_nodes.Size(), + "Unexpected size mismatch for nodes!"); + new_mesh.SetNodes(new_nodes); + } + orig_mesh = std::make_unique(std::move(new_mesh)); // Call move constructor +} + +void ReorderMeshElements(mfem::Mesh &mesh, bool print) { mfem::Array ordering; @@ -1713,64 +2094,9 @@ void ReorderMesh(mfem::Mesh &mesh, bool print) mesh.ReorderElements(ordering); } -std::unique_ptr GetMeshPartitioning(const mfem::Mesh &mesh, int size, - const std::string &partition, bool print) -{ - MFEM_VERIFY(size <= mesh.GetNE(), "Mesh partitioning must have parts <= mesh elements (" - << size << " vs. " << mesh.GetNE() << ")!"); - if (partition.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(partition); - 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::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 +2116,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 +2148,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 +2185,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 +2204,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 +2286,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 +2294,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 +2303,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 +2341,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); } } } @@ -2101,19 +2377,18 @@ 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(); // 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++) @@ -2123,8 +2398,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++; } } @@ -2134,126 +2409,120 @@ 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, 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) diff --git a/palace/utils/geodata.hpp b/palace/utils/geodata.hpp index 9c9ef3c38..285ebd52f 100644 --- a/palace/utils/geodata.hpp +++ b/palace/utils/geodata.hpp @@ -24,8 +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, - bool clean_elem, bool add_bdr, bool add_subdomain); +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..9894835ae 100644 --- a/scripts/schema/config/model.json +++ b/scripts/schema/config/model.json @@ -9,9 +9,14 @@ "Mesh": { "type": "string" }, "L0": { "type": "number", "exclusiveMinimum": 0.0 }, "Lc": { "type": "number", "exclusiveMinimum": 0.0 }, - "Partition": { "type": "string" }, - "ReorientTetMesh": { "type": "boolean" }, "RemoveCurvature": { "type": "boolean" }, + "MakeSimplex": { "type": "boolean" }, + "MakeHexahedral": { "type": "boolean" }, + "ReorderElements": { "type": "boolean" }, + "CleanUnusedElements": { "type": "boolean" }, + "AddInterfaceBoundaryElements": { "type": "boolean" }, + "ReorientTetMesh": { "type": "boolean" }, + "Partitioning": { "type": "string" }, "Refinement": { "type": "object", @@ -29,6 +34,7 @@ "SaveAdaptIterations": {"type": "boolean"}, "SaveAdaptMesh": {"type": "boolean"}, "UniformLevels": { "type": "integer", "minimum": 0 }, + "SerialUniformLevels": { "type": "integer", "minimum": 0 }, "Boxes": { "type": "array",