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