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