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 35254 : MultiAppGeneralFieldNearestLocationTransfer::validParams()
33 : {
34 35254 : InputParameters params = MultiAppGeneralFieldTransfer::validParams();
35 35254 : 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 105762 : params.addParam<unsigned int>("num_nearest_points",
39 70508 : 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 35254 : params.set<Real>("extrapolation_constant") = GeneralFieldTransfer::BetterOutOfMeshValue;
46 35254 : 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 35254 : 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 70508 : MooseEnum("nodes centroids variable_default", "variable_default")};
58 35254 : 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 105762 : params.addParam<bool>("group_subapps",
64 70508 : false,
65 : "Whether to group source locations and values from all subapps "
66 : "when considering a nearest-position algorithm");
67 :
68 70508 : return params;
69 105762 : }
70 :
71 3362 : MultiAppGeneralFieldNearestLocationTransfer::MultiAppGeneralFieldNearestLocationTransfer(
72 3362 : const InputParameters & parameters)
73 : : MultiAppGeneralFieldTransfer(parameters),
74 : SolutionInvalidInterface(this),
75 3362 : _num_nearest_points(getParam<unsigned int>("num_nearest_points")),
76 6724 : _group_subapps(getParam<bool>("group_subapps"))
77 : {
78 3362 : 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 3578 : if (_nearest_positions_obj &&
84 3578 : (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 3358 : 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 3434 : else if (_group_subapps &&
93 76 : (_from_mesh_division_behavior == MeshDivisionTransferUse::MATCH_DIVISION_INDEX ||
94 76 : _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 3358 : 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 3358 : }
108 :
109 : void
110 3258 : MultiAppGeneralFieldNearestLocationTransfer::initialSetup()
111 : {
112 3258 : MultiAppGeneralFieldTransfer::initialSetup();
113 :
114 : // Handle the source types ahead of time
115 3226 : const auto & source_types = getParam<std::vector<MooseEnum>>("source_type");
116 3226 : _source_is_nodes.resize(_from_var_names.size());
117 3226 : _use_zero_dof_for_value.resize(_from_var_names.size());
118 3226 : 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 6494 : for (const auto var_index : index_range(_from_var_names))
122 : {
123 : // Local app does not own any of the source problems
124 3274 : if (_from_problems.empty())
125 3 : break;
126 :
127 : // Get some info on the source variable
128 3271 : FEProblemBase & from_problem = *_from_problems[0];
129 : MooseVariableFieldBase & from_var =
130 3271 : from_problem.getVariable(0,
131 3271 : _from_var_names[var_index],
132 : Moose::VarKindType::VAR_ANY,
133 : Moose::VarFieldType::VAR_FIELD_ANY);
134 3271 : System & from_sys = from_var.sys().system();
135 3271 : const auto & fe_type = from_sys.variable_type(from_var.number());
136 :
137 : // Select where to get the origin values
138 3271 : if (source_types[var_index] == "nodes")
139 48 : _source_is_nodes[var_index] = true;
140 3223 : else if (source_types[var_index] == "centroids")
141 48 : _source_is_nodes[var_index] = false;
142 : else
143 : {
144 : // "Nodal" variables are continuous for sure at nodes
145 3175 : if (from_var.isNodal())
146 2391 : _source_is_nodes[var_index] = true;
147 : // for everyone else, we know they are continuous at centroids
148 : else
149 784 : _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 5710 : if ((_source_is_nodes[var_index] && fe_type.family == LAGRANGE &&
159 6746 : !from_problem.mesh().hasSecondOrderElements()) ||
160 4307 : (!_source_is_nodes[var_index] && fe_type.order == CONSTANT))
161 3019 : _use_zero_dof_for_value[var_index] = true;
162 : else
163 252 : _use_zero_dof_for_value[var_index] = false;
164 :
165 : // Check with the source variable that it is not discontinuous at the source
166 3271 : if (_source_is_nodes[var_index])
167 4878 : if (from_var.getContinuity() == DISCONTINUOUS ||
168 2439 : 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 3271 : 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 6536 : _to_problems[0]->getVariable(0,
178 3268 : _to_var_names[var_index],
179 : Moose::VarKindType::VAR_ANY,
180 : Moose::VarFieldType::VAR_FIELD_ANY);
181 3268 : System & to_sys = to_var.sys().system();
182 3268 : const auto & to_fe_type = to_sys.variable_type(to_var.number());
183 3268 : if (_source_is_nodes[var_index])
184 : {
185 2436 : if (to_fe_type.order == CONSTANT)
186 108 : 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 832 : else if (to_var.isNodal())
193 72 : 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 3226 : if (!_from_mesh_divisions.empty())
202 384 : for (const auto mesh_div : _from_mesh_divisions)
203 240 : 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 3226 : }
208 :
209 : void
210 51420 : MultiAppGeneralFieldNearestLocationTransfer::prepareEvaluationOfInterpValues(
211 : const unsigned int var_index)
212 : {
213 51420 : _local_kdtrees.clear();
214 51420 : _local_points.clear();
215 51420 : _local_values.clear();
216 51420 : buildKDTrees(var_index);
217 51420 : }
218 :
219 : void
220 51420 : MultiAppGeneralFieldNearestLocationTransfer::buildKDTrees(const unsigned int var_index)
221 : {
222 51420 : computeNumSources();
223 51420 : const auto num_apps_per_tree = getNumAppsPerTree();
224 51420 : _local_kdtrees.resize(_num_sources);
225 51420 : _local_points.resize(_num_sources);
226 51420 : _local_values.resize(_num_sources);
227 51420 : 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 107048 : 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 111496 : for (const auto app_i : make_range(num_apps_per_tree))
235 : {
236 : // Get the current app index
237 55868 : const auto i_from = getAppIndex(i_source, app_i);
238 : // Current position index, if using nearest positions (not used for use_nearest_app)
239 55868 : const auto i_pos = _group_subapps ? i_source : (i_source % getNumDivisions());
240 :
241 : // Get access to the variable and some variable information
242 55868 : FEProblemBase & from_problem = *_from_problems[i_from];
243 55868 : auto & from_mesh = from_problem.mesh(_displaced_source_mesh);
244 : MooseVariableFieldBase & from_var =
245 55868 : from_problem.getVariable(0,
246 55868 : _from_var_names[var_index],
247 : Moose::VarKindType::VAR_ANY,
248 : Moose::VarFieldType::VAR_FIELD_ANY);
249 55868 : System & from_sys = from_var.sys().system();
250 55868 : const unsigned int from_var_num = from_sys.variable_number(getFromVarName(var_index));
251 :
252 : // Build KDTree from nodes and nodal values
253 55868 : if (_source_is_nodes[var_index])
254 : {
255 9119456 : for (const auto & node : from_mesh.getMesh().local_node_ptr_range())
256 : {
257 4533646 : if (node->n_dofs(from_sys.number(), from_var_num) < 1)
258 289704 : continue;
259 :
260 4412134 : if (!_from_blocks.empty() && !inBlocks(_from_blocks, from_mesh, node))
261 5152 : continue;
262 :
263 4406982 : if (!_from_boundaries.empty() && !onBoundaries(_from_boundaries, from_mesh, node))
264 34544 : 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 4372438 : if (!_from_mesh_divisions.empty())
270 : {
271 87080 : const auto tree_division_index = i_source % getNumDivisions();
272 87080 : const auto node_div_index = _from_mesh_divisions[i_from]->divisionIndex(*node);
273 :
274 : // Spatial restriction is always active
275 87080 : if (node_div_index == MooseMeshDivision::INVALID_DIVISION_INDEX)
276 56552 : 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 30528 : else if ((_from_mesh_division_behavior ==
280 3328 : MeshDivisionTransferUse::MATCH_DIVISION_INDEX ||
281 3328 : _to_mesh_division_behavior == MeshDivisionTransferUse::MATCH_DIVISION_INDEX ||
282 3328 : _from_mesh_division_behavior ==
283 33856 : MeshDivisionTransferUse::MATCH_SUBAPP_INDEX) &&
284 : tree_division_index != node_div_index)
285 26880 : continue;
286 : }
287 :
288 : // Transformed node is in the reference space, as is the _nearest_positions_obj
289 4289006 : 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 4359598 : if (!_use_nearest_app && _nearest_positions_obj &&
298 70592 : !closestToPosition(i_pos, transformed_node))
299 45064 : continue;
300 :
301 4243942 : _local_points[i_source].push_back(transformed_node);
302 4243942 : const auto dof = node->dof_number(from_sys.number(), from_var_num, 0);
303 4243942 : _local_values[i_source].push_back((*from_sys.solution)(dof));
304 4243942 : if (!_use_zero_dof_for_value[var_index])
305 53920 : flagInvalidSolution(
306 : "Nearest-location is not implemented for this source variable type on "
307 : "this mesh. Returning value at dof 0");
308 52164 : }
309 : }
310 : // Build KDTree from centroids and value at centroids
311 : else
312 : {
313 7408 : for (auto & elem : as_range(from_mesh.getMesh().local_elements_begin(),
314 376616 : from_mesh.getMesh().local_elements_end()))
315 : {
316 182752 : if (elem->n_dofs(from_sys.number(), from_var_num) < 1)
317 99608 : continue;
318 :
319 180736 : if (!_from_blocks.empty() && !inBlocks(_from_blocks, elem))
320 3560 : continue;
321 :
322 177176 : if (!_from_boundaries.empty() && !onBoundaries(_from_boundaries, from_mesh, elem))
323 2048 : continue;
324 :
325 : // Handle the various source mesh divisions behaviors
326 175128 : if (!_from_mesh_divisions.empty())
327 : {
328 67424 : const auto tree_division_index = i_source % getNumDivisions();
329 67424 : const auto elem_div_index = _from_mesh_divisions[i_from]->divisionIndex(*elem);
330 :
331 : // Spatial restriction is always active
332 67424 : if (elem_div_index == MooseMeshDivision::INVALID_DIVISION_INDEX)
333 39744 : 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 27680 : else if ((_from_mesh_division_behavior ==
337 2880 : MeshDivisionTransferUse::MATCH_DIVISION_INDEX ||
338 2880 : _from_mesh_division_behavior ==
339 30560 : MeshDivisionTransferUse::MATCH_SUBAPP_INDEX) &&
340 : tree_division_index != elem_div_index)
341 24240 : continue;
342 : }
343 :
344 111144 : 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 111144 : (*_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 155944 : if (!_use_nearest_app && _nearest_positions_obj &&
353 44800 : !closestToPosition(i_pos, transformed_vertex_average))
354 28000 : continue;
355 :
356 83144 : _local_points[i_source].push_back(transformed_vertex_average);
357 83144 : if (_use_zero_dof_for_value[var_index])
358 : {
359 79704 : auto dof = elem->dof_number(from_sys.number(), from_var_num, 0);
360 79704 : _local_values[i_source].push_back((*from_sys.solution)(dof));
361 : }
362 : else
363 3440 : _local_values[i_source].push_back(
364 3440 : from_sys.point_value(from_var_num, vertex_average, elem));
365 3704 : }
366 : }
367 55868 : 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 55628 : std::make_shared<KDTree>(_local_points[i_source], max_leaf_size);
373 55628 : _local_kdtrees[i_source] = _kd_tree;
374 55628 : }
375 51420 : }
376 :
377 : void
378 79120 : MultiAppGeneralFieldNearestLocationTransfer::evaluateInterpValues(
379 : const std::vector<std::pair<Point, unsigned int>> & incoming_points,
380 : std::vector<std::pair<Real, Real>> & outgoing_vals)
381 : {
382 79120 : evaluateInterpValuesNearestNode(incoming_points, outgoing_vals);
383 79120 : }
384 :
385 : void
386 79120 : MultiAppGeneralFieldNearestLocationTransfer::evaluateInterpValuesNearestNode(
387 : const std::vector<std::pair<Point, unsigned int>> & incoming_points,
388 : std::vector<std::pair<Real, Real>> & outgoing_vals)
389 : {
390 79120 : dof_id_type i_pt = 0;
391 6198407 : for (const auto & [pt, mesh_div] : incoming_points)
392 : {
393 : // Reset distance
394 6119287 : outgoing_vals[i_pt].second = std::numeric_limits<Real>::max();
395 6119287 : bool point_found = false;
396 :
397 : // Loop on all sources
398 12553132 : 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 6433845 : if (!checkRestrictionsForSource(pt, mesh_div, i_from))
402 206880 : continue;
403 :
404 : // TODO: Pre-allocate these two work arrays. They will be regularly resized by the
405 : // searches
406 6226965 : std::vector<std::size_t> return_index(_num_nearest_points);
407 6226965 : 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 6226965 : if (_local_kdtrees[i_from]->numberCandidatePoints())
411 : {
412 6213028 : 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 6213028 : _local_kdtrees[i_from]->neighborSearch(
416 : pt, _num_nearest_points, return_index, return_dist_sqr);
417 6213028 : Real val_sum = 0, dist_sum = 0;
418 12426056 : for (const auto index : return_index)
419 : {
420 6213028 : val_sum += _local_values[i_from][index];
421 6213028 : 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 6213028 : const auto new_distance = dist_sum / return_dist_sqr.size();
426 6213028 : if (new_distance < outgoing_vals[i_pt].second)
427 6149277 : outgoing_vals[i_pt] = {val_sum / return_index.size(), new_distance};
428 : }
429 6226965 : }
430 :
431 : // none of the source problem meshes were within the restrictions set
432 6119287 : if (!point_found)
433 5143 : outgoing_vals[i_pt] = {GeneralFieldTransfer::BetterOutOfMeshValue,
434 10286 : GeneralFieldTransfer::BetterOutOfMeshValue};
435 6114144 : 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 14084 : unsigned int num_equidistant_problems = 0;
442 :
443 45776 : 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 31692 : if (!checkRestrictionsForSource(pt, mesh_div, i_from))
447 15660 : continue;
448 :
449 : // TODO: Pre-allocate these two work arrays. They will be regularly resized by the searches
450 16032 : std::vector<std::size_t> return_index(_num_nearest_points + 1);
451 16032 : std::vector<Real> return_dist_sqr(_num_nearest_points + 1);
452 :
453 : // NOTE: app_index is not valid if _group_subapps = true
454 16032 : const auto app_index = i_from / getNumDivisions();
455 16032 : const auto num_search = _num_nearest_points + 1;
456 :
457 16032 : if (_local_kdtrees[i_from]->numberCandidatePoints())
458 : {
459 14568 : _local_kdtrees[i_from]->neighborSearch(pt, num_search, return_index, return_dist_sqr);
460 14568 : 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 14568 : const Point local_pt = getPointInLocalSourceFrame(app_index, pt);
465 :
466 19320 : if (!_nearest_positions_obj &&
467 4752 : _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 14568 : std::vector<std::pair<Real, std::size_t>> zipped_nearest_points;
474 43704 : for (const auto i : make_range(num_found))
475 29136 : zipped_nearest_points.push_back(std::make_pair(return_dist_sqr[i], return_index[i]));
476 14568 : 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 29136 : if (num_found > 1 && num_found == num_search &&
484 14568 : MooseUtils::absoluteFuzzyEqual(zipped_nearest_points[num_found - 1].first,
485 14568 : 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 14568 : Real dist_sum = 0;
496 29136 : for (const auto i : make_range(num_search - 1))
497 : {
498 14568 : auto index = zipped_nearest_points[i].second;
499 14568 : 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 14568 : if (MooseUtils::absoluteFuzzyEqual(dist_sum / return_dist_sqr.size(),
504 14568 : 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 14568 : }
516 16032 : }
517 : }
518 :
519 : // Move to next point
520 6119287 : i_pt++;
521 : }
522 79120 : }
523 :
524 : bool
525 16000 : 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 44704 : for (const auto & i_node : make_range(elem->n_nodes()))
532 : {
533 40672 : const auto & node = elem->node_ptr(i_node);
534 40672 : const auto & node_blocks = mesh.getNodeBlockIds(*node);
535 40672 : std::set<SubdomainID> u;
536 40672 : std::set_intersection(blocks.begin(),
537 : blocks.end(),
538 : node_blocks.begin(),
539 : node_blocks.end(),
540 : std::inserter(u, u.begin()));
541 40672 : if (!u.empty())
542 11968 : return true;
543 40672 : }
544 4032 : return false;
545 : }
546 :
547 : void
548 51420 : 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 102708 : if (!_from_mesh_divisions.empty() ||
553 51288 : (!_use_nearest_app && _nearest_positions_obj && !_group_subapps))
554 418 : _num_sources = _from_problems.size() * getNumDivisions();
555 : // If we group apps, then we only use one tree per division (nearest-position region)
556 51002 : else if (_nearest_positions_obj && _group_subapps)
557 290 : _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 50712 : _num_sources = _from_problems.size();
564 51420 : }
565 :
566 : unsigned int
567 6449171 : 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 6449171 : if (_use_nearest_app)
572 16092 : return kdtree_index;
573 : // We are looping over all the apps that are grouped together
574 6433079 : else if (_group_subapps)
575 25315 : 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 6407764 : return kdtree_index / getNumDivisions();
580 : }
581 :
582 : unsigned int
583 51420 : MultiAppGeneralFieldNearestLocationTransfer::getNumAppsPerTree() const
584 : {
585 51420 : if (_use_nearest_app)
586 286 : return 1;
587 51134 : else if (_group_subapps)
588 290 : return _from_meshes.size();
589 : else
590 50844 : return 1;
591 : }
592 :
593 : unsigned int
594 13275189 : 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 13275189 : if (_use_nearest_app)
600 18288 : return _from_meshes.size();
601 : // Each nearest-position region is a division
602 13256901 : else if (_nearest_positions_obj && !_group_subapps)
603 110522 : return _nearest_positions_obj->getNumPositions(_fe_problem.getCurrentExecuteOnFlag() ==
604 110522 : EXEC_INITIAL);
605 : // Assume all mesh divisions (on each sub-app) has the same number of divisions. This is checked
606 13146379 : else if (!_from_mesh_divisions.empty())
607 747124 : return _from_mesh_divisions[0]->getNumDivisions();
608 : // Grouping subapps or no special mode, we do not subdivide
609 : else
610 12399255 : return 1;
611 : }
612 :
613 : Point
614 14568 : MultiAppGeneralFieldNearestLocationTransfer::getPointInLocalSourceFrame(unsigned int i_from,
615 : const Point & pt) const
616 : {
617 :
618 19320 : if (!_nearest_positions_obj &&
619 4752 : (!_from_transforms[getGlobalSourceAppIndex(i_from)]->hasCoordinateSystemTypeChange() ||
620 0 : _skip_coordinate_collapsing))
621 4752 : return _from_transforms[getGlobalSourceAppIndex(i_from)]->mapBack(pt);
622 9816 : else if (!_nearest_positions_obj || !_group_subapps)
623 4392 : return pt - _from_positions[i_from];
624 : else
625 5424 : return pt;
626 : }
627 :
628 : bool
629 6465537 : 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 6465537 : 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 138581 : auto position_index = i_from; // if _group_subapps
639 138581 : if (_use_nearest_app)
640 26297 : position_index = getGlobalSourceAppIndex(i_from);
641 112284 : else if (!_group_subapps)
642 52456 : 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 138581 : if (!closestToPosition(position_index, pt))
647 72234 : return false;
648 : }
649 :
650 : // Application index depends on which source/grouping mode we are using
651 6393303 : const unsigned int app_index = getAppIndex(i_from, i_from / getNumDivisions());
652 :
653 : // Check mesh restriction before anything
654 6393303 : 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 6393303 : 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 195716 : 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 195716 : if (_from_mesh_division_behavior == MeshDivisionTransferUse::MATCH_SUBAPP_INDEX &&
678 : mesh_div != kd_div_index)
679 16104 : 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 179612 : else if ((_from_mesh_division_behavior == MeshDivisionTransferUse::MATCH_DIVISION_INDEX ||
684 179612 : _to_mesh_division_behavior == MeshDivisionTransferUse::MATCH_DIVISION_INDEX) &&
685 : mesh_div != kd_div_index)
686 133152 : 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 6245497 : if (_to_mesh_division_behavior == MeshDivisionTransferUse::MATCH_SUBAPP_INDEX &&
692 1450 : mesh_div != getGlobalSourceAppIndex(app_index))
693 1050 : return false;
694 :
695 6242997 : return true;
696 : }
|
