From 61b9b89fbc701499d7e8d9a646b6ab393240de54 Mon Sep 17 00:00:00 2001 From: Leila Ghaffari Date: Tue, 13 Aug 2024 11:50:38 -0700 Subject: [PATCH] Fix format --- palace/fem/fespace.cpp | 11 ++++- palace/fem/libceed/restriction.cpp | 32 +++++++-------- palace/utils/configfile.cpp | 33 ++++++++++----- palace/utils/geodata.cpp | 66 ++++++++++++++++++++++-------- 4 files changed, 96 insertions(+), 46 deletions(-) diff --git a/palace/fem/fespace.cpp b/palace/fem/fespace.cpp index edb56e09d..d6153fa6d 100644 --- a/palace/fem/fespace.cpp +++ b/palace/fem/fespace.cpp @@ -138,7 +138,11 @@ CeedBasis FiniteElementSpace::BuildCeedBasis(const mfem::FiniteElementSpace &fes mfem::IsoparametricTransformation T; const mfem::FiniteElement *fe_nodal; fe_nodal = fespace.GetMesh()->GetNodalFESpace()->FEColl()->FiniteElementForGeometry(geom); - if (!fe_nodal) fe_nodal = fespace.GetMesh()->GetNodalFESpace()->FEColl()->TraceFiniteElementForGeometry(geom); + if (!fe_nodal) + { + fe_nodal = + fespace.GetMesh()->GetNodalFESpace()->FEColl()->TraceFiniteElementForGeometry(geom); + } T.SetFE(fe_nodal); const int q_order = fem::DefaultIntegrationOrder::Get(T); const mfem::IntegrationRule &ir = mfem::IntRules.Get(geom, q_order); @@ -147,7 +151,10 @@ CeedBasis FiniteElementSpace::BuildCeedBasis(const mfem::FiniteElementSpace &fes CeedBasis val; const mfem::FiniteElement *fe; fe = fespace.FEColl()->FiniteElementForGeometry(geom); - if (!fe) fe = fespace.FEColl()->TraceFiniteElementForGeometry(geom); + if (!fe) + { + fe = fespace.FEColl()->TraceFiniteElementForGeometry(geom); + } const int vdim = fespace.GetVDim(); ceed::InitBasis(*fe, ir, vdim, ceed, &val); return val; diff --git a/palace/fem/libceed/restriction.cpp b/palace/fem/libceed/restriction.cpp index 6eaee312d..f63450aeb 100644 --- a/palace/fem/libceed/restriction.cpp +++ b/palace/fem/libceed/restriction.cpp @@ -167,16 +167,15 @@ void InitNativeRestr(const mfem::FiniteElementSpace &fespace, fespace.GetMesh()->GetBdrElementAdjacentElement(e, elem_id, face_info); mfem::Geometry::Type face_geom = fespace.GetMesh()->GetBdrElementGeometry(e); face_info = fespace.GetMesh()->EncodeFaceInfo( - fespace.GetMesh()->DecodeFaceInfoLocalIndex(face_info), - mfem::Geometry::GetInverseOrientation( - face_geom, fespace.GetMesh()->DecodeFaceInfoOrientation(face_info)) - ); + fespace.GetMesh()->DecodeFaceInfoLocalIndex(face_info), + mfem::Geometry::GetInverseOrientation( + face_geom, fespace.GetMesh()->DecodeFaceInfoOrientation(face_info))); mfem::IntegrationPointTransformation Loc1; - fespace.GetMesh()->GetLocalFaceTransformation(fespace.GetMesh()->GetBdrElementType(e), - fespace.GetMesh()->GetElementType(elem_id), - Loc1.Transf, face_info); + fespace.GetMesh()->GetLocalFaceTransformation( + fespace.GetMesh()->GetBdrElementType(e), + fespace.GetMesh()->GetElementType(elem_id), Loc1.Transf, face_info); const mfem::FiniteElement *face_el = fespace.GetTraceElement(elem_id, face_geom); - MFEM_VERIFY(dynamic_cast(face_el), + MFEM_VERIFY(dynamic_cast(face_el), "Mesh requires nodal Finite Element."); mfem::IntegrationRule face_ir(face_el->GetDof()); Loc1.Transf.ElementNo = elem_id; @@ -188,7 +187,7 @@ void InitNativeRestr(const mfem::FiniteElementSpace &fespace, face_tr.ElementType = mfem::ElementTransformation::BDR_ELEMENT; face_tr.mesh = fespace.GetMesh(); face_tr.Attribute = fespace.GetMesh()->GetBdrAttribute(e); - mfem::DenseMatrix &face_pm = face_tr.GetPointMat(); // dim x dof + mfem::DenseMatrix &face_pm = face_tr.GetPointMat(); // dim x dof face_tr.Reset(); fespace.GetMesh()->GetNodes()->GetVectorValues(Loc1.Transf, face_ir, face_pm); @@ -196,15 +195,18 @@ void InitNativeRestr(const mfem::FiniteElementSpace &fespace, mfem::DenseMatrix elem_pm; const mfem::FiniteElement *fe_elem = fespace.GetFE(elem_id); const mfem::IntegrationRule &fe_elem_nodes = fe_elem->GetNodes(); - mfem::ElementTransformation *T = fespace.GetMesh()->GetElementTransformation(elem_id); + mfem::ElementTransformation *T = + fespace.GetMesh()->GetElementTransformation(elem_id); T->Transform(fe_elem_nodes, elem_pm); // Find the dofs + // TODO: Instead of comparing the point matrices for each element, we might + // be able to configure the face info once using the reference element. mfem::real_t tol = 1E-5; mfem::Array elem_dofs; fespace.GetElementDofs(elem_id, elem_dofs, dof_trans); dofs.SetSize(P); - for (int l=0; l< P; l++) + for (int l = 0; l < P; l++) { double norm_f = 0.0; for (int m = 0; m < face_pm.Height(); m++) @@ -212,7 +214,7 @@ void InitNativeRestr(const mfem::FiniteElementSpace &fespace, norm_f += face_pm(m, l) * face_pm(m, l); } norm_f = std::sqrt(norm_f); - for (int m=0; m< elem_pm.Width(); m++) + for (int m = 0; m < elem_pm.Width(); m++) { double norm_e = 0.0; for (int n = 0; n < elem_pm.Height(); n++) @@ -230,11 +232,10 @@ void InitNativeRestr(const mfem::FiniteElementSpace &fespace, else { double relative_tol = tol * std::max(norm_f, norm_e); - bool same_coord = false; double diff = 0.0; - for (int o = 0; o < elem_pm.Height(); o++) + for (int n = 0; n < elem_pm.Height(); n++) { - diff += std::fabs(elem_pm(o, m) - face_pm(o, l)); + diff += std::fabs(elem_pm(n, m) - face_pm(n, l)); } if (diff <= relative_tol) { @@ -319,7 +320,6 @@ void InitNativeRestr(const mfem::FiniteElementSpace &fespace, } } } - } use_el_orients += use_el_orients_loc; } diff --git a/palace/utils/configfile.cpp b/palace/utils/configfile.cpp index 226c7558c..64f69c913 100644 --- a/palace/utils/configfile.cpp +++ b/palace/utils/configfile.cpp @@ -1033,9 +1033,8 @@ void PeriodicBoundaryData::SetUp(json &boundaries) { return; } - MFEM_VERIFY( - periodic->is_array(), - "\"Periodic\" should specify an array in the configuration file!"); + MFEM_VERIFY(periodic->is_array(), + "\"Periodic\" should specify an array in the configuration file!"); for (auto it = periodic->begin(); it != periodic->end(); ++it) { MFEM_VERIFY(it->find("DonorAttributes") != it->end(), @@ -1046,9 +1045,10 @@ void PeriodicBoundaryData::SetUp(json &boundaries) "configuration file!"); PeriodicData &data = vecdata.emplace_back(); data.donor_attributes = it->at("DonorAttributes").get>(); // Required - data.receiver_attributes = it->at("ReceiverAttributes").get>(); // Required - auto direction = it->at("Direction").get(); // Required (TODO: the user can provide the vectors) - auto distance = it->at("Distance").get(); // Required + data.receiver_attributes = + it->at("ReceiverAttributes").get>(); // Required + auto direction = it->at("Direction").get(); // Required + auto distance = it->at("Distance").get(); // Required for (auto &c : direction) { c = std::tolower(c); @@ -1059,9 +1059,18 @@ void PeriodicBoundaryData::SetUp(json &boundaries) const bool xfound = xpos != std::string::npos; const bool yfound = ypos != std::string::npos; const bool zfound = zpos != std::string::npos; - if (xfound) data.translation[0] = distance; - if (yfound) data.translation[1] = distance; - if (zfound) data.translation[2] = distance; + if (xfound) + { + data.translation[0] = distance; + } + if (yfound) + { + data.translation[1] = distance; + } + if (zfound) + { + data.translation[2] = distance; + } // Cleanup it->erase("DonorAttributes"); @@ -1429,8 +1438,10 @@ void BoundaryData::SetUp(json &config) } for (const auto &data : periodic) { - donor_attributes.insert(donor_attributes.end(), data.donor_attributes.begin(), data.donor_attributes.end()); - receiver_attributes.insert(receiver_attributes.end(), data.receiver_attributes.begin(), data.receiver_attributes.end()); + donor_attributes.insert(donor_attributes.end(), data.donor_attributes.begin(), + data.donor_attributes.end()); + receiver_attributes.insert(receiver_attributes.end(), data.receiver_attributes.begin(), + data.receiver_attributes.end()); translation.insert(translation.end(), data.translation.begin(), data.translation.end()); } for (const auto &[idx, data] : waveport) diff --git a/palace/utils/geodata.cpp b/palace/utils/geodata.cpp index e0e9d0ad5..64b6a12a9 100644 --- a/palace/utils/geodata.cpp +++ b/palace/utils/geodata.cpp @@ -38,7 +38,8 @@ namespace constexpr auto MSH_FLT_PRECISION = std::numeric_limits::max_digits10; // Load the serial mesh from disk. -std::unique_ptr LoadMesh(const std::string &, bool, const config::BoundaryData &, double); +std::unique_ptr LoadMesh(const std::string &, bool, + const config::BoundaryData &, double); // 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 @@ -96,7 +97,8 @@ std::unique_ptr ReadMesh(MPI_Comm comm, const IoData &iodata) if ((use_mesh_partitioner && Mpi::Root(comm)) || (!use_mesh_partitioner && Mpi::Root(node_comm))) { - smesh = LoadMesh(iodata.model.mesh, iodata.model.remove_curvature, iodata.boundaries, iodata.model.L0); + smesh = LoadMesh(iodata.model.mesh, iodata.model.remove_curvature, iodata.boundaries, + iodata.model.L0); MFEM_VERIFY(!(smesh->Nonconforming() && use_mesh_partitioner), "Cannot use mesh partitioner on a nonconforming mesh!"); } @@ -228,7 +230,8 @@ std::unique_ptr ReadMesh(MPI_Comm comm, const IoData &iodata) } int width = 1 + static_cast(std::log10(Mpi::Size(comm) - 1)); std::unique_ptr gsmesh = - LoadMesh(iodata.model.mesh, iodata.model.remove_curvature, iodata.boundaries, iodata.model.L0); + LoadMesh(iodata.model.mesh, iodata.model.remove_curvature, iodata.boundaries, + iodata.model.L0); std::unique_ptr gpartitioning = GetMeshPartitioning(*gsmesh, Mpi::Size(comm)); mfem::ParMesh gpmesh(comm, *gsmesh, gpartitioning.get(), 0); { @@ -1664,7 +1667,8 @@ template void AttrToMarker(int, const std::vector &, mfem::Array &, bo namespace { -std::unique_ptr LoadMesh(const std::string &mesh_file, bool remove_curvature, const config::BoundaryData &boundaries, double L0) +std::unique_ptr LoadMesh(const std::string &mesh_file, bool remove_curvature, + const config::BoundaryData &boundaries, double L0) { // 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 @@ -1721,18 +1725,27 @@ std::unique_ptr LoadMesh(const std::string &mesh_file, bool remove_c { mesh->EnsureNodes(); } - if (!boundaries.periodic.empty()) { + if (!boundaries.periodic.empty()) + { mfem::real_t tol = 1E-5 / L0; auto periodic_mesh = std::move(mesh); for (auto &data : boundaries.periodic) { std::vector translation; mfem::Vector translation_vec(data.translation.size()); - std::copy(data.translation.begin(), data.translation.end(), translation_vec.GetData()); - for (int i = 0; i < translation_vec.Size(); ++i) translation_vec[i] /= L0; + std::copy(data.translation.begin(), data.translation.end(), + translation_vec.GetData()); + for (int i = 0; i < translation_vec.Size(); ++i) + { + translation_vec[i] /= L0; + } translation.push_back(translation_vec); - auto p_mesh = std::make_unique(mfem::Mesh::MakePeriodic(*periodic_mesh, periodic_mesh->CreatePeriodicVertexMapping(translation, tol))); - if (p_mesh) periodic_mesh = std::move(p_mesh); + auto p_mesh = std::make_unique(mfem::Mesh::MakePeriodic( + *periodic_mesh, periodic_mesh->CreatePeriodicVertexMapping(translation, tol))); + if (p_mesh) + { + periodic_mesh = std::move(p_mesh); + } } mesh = std::move(periodic_mesh); } @@ -2172,7 +2185,8 @@ std::map> CheckMesh(std::unique_ptr &orig_me int new_ne = orig_mesh->GetNE(); int new_nbe = orig_mesh->GetNBE(); std::vector elem_delete(orig_mesh->GetNE(), false), - bdr_elem_delete(orig_mesh->GetNBE(), false), bdr_elem_periodic(orig_mesh->GetNBE(), false); + bdr_elem_delete(orig_mesh->GetNBE(), false), + bdr_elem_periodic(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++) @@ -2181,14 +2195,27 @@ std::map> CheckMesh(std::unique_ptr &orig_me bool periodic_be = false; for (auto &data : iodata.boundaries.periodic) { - for (auto d_attr : data.donor_attributes) if (d_attr == attr) periodic_be = true; - for (auto r_attr : data.receiver_attributes) if (r_attr == attr) periodic_be = true; + for (auto d_attr : data.donor_attributes) + { + if (d_attr == attr) + { + periodic_be = true; + } + } + for (auto r_attr : data.receiver_attributes) + { + if (r_attr == attr) + { + periodic_be = true; + } + } } int face, o, e1, e2; orig_mesh->GetBdrElementFace(be, &face, &o); orig_mesh->GetFaceElements(face, &e1, &e2); - // If there are two elements associated with the face on the boundary, then the mesh is periodic. + // If there are two elements associated with the face on the boundary, then the mesh is + // periodic. if ((e1 >= 0) && (e2 >= 0)) { auto Normal = [&](int e) @@ -2205,12 +2232,16 @@ std::map> CheckMesh(std::unique_ptr &orig_me }; auto normal1 = Normal(e1); auto normal2 = Normal(e2); - if ((normal1 * normal1 == 0) && (normal2 * normal2 == 0)) periodic_be = true; - if (normal1 * normal2 > 0) periodic_be = true; - + if ((normal1 * normal1 == 0) && (normal2 * normal2 == 0)) + { + periodic_be = true; + } + if (normal1 * normal2 > 0) + { + periodic_be = true; + } int f, o; orig_mesh->GetBdrElementFace(be, &f, &o); - face_to_be[f] = be; if (periodic_be) { bdr_elem_periodic[be] = true; @@ -2220,6 +2251,7 @@ std::map> CheckMesh(std::unique_ptr &orig_me MFEM_VERIFY(face_to_be.find(f) == face_to_be.end(), "Mesh should not define boundary elements multiple times!"); } + face_to_be[f] = be; } }