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 "PolycrystalVoronoi.h"
11 : #include "IndirectSort.h"
12 : #include "MooseRandom.h"
13 : #include "MooseMesh.h"
14 : #include "MooseVariable.h"
15 : #include "NonlinearSystemBase.h"
16 : #include "DelimitedFileReader.h"
17 :
18 : registerMooseObject("PhaseFieldApp", PolycrystalVoronoi);
19 :
20 : InputParameters
21 1238 : PolycrystalVoronoi::validParams()
22 : {
23 1238 : InputParameters params = PolycrystalUserObjectBase::validParams();
24 1238 : params.addClassDescription(
25 : "Random Voronoi tessellation polycrystal (used by PolycrystalVoronoiAction)");
26 2476 : params.addParam<unsigned int>(
27 2476 : "grain_num", 0, "Number of grains being represented by the order parameters");
28 2476 : params.addParam<unsigned int>("rand_seed", 0, "The random seed");
29 2476 : params.addParam<bool>(
30 2476 : "columnar_3D", false, "3D microstructure will be columnar in the z-direction?");
31 2476 : params.addParam<bool>(
32 2476 : "use_kdtree", false, "Whether or not to use a KD tree to speedup grain search");
33 3714 : params.addRangeCheckedParam<unsigned int>("point_patch_size",
34 2476 : 1,
35 : "point_patch_size > 0",
36 : "How many nearest points KDTree should return");
37 3714 : params.addRangeCheckedParam<unsigned int>("grain_patch_size",
38 2476 : 10,
39 : "grain_patch_size > 0",
40 : "How many nearest grains KDTree should return");
41 2476 : params.addParam<FileName>(
42 : "file_name",
43 : "",
44 : "File containing grain centroids, if file_name is provided, the centroids "
45 : "from the file will be used.");
46 2476 : params.addParam<Real>("int_width", 0.0, "Width of diffuse interfaces");
47 1238 : return params;
48 0 : }
49 :
50 619 : PolycrystalVoronoi::PolycrystalVoronoi(const InputParameters & parameters)
51 : : PolycrystalUserObjectBase(parameters),
52 619 : _grain_num(getParam<unsigned int>("grain_num")),
53 1238 : _columnar_3D(getParam<bool>("columnar_3D")),
54 1238 : _rand_seed(getParam<unsigned int>("rand_seed")),
55 1238 : _int_width(getParam<Real>("int_width")),
56 1238 : _file_name(getParam<FileName>("file_name")),
57 : _kd_tree(nullptr),
58 1238 : _use_kdtree(getParam<bool>("use_kdtree")),
59 1238 : _point_patch_size(getParam<unsigned int>("point_patch_size")),
60 1857 : _grain_patch_size(getParam<unsigned int>("grain_patch_size"))
61 : {
62 619 : if (_file_name == "" && _grain_num == 0)
63 0 : mooseError("grain_num must be provided if the grain centroids are not read from a file");
64 :
65 619 : if (_file_name != "" && _grain_num > 0)
66 0 : mooseWarning("grain_num is ignored and will be determined from the file.");
67 619 : }
68 :
69 : void
70 10552196 : PolycrystalVoronoi::getGrainsBasedOnPoint(const Point & point,
71 : std::vector<unsigned int> & grains) const
72 : {
73 10552196 : grains.resize(0);
74 10552196 : Real d_min = _range.norm();
75 : Real distance;
76 : auto n_grains = _centerpoints.size();
77 : auto min_index = n_grains;
78 :
79 10552196 : if (_use_kdtree)
80 : {
81 : mooseAssert(_kd_tree, "KD tree is not constructed yet");
82 : mooseAssert(_grain_gtl_ids.size() == _new_points.size(),
83 : "The number of grain global IDs does not match that of new center points");
84 :
85 770830 : std::vector<std::size_t> return_index(_point_patch_size);
86 770830 : std::vector<Real> return_dist_sqr(_point_patch_size);
87 :
88 770830 : _kd_tree->neighborSearch(point, _point_patch_size, return_index, return_dist_sqr);
89 :
90 770830 : min_index = _grain_gtl_ids[return_index[0]];
91 :
92 770830 : d_min = return_dist_sqr[0];
93 :
94 : // By default, _point_patch_size is one. A larger _point_patch_size
95 : // may be useful if nearest nodes are not unique, and
96 : // we want to select the node that has the smallest ID
97 1261546 : for (unsigned int i = 1; i < _point_patch_size; i++)
98 : {
99 490716 : if (d_min > return_dist_sqr[i])
100 : {
101 0 : min_index = _grain_gtl_ids[return_index[i]];
102 : d_min = return_dist_sqr[i];
103 : }
104 490716 : else if (d_min == return_dist_sqr[i])
105 : {
106 9600 : min_index = min_index > _grain_gtl_ids[return_index[i]] ? _grain_gtl_ids[return_index[i]]
107 : : min_index;
108 : }
109 : }
110 770830 : }
111 : else
112 : {
113 : // Find the closest centerpoint to the current point
114 117229003 : for (MooseIndex(_centerpoints) grain = 0; grain < n_grains; ++grain)
115 : {
116 107447637 : distance = _mesh.minPeriodicDistance(_vars[0]->number(), _centerpoints[grain], point);
117 107447637 : if (distance < d_min)
118 : {
119 : d_min = distance;
120 : min_index = grain;
121 : }
122 : }
123 : }
124 :
125 : mooseAssert(min_index < n_grains, "Couldn't find closest Voronoi cell");
126 10552196 : Point current_grain = _centerpoints[min_index];
127 10552196 : grains.push_back(min_index); // closest centerpoint always gets included
128 :
129 10552196 : if (_int_width > 0.0)
130 : {
131 2188650 : if (_use_kdtree)
132 : {
133 : mooseAssert(_grain_patch_size < _grain_gtl_ids.size(),
134 : "Number of neighboring grains should not exceed number of global grains");
135 :
136 607258 : std::vector<std::size_t> return_index(_grain_patch_size);
137 607258 : _kd_tree->neighborSearch(current_grain, _grain_patch_size, return_index);
138 :
139 : std::set<dof_id_type> neighbor_grains;
140 6679838 : for (unsigned int i = 0; i < _grain_patch_size; i++)
141 6072580 : if (_grain_gtl_ids[return_index[i]] != min_index)
142 5465322 : neighbor_grains.insert(_grain_gtl_ids[return_index[i]]);
143 :
144 3036290 : for (auto it = neighbor_grains.begin(); it != neighbor_grains.end(); ++it)
145 2429032 : if ((*it) != min_index)
146 : {
147 2429032 : Point next_grain = _centerpoints[*it];
148 2429032 : Point N = findNormalVector(point, current_grain, next_grain);
149 2429032 : Point cntr = findCenterPoint(point, current_grain, next_grain);
150 : distance = N * (cntr - point);
151 2429032 : if (distance < _int_width)
152 250890 : grains.push_back(*it); // also include all grains with nearby boundaries
153 : }
154 607258 : }
155 : else
156 : {
157 9191882 : for (MooseIndex(_centerpoints) grain = 0; grain < n_grains; ++grain)
158 7610490 : if (grain != min_index)
159 : {
160 6029098 : Point next_grain = _centerpoints[grain];
161 6029098 : Point N = findNormalVector(point, current_grain, next_grain);
162 6029098 : Point cntr = findCenterPoint(point, current_grain, next_grain);
163 : distance = N * (cntr - point);
164 6029098 : if (distance < _int_width)
165 657105 : grains.push_back(grain); // also include all grains with nearby boundaries
166 : }
167 : }
168 : }
169 10552196 : }
170 :
171 : Real
172 10171249 : PolycrystalVoronoi::getVariableValue(unsigned int op_index, const Point & p) const
173 : {
174 : std::vector<unsigned int> grain_ids;
175 10171249 : getGrainsBasedOnPoint(p, grain_ids);
176 :
177 : // Now see if any of those grains are represented by the passed in order parameter
178 10171249 : unsigned int active_grain_on_op = invalid_id;
179 18891730 : for (auto grain_id : grain_ids)
180 10820721 : if (op_index == _grain_to_op.at(grain_id))
181 : {
182 2100240 : active_grain_on_op = grain_id;
183 2100240 : break;
184 : }
185 :
186 : Real profile_val = 0.0;
187 10171249 : if (active_grain_on_op != invalid_id)
188 2100240 : profile_val = computeDiffuseInterface(p, active_grain_on_op, grain_ids);
189 :
190 10171249 : return profile_val;
191 10171249 : }
192 :
193 : void
194 362 : PolycrystalVoronoi::precomputeGrainStructure()
195 : {
196 : // Set up domain bounds with mesh tools
197 1448 : for (const auto i : make_range(Moose::dim))
198 : {
199 1086 : _bottom_left(i) = _mesh.getMinInDimension(i);
200 1086 : _top_right(i) = _mesh.getMaxInDimension(i);
201 : }
202 362 : _range = _top_right - _bottom_left;
203 :
204 362 : if (!_file_name.empty())
205 : {
206 10 : MooseUtils::DelimitedFileReader txt_reader(_file_name, &_communicator);
207 :
208 10 : txt_reader.read();
209 10 : std::vector<std::string> col_names = txt_reader.getNames();
210 10 : std::vector<std::vector<Real>> data = txt_reader.getData();
211 10 : _grain_num = data[0].size();
212 10 : _centerpoints.resize(_grain_num);
213 :
214 30 : for (unsigned int i = 0; i < col_names.size(); ++i)
215 : {
216 : // Check vector lengths
217 20 : if (data[i].size() != _grain_num)
218 0 : paramError("Columns in ", _file_name, " do not have uniform lengths.");
219 : }
220 :
221 280 : for (unsigned int grain = 0; grain < _grain_num; ++grain)
222 : {
223 270 : _centerpoints[grain](0) = data[0][grain];
224 270 : _centerpoints[grain](1) = data[1][grain];
225 270 : if (col_names.size() > 2)
226 0 : _centerpoints[grain](2) = data[2][grain];
227 : else
228 270 : _centerpoints[grain](2) = 0.0;
229 : }
230 10 : }
231 : else
232 : {
233 352 : MooseRandom::seed(_rand_seed);
234 :
235 : // Randomly generate the centers of the individual grains represented by the Voronoi
236 : // tessellation
237 352 : _centerpoints.resize(_grain_num);
238 352 : std::vector<Real> distances(_grain_num);
239 :
240 3651 : for (auto grain = decltype(_grain_num)(0); grain < _grain_num; grain++)
241 : {
242 13196 : for (const auto i : make_range(Moose::dim))
243 9897 : _centerpoints[grain](i) = _bottom_left(i) + _range(i) * MooseRandom::rand();
244 3299 : if (_columnar_3D)
245 32 : _centerpoints[grain](2) = _bottom_left(2) + _range(2) * 0.5;
246 : }
247 352 : }
248 :
249 : // Build a KDTree that is used to speedup point search
250 362 : buildSearchTree();
251 362 : }
252 :
253 : void
254 391 : PolycrystalVoronoi::buildSearchTree()
255 : {
256 391 : if (!_use_kdtree)
257 378 : return;
258 :
259 : auto midplane = _bottom_left + _range / 2.0;
260 13 : dof_id_type local_grain_id = 0;
261 13 : _grain_gtl_ids.clear();
262 13 : _grain_gtl_ids.reserve(_centerpoints.size() * std::pow(2.0, _mesh.dimension()));
263 13 : _new_points.clear();
264 : // Use more memory. Factor is 2^dim
265 13 : _new_points.reserve(_centerpoints.size() * std::pow(2.0, _mesh.dimension()));
266 : // Domain will be extended along the periodic directions
267 : // For each direction, a half domain is constructed
268 13 : std::vector<std::vector<Real>> xyzs(LIBMESH_DIM);
269 70 : for (auto & point : _centerpoints)
270 : {
271 : // Cear up
272 228 : for (unsigned int i = 0; i < LIBMESH_DIM; i++)
273 171 : xyzs[i].clear();
274 :
275 : // Have at least the original point
276 228 : for (unsigned int i = 0; i < LIBMESH_DIM; i++)
277 171 : xyzs[i].push_back(point(i));
278 :
279 : // Add new coords when there exists periodic boundary conditions
280 : // We extend half domain
281 171 : for (unsigned int i = 0; i < _mesh.dimension(); i++)
282 114 : if (_mesh.isTranslatedPeriodic(_vars[0]->number(), i))
283 114 : xyzs[i].push_back(point(i) <= midplane(i) ? point(i) + _range(i) : point(i) - _range(i));
284 :
285 : // Construct all combinations
286 171 : for (auto x : xyzs[0])
287 342 : for (auto y : xyzs[1])
288 456 : for (auto z : xyzs[2])
289 : {
290 228 : _new_points.push_back(Point(x, y, z));
291 228 : _grain_gtl_ids.push_back(local_grain_id);
292 : }
293 :
294 57 : local_grain_id++;
295 : }
296 :
297 : // Build a KDTree that is used to speedup point search
298 : // We should not destroy _new_points after the tree is contributed
299 : // KDTree use a reference to "_new_points"
300 26 : _kd_tree = std::make_unique<KDTree>(_new_points, _mesh.getMaxLeafSize());
301 13 : }
302 :
303 : Real
304 2100240 : PolycrystalVoronoi::computeDiffuseInterface(const Point & point,
305 : const unsigned int & gr_index,
306 : const std::vector<unsigned int> & grain_ids) const
307 : {
308 : Real val = 1.0;
309 2100240 : Point current_grain = _centerpoints[gr_index];
310 4637356 : for (auto i : grain_ids)
311 2537116 : if (i != gr_index)
312 : {
313 436876 : Point next_grain = _centerpoints[i];
314 436876 : Point N = findNormalVector(point, current_grain, next_grain);
315 436876 : Point cntr = findCenterPoint(point, current_grain, next_grain);
316 436876 : for (unsigned int vcomp = 0; vcomp < 3; ++vcomp)
317 436876 : if (N(vcomp) != 0.0)
318 : {
319 436876 : Real L = findLinePoint(point, N, cntr, vcomp);
320 436876 : val *= 0.5 * (1.0 - std::tanh(2.0 * (point(vcomp) - L) * N(vcomp) / _int_width));
321 436876 : break;
322 : }
323 : }
324 2100240 : return val;
325 : }
326 :
327 : Point
328 8895006 : PolycrystalVoronoi::findNormalVector(const Point & point, const Point & p1, const Point & p2) const
329 : {
330 8895006 : Point pa = point + _mesh.minPeriodicVector(_vars[0]->number(), point, p1);
331 8895006 : Point pb = point + _mesh.minPeriodicVector(_vars[0]->number(), point, p2);
332 : Point N = pb - pa;
333 8895006 : return N / N.norm();
334 : }
335 :
336 : Point
337 8895006 : PolycrystalVoronoi::findCenterPoint(const Point & point, const Point & p1, const Point & p2) const
338 : {
339 8895006 : Point pa = point + _mesh.minPeriodicVector(_vars[0]->number(), point, p1);
340 8895006 : Point pb = point + _mesh.minPeriodicVector(_vars[0]->number(), point, p2);
341 8895006 : return 0.5 * (pa + pb);
342 : }
343 :
344 : Real
345 436876 : PolycrystalVoronoi::findLinePoint(const Point & point,
346 : const Point & N,
347 : const Point & cntr,
348 : const unsigned int vcomp) const
349 : {
350 436876 : const Real l_sum = N((vcomp + 1) % 3) * (point((vcomp + 1) % 3) - cntr((vcomp + 1) % 3)) +
351 436876 : N((vcomp + 2) % 3) * (point((vcomp + 2) % 3) - cntr((vcomp + 2) % 3));
352 :
353 436876 : return cntr(vcomp) - l_sum / N(vcomp);
354 : }
|
