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 "PropertyReadFile.h"
11 : #include "MooseRandom.h"
12 : #include "MooseMesh.h"
13 :
14 : #include <fstream>
15 :
16 : registerMooseObject("MooseApp", PropertyReadFile);
17 : registerMooseObjectRenamed("MooseApp",
18 : ElementPropertyReadFile,
19 : "06/30/2021 24:00",
20 : PropertyReadFile);
21 :
22 : InputParameters
23 29242 : PropertyReadFile::validParams()
24 : {
25 29242 : InputParameters params = GeneralUserObject::validParams();
26 29242 : params.addClassDescription("User Object to read property data from an external file and assign "
27 : "it to elements / nodes / subdomains etc.");
28 :
29 : // Data source file(s)
30 29242 : params.addParam<std::vector<FileName>>(
31 : "prop_file_name",
32 : "Name(s) of the property file name. If specifying multiple files, the next file will be read "
33 : "at initialization right before every execution");
34 29242 : params.addParam<FileName>("prop_file_names_csv",
35 : "Name(s) of a CSV file containing the list of the property file names "
36 : "as its first column, with no header.");
37 :
38 : // Data organization
39 29242 : params.addRequiredParam<unsigned int>("nprop", "Number of tabulated property values");
40 87726 : params.addParam<unsigned int>(
41 58484 : "nvoronoi", 0, "Number of voronoi tesselations/grains/nearest neighbor regions");
42 29242 : params.addDeprecatedParam<unsigned int>(
43 : "ngrain", "Number of grains.", "ngrain is deprecated, use nvoronoi instead");
44 29242 : params.addParam<unsigned int>("nblock", 0, "Number of blocks");
45 87726 : params.addRequiredParam<MooseEnum>("read_type",
46 58484 : MooseEnum("element voronoi block node grain"),
47 : "Type of property distribution: "
48 : "element:by element "
49 : "node: by node "
50 : "voronoi:nearest neighbor / voronoi grain structure "
51 : "block:by mesh block "
52 : "grain: deprecated, use voronoi");
53 87726 : params.addParam<bool>("use_random_voronoi",
54 58484 : false,
55 : "Wether to generate random positions for the Voronoi tesselation");
56 29242 : params.addParam<unsigned int>("rand_seed", 2000, "random seed");
57 87726 : params.addParam<MooseEnum>(
58 : "rve_type",
59 58484 : MooseEnum("periodic none", "none"),
60 : "Periodic or non-periodic grain distribution: Default is non-periodic");
61 87726 : params.addParam<bool>(
62 58484 : "use_zero_based_block_indexing", true, "Are the blocks numbered starting at zero?");
63 :
64 : // Data reading schedule
65 87726 : params.addParam<bool>(
66 : "load_first_file_on_construction",
67 58484 : true,
68 : "Whether to read the first CSV file on construction or on the first execution");
69 :
70 : // Set an execution schedule to what makes sense currently
71 : // We do not allow INITIAL because we read the file at construction
72 29242 : ExecFlagEnum & exec_enum = params.set<ExecFlagEnum>("execute_on", true);
73 29242 : exec_enum.removeAvailableFlags(EXEC_INITIAL, EXEC_FINAL, EXEC_LINEAR, EXEC_NONLINEAR);
74 29242 : params.setDocString("execute_on", exec_enum.getDocString());
75 : // we must read the files as early as possible
76 29242 : params.set<bool>("force_preaux") = true;
77 :
78 29242 : return params;
79 0 : }
80 :
81 352 : PropertyReadFile::PropertyReadFile(const InputParameters & parameters)
82 : : GeneralUserObject(parameters),
83 352 : _prop_file_names(getFileNames()),
84 344 : _current_file_index(declareRestartableData<unsigned int>("file_index", 0)),
85 : // index of files must be capped if restarting after having read all files
86 344 : _reader(
87 344 : _prop_file_names[std::min(_current_file_index, (unsigned int)_prop_file_names.size() - 1)]),
88 344 : _read_type(getParam<MooseEnum>("read_type").getEnum<ReadTypeEnum>()),
89 344 : _use_random_tesselation(getParam<bool>("use_random_voronoi")),
90 344 : _rand_seed(getParam<unsigned int>("rand_seed")),
91 344 : _rve_type(getParam<MooseEnum>("rve_type")),
92 344 : _block_zero(getParam<bool>("use_zero_based_block_indexing")),
93 688 : _ngrain(isParamValid("ngrain") ? getParam<unsigned int>("ngrain")
94 344 : : getParam<unsigned int>("nvoronoi")),
95 344 : _mesh(_fe_problem.mesh()),
96 344 : _nprop(getParam<unsigned int>("nprop")),
97 688 : _nvoronoi(isParamValid("ngrain") ? getParam<unsigned int>("ngrain")
98 344 : : getParam<unsigned int>("nvoronoi")),
99 344 : _nblock(getParam<unsigned int>("nblock")),
100 344 : _initialize_called_once(declareRestartableData<bool>("initialize_called", false)),
101 1040 : _load_on_construction(getParam<bool>("load_first_file_on_construction"))
102 : {
103 344 : if (!_use_random_tesselation && parameters.isParamSetByUser("rand_seed"))
104 0 : paramError("rand_seed",
105 : "Random seeds should only be provided if random tesselation is desired");
106 :
107 344 : Point mesh_min, mesh_max;
108 1376 : for (unsigned int i : make_range(Moose::dim))
109 : {
110 1032 : mesh_min(i) = _mesh.getMinInDimension(i);
111 1032 : mesh_max(i) = _mesh.getMaxInDimension(i);
112 : }
113 344 : _bounding_box = MooseUtils::buildBoundingBox(mesh_min, mesh_max);
114 :
115 344 : if (_load_on_construction)
116 344 : readData();
117 328 : }
118 :
119 : std::vector<FileName>
120 352 : PropertyReadFile::getFileNames()
121 : {
122 352 : std::vector<FileName> prop_file_names;
123 : // Retrieve the property file names
124 352 : if (isParamValid("prop_file_name") && !isParamValid("prop_file_names_csv"))
125 332 : prop_file_names = getParam<std::vector<FileName>>("prop_file_name");
126 20 : else if (isParamValid("prop_file_names_csv") && !isParamValid("prop_file_name"))
127 : {
128 : MooseUtils::DelimitedFileOfStringReader file_name_reader(
129 16 : getParam<FileName>("prop_file_names_csv"));
130 16 : file_name_reader.setHeaderFlag(MooseUtils::DelimitedFileOfStringReader::HeaderFlag::OFF);
131 16 : file_name_reader.read();
132 16 : if (file_name_reader.getData().size())
133 : {
134 12 : const auto fdata = file_name_reader.getData(0);
135 36 : for (const auto & fname : fdata)
136 24 : prop_file_names.push_back(fname);
137 12 : }
138 16 : }
139 : else
140 4 : paramError("prop_file_name",
141 : "The names of the property files must be provided with either 'prop_file_name' "
142 : "or 'prop_file_names_csv'. Providing both or none is not supported.");
143 348 : if (prop_file_names.empty())
144 4 : paramError("prop_file_name", "A property file should have been specified.");
145 344 : return prop_file_names;
146 0 : }
147 :
148 : void
149 284 : PropertyReadFile::initialize()
150 : {
151 : // Since we read at construction, no need to re-read the file on first initialize
152 284 : if (!_initialize_called_once)
153 : {
154 262 : _initialize_called_once = true;
155 262 : return;
156 : }
157 :
158 : // Set then read new file, only if we have not reached the last file
159 22 : if (_current_file_index < _prop_file_names.size())
160 : {
161 22 : _reader.setFileName(_prop_file_names[_current_file_index]);
162 22 : readData();
163 : }
164 0 : else if (_prop_file_names.size() > 1 && _current_file_index == _prop_file_names.size())
165 0 : mooseInfo("Last file specified has been read. The file will no longer be updated.");
166 : }
167 :
168 : void
169 366 : PropertyReadFile::readData()
170 : {
171 366 : if (_read_type == ReadTypeEnum::ELEMENT && _mesh.getMesh().allow_renumbering())
172 0 : mooseWarning("CSV data is sorted by element, but mesh element renumbering is on, be careful!");
173 366 : if (_read_type == ReadTypeEnum::NODE && _mesh.getMesh().allow_renumbering())
174 0 : mooseWarning("CSV data is sorted by node, but mesh node renumbering is on, be careful!");
175 :
176 366 : _reader.setFormatFlag(MooseUtils::DelimitedFileReader::FormatFlag::ROWS);
177 366 : _reader.read();
178 :
179 366 : unsigned int nobjects = 0;
180 :
181 366 : switch (_read_type)
182 : {
183 130 : case ReadTypeEnum::ELEMENT:
184 130 : nobjects = _mesh.nElem();
185 130 : break;
186 :
187 80 : case ReadTypeEnum::NODE:
188 80 : nobjects = _mesh.nNodes();
189 80 : break;
190 :
191 80 : case ReadTypeEnum::VORONOI:
192 80 : if (_nvoronoi <= 0)
193 4 : paramError("nvoronoi", "Provide non-zero number of voronoi tesselations/grains.");
194 76 : nobjects = _nvoronoi;
195 76 : break;
196 :
197 : // TODO Delete after Grizzly update see idaholab/Grizzly#182
198 0 : case ReadTypeEnum::GRAIN:
199 0 : if (_nvoronoi <= 0)
200 0 : paramError("ngrain", "Provide non-zero number of voronoi tesselations/grains.");
201 0 : nobjects = _nvoronoi;
202 0 : break;
203 :
204 76 : case ReadTypeEnum::BLOCK:
205 76 : if (_nblock <= 0)
206 4 : paramError("nblock", "Provide non-zero number of blocks.");
207 72 : nobjects = _nblock;
208 72 : break;
209 : }
210 :
211 : // make sure the data from file has enough rows and columns
212 358 : if (_reader.getData().size() < nobjects)
213 4 : mooseError("Data in ",
214 4 : _prop_file_names[_current_file_index],
215 : " does not have enough rows for ",
216 : nobjects,
217 : " objects.");
218 354 : if (_reader.getData().size() > nobjects)
219 4 : mooseWarning("Data size in ",
220 4 : _prop_file_names[_current_file_index],
221 : " is larger than ",
222 : nobjects,
223 : " objects, some data will not be used.");
224 4798 : for (unsigned int i = 0; i < nobjects; i++)
225 4448 : if (_reader.getData(i).size() < _nprop)
226 0 : mooseError("Row ",
227 : i,
228 : " in ",
229 0 : _prop_file_names[_current_file_index],
230 : " has number of data less than ",
231 0 : _nprop);
232 :
233 350 : if (_read_type == ReadTypeEnum::VORONOI || _read_type == ReadTypeEnum::GRAIN)
234 76 : initVoronoiCenterPoints();
235 :
236 350 : _current_file_index++;
237 350 : }
238 :
239 : void
240 76 : PropertyReadFile::initVoronoiCenterPoints()
241 : {
242 76 : _center.resize(_nvoronoi);
243 :
244 : // Generate a random tesselation
245 76 : if (_use_random_tesselation)
246 : {
247 10 : MooseRandom::seed(_rand_seed);
248 40 : for (const auto i : make_range(_nvoronoi))
249 120 : for (const auto j : make_range(Moose::dim))
250 180 : _center[i](j) = _bounding_box.min()(j) +
251 180 : MooseRandom::rand() * (_bounding_box.max() - _bounding_box.min())(j);
252 : }
253 : // Read tesselation from file
254 : else
255 : {
256 264 : for (const auto i : make_range(_nvoronoi))
257 792 : for (const auto j : make_range(Moose::dim))
258 594 : _center[i](j) = _reader.getData(i)[j];
259 : }
260 76 : }
261 :
262 : Real
263 1664 : PropertyReadFile::getData(const Elem * const elem, const unsigned int prop_num) const
264 : {
265 1664 : if (prop_num > _nprop)
266 0 : paramError(
267 0 : "nprop", "Property number ", prop_num, " greater than total number of properties ", _nprop);
268 :
269 1664 : Real data = 0.0;
270 1664 : switch (_read_type)
271 : {
272 1216 : case ReadTypeEnum::ELEMENT:
273 1216 : data = getElementData(elem, prop_num);
274 1216 : break;
275 :
276 0 : case ReadTypeEnum::VORONOI:
277 0 : data = getVoronoiData(elem->vertex_average(), prop_num);
278 0 : break;
279 :
280 : // TODO: Delete after Grizzly update
281 0 : case ReadTypeEnum::GRAIN:
282 0 : data = getVoronoiData(elem->vertex_average(), prop_num);
283 0 : break;
284 :
285 448 : case ReadTypeEnum::BLOCK:
286 448 : data = getBlockData(elem, prop_num);
287 448 : break;
288 :
289 0 : case ReadTypeEnum::NODE:
290 0 : mooseError(
291 : "PropertyReadFile has data sorted by node, it should note be retrieved by element");
292 : break;
293 : }
294 1664 : return data;
295 : }
296 :
297 : Real
298 1216 : PropertyReadFile::getElementData(const Elem * elem, unsigned int prop_num) const
299 : {
300 1216 : unsigned int jelem = elem->id();
301 1216 : if (jelem >= _mesh.nElem())
302 0 : mooseError("Element ID ",
303 : jelem,
304 : " greater than than total number of element in mesh: ",
305 0 : _mesh.nElem(),
306 : ". Elements should be numbered consecutively, and ids should start from 0.");
307 1216 : return _reader.getData(jelem)[prop_num];
308 : }
309 :
310 : Real
311 448 : PropertyReadFile::getNodeData(const Node * const node, const unsigned int prop_num) const
312 : {
313 448 : unsigned int jnode = node->id();
314 448 : if (jnode >= _mesh.nNodes())
315 0 : mooseError("Node ID ",
316 : jnode,
317 : " greater than than total number of nodes in mesh: ",
318 0 : _mesh.nNodes(),
319 : ". Nodes should be numbered consecutively, with ids starting from 0.");
320 448 : return _reader.getData(jnode)[prop_num];
321 : }
322 :
323 : Real
324 448 : PropertyReadFile::getBlockData(const Elem * elem, unsigned int prop_num) const
325 : {
326 448 : unsigned int offset = 0;
327 448 : if (!_block_zero)
328 0 : offset = 1;
329 :
330 448 : unsigned int elem_subdomain_id = elem->subdomain_id();
331 448 : if (elem_subdomain_id >= _nblock + offset)
332 0 : paramError("nblock",
333 : "Element block id ",
334 : elem_subdomain_id,
335 : " greater than than total number of blocks in mesh: ",
336 0 : _nblock,
337 : ". Blocks should be numbered consecutively, starting from 0.");
338 448 : return _reader.getData(elem_subdomain_id - offset)[prop_num];
339 : }
340 :
341 : Real
342 1344 : PropertyReadFile::getVoronoiData(const Point & point, const unsigned int prop_num) const
343 : {
344 1344 : Real min_dist = std::numeric_limits<Real>::max();
345 1344 : unsigned int ivoronoi = 0;
346 :
347 5376 : for (const auto i : make_range(_nvoronoi))
348 : {
349 4032 : Real dist = 0.0;
350 4032 : switch (_rve_type)
351 : {
352 1344 : case 0:
353 : // Calculates minimum periodic distance when "periodic" is specified
354 : // for rve_type
355 1344 : dist = minPeriodicDistance(_center[i], point);
356 1344 : break;
357 :
358 2688 : default:
359 : // Calculates minimum distance when nothing is specified
360 : // for rve_type
361 2688 : Point dist_vec = _center[i] - point;
362 2688 : dist = dist_vec.norm();
363 : }
364 :
365 4032 : if (dist < min_dist)
366 : {
367 2415 : min_dist = dist;
368 2415 : ivoronoi = i;
369 : }
370 : }
371 1344 : return _reader.getData(ivoronoi)[prop_num];
372 : }
373 :
374 : // TODO: this should probably use the built-in min periodic distance!
375 : Real
376 1344 : PropertyReadFile::minPeriodicDistance(const Point & c, const Point & p) const
377 : {
378 1344 : Point dist_vec = c - p;
379 1344 : Real min_dist = dist_vec.norm();
380 :
381 1344 : Real fac[3] = {-1.0, 0.0, 1.0};
382 5376 : for (unsigned int i = 0; i < 3; i++)
383 16128 : for (unsigned int j = 0; j < 3; j++)
384 48384 : for (unsigned int k = 0; k < 3; k++)
385 : {
386 36288 : Point p1;
387 36288 : p1(0) = p(0) + fac[i] * (_bounding_box.max() - _bounding_box.min())(0);
388 36288 : p1(1) = p(1) + fac[j] * (_bounding_box.max() - _bounding_box.min())(1);
389 36288 : p1(2) = p(2) + fac[k] * (_bounding_box.max() - _bounding_box.min())(2);
390 :
391 36288 : dist_vec = c - p1;
392 36288 : Real dist = dist_vec.norm();
393 :
394 36288 : if (dist < min_dist)
395 1246 : min_dist = dist;
396 : }
397 :
398 1344 : return min_dist;
399 : }
|
