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 "NearestNodeLocator.h"
11 : #include "MooseMesh.h"
12 : #include "SubProblem.h"
13 : #include "SecondaryNeighborhoodThread.h"
14 : #include "NearestNodeThread.h"
15 : #include "Moose.h"
16 : #include "KDTree.h"
17 : #include "Conversion.h"
18 : #include "MooseApp.h"
19 :
20 : // libMesh
21 : #include "libmesh/boundary_info.h"
22 : #include "libmesh/elem.h"
23 : #include "libmesh/plane.h"
24 : #include "libmesh/mesh_tools.h"
25 :
26 1812 : NearestNodeLocator::NearestNodeLocator(SubProblem & subproblem,
27 : MooseMesh & mesh,
28 : BoundaryID boundary1,
29 1812 : BoundaryID boundary2)
30 : : Restartable(subproblem.getMooseApp(),
31 3624 : Moose::stringify(boundary1) + Moose::stringify(boundary2),
32 : "NearestNodeLocator",
33 : 0),
34 1812 : PerfGraphInterface(subproblem.getMooseApp().perfGraph(),
35 3624 : "NearestNodeLocator_" + Moose::stringify(boundary1) + "_" +
36 1812 : Moose::stringify(boundary2)),
37 1812 : _subproblem(subproblem),
38 1812 : _mesh(mesh),
39 1812 : _boundary1(boundary1),
40 1812 : _boundary2(boundary2),
41 1812 : _first(true),
42 1812 : _reinit_iteration(true),
43 16308 : _patch_update_strategy(_mesh.getPatchUpdateStrategy())
44 : {
45 : /*
46 : //sanity check on boundary ids
47 : const std::set<BoundaryID>& bids=_mesh.getBoundaryIDs();
48 : std::set<BoundaryID>::const_iterator sit;
49 : sit=bids.find(_boundary1);
50 : if (sit == bids.end())
51 : mooseError("NearestNodeLocator being created for boundaries ", _boundary1, " and ", _boundary2,
52 : ", but boundary ", _boundary1, " does not exist");
53 : sit=bids.find(_boundary2);
54 : if (sit == bids.end())
55 : mooseError("NearestNodeLocator being created for boundaries ", _boundary1, " and ", _boundary2,
56 : ", but boundary ", _boundary2, " does not exist");
57 : */
58 :
59 : // Request the nodeToElem map upfront
60 1812 : _mesh.nodeToElemMap();
61 1812 : }
62 :
63 3612 : NearestNodeLocator::~NearestNodeLocator() = default;
64 :
65 : void
66 172969 : NearestNodeLocator::findNodes()
67 : {
68 864845 : TIME_SECTION("findNodes", 3, "Finding Nearest Nodes");
69 :
70 : /**
71 : * If this is the first time through we're going to build up a "neighborhood" of nodes
72 : * surrounding each of the secondary nodes. This will speed searching later.
73 : */
74 172969 : const std::map<dof_id_type, std::vector<dof_id_type>> & node_to_elem_map = _mesh.nodeToElemMap();
75 :
76 172969 : if (_first || (_reinit_iteration && _patch_update_strategy == Moose::Iteration))
77 : {
78 11448 : _first = false;
79 :
80 : // Trial secondary nodes are all the nodes on the secondary side
81 : // We only keep the ones that are either on this processor or are likely
82 : // to interact with elements on this processor (ie nodes owned by this processor
83 : // are in the "neighborhood" of the secondary node
84 11448 : std::vector<dof_id_type> trial_secondary_nodes;
85 11448 : std::vector<dof_id_type> trial_primary_nodes;
86 :
87 : // Build a bounding box. No reason to consider nodes outside of our inflated BB
88 11448 : std::unique_ptr<BoundingBox> my_inflated_box = nullptr;
89 :
90 11448 : const std::vector<Real> & inflation = _mesh.getGhostedBoundaryInflation();
91 :
92 : // This means there was a user specified inflation... so we can build a BB
93 11448 : if (inflation.size() > 0)
94 : {
95 0 : BoundingBox my_box = MeshTools::create_local_bounding_box(_mesh);
96 :
97 0 : Point distance;
98 0 : for (unsigned int i = 0; i < inflation.size(); ++i)
99 0 : distance(i) = inflation[i];
100 :
101 : my_inflated_box =
102 0 : std::make_unique<BoundingBox>(my_box.first - distance, my_box.second + distance);
103 : }
104 :
105 : // Data structures to hold the boundary nodes
106 11448 : ConstBndNodeRange & bnd_nodes = *_mesh.getBoundaryNodeRange();
107 2312825 : for (const auto & bnode : bnd_nodes)
108 : {
109 2301377 : BoundaryID boundary_id = bnode->_bnd_id;
110 2301377 : dof_id_type node_id = bnode->_node->id();
111 :
112 : // If we have a BB only consider saving this node if it's in our inflated BB
113 2301377 : if (!my_inflated_box || (my_inflated_box->contains_point(*bnode->_node)))
114 : {
115 2301377 : if (boundary_id == _boundary1)
116 1011695 : trial_primary_nodes.push_back(node_id);
117 1289682 : else if (boundary_id == _boundary2)
118 66227 : trial_secondary_nodes.push_back(node_id);
119 : }
120 : }
121 :
122 : // Convert trial primary nodes to a vector of Points. This will be used to
123 : // construct the Kdtree.
124 11448 : std::vector<Point> primary_points(trial_primary_nodes.size());
125 1023143 : for (unsigned int i = 0; i < trial_primary_nodes.size(); ++i)
126 : {
127 1011695 : const Node & node = _mesh.nodeRef(trial_primary_nodes[i]);
128 1011695 : primary_points[i] = node;
129 : }
130 :
131 : // Create object kd_tree of class KDTree using the coordinates of trial
132 : // primary nodes.
133 11448 : KDTree kd_tree(primary_points, _mesh.getMaxLeafSize());
134 :
135 : NodeIdRange trial_secondary_node_range(
136 11448 : trial_secondary_nodes.begin(), trial_secondary_nodes.end(), 1);
137 :
138 : SecondaryNeighborhoodThread snt(
139 11448 : _mesh, trial_primary_nodes, node_to_elem_map, _mesh.getPatchSize(), kd_tree);
140 :
141 11448 : Threads::parallel_reduce(trial_secondary_node_range, snt);
142 :
143 11448 : _secondary_nodes = snt._secondary_nodes;
144 11448 : _neighbor_nodes = snt._neighbor_nodes;
145 :
146 : // If 'iteration' patch update strategy is used, a second neighborhood
147 : // search using the ghosting_patch_size, which is larger than the regular
148 : // patch_size used for contact search, is conducted. The ghosted element set
149 : // given by this search is used for ghosting the elements connected to the
150 : // secondary and neighboring primary nodes.
151 11448 : if (_patch_update_strategy == Moose::Iteration)
152 : {
153 : SecondaryNeighborhoodThread snt_ghosting(
154 8658 : _mesh, trial_primary_nodes, node_to_elem_map, _mesh.getGhostingPatchSize(), kd_tree);
155 :
156 8658 : Threads::parallel_reduce(trial_secondary_node_range, snt_ghosting);
157 :
158 771822 : for (const auto & dof : snt_ghosting._ghosted_elems)
159 763164 : _subproblem.addGhostedElem(dof);
160 8658 : }
161 : else
162 : {
163 107025 : for (const auto & dof : snt._ghosted_elems)
164 104235 : _subproblem.addGhostedElem(dof);
165 : }
166 :
167 : // Cache the secondary_node_range so we don't have to build it each time
168 : _secondary_node_range =
169 11448 : std::make_unique<NodeIdRange>(_secondary_nodes.begin(), _secondary_nodes.end(), 1);
170 11448 : }
171 :
172 172969 : _nearest_node_info.clear();
173 :
174 172969 : NearestNodeThread nnt(_mesh, _neighbor_nodes);
175 :
176 172969 : Threads::parallel_reduce(*_secondary_node_range, nnt);
177 :
178 172969 : _max_patch_percentage = nnt._max_patch_percentage;
179 :
180 172969 : _nearest_node_info = nnt._nearest_node_info;
181 :
182 172969 : if (_patch_update_strategy == Moose::Iteration)
183 : {
184 : // Get the set of elements that are currently being ghosted
185 8658 : std::set<dof_id_type> ghost = _subproblem.ghostedElems();
186 :
187 34257 : for (const auto & node_id : *_secondary_node_range)
188 : {
189 25599 : const Node * nearest_node = _nearest_node_info[node_id]._nearest_node;
190 :
191 : // Check if the elements attached to the nearest node are within the ghosted
192 : // set of elements. If not produce an error.
193 25599 : auto node_to_elem_pair = node_to_elem_map.find(nearest_node->id());
194 :
195 25599 : if (node_to_elem_pair != node_to_elem_map.end())
196 : {
197 25599 : const std::vector<dof_id_type> & elems_connected_to_node = node_to_elem_pair->second;
198 76797 : for (const auto & dof : elems_connected_to_node)
199 85932 : if (std::find(ghost.begin(), ghost.end(), dof) == ghost.end() &&
200 34734 : _mesh.elemPtr(dof)->processor_id() != _mesh.processor_id())
201 0 : mooseError("Error in NearestNodeLocator: The nearest neighbor lies outside the "
202 : "ghosted set of elements. Increase the ghosting_patch_size parameter in the "
203 : "mesh block and try again.");
204 : }
205 : }
206 8658 : }
207 172969 : }
208 :
209 : void
210 1002 : NearestNodeLocator::reinit()
211 : {
212 5010 : TIME_SECTION("reinit", 3, "Reinitializing Nearest Node Search");
213 :
214 : // Reset all data
215 1002 : _secondary_node_range.reset();
216 1002 : _nearest_node_info.clear();
217 :
218 1002 : _first = true;
219 :
220 1002 : _secondary_nodes.clear();
221 1002 : _neighbor_nodes.clear();
222 :
223 1002 : _new_ghosted_elems.clear();
224 :
225 : // After a call from system reinit, mesh has been updated with initial adaptivity.
226 : // Moose::Iteration relies on data generated for ghosting (i.e. trial_primary_nodes)
227 1002 : _reinit_iteration = true;
228 :
229 : // Redo the search
230 1002 : findNodes();
231 :
232 1002 : _reinit_iteration = false;
233 1002 : }
234 :
235 : Real
236 523282 : NearestNodeLocator::distance(dof_id_type node_id)
237 : {
238 523282 : return _nearest_node_info[node_id]._distance;
239 : }
240 :
241 : const Node *
242 566890 : NearestNodeLocator::nearestNode(dof_id_type node_id)
243 : {
244 566890 : const Node * returnval = _nearest_node_info[node_id]._nearest_node;
245 : libmesh_assert(_mesh.getMesh().get_boundary_info().has_boundary_id(returnval, _boundary1));
246 566890 : return returnval;
247 : }
248 :
249 : void
250 0 : NearestNodeLocator::updatePatch(std::vector<dof_id_type> & secondary_nodes)
251 : {
252 0 : TIME_SECTION("updatePatch", 3, "Updating Nearest Node Search Patch");
253 :
254 0 : std::vector<dof_id_type> trial_primary_nodes;
255 :
256 : // Build a bounding box. No reason to consider nodes outside of our inflated BB
257 0 : std::unique_ptr<BoundingBox> my_inflated_box = nullptr;
258 :
259 0 : const std::vector<Real> & inflation = _mesh.getGhostedBoundaryInflation();
260 :
261 : // This means there was a user specified inflation... so we can build a BB
262 0 : if (inflation.size() > 0)
263 : {
264 0 : BoundingBox my_box = MeshTools::create_local_bounding_box(_mesh);
265 :
266 0 : Point distance;
267 0 : for (unsigned int i = 0; i < inflation.size(); ++i)
268 0 : distance(i) = inflation[i];
269 :
270 : my_inflated_box =
271 0 : std::make_unique<BoundingBox>(my_box.first - distance, my_box.second + distance);
272 : }
273 :
274 : // Data structures to hold the boundary nodes
275 0 : ConstBndNodeRange & bnd_nodes = *_mesh.getBoundaryNodeRange();
276 0 : for (const auto & bnode : bnd_nodes)
277 : {
278 0 : BoundaryID boundary_id = bnode->_bnd_id;
279 0 : dof_id_type node_id = bnode->_node->id();
280 :
281 : // If we have a BB only consider saving this node if it's in our inflated BB
282 0 : if (!my_inflated_box || (my_inflated_box->contains_point(*bnode->_node)))
283 : {
284 0 : if (boundary_id == _boundary1)
285 0 : trial_primary_nodes.push_back(node_id);
286 : }
287 : }
288 :
289 : // Convert trial primary nodes to a vector of Points. This will be used to construct the KDTree.
290 0 : std::vector<Point> primary_points(trial_primary_nodes.size());
291 0 : for (unsigned int i = 0; i < trial_primary_nodes.size(); ++i)
292 : {
293 0 : const Node & node = _mesh.nodeRef(trial_primary_nodes[i]);
294 0 : primary_points[i] = node;
295 : }
296 :
297 0 : const std::map<dof_id_type, std::vector<dof_id_type>> & node_to_elem_map = _mesh.nodeToElemMap();
298 :
299 : // Create object kd_tree of class KDTree using the coordinates of trial
300 : // primary nodes.
301 0 : KDTree kd_tree(primary_points, _mesh.getMaxLeafSize());
302 :
303 0 : NodeIdRange secondary_node_range(secondary_nodes.begin(), secondary_nodes.end(), 1);
304 :
305 : SecondaryNeighborhoodThread snt(
306 0 : _mesh, trial_primary_nodes, node_to_elem_map, _mesh.getPatchSize(), kd_tree);
307 :
308 0 : Threads::parallel_reduce(secondary_node_range, snt);
309 :
310 : // Calculate new ghosting patch for the secondary_node_range
311 : SecondaryNeighborhoodThread snt_ghosting(
312 0 : _mesh, trial_primary_nodes, node_to_elem_map, _mesh.getGhostingPatchSize(), kd_tree);
313 :
314 0 : Threads::parallel_reduce(secondary_node_range, snt_ghosting);
315 :
316 : // Add the new set of elements that need to be ghosted into _new_ghosted_elems
317 0 : for (const auto & dof : snt_ghosting._ghosted_elems)
318 0 : _new_ghosted_elems.push_back(dof);
319 :
320 0 : std::vector<dof_id_type> tracked_secondary_nodes = snt._secondary_nodes;
321 :
322 : // Update the neighbor nodes (patch) for these tracked secondary nodes
323 0 : for (const auto & node_id : tracked_secondary_nodes)
324 0 : _neighbor_nodes[node_id] = snt._neighbor_nodes[node_id];
325 :
326 : NodeIdRange tracked_secondary_node_range(
327 0 : tracked_secondary_nodes.begin(), tracked_secondary_nodes.end(), 1);
328 :
329 0 : NearestNodeThread nnt(_mesh, snt._neighbor_nodes);
330 :
331 0 : Threads::parallel_reduce(tracked_secondary_node_range, nnt);
332 :
333 0 : _max_patch_percentage = nnt._max_patch_percentage;
334 :
335 : // Get the set of elements that are currently being ghosted
336 0 : std::set<dof_id_type> ghost = _subproblem.ghostedElems();
337 :
338 : // Update the nearest node information corresponding to these tracked secondary nodes
339 0 : for (const auto & node_id : tracked_secondary_node_range)
340 : {
341 0 : _nearest_node_info[node_id] = nnt._nearest_node_info[node_id];
342 :
343 : // Check if the elements attached to the nearest node are within the ghosted
344 : // set of elements. If not produce an error.
345 0 : const Node * nearest_node = nnt._nearest_node_info[node_id]._nearest_node;
346 :
347 0 : auto node_to_elem_pair = node_to_elem_map.find(nearest_node->id());
348 :
349 0 : if (node_to_elem_pair != node_to_elem_map.end())
350 : {
351 0 : const std::vector<dof_id_type> & elems_connected_to_node = node_to_elem_pair->second;
352 0 : for (const auto & dof : elems_connected_to_node)
353 0 : if (std::find(ghost.begin(), ghost.end(), dof) == ghost.end() &&
354 0 : _mesh.elemPtr(dof)->processor_id() != _mesh.processor_id())
355 0 : mooseError("Error in NearestNodeLocator: The nearest neighbor lies outside the ghosted "
356 : "set of elements. Increase the ghosting_patch_size parameter in the mesh "
357 : "block and try again.");
358 : }
359 : }
360 0 : }
361 :
362 : void
363 662 : NearestNodeLocator::updateGhostedElems()
364 : {
365 3310 : TIME_SECTION("updateGhostedElems", 5, "Updating Nearest Node Search Because of Ghosting");
366 :
367 : // When 'iteration' patch update strategy is used, add the elements in
368 : // _new_ghosted_elems, which were accumulated in the nonlinear iterations
369 : // during the previous time step, to the list of ghosted elements. Also clear
370 : // the _new_ghosted_elems array for storing the ghosted elements from the
371 : // nonlinear iterations in the current time step.
372 :
373 662 : for (const auto & dof : _new_ghosted_elems)
374 0 : _subproblem.addGhostedElem(dof);
375 :
376 662 : _new_ghosted_elems.clear();
377 662 : }
378 : //===================================================================
379 1202737 : NearestNodeLocator::NearestNodeInfo::NearestNodeInfo()
380 1202737 : : _nearest_node(nullptr), _distance(std::numeric_limits<Real>::max())
381 : {
382 1202737 : }
|
