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 9746 : dataStore(std::ostream & stream, GrainTracker::PartialFeatureData & feature, void * context)
29 : {
30 9746 : storeHelper(stream, feature.boundary_intersection, context);
31 9746 : storeHelper(stream, feature.id, context);
32 9746 : storeHelper(stream, feature.centroid, context);
33 9746 : storeHelper(stream, feature.status, context);
34 9746 : }
35 :
36 : template <>
37 : void
38 16781 : dataLoad(std::istream & stream, GrainTracker::PartialFeatureData & feature, void * context)
39 : {
40 16781 : loadHelper(stream, feature.boundary_intersection, context);
41 16781 : loadHelper(stream, feature.id, context);
42 16781 : loadHelper(stream, feature.centroid, context);
43 16781 : loadHelper(stream, feature.status, context);
44 16781 : }
45 :
46 : registerMooseObject("PhaseFieldApp", GrainTracker);
47 :
48 : InputParameters
49 1462 : GrainTracker::validParams()
50 : {
51 1462 : InputParameters params = FeatureFloodCount::validParams();
52 1462 : params += GrainTrackerInterface::validParams();
53 :
54 : // FeatureFloodCount adds a relationship manager, but we need to extend that for GrainTracker
55 : params.clearRelationshipManagers();
56 :
57 2924 : params.addRelationshipManager(
58 : "ElementSideNeighborLayers",
59 : Moose::RelationshipManagerType::GEOMETRIC,
60 :
61 2193 : [](const InputParameters & obj_params, InputParameters & rm_params)
62 2193 : { rm_params.set<unsigned short>("layers") = obj_params.get<unsigned short>("halo_level"); }
63 :
64 : );
65 :
66 2924 : params.addRelationshipManager("ElementSideNeighborLayers",
67 : Moose::RelationshipManagerType::ALGEBRAIC);
68 :
69 : // The GrainTracker requires non-volatile storage for tracking grains across invocations.
70 1462 : params.set<bool>("restartable_required") = true;
71 :
72 2924 : params.addParam<Real>("bound_value",
73 2924 : 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 1462 : params.addClassDescription("Grain Tracker object for running reduced order parameter simulations "
79 : "without grain coalescence.");
80 :
81 1462 : return params;
82 0 : }
83 :
84 731 : GrainTracker::GrainTracker(const InputParameters & parameters)
85 : : FeatureFloodCount(parameters),
86 : GrainTrackerInterface(),
87 731 : _tracking_step(getParam<int>("tracking_step")),
88 1462 : _halo_level(getParam<unsigned short>("halo_level")),
89 1462 : _max_remap_recursion_depth(getParam<unsigned short>("max_remap_recursion_depth")),
90 1462 : _n_reserve_ops(getParam<unsigned short>("reserve_op")),
91 731 : _reserve_op_index(_n_reserve_ops <= _n_vars ? _n_vars - _n_reserve_ops : 0),
92 1462 : _reserve_op_threshold(getParam<Real>("reserve_op_threshold")),
93 1462 : _bound_value(getParam<Real>("bound_value")),
94 1462 : _remap(getParam<bool>("remap_grains")),
95 1462 : _tolerate_failure(getParam<bool>("tolerate_failure")),
96 731 : _poly_ic_uo(parameters.isParamValid("polycrystal_ic_uo")
97 1130 : ? &getUserObject<PolycrystalUserObjectBase>("polycrystal_ic_uo")
98 : : nullptr),
99 1462 : _verbosity_level(getParam<short>("verbosity_level")),
100 1462 : _first_time(declareRestartableData<bool>("first_time", true)),
101 1462 : _error_on_grain_creation(getParam<bool>("error_on_grain_creation")),
102 731 : _reserve_grain_first_index(0),
103 731 : _old_max_grain_id(0),
104 1462 : _max_curr_grain_id(declareRestartableData<unsigned int>("max_curr_grain_id", invalid_id)),
105 1462 : _is_transient(_subproblem.isTransient())
106 : {
107 731 : 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 731 : 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 731 : }
115 :
116 1454 : GrainTracker::~GrainTracker() {}
117 :
118 : Real
119 6032498 : GrainTracker::getEntityValue(dof_id_type entity_id,
120 : FieldType field_type,
121 : std::size_t var_index) const
122 : {
123 6032498 : if (_t_step < _tracking_step)
124 : return 0;
125 :
126 6032498 : return FeatureFloodCount::getEntityValue(entity_id, field_type, var_index);
127 : }
128 :
129 : const std::vector<unsigned int> &
130 175102504 : GrainTracker::getVarToFeatureVector(dof_id_type elem_id) const
131 : {
132 175102504 : return FeatureFloodCount::getVarToFeatureVector(elem_id);
133 : }
134 :
135 : unsigned int
136 476 : GrainTracker::getFeatureVar(unsigned int feature_id) const
137 : {
138 476 : 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 1093 : GrainTracker::getTotalFeatureCount() const
150 : {
151 : // Note: This value is parallel consistent, see assignGrains()/trackGrains()
152 1093 : return _max_curr_grain_id == invalid_id ? 0 : _max_curr_grain_id + 1;
153 : }
154 :
155 : Point
156 14723084 : GrainTracker::getGrainCentroid(unsigned int grain_id) const
157 : {
158 : mooseAssert(grain_id < _feature_id_to_local_index.size(), "Grain ID out of bounds");
159 14723084 : auto grain_index = _feature_id_to_local_index[grain_id];
160 :
161 14723084 : 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 14723084 : return _feature_sets[_feature_id_to_local_index[grain_id]]._centroid;
167 : }
168 :
169 : // Inactive grain
170 : return Point();
171 : }
172 :
173 : bool
174 476 : 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 476 : auto feature_index = _feature_id_to_local_index[feature_id];
180 476 : if (feature_index != invalid_size_t)
181 : {
182 : mooseAssert(feature_index < _feature_sets.size(), "Grain index out of bounds");
183 398 : return _feature_sets[feature_index]._boundary_intersection != BoundaryIntersection::NONE;
184 : }
185 :
186 : return false;
187 : }
188 :
189 : bool
190 476 : 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 476 : auto feature_index = _feature_id_to_local_index[feature_id];
196 476 : if (feature_index != invalid_size_t)
197 : {
198 : mooseAssert(feature_index < _feature_sets.size(), "Grain index out of bounds");
199 398 : return ((_feature_sets[feature_index]._boundary_intersection &
200 398 : BoundaryIntersection::SPECIFIED_BOUNDARY) == BoundaryIntersection::SPECIFIED_BOUNDARY);
201 : }
202 :
203 : return false;
204 : }
205 :
206 : bool
207 476 : 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 476 : auto feature_index = _feature_id_to_local_index[feature_id];
213 476 : if (feature_index != invalid_size_t)
214 : {
215 : mooseAssert(feature_index < _feature_sets.size(), "Grain index out of bounds");
216 398 : bool primary = ((_feature_sets[feature_index]._boundary_intersection &
217 : BoundaryIntersection::PRIMARY_PERCOLATION_BOUNDARY) ==
218 398 : BoundaryIntersection::PRIMARY_PERCOLATION_BOUNDARY);
219 : bool secondary = ((_feature_sets[feature_index]._boundary_intersection &
220 : BoundaryIntersection::SECONDARY_PERCOLATION_BOUNDARY) ==
221 398 : BoundaryIntersection::SECONDARY_PERCOLATION_BOUNDARY);
222 398 : return (primary && secondary);
223 : }
224 :
225 : return false;
226 : }
227 :
228 : void
229 2616 : GrainTracker::initialize()
230 : {
231 : // Don't track grains if the current simulation step is before the specified tracking step
232 2616 : 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 2616 : if (!_first_time)
241 2169 : _feature_sets_old.swap(_feature_sets);
242 :
243 2616 : FeatureFloodCount::initialize();
244 : }
245 :
246 : void
247 915 : GrainTracker::meshChanged()
248 : {
249 : // Update the element ID ranges for use when computing halo maps
250 915 : if (_compute_halo_maps && _mesh.isDistributedMesh())
251 : {
252 33 : _all_ranges.clear();
253 :
254 33 : auto range = std::make_pair(std::numeric_limits<dof_id_type>::max(),
255 : std::numeric_limits<dof_id_type>::min());
256 60482 : for (const auto & current_elem : _mesh.getMesh().active_local_element_ptr_range())
257 : {
258 30208 : auto id = current_elem->id();
259 30208 : if (id < range.first)
260 33 : range.first = id;
261 30175 : else if (id > range.second)
262 30175 : range.second = id;
263 33 : }
264 :
265 33 : _communicator.gather(0, range, _all_ranges);
266 : }
267 :
268 915 : FeatureFloodCount::meshChanged();
269 915 : }
270 :
271 : void
272 2616 : GrainTracker::execute()
273 : {
274 : // Don't track grains if the current simulation step is before the specified tracking step
275 2616 : if (_t_step < _tracking_step)
276 : return;
277 :
278 2616 : if (_poly_ic_uo && _first_time)
279 : return;
280 :
281 2379 : FeatureFloodCount::execute();
282 : }
283 :
284 : Real
285 4567797 : 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 4567797 : if (var_index >= _reserve_op_index)
291 142440 : return _reserve_op_threshold;
292 : else
293 4425357 : return _step_threshold;
294 : }
295 :
296 : void
297 237 : GrainTracker::prepopulateState(const FeatureFloodCount & ffc_object)
298 : {
299 : mooseAssert(_first_time, "This method should only be called on the first invocation");
300 :
301 237 : _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 237 : if (_is_primary)
309 : {
310 : const auto & features = ffc_object.getFeatures();
311 2122 : for (auto & feature : features)
312 2001 : _feature_sets.emplace_back(feature.duplicate());
313 :
314 121 : _feature_count = _feature_sets.size();
315 : }
316 : else
317 : {
318 : const auto & features = ffc_object.getFeatures();
319 : _partial_feature_sets[0].clear();
320 920 : for (auto & feature : features)
321 804 : _partial_feature_sets[0].emplace_back(feature.duplicate());
322 : }
323 :
324 : // Make sure that feature count is communicated to all ranks
325 237 : _communicator.broadcast(_feature_count);
326 237 : }
327 :
328 : void
329 2616 : GrainTracker::finalize()
330 : {
331 : // Don't track grains if the current simulation step is before the specified tracking step
332 2616 : if (_t_step < _tracking_step)
333 : return;
334 :
335 5232 : TIME_SECTION("finalize", 3, "Finalizing GrainTracker");
336 :
337 : // Expand the depth of the halos around all grains
338 2616 : auto num_halo_layers = _halo_level >= 1
339 2616 : ? _halo_level - 1
340 : : 0; // The first level of halos already exists so subtract one
341 :
342 2616 : if (_poly_ic_uo && _first_time)
343 237 : prepopulateState(*_poly_ic_uo);
344 : else
345 : {
346 2379 : expandEdgeHalos(num_halo_layers);
347 :
348 : // Build up the grain map on the root processor
349 2379 : communicateAndMerge();
350 : }
351 :
352 : /**
353 : * Assign or Track Grains
354 : */
355 2616 : if (_first_time)
356 447 : assignGrains();
357 : else
358 2169 : trackGrains();
359 :
360 2616 : if (_verbosity_level > 1)
361 15 : _console << "Finished inside of trackGrains" << std::endl;
362 :
363 : /**
364 : * Broadcast essential data
365 : */
366 2616 : broadcastAndUpdateGrainData();
367 :
368 : /**
369 : * Remap Grains
370 : */
371 2616 : if (_remap)
372 2607 : remapGrains();
373 :
374 2616 : updateFieldInfo();
375 2616 : 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 2616 : _first_time = false;
380 :
381 : // TODO: Release non essential memory
382 2616 : if (_verbosity_level > 0)
383 2616 : _console << "Finished inside of GrainTracker\n" << std::endl;
384 : }
385 :
386 : void
387 2616 : GrainTracker::broadcastAndUpdateGrainData()
388 : {
389 5232 : TIME_SECTION("broadcastAndUpdateGrainData", 3, "Broadcasting and Updating Grain Data");
390 :
391 : std::vector<PartialFeatureData> root_feature_data;
392 2616 : std::vector<std::string> send_buffer(1), recv_buffer;
393 :
394 2616 : if (_is_primary)
395 : {
396 1033 : root_feature_data.reserve(_feature_sets.size());
397 :
398 : // Populate a subset of the information in a small data structure
399 1033 : std::transform(_feature_sets.begin(),
400 1033 : _feature_sets.end(),
401 : std::back_inserter(root_feature_data),
402 : [](FeatureData & feature)
403 : {
404 : PartialFeatureData partial_feature;
405 9746 : partial_feature.boundary_intersection = feature._boundary_intersection;
406 9746 : partial_feature.id = feature._id;
407 9746 : partial_feature.centroid = feature._centroid;
408 9746 : partial_feature.status = feature._status;
409 : return partial_feature;
410 : });
411 :
412 1033 : std::ostringstream oss;
413 1033 : dataStore(oss, root_feature_data, this);
414 1033 : send_buffer[0].assign(oss.str());
415 1033 : }
416 :
417 : // Broadcast the data to all ranks
418 2616 : _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 2616 : if (!_is_primary)
426 : {
427 1583 : std::istringstream iss;
428 : iss.str(recv_buffer[0]);
429 1583 : iss.clear();
430 :
431 1583 : dataLoad(iss, root_feature_data, this);
432 :
433 18364 : for (const auto & partial_data : root_feature_data)
434 : {
435 : // See if this processor has a record of this grain
436 16781 : if (partial_data.id < _feature_id_to_local_index.size() &&
437 16781 : _feature_id_to_local_index[partial_data.id] != invalid_size_t)
438 : {
439 8113 : auto & grain = _feature_sets[_feature_id_to_local_index[partial_data.id]];
440 8113 : grain._boundary_intersection = partial_data.boundary_intersection;
441 8113 : grain._centroid = partial_data.centroid;
442 8113 : if (partial_data.status == Status::INACTIVE)
443 0 : grain._status = Status::INACTIVE;
444 : }
445 : }
446 1583 : }
447 5232 : }
448 :
449 : void
450 447 : 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 447 : if (_is_primary)
460 : {
461 : // Find the largest grain ID, this requires sorting if the ID is not already set
462 241 : sortAndLabel();
463 :
464 241 : if (_feature_sets.empty())
465 : {
466 6 : _max_curr_grain_id = invalid_id;
467 6 : _reserve_grain_first_index = 0;
468 : }
469 : else
470 : {
471 235 : _max_curr_grain_id = _feature_sets.back()._id;
472 235 : _reserve_grain_first_index = _max_curr_grain_id + 1;
473 : }
474 :
475 2920 : for (auto & grain : _feature_sets)
476 2679 : 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 447 : scatterAndUpdateRanks();
486 :
487 : // Build up an id to index map
488 447 : _communicator.broadcast(_max_curr_grain_id);
489 447 : buildFeatureIdToLocalIndices(_max_curr_grain_id);
490 :
491 : // Now trigger the newGrainCreated() callback on all ranks
492 447 : if (_max_curr_grain_id != invalid_id)
493 4997 : for (unsigned int new_id = 0; new_id <= _max_curr_grain_id; ++new_id)
494 4558 : newGrainCreated(new_id);
495 447 : }
496 :
497 : void
498 2169 : GrainTracker::trackGrains()
499 : {
500 4338 : 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 2169 : 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 2169 : if (_is_primary)
512 : {
513 : // Reset Status on active unique grains
514 792 : std::vector<unsigned int> map_sizes(_maps_size);
515 7915 : for (auto & grain : _feature_sets_old)
516 : {
517 7123 : if (grain._status != Status::INACTIVE)
518 : {
519 7123 : grain._status = Status::CLEAR;
520 7123 : map_sizes[grain._var_index]++;
521 : }
522 : }
523 :
524 : // Print out stats on overall tracking changes per var_index
525 792 : if (_verbosity_level > 0)
526 : {
527 792 : _console << "\nGrain Tracker Status:";
528 5761 : for (const auto map_num : make_range(_maps_size))
529 : {
530 4969 : _console << "\nGrains active index " << map_num << ": " << map_sizes[map_num] << " -> "
531 4969 : << _feature_counts_per_map[map_num];
532 4969 : if (map_sizes[map_num] > _feature_counts_per_map[map_num])
533 74 : _console << "--";
534 4895 : else if (map_sizes[map_num] < _feature_counts_per_map[map_num])
535 18 : _console << "++";
536 : }
537 792 : _console << '\n' << std::endl;
538 : }
539 :
540 : // Before we track grains, lets sort them so that we get parallel consistent answers
541 792 : 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 792 : std::vector<std::size_t> new_grain_index_to_existing_grain_index(_feature_sets.size(),
550 792 : invalid_size_t);
551 :
552 7915 : for (const auto old_grain_index : index_range(_feature_sets_old))
553 : {
554 : auto & old_grain = _feature_sets_old[old_grain_index];
555 :
556 7123 : 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 7123 : std::lower_bound(_feature_sets.begin(), _feature_sets.end(), old_grain._var_index,
570 : [](const FeatureData & item, std::size_t var_index)
571 : {
572 26893 : 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 20610 : for (MooseIndex(_feature_sets)
579 7123 : new_grain_index = std::distance(_feature_sets.begin(), start_it);
580 20610 : new_grain_index < _feature_sets.size() &&
581 19661 : _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 13487 : if (new_grain.boundingBoxesIntersect(old_grain))
591 : {
592 : any_boxes_intersect = true;
593 10074 : Real curr_centroid_diff = centroidRegionDistance(old_grain._bboxes, new_grain._bboxes);
594 10074 : 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 7123 : 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 7123 : 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 7057 : auto curr_index = new_grain_index_to_existing_grain_index[closest_match_index];
616 7057 : if (curr_index != invalid_size_t)
617 : {
618 : // The new feature being competed for
619 11 : 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 11 : auto centroid_diff1 = centroidRegionDistance(new_grain._bboxes, old_grain._bboxes);
625 11 : auto centroid_diff2 = centroidRegionDistance(new_grain._bboxes, other_old_grain._bboxes);
626 :
627 11 : auto & inactive_grain = (centroid_diff1 < centroid_diff2) ? other_old_grain : old_grain;
628 :
629 11 : inactive_grain._status = Status::INACTIVE;
630 11 : if (_verbosity_level > 0)
631 : {
632 13 : _console << COLOR_GREEN << "Marking Grain " << inactive_grain._id
633 11 : << " as INACTIVE (variable index: " << inactive_grain._var_index << ")\n"
634 13 : << COLOR_DEFAULT;
635 11 : 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 11 : 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 7046 : new_grain_index_to_existing_grain_index[closest_match_index] = old_grain_index;
649 : }
650 : }
651 :
652 : // Mark all resolved grain matches
653 7859 : for (const auto new_index : index_range(new_grain_index_to_existing_grain_index))
654 : {
655 7067 : 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 7067 : if (curr_index == invalid_size_t)
659 21 : continue;
660 :
661 : mooseAssert(_feature_sets_old[curr_index]._id != invalid_id,
662 : "Invalid ID in old grain structure");
663 :
664 7046 : _feature_sets[new_index]._id = _feature_sets_old[curr_index]._id; // Transfer ID
665 7046 : _feature_sets[new_index]._status = Status::MARKED; // Mark the status in the new set
666 7046 : _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 7859 : for (const auto grain_num : index_range(_feature_sets))
685 : {
686 7067 : auto & grain = _feature_sets[grain_num];
687 :
688 : // New Grain
689 7067 : 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 21 : std::lower_bound(_feature_sets.begin(), _feature_sets.end(), grain._var_index,
712 : [](const FeatureData & item, std::size_t var_index)
713 : {
714 48 : return item._var_index < var_index;
715 : });
716 : // clang-format on
717 :
718 : // Loop over matching variable indices
719 57 : for (MooseIndex(_feature_sets)
720 21 : new_grain_index = std::distance(_feature_sets.begin(), start_it);
721 57 : new_grain_index < _feature_sets.size() &&
722 45 : _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 15 : if (grain_num != new_grain_index && // Make sure indices aren't pointing at the same grain
729 30 : other_grain._status == Status::MARKED && // and that the other grain is indeed marked
730 60 : other_grain.boundingBoxesIntersect(grain) && // and the bboxes intersect
731 9 : 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 9 : grain._id = other_grain._id; // Set the duplicate ID
735 9 : grain._status = Status::MARKED; // Mark it
736 :
737 9 : if (_verbosity_level > 0)
738 11 : _console << COLOR_YELLOW << "Split Grain Detected #" << grain._id
739 9 : << " (variable index: " << grain._var_index << ")\n"
740 11 : << COLOR_DEFAULT;
741 9 : if (_verbosity_level > 1)
742 0 : _console << grain << other_grain;
743 : }
744 : }
745 :
746 21 : 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 9 : 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 60 : for (const auto old_grain_index : index_range(_feature_sets_old))
775 : {
776 : auto & old_grain = _feature_sets_old[old_grain_index];
777 :
778 51 : 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 18 : if (grain._var_index == old_grain._var_index &&
787 60 : grain.boundingBoxesIntersect(old_grain) && grain.halosIntersect(old_grain))
788 9 : old_grain_indices.push_back(old_grain_index);
789 : }
790 :
791 9 : if (old_grain_indices.size() == 1)
792 : {
793 9 : grain._id = _feature_sets_old[old_grain_indices[0]]._id;
794 9 : grain._status = Status::MARKED;
795 :
796 9 : if (_verbosity_level > 0)
797 11 : _console << COLOR_YELLOW << "Split Grain Detected #" << grain._id
798 9 : << " (variable index: " << grain._var_index << ")\n"
799 11 : << 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 9 : }
807 :
808 : // Must be a nucleating grain (status is still not set)
809 21 : if (grain._status == Status::CLEAR)
810 : {
811 12 : auto new_index = getNextUniqueID();
812 12 : grain._id = new_index; // Set the ID
813 12 : grain._status = Status::MARKED; // Mark it
814 :
815 12 : if (_verbosity_level > 0)
816 14 : _console << COLOR_YELLOW << "Nucleating Grain Detected "
817 12 : << " (variable index: " << grain._var_index << ")\n"
818 14 : << COLOR_DEFAULT;
819 12 : if (_verbosity_level > 1)
820 0 : _console << grain;
821 : }
822 : }
823 : }
824 :
825 : // Case 2 (inactive grains in _feature_sets_old)
826 7915 : for (auto & grain : _feature_sets_old)
827 : {
828 7123 : if (grain._status == Status::CLEAR)
829 : {
830 66 : grain._status = Status::INACTIVE;
831 66 : if (_verbosity_level > 0)
832 : {
833 79 : _console << COLOR_GREEN << "Marking Grain " << grain._id
834 66 : << " as INACTIVE (variable index: " << grain._var_index << ")\n"
835 79 : << COLOR_DEFAULT;
836 66 : if (_verbosity_level > 1)
837 0 : _console << grain;
838 : }
839 : }
840 : }
841 792 : } // is_primary
842 :
843 : /*************************************************************
844 : ****************** COLLECTIVE WORK SECTION ******************
845 : *************************************************************/
846 :
847 : // Make IDs on all non-primary ranks consistent
848 2169 : scatterAndUpdateRanks();
849 :
850 : // Build up an id to index map
851 2169 : _communicator.broadcast(_max_curr_grain_id);
852 2169 : buildFeatureIdToLocalIndices(_max_curr_grain_id);
853 :
854 : /**
855 : * Trigger callback for new grains
856 : */
857 2169 : if (_old_max_grain_id < _max_curr_grain_id)
858 : {
859 21 : 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 13 : if (new_id >= _reserve_grain_first_index + _n_reserve_ops)
863 : {
864 : // See if we've been instructed to terminate with an error
865 9 : 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 9 : newGrainCreated(new_id);
870 : }
871 : }
872 : }
873 2169 : }
874 :
875 : void
876 4567 : GrainTracker::newGrainCreated(unsigned int new_grain_id)
877 : {
878 4567 : 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 6 : 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 6 : const auto & grain = _feature_sets[grain_index];
886 7 : _console << COLOR_YELLOW
887 6 : << "\n*****************************************************************************"
888 6 : << "\nCouldn't find a matching grain while working on variable index: "
889 6 : << grain._var_index << "\nCreating new unique grain: " << new_grain_id << '\n'
890 6 : << grain
891 6 : << "\n*****************************************************************************\n"
892 7 : << COLOR_DEFAULT;
893 : }
894 4567 : }
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 2607 : GrainTracker::remapGrains()
910 : {
911 : // Don't remap grains if the current simulation step is before the specified tracking step
912 2607 : if (_t_step < _tracking_step)
913 0 : return;
914 :
915 5214 : TIME_SECTION("remapGrains", 3, "Remapping Grains");
916 :
917 2607 : 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 2607 : 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 10487 : for (auto & grain : _feature_sets)
944 : {
945 : // Unmark the grain so it can be used in the remap loop
946 9462 : grain._status = Status::CLEAR;
947 :
948 9462 : 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 10487 : for (const auto i : index_range(_feature_sets))
953 : {
954 9462 : auto & grain1 = _feature_sets[i];
955 :
956 147698 : for (const auto j : index_range(_feature_sets))
957 : {
958 138236 : auto & grain2 = _feature_sets[j];
959 138236 : if (i == j)
960 9462 : continue;
961 :
962 : // The first condition below is there to prevent symmetric checks (duplicate values)
963 128774 : if (i < j && grain1._id == grain2._id)
964 : {
965 10 : split_pairs.push_front(std::make_pair(i, j));
966 10 : 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 1163 : do
993 : {
994 : grains_remapped = false;
995 : notify_ids.clear();
996 :
997 12185 : for (auto & grain1 : _feature_sets)
998 : {
999 : // We need to remap any grains represented on any variable index above the cuttoff
1000 11022 : if (grain1._var_index >= _reserve_op_index)
1001 : {
1002 12 : if (_verbosity_level > 0)
1003 14 : _console << COLOR_YELLOW << "\nGrain #" << grain1._id
1004 12 : << " detected on a reserved order parameter #" << grain1._var_index
1005 12 : << ", remapping to another variable\n"
1006 14 : << COLOR_DEFAULT;
1007 :
1008 12 : for (const auto max : make_range(0, _max_remap_recursion_depth + 1))
1009 12 : if (max < _max_remap_recursion_depth)
1010 : {
1011 12 : 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 168134 : for (auto & grain2 : _feature_sets)
1031 : {
1032 : // Don't compare a grain with itself and don't try to remap inactive grains
1033 157112 : if (&grain1 == &grain2)
1034 11022 : continue;
1035 :
1036 159994 : if (grain1._var_index == grain2._var_index && // grains represented by same variable?
1037 27788 : grain1._id != grain2._id && // are they part of different grains?
1038 165577 : grain1.boundingBoxesIntersect(grain2) && // do bboxes intersect (coarse level)?
1039 5603 : grain1.halosIntersect(grain2)) // do they actually overlap (fine level)?
1040 : {
1041 212 : if (_verbosity_level > 0)
1042 256 : _console << COLOR_YELLOW << "Grain #" << grain1._id << " intersects Grain #"
1043 212 : << grain2._id << " (variable index: " << grain1._var_index << ")\n"
1044 256 : << COLOR_DEFAULT;
1045 :
1046 404 : for (const auto max : make_range(0, _max_remap_recursion_depth + 1))
1047 : {
1048 383 : if (max < _max_remap_recursion_depth)
1049 : {
1050 362 : 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 1025 : 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 10487 : 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 9462 : auto old_var_index = grain_id_to_existing_var_index[grain._id];
1102 :
1103 9462 : if (old_var_index != grain._var_index)
1104 : {
1105 : mooseAssert(static_cast<bool>(grain._status & Status::DIRTY), "grain status is incorrect");
1106 :
1107 257 : grain_id_to_new_var.emplace_hint(
1108 : grain_id_to_new_var.end(),
1109 257 : 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 257 : grain._var_index = old_var_index;
1120 : // Clear the DIRTY status as well for consistency
1121 : grain._status &= ~Status::DIRTY;
1122 : }
1123 : }
1124 :
1125 1025 : if (!grain_id_to_new_var.empty())
1126 : {
1127 141 : 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 2607 : _communicator.broadcast(grain_id_to_new_var);
1141 :
1142 : // Perform swaps if any occurred
1143 2607 : if (!grain_id_to_new_var.empty())
1144 : {
1145 : // Cache for holding values during swaps
1146 490 : std::vector<std::map<Node *, CacheValues>> cache(_n_vars);
1147 :
1148 : // Perform the actual swaps on all processors
1149 4050 : for (auto & grain : _feature_sets)
1150 : {
1151 : // See if this grain was remapped
1152 3560 : auto new_var_it = grain_id_to_new_var.find(grain._id);
1153 3560 : if (new_var_it != grain_id_to_new_var.end())
1154 557 : swapSolutionValues(grain, new_var_it->second, cache, RemapCacheMode::FILL);
1155 : }
1156 :
1157 4050 : for (auto & grain : _feature_sets)
1158 : {
1159 : // See if this grain was remapped
1160 3560 : auto new_var_it = grain_id_to_new_var.find(grain._id);
1161 3560 : if (new_var_it != grain_id_to_new_var.end())
1162 557 : swapSolutionValues(grain, new_var_it->second, cache, RemapCacheMode::USE);
1163 : }
1164 :
1165 490 : _sys.solution().close();
1166 490 : _sys.solutionOld().close();
1167 490 : _sys.solutionOlder().close();
1168 :
1169 490 : _sys.system().update();
1170 :
1171 490 : if (_verbosity_level > 1)
1172 0 : _console << "Swaps complete" << std::endl;
1173 490 : }
1174 : }
1175 :
1176 : void
1177 845 : 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 13217 : for (const auto i : index_range(_feature_sets))
1204 : {
1205 12372 : auto & other_grain = _feature_sets[i];
1206 :
1207 12372 : if (other_grain._var_index == grain._var_index || other_grain._var_index >= _reserve_op_index)
1208 2118 : continue;
1209 :
1210 : auto target_var_index = other_grain._var_index;
1211 : auto target_grain_index = i;
1212 10254 : auto target_grain_id = other_grain._id;
1213 :
1214 10254 : Real curr_bbox_diff = boundingRegionDistance(grain._bboxes, other_grain._bboxes);
1215 :
1216 : GrainDistance grain_distance_obj(
1217 10254 : 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 10254 : if (curr_bbox_diff == -1.0 && !min_distances[target_var_index].empty())
1221 : {
1222 3886 : Real last_distance = min_distances[target_var_index].begin()->_distance;
1223 3886 : if (last_distance < 0)
1224 3207 : 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 11049 : while (insert_it != min_distances[target_var_index].end() && !(grain_distance_obj < *insert_it))
1230 : ++insert_it;
1231 10254 : 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 7521 : 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 6676 : if (grain._var_index == var_index)
1243 833 : continue;
1244 :
1245 5843 : if (min_distances[var_index].empty())
1246 189 : min_distances[var_index].emplace_front(std::numeric_limits<Real>::max(), var_index);
1247 : }
1248 845 : }
1249 :
1250 : bool
1251 1186 : GrainTracker::attemptGrainRenumber(FeatureData & grain, unsigned int depth, unsigned int max_depth)
1252 : {
1253 : // End the recursion of our breadth first search
1254 1186 : if (depth > max_depth)
1255 : return false;
1256 :
1257 845 : std::size_t curr_var_index = grain._var_index;
1258 :
1259 : std::vector<std::map<Node *, CacheValues>> cache;
1260 :
1261 845 : 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 845 : 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 845 : 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 16426 : if (lhs.empty())
1287 : return false;
1288 14854 : else if (rhs.empty())
1289 : return true;
1290 : else
1291 12351 : return lhs.begin()->_distance > rhs.begin()->_distance;
1292 : });
1293 : // clang-format on
1294 :
1295 1567 : for (auto & list_ref : min_distances)
1296 : {
1297 : const auto target_it = list_ref.begin();
1298 1567 : if (target_it == list_ref.end())
1299 : continue;
1300 :
1301 : // If the distance is positive we can just remap and be done
1302 1567 : if (target_it->_distance > 0)
1303 : {
1304 140 : if (_verbosity_level > 0)
1305 : {
1306 167 : _console << COLOR_GREEN << "- Depth " << depth << ": Remapping grain #" << grain._id
1307 140 : << " from variable index " << curr_var_index << " to " << target_it->_var_index;
1308 140 : if (target_it->_distance == std::numeric_limits<Real>::max())
1309 118 : _console << " which currently contains zero grains.\n\n" << COLOR_DEFAULT;
1310 : else
1311 40 : _console << " whose closest grain (#" << target_it->_grain_id << ") is at a distance of "
1312 40 : << std::sqrt(target_it->_distance) << "\n\n"
1313 49 : << COLOR_DEFAULT;
1314 : }
1315 :
1316 : grain._status |= Status::DIRTY;
1317 140 : grain._var_index = target_it->_var_index;
1318 140 : 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 1427 : std::ostringstream oss;
1328 2892 : while (!intersection_hit && next_target_it != list_ref.end())
1329 : {
1330 1465 : 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 1465 : FeatureData & next_target_grain = _feature_sets[next_target_it->_grain_index];
1336 :
1337 : // If any grains touch we're done here
1338 1465 : if (grain.halosIntersect(next_target_grain))
1339 : intersection_hit = true;
1340 : else
1341 : {
1342 107 : if (num_close_targets > 0)
1343 9 : oss << ", "; // delimiter
1344 107 : oss << "#" << next_target_it->_grain_id;
1345 : }
1346 :
1347 : ++next_target_it;
1348 1465 : ++num_close_targets;
1349 : }
1350 :
1351 1427 : if (!intersection_hit)
1352 : {
1353 69 : if (_verbosity_level > 0)
1354 : {
1355 83 : _console << COLOR_GREEN << "- Depth " << depth << ": Remapping grain #" << grain._id
1356 69 : << " from variable index " << curr_var_index << " to " << target_it->_var_index;
1357 :
1358 69 : if (num_close_targets == 1)
1359 180 : _console << " whose closest grain (" << oss.str()
1360 60 : << ") is inside our bounding box but whose halo is not touching.\n\n"
1361 72 : << COLOR_DEFAULT;
1362 : else
1363 27 : _console << " whose closest grains (" << oss.str()
1364 9 : << ") are inside our bounding box but whose halos are not touching.\n\n"
1365 11 : << COLOR_DEFAULT;
1366 : }
1367 :
1368 : grain._status |= Status::DIRTY;
1369 69 : grain._var_index = target_it->_var_index;
1370 : return true;
1371 : }
1372 1427 : }
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 1358 : 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 1358 : 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 1067 : 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 770 : grain._var_index = target_it->_var_index;
1399 : grain._status |= Status::MARKED;
1400 770 : if (attemptGrainRenumber(target_grain, depth + 1, max_depth))
1401 : {
1402 : // SUCCESS!
1403 48 : if (_verbosity_level > 0)
1404 60 : _console << COLOR_GREEN << "- Depth " << depth << ": Remapping grain #" << grain._id
1405 48 : << " from variable index " << curr_var_index << " to " << target_it->_var_index
1406 48 : << "\n\n"
1407 60 : << COLOR_DEFAULT;
1408 :
1409 : // Now we need to mark the grain as DIRTY since the recursion succeeded.
1410 : grain._status |= Status::DIRTY;
1411 48 : return true;
1412 : }
1413 : else
1414 : // FAILURE, We need to set our var index back after failed recursive step
1415 722 : 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 845 : }
1424 :
1425 : void
1426 1114 : 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 1114 : MeshBase & mesh = _mesh.getMesh();
1432 :
1433 : // Remap the grain
1434 : std::set<Node *> updated_nodes_tmp; // Used only in the elemental case
1435 54736 : for (auto entity : grain._local_ids)
1436 : {
1437 53622 : if (_is_elemental)
1438 : {
1439 53622 : Elem * elem = mesh.query_elem_ptr(entity);
1440 53622 : if (!elem)
1441 0 : continue;
1442 :
1443 268110 : for (unsigned int i = 0; i < elem->n_nodes(); ++i)
1444 : {
1445 214488 : Node * curr_node = elem->node_ptr(i);
1446 214488 : 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 69744 : updated_nodes_tmp.insert(curr_node);
1450 69744 : 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 1114 : if (cache_mode == RemapCacheMode::USE || cache_mode == RemapCacheMode::BYPASS)
1461 557 : grain._var_index = new_var_index;
1462 1114 : }
1463 :
1464 : void
1465 69744 : 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 69744 : 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 58810 : _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 58810 : if (cache_mode == RemapCacheMode::FILL || cache_mode == RemapCacheMode::BYPASS)
1480 : {
1481 29405 : current = _vars[curr_var_index]->dofValues()[0];
1482 29405 : if (_is_transient)
1483 : {
1484 29405 : old = _vars[curr_var_index]->dofValuesOld()[0];
1485 29405 : 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 29405 : current = cache_it->second.current;
1493 29405 : old = cache_it->second.old;
1494 29405 : older = cache_it->second.older;
1495 : }
1496 :
1497 : // Cache the value or use it!
1498 58810 : if (cache_mode == RemapCacheMode::FILL)
1499 : {
1500 29405 : cache[curr_var_index][curr_node].current = current;
1501 29405 : cache[curr_var_index][curr_node].old = old;
1502 29405 : cache[curr_var_index][curr_node].older = older;
1503 : }
1504 : else // USE or BYPASS
1505 : {
1506 29405 : const auto & dof_index = _vars[new_var_index]->nodalDofIndex();
1507 :
1508 : // Transfer this solution from the old to the current
1509 29405 : _sys.solution().set(dof_index, current);
1510 29405 : if (_is_transient)
1511 : {
1512 29405 : _sys.solutionOld().set(dof_index, old);
1513 29405 : _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 58810 : if (cache_mode == RemapCacheMode::FILL || cache_mode == RemapCacheMode::BYPASS)
1530 : {
1531 29405 : const auto & dof_index = _vars[curr_var_index]->nodalDofIndex();
1532 :
1533 : // Set the DOF for the current variable to zero
1534 29405 : _sys.solution().set(dof_index, -_bound_value);
1535 29405 : if (_is_transient)
1536 : {
1537 29405 : _sys.solutionOld().set(dof_index, -_bound_value);
1538 29405 : _sys.solutionOlder().set(dof_index, -_bound_value);
1539 : }
1540 : }
1541 : }
1542 69744 : }
1543 :
1544 : void
1545 2616 : GrainTracker::updateFieldInfo()
1546 : {
1547 5232 : TIME_SECTION("updateFieldInfo", 3, "Updating Field Info");
1548 :
1549 20480 : 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 20626 : for (const auto & grain : _feature_sets)
1555 : {
1556 18010 : std::size_t curr_var = grain._var_index;
1557 18010 : std::size_t map_index = (_single_map_mode || _condense_map_info) ? 0 : curr_var;
1558 :
1559 1128110 : 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 1110100 : if (_is_elemental)
1564 : {
1565 1095132 : const Elem * elem = _mesh.elemPtr(entity);
1566 1095132 : std::vector<Point> centroid(1, elem->vertex_average());
1567 1095132 : if (_poly_ic_uo && _first_time)
1568 : {
1569 182961 : entity_value = _poly_ic_uo->getVariableValue(grain._var_index, centroid[0]);
1570 : }
1571 : else
1572 : {
1573 912171 : _fe_problem.reinitElemPhys(elem, centroid, 0);
1574 912171 : entity_value = _vars[curr_var]->sln()[0];
1575 : }
1576 1095132 : }
1577 : else
1578 : {
1579 14968 : auto node_ptr = _mesh.nodePtr(entity);
1580 14968 : entity_value = _vars[curr_var]->getNodalValue(*node_ptr);
1581 : }
1582 :
1583 1110100 : if (entity_value != std::numeric_limits<Real>::lowest() &&
1584 410556 : (tmp_map.find(entity) == tmp_map.end() || entity_value > tmp_map[entity]))
1585 : {
1586 : mooseAssert(grain._id != invalid_id, "Missing Grain ID");
1587 874757 : _feature_maps[map_index][entity] = grain._id;
1588 :
1589 874757 : if (_var_index_mode)
1590 874757 : _var_index_maps[map_index][entity] = grain._var_index;
1591 :
1592 874757 : tmp_map[entity] = entity_value;
1593 : }
1594 :
1595 1110100 : if (_compute_var_to_feature_map)
1596 : {
1597 204010 : auto insert_pair = moose_try_emplace(
1598 204010 : _entity_var_to_features, entity, std::vector<unsigned int>(_n_vars, invalid_id));
1599 : auto & vec_ref = insert_pair.first->second;
1600 :
1601 204010 : if (insert_pair.second)
1602 : {
1603 : // insert the reserve op numbers (if appropriate)
1604 146171 : for (const auto reserve_index : make_range(_n_reserve_ops))
1605 0 : vec_ref[reserve_index] = _reserve_grain_first_index + reserve_index;
1606 : }
1607 204010 : vec_ref[grain._var_index] = grain._id;
1608 : }
1609 : }
1610 :
1611 18010 : if (_compute_halo_maps)
1612 1166115 : for (auto entity : grain._halo_ids)
1613 1151167 : _halo_ids[grain._var_index][entity] = grain._var_index;
1614 :
1615 533661 : for (auto entity : grain._ghosted_ids)
1616 515651 : _ghosted_entity_ids[entity] = 1;
1617 : }
1618 :
1619 2616 : communicateHaloMap();
1620 2616 : }
1621 :
1622 : void
1623 2616 : GrainTracker::communicateHaloMap()
1624 : {
1625 2616 : 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 2036 : const bool isDistributedMesh = _mesh.isDistributedMesh();
1635 :
1636 2036 : if (_is_primary)
1637 : {
1638 631 : std::vector<std::vector<std::pair<std::size_t, dof_id_type>>> root_halo_ids(_n_procs);
1639 631 : counts.resize(_n_procs);
1640 :
1641 : // Loop over the _halo_ids "field" and build minimal lists for all of the other ranks
1642 5243 : for (const auto var_index : index_range(_halo_ids))
1643 : {
1644 390771 : for (const auto & entity_pair : _halo_ids[var_index])
1645 : {
1646 386159 : auto entity_id = entity_pair.first;
1647 386159 : if (isDistributedMesh)
1648 : {
1649 : // Check to see which contiguous range this entity ID falls into
1650 : auto range_it =
1651 30301 : 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 30301 : root_halo_ids[proc_id].push_back(std::make_pair(var_index, entity_id));
1664 : }
1665 : else
1666 : {
1667 : DofObject * halo_entity;
1668 355858 : if (_is_elemental)
1669 355858 : halo_entity = _mesh.queryElemPtr(entity_id);
1670 : else
1671 0 : halo_entity = _mesh.queryNodePtr(entity_id);
1672 :
1673 : if (halo_entity)
1674 355858 : root_halo_ids[halo_entity->processor_id()].push_back(
1675 355858 : std::make_pair(var_index, entity_id));
1676 : }
1677 : }
1678 : }
1679 :
1680 : // Build up the counts vector for MPI scatter
1681 2667 : 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 2036 : counts[counter] = vector_ref.size();
1685 2036 : counter++;
1686 : }
1687 631 : }
1688 :
1689 2036 : _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 388195 : for (const auto & halo_pair : local_halo_ids)
1693 386159 : _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 16984 : for (const auto & grain : _feature_sets)
1702 782716 : for (auto local_id : grain._local_ids)
1703 767768 : _halo_ids[grain._var_index].erase(local_id);
1704 2036 : }
1705 2616 : }
1706 :
1707 : Real
1708 10096 : 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 20930 : for (const auto & bbox1 : bboxes1)
1716 : {
1717 : const auto centroid_point1 = (bbox1.max() + bbox1.min()) / 2.0;
1718 :
1719 23509 : 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 12675 : auto curr_distance = _mesh.minPeriodicDistance(_var_number, centroid_point1, centroid_point2);
1725 :
1726 12675 : if (curr_distance < min_distance)
1727 : min_distance = curr_distance;
1728 : }
1729 : }
1730 :
1731 10096 : return min_distance;
1732 : }
1733 :
1734 : Real
1735 10254 : GrainTracker::boundingRegionDistance(std::vector<BoundingBox> & bboxes1,
1736 : std::vector<BoundingBox> & bboxes2) const
1737 : {
1738 : /**
1739 : * The region that each grain covers is represented by a bounding box large enough to encompassing
1740 : * all the points within that grain. When using periodic boundaries, we may have several discrete
1741 : * "pieces" of a grain each represented by a bounding box. The distance between any two grains
1742 : * is defined as the minimum distance between any pair of boxes, one selected from each grain.
1743 : */
1744 : auto min_distance = std::numeric_limits<Real>::max();
1745 11719 : for (const auto & bbox1 : bboxes1)
1746 : {
1747 11719 : for (const auto & bbox2 : bboxes2)
1748 : {
1749 : // AABB squared distance
1750 : Real curr_distance = 0.0;
1751 : bool boxes_overlap = true;
1752 41016 : for (unsigned int dim = 0; dim < LIBMESH_DIM; ++dim)
1753 : {
1754 : const auto & min1 = bbox1.min()(dim);
1755 : const auto & max1 = bbox1.max()(dim);
1756 : const auto & min2 = bbox2.min()(dim);
1757 : const auto & max2 = bbox2.max()(dim);
1758 :
1759 30762 : if (min1 > max2)
1760 : {
1761 971 : const auto delta = max2 - min1;
1762 971 : curr_distance += delta * delta;
1763 : boxes_overlap = false;
1764 : }
1765 29791 : else if (min2 > max1)
1766 : {
1767 506 : const auto delta = max1 - min2;
1768 506 : curr_distance += delta * delta;
1769 : boxes_overlap = false;
1770 : }
1771 : }
1772 :
1773 10254 : if (boxes_overlap)
1774 : return -1.0; /* all overlaps are treated the same */
1775 :
1776 1465 : if (curr_distance < min_distance)
1777 : min_distance = curr_distance;
1778 : }
1779 : }
1780 :
1781 : return min_distance;
1782 : }
1783 :
1784 : unsigned int
1785 12 : GrainTracker::getNextUniqueID()
1786 : {
1787 : /**
1788 : * Get the next unique grain ID but make sure to respect
1789 : * reserve ids. Note, that the first valid ID for a new
1790 : * grain is _reserve_grain_first_index + _n_reserve_ops because
1791 : * _reserve_grain_first_index IS a valid index. It does not
1792 : * point to the last valid index of the non-reserved grains.
1793 : */
1794 12 : _max_curr_grain_id = std::max(_max_curr_grain_id == invalid_id ? 0 : _max_curr_grain_id + 1,
1795 12 : _reserve_grain_first_index + _n_reserve_ops /* no +1 here!*/);
1796 :
1797 12 : return _max_curr_grain_id;
1798 : }
1799 :
1800 : /*************************************************
1801 : ************** Helper Structures ****************
1802 : ************************************************/
1803 189 : GrainDistance::GrainDistance(Real distance, std::size_t var_index)
1804 : : GrainDistance(distance,
1805 : var_index,
1806 : std::numeric_limits<std::size_t>::max(),
1807 189 : std::numeric_limits<unsigned int>::max())
1808 : {
1809 189 : }
1810 :
1811 10443 : GrainDistance::GrainDistance(Real distance,
1812 : std::size_t var_index,
1813 : std::size_t grain_index,
1814 10443 : unsigned int grain_id)
1815 10443 : : _distance(distance), _var_index(var_index), _grain_index(grain_index), _grain_id(grain_id)
1816 : {
1817 10443 : }
1818 :
1819 : bool
1820 4684 : GrainDistance::operator<(const GrainDistance & rhs) const
1821 : {
1822 4684 : return _distance < rhs._distance;
1823 : }
|
