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