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 "DenseMatrix.h"
11 : #include "PolycrystalUserObjectBase.h"
12 : #include "NonlinearSystemBase.h"
13 : #include "MooseMesh.h"
14 : #include "MooseVariable.h"
15 :
16 : #include <vector>
17 : #include <map>
18 : #include <algorithm>
19 :
20 : InputParameters
21 1130 : PolycrystalUserObjectBase::validParams()
22 : {
23 1130 : InputParameters params = FeatureFloodCount::validParams();
24 1130 : params.addClassDescription("This object provides the base capability for creating proper reduced "
25 : "order parameter polycrystal initial conditions.");
26 2260 : params.addRequiredCoupledVarWithAutoBuild(
27 : "variable", "var_name_base", "op_num", "Array of coupled variables");
28 2260 : params.addParam<bool>("output_adjacency_matrix",
29 2260 : false,
30 : "Output the Grain Adjacency Matrix used in the coloring algorithms. "
31 : "Additionally, the grain to OP assignments will be printed");
32 1130 : params.addParam<MooseEnum>("coloring_algorithm",
33 2260 : PolycrystalUserObjectBase::coloringAlgorithms(),
34 2260 : PolycrystalUserObjectBase::coloringAlgorithmDescriptions());
35 :
36 : // FeatureFloodCount adds a relationship manager, but we need to extend that for PolycrystalIC
37 : params.clearRelationshipManagers();
38 :
39 2260 : params.addRelationshipManager(
40 : "ElementSideNeighborLayers",
41 : Moose::RelationshipManagerType::GEOMETRIC,
42 :
43 1695 : [](const InputParameters & /*obj_params*/, InputParameters & rm_params)
44 1695 : { rm_params.set<unsigned short>("layers") = 2; }
45 :
46 : );
47 :
48 2260 : params.addRelationshipManager("ElementSideNeighborLayers",
49 : Moose::RelationshipManagerType::ALGEBRAIC);
50 :
51 : // Hide the output of the IC objects by default, it doesn't change over time
52 3390 : params.set<std::vector<OutputName>>("outputs") = {"none"};
53 :
54 : /// Run this user object more than once on the initial condition to handle initial adaptivity
55 1130 : params.set<bool>("allow_duplicate_execution_on_initial") = true;
56 :
57 : // This object should only be executed _before_ the initial condition
58 1130 : ExecFlagEnum execute_options = MooseUtils::getDefaultExecFlagEnum();
59 1130 : execute_options = EXEC_INITIAL;
60 2260 : params.set<ExecFlagEnum>("execute_on") = execute_options;
61 :
62 1130 : return params;
63 0 : }
64 :
65 565 : PolycrystalUserObjectBase::PolycrystalUserObjectBase(const InputParameters & parameters)
66 : : FeatureFloodCount(parameters),
67 565 : _dim(_mesh.dimension()),
68 565 : _op_num(_vars.size()),
69 1130 : _coloring_algorithm(getParam<MooseEnum>("coloring_algorithm")),
70 565 : _colors_assigned(false),
71 1130 : _output_adjacency_matrix(getParam<bool>("output_adjacency_matrix")),
72 1130 : _num_chunks(FeatureFloodCount::invalid_proc_id)
73 : {
74 : mooseAssert(_single_map_mode, "Do not turn off single_map_mode with this class");
75 565 : }
76 :
77 : void
78 510 : PolycrystalUserObjectBase::initialSetup()
79 : {
80 : /**
81 : * For polycrystal ICs we need to assume that each of the variables has the same periodicity.
82 : * Since BCs are handled elsewhere in the system, we'll have to check this case explicitly.
83 : */
84 510 : if (_op_num < 1)
85 0 : mooseError("No coupled variables found");
86 :
87 510 : const auto first_variable_value = _mesh.queryPeriodicDimensions(*_vars[0]);
88 3292 : for (unsigned int i = 1; i < _vars.size(); ++i)
89 2782 : if (_mesh.queryPeriodicDimensions(*_vars[i]) != first_variable_value)
90 0 : mooseError("Coupled polycrystal variables differ in periodicity");
91 :
92 510 : FeatureFloodCount::initialSetup();
93 510 : }
94 :
95 : void
96 595 : PolycrystalUserObjectBase::initialize()
97 : {
98 595 : if (_colors_assigned && !_fe_problem.hasInitialAdaptivity())
99 : return;
100 :
101 : _entity_to_grain_cache.clear();
102 :
103 424 : FeatureFloodCount::initialize();
104 : }
105 :
106 : void
107 595 : PolycrystalUserObjectBase::execute()
108 : {
109 595 : if (!_colors_assigned)
110 406 : precomputeGrainStructure();
111 : // No need to rerun the object if the mesh hasn't changed
112 378 : else if (!_fe_problem.hasInitialAdaptivity())
113 171 : return;
114 :
115 848 : TIME_SECTION("execute", 2, "Computing Polycrystal Initial Condition");
116 :
117 : /**
118 : * We need one map per grain when creating the initial condition to support overlapping features.
119 : * Luckily, this is a fairly sparse structure.
120 : */
121 424 : _entities_visited.resize(getNumGrains());
122 :
123 : /**
124 : * This loop is similar to the one found in the base class however, there are two key differences
125 : * between building up the initial condition and discovering features based on solution variables:
126 : *
127 : * 1) When building up the initial condition, we aren't inspecting the actual variable values so
128 : * we don't need to loop over all of the coupled variables.
129 : * 2) We want to discover all features on a single pass since there may be thousands of features
130 : * in a simulation. However, we can only actively flood a single feature at a time. To make
131 : * sure that we pick up all features that might start on a given entity, we'll keep retrying
132 : * the flood routine on the same entity as long as new discoveries are being made. We know
133 : * this information from the return value of flood.
134 : */
135 478045 : for (const auto & current_elem : _fe_problem.getNonlinearEvaluableElementRange())
136 : {
137 : // Loop over elements or nodes
138 477621 : if (_is_elemental)
139 485253 : while (flood(current_elem, invalid_size_t))
140 : ;
141 : else
142 : {
143 0 : auto n_nodes = current_elem->n_vertices();
144 0 : for (auto i = decltype(n_nodes)(0); i < n_nodes; ++i)
145 : {
146 0 : const Node * current_node = current_elem->node_ptr(i);
147 :
148 0 : while (flood(current_node, invalid_size_t))
149 : ;
150 : }
151 : }
152 : }
153 424 : }
154 :
155 : processor_id_type
156 424 : PolycrystalUserObjectBase::numberOfDistributedMergeHelpers() const
157 : {
158 : mooseAssert(_num_chunks != FeatureFloodCount::invalid_proc_id,
159 : "prepareDataForTransfer() hasn't been called yet");
160 :
161 424 : return _num_chunks;
162 : }
163 :
164 : void
165 424 : PolycrystalUserObjectBase::prepareDataForTransfer()
166 : {
167 424 : FeatureFloodCount::prepareDataForTransfer();
168 :
169 : /**
170 : * With this class, all of the partial features are crammed into a single "outer" entry in our
171 : * data structure because we don't know the variable assignments during the initial condition
172 : * (that is the whole point of this class!). However, what we do know is _the_ feature number
173 : * that each piece belongs to, which is even more useful for merging. However, with extensive
174 : * testing on larger problems, even ordering optimally ordering these pieces is far too much
175 : * work for a single core to process. In order create a scalable algorithm, we'll first order
176 : * our data, then break it into complete chunks for several processors to work on concurrently.
177 : *
178 : * After sorting the data structure, it'll look something like this on some rank:
179 : * _partial_feature_sets (showing the unique_id and fake variable number):
180 : * 0 0 0 0 1 1 1 3 4 4 4 4 4 6 6 ...
181 : *
182 : * We may very well have gaps in our numbering. This is a local view on one rank.
183 : * We'd like to transform the data into something like the following
184 : * 0 0 0 0 1 1 1 <- First 3 items here (even when we don't have all features)
185 : * 1 3 4 4 4 4 4 <- Next 3 items here
186 : * 2 6 6 ... <- Next 3 or remainder (linear partitioning)
187 : *
188 : * The way to break up this work is to simply break it into min(n_features, n_procs) chunks.
189 : * e.g. If we have 70 features and 8 cores, we'll partition the work into 8 chunks. In the
190 : * odd case where we have more processors than features (e.g. 100 processors merging 10 features),
191 : * we'll use end up using a subset of the available cores.
192 : *
193 : * To get all this started we just need the number of features, but we don't normally know
194 : * that until we've merged everything together... sigh... Wait! We can fall back on our
195 : * sorted data structure though and figure out before we sort and merge. We'll need
196 : * one more parallel communication, that won't hurt anything, right?!
197 : */
198 424 : _partial_feature_sets[0].sort();
199 :
200 : // Get the largest ID seen on any rank
201 424 : auto largest_id = _partial_feature_sets[0].back()._id;
202 424 : _communicator.max(largest_id);
203 : mooseAssert(largest_id != invalid_size_t, "Largest ID should not be invalid");
204 :
205 : /**
206 : * With this class there are no guarentees that our IDs have a contiguous zero-based numbering.
207 : * However, for many of the common derived classes they do (generated grain structures).
208 : * If we have holes in our numbering, we might not get an even partition, but it shouldn't break.
209 : * We just need the best guess at a total number (before we can actually count) which should
210 : * be bounded by the largest_id + 1.
211 : */
212 424 : auto total_items = largest_id + 1;
213 :
214 424 : _num_chunks = std::min(_app.n_processors(), total_items);
215 :
216 : /**
217 : * Here we are resizing our data structures that we normally size upon construction. This is to
218 : * support the parallel merge capability that's in the FeatureFloodCount class. We'll need to undo
219 : * this later, there are a few assumptions built on the sizes of these data structures.
220 : *
221 : * See FeatureFloodCount::consolidateMergedFeatures for the "un-sizing" of these structures.
222 : */
223 424 : _partial_feature_sets.resize(_num_chunks);
224 424 : _feature_counts_per_map.resize(_num_chunks);
225 :
226 8056 : for (auto it = _partial_feature_sets[0].begin(); it != _partial_feature_sets[0].end();
227 : /* No increment on it*/)
228 : {
229 7632 : auto chunk = MooseUtils::linearPartitionChunk(total_items, _num_chunks, it->_id);
230 :
231 7632 : if (chunk)
232 : {
233 2501 : _partial_feature_sets[chunk].emplace_back(std::move(*it));
234 : it = _partial_feature_sets[0].erase(it); // it is incremented here!
235 : }
236 : else
237 : ++it;
238 : }
239 424 : }
240 :
241 : void
242 595 : PolycrystalUserObjectBase::finalize()
243 : {
244 595 : if (_colors_assigned && !_fe_problem.hasInitialAdaptivity())
245 171 : return;
246 :
247 848 : TIME_SECTION("finalize", 2, "Finalizing Polycrystal Initial Condition");
248 :
249 : // TODO: Possibly retrieve the halo thickness from the active GrainTracker object?
250 : constexpr unsigned int halo_thickness = 2;
251 :
252 424 : expandEdgeHalos(halo_thickness - 1);
253 :
254 424 : FeatureFloodCount::finalize();
255 :
256 424 : if (!_colors_assigned)
257 : {
258 : // Resize the color assignment vector here. All ranks need a copy of this
259 406 : _grain_idx_to_op.resize(_feature_count, PolycrystalUserObjectBase::INVALID_COLOR);
260 406 : if (_is_primary)
261 : {
262 273 : buildGrainAdjacencyMatrix();
263 :
264 273 : assignOpsToGrains();
265 :
266 269 : if (_output_adjacency_matrix)
267 103 : printGrainAdjacencyMatrix();
268 : }
269 :
270 : // Communicate the coloring map with all ranks
271 402 : _communicator.broadcast(_grain_to_op);
272 :
273 : /**
274 : * All ranks: Update the variable indices based on the graph coloring algorithm.
275 : */
276 5890 : for (auto & feature : _feature_sets)
277 5488 : feature._var_index = _grain_to_op.at(feature._id);
278 : }
279 :
280 420 : _colors_assigned = true;
281 420 : }
282 :
283 : void
284 424 : PolycrystalUserObjectBase::mergeSets()
285 : {
286 : // When working with _distribute_merge_work all of the maps will be empty except for one
287 1508 : for (const auto map_num : index_range(_partial_feature_sets))
288 : {
289 : /**
290 : * With initial conditions we know the grain IDs of every grain (even partial grains). We can
291 : * use this information to put all mergeable features adjacent to one and other in the list so
292 : * that merging is simply O(n).
293 : */
294 1084 : _partial_feature_sets[map_num].sort();
295 :
296 : auto it1 = _partial_feature_sets[map_num].begin();
297 : auto it_end = _partial_feature_sets[map_num].end();
298 8292 : while (it1 != it_end)
299 : {
300 : auto it2 = it1;
301 7632 : if (++it2 == it_end)
302 : break;
303 :
304 7208 : if (areFeaturesMergeable(*it1, *it2))
305 : {
306 4144 : it1->merge(std::move(*it2));
307 : _partial_feature_sets[map_num].erase(it2);
308 : }
309 : else
310 : ++it1; // Only increment if we have a mismatch
311 : }
312 : }
313 424 : }
314 :
315 : void
316 136 : PolycrystalUserObjectBase::restoreOriginalDataStructures(std::vector<std::list<FeatureData>> & orig)
317 : {
318 : // Move all the data back into the first list
319 : auto & master_list = orig[0];
320 660 : for (MooseIndex(_maps_size) map_num = 1; map_num < orig.size(); ++map_num)
321 : {
322 524 : master_list.splice(master_list.end(), orig[map_num]);
323 : orig[map_num].clear();
324 : }
325 :
326 136 : orig.resize(1);
327 136 : }
328 :
329 : bool
330 1252021 : PolycrystalUserObjectBase::isNewFeatureOrConnectedRegion(const DofObject * dof_object,
331 : std::size_t & current_index,
332 : FeatureData *& feature,
333 : Status & status,
334 : unsigned int & new_id)
335 : {
336 : mooseAssert(_t_step == 0, "PolyIC only works if we begin in the initial condition");
337 :
338 : // Retrieve the id of the current entity
339 1252021 : auto entity_id = dof_object->id();
340 : auto grains_it = _entity_to_grain_cache.lower_bound(entity_id);
341 :
342 1252021 : if (grains_it == _entity_to_grain_cache.end() || grains_it->first != entity_id)
343 : {
344 : std::vector<unsigned int> grain_ids;
345 :
346 386897 : if (_is_elemental)
347 386897 : getGrainsBasedOnElem(*static_cast<const Elem *>(dof_object), grain_ids);
348 : else
349 0 : getGrainsBasedOnPoint(*static_cast<const Node *>(dof_object), grain_ids);
350 :
351 386897 : grains_it = _entity_to_grain_cache.emplace_hint(grains_it, entity_id, std::move(grain_ids));
352 386897 : }
353 :
354 : /**
355 : * When building the IC, we can't use the _entities_visited data structure the same way as we do
356 : * for the base class. We need to discover multiple overlapping grains in a single pass. However
357 : * we don't know what grain we are working on when we enter the flood routine (when that check is
358 : * normally made). Only after we've made the callback to the child class do we know which grains
359 : * we are operating on (at least until we've triggered the recursion). We need to see if there
360 : * is at least one active grain where we haven't already visited the current entity before
361 : * continuing.
362 : */
363 1252021 : auto saved_grain_id = invalid_id;
364 1252021 : if (current_index == invalid_size_t)
365 : {
366 992450 : for (auto grain_id : grains_it->second)
367 : {
368 : auto map_it = _grain_to_op.find(grain_id);
369 : mooseAssert(!_colors_assigned || map_it != _grain_to_op.end(), "grain_id missing");
370 514829 : auto map_num = _colors_assigned ? map_it->second : grain_id;
371 :
372 514829 : if (_entities_visited[map_num].find(entity_id) == _entities_visited[map_num].end())
373 : {
374 : saved_grain_id = grain_id;
375 7632 : current_index = map_num;
376 : break;
377 : }
378 : }
379 :
380 485253 : if (current_index == invalid_size_t)
381 : return false;
382 : }
383 766768 : else if (_entities_visited[current_index].find(entity_id) !=
384 : _entities_visited[current_index].end())
385 : return false;
386 :
387 774400 : if (!feature)
388 : {
389 7632 : new_id = saved_grain_id;
390 : status &= ~Status::INACTIVE;
391 :
392 7632 : return true;
393 : }
394 : else
395 : {
396 : const auto & grain_ids = grains_it->second;
397 766768 : if (std::find(grain_ids.begin(), grain_ids.end(), feature->_id) != grain_ids.end())
398 : return true;
399 :
400 : /**
401 : * If we get here the current entity is not part of the active feature, however we now want to
402 : * look at neighbors.
403 : *
404 : */
405 267917 : if (_is_elemental)
406 : {
407 267917 : Elem * elem = _mesh.queryElemPtr(entity_id);
408 : mooseAssert(elem, "Element is nullptr");
409 :
410 : std::vector<const Elem *> all_active_neighbors;
411 267917 : MeshBase & mesh = _mesh.getMesh();
412 :
413 1349805 : for (auto i = decltype(elem->n_neighbors())(0); i < elem->n_neighbors(); ++i)
414 : {
415 : const Elem * neighbor_ancestor = nullptr;
416 :
417 : /**
418 : * Retrieve only the active neighbors for each side of this element, append them to the list
419 : * of active neighbors
420 : */
421 : neighbor_ancestor = elem->neighbor_ptr(i);
422 :
423 1081888 : if (neighbor_ancestor)
424 : {
425 : /**
426 : * In general, {evaluable elements} >= {local elements} U {algebraic ghosting elements}.
427 : * That is, the number of evaluable elements does NOT necessarily equal to the number of
428 : * local and algebraic ghosting elements. The neighbors of evaluable elements can be
429 : * remote even though we have two layers of geometric ghosting elements.
430 : */
431 1052672 : if (neighbor_ancestor->is_remote())
432 0 : continue;
433 :
434 1052672 : neighbor_ancestor->active_family_tree_by_neighbor(all_active_neighbors, elem, false);
435 : }
436 : else // if (expand_halos_only /*&& feature->_periodic_nodes.empty()*/)
437 : {
438 29216 : neighbor_ancestor = elem->topological_neighbor(i, mesh, *_point_locator, _pbs);
439 :
440 : /**
441 : * If the current element (passed into this method) doesn't have a connected neighbor but
442 : * does have a topological neighbor, this might be a new disjoint region that we'll
443 : * need to represent with a separate bounding box. To find out for sure, we'll need
444 : * see if the new neighbors are present in any of the halo or disjoint halo sets. If
445 : * they are not present, this is a new region.
446 : */
447 29216 : if (neighbor_ancestor)
448 : {
449 : /**
450 : * In general, {evaluable elements} >= {local elements} U {algebraic ghosting elements}.
451 : * That is, the number of evaluable elements does NOT necessarily equal to the number of
452 : * local and algebraic ghosting elements. The neighbors of evaluable elements can be
453 : * remote even though we have two layers of geometric ghosting elements.
454 : */
455 10682 : if (neighbor_ancestor->is_remote())
456 0 : continue;
457 :
458 10682 : neighbor_ancestor->active_family_tree_by_topological_neighbor(
459 10682 : all_active_neighbors, elem, mesh, *_point_locator, _pbs, false);
460 : }
461 : }
462 : }
463 :
464 1036685 : for (const auto neighbor : all_active_neighbors)
465 : {
466 : // Retrieve the id of the current entity
467 876251 : auto neighbor_id = neighbor->id();
468 : auto neighbor_it = _entity_to_grain_cache.lower_bound(neighbor_id);
469 :
470 876251 : if (neighbor_it == _entity_to_grain_cache.end() || neighbor_it->first != neighbor_id)
471 : {
472 : std::vector<unsigned int> more_grain_ids;
473 :
474 92211 : getGrainsBasedOnElem(*static_cast<const Elem *>(neighbor), more_grain_ids);
475 :
476 92211 : neighbor_it = _entity_to_grain_cache.emplace_hint(
477 : neighbor_it, neighbor_id, std::move(more_grain_ids));
478 92211 : }
479 :
480 : const auto & more_grain_ids = neighbor_it->second;
481 876251 : if (std::find(more_grain_ids.begin(), more_grain_ids.end(), feature->_id) !=
482 : more_grain_ids.end())
483 107483 : return true;
484 : }
485 267917 : }
486 :
487 160434 : return false;
488 : }
489 : }
490 :
491 : bool
492 7208 : PolycrystalUserObjectBase::areFeaturesMergeable(const FeatureData & f1,
493 : const FeatureData & f2) const
494 : {
495 7208 : if (f1._id != f2._id)
496 3064 : return false;
497 :
498 : mooseAssert(f1._var_index == f2._var_index, "Feature should be mergeable but aren't");
499 : return true;
500 : }
501 :
502 : void
503 273 : PolycrystalUserObjectBase::buildGrainAdjacencyMatrix()
504 : {
505 : mooseAssert(_is_primary, "This routine should only be called on the primary rank");
506 :
507 546 : _adjacency_matrix = std::make_unique<DenseMatrix<Real>>(_feature_count, _feature_count);
508 3641 : for (MooseIndex(_feature_sets) i = 0; i < _feature_sets.size(); ++i)
509 : {
510 38921 : for (MooseIndex(_feature_sets) j = i + 1; j < _feature_sets.size(); ++j)
511 : {
512 54176 : if (_feature_sets[i].boundingBoxesIntersect(_feature_sets[j]) &&
513 18623 : _feature_sets[i].halosIntersect(_feature_sets[j]))
514 : {
515 : // Our grain adjacency matrix is symmetrical
516 12398 : (*_adjacency_matrix)(i, j) = 1;
517 12398 : (*_adjacency_matrix)(j, i) = 1;
518 : }
519 : }
520 : }
521 273 : }
522 :
523 : void
524 273 : PolycrystalUserObjectBase::assignOpsToGrains()
525 : {
526 : mooseAssert(_is_primary, "This routine should only be called on the primary rank");
527 :
528 546 : TIME_SECTION("assignOpsToGrains", 2, "Assigning OPs to grains");
529 :
530 : // Use a simple backtracking coloring algorithm
531 273 : if (_coloring_algorithm == "bt")
532 : {
533 170 : paramInfo("coloring_algorithm",
534 : "The backtracking algorithm has exponential complexity. If you are using very few "
535 : "order parameters,\nor you have several hundred grains or more, you should use one "
536 : "of the PETSc coloring algorithms such as \"jp\".");
537 :
538 170 : if (!colorGraph(0))
539 2 : paramError("op_num",
540 : "Unable to find a valid grain to op coloring, Make sure you have created enough "
541 : "variables to hold a\nvalid polycrystal initial condition (no grains represented "
542 : "by the same variable should be allowed to\ntouch, ~8 for 2D, ~25 for 3D)?");
543 : }
544 : else // PETSc Coloring algorithms
545 : {
546 : const std::string & ca_str = _coloring_algorithm;
547 : Real * am_data = _adjacency_matrix->get_values().data();
548 :
549 : try
550 : {
551 103 : Moose::PetscSupport::colorAdjacencyMatrix(
552 103 : am_data, _feature_count, _vars.size(), _grain_idx_to_op, ca_str.c_str());
553 : }
554 2 : catch (std::runtime_error & e)
555 : {
556 2 : paramError("op_num",
557 : "Unable to find a valid grain to op coloring, Make sure you have created enough "
558 : "variables to hold a\nvalid polycrystal initial condition (no grains represented "
559 : "by the same variable should be allowed to\ntouch, ~8 for 2D, ~25 for 3D)?");
560 0 : }
561 : }
562 :
563 : /**
564 : * Now we have a vector giving us a coloring based on the indices in our features, but we need to
565 : * build a map in case our features have non-contiguous IDs.
566 : */
567 : mooseAssert(_grain_to_op.empty(), "grain_to_op data structure should be empty here");
568 3497 : for (MooseIndex(_grain_idx_to_op) i = 0; i < _grain_idx_to_op.size(); ++i)
569 3228 : _grain_to_op.emplace_hint(_grain_to_op.end(), _feature_sets[i]._id, _grain_idx_to_op[i]);
570 269 : }
571 :
572 : bool
573 58225596 : PolycrystalUserObjectBase::colorGraph(unsigned int vertex)
574 : {
575 : // Base case: All grains are assigned
576 58225596 : if (vertex == _feature_count)
577 : return true;
578 :
579 : // Consider this grain and try different ops
580 582064384 : for (unsigned int color_idx = 0; color_idx < _op_num; ++color_idx)
581 : {
582 : // We'll try to spread these colors around a bit rather than
583 : // packing them all on the first few colors if we have several colors.
584 523841156 : unsigned int color = (vertex + color_idx) % _op_num;
585 :
586 523841156 : if (isGraphValid(vertex, color))
587 : {
588 58225426 : _grain_idx_to_op[vertex] = color;
589 :
590 58225426 : if (colorGraph(vertex + 1))
591 : return true;
592 :
593 : // Backtrack...
594 58223226 : _grain_idx_to_op[vertex] = PolycrystalUserObjectBase::INVALID_COLOR;
595 : }
596 : }
597 :
598 : return false;
599 : }
600 :
601 : bool
602 523841156 : PolycrystalUserObjectBase::isGraphValid(unsigned int vertex, unsigned int color)
603 : {
604 : // See if the proposed color is valid based on the current neighbor colors
605 8050911539 : for (unsigned int neighbor = 0; neighbor < _feature_count; ++neighbor)
606 7992686113 : if ((*_adjacency_matrix)(vertex, neighbor) && color == _grain_idx_to_op[neighbor])
607 : return false;
608 : return true;
609 : }
610 :
611 : void
612 103 : PolycrystalUserObjectBase::printGrainAdjacencyMatrix() const
613 : {
614 103 : _console << "Grain Adjacency Matrix:\n";
615 1856 : for (unsigned int i = 0; i < _adjacency_matrix->m(); i++)
616 : {
617 47808 : for (unsigned int j = 0; j < _adjacency_matrix->n(); j++)
618 46055 : _console << _adjacency_matrix->el(i, j) << " ";
619 1753 : _console << '\n';
620 : }
621 :
622 103 : _console << "Grain to OP assignments:\n";
623 1856 : for (auto op : _grain_idx_to_op)
624 1753 : _console << op << " ";
625 103 : _console << '\n' << std::endl;
626 103 : }
627 :
628 : MooseEnum
629 1130 : PolycrystalUserObjectBase::coloringAlgorithms()
630 : {
631 2260 : return MooseEnum("jp power greedy bt", "jp");
632 : }
633 :
634 : std::string
635 1130 : PolycrystalUserObjectBase::coloringAlgorithmDescriptions()
636 : {
637 : return "The grain neighbor graph coloring algorithm to use: \"jp\" (DEFAULT) Jones and "
638 : "Plassmann, an efficient coloring algorithm, \"power\" an alternative stochastic "
639 : "algorithm, \"greedy\", a greedy assignment algorithm with stochastic updates to "
640 : "guarantee a valid coloring, \"bt\", a back tracking algorithm that produces good "
641 : "distributions but may experience exponential run time in the worst case scenario "
642 1130 : "(works well on medium to large 2D problems)";
643 : }
644 :
645 : const unsigned int PolycrystalUserObjectBase::INVALID_COLOR =
646 : std::numeric_limits<unsigned int>::max();
647 : const unsigned int PolycrystalUserObjectBase::HALO_THICKNESS = 4;
|
