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 880 : PolycrystalVoronoi::validParams()
22 : {
23 880 : InputParameters params = PolycrystalUserObjectBase::validParams();
24 880 : params.addClassDescription(
25 : "Random Voronoi tessellation polycrystal (used by PolycrystalVoronoiAction)");
26 1760 : params.addParam<unsigned int>(
27 1760 : "grain_num", 0, "Number of grains being represented by the order parameters");
28 1760 : params.addParam<unsigned int>("rand_seed", 0, "The random seed");
29 1760 : params.addParam<bool>(
30 1760 : "columnar_3D", false, "3D microstructure will be columnar in the z-direction?");
31 1760 : params.addParam<bool>(
32 1760 : "use_kdtree", false, "Whether or not to use a KD tree to speedup grain search");
33 2640 : params.addRangeCheckedParam<unsigned int>("point_patch_size",
34 1760 : 1,
35 : "point_patch_size > 0",
36 : "How many nearest points KDTree should return");
37 2640 : params.addRangeCheckedParam<unsigned int>("grain_patch_size",
38 1760 : 10,
39 : "grain_patch_size > 0",
40 : "How many nearest grains KDTree should return");
41 1760 : 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 1760 : params.addParam<Real>("int_width", 0.0, "Width of diffuse interfaces");
47 880 : return params;
48 0 : }
49 :
50 440 : PolycrystalVoronoi::PolycrystalVoronoi(const InputParameters & parameters)
51 : : PolycrystalUserObjectBase(parameters),
52 440 : _grain_num(getParam<unsigned int>("grain_num")),
53 880 : _columnar_3D(getParam<bool>("columnar_3D")),
54 880 : _rand_seed(getParam<unsigned int>("rand_seed")),
55 880 : _int_width(getParam<Real>("int_width")),
56 880 : _file_name(getParam<FileName>("file_name")),
57 : _kd_tree(nullptr),
58 880 : _use_kdtree(getParam<bool>("use_kdtree")),
59 880 : _point_patch_size(getParam<unsigned int>("point_patch_size")),
60 1320 : _grain_patch_size(getParam<unsigned int>("grain_patch_size"))
61 : {
62 440 : 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 440 : if (_file_name != "" && _grain_num > 0)
66 0 : mooseWarning("grain_num is ignored and will be determined from the file.");
67 440 : }
68 :
69 : void
70 9431780 : PolycrystalVoronoi::getGrainsBasedOnPoint(const Point & point,
71 : std::vector<unsigned int> & grains) const
72 : {
73 9431780 : grains.resize(0);
74 9431780 : Real d_min = _range.norm();
75 : Real distance;
76 : auto n_grains = _centerpoints.size();
77 : auto min_index = n_grains;
78 :
79 9431780 : 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 743444 : std::vector<std::size_t> return_index(_point_patch_size);
86 743444 : std::vector<Real> return_dist_sqr(_point_patch_size);
87 :
88 743444 : _kd_tree->neighborSearch(point, _point_patch_size, return_index, return_dist_sqr);
89 :
90 743444 : min_index = _grain_gtl_ids[return_index[0]];
91 :
92 743444 : 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 1152002 : for (unsigned int i = 1; i < _point_patch_size; i++)
98 : {
99 408558 : 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 408558 : else if (d_min == return_dist_sqr[i])
105 : {
106 8000 : min_index = min_index > _grain_gtl_ids[return_index[i]] ? _grain_gtl_ids[return_index[i]]
107 : : min_index;
108 : }
109 : }
110 743444 : }
111 : else
112 : {
113 : // Find the closest centerpoint to the current point
114 105737349 : for (MooseIndex(_centerpoints) grain = 0; grain < n_grains; ++grain)
115 : {
116 97049013 : distance = _mesh.minPeriodicDistance(*_vars[0], _centerpoints[grain], point);
117 97049013 : 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 9431780 : Point current_grain = _centerpoints[min_index];
127 9431780 : grains.push_back(min_index); // closest centerpoint always gets included
128 :
129 9431780 : if (_int_width > 0.0)
130 : {
131 2034794 : 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 8309350 : for (MooseIndex(_centerpoints) grain = 0; grain < n_grains; ++grain)
158 6881814 : if (grain != min_index)
159 : {
160 5454278 : Point next_grain = _centerpoints[grain];
161 5454278 : Point N = findNormalVector(point, current_grain, next_grain);
162 5454278 : Point cntr = findCenterPoint(point, current_grain, next_grain);
163 : distance = N * (cntr - point);
164 5454278 : if (distance < _int_width)
165 592200 : grains.push_back(grain); // also include all grains with nearby boundaries
166 : }
167 : }
168 : }
169 9431780 : }
170 :
171 : Real
172 9090673 : PolycrystalVoronoi::getVariableValue(unsigned int op_index, const Point & p) const
173 : {
174 : std::vector<unsigned int> grain_ids;
175 9090673 : getGrainsBasedOnPoint(p, grain_ids);
176 :
177 : // Now see if any of those grains are represented by the passed in order parameter
178 9090673 : unsigned int active_grain_on_op = invalid_id;
179 16903930 : for (auto grain_id : grain_ids)
180 9693953 : if (op_index == _grain_to_op.at(grain_id))
181 : {
182 1880696 : active_grain_on_op = grain_id;
183 1880696 : break;
184 : }
185 :
186 : Real profile_val = 0.0;
187 9090673 : if (active_grain_on_op != invalid_id)
188 1880696 : profile_val = computeDiffuseInterface(p, active_grain_on_op, grain_ids);
189 :
190 9090673 : return profile_val;
191 9090673 : }
192 :
193 : void
194 292 : PolycrystalVoronoi::precomputeGrainStructure()
195 : {
196 : // Set up domain bounds with mesh tools
197 1168 : for (const auto i : make_range(Moose::dim))
198 : {
199 876 : _bottom_left(i) = _mesh.getMinInDimension(i);
200 876 : _top_right(i) = _mesh.getMaxInDimension(i);
201 : }
202 292 : _range = _top_right - _bottom_left;
203 :
204 292 : 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 282 : MooseRandom::seed(_rand_seed);
234 :
235 : // Randomly generate the centers of the individual grains represented by the Voronoi
236 : // tessellation
237 282 : _centerpoints.resize(_grain_num);
238 282 : std::vector<Real> distances(_grain_num);
239 :
240 2944 : for (auto grain = decltype(_grain_num)(0); grain < _grain_num; grain++)
241 : {
242 10648 : for (const auto i : make_range(Moose::dim))
243 7986 : _centerpoints[grain](i) = _bottom_left(i) + _range(i) * MooseRandom::rand();
244 2662 : if (_columnar_3D)
245 24 : _centerpoints[grain](2) = _bottom_left(2) + _range(2) * 0.5;
246 : }
247 282 : }
248 :
249 : // Build a KDTree that is used to speedup point search
250 292 : buildSearchTree();
251 292 : }
252 :
253 : void
254 315 : PolycrystalVoronoi::buildSearchTree()
255 : {
256 315 : if (!_use_kdtree)
257 304 : return;
258 :
259 : auto midplane = _bottom_left + _range / 2.0;
260 11 : dof_id_type local_grain_id = 0;
261 11 : _grain_gtl_ids.clear();
262 11 : _grain_gtl_ids.reserve(_centerpoints.size() * std::pow(2.0, _mesh.dimension()));
263 11 : _new_points.clear();
264 : // Use more memory. Factor is 2^dim
265 11 : _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 11 : std::vector<std::vector<Real>> xyzs(LIBMESH_DIM);
269 11 : const auto & periodic_dims = _mesh.queryPeriodicDimensions(*_vars[0]);
270 60 : for (auto & point : _centerpoints)
271 : {
272 : // Cear up
273 196 : for (unsigned int i = 0; i < LIBMESH_DIM; i++)
274 147 : xyzs[i].clear();
275 :
276 : // Have at least the original point
277 196 : for (unsigned int i = 0; i < LIBMESH_DIM; i++)
278 147 : xyzs[i].push_back(point(i));
279 :
280 : // Add new coords when there exists periodic boundary conditions
281 : // We extend half domain
282 147 : for (unsigned int i = 0; i < _mesh.dimension(); i++)
283 98 : if (periodic_dims[i])
284 98 : xyzs[i].push_back(point(i) <= midplane(i) ? point(i) + _range(i) : point(i) - _range(i));
285 :
286 : // Construct all combinations
287 147 : for (auto x : xyzs[0])
288 294 : for (auto y : xyzs[1])
289 392 : for (auto z : xyzs[2])
290 : {
291 196 : _new_points.push_back(Point(x, y, z));
292 196 : _grain_gtl_ids.push_back(local_grain_id);
293 : }
294 :
295 49 : local_grain_id++;
296 : }
297 :
298 : // Build a KDTree that is used to speedup point search
299 : // We should not destroy _new_points after the tree is contributed
300 : // KDTree use a reference to "_new_points"
301 22 : _kd_tree = std::make_unique<KDTree>(_new_points, _mesh.getMaxLeafSize());
302 11 : }
303 :
304 : Real
305 1880696 : PolycrystalVoronoi::computeDiffuseInterface(const Point & point,
306 : const unsigned int & gr_index,
307 : const std::vector<unsigned int> & grain_ids) const
308 : {
309 : Real val = 1.0;
310 1880696 : Point current_grain = _centerpoints[gr_index];
311 4166266 : for (auto i : grain_ids)
312 2285570 : if (i != gr_index)
313 : {
314 404874 : Point next_grain = _centerpoints[i];
315 404874 : Point N = findNormalVector(point, current_grain, next_grain);
316 404874 : Point cntr = findCenterPoint(point, current_grain, next_grain);
317 404874 : for (unsigned int vcomp = 0; vcomp < 3; ++vcomp)
318 404874 : if (N(vcomp) != 0.0)
319 : {
320 404874 : Real L = findLinePoint(point, N, cntr, vcomp);
321 404874 : val *= 0.5 * (1.0 - std::tanh(2.0 * (point(vcomp) - L) * N(vcomp) / _int_width));
322 404874 : break;
323 : }
324 : }
325 1880696 : return val;
326 : }
327 :
328 : Point
329 8288184 : PolycrystalVoronoi::findNormalVector(const Point & point, const Point & p1, const Point & p2) const
330 : {
331 8288184 : Point pa = point + _mesh.minPeriodicVector(*_vars[0], point, p1);
332 8288184 : Point pb = point + _mesh.minPeriodicVector(*_vars[0], point, p2);
333 : Point N = pb - pa;
334 8288184 : return N / N.norm();
335 : }
336 :
337 : Point
338 8288184 : PolycrystalVoronoi::findCenterPoint(const Point & point, const Point & p1, const Point & p2) const
339 : {
340 8288184 : Point pa = point + _mesh.minPeriodicVector(*_vars[0], point, p1);
341 8288184 : Point pb = point + _mesh.minPeriodicVector(*_vars[0], point, p2);
342 8288184 : return 0.5 * (pa + pb);
343 : }
344 :
345 : Real
346 404874 : PolycrystalVoronoi::findLinePoint(const Point & point,
347 : const Point & N,
348 : const Point & cntr,
349 : const unsigned int vcomp) const
350 : {
351 404874 : const Real l_sum = N((vcomp + 1) % 3) * (point((vcomp + 1) % 3) - cntr((vcomp + 1) % 3)) +
352 404874 : N((vcomp + 2) % 3) * (point((vcomp + 2) % 3) - cntr((vcomp + 2) % 3));
353 :
354 404874 : return cntr(vcomp) - l_sum / N(vcomp);
355 : }
|
