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 "ActivateElementsUserObjectBase.h"
11 : #include "DisplacedProblem.h"
12 :
13 : #include "libmesh/quadrature.h"
14 : #include "libmesh/parallel_algebra.h"
15 : #include "libmesh/parallel.h"
16 : #include "libmesh/point.h"
17 : #include "libmesh/dof_map.h"
18 :
19 : #include "libmesh/parallel_ghost_sync.h"
20 : #include "libmesh/mesh_communication.h"
21 :
22 : InputParameters
23 28530 : ActivateElementsUserObjectBase::validParams()
24 : {
25 28530 : InputParameters params = ElementUserObject::validParams();
26 28530 : params.addClassDescription("Determine activated elements.");
27 28530 : params.addRequiredParam<subdomain_id_type>("active_subdomain_id", "The active subdomain ID.");
28 28530 : params.addParam<subdomain_id_type>(
29 : "inactive_subdomain_id",
30 : libMesh::invalid_uint,
31 : "The inactivate subdomain ID, i.e., the subdomain that you want to keep the same.");
32 28530 : params.addRequiredParam<std::vector<BoundaryName>>("expand_boundary_name",
33 : "The expanded boundary name.");
34 28530 : params.registerBase("MeshModifier");
35 28530 : return params;
36 0 : }
37 :
38 0 : ActivateElementsUserObjectBase::ActivateElementsUserObjectBase(const InputParameters & parameters)
39 : : ElementUserObject(parameters),
40 0 : _active_subdomain_id(declareRestartableData<subdomain_id_type>(
41 : "active_subdomain_id", getParam<subdomain_id_type>("active_subdomain_id"))),
42 0 : _inactive_subdomain_id(declareRestartableData<subdomain_id_type>(
43 : "inactive_subdomain_id", getParam<subdomain_id_type>("inactive_subdomain_id"))),
44 0 : _expand_boundary_name(getParam<std::vector<BoundaryName>>("expand_boundary_name"))
45 : {
46 0 : setNewBoundayName();
47 0 : }
48 :
49 : void
50 0 : ActivateElementsUserObjectBase::setNewBoundayName()
51 : {
52 : mooseAssert(!_expand_boundary_name.empty(), "Expanded boundary name is empty");
53 : // add the new boundary and get its boundary id
54 0 : _boundary_ids = _mesh.getBoundaryIDs(_expand_boundary_name, true);
55 : mooseAssert(!_boundary_ids.empty(), "Boundary ID is empty.");
56 0 : _mesh.setBoundaryName(_boundary_ids[0], _expand_boundary_name[0]);
57 0 : _mesh.getMesh().get_boundary_info().sideset_name(_boundary_ids[0]) = _expand_boundary_name[0];
58 0 : _mesh.getMesh().get_boundary_info().nodeset_name(_boundary_ids[0]) = _expand_boundary_name[0];
59 :
60 0 : auto displaced_problem = _fe_problem.getDisplacedProblem();
61 0 : if (displaced_problem)
62 : {
63 0 : _disp_boundary_ids = displaced_problem->mesh().getBoundaryIDs(_expand_boundary_name, true);
64 : mooseAssert(!_disp_boundary_ids.empty(), "Boundary ID in the displaced mesh is empty");
65 :
66 0 : displaced_problem->mesh().setBoundaryName(_disp_boundary_ids[0], _expand_boundary_name[0]);
67 0 : displaced_problem->mesh().getMesh().get_boundary_info().sideset_name(_disp_boundary_ids[0]) =
68 0 : _expand_boundary_name[0];
69 0 : displaced_problem->mesh().getMesh().get_boundary_info().nodeset_name(_disp_boundary_ids[0]) =
70 0 : _expand_boundary_name[0];
71 : }
72 0 : }
73 :
74 : void
75 0 : ActivateElementsUserObjectBase::execute()
76 : {
77 0 : if (isElementActivated() && _current_elem->subdomain_id() != _active_subdomain_id &&
78 0 : _current_elem->subdomain_id() != _inactive_subdomain_id)
79 : {
80 : /*
81 : _current_elem subdomain id is not assignable
82 : create a copy of this element from MooseMesh
83 : */
84 0 : dof_id_type ele_id = _current_elem->id();
85 0 : Elem * ele = _mesh.elemPtr(ele_id);
86 :
87 : // Add element to the activate subdomain
88 0 : ele->subdomain_id() = _active_subdomain_id;
89 :
90 : // Reassign element in the reference mesh while using a displaced mesh
91 0 : auto displaced_problem = _fe_problem.getDisplacedProblem();
92 0 : if (displaced_problem)
93 : {
94 0 : Elem * disp_ele = displaced_problem->mesh().elemPtr(ele_id);
95 0 : disp_ele->subdomain_id() = _active_subdomain_id;
96 : }
97 :
98 : // Save the newly activated element id and node for updating boundary info later
99 0 : _newly_activated_elem.insert(ele_id);
100 0 : for (unsigned int i = 0; i < ele->n_nodes(); ++i)
101 0 : _newly_activated_node.insert(ele->node_id(i));
102 0 : }
103 0 : }
104 :
105 : void
106 0 : ActivateElementsUserObjectBase::finalize()
107 : {
108 : /*
109 : Synchronize ghost element subdomain ID
110 : Note: this needs to be done before updating boundary info because
111 : updating boundary requires the updated element subdomain ids
112 : */
113 0 : libMesh::SyncSubdomainIds sync(_mesh.getMesh());
114 0 : Parallel::sync_dofobject_data_by_id(_mesh.getMesh().comm(),
115 0 : _mesh.getMesh().elements_begin(),
116 0 : _mesh.getMesh().elements_end(),
117 : sync);
118 : // Update boundary info
119 0 : updateBoundaryInfo(_mesh);
120 :
121 : // Similarly for the displaced mesh
122 0 : auto displaced_problem = _fe_problem.getDisplacedProblem();
123 0 : if (displaced_problem)
124 : {
125 0 : libMesh::SyncSubdomainIds sync_mesh(displaced_problem->mesh().getMesh());
126 0 : Parallel::sync_dofobject_data_by_id(displaced_problem->mesh().getMesh().comm(),
127 0 : displaced_problem->mesh().getMesh().elements_begin(),
128 0 : displaced_problem->mesh().getMesh().elements_end(),
129 : sync_mesh);
130 0 : updateBoundaryInfo(displaced_problem->mesh());
131 : }
132 :
133 : // Reinit equation systems
134 0 : _fe_problem.meshChanged(
135 : /*intermediate_change=*/false, /*contract_mesh=*/false, /*clean_refinement_flags=*/false);
136 :
137 : // Get storage ranges for the newly activated elements and boundary nodes
138 0 : ConstElemRange & elem_range = *this->getNewlyActivatedElementRange();
139 0 : ConstBndNodeRange & bnd_node_range = *this->getNewlyActivatedBndNodeRange();
140 :
141 : // Apply initial condition for the newly activated elements
142 0 : initSolutions(elem_range, bnd_node_range);
143 :
144 : // Initialize stateful material properties for the newly activated elements
145 0 : _fe_problem.initElementStatefulProps(elem_range, false);
146 :
147 : // Clear the list
148 0 : _newly_activated_elem.clear();
149 0 : _newly_activated_node.clear();
150 :
151 0 : _node_to_remove_from_bnd.clear();
152 0 : }
153 :
154 : void
155 0 : ActivateElementsUserObjectBase::getNodesToRemoveFromBnd(std::set<dof_id_type> & remove_set,
156 : std::set<dof_id_type> & add_set)
157 : {
158 : // get the difference between the remove_set and the add_set,
159 : // save the difference in _node_to_remove_from_bnd
160 0 : int sz = remove_set.size() + add_set.size();
161 0 : std::vector<dof_id_type> v(sz);
162 0 : std::vector<dof_id_type>::iterator it = std::set_difference(
163 : remove_set.begin(), remove_set.end(), add_set.begin(), add_set.end(), v.begin());
164 0 : v.resize(it - v.begin());
165 0 : _node_to_remove_from_bnd.clear();
166 0 : for (auto id : v)
167 0 : _node_to_remove_from_bnd.insert(id);
168 0 : }
169 :
170 : void
171 0 : ActivateElementsUserObjectBase::insertNodeIdsOnSide(const Elem * ele,
172 : const unsigned short int side,
173 : std::set<dof_id_type> & node_ids)
174 : {
175 0 : for (unsigned int i = 0; i < ele->side_ptr(side)->n_nodes(); ++i)
176 0 : node_ids.insert(ele->side_ptr(side)->node_id(i));
177 0 : }
178 :
179 : void
180 0 : ActivateElementsUserObjectBase::updateBoundaryInfo(MooseMesh & mesh)
181 : {
182 : // save the removed ghost sides and associated nodes to sync across processors
183 : std::unordered_map<processor_id_type, std::vector<std::pair<dof_id_type, unsigned int>>>
184 0 : ghost_sides_to_remove;
185 0 : std::unordered_map<processor_id_type, std::vector<dof_id_type>> ghost_nodes_to_remove;
186 :
187 : // save nodes are added and removed
188 0 : std::set<dof_id_type> add_nodes, remove_nodes;
189 :
190 0 : for (auto ele_id : _newly_activated_elem)
191 : {
192 0 : Elem * ele = mesh.elemPtr(ele_id);
193 0 : for (auto s : ele->side_index_range())
194 : {
195 0 : Elem * neighbor_ele = ele->neighbor_ptr(s);
196 0 : if (neighbor_ele == nullptr)
197 : {
198 : // add this side to boundary
199 0 : mesh.getMesh().get_boundary_info().add_side(ele, s, _boundary_ids[0]);
200 0 : insertNodeIdsOnSide(ele, s, add_nodes);
201 : }
202 : else
203 : {
204 0 : if (neighbor_ele->subdomain_id() != _active_subdomain_id &&
205 0 : neighbor_ele->subdomain_id() != _inactive_subdomain_id)
206 : {
207 : // add this side to boundary
208 0 : mesh.getMesh().get_boundary_info().add_side(ele, s, _boundary_ids[0]);
209 0 : insertNodeIdsOnSide(ele, s, add_nodes);
210 : }
211 : else
212 : {
213 : // remove this side from the boundary
214 0 : mesh.getMesh().get_boundary_info().remove_side(ele, s);
215 0 : insertNodeIdsOnSide(ele, s, remove_nodes);
216 :
217 : // remove the neighbor side from the boundary
218 0 : unsigned int neighbor_s = neighbor_ele->which_neighbor_am_i(ele);
219 0 : mesh.getMesh().get_boundary_info().remove_side(neighbor_ele, neighbor_s);
220 0 : insertNodeIdsOnSide(neighbor_ele, neighbor_s, remove_nodes);
221 :
222 0 : if (neighbor_ele->processor_id() != this->processor_id())
223 0 : ghost_sides_to_remove[neighbor_ele->processor_id()].emplace_back(neighbor_ele->id(),
224 : neighbor_s);
225 : }
226 : }
227 : }
228 : }
229 : // make sure to remove nodes that are not in the add list
230 0 : getNodesToRemoveFromBnd(remove_nodes, add_nodes);
231 0 : for (auto node_id : _node_to_remove_from_bnd)
232 0 : mesh.getMesh().get_boundary_info().remove_node(mesh.nodePtr(node_id), _boundary_ids[0]);
233 :
234 : // synchronize boundary information across processors
235 0 : push_boundary_side_info(mesh, ghost_sides_to_remove);
236 0 : push_boundary_node_info(mesh, ghost_nodes_to_remove);
237 0 : mesh.getMesh().get_boundary_info().parallel_sync_side_ids();
238 0 : mesh.getMesh().get_boundary_info().parallel_sync_node_ids();
239 0 : mesh.update();
240 0 : }
241 :
242 : void
243 0 : ActivateElementsUserObjectBase::push_boundary_side_info(
244 : MooseMesh & mesh,
245 : std::unordered_map<processor_id_type, std::vector<std::pair<dof_id_type, unsigned int>>> &
246 : elems_to_push)
247 : {
248 : auto elem_action_functor =
249 0 : [&mesh, this](processor_id_type,
250 0 : const std::vector<std::pair<dof_id_type, unsigned int>> & received_elem)
251 : {
252 : // remove the side
253 0 : for (const auto & pr : received_elem)
254 0 : mesh.getMesh().get_boundary_info().remove_side(
255 0 : mesh.getMesh().elem_ptr(pr.first), pr.second, this->getExpandedBoundaryID());
256 0 : };
257 :
258 0 : Parallel::push_parallel_vector_data(
259 0 : mesh.getMesh().get_boundary_info().comm(), elems_to_push, elem_action_functor);
260 0 : }
261 :
262 : void
263 0 : ActivateElementsUserObjectBase::push_boundary_node_info(
264 : MooseMesh & mesh,
265 : std::unordered_map<processor_id_type, std::vector<dof_id_type>> & nodes_to_push)
266 : {
267 : auto node_action_functor =
268 0 : [&mesh, this](processor_id_type, const std::vector<dof_id_type> & received_nodes)
269 : {
270 0 : for (const auto & pr : received_nodes)
271 : {
272 : // remove the node
273 0 : mesh.getMesh().get_boundary_info().remove_node(mesh.getMesh().node_ptr(pr),
274 0 : this->getExpandedBoundaryID());
275 : }
276 0 : };
277 :
278 0 : Parallel::push_parallel_vector_data(
279 0 : mesh.getMesh().get_boundary_info().comm(), nodes_to_push, node_action_functor);
280 0 : }
281 :
282 : ConstElemRange *
283 0 : ActivateElementsUserObjectBase::getNewlyActivatedElementRange()
284 : {
285 : // deletes the object first
286 0 : _activated_elem_range.reset();
287 :
288 : // create a vector of the newly activated elements
289 0 : std::vector<Elem *> elems;
290 0 : for (auto elem_id : _newly_activated_elem)
291 0 : elems.push_back(_mesh.elemPtr(elem_id));
292 :
293 : // Make some fake element iterators defining this vector of
294 : // elements
295 0 : Elem * const * elempp = const_cast<Elem * const *>(elems.data());
296 0 : Elem * const * elemend = elempp + elems.size();
297 :
298 : const auto elems_begin =
299 0 : MeshBase::const_element_iterator(elempp, elemend, Predicates::NotNull<Elem * const *>());
300 :
301 : const auto elems_end =
302 0 : MeshBase::const_element_iterator(elemend, elemend, Predicates::NotNull<Elem * const *>());
303 0 : if (!_activated_elem_range)
304 0 : _activated_elem_range = std::make_unique<ConstElemRange>(elems_begin, elems_end);
305 :
306 0 : return _activated_elem_range.get();
307 0 : }
308 :
309 : ConstBndNodeRange *
310 0 : ActivateElementsUserObjectBase::getNewlyActivatedBndNodeRange()
311 : {
312 : // deletes the object first
313 0 : _activated_bnd_node_range.reset();
314 :
315 : // create a vector of the newly activated nodes
316 0 : std::vector<const BndNode *> nodes;
317 0 : std::set<const BndNode *> set_nodes;
318 0 : ConstBndNodeRange & bnd_nodes = *_mesh.getBoundaryNodeRange();
319 0 : for (auto & bnode : bnd_nodes)
320 : {
321 0 : dof_id_type bnode_id = bnode->_node->id();
322 0 : auto it = _newly_activated_node.find(bnode_id);
323 0 : if (it != _newly_activated_node.end())
324 0 : set_nodes.insert(bnode);
325 : }
326 :
327 0 : nodes.assign(set_nodes.begin(), set_nodes.end());
328 :
329 : // Make some fake element iterators defining this vector of
330 : // nodes
331 0 : BndNode * const * nodepp = const_cast<BndNode * const *>(nodes.data());
332 0 : BndNode * const * nodeend = nodepp + nodes.size();
333 :
334 : const auto nodes_begin =
335 0 : MooseMesh::const_bnd_node_iterator(nodepp, nodeend, Predicates::NotNull<BndNode * const *>());
336 :
337 : const auto nodes_end = MooseMesh::const_bnd_node_iterator(
338 0 : nodeend, nodeend, Predicates::NotNull<BndNode * const *>());
339 :
340 0 : if (!_activated_bnd_node_range)
341 0 : _activated_bnd_node_range = std::make_unique<ConstBndNodeRange>(nodes_begin, nodes_end);
342 :
343 0 : return _activated_bnd_node_range.get();
344 0 : }
345 :
346 : ConstNodeRange *
347 0 : ActivateElementsUserObjectBase::getNewlyActivatedNodeRange()
348 : {
349 : // deletes the object first
350 0 : _activated_node_range.reset();
351 :
352 : // create a vector of the newly activated nodes
353 0 : std::vector<const Node *> nodes;
354 0 : for (auto elem_id : _newly_activated_elem)
355 : {
356 0 : const Node * const * elem_nodes = _mesh.elemPtr(elem_id)->get_nodes();
357 0 : unsigned int n_nodes = _mesh.elemPtr(elem_id)->n_nodes();
358 0 : for (unsigned int n = 0; n < n_nodes; ++n)
359 : {
360 : // check if all the elements connected to this node are newly activated
361 0 : const Node * nd = elem_nodes[n];
362 0 : if (isNewlyActivated(nd))
363 0 : nodes.push_back(nd);
364 : }
365 : }
366 :
367 : // Make some fake node iterators defining this vector of
368 : // nodes
369 0 : Node * const * nodepp = const_cast<Node * const *>(nodes.data());
370 0 : Node * const * nodeend = nodepp + nodes.size();
371 :
372 : const auto nodes_begin =
373 0 : MeshBase::const_node_iterator(nodepp, nodeend, Predicates::NotNull<Node * const *>());
374 :
375 : const auto nodes_end =
376 0 : MeshBase::const_node_iterator(nodeend, nodeend, Predicates::NotNull<Node * const *>());
377 :
378 0 : if (!_activated_node_range)
379 0 : _activated_node_range = std::make_unique<ConstNodeRange>(nodes_begin, nodes_end);
380 :
381 0 : return _activated_node_range.get();
382 0 : }
383 :
384 : bool
385 0 : ActivateElementsUserObjectBase::isNewlyActivated(const Node * nd)
386 : {
387 0 : const auto & node_to_elem_map = _mesh.nodeToElemMap();
388 0 : auto node_to_elem_pair = node_to_elem_map.find(nd->id());
389 0 : if (node_to_elem_pair != node_to_elem_map.end())
390 : {
391 0 : const std::vector<dof_id_type> & connected_ele_ids = node_to_elem_pair->second;
392 0 : for (auto connected_ele_id : connected_ele_ids)
393 : {
394 : // check the connected elements
395 0 : if (_mesh.elemPtr(connected_ele_id)->subdomain_id() == _inactive_subdomain_id)
396 0 : return false;
397 0 : if (_mesh.elemPtr(connected_ele_id)->subdomain_id() == _active_subdomain_id &&
398 0 : std::find(_newly_activated_elem.begin(), _newly_activated_elem.end(), connected_ele_id) ==
399 0 : _newly_activated_elem.end())
400 0 : return false;
401 : }
402 : }
403 0 : return true;
404 : }
405 :
406 : void
407 0 : ActivateElementsUserObjectBase::initSolutions(ConstElemRange & elem_range,
408 : ConstBndNodeRange & bnd_node_range)
409 : {
410 : // project initial condition to the current solution
411 0 : _fe_problem.projectInitialConditionOnCustomRange(elem_range, bnd_node_range);
412 :
413 0 : auto & nl = _fe_problem.getNonlinearSystemBase(_sys.number());
414 0 : NumericVector<Number> & current_solution = *nl.system().current_local_solution;
415 0 : NumericVector<Number> & old_solution = nl.solutionOld();
416 0 : NumericVector<Number> & older_solution = nl.solutionOlder();
417 :
418 : NumericVector<Number> & current_aux_solution =
419 0 : *_fe_problem.getAuxiliarySystem().system().current_local_solution;
420 0 : NumericVector<Number> & old_aux_solution = _fe_problem.getAuxiliarySystem().solutionOld();
421 0 : NumericVector<Number> & older_aux_solution = _fe_problem.getAuxiliarySystem().solutionOlder();
422 :
423 0 : DofMap & dof_map = nl.dofMap();
424 0 : DofMap & dof_map_aux = _fe_problem.getAuxiliarySystem().dofMap();
425 :
426 0 : std::set<dof_id_type> dofs, dofs_aux;
427 : // get dofs for the newly added elements
428 0 : for (auto & elem : elem_range)
429 : {
430 0 : std::vector<dof_id_type> di, di_aux;
431 0 : dof_map.dof_indices(elem, di);
432 0 : dof_map_aux.dof_indices(elem, di_aux);
433 0 : for (unsigned int i = 0; i < di.size(); ++i)
434 0 : dofs.insert(di[i]);
435 0 : for (unsigned int i = 0; i < di_aux.size(); ++i)
436 0 : dofs_aux.insert(di_aux[i]);
437 :
438 0 : di.clear();
439 0 : di_aux.clear();
440 0 : }
441 :
442 : // update solutions
443 0 : for (auto dof : dofs)
444 : {
445 0 : old_solution.set(dof, current_solution(dof));
446 0 : older_solution.set(dof, current_solution(dof));
447 : }
448 : // update aux solutions
449 0 : for (auto dof_aux : dofs_aux)
450 : {
451 0 : old_aux_solution.set(dof_aux, current_aux_solution(dof_aux));
452 0 : older_aux_solution.set(dof_aux, current_aux_solution(dof_aux));
453 : }
454 :
455 0 : dofs.clear();
456 0 : dofs_aux.clear();
457 :
458 0 : current_solution.close();
459 0 : old_solution.close();
460 0 : older_solution.close();
461 :
462 0 : current_aux_solution.close();
463 0 : old_aux_solution.close();
464 0 : older_aux_solution.close();
465 :
466 0 : _fe_problem.restoreSolutions();
467 0 : }
|
