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 10594 : ElementSubdomainModifierBase::validParams()
24 : {
25 10594 : 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 21188 : params.set<bool>("use_displaced_mesh") = false;
32 21188 : params.suppressParameter<bool>("use_displaced_mesh");
33 :
34 42376 : 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 42376 : 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 52970 : 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 31782 : params.addParam<bool>(
61 : "old_subdomain_reinitialized",
62 21188 : 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 42376 : params.addParam<std::vector<VariableName>>(
69 : "reinitialize_variables", {}, "Which variables to reinitialize when subdomain changes.");
70 42376 : MooseEnum reinit_strategy("IC POLYNOMIAL_NEIGHBOR POLYNOMIAL_WHOLE POLYNOMIAL_NEARBY NONE", "IC");
71 52970 : 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 42376 : 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 31782 : params.addParam<int>(
84 : "nearby_kd_tree_leaf_max_size",
85 21188 : 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 42376 : params.addParam<double>(
89 : "nearby_distance_threshold",
90 21188 : -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 42376 : 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 31782 : params.addParam<bool>("skip_restore_subdomain_changes",
102 21188 : false,
103 : "Skip restoring the subdomain changes if the timestep is not advanced.");
104 :
105 10594 : params.registerBase("MeshModifier");
106 :
107 21188 : return params;
108 31782 : }
109 :
110 735 : ElementSubdomainModifierBase::ElementSubdomainModifierBase(const InputParameters & parameters)
111 : : ElementUserObject(parameters),
112 735 : _displaced_problem(_fe_problem.getDisplacedProblem().get()),
113 735 : _displaced_mesh(_displaced_problem ? &_displaced_problem->mesh() : nullptr),
114 735 : _nl_sys(_fe_problem.getNonlinearSystemBase(systemNumber())),
115 735 : _aux_sys(_fe_problem.getAuxiliarySystem()),
116 1470 : _t_step_old(declareRestartableData<int>("t_step_old", 0)),
117 735 : _restep(false),
118 1470 : _skip_restore_subdomain_changes(getParam<bool>("skip_restore_subdomain_changes")),
119 1470 : _old_subdomain_reinitialized(getParam<bool>("old_subdomain_reinitialized")),
120 1470 : _pr_names(getParam<std::vector<UserObjectName>>("polynomial_fitters")),
121 1470 : _reinit_vars(getParam<std::vector<VariableName>>("reinitialize_variables")),
122 1470 : _leaf_max_size(getParam<int>("nearby_kd_tree_leaf_max_size")),
123 7350 : _nearby_distance_threshold(getParam<double>("nearby_distance_threshold"))
124 : {
125 735 : }
126 :
127 : void
128 726 : ElementSubdomainModifierBase::initialSetup()
129 : {
130 : const std::vector<SubdomainName> subdomain_names_to_reinitialize =
131 1452 : getParam<std::vector<SubdomainName>>("reinitialize_subdomains");
132 :
133 726 : if (std::find(subdomain_names_to_reinitialize.begin(),
134 : subdomain_names_to_reinitialize.end(),
135 1452 : "ANY_BLOCK_ID") != subdomain_names_to_reinitialize.end())
136 429 : _subdomain_ids_to_reinitialize.push_back(Moose::ANY_BLOCK_ID);
137 : else
138 297 : _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 726 : _subdomain_ids_to_reinitialize.end());
142 :
143 1192 : if (_old_subdomain_reinitialized == false &&
144 233 : (std::find(_subdomain_ids_to_reinitialize.begin(),
145 : _subdomain_ids_to_reinitialize.end(),
146 1192 : Moose::ANY_BLOCK_ID) != _subdomain_ids_to_reinitialize.end() ||
147 233 : 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 726 : 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 1452 : auto bnd_names = getParam<std::vector<BoundaryName>>("moving_boundaries");
161 726 : auto bnd_ids = _mesh.getBoundaryIDs(bnd_names, true);
162 : const auto bnd_subdomains =
163 1452 : getParam<std::vector<std::vector<SubdomainName>>>("moving_boundary_subdomain_pairs");
164 :
165 726 : if (bnd_names.size() == 1 && bnd_subdomains.size() > 1)
166 : {
167 223 : bnd_names.insert(bnd_names.end(), bnd_subdomains.size() - 1, bnd_names[0]);
168 223 : bnd_ids.insert(bnd_ids.end(), bnd_subdomains.size() - 1, bnd_ids[0]);
169 : }
170 503 : 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 1409 : for (auto i : index_range(bnd_names))
180 : {
181 683 : _moving_boundary_names[bnd_ids[i]] = bnd_names[i];
182 :
183 683 : if (bnd_subdomains[i].size() == 2)
184 408 : _moving_boundaries[{_mesh.getSubdomainID(bnd_subdomains[i][0]),
185 816 : _mesh.getSubdomainID(bnd_subdomains[i][1])}] = bnd_ids[i];
186 275 : else if (bnd_subdomains[i].size() == 1)
187 275 : _moving_boundaries[{_mesh.getSubdomainID(bnd_subdomains[i][0]), Moose::INVALID_BLOCK_ID}] =
188 275 : 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 1149 : for (const auto & var_name : _reinit_vars)
198 423 : 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 1452 : const auto reinit_strategy_in = getParam<std::vector<MooseEnum>>("reinitialization_strategy");
205 1452 : const auto restore_overridden_dofs_in = getParam<std::vector<bool>>("restore_overridden_dofs");
206 :
207 726 : if (std::any_of(reinit_strategy_in.begin(),
208 : reinit_strategy_in.end(),
209 2381 : [](const MooseEnum & val) { return val == "POLYNOMIAL_NEARBY"; }) &&
210 804 : !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 726 : if (std::all_of(reinit_strategy_in.begin(),
215 : reinit_strategy_in.end(),
216 3055 : [](const MooseEnum & val) { return val != "POLYNOMIAL_NEARBY"; }) &&
217 2826 : (isParamSetByUser("nearby_distance_threshold") ||
218 2817 : isParamSetByUser("nearby_kd_tree_leaf_max_size")))
219 3 : mooseWarning("The 'nearby_distance_threshold' and 'nearby_kd_tree_leaf_max_size' parameters "
220 : "will be ignored because no 'reinitialization_strategy' is set to "
221 : "POLYNOMIAL_NEARBY.");
222 :
223 723 : if (reinit_strategy_in.size() == 1)
224 639 : _reinit_strategy.resize(_reinit_vars.size(), reinit_strategy_in[0].getEnum<ReinitStrategy>());
225 84 : else if (reinit_strategy_in.size() == _reinit_vars.size())
226 327 : for (const auto & e : reinit_strategy_in)
227 246 : _reinit_strategy.push_back(e.getEnum<ReinitStrategy>());
228 : else
229 6 : paramError(
230 : "reinitialization_strategy",
231 : "The 'reinitialization_strategy' parameter must have either a single value or a number "
232 : "of values equal to the number of 'reinitialize_variables'. "
233 : "Got ",
234 : reinit_strategy_in.size(),
235 : " strategies for ",
236 : _reinit_vars.size(),
237 : " variables.");
238 :
239 720 : if (restore_overridden_dofs_in.size() == 1)
240 : {
241 188 : if (restore_overridden_dofs_in[0])
242 : _vars_to_restore_overridden_dofs =
243 188 : _reinit_vars; // Restore overridden DOFs for all reinitialized variables
244 : }
245 532 : else if (restore_overridden_dofs_in.size() == _reinit_vars.size())
246 : {
247 568 : for (auto i : index_range(_reinit_vars))
248 52 : if (restore_overridden_dofs_in[i])
249 26 : _vars_to_restore_overridden_dofs.push_back(_reinit_vars[i]);
250 : }
251 : else
252 : {
253 16 : if (!restore_overridden_dofs_in.empty())
254 6 : paramError(
255 : "restore_overridden_dofs",
256 : "The 'restore_overridden_dofs' parameter must have either a single value or a number "
257 : "of values equal to the number of 'reinitialize_variables'. "
258 : "Got ",
259 : restore_overridden_dofs_in.size(),
260 : " restore_overridden_dofs for ",
261 : _reinit_vars.size(),
262 : " variables.");
263 : }
264 :
265 : // For all the other variables, we set the reinitialization strategy to IC
266 2644 : for (const auto & var_name : _fe_problem.getVariableNames())
267 1927 : if (std::find(_reinit_vars.begin(), _reinit_vars.end(), var_name) == _reinit_vars.end())
268 : {
269 1531 : _reinit_vars.push_back(var_name);
270 1531 : _reinit_strategy.push_back(ReinitStrategy::IC);
271 717 : }
272 :
273 : // Map variable names to the index of the nodal patch recovery user object
274 717 : _pr.resize(_pr_names.size());
275 717 : std::size_t pr_count = 0;
276 2638 : for (auto i : index_range(_reinit_vars))
277 1924 : if (_reinit_strategy[i] == ReinitStrategy::POLYNOMIAL_NEIGHBOR ||
278 3585 : _reinit_strategy[i] == ReinitStrategy::POLYNOMIAL_WHOLE ||
279 1661 : _reinit_strategy[i] == ReinitStrategy::POLYNOMIAL_NEARBY)
280 : {
281 289 : _var_name_to_pr_idx[_reinit_vars[i]] = pr_count;
282 289 : if (pr_count >= _pr_names.size())
283 6 : paramError("polynomial_fitters",
284 : "The number of polynomial fitters (",
285 : _pr_names.size(),
286 : ") is less than the number of variables to reinitialize with polynomial "
287 : "extrapolation.");
288 572 : _pr[pr_count] =
289 286 : &_fe_problem.getUserObject<NodalPatchRecoveryVariable>(_pr_names[pr_count], /*tid=*/0);
290 : mooseAssert(_pr[pr_count]->variableName() == _reinit_vars[i],
291 : "The patch recovery UserObject's variable name must match the variable being "
292 : "reinitialized in ElementSubdomainModifierBase.");
293 286 : _depend_uo.insert(_pr_names[pr_count]);
294 286 : pr_count++;
295 : }
296 714 : if (_pr_names.size() != pr_count)
297 6 : paramError("polynomial_fitters",
298 : "Mismatch between number of reinitialization strategies using polynomial "
299 : "extrapolation and polynomial fitters (expected: ",
300 : pr_count,
301 : ", given: ",
302 : _pr_names.size(),
303 : ").");
304 711 : }
305 :
306 : void
307 3258 : ElementSubdomainModifierBase::timestepSetup()
308 : {
309 3258 : if (_t_step == _t_step_old && !_skip_restore_subdomain_changes)
310 : {
311 30 : mooseInfoRepeated(name(), ": Restoring element subdomain changes.");
312 :
313 : // Reverse the subdomain changes
314 30 : auto moved_elem_reversed = _moved_elems;
315 625 : for (auto & [elem_id, subdomain] : moved_elem_reversed)
316 595 : std::swap(subdomain.first, subdomain.second);
317 :
318 30 : _restep = true;
319 30 : modify(moved_elem_reversed);
320 30 : _restep = false;
321 30 : }
322 :
323 3258 : _t_step_old = _t_step;
324 3258 : }
325 :
326 : void
327 3562 : ElementSubdomainModifierBase::modify(
328 : const std::unordered_map<dof_id_type, std::pair<SubdomainID, SubdomainID>> & moved_elems)
329 : {
330 3562 : if (!_restep)
331 : // Cache the moved elements for potential restore
332 3532 : _moved_elems = moved_elems;
333 :
334 : // If nothing need to change, just return.
335 : // This will skip all mesh changes, and so no adaptivity mesh files will be written.
336 3562 : auto n_moved_elem = moved_elems.size();
337 3562 : gatherSum(n_moved_elem);
338 3562 : if (n_moved_elem == 0)
339 1496 : return;
340 :
341 : // Create moving boundaries on the undisplaced and displaced meshes
342 : //
343 : // Note: We do this _everytime_ because previous execution might have removed the sidesets and
344 : // nodesets. Most of the moving boundary algorithms below assume that the moving sidesets and
345 : // nodesets already exist on the mesh.
346 2066 : createMovingBoundaries(_mesh);
347 2066 : if (_displaced_mesh)
348 114 : createMovingBoundaries(*_displaced_mesh);
349 :
350 : // This has to be done *before* subdomain changes are applied
351 2066 : findReinitializedElemsAndNodes(moved_elems);
352 :
353 : // Apply cached subdomain changes
354 2066 : applySubdomainChanges(moved_elems, _mesh);
355 2066 : if (_displaced_mesh)
356 114 : applySubdomainChanges(moved_elems, *_displaced_mesh);
357 :
358 : // Update moving boundaries
359 2066 : gatherMovingBoundaryChanges(moved_elems);
360 2066 : applyMovingBoundaryChanges(_mesh);
361 2066 : if (_displaced_mesh)
362 114 : applyMovingBoundaryChanges(*_displaced_mesh);
363 :
364 : // Some variable reinitialization strategies require patch elements to be gathered
365 : // This has to be done *before* reinitializing the equation systems because we need to find
366 : // currently evaluable elements
367 2066 : if (!_restep)
368 : {
369 2037 : _evaluable_elems.clear();
370 2037 : _patch_elem_ids.clear();
371 7457 : for (auto i : index_range(_reinit_vars))
372 5420 : prepareVariableForReinitialization(_reinit_vars[i], _reinit_strategy[i]);
373 : }
374 :
375 : // Reinit equation systems
376 2066 : _fe_problem.meshChanged(
377 : /*intermediate_change=*/false, /*contract_mesh=*/false, /*clean_refinement_flags=*/false);
378 :
379 : // Initialize solution and stateful material properties
380 2066 : if (!_restep)
381 : {
382 2037 : applyIC();
383 2037 : if (_fe_problem.getMaterialWarehouse().hasActiveObjects(0))
384 577 : initElementStatefulProps();
385 : }
386 : }
387 :
388 : void
389 2180 : ElementSubdomainModifierBase::createMovingBoundaries(MooseMesh & mesh)
390 : {
391 2180 : auto & bnd_info = mesh.getMesh().get_boundary_info();
392 3813 : for (const auto & [bnd_id, bnd_name] : _moving_boundary_names)
393 : {
394 1633 : bnd_info.sideset_name(bnd_id) = bnd_name;
395 1633 : bnd_info.nodeset_name(bnd_id) = bnd_name;
396 : }
397 2180 : }
398 :
399 : void
400 2180 : ElementSubdomainModifierBase::applySubdomainChanges(
401 : const std::unordered_map<dof_id_type, std::pair<SubdomainID, SubdomainID>> & moved_elems,
402 : MooseMesh & mesh)
403 : {
404 33055 : for (const auto & [elem_id, subdomain] : moved_elems)
405 : {
406 : // Change the element's subdomain ID
407 30875 : auto elem = mesh.elemPtr(elem_id);
408 30875 : const auto & [from, to] = subdomain;
409 : mooseAssert(elem->subdomain_id() == from, "Inconsistent element subdomain ID.");
410 30875 : elem->subdomain_id() = to;
411 :
412 : // Change the ancestors' (if any) subdomain ID
413 30875 : setAncestorsSubdomainIDs(elem, to);
414 : }
415 :
416 : // Synchronize ghost element subdomain changes
417 2180 : libMesh::SyncSubdomainIds sync(mesh.getMesh());
418 2180 : Parallel::sync_dofobject_data_by_id(
419 4360 : mesh.getMesh().comm(), mesh.getMesh().elements_begin(), mesh.getMesh().elements_end(), sync);
420 2180 : }
421 :
422 : void
423 30875 : ElementSubdomainModifierBase::setAncestorsSubdomainIDs(Elem * elem, const SubdomainID subdomain_id)
424 : {
425 30875 : auto curr_elem = elem;
426 :
427 38863 : for (unsigned int i = curr_elem->level(); i > 0; --i)
428 : {
429 : // Change the parent's subdomain
430 7988 : curr_elem = curr_elem->parent();
431 7988 : curr_elem->subdomain_id() = subdomain_id;
432 : }
433 30875 : }
434 :
435 : void
436 2066 : ElementSubdomainModifierBase::gatherMovingBoundaryChanges(
437 : const std::unordered_map<dof_id_type, std::pair<SubdomainID, SubdomainID>> & moved_elems)
438 : {
439 : // Clear moving boundary changes from last execution
440 2066 : _add_element_sides.clear();
441 2066 : _add_neighbor_sides.clear();
442 2066 : _remove_element_sides.clear();
443 2066 : _remove_neighbor_sides.clear();
444 :
445 2066 : const auto & sidesets = _mesh.getMesh().get_boundary_info().get_sideset_map();
446 :
447 32337 : for (const auto & [elem_id, subdomain_assignment] : moved_elems)
448 : {
449 30271 : auto elem = _mesh.elemPtr(elem_id);
450 :
451 : // The existing moving boundaries on the element side should be removed
452 41475 : for (auto itr = sidesets.lower_bound(elem); itr != sidesets.upper_bound(elem); itr++)
453 11204 : if (_moving_boundary_names.count(itr->second.second))
454 3392 : _remove_element_sides[elem->id()].emplace(itr->second.first, itr->second.second);
455 :
456 157947 : for (auto side : elem->side_index_range())
457 : {
458 127676 : auto neigh = elem->neighbor_ptr(side);
459 :
460 : // Don't mess with remote element neighbor
461 127676 : if (neigh && neigh == libMesh::remote_elem)
462 0 : continue;
463 : // If neighbor doesn't exist
464 127676 : else if (!neigh)
465 8808 : gatherMovingBoundaryChangesHelper(elem, side, nullptr, 0);
466 : // If neighbor exists
467 : else
468 : {
469 118868 : auto neigh_side = neigh->which_neighbor_am_i(elem);
470 :
471 118868 : if (neigh->active())
472 116118 : gatherMovingBoundaryChangesHelper(elem, side, neigh, neigh_side);
473 : else
474 : {
475 : // Find the active neighbors of the element
476 2750 : std::vector<const Elem *> active_neighs;
477 : // Neighbor has active children, they are neighbors of the element along that side
478 : mooseAssert(!neigh->subactive(),
479 : "The case where the active neighbor is an ancestor of this neighbor is not "
480 : "handled at this time.");
481 2750 : neigh->active_family_tree_by_neighbor(active_neighs, elem);
482 :
483 8250 : for (auto active_neigh : active_neighs)
484 5500 : gatherMovingBoundaryChangesHelper(elem, side, active_neigh, neigh_side);
485 2750 : }
486 : }
487 : }
488 : }
489 2066 : }
490 :
491 : void
492 130426 : ElementSubdomainModifierBase::gatherMovingBoundaryChangesHelper(const Elem * elem,
493 : unsigned short side,
494 : const Elem * neigh,
495 : unsigned short neigh_side)
496 : {
497 130426 : const auto & sidesets = _mesh.getMesh().get_boundary_info().get_sideset_map();
498 :
499 : // Detect element side change
500 130426 : SubdomainPair subdomain_pair = {elem->subdomain_id(),
501 130426 : neigh ? neigh->subdomain_id() : Moose::INVALID_BLOCK_ID};
502 130426 : if (_moving_boundaries.count(subdomain_pair))
503 10578 : _add_element_sides[elem->id()].emplace(side, _moving_boundaries.at(subdomain_pair));
504 :
505 130426 : if (neigh)
506 : {
507 : // The existing moving boundaries on the neighbor side should be removed
508 164728 : for (auto itr = sidesets.lower_bound(neigh); itr != sidesets.upper_bound(neigh); itr++)
509 43110 : if (itr->second.first == neigh_side && _moving_boundary_names.count(itr->second.second))
510 10340 : _remove_neighbor_sides[neigh->id()].emplace(itr->second.first, itr->second.second);
511 :
512 : // Detect neighbor side change (by reversing the subdomain pair)
513 121618 : subdomain_pair = {subdomain_pair.second, subdomain_pair.first};
514 121618 : if (_moving_boundaries.count(subdomain_pair))
515 2812 : _add_neighbor_sides[neigh->id()].emplace(neigh_side, _moving_boundaries.at(subdomain_pair));
516 : }
517 130426 : }
518 :
519 : void
520 2180 : ElementSubdomainModifierBase::applyMovingBoundaryChanges(MooseMesh & mesh)
521 : {
522 2180 : auto & bnd_info = mesh.getMesh().get_boundary_info();
523 :
524 : // Remove all boundary nodes from the previous moving boundaries
525 2180 : auto nodesets = bnd_info.get_nodeset_map();
526 202383 : for (const auto & [node_id, bnd] : nodesets)
527 200203 : if (_moving_boundary_names.count(bnd))
528 33203 : bnd_info.remove_node(node_id, bnd);
529 :
530 : // Keep track of ghost element changes
531 : std::unordered_map<processor_id_type,
532 : std::vector<std::tuple<dof_id_type, unsigned short, BoundaryID>>>
533 2180 : add_ghost_sides, remove_ghost_sides;
534 :
535 : // Remove element sides from moving boundaries
536 5498 : for (const auto & [elem_id, sides] : _remove_element_sides)
537 6842 : for (const auto & [side, bnd] : sides)
538 3524 : bnd_info.remove_side(mesh.elemPtr(elem_id), side, bnd);
539 :
540 : // Remove neighbor sides from moving boundaries
541 11272 : for (const auto & [elem_id, sides] : _remove_neighbor_sides)
542 : {
543 9092 : auto elem = mesh.elemPtr(elem_id);
544 19432 : for (const auto & [side, bnd] : sides)
545 : {
546 10340 : bnd_info.remove_side(elem, side, bnd);
547 : // Keep track of changes to ghosted elements
548 10340 : if (elem->processor_id() != processor_id())
549 128 : remove_ghost_sides[elem->processor_id()].push_back({elem_id, side, bnd});
550 : }
551 : }
552 :
553 2180 : Parallel::push_parallel_vector_data(
554 2180 : bnd_info.comm(),
555 : remove_ghost_sides,
556 2180 : [&mesh,
557 : &bnd_info](processor_id_type,
558 : const std::vector<std::tuple<dof_id_type, unsigned short, BoundaryID>> & received)
559 : {
560 204 : for (const auto & [elem_id, side, bnd] : received)
561 128 : bnd_info.remove_side(mesh.elemPtr(elem_id), side, bnd);
562 76 : });
563 :
564 : // Add element sides to moving boundaries
565 9808 : for (const auto & [elem_id, sides] : _add_element_sides)
566 18172 : for (const auto & [side, bnd] : sides)
567 10544 : bnd_info.add_side(mesh.elemPtr(elem_id), side, bnd);
568 :
569 : // Add neighbor sides to moving boundaries
570 4343 : for (const auto & [elem_id, sides] : _add_neighbor_sides)
571 : {
572 2163 : auto elem = mesh.elemPtr(elem_id);
573 4419 : for (const auto & [side, bnd] : sides)
574 : {
575 2256 : bnd_info.add_side(elem, side, bnd);
576 : // Keep track of changes to ghosted elements
577 2256 : if (elem->processor_id() != processor_id())
578 13 : add_ghost_sides[elem->processor_id()].push_back({elem_id, side, bnd});
579 : }
580 : }
581 :
582 2180 : Parallel::push_parallel_vector_data(
583 2180 : bnd_info.comm(),
584 : add_ghost_sides,
585 2180 : [&mesh,
586 : &bnd_info](processor_id_type,
587 : const std::vector<std::tuple<dof_id_type, unsigned short, BoundaryID>> & received)
588 : {
589 22 : for (const auto & [elem_id, side, bnd] : received)
590 13 : bnd_info.add_side(mesh.elemPtr(elem_id), side, bnd);
591 9 : });
592 :
593 2180 : bnd_info.parallel_sync_side_ids();
594 2180 : bnd_info.parallel_sync_node_ids();
595 2180 : mesh.update();
596 2180 : }
597 :
598 : void
599 5420 : ElementSubdomainModifierBase::prepareVariableForReinitialization(const VariableName & var_name,
600 : ReinitStrategy reinit_strategy)
601 : {
602 5420 : switch (reinit_strategy)
603 : {
604 4683 : case ReinitStrategy::NONE:
605 : case ReinitStrategy::IC:
606 : // No additional preparation needed for IC and NONE strategies
607 4683 : break;
608 737 : case ReinitStrategy::POLYNOMIAL_NEIGHBOR:
609 : case ReinitStrategy::POLYNOMIAL_WHOLE:
610 : case ReinitStrategy::POLYNOMIAL_NEARBY:
611 : {
612 737 : if (_var_name_to_pr_idx.find(var_name) == _var_name_to_pr_idx.end())
613 0 : return;
614 737 : const int pr_idx = _var_name_to_pr_idx[var_name];
615 : // The patch elements might be different for each variable
616 737 : gatherPatchElements(var_name, reinit_strategy);
617 :
618 : // Notify the patch recovery user object about the patch elements
619 737 : _pr[pr_idx]->sync(_patch_elem_ids[var_name]);
620 :
621 737 : break;
622 : }
623 0 : default:
624 0 : mooseError("Unknown reinitialization strategy");
625 : break;
626 : }
627 : }
628 :
629 : void
630 3577 : ElementSubdomainModifierBase::meshChanged()
631 : {
632 : // Clear cached ranges
633 3577 : _reinitialized_elem_range.reset();
634 3577 : _reinitialized_bnd_node_range.reset();
635 3577 : _reinitialized_node_range.reset();
636 :
637 3577 : updateAMRMovingBoundary(_mesh);
638 3577 : if (_displaced_mesh)
639 200 : updateAMRMovingBoundary(*_displaced_mesh);
640 3577 : }
641 :
642 : void
643 3777 : ElementSubdomainModifierBase::updateAMRMovingBoundary(MooseMesh & mesh)
644 : {
645 3777 : auto & bnd_info = mesh.getMesh().get_boundary_info();
646 3777 : auto sidesets = bnd_info.get_sideset_map();
647 267985 : for (const auto & i : sidesets)
648 : {
649 264208 : auto elem = i.first;
650 264208 : auto side = i.second.first;
651 264208 : auto bnd = i.second.second;
652 264208 : if (_moving_boundary_names.count(bnd) && !elem->active())
653 : {
654 5294 : bnd_info.remove_side(elem, side, bnd);
655 :
656 5294 : std::vector<const Elem *> elem_family;
657 5294 : elem->active_family_tree_by_side(elem_family, side);
658 16909 : for (auto felem : elem_family)
659 11615 : bnd_info.add_side(felem, side, bnd);
660 5294 : }
661 : }
662 :
663 3777 : bnd_info.parallel_sync_side_ids();
664 3777 : bnd_info.parallel_sync_node_ids();
665 3777 : }
666 :
667 : void
668 2066 : ElementSubdomainModifierBase::findReinitializedElemsAndNodes(
669 : const std::unordered_map<dof_id_type, std::pair<SubdomainID, SubdomainID>> & moved_elems)
670 : {
671 : // Clear cached element reinitialization data
672 2066 : _reinitialized_elems.clear();
673 2066 : _reinitialized_nodes.clear();
674 :
675 : // One more algorithm:
676 : // (1) Loop over moved elements
677 : // (2) If neighbor element processor ID is not the same as current processor ID (ghost element),
678 : // push the moved element ID to the neighbor processor
679 :
680 2066 : std::unordered_map<processor_id_type, std::unordered_set<dof_id_type>> push_data_set;
681 2066 : std::unordered_map<processor_id_type, std::vector<dof_id_type>> push_data;
682 :
683 32337 : for (const auto & [elem_id, subdomain] : moved_elems)
684 : {
685 : mooseAssert(_mesh.elemPtr(elem_id)->active(), "Moved elements should be active");
686 : // Default: any element that changes subdomain is reinitialized
687 30271 : if (std::find(_subdomain_ids_to_reinitialize.begin(),
688 : _subdomain_ids_to_reinitialize.end(),
689 60542 : Moose::ANY_BLOCK_ID) != _subdomain_ids_to_reinitialize.end())
690 21763 : _reinitialized_elems.insert(elem_id);
691 : else // Reinitialize if new subdomain is in list of subdomains to be reinitialized
692 : {
693 8508 : const auto & [from, to] = subdomain;
694 8508 : if (subdomainIsReinitialized(to) && _old_subdomain_reinitialized)
695 78 : _reinitialized_elems.insert(elem_id);
696 : // Only reinitialize if original subdomain is not in list of subdomains
697 13392 : else if (subdomainIsReinitialized(to) && !_old_subdomain_reinitialized &&
698 4962 : !subdomainIsReinitialized(from))
699 4104 : _reinitialized_elems.insert(elem_id);
700 : else // New subdomain is not in list of subdomains
701 4326 : continue;
702 : }
703 25945 : const auto & elem = _mesh.elemPtr(elem_id);
704 :
705 : // (1) Loop over nodes of moved elements
706 : // (2) node to element map is used to find neighbor elements
707 : // (3) If neighbor element processor ID is not the same as current processor ID (means that the
708 : // current element is ghosted element to the neighbor processor), push the moved element (or
709 : // reinitialized, or newly-activated) ID to the neighbor processor
710 144512 : for (const auto & node : elem->node_ref_range())
711 757887 : for (const auto & neigh_id : _mesh.nodeToElemMap().at(node.id()))
712 639320 : if (neigh_id != elem_id) // Don't check the element itself
713 : {
714 520753 : const auto neigh_elem = _mesh.elemPtr(neigh_id);
715 520753 : if (neigh_elem->processor_id() != processor_id())
716 18365 : push_data_set[neigh_elem->processor_id()].insert(elem_id);
717 : }
718 :
719 144512 : for (unsigned int i = 0; i < elem->n_nodes(); ++i)
720 118567 : if (nodeIsNewlyReinitialized(elem->node_id(i)))
721 15195 : _reinitialized_nodes.insert(elem->node_id(i));
722 : }
723 :
724 2682 : for (auto & [pid, s] : push_data_set)
725 1848 : push_data[pid] = {s.begin(), s.end()};
726 :
727 2066 : _semi_local_reinitialized_elems = _reinitialized_elems;
728 :
729 : auto push_receiver =
730 616 : [this](const processor_id_type, const std::vector<dof_id_type> & received_data)
731 : {
732 2688 : for (const auto & id : received_data)
733 2072 : _semi_local_reinitialized_elems.insert(id);
734 616 : };
735 :
736 2066 : Parallel::push_parallel_vector_data(_mesh.comm(), push_data, push_receiver);
737 2066 : }
738 :
739 : bool
740 197614 : ElementSubdomainModifierBase::subdomainIsReinitialized(SubdomainID id) const
741 : {
742 : // Default: any element that changes subdomain is reinitialized
743 197614 : if (std::find(_subdomain_ids_to_reinitialize.begin(),
744 : _subdomain_ids_to_reinitialize.end(),
745 395228 : Moose::ANY_BLOCK_ID) != _subdomain_ids_to_reinitialize.end())
746 96724 : return true;
747 :
748 : // Is subdomain in list of subdomains to be reinitialized
749 100890 : return std::find(_subdomain_ids_to_reinitialize.begin(),
750 : _subdomain_ids_to_reinitialize.end(),
751 201780 : id) != _subdomain_ids_to_reinitialize.end();
752 : }
753 :
754 : bool
755 118567 : ElementSubdomainModifierBase::nodeIsNewlyReinitialized(dof_id_type node_id) const
756 : {
757 : // If any of the node neighbor elements has reinitialized, then the node is NOT newly
758 : // reinitialized.
759 190909 : for (auto neighbor_elem_id : _mesh.nodeToElemMap().at(node_id))
760 175714 : if (subdomainIsReinitialized(_mesh.elemPtr(neighbor_elem_id)->subdomain_id()))
761 103372 : return false;
762 15195 : return true;
763 : }
764 :
765 : void
766 2037 : ElementSubdomainModifierBase::applyIC()
767 : {
768 : // Before reinitializing variables, some DOFs may be overwritten.
769 : // By default, these overwritten DOF values are NOT restored.
770 : // If the user sets `restore_overridden_dofs` to true, we first save the current
771 : // values of these DOFs, then restore them after reinitialization.
772 2944 : for (const auto & var_name : _vars_to_restore_overridden_dofs)
773 907 : storeOverriddenDofValues(var_name);
774 :
775 : // ReinitStrategy::IC
776 2037 : std::set<VariableName> ic_vars;
777 7457 : for (auto i : index_range(_reinit_vars))
778 5420 : if (_reinit_strategy[i] == ReinitStrategy::IC)
779 4649 : ic_vars.insert(_reinit_vars[i]);
780 2037 : if (!ic_vars.empty())
781 2037 : _fe_problem.projectInitialConditionOnCustomRange(
782 : reinitializedElemRange(), reinitializedBndNodeRange(), ic_vars);
783 :
784 : // ReinitStrategy::POLYNOMIAL_NEIGHBOR, POLYNOMIAL_WHOLE, POLYNOMIAL_NEARBY
785 2774 : for (const auto & [var, patch] : _patch_elem_ids)
786 737 : extrapolatePolynomial(var);
787 :
788 : // See the comment above, now we restore the values of the dofs that were overridden
789 2944 : for (const auto & var_name : _vars_to_restore_overridden_dofs)
790 907 : restoreOverriddenDofValues(var_name);
791 :
792 : mooseAssert(_fe_problem.numSolverSystems() <= 1,
793 : "This code was written for a single nonlinear system");
794 : // Set old and older solutions on the reinitialized dofs to the reinitialized values
795 : // note: from current -> old -> older
796 2037 : setOldAndOlderSolutions(_fe_problem.getNonlinearSystemBase(_sys.number()),
797 : reinitializedElemRange(),
798 : reinitializedBndNodeRange());
799 4074 : setOldAndOlderSolutions(
800 2037 : _fe_problem.getAuxiliarySystem(), reinitializedElemRange(), reinitializedBndNodeRange());
801 :
802 : // Note: Need method to handle solve failures at timesteps where subdomain changes. The old
803 : // solutions are now set to the reinitialized values. Does this impact restoring solutions
804 2037 : }
805 :
806 : void
807 907 : ElementSubdomainModifierBase::storeOverriddenDofValues(const VariableName & var_name)
808 : {
809 907 : const auto & sys = _fe_problem.getSystem(var_name);
810 907 : const auto & current_solution = *sys.current_local_solution;
811 907 : const auto & dof_map = sys.get_dof_map();
812 907 : const auto & var = _fe_problem.getStandardVariable(0, var_name);
813 907 : const auto var_num = var.number();
814 :
815 : // Get the DOFs on the reinitialized elements
816 : // Here we should loop over both ghosted and local reinitialized elements.
817 : // The ghosted elements here can take care of DoFs that is belong to the reinitialized
818 : // elements but are not on the current processor.
819 907 : std::set<dof_id_type> reinitialized_dofs;
820 5963 : for (const auto & elem_id : _semi_local_reinitialized_elems)
821 : {
822 5056 : const auto & elem = _mesh.elemPtr(elem_id);
823 5056 : std::vector<dof_id_type> elem_dofs;
824 5056 : dof_map.dof_indices(elem, elem_dofs, var_num);
825 5056 : reinitialized_dofs.insert(elem_dofs.begin(), elem_dofs.end());
826 5056 : }
827 :
828 : // Get existing DOFs on the active elements excluding reinitialized elements
829 907 : std::set<dof_id_type> existing_dofs;
830 46085 : for (const auto * elem : *_mesh.getActiveLocalElementRange())
831 : {
832 45178 : if (_reinitialized_elems.count(elem->id()))
833 4422 : continue; // Skip reinitialized elements
834 40756 : std::vector<dof_id_type> elem_dofs;
835 40756 : dof_map.dof_indices(elem, elem_dofs, var_num);
836 40756 : existing_dofs.insert(elem_dofs.begin(), elem_dofs.end());
837 40756 : }
838 :
839 : // Get the DOFs on the nodes that are overridden on reinitialized elements
840 907 : std::vector<dof_id_type> overridden_dofs;
841 907 : std::set_intersection(reinitialized_dofs.begin(),
842 : reinitialized_dofs.end(),
843 : existing_dofs.begin(),
844 : existing_dofs.end(),
845 : std::back_inserter(overridden_dofs));
846 :
847 : // Values before overriding (to be restored later)
848 907 : std::vector<Number> values;
849 5102 : for (auto dof : overridden_dofs)
850 4195 : values.push_back(current_solution(dof));
851 :
852 907 : _overridden_values_on_reinit_elems[var_name] = {overridden_dofs, values};
853 907 : }
854 :
855 : void
856 907 : ElementSubdomainModifierBase::restoreOverriddenDofValues(const VariableName & var_name)
857 : {
858 907 : const auto sn = _fe_problem.systemNumForVariable(var_name);
859 907 : auto & sys = _fe_problem.getSystemBase(sn);
860 907 : auto & sol = sys.solution();
861 907 : const auto & dof_map = sys.dofMap();
862 907 : const auto & [dof_ids, values] = _overridden_values_on_reinit_elems[var_name];
863 :
864 907 : std::unordered_map<processor_id_type, std::vector<std::pair<dof_id_type, Number>>> push_data;
865 :
866 5102 : for (const int i : index_range(dof_ids))
867 : {
868 4195 : if (dof_map.dof_owner(dof_ids[i]) == processor_id())
869 4089 : sol.set(dof_ids[i], values[i]);
870 : else
871 106 : push_data[dof_map.dof_owner(dof_ids[i])].emplace_back(dof_ids[i], values[i]);
872 : }
873 :
874 83 : auto push_receiver = [&](const processor_id_type,
875 : const std::vector<std::pair<dof_id_type, Number>> & received_data)
876 : {
877 189 : for (const auto & [id, value] : received_data)
878 106 : sol.set(id, value);
879 83 : };
880 :
881 907 : Parallel::push_parallel_vector_data(_mesh.comm(), push_data, push_receiver);
882 :
883 907 : sol.close();
884 907 : sol.localize(*sys.system().current_local_solution, sys.dofMap().get_send_list());
885 907 : }
886 :
887 : void
888 577 : ElementSubdomainModifierBase::initElementStatefulProps()
889 : {
890 577 : _fe_problem.initElementStatefulProps(reinitializedElemRange(), /*threaded=*/true);
891 577 : }
892 :
893 : ConstElemRange &
894 9462 : ElementSubdomainModifierBase::reinitializedElemRange()
895 : {
896 9462 : if (_reinitialized_elem_range)
897 7425 : return *_reinitialized_elem_range.get();
898 :
899 : // Create a vector of the newly reinitialized elements
900 2037 : std::vector<Elem *> elems;
901 27814 : for (auto elem_id : _reinitialized_elems)
902 25777 : elems.push_back(_mesh.elemPtr(elem_id));
903 :
904 : // Make some fake element iterators defining this vector of elements
905 2037 : Elem * const * elem_itr_begin = const_cast<Elem * const *>(elems.data());
906 2037 : Elem * const * elem_itr_end = elem_itr_begin + elems.size();
907 :
908 : const auto elems_begin = MeshBase::const_element_iterator(
909 2037 : elem_itr_begin, elem_itr_end, Predicates::NotNull<Elem * const *>());
910 : const auto elems_end = MeshBase::const_element_iterator(
911 2037 : elem_itr_end, elem_itr_end, Predicates::NotNull<Elem * const *>());
912 :
913 2037 : _reinitialized_elem_range = std::make_unique<ConstElemRange>(elems_begin, elems_end);
914 :
915 2037 : return reinitializedElemRange();
916 2037 : }
917 :
918 : ConstBndNodeRange &
919 8148 : ElementSubdomainModifierBase::reinitializedBndNodeRange()
920 : {
921 8148 : if (_reinitialized_bnd_node_range)
922 6111 : return *_reinitialized_bnd_node_range.get();
923 :
924 : // Create a vector of the newly reinitialized boundary nodes
925 2037 : std::vector<const BndNode *> nodes;
926 2037 : auto bnd_nodes = _mesh.getBoundaryNodeRange();
927 191314 : for (auto bnd_node : *bnd_nodes)
928 189277 : if (bnd_node->_node)
929 189277 : if (_reinitialized_nodes.count(bnd_node->_node->id()))
930 4412 : nodes.push_back(bnd_node);
931 :
932 2037 : BndNode * const * bnd_node_itr_begin = const_cast<BndNode * const *>(nodes.data());
933 2037 : BndNode * const * bnd_node_itr_end = bnd_node_itr_begin + nodes.size();
934 :
935 : const auto bnd_nodes_begin = MooseMesh::const_bnd_node_iterator(
936 2037 : bnd_node_itr_begin, bnd_node_itr_end, Predicates::NotNull<const BndNode * const *>());
937 : const auto bnd_nodes_end = MooseMesh::const_bnd_node_iterator(
938 2037 : bnd_node_itr_end, bnd_node_itr_end, Predicates::NotNull<const BndNode * const *>());
939 :
940 : _reinitialized_bnd_node_range =
941 2037 : std::make_unique<ConstBndNodeRange>(bnd_nodes_begin, bnd_nodes_end);
942 :
943 2037 : return reinitializedBndNodeRange();
944 2037 : }
945 :
946 : ConstNodeRange &
947 0 : ElementSubdomainModifierBase::reinitializedNodeRange()
948 : {
949 0 : if (_reinitialized_node_range)
950 0 : return *_reinitialized_node_range.get();
951 :
952 : // Create a vector of the newly reinitialized nodes
953 0 : std::vector<const Node *> nodes;
954 :
955 0 : for (auto node_id : _reinitialized_nodes)
956 0 : nodes.push_back(_mesh.nodePtr(node_id)); // displaced mesh shares the same node object
957 :
958 0 : Node * const * node_itr_begin = const_cast<Node * const *>(nodes.data());
959 0 : Node * const * node_itr_end = node_itr_begin + nodes.size();
960 :
961 : const auto nodes_begin = MeshBase::const_node_iterator(
962 0 : node_itr_begin, node_itr_end, Predicates::NotNull<const Node * const *>());
963 : const auto nodes_end = MeshBase::const_node_iterator(
964 0 : node_itr_end, node_itr_end, Predicates::NotNull<const Node * const *>());
965 :
966 0 : _reinitialized_node_range = std::make_unique<ConstNodeRange>(nodes_begin, nodes_end);
967 :
968 0 : return *_reinitialized_node_range.get();
969 0 : }
970 :
971 : void
972 4074 : ElementSubdomainModifierBase::setOldAndOlderSolutions(SystemBase & sys,
973 : ConstElemRange & elem_range,
974 : ConstBndNodeRange & bnd_node_range)
975 : {
976 12898 : for (auto bnd : bnd_node_range)
977 : {
978 8824 : const Node * bnode = bnd->_node;
979 8824 : if (!bnode)
980 0 : continue;
981 : }
982 :
983 : // Don't do anything if this is a steady simulation
984 4074 : if (!sys.hasSolutionState(1))
985 22 : return;
986 :
987 4052 : NumericVector<Number> & current_solution = *sys.system().current_local_solution;
988 4052 : NumericVector<Number> & old_solution = sys.solutionOld();
989 4052 : NumericVector<Number> * older_solution = sys.hasSolutionState(2) ? &sys.solutionOlder() : nullptr;
990 :
991 : // Get dofs for the reinitialized elements and nodes
992 4052 : DofMap & dof_map = sys.dofMap();
993 4052 : std::vector<dof_id_type> dofs;
994 :
995 55270 : for (auto & elem : elem_range)
996 : {
997 51218 : std::vector<dof_id_type> elem_dofs;
998 51218 : dof_map.dof_indices(elem, elem_dofs);
999 51218 : dofs.insert(dofs.end(), elem_dofs.begin(), elem_dofs.end());
1000 51218 : }
1001 :
1002 12876 : for (auto & bnd_node : bnd_node_range)
1003 : {
1004 8824 : std::vector<dof_id_type> bnd_node_dofs;
1005 8824 : dof_map.dof_indices(bnd_node->_node, bnd_node_dofs);
1006 8824 : dofs.insert(dofs.end(), bnd_node_dofs.begin(), bnd_node_dofs.end());
1007 8824 : }
1008 :
1009 : // Set the old and older solutions to match the reinitialization
1010 226585 : for (auto dof : dofs)
1011 : {
1012 222533 : old_solution.set(dof, current_solution(dof));
1013 222533 : if (older_solution)
1014 738 : older_solution->set(dof, current_solution(dof));
1015 : }
1016 :
1017 4052 : old_solution.close();
1018 4052 : if (older_solution)
1019 34 : older_solution->close();
1020 4052 : }
1021 :
1022 : void
1023 737 : ElementSubdomainModifierBase::gatherPatchElements(const VariableName & var_name,
1024 : ReinitStrategy reinit_strategy)
1025 : {
1026 737 : _patch_elem_ids[var_name].clear();
1027 :
1028 : // First collect all elements who own dofs in the current dofmap
1029 737 : auto & sys = _fe_problem.getSystem(var_name);
1030 :
1031 : // Cache evaluable elements for the system if not already done
1032 737 : if (!_evaluable_elems.count(sys.number()))
1033 : {
1034 465 : auto & [candidate_elems, candidate_elem_ids] = _evaluable_elems[sys.number()];
1035 465 : const auto & dof_map = sys.get_dof_map();
1036 465 : std::vector<dof_id_type> elem_dofs;
1037 465 : auto vn = sys.variable_number(static_cast<std::string>(var_name));
1038 29718 : for (const auto elem : *_mesh.getActiveLocalElementRange())
1039 : {
1040 29253 : if (std::find(_reinitialized_elems.begin(), _reinitialized_elems.end(), elem->id()) !=
1041 58506 : _reinitialized_elems.end())
1042 2771 : continue; // Skip elements that were reinitialized
1043 :
1044 26482 : dof_map.dof_indices(elem, elem_dofs, vn);
1045 26482 : if (!elem_dofs.empty())
1046 : {
1047 9790 : candidate_elems.insert(elem);
1048 9790 : candidate_elem_ids.push_back(elem->id());
1049 : }
1050 : }
1051 465 : }
1052 737 : auto & [candidate_elems, candidate_elem_ids] = _evaluable_elems[sys.number()];
1053 :
1054 : // Now we gather patch elements based on the reinit strategy
1055 737 : auto & patch_elems = _patch_elem_ids[var_name];
1056 :
1057 737 : switch (reinit_strategy)
1058 : {
1059 567 : case ReinitStrategy::POLYNOMIAL_NEIGHBOR:
1060 : {
1061 11104 : auto has_neighbor_in_reinit_elems = [&](const Elem * elem) -> bool
1062 : {
1063 61482 : for (const auto & node : elem->node_ref_range())
1064 287881 : for (const auto & neigh_id : _mesh.nodeToElemMap().at(node.id()))
1065 : // here we need to use _global_reinitialized_elems gathering from all processors
1066 237503 : if (_semi_local_reinitialized_elems.count(neigh_id))
1067 2704 : return true;
1068 8400 : return false;
1069 567 : };
1070 : // Loop over all candidate elements, for each element, if any of its point neighbor belongs
1071 : // to the reinitialized elements, we will include that element in the patch element set.
1072 11671 : for (const auto * elem : candidate_elems)
1073 11104 : if (has_neighbor_in_reinit_elems(elem))
1074 2704 : patch_elems.push_back(elem->id());
1075 567 : break;
1076 : }
1077 102 : case ReinitStrategy::POLYNOMIAL_WHOLE:
1078 : {
1079 : // This is simple: all candidate elements are patch elements
1080 102 : patch_elems = candidate_elem_ids;
1081 102 : break;
1082 : }
1083 68 : case ReinitStrategy::POLYNOMIAL_NEARBY:
1084 : {
1085 68 : std::vector<Point> kd_points;
1086 68 : std::vector<dof_id_type> global_candidate_elem_ids;
1087 :
1088 68 : if (_mesh.isDistributedMesh())
1089 : {
1090 12 : std::vector<std::pair<Point, dof_id_type>> pts_ids(candidate_elem_ids.size());
1091 116 : for (std::size_t i = 0; i < candidate_elem_ids.size(); ++i)
1092 208 : pts_ids[i] = {_mesh.elemPtr(candidate_elem_ids[i])->vertex_average(),
1093 208 : candidate_elem_ids[i]};
1094 12 : _mesh.comm().allgather(pts_ids);
1095 220 : for (const auto & [pt, id] : pts_ids)
1096 : {
1097 208 : kd_points.push_back(pt);
1098 208 : global_candidate_elem_ids.push_back(id);
1099 : }
1100 12 : }
1101 : else
1102 : {
1103 56 : _mesh.comm().allgather(candidate_elem_ids);
1104 56 : global_candidate_elem_ids = candidate_elem_ids;
1105 1036 : for (const auto & id : candidate_elem_ids)
1106 980 : kd_points.push_back(_mesh.elemPtr(id)->vertex_average());
1107 : }
1108 :
1109 68 : const auto kd_tree = std::make_unique<KDTree>(kd_points, _leaf_max_size);
1110 :
1111 68 : std::vector<nanoflann::ResultItem<std::size_t, Real>> query_result;
1112 322 : for (const auto & elem_id : _reinitialized_elems)
1113 : {
1114 254 : const Point & centroid = _mesh.elemPtr(elem_id)->vertex_average();
1115 254 : kd_tree->radiusSearch(centroid, _nearby_distance_threshold, query_result);
1116 1344 : for (const auto & [qid, dist] : query_result)
1117 1090 : patch_elems.push_back(global_candidate_elem_ids[qid]);
1118 : }
1119 68 : break;
1120 68 : }
1121 0 : case ReinitStrategy::NONE:
1122 : // do nothing
1123 0 : break;
1124 0 : default:
1125 0 : mooseError("Unknown reinitialization strategy");
1126 : break;
1127 : }
1128 :
1129 : // every processor should have the same patch elements to do the polynomial extrapolation,
1130 : // so we gather them across all processors
1131 737 : _mesh.comm().allgather(patch_elems);
1132 :
1133 : // Remove duplicates from the patch elements (espcially important for POLYNOMIAL_NEARBY)
1134 737 : std::sort(patch_elems.begin(), patch_elems.end());
1135 737 : patch_elems.erase(std::unique(patch_elems.begin(), patch_elems.end()), patch_elems.end());
1136 737 : }
1137 :
1138 : void
1139 737 : ElementSubdomainModifierBase::extrapolatePolynomial(const VariableName & var_name)
1140 : {
1141 : const auto & coef =
1142 737 : _pr[_var_name_to_pr_idx[var_name]]->getCachedCoefficients(_patch_elem_ids[var_name]);
1143 :
1144 737 : const unsigned dim = _mesh.dimension();
1145 :
1146 737 : libMesh::Parameters function_parameters;
1147 :
1148 737 : const auto & multi_index = _pr[_var_name_to_pr_idx[var_name]]->multiIndex();
1149 :
1150 1474 : function_parameters.set<std::vector<std::vector<unsigned int>>>("multi_index") = multi_index;
1151 :
1152 737 : std::vector<Real> coef_vec(coef.size());
1153 3109 : for (auto i = 0; i < coef.size(); ++i)
1154 2372 : coef_vec[i] = coef(i);
1155 :
1156 737 : function_parameters.set<std::vector<Real>>("multi_index_coefficients") = coef_vec;
1157 1474 : function_parameters.set<unsigned int>("dimension_for_projection") = dim;
1158 :
1159 : // Define projection function
1160 9271 : auto poly_func = [](const Point & p,
1161 : const libMesh::Parameters & parameters,
1162 : const std::string &,
1163 : const std::string &) -> libMesh::Number
1164 : {
1165 : const auto & multi_index =
1166 9271 : parameters.get<std::vector<std::vector<unsigned int>>>("multi_index");
1167 9271 : const auto & coeffs = parameters.get<std::vector<Real>>("multi_index_coefficients");
1168 :
1169 9271 : Real val = 0.0;
1170 :
1171 57321 : for (unsigned int r = 0; r < multi_index.size(); r++)
1172 : {
1173 48050 : Real monomial = 1.0;
1174 173060 : for (unsigned int d = 0; d < multi_index[r].size(); d++)
1175 : {
1176 125010 : const auto power = multi_index[r][d];
1177 125010 : if (power == 0)
1178 77558 : continue;
1179 :
1180 47452 : monomial *= std::pow(p(d), power);
1181 : }
1182 48050 : val += coeffs[r] * monomial;
1183 : }
1184 :
1185 9271 : return val;
1186 : };
1187 :
1188 : // Define gradient
1189 0 : auto poly_func_grad = [](const Point & p,
1190 : const libMesh::Parameters & parameters,
1191 : const std::string &,
1192 : const std::string &) -> libMesh::Gradient
1193 : {
1194 0 : const unsigned int dim = parameters.get<unsigned int>("dimension_for_projection");
1195 :
1196 : const auto & multi_index =
1197 0 : parameters.get<std::vector<std::vector<unsigned int>>>("multi_index");
1198 0 : const auto & coeffs = parameters.get<std::vector<Real>>("multi_index_coefficients");
1199 :
1200 0 : libMesh::Gradient grad; // Zero-initialized
1201 :
1202 0 : for (unsigned int r = 0; r < multi_index.size(); ++r)
1203 : {
1204 0 : const auto & powers = multi_index[r];
1205 0 : const Real coef = coeffs[r];
1206 :
1207 0 : for (unsigned int d = 0; d < dim; ++d) // Loop over dimension
1208 : {
1209 0 : const auto power_d = powers[d];
1210 0 : if (power_d == 0)
1211 0 : continue;
1212 :
1213 : // Compute partial derivative in direction d
1214 0 : Real partial = coef * power_d;
1215 :
1216 0 : for (unsigned int i = 0; i < powers.size(); ++i)
1217 : {
1218 0 : if (i == d)
1219 : {
1220 0 : if (powers[i] > 1)
1221 0 : partial *= std::pow(p(i), powers[i] - 1); // reduce power by 1
1222 : }
1223 : else
1224 : {
1225 0 : if (powers[i] > 0)
1226 0 : partial *= std::pow(p(i), powers[i]); // full power
1227 : }
1228 : }
1229 :
1230 0 : grad(d) += partial;
1231 : }
1232 : }
1233 :
1234 0 : return grad;
1235 : };
1236 :
1237 2211 : std::vector<VariableName> var_names = {var_name};
1238 737 : _fe_problem.projectFunctionOnCustomRange(
1239 : reinitializedElemRange(), poly_func, poly_func_grad, function_parameters, var_names);
1240 1474 : }
|
