From 77b2be78062488f0c6976345e6afe447a681287f Mon Sep 17 00:00:00 2001 From: Gary Hu Date: Thu, 2 Oct 2025 16:40:09 -0500 Subject: [PATCH] Allow specifying element range in project_vector --- include/systems/system.h | 22 ++++++---- src/systems/system_projection.C | 73 +++++++++++++++++++++------------ 2 files changed, 61 insertions(+), 34 deletions(-) diff --git a/include/systems/system.h b/include/systems/system.h index 05fd94833ff..1e5debcbf01 100644 --- a/include/systems/system.h +++ b/include/systems/system.h @@ -40,6 +40,7 @@ #include #include #include +#include // This define may be useful in --disable-optional builds when it is // possible that libmesh will not have any solvers available. @@ -561,7 +562,8 @@ class System : public ReferenceCountedObject, void project_vector (NumericVector & new_vector, FunctionBase * f, FunctionBase * g = nullptr, - int is_adjoint = -1) const; + int is_adjoint = -1, + std::optional active_local_range = std::nullopt) const; /** * Projects arbitrary functions onto a vector of degree of freedom @@ -577,7 +579,8 @@ class System : public ReferenceCountedObject, void project_vector (NumericVector & new_vector, FEMFunctionBase * f, FEMFunctionBase * g = nullptr, - int is_adjoint = -1) const; + int is_adjoint = -1, + std::optional active_local_range = std::nullopt) const; /** * Projects arbitrary functions onto a vector of degree of freedom @@ -594,7 +597,8 @@ class System : public ReferenceCountedObject, GradientFunctionPointer gptr, const Parameters & parameters, NumericVector & new_vector, - int is_adjoint = -1) const; + int is_adjoint = -1, + std::optional active_local_range = std::nullopt) const; /** * Projects arbitrary boundary functions onto a vector of degree of @@ -651,7 +655,8 @@ class System : public ReferenceCountedObject, NumericVector & new_vector, FunctionBase * f, FunctionBase * g = nullptr, - int is_adjoint = -1) const; + int is_adjoint = -1, + std::optional active_local_range = std::nullopt) const; /** * Projects arbitrary boundary functions onto a vector of degree of @@ -674,7 +679,8 @@ class System : public ReferenceCountedObject, GradientFunctionPointer gptr, const Parameters & parameters, NumericVector & new_vector, - int is_adjoint = -1) const; + int is_adjoint = -1, + std::optional active_local_range = std::nullopt) const; /** * \returns The system number. @@ -2005,7 +2011,8 @@ class System : public ReferenceCountedObject, * primal constraints if is_adjoint is non-negative. */ void project_vector (NumericVector &, - int is_adjoint = -1) const; + int is_adjoint = -1, + std::optional active_local_range = std::nullopt) const; /** * Projects the vector defined on the old mesh onto the @@ -2017,7 +2024,8 @@ class System : public ReferenceCountedObject, */ void project_vector (const NumericVector &, NumericVector &, - int is_adjoint = -1) const; + int is_adjoint = -1, + std::optional active_local_range = std::nullopt) const; /* * If we have e.g. a element space constrained by spline values, we diff --git a/src/systems/system_projection.C b/src/systems/system_projection.C index 3970fd4a446..9b13503741f 100644 --- a/src/systems/system_projection.C +++ b/src/systems/system_projection.C @@ -245,7 +245,8 @@ public: // ------------------------------------------------------------ // System implementation void System::project_vector (NumericVector & vector, - int is_adjoint) const + int is_adjoint, + std::optional active_local_range) const { // Create a copy of the vector, which currently // contains the old data. @@ -253,7 +254,7 @@ void System::project_vector (NumericVector & vector, old_vector (vector.clone()); // Project the old vector to the new vector - this->project_vector (*old_vector, vector, is_adjoint); + this->project_vector (*old_vector, vector, is_adjoint, active_local_range); } @@ -264,7 +265,8 @@ void System::project_vector (NumericVector & vector, */ void System::project_vector (const NumericVector & old_v, NumericVector & new_v, - int is_adjoint) const + int is_adjoint, + std::optional active_local_range) const { LOG_SCOPE ("project_vector(old,new)", "System"); @@ -285,9 +287,12 @@ void System::project_vector (const NumericVector & old_v, std::unique_ptr> local_old_vector_built; const NumericVector * old_vector_ptr = nullptr; - ConstElemRange active_local_elem_range - (this->get_mesh().active_local_elements_begin(), - this->get_mesh().active_local_elements_end()); + if (!active_local_range) + { + active_local_range.emplace + (this->get_mesh().active_local_elements_begin(), + this->get_mesh().active_local_elements_end()); + } // If the old vector was uniprocessor, make the new // vector uniprocessor @@ -304,7 +309,7 @@ void System::project_vector (const NumericVector & old_v, { // Build a send list for efficient localization BuildProjectionList projection_list(*this); - Threads::parallel_reduce (active_local_elem_range, + Threads::parallel_reduce (active_local_range.value(), projection_list); // Create a sorted, unique send_list @@ -328,7 +333,7 @@ void System::project_vector (const NumericVector & old_v, { // Build a send list for efficient localization BuildProjectionList projection_list(*this); - Threads::parallel_reduce (active_local_elem_range, + Threads::parallel_reduce (active_local_range.value(), projection_list); // Create a sorted, unique send_list @@ -389,7 +394,7 @@ void System::project_vector (const NumericVector & old_v, g(*this, old_vector, ®ular_vars); FEMProjector projector(*this, f, &g, setter, regular_vars); - projector.project(active_local_elem_range); + projector.project(active_local_range.value()); } if (!vector_vars.empty()) @@ -403,7 +408,7 @@ void System::project_vector (const NumericVector & old_v, OldSolutionValue g_vector(*this, old_vector, &vector_vars); FEMVectorProjector vector_projector(*this, f_vector, &g_vector, setter, vector_vars); - vector_projector.project(active_local_elem_range); + vector_projector.project(active_local_range.value()); } // Copy the SCALAR dofs from old_vector to new_vector @@ -1068,11 +1073,12 @@ void System::project_vector (ValueFunctionPointer fptr, GradientFunctionPointer gptr, const Parameters & function_parameters, NumericVector & new_vector, - int is_adjoint) const + int is_adjoint, + std::optional active_local_range) const { WrappedFunction f(*this, fptr, &function_parameters); WrappedFunction g(*this, gptr, &function_parameters); - this->project_vector(new_vector, &f, &g, is_adjoint); + this->project_vector(new_vector, &f, &g, is_adjoint, active_local_range); } /** @@ -1082,7 +1088,8 @@ void System::project_vector (ValueFunctionPointer fptr, void System::project_vector (NumericVector & new_vector, FunctionBase * f, FunctionBase * g, - int is_adjoint) const + int is_adjoint, + std::optional active_local_range) const { LOG_SCOPE ("project_vector(FunctionBase)", "System"); @@ -1094,10 +1101,10 @@ void System::project_vector (NumericVector & new_vector, { WrappedFunctor g_fem(*g); - this->project_vector(new_vector, &f_fem, &g_fem, is_adjoint); + this->project_vector(new_vector, &f_fem, &g_fem, is_adjoint, active_local_range); } else - this->project_vector(new_vector, &f_fem, nullptr, is_adjoint); + this->project_vector(new_vector, &f_fem, nullptr, is_adjoint, active_local_range); } @@ -1108,15 +1115,19 @@ void System::project_vector (NumericVector & new_vector, void System::project_vector (NumericVector & new_vector, FEMFunctionBase * f, FEMFunctionBase * g, - int is_adjoint) const + int is_adjoint, + std::optional active_local_range) const { LOG_SCOPE ("project_fem_vector()", "System"); libmesh_assert (f); - ConstElemRange active_local_range - (this->get_mesh().active_local_elements_begin(), - this->get_mesh().active_local_elements_end() ); + if (!active_local_range) + { + active_local_range.emplace + (this->get_mesh().active_local_elements_begin(), + this->get_mesh().active_local_elements_end()); + } VectorSetAction setter(new_vector); @@ -1137,12 +1148,12 @@ void System::project_vector (NumericVector & new_vector, FEMFunctionWrapper gw(*g); FEMProjector projector(*this, fw, &gw, setter, vars); - projector.project(active_local_range); + projector.project(active_local_range.value()); } else { FEMProjector projector(*this, fw, nullptr, setter, vars); - projector.project(active_local_range); + projector.project(active_local_range.value()); } // Also, load values into the SCALAR dofs @@ -1197,7 +1208,7 @@ void System::project_vector (NumericVector & new_vector, { // Look for a spline node: a NodeElem with a rational variable // on it. - for (auto & elem : active_local_range) + for (auto & elem : active_local_range.value()) if (elem->type() == NODEELEM) for (auto rational_var : rational_vars) if (rational_var->active_on_subdomain(elem->subdomain_id())) @@ -1274,12 +1285,13 @@ void System::boundary_project_vector (const std::set & b, GradientFunctionPointer gptr, const Parameters & function_parameters, NumericVector & new_vector, - int is_adjoint) const + int is_adjoint, + std::optional active_local_range) const { WrappedFunction f(*this, fptr, &function_parameters); WrappedFunction g(*this, gptr, &function_parameters); this->boundary_project_vector(b, variables, new_vector, &f, &g, - is_adjoint); + is_adjoint, active_local_range); } /** @@ -1291,13 +1303,20 @@ void System::boundary_project_vector (const std::set & b, NumericVector & new_vector, FunctionBase * f, FunctionBase * g, - int is_adjoint) const + int is_adjoint, + std::optional active_local_range) const { LOG_SCOPE ("boundary_project_vector()", "System"); + if (!active_local_range) + { + active_local_range.emplace + (this->get_mesh().active_local_elements_begin(), + this->get_mesh().active_local_elements_end()); + } + Threads::parallel_for - (ConstElemRange (this->get_mesh().active_local_elements_begin(), - this->get_mesh().active_local_elements_end() ), + (active_local_range.value(), BoundaryProjectSolution(b, variables, *this, f, g, this->get_equation_systems().parameters, new_vector)