diff --git a/palace/utils/configfile.cpp b/palace/utils/configfile.cpp index d148e32f0..37d368940 100644 --- a/palace/utils/configfile.cpp +++ b/palace/utils/configfile.cpp @@ -289,10 +289,61 @@ void RefinementData::SetUp(json &model) } } + auto adapt = refinement->find("Adaptation"); + if (adapt != refinement->end()) + { + MFEM_ABORT("Placeholder: Not currently supported!"); + + // Load Values + adaptation.tolerance = adapt->value("Tol", adaptation.tolerance); + adaptation.max_its = adapt->value("MaxIts", adaptation.max_its); + adaptation.min_its = adapt->value("MinIts", adaptation.min_its); + adaptation.update_fraction = adapt->value("UpdateFraction", adaptation.update_fraction); + adaptation.max_nc_levels = adapt->value("MaxNCLevels", adaptation.max_nc_levels); + adaptation.dof_limit = adapt->value("DOFLimit", adaptation.dof_limit); + adaptation.coarsening_fraction = + adapt->value("CoarseningFraction", adaptation.coarsening_fraction); + adaptation.save_step = adapt->value("SaveStep", adaptation.save_step); + adaptation.nonconformal = adapt->value("Nonconformal", adaptation.nonconformal); + adaptation.maximum_imbalance = + adapt->value("MaximumImbalance", adaptation.maximum_imbalance); + + // Perform Checks + MFEM_VERIFY(adaptation.tolerance > 0, "\"Tol\" must be strictly positive"); + MFEM_VERIFY(adaptation.max_its >= 0, "\"MaxIts\" must be non-negative"); + MFEM_VERIFY(adaptation.min_its >= 0, "\"MinIts\" must be non-negative"); + MFEM_VERIFY(adaptation.min_its <= adaptation.max_its, + "\"MinIts\" must be smaller than \"MaxIts\": " << adaptation.min_its << "," + << adaptation.max_its); + MFEM_VERIFY(adaptation.update_fraction > 0 && adaptation.update_fraction < 1, + "\"UpdateFraction\" must be in (0,1)"); + MFEM_VERIFY(adaptation.coarsening_fraction >= 0 && adaptation.coarsening_fraction < 1, + "\"CoarseningFraction\" must be in [0, 1)"); + MFEM_VERIFY(adaptation.max_nc_levels >= 0, "\"MaxNCLevels\" must non-negative"); + MFEM_VERIFY(adaptation.dof_limit >= 0, "\"DOFLimit\" must be non-negative"); + MFEM_VERIFY(adaptation.save_step >= 0, "\"SaveStep\" must be non-negative"); + MFEM_VERIFY(adaptation.maximum_imbalance >= 1, + "\"MaximumImbalance\" must be greater than or equal to 1"); + + // Cleanup + const auto fields = { + "Tol", "MaxIts", "MinIts", "UpdateFraction", "CoarseningFraction", + "MaxNCLevels", "DOFLimit", "SaveStep", "Nonconformal", "MaximumImbalance"}; + for (const auto &f : fields) + { + adapt->erase(f); + } + + MFEM_VERIFY(adapt->empty(), + "Found an unsupported configuration file keyword under \"Adaptation\"!\n" + << adapt->dump(2)); + } + // Cleanup refinement->erase("UniformLevels"); refinement->erase("Boxes"); refinement->erase("Spheres"); + refinement->erase("Adaptation"); MFEM_VERIFY(refinement->empty(), "Found an unsupported configuration file keyword under \"Refinement\"!\n" << refinement->dump(2)); diff --git a/palace/utils/configfile.hpp b/palace/utils/configfile.hpp index 822b6083c..55242248b 100644 --- a/palace/utils/configfile.hpp +++ b/palace/utils/configfile.hpp @@ -107,11 +107,43 @@ struct SphereRefinementData std::vector center = {}; }; +// Stores data specifying the adaptive mesh refinement algorithm. +struct AdaptiveRefinementData +{ + // Non-dimensional tolerance used to specify convergence of the AMR. + double tolerance = 1e-2; + // Maximum number of iterations to perform during the AMR. + int max_its = 0; + // Minimum number of iterations to perform during the AMR. + int min_its = 0; + // Dörfler update fraction. The set of marked elements is the minimum set + // that contains update_fraction of the total error. + double update_fraction = 0.4; + // Whether or not to perform coarsening during the AMR. + double coarsening_fraction = 0.0; + // Maximum difference in non-conformal refinements between two adjacent + // elements. Default = 0 implies there is no constraint on local non-conformity. + int max_nc_levels = 0; + // If a refinement results in a greater number of DOFs than this value, no + // future refinement will be allowed unless coarsening is allowed to occur. + int dof_limit = 0; + // Frequency with which to store the post processing results for a given + // adaptation, e.g. save_step = 3 means save every third adaptation. + int save_step = 0; + // Whether or not to perform nonconformal adaptation. + bool nonconformal = true; + // Maximum allowable ratio of number of elements across processors before + // rebalancing is performed. + double maximum_imbalance = 1.5; +}; + struct RefinementData { public: // Parallel uniform mesh refinement levels. int uniform_ref_levels = 0; + // Adaptive refinement configuration data. + AdaptiveRefinementData adaptation; private: // Refinement data for mesh regions. diff --git a/palace/utils/geodata.cpp b/palace/utils/geodata.cpp index 0e1014654..a1e0420ca 100644 --- a/palace/utils/geodata.cpp +++ b/palace/utils/geodata.cpp @@ -22,34 +22,51 @@ namespace // Floating point precision for mesh IO. This precision is important, make sure nothing is // lost! -const auto MSH_FLT_PRECISION = std::numeric_limits::max_digits10; +constexpr auto MSH_FLT_PRECISION = std::numeric_limits::max_digits10; // Load the serial mesh from disk. -std::unique_ptr LoadMesh(const std::string &); +std::unique_ptr LoadMesh(const std::string &path); // Optionally reorder mesh elements based on MFEM's internal reordeing tools for improved // cache usage. -void ReorderMesh(mfem::Mesh &); +void ReorderMesh(mfem::Mesh &mesh); // Generate element-based mesh partitioning, using either a provided file or METIS. -std::unique_ptr GetMeshPartitioning(mfem::Mesh &, int, const std::string &); +std::unique_ptr GetMeshPartitioning(mfem::Mesh &mesh, int size, + const std::string &partition = ""); + +// Given a serial CONFORMAL mesh on the root processor, and element +// partitioning, create a parallel mesh over the given communicator. +std::unique_ptr +DistributeMesh(MPI_Comm comm, std::unique_ptr &mesh, + const std::unique_ptr &partitioning = nullptr); // Cleanup 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 std::unique_ptr &, const IoData &, - bool, bool, bool); - -// Given a serial mesh on the root processor and element partitioning, create a parallel -// mesh oer the given communicator. -std::unique_ptr DistributeMesh(MPI_Comm, std::unique_ptr &, - std::unique_ptr &); +std::map> CheckMesh(mfem::Mesh &orig_mesh, + const std::unique_ptr &partitioning, + const IoData &iodata, bool clean_elem, + bool add_bdr, bool add_subdomain); // Get list of domain and boundary attribute markers used in configuration file for mesh // cleaning. -void GetUsedAttributeMarkers(const IoData &, int, int, mfem::Array &, - mfem::Array &); +void GetUsedAttributeMarkers(const IoData &iodata, int n_mat, int n_bdr, + mfem::Array &mat_marker, mfem::Array &bdr_marker); + +// Simplified helper for describing the element types in a mesh +struct ElementTypeInfo +{ + bool has_simplices; + bool has_tensors; + bool has_wedges; + bool has_pyramids; +}; +ElementTypeInfo CheckElements(mfem::Mesh &mesh) +{ + auto meshgen = mesh.MeshGenerator(); + return {bool(meshgen & 1), bool(meshgen & 2), bool(meshgen & 4), bool(meshgen & 8)}; +} } // namespace @@ -60,16 +77,24 @@ std::unique_ptr ReadMesh(MPI_Comm comm, const IoData &iodata, boo bool clean, bool add_bdr, bool unassembled, Timer &timer) { - // On root, read the serial mesh (converting format if necessary), and do all necessary - // serial preprocessing. When finished, distribute the mesh to all processes. Count disk - // I/O time separately for the mesh read from file. - std::unique_ptr smesh; + // 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. + // Count disk I/O time separately for the mesh read from file. auto t0 = timer.Now(); - if (Mpi::Root(comm)) + + // If not adapting, or performing conformal adaptation, can use the mesh partitioner. + std::unique_ptr smesh; + std::unique_ptr partitioning; + const auto use_amr = iodata.model.refinement.adaptation.max_its > 0; + const bool use_mesh_partitioner = + !use_amr || !iodata.model.refinement.adaptation.nonconformal; + if (Mpi::Root(comm) || !use_mesh_partitioner) { + // If using the mesh partitioner, only the root node needs to load the mesh. + smesh = LoadMesh(iodata.model.mesh); + // Optionally reorder elements (and vertices) based on spatial location after loading // the serial mesh. - smesh = LoadMesh(iodata.model.mesh); if (reorder) { ReorderMesh(*smesh); @@ -78,24 +103,40 @@ std::unique_ptr ReadMesh(MPI_Comm comm, const IoData &iodata, boo Mpi::Barrier(comm); timer.io_time += timer.Now() - t0; - std::unique_ptr partitioning; - if (Mpi::Root(comm)) + if (Mpi::Root(comm) || !use_mesh_partitioner) { - // Generate the parallel mesh partitioning on the root process. + const auto element_types = CheckElements(*smesh); + + MFEM_VERIFY(!use_amr || !element_types.has_tensors || + iodata.model.refinement.adaptation.nonconformal, + "If there are tensor elements, AMR must be nonconformal"); + MFEM_VERIFY(!use_amr || !element_types.has_pyramids || + iodata.model.refinement.adaptation.nonconformal, + "If there are pyramid elements, AMR must be nonconformal"); + MFEM_VERIFY(!use_amr || !element_types.has_wedges || + iodata.model.refinement.adaptation.nonconformal, + "If there are wedge elements, AMR must be nonconformal"); + partitioning = GetMeshPartitioning(*smesh, Mpi::Size(comm), iodata.model.partition); // Clean up unused domain elements from the mesh, add new boundary elements for material // interfaces if not present, and optionally (when running unassembled) add subdomain // interface boundary elements. - std::map> attr_map = - CheckMesh(smesh, partitioning, iodata, clean, add_bdr, unassembled); + static_cast(CheckMesh(*smesh, partitioning, iodata, clean, add_bdr, unassembled)); } - // Construct the parallel mesh data structure by distributing the serial mesh from the - // root process. The serial mesh and partitioning are deleted inside. - std::unique_ptr mesh = DistributeMesh(comm, smesh, partitioning); + std::unique_ptr pmesh; + if (iodata.model.refinement.adaptation.nonconformal && use_amr) + { + smesh->EnsureNCMesh(true); + pmesh = std::make_unique(comm, *smesh, partitioning.get()); + } + else + { + pmesh = DistributeMesh(comm, smesh, partitioning); + } -#if 0 + if constexpr (false) { std::string tmp = iodata.problem.output; if (tmp.back() != '/') @@ -107,28 +148,29 @@ std::unique_ptr ReadMesh(MPI_Comm comm, const IoData &iodata, boo { std::filesystem::create_directories(tmp); } - int width = 1 + static_cast(std::log10(Mpi::Size(comm)-1)); + int width = 1 + static_cast(std::log10(Mpi::Size(comm) - 1)); std::unique_ptr gsmesh = LoadMesh(iodata.model.mesh); std::unique_ptr gpartitioning = GetMeshPartitioning(*gsmesh, Mpi::Size(comm)); mfem::ParMesh gpmesh(comm, *gsmesh, gpartitioning.get(), 0); { - std::string pfile = mfem::MakeParFilename(tmp + "part.", Mpi::Rank(comm), ".mesh", width); + std::string pfile = + mfem::MakeParFilename(tmp + "part.", Mpi::Rank(comm), ".mesh", width); std::ofstream fo(pfile); // mfem::ofgzstream fo(pfile, true); // Use zlib compression if available fo.precision(MSH_FLT_PRECISION); gpmesh.ParPrint(fo); } { - std::string pfile = mfem::MakeParFilename(tmp + "final.", Mpi::Rank(comm), ".mesh", width); + std::string pfile = + mfem::MakeParFilename(tmp + "final.", Mpi::Rank(comm), ".mesh", width); std::ofstream fo(pfile); // mfem::ofgzstream fo(pfile, true); // Use zlib compression if available fo.precision(MSH_FLT_PRECISION); - mesh->ParPrint(fo); + pmesh->ParPrint(fo); } } -#endif - return mesh; + return pmesh; } void RefineMesh(const IoData &iodata, std::vector> &mesh) @@ -177,7 +219,7 @@ void RefineMesh(const IoData &iodata, std::vector { if (mesh.capacity() > 1) { - mesh.push_back(std::make_unique(*mesh.back())); + mesh.emplace_back(std::make_unique(*mesh.back())); } mesh.back()->UniformRefinement(); } @@ -296,7 +338,7 @@ void RefineMesh(const IoData &iodata, std::vector // (adds hanging nodes). if (mesh.capacity() > 1) { - mesh.push_back(std::make_unique(*mesh.back())); + mesh.emplace_back(std::make_unique(*mesh.back())); } mesh.back()->GeneralRefinement(refs, -1); region_ref_level++; @@ -670,7 +712,7 @@ std::unique_ptr LoadMesh(const std::string &path) { MFEM_ABORT("Unable to open mesh file \"" << path << "\"!"); } - std::unique_ptr mesh = std::make_unique(fi, 1, 1, true); + auto mesh = std::make_unique(fi, 1, 1, true); mesh->EnsureNodes(); return mesh; } @@ -679,25 +721,6 @@ void ReorderMesh(mfem::Mesh &mesh) { mfem::Array ordering; -#if 0 - // Gecko reordering. - Mpi::Print(\n); - 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, eriod, seed, true); - if (cost < best_cost) - { - ordering = tentative; - best_cost = cost; - } - } - Mpi::Print("Final cost: {:e}\n", best_cost); -#endif - // (Faster) Hilbert reordering. mesh.GetHilbertElementOrdering(ordering); mesh.ReorderElements(ordering); @@ -750,7 +773,7 @@ std::unique_ptr GetMeshPartitioning(mfem::Mesh &mesh, int size, return partitioning; } -std::map> CheckMesh(std::unique_ptr &orig_mesh, +std::map> CheckMesh(mfem::Mesh &orig_mesh, const std::unique_ptr &partitioning, const IoData &iodata, bool clean_elem, bool add_bdr, bool add_subdomain) @@ -761,20 +784,22 @@ std::map> CheckMesh(std::unique_ptr &orig_me // interfaces if these elements do not yet exist. // - If desired, create a new mesh which has removed all domain elements which do not have // an associated material property specified in the input file. - MFEM_VERIFY(orig_mesh->Dimension() == 3 && !orig_mesh->Nonconforming(), + MFEM_VERIFY(orig_mesh.Dimension() == 3 && !orig_mesh.Nonconforming(), "Nonconforming or 2D meshes have not been tested yet!"); + MFEM_VERIFY(dynamic_cast(&orig_mesh) == nullptr, + "This function does not work for ParMesh"); mfem::Array mat_marker, bdr_marker; - GetUsedAttributeMarkers(iodata, orig_mesh->attributes.Max(), - orig_mesh->bdr_attributes.Max(), mat_marker, bdr_marker); + GetUsedAttributeMarkers(iodata, orig_mesh.attributes.Max(), + orig_mesh.bdr_attributes.Max(), mat_marker, bdr_marker); bool warn = false; - for (int be = 0; be < orig_mesh->GetNBE(); be++) + for (int be = 0; be < orig_mesh.GetNBE(); be++) { - int attr = orig_mesh->GetBdrAttribute(be); + int attr = orig_mesh.GetBdrAttribute(be); if (!bdr_marker[attr - 1]) { int f, o, e1, e2; - orig_mesh->GetBdrElementFace(be, &f, &o); - orig_mesh->GetFaceElements(f, &e1, &e2); + orig_mesh.GetBdrElementFace(be, &f, &o); + orig_mesh.GetFaceElements(f, &e1, &e2); if (e1 < 0 || e2 < 0) // Internal boundary elements are allowed to have no BC { warn = true; @@ -797,33 +822,33 @@ std::map> CheckMesh(std::unique_ptr &orig_me } // Count deleted or added domain and boundary elements. - int new_ne = orig_mesh->GetNE(); - int new_nbdr = orig_mesh->GetNBE(); + 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); - for (int be = 0; be < orig_mesh->GetNBE(); be++) + elem_delete.SetSize(orig_mesh.GetNE(), false); + bdr_delete.SetSize(orig_mesh.GetNBE(), false); + orig_bdr_faces.SetSize(orig_mesh.GetNumFaces(), -1); + for (int be = 0; be < orig_mesh.GetNBE(); be++) { int f, o; - orig_mesh->GetBdrElementFace(be, &f, &o); + orig_mesh.GetBdrElementFace(be, &f, &o); MFEM_VERIFY(orig_bdr_faces[f] < 0, "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); + add_bdr_faces.SetSize(orig_mesh.GetNumFaces(), -1); } if (clean_elem) { // Delete domain and boundary elements which have no associated material or BC attribute // from the mesh. - for (int e = 0; e < orig_mesh->GetNE(); e++) + for (int e = 0; e < orig_mesh.GetNE(); e++) { - int attr = orig_mesh->GetAttribute(e); + int attr = orig_mesh.GetAttribute(e); if (!mat_marker[attr - 1]) { elem_delete[e] = true; @@ -833,13 +858,13 @@ 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 f = 0; f < orig_mesh.GetNumFaces(); f++) { const int &be = orig_bdr_faces[f]; if (be >= 0) { int e1, e2; - orig_mesh->GetFaceElements(f, &e1, &e2); + orig_mesh.GetFaceElements(f, &e1, &e2); if ((e1 < 0 || elem_delete[e1]) && (e2 < 0 || elem_delete[e2])) { // Mpi::Print("Deleting an unattached boundary element!\n"); @@ -848,15 +873,15 @@ std::map> CheckMesh(std::unique_ptr &orig_me } } } - if (new_ne < orig_mesh->GetNE()) + if (new_ne < orig_mesh.GetNE()) { Mpi::Print("Removed {:d} unmarked domain elements from the mesh\n", - orig_mesh->GetNE() - new_ne); + orig_mesh.GetNE() - new_ne); } - if (new_nbdr < orig_mesh->GetNBE()) + if (new_nbdr < orig_mesh.GetNBE()) { Mpi::Print("Removed {:d} unattached boundary elements from the mesh\n", - orig_mesh->GetNBE() - new_nbdr); + orig_mesh.GetNBE() - new_nbdr); } } int new_ne_step1 = new_ne; @@ -866,16 +891,17 @@ std::map> CheckMesh(std::unique_ptr &orig_me { // 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!"); + 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; - for (int f = 0; f < orig_mesh->GetNumFaces(); f++) + 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); + 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)) @@ -884,7 +910,7 @@ std::map> CheckMesh(std::unique_ptr &orig_me add_bdr_faces[f] = 1; add_bdr_ext++; } - else if (orig_mesh->GetAttribute(e1) != orig_mesh->GetAttribute(e2)) + 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"); @@ -910,19 +936,20 @@ std::map> CheckMesh(std::unique_ptr &orig_me if (add_subdomain) { - // Add new boundary elements at interfaces between elements beloning to different + // 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++) + 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); + + 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]) @@ -949,61 +976,54 @@ std::map> CheckMesh(std::unique_ptr &orig_me // 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()) + new_ne_step2 == orig_mesh.GetNE() && new_nbdr == new_nbdr_step1 && + new_nbdr_step1 == new_nbdr_step2 && new_nbdr_step2 == orig_mesh.GetNBE()) { return new_attr_map; } - std::unique_ptr new_mesh = - std::make_unique(orig_mesh->Dimension(), orig_mesh->GetNV(), new_ne, - new_nbdr, orig_mesh->SpaceDimension()); + + auto new_mesh = mfem::Mesh(orig_mesh.Dimension(), orig_mesh.GetNV(), new_ne, new_nbdr, + orig_mesh.SpaceDimension()); // Copy vertices and non-deleted domain and boundary elements. - for (int v = 0; v < orig_mesh->GetNV(); v++) + for (int v = 0; v < orig_mesh.GetNV(); v++) { - new_mesh->AddVertex(orig_mesh->GetVertex(v)); + new_mesh.AddVertex(orig_mesh.GetVertex(v)); } - for (int e = 0; e < orig_mesh->GetNE(); e++) + for (int e = 0; e < orig_mesh.GetNE(); e++) { if (!elem_delete[e]) { - mfem::Element *ne = orig_mesh->GetElement(e)->Duplicate(new_mesh.get()); - new_mesh->AddElement(ne); + mfem::Element *ne = orig_mesh.GetElement(e)->Duplicate(&new_mesh); + new_mesh.AddElement(ne); } } - for (int be = 0; be < orig_mesh->GetNBE(); be++) + for (int be = 0; be < orig_mesh.GetNBE(); be++) { if (!bdr_delete[be]) { - mfem::Element *ne = orig_mesh->GetBdrElement(be)->Duplicate(new_mesh.get()); - new_mesh->AddBdrElement(ne); + mfem::Element *ne = orig_mesh.GetBdrElement(be)->Duplicate(&new_mesh); + new_mesh.AddBdrElement(ne); } } // Add new boundary elements. if (add_bdr || add_subdomain) { + auto FlipVertices = [](mfem::Element *e) { mfem::Array v; e->GetVertices(v); - int start = 0, end = v.Size() - 1; - while (start < end) - { - int t = v[start]; - v[start] = v[end]; - v[end] = t; - start++; - end--; - } + std::reverse(v.begin(), v.end()); e->SetVertices(v.HostRead()); }; // 1-based, some boundary attributes may be empty since they were removed from the // original mesh, but to keep indices the same as config file we don't compact the // list. - int max_bdr_attr = orig_mesh->bdr_attributes.Max(); - for (int f = 0; f < orig_mesh->GetNumFaces(); f++) + int max_bdr_attr = orig_mesh.bdr_attributes.Max(); + for (int f = 0; f < orig_mesh.GetNumFaces(); f++) { if (add_bdr_faces[f] > 0) { @@ -1013,23 +1033,23 @@ std::map> CheckMesh(std::unique_ptr &orig_me // inverse so that the attributes of e1 and e2 can be easily referenced using the // new attribute. Since attributes are in 1-based indexing, a, b > 0. int e1, e2, a = 0, b = 0; - orig_mesh->GetFaceElements(f, &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) { - a = std::max(orig_mesh->GetAttribute(e1), orig_mesh->GetAttribute(e2)); - b = (a == orig_mesh->GetAttribute(e1)) ? orig_mesh->GetAttribute(e2) - : orig_mesh->GetAttribute(e1); + a = std::max(orig_mesh.GetAttribute(e1), orig_mesh.GetAttribute(e2)); + b = (a == orig_mesh.GetAttribute(e1)) ? orig_mesh.GetAttribute(e2) + : orig_mesh.GetAttribute(e1); } else if (!no_e1) { - a = orig_mesh->GetAttribute(e1); + a = orig_mesh.GetAttribute(e1); b = 0; } else if (!no_e2) { - a = orig_mesh->GetAttribute(e2); + a = orig_mesh.GetAttribute(e2); b = 0; } MFEM_VERIFY(a + b > 0, "Invalid new boundary element attribute!"); @@ -1040,16 +1060,16 @@ std::map> CheckMesh(std::unique_ptr &orig_me } // Add the boundary elements with the new boundary attribute. - mfem::Element *ne = orig_mesh->GetFace(f)->Duplicate(new_mesh.get()); + mfem::Element *ne = orig_mesh.GetFace(f)->Duplicate(&new_mesh); ne->SetAttribute(new_attr); - new_mesh->AddBdrElement(ne); + new_mesh.AddBdrElement(ne); if (add_bdr_faces[f] > 1) { // Flip order of vertices to reverse normal direction of second added element. - ne = orig_mesh->GetFace(f)->Duplicate(new_mesh.get()); + ne = orig_mesh.GetFace(f)->Duplicate(&new_mesh); FlipVertices(ne); ne->SetAttribute(new_attr); - new_mesh->AddBdrElement(ne); + new_mesh.AddBdrElement(ne); // Mpi::Print("Adding two BE with attr {:d} from elements {:d} and {:d}\n", // new_attr, a, b); } @@ -1061,22 +1081,22 @@ std::map> CheckMesh(std::unique_ptr &orig_me // projecting nodes onto the new mesh for the non-trimmed vdofs (accounts for new // boundary elements too since no new dofs are added). See the MFEM trimmer miniapp for // reference. - new_mesh->FinalizeTopology(); - new_mesh->Finalize(); - new_mesh->RemoveUnusedVertices(); - if (orig_mesh->GetNodes()) + new_mesh.FinalizeTopology(); + new_mesh.Finalize(); + new_mesh.RemoveUnusedVertices(); + if (orig_mesh.GetNodes()) { - const mfem::GridFunction *nodes = orig_mesh->GetNodes(); + 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; - new_mesh->SetCurvature(order, discont, sdim, ordering); - mfem::GridFunction *new_nodes = new_mesh->GetNodes(); + 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 @@ -1084,7 +1104,7 @@ std::map> CheckMesh(std::unique_ptr &orig_me mfem::Array vdofs, new_vdofs; mfem::Vector loc_vec; int te = 0; - for (int e = 0; e < orig_mesh->GetNE(); e++) + for (int e = 0; e < orig_mesh.GetNE(); e++) { if (!elem_delete[e]) { @@ -1096,32 +1116,27 @@ std::map> CheckMesh(std::unique_ptr &orig_me } } } + // Reperform Finalization, as the true argument allows for element vertex reordering, + // to ensure that the longest edge is first. This prepares for conformal adaptation. + new_mesh.Finalize(true); orig_mesh = std::move(new_mesh); return new_attr_map; } std::unique_ptr DistributeMesh(MPI_Comm comm, std::unique_ptr &smesh, - std::unique_ptr &partitioning) + const std::unique_ptr &partitioning) { // Take a serial mesh and partitioning on the root process and construct the global // parallel mesh. For now, prefer the MPI-based version. -#if 0 + if constexpr (false) { // Write each processor's component to file. - std::string tmp = iodata.problem.output; - if (tmp.back() != '/') - { - tmp += '/'; - } - tmp += "tmp/"; + std::string tmp = (std::filesystem::temp_directory_path() += "/"); + 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), partitioning.get()); for (int i = 0; i < Mpi::Size(comm); i++) { @@ -1160,7 +1175,7 @@ std::unique_ptr DistributeMesh(MPI_Comm comm, } return pmesh; } -#endif + else { // Send each processor's component as a byte string. std::vector so; diff --git a/scripts/schema/config/model.json b/scripts/schema/config/model.json index 00e693fba..beba3d363 100644 --- a/scripts/schema/config/model.json +++ b/scripts/schema/config/model.json @@ -19,6 +19,25 @@ "properties": { "UniformLevels": { "type": "integer", "minimum": 0 }, + "Adaptation": + { + "type": "object", + "additionalProperties": false, + "required": ["MaxIts"], + "properties": + { + "Tol": {"type": "number", "exclusiveMinimum": 0.0}, + "MaxIts": {"type": "integer", "inclusiveMinimum": 0}, + "MinIts": {"type": "integer", "inclusiveMinimum": 0}, + "UpdateFraction": {"type": "number", "exclusiveMinimum": 0.0, "exclusiveMaximum": 1.0}, + "CoarseningFraction": {"type": "number", "exclusiveMaximum": 1.0, "inclusiveMinimum": 0.0}, + "DOFLimit": {"type": "number", "inclusiveMinimum": 0}, + "SaveStep": {"type": "integer", "inclusiveMinimum": 0}, + "MaxNCLevels": {"type": "integer", "inclusiveMinimum": 0}, + "Nonconformal": {"type": "boolean"}, + "MaximumImbalance": {"type": "number", "inclusiveMinimum": 1.0} + } + }, "Boxes": { "type": "array",