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 "GrainTracker.h"
11 :
12 : // MOOSE includes
13 : #include "PolycrystalUserObjectBase.h"
14 : #include "GeneratedMesh.h"
15 : #include "MooseMesh.h"
16 : #include "MooseVariable.h"
17 : #include "NonlinearSystem.h"
18 :
19 : #include "libmesh/periodic_boundary_base.h"
20 :
21 : // C++ includes
22 : #include <algorithm>
23 : #include <limits>
24 : #include <numeric>
25 :
26 : template <>
27 : void
28 8205 : dataStore(std::ostream & stream, GrainTracker::PartialFeatureData & feature, void * context)
29 : {
30 8205 : storeHelper(stream, feature.boundary_intersection, context);
31 8205 : storeHelper(stream, feature.id, context);
32 8205 : storeHelper(stream, feature.centroid, context);
33 8205 : storeHelper(stream, feature.status, context);
34 8205 : }
35 :
36 : template <>
37 : void
38 13130 : dataLoad(std::istream & stream, GrainTracker::PartialFeatureData & feature, void * context)
39 : {
40 13130 : loadHelper(stream, feature.boundary_intersection, context);
41 13130 : loadHelper(stream, feature.id, context);
42 13130 : loadHelper(stream, feature.centroid, context);
43 13130 : loadHelper(stream, feature.status, context);
44 13130 : }
45 :
46 : registerMooseObject("PhaseFieldApp", GrainTracker);
47 :
48 : InputParameters
49 1038 : GrainTracker::validParams()
50 : {
51 1038 : InputParameters params = FeatureFloodCount::validParams();
52 1038 : params += GrainTrackerInterface::validParams();
53 :
54 : // FeatureFloodCount adds a relationship manager, but we need to extend that for GrainTracker
55 : params.clearRelationshipManagers();
56 :
57 2076 : params.addRelationshipManager(
58 : "ElementSideNeighborLayers",
59 : Moose::RelationshipManagerType::GEOMETRIC,
60 :
61 1557 : [](const InputParameters & obj_params, InputParameters & rm_params)
62 1557 : { rm_params.set<unsigned short>("layers") = obj_params.get<unsigned short>("halo_level"); }
63 :
64 : );
65 :
66 2076 : params.addRelationshipManager("ElementSideNeighborLayers",
67 : Moose::RelationshipManagerType::ALGEBRAIC);
68 :
69 : // The GrainTracker requires non-volatile storage for tracking grains across invocations.
70 1038 : params.set<bool>("restartable_required") = true;
71 :
72 2076 : params.addParam<Real>("bound_value",
73 2076 : 0.0,
74 : "Absolute value of the lower bound for the variable value that represents "
75 : "a region not assigned to the grain. Must be positive, but the actual "
76 : "value used is -bound_value.");
77 :
78 1038 : params.addClassDescription("Grain Tracker object for running reduced order parameter simulations "
79 : "without grain coalescence.");
80 :
81 1038 : return params;
82 0 : }
83 :
84 519 : GrainTracker::GrainTracker(const InputParameters & parameters)
85 : : FeatureFloodCount(parameters),
86 : GrainTrackerInterface(),
87 519 : _tracking_step(getParam<int>("tracking_step")),
88 1038 : _halo_level(getParam<unsigned short>("halo_level")),
89 1038 : _max_remap_recursion_depth(getParam<unsigned short>("max_remap_recursion_depth")),
90 1038 : _n_reserve_ops(getParam<unsigned short>("reserve_op")),
91 519 : _reserve_op_index(_n_reserve_ops <= _n_vars ? _n_vars - _n_reserve_ops : 0),
92 1038 : _reserve_op_threshold(getParam<Real>("reserve_op_threshold")),
93 1038 : _bound_value(getParam<Real>("bound_value")),
94 1038 : _remap(getParam<bool>("remap_grains")),
95 1038 : _tolerate_failure(getParam<bool>("tolerate_failure")),
96 519 : _poly_ic_uo(parameters.isParamValid("polycrystal_ic_uo")
97 798 : ? &getUserObject<PolycrystalUserObjectBase>("polycrystal_ic_uo")
98 : : nullptr),
99 1038 : _verbosity_level(getParam<short>("verbosity_level")),
100 1038 : _first_time(declareRestartableData<bool>("first_time", true)),
101 1038 : _error_on_grain_creation(getParam<bool>("error_on_grain_creation")),
102 519 : _reserve_grain_first_index(0),
103 519 : _old_max_grain_id(0),
104 1038 : _max_curr_grain_id(declareRestartableData<unsigned int>("max_curr_grain_id", invalid_id)),
105 1038 : _is_transient(_subproblem.isTransient())
106 : {
107 519 : if (_tolerate_failure)
108 12 : paramInfo("tolerate_failure",
109 : "Tolerate failure has been set to true. Non-physical simulation results "
110 : "are possible, you will be notified in the event of a failed remapping operation.");
111 :
112 519 : if (_tracking_step > 0 && _poly_ic_uo)
113 0 : mooseError("Can't start tracking after the initial condition when using a polycrystal_ic_uo");
114 519 : }
115 :
116 515 : GrainTracker::~GrainTracker() {}
117 :
118 : Real
119 5480451 : GrainTracker::getEntityValue(dof_id_type entity_id,
120 : FieldType field_type,
121 : std::size_t var_index) const
122 : {
123 5480451 : if (_t_step < _tracking_step)
124 : return 0;
125 :
126 5480451 : return FeatureFloodCount::getEntityValue(entity_id, field_type, var_index);
127 : }
128 :
129 : const std::vector<unsigned int> &
130 151355083 : GrainTracker::getVarToFeatureVector(dof_id_type elem_id) const
131 : {
132 151355083 : return FeatureFloodCount::getVarToFeatureVector(elem_id);
133 : }
134 :
135 : unsigned int
136 378 : GrainTracker::getFeatureVar(unsigned int feature_id) const
137 : {
138 378 : return FeatureFloodCount::getFeatureVar(feature_id);
139 : }
140 :
141 : std::size_t
142 0 : GrainTracker::getNumberActiveGrains() const
143 : {
144 : // Note: This value is parallel consistent, see FeatureFloodCount::communicateAndMerge()
145 0 : return _feature_count;
146 : }
147 :
148 : std::size_t
149 853 : GrainTracker::getTotalFeatureCount() const
150 : {
151 : // Note: This value is parallel consistent, see assignGrains()/trackGrains()
152 853 : return _max_curr_grain_id == invalid_id ? 0 : _max_curr_grain_id + 1;
153 : }
154 :
155 : Point
156 12586316 : GrainTracker::getGrainCentroid(unsigned int grain_id) const
157 : {
158 : mooseAssert(grain_id < _feature_id_to_local_index.size(), "Grain ID out of bounds");
159 12586316 : auto grain_index = _feature_id_to_local_index[grain_id];
160 :
161 12586316 : if (grain_index != invalid_size_t)
162 : {
163 : mooseAssert(_feature_id_to_local_index[grain_id] < _feature_sets.size(),
164 : "Grain index out of bounds");
165 : // Note: This value is parallel consistent, see GrainTracker::broadcastAndUpdateGrainData()
166 12586316 : return _feature_sets[_feature_id_to_local_index[grain_id]]._centroid;
167 : }
168 :
169 : // Inactive grain
170 : return Point();
171 : }
172 :
173 : bool
174 378 : GrainTracker::doesFeatureIntersectBoundary(unsigned int feature_id) const
175 : {
176 : // TODO: This data structure may need to be turned into a Multimap
177 : mooseAssert(feature_id < _feature_id_to_local_index.size(), "Grain ID out of bounds");
178 :
179 378 : auto feature_index = _feature_id_to_local_index[feature_id];
180 378 : if (feature_index != invalid_size_t)
181 : {
182 : mooseAssert(feature_index < _feature_sets.size(), "Grain index out of bounds");
183 313 : return _feature_sets[feature_index]._boundary_intersection != BoundaryIntersection::NONE;
184 : }
185 :
186 : return false;
187 : }
188 :
189 : bool
190 378 : GrainTracker::doesFeatureIntersectSpecifiedBoundary(unsigned int feature_id) const
191 : {
192 : // TODO: This data structure may need to be turned into a Multimap
193 : mooseAssert(feature_id < _feature_id_to_local_index.size(), "Grain ID out of bounds");
194 :
195 378 : auto feature_index = _feature_id_to_local_index[feature_id];
196 378 : if (feature_index != invalid_size_t)
197 : {
198 : mooseAssert(feature_index < _feature_sets.size(), "Grain index out of bounds");
199 313 : return ((_feature_sets[feature_index]._boundary_intersection &
200 313 : BoundaryIntersection::SPECIFIED_BOUNDARY) == BoundaryIntersection::SPECIFIED_BOUNDARY);
201 : }
202 :
203 : return false;
204 : }
205 :
206 : bool
207 378 : GrainTracker::isFeaturePercolated(unsigned int feature_id) const
208 : {
209 : // TODO: This data structure may need to be turned into a Multimap
210 : mooseAssert(feature_id < _feature_id_to_local_index.size(), "Grain ID out of bounds");
211 :
212 378 : auto feature_index = _feature_id_to_local_index[feature_id];
213 378 : if (feature_index != invalid_size_t)
214 : {
215 : mooseAssert(feature_index < _feature_sets.size(), "Grain index out of bounds");
216 313 : bool primary = ((_feature_sets[feature_index]._boundary_intersection &
217 : BoundaryIntersection::PRIMARY_PERCOLATION_BOUNDARY) ==
218 313 : BoundaryIntersection::PRIMARY_PERCOLATION_BOUNDARY);
219 : bool secondary = ((_feature_sets[feature_index]._boundary_intersection &
220 : BoundaryIntersection::SECONDARY_PERCOLATION_BOUNDARY) ==
221 313 : BoundaryIntersection::SECONDARY_PERCOLATION_BOUNDARY);
222 313 : return (primary && secondary);
223 : }
224 :
225 : return false;
226 : }
227 :
228 : void
229 2098 : GrainTracker::initialize()
230 : {
231 : // Don't track grains if the current simulation step is before the specified tracking step
232 2098 : if (_t_step < _tracking_step)
233 : return;
234 :
235 : /**
236 : * If we are passed the first time, we need to save the existing grains before beginning the
237 : * tracking on the current step. We'll do that with a swap since the _feature_sets contents will
238 : * be cleared anyway.
239 : */
240 2098 : if (!_first_time)
241 1738 : _feature_sets_old.swap(_feature_sets);
242 :
243 2098 : FeatureFloodCount::initialize();
244 : }
245 :
246 : void
247 675 : GrainTracker::meshChanged()
248 : {
249 : // Update the element ID ranges for use when computing halo maps
250 675 : if (_compute_halo_maps && _mesh.isDistributedMesh())
251 : {
252 21 : _all_ranges.clear();
253 :
254 21 : auto range = std::make_pair(std::numeric_limits<dof_id_type>::max(),
255 : std::numeric_limits<dof_id_type>::min());
256 45354 : for (const auto & current_elem : _mesh.getMesh().active_local_element_ptr_range())
257 : {
258 22656 : auto id = current_elem->id();
259 22656 : if (id < range.first)
260 21 : range.first = id;
261 22635 : else if (id > range.second)
262 22635 : range.second = id;
263 21 : }
264 :
265 21 : _communicator.gather(0, range, _all_ranges);
266 : }
267 :
268 675 : FeatureFloodCount::meshChanged();
269 675 : }
270 :
271 : void
272 2098 : GrainTracker::execute()
273 : {
274 : // Don't track grains if the current simulation step is before the specified tracking step
275 2098 : if (_t_step < _tracking_step)
276 : return;
277 :
278 2098 : if (_poly_ic_uo && _first_time)
279 : return;
280 :
281 1909 : FeatureFloodCount::execute();
282 : }
283 :
284 : Real
285 3749247 : GrainTracker::getThreshold(std::size_t var_index) const
286 : {
287 : // If we are inspecting a reserve op parameter, we need to make sure
288 : // that there is an entity above the reserve_op threshold before
289 : // starting the flood of the feature.
290 3749247 : if (var_index >= _reserve_op_index)
291 118988 : return _reserve_op_threshold;
292 : else
293 3630259 : return _step_threshold;
294 : }
295 :
296 : void
297 189 : GrainTracker::prepopulateState(const FeatureFloodCount & ffc_object)
298 : {
299 : mooseAssert(_first_time, "This method should only be called on the first invocation");
300 :
301 189 : _feature_sets.clear();
302 :
303 : /**
304 : * The minimum information needed to bootstrap the GrainTracker is as follows:
305 : * _feature_sets
306 : * _feature_count
307 : */
308 189 : if (_is_primary)
309 : {
310 : const auto & features = ffc_object.getFeatures();
311 1811 : for (auto & feature : features)
312 1708 : _feature_sets.emplace_back(feature.duplicate());
313 :
314 103 : _feature_count = _feature_sets.size();
315 : }
316 : else
317 : {
318 : const auto & features = ffc_object.getFeatures();
319 : _partial_feature_sets[0].clear();
320 637 : for (auto & feature : features)
321 551 : _partial_feature_sets[0].emplace_back(feature.duplicate());
322 : }
323 :
324 : // Make sure that feature count is communicated to all ranks
325 189 : _communicator.broadcast(_feature_count);
326 189 : }
327 :
328 : void
329 2098 : GrainTracker::finalize()
330 : {
331 : // Don't track grains if the current simulation step is before the specified tracking step
332 2098 : if (_t_step < _tracking_step)
333 0 : return;
334 :
335 4196 : TIME_SECTION("finalize", 3, "Finalizing GrainTracker");
336 :
337 : // Expand the depth of the halos around all grains
338 2098 : auto num_halo_layers = _halo_level >= 1
339 2098 : ? _halo_level - 1
340 : : 0; // The first level of halos already exists so subtract one
341 :
342 2098 : if (_poly_ic_uo && _first_time)
343 189 : prepopulateState(*_poly_ic_uo);
344 : else
345 : {
346 1909 : expandEdgeHalos(num_halo_layers);
347 :
348 : // Build up the grain map on the root processor
349 1909 : communicateAndMerge();
350 : }
351 :
352 : /**
353 : * Assign or Track Grains
354 : */
355 2098 : if (_first_time)
356 360 : assignGrains();
357 : else
358 1738 : trackGrains();
359 :
360 2098 : if (_verbosity_level > 1)
361 15 : _console << "Finished inside of trackGrains" << std::endl;
362 :
363 : /**
364 : * Broadcast essential data
365 : */
366 2098 : broadcastAndUpdateGrainData();
367 :
368 : /**
369 : * Remap Grains
370 : */
371 2098 : if (_remap)
372 2090 : remapGrains();
373 :
374 2098 : updateFieldInfo();
375 2098 : if (_verbosity_level > 1)
376 15 : _console << "Finished inside of updateFieldInfo" << std::endl;
377 :
378 : // Set the first time flag false here (after all methods of finalize() have completed)
379 2098 : _first_time = false;
380 :
381 : // TODO: Release non essential memory
382 2098 : if (_verbosity_level > 0)
383 2098 : _console << "Finished inside of GrainTracker\n" << std::endl;
384 2098 : }
385 :
386 : void
387 2098 : GrainTracker::broadcastAndUpdateGrainData()
388 : {
389 4196 : TIME_SECTION("broadcastAndUpdateGrainData", 3, "Broadcasting and Updating Grain Data");
390 :
391 : std::vector<PartialFeatureData> root_feature_data;
392 2098 : std::vector<std::string> send_buffer(1), recv_buffer;
393 :
394 2098 : if (_is_primary)
395 : {
396 873 : root_feature_data.reserve(_feature_sets.size());
397 :
398 : // Populate a subset of the information in a small data structure
399 873 : std::transform(_feature_sets.begin(),
400 873 : _feature_sets.end(),
401 : std::back_inserter(root_feature_data),
402 : [](FeatureData & feature)
403 : {
404 : PartialFeatureData partial_feature;
405 8205 : partial_feature.boundary_intersection = feature._boundary_intersection;
406 8205 : partial_feature.id = feature._id;
407 8205 : partial_feature.centroid = feature._centroid;
408 8205 : partial_feature.status = feature._status;
409 : return partial_feature;
410 : });
411 :
412 873 : std::ostringstream oss;
413 873 : dataStore(oss, root_feature_data, this);
414 873 : send_buffer[0].assign(oss.str());
415 873 : }
416 :
417 : // Broadcast the data to all ranks
418 2098 : _communicator.broadcast_packed_range((void *)(nullptr),
419 : send_buffer.begin(),
420 : send_buffer.end(),
421 : (void *)(nullptr),
422 : std::back_inserter(recv_buffer));
423 :
424 : // Unpack and update
425 2098 : if (!_is_primary)
426 : {
427 1225 : std::istringstream iss;
428 : iss.str(recv_buffer[0]);
429 1225 : iss.clear();
430 :
431 1225 : dataLoad(iss, root_feature_data, this);
432 :
433 14355 : for (const auto & partial_data : root_feature_data)
434 : {
435 : // See if this processor has a record of this grain
436 13130 : if (partial_data.id < _feature_id_to_local_index.size() &&
437 13130 : _feature_id_to_local_index[partial_data.id] != invalid_size_t)
438 : {
439 6176 : auto & grain = _feature_sets[_feature_id_to_local_index[partial_data.id]];
440 6176 : grain._boundary_intersection = partial_data.boundary_intersection;
441 6176 : grain._centroid = partial_data.centroid;
442 6176 : if (partial_data.status == Status::INACTIVE)
443 0 : grain._status = Status::INACTIVE;
444 : }
445 : }
446 1225 : }
447 2098 : }
448 :
449 : void
450 360 : GrainTracker::assignGrains()
451 : {
452 : mooseAssert(_first_time, "assignGrains may only be called on the first tracking step");
453 :
454 : /**
455 : * We need to assign grainIDs to get the simulation going. We'll use the default sorting that
456 : * doesn't require valid grainIDs (relies on _min_entity_id and _var_index). These will be the
457 : * unique grain numbers that we must track for remainder of the simulation.
458 : */
459 360 : if (_is_primary)
460 : {
461 : // Find the largest grain ID, this requires sorting if the ID is not already set
462 207 : sortAndLabel();
463 :
464 207 : if (_feature_sets.empty())
465 : {
466 5 : _max_curr_grain_id = invalid_id;
467 5 : _reserve_grain_first_index = 0;
468 : }
469 : else
470 : {
471 202 : _max_curr_grain_id = _feature_sets.back()._id;
472 202 : _reserve_grain_first_index = _max_curr_grain_id + 1;
473 : }
474 :
475 2491 : for (auto & grain : _feature_sets)
476 2284 : grain._status = Status::MARKED; // Mark the grain
477 :
478 : } // is_primary
479 :
480 : /*************************************************************
481 : ****************** COLLECTIVE WORK SECTION ******************
482 : *************************************************************/
483 :
484 : // Make IDs on all non-primary ranks consistent
485 360 : scatterAndUpdateRanks();
486 :
487 : // Build up an id to index map
488 360 : _communicator.broadcast(_max_curr_grain_id);
489 360 : buildFeatureIdToLocalIndices(_max_curr_grain_id);
490 :
491 : // Now trigger the newGrainCreated() callback on all ranks
492 360 : if (_max_curr_grain_id != invalid_id)
493 4003 : for (unsigned int new_id = 0; new_id <= _max_curr_grain_id; ++new_id)
494 3649 : newGrainCreated(new_id);
495 360 : }
496 :
497 : void
498 1738 : GrainTracker::trackGrains()
499 : {
500 3476 : TIME_SECTION("trackGrains", 3, "Tracking Grains");
501 :
502 : mooseAssert(!_first_time, "Track grains may only be called when _tracking_step > _t_step");
503 :
504 : // Used to track indices for which to trigger the new grain callback on (used on all ranks)
505 1738 : auto _old_max_grain_id = _max_curr_grain_id;
506 :
507 : /**
508 : * Only the primary rank does tracking, the remaining ranks
509 : * wait to receive local to global indices from the primary.
510 : */
511 1738 : if (_is_primary)
512 : {
513 : // Reset Status on active unique grains
514 666 : std::vector<unsigned int> map_sizes(_maps_size);
515 6633 : for (auto & grain : _feature_sets_old)
516 : {
517 5967 : if (grain._status != Status::INACTIVE)
518 : {
519 5967 : grain._status = Status::CLEAR;
520 5967 : map_sizes[grain._var_index]++;
521 : }
522 : }
523 :
524 : // Print out stats on overall tracking changes per var_index
525 666 : if (_verbosity_level > 0)
526 : {
527 666 : _console << "\nGrain Tracker Status:";
528 4798 : for (const auto map_num : make_range(_maps_size))
529 : {
530 4132 : _console << "\nGrains active index " << map_num << ": " << map_sizes[map_num] << " -> "
531 4132 : << _feature_counts_per_map[map_num];
532 4132 : if (map_sizes[map_num] > _feature_counts_per_map[map_num])
533 61 : _console << "--";
534 4071 : else if (map_sizes[map_num] < _feature_counts_per_map[map_num])
535 15 : _console << "++";
536 : }
537 666 : _console << '\n' << std::endl;
538 : }
539 :
540 : // Before we track grains, lets sort them so that we get parallel consistent answers
541 666 : std::sort(_feature_sets.begin(), _feature_sets.end());
542 :
543 : /**
544 : * To track grains across time steps, we will loop over our unique grains and link each one up
545 : * with one of our new unique grains. The criteria for doing this will be to find the unique
546 : * grain in the new list with a matching variable index whose centroid is closest to this
547 : * unique grain.
548 : */
549 666 : std::vector<std::size_t> new_grain_index_to_existing_grain_index(_feature_sets.size(),
550 666 : invalid_size_t);
551 :
552 6633 : for (const auto old_grain_index : index_range(_feature_sets_old))
553 : {
554 : auto & old_grain = _feature_sets_old[old_grain_index];
555 :
556 5967 : if (old_grain._status == Status::INACTIVE) // Don't try to find matches for inactive grains
557 : continue;
558 :
559 : std::size_t closest_match_index = invalid_size_t;
560 : Real min_centroid_diff = std::numeric_limits<Real>::max();
561 :
562 : /**
563 : * The _feature_sets vector is constructed by _var_index so we can avoid looping over all
564 : * indices. We can quickly jump to the first matching index to reduce the number of
565 : * comparisons and terminate our loop when our variable index stops matching.
566 : */
567 : // clang-format off
568 : auto start_it =
569 5967 : std::lower_bound(_feature_sets.begin(), _feature_sets.end(), old_grain._var_index,
570 : [](const FeatureData & item, std::size_t var_index)
571 : {
572 22539 : return item._var_index < var_index;
573 : });
574 : // clang-format on
575 :
576 : // We only need to examine grains that have matching variable indices
577 : bool any_boxes_intersect = false;
578 17374 : for (MooseIndex(_feature_sets)
579 5967 : new_grain_index = std::distance(_feature_sets.begin(), start_it);
580 17374 : new_grain_index < _feature_sets.size() &&
581 16562 : _feature_sets[new_grain_index]._var_index == old_grain._var_index;
582 : ++new_grain_index)
583 : {
584 : auto & new_grain = _feature_sets[new_grain_index];
585 :
586 : /**
587 : * Don't try to do any matching unless the bounding boxes at least overlap. This is to avoid
588 : * the corner case of having a grain split and a grain disappear during the same time step!
589 : */
590 11407 : if (new_grain.boundingBoxesIntersect(old_grain))
591 : {
592 : any_boxes_intersect = true;
593 8442 : Real curr_centroid_diff = centroidRegionDistance(old_grain._bboxes, new_grain._bboxes);
594 8442 : if (curr_centroid_diff <= min_centroid_diff)
595 : {
596 : closest_match_index = new_grain_index;
597 : min_centroid_diff = curr_centroid_diff;
598 : }
599 : }
600 : }
601 :
602 5967 : if (_verbosity_level > 2 && !any_boxes_intersect)
603 0 : _console << "\nNo intersecting bounding boxes found while trying to match grain "
604 0 : << old_grain;
605 :
606 : // found a match
607 5967 : if (closest_match_index != invalid_size_t)
608 : {
609 : /**
610 : * It's possible that multiple existing grains will map to a single new grain (indicated by
611 : * finding multiple matches when we are building this map). This will happen any time a
612 : * grain disappears during this time step. We need to figure out the rightful owner in this
613 : * case and inactivate the old grain.
614 : */
615 5912 : auto curr_index = new_grain_index_to_existing_grain_index[closest_match_index];
616 5912 : if (curr_index != invalid_size_t)
617 : {
618 : // The new feature being competed for
619 9 : auto & new_grain = _feature_sets[closest_match_index];
620 :
621 : // The other old grain competing to match up to the same new grain
622 : auto & other_old_grain = _feature_sets_old[curr_index];
623 :
624 9 : auto centroid_diff1 = centroidRegionDistance(new_grain._bboxes, old_grain._bboxes);
625 9 : auto centroid_diff2 = centroidRegionDistance(new_grain._bboxes, other_old_grain._bboxes);
626 :
627 9 : auto & inactive_grain = (centroid_diff1 < centroid_diff2) ? other_old_grain : old_grain;
628 :
629 9 : inactive_grain._status = Status::INACTIVE;
630 9 : if (_verbosity_level > 0)
631 : {
632 9 : _console << COLOR_GREEN << "Marking Grain " << inactive_grain._id
633 9 : << " as INACTIVE (variable index: " << inactive_grain._var_index << ")\n"
634 9 : << COLOR_DEFAULT;
635 9 : if (_verbosity_level > 1)
636 0 : _console << inactive_grain;
637 : }
638 :
639 : /**
640 : * If the grain we just marked inactive was the one whose index was in the new grain
641 : * to existing grain map (other_old_grain). Then we need to update the map to point
642 : * to the new match winner.
643 : */
644 9 : if (&inactive_grain == &other_old_grain)
645 0 : new_grain_index_to_existing_grain_index[closest_match_index] = old_grain_index;
646 : }
647 : else
648 5903 : new_grain_index_to_existing_grain_index[closest_match_index] = old_grain_index;
649 : }
650 : }
651 :
652 : // Mark all resolved grain matches
653 6587 : for (const auto new_index : index_range(new_grain_index_to_existing_grain_index))
654 : {
655 5921 : auto curr_index = new_grain_index_to_existing_grain_index[new_index];
656 :
657 : // This may be a new grain, we'll handle that case below
658 5921 : if (curr_index == invalid_size_t)
659 18 : continue;
660 :
661 : mooseAssert(_feature_sets_old[curr_index]._id != invalid_id,
662 : "Invalid ID in old grain structure");
663 :
664 5903 : _feature_sets[new_index]._id = _feature_sets_old[curr_index]._id; // Transfer ID
665 5903 : _feature_sets[new_index]._status = Status::MARKED; // Mark the status in the new set
666 5903 : _feature_sets_old[curr_index]._status = Status::MARKED; // Mark the status in the old set
667 : }
668 :
669 : /**
670 : * At this point we have should have only two cases left to handle:
671 : * Case 1: A grain in the new set who has an unset status (These are new grains, previously
672 : * untracked) This case is easy to understand. Since we are matching up grains by
673 : * looking at the old set and finding closest matches in the new set, any grain in
674 : * the new set that isn't matched up is simply new since some other grain satisfied
675 : * each and every request from the old set.
676 : *
677 : * Case 2: A grain in the old set who has an unset status (These are inactive grains that
678 : * haven't been marked) We can only fall into this case when the very last grain on
679 : * a given variable disappears during the current time step. In that case we never have
680 : * a matching _var_index in the comparison loop above so that old grain never competes
681 : * for any new grain which means it can't be marked inactive in the loop above.
682 : */
683 : // Case 1 (new grains in _feature_sets):
684 6587 : for (const auto grain_num : index_range(_feature_sets))
685 : {
686 5921 : auto & grain = _feature_sets[grain_num];
687 :
688 : // New Grain
689 5921 : if (grain._status == Status::CLEAR)
690 : {
691 : /**
692 : * Now we need to figure out what kind of "new" grain this is. Is it a nucleating grain that
693 : * we're just barely seeing for the first time or is it a "splitting" grain. A grain that
694 : * gets pinched into two or more pieces usually as it is being absorbed by other grains or
695 : * possibly due to external forces. We have to handle splitting grains this way so as to
696 : * no confuse them with regular grains that just happen to be in contact in this step.
697 : *
698 : * Splitting Grain: An grain that is unmatched by any old grain
699 : * on the same order parameter with touching halos.
700 : *
701 : * Nucleating Grain: A completely new grain appearing somewhere in the domain
702 : * not overlapping any other grain's halo.
703 : *
704 : * To figure out which case we are dealing with, we have to make another pass over all of
705 : * the existing grains with matching variable indices to see if any of them have overlapping
706 : * halos.
707 : */
708 :
709 : // clang-format off
710 : auto start_it =
711 18 : std::lower_bound(_feature_sets.begin(), _feature_sets.end(), grain._var_index,
712 : [](const FeatureData & item, std::size_t var_index)
713 : {
714 42 : return item._var_index < var_index;
715 : });
716 : // clang-format on
717 :
718 : // Loop over matching variable indices
719 50 : for (MooseIndex(_feature_sets)
720 18 : new_grain_index = std::distance(_feature_sets.begin(), start_it);
721 50 : new_grain_index < _feature_sets.size() &&
722 40 : _feature_sets[new_grain_index]._var_index == grain._var_index;
723 : ++new_grain_index)
724 : {
725 : auto & other_grain = _feature_sets[new_grain_index];
726 :
727 : // Splitting grain?
728 14 : if (grain_num != new_grain_index && // Make sure indices aren't pointing at the same grain
729 28 : other_grain._status == Status::MARKED && // and that the other grain is indeed marked
730 54 : other_grain.boundingBoxesIntersect(grain) && // and the bboxes intersect
731 8 : other_grain.halosIntersect(grain)) // and the halos also intersect
732 : // TODO: Inspect combined volume and see if it's "close" to the expected value
733 : {
734 8 : grain._id = other_grain._id; // Set the duplicate ID
735 8 : grain._status = Status::MARKED; // Mark it
736 :
737 8 : if (_verbosity_level > 0)
738 8 : _console << COLOR_YELLOW << "Split Grain Detected #" << grain._id
739 8 : << " (variable index: " << grain._var_index << ")\n"
740 8 : << COLOR_DEFAULT;
741 8 : if (_verbosity_level > 1)
742 0 : _console << grain << other_grain;
743 : }
744 : }
745 :
746 18 : if (grain._var_index < _reserve_op_index)
747 : {
748 : /**
749 : * The "try-harder loop":
750 : * OK so we still have an extra grain in the new set that isn't matched up against the
751 : * old set and since the order parameter isn't reserved. We aren't really expecting a new
752 : * grain. Let's try to make a few more attempts to see if this is a split grain even
753 : * though it failed to match the criteria above. This might happen if the halo front is
754 : * advancing too fast!
755 : *
756 : * In this loop we'll make an attempt to match up this new grain to the old halos. If
757 : * adaptivity is happening this could fail as elements in the new set may be at a
758 : * different level than in the old set. If we get multiple matches, we'll compare the
759 : * grain volumes (based on elements, not integrated to choose the closest).
760 : *
761 : * Future ideas:
762 : * Look at the volume fraction of the new grain and overlay it over the volume fraction
763 : * of the old grain (would require more saved information, or an aux field hanging around
764 : * (subject to projection problems).
765 : */
766 8 : if (_verbosity_level > 1)
767 0 : _console << COLOR_YELLOW
768 : << "Trying harder to detect a split grain while examining grain on variable "
769 0 : "index "
770 0 : << grain._var_index << '\n'
771 0 : << COLOR_DEFAULT;
772 :
773 : std::vector<std::size_t> old_grain_indices;
774 56 : for (const auto old_grain_index : index_range(_feature_sets_old))
775 : {
776 : auto & old_grain = _feature_sets_old[old_grain_index];
777 :
778 48 : if (old_grain._status == Status::INACTIVE)
779 0 : continue;
780 :
781 : /**
782 : * Note that the old grains we are looking at will already be marked from the earlier
783 : * tracking phase. We are trying to see if this unmatched grain is part of a larger
784 : * whole. To do that we'll look at the halos across the time step.
785 : */
786 17 : if (grain._var_index == old_grain._var_index &&
787 56 : grain.boundingBoxesIntersect(old_grain) && grain.halosIntersect(old_grain))
788 8 : old_grain_indices.push_back(old_grain_index);
789 : }
790 :
791 8 : if (old_grain_indices.size() == 1)
792 : {
793 8 : grain._id = _feature_sets_old[old_grain_indices[0]]._id;
794 8 : grain._status = Status::MARKED;
795 :
796 8 : if (_verbosity_level > 0)
797 8 : _console << COLOR_YELLOW << "Split Grain Detected #" << grain._id
798 8 : << " (variable index: " << grain._var_index << ")\n"
799 8 : << COLOR_DEFAULT;
800 : }
801 0 : else if (old_grain_indices.size() > 1)
802 : _console
803 0 : << COLOR_RED << "Split Grain Likely Detected #" << grain._id
804 0 : << " Need more information to find correct candidate - contact a developer!\n\n"
805 0 : << COLOR_DEFAULT;
806 8 : }
807 :
808 : // Must be a nucleating grain (status is still not set)
809 18 : if (grain._status == Status::CLEAR)
810 : {
811 10 : auto new_index = getNextUniqueID();
812 10 : grain._id = new_index; // Set the ID
813 10 : grain._status = Status::MARKED; // Mark it
814 :
815 10 : if (_verbosity_level > 0)
816 10 : _console << COLOR_YELLOW << "Nucleating Grain Detected "
817 10 : << " (variable index: " << grain._var_index << ")\n"
818 10 : << COLOR_DEFAULT;
819 10 : if (_verbosity_level > 1)
820 0 : _console << grain;
821 : }
822 : }
823 : }
824 :
825 : // Case 2 (inactive grains in _feature_sets_old)
826 6633 : for (auto & grain : _feature_sets_old)
827 : {
828 5967 : if (grain._status == Status::CLEAR)
829 : {
830 55 : grain._status = Status::INACTIVE;
831 55 : if (_verbosity_level > 0)
832 : {
833 55 : _console << COLOR_GREEN << "Marking Grain " << grain._id
834 55 : << " as INACTIVE (variable index: " << grain._var_index << ")\n"
835 55 : << COLOR_DEFAULT;
836 55 : if (_verbosity_level > 1)
837 0 : _console << grain;
838 : }
839 : }
840 : }
841 666 : } // is_primary
842 :
843 : /*************************************************************
844 : ****************** COLLECTIVE WORK SECTION ******************
845 : *************************************************************/
846 :
847 : // Make IDs on all non-primary ranks consistent
848 1738 : scatterAndUpdateRanks();
849 :
850 : // Build up an id to index map
851 1738 : _communicator.broadcast(_max_curr_grain_id);
852 1738 : buildFeatureIdToLocalIndices(_max_curr_grain_id);
853 :
854 : /**
855 : * Trigger callback for new grains
856 : */
857 1738 : if (_old_max_grain_id < _max_curr_grain_id)
858 : {
859 17 : for (auto new_id = _old_max_grain_id + 1; new_id <= _max_curr_grain_id; ++new_id)
860 : {
861 : // Don't trigger the callback on the reserve IDs
862 11 : if (new_id >= _reserve_grain_first_index + _n_reserve_ops)
863 : {
864 : // See if we've been instructed to terminate with an error
865 7 : if (!_first_time && _error_on_grain_creation)
866 0 : mooseError(
867 : "Error: New grain detected and \"error_on_new_grain_creation\" is set to true");
868 : else
869 7 : newGrainCreated(new_id);
870 : }
871 : }
872 : }
873 1738 : }
874 :
875 : void
876 3602 : GrainTracker::newGrainCreated(unsigned int new_grain_id)
877 : {
878 3602 : if (!_first_time && _is_primary)
879 : {
880 : mooseAssert(new_grain_id < _feature_id_to_local_index.size(), "new_grain_id is out of bounds");
881 5 : auto grain_index = _feature_id_to_local_index[new_grain_id];
882 : mooseAssert(grain_index != invalid_size_t && grain_index < _feature_sets.size(),
883 : "new_grain_id appears to be invalid");
884 :
885 5 : const auto & grain = _feature_sets[grain_index];
886 5 : _console << COLOR_YELLOW
887 5 : << "\n*****************************************************************************"
888 5 : << "\nCouldn't find a matching grain while working on variable index: "
889 5 : << grain._var_index << "\nCreating new unique grain: " << new_grain_id << '\n'
890 5 : << grain
891 5 : << "\n*****************************************************************************\n"
892 5 : << COLOR_DEFAULT;
893 : }
894 3602 : }
895 :
896 : std::vector<unsigned int>
897 0 : GrainTracker::getNewGrainIDs() const
898 : {
899 0 : std::vector<unsigned int> new_ids(_max_curr_grain_id - _old_max_grain_id);
900 0 : auto new_id = _old_max_grain_id + 1;
901 :
902 : // Generate the new ids
903 : std::iota(new_ids.begin(), new_ids.end(), new_id);
904 :
905 0 : return new_ids;
906 : }
907 :
908 : void
909 2090 : GrainTracker::remapGrains()
910 : {
911 : // Don't remap grains if the current simulation step is before the specified tracking step
912 2090 : if (_t_step < _tracking_step)
913 0 : return;
914 :
915 4180 : TIME_SECTION("remapGrains", 3, "Remapping Grains");
916 :
917 2090 : if (_verbosity_level > 1)
918 15 : _console << "Running remap Grains\n" << std::endl;
919 :
920 : /**
921 : * Map used for communicating remap indices to all ranks
922 : * This map isn't populated until after the remap loop.
923 : * It's declared here before we enter the root scope
924 : * since it's needed by all ranks during the broadcast.
925 : */
926 : std::map<unsigned int, std::size_t> grain_id_to_new_var;
927 :
928 : // Items are added to this list when split grains are found
929 : std::list<std::pair<std::size_t, std::size_t>> split_pairs;
930 :
931 : /**
932 : * The remapping algorithm is recursive. We will use the status variable in each FeatureData
933 : * to track which grains are currently being remapped so we don't have runaway recursion.
934 : * To begin we need to clear all of the active (MARKED) flags (CLEAR).
935 : *
936 : * Additionally we need to record each grain's variable index so that we can communicate
937 : * changes to the non-root ranks later in a single batch.
938 : */
939 2090 : if (_is_primary)
940 : {
941 : // Build the map to detect difference in _var_index mappings after the remap operation
942 : std::map<unsigned int, std::size_t> grain_id_to_existing_var_index;
943 8822 : for (auto & grain : _feature_sets)
944 : {
945 : // Unmark the grain so it can be used in the remap loop
946 7956 : grain._status = Status::CLEAR;
947 :
948 7956 : grain_id_to_existing_var_index[grain._id] = grain._var_index;
949 : }
950 :
951 : // Make sure that all split pieces of any grain are on the same OP
952 8822 : for (const auto i : index_range(_feature_sets))
953 : {
954 7956 : auto & grain1 = _feature_sets[i];
955 :
956 124014 : for (const auto j : index_range(_feature_sets))
957 : {
958 116058 : auto & grain2 = _feature_sets[j];
959 116058 : if (i == j)
960 7956 : continue;
961 :
962 : // The first condition below is there to prevent symmetric checks (duplicate values)
963 108102 : if (i < j && grain1._id == grain2._id)
964 : {
965 9 : split_pairs.push_front(std::make_pair(i, j));
966 9 : if (grain1._var_index != grain2._var_index)
967 : {
968 0 : if (_verbosity_level > 0)
969 0 : _console << COLOR_YELLOW << "Split Grain (#" << grain1._id
970 0 : << ") detected on unmatched OPs (" << grain1._var_index << ", "
971 0 : << grain2._var_index << ") attempting to remap to " << grain1._var_index
972 0 : << ".\n"
973 0 : << COLOR_DEFAULT;
974 :
975 : /**
976 : * We're not going to try very hard to look for a suitable remapping. Just set it to
977 : * what we want and hope it all works out. Make the GrainTracker great again!
978 : */
979 0 : grain1._var_index = grain2._var_index;
980 : grain1._status |= Status::DIRTY;
981 : }
982 : }
983 : }
984 : }
985 :
986 : /**
987 : * Loop over each grain and see if any grains represented by the same variable are "touching"
988 : */
989 : bool grains_remapped;
990 :
991 : std::set<unsigned int> notify_ids;
992 980 : do
993 : {
994 : grains_remapped = false;
995 : notify_ids.clear();
996 :
997 10228 : for (auto & grain1 : _feature_sets)
998 : {
999 : // We need to remap any grains represented on any variable index above the cuttoff
1000 9248 : if (grain1._var_index >= _reserve_op_index)
1001 : {
1002 10 : if (_verbosity_level > 0)
1003 10 : _console << COLOR_YELLOW << "\nGrain #" << grain1._id
1004 10 : << " detected on a reserved order parameter #" << grain1._var_index
1005 10 : << ", remapping to another variable\n"
1006 10 : << COLOR_DEFAULT;
1007 :
1008 10 : for (const auto max : make_range(0, _max_remap_recursion_depth + 1))
1009 10 : if (max < _max_remap_recursion_depth)
1010 : {
1011 10 : if (attemptGrainRenumber(grain1, 0, max))
1012 : break;
1013 : }
1014 0 : else if (!attemptGrainRenumber(grain1, 0, max))
1015 : {
1016 0 : _console << std::flush;
1017 0 : std::stringstream oss;
1018 : oss << "Unable to find any suitable order parameters for remapping while working "
1019 0 : << "with Grain #" << grain1._id << ", which is on a reserve order parameter.\n"
1020 : << "\n\nPossible Resolutions:\n"
1021 : << "\t- Add more order parameters to your simulation (8 for 2D, 28 for 3D)\n"
1022 : << "\t- Increase adaptivity or reduce your grain boundary widths\n"
1023 0 : << "\t- Make sure you are not starting with too many grains for the mesh size\n";
1024 0 : mooseError(oss.str());
1025 0 : }
1026 :
1027 : grains_remapped = true;
1028 : }
1029 :
1030 141004 : for (auto & grain2 : _feature_sets)
1031 : {
1032 : // Don't compare a grain with itself and don't try to remap inactive grains
1033 131756 : if (&grain1 == &grain2)
1034 9248 : continue;
1035 :
1036 134314 : if (grain1._var_index == grain2._var_index && // grains represented by same variable?
1037 23594 : grain1._id != grain2._id && // are they part of different grains?
1038 138985 : grain1.boundingBoxesIntersect(grain2) && // do bboxes intersect (coarse level)?
1039 4689 : grain1.halosIntersect(grain2)) // do they actually overlap (fine level)?
1040 : {
1041 180 : if (_verbosity_level > 0)
1042 180 : _console << COLOR_YELLOW << "Grain #" << grain1._id << " intersects Grain #"
1043 180 : << grain2._id << " (variable index: " << grain1._var_index << ")\n"
1044 180 : << COLOR_DEFAULT;
1045 :
1046 367 : for (const auto max : make_range(0, _max_remap_recursion_depth + 1))
1047 : {
1048 346 : if (max < _max_remap_recursion_depth)
1049 : {
1050 325 : if (attemptGrainRenumber(grain1, 0, max))
1051 : {
1052 : grains_remapped = true;
1053 : break;
1054 : }
1055 : }
1056 42 : else if (!attemptGrainRenumber(grain1, 0, max) &&
1057 21 : !attemptGrainRenumber(grain2, 0, max))
1058 : {
1059 15 : notify_ids.insert(grain1._id);
1060 15 : notify_ids.insert(grain2._id);
1061 : }
1062 : }
1063 : }
1064 : }
1065 : }
1066 : } while (grains_remapped);
1067 :
1068 866 : if (!notify_ids.empty())
1069 : {
1070 3 : _console << std::flush;
1071 3 : std::stringstream oss;
1072 : oss << "Unable to find any suitable order parameters for remapping while working "
1073 : << "with the following grain IDs:\n"
1074 6 : << Moose::stringify(notify_ids, ", ", "", true) << "\n\nPossible Resolutions:\n"
1075 : << "\t- Add more order parameters to your simulation (8 for 2D, 28 for 3D)\n"
1076 : << "\t- Increase adaptivity or reduce your grain boundary widths\n"
1077 6 : << "\t- Make sure you are not starting with too many grains for the mesh size\n";
1078 :
1079 3 : if (_tolerate_failure)
1080 3 : mooseWarning(oss.str());
1081 : else
1082 0 : mooseError(oss.str());
1083 3 : }
1084 :
1085 : // Verify that split grains are still intact
1086 : for (auto & split_pair : split_pairs)
1087 : if (_feature_sets[split_pair.first]._var_index != _feature_sets[split_pair.first]._var_index)
1088 : mooseError("Split grain remapped - This case is currently not handled");
1089 :
1090 : /**
1091 : * The remapping loop is complete but only on the primary process.
1092 : * Now we need to build the remap map and communicate it to the
1093 : * remaining processors.
1094 : */
1095 8822 : for (auto & grain : _feature_sets)
1096 : {
1097 : mooseAssert(grain_id_to_existing_var_index.find(grain._id) !=
1098 : grain_id_to_existing_var_index.end(),
1099 : "Missing unique ID");
1100 :
1101 7956 : auto old_var_index = grain_id_to_existing_var_index[grain._id];
1102 :
1103 7956 : if (old_var_index != grain._var_index)
1104 : {
1105 : mooseAssert(static_cast<bool>(grain._status & Status::DIRTY), "grain status is incorrect");
1106 :
1107 218 : grain_id_to_new_var.emplace_hint(
1108 : grain_id_to_new_var.end(),
1109 218 : std::pair<unsigned int, std::size_t>(grain._id, grain._var_index));
1110 :
1111 : /**
1112 : * Since the remapping algorithm only runs on the root process,
1113 : * the variable index on the primary's grains is inconsistent from
1114 : * the rest of the ranks. These are the grains with a status of
1115 : * DIRTY. As we build this map we will temporarily switch these
1116 : * variable indices back to the correct value so that all
1117 : * processors use the same algorithm to remap.
1118 : */
1119 218 : grain._var_index = old_var_index;
1120 : // Clear the DIRTY status as well for consistency
1121 : grain._status &= ~Status::DIRTY;
1122 : }
1123 : }
1124 :
1125 866 : if (!grain_id_to_new_var.empty())
1126 : {
1127 117 : if (_verbosity_level > 1)
1128 : {
1129 0 : _console << "Final remapping tally:\n";
1130 0 : for (const auto & remap_pair : grain_id_to_new_var)
1131 0 : _console << "Grain #" << remap_pair.first << " var_index "
1132 0 : << grain_id_to_existing_var_index[remap_pair.first] << " -> "
1133 0 : << remap_pair.second << '\n';
1134 0 : _console << "Communicating swaps with remaining processors..." << std::endl;
1135 : }
1136 : }
1137 : } // root processor
1138 :
1139 : // Communicate the std::map to all ranks
1140 2090 : _communicator.broadcast(grain_id_to_new_var);
1141 :
1142 : // Perform swaps if any occurred
1143 2090 : if (!grain_id_to_new_var.empty())
1144 : {
1145 : // Cache for holding values during swaps
1146 394 : std::vector<std::map<Node *, CacheValues>> cache(_n_vars);
1147 :
1148 : // Perform the actual swaps on all processors
1149 3267 : for (auto & grain : _feature_sets)
1150 : {
1151 : // See if this grain was remapped
1152 2873 : auto new_var_it = grain_id_to_new_var.find(grain._id);
1153 2873 : if (new_var_it != grain_id_to_new_var.end())
1154 457 : swapSolutionValues(grain, new_var_it->second, cache, RemapCacheMode::FILL);
1155 : }
1156 :
1157 3267 : for (auto & grain : _feature_sets)
1158 : {
1159 : // See if this grain was remapped
1160 2873 : auto new_var_it = grain_id_to_new_var.find(grain._id);
1161 2873 : if (new_var_it != grain_id_to_new_var.end())
1162 457 : swapSolutionValues(grain, new_var_it->second, cache, RemapCacheMode::USE);
1163 : }
1164 :
1165 394 : _sys.solution().close();
1166 394 : _sys.solutionOld().close();
1167 394 : _sys.solutionOlder().close();
1168 :
1169 394 : _sys.system().update();
1170 :
1171 394 : if (_verbosity_level > 1)
1172 0 : _console << "Swaps complete" << std::endl;
1173 394 : }
1174 2090 : }
1175 :
1176 : void
1177 798 : GrainTracker::computeMinDistancesFromGrain(FeatureData & grain,
1178 : std::vector<std::list<GrainDistance>> & min_distances)
1179 : {
1180 : /**
1181 : * In the diagram below assume we have 4 order parameters. The grain with the asterisk needs to
1182 : * be remapped. All order parameters are used in neighboring grains. For all "touching" grains,
1183 : * the value of the corresponding entry in min_distances will be a negative integer representing
1184 : * the number of immediate neighbors with that order parameter.
1185 : *
1186 : * Note: Only the first member of the pair (the distance) is shown in the array below.
1187 : * e.g. [-2.0, -max, -1.0, -2.0]
1188 : *
1189 : * After sorting, variable index 2 (value: -1.0) be at the end of the array and will be the first
1190 : * variable we attempt to renumber the current grain to.
1191 : *
1192 : * __ ___
1193 : * \ 0 / \
1194 : * 2 \___/ 1 \___
1195 : * / \ / \
1196 : * __/ 1 \___/ 2 \
1197 : * \ * / \ /
1198 : * 3 \___/ 3 \___/
1199 : * / \ /
1200 : * __/ 0 \___/
1201 : *
1202 : */
1203 12626 : for (const auto i : index_range(_feature_sets))
1204 : {
1205 11828 : auto & other_grain = _feature_sets[i];
1206 :
1207 11828 : if (other_grain._var_index == grain._var_index || other_grain._var_index >= _reserve_op_index)
1208 2012 : continue;
1209 :
1210 : auto target_var_index = other_grain._var_index;
1211 : auto target_grain_index = i;
1212 9816 : auto target_grain_id = other_grain._id;
1213 :
1214 9816 : Real curr_bbox_diff = boundingRegionDistance(grain._bboxes, other_grain._bboxes);
1215 :
1216 : GrainDistance grain_distance_obj(
1217 9816 : curr_bbox_diff, target_var_index, target_grain_index, target_grain_id);
1218 :
1219 : // To handle touching halos we penalize the top pick each time we see another
1220 9816 : if (curr_bbox_diff == -1.0 && !min_distances[target_var_index].empty())
1221 : {
1222 3762 : Real last_distance = min_distances[target_var_index].begin()->_distance;
1223 3762 : if (last_distance < 0)
1224 3095 : grain_distance_obj._distance += last_distance;
1225 : }
1226 :
1227 : // Insertion sort into a list
1228 : auto insert_it = min_distances[target_var_index].begin();
1229 10571 : while (insert_it != min_distances[target_var_index].end() && !(grain_distance_obj < *insert_it))
1230 : ++insert_it;
1231 9816 : min_distances[target_var_index].insert(insert_it, grain_distance_obj);
1232 : }
1233 :
1234 : /**
1235 : * See if we have any completely open OPs (excluding reserve order parameters) or the order
1236 : * parameter corresponding to this grain, we need to put them in the list or the grain tracker
1237 : * won't realize that those vars are available for remapping.
1238 : */
1239 7112 : for (const auto var_index : make_range(_reserve_op_index))
1240 : {
1241 : // Don't put an entry in for matching variable indices (i.e. we can't remap to ourselves)
1242 6314 : if (grain._var_index == var_index)
1243 788 : continue;
1244 :
1245 5526 : if (min_distances[var_index].empty())
1246 155 : min_distances[var_index].emplace_front(std::numeric_limits<Real>::max(), var_index);
1247 : }
1248 798 : }
1249 :
1250 : bool
1251 1091 : GrainTracker::attemptGrainRenumber(FeatureData & grain, unsigned int depth, unsigned int max_depth)
1252 : {
1253 : // End the recursion of our breadth first search
1254 1091 : if (depth > max_depth)
1255 : return false;
1256 :
1257 798 : std::size_t curr_var_index = grain._var_index;
1258 :
1259 : std::vector<std::map<Node *, CacheValues>> cache;
1260 :
1261 798 : std::vector<std::list<GrainDistance>> min_distances(_vars.size());
1262 :
1263 : /**
1264 : * We have two grains that are getting close represented by the same order parameter.
1265 : * We need to map to the variable whose closest grain to this one is furthest away by bounding
1266 : * region to bounding region distance.
1267 : */
1268 798 : computeMinDistancesFromGrain(grain, min_distances);
1269 :
1270 : /**
1271 : * We have a vector of the distances to the closest grains represented by each of our variables.
1272 : * We just need to pick a suitable grain to replace with. We will start with the maximum of this
1273 : * this list: (max of the mins), but will settle for next to largest and so forth as we make more
1274 : * attempts at remapping grains. This is a graph coloring problem so more work will be required
1275 : * to optimize this process.
1276 : *
1277 : * Note: We don't have an explicit check here to avoid remapping a variable to itself. This is
1278 : * unnecessary since the min_distance of a variable is explicitly set up above.
1279 : */
1280 : // clang-format off
1281 798 : std::sort(min_distances.begin(), min_distances.end(),
1282 : [](const std::list<GrainDistance> & lhs, const std::list<GrainDistance> & rhs)
1283 : {
1284 : // Sort lists in reverse order (largest distance first)
1285 : // These empty cases are here to make this comparison stable
1286 15408 : if (lhs.empty())
1287 : return false;
1288 13912 : else if (rhs.empty())
1289 : return true;
1290 : else
1291 11595 : return lhs.begin()->_distance > rhs.begin()->_distance;
1292 : });
1293 : // clang-format on
1294 :
1295 1469 : for (auto & list_ref : min_distances)
1296 : {
1297 : const auto target_it = list_ref.begin();
1298 1469 : if (target_it == list_ref.end())
1299 : continue;
1300 :
1301 : // If the distance is positive we can just remap and be done
1302 1469 : if (target_it->_distance > 0)
1303 : {
1304 117 : if (_verbosity_level > 0)
1305 : {
1306 117 : _console << COLOR_GREEN << "- Depth " << depth << ": Remapping grain #" << grain._id
1307 117 : << " from variable index " << curr_var_index << " to " << target_it->_var_index;
1308 117 : if (target_it->_distance == std::numeric_limits<Real>::max())
1309 82 : _console << " which currently contains zero grains.\n\n" << COLOR_DEFAULT;
1310 : else
1311 35 : _console << " whose closest grain (#" << target_it->_grain_id << ") is at a distance of "
1312 35 : << std::sqrt(target_it->_distance) << "\n\n"
1313 35 : << COLOR_DEFAULT;
1314 : }
1315 :
1316 : grain._status |= Status::DIRTY;
1317 117 : grain._var_index = target_it->_var_index;
1318 117 : return true;
1319 : }
1320 :
1321 : // If the distance isn't positive we just need to make sure that none of the grains represented
1322 : // by the target variable index would intersect this one if we were to remap
1323 : {
1324 : auto next_target_it = target_it;
1325 : bool intersection_hit = false;
1326 : unsigned short num_close_targets = 0;
1327 1352 : std::ostringstream oss;
1328 2740 : while (!intersection_hit && next_target_it != list_ref.end())
1329 : {
1330 1388 : if (next_target_it->_distance > 0)
1331 : break;
1332 :
1333 : mooseAssert(next_target_it->_grain_index < _feature_sets.size(),
1334 : "Error in indexing target grain in attemptGrainRenumber");
1335 1388 : FeatureData & next_target_grain = _feature_sets[next_target_it->_grain_index];
1336 :
1337 : // If any grains touch we're done here
1338 1388 : if (grain.halosIntersect(next_target_grain))
1339 : intersection_hit = true;
1340 : else
1341 : {
1342 94 : if (num_close_targets > 0)
1343 8 : oss << ", "; // delimiter
1344 94 : oss << "#" << next_target_it->_grain_id;
1345 : }
1346 :
1347 : ++next_target_it;
1348 1388 : ++num_close_targets;
1349 : }
1350 :
1351 1352 : if (!intersection_hit)
1352 : {
1353 58 : if (_verbosity_level > 0)
1354 : {
1355 58 : _console << COLOR_GREEN << "- Depth " << depth << ": Remapping grain #" << grain._id
1356 58 : << " from variable index " << curr_var_index << " to " << target_it->_var_index;
1357 :
1358 58 : if (num_close_targets == 1)
1359 150 : _console << " whose closest grain (" << oss.str()
1360 50 : << ") is inside our bounding box but whose halo is not touching.\n\n"
1361 50 : << COLOR_DEFAULT;
1362 : else
1363 24 : _console << " whose closest grains (" << oss.str()
1364 8 : << ") are inside our bounding box but whose halos are not touching.\n\n"
1365 8 : << COLOR_DEFAULT;
1366 : }
1367 :
1368 : grain._status |= Status::DIRTY;
1369 58 : grain._var_index = target_it->_var_index;
1370 : return true;
1371 : }
1372 1352 : }
1373 :
1374 : // If we reach this part of the loop, there is no simple renumbering that can be done.
1375 : mooseAssert(target_it->_grain_index < _feature_sets.size(),
1376 : "Error in indexing target grain in attemptGrainRenumber");
1377 1294 : FeatureData & target_grain = _feature_sets[target_it->_grain_index];
1378 :
1379 : /**
1380 : * If we get to this case and the best distance is less than -1, we are in big trouble.
1381 : * This means that grains represented by all of the remaining order parameters are
1382 : * overlapping this one in at least two places. We'd have to maintain multiple recursive
1383 : * chains, or just start over from scratch...
1384 : * Let's just return false and see if there is another remapping option.
1385 : */
1386 1294 : if (target_it->_distance < -1)
1387 : return false;
1388 :
1389 : // Make sure this grain isn't marked. If it is, we can't recurse here
1390 1011 : if ((target_grain._status & Status::MARKED) == Status::MARKED)
1391 : return false;
1392 :
1393 : /**
1394 : * Propose a new variable index for the current grain and recurse.
1395 : * We don't need to mark the status as DIRTY here since the recursion
1396 : * may fail. For now, we'll just add MARKED to the status.
1397 : */
1398 714 : grain._var_index = target_it->_var_index;
1399 : grain._status |= Status::MARKED;
1400 714 : if (attemptGrainRenumber(target_grain, depth + 1, max_depth))
1401 : {
1402 : // SUCCESS!
1403 43 : if (_verbosity_level > 0)
1404 43 : _console << COLOR_GREEN << "- Depth " << depth << ": Remapping grain #" << grain._id
1405 43 : << " from variable index " << curr_var_index << " to " << target_it->_var_index
1406 43 : << "\n\n"
1407 43 : << COLOR_DEFAULT;
1408 :
1409 : // Now we need to mark the grain as DIRTY since the recursion succeeded.
1410 : grain._status |= Status::DIRTY;
1411 43 : return true;
1412 : }
1413 : else
1414 : // FAILURE, We need to set our var index back after failed recursive step
1415 671 : grain._var_index = curr_var_index;
1416 :
1417 : // ALWAYS "unmark" (or clear the MARKED status) after recursion so it can be used by other remap
1418 : // operations
1419 : grain._status &= ~Status::MARKED;
1420 : }
1421 :
1422 : return false;
1423 798 : }
1424 :
1425 : void
1426 914 : GrainTracker::swapSolutionValues(FeatureData & grain,
1427 : std::size_t new_var_index,
1428 : std::vector<std::map<Node *, CacheValues>> & cache,
1429 : RemapCacheMode cache_mode)
1430 : {
1431 914 : MeshBase & mesh = _mesh.getMesh();
1432 :
1433 : // Remap the grain
1434 : std::set<Node *> updated_nodes_tmp; // Used only in the elemental case
1435 45612 : for (auto entity : grain._local_ids)
1436 : {
1437 44698 : if (_is_elemental)
1438 : {
1439 44698 : Elem * elem = mesh.query_elem_ptr(entity);
1440 44698 : if (!elem)
1441 0 : continue;
1442 :
1443 223490 : for (unsigned int i = 0; i < elem->n_nodes(); ++i)
1444 : {
1445 178792 : Node * curr_node = elem->node_ptr(i);
1446 178792 : if (updated_nodes_tmp.find(curr_node) == updated_nodes_tmp.end())
1447 : {
1448 : // cache this node so we don't attempt to remap it again within this loop
1449 58028 : updated_nodes_tmp.insert(curr_node);
1450 58028 : swapSolutionValuesHelper(curr_node, grain._var_index, new_var_index, cache, cache_mode);
1451 : }
1452 : }
1453 : }
1454 : else
1455 0 : swapSolutionValuesHelper(
1456 0 : mesh.query_node_ptr(entity), grain._var_index, new_var_index, cache, cache_mode);
1457 : }
1458 :
1459 : // Update the variable index in the unique grain datastructure after swaps are complete
1460 914 : if (cache_mode == RemapCacheMode::USE || cache_mode == RemapCacheMode::BYPASS)
1461 457 : grain._var_index = new_var_index;
1462 914 : }
1463 :
1464 : void
1465 58028 : GrainTracker::swapSolutionValuesHelper(Node * curr_node,
1466 : std::size_t curr_var_index,
1467 : std::size_t new_var_index,
1468 : std::vector<std::map<Node *, CacheValues>> & cache,
1469 : RemapCacheMode cache_mode)
1470 : {
1471 58028 : if (curr_node && curr_node->processor_id() == processor_id())
1472 : {
1473 : // Reinit the node so we can get and set values of the solution here
1474 49470 : _subproblem.reinitNode(curr_node, 0);
1475 :
1476 : // Local variables to hold values being transferred
1477 : Real current, old = 0, older = 0;
1478 : // Retrieve the value either from the old variable or cache
1479 49470 : if (cache_mode == RemapCacheMode::FILL || cache_mode == RemapCacheMode::BYPASS)
1480 : {
1481 24735 : current = _vars[curr_var_index]->dofValues()[0];
1482 24735 : if (_is_transient)
1483 : {
1484 24735 : old = _vars[curr_var_index]->dofValuesOld()[0];
1485 24735 : older = _vars[curr_var_index]->dofValuesOlder()[0];
1486 : }
1487 : }
1488 : else // USE
1489 : {
1490 : const auto cache_it = cache[curr_var_index].find(curr_node);
1491 : mooseAssert(cache_it != cache[curr_var_index].end(), "Error in cache");
1492 24735 : current = cache_it->second.current;
1493 24735 : old = cache_it->second.old;
1494 24735 : older = cache_it->second.older;
1495 : }
1496 :
1497 : // Cache the value or use it!
1498 49470 : if (cache_mode == RemapCacheMode::FILL)
1499 : {
1500 24735 : cache[curr_var_index][curr_node].current = current;
1501 24735 : cache[curr_var_index][curr_node].old = old;
1502 24735 : cache[curr_var_index][curr_node].older = older;
1503 : }
1504 : else // USE or BYPASS
1505 : {
1506 24735 : const auto & dof_index = _vars[new_var_index]->nodalDofIndex();
1507 :
1508 : // Transfer this solution from the old to the current
1509 24735 : _sys.solution().set(dof_index, current);
1510 24735 : if (_is_transient)
1511 : {
1512 24735 : _sys.solutionOld().set(dof_index, old);
1513 24735 : _sys.solutionOlder().set(dof_index, older);
1514 : }
1515 : }
1516 :
1517 : /**
1518 : * Finally zero out the old variable. When using the FILL/USE combination to
1519 : * read/write variables, it's important to zero the variable on the FILL
1520 : * stage and not the USE stage. The reason for this is handling swaps as
1521 : * illustrated in the following diagram
1522 : * ___ ___
1523 : * / \/ \ If adjacent grains (overlapping flood region) end up
1524 : * / 1 /\ 2 \ swapping variable indices and variables are zeroed on
1525 : * \ 2*\/ 1* / "USE", the overlap region will be incorrectly zeroed
1526 : * \___/\___/ by whichever variable is written to second.
1527 : *.
1528 : */
1529 49470 : if (cache_mode == RemapCacheMode::FILL || cache_mode == RemapCacheMode::BYPASS)
1530 : {
1531 24735 : const auto & dof_index = _vars[curr_var_index]->nodalDofIndex();
1532 :
1533 : // Set the DOF for the current variable to zero
1534 24735 : _sys.solution().set(dof_index, -_bound_value);
1535 24735 : if (_is_transient)
1536 : {
1537 24735 : _sys.solutionOld().set(dof_index, -_bound_value);
1538 24735 : _sys.solutionOlder().set(dof_index, -_bound_value);
1539 : }
1540 : }
1541 : }
1542 58028 : }
1543 :
1544 : void
1545 2098 : GrainTracker::updateFieldInfo()
1546 : {
1547 4196 : TIME_SECTION("updateFieldInfo", 3, "Updating Field Info");
1548 :
1549 16426 : for (const auto map_num : make_range(_maps_size))
1550 : _feature_maps[map_num].clear();
1551 :
1552 : std::map<dof_id_type, Real> tmp_map;
1553 :
1554 16568 : for (const auto & grain : _feature_sets)
1555 : {
1556 14470 : std::size_t curr_var = grain._var_index;
1557 14470 : std::size_t map_index = (_single_map_mode || _condense_map_info) ? 0 : curr_var;
1558 :
1559 945402 : for (auto entity : grain._local_ids)
1560 : {
1561 : // Highest variable value at this entity wins
1562 : Real entity_value = std::numeric_limits<Real>::lowest();
1563 930932 : if (_is_elemental)
1564 : {
1565 920680 : const Elem * elem = _mesh.elemPtr(entity);
1566 920680 : std::vector<Point> centroid(1, elem->vertex_average());
1567 920680 : if (_poly_ic_uo && _first_time)
1568 : {
1569 157353 : entity_value = _poly_ic_uo->getVariableValue(grain._var_index, centroid[0]);
1570 : }
1571 : else
1572 : {
1573 763327 : _fe_problem.reinitElemPhys(elem, centroid, 0);
1574 763327 : entity_value = _vars[curr_var]->sln()[0];
1575 : }
1576 920680 : }
1577 : else
1578 : {
1579 10252 : auto node_ptr = _mesh.nodePtr(entity);
1580 10252 : entity_value = _vars[curr_var]->getNodalValue(*node_ptr);
1581 : }
1582 :
1583 930932 : if (entity_value != std::numeric_limits<Real>::lowest() &&
1584 342254 : (tmp_map.find(entity) == tmp_map.end() || entity_value > tmp_map[entity]))
1585 : {
1586 : mooseAssert(grain._id != invalid_id, "Missing Grain ID");
1587 735105 : _feature_maps[map_index][entity] = grain._id;
1588 :
1589 735105 : if (_var_index_mode)
1590 735105 : _var_index_maps[map_index][entity] = grain._var_index;
1591 :
1592 735105 : tmp_map[entity] = entity_value;
1593 : }
1594 :
1595 930932 : if (_compute_var_to_feature_map)
1596 : {
1597 177651 : auto insert_pair = moose_try_emplace(
1598 177651 : _entity_var_to_features, entity, std::vector<unsigned int>(_n_vars, invalid_id));
1599 : auto & vec_ref = insert_pair.first->second;
1600 :
1601 177651 : if (insert_pair.second)
1602 : {
1603 : // insert the reserve op numbers (if appropriate)
1604 128458 : for (const auto reserve_index : make_range(_n_reserve_ops))
1605 0 : vec_ref[reserve_index] = _reserve_grain_first_index + reserve_index;
1606 : }
1607 177651 : vec_ref[grain._var_index] = grain._id;
1608 : }
1609 : }
1610 :
1611 14470 : if (_compute_halo_maps)
1612 933823 : for (auto entity : grain._halo_ids)
1613 921850 : _halo_ids[grain._var_index][entity] = grain._var_index;
1614 :
1615 413404 : for (auto entity : grain._ghosted_ids)
1616 398934 : _ghosted_entity_ids[entity] = 1;
1617 : }
1618 :
1619 2098 : communicateHaloMap();
1620 2098 : }
1621 :
1622 : void
1623 2098 : GrainTracker::communicateHaloMap()
1624 : {
1625 2098 : if (_compute_halo_maps)
1626 : {
1627 : // rank var_index entity_id
1628 : std::vector<std::pair<std::size_t, dof_id_type>> halo_ids_all;
1629 :
1630 : std::vector<int> counts;
1631 : std::vector<std::pair<std::size_t, dof_id_type>> local_halo_ids;
1632 : std::size_t counter = 0;
1633 :
1634 1639 : const bool isDistributedMesh = _mesh.isDistributedMesh();
1635 :
1636 1639 : if (_is_primary)
1637 : {
1638 526 : std::vector<std::vector<std::pair<std::size_t, dof_id_type>>> root_halo_ids(_n_procs);
1639 526 : counts.resize(_n_procs);
1640 :
1641 : // Loop over the _halo_ids "field" and build minimal lists for all of the other ranks
1642 4346 : for (const auto var_index : index_range(_halo_ids))
1643 : {
1644 322041 : for (const auto & entity_pair : _halo_ids[var_index])
1645 : {
1646 318221 : auto entity_id = entity_pair.first;
1647 318221 : if (isDistributedMesh)
1648 : {
1649 : // Check to see which contiguous range this entity ID falls into
1650 : auto range_it =
1651 25325 : std::lower_bound(_all_ranges.begin(),
1652 : _all_ranges.end(),
1653 : entity_id,
1654 : [](const std::pair<dof_id_type, dof_id_type> range,
1655 : dof_id_type entity_id) { return range.second < entity_id; });
1656 :
1657 : mooseAssert(range_it != _all_ranges.end(), "No range round?");
1658 :
1659 : // Recover the index from the iterator
1660 : auto proc_id = std::distance(_all_ranges.begin(), range_it);
1661 :
1662 : // Now add this halo entity to the map for the corresponding proc to scatter latter
1663 25325 : root_halo_ids[proc_id].push_back(std::make_pair(var_index, entity_id));
1664 : }
1665 : else
1666 : {
1667 : DofObject * halo_entity;
1668 292896 : if (_is_elemental)
1669 292896 : halo_entity = _mesh.queryElemPtr(entity_id);
1670 : else
1671 0 : halo_entity = _mesh.queryNodePtr(entity_id);
1672 :
1673 : if (halo_entity)
1674 292896 : root_halo_ids[halo_entity->processor_id()].push_back(
1675 292896 : std::make_pair(var_index, entity_id));
1676 : }
1677 : }
1678 : }
1679 :
1680 : // Build up the counts vector for MPI scatter
1681 2165 : for (const auto & vector_ref : root_halo_ids)
1682 : {
1683 : std::copy(vector_ref.begin(), vector_ref.end(), std::back_inserter(halo_ids_all));
1684 1639 : counts[counter] = vector_ref.size();
1685 1639 : counter++;
1686 : }
1687 526 : }
1688 :
1689 1639 : _communicator.scatter(halo_ids_all, counts, local_halo_ids);
1690 :
1691 : // Now add the contributions from the root process to the processor local maps
1692 319860 : for (const auto & halo_pair : local_halo_ids)
1693 318221 : _halo_ids[halo_pair.first].emplace(std::make_pair(halo_pair.second, halo_pair.first));
1694 :
1695 : /**
1696 : * Finally remove halo markings from interior regions. This step is necessary because we expand
1697 : * halos _before_ we do communication but that expansion can and will likely go into the
1698 : * interior of the grain (from a single processor's perspective). We could expand halos after
1699 : * merging, but that would likely be less scalable.
1700 : */
1701 13612 : for (const auto & grain : _feature_sets)
1702 645152 : for (auto local_id : grain._local_ids)
1703 633179 : _halo_ids[grain._var_index].erase(local_id);
1704 1639 : }
1705 2098 : }
1706 :
1707 : Real
1708 8460 : GrainTracker::centroidRegionDistance(std::vector<BoundingBox> & bboxes1,
1709 : std::vector<BoundingBox> & bboxes2) const
1710 : {
1711 : /**
1712 : * Find the minimum centroid distance between any to pieces of the grains.
1713 : */
1714 : auto min_distance = std::numeric_limits<Real>::max();
1715 17507 : for (const auto & bbox1 : bboxes1)
1716 : {
1717 : const auto centroid_point1 = (bbox1.max() + bbox1.min()) / 2.0;
1718 :
1719 19559 : for (const auto & bbox2 : bboxes2)
1720 : {
1721 : const auto centroid_point2 = (bbox2.max() + bbox2.min()) / 2.0;
1722 :
1723 : // Here we'll calculate a distance between the centroids
1724 : auto curr_distance =
1725 10512 : _mesh.minPeriodicDistance(*_fe_vars[0], centroid_point1, centroid_point2);
1726 :
1727 10512 : if (curr_distance < min_distance)
1728 : min_distance = curr_distance;
1729 : }
1730 : }
1731 :
1732 8460 : return min_distance;
1733 : }
1734 :
1735 : Real
1736 9816 : GrainTracker::boundingRegionDistance(std::vector<BoundingBox> & bboxes1,
1737 : std::vector<BoundingBox> & bboxes2) const
1738 : {
1739 : /**
1740 : * The region that each grain covers is represented by a bounding box large enough to encompassing
1741 : * all the points within that grain. When using periodic boundaries, we may have several discrete
1742 : * "pieces" of a grain each represented by a bounding box. The distance between any two grains
1743 : * is defined as the minimum distance between any pair of boxes, one selected from each grain.
1744 : */
1745 : auto min_distance = std::numeric_limits<Real>::max();
1746 11229 : for (const auto & bbox1 : bboxes1)
1747 : {
1748 11229 : for (const auto & bbox2 : bboxes2)
1749 : {
1750 : // AABB squared distance
1751 : Real curr_distance = 0.0;
1752 : bool boxes_overlap = true;
1753 39264 : for (unsigned int dim = 0; dim < LIBMESH_DIM; ++dim)
1754 : {
1755 : const auto & min1 = bbox1.min()(dim);
1756 : const auto & max1 = bbox1.max()(dim);
1757 : const auto & min2 = bbox2.min()(dim);
1758 : const auto & max2 = bbox2.max()(dim);
1759 :
1760 29448 : if (min1 > max2)
1761 : {
1762 955 : const auto delta = max2 - min1;
1763 955 : curr_distance += delta * delta;
1764 : boxes_overlap = false;
1765 : }
1766 28493 : else if (min2 > max1)
1767 : {
1768 470 : const auto delta = max1 - min2;
1769 470 : curr_distance += delta * delta;
1770 : boxes_overlap = false;
1771 : }
1772 : }
1773 :
1774 9816 : if (boxes_overlap)
1775 : return -1.0; /* all overlaps are treated the same */
1776 :
1777 1413 : if (curr_distance < min_distance)
1778 : min_distance = curr_distance;
1779 : }
1780 : }
1781 :
1782 : return min_distance;
1783 : }
1784 :
1785 : unsigned int
1786 10 : GrainTracker::getNextUniqueID()
1787 : {
1788 : /**
1789 : * Get the next unique grain ID but make sure to respect
1790 : * reserve ids. Note, that the first valid ID for a new
1791 : * grain is _reserve_grain_first_index + _n_reserve_ops because
1792 : * _reserve_grain_first_index IS a valid index. It does not
1793 : * point to the last valid index of the non-reserved grains.
1794 : */
1795 10 : _max_curr_grain_id = std::max(_max_curr_grain_id == invalid_id ? 0 : _max_curr_grain_id + 1,
1796 10 : _reserve_grain_first_index + _n_reserve_ops /* no +1 here!*/);
1797 :
1798 10 : return _max_curr_grain_id;
1799 : }
1800 :
1801 : /*************************************************
1802 : ************** Helper Structures ****************
1803 : ************************************************/
1804 155 : GrainDistance::GrainDistance(Real distance, std::size_t var_index)
1805 : : GrainDistance(distance,
1806 : var_index,
1807 : std::numeric_limits<std::size_t>::max(),
1808 155 : std::numeric_limits<unsigned int>::max())
1809 : {
1810 155 : }
1811 :
1812 9971 : GrainDistance::GrainDistance(Real distance,
1813 : std::size_t var_index,
1814 : std::size_t grain_index,
1815 9971 : unsigned int grain_id)
1816 9971 : : _distance(distance), _var_index(var_index), _grain_index(grain_index), _grain_id(grain_id)
1817 : {
1818 9971 : }
1819 :
1820 : bool
1821 4520 : GrainDistance::operator<(const GrainDistance & rhs) const
1822 : {
1823 4520 : return _distance < rhs._distance;
1824 : }
|
