www.mooseframework.org
PolycrystalICTools.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 "PolycrystalICTools.h"
9 
10 // MOOSE includes
11 #include "MooseMesh.h"
12 #include "MooseVariable.h"
13 
14 #include "libmesh/mesh_tools.h"
15 #include "libmesh/periodic_boundaries.h"
16 #include "libmesh/point_locator_base.h"
17 
18 namespace GraphColoring
19 {
20 const unsigned int INVALID_COLOR = std::numeric_limits<unsigned int>::max();
21 }
22 
23 namespace PolycrystalICTools
24 {
25 const unsigned int HALO_THICKNESS = 4;
26 }
27 
28 // Forward declarations
29 bool colorGraph(const PolycrystalICTools::AdjacencyMatrix<Real> & adjacency_matrix,
30  std::vector<unsigned int> & colors,
31  unsigned int n_vertices,
32  unsigned int n_ops,
33  unsigned int vertex);
34 bool isGraphValid(const PolycrystalICTools::AdjacencyMatrix<Real> & adjacency_matrix,
35  std::vector<unsigned int> & colors,
36  unsigned int n_vertices,
37  unsigned int vertex,
38  unsigned int color);
39 void visitElementalNeighbors(const Elem * elem,
40  const MeshBase & mesh,
41  const PointLocatorBase & point_locator,
42  const PeriodicBoundaries * pb,
43  std::set<dof_id_type> & halo_ids);
44 
45 std::vector<unsigned int>
46 PolycrystalICTools::assignPointsToVariables(const std::vector<Point> & centerpoints,
47  const Real op_num,
48  const MooseMesh & mesh,
49  const MooseVariable & var)
50 {
51  Real grain_num = centerpoints.size();
52 
53  std::vector<unsigned int> assigned_op(grain_num);
54  std::vector<int> min_op_ind(op_num);
55  std::vector<Real> min_op_dist(op_num);
56 
57  // Assign grains to specific order parameters in a way that maximizes the distance
58  for (unsigned int grain = 0; grain < grain_num; grain++)
59  {
60  // Determine the distance to the closest center assigned to each order parameter
61  if (grain >= op_num)
62  {
63  // We can set the array to the distances to the grains 0..op_num-1 (see assignment in the else
64  // case)
65  for (unsigned int i = 0; i < op_num; ++i)
66  {
67  min_op_dist[i] =
68  mesh.minPeriodicDistance(var.number(), centerpoints[grain], centerpoints[i]);
69  min_op_ind[assigned_op[i]] = i;
70  }
71 
72  // Now check if any of the extra grains are even closer
73  for (unsigned int i = op_num; i < grain; ++i)
74  {
75  Real dist = mesh.minPeriodicDistance(var.number(), centerpoints[grain], centerpoints[i]);
76  if (min_op_dist[assigned_op[i]] > dist)
77  {
78  min_op_dist[assigned_op[i]] = dist;
79  min_op_ind[assigned_op[i]] = i;
80  }
81  }
82  }
83  else
84  {
85  assigned_op[grain] = grain;
86  continue;
87  }
88 
89  // Assign the current center point to the order parameter that is furthest away.
90  unsigned int mx_ind = 0;
91  for (unsigned int i = 1; i < op_num; ++i) // Find index of max
92  if (min_op_dist[mx_ind] < min_op_dist[i])
93  mx_ind = i;
94 
95  assigned_op[grain] = mx_ind;
96  }
97 
98  return assigned_op;
99 }
100 
101 unsigned int
103  const std::vector<Point> & centerpoints,
104  const MooseMesh & mesh,
105  const MooseVariable & var,
106  const Real maxsize)
107 {
108  unsigned int grain_num = centerpoints.size();
109 
110  Real min_distance = maxsize;
111  unsigned int min_index = grain_num;
112  // Loops through all of the grain centers and finds the center that is closest to the point p
113  for (unsigned int grain = 0; grain < grain_num; grain++)
114  {
115  Real distance = mesh.minPeriodicDistance(var.number(), centerpoints[grain], p);
116 
117  if (min_distance > distance)
118  {
119  min_distance = distance;
120  min_index = grain;
121  }
122  }
123 
124  if (min_index >= grain_num)
125  mooseError("ERROR in PolycrystalVoronoiVoidIC: didn't find minimum values in grain_value_calc");
126 
127  return min_index;
128 }
129 
132  const std::map<dof_id_type, unsigned int> & entity_to_grain,
133  MooseMesh & mesh,
134  const PeriodicBoundaries * pb,
135  unsigned int n_grains,
136  bool is_elemental)
137 {
138  if (is_elemental)
139  return buildElementalGrainAdjacencyMatrix(entity_to_grain, mesh, pb, n_grains);
140  else
141  return buildNodalGrainAdjacencyMatrix(entity_to_grain, mesh, pb, n_grains);
142 }
143 
146  const std::map<dof_id_type, unsigned int> & element_to_grain,
147  MooseMesh & mesh,
148  const PeriodicBoundaries * pb,
149  unsigned int n_grains)
150 {
151  AdjacencyMatrix<Real> adjacency_matrix(n_grains);
152 
153  // We can't call this in the constructor, it appears that _mesh_ptr is always NULL there.
154  mesh.errorIfDistributedMesh("advanced_op_assignment = true");
155 
156  std::vector<const Elem *> all_active_neighbors;
157 
158  std::vector<std::set<dof_id_type>> local_ids(n_grains);
159  std::vector<std::set<dof_id_type>> halo_ids(n_grains);
160 
161  std::unique_ptr<PointLocatorBase> point_locator = mesh.getMesh().sub_point_locator();
162  const auto end = mesh.getMesh().active_elements_end();
163  for (auto el = mesh.getMesh().active_elements_begin(); el != end; ++el)
164  {
165  const Elem * elem = *el;
166  std::map<dof_id_type, unsigned int>::const_iterator grain_it =
167  element_to_grain.find(elem->id());
168  mooseAssert(grain_it != element_to_grain.end(), "Element not found in map");
169  unsigned int my_grain = grain_it->second;
170 
171  all_active_neighbors.clear();
172  // Loop over all neighbors (at the the same level as the current element)
173  for (unsigned int i = 0; i < elem->n_neighbors(); ++i)
174  {
175  const Elem * neighbor_ancestor = elem->topological_neighbor(i, mesh, *point_locator, pb);
176  if (neighbor_ancestor)
177  // Retrieve only the active neighbors for each side of this element, append them to the list
178  // of active neighbors
179  neighbor_ancestor->active_family_tree_by_topological_neighbor(
180  all_active_neighbors, elem, mesh, *point_locator, pb, false);
181  }
182 
183  // Loop over all active element neighbors
184  for (std::vector<const Elem *>::const_iterator neighbor_it = all_active_neighbors.begin();
185  neighbor_it != all_active_neighbors.end();
186  ++neighbor_it)
187  {
188  const Elem * neighbor = *neighbor_it;
189  std::map<dof_id_type, unsigned int>::const_iterator grain_it2 =
190  element_to_grain.find(neighbor->id());
191  mooseAssert(grain_it2 != element_to_grain.end(), "Element not found in map");
192  unsigned int their_grain = grain_it2->second;
193 
194  if (my_grain != their_grain)
195  {
204  // First add corresponding element and grain information
205  local_ids[my_grain].insert(elem->id());
206  local_ids[their_grain].insert(neighbor->id());
207 
208  // Now add opposing element and grain information
209  halo_ids[my_grain].insert(neighbor->id());
210  halo_ids[their_grain].insert(elem->id());
211  }
212  // adjacency_matrix[my_grain][their_grain] = 1;
213  }
214  }
215 
216  // Build up the halos
217  std::set<dof_id_type> set_difference;
218  for (unsigned int i = 0; i < n_grains; ++i)
219  {
220  std::set<dof_id_type> orig_halo_ids(halo_ids[i]);
221 
222  for (unsigned int halo_level = 0; halo_level < PolycrystalICTools::HALO_THICKNESS; ++halo_level)
223  {
224  for (std::set<dof_id_type>::iterator entity_it = orig_halo_ids.begin();
225  entity_it != orig_halo_ids.end();
226  ++entity_it)
227  {
228  if (true)
230  mesh.elemPtr(*entity_it), mesh.getMesh(), *point_locator, pb, halo_ids[i]);
231  else
232  mooseError("Unimplemented");
233  }
234 
235  set_difference.clear();
236  std::set_difference(
237  halo_ids[i].begin(),
238  halo_ids[i].end(),
239  local_ids[i].begin(),
240  local_ids[i].end(),
241  std::insert_iterator<std::set<dof_id_type>>(set_difference, set_difference.begin()));
242 
243  halo_ids[i].swap(set_difference);
244  }
245  }
246 
247  // Finally look at the halo intersections to build the connectivity graph
248  std::set<dof_id_type> set_intersection;
249  for (unsigned int i = 0; i < n_grains; ++i)
250  for (unsigned int j = i + 1; j < n_grains; ++j)
251  {
252  set_intersection.clear();
253  std::set_intersection(
254  halo_ids[i].begin(),
255  halo_ids[i].end(),
256  halo_ids[j].begin(),
257  halo_ids[j].end(),
258  std::insert_iterator<std::set<dof_id_type>>(set_intersection, set_intersection.begin()));
259 
260  if (!set_intersection.empty())
261  {
262  adjacency_matrix(i, j) = 1.;
263  adjacency_matrix(j, i) = 1.;
264  }
265  }
266 
267  return adjacency_matrix;
268 }
269 
272  const std::map<dof_id_type, unsigned int> & node_to_grain,
273  MooseMesh & mesh,
274  const PeriodicBoundaries * /*pb*/,
275  unsigned int n_grains)
276 {
277  // Build node to elem map
278  std::vector<std::vector<const Elem *>> nodes_to_elem_map;
279  MeshTools::build_nodes_to_elem_map(mesh.getMesh(), nodes_to_elem_map);
280 
281  AdjacencyMatrix<Real> adjacency_matrix(n_grains);
282 
283  const auto end = mesh.getMesh().active_nodes_end();
284  for (auto nl = mesh.getMesh().active_nodes_begin(); nl != end; ++nl)
285  {
286  const Node * node = *nl;
287  std::map<dof_id_type, unsigned int>::const_iterator grain_it = node_to_grain.find(node->id());
288  mooseAssert(grain_it != node_to_grain.end(), "Node not found in map");
289  unsigned int my_grain = grain_it->second;
290 
291  std::vector<const Node *> nodal_neighbors;
292  MeshTools::find_nodal_neighbors(mesh.getMesh(), *node, nodes_to_elem_map, nodal_neighbors);
293 
294  // Loop over all nodal neighbors
295  for (unsigned int i = 0; i < nodal_neighbors.size(); ++i)
296  {
297  const Node * neighbor_node = nodal_neighbors[i];
298  std::map<dof_id_type, unsigned int>::const_iterator grain_it2 =
299  node_to_grain.find(neighbor_node->id());
300  mooseAssert(grain_it2 != node_to_grain.end(), "Node not found in map");
301  unsigned int their_grain = grain_it2->second;
302 
303  if (my_grain != their_grain)
311  adjacency_matrix(my_grain, their_grain) = 1.;
312  }
313  }
314 
315  return adjacency_matrix;
316 }
317 
318 std::vector<unsigned int>
320  unsigned int n_grains,
321  unsigned int n_ops,
322  const MooseEnum & coloring_algorithm)
323 {
324  Moose::perf_log.push("assignOpsToGrains()", "PolycrystalICTools");
325 
326  std::vector<unsigned int> grain_to_op(n_grains, GraphColoring::INVALID_COLOR);
327 
328  // Use a simple backtracking coloring algorithm
329  if (coloring_algorithm == "bt")
330  {
331  if (!colorGraph(adjacency_matrix, grain_to_op, n_grains, n_ops, 0))
332  mooseError(
333  "Unable to find a valid Grain to op configuration, do you have enough op variables?");
334  }
335  else // PETSc Coloring algorithms
336  {
337 #ifdef LIBMESH_HAVE_PETSC
338  const std::string & ca_str = coloring_algorithm;
339  Real * am_data = adjacency_matrix.rawDataPtr();
340  Moose::PetscSupport::colorAdjacencyMatrix(
341  am_data, n_grains, n_ops, grain_to_op, ca_str.c_str());
342 #else
343  mooseError("Selected coloring algorithm requires PETSc");
344 #endif
345  }
346 
347  Moose::perf_log.pop("assignOpsToGrains()", "PolycrystalICTools");
348 
349  return grain_to_op;
350 }
351 
352 MooseEnum
354 {
355  return MooseEnum("legacy bt jp power greedy", "legacy");
356 }
357 
358 std::string
360 {
361  return "The grain neighbor graph coloring algorithm to use. \"legacy\" is the original "
362  "algorithm "
363  "which does not guarantee a valid coloring. \"bt\" is a simple backtracking algorithm "
364  "which will produce a valid coloring but has potential exponential run time. The "
365  "remaining algorithms require PETSc but are recommended for larger problems (See "
366  "MatColoringType)";
367 }
368 
372 void
373 visitElementalNeighbors(const Elem * elem,
374  const MeshBase & mesh,
375  const PointLocatorBase & point_locator,
376  const PeriodicBoundaries * pb,
377  std::set<dof_id_type> & halo_ids)
378 {
379  mooseAssert(elem, "Elem is NULL");
380 
381  std::vector<const Elem *> all_active_neighbors;
382 
383  // Loop over all neighbors (at the the same level as the current element)
384  for (unsigned int i = 0; i < elem->n_neighbors(); ++i)
385  {
386  const Elem * neighbor_ancestor = elem->topological_neighbor(i, mesh, point_locator, pb);
387  if (neighbor_ancestor)
388  // Retrieve only the active neighbors for each side of this element, append them to the list
389  // of active neighbors
390  neighbor_ancestor->active_family_tree_by_topological_neighbor(
391  all_active_neighbors, elem, mesh, point_locator, pb, false);
392  }
393 
394  // Loop over all active element neighbors
395  for (std::vector<const Elem *>::const_iterator neighbor_it = all_active_neighbors.begin();
396  neighbor_it != all_active_neighbors.end();
397  ++neighbor_it)
398  if (*neighbor_it)
399  halo_ids.insert((*neighbor_it)->id());
400 }
401 
405 bool
407  std::vector<unsigned int> & colors,
408  unsigned int n_vertices,
409  unsigned int n_colors,
410  unsigned int vertex)
411 {
412  // Base case: All grains are assigned
413  if (vertex == n_vertices)
414  return true;
415 
416  // Consider this grain and try different ops
417  for (unsigned int color_idx = 0; color_idx < n_colors; ++color_idx)
418  {
419  // We'll try to spread these colors around a bit rather than
420  // packing them all on the first few colors if we have several colors.
421  unsigned int color = (vertex + color_idx) % n_colors;
422 
423  if (isGraphValid(adjacency_matrix, colors, n_vertices, vertex, color))
424  {
425  colors[vertex] = color;
426 
427  if (colorGraph(adjacency_matrix, colors, n_vertices, n_colors, vertex + 1))
428  return true;
429 
430  // Backtrack...
431  colors[vertex] = GraphColoring::INVALID_COLOR;
432  }
433  }
434 
435  return false;
436 }
437 
438 bool
440  std::vector<unsigned int> & colors,
441  unsigned int n_vertices,
442  unsigned int vertex,
443  unsigned int color)
444 {
445  // See if the proposed color is valid based on the current neighbor colors
446  for (unsigned int neighbor = 0; neighbor < n_vertices; ++neighbor)
447  if (adjacency_matrix(vertex, neighbor) && color == colors[neighbor])
448  return false;
449  return true;
450 }
const unsigned int HALO_THICKNESS
AdjacencyMatrix< Real > buildGrainAdjacencyMatrix(const std::map< dof_id_type, unsigned int > &entity_to_grain, MooseMesh &mesh, const PeriodicBoundaries *pb, unsigned int n_grains, bool is_elemental)
bool isGraphValid(const PolycrystalICTools::AdjacencyMatrix< Real > &adjacency_matrix, std::vector< unsigned int > &colors, unsigned int n_vertices, unsigned int vertex, unsigned int color)
const unsigned int INVALID_COLOR
void visitElementalNeighbors(const Elem *elem, const MeshBase &mesh, const PointLocatorBase &point_locator, const PeriodicBoundaries *pb, std::set< dof_id_type > &halo_ids)
Utility routines.
bool colorGraph(const PolycrystalICTools::AdjacencyMatrix< Real > &adjacency_matrix, std::vector< unsigned int > &colors, unsigned int n_vertices, unsigned int n_ops, unsigned int vertex)
Backtracking graph coloring routines.
unsigned int assignPointToGrain(const Point &p, const std::vector< Point > &centerpoints, const MooseMesh &mesh, const MooseVariable &var, const Real maxsize)
std::vector< unsigned int > assignPointsToVariables(const std::vector< Point > &centerpoints, const Real op_num, const MooseMesh &mesh, const MooseVariable &var)
AdjacencyMatrix< Real > buildNodalGrainAdjacencyMatrix(const std::map< dof_id_type, unsigned int > &node_to_grain, MooseMesh &mesh, const PeriodicBoundaries *pb, unsigned int n_grains)
MooseEnum coloringAlgorithms()
AdjacencyMatrix< Real > buildElementalGrainAdjacencyMatrix(const std::map< dof_id_type, unsigned int > &element_to_grain, MooseMesh &mesh, const PeriodicBoundaries *pb, unsigned int n_grains)
Simple 2D block matrix indicating graph adjacency.
std::vector< unsigned int > assignOpsToGrains(AdjacencyMatrix< Real > &adjacency_matrix, unsigned int n_grains, unsigned int n_ops, const MooseEnum &coloring_algorithm)
std::string coloringAlgorithmDescriptions()