Line data Source code
1 : //* This file is part of the MOOSE framework
2 : //* https://mooseframework.inl.gov
3 : //*
4 : //* All rights reserved, see COPYRIGHT for full restrictions
5 : //* https://github.com/idaholab/moose/blob/master/COPYRIGHT
6 : //*
7 : //* Licensed under LGPL 2.1, please see LICENSE for details
8 : //* https://www.gnu.org/licenses/lgpl-2.1.html
9 :
10 : #include "ElementSubdomainModifierBase.h"
11 : #include "DisplacedProblem.h"
12 : #include "MaterialWarehouse.h"
13 :
14 : #include "libmesh/dof_map.h"
15 : #include "libmesh/remote_elem.h"
16 : #include "libmesh/parallel_ghost_sync.h"
17 : #include "libmesh/numeric_vector.h"
18 : #include "libmesh/parameters.h"
19 : #include <iterator>
20 : #include <unordered_set>
21 :
22 : InputParameters
23 44297 : ElementSubdomainModifierBase::validParams()
24 : {
25 44297 : InputParameters params = ElementUserObject::validParams();
26 :
27 : // ESMs only operate on the undisplaced mesh to gather subdomain and sideset information. This
28 : // information is used to modify both the undisplaced and the displaced meshes. It is the
29 : // developer's responsibility to make sure the element IDs and sideset info are consistent across
30 : // both meshes.
31 88594 : params.set<bool>("use_displaced_mesh") = false;
32 88594 : params.suppressParameter<bool>("use_displaced_mesh");
33 :
34 177188 : params.addParam<std::vector<BoundaryName>>(
35 : "moving_boundaries",
36 : {},
37 : "Moving boundaries between subdomains. These boundaries (both sidesets and nodesets) will be "
38 : "updated for elements that change subdomain. The subdomains that each moving "
39 : "boundary lies between shall be specified using the parameter "
40 : "'moving_boundary_subdomain_pairs'. If one boundary and multiple subdomain pairs are "
41 : "specified, then it is assumed that the pairs all apply to the boundary. A boundary will be "
42 : "created on the mesh if it does not already exist.");
43 177188 : params.addParam<std::vector<std::vector<SubdomainName>>>(
44 : "moving_boundary_subdomain_pairs",
45 : {},
46 : "The subdomain pairs associated with each moving boundary. For each pair of subdomains, only "
47 : "the element side from the first subdomain will be added to the moving boundary, i.e., the "
48 : "side normal is pointing from the first subdomain to the second subdomain. The pairs shall "
49 : "be delimited by ';'. If a pair only has one subdomain, the moving boundary is associated "
50 : "with the subdomain's external boundary, i.e., when the elements have no neighboring "
51 : "elements.");
52 :
53 221485 : params.addParam<std::vector<SubdomainName>>(
54 : "reinitialize_subdomains",
55 : {"ANY_BLOCK_ID"},
56 : "By default, any element which changes subdomain is reinitialized. If a list of subdomains "
57 : "(IDs or names) is provided, then only elements whose new subdomain is in the list will be "
58 : "reinitialized. If an empty list is set, then no elements will be reinitialized.");
59 :
60 132891 : params.addParam<bool>(
61 : "old_subdomain_reinitialized",
62 88594 : true,
63 : "This parameter must be set with a non-empty list in 'reinitialize_subdomains'. When set to "
64 : "the default true, the element's old subdomain is not considered when determining if an "
65 : "element should be reinitialized. If set to false, only elements whose old subdomain was not "
66 : "in 'reinitialize_subdomains' are reinitialized. ");
67 :
68 177188 : params.addParam<std::vector<VariableName>>(
69 : "reinitialize_variables", {}, "Which variables to reinitialize when subdomain changes.");
70 177188 : MooseEnum reinit_strategy("IC POLYNOMIAL_NEIGHBOR POLYNOMIAL_WHOLE POLYNOMIAL_NEARBY", "IC");
71 221485 : params.addParam<std::vector<MooseEnum>>(
72 : "reinitialization_strategy",
73 : {reinit_strategy},
74 : "The strategy used to reinitialize the solution when elements change subdomain. If multiple "
75 : "strategies are provided, each strategy will be applied to the corresponding variable. If "
76 : "only one strategy is provided, it will be applied to all variables.");
77 177188 : params.addParam<std::vector<UserObjectName>>(
78 : "polynomial_fitters",
79 : {},
80 : "List of NodalPatchRecovery UserObjects used for polynomial fitting during variable "
81 : "reinitialization. Required only if 'reinitialization_strategy' includes "
82 : "POLYNOMIAL_NEIGHBOR, POLYNOMIAL_WHOLE, or POLYNOMIAL_NEARBY.");
83 132891 : params.addParam<int>(
84 : "nearby_kd_tree_leaf_max_size",
85 88594 : 10,
86 : "Maximum number of elements in a leaf node of the K-D tree used to search for nearby "
87 : "elements. Only needed if 'reinitialization_strategy' is set to POLYNOMIAL_NEARBY.");
88 177188 : params.addParam<double>(
89 : "nearby_distance_threshold",
90 88594 : -1.0,
91 : "Threshold for considering elements as 'nearby' in the K-D tree search. Only elements within "
92 : "this distance will be considered for polynomial fitting.");
93 177188 : params.addParam<std::vector<bool>>(
94 : "restore_overridden_dofs",
95 : {},
96 : "A list of boolean flags, one for each variable in 'reinitialize_variables', specifying "
97 : "whether overridden DOF values should be restored after reinitialization for each variable. "
98 : "This is useful when the solved values on these DOFs should be preserved. If the list is "
99 : "empty, overridden DOF values will NOT be restored for any variable by default.");
100 :
101 132891 : params.addParam<bool>("skip_restore_subdomain_changes",
102 88594 : false,
103 : "Skip restoring the subdomain changes if the timestep is not advanced.");
104 :
105 44297 : params.registerBase("MeshModifier");
106 :
107 88594 : return params;
108 132891 : }
109 :
110 779 : ElementSubdomainModifierBase::ElementSubdomainModifierBase(const InputParameters & parameters)
111 : : ElementUserObject(parameters),
112 779 : _displaced_problem(_fe_problem.getDisplacedProblem().get()),
113 779 : _displaced_mesh(_displaced_problem ? &_displaced_problem->mesh() : nullptr),
114 779 : _nl_sys(_fe_problem.getNonlinearSystemBase(systemNumber())),
115 779 : _aux_sys(_fe_problem.getAuxiliarySystem()),
116 1558 : _t_step_old(declareRestartableData<int>("t_step_old", 0)),
117 779 : _restep(false),
118 1558 : _skip_restore_subdomain_changes(getParam<bool>("skip_restore_subdomain_changes")),
119 1558 : _old_subdomain_reinitialized(getParam<bool>("old_subdomain_reinitialized")),
120 1558 : _pr_names(getParam<std::vector<UserObjectName>>("polynomial_fitters")),
121 1558 : _reinit_vars(getParam<std::vector<VariableName>>("reinitialize_variables")),
122 1558 : _leaf_max_size(getParam<int>("nearby_kd_tree_leaf_max_size")),
123 7790 : _nearby_distance_threshold(getParam<double>("nearby_distance_threshold"))
124 : {
125 779 : }
126 :
127 : void
128 770 : ElementSubdomainModifierBase::initialSetup()
129 : {
130 : const std::vector<SubdomainName> subdomain_names_to_reinitialize =
131 1540 : getParam<std::vector<SubdomainName>>("reinitialize_subdomains");
132 :
133 770 : if (std::find(subdomain_names_to_reinitialize.begin(),
134 : subdomain_names_to_reinitialize.end(),
135 1540 : "ANY_BLOCK_ID") != subdomain_names_to_reinitialize.end())
136 465 : _subdomain_ids_to_reinitialize.push_back(Moose::ANY_BLOCK_ID);
137 : else
138 305 : _subdomain_ids_to_reinitialize = _mesh.getSubdomainIDs(subdomain_names_to_reinitialize);
139 :
140 : std::set<SubdomainID> set_subdomain_ids_to_reinitialize(_subdomain_ids_to_reinitialize.begin(),
141 770 : _subdomain_ids_to_reinitialize.end());
142 :
143 1242 : if (_old_subdomain_reinitialized == false &&
144 236 : (std::find(_subdomain_ids_to_reinitialize.begin(),
145 : _subdomain_ids_to_reinitialize.end(),
146 1242 : Moose::ANY_BLOCK_ID) != _subdomain_ids_to_reinitialize.end() ||
147 236 : set_subdomain_ids_to_reinitialize == _mesh.meshSubdomains()))
148 0 : paramError("old_subdomain_reinitialized",
149 : "'old_subdomain_reinitialized' can only be set to false if "
150 : "reinitialize_subdomains does "
151 : "not cover the whole model, otherwise no elements will be reinitialized as it is "
152 : "impossible for an element's old subdomain to not be in the list.");
153 770 : else if (_old_subdomain_reinitialized == false && _subdomain_ids_to_reinitialize.empty())
154 0 : paramError("old_subdomain_reinitialized",
155 : "'old_subdomain_reinitialized' can only be set to false if "
156 : "reinitialize_subdomains is set to a non-empty list of subdomains, otherwise no "
157 : "elements will be reinitialized, as it is impossible for an element's new subdomain "
158 : "to be in the list.");
159 :
160 1540 : auto bnd_names = getParam<std::vector<BoundaryName>>("moving_boundaries");
161 770 : auto bnd_ids = _mesh.getBoundaryIDs(bnd_names, true);
162 : const auto bnd_subdomains =
163 1540 : getParam<std::vector<std::vector<SubdomainName>>>("moving_boundary_subdomain_pairs");
164 :
165 770 : if (bnd_names.size() == 1 && bnd_subdomains.size() > 1)
166 : {
167 226 : bnd_names.insert(bnd_names.end(), bnd_subdomains.size() - 1, bnd_names[0]);
168 226 : bnd_ids.insert(bnd_ids.end(), bnd_subdomains.size() - 1, bnd_ids[0]);
169 : }
170 544 : else if (bnd_names.size() != bnd_subdomains.size())
171 0 : paramError("moving_boundary_subdomain_pairs",
172 : "Each moving boundary must correspond to a pair of subdomains. ",
173 : bnd_names.size(),
174 : " boundaries are specified by the parameter 'moving_boundaries', while ",
175 : bnd_subdomains.size(),
176 : " subdomain pairs are provided. Alternatively, if one boundary and multiple "
177 : "subdomain pairs are provided, then the subdomain pairs all apply to one boundary.");
178 :
179 1478 : for (auto i : index_range(bnd_names))
180 : {
181 708 : _moving_boundary_names[bnd_ids[i]] = bnd_names[i];
182 :
183 708 : if (bnd_subdomains[i].size() == 2)
184 426 : _moving_boundaries[{_mesh.getSubdomainID(bnd_subdomains[i][0]),
185 852 : _mesh.getSubdomainID(bnd_subdomains[i][1])}] = bnd_ids[i];
186 282 : else if (bnd_subdomains[i].size() == 1)
187 282 : _moving_boundaries[{_mesh.getSubdomainID(bnd_subdomains[i][0]), Moose::INVALID_BLOCK_ID}] =
188 282 : bnd_ids[i];
189 : else
190 0 : paramError("moving_boundary_subdomain_pairs",
191 : "Each subdomain pair must contain 1 or 2 subdomain names, but ",
192 0 : bnd_subdomains[i].size(),
193 : " are given.");
194 : }
195 :
196 : // Check if the variables to reinitialize exist in the system
197 1180 : for (const auto & var_name : _reinit_vars)
198 410 : if (!_nl_sys.hasVariable(var_name) && !_aux_sys.hasVariable(var_name))
199 0 : paramError("reinitialize_variables", "Variable ", var_name, " does not exist.");
200 :
201 : // Determine the reinitialization strategy for each variable.
202 : // (1) If they are of the same size, we perform a 1-to-1 mapping.
203 : // (2) If only one strategy or restore flag is provided, it applies to all variables.
204 1540 : const auto reinit_strategy_in = getParam<std::vector<MooseEnum>>("reinitialization_strategy");
205 1540 : const auto restore_overridden_dofs_in = getParam<std::vector<bool>>("restore_overridden_dofs");
206 :
207 770 : if (std::any_of(reinit_strategy_in.begin(),
208 : reinit_strategy_in.end(),
209 2498 : [](const MooseEnum & val) { return val == "POLYNOMIAL_NEARBY"; }) &&
210 854 : !isParamSetByUser("nearby_distance_threshold"))
211 0 : mooseError("The 'nearby_distance_threshold' parameter must be set when using the "
212 : "POLYNOMIAL_NEARBY reinitialization strategy.");
213 :
214 770 : if (reinit_strategy_in.size() == 1)
215 692 : _reinit_strategy.resize(_reinit_vars.size(), reinit_strategy_in[0].getEnum<ReinitStrategy>());
216 78 : else if (reinit_strategy_in.size() == _reinit_vars.size())
217 300 : for (const auto & e : reinit_strategy_in)
218 226 : _reinit_strategy.push_back(e.getEnum<ReinitStrategy>());
219 : else
220 8 : paramError(
221 : "reinitialization_strategy",
222 : "The 'reinitialization_strategy' parameter must have either a single value or a number "
223 : "of values equal to the number of 'reinitialize_variables'. "
224 : "Got ",
225 : reinit_strategy_in.size(),
226 : " strategies for ",
227 : _reinit_vars.size(),
228 : " variables.");
229 :
230 766 : if (restore_overridden_dofs_in.size() == 1)
231 : {
232 190 : if (restore_overridden_dofs_in[0])
233 : _vars_to_restore_overridden_dofs =
234 190 : _reinit_vars; // Restore overridden DOFs for all reinitialized variables
235 : }
236 576 : else if (restore_overridden_dofs_in.size() == _reinit_vars.size())
237 : {
238 614 : for (auto i : index_range(_reinit_vars))
239 56 : if (restore_overridden_dofs_in[i])
240 28 : _vars_to_restore_overridden_dofs.push_back(_reinit_vars[i]);
241 : }
242 : else
243 : {
244 18 : if (!restore_overridden_dofs_in.empty())
245 8 : paramError(
246 : "restore_overridden_dofs",
247 : "The 'restore_overridden_dofs' parameter must have either a single value or a number "
248 : "of values equal to the number of 'reinitialize_variables'. "
249 : "Got ",
250 : restore_overridden_dofs_in.size(),
251 : " restore_overridden_dofs for ",
252 : _reinit_vars.size(),
253 : " variables.");
254 : }
255 :
256 : // For all the other variables, we set the reinitialization strategy to IC
257 2795 : for (const auto & var_name : _fe_problem.getVariableNames())
258 2033 : if (std::find(_reinit_vars.begin(), _reinit_vars.end(), var_name) == _reinit_vars.end())
259 : {
260 1647 : _reinit_vars.push_back(var_name);
261 1647 : _reinit_strategy.push_back(ReinitStrategy::IC);
262 762 : }
263 :
264 : // Map variable names to the index of the nodal patch recovery user object
265 762 : _pr.resize(_pr_names.size());
266 762 : std::size_t pr_count = 0;
267 2787 : for (auto i : index_range(_reinit_vars))
268 2029 : if (_reinit_strategy[i] == ReinitStrategy::POLYNOMIAL_NEIGHBOR ||
269 3788 : _reinit_strategy[i] == ReinitStrategy::POLYNOMIAL_WHOLE ||
270 1759 : _reinit_strategy[i] == ReinitStrategy::POLYNOMIAL_NEARBY)
271 : {
272 298 : _var_name_to_pr_idx[_reinit_vars[i]] = pr_count;
273 298 : if (pr_count >= _pr_names.size())
274 8 : paramError("polynomial_fitters",
275 : "The number of polynomial fitters (",
276 : _pr_names.size(),
277 : ") is less than the number of variables to reinitialize with polynomial "
278 : "extrapolation.");
279 588 : _pr[pr_count] =
280 294 : &_fe_problem.getUserObject<NodalPatchRecoveryVariable>(_pr_names[pr_count], /*tid=*/0);
281 : mooseAssert(_pr[pr_count]->variableName() == _reinit_vars[i],
282 : "The patch recovery UserObject's variable name must match the variable being "
283 : "reinitialized in ElementSubdomainModifierBase.");
284 294 : _depend_uo.insert(_pr_names[pr_count]);
285 294 : pr_count++;
286 : }
287 758 : if (_pr_names.size() != pr_count)
288 8 : paramError("polynomial_fitters",
289 : "Mismatch between number of reinitialization strategies using polynomial "
290 : "extrapolation and polynomial fitters (expected: ",
291 : pr_count,
292 : ", given: ",
293 : _pr_names.size(),
294 : ").");
295 754 : }
296 :
297 : void
298 3497 : ElementSubdomainModifierBase::timestepSetup()
299 : {
300 3497 : if (_t_step == _t_step_old && !_skip_restore_subdomain_changes)
301 : {
302 29 : mooseInfoRepeated(name(), ": Restoring element subdomain changes.");
303 :
304 : // Reverse the subdomain changes
305 29 : auto moved_elem_reversed = _moved_elems;
306 617 : for (auto & [elem_id, subdomain] : moved_elem_reversed)
307 588 : std::swap(subdomain.first, subdomain.second);
308 :
309 29 : _restep = true;
310 29 : modify(moved_elem_reversed);
311 29 : _restep = false;
312 29 : }
313 :
314 3497 : _t_step_old = _t_step;
315 3497 : }
316 :
317 : void
318 3842 : ElementSubdomainModifierBase::modify(
319 : const std::unordered_map<dof_id_type, std::pair<SubdomainID, SubdomainID>> & moved_elems)
320 : {
321 3842 : if (!_restep)
322 : // Cache the moved elements for potential restore
323 3813 : _moved_elems = moved_elems;
324 :
325 : // If nothing need to change, just return.
326 : // This will skip all mesh changes, and so no adaptivity mesh files will be written.
327 3842 : auto n_moved_elem = moved_elems.size();
328 3842 : gatherSum(n_moved_elem);
329 3842 : if (n_moved_elem == 0)
330 1623 : return;
331 :
332 : // Create moving boundaries on the undisplaced and displaced meshes
333 : //
334 : // Note: We do this _everytime_ because previous execution might have removed the sidesets and
335 : // nodesets. Most of the moving boundary algorithms below assume that the moving sidesets and
336 : // nodesets already exist on the mesh.
337 2219 : createMovingBoundaries(_mesh);
338 2219 : if (_displaced_mesh)
339 124 : createMovingBoundaries(*_displaced_mesh);
340 :
341 : // This has to be done *before* subdomain changes are applied
342 2219 : findReinitializedElemsAndNodes(moved_elems);
343 :
344 : // Apply cached subdomain changes
345 2219 : applySubdomainChanges(moved_elems, _mesh);
346 2219 : if (_displaced_mesh)
347 124 : applySubdomainChanges(moved_elems, *_displaced_mesh);
348 :
349 : // Update moving boundaries
350 2219 : gatherMovingBoundaryChanges(moved_elems);
351 2219 : applyMovingBoundaryChanges(_mesh);
352 2219 : if (_displaced_mesh)
353 124 : applyMovingBoundaryChanges(*_displaced_mesh);
354 :
355 : // Some variable reinitialization strategies require patch elements to be gathered
356 : // This has to be done *before* reinitializing the equation systems because we need to find
357 : // currently evaluable elements
358 2219 : if (!_restep)
359 : {
360 2191 : _evaluable_elems.clear();
361 2191 : _patch_elem_ids.clear();
362 7989 : for (auto i : index_range(_reinit_vars))
363 5798 : prepareVariableForReinitialization(_reinit_vars[i], _reinit_strategy[i]);
364 : }
365 :
366 : // Reinit equation systems
367 2219 : _fe_problem.meshChanged(
368 : /*intermediate_change=*/false, /*contract_mesh=*/false, /*clean_refinement_flags=*/false);
369 :
370 : // Initialize solution and stateful material properties
371 2219 : if (!_restep)
372 : {
373 2191 : applyIC();
374 2191 : if (_fe_problem.getMaterialWarehouse().hasActiveObjects(0))
375 595 : initElementStatefulProps();
376 : }
377 : }
378 :
379 : void
380 2343 : ElementSubdomainModifierBase::createMovingBoundaries(MooseMesh & mesh)
381 : {
382 2343 : auto & bnd_info = mesh.getMesh().get_boundary_info();
383 4085 : for (const auto & [bnd_id, bnd_name] : _moving_boundary_names)
384 : {
385 1742 : bnd_info.sideset_name(bnd_id) = bnd_name;
386 1742 : bnd_info.nodeset_name(bnd_id) = bnd_name;
387 : }
388 2343 : }
389 :
390 : void
391 2343 : ElementSubdomainModifierBase::applySubdomainChanges(
392 : const std::unordered_map<dof_id_type, std::pair<SubdomainID, SubdomainID>> & moved_elems,
393 : MooseMesh & mesh)
394 : {
395 37028 : for (const auto & [elem_id, subdomain] : moved_elems)
396 : {
397 : // Change the element's subdomain ID
398 34685 : auto elem = mesh.elemPtr(elem_id);
399 34685 : const auto & [from, to] = subdomain;
400 : mooseAssert(elem->subdomain_id() == from, "Inconsistent element subdomain ID.");
401 34685 : elem->subdomain_id() = to;
402 :
403 : // Change the ancestors' (if any) subdomain ID
404 34685 : setAncestorsSubdomainIDs(elem, to);
405 : }
406 :
407 : // Synchronize ghost element subdomain changes
408 2343 : libMesh::SyncSubdomainIds sync(mesh.getMesh());
409 2343 : Parallel::sync_dofobject_data_by_id(
410 4686 : mesh.getMesh().comm(), mesh.getMesh().elements_begin(), mesh.getMesh().elements_end(), sync);
411 2343 : }
412 :
413 : void
414 34685 : ElementSubdomainModifierBase::setAncestorsSubdomainIDs(Elem * elem, const SubdomainID subdomain_id)
415 : {
416 34685 : auto curr_elem = elem;
417 :
418 43733 : for (unsigned int i = curr_elem->level(); i > 0; --i)
419 : {
420 : // Change the parent's subdomain
421 9048 : curr_elem = curr_elem->parent();
422 9048 : curr_elem->subdomain_id() = subdomain_id;
423 : }
424 34685 : }
425 :
426 : void
427 2219 : ElementSubdomainModifierBase::gatherMovingBoundaryChanges(
428 : const std::unordered_map<dof_id_type, std::pair<SubdomainID, SubdomainID>> & moved_elems)
429 : {
430 : // Clear moving boundary changes from last execution
431 2219 : _add_element_sides.clear();
432 2219 : _add_neighbor_sides.clear();
433 2219 : _remove_element_sides.clear();
434 2219 : _remove_neighbor_sides.clear();
435 :
436 2219 : const auto & sidesets = _mesh.getMesh().get_boundary_info().get_sideset_map();
437 :
438 36226 : for (const auto & [elem_id, subdomain_assignment] : moved_elems)
439 : {
440 34007 : auto elem = _mesh.elemPtr(elem_id);
441 :
442 : // The existing moving boundaries on the element side should be removed
443 46537 : for (auto itr = sidesets.lower_bound(elem); itr != sidesets.upper_bound(elem); itr++)
444 12530 : if (_moving_boundary_names.count(itr->second.second))
445 3792 : _remove_element_sides[elem->id()].emplace(itr->second.first, itr->second.second);
446 :
447 177483 : for (auto side : elem->side_index_range())
448 : {
449 143476 : auto neigh = elem->neighbor_ptr(side);
450 :
451 : // Don't mess with remote element neighbor
452 143476 : if (neigh && neigh == libMesh::remote_elem)
453 0 : continue;
454 : // If neighbor doesn't exist
455 143476 : else if (!neigh)
456 9863 : gatherMovingBoundaryChangesHelper(elem, side, nullptr, 0);
457 : // If neighbor exists
458 : else
459 : {
460 133613 : auto neigh_side = neigh->which_neighbor_am_i(elem);
461 :
462 133613 : if (neigh->active())
463 130513 : gatherMovingBoundaryChangesHelper(elem, side, neigh, neigh_side);
464 : else
465 : {
466 : // Find the active neighbors of the element
467 3100 : std::vector<const Elem *> active_neighs;
468 : // Neighbor has active children, they are neighbors of the element along that side
469 : mooseAssert(!neigh->subactive(),
470 : "The case where the active neighbor is an ancestor of this neighbor is not "
471 : "handled at this time.");
472 3100 : neigh->active_family_tree_by_neighbor(active_neighs, elem);
473 :
474 9300 : for (auto active_neigh : active_neighs)
475 6200 : gatherMovingBoundaryChangesHelper(elem, side, active_neigh, neigh_side);
476 3100 : }
477 : }
478 : }
479 : }
480 2219 : }
481 :
482 : void
483 146576 : ElementSubdomainModifierBase::gatherMovingBoundaryChangesHelper(const Elem * elem,
484 : unsigned short side,
485 : const Elem * neigh,
486 : unsigned short neigh_side)
487 : {
488 146576 : const auto & sidesets = _mesh.getMesh().get_boundary_info().get_sideset_map();
489 :
490 : // Detect element side change
491 146576 : SubdomainPair subdomain_pair = {elem->subdomain_id(),
492 146576 : neigh ? neigh->subdomain_id() : Moose::INVALID_BLOCK_ID};
493 146576 : if (_moving_boundaries.count(subdomain_pair))
494 11787 : _add_element_sides[elem->id()].emplace(side, _moving_boundaries.at(subdomain_pair));
495 :
496 146576 : if (neigh)
497 : {
498 : // The existing moving boundaries on the neighbor side should be removed
499 184957 : for (auto itr = sidesets.lower_bound(neigh); itr != sidesets.upper_bound(neigh); itr++)
500 48244 : if (itr->second.first == neigh_side && _moving_boundary_names.count(itr->second.second))
501 11633 : _remove_neighbor_sides[neigh->id()].emplace(itr->second.first, itr->second.second);
502 :
503 : // Detect neighbor side change (by reversing the subdomain pair)
504 136713 : subdomain_pair = {subdomain_pair.second, subdomain_pair.first};
505 136713 : if (_moving_boundaries.count(subdomain_pair))
506 3152 : _add_neighbor_sides[neigh->id()].emplace(neigh_side, _moving_boundaries.at(subdomain_pair));
507 : }
508 146576 : }
509 :
510 : void
511 2343 : ElementSubdomainModifierBase::applyMovingBoundaryChanges(MooseMesh & mesh)
512 : {
513 2343 : auto & bnd_info = mesh.getMesh().get_boundary_info();
514 :
515 : // Remove all boundary nodes from the previous moving boundaries
516 2343 : auto nodesets = bnd_info.get_nodeset_map();
517 214551 : for (const auto & [node_id, bnd] : nodesets)
518 212208 : if (_moving_boundary_names.count(bnd))
519 36302 : bnd_info.remove_node(node_id, bnd);
520 :
521 : // Keep track of ghost element changes
522 : std::unordered_map<processor_id_type,
523 : std::vector<std::tuple<dof_id_type, unsigned short, BoundaryID>>>
524 2343 : add_ghost_sides, remove_ghost_sides;
525 :
526 : // Remove element sides from moving boundaries
527 6057 : for (const auto & [elem_id, sides] : _remove_element_sides)
528 7654 : for (const auto & [side, bnd] : sides)
529 3940 : bnd_info.remove_side(mesh.elemPtr(elem_id), side, bnd);
530 :
531 : // Remove neighbor sides from moving boundaries
532 12573 : for (const auto & [elem_id, sides] : _remove_neighbor_sides)
533 : {
534 10230 : auto elem = mesh.elemPtr(elem_id);
535 21863 : for (const auto & [side, bnd] : sides)
536 : {
537 11633 : bnd_info.remove_side(elem, side, bnd);
538 : // Keep track of changes to ghosted elements
539 11633 : if (elem->processor_id() != processor_id())
540 125 : remove_ghost_sides[elem->processor_id()].push_back({elem_id, side, bnd});
541 : }
542 : }
543 :
544 2343 : Parallel::push_parallel_vector_data(
545 2343 : bnd_info.comm(),
546 : remove_ghost_sides,
547 2343 : [&mesh,
548 : &bnd_info](processor_id_type,
549 : const std::vector<std::tuple<dof_id_type, unsigned short, BoundaryID>> & received)
550 : {
551 198 : for (const auto & [elem_id, side, bnd] : received)
552 125 : bnd_info.remove_side(mesh.elemPtr(elem_id), side, bnd);
553 73 : });
554 :
555 : // Add element sides to moving boundaries
556 10843 : for (const auto & [elem_id, sides] : _add_element_sides)
557 20247 : for (const auto & [side, bnd] : sides)
558 11747 : bnd_info.add_side(mesh.elemPtr(elem_id), side, bnd);
559 :
560 : // Add neighbor sides to moving boundaries
561 4767 : for (const auto & [elem_id, sides] : _add_neighbor_sides)
562 : {
563 2424 : auto elem = mesh.elemPtr(elem_id);
564 4950 : for (const auto & [side, bnd] : sides)
565 : {
566 2526 : bnd_info.add_side(elem, side, bnd);
567 : // Keep track of changes to ghosted elements
568 2526 : if (elem->processor_id() != processor_id())
569 13 : add_ghost_sides[elem->processor_id()].push_back({elem_id, side, bnd});
570 : }
571 : }
572 :
573 2343 : Parallel::push_parallel_vector_data(
574 2343 : bnd_info.comm(),
575 : add_ghost_sides,
576 2343 : [&mesh,
577 : &bnd_info](processor_id_type,
578 : const std::vector<std::tuple<dof_id_type, unsigned short, BoundaryID>> & received)
579 : {
580 22 : for (const auto & [elem_id, side, bnd] : received)
581 13 : bnd_info.add_side(mesh.elemPtr(elem_id), side, bnd);
582 9 : });
583 :
584 2343 : bnd_info.parallel_sync_side_ids();
585 2343 : bnd_info.parallel_sync_node_ids();
586 2343 : mesh.update();
587 2343 : }
588 :
589 : void
590 5798 : ElementSubdomainModifierBase::prepareVariableForReinitialization(const VariableName & var_name,
591 : ReinitStrategy reinit_strategy)
592 : {
593 5798 : switch (reinit_strategy)
594 : {
595 5033 : case ReinitStrategy::IC:
596 : // No additional preparation needed for IC
597 5033 : break;
598 765 : case ReinitStrategy::POLYNOMIAL_NEIGHBOR:
599 : case ReinitStrategy::POLYNOMIAL_WHOLE:
600 : case ReinitStrategy::POLYNOMIAL_NEARBY:
601 : {
602 765 : if (_var_name_to_pr_idx.find(var_name) == _var_name_to_pr_idx.end())
603 0 : return;
604 765 : const int pr_idx = _var_name_to_pr_idx[var_name];
605 : // The patch elements might be different for each variable
606 765 : gatherPatchElements(var_name, reinit_strategy);
607 :
608 : // Notify the patch recovery user object about the patch elements
609 765 : _pr[pr_idx]->sync(_patch_elem_ids[var_name]);
610 :
611 765 : break;
612 : }
613 0 : default:
614 0 : mooseError("Unknown reinitialization strategy");
615 : break;
616 : }
617 : }
618 :
619 : void
620 3839 : ElementSubdomainModifierBase::meshChanged()
621 : {
622 : // Clear cached ranges
623 3839 : _reinitialized_elem_range.reset();
624 3839 : _reinitialized_bnd_node_range.reset();
625 3839 : _reinitialized_node_range.reset();
626 :
627 3839 : updateAMRMovingBoundary(_mesh);
628 3839 : if (_displaced_mesh)
629 216 : updateAMRMovingBoundary(*_displaced_mesh);
630 3839 : }
631 :
632 : void
633 4055 : ElementSubdomainModifierBase::updateAMRMovingBoundary(MooseMesh & mesh)
634 : {
635 4055 : auto & bnd_info = mesh.getMesh().get_boundary_info();
636 4055 : auto sidesets = bnd_info.get_sideset_map();
637 290494 : for (const auto & i : sidesets)
638 : {
639 286439 : auto elem = i.first;
640 286439 : auto side = i.second.first;
641 286439 : auto bnd = i.second.second;
642 286439 : if (_moving_boundary_names.count(bnd) && !elem->active())
643 : {
644 5834 : bnd_info.remove_side(elem, side, bnd);
645 :
646 5834 : std::vector<const Elem *> elem_family;
647 5834 : elem->active_family_tree_by_side(elem_family, side);
648 18649 : for (auto felem : elem_family)
649 12815 : bnd_info.add_side(felem, side, bnd);
650 5834 : }
651 : }
652 :
653 4055 : bnd_info.parallel_sync_side_ids();
654 4055 : bnd_info.parallel_sync_node_ids();
655 4055 : }
656 :
657 : void
658 2219 : ElementSubdomainModifierBase::findReinitializedElemsAndNodes(
659 : const std::unordered_map<dof_id_type, std::pair<SubdomainID, SubdomainID>> & moved_elems)
660 : {
661 : // Clear cached element reinitialization data
662 2219 : _reinitialized_elems.clear();
663 2219 : _reinitialized_nodes.clear();
664 :
665 : // One more algorithm:
666 : // (1) Loop over moved elements
667 : // (2) If neighbor element processor ID is not the same as current processor ID (ghost element),
668 : // push the moved element ID to the neighbor processor
669 :
670 2219 : std::unordered_map<processor_id_type, std::unordered_set<dof_id_type>> push_data_set;
671 2219 : std::unordered_map<processor_id_type, std::vector<dof_id_type>> push_data;
672 :
673 36226 : for (const auto & [elem_id, subdomain] : moved_elems)
674 : {
675 : mooseAssert(_mesh.elemPtr(elem_id)->active(), "Moved elements should be active");
676 : // Default: any element that changes subdomain is reinitialized
677 34007 : if (std::find(_subdomain_ids_to_reinitialize.begin(),
678 : _subdomain_ids_to_reinitialize.end(),
679 68014 : Moose::ANY_BLOCK_ID) != _subdomain_ids_to_reinitialize.end())
680 24691 : _reinitialized_elems.insert(elem_id);
681 : else // Reinitialize if new subdomain is in list of subdomains to be reinitialized
682 : {
683 9316 : const auto & [from, to] = subdomain;
684 9316 : if (subdomainIsReinitialized(to) && _old_subdomain_reinitialized)
685 88 : _reinitialized_elems.insert(elem_id);
686 : // Only reinitialize if original subdomain is not in list of subdomains
687 14631 : else if (subdomainIsReinitialized(to) && !_old_subdomain_reinitialized &&
688 5403 : !subdomainIsReinitialized(from))
689 4440 : _reinitialized_elems.insert(elem_id);
690 : else // New subdomain is not in list of subdomains
691 4788 : continue;
692 : }
693 29219 : const auto & elem = _mesh.elemPtr(elem_id);
694 :
695 : // (1) Loop over nodes of moved elements
696 : // (2) node to element map is used to find neighbor elements
697 : // (3) If neighbor element processor ID is not the same as current processor ID (means that the
698 : // current element is ghosted element to the neighbor processor), push the moved element (or
699 : // reinitialized, or newly-activated) ID to the neighbor processor
700 162825 : for (const auto & node : elem->node_ref_range())
701 854651 : for (const auto & neigh_id : _mesh.nodeToElemMap().at(node.id()))
702 721045 : if (neigh_id != elem_id) // Don't check the element itself
703 : {
704 587439 : const auto neigh_elem = _mesh.elemPtr(neigh_id);
705 587439 : if (neigh_elem->processor_id() != processor_id())
706 18295 : push_data_set[neigh_elem->processor_id()].insert(elem_id);
707 : }
708 :
709 162825 : for (unsigned int i = 0; i < elem->n_nodes(); ++i)
710 133606 : if (nodeIsNewlyReinitialized(elem->node_id(i)))
711 16572 : _reinitialized_nodes.insert(elem->node_id(i));
712 : }
713 :
714 2826 : for (auto & [pid, s] : push_data_set)
715 1821 : push_data[pid] = {s.begin(), s.end()};
716 :
717 2219 : _semi_local_reinitialized_elems = _reinitialized_elems;
718 :
719 : auto push_receiver =
720 607 : [this](const processor_id_type, const std::vector<dof_id_type> & received_data)
721 : {
722 2659 : for (const auto & id : received_data)
723 2052 : _semi_local_reinitialized_elems.insert(id);
724 607 : };
725 :
726 2219 : Parallel::push_parallel_vector_data(_mesh.comm(), push_data, push_receiver);
727 2219 : }
728 :
729 : bool
730 220368 : ElementSubdomainModifierBase::subdomainIsReinitialized(SubdomainID id) const
731 : {
732 : // Default: any element that changes subdomain is reinitialized
733 220368 : if (std::find(_subdomain_ids_to_reinitialize.begin(),
734 : _subdomain_ids_to_reinitialize.end(),
735 440736 : Moose::ANY_BLOCK_ID) != _subdomain_ids_to_reinitialize.end())
736 109792 : return true;
737 :
738 : // Is subdomain in list of subdomains to be reinitialized
739 110576 : return std::find(_subdomain_ids_to_reinitialize.begin(),
740 : _subdomain_ids_to_reinitialize.end(),
741 221152 : id) != _subdomain_ids_to_reinitialize.end();
742 : }
743 :
744 : bool
745 133606 : ElementSubdomainModifierBase::nodeIsNewlyReinitialized(dof_id_type node_id) const
746 : {
747 : // If any of the node neighbor elements has reinitialized, then the node is NOT newly
748 : // reinitialized.
749 212993 : for (auto neighbor_elem_id : _mesh.nodeToElemMap().at(node_id))
750 196421 : if (subdomainIsReinitialized(_mesh.elemPtr(neighbor_elem_id)->subdomain_id()))
751 117034 : return false;
752 16572 : return true;
753 : }
754 :
755 : void
756 2191 : ElementSubdomainModifierBase::applyIC()
757 : {
758 : // Before reinitializing variables, some DOFs may be overwritten.
759 : // By default, these overwritten DOF values are NOT restored.
760 : // If the user sets `restore_overridden_dofs` to true, we first save the current
761 : // values of these DOFs, then restore them after reinitialization.
762 3067 : for (const auto & var_name : _vars_to_restore_overridden_dofs)
763 876 : storeOverriddenDofValues(var_name);
764 :
765 : // ReinitStrategy::IC
766 2191 : std::set<VariableName> ic_vars;
767 7989 : for (auto i : index_range(_reinit_vars))
768 5798 : if (_reinit_strategy[i] == ReinitStrategy::IC)
769 5033 : ic_vars.insert(_reinit_vars[i]);
770 2191 : if (!ic_vars.empty())
771 2191 : _fe_problem.projectInitialConditionOnCustomRange(
772 : reinitializedElemRange(), reinitializedBndNodeRange(), ic_vars);
773 :
774 : // ReinitStrategy::POLYNOMIAL_NEIGHBOR, POLYNOMIAL_WHOLE, POLYNOMIAL_NEARBY
775 2956 : for (const auto & [var, patch] : _patch_elem_ids)
776 765 : extrapolatePolynomial(var);
777 :
778 : // See the comment above, now we restore the values of the dofs that were overridden
779 3067 : for (const auto & var_name : _vars_to_restore_overridden_dofs)
780 876 : restoreOverriddenDofValues(var_name);
781 :
782 : mooseAssert(_fe_problem.numSolverSystems() <= 1,
783 : "This code was written for a single nonlinear system");
784 : // Set old and older solutions on the reinitialized dofs to the reinitialized values
785 : // note: from current -> old -> older
786 2191 : setOldAndOlderSolutions(_fe_problem.getNonlinearSystemBase(_sys.number()),
787 : reinitializedElemRange(),
788 : reinitializedBndNodeRange());
789 4382 : setOldAndOlderSolutions(
790 2191 : _fe_problem.getAuxiliarySystem(), reinitializedElemRange(), reinitializedBndNodeRange());
791 :
792 : // Note: Need method to handle solve failures at timesteps where subdomain changes. The old
793 : // solutions are now set to the reinitialized values. Does this impact restoring solutions
794 2191 : }
795 :
796 : void
797 876 : ElementSubdomainModifierBase::storeOverriddenDofValues(const VariableName & var_name)
798 : {
799 876 : const auto & sys = _fe_problem.getSystem(var_name);
800 876 : const auto & current_solution = *sys.current_local_solution;
801 876 : const auto & dof_map = sys.get_dof_map();
802 876 : const auto & var = _fe_problem.getStandardVariable(0, var_name);
803 876 : const auto var_num = var.number();
804 :
805 : // Get the DOFs on the reinitialized elements
806 : // Here we should loop over both ghosted and local reinitialized elements.
807 : // The ghosted elements here can take care of DoFs that is belong to the reinitialized
808 : // elements but are not on the current processor.
809 876 : std::set<dof_id_type> reinitialized_dofs;
810 5964 : for (const auto & elem_id : _semi_local_reinitialized_elems)
811 : {
812 5088 : const auto & elem = _mesh.elemPtr(elem_id);
813 5088 : std::vector<dof_id_type> elem_dofs;
814 5088 : dof_map.dof_indices(elem, elem_dofs, var_num);
815 5088 : reinitialized_dofs.insert(elem_dofs.begin(), elem_dofs.end());
816 5088 : }
817 :
818 : // Get existing DOFs on the active elements excluding reinitialized elements
819 876 : std::set<dof_id_type> existing_dofs;
820 47328 : for (const auto * elem : *_mesh.getActiveLocalElementRange())
821 : {
822 46452 : if (_reinitialized_elems.count(elem->id()))
823 4514 : continue; // Skip reinitialized elements
824 41938 : std::vector<dof_id_type> elem_dofs;
825 41938 : dof_map.dof_indices(elem, elem_dofs, var_num);
826 41938 : existing_dofs.insert(elem_dofs.begin(), elem_dofs.end());
827 41938 : }
828 :
829 : // Get the DOFs on the nodes that are overridden on reinitialized elements
830 876 : std::vector<dof_id_type> overridden_dofs;
831 876 : std::set_intersection(reinitialized_dofs.begin(),
832 : reinitialized_dofs.end(),
833 : existing_dofs.begin(),
834 : existing_dofs.end(),
835 : std::back_inserter(overridden_dofs));
836 :
837 : // Values before overriding (to be restored later)
838 876 : std::vector<Number> values;
839 5180 : for (auto dof : overridden_dofs)
840 4304 : values.push_back(current_solution(dof));
841 :
842 876 : _overridden_values_on_reinit_elems[var_name] = {overridden_dofs, values};
843 876 : }
844 :
845 : void
846 876 : ElementSubdomainModifierBase::restoreOverriddenDofValues(const VariableName & var_name)
847 : {
848 876 : const auto sn = _fe_problem.systemNumForVariable(var_name);
849 876 : auto & sys = _fe_problem.getSystemBase(sn);
850 876 : auto & sol = sys.solution();
851 876 : const auto & dof_map = sys.dofMap();
852 876 : const auto & [dof_ids, values] = _overridden_values_on_reinit_elems[var_name];
853 :
854 876 : std::unordered_map<processor_id_type, std::vector<std::pair<dof_id_type, Number>>> push_data;
855 :
856 5180 : for (const int i : index_range(dof_ids))
857 : {
858 4304 : if (dof_map.dof_owner(dof_ids[i]) == processor_id())
859 4207 : sol.set(dof_ids[i], values[i]);
860 : else
861 97 : push_data[dof_map.dof_owner(dof_ids[i])].emplace_back(dof_ids[i], values[i]);
862 : }
863 :
864 74 : auto push_receiver = [&](const processor_id_type,
865 : const std::vector<std::pair<dof_id_type, Number>> & received_data)
866 : {
867 171 : for (const auto & [id, value] : received_data)
868 97 : sol.set(id, value);
869 74 : };
870 :
871 876 : Parallel::push_parallel_vector_data(_mesh.comm(), push_data, push_receiver);
872 :
873 876 : sol.close();
874 876 : sol.localize(*sys.system().current_local_solution, sys.dofMap().get_send_list());
875 876 : }
876 :
877 : void
878 595 : ElementSubdomainModifierBase::initElementStatefulProps()
879 : {
880 595 : _fe_problem.initElementStatefulProps(reinitializedElemRange(), /*threaded=*/true);
881 595 : }
882 :
883 : ConstElemRange &
884 10124 : ElementSubdomainModifierBase::reinitializedElemRange()
885 : {
886 10124 : if (_reinitialized_elem_range)
887 7933 : return *_reinitialized_elem_range.get();
888 :
889 : // Create a vector of the newly reinitialized elements
890 2191 : std::vector<Elem *> elems;
891 31242 : for (auto elem_id : _reinitialized_elems)
892 29051 : elems.push_back(_mesh.elemPtr(elem_id));
893 :
894 : // Make some fake element iterators defining this vector of elements
895 2191 : Elem * const * elem_itr_begin = const_cast<Elem * const *>(elems.data());
896 2191 : Elem * const * elem_itr_end = elem_itr_begin + elems.size();
897 :
898 : const auto elems_begin = MeshBase::const_element_iterator(
899 2191 : elem_itr_begin, elem_itr_end, Predicates::NotNull<Elem * const *>());
900 : const auto elems_end = MeshBase::const_element_iterator(
901 2191 : elem_itr_end, elem_itr_end, Predicates::NotNull<Elem * const *>());
902 :
903 2191 : _reinitialized_elem_range = std::make_unique<ConstElemRange>(elems_begin, elems_end);
904 :
905 2191 : return reinitializedElemRange();
906 2191 : }
907 :
908 : ConstBndNodeRange &
909 8764 : ElementSubdomainModifierBase::reinitializedBndNodeRange()
910 : {
911 8764 : if (_reinitialized_bnd_node_range)
912 6573 : return *_reinitialized_bnd_node_range.get();
913 :
914 : // Create a vector of the newly reinitialized boundary nodes
915 2191 : std::vector<const BndNode *> nodes;
916 2191 : auto bnd_nodes = _mesh.getBoundaryNodeRange();
917 204758 : for (auto bnd_node : *bnd_nodes)
918 202567 : if (bnd_node->_node)
919 202567 : if (_reinitialized_nodes.count(bnd_node->_node->id()))
920 4657 : nodes.push_back(bnd_node);
921 :
922 2191 : BndNode * const * bnd_node_itr_begin = const_cast<BndNode * const *>(nodes.data());
923 2191 : BndNode * const * bnd_node_itr_end = bnd_node_itr_begin + nodes.size();
924 :
925 : const auto bnd_nodes_begin = MooseMesh::const_bnd_node_iterator(
926 2191 : bnd_node_itr_begin, bnd_node_itr_end, Predicates::NotNull<const BndNode * const *>());
927 : const auto bnd_nodes_end = MooseMesh::const_bnd_node_iterator(
928 2191 : bnd_node_itr_end, bnd_node_itr_end, Predicates::NotNull<const BndNode * const *>());
929 :
930 : _reinitialized_bnd_node_range =
931 2191 : std::make_unique<ConstBndNodeRange>(bnd_nodes_begin, bnd_nodes_end);
932 :
933 2191 : return reinitializedBndNodeRange();
934 2191 : }
935 :
936 : ConstNodeRange &
937 0 : ElementSubdomainModifierBase::reinitializedNodeRange()
938 : {
939 0 : if (_reinitialized_node_range)
940 0 : return *_reinitialized_node_range.get();
941 :
942 : // Create a vector of the newly reinitialized nodes
943 0 : std::vector<const Node *> nodes;
944 :
945 0 : for (auto node_id : _reinitialized_nodes)
946 0 : nodes.push_back(_mesh.nodePtr(node_id)); // displaced mesh shares the same node object
947 :
948 0 : Node * const * node_itr_begin = const_cast<Node * const *>(nodes.data());
949 0 : Node * const * node_itr_end = node_itr_begin + nodes.size();
950 :
951 : const auto nodes_begin = MeshBase::const_node_iterator(
952 0 : node_itr_begin, node_itr_end, Predicates::NotNull<const Node * const *>());
953 : const auto nodes_end = MeshBase::const_node_iterator(
954 0 : node_itr_end, node_itr_end, Predicates::NotNull<const Node * const *>());
955 :
956 0 : _reinitialized_node_range = std::make_unique<ConstNodeRange>(nodes_begin, nodes_end);
957 :
958 0 : return *_reinitialized_node_range.get();
959 0 : }
960 :
961 : void
962 4382 : ElementSubdomainModifierBase::setOldAndOlderSolutions(SystemBase & sys,
963 : ConstElemRange & elem_range,
964 : ConstBndNodeRange & bnd_node_range)
965 : {
966 13696 : for (auto bnd : bnd_node_range)
967 : {
968 9314 : const Node * bnode = bnd->_node;
969 9314 : if (!bnode)
970 0 : continue;
971 : }
972 :
973 : // Don't do anything if this is a steady simulation
974 4382 : if (!sys.hasSolutionState(1))
975 24 : return;
976 :
977 4358 : NumericVector<Number> & current_solution = *sys.system().current_local_solution;
978 4358 : NumericVector<Number> & old_solution = sys.solutionOld();
979 4358 : NumericVector<Number> * older_solution = sys.hasSolutionState(2) ? &sys.solutionOlder() : nullptr;
980 :
981 : // Get dofs for the reinitialized elements and nodes
982 4358 : DofMap & dof_map = sys.dofMap();
983 4358 : std::vector<dof_id_type> dofs;
984 :
985 62082 : for (auto & elem : elem_range)
986 : {
987 57724 : std::vector<dof_id_type> elem_dofs;
988 57724 : dof_map.dof_indices(elem, elem_dofs);
989 57724 : dofs.insert(dofs.end(), elem_dofs.begin(), elem_dofs.end());
990 57724 : }
991 :
992 13672 : for (auto & bnd_node : bnd_node_range)
993 : {
994 9314 : std::vector<dof_id_type> bnd_node_dofs;
995 9314 : dof_map.dof_indices(bnd_node->_node, bnd_node_dofs);
996 9314 : dofs.insert(dofs.end(), bnd_node_dofs.begin(), bnd_node_dofs.end());
997 9314 : }
998 :
999 : // Set the old and older solutions to match the reinitialization
1000 253213 : for (auto dof : dofs)
1001 : {
1002 248855 : old_solution.set(dof, current_solution(dof));
1003 248855 : if (older_solution)
1004 824 : older_solution->set(dof, current_solution(dof));
1005 : }
1006 :
1007 4358 : old_solution.close();
1008 4358 : if (older_solution)
1009 37 : older_solution->close();
1010 4358 : }
1011 :
1012 : void
1013 765 : ElementSubdomainModifierBase::gatherPatchElements(const VariableName & var_name,
1014 : ReinitStrategy reinit_strategy)
1015 : {
1016 765 : _patch_elem_ids[var_name].clear();
1017 :
1018 : // First collect all elements who own dofs in the current dofmap
1019 765 : auto & sys = _fe_problem.getSystem(var_name);
1020 :
1021 : // Cache evaluable elements for the system if not already done
1022 765 : if (!_evaluable_elems.count(sys.number()))
1023 : {
1024 469 : auto & [candidate_elems, candidate_elem_ids] = _evaluable_elems[sys.number()];
1025 469 : const auto & dof_map = sys.get_dof_map();
1026 469 : std::vector<dof_id_type> elem_dofs;
1027 469 : auto vn = sys.variable_number(static_cast<std::string>(var_name));
1028 31829 : for (const auto elem : *_mesh.getActiveLocalElementRange())
1029 : {
1030 31360 : if (std::find(_reinitialized_elems.begin(), _reinitialized_elems.end(), elem->id()) !=
1031 62720 : _reinitialized_elems.end())
1032 2952 : continue; // Skip elements that were reinitialized
1033 :
1034 28408 : dof_map.dof_indices(elem, elem_dofs, vn);
1035 28408 : if (!elem_dofs.empty())
1036 : {
1037 10456 : candidate_elems.insert(elem);
1038 10456 : candidate_elem_ids.push_back(elem->id());
1039 : }
1040 : }
1041 469 : }
1042 765 : auto & [candidate_elems, candidate_elem_ids] = _evaluable_elems[sys.number()];
1043 :
1044 : // Now we gather patch elements based on the reinit strategy
1045 765 : auto & patch_elems = _patch_elem_ids[var_name];
1046 :
1047 765 : switch (reinit_strategy)
1048 : {
1049 580 : case ReinitStrategy::POLYNOMIAL_NEIGHBOR:
1050 : {
1051 11926 : auto has_neighbor_in_reinit_elems = [&](const Elem * elem) -> bool
1052 : {
1053 66605 : for (const auto & node : elem->node_ref_range())
1054 314560 : for (const auto & neigh_id : _mesh.nodeToElemMap().at(node.id()))
1055 : // here we need to use _global_reinitialized_elems gathering from all processors
1056 259881 : if (_semi_local_reinitialized_elems.count(neigh_id))
1057 2912 : return true;
1058 9014 : return false;
1059 580 : };
1060 : // Loop over all candidate elements, for each element, if any of its point neighbor belongs
1061 : // to the reinitialized elements, we will include that element in the patch element set.
1062 12506 : for (const auto * elem : candidate_elems)
1063 11926 : if (has_neighbor_in_reinit_elems(elem))
1064 2912 : patch_elems.push_back(elem->id());
1065 580 : break;
1066 : }
1067 111 : case ReinitStrategy::POLYNOMIAL_WHOLE:
1068 : {
1069 : // This is simple: all candidate elements are patch elements
1070 111 : patch_elems = candidate_elem_ids;
1071 111 : break;
1072 : }
1073 74 : case ReinitStrategy::POLYNOMIAL_NEARBY:
1074 : {
1075 74 : std::vector<Point> kd_points;
1076 74 : std::vector<dof_id_type> global_candidate_elem_ids;
1077 :
1078 74 : if (_mesh.isDistributedMesh())
1079 : {
1080 12 : std::vector<std::pair<Point, dof_id_type>> pts_ids(candidate_elem_ids.size());
1081 116 : for (std::size_t i = 0; i < candidate_elem_ids.size(); ++i)
1082 208 : pts_ids[i] = {_mesh.elemPtr(candidate_elem_ids[i])->vertex_average(),
1083 208 : candidate_elem_ids[i]};
1084 12 : _mesh.comm().allgather(pts_ids);
1085 220 : for (const auto & [pt, id] : pts_ids)
1086 : {
1087 208 : kd_points.push_back(pt);
1088 208 : global_candidate_elem_ids.push_back(id);
1089 : }
1090 12 : }
1091 : else
1092 : {
1093 62 : _mesh.comm().allgather(candidate_elem_ids);
1094 62 : global_candidate_elem_ids = candidate_elem_ids;
1095 1146 : for (const auto & id : candidate_elem_ids)
1096 1084 : kd_points.push_back(_mesh.elemPtr(id)->vertex_average());
1097 : }
1098 :
1099 74 : const auto kd_tree = std::make_unique<KDTree>(kd_points, _leaf_max_size);
1100 :
1101 74 : std::vector<nanoflann::ResultItem<std::size_t, Real>> query_result;
1102 358 : for (const auto & elem_id : _reinitialized_elems)
1103 : {
1104 284 : const Point & centroid = _mesh.elemPtr(elem_id)->vertex_average();
1105 284 : kd_tree->radiusSearch(centroid, _nearby_distance_threshold, query_result);
1106 1504 : for (const auto & [qid, dist] : query_result)
1107 1220 : patch_elems.push_back(global_candidate_elem_ids[qid]);
1108 : }
1109 74 : break;
1110 74 : }
1111 0 : default:
1112 0 : mooseError("Unknown reinitialization strategy");
1113 : break;
1114 : }
1115 :
1116 : // every processor should have the same patch elements to do the polynomial extrapolation,
1117 : // so we gather them across all processors
1118 765 : _mesh.comm().allgather(patch_elems);
1119 :
1120 : // Remove duplicates from the patch elements (espcially important for POLYNOMIAL_NEARBY)
1121 765 : std::sort(patch_elems.begin(), patch_elems.end());
1122 765 : patch_elems.erase(std::unique(patch_elems.begin(), patch_elems.end()), patch_elems.end());
1123 765 : }
1124 :
1125 : void
1126 765 : ElementSubdomainModifierBase::extrapolatePolynomial(const VariableName & var_name)
1127 : {
1128 : const auto & coef =
1129 765 : _pr[_var_name_to_pr_idx[var_name]]->getCachedCoefficients(_patch_elem_ids[var_name]);
1130 :
1131 765 : const unsigned dim = _mesh.dimension();
1132 :
1133 765 : libMesh::Parameters function_parameters;
1134 :
1135 765 : const auto & multi_index = _pr[_var_name_to_pr_idx[var_name]]->multiIndex();
1136 :
1137 1530 : function_parameters.set<std::vector<std::vector<unsigned int>>>("multi_index") = multi_index;
1138 :
1139 765 : std::vector<Real> coef_vec(coef.size());
1140 3235 : for (auto i = 0; i < coef.size(); ++i)
1141 2470 : coef_vec[i] = coef(i);
1142 :
1143 765 : function_parameters.set<std::vector<Real>>("multi_index_coefficients") = coef_vec;
1144 1530 : function_parameters.set<unsigned int>("dimension_for_projection") = dim;
1145 :
1146 : // Define projection function
1147 66566 : auto poly_func = [](const Point & p,
1148 : const libMesh::Parameters & parameters,
1149 : const std::string &,
1150 : const std::string &) -> libMesh::Number
1151 : {
1152 : const auto & multi_index =
1153 66566 : parameters.get<std::vector<std::vector<unsigned int>>>("multi_index");
1154 66566 : const auto & coeffs = parameters.get<std::vector<Real>>("multi_index_coefficients");
1155 :
1156 66566 : Real val = 0.0;
1157 :
1158 338357 : for (unsigned int r = 0; r < multi_index.size(); r++)
1159 : {
1160 271791 : Real monomial = 1.0;
1161 918363 : for (unsigned int d = 0; d < multi_index[r].size(); d++)
1162 : {
1163 646572 : const auto power = multi_index[r][d];
1164 646572 : if (power == 0)
1165 410450 : continue;
1166 :
1167 236122 : monomial *= std::pow(p(d), power);
1168 : }
1169 271791 : val += coeffs[r] * monomial;
1170 : }
1171 :
1172 66566 : return val;
1173 : };
1174 :
1175 : // Define gradient
1176 0 : auto poly_func_grad = [](const Point & p,
1177 : const libMesh::Parameters & parameters,
1178 : const std::string &,
1179 : const std::string &) -> libMesh::Gradient
1180 : {
1181 0 : const unsigned int dim = parameters.get<unsigned int>("dimension_for_projection");
1182 :
1183 : const auto & multi_index =
1184 0 : parameters.get<std::vector<std::vector<unsigned int>>>("multi_index");
1185 0 : const auto & coeffs = parameters.get<std::vector<Real>>("multi_index_coefficients");
1186 :
1187 0 : libMesh::Gradient grad; // Zero-initialized
1188 :
1189 0 : for (unsigned int r = 0; r < multi_index.size(); ++r)
1190 : {
1191 0 : const auto & powers = multi_index[r];
1192 0 : const Real coef = coeffs[r];
1193 :
1194 0 : for (unsigned int d = 0; d < dim; ++d) // Loop over dimension
1195 : {
1196 0 : const auto power_d = powers[d];
1197 0 : if (power_d == 0)
1198 0 : continue;
1199 :
1200 : // Compute partial derivative in direction d
1201 0 : Real partial = coef * power_d;
1202 :
1203 0 : for (unsigned int i = 0; i < powers.size(); ++i)
1204 : {
1205 0 : if (i == d)
1206 : {
1207 0 : if (powers[i] > 1)
1208 0 : partial *= std::pow(p(i), powers[i] - 1); // reduce power by 1
1209 : }
1210 : else
1211 : {
1212 0 : if (powers[i] > 0)
1213 0 : partial *= std::pow(p(i), powers[i]); // full power
1214 : }
1215 : }
1216 :
1217 0 : grad(d) += partial;
1218 : }
1219 : }
1220 :
1221 0 : return grad;
1222 : };
1223 :
1224 765 : _fe_problem.projectFunctionOnCustomRange(
1225 : reinitializedElemRange(), poly_func, poly_func_grad, function_parameters, var_name);
1226 765 : }
|
