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/parallel_algebra.h"
15 : #include "libmesh/parallel.h"
16 : #include "libmesh/dof_map.h"
17 : #include "libmesh/remote_elem.h"
18 : #include "libmesh/parallel_ghost_sync.h"
19 : #include "libmesh/petsc_vector.h"
20 :
21 : InputParameters
22 43678 : ElementSubdomainModifierBase::validParams()
23 : {
24 43678 : InputParameters params = ElementUserObject::validParams();
25 :
26 : // ESMs only operate on the undisplaced mesh to gather subdomain and sideset information. This
27 : // information is used to modify both the undisplaced and the displaced meshes. It is the
28 : // developer's responsibility to make sure the element IDs and sideset info are consistent across
29 : // both meshes.
30 43678 : params.set<bool>("use_displaced_mesh") = false;
31 43678 : params.suppressParameter<bool>("use_displaced_mesh");
32 :
33 43678 : params.addDeprecatedParam<BoundaryName>(
34 : "moving_boundary_name",
35 : "Name of the moving boundary",
36 : "This has been replaced by 'moving_boundaries' and 'moving_boundary_subdomain_pairs'.");
37 43678 : params.addDeprecatedParam<BoundaryName>(
38 : "complement_moving_boundary_name",
39 : "Name of the moving boundary on the complement subdomain(s)",
40 : "This has been replaced by 'moving_boundaries' and 'moving_boundary_subdomain_pairs'.");
41 131034 : params.addDeprecatedParam<bool>(
42 : "apply_initial_conditions",
43 87356 : true,
44 : "Whether to apply initial conditions on the moved nodes and elements",
45 : "This has been replaced by 'reinitialize_subdomains'");
46 :
47 43678 : params.addParam<std::vector<BoundaryName>>(
48 : "moving_boundaries",
49 : {},
50 : "Moving boundaries between subdomains. These boundaries (both sidesets and nodesets) will be "
51 : "updated for elements that change subdomain. The subdomains that each moving "
52 : "boundary lies between shall be specified using the parameter "
53 : "'moving_boundary_subdomain_pairs'. If one boundary and multiple subdomain pairs are "
54 : "specified, then it is assumed that the pairs all apply to the boundary. A boundary will be "
55 : "created on the mesh if it does not already exist.");
56 43678 : params.addParam<std::vector<std::vector<SubdomainName>>>(
57 : "moving_boundary_subdomain_pairs",
58 : {},
59 : "The subdomain pairs associated with each moving boundary. For each pair of subdomains, only "
60 : "the element side from the first subdomain will be added to the moving boundary, i.e., the "
61 : "side normal is pointing from the first subdomain to the second subdomain. The pairs shall "
62 : "be delimited by ';'. If a pair only has one subdomain, the moving boundary is associated "
63 : "with the subdomain's external boundary, i.e., when the elements have no neighboring "
64 : "elements.");
65 :
66 131034 : params.addParam<std::vector<SubdomainName>>(
67 : "reinitialize_subdomains",
68 : {"ANY_BLOCK_ID"},
69 : "By default, any element which changes subdomain is reinitialized. If a list of subdomains "
70 : "(IDs or names) is provided, then only elements whose new subdomain is in the list will be "
71 : "reinitialized. If an empty list is set, then no elements will be reinitialized.");
72 131034 : params.addParam<bool>(
73 : "old_subdomain_reinitialized",
74 87356 : true,
75 : "This parameter must be set with a non-empty list in 'reinitialize_subdomains'. When set to "
76 : "the default true, the element's old subdomain is not considered when determining if an "
77 : "element should be reinitialized. If set to false, only elements whose old subdomain was not "
78 : "in 'reinitialize_subdomains' are reinitialized. ");
79 :
80 43678 : params.registerBase("MeshModifier");
81 :
82 43678 : return params;
83 43678 : }
84 :
85 460 : ElementSubdomainModifierBase::ElementSubdomainModifierBase(const InputParameters & parameters)
86 : : ElementUserObject(parameters),
87 460 : _displaced_problem(_fe_problem.getDisplacedProblem().get()),
88 460 : _displaced_mesh(_displaced_problem ? &_displaced_problem->mesh() : nullptr),
89 920 : _old_subdomain_reinitialized(getParam<bool>("old_subdomain_reinitialized"))
90 : {
91 1380 : if (isParamSetByUser("moving_boundary_name") ||
92 920 : isParamSetByUser("complement_moving_boundary_name"))
93 0 : mooseError(
94 : "'moving_boundary_name' and 'complement_moving_boundary_name' have been replaced by "
95 : "'moving_boundaries' and 'moving_boundary_subdomain_pairs'. See the documentation in "
96 : "https://mooseframework.inl.gov/blackbear/source/userobjects/"
97 : "ElementSubdomainModifier.html for more information. "
98 : "Additionally, SidesetAroundSubdomainUpdater can now be used to update boundaries "
99 : "that are defined around a subdomain, or between two subdomains.");
100 460 : }
101 :
102 : void
103 455 : ElementSubdomainModifierBase::initialSetup()
104 : {
105 : // When 'apply_initial_conditions' is fully deprecated, change this to 'const' vector
106 : std::vector<SubdomainName> subdomain_names_to_reinitialize =
107 455 : getParam<std::vector<SubdomainName>>("reinitialize_subdomains");
108 :
109 : // When 'apply_initial_conditions' is fully deprecated, remove this block
110 455 : if (isParamSetByUser("apply_initial_conditions") && getParam<bool>("apply_initial_conditions"))
111 0 : subdomain_names_to_reinitialize = {"ANY_BLOCK_ID"};
112 455 : else if (isParamSetByUser("apply_initial_conditions"))
113 0 : subdomain_names_to_reinitialize = {};
114 :
115 455 : if (std::find(subdomain_names_to_reinitialize.begin(),
116 : subdomain_names_to_reinitialize.end(),
117 910 : "ANY_BLOCK_ID") != subdomain_names_to_reinitialize.end())
118 364 : _subdomain_ids_to_reinitialize.push_back(Moose::ANY_BLOCK_ID);
119 : else
120 91 : _subdomain_ids_to_reinitialize = _mesh.getSubdomainIDs(subdomain_names_to_reinitialize);
121 :
122 : std::set<SubdomainID> set_subdomain_ids_to_reinitialize(_subdomain_ids_to_reinitialize.begin(),
123 455 : _subdomain_ids_to_reinitialize.end());
124 :
125 507 : if (_old_subdomain_reinitialized == false &&
126 26 : (std::find(_subdomain_ids_to_reinitialize.begin(),
127 : _subdomain_ids_to_reinitialize.end(),
128 507 : Moose::ANY_BLOCK_ID) != _subdomain_ids_to_reinitialize.end() ||
129 26 : set_subdomain_ids_to_reinitialize == _mesh.meshSubdomains()))
130 0 : paramError("old_subdomain_reinitialized",
131 : "'old_subdomain_reinitialized' can only be set to false if "
132 : "reinitialize_subdomains does "
133 : "not cover the whole model, otherwise no elements will be reinitialized as it is "
134 : "impossible for an element's old subdomain to not be in the list.");
135 455 : else if (_old_subdomain_reinitialized == false && _subdomain_ids_to_reinitialize.empty())
136 0 : paramError("old_subdomain_reinitialized",
137 : "'old_subdomain_reinitialized' can only be set to false if "
138 : "reinitialize_subdomains is set to a non-empty list of subdomains, otherwise no "
139 : "elements will be reinitialized, as it is impossible for an element's new subdomain "
140 : "to be in the list.");
141 :
142 455 : auto bnd_names = getParam<std::vector<BoundaryName>>("moving_boundaries");
143 455 : auto bnd_ids = _mesh.getBoundaryIDs(bnd_names, true);
144 : const auto bnd_subdomains =
145 455 : getParam<std::vector<std::vector<SubdomainName>>>("moving_boundary_subdomain_pairs");
146 :
147 455 : if (bnd_names.size() == 1 && bnd_subdomains.size() > 1)
148 : {
149 13 : bnd_names.insert(bnd_names.end(), bnd_subdomains.size() - 1, bnd_names[0]);
150 13 : bnd_ids.insert(bnd_ids.end(), bnd_subdomains.size() - 1, bnd_ids[0]);
151 : }
152 442 : else if (bnd_names.size() != bnd_subdomains.size())
153 0 : paramError("moving_boundary_subdomain_pairs",
154 : "Each moving boundary must correspond to a pair of subdomains. ",
155 : bnd_names.size(),
156 : " boundaries are specified by the parameter 'moving_boundaries', while ",
157 : bnd_subdomains.size(),
158 : " subdomain pairs are provided. Alternatively, if one boundary and multiple "
159 : "subdomain pairs are provided, then the subdomain pairs all apply to one boundary.");
160 :
161 616 : for (auto i : index_range(bnd_names))
162 : {
163 161 : _moving_boundary_names[bnd_ids[i]] = bnd_names[i];
164 :
165 161 : if (bnd_subdomains[i].size() == 2)
166 148 : _moving_boundaries[{_mesh.getSubdomainID(bnd_subdomains[i][0]),
167 296 : _mesh.getSubdomainID(bnd_subdomains[i][1])}] = bnd_ids[i];
168 13 : else if (bnd_subdomains[i].size() == 1)
169 13 : _moving_boundaries[{_mesh.getSubdomainID(bnd_subdomains[i][0]), Moose::INVALID_BLOCK_ID}] =
170 13 : bnd_ids[i];
171 : else
172 0 : paramError("moving_boundary_subdomain_pairs",
173 : "Each subdomain pair must contain 1 or 2 subdomain names, but ",
174 0 : bnd_subdomains[i].size(),
175 : " are given.");
176 : }
177 455 : }
178 :
179 : void
180 2663 : ElementSubdomainModifierBase::modify(
181 : const std::unordered_map<dof_id_type, std::pair<SubdomainID, SubdomainID>> & moved_elems)
182 : {
183 :
184 : // If nothing need to change, just return.
185 : // This will skip all mesh changes, and so no adaptivity mesh files will be written.
186 2663 : auto n_moved_elem = moved_elems.size();
187 2663 : gatherSum(n_moved_elem);
188 2663 : if (n_moved_elem == 0)
189 1291 : return;
190 :
191 : // Create moving boundaries on the undisplaced and displaced meshes
192 : //
193 : // Note: We do this _everytime_ because previous execution might have removed the sidesets and
194 : // nodesets. Most of the moving boundary algorithms below assume that the moving sidesets and
195 : // nodesets already exist on the mesh.
196 1372 : createMovingBoundaries(_mesh);
197 1372 : if (_displaced_mesh)
198 44 : createMovingBoundaries(*_displaced_mesh);
199 :
200 : // This has to be done _before_ subdomain changes are applied
201 1372 : findReinitializedElemsAndNodes(moved_elems);
202 :
203 : // Apply cached subdomain changes
204 1372 : applySubdomainChanges(moved_elems, _mesh);
205 1372 : if (_displaced_mesh)
206 44 : applySubdomainChanges(moved_elems, *_displaced_mesh);
207 :
208 : // Update moving boundaries
209 1372 : gatherMovingBoundaryChanges(moved_elems);
210 1372 : applyMovingBoundaryChanges(_mesh);
211 1372 : if (_displaced_mesh)
212 44 : applyMovingBoundaryChanges(*_displaced_mesh);
213 :
214 : // Reinit equation systems
215 1372 : _fe_problem.meshChanged(
216 : /*intermediate_change=*/false, /*contract_mesh=*/false, /*clean_refinement_flags=*/false);
217 :
218 : // Initialize solution and stateful material properties
219 1372 : applyIC(/*displaced=*/false);
220 1372 : if (_fe_problem.getMaterialWarehouse().hasActiveObjects(0))
221 48 : initElementStatefulProps(/*displaced=*/false);
222 :
223 1372 : if (_displaced_mesh)
224 : {
225 44 : applyIC(/*displaced=*/true);
226 44 : if (_fe_problem.getMaterialWarehouse().hasActiveObjects(0))
227 0 : initElementStatefulProps(/*displaced=*/true);
228 : }
229 : }
230 :
231 : void
232 1416 : ElementSubdomainModifierBase::createMovingBoundaries(MooseMesh & mesh)
233 : {
234 1416 : auto & bnd_info = mesh.getMesh().get_boundary_info();
235 2112 : for (const auto & [bnd_id, bnd_name] : _moving_boundary_names)
236 : {
237 696 : bnd_info.sideset_name(bnd_id) = bnd_name;
238 696 : bnd_info.nodeset_name(bnd_id) = bnd_name;
239 : }
240 1416 : }
241 :
242 : void
243 1416 : ElementSubdomainModifierBase::applySubdomainChanges(
244 : const std::unordered_map<dof_id_type, std::pair<SubdomainID, SubdomainID>> & moved_elems,
245 : MooseMesh & mesh)
246 : {
247 28053 : for (const auto & [elem_id, subdomain] : moved_elems)
248 : {
249 : // Change the element's subdomain ID
250 26637 : auto elem = mesh.elemPtr(elem_id);
251 26637 : const auto & [from, to] = subdomain;
252 : mooseAssert(elem->subdomain_id() == from, "Inconsistent element subdomain ID.");
253 26637 : elem->subdomain_id() = to;
254 :
255 : // Change the ancestors' (if any) subdomain ID
256 26637 : setAncestorsSubdomainIDs(elem, to);
257 : }
258 :
259 : // Synchronize ghost element subdomain changes
260 1416 : libMesh::SyncSubdomainIds sync(mesh.getMesh());
261 1416 : Parallel::sync_dofobject_data_by_id(
262 2832 : mesh.getMesh().comm(), mesh.getMesh().elements_begin(), mesh.getMesh().elements_end(), sync);
263 1416 : }
264 :
265 : void
266 26637 : ElementSubdomainModifierBase::setAncestorsSubdomainIDs(Elem * elem, const SubdomainID subdomain_id)
267 : {
268 26637 : auto curr_elem = elem;
269 :
270 34697 : for (unsigned int i = curr_elem->level(); i > 0; --i)
271 : {
272 : // Change the parent's subdomain
273 8060 : curr_elem = curr_elem->parent();
274 8060 : curr_elem->subdomain_id() = subdomain_id;
275 : }
276 26637 : }
277 :
278 : void
279 1372 : ElementSubdomainModifierBase::gatherMovingBoundaryChanges(
280 : const std::unordered_map<dof_id_type, std::pair<SubdomainID, SubdomainID>> & moved_elems)
281 : {
282 : // Clear moving boundary changes from last execution
283 1372 : _add_element_sides.clear();
284 1372 : _add_neighbor_sides.clear();
285 1372 : _remove_element_sides.clear();
286 1372 : _remove_neighbor_sides.clear();
287 :
288 1372 : const auto & sidesets = _mesh.getMesh().get_boundary_info().get_sideset_map();
289 :
290 27561 : for (const auto & [elem_id, subdomain_assignment] : moved_elems)
291 : {
292 26189 : auto elem = _mesh.elemPtr(elem_id);
293 :
294 : // The existing moving boundaries on the element side should be removed
295 35235 : for (auto itr = sidesets.lower_bound(elem); itr != sidesets.upper_bound(elem); itr++)
296 9046 : if (_moving_boundary_names.count(itr->second.second))
297 3008 : _remove_element_sides[elem->id()].emplace(itr->second.first, itr->second.second);
298 :
299 135313 : for (auto side : elem->side_index_range())
300 : {
301 109124 : auto neigh = elem->neighbor_ptr(side);
302 :
303 : // Don't mess with remote element neighbor
304 109124 : if (neigh && neigh == libMesh::remote_elem)
305 0 : continue;
306 : // If neighbor doesn't exist
307 109124 : else if (!neigh)
308 6730 : gatherMovingBoundaryChangesHelper(elem, side, nullptr, 0);
309 : // If neighbor exists
310 : else
311 : {
312 102394 : auto neigh_side = neigh->which_neighbor_am_i(elem);
313 :
314 102394 : if (neigh->active())
315 99624 : gatherMovingBoundaryChangesHelper(elem, side, neigh, neigh_side);
316 : else
317 : {
318 : // Find the active neighbors of the element
319 2770 : std::vector<const Elem *> active_neighs;
320 : // Neighbor has active children, they are neighbors of the element along that side
321 : mooseAssert(!neigh->subactive(),
322 : "The case where the active neighbor is an ancestor of this neighbor is not "
323 : "handled at this time.");
324 2770 : neigh->active_family_tree_by_neighbor(active_neighs, elem);
325 :
326 8310 : for (auto active_neigh : active_neighs)
327 5540 : gatherMovingBoundaryChangesHelper(elem, side, active_neigh, neigh_side);
328 2770 : }
329 : }
330 : }
331 : }
332 1372 : }
333 :
334 : void
335 111894 : ElementSubdomainModifierBase::gatherMovingBoundaryChangesHelper(const Elem * elem,
336 : unsigned short side,
337 : const Elem * neigh,
338 : unsigned short neigh_side)
339 : {
340 111894 : const auto & sidesets = _mesh.getMesh().get_boundary_info().get_sideset_map();
341 :
342 : // Detect element side change
343 111894 : SubdomainPair subdomain_pair = {elem->subdomain_id(),
344 111894 : neigh ? neigh->subdomain_id() : Moose::INVALID_BLOCK_ID};
345 111894 : if (_moving_boundaries.count(subdomain_pair))
346 8158 : _add_element_sides[elem->id()].emplace(side, _moving_boundaries.at(subdomain_pair));
347 :
348 111894 : if (neigh)
349 : {
350 : // The existing moving boundaries on the neighbor side should be removed
351 139659 : for (auto itr = sidesets.lower_bound(neigh); itr != sidesets.upper_bound(neigh); itr++)
352 34495 : if (itr->second.first == neigh_side && _moving_boundary_names.count(itr->second.second))
353 9536 : _remove_neighbor_sides[neigh->id()].emplace(itr->second.first, itr->second.second);
354 :
355 : // Detect neighbor side change (by reversing the subdomain pair)
356 105164 : subdomain_pair = {subdomain_pair.second, subdomain_pair.first};
357 105164 : if (_moving_boundaries.count(subdomain_pair))
358 2752 : _add_neighbor_sides[neigh->id()].emplace(neigh_side, _moving_boundaries.at(subdomain_pair));
359 : }
360 111894 : }
361 :
362 : void
363 1416 : ElementSubdomainModifierBase::applyMovingBoundaryChanges(MooseMesh & mesh)
364 : {
365 1416 : auto & bnd_info = mesh.getMesh().get_boundary_info();
366 :
367 : // Remove all boundary nodes from the previous moving boundaries
368 1416 : auto nodesets = bnd_info.get_nodeset_map();
369 152252 : for (const auto & [node_id, bnd] : nodesets)
370 150836 : if (_moving_boundary_names.count(bnd))
371 26635 : bnd_info.remove_node(node_id, bnd);
372 :
373 : // Keep track of ghost element changes
374 : std::unordered_map<processor_id_type,
375 : std::vector<std::tuple<dof_id_type, unsigned short, BoundaryID>>>
376 1416 : add_ghost_sides, remove_ghost_sides;
377 :
378 : // Remove element sides from moving boundaries
379 4344 : for (const auto & [elem_id, sides] : _remove_element_sides)
380 5936 : for (const auto & [side, bnd] : sides)
381 3008 : bnd_info.remove_side(mesh.elemPtr(elem_id), side, bnd);
382 :
383 : // Remove neighbor sides from moving boundaries
384 9947 : for (const auto & [elem_id, sides] : _remove_neighbor_sides)
385 : {
386 8531 : auto elem = mesh.elemPtr(elem_id);
387 18067 : for (const auto & [side, bnd] : sides)
388 : {
389 9536 : bnd_info.remove_side(elem, side, bnd);
390 : // Keep track of changes to ghosted elements
391 9536 : if (elem->processor_id() != processor_id())
392 83 : remove_ghost_sides[elem->processor_id()].push_back({elem_id, side, bnd});
393 : }
394 : }
395 :
396 : // Add element sides to moving boundaries
397 7450 : for (const auto & [elem_id, sides] : _add_element_sides)
398 14150 : for (const auto & [side, bnd] : sides)
399 8116 : bnd_info.add_side(mesh.elemPtr(elem_id), side, bnd);
400 :
401 : // Add neighbor sides to moving boundaries
402 3534 : for (const auto & [elem_id, sides] : _add_neighbor_sides)
403 : {
404 2118 : auto elem = mesh.elemPtr(elem_id);
405 4314 : for (const auto & [side, bnd] : sides)
406 : {
407 2196 : bnd_info.add_side(elem, side, bnd);
408 : // Keep track of changes to ghosted elements
409 2196 : if (elem->processor_id() != processor_id())
410 13 : add_ghost_sides[elem->processor_id()].push_back({elem_id, side, bnd});
411 : }
412 : }
413 :
414 1416 : Parallel::push_parallel_vector_data(
415 1416 : bnd_info.comm(),
416 : add_ghost_sides,
417 1416 : [&mesh,
418 : &bnd_info](processor_id_type,
419 13 : const std::vector<std::tuple<dof_id_type, unsigned short, BoundaryID>> & received)
420 : {
421 22 : for (const auto & [elem_id, side, bnd] : received)
422 13 : bnd_info.add_side(mesh.elemPtr(elem_id), side, bnd);
423 9 : });
424 :
425 1416 : Parallel::push_parallel_vector_data(
426 1416 : bnd_info.comm(),
427 : remove_ghost_sides,
428 1416 : [&mesh,
429 : &bnd_info](processor_id_type,
430 83 : const std::vector<std::tuple<dof_id_type, unsigned short, BoundaryID>> & received)
431 : {
432 114 : for (const auto & [elem_id, side, bnd] : received)
433 83 : bnd_info.remove_side(mesh.elemPtr(elem_id), side, bnd);
434 31 : });
435 :
436 1416 : bnd_info.parallel_sync_side_ids();
437 1416 : bnd_info.parallel_sync_node_ids();
438 1416 : mesh.update();
439 1416 : }
440 :
441 : void
442 2668 : ElementSubdomainModifierBase::meshChanged()
443 : {
444 : // Clear cached ranges
445 2668 : _reinitialized_elem_range.reset();
446 2668 : _reinitialized_displaced_elem_range.reset();
447 2668 : _reinitialized_bnd_node_range.reset();
448 2668 : _reinitialized_displaced_bnd_node_range.reset();
449 :
450 2668 : updateAMRMovingBoundary(_mesh);
451 2668 : if (_displaced_mesh)
452 48 : updateAMRMovingBoundary(*_displaced_mesh);
453 2668 : }
454 :
455 : void
456 2716 : ElementSubdomainModifierBase::updateAMRMovingBoundary(MooseMesh & mesh)
457 : {
458 2716 : auto & bnd_info = mesh.getMesh().get_boundary_info();
459 2716 : auto sidesets = bnd_info.get_sideset_map();
460 225664 : for (const auto & i : sidesets)
461 : {
462 222948 : auto elem = i.first;
463 222948 : auto side = i.second.first;
464 222948 : auto bnd = i.second.second;
465 222948 : if (_moving_boundary_names.count(bnd) && !elem->active())
466 : {
467 5334 : bnd_info.remove_side(elem, side, bnd);
468 :
469 5334 : std::vector<const Elem *> elem_family;
470 5334 : elem->active_family_tree_by_side(elem_family, side);
471 17029 : for (auto felem : elem_family)
472 11695 : bnd_info.add_side(felem, side, bnd);
473 5334 : }
474 : }
475 :
476 2716 : bnd_info.parallel_sync_side_ids();
477 2716 : bnd_info.parallel_sync_node_ids();
478 2716 : }
479 :
480 : void
481 1372 : ElementSubdomainModifierBase::findReinitializedElemsAndNodes(
482 : const std::unordered_map<dof_id_type, std::pair<SubdomainID, SubdomainID>> & moved_elems)
483 : {
484 : // Clear cached element reinitialization data
485 1372 : _reinitialized_elems.clear();
486 1372 : _reinitialized_nodes.clear();
487 :
488 27561 : for (const auto & [elem_id, subdomain] : moved_elems)
489 : {
490 : mooseAssert(_mesh.elemPtr(elem_id)->active(), "Moved elements should be active");
491 : // Default: any element that changes subdomain is reinitialized
492 26189 : if (std::find(_subdomain_ids_to_reinitialize.begin(),
493 : _subdomain_ids_to_reinitialize.end(),
494 52378 : Moose::ANY_BLOCK_ID) != _subdomain_ids_to_reinitialize.end())
495 21237 : _reinitialized_elems.insert(elem_id);
496 : else // Reinitialize if new subdomain is in list of subdomains to be reinitialized
497 : {
498 4952 : const auto & [from, to] = subdomain;
499 4952 : if (subdomainIsReinitialized(to) && _old_subdomain_reinitialized)
500 80 : _reinitialized_elems.insert(elem_id);
501 : // Only reinitialize if original subdomain is not in list of subdomains
502 6832 : else if (subdomainIsReinitialized(to) && !_old_subdomain_reinitialized &&
503 1960 : !subdomainIsReinitialized(from))
504 1120 : _reinitialized_elems.insert(elem_id);
505 : else // New subdomain is not in list of subdomains
506 3752 : continue;
507 : }
508 :
509 22437 : const auto elem = _mesh.elemPtr(elem_id);
510 120921 : for (unsigned int i = 0; i < elem->n_nodes(); ++i)
511 98484 : if (nodeIsNewlyReinitialized(elem->node_id(i)))
512 3936 : _reinitialized_nodes.insert(elem->node_id(i));
513 : }
514 1372 : }
515 :
516 : bool
517 121228 : ElementSubdomainModifierBase::subdomainIsReinitialized(SubdomainID id) const
518 : {
519 : // Default: any element that changes subdomain is reinitialized
520 121228 : if (std::find(_subdomain_ids_to_reinitialize.begin(),
521 : _subdomain_ids_to_reinitialize.end(),
522 242456 : Moose::ANY_BLOCK_ID) != _subdomain_ids_to_reinitialize.end())
523 93684 : return true;
524 :
525 : // Is subdomain in list of subdomains to be reinitialized
526 27544 : return std::find(_subdomain_ids_to_reinitialize.begin(),
527 : _subdomain_ids_to_reinitialize.end(),
528 55088 : id) != _subdomain_ids_to_reinitialize.end();
529 : }
530 :
531 : bool
532 98484 : ElementSubdomainModifierBase::nodeIsNewlyReinitialized(dof_id_type node_id) const
533 : {
534 : // If any of the node neighbors are already reinitialized, then the node is NOT newly
535 : // reinitialized.
536 113380 : for (auto neighbor_elem_id : _mesh.nodeToElemMap().at(node_id))
537 109444 : if (subdomainIsReinitialized(_mesh.elemPtr(neighbor_elem_id)->subdomain_id()))
538 94548 : return false;
539 3936 : return true;
540 : }
541 :
542 : void
543 1416 : ElementSubdomainModifierBase::applyIC(bool displaced)
544 : {
545 1416 : _fe_problem.projectInitialConditionOnCustomRange(reinitializedElemRange(displaced),
546 : reinitializedBndNodeRange(displaced));
547 :
548 : mooseAssert(_fe_problem.numSolverSystems() < 2,
549 : "This code was written for a single nonlinear system");
550 : // Set old and older solutions on the reinitialized dofs to the reinitialized values
551 1416 : setOldAndOlderSolutions(_fe_problem.getNonlinearSystemBase(_sys.number()),
552 : reinitializedElemRange(displaced),
553 : reinitializedBndNodeRange(displaced));
554 1416 : setOldAndOlderSolutions(_fe_problem.getAuxiliarySystem(),
555 : reinitializedElemRange(displaced),
556 : reinitializedBndNodeRange(displaced));
557 :
558 : // Note: Need method to handle solve failures at timesteps where subdomain changes. The old
559 : // solutions are now set to the reinitialized values. Does this impact restoring solutions
560 1416 : }
561 :
562 : void
563 48 : ElementSubdomainModifierBase::initElementStatefulProps(bool displaced)
564 : {
565 48 : _fe_problem.initElementStatefulProps(reinitializedElemRange(displaced), /*threaded=*/true);
566 48 : }
567 :
568 : ConstElemRange &
569 5712 : ElementSubdomainModifierBase::reinitializedElemRange(bool displaced)
570 : {
571 5712 : if (!displaced && _reinitialized_elem_range)
572 4164 : return *_reinitialized_elem_range.get();
573 :
574 1548 : if (displaced && _reinitialized_displaced_elem_range)
575 132 : return *_reinitialized_displaced_elem_range.get();
576 :
577 : // Create a vector of the newly reinitialized elements
578 1416 : std::vector<Elem *> elems;
579 24301 : for (auto elem_id : _reinitialized_elems)
580 22885 : elems.push_back(displaced ? _displaced_mesh->elemPtr(elem_id) : _mesh.elemPtr(elem_id));
581 :
582 : // Make some fake element iterators defining this vector of elements
583 1416 : Elem * const * elem_itr_begin = const_cast<Elem * const *>(elems.data());
584 1416 : Elem * const * elem_itr_end = elem_itr_begin + elems.size();
585 :
586 : const auto elems_begin = MeshBase::const_element_iterator(
587 1416 : elem_itr_begin, elem_itr_end, Predicates::NotNull<Elem * const *>());
588 : const auto elems_end = MeshBase::const_element_iterator(
589 1416 : elem_itr_end, elem_itr_end, Predicates::NotNull<Elem * const *>());
590 :
591 1416 : if (!displaced)
592 1372 : _reinitialized_elem_range = std::make_unique<ConstElemRange>(elems_begin, elems_end);
593 : else
594 44 : _reinitialized_displaced_elem_range = std::make_unique<ConstElemRange>(elems_begin, elems_end);
595 :
596 1416 : return reinitializedElemRange(displaced);
597 1416 : }
598 :
599 : ConstBndNodeRange &
600 5664 : ElementSubdomainModifierBase::reinitializedBndNodeRange(bool displaced)
601 : {
602 5664 : if (!displaced && _reinitialized_bnd_node_range)
603 4116 : return *_reinitialized_bnd_node_range.get();
604 :
605 1548 : if (displaced && _reinitialized_displaced_bnd_node_range)
606 132 : return *_reinitialized_displaced_bnd_node_range.get();
607 :
608 : // Create a vector of the newly reinitialized boundary nodes
609 1416 : std::vector<const BndNode *> nodes;
610 : auto bnd_nodes =
611 1416 : displaced ? _displaced_mesh->getBoundaryNodeRange() : _mesh.getBoundaryNodeRange();
612 149367 : for (auto bnd_node : *bnd_nodes)
613 147951 : if (bnd_node->_node)
614 147951 : if (_reinitialized_nodes.count(bnd_node->_node->id()))
615 288 : nodes.push_back(bnd_node);
616 :
617 : // Make some fake node iterators defining this vector of nodes
618 1416 : BndNode * const * bnd_node_itr_begin = const_cast<BndNode * const *>(nodes.data());
619 1416 : BndNode * const * bnd_node_itr_end = bnd_node_itr_begin + nodes.size();
620 :
621 : const auto bnd_nodes_begin = MooseMesh::const_bnd_node_iterator(
622 1416 : bnd_node_itr_begin, bnd_node_itr_end, Predicates::NotNull<const BndNode * const *>());
623 : const auto bnd_nodes_end = MooseMesh::const_bnd_node_iterator(
624 1416 : bnd_node_itr_end, bnd_node_itr_end, Predicates::NotNull<const BndNode * const *>());
625 :
626 1416 : if (!displaced)
627 : _reinitialized_bnd_node_range =
628 1372 : std::make_unique<ConstBndNodeRange>(bnd_nodes_begin, bnd_nodes_end);
629 : else
630 : _reinitialized_displaced_bnd_node_range =
631 44 : std::make_unique<ConstBndNodeRange>(bnd_nodes_begin, bnd_nodes_end);
632 :
633 1416 : return reinitializedBndNodeRange(displaced);
634 1416 : }
635 :
636 : void
637 2832 : ElementSubdomainModifierBase::setOldAndOlderSolutions(SystemBase & sys,
638 : ConstElemRange & elem_range,
639 : ConstBndNodeRange & bnd_node_range)
640 : {
641 : // Don't do anything if this is a steady simulation
642 2832 : if (!sys.hasSolutionState(1))
643 22 : return;
644 :
645 2810 : NumericVector<Number> & current_solution = *sys.system().current_local_solution;
646 2810 : NumericVector<Number> & old_solution = sys.solutionOld();
647 2810 : NumericVector<Number> * older_solution = sys.hasSolutionState(2) ? &sys.solutionOlder() : nullptr;
648 :
649 : // Get dofs for the reinitialized elements and nodes
650 2810 : DofMap & dof_map = sys.dofMap();
651 2810 : std::vector<dof_id_type> dofs;
652 :
653 48244 : for (auto & elem : elem_range)
654 : {
655 45434 : std::vector<dof_id_type> elem_dofs;
656 45434 : dof_map.dof_indices(elem, elem_dofs);
657 45434 : dofs.insert(dofs.end(), elem_dofs.begin(), elem_dofs.end());
658 45434 : }
659 :
660 3386 : for (auto & bnd_node : bnd_node_range)
661 : {
662 576 : std::vector<dof_id_type> bnd_node_dofs;
663 576 : dof_map.dof_indices(bnd_node->_node, bnd_node_dofs);
664 576 : dofs.insert(dofs.end(), bnd_node_dofs.begin(), bnd_node_dofs.end());
665 576 : }
666 :
667 : // Set the old and older solutions to match the reinitialization
668 171782 : for (auto dof : dofs)
669 : {
670 168972 : old_solution.set(dof, current_solution(dof));
671 168972 : if (older_solution)
672 0 : older_solution->set(dof, current_solution(dof));
673 : }
674 :
675 2810 : old_solution.close();
676 2810 : if (older_solution)
677 0 : older_solution->close();
678 2810 : }
|
