www.mooseframework.org
FeatureFloodCount.C
Go to the documentation of this file.
1 //* This file is part of the MOOSE framework
2 //* https://www.mooseframework.org
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 "FeatureFloodCount.h"
11 #include "IndirectSort.h"
12 #include "MooseMesh.h"
13 #include "MooseUtils.h"
14 #include "MooseVariable.h"
15 #include "SubProblem.h"
16 
17 #include "Assembly.h"
18 #include "FEProblem.h"
19 #include "NonlinearSystem.h"
20 
21 #include "libmesh/dof_map.h"
22 #include "libmesh/mesh_tools.h"
23 #include "libmesh/periodic_boundaries.h"
24 #include "libmesh/point_locator_base.h"
25 #include "libmesh/remote_elem.h"
26 
27 #include <algorithm>
28 #include <limits>
29 
30 template <>
31 void
32 dataStore(std::ostream & stream, FeatureFloodCount::FeatureData & feature, void * context)
33 {
38  storeHelper(stream, feature._ghosted_ids, context);
39  storeHelper(stream, feature._halo_ids, context);
40  storeHelper(stream, feature._disjoint_halo_ids, context);
41  storeHelper(stream, feature._periodic_nodes, context);
42  storeHelper(stream, feature._var_index, context);
43  storeHelper(stream, feature._id, context);
44  storeHelper(stream, feature._bboxes, context);
45  storeHelper(stream, feature._orig_ids, context);
46  storeHelper(stream, feature._min_entity_id, context);
47  storeHelper(stream, feature._vol_count, context);
48  storeHelper(stream, feature._centroid, context);
49  storeHelper(stream, feature._status, context);
50  storeHelper(stream, feature._boundary_intersection, context);
51 }
52 
53 template <>
54 void
55 dataStore(std::ostream & stream, BoundingBox & bbox, void * context)
56 {
57  storeHelper(stream, bbox.min(), context);
58  storeHelper(stream, bbox.max(), context);
59 }
60 
61 template <>
62 void
63 dataLoad(std::istream & stream, FeatureFloodCount::FeatureData & feature, void * context)
64 {
69  loadHelper(stream, feature._ghosted_ids, context);
70  loadHelper(stream, feature._halo_ids, context);
71  loadHelper(stream, feature._disjoint_halo_ids, context);
72  loadHelper(stream, feature._periodic_nodes, context);
73  loadHelper(stream, feature._var_index, context);
74  loadHelper(stream, feature._id, context);
75  loadHelper(stream, feature._bboxes, context);
76  loadHelper(stream, feature._orig_ids, context);
77  loadHelper(stream, feature._min_entity_id, context);
78  loadHelper(stream, feature._vol_count, context);
79  loadHelper(stream, feature._centroid, context);
80  loadHelper(stream, feature._status, context);
81  loadHelper(stream, feature._boundary_intersection, context);
82 }
83 
84 template <>
85 void
86 dataLoad(std::istream & stream, BoundingBox & bbox, void * context)
87 {
88  loadHelper(stream, bbox.min(), context);
89  loadHelper(stream, bbox.max(), context);
90 }
91 
92 // Utility routines
93 void updateBBoxExtremesHelper(BoundingBox & bbox, const Point & node);
94 void updateBBoxExtremesHelper(BoundingBox & bbox, const Elem & elem);
95 
96 registerMooseObject("PhaseFieldApp", FeatureFloodCount);
97 
100 {
103 
104  params.addRequiredCoupledVar(
105  "variable",
106  "The variable(s) for which to find connected regions of interests, i.e. \"features\".");
107  params.addParam<Real>(
108  "threshold", 0.5, "The threshold value for which a new feature may be started");
109  params.addParam<Real>(
110  "connecting_threshold",
111  "The threshold for which an existing feature may be extended (defaults to \"threshold\")");
112  params.addParam<bool>("use_single_map",
113  true,
114  "Determine whether information is tracked per "
115  "coupled variable or consolidated into one "
116  "(default: true)");
117  params.addParam<bool>(
118  "condense_map_info",
119  false,
120  "Determines whether we condense all the node values when in multimap mode (default: false)");
121  params.addParam<bool>("use_global_numbering",
122  true,
123  "Determine whether or not global numbers are "
124  "used to label features on multiple maps "
125  "(default: true)");
126  params.addParam<bool>("enable_var_coloring",
127  false,
128  "Instruct the Postprocessor to populate the variable index map.");
129  params.addParam<bool>(
130  "compute_halo_maps",
131  false,
132  "Instruct the Postprocessor to communicate proper halo information to all ranks");
133  params.addParam<bool>("compute_var_to_feature_map",
134  false,
135  "Instruct the Postprocessor to compute the active vars to features map");
136  params.addParam<bool>(
137  "use_less_than_threshold_comparison",
138  true,
139  "Controls whether features are defined to be less than or greater than the threshold value.");
140 
141  params.addParam<std::vector<BoundaryName>>(
142  "primary_percolation_boundaries",
143  "A list of boundaries used in conjunction with the corresponding "
144  "\"secondary_percolation_boundaries\" parameter for determining if a feature creates a path "
145  "connecting any pair of boundaries");
146  params.addParam<std::vector<BoundaryName>>(
147  "secondary_percolation_boundaries",
148  "Paired boundaries with \"primaryary_percolation_boundaries\" parameter");
149  params.addParam<std::vector<BoundaryName>>(
150  "specified_boundaries",
151  "An optional list of boundaries; if supplied, each feature is checked to determine whether "
152  "it intersects any of the specified boundaries in this list.");
153 
161  params.set<bool>("use_displaced_mesh") = false;
162 
163  // The FeatureFloodCount object does not require that any state (restartable information) is
164  // maintained. This Boolean is set to false so that we don't ask MOOSE to save a potentially
165  // large data structure for no reason. It is set for true in at least one derived class
166  // (GrainTracker).
167  params.addPrivateParam<bool>("restartable_required", false);
168 
169  params.addParamNamesToGroup(
170  "use_single_map condense_map_info use_global_numbering primary_percolation_boundaries",
171  "Advanced");
172 
173  MooseEnum flood_type("NODAL ELEMENTAL", "ELEMENTAL");
174  params.addParam<MooseEnum>("flood_entity_type",
175  flood_type,
176  "Determines whether the flood algorithm runs on nodes or elements");
177 
178  params.addClassDescription("The object is able to find and count \"connected components\" in any "
179  "solution field or number of solution fields. A primary example would "
180  "be to count \"bubbles\".");
181 
182  params.addRelationshipManager("ElementSideNeighborLayers",
185 
186  return params;
187 }
188 
190  : GeneralPostprocessor(parameters),
191  Coupleable(this, false),
193  BoundaryRestrictable(this, false),
194  _fe_vars(getCoupledMooseVars()),
195  _vars(getCoupledStandardMooseVars()),
196  _dof_map(_vars[0]->dofMap()),
197  _threshold(getParam<Real>("threshold")),
198  _connecting_threshold(isParamValid("connecting_threshold")
199  ? getParam<Real>("connecting_threshold")
200  : getParam<Real>("threshold")),
201  _mesh(_subproblem.mesh()),
202  _var_number(_fe_vars[0]->number()),
203  _single_map_mode(getParam<bool>("use_single_map")),
204  _condense_map_info(getParam<bool>("condense_map_info")),
205  _global_numbering(getParam<bool>("use_global_numbering")),
206  _var_index_mode(getParam<bool>("enable_var_coloring")),
207  _compute_halo_maps(getParam<bool>("compute_halo_maps")),
208  _compute_var_to_feature_map(getParam<bool>("compute_var_to_feature_map")),
209  _use_less_than_threshold_comparison(getParam<bool>("use_less_than_threshold_comparison")),
210  _n_vars(_fe_vars.size()),
211  _maps_size(_single_map_mode ? 1 : _fe_vars.size()),
212  _n_procs(_app.n_processors()),
213  _feature_counts_per_map(_maps_size),
214  _feature_count(0),
215  _partial_feature_sets(_maps_size),
216  _feature_sets(getParam<bool>("restartable_required")
217  ? declareRestartableData<std::vector<FeatureData>>("feature_sets")
218  : _volatile_feature_sets),
219  _feature_maps(_maps_size),
220  _pbs(nullptr),
221  _element_average_value(parameters.isParamValid("elem_avg_value")
222  ? getPostprocessorValue("elem_avg_value")
223  : _real_zero),
224  _halo_ids(_maps_size),
225  _is_elemental(getParam<MooseEnum>("flood_entity_type") == "ELEMENTAL"),
226  _is_primary(processor_id() == 0)
227 {
228  if (_var_index_mode)
229  _var_index_maps.resize(_maps_size);
230 
232 
234 
235  if (_subproblem.isTransient())
236  {
237  // tell MOOSE that we are going to need old and older DoF values
238  for (auto & var : _vars)
239  {
240  var->dofValuesOld();
241  var->dofValuesOlder();
242  }
243  }
244 
245  if (parameters.isParamValid("primary_percolation_boundaries"))
247  parameters.get<std::vector<BoundaryName>>("primary_percolation_boundaries"));
248  if (parameters.isParamValid("secondary_percolation_boundaries"))
250  parameters.get<std::vector<BoundaryName>>("secondary_percolation_boundaries"));
251 
252  if (_primary_perc_bnds.empty() != _secondary_perc_bnds.empty())
253  paramError("primary_percolation_boundaries",
254  "primary_percolation_boundaries and secondary_percolation_boundaries must both be "
255  "supplied when checking for percolation");
256 
257  if (parameters.isParamValid("specified_boundaries"))
259  _mesh.getBoundaryIDs(parameters.get<std::vector<BoundaryName>>("specified_boundaries"));
260 }
261 
262 void
264 {
265  // We need one map per coupled variable for normal runs to support overlapping features
266  _entities_visited.resize(_vars.size());
267 
268  // Get a pointer to the PeriodicBoundaries buried in libMesh
269  _pbs = _fe_problem.getNonlinearSystemBase(_sys.number()).dofMap().get_periodic_boundaries();
270 
271  meshChanged();
272 
281 }
282 
283 void
285 {
286  // Clear the feature marking maps and region counters and other data structures
287  for (const auto map_num : make_range(_maps_size))
288  {
289  _feature_maps[map_num].clear();
290  _partial_feature_sets[map_num].clear();
291 
292  if (_var_index_mode)
293  _var_index_maps[map_num].clear();
294 
295  _halo_ids[map_num].clear();
296  }
297 
298  _feature_sets.clear();
299 
300  // Calculate the thresholds for this iteration
303 
304  _ghosted_entity_ids.clear();
305 
306  // Reset the feature count and max local size
307  _feature_count = 0;
308 
309  _entity_var_to_features.clear();
310 
311  for (auto & map_ref : _entities_visited)
312  map_ref.clear();
313 }
314 
315 void
317 {
318 }
319 
320 void
322 {
323  _point_locator = _mesh.getMesh().sub_point_locator();
324 
326 
327  // Build a new node to element map
328  _nodes_to_elem_map.clear();
329  MeshTools::build_nodes_to_elem_map(_mesh.getMesh(), _nodes_to_elem_map);
330 
336  _all_boundary_entity_ids.clear();
337  if (_is_elemental)
338  for (auto elem_it = _mesh.bndElemsBegin(), elem_end = _mesh.bndElemsEnd(); elem_it != elem_end;
339  ++elem_it)
340  _all_boundary_entity_ids.insert((*elem_it)->_elem->id());
341 }
342 
343 void
345 {
346  TIME_SECTION("execute", 3, "Flooding Features");
347 
348  // Iterate only over boundaries if restricted
350  {
351  mooseInfo("Using EXPERIMENTAL boundary restricted FeatureFloodCount object!\n");
352 
353  // Set the boundary range pointer for use during flooding
355 
356  auto rank = processor_id();
357 
358  for (const auto & belem : *_bnd_elem_range)
359  {
360  const Elem * elem = belem->_elem;
361  BoundaryID boundary_id = belem->_bnd_id;
362 
363  if (elem->processor_id() == rank)
364  {
365  if (hasBoundary(boundary_id))
366  for (const auto var_num : index_range(_vars))
367  flood(elem, var_num);
368  }
369  }
370  }
371  else // Normal volumetric operation
372  {
373  for (const auto & current_elem : _mesh.getMesh().active_local_element_ptr_range())
374  {
375  // Loop over elements or nodes
376  if (_is_elemental)
377  {
378  for (const auto var_num : index_range(_vars))
379  flood(current_elem, var_num);
380  }
381  else
382  {
383  auto n_nodes = current_elem->n_vertices();
384  for (const auto i : make_range(n_nodes))
385  {
386  const Node * current_node = current_elem->node_ptr(i);
387 
388  for (const auto var_num : index_range(_vars))
389  flood(current_node, var_num);
390  }
391  }
392  }
393  }
394 }
395 
398 {
399  return _app.n_processors() >= _maps_size ? _maps_size : 1;
400 }
401 
402 void
404 {
405  TIME_SECTION("communicateAndMerge", 3, "Communicating and Merging");
406 
407  // First we need to transform the raw data into a usable data structure
409 
417  std::vector<std::string> send_buffers(1);
418 
425  std::vector<std::string> recv_buffers, deserialize_buffers;
426 
433  const auto n_merging_procs = numberOfDistributedMergeHelpers();
434 
435  if (n_merging_procs > 1)
436  {
437  auto rank = processor_id();
438  bool is_merging_processor = rank < n_merging_procs;
439 
440  if (is_merging_processor)
441  recv_buffers.reserve(_app.n_processors());
442 
443  for (const auto i : make_range(n_merging_procs))
444  {
445  serialize(send_buffers[0], i);
446 
452  (void *)(nullptr),
453  send_buffers.begin(),
454  send_buffers.end(),
455  std::back_inserter(recv_buffers));
456 
463  if (rank == i)
464  recv_buffers.swap(deserialize_buffers);
465  else
466  recv_buffers.clear();
467  }
468 
469  // Setup a new communicator for doing merging communication operations
470  Parallel::Communicator merge_comm;
471 
472  _communicator.split(is_merging_processor ? 0 : MPI_UNDEFINED, rank, merge_comm);
473 
474  if (is_merging_processor)
475  {
483  std::vector<std::list<FeatureData>> tmp_data(_partial_feature_sets.size());
484  tmp_data.swap(_partial_feature_sets);
485 
486  deserialize(deserialize_buffers, processor_id());
487 
488  send_buffers[0].clear();
489  recv_buffers.clear();
490  deserialize_buffers.clear();
491 
492  // Merge one variable's worth of data
493  mergeSets();
494 
495  // Now we need to serialize again to send to the primary (only the processors who did work)
496  serialize(send_buffers[0]);
497 
498  // Free up as much memory as possible here before we do global communication
500 
505  merge_comm.gather_packed_range(0,
506  (void *)(nullptr),
507  send_buffers.begin(),
508  send_buffers.end(),
509  std::back_inserter(recv_buffers));
510 
511  if (_is_primary)
512  {
513  // The root process now needs to deserialize all of the data
514  deserialize(recv_buffers);
515 
516  send_buffers[0].clear();
517  recv_buffers.clear();
518 
519  consolidateMergedFeatures(&tmp_data);
520  }
521  else
522  // Restore our original data on non-zero ranks
523  tmp_data.swap(_partial_feature_sets);
524  }
525  }
526 
527  // Serialized merging (primary does all the work)
528  else
529  {
530  if (_is_primary)
531  recv_buffers.reserve(_app.n_processors());
532 
533  serialize(send_buffers[0]);
534 
535  // Free up as much memory as possible here before we do global communication
537 
543  (void *)(nullptr),
544  send_buffers.begin(),
545  send_buffers.end(),
546  std::back_inserter(recv_buffers));
547 
548  if (_is_primary)
549  {
550  // The root process now needs to deserialize all of the data
551  deserialize(recv_buffers);
552  recv_buffers.clear();
553 
554  mergeSets();
555 
557  }
558  }
559 
560  if (!_is_primary)
562 
563  // Make sure that feature count is communicated to all ranks
565 }
566 
567 void
569 {
570  mooseAssert(_is_primary, "sortAndLabel can only be called on the primary");
571 
577  std::sort(_feature_sets.begin(), _feature_sets.end());
578 
579 #ifndef NDEBUG
580 
585  unsigned int feature_offset = 0;
586  for (const auto map_num : make_range(_maps_size))
587  {
588  // Skip empty map checks
589  if (_feature_counts_per_map[map_num] == 0)
590  continue;
591 
592  // Check the begin and end of the current range
593  auto range_front = feature_offset;
594  auto range_back = feature_offset + _feature_counts_per_map[map_num] - 1;
595 
596  mooseAssert(range_front <= range_back && range_back < _feature_count,
597  "Indexing error in feature sets");
598 
599  if (!_single_map_mode && (_feature_sets[range_front]._var_index != map_num ||
600  _feature_sets[range_back]._var_index != map_num))
601  mooseError("Error in _feature_sets sorting, map index: ", map_num);
602 
603  feature_offset += _feature_counts_per_map[map_num];
604  }
605 #endif
606 
607  // Label the features with an ID based on the sorting (processor number independent value)
608  for (const auto i : index_range(_feature_sets))
609  if (_feature_sets[i]._id == invalid_id)
610  _feature_sets[i]._id = i;
611 }
612 
613 void
614 FeatureFloodCount::buildLocalToGlobalIndices(std::vector<std::size_t> & local_to_global_all,
615  std::vector<int> & counts) const
616 {
617  mooseAssert(_is_primary, "This method must only be called on the root processor");
618 
619  counts.assign(_n_procs, 0);
620  // Now size the individual counts vectors based on the largest index seen per processor
621  for (const auto & feature : _feature_sets)
622  for (const auto & local_index_pair : feature._orig_ids)
623  {
624  // local_index_pair.first = ranks, local_index_pair.second = local_index
625  mooseAssert(local_index_pair.first < _n_procs, "Processor ID is out of range");
626  if (local_index_pair.second >= static_cast<std::size_t>(counts[local_index_pair.first]))
627  counts[local_index_pair.first] = local_index_pair.second + 1;
628  }
629 
630  // Build the offsets vector
631  unsigned int globalsize = 0;
632  std::vector<int> offsets(_n_procs); // Type is signed for use with the MPI API
633  for (const auto i : index_range(offsets))
634  {
635  offsets[i] = globalsize;
636  globalsize += counts[i];
637  }
638 
639  // Finally populate the primary vector
640  local_to_global_all.resize(globalsize, FeatureFloodCount::invalid_size_t);
641  for (const auto & feature : _feature_sets)
642  {
643  // Get the local indices from the feature and build a map
644  for (const auto & local_index_pair : feature._orig_ids)
645  {
646  auto rank = local_index_pair.first;
647  mooseAssert(rank < _n_procs, rank << ", " << _n_procs);
648 
649  auto local_index = local_index_pair.second;
650  auto stacked_local_index = offsets[rank] + local_index;
651 
652  mooseAssert(stacked_local_index < globalsize,
653  "Global index: " << stacked_local_index << " is out of range");
654  local_to_global_all[stacked_local_index] = feature._id;
655  }
656  }
657 }
658 
659 void
661 {
662  _feature_id_to_local_index.assign(max_id + 1, invalid_size_t);
663  for (const auto feature_index : index_range(_feature_sets))
664  {
665  if (_feature_sets[feature_index]._status != Status::INACTIVE)
666  {
667  mooseAssert(_feature_sets[feature_index]._id <= max_id,
668  "Feature ID out of range(" << _feature_sets[feature_index]._id << ')');
669  _feature_id_to_local_index[_feature_sets[feature_index]._id] = feature_index;
670  }
671  }
672 }
673 
674 void
676 {
677  TIME_SECTION("finalize", 3, "Finalizing Feature Identification");
678 
679  // Gather all information on processor zero and merge
681 
682  // Sort and label the features
683  if (_is_primary)
684  sortAndLabel();
685 
686  // Send out the local to global mappings
688 
689  // Populate _feature_maps and _var_index_maps
690  updateFieldInfo();
691 }
692 
693 const std::vector<unsigned int> &
695 {
696  mooseDoOnce(if (!_compute_var_to_feature_map) mooseError(
697  "Please set \"compute_var_to_feature_map = true\" to use this interface method"));
698 
699  const auto pos = _entity_var_to_features.find(elem_id);
700  if (pos != _entity_var_to_features.end())
701  {
702  mooseAssert(pos->second.size() == _n_vars, "Variable to feature vector not sized properly");
703  return pos->second;
704  }
705  else
706  return _empty_var_to_features;
707 }
708 
709 void
711 {
712  // local to global map (one per processor)
713  std::vector<int> counts;
714  std::vector<std::size_t> local_to_global_all;
715  if (_is_primary)
716  buildLocalToGlobalIndices(local_to_global_all, counts);
717 
718  // Scatter local_to_global indices to all processors and store in class member variable
719  _communicator.scatter(local_to_global_all, counts, _local_to_global_feature_map);
720 
721  std::size_t largest_global_index = std::numeric_limits<std::size_t>::lowest();
722  if (!_is_primary)
723  {
725 
732  for (auto & list_ref : _partial_feature_sets)
733  {
734  for (auto & feature : list_ref)
735  {
736  mooseAssert(feature._orig_ids.size() == 1, "feature._orig_ids length doesn't make sense");
737 
738  auto global_index = FeatureFloodCount::invalid_size_t;
739  auto local_index = feature._orig_ids.begin()->second;
740 
741  if (local_index < _local_to_global_feature_map.size())
742  global_index = _local_to_global_feature_map[local_index];
743 
744  if (global_index != FeatureFloodCount::invalid_size_t)
745  {
746  if (global_index > largest_global_index)
747  largest_global_index = global_index;
748 
749  // Set the correct global index
750  feature._id = global_index;
751 
759  feature._status &= ~Status::INACTIVE;
760 
761  // Move the feature into the correct place
762  _feature_sets[local_index] = std::move(feature);
763  }
764  }
765  }
766  }
767  else
768  {
769  for (auto global_index : local_to_global_all)
770  if (global_index != FeatureFloodCount::invalid_size_t && global_index > largest_global_index)
771  largest_global_index = global_index;
772  }
773 
774  // communicate the boundary intersection state
775  std::vector<std::pair<unsigned int, int>> intersection_state;
776  for (auto & feature : _feature_sets)
777  intersection_state.emplace_back(feature._id, static_cast<int>(feature._boundary_intersection));
778 
779  // gather on root
780  _communicator.gather(0, intersection_state);
781 
782  // consolidate
783  std::map<unsigned int, int> consolidated_intersection_state;
784  if (_is_primary)
785  for (const auto & [id, state] : intersection_state)
786  consolidated_intersection_state[id] |= state;
787 
788  // broadcast result
789  _communicator.broadcast(consolidated_intersection_state, 0);
790 
791  // apply broadcast changes
792  for (auto & feature : _feature_sets)
793  feature._boundary_intersection |=
794  static_cast<BoundaryIntersection>(consolidated_intersection_state[feature._id]);
795 
796  buildFeatureIdToLocalIndices(largest_global_index);
797 }
798 
799 Real
801 {
802  return static_cast<Real>(_feature_count);
803 }
804 
805 std::size_t
807 {
808  // Note: This value is parallel consistent, see FeatureFloodCount::communicateAndMerge()
809  return _feature_count;
810 }
811 
812 std::size_t
814 {
820  return _feature_count;
821 }
822 
823 unsigned int
824 FeatureFloodCount::getFeatureVar(unsigned int feature_id) const
825 {
826  // Some processors don't contain the largest feature id, in that case we just return invalid_id
827  if (feature_id >= _feature_id_to_local_index.size())
828  return invalid_id;
829 
830  auto local_index = _feature_id_to_local_index[feature_id];
831  if (local_index != invalid_size_t)
832  {
833  mooseAssert(local_index < _feature_sets.size(), "local_index out of bounds");
834  return _feature_sets[local_index]._status != Status::INACTIVE
835  ? _feature_sets[local_index]._var_index
836  : invalid_id;
837  }
838 
839  return invalid_id;
840 }
841 
842 bool
844 {
845  // Some processors don't contain the largest feature id, in that case we just return invalid_id
846  if (feature_id >= _feature_id_to_local_index.size())
847  return false;
848 
849  auto local_index = _feature_id_to_local_index[feature_id];
850 
851  if (local_index != invalid_size_t)
852  {
853  mooseAssert(local_index < _feature_sets.size(), "local_index out of bounds");
854  return _feature_sets[local_index]._status != Status::INACTIVE
855  ? _feature_sets[local_index]._boundary_intersection != BoundaryIntersection::NONE
856  : false;
857  }
858 
859  return false;
860 }
861 
862 bool
864 {
865  // Some processors don't contain the largest feature id, in that case we just return invalid_id
866  if (feature_id >= _feature_id_to_local_index.size())
867  return false;
868 
869  auto local_index = _feature_id_to_local_index[feature_id];
870 
871  if (local_index != invalid_size_t)
872  {
873  mooseAssert(local_index < _feature_sets.size(), "local_index out of bounds");
874  return _feature_sets[local_index]._status != Status::INACTIVE
875  ? ((_feature_sets[local_index]._boundary_intersection &
878  : false;
879  }
880 
881  return false;
882 }
883 
884 bool
885 FeatureFloodCount::isFeaturePercolated(unsigned int feature_id) const
886 {
887  // Some processors don't contain the largest feature id, in that case we just return invalid_id
888  if (feature_id >= _feature_id_to_local_index.size())
889  return false;
890 
891  auto local_index = _feature_id_to_local_index[feature_id];
892 
893  if (local_index != invalid_size_t)
894  {
895  mooseAssert(local_index < _feature_sets.size(), "local_index out of bounds");
896  bool primary = ((_feature_sets[local_index]._boundary_intersection &
899  bool secondary = ((_feature_sets[local_index]._boundary_intersection &
902  return _feature_sets[local_index]._status != Status::INACTIVE ? (primary && secondary) : false;
903  }
904 
905  return false;
906 }
907 
908 Point
909 FeatureFloodCount::featureCentroid(unsigned int feature_id) const
910 {
911  if (feature_id >= _feature_id_to_local_index.size())
912  return invalid_id;
913 
914  auto local_index = _feature_id_to_local_index[feature_id];
915 
916  Real invalid_coord = std::numeric_limits<Real>::max();
917  Point p(invalid_coord, invalid_coord, invalid_coord);
918  if (local_index != invalid_size_t)
919  {
920  mooseAssert(local_index < _feature_sets.size(), "local_index out of bounds");
921  p = _feature_sets[local_index]._centroid;
922  }
923  return p;
924 }
925 
926 Real
928  FieldType field_type,
929  std::size_t var_index) const
930 {
931  auto use_default = false;
932  if (var_index == invalid_size_t)
933  {
934  use_default = true;
935  var_index = 0;
936  }
937 
938  mooseAssert(var_index < _maps_size, "Index out of range");
939 
940  switch (field_type)
941  {
943  {
944  const auto entity_it = _feature_maps[var_index].find(entity_id);
945 
946  if (entity_it != _feature_maps[var_index].end())
947  return entity_it->second; // + _region_offsets[var_index];
948  else
949  return -1;
950  }
951 
953  {
954  mooseAssert(
956  "\"enable_var_coloring\" must be set to true to pull back the VARIABLE_COLORING field");
957 
958  const auto entity_it = _var_index_maps[var_index].find(entity_id);
959 
960  if (entity_it != _var_index_maps[var_index].end())
961  return entity_it->second;
962  else
963  return -1;
964  }
965 
967  {
968  const auto entity_it = _ghosted_entity_ids.find(entity_id);
969 
970  if (entity_it != _ghosted_entity_ids.end())
971  return entity_it->second;
972  else
973  return -1;
974  }
975 
976  case FieldType::HALOS:
977  {
978  if (!use_default)
979  {
980  const auto entity_it = _halo_ids[var_index].find(entity_id);
981  if (entity_it != _halo_ids[var_index].end())
982  return entity_it->second;
983  }
984  else
985  {
986  // Showing halos in reverse order for backwards compatibility
987  for (auto map_num = _maps_size;
988  map_num-- /* don't compare greater than zero for unsigned */;)
989  {
990  const auto entity_it = _halo_ids[map_num].find(entity_id);
991 
992  if (entity_it != _halo_ids[map_num].end())
993  return entity_it->second;
994  }
995  }
996  return -1;
997  }
998 
999  case FieldType::CENTROID:
1000  {
1001  if (_periodic_node_map.size())
1002  mooseDoOnce(mooseWarning(
1003  "Centroids are not correct when using periodic boundaries, contact the MOOSE team"));
1004 
1005  // If this element contains the centroid of one of features, return one
1006  const auto * elem_ptr = _mesh.elemPtr(entity_id);
1007 
1008  for (const auto & feature : _feature_sets)
1009  {
1010  if (feature._status == Status::INACTIVE)
1011  continue;
1012 
1013  if (elem_ptr->contains_point(feature._centroid))
1014  return 1;
1015  }
1016 
1017  return 0;
1018  }
1019 
1021  {
1022  auto ids = getVarToFeatureVector(entity_id);
1023  if (ids.size() != 0)
1025  return 0;
1026  }
1027 
1028  default:
1029  return 0;
1030  }
1031 }
1032 
1033 void
1035 {
1036  TIME_SECTION("prepareDataForTransfer", 3, "Preparing Data For Transfer");
1037 
1038  MeshBase & mesh = _mesh.getMesh();
1039 
1040  FeatureData::container_type local_ids_no_ghost, set_difference;
1041 
1042  for (auto & list_ref : _partial_feature_sets)
1043  {
1044  for (auto & feature : list_ref)
1045  {
1046  // See if the feature intersects a boundary or perhaps one of the percolation boundaries.
1047  updateBoundaryIntersections(feature);
1048 
1049  // Periodic node ids
1050  appendPeriodicNeighborNodes(feature);
1051 
1057  FeatureFloodCount::sort(feature._ghosted_ids);
1058  FeatureFloodCount::sort(feature._local_ids);
1059  FeatureFloodCount::sort(feature._halo_ids);
1060  FeatureFloodCount::sort(feature._disjoint_halo_ids);
1061  FeatureFloodCount::sort(feature._periodic_nodes);
1062 
1063  // Now extend the bounding box by the halo region
1064  if (_is_elemental)
1065  feature.updateBBoxExtremes(mesh);
1066  else
1067  {
1068  for (auto & halo_id : feature._halo_ids)
1069  updateBBoxExtremesHelper(feature._bboxes[0], mesh.point(halo_id));
1070  }
1071 
1072  mooseAssert(!feature._local_ids.empty(), "local entity ids cannot be empty");
1073 
1078  feature._min_entity_id = *feature._local_ids.begin();
1079  }
1080  }
1081 }
1082 
1083 void
1084 FeatureFloodCount::serialize(std::string & serialized_buffer, unsigned int var_num)
1085 {
1086  // stream for serializing the _partial_feature_sets data structure to a byte stream
1087  std::ostringstream oss;
1088 
1089  mooseAssert(var_num == invalid_id || var_num < _partial_feature_sets.size(),
1090  "var_num out of range");
1091 
1092  // Serialize everything
1093  if (var_num == invalid_id)
1094  dataStore(oss, _partial_feature_sets, this);
1095  else
1096  dataStore(oss, _partial_feature_sets[var_num], this);
1097 
1098  // Populate the passed in string pointer with the string stream's buffer contents
1099  serialized_buffer.assign(oss.str());
1100 }
1101 
1102 void
1103 FeatureFloodCount::deserialize(std::vector<std::string> & serialized_buffers, unsigned int var_num)
1104 {
1105  // The input string stream used for deserialization
1106  std::istringstream iss;
1107 
1108  auto rank = processor_id();
1109 
1110  for (const auto proc_id : index_range(serialized_buffers))
1111  {
1123  if (var_num == invalid_id && proc_id == rank)
1124  continue;
1125 
1126  iss.str(serialized_buffers[proc_id]); // populate the stream with a new buffer
1127  iss.clear(); // reset the string stream state
1128 
1129  // Load the gathered data into the data structure.
1130  if (var_num == invalid_id)
1131  dataLoad(iss, _partial_feature_sets, this);
1132  else
1133  dataLoad(iss, _partial_feature_sets[var_num], this);
1134  }
1135 }
1136 
1137 void
1139 {
1140  TIME_SECTION("mergeSets", 3, "Merging Sets");
1141 
1142  // When working with _distribute_merge_work all of the maps will be empty except for one
1143  for (const auto map_num : make_range(_maps_size))
1144  {
1145  for (auto it1 = _partial_feature_sets[map_num].begin();
1146  it1 != _partial_feature_sets[map_num].end();
1147  /* No increment on it1 */)
1148  {
1149  bool merge_occured = false;
1150  for (auto it2 = _partial_feature_sets[map_num].begin();
1151  it2 != _partial_feature_sets[map_num].end();
1152  ++it2)
1153  {
1154  if (it1 != it2 && areFeaturesMergeable(*it1, *it2))
1155  {
1156  it2->merge(std::move(*it1));
1157 
1162  _partial_feature_sets[map_num].emplace_back(std::move(*it2));
1163 
1173  _partial_feature_sets[map_num].erase(it2);
1174  it1 = _partial_feature_sets[map_num].erase(it1); // it1 is incremented here!
1175 
1176  // A merge occurred, this is used to determine whether or not we increment the outer
1177  // iterator
1178  merge_occured = true;
1179 
1180  // We need to start the list comparison over for the new it1 so break here
1181  break;
1182  }
1183  } // it2 loop
1184 
1185  if (!merge_occured) // No merges so we need to manually increment the outer iterator
1186  ++it1;
1187 
1188  } // it1 loop
1189  } // map loop
1190 }
1191 
1192 void
1193 FeatureFloodCount::consolidateMergedFeatures(std::vector<std::list<FeatureData>> * saved_data)
1194 {
1195  TIME_SECTION("consolidateMergedFeatures", 3, "Consolidating Merged Features");
1196 
1205  mooseAssert(_is_primary,
1206  "cosolidateMergedFeatures() may only be called on the primary processor");
1207  mooseAssert(saved_data == nullptr || saved_data->size() == _partial_feature_sets.size(),
1208  "Data structure size mismatch");
1209 
1210  // Offset where the current set of features with the same variable id starts in the flat vector
1211  unsigned int feature_offset = 0;
1212  // Set the member feature count to zero and start counting the actual features
1213  _feature_count = 0;
1214  for (const auto map_num : index_range(_partial_feature_sets))
1215  {
1216  for (auto & feature : _partial_feature_sets[map_num])
1217  {
1218  if (saved_data)
1219  {
1220  for (auto it = (*saved_data)[map_num].begin(); it != (*saved_data)[map_num].end();
1221  /* no increment */)
1222  {
1223  if (feature.canConsolidate(*it))
1224  {
1225  feature.consolidate(std::move(*it));
1226  it = (*saved_data)[map_num].erase(it); // increment
1227  }
1228  else
1229  ++it;
1230  }
1231  }
1232 
1233  // If after merging we still have an inactive feature, discard it
1234  if (feature._status == Status::CLEAR)
1235  {
1236  // First we need to calculate the centroid now that we are doing merging all partial
1237  // features
1238  if (feature._vol_count != 0)
1239  feature._centroid /= feature._vol_count;
1240 
1241  _feature_sets.emplace_back(std::move(feature));
1242  ++_feature_count;
1243  }
1244  }
1245 
1246  // Record the feature numbers just for the current map
1247  _feature_counts_per_map[map_num] = _feature_count - feature_offset;
1248 
1249  // Now update the running feature count so we can calculate the next map's contribution
1250  feature_offset = _feature_count;
1251 
1252  // Clean up the "moved" objects
1253  _partial_feature_sets[map_num].clear();
1254  if (saved_data)
1255  (*saved_data)[map_num].clear();
1256  }
1257 
1258  // We may have resided our data structures for the communicateAndMerge step. We'll restore the
1259  // original size here just in case we need to loop over the assumed size (i.e. _maps_size)
1260  // elsewhere in this or derived objects.
1261  if (_partial_feature_sets.size() != _maps_size)
1262  {
1264 
1267  }
1268 
1273 }
1274 
1275 bool
1277 {
1278  return f1.mergeable(f2);
1279 }
1280 
1281 void
1283 {
1284  for (const auto i : index_range(_feature_sets))
1285  {
1286  auto & feature = _feature_sets[i];
1287 
1288  // If the developer has requested _condense_map_info we'll make sure we only update the zeroth
1289  // map
1290  auto map_index = (_single_map_mode || _condense_map_info) ? decltype(feature._var_index)(0)
1291  : feature._var_index;
1292 
1293  // Loop over the entity ids of this feature and update our local map
1294  for (auto entity : feature._local_ids)
1295  {
1296  _feature_maps[map_index][entity] = static_cast<int>(feature._id);
1297 
1298  if (_var_index_mode)
1299  _var_index_maps[map_index][entity] = feature._var_index;
1300 
1301  // Fill in the data structure that keeps track of all features per elem
1303  {
1304  auto insert_pair = moose_try_emplace(
1305  _entity_var_to_features, entity, std::vector<unsigned int>(_n_vars, invalid_id));
1306  auto & vec_ref = insert_pair.first->second;
1307  vec_ref[feature._var_index] = feature._id;
1308  }
1309  }
1310 
1311  if (_compute_halo_maps)
1312  // Loop over the halo ids to update cells with halo information
1313  for (auto entity : feature._halo_ids)
1314  _halo_ids[map_index][entity] = static_cast<int>(feature._id);
1315 
1316  // Loop over the ghosted ids to update cells with ghost information
1317  for (auto entity : feature._ghosted_ids)
1318  _ghosted_entity_ids[entity] = 1;
1319 
1320  // TODO: Fixme
1321  if (!_global_numbering)
1322  mooseError("Local numbering currently disabled");
1323  }
1324 }
1325 
1326 bool
1327 FeatureFloodCount::flood(const DofObject * dof_object, std::size_t current_index)
1328 
1329 {
1330  // if (dof_object == nullptr || dof_object == libMesh::remote_elem)
1331  // return false;
1332  mooseAssert(dof_object, "DOF object is nullptr");
1333  mooseAssert(_entity_queue.empty(), "Entity queue is not empty when starting a feature");
1334 
1335  // Kick off the exploration of a new feature
1336  _entity_queue.push_front(dof_object);
1337 
1338  bool return_value = false;
1339  FeatureData * feature = nullptr;
1340  while (!_entity_queue.empty())
1341  {
1342  const DofObject * curr_dof_object = _entity_queue.back();
1343  const Elem * elem = _is_elemental ? static_cast<const Elem *>(curr_dof_object) : nullptr;
1344  _entity_queue.pop_back();
1345 
1346  // Retrieve the id of the current entity
1347  auto entity_id = curr_dof_object->id();
1348 
1349  // Has this entity already been marked? - if so move along
1350  if (current_index != invalid_size_t &&
1351  _entities_visited[current_index].find(entity_id) != _entities_visited[current_index].end())
1352  continue;
1353 
1354  // Are we outside of the range we should be working in?
1355  if (_is_elemental && !_dof_map.is_evaluable(*elem))
1356  continue;
1357 
1358  // See if the current entity either starts a new feature or continues an existing feature
1359  auto new_id = invalid_id; // Writable reference to hold an optional id;
1360  Status status =
1361  Status::INACTIVE; // Status is inactive until we find an entity above the starting threshold
1362 
1363  // Make sure that the Assembly object has the right element and subdomain information set
1364  // since we are moving through the mesh in a manual fashion.
1365  if (_is_elemental)
1367 
1368  if (!isNewFeatureOrConnectedRegion(curr_dof_object, current_index, feature, status, new_id))
1369  {
1370  // If we have an active feature, we just found a halo entity
1371  if (feature)
1372  feature->_halo_ids.insert(feature->_halo_ids.end(), entity_id);
1373  continue;
1374  }
1375 
1376  mooseAssert(current_index != invalid_size_t, "current_index is invalid");
1377 
1386  return_value = true;
1387  _entities_visited[current_index].insert(entity_id);
1388 
1389  auto map_num = _single_map_mode ? decltype(current_index)(0) : current_index;
1390 
1391  // New Feature (we need to create it and add it to our data structure)
1392  if (!feature)
1393  {
1394  _partial_feature_sets[map_num].emplace_back(
1395  current_index, _feature_count++, processor_id(), status);
1396 
1397  // Get a handle to the feature we will update (always the last feature in the data structure)
1398  feature = &_partial_feature_sets[map_num].back();
1399 
1400  // If new_id is valid, we'll set it in the feature here.
1401  if (new_id != invalid_id)
1402  feature->_id = new_id;
1403  }
1404 
1405  // Insert the current entity into the local ids data structure
1406  feature->_local_ids.insert(feature->_local_ids.end(), entity_id);
1407 
1413  if (_is_elemental && processor_id() == curr_dof_object->processor_id())
1414  {
1415  // Keep track of how many elements participate in the centroid averaging
1416  feature->_vol_count++;
1417 
1418  // Sum the centroid values for now, we'll average them later
1419  feature->_centroid += elem->vertex_average();
1420 
1421  // // Does the volume intersect the boundary?
1422  // if (_all_boundary_entity_ids.find(elem->id()) != _all_boundary_entity_ids.end())
1423  // feature->_intersects_boundary = true;
1424  }
1425 
1426  if (_is_elemental)
1428  feature,
1429  /*expand_halos_only =*/false,
1430  /*disjoint_only =*/false);
1431  else
1432  visitNodalNeighbors(static_cast<const Node *>(curr_dof_object),
1433  feature,
1434  /*expand_halos_only =*/false);
1435  }
1436 
1437  return return_value;
1438 }
1439 
1440 Real FeatureFloodCount::getThreshold(std::size_t /*current_index*/) const
1441 {
1442  return _step_threshold;
1443 }
1444 
1445 Real FeatureFloodCount::getConnectingThreshold(std::size_t /*current_index*/) const
1446 {
1448 }
1449 
1450 bool
1451 FeatureFloodCount::compareValueWithThreshold(Real entity_value, Real threshold) const
1452 {
1453  return ((_use_less_than_threshold_comparison && (entity_value >= threshold)) ||
1454  (!_use_less_than_threshold_comparison && (entity_value <= threshold)));
1455 }
1456 
1457 bool
1459  std::size_t & current_index,
1460  FeatureData *& feature,
1461  Status & status,
1462  unsigned int & /*new_id*/)
1463 {
1464  // Get the value of the current variable for the current entity
1465  Real entity_value;
1466  if (_is_elemental)
1467  {
1468  const Elem * elem = static_cast<const Elem *>(dof_object);
1469  std::vector<Point> centroid(1, elem->vertex_average());
1470  _subproblem.reinitElemPhys(elem, centroid, 0);
1471  entity_value = _vars[current_index]->sln()[0];
1472  }
1473  else
1474  entity_value = _vars[current_index]->getNodalValue(*static_cast<const Node *>(dof_object));
1475 
1476  // If the value compares against our starting threshold, this is definitely part of a feature
1477  // we'll keep
1478  if (compareValueWithThreshold(entity_value, getThreshold(current_index)))
1479  {
1480  Status * status_ptr = &status;
1481 
1482  if (feature)
1483  status_ptr = &feature->_status;
1484 
1485  // Update an existing feature's status or clear the flag on the passed in status
1486  *status_ptr &= ~Status::INACTIVE;
1487  return true;
1488  }
1489 
1494  return compareValueWithThreshold(entity_value, getConnectingThreshold(current_index));
1495 }
1496 
1497 void
1499 {
1500  const auto & node_to_elem_map = _mesh.nodeToActiveSemilocalElemMap();
1501  FeatureData::container_type expanded_local_ids;
1502  auto my_processor_id = processor_id();
1503 
1510  for (auto & list_ref : _partial_feature_sets)
1511  {
1512  for (auto & feature : list_ref)
1513  {
1514  expanded_local_ids.clear();
1515 
1516  for (auto entity : feature._local_ids)
1517  {
1518  const Elem * elem = _mesh.elemPtr(entity);
1519  mooseAssert(elem, "elem pointer is NULL");
1520 
1521  // Get the nodes on a current element so that we can add in point neighbors
1522  auto n_nodes = elem->n_vertices();
1523  for (const auto i : make_range(n_nodes))
1524  {
1525  const Node * current_node = elem->node_ptr(i);
1526 
1527  auto elem_vector_it = node_to_elem_map.find(current_node->id());
1528  if (elem_vector_it == node_to_elem_map.end())
1529  mooseError("Error in node to elem map");
1530 
1531  const auto & elem_vector = elem_vector_it->second;
1532 
1533  std::copy(elem_vector.begin(),
1534  elem_vector.end(),
1535  std::insert_iterator<FeatureData::container_type>(expanded_local_ids,
1536  expanded_local_ids.end()));
1537 
1538  // Now see which elements need to go into the ghosted set
1539  for (auto entity : elem_vector)
1540  {
1541  const Elem * neighbor = _mesh.elemPtr(entity);
1542  mooseAssert(neighbor, "neighbor pointer is NULL");
1543 
1544  if (neighbor->processor_id() != my_processor_id)
1545  feature._ghosted_ids.insert(feature._ghosted_ids.end(), elem->id());
1546  }
1547  }
1548  }
1549 
1550  // Replace the existing local ids with the expanded local ids
1551  feature._local_ids.swap(expanded_local_ids);
1552 
1553  // Copy the expanded local_ids into the halo_ids container
1554  feature._halo_ids = feature._local_ids;
1555  }
1556  }
1557 }
1558 
1559 void
1560 FeatureFloodCount::expandEdgeHalos(unsigned int num_layers_to_expand)
1561 {
1562  if (num_layers_to_expand == 0)
1563  return;
1564 
1565  TIME_SECTION("expandEdgeHalos", 3, "Expanding Edge Halos");
1566 
1567  for (auto & list_ref : _partial_feature_sets)
1568  {
1569  for (auto & feature : list_ref)
1570  {
1571  for (unsigned short halo_level = 0; halo_level < num_layers_to_expand; ++halo_level)
1572  {
1577  FeatureData::container_type orig_halo_ids(feature._halo_ids);
1578  for (auto entity : orig_halo_ids)
1579  {
1580  if (_is_elemental)
1582  &feature,
1583  /*expand_halos_only =*/true,
1584  /*disjoint_only =*/false);
1585  else
1587  &feature,
1588  /*expand_halos_only =*/true);
1589  }
1590 
1595  FeatureData::container_type disjoint_orig_halo_ids(feature._disjoint_halo_ids);
1596  for (auto entity : disjoint_orig_halo_ids)
1597  {
1598  if (_is_elemental)
1600 
1601  &feature,
1602  /*expand_halos_only =*/true,
1603  /*disjoint_only =*/true);
1604  else
1606 
1607  &feature,
1608  /*expand_halos_only =*/true);
1609  }
1610  }
1611  }
1612  }
1613 }
1614 
1615 void
1617  FeatureData * feature,
1618  bool expand_halos_only,
1619  bool disjoint_only)
1620 {
1621  mooseAssert(elem, "Elem is NULL");
1622 
1623  std::vector<const Elem *> all_active_neighbors;
1624  MeshBase & mesh = _mesh.getMesh();
1625 
1626  // Loop over all neighbors (at the the same level as the current element)
1627  for (const auto i : make_range(elem->n_neighbors()))
1628  {
1629  const Elem * neighbor_ancestor = nullptr;
1630  bool topological_neighbor = false;
1631 
1636  neighbor_ancestor = elem->neighbor_ptr(i);
1637 
1638  if (neighbor_ancestor)
1639  {
1646  if (neighbor_ancestor->is_remote())
1647  continue;
1648 
1649  neighbor_ancestor->active_family_tree_by_neighbor(all_active_neighbors, elem, false);
1650  }
1651  else
1652  {
1653  neighbor_ancestor = elem->topological_neighbor(i, mesh, *_point_locator, _pbs);
1654 
1662  if (neighbor_ancestor)
1663  {
1670  if (neighbor_ancestor->is_remote())
1671  continue;
1672 
1673  neighbor_ancestor->active_family_tree_by_topological_neighbor(
1674  all_active_neighbors, elem, mesh, *_point_locator, _pbs, false);
1675 
1676  topological_neighbor = true;
1677  }
1678  else
1679  {
1685  updateBBoxExtremesHelper(feature->_bboxes[0], *elem);
1686  }
1687  }
1688 
1689  visitNeighborsHelper(elem,
1690  all_active_neighbors,
1691  feature,
1692  expand_halos_only,
1693  topological_neighbor,
1694  disjoint_only);
1695 
1696  all_active_neighbors.clear();
1697  }
1698 }
1699 
1700 void
1702  FeatureData * feature,
1703  bool expand_halos_only)
1704 {
1705  mooseAssert(node, "Node is NULL");
1706 
1707  std::vector<const Node *> all_active_neighbors;
1708  MeshTools::find_nodal_neighbors(_mesh.getMesh(), *node, _nodes_to_elem_map, all_active_neighbors);
1709 
1710  visitNeighborsHelper(node, all_active_neighbors, feature, expand_halos_only, false, false);
1711 }
1712 
1713 template <typename T>
1714 void
1716  std::vector<const T *> neighbor_entities,
1717  FeatureData * feature,
1718  bool expand_halos_only,
1719  bool topological_neighbor,
1720  bool disjoint_only)
1721 {
1722  // Loop over all active element neighbors
1723  for (const auto neighbor : neighbor_entities)
1724  {
1725  if (neighbor && (!_is_boundary_restricted || isBoundaryEntity(neighbor)))
1726  {
1727  if (expand_halos_only)
1728  {
1729  auto entity_id = neighbor->id();
1730 
1731  if (topological_neighbor || disjoint_only)
1732  feature->_disjoint_halo_ids.insert(feature->_disjoint_halo_ids.end(), entity_id);
1733  else if (!FeatureFloodCount::contains(feature->_local_ids, entity_id))
1734  feature->_halo_ids.insert(feature->_halo_ids.end(), entity_id);
1735  }
1736  else
1737  {
1738  auto my_processor_id = processor_id();
1739 
1740  if (!topological_neighbor && neighbor->processor_id() != my_processor_id)
1741  feature->_ghosted_ids.insert(feature->_ghosted_ids.end(), curr_entity->id());
1742 
1752  if (curr_entity->processor_id() == my_processor_id ||
1753  neighbor->processor_id() == my_processor_id)
1754  {
1761  if (topological_neighbor || disjoint_only)
1762  feature->_disjoint_halo_ids.insert(feature->_disjoint_halo_ids.end(), neighbor->id());
1763  else
1764  _entity_queue.push_front(neighbor);
1765  }
1766  }
1767  }
1768  }
1769 }
1770 
1771 void
1773 {
1774  if (_is_elemental)
1775  {
1776  for (auto entity : feature._local_ids)
1777  {
1778  // See if this feature is on a boundary if we haven't already figured that out
1781  {
1782  Elem * elem = _mesh.elemPtr(entity);
1783  if (elem && elem->on_boundary())
1785  }
1786 
1787  // Now see if the feature touches the primary and/or secondary boundary IDs if we haven't
1788  // figured that out already
1791  {
1792  for (auto primary_id : _primary_perc_bnds)
1793  if (_mesh.isBoundaryElem(entity, primary_id))
1795  }
1796 
1799  {
1800  for (auto secondary_id : _secondary_perc_bnds)
1801  if (_mesh.isBoundaryElem(entity, secondary_id))
1803  }
1804 
1805  // See if the feature contacts any of the user-specified boundaries if we haven't
1806  // done so already
1809  {
1810  for (auto specified_id : _specified_bnds)
1811  if (_mesh.isBoundaryElem(entity, specified_id))
1813  }
1814  }
1815  }
1816 }
1817 
1818 void
1820 {
1821  if (_is_elemental)
1822  {
1823  for (auto entity : feature._local_ids)
1824  {
1825  Elem * elem = _mesh.elemPtr(entity);
1826 
1827  for (const auto node_n : make_range(elem->n_nodes()))
1828  {
1829  auto iters = _periodic_node_map.equal_range(elem->node_id(node_n));
1830 
1831  for (auto it = iters.first; it != iters.second; ++it)
1832  {
1833  feature._periodic_nodes.insert(feature._periodic_nodes.end(), it->first);
1834  feature._periodic_nodes.insert(feature._periodic_nodes.end(), it->second);
1835  }
1836  }
1837  }
1838  }
1839  else
1840  {
1841  for (auto entity : feature._local_ids)
1842  {
1843  auto iters = _periodic_node_map.equal_range(entity);
1844 
1845  for (auto it = iters.first; it != iters.second; ++it)
1846  {
1847  feature._periodic_nodes.insert(feature._periodic_nodes.end(), it->first);
1848  feature._periodic_nodes.insert(feature._periodic_nodes.end(), it->second);
1849  }
1850  }
1851  }
1852 
1853  // TODO: Remove duplicates
1854 }
1855 
1856 template <typename T>
1857 bool
1859 {
1860  mooseAssert(_bnd_elem_range, "Boundary Element Range is nullptr");
1861 
1862  if (entity)
1863  for (const auto & belem : *_bnd_elem_range)
1864  // Only works for Elements
1865  if (belem->_elem->id() == entity->id() && hasBoundary(belem->_bnd_id))
1866  return true;
1867 
1868  return false;
1869 }
1870 
1871 void
1873 {
1874  // First update the primary bounding box (all topologically connected)
1875  for (auto & halo_id : _halo_ids)
1876  updateBBoxExtremesHelper(_bboxes[0], mesh.elem_ref(halo_id));
1877  for (auto & ghost_id : _ghosted_ids)
1878  updateBBoxExtremesHelper(_bboxes[0], mesh.elem_ref(ghost_id));
1879 
1880  FeatureData::container_type all_primary_region_ids;
1881  FeatureFloodCount::reserve(all_primary_region_ids, _local_ids.size() + _halo_ids.size());
1882  std::set_union(_local_ids.begin(),
1883  _local_ids.end(),
1884  _halo_ids.begin(),
1885  _halo_ids.end(),
1886  std::insert_iterator<FeatureData::container_type>(all_primary_region_ids,
1887  all_primary_region_ids.begin()));
1888 
1889  // Remove all of the IDs that are in the primary region
1890  std::list<dof_id_type> disjoint_elem_id_list;
1891  std::set_difference(_disjoint_halo_ids.begin(),
1892  _disjoint_halo_ids.end(),
1893  all_primary_region_ids.begin(),
1894  all_primary_region_ids.end(),
1895  std::insert_iterator<std::list<dof_id_type>>(disjoint_elem_id_list,
1896  disjoint_elem_id_list.begin()));
1897 
1904  constexpr auto TOLERANCE = libMesh::TOLERANCE * libMesh::TOLERANCE;
1905  for (auto elem_id : disjoint_elem_id_list)
1906  {
1907  BoundingBox elem_bbox;
1908  updateBBoxExtremesHelper(elem_bbox, *mesh.elem_ptr(elem_id));
1909 
1910  bool found_match = false;
1911  for (auto & bbox : _bboxes)
1912  {
1913  if (bbox.intersects(elem_bbox, TOLERANCE))
1914  {
1915  updateBBoxExtremes(bbox, elem_bbox);
1916  found_match = true;
1917  break;
1918  }
1919  }
1920 
1921  if (!found_match)
1922  _bboxes.push_back(elem_bbox);
1923  }
1924 
1925  mergeBBoxes(_bboxes, false);
1926 
1931  FeatureData::container_type set_union;
1932  FeatureFloodCount::reserve(set_union, _halo_ids.size() + _disjoint_halo_ids.size());
1933  std::set_union(_halo_ids.begin(),
1934  _halo_ids.end(),
1935  _disjoint_halo_ids.begin(),
1936  _disjoint_halo_ids.end(),
1937  std::insert_iterator<FeatureData::container_type>(set_union, set_union.begin()));
1938  _halo_ids.swap(set_union);
1939  _disjoint_halo_ids.clear();
1940 }
1941 
1942 void
1943 FeatureFloodCount::FeatureData::updateBBoxExtremes(BoundingBox & bbox, const BoundingBox & rhs_bbox)
1944 {
1945  for (const auto i : make_range(Moose::dim))
1946  {
1947  bbox.min()(i) = std::min(bbox.min()(i), rhs_bbox.min()(i));
1948  bbox.max()(i) = std::max(bbox.max()(i), rhs_bbox.max()(i));
1949  }
1950 }
1951 
1952 bool
1954 {
1955  // See if any of the bounding boxes in either FeatureData object intersect
1956  for (const auto & bbox_lhs : _bboxes)
1957  for (const auto & bbox_rhs : rhs._bboxes)
1958  if (bbox_lhs.intersects(bbox_rhs, libMesh::TOLERANCE * libMesh::TOLERANCE))
1959  return true;
1960 
1961  return false;
1962 }
1963 
1964 bool
1966 {
1968 }
1969 
1970 bool
1972 {
1973  return MooseUtils::setsIntersect(_periodic_nodes, rhs._periodic_nodes);
1974 }
1975 
1976 bool
1978 {
1979  return MooseUtils::setsIntersect(_ghosted_ids, rhs._ghosted_ids);
1980 }
1981 
1982 bool
1984 {
1985  return (_var_index == rhs._var_index && // the sets have matching variable indices and
1986  ((boundingBoxesIntersect(rhs) && // (if the feature's bboxes intersect and
1987  ghostedIntersect(rhs)) // the ghosted entities also intersect)
1988  || // or
1989  periodicBoundariesIntersect(rhs))); // periodic node sets intersect)
1990 }
1991 
1992 bool
1994 {
1995  for (const auto & orig_id_pair1 : _orig_ids)
1996  for (const auto & orig_id_pair2 : rhs._orig_ids)
1997  if (orig_id_pair1 == orig_id_pair2)
1998  return true;
1999 
2000  return false;
2001 }
2002 
2003 void
2005 {
2006  mooseAssert(_var_index == rhs._var_index, "Mismatched variable index in merge");
2007  mooseAssert(_id == rhs._id, "Mismatched auxiliary id in merge");
2008 
2009  FeatureData::container_type set_union;
2010 
2011  FeatureFloodCount::reserve(set_union, _local_ids.size() + rhs._local_ids.size());
2012  std::set_union(_local_ids.begin(),
2013  _local_ids.end(),
2014  rhs._local_ids.begin(),
2015  rhs._local_ids.end(),
2016  std::insert_iterator<FeatureData::container_type>(set_union, set_union.begin()));
2017  _local_ids.swap(set_union);
2018 
2019  set_union.clear();
2020  FeatureFloodCount::reserve(set_union, _periodic_nodes.size() + rhs._periodic_nodes.size());
2021  std::set_union(_periodic_nodes.begin(),
2022  _periodic_nodes.end(),
2023  rhs._periodic_nodes.begin(),
2024  rhs._periodic_nodes.end(),
2025  std::insert_iterator<FeatureData::container_type>(set_union, set_union.begin()));
2026  _periodic_nodes.swap(set_union);
2027 
2028  set_union.clear();
2029  FeatureFloodCount::reserve(set_union, _ghosted_ids.size() + rhs._ghosted_ids.size());
2030  std::set_union(_ghosted_ids.begin(),
2031  _ghosted_ids.end(),
2032  rhs._ghosted_ids.begin(),
2033  rhs._ghosted_ids.end(),
2034  std::insert_iterator<FeatureData::container_type>(set_union, set_union.begin()));
2035 
2043  bool physical_intersection = (_ghosted_ids.size() + rhs._ghosted_ids.size() > set_union.size());
2044  _ghosted_ids.swap(set_union);
2045 
2050  _bboxes.reserve(_bboxes.size() + rhs._bboxes.size());
2051  std::copy(rhs._bboxes.begin(), rhs._bboxes.end(), std::back_inserter(_bboxes));
2052 
2053  mergeBBoxes(_bboxes, physical_intersection);
2054 
2055  set_union.clear();
2056  FeatureFloodCount::reserve(set_union, _disjoint_halo_ids.size() + rhs._disjoint_halo_ids.size());
2057  std::set_union(_disjoint_halo_ids.begin(),
2058  _disjoint_halo_ids.end(),
2059  rhs._disjoint_halo_ids.begin(),
2060  rhs._disjoint_halo_ids.end(),
2061  std::insert_iterator<FeatureData::container_type>(set_union, set_union.begin()));
2062  _disjoint_halo_ids.swap(set_union);
2063 
2064  set_union.clear();
2065  FeatureFloodCount::reserve(set_union, _halo_ids.size() + rhs._halo_ids.size());
2066  std::set_union(_halo_ids.begin(),
2067  _halo_ids.end(),
2068  rhs._halo_ids.begin(),
2069  rhs._halo_ids.end(),
2070  std::insert_iterator<FeatureData::container_type>(set_union, set_union.begin()));
2071  _halo_ids.swap(set_union);
2072 
2073  // Keep track of the original ids so we can notify other processors of the local to global mapping
2074  _orig_ids.splice(_orig_ids.end(), std::move(rhs._orig_ids));
2075 
2076  // Update the min feature id
2077  _min_entity_id = std::min(_min_entity_id, rhs._min_entity_id);
2078 
2084  mooseAssert((_status & Status::MARKED & Status::DIRTY) == Status::CLEAR,
2085  "Flags in invalid state");
2086 
2087  // Logical AND here to combine flags (INACTIVE & INACTIVE == INACTIVE, all other combos are CLEAR)
2088  _status &= rhs._status;
2089 
2090  // Logical OR here to make sure we maintain boundary intersection attribute when joining
2091  _boundary_intersection |= rhs._boundary_intersection;
2092 
2093  _vol_count += rhs._vol_count;
2094  _centroid += rhs._centroid;
2095 }
2096 
2097 void
2099 {
2100  mooseAssert(_var_index == rhs._var_index, "Mismatched variable index in merge");
2101  mooseAssert(_id == rhs._id, "Mismatched auxiliary id in merge");
2102 
2103  FeatureData::container_type set_union;
2104  FeatureFloodCount::reserve(_local_ids, _local_ids.size() + rhs._local_ids.size());
2105  std::set_union(_local_ids.begin(),
2106  _local_ids.end(),
2107  rhs._local_ids.begin(),
2108  rhs._local_ids.end(),
2109  std::insert_iterator<FeatureData::container_type>(set_union, set_union.begin()));
2110  _local_ids.swap(set_union);
2111 
2112  mooseAssert((_status & Status::MARKED & Status::DIRTY) == Status::CLEAR,
2113  "Flags in invalid state");
2114 }
2115 
2116 void
2118 {
2119  _local_ids.clear();
2120  _periodic_nodes.clear();
2121  _halo_ids.clear();
2122  _disjoint_halo_ids.clear();
2123  _ghosted_ids.clear();
2124  _bboxes.clear();
2125  _orig_ids.clear();
2126 }
2127 
2128 void
2129 FeatureFloodCount::FeatureData::mergeBBoxes(std::vector<BoundingBox> & bboxes,
2130  bool physical_intersection)
2131 {
2140  std::list<BoundingBox> box_list(bboxes.begin(), bboxes.end());
2141 
2142  auto box_expanded = false;
2143  for (auto it1 = box_list.begin(); it1 != box_list.end(); /* No increment on it1 */)
2144  {
2145  auto merge_occured = false;
2146  for (auto it2 = box_list.begin(); it2 != box_list.end(); ++it2)
2147  {
2148  if (it1 != it2 && it1->intersects(*it2, TOLERANCE))
2149  {
2150  updateBBoxExtremes(*it2, *it1);
2151  box_list.emplace_back(std::move(*it2));
2152 
2153  box_list.erase(it2);
2154  it1 = box_list.erase(it1); // it1 is incremented here!
2155 
2156  box_expanded = true;
2157  merge_occured = true;
2158 
2159  break;
2160  }
2161  }
2162 
2163  if (!merge_occured)
2164  ++it1;
2165  }
2166 
2170  bboxes.resize(box_list.size());
2171  std::copy(box_list.begin(), box_list.end(), bboxes.begin());
2172 
2173  // Error check
2174  if (physical_intersection && !box_expanded)
2175  {
2176  std::ostringstream oss;
2177  oss << "LHS BBoxes:\n";
2178  for (const auto i : index_range(_bboxes))
2179  oss << "Max: " << _bboxes[i].max() << " Min: " << _bboxes[i].min() << '\n';
2180 
2181  ::mooseError("No Bounding Boxes Expanded - This is a catastrophic error!\n", oss.str());
2182  }
2183 }
2184 
2185 std::ostream &
2186 operator<<(std::ostream & out, const FeatureFloodCount::FeatureData & feature)
2187 {
2188  static const bool debug = true;
2189 
2190  out << "Grain ID: ";
2191  if (feature._id != FeatureFloodCount::invalid_id)
2192  out << feature._id;
2193  else
2194  out << "invalid";
2195 
2196  if (debug)
2197  {
2198  out << "\nGhosted Entities: ";
2199  for (auto ghosted_id : feature._ghosted_ids)
2200  out << ghosted_id << " ";
2201 
2202  out << "\nLocal Entities: ";
2203  for (auto local_id : feature._local_ids)
2204  out << local_id << " ";
2205 
2206  out << "\nHalo Entities: ";
2207  for (auto halo_id : feature._halo_ids)
2208  out << halo_id << " ";
2209 
2210  out << "\nPeriodic Node IDs: ";
2211  for (auto periodic_node : feature._periodic_nodes)
2212  out << periodic_node << " ";
2213  }
2214 
2215  out << "\nBBoxes:";
2216  for (const auto & bbox : feature._bboxes)
2217  {
2218  out << "\nMax: " << bbox.max() << " Min: " << bbox.min();
2219  }
2220 
2221  out << "\nStatus: ";
2223  out << "CLEAR";
2224  if (static_cast<bool>(feature._status & FeatureFloodCount::Status::MARKED))
2225  out << " MARKED";
2226  if (static_cast<bool>(feature._status & FeatureFloodCount::Status::DIRTY))
2227  out << " DIRTY";
2228  if (static_cast<bool>(feature._status & FeatureFloodCount::Status::INACTIVE))
2229  out << " INACTIVE";
2230 
2231  if (debug)
2232  {
2233  out << "\nOrig IDs (rank, index): ";
2234  for (const auto & orig_pair : feature._orig_ids)
2235  out << '(' << orig_pair.first << ", " << orig_pair.second << ") ";
2236  out << "\nVar_index: " << feature._var_index;
2237  out << "\nMin Entity ID: " << feature._min_entity_id;
2238  }
2239  out << "\n" << std::endl;
2240 
2241  return out;
2242 }
2243 
2244 /*****************************************************************************************
2245  ******************************* Utility Routines ****************************************
2246  *****************************************************************************************
2247  */
2248 void
2249 updateBBoxExtremesHelper(BoundingBox & bbox, const Point & node)
2250 {
2251  for (const auto i : make_range(Moose::dim))
2252  {
2253  bbox.min()(i) = std::min(bbox.min()(i), node(i));
2254  bbox.max()(i) = std::max(bbox.max()(i), node(i));
2255  }
2256 }
2257 
2258 void
2259 updateBBoxExtremesHelper(BoundingBox & bbox, const Elem & elem)
2260 {
2261  for (const auto node_n : make_range(elem.n_nodes()))
2262  updateBBoxExtremesHelper(bbox, *(elem.node_ptr(node_n)));
2263 }
2264 
2265 // Constants
2266 const std::size_t FeatureFloodCount::invalid_size_t = std::numeric_limits<std::size_t>::max();
2267 const unsigned int FeatureFloodCount::invalid_id = std::numeric_limits<unsigned int>::max();
2269  std::numeric_limits<processor_id_type>::max();
virtual bnd_elem_iterator bndElemsEnd()
FeatureFloodCount(const InputParameters &parameters)
void gather_packed_range(const unsigned int root_id, Context *context, Iter range_begin, const Iter range_end, OutputIter out, std::size_t approx_buffer_size=1000000) const
std::vector< BoundaryID > _primary_perc_bnds
Point _centroid
The centroid of the feature (average of coordinates from entities participating in the volume calcula...
std::size_t getNumberActiveFeatures() const
Return the number of active features.
bool halosIntersect(const FeatureData &rhs) const
Determine if one of this FeaturesData&#39;s member sets intersects the other FeatureData&#39;s corresponding ...
std::ostream & operator<<(std::ostream &out, const FeatureFloodCount::FeatureData &feature)
std::multimap< dof_id_type, dof_id_type > _periodic_node_map
The data structure which is a list of nodes that are constrained to other nodes based on the imposed ...
void visitNodalNeighbors(const Node *node, FeatureData *feature, bool expand_halos_only)
These two routines are utility routines used by the flood routine and by derived classes for visiting...
std::vector< BoundingBox > _bboxes
The vector of bounding boxes completely enclosing this feature (multiple used with periodic constrain...
const PostprocessorValue & _element_average_value
Average value of the domain which can optionally be used to find features in a field.
void expandEdgeHalos(unsigned int num_layers_to_expand)
This method expands the existing halo set by some width determined by the passed in value...
std::vector< BoundaryID > _specified_bnds
std::vector< BoundaryID > _secondary_perc_bnds
void expandPointHalos()
This method takes all of the partial features and expands the local, ghosted, and halo sets around th...
container_type _halo_ids
Holds the ids surrounding the feature.
void scatter(const std::vector< T, A > &data, T &recv, const unsigned int root_id=0) const
const std::size_t _n_vars
static const std::size_t invalid_size_t
virtual void updateFieldInfo()
This method is used to populate any of the data structures used for storing field data (nodal or elem...
static InputParameters validParams()
virtual bool boundaryRestricted() const
std::vector< dof_id_type > container_type
The primary underlying container type used to hold the data in each FeatureData.
void appendPeriodicNeighborNodes(FeatureData &feature) const
This routine adds the periodic node information to our data structure prior to packing the data this ...
virtual Elem * elemPtr(const dof_id_type i)
void addParam(const std::string &name, const std::initializer_list< typename T::value_type > &value, const std::string &doc_string)
virtual const std::vector< unsigned int > & getVarToFeatureVector(dof_id_type elem_id) const
Returns a list of active unique feature ids for a particular element.
void serialize(std::string &serialized_buffer, unsigned int var_num=invalid_id)
This routines packs the _partial_feature_sets data into a structure suitable for parallel communicati...
void addPrivateParam(const std::string &name, const T &value)
const std::map< dof_id_type, std::vector< dof_id_type > > & nodeToActiveSemilocalElemMap()
std::vector< std::pair< R1, R2 > > get(const std::string &param1, const std::string &param2) const
static constexpr Real TOLERANCE
const bool _condense_map_info
void mooseInfo(Args &&... args) const
void gather(const unsigned int root_id, const T &send_data, std::vector< T, A > &recv) const
container_type _ghosted_ids
Holds the ghosted ids for a feature (the ids which will be used for stitching.
Status
This enumeration is used to indicate status of the grains in the _unique_grains data structure...
std::vector< std::set< dof_id_type > > _entities_visited
This variable keeps track of which nodes have been visited during execution.
registerMooseObject("PhaseFieldApp", FeatureFloodCount)
void updateBBoxExtremes(MeshBase &mesh)
Update the minimum and maximum coordinates of a bounding box given a Point, Elem or BBox parameter...
bool mergeable(const FeatureData &rhs) const
The routine called to see if two features are mergeable:
bool boundingBoxesIntersect(const FeatureData &rhs) const
Determines if any of this FeatureData&#39;s bounding boxes overlap with the other FeatureData&#39;s bounding ...
T & set(const std::string &name, bool quiet_mode=false)
static void sort(std::set< T > &)
virtual void finalize() override
bool setsIntersect(InputIterator first1, InputIterator last1, InputIterator first2, InputIterator last2)
container_type _disjoint_halo_ids
Holds halo ids that extend onto a non-topologically connected surface.
std::vector< FeatureData > & _feature_sets
The data structure used to hold the globally unique features.
MeshBase & mesh
bool isBoundaryEntity(const T *entity) const
Returns a Boolean indicating whether the entity is on one of the desired boundaries.
const bool _is_primary
Convenience variable for testing primary rank.
static constexpr std::size_t dim
virtual bnd_elem_iterator bndElemsBegin()
std::map< dof_id_type, std::vector< unsigned int > > _entity_var_to_features
virtual std::size_t getTotalFeatureCount() const
Returns the total feature count (active and inactive ids, useful for sizing vectors) ...
void communicateAndMerge()
This routine handles all of the serialization, communication and deserialization of the data structur...
virtual void initialize() override
const Parallel::Communicator & _communicator
void addRelationshipManager(const std::string &name, Moose::RelationshipManagerType rm_type, Moose::RelationshipManagerInputParameterCallback input_parameter_callback=nullptr)
std::pair< typename M::iterator, bool > moose_try_emplace(M &m, const typename M::key_type &k, Args &&... args)
virtual bool areFeaturesMergeable(const FeatureData &f1, const FeatureData &f2) const
Method for determining whether two features are mergeable.
void mooseWarning(Args &&... args) const
std::map< dof_id_type, int > _ghosted_entity_ids
The map for holding reconstructed ghosted element information.
auto max(const L &left, const R &right)
virtual bool isFeaturePercolated(unsigned int feature_id) const
Returns a Boolean indicating whether this feature is percolated (e.g.
virtual void buildLocalToGlobalIndices(std::vector< std::size_t > &local_to_global_all, std::vector< int > &counts) const
This routine populates a stacked vector of local to global indices per rank and the associated count ...
SubProblem & _subproblem
void storeHelper(std::ostream &stream, P &data, void *context)
void sortAndLabel()
Sort and assign ids to features based on their position in the container after sorting.
MPI_Status status
static InputParameters validParams()
virtual void consolidateMergedFeatures(std::vector< std::list< FeatureData >> *saved_data=nullptr)
This method consolidates all of the merged information from _partial_feature_sets into the _feature_s...
uint8_t processor_id_type
static bool contains(std::set< T > &container, const T &item)
virtual unsigned int getFeatureVar(unsigned int feature_id) const
Returns the variable representing the passed in feature.
processor_id_type n_processors() const
const dof_id_type n_nodes
dof_id_type _min_entity_id
The minimum entity seen in the _local_ids, used for sorting features.
std::vector< std::map< dof_id_type, int > > _halo_ids
The data structure for looking up halos around features.
std::vector< MooseVariable * > _vars
The vector of coupled in variables cast to MooseVariable.
const Real _connecting_threshold
The threshold above (or below) which neighboring entities are flooded (where regions can be extended ...
virtual void execute() override
MeshBase & getMesh()
void dataStore(std::ostream &stream, FeatureFloodCount::FeatureData &feature, void *context)
std::size_t _var_index
The Moose variable where this feature was found (often the "order parameter")
virtual void clearDataStructures()
Helper routine for clearing up data structures during initialize and prior to parallel communication...
boundary_id_type BoundaryID
virtual const Node * nodePtr(const dof_id_type i) const
virtual bool isTransient() const=0
BoundaryIntersection
This enumeration is used to inidacate status of boundary intersections.
unsigned long _var_number
This variable is used to build the periodic node map.
This object will mark nodes or elements of continuous regions all with a unique number for the purpos...
virtual processor_id_type numberOfDistributedMergeHelpers() const
Returns a number indicating the number of merge helpers when running in parallel based on certain imp...
container_type _periodic_nodes
Holds the nodes that belong to the feature on a periodic boundary.
void visitElementalNeighbors(const Elem *elem, FeatureData *feature, bool expand_halos_only, bool disjoint_only)
void split(int color, int key, Communicator &target) const
const std::size_t _maps_size
Convenience variable holding the size of all the datastructures size by the number of maps...
const bool _use_less_than_threshold_comparison
Use less-than when comparing values against the threshold value.
bool ghostedIntersect(const FeatureData &rhs) const
const bool _global_numbering
This variable is used to indicate whether or not we identify features with unique numbers on multiple...
void dataLoad(std::istream &stream, FeatureFloodCount::FeatureData &feature, void *context)
const DofMap & _dof_map
Reference to the dof_map containing the coupled variables.
static const unsigned int invalid_id
NonlinearSystemBase & getNonlinearSystemBase(const unsigned int sys_num)
void paramError(const std::string &param, Args... args) const
unsigned int number() const
bool canConsolidate(const FeatureData &rhs) const
This routine indicates whether two features can be consolidated, that is, one feature is reasonably e...
std::vector< std::map< dof_id_type, int > > _feature_maps
The feature maps contain the raw flooded node information and eventually the unique grain numbers...
virtual Real getConnectingThreshold(std::size_t current_index) const
Return the "connecting" comparison threshold to use when inspecting an entity during the flood stage...
const bool _single_map_mode
This variable is used to indicate whether or not multiple maps are used during flooding.
unsigned int _feature_count
The number of features seen by this object (same as summing _feature_counts_per_map) ...
virtual Real getValue() const override
virtual void setCurrentSubdomainID(const Elem *elem, const THREAD_ID tid) override
void broadcast(T &data, const unsigned int root_id=0, const bool identical_sizes=false) const
static InputParameters validParams()
virtual bool doesFeatureIntersectSpecifiedBoundary(unsigned int feature_id) const
Returns a Boolean indicating whether this feature intersects boundaries in a user-supplied list...
const Real _threshold
The threshold above (or below) where an entity may begin a new region (feature)
void consolidate(FeatureData &&rhs)
Consolidates features, i.e.
bool hasBoundary(const BoundaryName &name) const
void addRequiredCoupledVar(const std::string &name, const std::string &doc_string)
std::unique_ptr< PointLocatorBase > _point_locator
virtual void mergeSets()
This routine is called on the primary rank only and stitches together the partial feature pieces seen...
bool _is_boundary_restricted
Indicates that this object should only run on one or more boundaries.
std::vector< MooseVariableFEBase * > _fe_vars
The vector of coupled in variables.
virtual bool isNewFeatureOrConnectedRegion(const DofObject *dof_object, std::size_t &current_index, FeatureData *&feature, Status &status, unsigned int &new_id)
Method called during the recursive flood routine that should return whether or not the current entity...
std::vector< std::list< FeatureData > > _partial_feature_sets
The data structure used to hold partial and communicated feature data, during the discovery and mergi...
static void reserve(std::set< T > &, std::size_t)
virtual void reinitElemPhys(const Elem *elem, const std::vector< Point > &phys_points_in_elem, const THREAD_ID tid)=0
void merge(FeatureData &&rhs)
Merges another Feature Data into this one.
virtual void restoreOriginalDataStructures(std::vector< std::list< FeatureData >> &)
const bool _compute_halo_maps
Indicates whether or not to communicate halo map information with all ranks.
PeriodicBoundaries * _pbs
A pointer to the periodic boundary constraints object.
void addMooseVariableDependency(MooseVariableFieldBase *var)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
std::unordered_map< dof_id_type, std::vector< const Elem * > > _nodes_to_elem_map
The data structure used to find neighboring elements give a node ID.
const bool _is_elemental
Determines if the flood counter is elements or not (nodes)
std::vector< unsigned int > _empty_var_to_features
void visitNeighborsHelper(const T *curr_entity, std::vector< const T *> neighbor_entities, FeatureData *feature, bool expand_halos_only, bool topological_neighbor, bool disjoint_only)
The actual logic for visiting neighbors is abstracted out here.
bool flood(const DofObject *dof_object, std::size_t current_index)
This method will check if the current entity is above the supplied threshold and "mark" it...
container_type _local_ids
Holds the local ids in the interior of a feature.
FEProblemBase & _fe_problem
OStreamProxy out
const processor_id_type _n_procs
Convenience variable holding the number of processors in this simulation.
IntRange< T > make_range(T beg, T end)
void updateBoundaryIntersections(FeatureData &feature) const
Update the feature&#39;s attributes to indicate boundary intersections.
virtual void meshChanged() override
std::vector< unsigned int > _feature_counts_per_map
The number of features seen by this object per map.
void mooseError(Args &&... args) const
std::vector< std::size_t > _feature_id_to_local_index
The vector recording the grain_id to local index (several indices will contain invalid_size_t) ...
ConstBndElemRange * _bnd_elem_range
Boundary element range pointer.
static const processor_id_type invalid_proc_id
void addClassDescription(const std::string &doc_string)
const InputParameters & parameters() const
std::list< std::pair< processor_id_type, unsigned int > > _orig_ids
Original processor/local ids.
std::vector< BoundaryID > getBoundaryIDs(const Elem *const elem, const unsigned short int side) const
virtual void initialSetup() override
BoundaryIntersection _boundary_intersection
Enumeration indicating boundary intersection status.
void updateBBoxExtremesHelper(BoundingBox &bbox, const Point &node)
Status _status
The status of a feature (used mostly in derived classes like the GrainTracker)
StoredRange< MooseMesh::const_bnd_elem_iterator, const BndElement *> * getBoundaryElementRange()
bool compareValueWithThreshold(Real entity_value, Real threshold) const
This method is used to determine whether the current entity value is part of a feature or not...
bool isBoundaryElem(dof_id_type elem_id) const
void deserialize(std::vector< std::string > &serialized_buffers, unsigned int var_num=invalid_id)
This routine takes the vector of byte buffers (one for each processor), deserializes them into a seri...
void mergeBBoxes(std::vector< BoundingBox > &bboxes, bool physical_intersection)
Located the overlapping bounding box between this Feature and the other Feature and expands that over...
MooseMesh & _mesh
A reference to the mesh.
processor_id_type processor_id() const
virtual Real getEntityValue(dof_id_type entity_id, FieldType field_type, std::size_t var_index=0) const
void buildFeatureIdToLocalIndices(unsigned int max_id)
This method builds a lookup map for retrieving the right local feature (by index) given a global inde...
void scatterAndUpdateRanks()
Calls buildLocalToGlobalIndices to build the individual local to global indicies for each rank and sc...
virtual bool doesFeatureIntersectBoundary(unsigned int feature_id) const
Returns a Boolean indicating whether this feature intersects any boundary.
std::deque< const DofObject * > _entity_queue
The data structure for maintaining entities to flood during discovery.
bool periodicBoundariesIntersect(const FeatureData &rhs) const
void loadHelper(std::istream &stream, P &data, void *context)
const bool _compute_var_to_feature_map
Indicates whether or not the var to feature map is populated.
auto index_range(const T &sizable)
void buildPeriodicNodeMap(std::multimap< dof_id_type, dof_id_type > &periodic_node_map, unsigned int var_number, PeriodicBoundaries *pbs) const
std::unordered_set< dof_id_type > _all_boundary_entity_ids
The set of entities on the boundary of the domain used for determining if features intersect any boun...
const bool _var_index_mode
This variable is used to indicate whether the maps will contain unique region information or just the...
std::size_t _vol_count
The count of entities contributing to the volume calculation.
unsigned int _id
An ID for this feature.
uint8_t dof_id_type
std::vector< std::map< dof_id_type, int > > _var_index_maps
This map keeps track of which variables own which nodes.
virtual Point featureCentroid(unsigned int feature_id) const
Returns the centroid of the designated feature (only supported without periodic boundaries) ...
void addParamNamesToGroup(const std::string &space_delim_names, const std::string group_name)
virtual Real getThreshold(std::size_t current_index) const
Return the starting comparison threshold to use when inspecting an entity during the flood stage...
std::vector< std::size_t > _local_to_global_feature_map
The vector recording the local to global feature indices.
bool isParamValid(const std::string &name) const
virtual void prepareDataForTransfer()
This routine uses the local flooded data to build up the local feature data structures (_partial feat...