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