www.mooseframework.org
FeatureFloodCount.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
3 /* */
4 /* All contents are licensed under LGPL V2.1 */
5 /* See LICENSE for full restrictions */
6 /****************************************************************/
7 
8 #include "FeatureFloodCount.h"
9 #include "IndirectSort.h"
10 #include "MooseMesh.h"
11 #include "MooseUtils.h"
12 #include "MooseVariable.h"
13 #include "SubProblem.h"
14 
15 #include "Assembly.h"
16 #include "FEProblem.h"
17 #include "NonlinearSystem.h"
18 
19 #include "libmesh/dof_map.h"
20 #include "libmesh/mesh_tools.h"
21 #include "libmesh/periodic_boundaries.h"
22 #include "libmesh/point_locator_base.h"
23 
24 #include <algorithm>
25 #include <limits>
26 
27 template <>
28 void
29 dataStore(std::ostream & stream, FeatureFloodCount::FeatureData & feature, void * context)
30 {
35  storeHelper(stream, feature._ghosted_ids, context);
36  storeHelper(stream, feature._halo_ids, context);
37  storeHelper(stream, feature._disjoint_halo_ids, context);
38  storeHelper(stream, feature._periodic_nodes, context);
39  storeHelper(stream, feature._var_index, context);
40  storeHelper(stream, feature._id, context);
41  storeHelper(stream, feature._bboxes, context);
42  storeHelper(stream, feature._orig_ids, context);
43  storeHelper(stream, feature._min_entity_id, context);
44  storeHelper(stream, feature._vol_count, context);
45  storeHelper(stream, feature._centroid, context);
46  storeHelper(stream, feature._status, context);
47  storeHelper(stream, feature._intersects_boundary, context);
48 }
49 
50 template <>
51 void
52 dataStore(std::ostream & stream, MeshTools::BoundingBox & bbox, void * context)
53 {
54  storeHelper(stream, bbox.min(), context);
55  storeHelper(stream, bbox.max(), context);
56 }
57 
58 template <>
59 void
60 dataLoad(std::istream & stream, FeatureFloodCount::FeatureData & feature, void * context)
61 {
66  loadHelper(stream, feature._ghosted_ids, context);
67  loadHelper(stream, feature._halo_ids, context);
68  loadHelper(stream, feature._disjoint_halo_ids, context);
69  loadHelper(stream, feature._periodic_nodes, context);
70  loadHelper(stream, feature._var_index, context);
71  loadHelper(stream, feature._id, context);
72  loadHelper(stream, feature._bboxes, context);
73  loadHelper(stream, feature._orig_ids, context);
74  loadHelper(stream, feature._min_entity_id, context);
75  loadHelper(stream, feature._vol_count, context);
76  loadHelper(stream, feature._centroid, context);
77  loadHelper(stream, feature._status, context);
78  loadHelper(stream, feature._intersects_boundary, context);
79 }
80 
81 template <>
82 void
83 dataLoad(std::istream & stream, MeshTools::BoundingBox & bbox, void * context)
84 {
85  loadHelper(stream, bbox.min(), context);
86  loadHelper(stream, bbox.max(), context);
87 }
88 
89 // Utility routines
90 void updateBBoxExtremesHelper(MeshTools::BoundingBox & bbox, const Point & node);
91 void updateBBoxExtremesHelper(MeshTools::BoundingBox & bbox, const Elem & elem);
92 bool areElemListsMergeable(const std::list<dof_id_type> & elem_list1,
93  const std::list<dof_id_type> & elem_list2,
94  MeshBase & mesh);
95 
96 template <>
97 InputParameters
99 {
100  InputParameters params = validParams<GeneralPostprocessor>();
101  params.addRequiredCoupledVar(
102  "variable",
103  "The variable(s) for which to find connected regions of interests, i.e. \"features\".");
104  params.addParam<Real>(
105  "threshold", 0.5, "The threshold value for which a new feature may be started");
106  params.addParam<Real>(
107  "connecting_threshold",
108  "The threshold for which an existing feature may be extended (defaults to \"threshold\")");
109  params.addParam<bool>("use_single_map",
110  true,
111  "Determine whether information is tracked per "
112  "coupled variable or consolidated into one "
113  "(default: true)");
114  params.addParam<bool>(
115  "condense_map_info",
116  false,
117  "Determines whether we condense all the node values when in multimap mode (default: false)");
118  params.addParam<bool>("use_global_numbering",
119  true,
120  "Determine whether or not global numbers are "
121  "used to label features on multiple maps "
122  "(default: true)");
123  params.addParam<bool>("enable_var_coloring",
124  false,
125  "Instruct the Postprocessor to populate the variable index map.");
126  params.addParam<bool>(
127  "compute_halo_maps",
128  false,
129  "Instruct the Postprocessor to communicate proper halo information to all ranks");
130  params.addParam<bool>("compute_var_to_feature_map",
131  false,
132  "Instruct the Postprocessor to compute the active vars to features map");
133  params.addParam<bool>(
134  "use_less_than_threshold_comparison",
135  true,
136  "Controls whether features are defined to be less than or greater than the threshold value.");
137 
145  params.set<bool>("use_displaced_mesh") = false;
146 
147  params.addParamNamesToGroup("use_single_map condense_map_info use_global_numbering", "Advanced");
148 
149  MooseEnum flood_type("NODAL ELEMENTAL", "ELEMENTAL");
150  params.addParam<MooseEnum>("flood_entity_type",
151  flood_type,
152  "Determines whether the flood algorithm runs on nodes or elements");
153  return params;
154 }
155 
156 FeatureFloodCount::FeatureFloodCount(const InputParameters & parameters)
157  : GeneralPostprocessor(parameters),
158  Coupleable(this, false),
159  MooseVariableDependencyInterface(),
160  ZeroInterface(parameters),
161  _vars(getCoupledMooseVars()),
162  _threshold(getParam<Real>("threshold")),
163  _connecting_threshold(isParamValid("connecting_threshold")
164  ? getParam<Real>("connecting_threshold")
165  : getParam<Real>("threshold")),
166  _mesh(_subproblem.mesh()),
167  _var_number(_vars[0]->number()),
168  _single_map_mode(getParam<bool>("use_single_map")),
169  _condense_map_info(getParam<bool>("condense_map_info")),
170  _global_numbering(getParam<bool>("use_global_numbering")),
171  _var_index_mode(getParam<bool>("enable_var_coloring")),
172  _compute_halo_maps(getParam<bool>("compute_halo_maps")),
173  _compute_var_to_feature_map(getParam<bool>("compute_var_to_feature_map")),
174  _use_less_than_threshold_comparison(getParam<bool>("use_less_than_threshold_comparison")),
175  _n_vars(_vars.size()),
176  _maps_size(_single_map_mode ? 1 : _vars.size()),
177  _n_procs(_app.n_processors()),
178  _feature_counts_per_map(_maps_size),
179  _feature_count(0),
180  _partial_feature_sets(_maps_size),
181  _feature_maps(_maps_size),
182  _pbs(nullptr),
183  _element_average_value(parameters.isParamValid("elem_avg_value")
184  ? getPostprocessorValue("elem_avg_value")
185  : _real_zero),
186  _halo_ids(_maps_size),
187  _is_elemental(getParam<MooseEnum>("flood_entity_type") == "ELEMENTAL"),
188  _is_master(processor_id() == 0)
189 {
190  if (_var_index_mode)
191  _var_index_maps.resize(_maps_size);
192 
193  addMooseVariableDependency(_vars);
194 }
195 
197 
198 void
200 {
201  // We need one map per coupled variable for normal runs to support overlapping features
202  _entities_visited.resize(_vars.size());
203 
204  // Get a pointer to the PeriodicBoundaries buried in libMesh
205  _pbs = _fe_problem.getNonlinearSystemBase().dofMap().get_periodic_boundaries();
206 
207  meshChanged();
208 
217 }
218 
219 void
221 {
222  // Clear the feature marking maps and region counters and other data structures
223  for (auto map_num = decltype(_maps_size)(0); map_num < _maps_size; ++map_num)
224  {
225  _feature_maps[map_num].clear();
226  _partial_feature_sets[map_num].clear();
227 
228  if (_var_index_mode)
229  _var_index_maps[map_num].clear();
230 
231  _halo_ids[map_num].clear();
232  }
233 
234  _feature_sets.clear();
235 
236  // Calculate the thresholds for this iteration
239 
240  _ghosted_entity_ids.clear();
241 
242  // Reset the feature count and max local size
243  _feature_count = 0;
244 
245  _entity_var_to_features.clear();
246 
247  for (auto & map_ref : _entities_visited)
248  map_ref.clear();
249 }
250 
251 void
253 {
254 }
255 
256 void
258 {
259  _point_locator = _mesh.getMesh().sub_point_locator();
260 
261  _mesh.buildPeriodicNodeMap(_periodic_node_map, _var_number, _pbs);
262 
263  // Build a new node to element map
264  _nodes_to_elem_map.clear();
265  MeshTools::build_nodes_to_elem_map(_mesh.getMesh(), _nodes_to_elem_map);
266 
272  _all_boundary_entity_ids.clear();
273  if (_is_elemental)
274  for (auto elem_it = _mesh.bndElemsBegin(), elem_end = _mesh.bndElemsEnd(); elem_it != elem_end;
275  ++elem_it)
276  _all_boundary_entity_ids.insert((*elem_it)->_elem->id());
277 }
278 
279 void
281 {
282  const auto end = _mesh.getMesh().active_local_elements_end();
283  for (auto el = _mesh.getMesh().active_local_elements_begin(); el != end; ++el)
284  {
285  const Elem * current_elem = *el;
286 
287  // Loop over elements or nodes
288  if (_is_elemental)
289  {
290  for (auto var_num = beginIndex(_vars); var_num < _vars.size(); ++var_num)
291  flood(current_elem, var_num, nullptr /* Designates inactive feature */);
292  }
293  else
294  {
295  auto n_nodes = current_elem->n_vertices();
296  for (auto i = decltype(n_nodes)(0); i < n_nodes; ++i)
297  {
298  const Node * current_node = current_elem->get_node(i);
299 
300  for (auto var_num = beginIndex(_vars); var_num < _vars.size(); ++var_num)
301  flood(current_node, var_num, nullptr /* Designates inactive feature */);
302  }
303  }
304  }
305 }
306 
307 void
309 {
310  // First we need to transform the raw data into a usable data structure
312 
320  std::vector<std::string> send_buffers(1);
321 
328  std::vector<std::string> recv_buffers;
329  if (_is_master)
330  recv_buffers.reserve(_app.n_processors());
331 
332  serialize(send_buffers[0]);
333 
334  // Free up as much memory as possible here before we do global communication
336 
341  _communicator.gather_packed_range(0,
342  (void *)(nullptr),
343  send_buffers.begin(),
344  send_buffers.end(),
345  std::back_inserter(recv_buffers));
346 
347  if (_is_master)
348  {
349  // The root process now needs to deserialize and merge all of the data
350  deserialize(recv_buffers);
351  recv_buffers.clear();
352 
353  mergeSets();
354  }
355 
356  // Make sure that feature count is communicated to all ranks
357  _communicator.broadcast(_feature_count);
358 }
359 
360 void
362 {
363  mooseAssert(_is_master, "sortAndLabel can only be called on the master");
364 
370  std::sort(_feature_sets.begin(), _feature_sets.end());
371 
372 #ifndef NDEBUG
373 
378  unsigned int feature_offset = 0;
379  for (auto map_num = beginIndex(_feature_counts_per_map); map_num < _maps_size; ++map_num)
380  {
381  // Skip empty map checks
382  if (_feature_counts_per_map[map_num] == 0)
383  continue;
384 
385  // Check the begin and end of the current range
386  auto range_front = feature_offset;
387  auto range_back = feature_offset + _feature_counts_per_map[map_num] - 1;
388 
389  mooseAssert(range_front <= range_back && range_back < _feature_count,
390  "Indexing error in feature sets");
391 
392  if (!_single_map_mode && (_feature_sets[range_front]._var_index != map_num ||
393  _feature_sets[range_back]._var_index != map_num))
394  mooseError("Error in _feature_sets sorting, map index: ", map_num);
395 
396  feature_offset += _feature_counts_per_map[map_num];
397  }
398 #endif
399 
400  // Label the features with an ID based on the sorting (processor number independent value)
401  for (auto i = beginIndex(_feature_sets); i < _feature_sets.size(); ++i)
402  if (_feature_sets[i]._id == invalid_id)
403  _feature_sets[i]._id = i;
404 }
405 
406 void
407 FeatureFloodCount::buildLocalToGlobalIndices(std::vector<std::size_t> & local_to_global_all,
408  std::vector<int> & counts) const
409 {
410  mooseAssert(_is_master, "This method must only be called on the root processor");
411 
412  counts.assign(_n_procs, 0);
413  // Now size the individual counts vectors based on the largest index seen per processor
414  for (const auto & feature : _feature_sets)
415  for (const auto & local_index_pair : feature._orig_ids)
416  {
417  // local_index_pair.first = ranks, local_index_pair.second = local_index
418  mooseAssert(local_index_pair.first < _n_procs, "Processor ID is out of range");
419  if (local_index_pair.second >= static_cast<std::size_t>(counts[local_index_pair.first]))
420  counts[local_index_pair.first] = local_index_pair.second + 1;
421  }
422 
423  // Build the offsets vector
424  unsigned int globalsize = 0;
425  std::vector<int> offsets(_n_procs); // Type is signed for use with the MPI API
426  for (auto i = beginIndex(offsets); i < offsets.size(); ++i)
427  {
428  offsets[i] = globalsize;
429  globalsize += counts[i];
430  }
431 
432  // Finally populate the master vector
433  local_to_global_all.resize(globalsize, FeatureFloodCount::invalid_size_t);
434  for (const auto & feature : _feature_sets)
435  {
436  // Get the local indices from the feature and build a map
437  for (const auto & local_index_pair : feature._orig_ids)
438  {
439  auto rank = local_index_pair.first;
440  mooseAssert(rank < _n_procs, rank << ", " << _n_procs);
441 
442  auto local_index = local_index_pair.second;
443  auto stacked_local_index = offsets[rank] + local_index;
444 
445  mooseAssert(stacked_local_index < globalsize,
446  "Global index: " << stacked_local_index << " is out of range");
447  local_to_global_all[stacked_local_index] = feature._id;
448  }
449  }
450 }
451 
452 void
454 {
455  _feature_id_to_local_index.assign(max_id + 1, invalid_size_t);
456  for (auto feature_index = beginIndex(_feature_sets); feature_index < _feature_sets.size();
457  ++feature_index)
458  {
459  if (_feature_sets[feature_index]._status != Status::INACTIVE)
460  {
461  mooseAssert(_feature_sets[feature_index]._id <= max_id,
462  "Feature ID out of range(" << _feature_sets[feature_index]._id << ')');
463  _feature_id_to_local_index[_feature_sets[feature_index]._id] = feature_index;
464  }
465  }
466 }
467 
468 void
470 {
471  // Gather all information on processor zero and merge
473 
474  // Sort and label the features
475  if (_is_master)
476  sortAndLabel();
477 
478  // Send out the local to global mappings
480 
481  // Populate _feature_maps and _var_index_maps
482  updateFieldInfo();
483 }
484 
485 const std::vector<unsigned int> &
487 {
488  mooseDoOnce(if (!_compute_var_to_feature_map) mooseError(
489  "Please set \"compute_var_to_feature_map = true\" to use this interface method"));
490 
491  const auto pos = _entity_var_to_features.find(elem_id);
492  if (pos != _entity_var_to_features.end())
493  {
494  mooseAssert(pos->second.size() == _n_vars, "Variable to feature vector not sized properly");
495  return pos->second;
496  }
497  else
498  return _empty_var_to_features;
499 }
500 
501 void
503 {
504  // local to global map (one per processor)
505  std::vector<int> counts;
506  std::vector<std::size_t> local_to_global_all;
507  if (_is_master)
508  buildLocalToGlobalIndices(local_to_global_all, counts);
509 
510  // Scatter local_to_global indices to all processors and store in class member variable
511  _communicator.scatter(local_to_global_all, counts, _local_to_global_feature_map);
512 
513  std::size_t largest_global_index = std::numeric_limits<std::size_t>::lowest();
514  if (!_is_master)
515  {
517 
524  for (auto & list_ref : _partial_feature_sets)
525  {
526  for (auto & feature : list_ref)
527  {
528  mooseAssert(feature._orig_ids.size() == 1, "feature._orig_ids length doesn't make sense");
529 
530  auto global_index = FeatureFloodCount::invalid_size_t;
531  auto local_index = feature._orig_ids.begin()->second;
532 
533  if (local_index < _local_to_global_feature_map.size())
534  global_index = _local_to_global_feature_map[local_index];
535 
536  if (global_index != FeatureFloodCount::invalid_size_t)
537  {
538  if (global_index > largest_global_index)
539  largest_global_index = global_index;
540 
541  // Set the correct global index
542  feature._id = global_index;
543 
551  feature._status &= ~Status::INACTIVE;
552 
553  // Move the feature into the correct place
554  _feature_sets[local_index] = std::move(feature);
555  }
556  }
557  }
558  }
559  else
560  {
561  for (auto global_index : local_to_global_all)
562  if (global_index != FeatureFloodCount::invalid_size_t && global_index > largest_global_index)
563  largest_global_index = global_index;
564  }
565 
566  buildFeatureIdToLocalIndices(largest_global_index);
567 }
568 
569 Real
571 {
572  return static_cast<Real>(_feature_count);
573 }
574 
575 std::size_t
577 {
583  return _feature_count;
584 }
585 
586 unsigned int
587 FeatureFloodCount::getFeatureVar(unsigned int feature_id) const
588 {
589  // Some processors don't contain the largest feature id, in that case we just return invalid_id
590  if (feature_id >= _feature_id_to_local_index.size())
591  return invalid_id;
592 
593  auto local_index = _feature_id_to_local_index[feature_id];
594  if (local_index != invalid_size_t)
595  {
596  mooseAssert(local_index < _feature_sets.size(), "local_index out of bounds");
597  return _feature_sets[local_index]._status != Status::INACTIVE
598  ? _feature_sets[feature_id]._var_index
599  : invalid_id;
600  }
601 
602  return invalid_id;
603 }
604 
605 bool
607 {
608  // TODO: This information is not parallel consistent when using FeatureFloodCounter
609 
610  // Some processors don't contain the largest feature id, in that case we just return invalid_id
611  if (feature_id >= _feature_id_to_local_index.size())
612  return false;
613 
614  auto local_index = _feature_id_to_local_index[feature_id];
615 
616  if (local_index != invalid_size_t)
617  {
618  mooseAssert(local_index < _feature_sets.size(), "local_index out of bounds");
619  return _feature_sets[local_index]._intersects_boundary;
620  }
621 
622  return false;
623 }
624 
625 Real
626 FeatureFloodCount::getEntityValue(dof_id_type entity_id,
627  FieldType field_type,
628  std::size_t var_index) const
629 {
630  auto use_default = false;
631  if (var_index == invalid_size_t)
632  {
633  use_default = true;
634  var_index = 0;
635  }
636 
637  mooseAssert(var_index < _maps_size, "Index out of range");
638 
639  switch (field_type)
640  {
642  {
643  const auto entity_it = _feature_maps[var_index].find(entity_id);
644 
645  if (entity_it != _feature_maps[var_index].end())
646  return entity_it->second; // + _region_offsets[var_index];
647  else
648  return -1;
649  }
650 
652  {
653  const auto entity_it = _var_index_maps[var_index].find(entity_id);
654 
655  if (entity_it != _var_index_maps[var_index].end())
656  return entity_it->second;
657  else
658  return -1;
659  }
660 
662  {
663  const auto entity_it = _ghosted_entity_ids.find(entity_id);
664 
665  if (entity_it != _ghosted_entity_ids.end())
666  return entity_it->second;
667  else
668  return -1;
669  }
670 
671  case FieldType::HALOS:
672  {
673  if (!use_default)
674  {
675  const auto entity_it = _halo_ids[var_index].find(entity_id);
676  if (entity_it != _halo_ids[var_index].end())
677  return entity_it->second;
678  }
679  else
680  {
681  // Showing halos in reverse order for backwards compatibility
682  for (auto map_num = _maps_size;
683  map_num-- /* don't compare greater than zero for unsigned */;)
684  {
685  const auto entity_it = _halo_ids[map_num].find(entity_id);
686 
687  if (entity_it != _halo_ids[map_num].end())
688  return entity_it->second;
689  }
690  }
691  return -1;
692  }
693 
694  case FieldType::CENTROID:
695  {
696  if (_periodic_node_map.size())
697  mooseDoOnce(mooseWarning(
698  "Centroids are not correct when using periodic boundaries, contact the MOOSE team"));
699 
700  // If this element contains the centroid of one of features, return one
701  const auto * elem_ptr = _mesh.elemPtr(entity_id);
702 
703  for (const auto & feature : _feature_sets)
704  {
705  if (feature._status == Status::INACTIVE)
706  continue;
707 
708  if (elem_ptr->contains_point(feature._centroid))
709  return 1;
710  }
711 
712  return 0;
713  }
714 
715  default:
716  return 0;
717  }
718 }
719 
720 void
722 {
723  MeshBase & mesh = _mesh.getMesh();
724 
725  std::set<dof_id_type> local_ids_no_ghost, set_difference;
726 
727  for (auto & list_ref : _partial_feature_sets)
728  {
729  for (auto & feature : list_ref)
730  {
731  // Now extend the bounding box by the halo region
732  if (_is_elemental)
733  feature.updateBBoxExtremes(mesh);
734  else
735  {
736  for (auto & halo_id : feature._halo_ids)
737  updateBBoxExtremesHelper(feature._bboxes[0], mesh.node(halo_id));
738  }
739 
746  std::set_difference(feature._local_ids.begin(),
747  feature._local_ids.end(),
748  feature._ghosted_ids.begin(),
749  feature._ghosted_ids.end(),
750  std::insert_iterator<std::set<dof_id_type>>(local_ids_no_ghost,
751  local_ids_no_ghost.begin()));
752 
753  std::set_difference(
754  feature._halo_ids.begin(),
755  feature._halo_ids.end(),
756  local_ids_no_ghost.begin(),
757  local_ids_no_ghost.end(),
758  std::insert_iterator<std::set<dof_id_type>>(set_difference, set_difference.begin()));
759  feature._halo_ids.swap(set_difference);
760  local_ids_no_ghost.clear();
761  set_difference.clear();
762 
763  mooseAssert(!feature._local_ids.empty(), "local entity ids cannot be empty");
764 
769  feature._min_entity_id = *feature._local_ids.begin();
770 
771  // Periodic node ids
773  }
774  }
775 }
776 
777 void
778 FeatureFloodCount::serialize(std::string & serialized_buffer)
779 {
780  // stream for serializing the _partial_feature_sets data structure to a byte stream
781  std::ostringstream oss;
782 
787  dataStore(oss, _partial_feature_sets, this);
788 
789  // Populate the passed in string pointer with the string stream's buffer contents
790  serialized_buffer.assign(oss.str());
791 }
792 
800 void
801 FeatureFloodCount::deserialize(std::vector<std::string> & serialized_buffers)
802 {
803  // The input string stream used for deserialization
804  std::istringstream iss;
805 
806  mooseAssert(serialized_buffers.size() == _app.n_processors(),
807  "Unexpected size of serialized_buffers: " << serialized_buffers.size());
808  auto rank = processor_id();
809  for (auto proc_id = beginIndex(serialized_buffers); proc_id < serialized_buffers.size();
810  ++proc_id)
811  {
816  if (proc_id == rank)
817  continue;
818 
819  iss.str(serialized_buffers[proc_id]); // populate the stream with a new buffer
820  iss.clear(); // reset the string stream state
821 
822  // Load the communicated data into all of the other processors' slots
823  dataLoad(iss, _partial_feature_sets, this);
824  }
825 }
826 
827 void
829 {
830  Moose::perf_log.push("mergeSets()", "FeatureFloodCount");
831 
832  // Since we gathered only on the root process, we only need to merge sets on the root process.
833  mooseAssert(_is_master, "mergeSets() should only be called on the root process");
834 
835  // Local variable used for sizing structures, it will be >= the actual number of features
836  for (auto map_num = decltype(_maps_size)(0); map_num < _maps_size; ++map_num)
837  {
838  for (auto it1 = _partial_feature_sets[map_num].begin();
839  it1 != _partial_feature_sets[map_num].end();
840  /* No increment on it1 */)
841  {
842  bool merge_occured = false;
843  for (auto it2 = _partial_feature_sets[map_num].begin();
844  it2 != _partial_feature_sets[map_num].end();
845  ++it2)
846  {
847  if (it1 != it2 && areFeaturesMergeable(*it1, *it2))
848  {
849  it2->merge(std::move(*it1));
850 
855  _partial_feature_sets[map_num].emplace_back(std::move(*it2));
856 
866  _partial_feature_sets[map_num].erase(it2);
867  it1 = _partial_feature_sets[map_num].erase(it1); // it1 is incremented here!
868 
869  // A merge occurred, this is used to determine whether or not we increment the outer
870  // iterator
871  merge_occured = true;
872 
873  // We need to start the list comparison over for the new it1 so break here
874  break;
875  }
876  } // it2 loop
877 
878  if (!merge_occured) // No merges so we need to manually increment the outer iterator
879  ++it1;
880 
881  } // it1 loop
882  } // map loop
883 
891  // Offset where the current set of features with the same variable id starts in the flat vector
892  unsigned int feature_offset = 0;
893  // Set the member feature count to zero and start counting the actual features
894  _feature_count = 0;
895 
896  for (auto map_num = decltype(_maps_size)(0); map_num < _maps_size; ++map_num)
897  {
898  std::set<dof_id_type> set_difference;
899  for (auto & feature : _partial_feature_sets[map_num])
900  {
901  // If after merging we still have an inactive feature, discard it
902  if (feature._status == Status::CLEAR)
903  {
904  // First we need to calculate the centroid now that we are doing merging all partial
905  // features
906  if (feature._vol_count != 0)
907  feature._centroid /= feature._vol_count;
908 
909  _feature_sets.emplace_back(std::move(feature));
910  ++_feature_count;
911  }
912  }
913 
914  // Record the feature numbers just for the current map
915  _feature_counts_per_map[map_num] = _feature_count - feature_offset;
916 
917  // Now update the running feature count so we can calculate the next map's contribution
918  feature_offset = _feature_count;
919 
920  // Clean up the "moved" objects
921  _partial_feature_sets[map_num].clear();
922  }
923 
929  Moose::perf_log.pop("mergeSets()", "FeatureFloodCount");
930 }
931 
932 bool
934 {
935  return f1.mergeable(f2);
936 }
937 
938 void
940 {
941  for (auto i = beginIndex(_feature_sets); i < _feature_sets.size(); ++i)
942  {
943  auto & feature = _feature_sets[i];
944 
945  // If the developer has requested _condense_map_info we'll make sure we only update the zeroth
946  // map
947  auto map_index = (_single_map_mode || _condense_map_info) ? decltype(feature._var_index)(0)
948  : feature._var_index;
949 
950  // Loop over the entity ids of this feature and update our local map
951  for (auto entity : feature._local_ids)
952  {
953  _feature_maps[map_index][entity] = static_cast<int>(feature._id);
954 
955  if (_var_index_mode)
956  _var_index_maps[map_index][entity] = feature._var_index;
957 
958  // Fill in the data structure that keeps track of all features per elem
960  {
961  auto insert_pair = moose_try_emplace(
962  _entity_var_to_features, entity, std::vector<unsigned int>(_n_vars, invalid_id));
963  auto & vec_ref = insert_pair.first->second;
964  vec_ref[feature._var_index] = feature._id;
965  }
966  }
967 
968  if (_compute_halo_maps)
969  // Loop over the halo ids to update cells with halo information
970  for (auto entity : feature._halo_ids)
971  _halo_ids[map_index][entity] = static_cast<int>(feature._id);
972 
973  // Loop over the ghosted ids to update cells with ghost information
974  for (auto entity : feature._ghosted_ids)
975  _ghosted_entity_ids[entity] = 1;
976 
977  // TODO: Fixme
978  if (!_global_numbering)
979  mooseError("Local numbering currently disabled");
980  }
981 }
982 
983 bool
984 FeatureFloodCount::flood(const DofObject * dof_object,
985  std::size_t current_index,
986  FeatureData * feature)
987 {
988  if (dof_object == nullptr)
989  return false;
990 
991  // Retrieve the id of the current entity
992  auto entity_id = dof_object->id();
993 
994  // Has this entity already been marked? - if so move along
995  if (current_index != invalid_size_t &&
996  _entities_visited[current_index].find(entity_id) != _entities_visited[current_index].end())
997  return false;
998 
999  // See if the current entity either starts a new feature or continues an existing feature
1000  auto new_id = invalid_id; // Writable reference to hold an optional id;
1001  Status status =
1002  Status::INACTIVE; // Status is inactive until we find an entity above the starting threshold
1003  if (!isNewFeatureOrConnectedRegion(dof_object, current_index, feature, status, new_id))
1004  return false;
1005 
1006  mooseAssert(current_index != invalid_size_t, "current_index is invalid");
1007 
1016  _entities_visited[current_index].emplace(entity_id);
1017 
1018  auto map_num = _single_map_mode ? decltype(current_index)(0) : current_index;
1019 
1020  // New Feature (we need to create it and add it to our data structure)
1021  if (!feature)
1022  {
1023  _partial_feature_sets[map_num].emplace_back(
1024  current_index, _feature_count++, processor_id(), status);
1025 
1026  // Get a handle to the feature we will update (always the last feature in the data structure)
1027  feature = &_partial_feature_sets[map_num].back();
1028 
1029  // If new_id is valid, we'll set it in the feature here.
1030  if (new_id != invalid_id)
1031  feature->_id = new_id;
1032  }
1033 
1034  // Insert the current entity into the local ids map
1035  feature->_local_ids.insert(entity_id);
1036 
1042  if (_is_elemental && processor_id() == dof_object->processor_id())
1043  {
1044  const Elem * elem = static_cast<const Elem *>(dof_object);
1045 
1046  // Keep track of how many elements participate in the centroid averaging
1047  feature->_vol_count++;
1048 
1049  // Sum the centroid values for now, we'll average them later
1050  feature->_centroid += elem->centroid();
1051 
1052  // Does the volume intersect the boundary?
1053  if (_all_boundary_entity_ids.find(elem->id()) != _all_boundary_entity_ids.end())
1054  feature->_intersects_boundary = true;
1055  }
1056 
1057  if (_is_elemental)
1058  visitElementalNeighbors(static_cast<const Elem *>(dof_object),
1059  current_index,
1060  feature,
1061  /*expand_halos_only =*/false,
1062  /*disjoint_only =*/false);
1063  else
1064  visitNodalNeighbors(static_cast<const Node *>(dof_object),
1065  current_index,
1066  feature,
1067  /*expand_halos_only =*/false);
1068 
1069  return true;
1070 }
1071 
1072 Real FeatureFloodCount::getThreshold(std::size_t /*current_index*/) const
1073 {
1074  return _step_threshold;
1075 }
1076 
1077 Real FeatureFloodCount::getConnectingThreshold(std::size_t /*current_index*/) const
1078 {
1080 }
1081 
1082 bool
1083 FeatureFloodCount::compareValueWithThreshold(Real entity_value, Real threshold) const
1084 {
1085  return ((_use_less_than_threshold_comparison && (entity_value >= threshold)) ||
1086  (!_use_less_than_threshold_comparison && (entity_value <= threshold)));
1087 }
1088 
1089 bool
1091  std::size_t & current_index,
1092  FeatureData *& feature,
1093  Status & status,
1094  unsigned int & /*new_id*/)
1095 {
1096  // Get the value of the current variable for the current entity
1097  Real entity_value;
1098  if (_is_elemental)
1099  {
1100  const Elem * elem = static_cast<const Elem *>(dof_object);
1101  std::vector<Point> centroid(1, elem->centroid());
1102  _subproblem.reinitElemPhys(elem, centroid, 0);
1103  entity_value = _vars[current_index]->sln()[0];
1104  }
1105  else
1106  entity_value = _vars[current_index]->getNodalValue(*static_cast<const Node *>(dof_object));
1107 
1108  // If the value compares against our starting threshold, this is definitely part of a feature
1109  // we'll keep
1110  if (compareValueWithThreshold(entity_value, getThreshold(current_index)))
1111  {
1112  Status * status_ptr = &status;
1113 
1114  if (feature)
1115  status_ptr = &feature->_status;
1116 
1117  // Update an existing feature's status or clear the flag on the passed in status
1118  *status_ptr &= ~Status::INACTIVE;
1119  return true;
1120  }
1121 
1126  return compareValueWithThreshold(entity_value, getConnectingThreshold(current_index));
1127 }
1128 
1129 void
1131 {
1132  const auto & node_to_elem_map = _mesh.nodeToActiveSemilocalElemMap();
1133  decltype(FeatureData::_local_ids) expanded_local_ids;
1134  auto my_processor_id = processor_id();
1135 
1142  for (auto & list_ref : _partial_feature_sets)
1143  {
1144  for (auto & feature : list_ref)
1145  {
1146  expanded_local_ids.clear();
1147 
1148  for (auto entity : feature._local_ids)
1149  {
1150  const Elem * elem = _mesh.elemPtr(entity);
1151  mooseAssert(elem, "elem pointer is NULL");
1152 
1153  // Get the nodes on a current element so that we can add in point neighbors
1154  auto n_nodes = elem->n_vertices();
1155  for (auto i = decltype(n_nodes)(0); i < n_nodes; ++i)
1156  {
1157  const Node * current_node = elem->get_node(i);
1158 
1159  auto elem_vector_it = node_to_elem_map.find(current_node->id());
1160  if (elem_vector_it == node_to_elem_map.end())
1161  mooseError("Error in node to elem map");
1162 
1163  const auto & elem_vector = elem_vector_it->second;
1164 
1165  expanded_local_ids.insert(elem_vector.begin(), elem_vector.end());
1166 
1167  // Now see which elements need to go into the ghosted set
1168  for (auto entity : elem_vector)
1169  {
1170  const Elem * neighbor = _mesh.elemPtr(entity);
1171  mooseAssert(neighbor, "neighbor pointer is NULL");
1172 
1173  if (neighbor->processor_id() != my_processor_id)
1174  feature._ghosted_ids.insert(elem->id());
1175  }
1176  }
1177  }
1178 
1179  // Replace the existing local ids with the expanded local ids
1180  feature._local_ids.swap(expanded_local_ids);
1181 
1182  // Copy the expanded local_ids into the halo_ids container
1183  feature._halo_ids = feature._local_ids;
1184  }
1185  }
1186 }
1187 
1188 void
1189 FeatureFloodCount::expandEdgeHalos(unsigned int num_layers_to_expand)
1190 {
1191  if (num_layers_to_expand == 0)
1192  return;
1193 
1194  for (auto & list_ref : _partial_feature_sets)
1195  {
1196  for (auto & feature : list_ref)
1197  {
1198  for (auto halo_level = decltype(num_layers_to_expand)(0); halo_level < num_layers_to_expand;
1199  ++halo_level)
1200  {
1205  std::set<dof_id_type> orig_halo_ids(feature._halo_ids);
1206  for (auto entity : orig_halo_ids)
1207  {
1208  if (_is_elemental)
1209  visitElementalNeighbors(_mesh.elemPtr(entity),
1210  feature._var_index,
1211  &feature,
1212  /*expand_halos_only =*/true,
1213  /*disjoint_only =*/false);
1214  else
1215  visitNodalNeighbors(_mesh.nodePtr(entity),
1216  feature._var_index,
1217  &feature,
1218  /*expand_halos_only =*/true);
1219  }
1220 
1225  std::set<dof_id_type> disjoint_orig_halo_ids(feature._disjoint_halo_ids);
1226  for (auto entity : disjoint_orig_halo_ids)
1227  {
1228  if (_is_elemental)
1229  visitElementalNeighbors(_mesh.elemPtr(entity),
1230  feature._var_index,
1231  &feature,
1232  /*expand_halos_only =*/true,
1233  /*disjoint_only =*/true);
1234  else
1235  visitNodalNeighbors(_mesh.nodePtr(entity),
1236  feature._var_index,
1237  &feature,
1238  /*expand_halos_only =*/true);
1239  }
1240  }
1241  }
1242  }
1243 }
1244 
1245 void
1247  std::size_t current_index,
1248  FeatureData * feature,
1249  bool expand_halos_only,
1250  bool disjoint_only)
1251 {
1252  mooseAssert(elem, "Elem is NULL");
1253 
1254  std::vector<const Elem *> all_active_neighbors;
1255  MeshBase & mesh = _mesh.getMesh();
1256 
1257  // Loop over all neighbors (at the the same level as the current element)
1258  // Loop over all neighbors (at the the same level as the current element)
1259  for (auto i = decltype(elem->n_neighbors())(0); i < elem->n_neighbors(); ++i)
1260  {
1261  const Elem * neighbor_ancestor = nullptr;
1262  bool topological_neighbor = false;
1263 
1268  neighbor_ancestor = elem->neighbor(i);
1269  if (neighbor_ancestor)
1270  neighbor_ancestor->active_family_tree_by_neighbor(all_active_neighbors, elem, false);
1271  else // if (expand_halos_only /*&& feature->_periodic_nodes.empty()*/)
1272  {
1273  neighbor_ancestor = elem->topological_neighbor(i, mesh, *_point_locator, _pbs);
1274 
1282  if (neighbor_ancestor)
1283  {
1284  neighbor_ancestor->active_family_tree_by_topological_neighbor(
1285  all_active_neighbors, elem, mesh, *_point_locator, _pbs, false);
1286 
1287  topological_neighbor = true;
1288 
1289  // for (const auto neighbor : all_active_neighbors)
1290  // feature->_disjoint_halo_ids.insert(neighbor->id());
1291  }
1292  }
1293 
1294  visitNeighborsHelper(elem,
1295  all_active_neighbors,
1296  current_index,
1297  feature,
1298  expand_halos_only,
1299  topological_neighbor,
1300  disjoint_only);
1301 
1302  all_active_neighbors.clear();
1303  }
1304 }
1305 
1306 void
1308  std::size_t current_index,
1309  FeatureData * feature,
1310  bool expand_halos_only)
1311 {
1312  mooseAssert(node, "Node is NULL");
1313 
1314  std::vector<const Node *> all_active_neighbors;
1315  MeshTools::find_nodal_neighbors(_mesh.getMesh(), *node, _nodes_to_elem_map, all_active_neighbors);
1316 
1318  node, all_active_neighbors, current_index, feature, expand_halos_only, false, false);
1319 }
1320 
1321 template <typename T>
1322 void
1324  std::vector<const T *> neighbor_entities,
1325  std::size_t current_index,
1326  FeatureData * feature,
1327  bool expand_halos_only,
1328  bool topological_neighbor,
1329  bool disjoint_only)
1330 {
1331  // Loop over all active element neighbors
1332  for (const auto neighbor : neighbor_entities)
1333  {
1334  if (neighbor)
1335  {
1336  if (expand_halos_only)
1337  {
1338  if (topological_neighbor || disjoint_only)
1339  feature->_disjoint_halo_ids.insert(neighbor->id());
1340  else
1341  feature->_halo_ids.insert(neighbor->id());
1342  }
1343  else
1344  {
1345  auto my_processor_id = processor_id();
1346 
1347  if (!topological_neighbor && neighbor->processor_id() != my_processor_id)
1348  feature->_ghosted_ids.insert(curr_entity->id());
1349 
1359  if (curr_entity->processor_id() == my_processor_id)
1360  {
1367  if (topological_neighbor || disjoint_only)
1368  feature->_disjoint_halo_ids.insert(neighbor->id());
1369  else
1370  {
1371  feature->_halo_ids.insert(neighbor->id());
1372 
1373  flood(neighbor, current_index, feature);
1374  }
1375  }
1376  }
1377  }
1378  }
1379 }
1380 
1381 void
1383 {
1384  if (_is_elemental)
1385  {
1386  for (auto entity : feature._local_ids)
1387  {
1388  Elem * elem = _mesh.elemPtr(entity);
1389 
1390  for (auto node_n = decltype(elem->n_nodes())(0); node_n < elem->n_nodes(); ++node_n)
1391  {
1392  auto iters = _periodic_node_map.equal_range(elem->node(node_n));
1393 
1394  for (auto it = iters.first; it != iters.second; ++it)
1395  {
1396  feature._periodic_nodes.insert(it->first);
1397  feature._periodic_nodes.insert(it->second);
1398  }
1399  }
1400  }
1401  }
1402  else
1403  {
1404  for (auto entity : feature._local_ids)
1405  {
1406  auto iters = _periodic_node_map.equal_range(entity);
1407 
1408  for (auto it = iters.first; it != iters.second; ++it)
1409  {
1410  feature._periodic_nodes.insert(it->first);
1411  feature._periodic_nodes.insert(it->second);
1412  }
1413  }
1414  }
1415 }
1416 
1417 void
1419 {
1420  // First update the primary bounding box (all topologically connected)
1421  for (auto & halo_id : _halo_ids)
1422  updateBBoxExtremesHelper(_bboxes[0], *mesh.elem(halo_id));
1423 
1424  // Remove all of the IDs that are in the primary region
1425  std::list<dof_id_type> disjoint_elem_id_list;
1426  std::set_difference(_disjoint_halo_ids.begin(),
1427  _disjoint_halo_ids.end(),
1428  _halo_ids.begin(),
1429  _halo_ids.end(),
1430  std::insert_iterator<std::list<dof_id_type>>(disjoint_elem_id_list,
1431  disjoint_elem_id_list.begin()));
1432 
1433  if (!disjoint_elem_id_list.empty())
1434  {
1443  std::list<std::list<dof_id_type>> disjoint_regions;
1444  for (auto elem_id : _disjoint_halo_ids)
1445  {
1446  disjoint_regions.emplace_back(std::list<dof_id_type>({elem_id}));
1447  }
1448 
1449  for (auto it1 = disjoint_regions.begin(); it1 != disjoint_regions.end(); /* No increment */)
1450  {
1451  bool merge_occured = false;
1452  for (auto it2 = disjoint_regions.begin(); it2 != disjoint_regions.end(); ++it2)
1453  {
1454  if (it1 != it2 && areElemListsMergeable(*it1, *it2, mesh))
1455  {
1456  it2->splice(it2->begin(), *it1);
1457 
1458  disjoint_regions.emplace_back(std::move(*it2));
1459  disjoint_regions.erase(it2);
1460  it1 = disjoint_regions.erase(it1);
1461 
1462  merge_occured = true;
1463 
1464  break;
1465  }
1466  }
1467 
1468  if (!merge_occured)
1469  ++it1;
1470  }
1471 
1472  // Finally create new bounding boxes for each disjoint region
1473  auto num_regions = disjoint_regions.size();
1474  // We have num_regions *new* bounding boxes plus the existing bounding box
1475  _bboxes.resize(num_regions + 1);
1476 
1477  decltype(num_regions) region = 1;
1478  for (const auto list_ref : disjoint_regions)
1479  {
1480  for (const auto elem_id : list_ref)
1481  updateBBoxExtremesHelper(_bboxes[region], *mesh.elem_ptr(elem_id));
1482 
1483  _halo_ids.insert(_disjoint_halo_ids.begin(), _disjoint_halo_ids.end());
1484  _disjoint_halo_ids.clear();
1485  ++region;
1486  }
1487  }
1488 }
1489 
1490 void
1492  const MeshTools::BoundingBox & rhs_bbox)
1493 {
1494  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
1495  {
1496  bbox.min()(i) = std::min(bbox.min()(i), rhs_bbox.min()(i));
1497  bbox.max()(i) = std::max(bbox.max()(i), rhs_bbox.max()(i));
1498  }
1499 }
1500 
1501 bool
1503 {
1504  // See if any of the bounding boxes in either FeatureData object intersect
1505  for (const auto & bbox_lhs : _bboxes)
1506  for (const auto & bbox_rhs : rhs._bboxes)
1507  if (bbox_lhs.intersects(bbox_rhs))
1508  return true;
1509 
1510  return false;
1511 }
1512 
1513 bool
1515 {
1516  return setsIntersect(
1517  _halo_ids.begin(), _halo_ids.end(), rhs._halo_ids.begin(), rhs._halo_ids.end());
1518 }
1519 
1520 bool
1522 {
1523  return setsIntersect(_periodic_nodes.begin(),
1524  _periodic_nodes.end(),
1525  rhs._periodic_nodes.begin(),
1526  rhs._periodic_nodes.end());
1527 }
1528 
1529 bool
1531 {
1532  return setsIntersect(
1533  _ghosted_ids.begin(), _ghosted_ids.end(), rhs._ghosted_ids.begin(), rhs._ghosted_ids.end());
1534 }
1535 
1536 bool
1538 {
1539  return (_var_index == rhs._var_index && // the sets have matching variable indices and
1540  ((boundingBoxesIntersect(rhs) && // (if the feature's bboxes intersect and
1541  ghostedIntersect(rhs)) // the ghosted entities also intersect)
1542  || // or
1543  periodicBoundariesIntersect(rhs))); // periodic node sets intersect)
1544 }
1545 
1546 void
1548 {
1549  mooseAssert(_var_index == rhs._var_index, "Mismatched variable index in merge");
1550  mooseAssert(_id == rhs._id, "Mismatched auxiliary id in merge");
1551 
1552  std::set<dof_id_type> set_union;
1553 
1561  std::set_union(_periodic_nodes.begin(),
1562  _periodic_nodes.end(),
1563  rhs._periodic_nodes.begin(),
1564  rhs._periodic_nodes.end(),
1565  std::insert_iterator<std::set<dof_id_type>>(set_union, set_union.begin()));
1566  _periodic_nodes.swap(set_union);
1567 
1568  set_union.clear();
1569  std::set_union(_local_ids.begin(),
1570  _local_ids.end(),
1571  rhs._local_ids.begin(),
1572  rhs._local_ids.end(),
1573  std::insert_iterator<std::set<dof_id_type>>(set_union, set_union.begin()));
1574  _local_ids.swap(set_union);
1575 
1576  set_union.clear();
1577  std::set_union(_ghosted_ids.begin(),
1578  _ghosted_ids.end(),
1579  rhs._ghosted_ids.begin(),
1580  rhs._ghosted_ids.end(),
1581  std::insert_iterator<std::set<dof_id_type>>(set_union, set_union.begin()));
1582 
1583  // Was there overlap in the physical region?
1584  bool physical_intersection = (_ghosted_ids.size() + rhs._ghosted_ids.size() > set_union.size());
1585  _ghosted_ids.swap(set_union);
1586 
1591  if (physical_intersection)
1592  expandBBox(rhs);
1593  else
1594  {
1595  // set_union.clear();
1596  // std::set_difference(_halo_ids.begin(),
1597  // _halo_ids.end(),
1598  // _disjoint_halo_ids.begin(),
1599  // _disjoint_halo_ids.end(),
1600  // std::insert_iterator<std::set<dof_id_type>>(set_union,
1601  // set_union.begin()));
1602  // _halo_ids.swap(set_union);
1603  //
1604  // set_union.clear();
1605  // std::set_difference(rhs._halo_ids.begin(),
1606  // rhs._halo_ids.end(),
1607  // rhs._disjoint_halo_ids.begin(),
1608  // rhs._disjoint_halo_ids.end(),
1609  // std::insert_iterator<std::set<dof_id_type>>(set_union,
1610  // set_union.begin()));
1611  // rhs._halo_ids.swap(set_union);
1612 
1613  std::move(rhs._bboxes.begin(), rhs._bboxes.end(), std::back_inserter(_bboxes));
1614  }
1615 
1616  set_union.clear();
1617  std::set_union(_disjoint_halo_ids.begin(),
1618  _disjoint_halo_ids.end(),
1619  rhs._disjoint_halo_ids.begin(),
1620  rhs._disjoint_halo_ids.end(),
1621  std::insert_iterator<std::set<dof_id_type>>(set_union, set_union.begin()));
1622  _disjoint_halo_ids.swap(set_union);
1623 
1624  set_union.clear();
1625  std::set_union(_halo_ids.begin(),
1626  _halo_ids.end(),
1627  rhs._halo_ids.begin(),
1628  rhs._halo_ids.end(),
1629  std::insert_iterator<std::set<dof_id_type>>(set_union, set_union.begin()));
1630  _halo_ids.swap(set_union);
1631 
1632  // Keep track of the original ids so we can notify other processors of the local to global mapping
1633  _orig_ids.splice(_orig_ids.end(), std::move(rhs._orig_ids));
1634 
1635  // Update the min feature id
1636  _min_entity_id = std::min(_min_entity_id, rhs._min_entity_id);
1637 
1643  mooseAssert((_status & Status::MARKED & Status::DIRTY) == Status::CLEAR,
1644  "Flags in invalid state");
1645 
1646  // Logical AND here to combine flags (INACTIVE & INACTIVE == INACTIVE, all other combos are CLEAR)
1647  _status &= rhs._status;
1648 
1649  _vol_count += rhs._vol_count;
1650  _centroid += rhs._centroid;
1651 }
1652 
1653 void
1655 {
1656  _local_ids.clear();
1657  _periodic_nodes.clear();
1658  _halo_ids.clear();
1659  _disjoint_halo_ids.clear();
1660  _ghosted_ids.clear();
1661  _bboxes.clear();
1662  _orig_ids.clear();
1663 }
1664 
1665 void
1667 {
1668  std::vector<bool> intersected_boxes(rhs._bboxes.size(), false);
1669 
1670  auto box_expanded = false;
1671  for (auto & bbox : _bboxes)
1672  for (auto j = beginIndex(rhs._bboxes); j < rhs._bboxes.size(); ++j)
1673  if (bbox.intersects(rhs._bboxes[j]))
1674  {
1675  updateBBoxExtremes(bbox, rhs._bboxes[j]);
1676  intersected_boxes[j] = true;
1677  box_expanded = true;
1678  }
1679 
1680  // Any bounding box in the rhs vector that doesn't intersect
1681  // needs to be appended to the lhs vector
1682  for (auto j = beginIndex(intersected_boxes); j < intersected_boxes.size(); ++j)
1683  if (!intersected_boxes[j])
1684  _bboxes.push_back(rhs._bboxes[j]);
1685 
1686  // Error check
1687  if (!box_expanded)
1688  {
1689  std::ostringstream oss;
1690  oss << "LHS BBoxes:\n";
1691  for (auto i = beginIndex(_bboxes); i < _bboxes.size(); ++i)
1692  oss << "Max: " << _bboxes[i].max() << " Min: " << _bboxes[i].min() << '\n';
1693 
1694  oss << "RHS BBoxes:\n";
1695  for (auto i = beginIndex(rhs._bboxes); i < rhs._bboxes.size(); ++i)
1696  oss << "Max: " << rhs._bboxes[i].max() << " Min: " << rhs._bboxes[i].min() << '\n';
1697 
1698  ::mooseError("No Bounding Boxes Expanded - This is a catastrophic error!\n", oss.str());
1699  }
1700 }
1701 
1702 std::ostream &
1703 operator<<(std::ostream & out, const FeatureFloodCount::FeatureData & feature)
1704 {
1705  static const bool debug = true;
1706 
1707  out << "Grain ID: ";
1708  if (feature._id != FeatureFloodCount::invalid_id)
1709  out << feature._id;
1710  else
1711  out << "invalid";
1712 
1713  if (debug)
1714  {
1715  out << "\nGhosted Entities: ";
1716  for (auto ghosted_id : feature._ghosted_ids)
1717  out << ghosted_id << " ";
1718 
1719  out << "\nLocal Entities: ";
1720  for (auto local_id : feature._local_ids)
1721  out << local_id << " ";
1722 
1723  out << "\nHalo Entities: ";
1724  for (auto halo_id : feature._halo_ids)
1725  out << halo_id << " ";
1726 
1727  out << "\nPeriodic Node IDs: ";
1728  for (auto periodic_node : feature._periodic_nodes)
1729  out << periodic_node << " ";
1730  }
1731 
1732  out << "\nBBoxes:";
1733  Real volume = 0;
1734  for (const auto & bbox : feature._bboxes)
1735  {
1736  out << "\nMax: " << bbox.max() << " Min: " << bbox.min();
1737  volume += (bbox.max()(0) - bbox.min()(0)) * (bbox.max()(1) - bbox.min()(1)) *
1738  (MooseUtils::absoluteFuzzyEqual(bbox.max()(2), bbox.min()(2))
1739  ? 1
1740  : bbox.max()(2) - bbox.min()(2));
1741  }
1742 
1743  out << "\nStatus: ";
1745  out << "CLEAR";
1746  if (static_cast<bool>(feature._status & FeatureFloodCount::Status::MARKED))
1747  out << " MARKED";
1748  if (static_cast<bool>(feature._status & FeatureFloodCount::Status::DIRTY))
1749  out << " DIRTY";
1750  if (static_cast<bool>(feature._status & FeatureFloodCount::Status::INACTIVE))
1751  out << " INACTIVE";
1752 
1753  if (debug)
1754  {
1755  out << "\nOrig IDs (rank, index): ";
1756  for (const auto & orig_pair : feature._orig_ids)
1757  out << '(' << orig_pair.first << ", " << orig_pair.second << ") ";
1758  out << "\nVar_index: " << feature._var_index;
1759  out << "\nMin Entity ID: " << feature._min_entity_id;
1760  }
1761  out << "\n\n";
1762 
1763  return out;
1764 }
1765 
1766 /*****************************************************************************************
1767  ******************************* Utility Routines ****************************************
1768  *****************************************************************************************
1769  */
1770 void
1771 updateBBoxExtremesHelper(MeshTools::BoundingBox & bbox, const Point & node)
1772 {
1773  for (unsigned int i = 0; i < LIBMESH_DIM; ++i)
1774  {
1775  bbox.min()(i) = std::min(bbox.min()(i), node(i));
1776  bbox.max()(i) = std::max(bbox.max()(i), node(i));
1777  }
1778 }
1779 
1780 void
1781 updateBBoxExtremesHelper(MeshTools::BoundingBox & bbox, const Elem & elem)
1782 {
1783  for (auto node_n = decltype(elem.n_nodes())(0); node_n < elem.n_nodes(); ++node_n)
1784  updateBBoxExtremesHelper(bbox, *(elem.get_node(node_n)));
1785 }
1786 
1787 bool
1788 areElemListsMergeable(const std::list<dof_id_type> & elem_list1,
1789  const std::list<dof_id_type> & elem_list2,
1790  MeshBase & mesh)
1791 {
1792  for (const auto elem_id1 : elem_list1)
1793  {
1794  const auto * elem1 = mesh.elem_ptr(elem_id1);
1795  for (const auto elem_id2 : elem_list2)
1796  {
1797  const auto * elem2 = mesh.elem_ptr(elem_id2);
1798  if (elem1->has_neighbor(elem2))
1799  return true;
1800  }
1801  }
1802  return false;
1803 }
1804 
1805 // Constants
1806 const std::size_t FeatureFloodCount::invalid_size_t = std::numeric_limits<std::size_t>::max();
1807 const unsigned int FeatureFloodCount::invalid_id = std::numeric_limits<unsigned int>::max();
FeatureFloodCount(const InputParameters &parameters)
void visitNeighborsHelper(const T *curr_entity, std::vector< const T * > neighbor_entities, std::size_t current_index, FeatureData *feature, bool expand_halos_only, bool topological_neighbor, bool disjoint_only)
The actual logic for visiting neighbors is abstracted out here.
virtual Real getEntityValue(dof_id_type entity_id, FieldType field_type, std::size_t var_index=0) const
Point _centroid
The centroid of the feature (average of coordinates from entities participating in the volume calcula...
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 ...
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...
void expandPointHalos()
This method takes all of the partial features and expands the local, ghosted, and halo sets around th...
bool boundingBoxesIntersect(const FeatureData &rhs) const
Determines if any of this FeatureData&#39;s bounding boxes overlap with the other FeatureData&#39;s bounding ...
virtual Real getValue() override
bool halosIntersect(const FeatureData &rhs) const
Determine if one of this FeaturesData&#39;s member sets intersects the other FeatureData&#39;s corresponding ...
const std::size_t _n_vars
void visitElementalNeighbors(const Elem *elem, std::size_t current_index, FeatureData *feature, bool expand_halos_only, bool disjoint_only)
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...
void serialize(std::string &serialized_buffer)
These routines packs/unpack the _feature_map data into a structure suitable for parallel communicatio...
std::set< dof_id_type > _local_ids
Holds the local ids in the interior of a feature.
virtual bool areFeaturesMergeable(const FeatureData &f1, const FeatureData &f2) const
Method for determining whether two features are mergeable.
bool areElemListsMergeable(const std::list< dof_id_type > &elem_list1, const std::list< dof_id_type > &elem_list2, MeshBase &mesh)
const bool _condense_map_info
std::set< dof_id_type > _periodic_nodes
Holds the nodes that belong to the feature on a periodic boundary.
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.
void updateBBoxExtremes(MeshBase &mesh)
Update the minimum and maximum coordinates of a bounding box given a Point, Elem or BBox parameter...
void deserialize(std::vector< std::string > &serialized_buffers)
This routine takes the vector of byte buffers (one for each processor), deserializes them into a seri...
virtual void finalize() override
bool ghostedIntersect(const FeatureData &rhs) const
std::vector< MeshTools::BoundingBox > _bboxes
The vector of bounding boxes completely enclosing this feature (multiple used with periodic constrain...
std::map< dof_id_type, std::vector< unsigned int > > _entity_var_to_features
std::set< dof_id_type > _disjoint_halo_ids
Holds halo ids that extend onto a non-topologically connected surface.
void communicateAndMerge()
This routine handles all of the serialization, communication and deserialization of the data structur...
virtual void initialize() override
bool mergeable(const FeatureData &rhs) const
The routine called to see if two features are mergeable:
virtual Real getThreshold(std::size_t current_index) const
Return the starting comparison threshold to use when inspecting an entity during the flood stage...
void expandBBox(const FeatureData &rhs)
Located the overlapping bounding box between this Feature and the other Feature and expands that over...
virtual unsigned int getFeatureVar(unsigned int feature_id) const
Returns the variable representing the passed in feature.
void visitNodalNeighbors(const Node *node, std::size_t current_index, FeatureData *feature, bool expand_halos_only)
These two routines are utility routines used by the flood routine and by derived classes for visiting...
std::map< dof_id_type, int > _ghosted_entity_ids
The map for holding reconstructed ghosted element information.
void appendPeriodicNeighborNodes(FeatureData &feature) const
This routine adds the periodic node information to our data structure prior to packing the data this ...
void sortAndLabel()
Sort and assign ids to features based on their position in the container after sorting.
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.
const Real _connecting_threshold
The threshold above (or below) which neighboring entities are flooded (where regions can be extended ...
virtual std::size_t getTotalFeatureCount() const
Returns the total feature count (active and inactive ids, useful for sizing vectors) ...
virtual void execute() override
bool _is_elemental
Determines if the flood counter is elements or not (nodes)
void dataStore(std::ostream &stream, FeatureFloodCount::FeatureData &feature, void *context)
bool _is_master
Convenience variable for testing master rank.
std::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...
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...
std::vector< FeatureData > _feature_sets
The data structure used to hold the globally unique features.
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 ...
unsigned long _var_number
This variable is used to build the periodic node map.
std::vector< std::vector< const Elem * > > _nodes_to_elem_map
The data structure used to find neighboring elements give a node ID.
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.
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)
static const unsigned int invalid_id
InputParameters validParams< FeatureFloodCount >()
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...
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) ...
bool periodicBoundariesIntersect(const FeatureData &rhs) const
const Real _threshold
The threshold above (or below) where an entity may begin a new region (feature)
void updateBBoxExtremesHelper(MeshTools::BoundingBox &bbox, const Point &node)
std::unique_ptr< PointLocatorBase > _point_locator
void mergeSets()
This routine is called on the master rank only and stitches together the partial feature pieces seen ...
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.
void merge(FeatureData &&rhs)
Merges another Feature Data into this one.
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.
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...
virtual Real getConnectingThreshold(std::size_t current_index) const
Return the "connecting" comparison threshold to use when inspecting an entity during the flood stage...
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.
std::vector< unsigned int > _empty_var_to_features
const processor_id_type _n_procs
Convenience variable holding the number of processors in this simulation.
virtual void meshChanged() override
std::vector< unsigned int > _feature_counts_per_map
The number of features seen by this object per map.
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) ...
std::set< dof_id_type > _halo_ids
Holds the ids surrounding the feature.
std::list< std::pair< processor_id_type, unsigned int > > _orig_ids
Original processor/local ids.
virtual bool doesFeatureIntersectBoundary(unsigned int feature_id) const
Returns a Boolean indicating whether this feature intersects any boundary.
virtual void initialSetup() override
Status _status
The status of a feature (used mostly in derived classes like the GrainTracker)
bool flood(const DofObject *dof_object, std::size_t current_index, FeatureData *feature)
This method will "mark" all entities on neighboring elements that are above the supplied threshold...
bool _intersects_boundary
Flag indicating whether this feature intersects a boundary.
MooseMesh & _mesh
A reference to the mesh.
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...
const bool _compute_var_to_feature_map
Indicates whether or not the var to feature map is populated.
std::set< dof_id_type > _ghosted_ids
Holds the ghosted ids for a feature (the ids which will be used for stitching.
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.
std::vector< std::map< dof_id_type, int > > _var_index_maps
This map keeps track of which variables own which nodes.
std::vector< std::size_t > _local_to_global_feature_map
The vector recording the local to global feature indices.
void prepareDataForTransfer()
This routine uses the local flooded data to build up the local feature data structures (_feature_sets...
static bool setsIntersect(InputIterator first1, InputIterator last1, InputIterator first2, InputIterator last2)
This method detects whether two sets intersect without building a result set.