Skip to content

Commit 77b2be7

Browse files
hugary1995hugary1995
committed
Allow specifying element range in project_vector
1 parent 2fd6181 commit 77b2be7

File tree

2 files changed

+61
-34
lines changed

2 files changed

+61
-34
lines changed

‎include/systems/system.h‎

Lines changed: 15 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -40,6 +40,7 @@
4040
#include <string_view>
4141
#include <vector>
4242
#include <memory>
43+
#include <optional>
4344

4445
// This define may be useful in --disable-optional builds when it is
4546
// possible that libmesh will not have any solvers available.
@@ -561,7 +562,8 @@ class System : public ReferenceCountedObject<System>,
561562
void project_vector (NumericVector<Number> & new_vector,
562563
FunctionBase<Number> * f,
563564
FunctionBase<Gradient> * g = nullptr,
564-
int is_adjoint = -1) const;
565+
int is_adjoint = -1,
566+
std::optional<ConstElemRange> active_local_range = std::nullopt) const;
565567

566568
/**
567569
* Projects arbitrary functions onto a vector of degree of freedom
@@ -577,7 +579,8 @@ class System : public ReferenceCountedObject<System>,
577579
void project_vector (NumericVector<Number> & new_vector,
578580
FEMFunctionBase<Number> * f,
579581
FEMFunctionBase<Gradient> * g = nullptr,
580-
int is_adjoint = -1) const;
582+
int is_adjoint = -1,
583+
std::optional<ConstElemRange> active_local_range = std::nullopt) const;
581584

582585
/**
583586
* Projects arbitrary functions onto a vector of degree of freedom
@@ -594,7 +597,8 @@ class System : public ReferenceCountedObject<System>,
594597
GradientFunctionPointer gptr,
595598
const Parameters & parameters,
596599
NumericVector<Number> & new_vector,
597-
int is_adjoint = -1) const;
600+
int is_adjoint = -1,
601+
std::optional<ConstElemRange> active_local_range = std::nullopt) const;
598602

599603
/**
600604
* Projects arbitrary boundary functions onto a vector of degree of
@@ -651,7 +655,8 @@ class System : public ReferenceCountedObject<System>,
651655
NumericVector<Number> & new_vector,
652656
FunctionBase<Number> * f,
653657
FunctionBase<Gradient> * g = nullptr,
654-
int is_adjoint = -1) const;
658+
int is_adjoint = -1,
659+
std::optional<ConstElemRange> active_local_range = std::nullopt) const;
655660

656661
/**
657662
* Projects arbitrary boundary functions onto a vector of degree of
@@ -674,7 +679,8 @@ class System : public ReferenceCountedObject<System>,
674679
GradientFunctionPointer gptr,
675680
const Parameters & parameters,
676681
NumericVector<Number> & new_vector,
677-
int is_adjoint = -1) const;
682+
int is_adjoint = -1,
683+
std::optional<ConstElemRange> active_local_range = std::nullopt) const;
678684

679685
/**
680686
* \returns The system number.
@@ -2005,7 +2011,8 @@ class System : public ReferenceCountedObject<System>,
20052011
* primal constraints if is_adjoint is non-negative.
20062012
*/
20072013
void project_vector (NumericVector<Number> &,
2008-
int is_adjoint = -1) const;
2014+
int is_adjoint = -1,
2015+
std::optional<ConstElemRange> active_local_range = std::nullopt) const;
20092016

20102017
/**
20112018
* Projects the vector defined on the old mesh onto the
@@ -2017,7 +2024,8 @@ class System : public ReferenceCountedObject<System>,
20172024
*/
20182025
void project_vector (const NumericVector<Number> &,
20192026
NumericVector<Number> &,
2020-
int is_adjoint = -1) const;
2027+
int is_adjoint = -1,
2028+
std::optional<ConstElemRange> active_local_range = std::nullopt) const;
20212029

20222030
/*
20232031
* If we have e.g. a element space constrained by spline values, we

‎src/systems/system_projection.C‎

Lines changed: 46 additions & 27 deletions
Original file line numberDiff line numberDiff line change
@@ -245,15 +245,16 @@ public:
245245
// ------------------------------------------------------------
246246
// System implementation
247247
void System::project_vector (NumericVector<Number> & vector,
248-
int is_adjoint) const
248+
int is_adjoint,
249+
std::optional<ConstElemRange> active_local_range) const
249250
{
250251
// Create a copy of the vector, which currently
251252
// contains the old data.
252253
std::unique_ptr<NumericVector<Number>>
253254
old_vector (vector.clone());
254255

255256
// Project the old vector to the new vector
256-
this->project_vector (*old_vector, vector, is_adjoint);
257+
this->project_vector (*old_vector, vector, is_adjoint, active_local_range);
257258
}
258259

259260

@@ -264,7 +265,8 @@ void System::project_vector (NumericVector<Number> & vector,
264265
*/
265266
void System::project_vector (const NumericVector<Number> & old_v,
266267
NumericVector<Number> & new_v,
267-
int is_adjoint) const
268+
int is_adjoint,
269+
std::optional<ConstElemRange> active_local_range) const
268270
{
269271
LOG_SCOPE ("project_vector(old,new)", "System");
270272

@@ -285,9 +287,12 @@ void System::project_vector (const NumericVector<Number> & old_v,
285287
std::unique_ptr<NumericVector<Number>> local_old_vector_built;
286288
const NumericVector<Number> * old_vector_ptr = nullptr;
287289

288-
ConstElemRange active_local_elem_range
289-
(this->get_mesh().active_local_elements_begin(),
290-
this->get_mesh().active_local_elements_end());
290+
if (!active_local_range)
291+
{
292+
active_local_range.emplace
293+
(this->get_mesh().active_local_elements_begin(),
294+
this->get_mesh().active_local_elements_end());
295+
}
291296

292297
// If the old vector was uniprocessor, make the new
293298
// vector uniprocessor
@@ -304,7 +309,7 @@ void System::project_vector (const NumericVector<Number> & old_v,
304309
{
305310
// Build a send list for efficient localization
306311
BuildProjectionList projection_list(*this);
307-
Threads::parallel_reduce (active_local_elem_range,
312+
Threads::parallel_reduce (active_local_range.value(),
308313
projection_list);
309314

310315
// Create a sorted, unique send_list
@@ -328,7 +333,7 @@ void System::project_vector (const NumericVector<Number> & old_v,
328333
{
329334
// Build a send list for efficient localization
330335
BuildProjectionList projection_list(*this);
331-
Threads::parallel_reduce (active_local_elem_range,
336+
Threads::parallel_reduce (active_local_range.value(),
332337
projection_list);
333338

334339
// Create a sorted, unique send_list
@@ -389,7 +394,7 @@ void System::project_vector (const NumericVector<Number> & old_v,
389394
g(*this, old_vector, &regular_vars);
390395

391396
FEMProjector projector(*this, f, &g, setter, regular_vars);
392-
projector.project(active_local_elem_range);
397+
projector.project(active_local_range.value());
393398
}
394399

395400
if (!vector_vars.empty())
@@ -403,7 +408,7 @@ void System::project_vector (const NumericVector<Number> & old_v,
403408
OldSolutionValue<Tensor, &FEMContext::point_gradient> g_vector(*this, old_vector, &vector_vars);
404409

405410
FEMVectorProjector vector_projector(*this, f_vector, &g_vector, setter, vector_vars);
406-
vector_projector.project(active_local_elem_range);
411+
vector_projector.project(active_local_range.value());
407412
}
408413

409414
// Copy the SCALAR dofs from old_vector to new_vector
@@ -1068,11 +1073,12 @@ void System::project_vector (ValueFunctionPointer fptr,
10681073
GradientFunctionPointer gptr,
10691074
const Parameters & function_parameters,
10701075
NumericVector<Number> & new_vector,
1071-
int is_adjoint) const
1076+
int is_adjoint,
1077+
std::optional<ConstElemRange> active_local_range) const
10721078
{
10731079
WrappedFunction<Number> f(*this, fptr, &function_parameters);
10741080
WrappedFunction<Gradient> g(*this, gptr, &function_parameters);
1075-
this->project_vector(new_vector, &f, &g, is_adjoint);
1081+
this->project_vector(new_vector, &f, &g, is_adjoint, active_local_range);
10761082
}
10771083

10781084
/**
@@ -1082,7 +1088,8 @@ void System::project_vector (ValueFunctionPointer fptr,
10821088
void System::project_vector (NumericVector<Number> & new_vector,
10831089
FunctionBase<Number> * f,
10841090
FunctionBase<Gradient> * g,
1085-
int is_adjoint) const
1091+
int is_adjoint,
1092+
std::optional<ConstElemRange> active_local_range) const
10861093
{
10871094
LOG_SCOPE ("project_vector(FunctionBase)", "System");
10881095

@@ -1094,10 +1101,10 @@ void System::project_vector (NumericVector<Number> & new_vector,
10941101
{
10951102
WrappedFunctor<Gradient> g_fem(*g);
10961103

1097-
this->project_vector(new_vector, &f_fem, &g_fem, is_adjoint);
1104+
this->project_vector(new_vector, &f_fem, &g_fem, is_adjoint, active_local_range);
10981105
}
10991106
else
1100-
this->project_vector(new_vector, &f_fem, nullptr, is_adjoint);
1107+
this->project_vector(new_vector, &f_fem, nullptr, is_adjoint, active_local_range);
11011108
}
11021109

11031110

@@ -1108,15 +1115,19 @@ void System::project_vector (NumericVector<Number> & new_vector,
11081115
void System::project_vector (NumericVector<Number> & new_vector,
11091116
FEMFunctionBase<Number> * f,
11101117
FEMFunctionBase<Gradient> * g,
1111-
int is_adjoint) const
1118+
int is_adjoint,
1119+
std::optional<ConstElemRange> active_local_range) const
11121120
{
11131121
LOG_SCOPE ("project_fem_vector()", "System");
11141122

11151123
libmesh_assert (f);
11161124

1117-
ConstElemRange active_local_range
1118-
(this->get_mesh().active_local_elements_begin(),
1119-
this->get_mesh().active_local_elements_end() );
1125+
if (!active_local_range)
1126+
{
1127+
active_local_range.emplace
1128+
(this->get_mesh().active_local_elements_begin(),
1129+
this->get_mesh().active_local_elements_end());
1130+
}
11201131

11211132
VectorSetAction<Number> setter(new_vector);
11221133

@@ -1137,12 +1148,12 @@ void System::project_vector (NumericVector<Number> & new_vector,
11371148
FEMFunctionWrapper<Gradient> gw(*g);
11381149

11391150
FEMProjector projector(*this, fw, &gw, setter, vars);
1140-
projector.project(active_local_range);
1151+
projector.project(active_local_range.value());
11411152
}
11421153
else
11431154
{
11441155
FEMProjector projector(*this, fw, nullptr, setter, vars);
1145-
projector.project(active_local_range);
1156+
projector.project(active_local_range.value());
11461157
}
11471158

11481159
// Also, load values into the SCALAR dofs
@@ -1197,7 +1208,7 @@ void System::project_vector (NumericVector<Number> & new_vector,
11971208
{
11981209
// Look for a spline node: a NodeElem with a rational variable
11991210
// on it.
1200-
for (auto & elem : active_local_range)
1211+
for (auto & elem : active_local_range.value())
12011212
if (elem->type() == NODEELEM)
12021213
for (auto rational_var : rational_vars)
12031214
if (rational_var->active_on_subdomain(elem->subdomain_id()))
@@ -1274,12 +1285,13 @@ void System::boundary_project_vector (const std::set<boundary_id_type> & b,
12741285
GradientFunctionPointer gptr,
12751286
const Parameters & function_parameters,
12761287
NumericVector<Number> & new_vector,
1277-
int is_adjoint) const
1288+
int is_adjoint,
1289+
std::optional<ConstElemRange> active_local_range) const
12781290
{
12791291
WrappedFunction<Number> f(*this, fptr, &function_parameters);
12801292
WrappedFunction<Gradient> g(*this, gptr, &function_parameters);
12811293
this->boundary_project_vector(b, variables, new_vector, &f, &g,
1282-
is_adjoint);
1294+
is_adjoint, active_local_range);
12831295
}
12841296

12851297
/**
@@ -1291,13 +1303,20 @@ void System::boundary_project_vector (const std::set<boundary_id_type> & b,
12911303
NumericVector<Number> & new_vector,
12921304
FunctionBase<Number> * f,
12931305
FunctionBase<Gradient> * g,
1294-
int is_adjoint) const
1306+
int is_adjoint,
1307+
std::optional<ConstElemRange> active_local_range) const
12951308
{
12961309
LOG_SCOPE ("boundary_project_vector()", "System");
12971310

1311+
if (!active_local_range)
1312+
{
1313+
active_local_range.emplace
1314+
(this->get_mesh().active_local_elements_begin(),
1315+
this->get_mesh().active_local_elements_end());
1316+
}
1317+
12981318
Threads::parallel_for
1299-
(ConstElemRange (this->get_mesh().active_local_elements_begin(),
1300-
this->get_mesh().active_local_elements_end() ),
1319+
(active_local_range.value(),
13011320
BoundaryProjectSolution(b, variables, *this, f, g,
13021321
this->get_equation_systems().parameters,
13031322
new_vector)

0 commit comments

Comments
 (0)