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 "MultiAppGeneralFieldNearestLocationTransfer.h"
11 :
12 : // MOOSE includes
13 : #include "FEProblem.h"
14 : #include "MooseMesh.h"
15 : #include "MooseTypes.h"
16 : #include "MooseVariableFE.h"
17 : #include "SystemBase.h"
18 : #include "Positions.h"
19 : #include "MooseAppCoordTransform.h"
20 :
21 : #include "libmesh/system.h"
22 :
23 : using namespace libMesh;
24 :
25 : registerMooseObject("MooseApp", MultiAppGeneralFieldNearestLocationTransfer);
26 : registerMooseObjectRenamed("MooseApp",
27 : MultiAppGeneralFieldNearestNodeTransfer,
28 : "12/31/2024 24:00",
29 : MultiAppGeneralFieldNearestLocationTransfer);
30 :
31 : InputParameters
32 37332 : MultiAppGeneralFieldNearestLocationTransfer::validParams()
33 : {
34 37332 : InputParameters params = MultiAppGeneralFieldTransfer::validParams();
35 74664 : params.addClassDescription(
36 : "Transfers field data at the MultiApp position by finding the value at the nearest "
37 : "neighbor(s) in the origin application.");
38 74664 : params.addParam<unsigned int>("num_nearest_points",
39 74664 : 1,
40 : "Number of nearest source (from) points will be chosen to "
41 : "construct a value for the target point. All points will be "
42 : "selected from the same origin mesh!");
43 :
44 : // Nearest node is historically more an extrapolation transfer
45 111996 : params.set<Real>("extrapolation_constant") = GeneralFieldTransfer::BetterOutOfMeshValue;
46 74664 : params.suppressParameter<Real>("extrapolation_constant");
47 : // We dont keep track of both point distance to app and to nearest node, so we cannot guarantee
48 : // that the nearest app (among the local apps, not globally) will hold the nearest location.
49 : // However, if the user knows this is true, we can use this heuristic to reduce the number of apps
50 : // that are requested to provide a candidate value. If the user is wrong, then the nearest
51 : // location is used, which can be from the non-nearest app.
52 186660 : params.renameParam("use_nearest_app", "assume_nearest_app_holds_nearest_location", "");
53 :
54 : // the default of node/centroid switching based on the variable is causing lots of mistakes and
55 : // bad results
56 : std::vector<MooseEnum> source_types = {
57 111996 : MooseEnum("nodes centroids variable_default", "variable_default")};
58 149328 : params.addParam<std::vector<MooseEnum>>(
59 : "source_type", source_types, "Where to get the source values from for each source variable");
60 :
61 : // choose whether to include data from multiple apps when performing nearest-position/
62 : // mesh-divisions based transfers
63 74664 : params.addParam<bool>("group_subapps",
64 74664 : false,
65 : "Whether to group source locations and values from all subapps "
66 : "when considering a nearest-position algorithm");
67 :
68 74664 : return params;
69 186660 : }
70 :
71 3675 : MultiAppGeneralFieldNearestLocationTransfer::MultiAppGeneralFieldNearestLocationTransfer(
72 3675 : const InputParameters & parameters)
73 : : MultiAppGeneralFieldTransfer(parameters),
74 7350 : _num_nearest_points(getParam<unsigned int>("num_nearest_points")),
75 14700 : _group_subapps(getParam<bool>("group_subapps"))
76 : {
77 3675 : if (_source_app_must_contain_point && _nearest_positions_obj)
78 8 : paramError("use_nearest_position",
79 : "We do not support using both nearest positions matching and checking if target "
80 : "points are within an app domain because the KDTrees for nearest-positions matching "
81 : "are (currently) built with data from multiple applications.");
82 3909 : if (_nearest_positions_obj &&
83 4861 : (isParamValid("from_mesh_divisions") || isParamValid("to_mesh_divisions")))
84 0 : paramError("use_nearest_position", "Cannot use nearest positions with mesh divisions");
85 :
86 : // Parameter checks on grouping subapp values
87 3671 : if (_group_subapps && _from_mesh_divisions.empty() && !_nearest_positions_obj)
88 0 : paramError(
89 : "group_subapps",
90 : "This option is only available for using mesh divisions or nearest positions regions");
91 3753 : else if (_group_subapps &&
92 82 : (_from_mesh_division_behavior == MeshDivisionTransferUse::MATCH_DIVISION_INDEX ||
93 82 : _from_mesh_division_behavior == MeshDivisionTransferUse::MATCH_SUBAPP_INDEX))
94 0 : paramError("group_subapps",
95 : "Cannot group subapps when considering nearest-location data as we would lose "
96 : "track of the division index of the source locations");
97 3671 : else if (_group_subapps && _use_nearest_app)
98 0 : paramError(
99 : "group_subapps",
100 : "When using the 'nearest child application' data, the source data (positions and values) "
101 : "are grouped on a per-application basis, so it cannot be agglomerated over all child "
102 : "applications.\nNote that the option to use nearest applications for source restrictions, "
103 : "but further split each child application's domain by regions closest to each position "
104 : "(here the the child application's centroid), which could be conceived when "
105 : "'group_subapps' = false, is also not available.");
106 3671 : }
107 :
108 : void
109 3523 : MultiAppGeneralFieldNearestLocationTransfer::initialSetup()
110 : {
111 3523 : MultiAppGeneralFieldTransfer::initialSetup();
112 :
113 : // Handle the source types ahead of time
114 6950 : const auto & source_types = getParam<std::vector<MooseEnum>>("source_type");
115 3475 : _source_is_nodes.resize(_from_var_names.size());
116 3475 : _use_zero_dof_for_value.resize(_from_var_names.size());
117 3475 : if (source_types.size() != _from_var_names.size())
118 0 : mooseError("Not enough source types specified for this number of variables. Source types must "
119 : "be specified for transfers with multiple variables");
120 6996 : for (const auto var_index : index_range(_from_var_names))
121 : {
122 : // Local app does not own any of the source problems
123 3527 : if (_from_problems.empty())
124 3 : break;
125 :
126 : // Get some info on the source variable
127 3524 : FEProblemBase & from_problem = *_from_problems[0];
128 : MooseVariableFieldBase & from_var =
129 3524 : from_problem.getVariable(0,
130 3524 : _from_var_names[var_index],
131 : Moose::VarKindType::VAR_ANY,
132 : Moose::VarFieldType::VAR_FIELD_ANY);
133 3524 : System & from_sys = from_var.sys().system();
134 3524 : const auto & fe_type = from_sys.variable_type(from_var.number());
135 :
136 : // Select where to get the origin values
137 3524 : if (source_types[var_index] == "nodes")
138 52 : _source_is_nodes[var_index] = true;
139 3472 : else if (source_types[var_index] == "centroids")
140 52 : _source_is_nodes[var_index] = false;
141 : else
142 : {
143 : // "Nodal" variables are continuous for sure at nodes
144 3420 : if (from_var.isNodal())
145 2572 : _source_is_nodes[var_index] = true;
146 : // for everyone else, we know they are continuous at centroids
147 : else
148 848 : _source_is_nodes[var_index] = false;
149 : }
150 :
151 : // Some variables can be sampled directly at their 0 dofs
152 : // - lagrange at nodes on a first order mesh
153 : // - anything constant and elemental obviously has the 0-dof value at the centroid (or
154 : // vertex-average). However, higher order elemental, even monomial, do not hold the centroid
155 : // value at dof index 0 For example: pyramid has dof 0 at the center of the base, prism has dof
156 : // 0 on an edge etc
157 6148 : if ((_source_is_nodes[var_index] && fe_type.family == LAGRANGE &&
158 7269 : !from_problem.mesh().hasSecondOrderElements()) ||
159 4645 : (!_source_is_nodes[var_index] && fe_type.order == CONSTANT))
160 3251 : _use_zero_dof_for_value[var_index] = true;
161 : else
162 273 : _use_zero_dof_for_value[var_index] = false;
163 :
164 : // Check with the source variable that it is not discontinuous at the source
165 3524 : if (_source_is_nodes[var_index])
166 5248 : if (from_var.getContinuity() == DISCONTINUOUS ||
167 2624 : from_var.getContinuity() == SIDE_DISCONTINUOUS)
168 0 : mooseError("Source variable cannot be sampled at nodes as it is discontinuous");
169 :
170 : // Local app does not own any of the target problems
171 3524 : if (_to_problems.empty())
172 3 : break;
173 :
174 : // Check with the target variable that we are not doing awful projections
175 : MooseVariableFieldBase & to_var =
176 7042 : _to_problems[0]->getVariable(0,
177 3521 : _to_var_names[var_index],
178 : Moose::VarKindType::VAR_ANY,
179 : Moose::VarFieldType::VAR_FIELD_ANY);
180 3521 : System & to_sys = to_var.sys().system();
181 3521 : const auto & to_fe_type = to_sys.variable_type(to_var.number());
182 3521 : if (_source_is_nodes[var_index])
183 : {
184 2621 : if (to_fe_type.order == CONSTANT)
185 117 : mooseWarning(
186 : "Transfer is projecting from nearest-nodes to centroids. This is likely causing "
187 : "floating point indetermination in the results because multiple nodes are 'nearest' to "
188 : "a centroid. Please consider using a ProjectionAux to build an elemental source "
189 : "variable (for example constant monomial) before transferring");
190 : }
191 900 : else if (to_var.isNodal())
192 78 : mooseWarning(
193 : "Transfer is projecting from nearest-centroids to nodes. This is likely causing "
194 : "floating point indetermination in the results because multiple centroids are "
195 : "'nearest' to a node. Please consider using a ProjectionAux to build a nodal source "
196 : "variable (for example linear Lagrange) before transferring");
197 : }
198 :
199 : // We need to improve the indexing if we are to allow this
200 3475 : if (!_from_mesh_divisions.empty())
201 420 : for (const auto mesh_div : _from_mesh_divisions)
202 264 : if (mesh_div->getNumDivisions() != _from_mesh_divisions[0]->getNumDivisions())
203 0 : paramError("from_mesh_division",
204 : "This transfer has only been implemented with a uniform number of source mesh "
205 : "divisions across all source applications");
206 3475 : }
207 :
208 : void
209 52520 : MultiAppGeneralFieldNearestLocationTransfer::prepareEvaluationOfInterpValues(
210 : const unsigned int var_index)
211 : {
212 52520 : _local_kdtrees.clear();
213 52520 : _local_points.clear();
214 52520 : _local_values.clear();
215 52520 : buildKDTrees(var_index);
216 52520 : }
217 :
218 : void
219 52520 : MultiAppGeneralFieldNearestLocationTransfer::buildKDTrees(const unsigned int var_index)
220 : {
221 52520 : computeNumSources();
222 52520 : const auto num_apps_per_tree = getNumAppsPerTree();
223 52520 : _local_kdtrees.resize(_num_sources);
224 52520 : _local_points.resize(_num_sources);
225 52520 : _local_values.resize(_num_sources);
226 52520 : unsigned int max_leaf_size = 0;
227 :
228 : // Construct a local KDTree for each source. A source can be a single app or multiple apps
229 : // combined (option for nearest-position / mesh-divisions)
230 108967 : for (const auto i_source : make_range(_num_sources))
231 : {
232 : // Nest a loop on apps in case multiple apps contribute to the same KD-Tree source
233 112942 : for (const auto app_i : make_range(num_apps_per_tree))
234 : {
235 : // Get the current app index
236 56495 : const auto i_from = getAppIndex(i_source, app_i);
237 : // Current position index, if using nearest positions (not used for use_nearest_app)
238 56495 : const auto i_pos = _group_subapps ? i_source : (i_source % getNumDivisions());
239 :
240 : // Get access to the variable and some variable information
241 56495 : FEProblemBase & from_problem = *_from_problems[i_from];
242 56495 : auto & from_mesh = from_problem.mesh(_displaced_source_mesh);
243 : MooseVariableFieldBase & from_var =
244 56495 : from_problem.getVariable(0,
245 56495 : _from_var_names[var_index],
246 : Moose::VarKindType::VAR_ANY,
247 : Moose::VarFieldType::VAR_FIELD_ANY);
248 56495 : System & from_sys = from_var.sys().system();
249 56495 : const unsigned int from_var_num = from_sys.variable_number(getFromVarName(var_index));
250 :
251 : // Build KDTree from nodes and nodal values
252 56495 : if (_source_is_nodes[var_index])
253 : {
254 4953021 : for (const auto & node : from_mesh.getMesh().local_node_ptr_range())
255 : {
256 4899277 : if (node->n_dofs(from_sys.number(), from_var_num) < 1)
257 310284 : continue;
258 :
259 4764126 : if (!_from_blocks.empty() && !inBlocks(_from_blocks, from_mesh, node))
260 5796 : continue;
261 :
262 4758330 : if (!_from_boundaries.empty() && !onBoundaries(_from_boundaries, from_mesh, node))
263 38960 : continue;
264 :
265 : // Handle the various source mesh divisions behaviors
266 : // NOTE: This could be more efficient, as instead of rejecting points in the
267 : // wrong division, we could just be adding them to the tree for the right division
268 4719370 : if (!_from_mesh_divisions.empty())
269 : {
270 97965 : const auto tree_division_index = i_source % getNumDivisions();
271 97965 : const auto node_div_index = _from_mesh_divisions[i_from]->divisionIndex(*node);
272 :
273 : // Spatial restriction is always active
274 97965 : if (node_div_index == MooseMeshDivision::INVALID_DIVISION_INDEX)
275 63621 : continue;
276 : // We fill one tree per division index for matching subapp index or division index. We
277 : // only accept source data from the division index
278 34344 : else if ((_from_mesh_division_behavior ==
279 3744 : MeshDivisionTransferUse::MATCH_DIVISION_INDEX ||
280 3744 : _to_mesh_division_behavior == MeshDivisionTransferUse::MATCH_DIVISION_INDEX ||
281 3744 : _from_mesh_division_behavior ==
282 38088 : MeshDivisionTransferUse::MATCH_SUBAPP_INDEX) &&
283 : tree_division_index != node_div_index)
284 30240 : continue;
285 : }
286 :
287 : // Transformed node is in the reference space, as is the _nearest_positions_obj
288 4625509 : const auto transformed_node = (*_from_transforms[getGlobalSourceAppIndex(i_from)])(*node);
289 :
290 : // Only add to the KDTree nodes that are closest to the 'position'
291 : // When querying values at a target point, the KDTree associated to the closest
292 : // position to the target point is queried
293 : // We do not need to check the positions when using nearest app as we will assume
294 : // (somewhat incorrectly) that all the points in each subapp are closer to that subapp
295 : // than to any other
296 4677377 : if (!_use_nearest_app && _nearest_positions_obj &&
297 51868 : !closestToPosition(i_pos, transformed_node))
298 36516 : continue;
299 :
300 4588993 : _local_points[i_source].push_back(transformed_node);
301 4588993 : const auto dof = node->dof_number(from_sys.number(), from_var_num, 0);
302 4588993 : _local_values[i_source].push_back((*from_sys.solution)(dof));
303 4588993 : if (!_use_zero_dof_for_value[var_index])
304 61363 : flagInvalidSolution(
305 : "Nearest-location is not implemented for this source variable type on "
306 : "this mesh. Returning value at dof 0");
307 53744 : }
308 : }
309 : // Build KDTree from centroids and value at centroids
310 : else
311 : {
312 5502 : for (auto & elem : as_range(from_mesh.getMesh().local_elements_begin(),
313 185829 : from_mesh.getMesh().local_elements_end()))
314 : {
315 177576 : if (elem->n_dofs(from_sys.number(), from_var_num) < 1)
316 101169 : continue;
317 :
318 177198 : if (!_from_blocks.empty() && !inBlocks(_from_blocks, elem))
319 4005 : continue;
320 :
321 173193 : if (!_from_boundaries.empty() && !onBoundaries(_from_boundaries, from_mesh, elem))
322 2304 : continue;
323 :
324 : // Handle the various source mesh divisions behaviors
325 170889 : if (!_from_mesh_divisions.empty())
326 : {
327 75852 : const auto tree_division_index = i_source % getNumDivisions();
328 75852 : const auto elem_div_index = _from_mesh_divisions[i_from]->divisionIndex(*elem);
329 :
330 : // Spatial restriction is always active
331 75852 : if (elem_div_index == MooseMeshDivision::INVALID_DIVISION_INDEX)
332 44712 : continue;
333 : // We fill one tree per division index for matching subapp index or division index. We
334 : // only accept source data from the division index
335 31140 : else if ((_from_mesh_division_behavior ==
336 3240 : MeshDivisionTransferUse::MATCH_DIVISION_INDEX ||
337 3240 : _from_mesh_division_behavior ==
338 34380 : MeshDivisionTransferUse::MATCH_SUBAPP_INDEX) &&
339 : tree_division_index != elem_div_index)
340 27270 : continue;
341 : }
342 :
343 98907 : const auto vertex_average = elem->vertex_average();
344 : // Transformed element vertex average is in the reference space, as is the
345 : // _nearest_positions_obj
346 : const auto transformed_vertex_average =
347 98907 : (*_from_transforms[getGlobalSourceAppIndex(i_from)])(vertex_average);
348 :
349 : // We do not need to check the positions when using nearest app as we will assume
350 : // (somewhat incorrectly) that all the points in each subapp are closer to that subapp
351 131307 : if (!_use_nearest_app && _nearest_positions_obj &&
352 32400 : !closestToPosition(i_pos, transformed_vertex_average))
353 22500 : continue;
354 :
355 76407 : _local_points[i_source].push_back(transformed_vertex_average);
356 76407 : if (_use_zero_dof_for_value[var_index])
357 : {
358 72537 : auto dof = elem->dof_number(from_sys.number(), from_var_num, 0);
359 72537 : _local_values[i_source].push_back((*from_sys.solution)(dof));
360 : }
361 : else
362 3870 : _local_values[i_source].push_back(
363 3870 : from_sys.point_value(from_var_num, vertex_average, elem));
364 2751 : }
365 : }
366 56495 : max_leaf_size = std::max(max_leaf_size, from_mesh.getMaxLeafSize());
367 : }
368 :
369 : // Make a KDTree from the accumulated points data
370 : std::shared_ptr<KDTree> _kd_tree =
371 56447 : std::make_shared<KDTree>(_local_points[i_source], max_leaf_size);
372 56447 : _local_kdtrees[i_source] = _kd_tree;
373 56447 : }
374 52520 : }
375 :
376 : void
377 78836 : MultiAppGeneralFieldNearestLocationTransfer::evaluateInterpValues(
378 : const std::vector<std::pair<Point, unsigned int>> & incoming_points,
379 : std::vector<std::pair<Real, Real>> & outgoing_vals)
380 : {
381 78836 : evaluateInterpValuesNearestNode(incoming_points, outgoing_vals);
382 78836 : }
383 :
384 : void
385 78836 : MultiAppGeneralFieldNearestLocationTransfer::evaluateInterpValuesNearestNode(
386 : const std::vector<std::pair<Point, unsigned int>> & incoming_points,
387 : std::vector<std::pair<Real, Real>> & outgoing_vals)
388 : {
389 78836 : dof_id_type i_pt = 0;
390 6518166 : for (const auto & [pt, mesh_div] : incoming_points)
391 : {
392 : // Reset distance
393 6439330 : outgoing_vals[i_pt].second = std::numeric_limits<Real>::max();
394 6439330 : bool point_found = false;
395 :
396 : // Loop on all sources
397 13188086 : for (const auto i_from : make_range(_num_sources))
398 : {
399 : // Examine all restrictions for the point. This source (KDTree+values) could be ruled out
400 6748756 : if (!checkRestrictionsForSource(pt, mesh_div, i_from))
401 190872 : continue;
402 :
403 : // TODO: Pre-allocate these two work arrays. They will be regularly resized by the
404 : // searches
405 13115768 : std::vector<std::size_t> return_index(_num_nearest_points);
406 6557884 : std::vector<Real> return_dist_sqr(_num_nearest_points);
407 :
408 : // KD Tree can be empty if no points are within block/boundary/bounding box restrictions
409 6557884 : if (_local_kdtrees[i_from]->numberCandidatePoints())
410 : {
411 6549788 : point_found = true;
412 : // Note that we do not need to use the transformed_pt (in the source app frame)
413 : // because the KDTree has been created in the reference frame
414 6549788 : _local_kdtrees[i_from]->neighborSearch(
415 : pt, _num_nearest_points, return_index, return_dist_sqr);
416 6549788 : Real val_sum = 0, dist_sum = 0;
417 13099576 : for (const auto index : return_index)
418 : {
419 6549788 : val_sum += _local_values[i_from][index];
420 6549788 : dist_sum += (_local_points[i_from][index] - pt).norm();
421 : }
422 :
423 : // If the new value found is closer than for other sources, use it
424 6549788 : const auto new_distance = dist_sum / return_dist_sqr.size();
425 6549788 : if (new_distance < outgoing_vals[i_pt].second)
426 6475022 : outgoing_vals[i_pt] = {val_sum / return_index.size(), new_distance};
427 : }
428 6557884 : }
429 :
430 : // none of the source problem meshes were within the restrictions set
431 6439330 : if (!point_found)
432 4769 : outgoing_vals[i_pt] = {GeneralFieldTransfer::BetterOutOfMeshValue,
433 9538 : GeneralFieldTransfer::BetterOutOfMeshValue};
434 6434561 : else if (_search_value_conflicts)
435 : {
436 : // Local source meshes conflicts can happen two ways:
437 : // - two (or more) problems' nearest nodes are equidistant to the target point
438 : // - within a problem, the num_nearest_points'th closest node is as close to the target
439 : // point as the num_nearest_points + 1'th closest node
440 14605 : unsigned int num_equidistant_problems = 0;
441 :
442 47481 : for (const auto i_from : make_range(_num_sources))
443 : {
444 : // Examine all restrictions for the point. This source (KDTree+values) could be ruled out
445 32876 : if (!checkRestrictionsForSource(pt, mesh_div, i_from))
446 16323 : continue;
447 :
448 : // TODO: Pre-allocate these two work arrays. They will be regularly resized by the searches
449 33106 : std::vector<std::size_t> return_index(_num_nearest_points + 1);
450 16553 : std::vector<Real> return_dist_sqr(_num_nearest_points + 1);
451 :
452 : // NOTE: app_index is not valid if _group_subapps = true
453 16553 : const auto app_index = i_from / getNumDivisions();
454 16553 : const auto num_search = _num_nearest_points + 1;
455 :
456 16553 : if (_local_kdtrees[i_from]->numberCandidatePoints())
457 : {
458 15089 : _local_kdtrees[i_from]->neighborSearch(pt, num_search, return_index, return_dist_sqr);
459 15089 : auto num_found = return_dist_sqr.size();
460 :
461 : // Local coordinates only accessible when not using nearest-position
462 : // as we did not keep the index of the source app, only the position index
463 15089 : const Point local_pt = getPointInLocalSourceFrame(app_index, pt);
464 :
465 20141 : if (!_nearest_positions_obj &&
466 5052 : _from_transforms[getGlobalSourceAppIndex(app_index)]->hasCoordinateSystemTypeChange())
467 0 : if (!_skip_coordinate_collapsing)
468 0 : mooseInfo("Search value conflict cannot find the origin point due to the "
469 : "non-uniqueness of the coordinate collapsing reverse mapping");
470 :
471 : // Look for too many equidistant nodes within a problem. First zip then sort by distance
472 15089 : std::vector<std::pair<Real, std::size_t>> zipped_nearest_points;
473 45267 : for (const auto i : make_range(num_found))
474 30178 : zipped_nearest_points.push_back(std::make_pair(return_dist_sqr[i], return_index[i]));
475 15089 : std::sort(zipped_nearest_points.begin(), zipped_nearest_points.end());
476 :
477 : // If two furthest are equally far from target point, then we have an indetermination in
478 : // what is sent in this communication round from this process. However, it may not
479 : // materialize to an actual conflict, as values sent from another process for the
480 : // desired target point could be closer (nearest). There is no way to know at this point
481 : // in the communication that a closer value exists somewhere else
482 30178 : if (num_found > 1 && num_found == num_search &&
483 15089 : MooseUtils::absoluteFuzzyEqual(zipped_nearest_points[num_found - 1].first,
484 15089 : zipped_nearest_points[num_found - 2].first))
485 : {
486 352 : if (_nearest_positions_obj)
487 16 : registerConflict(app_index, 0, pt, outgoing_vals[i_pt].second, true);
488 : else
489 336 : registerConflict(app_index, 0, local_pt, outgoing_vals[i_pt].second, true);
490 : }
491 :
492 : // Recompute the distance for this problem. If it matches the cached value more than
493 : // once it means multiple problems provide equidistant values for this point
494 15089 : Real dist_sum = 0;
495 30178 : for (const auto i : make_range(num_search - 1))
496 : {
497 15089 : auto index = zipped_nearest_points[i].second;
498 15089 : dist_sum += (_local_points[i_from][index] - pt).norm();
499 : }
500 :
501 : // Compare to the selected value found after looking at all the problems
502 15089 : if (MooseUtils::absoluteFuzzyEqual(dist_sum / return_dist_sqr.size(),
503 15089 : outgoing_vals[i_pt].second))
504 : {
505 88 : num_equidistant_problems++;
506 88 : if (num_equidistant_problems > 1)
507 : {
508 0 : if (_nearest_positions_obj)
509 0 : registerConflict(app_index, 0, pt, outgoing_vals[i_pt].second, true);
510 : else
511 0 : registerConflict(app_index, 0, local_pt, outgoing_vals[i_pt].second, true);
512 : }
513 : }
514 15089 : }
515 16553 : }
516 : }
517 :
518 : // Move to next point
519 6439330 : i_pt++;
520 : }
521 78836 : }
522 :
523 : bool
524 18000 : MultiAppGeneralFieldNearestLocationTransfer::inBlocks(const std::set<SubdomainID> & blocks,
525 : const MooseMesh & mesh,
526 : const Elem * elem) const
527 : {
528 : // We need to override the definition of block restriction for an element
529 : // because we have to consider whether each node of an element is adjacent to a block
530 50292 : for (const auto & i_node : make_range(elem->n_nodes()))
531 : {
532 45756 : const auto & node = elem->node_ptr(i_node);
533 45756 : const auto & node_blocks = mesh.getNodeBlockIds(*node);
534 45756 : std::set<SubdomainID> u;
535 45756 : std::set_intersection(blocks.begin(),
536 : blocks.end(),
537 : node_blocks.begin(),
538 : node_blocks.end(),
539 : std::inserter(u, u.begin()));
540 45756 : if (!u.empty())
541 13464 : return true;
542 45756 : }
543 4536 : return false;
544 : }
545 :
546 : void
547 52520 : MultiAppGeneralFieldNearestLocationTransfer::computeNumSources()
548 : {
549 : // Number of source = number of KDTrees.
550 : // Using mesh divisions or nearest-positions, for every app we use 1 tree per division
551 104896 : if (!_from_mesh_divisions.empty() ||
552 52376 : (!_use_nearest_app && _nearest_positions_obj && !_group_subapps))
553 216 : _num_sources = _from_problems.size() * getNumDivisions();
554 : // If we group apps, then we only use one tree per division (nearest-position region)
555 52304 : else if (_nearest_positions_obj && _group_subapps)
556 76 : _num_sources = _nearest_positions_obj->getNumPositions(_fe_problem.getCurrentExecuteOnFlag() ==
557 : EXEC_INITIAL);
558 : // Regular case: 1 KDTree per app
559 : // Also if use_nearest_app = true, the number of problems is better than the number of positions,
560 : // because some of the positions are positions of child applications that are not local
561 : else
562 52228 : _num_sources = _from_problems.size();
563 52520 : }
564 :
565 : unsigned int
566 6798308 : MultiAppGeneralFieldNearestLocationTransfer::getAppIndex(
567 : unsigned int kdtree_index, unsigned int nested_loop_on_app_index) const
568 : {
569 : // Each app is mapped to a single KD Tree
570 6798308 : if (_use_nearest_app)
571 6489 : return kdtree_index;
572 : // We are looping over all the apps that are grouped together
573 6791819 : else if (_group_subapps)
574 11965 : return nested_loop_on_app_index;
575 : // There are num_divisions trees for each app, inner ordering is divisions, so dividing by the
576 : // number of divisions gets us the index of the application
577 : else
578 6779854 : return kdtree_index / getNumDivisions();
579 : }
580 :
581 : unsigned int
582 52520 : MultiAppGeneralFieldNearestLocationTransfer::getNumAppsPerTree() const
583 : {
584 52520 : if (_use_nearest_app)
585 72 : return 1;
586 52448 : else if (_group_subapps)
587 76 : return _from_meshes.size();
588 : else
589 52372 : return 1;
590 : }
591 :
592 : unsigned int
593 14009368 : MultiAppGeneralFieldNearestLocationTransfer::getNumDivisions() const
594 : {
595 : // This is not used currently, but conceptually it is better to only divide the domain with the
596 : // local of local applications rather than the global number of positions (# global applications
597 : // here)
598 14009368 : if (_use_nearest_app)
599 8685 : return _from_meshes.size();
600 : // Each nearest-position region is a division
601 14000683 : else if (_nearest_positions_obj && !_group_subapps)
602 49356 : return _nearest_positions_obj->getNumPositions(_fe_problem.getCurrentExecuteOnFlag() ==
603 49356 : EXEC_INITIAL);
604 : // Assume all mesh divisions (on each sub-app) has the same number of divisions. This is checked
605 13951327 : else if (!_from_mesh_divisions.empty())
606 834753 : return _from_mesh_divisions[0]->getNumDivisions();
607 : // Grouping subapps or no special mode, we do not subdivide
608 : else
609 13116574 : return 1;
610 : }
611 :
612 : Point
613 15089 : MultiAppGeneralFieldNearestLocationTransfer::getPointInLocalSourceFrame(unsigned int i_from,
614 : const Point & pt) const
615 : {
616 :
617 20141 : if (!_nearest_positions_obj &&
618 5052 : (!_from_transforms[getGlobalSourceAppIndex(i_from)]->hasCoordinateSystemTypeChange() ||
619 0 : _skip_coordinate_collapsing))
620 5052 : return _from_transforms[getGlobalSourceAppIndex(i_from)]->mapBack(pt);
621 10037 : else if (!_nearest_positions_obj || !_group_subapps)
622 4392 : return pt - _from_positions[i_from];
623 : else
624 5645 : return pt;
625 : }
626 :
627 : bool
628 6781632 : MultiAppGeneralFieldNearestLocationTransfer::checkRestrictionsForSource(
629 : const Point & pt, const unsigned int mesh_div, const unsigned int i_from) const
630 : {
631 : // Only use the KDTree from the closest position if in "nearest-position" mode
632 6781632 : if (_nearest_positions_obj)
633 : {
634 : // See computeNumSources for the number of sources. i_from is the index in the source loop
635 : // i_from is local if looping on _from_problems as sources, positions are indexed globally
636 : // i_from is already indexing in positions if using group_subapps
637 69207 : auto position_index = i_from; // if _group_subapps
638 69207 : if (_use_nearest_app)
639 11367 : position_index = getGlobalSourceAppIndex(i_from);
640 57840 : else if (!_group_subapps)
641 22596 : position_index = i_from % getNumDivisions();
642 :
643 : // NOTE: if two positions are equi-distant to the point, this will chose one
644 : // This problem is detected if using search_value_conflicts in this call
645 69207 : if (!closestToPosition(position_index, pt))
646 39819 : return false;
647 : }
648 :
649 : // Application index depends on which source/grouping mode we are using
650 6741813 : const unsigned int app_index = getAppIndex(i_from, i_from / getNumDivisions());
651 :
652 : // Check mesh restriction before anything
653 6741813 : if (_source_app_must_contain_point)
654 : {
655 : // We have to be careful that getPointInLocalSourceFrame returns in the reference frame
656 0 : if (_nearest_positions_obj)
657 0 : mooseError("Nearest-positions + source_app_must_contain_point not implemented");
658 : // Transform the point to place it in the local coordinate system
659 0 : const auto local_pt = getPointInLocalSourceFrame(app_index, pt);
660 0 : if (!inMesh(_from_point_locators[app_index].get(), local_pt))
661 0 : return false;
662 : }
663 :
664 : // Check the mesh division. We have handled the restriction of the source locations when
665 : // building the nearest-neighbor trees. We only need to check that we meet the required
666 : // source division index.
667 6741813 : if (!_from_mesh_divisions.empty())
668 : {
669 : mooseAssert(mesh_div != MooseMeshDivision::INVALID_DIVISION_INDEX,
670 : "We should not be receiving point requests with an invalid "
671 : "source mesh division index");
672 218280 : const unsigned int kd_div_index = i_from % getNumDivisions();
673 :
674 : // If matching source mesh divisions to target apps, we check that the index of the target
675 : // application, which was passed in the point request, is equal to the current mesh division
676 218280 : if (_from_mesh_division_behavior == MeshDivisionTransferUse::MATCH_SUBAPP_INDEX &&
677 : mesh_div != kd_div_index)
678 17568 : return false;
679 : // If matching source mesh divisions to target mesh divisions, we check that the index of the
680 : // target mesh division, which was passed in the point request, is equal to the current mesh
681 : // division
682 200712 : else if ((_from_mesh_division_behavior == MeshDivisionTransferUse::MATCH_DIVISION_INDEX ||
683 200712 : _to_mesh_division_behavior == MeshDivisionTransferUse::MATCH_DIVISION_INDEX) &&
684 : mesh_div != kd_div_index)
685 148608 : return false;
686 : }
687 :
688 : // If matching target apps to source mesh divisions, we check that the global index of the
689 : // application is equal to the target mesh division index, which was passed in the point request
690 6577287 : if (_to_mesh_division_behavior == MeshDivisionTransferUse::MATCH_SUBAPP_INDEX &&
691 1650 : mesh_div != getGlobalSourceAppIndex(app_index))
692 1200 : return false;
693 :
694 6574437 : return true;
695 : }
|
