diff --git a/src/serac/numerics/functional/domain.cpp b/src/serac/numerics/functional/domain.cpp index dca0342319..33f77ba516 100644 --- a/src/serac/numerics/functional/domain.cpp +++ b/src/serac/numerics/functional/domain.cpp @@ -535,6 +535,49 @@ void Domain::insert_restriction(const serac::fes_t* fes, FunctionSpace space) const BlockElementRestriction& Domain::get_restriction(FunctionSpace space) { return restriction_operators.at(space); }; +void Domain::insert_shared_interior_face_list() +{ + // Weights only need to be computed for Domain of InteriorFaces type + SLIC_ERROR_ROOT_IF(type_ != Domain::Type::InteriorFaces, "This method is only for interior face domains"); + + // make a list if we don't already have one + if (shared_interior_face_ids_.empty()) { + if (dim_ == 1) { + int i = 0; + for (int f : mfem_edge_ids_) { + mfem::Mesh::FaceInformation info = mesh_.GetFaceInformation(f); + + if (info.IsShared()) { + shared_interior_face_ids_.push_back(i); + } + + ++i; + } + } else if (dim_ == 2) { + int i = 0; + for (int f : mfem_tri_ids_) { + mfem::Mesh::FaceInformation info = mesh_.GetFaceInformation(f); + + if (info.IsShared()) { + shared_interior_face_ids_.push_back(i); + } + + ++i; + } + + for (int f : mfem_quad_ids_) { + mfem::Mesh::FaceInformation info = mesh_.GetFaceInformation(f); + + if (info.IsShared()) { + shared_interior_face_ids_.push_back(i); + } + + ++i; + } + } + } +} + /////////////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////////// /////////////////////////////////////////////////////////////////////////////////////// @@ -604,6 +647,8 @@ Domain InteriorFaces(const mesh_t& mesh) } } + output.insert_shared_interior_face_list(); + return output; } diff --git a/src/serac/numerics/functional/domain.hpp b/src/serac/numerics/functional/domain.hpp index cf5171b7ba..fb29cfec9e 100644 --- a/src/serac/numerics/functional/domain.hpp +++ b/src/serac/numerics/functional/domain.hpp @@ -96,6 +96,9 @@ struct Domain { */ std::map restriction_operators; + /// @brief Ids of interior faces that lie on the boundary shared by two processors + std::vector shared_interior_face_ids_; + /** * @brief empty Domain constructor, with connectivity info to be populated later */ @@ -244,6 +247,12 @@ struct Domain { /// thereby populate the element lists. void addElements(const std::vector& geom_id, const std::vector& elem_id, mfem::Geometry::Type element_geometry); + + /** + * @brief Find the list of interior faces shared by two processors to make sure + * these faces are only integrated once. + */ + void insert_shared_interior_face_list(); }; /// @brief constructs a domain from all the elements in a mesh diff --git a/src/serac/numerics/functional/element_restriction.cpp b/src/serac/numerics/functional/element_restriction.cpp index 998a4d9e9b..8a5741933f 100644 --- a/src/serac/numerics/functional/element_restriction.cpp +++ b/src/serac/numerics/functional/element_restriction.cpp @@ -4,9 +4,6 @@ #include "serac/numerics/functional/geometry.hpp" -// TODO REMOVE AFTER DEBUGGING -#include "serac/infrastructure/mpi_fstream.hpp" - std::vector > lexicographic_permutations(int p) { // p == 0 is admissible for L2 spaces, but lexicographic permutations @@ -474,15 +471,6 @@ axom::Array GetFaceDofs(const serac::f std::vector > local_face_dofs = geom_local_face_dofs(p); std::vector > lex_perm = lexicographic_permutations(p); - // sam: this function contains several comments that relate to an investigation into adopting - // mfem's FaceNbrData pattern (which ended up being too invasive to implement initially). - // I leave some of the commented code here so Julian's recommendations are not lost. - // - // int components_per_node = fes->GetVDim(); - // bool by_vdim = fes->GetOrdering() == mfem::Ordering::byVDIM; - // fes->ExchangeFaceNbrData(); - // int LSize = fes->GetProlongationMatrix()->Height(); - uint64_t n = 0; for (int f : mfem_face_ids) { @@ -499,24 +487,6 @@ axom::Array GetFaceDofs(const serac::f mfem::Array elem_ids; face_to_elem->GetRow(f, elem_ids); -#if 0 - mpi::out << f << " elem ids: "; - for (auto elem : elem_ids) { - mpi::out << elem << " "; - } - mpi::out << std::endl; - - mfem::Array test_dofs; - fes->GetFaceDofs(f, test_dofs); - - mpi::out << " dofs (sz = " << test_dofs.Size() << "): "; - for (int z = 0; z < test_dofs.Size(); z++) { - mpi::out << test_dofs[z] << " "; - } - mpi::out << std::endl; - mpi::out << std::endl; -#endif - for (auto elem : elem_ids) { // 2a. get the list of faces (and their orientations) that belong to that element ... mfem::Array elem_side_ids, orientations; @@ -537,6 +507,7 @@ axom::Array GetFaceDofs(const serac::f } // 2b. ... and find `i` such that `elem_side_ids[i] == f` + // This generates a local_face_index that coincides with ones from GetFaceInformation int i; for (i = 0; i < elem_side_ids.Size(); i++) { if (elem_side_ids[i] == f) break; @@ -561,34 +532,65 @@ axom::Array GetFaceDofs(const serac::f face_dofs.push_back(uint64_t(elem_dof_ids[local_face_dofs[uint32_t(elem_geom)](i, k)])); } -#if 0 - // (5. optional) add remaining dofs that were omitted on shared faces + // 5. add remaining dofs that were omitted on shared faces + // FaceNbrData includes all dofs in the "volumetric" element adjacent to the shared faces + // We need to again find our face on the ghost element and extract only the face dofs + // This may not yet work for NC faces if (info.IsShared()) { + // FaceNbrData should be already exchanged in Functional when invoked by ParGridFunction + // which also invokes ExchangeFaceNbrData on subsequent ParFiniteElementSpace and ParMesh + int components_per_node = fes->GetVDim(); + bool by_vdim = fes->GetOrdering() == mfem::Ordering::byVDIM; + int LSize = fes->GetProlongationMatrix()->Height(); + + // Additional safeguard to make sure ByNODES ordering in vector L2 field raise an error + SLIC_ERROR_IF(!by_vdim && components_per_node > 1, "Unsupported: L2 vector field ordered by nodes"); mfem::Array shared_elem_vdof_ids; + mfem::Array shared_elem_dof_ids; // on my processor - // VVVVVVV + // VVVVVVV // dofs for the face [0 1 3 4 | 5 8 12 13] // ^^^^^^^^^ // on the other processor int other_element_id = info.element[1].index; - fes->GetFaceNbrElementVDofs(other_element_id, shared_elem_vdof_ids); // indices into vector from FaceNbrData + fes->GetFaceNbrElementVDofs(other_element_id, shared_elem_vdof_ids); // indices into vector from FaceNbrData - mpi::out << " this face is also shared so it has additional dofs: "; - int dofs_per_face = shared_elem_vdof_ids.Size() / components_per_node; - int stride = (by_vdim) ? components_per_node : 1; - - // sam: is this right? // byVDIM == x y z x y z x y z x y z // byNODES == x x x x y y y y z z z z - for (int k = 0; k < dofs_per_face; k++) { - face_dofs.push_back(uint64_t(fes->VDofToDof(shared_elem_vdof_ids[k * stride]) + LSize)); - mpi::out << face_dofs.back().index() << " "; + int dofs_per_element = shared_elem_vdof_ids.Size() / components_per_node; + int stride = (by_vdim) ? components_per_node : 1; + + shared_elem_dof_ids.SetSize(dofs_per_element); + for (int k = 0; k < dofs_per_element; ++k) { + shared_elem_dof_ids[k] = fes->VDofToDof(shared_elem_vdof_ids[k * stride]); // Change vdof to dof indices + } + + // Find the dofs on the shared face on the ghost element side + // Neighbor element local_face_index and orientation is already available from GetFaceInformation + mfem::Geometry::Type ghost_geom = fes->GetParMesh()->face_nbr_elements[other_element_id]->GetGeometryType(); + int j = info.element[1].local_face_id; + int ghost_orientation; + if (mesh->Dimension() == 2) { + // In 2D, orientations[i] sometimes is inverted as compared to the ones from GetFaceInformation + // In this case, we stay with the orientations constructed previously by GetElementEdges, which is strictly + // CCW So we set ghost orientation to be exactly the opposite of the original orientation (orientations[i]) + // side note: even if you use info.element[1].orientation, looks like it's still fine. + // The only difference is the order those face dof indices gets registered in std::vector face_dofs + ghost_orientation = (orientation == 0) ? 1 : 0; + } else { + // In 3D, I think it's consistently CCW. So we directly use the orientation from GetFaceInformation + // This orientation is already consistent with the permutation in {0, 1, ... , n} form + ghost_orientation = info.element[1].orientation; + } + + // extract only the dofs that correspond to side `j` + for (auto k : face_perm(ghost_orientation)) { + face_dofs.push_back(uint64_t(shared_elem_dof_ids[local_face_dofs[uint32_t(ghost_geom)](j, k)] + + LSize / components_per_node)); } - mpi::out << std::endl; } -#endif } // H1 and Hcurl spaces are more straight-forward, since diff --git a/src/serac/numerics/functional/functional.hpp b/src/serac/numerics/functional/functional.hpp index f6d03cb5ae..c9fdd7e8ef 100644 --- a/src/serac/numerics/functional/functional.hpp +++ b/src/serac/numerics/functional/functional.hpp @@ -106,25 +106,6 @@ inline void check_for_unsupported_elements(const mfem::Mesh& mesh) } } -/** - * @brief function for verifying that DG spaces aren't used on interior face integrals over meshes that contain "shared" - * faces - * - * sam: I would like to support these "shared" faces, but apparently mfem handles them in a fundamentally different way - * than the "finite element operator decomposition" pattern used by everything else (see: "ExchangeFaceNbrData") - */ -inline void check_interior_face_compatibility(const mfem::Mesh& mesh, const FunctionSpace space) -{ - if (space.family == Family::L2) { - const mfem::ParMesh* pmesh = dynamic_cast(&mesh); - if (pmesh) { - SLIC_ERROR_IF( - pmesh->GetNSharedFaces() > 0, - "interior face integrals involving DG function spaces don't currently support meshes with shared faces"); - } - } -} - /** * @brief create an mfem::ParFiniteElementSpace from one of serac's * tag types: H1, Hcurl, L2 @@ -166,6 +147,89 @@ generateParFiniteElementSpace(mfem::ParMesh* mesh) return std::pair(std::move(fes), std::move(fec)); } +/** + * @brief helper function to locally cast away const on FE space so we can update face neighbor + * data with ExchangeFaceNbrData. This is ok because : 1) the original trial FE space is declared + * without const; 2) we constrained the non-constness locally; 3) the locally owned data associated + * with the trial function space is NOT altered and ONLY ghost data is updated. + */ +inline void updateFaceNbrData(const mfem::ParFiniteElementSpace* const_trial_space, mfem::ParGridFunction& trial_pgf, + mfem::Vector& trial_tdof_vals) +{ + mfem::ParFiniteElementSpace* nonconst_trial_space = const_cast(const_trial_space); + + // Link the ParGridFunction to external FE space and data by reference + trial_pgf.MakeRef(nonconst_trial_space, trial_tdof_vals.GetData()); + + // Exchange face nbr data to set the vector in ParGridFunction with the right size + // This should automatically invoke ExchangeFaceNbrData on the subsequent ParFESpace and ParMesh + trial_pgf.ExchangeFaceNbrData(); +} + +/** + * @brief helper functional to reorder the ordering of FaceNbrData for L2 space to byVDIM and append this vector + * to the end of local dof vector, which will result in a vector in form [ --- L --- | --- FND --- ] + */ +inline void appendFaceNbrData(const mfem::ParFiniteElementSpace* trial_space, const mfem::ParGridFunction& trial_pgf, + const int LSize, mfem::Vector& input_L) +{ + int num_ghost_elem = trial_space->GetParMesh()->GetNFaceNeighborElements(); + int components_per_node = trial_space->GetVDim(); + const mfem::Vector& face_neighbor_data = trial_pgf.FaceNbrData(); + + int offset = 0; + for (int n = 0; n < num_ghost_elem; ++n) { + mfem::Array old_shared_elem_vdof_ids; + trial_space->GetFaceNbrElementVDofs(n, old_shared_elem_vdof_ids); + int dofs_per_element = old_shared_elem_vdof_ids.Size() / components_per_node; + + // Put the entries into their new index + for (int m = 0; m < old_shared_elem_vdof_ids.Size(); ++m) { + // -------------- dofs before -------------- ---- component ----- + int new_id = (m % dofs_per_element) * components_per_node + m / dofs_per_element + offset; + input_L[LSize + new_id] = face_neighbor_data(old_shared_elem_vdof_ids[m]); + } + // Increase offset for next ghost element + offset += old_shared_elem_vdof_ids.Size(); + } +} + +/** + * @brief helper functional to reorder the face_nbr_glob_dof_map for L2 space to byVDIM + */ +inline void rearrangeFaceNbrDofGlobalIndex(const mfem::ParFiniteElementSpace* trial_space, + mfem::Array& face_nbr_glob_vdof_map) +{ + int components_per_node = trial_space->GetVDim(); + if (components_per_node == 1) { + face_nbr_glob_vdof_map = trial_space->face_nbr_glob_dof_map; + return; + } + + if (components_per_node > 1 && trial_space->GetOrdering() == mfem::Ordering::byNODES) { + SLIC_ERROR_ROOT("Unsupported: L2 vector field ordered by nodes"); + } + + face_nbr_glob_vdof_map.SetSize(trial_space->face_nbr_glob_dof_map.Size()); + int num_ghost_elem = trial_space->GetParMesh()->GetNFaceNeighborElements(); + const HYPRE_BigInt* face_nbr_glob_dof_index = trial_space->face_nbr_glob_dof_map; + + int offset = 0; + for (int n = 0; n < num_ghost_elem; ++n) { + mfem::Array old_shared_elem_vdof_ids; + trial_space->GetFaceNbrElementVDofs(n, old_shared_elem_vdof_ids); + int dofs_per_element = old_shared_elem_vdof_ids.Size() / components_per_node; + + // Put the entries into their new index + for (int m = 0; m < old_shared_elem_vdof_ids.Size(); ++m) { + int new_id = (m % dofs_per_element) * components_per_node + m / dofs_per_element + offset; + face_nbr_glob_vdof_map[new_id] = face_nbr_glob_dof_index[old_shared_elem_vdof_ids[m]]; + } + // Increase offset for next ghost element + offset += old_shared_elem_vdof_ids.Size(); + } +} + /// @cond template class Functional; @@ -247,23 +311,40 @@ class Functional { { SERAC_MARK_FUNCTION; - for (uint32_t i = 0; i < num_trial_spaces; i++) { - P_trial_[i] = trial_space_[i]->GetProlongationMatrix(); - - input_L_[i].SetSize(P_trial_[i]->Height(), mfem::Device::GetMemoryType()); - - // create the necessary number of empty mfem::Vectors, to be resized later - input_E_.push_back({}); - input_E_buffer_.push_back({}); - } - test_function_space_ = {test::family, test::order, test::components}; std::array trial_families = {trials::family...}; std::array trial_orders = {trials::order...}; std::array trial_components = {trials::components...}; + for (uint32_t i = 0; i < num_trial_spaces; i++) { trial_function_spaces_[i] = {trial_families[i], trial_orders[i], trial_components[i]}; + + P_trial_[i] = trial_space_[i]->GetProlongationMatrix(); + + // For L2 spaces, configure a ParGridFunction with external data to exchange face neighbor data. + // This is for DG method where interior faces on the processor boundary need to access data on + // the neighrbor element + if (trial_function_spaces_[i].family == Family::L2) { + // Set the vector size to store tdof data + input_ldof_values_[i].SetSize(P_trial_[i]->Height(), mem_type); + + // Initialize face neighbor data vector + updateFaceNbrData(trial_space_[i], trial_pargrid_functions_[i], input_ldof_values_[i]); + + // Set the local input vector size to store local + ghost data + input_L_[i].SetSize(input_ldof_values_[i].Size() + trial_pargrid_functions_[i].FaceNbrData().Size(), mem_type); + + // Rearrange the global index of face neighbor vdofs to the correct ordering. + // This only needs to be done once + rearrangeFaceNbrDofGlobalIndex(trial_space_[i], face_nbr_glob_vdof_maps_[i]); + } else { + input_L_[i].SetSize(P_trial_[i]->Height(), mem_type); + } + + // create the necessary number of empty mfem::Vectors, to be resized later + input_E_.push_back({}); + input_E_buffer_.push_back({}); } // for (auto type : {Domain::Type::Elements, Domain::Type::BoundaryElements, Domain::Type::InteriorFaces}) { @@ -272,7 +353,16 @@ class Functional { P_test_ = test_space_->GetProlongationMatrix(); - output_L_.SetSize(P_test_->Height(), mem_type); + if (test_function_space_.family == Family::L2) { + // We only need these variables to set the output_L_ vector to the right size + mfem::ParGridFunction X_pgf; + mfem::Vector X_ldof_values(P_test_->Height(), mem_type); + + updateFaceNbrData(test_space_, X_pgf, X_ldof_values); + output_L_.SetSize(X_ldof_values.Size() + X_pgf.FaceNbrData().Size(), mem_type); + } else { + output_L_.SetSize(P_test_->Height(), mem_type); + } output_T_.SetSize(test_fes->GetTrueVSize(), mem_type); @@ -357,10 +447,8 @@ class Functional { std::vector arg_vec = {args...}; for (uint32_t i : arg_vec) { domain.insert_restriction(trial_space_[i], trial_function_spaces_[i]); - check_interior_face_compatibility(domain.mesh_, trial_function_spaces_[i]); } domain.insert_restriction(test_space_, test_function_space_); - check_interior_face_compatibility(domain.mesh_, test_function_space_); using signature = test(decltype(serac::type(trial_spaces))...); integrals_.push_back( @@ -419,7 +507,29 @@ class Functional { */ void ActionOfGradient(const mfem::Vector& input_T, mfem::Vector& output_T, uint32_t which) const { - P_trial_[which]->Mult(input_T, input_L_[which]); + // Please refer to 'operator()' below for detailed comments of this implementation + if (trial_function_spaces_[which].family == Family::L2) { + // copy input_L[which] and facenbrdata[which] into a common array like: + // first part second part + // [ --- L --- | --- FND --- ] + P_trial_[which]->Mult(input_T, input_ldof_values_[which]); + input_L_[which].SetVector(input_ldof_values_[which], 0); + updateFaceNbrData(trial_space_[which], trial_pargrid_functions_[which], input_ldof_values_[which]); + + // Ensure ghost data ordering is consistent with local data + if (trial_space_[which]->GetVDim() == 1) { + input_L_[which].SetVector(trial_pargrid_functions_[which].FaceNbrData(), input_ldof_values_[which].Size()); + } else { + if (trial_space_[which]->GetOrdering() == mfem::Ordering::byVDIM) { + int local_size = input_ldof_values_[which].Size(); + appendFaceNbrData(trial_space_[which], trial_pargrid_functions_[which], local_size, input_L_[which]); + } else { + SLIC_ERROR_ROOT("Unsupported: L2 vector field ordered by nodes"); + } + } + } else { + P_trial_[which]->Mult(input_T, input_L_[which]); + } output_L_ = 0.0; @@ -444,7 +554,17 @@ class Functional { } // scatter-add to compute global residuals - P_test_->MultTranspose(output_L_, output_T); + if (test_function_space_.family == Family::L2) { + // Extract a subvector from output_L_ for local residuals excluding ghost dofs for L2 spaces + // The output vector is in form [ --- L --- | --- FND --- ] + mfem::Vector output_ldof_values(P_test_->Height(), mem_type); + for (int j = 0; j < output_ldof_values.Size(); ++j) { + output_ldof_values[j] = output_L_[j]; + } + P_test_->MultTranspose(output_ldof_values, output_T); + } else { + P_test_->MultTranspose(output_L_, output_T); + } } /** @@ -466,36 +586,54 @@ class Functional { // get the values for each local processor for (uint32_t i = 0; i < num_trial_spaces; i++) { -#if 0 if (trial_function_spaces_[i].family == Family::L2) { - - // make a PGF on the fly - // TODO: don't allocate/deallocate this data every invocation - mfem::ParGridFunction X; - X.MakeRef(trial_space_[i], input_L_[i].GetData()); - - // call exchange face nbr - X.ExchangeFaceNbrData(); - - // copy input_L[i] and facenbrdata [i] into a common array like: + // copy input_L[i] and facenbrdata[i] into a common array like: // first part second part // [ --- L --- | --- FND --- ] - mfem::Vector first_part; - first_part.MakeRef(input_L_[i], 0); - P_trial_[i]->Mult(*input_T[i], first_part); - - mfem::Vector second_part; - second_part.MakeRef(input_L_[i], first_part.Size()); - second_part = X.FaceNbrData(); - + P_trial_[i]->Mult(*input_T[i], input_ldof_values_[i]); + + // MakeRef from mfem::Vector will actually delete the data in the original vector, so we assign the + // values directly by SetVector, which invokes additional operations. This can be optimized later. + input_L_[i].SetVector(input_ldof_values_[i], 0); + + // Update again with new data + updateFaceNbrData(trial_space_[i], trial_pargrid_functions_[i], input_ldof_values_[i]); + + // Set the values of neighbor ghost dofs + // The entries in face_nbr_data vector is ALWAYS arranged byNODES per "volumetric" element. + // For example with two quadrilateral ghost elements we would have such face_nbr_data + // ----------- Ghost elem 1 ---------- ----------- Ghost elem 2 ---------- + // {X1 X1 X1 X1 Y1 Y1 Y1 Y1 Z1 Z1 Z1 Z1 X2 X2 X2 X2 Y2 Y2 Y2 Y2 Z2 Z2 Z2 Z2} + // + // This is because 1) mfem prepares face neighbor data element by element to communicate later (standard). + // 2) mfem uses GetElementVDofs to gather the data for communication which returns a set of indices ALWAYS + // ordered byNODES (???). If the ordering we set for ParFiniteElementSpace is byVDIM, we need the entries to be + // ----------- Ghost elem 1 ---------- ----------- Ghost elem 2 ---------- + // {X1 Y1 Z1 X1 Y1 Z1 X1 Y1 Z1 X1 Y1 Z1 X2 Y2 Z2 X2 Y2 Z2 X2 Y2 Z2 X2 Y2 Z2} + // so it's consistent with the local data. Therefore, we need to manually change this ordering. + if (trial_space_[i]->GetVDim() == 1) { + // For scalar fields, this weird ordering doesn't cause any problems. + input_L_[i].SetVector(trial_pargrid_functions_[i].FaceNbrData(), input_ldof_values_[i].Size()); + } else { + if (trial_space_[i]->GetOrdering() == mfem::Ordering::byVDIM) { + // For vector fields with byVDIM ordering, we manually set the entries in VDIM order. By the way this + // is really annoying because we have to get face neighbor vdof indices again in element restriction. + int local_size = input_ldof_values_[i].Size(); + appendFaceNbrData(trial_space_[i], trial_pargrid_functions_[i], local_size, input_L_[i]); + } else { + // For vector fields with byNODES ordering, we need to prepare the local input vector in the form + // { true_comp_X ghost_comp_X true_comp_Y ghost_comp_Y true_comp_Z ghost_comp_Z }. + // In this form, to ensure the correct mapping between dof and vdof, we need to change the ndofs in + // FiniteElementSpace to include ones from ghost element. Additionally, we need to prepare + // ghost_comp_X / ghost_comp_Y / ghost_comp_Z from element wise entries in face_nbr_data, + // eg. ghost_comp_X includes X component entries of all ghost element dofs. + // This requires significant renumbering of the dofs and therefore is not supported for now. + SLIC_ERROR_ROOT("Unsupported: L2 vector field ordered by nodes"); + } + } } else { - P_trial_[i]->Mult(*input_T[i], input_L_[i]); - } -#else - P_trial_[i]->Mult(*input_T[i], input_L_[i]); -#endif } output_L_ = 0.0; @@ -521,7 +659,17 @@ class Functional { } // scatter-add to compute global residuals - P_test_->MultTranspose(output_L_, output_T_); + if (test_function_space_.family == Family::L2) { + // Extract a subvector from output_L_ for local residuals excluding ghost dofs for L2 spaces + // The output vector is in form [ --- L --- | --- FND --- ] + mfem::Vector output_ldof_values(P_test_->Height(), mem_type); + for (int j = 0; j < output_ldof_values.Size(); ++j) { + output_ldof_values[j] = output_L_[j]; + } + P_test_->MultTranspose(output_ldof_values, output_T_); + } else { + P_test_->MultTranspose(output_L_, output_T_); + } if constexpr (wrt != NO_DIFFERENTIATION) { // if the user has indicated they'd like to evaluate and differentiate w.r.t. @@ -770,9 +918,51 @@ class Functional { auto* R = form_.test_space_->Dof_TrueDof_Matrix(); - auto* A_hypre = - new mfem::HypreParMatrix(test_space_->GetComm(), test_space_->GlobalVSize(), trial_space_->GlobalVSize(), - test_space_->GetDofOffsets(), trial_space_->GetDofOffsets(), &A_local); + // If either test_space_ or trial_space_ is L2, the local matrix (A_local) will includ ghost rows and columns like + // ------------ ------- + // | || | + // | diagonal || off-d | + // | block || block | + // | || | + // ------------ ------- + // ------------ ------- + // | ghost || XXX | + // | rows || XXX | + // ------------ ------- + // We only need + // ------------ ------- + // | || | + // | diagonal || off-d | + // | block || block | + // | || | + // ------------ ------- + // to compute local residual and the neighbor dof values will be communicated to multiply into off-digonal block. + // So we construct a HypreParMatrix that contains the off-diagonal block from the local sparse matrix + mfem::HypreParMatrix* A_hypre; + if (dynamic_cast(test_space_->FEColl()) || + dynamic_cast(trial_space_->FEColl())) { + int lrows = test_space_->GetVSize(); + int lcols = trial_space_->GetVSize(); + HYPRE_BigInt col_offset = trial_space_->GetMyDofOffset(); + mfem::Array glob_J(A_local.NumNonZeroElems()); + + int* J = A_local.GetJ(); + for (int j = 0; j < glob_J.Size(); ++j) { + if (J[j] < lcols) { + glob_J[j] = J[j] + col_offset; + } else { + glob_J[j] = form_.face_nbr_glob_vdof_maps_[which_argument][J[j] - lcols]; + } + } + A_hypre = new mfem::HypreParMatrix(test_space_->GetComm(), lrows, test_space_->GlobalVSize(), + trial_space_->GlobalVSize(), A_local.GetI(), glob_J, A_local.GetData(), + test_space_->GetDofOffsets(), trial_space_->GetDofOffsets()); + glob_J.DeleteAll(); + } else { + A_hypre = + new mfem::HypreParMatrix(test_space_->GetComm(), test_space_->GlobalVSize(), trial_space_->GlobalVSize(), + test_space_->GetDofOffsets(), trial_space_->GetDofOffsets(), &A_local); + } auto* P = trial_space_->Dof_TrueDof_Matrix(); @@ -823,6 +1013,15 @@ class Functional { std::array trial_function_spaces_; FunctionSpace test_function_space_; + /// @brief Manage global index of face neighbor vdofs + std::array, num_trial_spaces> face_nbr_glob_vdof_maps_; + + /// @brief Cache for access of neighbor element info for processor boundary faces + mutable std::array trial_pargrid_functions_; + + /// @brief Cache to store values of L2 ldofs on the local processor + mutable std::array input_ldof_values_; + /** * @brief Operator that converts true (global) DOF values to local (current rank) DOF values * for the test space diff --git a/src/serac/numerics/functional/functional_qoi.inl b/src/serac/numerics/functional/functional_qoi.inl index 76cfa97257..6244ff136b 100644 --- a/src/serac/numerics/functional/functional_qoi.inl +++ b/src/serac/numerics/functional/functional_qoi.inl @@ -88,29 +88,42 @@ class Functional { { SERAC_MARK_FUNCTION; + std::array trial_families = {trials::family...}; + std::array trial_orders = {trials::order...}; + std::array trial_components = {trials::components...}; + for (uint32_t i = 0; i < num_trial_spaces; i++) { + trial_function_spaces_[i] = {trial_families[i], trial_orders[i], trial_components[i]}; + P_trial_[i] = trial_space_[i]->GetProlongationMatrix(); - input_L_[i].SetSize(P_trial_[i]->Height(), mfem::Device::GetMemoryType()); + // For L2 spaces, configure a ParGridFunction with external data to exchange face neighbor data. + // This is for DG method where interior faces on the processor boundary need to access data on + // the neighrbor element + if (trial_function_spaces_[i].family == Family::L2) { + // Set the vector size to store tdof data + input_ldof_values_[i].SetSize(P_trial_[i]->Height(), mem_type); + + // Initialize face neighbor data vector + updateFaceNbrData(trial_space_[i], trial_pargrid_functions_[i], input_ldof_values_[i]); + + // Set the local input vector size to store local + ghost data + input_L_[i].SetSize(input_ldof_values_[i].Size() + trial_pargrid_functions_[i].FaceNbrData().Size(), mem_type); + } else { + input_L_[i].SetSize(P_trial_[i]->Height(), mem_type); + } // create the necessary number of empty mfem::Vectors, to be resized later input_E_.push_back({}); input_E_buffer_.push_back({}); } - std::array trial_families = {trials::family...}; - std::array trial_orders = {trials::order...}; - std::array trial_components = {trials::components...}; - for (uint32_t i = 0; i < num_trial_spaces; i++) { - trial_function_spaces_[i] = {trial_families[i], trial_orders[i], trial_components[i]}; - } - G_test_ = QoIElementRestriction(); P_test_ = QoIProlongation(trial_fes[0]->GetParMesh()->GetComm()); output_L_.SetSize(1, mem_type); - output_T_.SetSize(1, mfem::Device::GetMemoryType()); + output_T_.SetSize(1, mem_type); // gradient objects depend on some member variables in // Functional, so we initialize the gradient objects last @@ -202,7 +215,6 @@ class Functional { std::vector arg_vec = {args...}; for (uint32_t i : arg_vec) { domain.insert_restriction(trial_space_[i], trial_function_spaces_[i]); - check_interior_face_compatibility(domain.mesh_, trial_function_spaces_[i]); } using signature = test(decltype(serac::type(trial_spaces))...); @@ -262,7 +274,25 @@ class Functional { */ double ActionOfGradient(const mfem::Vector& input_T, uint32_t which) const { - P_trial_[which]->Mult(input_T, input_L_[which]); + // Please refer to 'ActionOfGradient' in 'functional.hpp' for detailed comments of this implementation + if (trial_function_spaces_[which].family == Family::L2) { + P_trial_[which]->Mult(input_T, input_ldof_values_[which]); + input_L_[which].SetVector(input_ldof_values_[which], 0); + updateFaceNbrData(trial_space_[which], trial_pargrid_functions_[which], input_ldof_values_[which]); + + if (trial_space_[which]->GetVDim() == 1) { + input_L_[which].SetVector(trial_pargrid_functions_[which].FaceNbrData(), input_ldof_values_[which].Size()); + } else { + if (trial_space_[which]->GetOrdering() == mfem::Ordering::byVDIM) { + int local_size = input_ldof_values_[which].Size(); + appendFaceNbrData(trial_space_[which], trial_pargrid_functions_[which], local_size, input_L_[which]); + } else { + SLIC_ERROR_ROOT("Unsupported: L2 vector field ordered by nodes"); + } + } + } else { + P_trial_[which]->Mult(input_T, input_L_[which]); + } output_L_ = 0.0; @@ -280,6 +310,13 @@ class Functional { integral.GradientMult(input_E_[which], output_E_, which); + // make sure shared interior faces are integrated only once + if (dom.type_ == Domain::Type::InteriorFaces) { + for (int shared_id : dom.shared_interior_face_ids_) { + output_E_[shared_id] *= 0.5; + } + } + // scatter-add to compute QoI value for the local processor G_test_.ScatterAdd(output_E_, output_L_); } @@ -306,9 +343,27 @@ class Functional { { const mfem::Vector* input_T[] = {&static_cast(args)...}; + // Please refer to 'operator()' in 'functional.hpp' for detailed comments of this implementation // get the values for each local processor for (uint32_t i = 0; i < num_trial_spaces; i++) { - P_trial_[i]->Mult(*input_T[i], input_L_[i]); + if (trial_function_spaces_[i].family == Family::L2) { + P_trial_[i]->Mult(*input_T[i], input_ldof_values_[i]); + input_L_[i].SetVector(input_ldof_values_[i], 0); + updateFaceNbrData(trial_space_[i], trial_pargrid_functions_[i], input_ldof_values_[i]); + + if (trial_space_[i]->GetVDim() == 1) { + input_L_[i].SetVector(trial_pargrid_functions_[i].FaceNbrData(), input_ldof_values_[i].Size()); + } else { + if (trial_space_[i]->GetOrdering() == mfem::Ordering::byVDIM) { + int local_size = input_ldof_values_[i].Size(); + appendFaceNbrData(trial_space_[i], trial_pargrid_functions_[i], local_size, input_L_[i]); + } else { + SLIC_ERROR_ROOT("Unsupported: L2 vector field ordered by nodes"); + } + } + } else { + P_trial_[i]->Mult(*input_T[i], input_L_[i]); + } } output_L_ = 0.0; @@ -329,6 +384,13 @@ class Functional { const bool update_qdata = false; integral.Mult(t, input_E_, output_E_, wrt, update_qdata); + // make sure shared interior faces are integrated only once + if (dom.type_ == Domain::Type::InteriorFaces) { + for (int shared_id : dom.shared_interior_face_ids_) { + output_E_[shared_id] *= 0.5; + } + } + // scatter-add to compute QoI value for the local processor G_test_.ScatterAdd(output_E_, output_L_); } @@ -427,6 +489,8 @@ class Functional { std::map> element_gradients[Domain::num_types]; + int lcols = form_.trial_space_[which_argument]->GetVSize(); + //////////////////////////////////////////////////////////////////////////////// for (auto& integral : form_.integrals_) { @@ -456,7 +520,11 @@ class Functional { for (uint32_t i = 0; i < cols_per_elem; i++) { int col = int(trial_vdofs[i].index()); - gradient_L_[col] += K_e(e, 0, i); + + // only add local DG dof gradients (other FE spaces satisfy col < lcols inherently) + if (col < lcols) { + gradient_L_[col] += K_e(e, 0, i); + } } } } @@ -492,6 +560,12 @@ class Functional { std::array trial_function_spaces_; + /// @brief Cache for access of neighbor element info for processor boundary faces + mutable std::array trial_pargrid_functions_; + + /// @brief Cache to store values of L2 ldofs on the local processor + mutable std::array input_ldof_values_; + /** * @brief Operator that converts true (global) DOF values to local (current rank) DOF values * for the test space diff --git a/src/serac/numerics/functional/geometric_factors.cpp b/src/serac/numerics/functional/geometric_factors.cpp index 972fef3d34..d41a7645b0 100644 --- a/src/serac/numerics/functional/geometric_factors.cpp +++ b/src/serac/numerics/functional/geometric_factors.cpp @@ -67,10 +67,10 @@ void compute_geometric_factors(mfem::Vector& positions_q, mfem::Vector& jacobian GeometricFactors::GeometricFactors(const Domain& domain, int q, mfem::Geometry::Type geom) { - // const mfem::ParGridFunction* nodes = static_cast< const mfem::ParGridFunction * >(domain.mesh_.GetNodes()); - // mfem::ParFiniteElementSpace * pfes = nodes->ParFESpace(); - const mfem::GridFunction* nodes = domain.mesh_.GetNodes(); - const mfem::FiniteElementSpace* fes = nodes->FESpace(); + const mfem::ParGridFunction* nodes = static_cast(domain.mesh_.GetNodes()); + mfem::ParFiniteElementSpace* fes = nodes->ParFESpace(); + // const mfem::GridFunction* nodes = domain.mesh_.GetNodes(); + // const mfem::FiniteElementSpace* fes = nodes->FESpace(); const std::vector& element_ids = domain.get_mfem_ids(geom); diff --git a/src/serac/numerics/functional/tests/CMakeLists.txt b/src/serac/numerics/functional/tests/CMakeLists.txt index a3f8b34eee..cf101153a4 100644 --- a/src/serac/numerics/functional/tests/CMakeLists.txt +++ b/src/serac/numerics/functional/tests/CMakeLists.txt @@ -20,7 +20,6 @@ set(functional_serial_test_sources simplex_basis_function_unit_tests.cpp element_restriction_tests.cpp dg_restriction_operators.cpp - functional_basic_dg.cpp bug_boundary_qoi.cpp domain_tests.cpp geometric_factors_tests.cpp @@ -41,12 +40,14 @@ set(functional_parallel_test_sources functional_with_domain.cpp functional_basic_h1_scalar.cpp functional_basic_h1_vector.cpp + functional_basic_dg.cpp functional_multiphysics.cpp functional_qoi.cpp functional_nonlinear.cpp functional_boundary_test.cpp functional_comparisons.cpp functional_comparison_L2.cpp + dg_ghost_face_index.cpp ) serac_add_tests(SOURCES ${functional_parallel_test_sources} diff --git a/src/serac/numerics/functional/tests/dg_ghost_face_index.cpp b/src/serac/numerics/functional/tests/dg_ghost_face_index.cpp new file mode 100644 index 0000000000..7aaff324c4 --- /dev/null +++ b/src/serac/numerics/functional/tests/dg_ghost_face_index.cpp @@ -0,0 +1,115 @@ +// Copyright (c) Lawrence Livermore National Security, LLC and +// other Serac Project Developers. See the top-level LICENSE file for +// details. +// +// SPDX-License-Identifier: (BSD-3-Clause) + +#include +#include + +#include "mfem.hpp" + +#include + +#include "serac/infrastructure/application_manager.hpp" +#include "serac/serac_config.hpp" +#include "serac/mesh_utils/mesh_utils_base.hpp" +#include "serac/numerics/stdfunction_operator.hpp" +#include "serac/numerics/functional/functional.hpp" +#include "serac/numerics/functional/tensor.hpp" + +using namespace serac; +using namespace serac::profiling; + +// This test initializes a DG field with nodal coordinates of the dofs, so that the +// discontinuous dof pairs across the interior faces have the value. For example +// on the interior face with the following dofs +// {1, 2} | {5, 6} +// | +// | +// {3, 4} | {7, 8} +// we have {1, 2} = {5, 6} and {3, 4} = {7, 8}. +// It then integrates the jump of dof values over all interior faces. +// If the ghost dof data is constructed correctly to align with locally owned data, +// then every entry in the residual vector should equal to zero. This is tested +// by the L2 norm of the residual equal to zero. +template +void L2_index_test(std::string meshfile) +{ + using test_space = L2; + using trial_space = L2; + + auto mesh = mesh::refineAndDistribute(buildMeshFromFile(meshfile), 1); + + auto [test_fespace, test_fec] = serac::generateParFiniteElementSpace(mesh.get()); + auto [trial_fespace, trial_fec] = serac::generateParFiniteElementSpace(mesh.get()); + + // Initialize the ParGridFunction by dof coordinates + mfem::ParGridFunction U_gf(trial_fespace.get()); + mfem::VectorFunctionCoefficient vcoef(dim, [](const mfem::Vector& X, mfem::Vector& F) { + int d = X.Size(); + F.SetSize(d); + for (int i = 0; i < d; ++i) { + F(i) = X(i); + } + }); + U_gf.ProjectCoefficient(vcoef); + + mfem::Vector U(trial_fespace->TrueVSize()); + U_gf.GetTrueDofs(U); + + // Construct the new functional object using the specified test and trial spaces + Functional residual(test_fespace.get(), {trial_fespace.get()}); + + Domain interior_faces = InteriorFaces(*mesh); + + // Define the integral of jumps over all interior faces + residual.AddInteriorFaceIntegral( + Dimension{}, DependsOn<0>{}, + [=](double /*t*/, auto X, auto velocity) { + // compute the surface normal + auto dX_dxi = get(X); + auto n = normalize(cross(dX_dxi)); + + // extract the velocity values from each side of the interface + // note: the orientation convention is such that the normal + // computed as above will point from from side 1->2 + auto [u_1, u_2] = velocity; + SLIC_INFO(axom::fmt::format("One size = {}, The other side = {}, Jump = {}", axom::fmt::streamed(u_1), + axom::fmt::streamed(u_2), axom::fmt::streamed(u_1 - u_2))); + + auto a = dot(u_2 - u_1, n); + + auto f_1 = u_1 * a; + auto f_2 = u_2 * a; + return serac::tuple{f_1, f_2}; + }, + interior_faces); + + double t = 0.0; + + auto value = residual(t, U); + EXPECT_NEAR(0., value.Norml2(), 1.e-12); +} + +TEST(index, L2_test_tris_and_quads_linear) +{ + L2_index_test<2, 1>(SERAC_REPO_DIR "/data/meshes/patch2D_tris_and_quads.mesh"); +} +TEST(index, L2_test_tris_and_quads_quadratic) +{ + L2_index_test<2, 2>(SERAC_REPO_DIR "/data/meshes/patch2D_tris_and_quads.mesh"); +} + +TEST(index, L2_test_tets_linear) { L2_index_test<3, 1>(SERAC_REPO_DIR "/data/meshes/patch3D_tets.mesh"); } +TEST(index, L2_test_tets_quadratic) { L2_index_test<3, 2>(SERAC_REPO_DIR "/data/meshes/patch3D_tets.mesh"); } + +TEST(index, L2_test_hexes_linear) { L2_index_test<3, 1>(SERAC_REPO_DIR "/data/meshes/patch3D_hexes.mesh"); } +TEST(index, L2_test_hexes_quadratic) { L2_index_test<3, 2>(SERAC_REPO_DIR "/data/meshes/patch3D_hexes.mesh"); } + +int main(int argc, char* argv[]) +{ + ::testing::InitGoogleTest(&argc, argv); + serac::ApplicationManager applicationManager(argc, argv); + return RUN_ALL_TESTS(); +} diff --git a/src/serac/numerics/functional/tests/dg_restriction_operators.cpp b/src/serac/numerics/functional/tests/dg_restriction_operators.cpp index 49def5842c..e3762202c1 100644 --- a/src/serac/numerics/functional/tests/dg_restriction_operators.cpp +++ b/src/serac/numerics/functional/tests/dg_restriction_operators.cpp @@ -7,6 +7,7 @@ #include #include "serac/mesh_utils/mesh_utils_base.hpp" +#include "serac/mesh_utils/mesh_utils.hpp" #include "serac/numerics/functional/domain.hpp" #include "serac/numerics/functional/element_restriction.hpp" #include "serac/numerics/functional/functional.hpp" @@ -286,13 +287,14 @@ void parametrized_test(int permutation) constexpr mfem::Geometry::Type face_geom = face_type(geom); constexpr int dim = dimension_of(geom); - mfem::Mesh mesh = generate_permuted_mesh(geom, permutation); + mfem::Mesh bmesh = generate_permuted_mesh(geom, permutation); + auto mesh = serac::mesh::refineAndDistribute(std::move(bmesh)); - Domain interior_faces = InteriorFaces(mesh); + Domain interior_faces = InteriorFaces(*mesh); // each one of these meshes should have two elements // and a single "face" that separates them - EXPECT_EQ(mesh.GetNE(), 2); + EXPECT_EQ(mesh->GetNE(), 2); std::vector face_ids; @@ -318,13 +320,13 @@ void parametrized_test(int permutation) auto Hcurl_fec = std::make_unique(polynomial_order, dim); auto L2_fec = std::make_unique(polynomial_order, dim, mfem::BasisType::GaussLobatto); - auto H1_fes = std::make_unique(&mesh, H1_fec.get()); - auto Hcurl_fes = std::make_unique(&mesh, Hcurl_fec.get()); - auto L2_fes = std::make_unique(&mesh, L2_fec.get()); + auto H1_fes = std::make_unique(mesh.get(), H1_fec.get()); + auto Hcurl_fes = std::make_unique(mesh.get(), Hcurl_fec.get()); + auto L2_fes = std::make_unique(mesh.get(), L2_fec.get()); - mfem::GridFunction H1_gf(H1_fes.get()); - mfem::GridFunction Hcurl_gf(Hcurl_fes.get()); - mfem::GridFunction L2_gf(L2_fes.get()); + mfem::ParGridFunction H1_gf(H1_fes.get()); + mfem::ParGridFunction Hcurl_gf(Hcurl_fes.get()); + mfem::ParGridFunction L2_gf(L2_fes.get()); mfem::FunctionCoefficient sfunc(scalar_func_mfem); mfem::VectorFunctionCoefficient vfunc(dim, vector_func_mfem); @@ -351,12 +353,7 @@ void parametrized_test(int permutation) using space = L2; - auto pmesh = mesh::refineAndDistribute(std::move(mesh), 0); - Domain pinterior_faces = InteriorFaces(*pmesh); - - auto L2_pfes = mfem::ParFiniteElementSpace(pmesh.get(), L2_fec.get()); - - Functional r({&L2_pfes}); + Functional r({L2_fes.get()}); r.AddInteriorFaceIntegral( Dimension{}, DependsOn<0>{}, @@ -369,10 +366,10 @@ void parametrized_test(int permutation) auto difference = rho1 - rho2; return difference * difference; }, - pinterior_faces); + interior_faces); - mfem::Vector U(L2_pfes.TrueVSize()); - mfem::ParGridFunction U_pgf(&L2_pfes); + mfem::Vector U(L2_fes->TrueVSize()); + mfem::ParGridFunction U_pgf(L2_fes.get()); U_pgf.ProjectCoefficient(sfunc); U_pgf.GetTrueDofs(U); diff --git a/src/serac/numerics/functional/tests/domain_tests.cpp b/src/serac/numerics/functional/tests/domain_tests.cpp index 35e9cb92eb..5ce35c8553 100644 --- a/src/serac/numerics/functional/tests/domain_tests.cpp +++ b/src/serac/numerics/functional/tests/domain_tests.cpp @@ -8,6 +8,7 @@ #include "serac/numerics/functional/domain.hpp" #include "serac/infrastructure/application_manager.hpp" +#include "serac/mesh_utils/mesh_utils.hpp" using namespace serac; @@ -31,16 +32,17 @@ mfem::Mesh import_mesh(std::string meshfile) TEST(domain, of_edges) { { - auto mesh = import_mesh("onehex.mesh"); - Domain d0 = Domain::ofEdges(mesh, std::function([](std::vector x) { - return (0.5 * (x[0][0] + x[1][0])) < 0.25; // x coordinate of edge midpoint - })); + auto bmesh = import_mesh("onehex.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(bmesh)); + Domain d0 = Domain::ofEdges(*mesh, std::function([](std::vector x) { + return (0.5 * (x[0][0] + x[1][0])) < 0.25; // x coordinate of edge midpoint + })); EXPECT_EQ(d0.edge_ids_.size(), 4); EXPECT_EQ(d0.dim_, 1); - Domain d1 = Domain::ofEdges(mesh, std::function([](std::vector x) { - return (0.5 * (x[0][1] + x[1][1])) < 0.25; // y coordinate of edge midpoint - })); + Domain d1 = Domain::ofEdges(*mesh, std::function([](std::vector x) { + return (0.5 * (x[0][1] + x[1][1])) < 0.25; // y coordinate of edge midpoint + })); EXPECT_EQ(d1.edge_ids_.size(), 4); EXPECT_EQ(d1.dim_, 1); @@ -58,16 +60,17 @@ TEST(domain, of_edges) } { - auto mesh = import_mesh("onetet.mesh"); - Domain d0 = Domain::ofEdges(mesh, std::function([](std::vector x) { - return (0.5 * (x[0][0] + x[1][0])) < 0.25; // x coordinate of edge midpoint - })); + auto bmesh = import_mesh("onetet.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(bmesh)); + Domain d0 = Domain::ofEdges(*mesh, std::function([](std::vector x) { + return (0.5 * (x[0][0] + x[1][0])) < 0.25; // x coordinate of edge midpoint + })); EXPECT_EQ(d0.edge_ids_.size(), 3); EXPECT_EQ(d0.dim_, 1); - Domain d1 = Domain::ofEdges(mesh, std::function([](std::vector x) { - return (0.5 * (x[0][1] + x[1][1])) < 0.25; // y coordinate of edge midpoint - })); + Domain d1 = Domain::ofEdges(*mesh, std::function([](std::vector x) { + return (0.5 * (x[0][1] + x[1][1])) < 0.25; // y coordinate of edge midpoint + })); EXPECT_EQ(d1.edge_ids_.size(), 3); EXPECT_EQ(d1.dim_, 1); @@ -86,17 +89,18 @@ TEST(domain, of_edges) { constexpr int dim = 2; - auto mesh = import_mesh("beam-quad.mesh"); - mesh.FinalizeQuadMesh(true); - Domain d0 = Domain::ofEdges(mesh, std::function([](std::vector x, int /* bdr_attr */) { - return (0.5 * (x[0][0] + x[1][0])) < 0.25; // x coordinate of edge midpoint - })); + auto bmesh = import_mesh("beam-quad.mesh"); + bmesh.FinalizeQuadMesh(true); + auto mesh = serac::mesh::refineAndDistribute(std::move(bmesh)); + Domain d0 = Domain::ofEdges(*mesh, std::function([](std::vector x, int /* bdr_attr */) { + return (0.5 * (x[0][0] + x[1][0])) < 0.25; // x coordinate of edge midpoint + })); EXPECT_EQ(d0.edge_ids_.size(), 1); EXPECT_EQ(d0.dim_, 1); - Domain d1 = Domain::ofEdges(mesh, std::function([](std::vector x, int /* bdr_attr */) { - return (0.5 * (x[0][1] + x[1][1])) < 0.25; // y coordinate of edge midpoint - })); + Domain d1 = Domain::ofEdges(*mesh, std::function([](std::vector x, int /* bdr_attr */) { + return (0.5 * (x[0][1] + x[1][1])) < 0.25; // y coordinate of edge midpoint + })); EXPECT_EQ(d1.edge_ids_.size(), 8); EXPECT_EQ(d1.dim_, 1); @@ -109,10 +113,10 @@ TEST(domain, of_edges) EXPECT_EQ(d3.dim_, 1); // check that by_attr compiles - Domain d4 = Domain::ofEdges(mesh, by_attr(3)); + Domain d4 = Domain::ofEdges(*mesh, by_attr(3)); EXPECT_EQ(d4.mfem_edge_ids_.size(), 16); - Domain d5 = Domain::ofBoundaryElements(mesh, [](std::vector, int) { return true; }); + Domain d5 = Domain::ofBoundaryElements(*mesh, [](std::vector, int) { return true; }); EXPECT_EQ(d5.edge_ids_.size(), 18); // 1x8 row of quads has 18 boundary edges } } @@ -121,16 +125,17 @@ TEST(domain, of_faces) { { constexpr int dim = 3; - auto mesh = import_mesh("onehex.mesh"); - Domain d0 = Domain::ofFaces(mesh, std::function([](std::vector vertices, int /*bdr_attr*/) { - return average(vertices)[0] < 0.25; // x coordinate of face center - })); + auto bmesh = import_mesh("onehex.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(bmesh)); + Domain d0 = Domain::ofFaces(*mesh, std::function([](std::vector vertices, int /*bdr_attr*/) { + return average(vertices)[0] < 0.25; // x coordinate of face center + })); EXPECT_EQ(d0.quad_ids_.size(), 1); EXPECT_EQ(d0.dim_, 2); - Domain d1 = Domain::ofFaces(mesh, std::function([](std::vector vertices, int /*bdr_attr*/) { - return average(vertices)[1] < 0.25; // y coordinate of face center - })); + Domain d1 = Domain::ofFaces(*mesh, std::function([](std::vector vertices, int /*bdr_attr*/) { + return average(vertices)[1] < 0.25; // y coordinate of face center + })); EXPECT_EQ(d1.quad_ids_.size(), 1); EXPECT_EQ(d1.dim_, 2); @@ -143,27 +148,28 @@ TEST(domain, of_faces) EXPECT_EQ(d3.dim_, 2); // check that by_attr compiles - Domain d4 = Domain::ofFaces(mesh, by_attr(3)); + Domain d4 = Domain::ofFaces(*mesh, by_attr(3)); - Domain d5 = Domain::ofBoundaryElements(mesh, [](std::vector, int) { return true; }); + Domain d5 = Domain::ofBoundaryElements(*mesh, [](std::vector, int) { return true; }); EXPECT_EQ(d5.quad_ids_.size(), 6); } { constexpr int dim = 3; - auto mesh = import_mesh("onetet.mesh"); - Domain d0 = Domain::ofFaces(mesh, std::function([](std::vector vertices, int /* bdr_attr */) { - // accept face if it contains a vertex whose x coordinate is less than 0.1 - for (auto v : vertices) { - if (v[0] < 0.1) return true; - } - return false; - })); + auto bmesh = import_mesh("onetet.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(bmesh)); + Domain d0 = Domain::ofFaces(*mesh, std::function([](std::vector vertices, int /* bdr_attr */) { + // accept face if it contains a vertex whose x coordinate is less than 0.1 + for (auto v : vertices) { + if (v[0] < 0.1) return true; + } + return false; + })); EXPECT_EQ(d0.tri_ids_.size(), 4); EXPECT_EQ(d0.dim_, 2); Domain d1 = Domain::ofFaces( - mesh, std::function([](std::vector x, int /* bdr_attr */) { return average(x)[1] < 0.1; })); + *mesh, std::function([](std::vector x, int /* bdr_attr */) { return average(x)[1] < 0.1; })); EXPECT_EQ(d1.tri_ids_.size(), 1); EXPECT_EQ(d1.dim_, 2); @@ -176,24 +182,25 @@ TEST(domain, of_faces) EXPECT_EQ(d3.dim_, 2); // check that by_attr compiles - Domain d4 = Domain::ofFaces(mesh, by_attr(3)); + Domain d4 = Domain::ofFaces(*mesh, by_attr(3)); - Domain d5 = Domain::ofBoundaryElements(mesh, [](std::vector, int) { return true; }); + Domain d5 = Domain::ofBoundaryElements(*mesh, [](std::vector, int) { return true; }); EXPECT_EQ(d5.tri_ids_.size(), 4); } { constexpr int dim = 2; - auto mesh = import_mesh("beam-quad.mesh"); - Domain d0 = Domain::ofFaces(mesh, std::function([](std::vector vertices, int /* attr */) { - return average(vertices)[0] < 2.25; // x coordinate of face center - })); + auto bmesh = import_mesh("beam-quad.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(bmesh)); + Domain d0 = Domain::ofFaces(*mesh, std::function([](std::vector vertices, int /* attr */) { + return average(vertices)[0] < 2.25; // x coordinate of face center + })); EXPECT_EQ(d0.quad_ids_.size(), 2); EXPECT_EQ(d0.dim_, 2); - Domain d1 = Domain::ofFaces(mesh, std::function([](std::vector vertices, int /* attr */) { - return average(vertices)[1] < 0.55; // y coordinate of face center - })); + Domain d1 = Domain::ofFaces(*mesh, std::function([](std::vector vertices, int /* attr */) { + return average(vertices)[1] < 0.55; // y coordinate of face center + })); EXPECT_EQ(d1.quad_ids_.size(), 8); EXPECT_EQ(d1.dim_, 2); @@ -206,7 +213,7 @@ TEST(domain, of_faces) EXPECT_EQ(d3.dim_, 2); // check that by_attr compiles - Domain d4 = Domain::ofFaces(mesh, by_attr(3)); + Domain d4 = Domain::ofFaces(*mesh, by_attr(3)); } } @@ -214,18 +221,19 @@ TEST(domain, of_elements) { { constexpr int dim = 3; - auto mesh = import_mesh("patch3D_tets_and_hexes.mesh"); - Domain d0 = Domain::ofElements(mesh, std::function([](std::vector vertices, int /*bdr_attr*/) { - return average(vertices)[0] < 0.7; // x coordinate of face center - })); + auto bmesh = import_mesh("patch3D_tets_and_hexes.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(bmesh)); + Domain d0 = Domain::ofElements(*mesh, std::function([](std::vector vertices, int /*bdr_attr*/) { + return average(vertices)[0] < 0.7; // x coordinate of face center + })); EXPECT_EQ(d0.tet_ids_.size(), 0); EXPECT_EQ(d0.hex_ids_.size(), 1); EXPECT_EQ(d0.dim_, 3); - Domain d1 = Domain::ofElements(mesh, std::function([](std::vector vertices, int /*bdr_attr*/) { - return average(vertices)[1] < 0.75; // y coordinate of face center - })); + Domain d1 = Domain::ofElements(*mesh, std::function([](std::vector vertices, int /*bdr_attr*/) { + return average(vertices)[1] < 0.75; // y coordinate of face center + })); EXPECT_EQ(d1.tet_ids_.size(), 6); EXPECT_EQ(d1.hex_ids_.size(), 1); EXPECT_EQ(d1.dim_, 3); @@ -241,20 +249,21 @@ TEST(domain, of_elements) EXPECT_EQ(d3.dim_, 3); // check that by_attr works - Domain d4 = Domain::ofElements(mesh, by_attr(3)); + Domain d4 = Domain::ofElements(*mesh, by_attr(3)); } { constexpr int dim = 2; - auto mesh = import_mesh("patch2D_tris_and_quads.mesh"); + auto bmesh = import_mesh("patch2D_tris_and_quads.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(bmesh)); Domain d0 = Domain::ofElements( - mesh, std::function([](std::vector vertices, int /* attr */) { return average(vertices)[0] < 0.45; })); + *mesh, std::function([](std::vector vertices, int /* attr */) { return average(vertices)[0] < 0.45; })); EXPECT_EQ(d0.tri_ids_.size(), 1); EXPECT_EQ(d0.quad_ids_.size(), 1); EXPECT_EQ(d0.dim_, 2); Domain d1 = Domain::ofElements( - mesh, std::function([](std::vector vertices, int /* attr */) { return average(vertices)[1] < 0.45; })); + *mesh, std::function([](std::vector vertices, int /* attr */) { return average(vertices)[1] < 0.45; })); EXPECT_EQ(d1.tri_ids_.size(), 1); EXPECT_EQ(d1.quad_ids_.size(), 1); EXPECT_EQ(d1.dim_, 2); @@ -270,7 +279,7 @@ TEST(domain, of_elements) EXPECT_EQ(d3.dim_, 2); // check that by_attr compiles - Domain d4 = Domain::ofElements(mesh, by_attr(3)); + Domain d4 = Domain::ofElements(*mesh, by_attr(3)); } } @@ -278,16 +287,17 @@ TEST(domain, entireDomain2d) { constexpr int dim = 2; constexpr int p = 1; - auto mesh = import_mesh("patch2D_tris_and_quads.mesh"); + auto bmesh = import_mesh("patch2D_tris_and_quads.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(bmesh)); - Domain d0 = EntireDomain(mesh); + Domain d0 = EntireDomain(*mesh); EXPECT_EQ(d0.dim_, 2); EXPECT_EQ(d0.tri_ids_.size(), 2); EXPECT_EQ(d0.quad_ids_.size(), 4); auto fec = mfem::H1_FECollection(p, dim); - auto fes = mfem::FiniteElementSpace(&mesh, &fec); + auto fes = mfem::ParFiniteElementSpace(mesh.get(), &fec); mfem::Array dof_indices = d0.dof_list(&fes); EXPECT_EQ(dof_indices.Size(), 8); @@ -297,16 +307,17 @@ TEST(domain, entireDomain3d) { constexpr int dim = 3; constexpr int p = 1; - auto mesh = import_mesh("patch3D_tets_and_hexes.mesh"); + auto bmesh = import_mesh("patch3D_tets_and_hexes.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(bmesh)); - Domain d0 = EntireDomain(mesh); + Domain d0 = EntireDomain(*mesh); EXPECT_EQ(d0.dim_, 3); EXPECT_EQ(d0.tet_ids_.size(), 12); EXPECT_EQ(d0.hex_ids_.size(), 7); auto fec = mfem::H1_FECollection(p, dim); - auto fes = mfem::FiniteElementSpace(&mesh, &fec); + auto fes = mfem::ParFiniteElementSpace(mesh.get(), &fec); mfem::Array dof_indices = d0.dof_list(&fes); EXPECT_EQ(dof_indices.Size(), 25); @@ -316,17 +327,18 @@ TEST(domain, of2dElementsFindsDofs) { constexpr int dim = 2; constexpr int p = 2; - auto mesh = import_mesh("patch2D_tris_and_quads.mesh"); + auto bmesh = import_mesh("patch2D_tris_and_quads.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(bmesh)); auto fec = mfem::H1_FECollection(p, dim); - auto fes = mfem::FiniteElementSpace(&mesh, &fec); + auto fes = mfem::ParFiniteElementSpace(mesh.get(), &fec); auto find_element_0 = [](std::vector vertices, int /* attr */) { auto centroid = average(vertices); return (centroid[0] < 0.5) && (centroid[1] < 0.25); }; - Domain d0 = Domain::ofElements(mesh, find_element_0); + Domain d0 = Domain::ofElements(*mesh, find_element_0); mfem::Array dof_indices = d0.dof_list(&fes); @@ -339,7 +351,7 @@ TEST(domain, of2dElementsFindsDofs) tensor target{{0.533, 0.424}}; return norm(centroid - target) < 1e-2; }; - Domain d1 = Domain::ofElements(mesh, find_element_4); + Domain d1 = Domain::ofElements(*mesh, find_element_4); Domain elements_0_and_4 = d0 | d1; @@ -348,7 +360,7 @@ TEST(domain, of2dElementsFindsDofs) /////////////////////////////////////// - Domain d2 = EntireDomain(mesh) - elements_0_and_4; + Domain d2 = EntireDomain(*mesh) - elements_0_and_4; dof_indices = d2.dof_list(&fes); @@ -359,10 +371,11 @@ TEST(domain, of3dElementsFindsDofs) { constexpr int dim = 3; constexpr int p = 2; - auto mesh = import_mesh("patch3D_tets_and_hexes.mesh"); + auto bmesh = import_mesh("patch3D_tets_and_hexes.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(bmesh)); auto fec = mfem::H1_FECollection(p, dim); - auto fes = mfem::FiniteElementSpace(&mesh, &fec); + auto fes = mfem::ParFiniteElementSpace(mesh.get(), &fec); auto find_element_0 = [](std::vector vertices, int /* attr */) { auto centroid = average(vertices); @@ -370,7 +383,7 @@ TEST(domain, of3dElementsFindsDofs) return norm(centroid - target) < 1e-2; }; - Domain d0 = Domain::ofElements(mesh, find_element_0); + Domain d0 = Domain::ofElements(*mesh, find_element_0); mfem::Array dof_indices = d0.dof_list(&fes); @@ -384,7 +397,7 @@ TEST(domain, of3dElementsFindsDofs) vec3 target{{3.275, 1.2, 0.725}}; return norm(centroid - target) < 1e-2; }; - Domain d1 = Domain::ofElements(mesh, find_element_1); + Domain d1 = Domain::ofElements(*mesh, find_element_1); Domain elements_0_and_1 = d0 | d1; @@ -395,7 +408,7 @@ TEST(domain, of3dElementsFindsDofs) ///////////////////////////////////////// - Domain d2 = EntireDomain(mesh) - elements_0_and_1; + Domain d2 = EntireDomain(*mesh) - elements_0_and_1; dof_indices = d2.dof_list(&fes); @@ -406,17 +419,18 @@ TEST(domain, of2dBoundaryElementsFindsDofs) { constexpr int dim = 2; constexpr int p = 2; - auto mesh = import_mesh("patch2D_tris_and_quads.mesh"); + auto bmesh = import_mesh("patch2D_tris_and_quads.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(bmesh)); auto find_right_boundary = [](std::vector vertices, int /* attr */) { return std::all_of(vertices.begin(), vertices.end(), [](vec2 X) { return X[0] > 1.0 - 1e-2; }); }; - Domain d0 = Domain::ofBoundaryElements(mesh, find_right_boundary); + Domain d0 = Domain::ofBoundaryElements(*mesh, find_right_boundary); EXPECT_EQ(d0.edge_ids_.size(), 1); auto fec = mfem::H1_FECollection(p, dim); - auto fes = mfem::FiniteElementSpace(&mesh, &fec); + auto fes = mfem::ParFiniteElementSpace(mesh.get(), &fec); mfem::Array dof_indices = d0.dof_list(&fes); @@ -426,7 +440,7 @@ TEST(domain, of2dBoundaryElementsFindsDofs) return std::all_of(vertices.begin(), vertices.end(), [](vec2 X) { return X[1] > 1.0 - 1e-2; }); }; - Domain d1 = Domain::ofBoundaryElements(mesh, find_top_boundary); + Domain d1 = Domain::ofBoundaryElements(*mesh, find_top_boundary); EXPECT_EQ(d1.edge_ids_.size(), 1); Domain d2 = d0 | d1; @@ -440,17 +454,18 @@ TEST(domain, of3dBoundaryElementsFindsDofs) { constexpr int dim = 3; constexpr int p = 2; - auto mesh = import_mesh("patch3D_tets.mesh"); + auto bmesh = import_mesh("patch3D_tets.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(bmesh)); auto find_xmax_boundary = [](std::vector vertices, int /* attr */) { return std::all_of(vertices.begin(), vertices.end(), [](vec3 X) { return X[0] > 1.0 - 1e-2; }); }; - Domain d0 = Domain::ofBoundaryElements(mesh, find_xmax_boundary); + Domain d0 = Domain::ofBoundaryElements(*mesh, find_xmax_boundary); EXPECT_EQ(d0.tri_ids_.size(), 2); auto fec = mfem::H1_FECollection(p, dim); - auto fes = mfem::FiniteElementSpace(&mesh, &fec); + auto fes = mfem::ParFiniteElementSpace(mesh.get(), &fec); mfem::Array dof_indices = d0.dof_list(&fes); @@ -460,7 +475,7 @@ TEST(domain, of3dBoundaryElementsFindsDofs) return std::all_of(vertices.begin(), vertices.end(), [](vec3 X) { return X[1] > 1.0 - 1e-2; }); }; - Domain d1 = Domain::ofBoundaryElements(mesh, find_ymax_boundary); + Domain d1 = Domain::ofBoundaryElements(*mesh, find_ymax_boundary); EXPECT_EQ(d1.tri_ids_.size(), 2); Domain d2 = d0 | d1; diff --git a/src/serac/numerics/functional/tests/element_restriction_tests.cpp b/src/serac/numerics/functional/tests/element_restriction_tests.cpp index 361ccfca14..2da708fb75 100644 --- a/src/serac/numerics/functional/tests/element_restriction_tests.cpp +++ b/src/serac/numerics/functional/tests/element_restriction_tests.cpp @@ -8,6 +8,8 @@ #include "serac/numerics/functional/domain.hpp" #include "serac/numerics/functional/element_restriction.hpp" +#include "serac/infrastructure/application_manager.hpp" +#include "serac/mesh_utils/mesh_utils.hpp" using namespace serac; @@ -26,18 +28,19 @@ TEST(patch_test_meshes, triangle_domains) { int p = 2; int dim = 2; - mfem::Mesh mesh(SERAC_REPO_DIR "/data/meshes/patch2D_tris.mesh"); + mfem::Mesh bmesh(SERAC_REPO_DIR "/data/meshes/patch2D_tris.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(bmesh)); auto H1_fec = std::make_unique(p, dim); - auto H1_fes = std::make_unique(&mesh, H1_fec.get()); + auto H1_fes = std::make_unique(mesh.get(), H1_fec.get()); auto Hcurl_fec = std::make_unique(p, dim); - auto Hcurl_fes = std::make_unique(&mesh, Hcurl_fec.get()); + auto Hcurl_fes = std::make_unique(mesh.get(), Hcurl_fec.get()); auto L2_fec = std::make_unique(p, dim, mfem::BasisType::GaussLobatto); - auto L2_fes = std::make_unique(&mesh, L2_fec.get()); + auto L2_fes = std::make_unique(mesh.get(), L2_fec.get()); - Domain whole = EntireDomain(mesh); + Domain whole = EntireDomain(*mesh); EXPECT_EQ(whole.mfem_edge_ids_.size(), 0); EXPECT_EQ(whole.mfem_tri_ids_.size(), 4); EXPECT_EQ(whole.mfem_quad_ids_.size(), 0); @@ -55,7 +58,7 @@ TEST(patch_test_meshes, triangle_domains) EXPECT_EQ(L2_BER.ESize(), 4 * 6); } - Domain boundary = EntireBoundary(mesh); + Domain boundary = EntireBoundary(*mesh); EXPECT_EQ(boundary.mfem_edge_ids_.size(), 4); EXPECT_EQ(boundary.mfem_tri_ids_.size(), 0); EXPECT_EQ(boundary.mfem_quad_ids_.size(), 0); @@ -73,7 +76,7 @@ TEST(patch_test_meshes, triangle_domains) EXPECT_EQ(L2_BER.ESize(), 4 * 3); } - Domain interior = InteriorFaces(mesh); + Domain interior = InteriorFaces(*mesh); EXPECT_EQ(interior.mfem_edge_ids_.size(), 4); EXPECT_EQ(interior.mfem_tri_ids_.size(), 0); EXPECT_EQ(interior.mfem_quad_ids_.size(), 0); @@ -96,17 +99,18 @@ TEST(patch_test_meshes, quadrilateral_domains) { int p = 2; int dim = 2; - mfem::Mesh mesh(SERAC_REPO_DIR "/data/meshes/patch2D_quads.mesh"); + mfem::Mesh bmesh(SERAC_REPO_DIR "/data/meshes/patch2D_quads.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(bmesh)); auto H1_fec = std::make_unique(p, dim); auto Hcurl_fec = std::make_unique(p, dim); auto L2_fec = std::make_unique(p, dim, mfem::BasisType::GaussLobatto); - auto H1_fes = std::make_unique(&mesh, H1_fec.get()); - auto Hcurl_fes = std::make_unique(&mesh, Hcurl_fec.get()); - auto L2_fes = std::make_unique(&mesh, L2_fec.get()); + auto H1_fes = std::make_unique(mesh.get(), H1_fec.get()); + auto Hcurl_fes = std::make_unique(mesh.get(), Hcurl_fec.get()); + auto L2_fes = std::make_unique(mesh.get(), L2_fec.get()); - Domain whole = EntireDomain(mesh); + Domain whole = EntireDomain(*mesh); EXPECT_EQ(whole.mfem_edge_ids_.size(), 0); EXPECT_EQ(whole.mfem_tri_ids_.size(), 0); EXPECT_EQ(whole.mfem_quad_ids_.size(), 5); @@ -124,7 +128,7 @@ TEST(patch_test_meshes, quadrilateral_domains) EXPECT_EQ(L2_BER.ESize(), 5 * 9); } - Domain boundary = EntireBoundary(mesh); + Domain boundary = EntireBoundary(*mesh); EXPECT_EQ(boundary.mfem_edge_ids_.size(), 4); EXPECT_EQ(boundary.mfem_tri_ids_.size(), 0); EXPECT_EQ(boundary.mfem_quad_ids_.size(), 0); @@ -142,7 +146,7 @@ TEST(patch_test_meshes, quadrilateral_domains) EXPECT_EQ(L2_BER.ESize(), 4 * 3); } - Domain interior = InteriorFaces(mesh); + Domain interior = InteriorFaces(*mesh); EXPECT_EQ(interior.mfem_edge_ids_.size(), 8); EXPECT_EQ(interior.mfem_tri_ids_.size(), 0); EXPECT_EQ(interior.mfem_quad_ids_.size(), 0); @@ -165,17 +169,18 @@ TEST(patch_test_meshes, triangle_and_quadrilateral_domains) { int p = 2; int dim = 2; - mfem::Mesh mesh(SERAC_REPO_DIR "/data/meshes/patch2D_tris_and_quads.mesh"); + mfem::Mesh bmesh(SERAC_REPO_DIR "/data/meshes/patch2D_tris_and_quads.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(bmesh)); auto H1_fec = std::make_unique(p, dim); auto Hcurl_fec = std::make_unique(p, dim); auto L2_fec = std::make_unique(p, dim, mfem::BasisType::GaussLobatto); - auto H1_fes = std::make_unique(&mesh, H1_fec.get()); - auto Hcurl_fes = std::make_unique(&mesh, Hcurl_fec.get()); - auto L2_fes = std::make_unique(&mesh, L2_fec.get()); + auto H1_fes = std::make_unique(mesh.get(), H1_fec.get()); + auto Hcurl_fes = std::make_unique(mesh.get(), Hcurl_fec.get()); + auto L2_fes = std::make_unique(mesh.get(), L2_fec.get()); - Domain whole = EntireDomain(mesh); + Domain whole = EntireDomain(*mesh); EXPECT_EQ(whole.mfem_edge_ids_.size(), 0); EXPECT_EQ(whole.mfem_tri_ids_.size(), 2); EXPECT_EQ(whole.mfem_quad_ids_.size(), 4); @@ -193,7 +198,7 @@ TEST(patch_test_meshes, triangle_and_quadrilateral_domains) EXPECT_EQ(L2_BER.ESize(), 2 * 6 + 4 * 9); } - Domain boundary = EntireBoundary(mesh); + Domain boundary = EntireBoundary(*mesh); EXPECT_EQ(boundary.mfem_edge_ids_.size(), 4); EXPECT_EQ(boundary.mfem_tri_ids_.size(), 0); EXPECT_EQ(boundary.mfem_quad_ids_.size(), 0); @@ -211,7 +216,7 @@ TEST(patch_test_meshes, triangle_and_quadrilateral_domains) EXPECT_EQ(L2_BER.ESize(), 4 * 3); } - Domain interior = InteriorFaces(mesh); + Domain interior = InteriorFaces(*mesh); EXPECT_EQ(interior.mfem_edge_ids_.size(), 9); EXPECT_EQ(interior.mfem_tri_ids_.size(), 0); EXPECT_EQ(interior.mfem_quad_ids_.size(), 0); @@ -234,17 +239,18 @@ TEST(patch_test_meshes, tetrahedron_domains) { int p = 2; int dim = 3; - mfem::Mesh mesh(SERAC_REPO_DIR "/data/meshes/patch3D_tets.mesh"); + mfem::Mesh bmesh(SERAC_REPO_DIR "/data/meshes/patch3D_tets.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(bmesh)); auto H1_fec = std::make_unique(p, dim); auto Hcurl_fec = std::make_unique(p, dim); auto L2_fec = std::make_unique(p, dim, mfem::BasisType::GaussLobatto); - auto H1_fes = std::make_unique(&mesh, H1_fec.get()); - auto Hcurl_fes = std::make_unique(&mesh, Hcurl_fec.get()); - auto L2_fes = std::make_unique(&mesh, L2_fec.get()); + auto H1_fes = std::make_unique(mesh.get(), H1_fec.get()); + auto Hcurl_fes = std::make_unique(mesh.get(), Hcurl_fec.get()); + auto L2_fes = std::make_unique(mesh.get(), L2_fec.get()); - Domain whole = EntireDomain(mesh); + Domain whole = EntireDomain(*mesh); EXPECT_EQ(whole.mfem_edge_ids_.size(), 0); EXPECT_EQ(whole.mfem_tri_ids_.size(), 0); EXPECT_EQ(whole.mfem_quad_ids_.size(), 0); @@ -262,7 +268,7 @@ TEST(patch_test_meshes, tetrahedron_domains) EXPECT_EQ(L2_BER.ESize(), 12 * 10); } - Domain boundary = EntireBoundary(mesh); + Domain boundary = EntireBoundary(*mesh); EXPECT_EQ(boundary.mfem_edge_ids_.size(), 0); EXPECT_EQ(boundary.mfem_tri_ids_.size(), 12); EXPECT_EQ(boundary.mfem_quad_ids_.size(), 0); @@ -280,7 +286,7 @@ TEST(patch_test_meshes, tetrahedron_domains) EXPECT_EQ(L2_BER.ESize(), 12 * 6); } - Domain interior = InteriorFaces(mesh); + Domain interior = InteriorFaces(*mesh); EXPECT_EQ(interior.mfem_edge_ids_.size(), 0); EXPECT_EQ(interior.mfem_tri_ids_.size(), 18); EXPECT_EQ(interior.mfem_quad_ids_.size(), 0); @@ -303,17 +309,18 @@ TEST(patch_test_meshes, hexahedron_domains) { int p = 2; int dim = 3; - mfem::Mesh mesh(SERAC_REPO_DIR "/data/meshes/patch3D_hexes.mesh"); + mfem::Mesh bmesh(SERAC_REPO_DIR "/data/meshes/patch3D_hexes.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(bmesh)); auto H1_fec = std::make_unique(p, dim); auto Hcurl_fec = std::make_unique(p, dim); auto L2_fec = std::make_unique(p, dim, mfem::BasisType::GaussLobatto); - auto H1_fes = std::make_unique(&mesh, H1_fec.get()); - auto Hcurl_fes = std::make_unique(&mesh, Hcurl_fec.get()); - auto L2_fes = std::make_unique(&mesh, L2_fec.get()); + auto H1_fes = std::make_unique(mesh.get(), H1_fec.get()); + auto Hcurl_fes = std::make_unique(mesh.get(), Hcurl_fec.get()); + auto L2_fes = std::make_unique(mesh.get(), L2_fec.get()); - Domain whole = EntireDomain(mesh); + Domain whole = EntireDomain(*mesh); EXPECT_EQ(whole.mfem_edge_ids_.size(), 0); EXPECT_EQ(whole.mfem_tri_ids_.size(), 0); EXPECT_EQ(whole.mfem_quad_ids_.size(), 0); @@ -331,7 +338,7 @@ TEST(patch_test_meshes, hexahedron_domains) EXPECT_EQ(L2_BER.ESize(), 7 * 27); } - Domain boundary = EntireBoundary(mesh); + Domain boundary = EntireBoundary(*mesh); EXPECT_EQ(boundary.mfem_edge_ids_.size(), 0); EXPECT_EQ(boundary.mfem_tri_ids_.size(), 0); EXPECT_EQ(boundary.mfem_quad_ids_.size(), 6); @@ -349,7 +356,7 @@ TEST(patch_test_meshes, hexahedron_domains) EXPECT_EQ(L2_BER.ESize(), 6 * 9); } - Domain interior = InteriorFaces(mesh); + Domain interior = InteriorFaces(*mesh); EXPECT_EQ(interior.mfem_edge_ids_.size(), 0); EXPECT_EQ(interior.mfem_tri_ids_.size(), 0); EXPECT_EQ(interior.mfem_quad_ids_.size(), 18); @@ -367,3 +374,10 @@ TEST(patch_test_meshes, hexahedron_domains) EXPECT_EQ(L2_BER.ESize(), 18 * (9 * 2)); } } + +int main(int argc, char* argv[]) +{ + ::testing::InitGoogleTest(&argc, argv); + serac::ApplicationManager applicationManager(argc, argv); + return RUN_ALL_TESTS(); +} diff --git a/src/serac/numerics/functional/tests/functional_basic_dg.cpp b/src/serac/numerics/functional/tests/functional_basic_dg.cpp index 6705c8713d..6c307e65eb 100644 --- a/src/serac/numerics/functional/tests/functional_basic_dg.cpp +++ b/src/serac/numerics/functional/tests/functional_basic_dg.cpp @@ -33,7 +33,7 @@ void L2_test(std::string meshfile) // int k = 0; // while (k == 0); - auto mesh = mesh::refineAndDistribute(buildMeshFromFile(meshfile), 0); + auto mesh = mesh::refineAndDistribute(buildMeshFromFile(meshfile), 1); auto fec = mfem::L2_FECollection(p, dim, mfem::BasisType::GaussLobatto); mfem::ParFiniteElementSpace fespace(mesh.get(), &fec, dim, serac::ordering); @@ -69,8 +69,8 @@ void L2_test(std::string meshfile) double t = 0.0; - auto value = residual(t, U); - // check_gradient(residual, t, U); + // auto value = residual(t, U); + check_gradient(residual, t, U); } TEST(basic, L2_test_tris_and_quads_linear) { L2_test<2, 1>(SERAC_REPO_DIR "/data/meshes/patch2D_tris_and_quads.mesh"); } @@ -195,11 +195,21 @@ TEST(basic, L2_mixed_scalar_test_tris_and_quads_linear) L2_scalar_valued_test<2, 1>(SERAC_REPO_DIR "/data/meshes/patch2D_tris_and_quads.mesh"); } +TEST(basic, L2_mixed_scalar_test_tris_and_quads_quadratic) +{ + L2_scalar_valued_test<2, 2>(SERAC_REPO_DIR "/data/meshes/patch2D_tris_and_quads.mesh"); +} + TEST(basic, L2_mixed_scalar_test_tets_and_hexes_linear) { L2_scalar_valued_test<3, 1>(SERAC_REPO_DIR "/data/meshes/patch3D_tets_and_hexes.mesh"); } +TEST(basic, L2_mixed_scalar_test_tets_and_hexes_quadratic) +{ + L2_scalar_valued_test<3, 2>(SERAC_REPO_DIR "/data/meshes/patch3D_tets_and_hexes.mesh"); +} + int main(int argc, char* argv[]) { ::testing::InitGoogleTest(&argc, argv); diff --git a/src/serac/numerics/functional/tests/geometric_factors_tests.cpp b/src/serac/numerics/functional/tests/geometric_factors_tests.cpp index 3a7a06f171..baea47021c 100644 --- a/src/serac/numerics/functional/tests/geometric_factors_tests.cpp +++ b/src/serac/numerics/functional/tests/geometric_factors_tests.cpp @@ -8,6 +8,7 @@ #include "serac/numerics/functional/geometric_factors.hpp" #include "serac/infrastructure/application_manager.hpp" +#include "serac/mesh_utils/mesh_utils.hpp" using namespace serac; @@ -30,11 +31,12 @@ mfem::Mesh import_mesh(std::string meshfile) TEST(geometric_factors, with_2D_domains) { - auto mesh = import_mesh("patch2D_tris_and_quads.mesh"); + auto bmesh = import_mesh("patch2D_tris_and_quads.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(bmesh)); // `d` will consist of one tri and one quad Domain d = Domain::ofElements( - mesh, std::function([](std::vector vertices, int /* attr */) { return average(vertices)[0] < 0.45; })); + *mesh, std::function([](std::vector vertices, int /* attr */) { return average(vertices)[0] < 0.45; })); int q = 2; @@ -59,12 +61,13 @@ TEST(geometric_factors, with_2D_domains) TEST(geometric_factors, with_3D_domains) { - auto mesh = import_mesh("patch3D_tets_and_hexes.mesh"); + auto bmesh = import_mesh("patch3D_tets_and_hexes.mesh"); + auto mesh = serac::mesh::refineAndDistribute(std::move(bmesh)); // `d` will consist of 6 tets and 1 hex - Domain d = Domain::ofElements(mesh, std::function([](std::vector vertices, int /*bdr_attr*/) { - return average(vertices)[1] < 0.75; // y coordinate of face center - })); + Domain d = Domain::ofElements(*mesh, std::function([](std::vector vertices, int /*bdr_attr*/) { + return average(vertices)[1] < 0.75; // y coordinate of face center + })); int q = 2; diff --git a/src/serac/numerics/functional/typedefs.hpp b/src/serac/numerics/functional/typedefs.hpp index dad72993f0..74eb0d4f9d 100644 --- a/src/serac/numerics/functional/typedefs.hpp +++ b/src/serac/numerics/functional/typedefs.hpp @@ -13,8 +13,8 @@ namespace serac { // interface change like that, so these typedefs mark the parts that would need to eventually change /// @cond -using mesh_t = mfem::Mesh; -using fes_t = mfem::FiniteElementSpace; +using mesh_t = mfem::ParMesh; +using fes_t = mfem::ParFiniteElementSpace; /// @endcond } // namespace serac diff --git a/src/serac/physics/state/finite_element_state.cpp b/src/serac/physics/state/finite_element_state.cpp index 0bb4d10d1c..29e81d7bda 100644 --- a/src/serac/physics/state/finite_element_state.cpp +++ b/src/serac/physics/state/finite_element_state.cpp @@ -75,7 +75,7 @@ void FiniteElementState::projectOnBoundary(mfem::VectorCoefficient& coef, const void FiniteElementState::project(mfem::Coefficient& coef, const Domain& domain) { - mfem::Array uniq_dof_ids = domain.dof_list(gridFunction().FESpace()); + mfem::Array uniq_dof_ids = domain.dof_list(gridFunction().ParFESpace()); mfem::ParGridFunction& grid_function = gridFunction(); grid_function.ProjectCoefficient(coef, uniq_dof_ids); setFromGridFunction(grid_function); @@ -83,7 +83,7 @@ void FiniteElementState::project(mfem::Coefficient& coef, const Domain& domain) void FiniteElementState::project(mfem::VectorCoefficient& coef, const Domain& domain) { - mfem::Array uniq_dof_ids = domain.dof_list(gridFunction().FESpace()); + mfem::Array uniq_dof_ids = domain.dof_list(gridFunction().ParFESpace()); mfem::ParGridFunction& grid_function = gridFunction(); grid_function.ProjectCoefficient(coef, uniq_dof_ids); setFromGridFunction(grid_function);