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 29294 : PropertyReadFile::validParams()
24 : {
25 29294 : InputParameters params = GeneralUserObject::validParams();
26 29294 : 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 29294 : 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 29294 : 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 29294 : params.addRequiredParam<unsigned int>("nprop", "Number of tabulated property values");
40 87882 : params.addParam<unsigned int>(
41 58588 : "nvoronoi", 0, "Number of voronoi tesselations/grains/nearest neighbor regions");
42 29294 : params.addDeprecatedParam<unsigned int>(
43 : "ngrain", "Number of grains.", "ngrain is deprecated, use nvoronoi instead");
44 29294 : params.addParam<unsigned int>("nblock", 0, "Number of blocks");
45 87882 : params.addRequiredParam<MooseEnum>("read_type",
46 58588 : 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 87882 : params.addParam<bool>("use_random_voronoi",
54 58588 : false,
55 : "Wether to generate random positions for the Voronoi tesselation");
56 29294 : params.addParam<unsigned int>("rand_seed", 2000, "random seed");
57 87882 : params.addParam<MooseEnum>(
58 : "rve_type",
59 58588 : MooseEnum("periodic none", "none"),
60 : "Periodic or non-periodic grain distribution: Default is non-periodic");
61 87882 : params.addParam<bool>(
62 58588 : "use_zero_based_block_indexing", true, "Are the blocks numbered starting at zero?");
63 :
64 : // Data reading schedule
65 87882 : params.addParam<bool>(
66 : "load_first_file_on_construction",
67 58588 : 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 29294 : ExecFlagEnum & exec_enum = params.set<ExecFlagEnum>("execute_on", true);
73 29294 : exec_enum.removeAvailableFlags(EXEC_INITIAL, EXEC_FINAL, EXEC_LINEAR, EXEC_NONLINEAR);
74 29294 : params.setDocString("execute_on", exec_enum.getDocString());
75 : // we must read the files as early as possible
76 29294 : params.set<bool>("force_preaux") = true;
77 :
78 29294 : return params;
79 0 : }
80 :
81 378 : PropertyReadFile::PropertyReadFile(const InputParameters & parameters)
82 : : GeneralUserObject(parameters),
83 378 : _prop_file_names(getFileNames()),
84 370 : _current_file_index(declareRestartableData<unsigned int>("file_index", 0)),
85 : // index of files must be capped if restarting after having read all files
86 370 : _reader(
87 370 : _prop_file_names[std::min(_current_file_index, (unsigned int)_prop_file_names.size() - 1)]),
88 370 : _read_type(getParam<MooseEnum>("read_type").getEnum<ReadTypeEnum>()),
89 370 : _use_random_tesselation(getParam<bool>("use_random_voronoi")),
90 370 : _rand_seed(getParam<unsigned int>("rand_seed")),
91 370 : _rve_type(getParam<MooseEnum>("rve_type")),
92 370 : _block_zero(getParam<bool>("use_zero_based_block_indexing")),
93 740 : _ngrain(isParamValid("ngrain") ? getParam<unsigned int>("ngrain")
94 370 : : getParam<unsigned int>("nvoronoi")),
95 370 : _mesh(_fe_problem.mesh()),
96 370 : _nprop(getParam<unsigned int>("nprop")),
97 740 : _nvoronoi(isParamValid("ngrain") ? getParam<unsigned int>("ngrain")
98 370 : : getParam<unsigned int>("nvoronoi")),
99 370 : _nblock(getParam<unsigned int>("nblock")),
100 370 : _initialize_called_once(declareRestartableData<bool>("initialize_called", false)),
101 1118 : _load_on_construction(getParam<bool>("load_first_file_on_construction"))
102 : {
103 370 : 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 370 : Point mesh_min, mesh_max;
108 1480 : for (unsigned int i : make_range(Moose::dim))
109 : {
110 1110 : mesh_min(i) = _mesh.getMinInDimension(i);
111 1110 : mesh_max(i) = _mesh.getMaxInDimension(i);
112 : }
113 370 : _bounding_box = MooseUtils::buildBoundingBox(mesh_min, mesh_max);
114 :
115 370 : if (_load_on_construction)
116 370 : readData();
117 354 : }
118 :
119 : std::vector<FileName>
120 378 : PropertyReadFile::getFileNames()
121 : {
122 378 : std::vector<FileName> prop_file_names;
123 : // Retrieve the property file names
124 378 : if (isParamValid("prop_file_name") && !isParamValid("prop_file_names_csv"))
125 357 : prop_file_names = getParam<std::vector<FileName>>("prop_file_name");
126 21 : else if (isParamValid("prop_file_names_csv") && !isParamValid("prop_file_name"))
127 : {
128 : MooseUtils::DelimitedFileOfStringReader file_name_reader(
129 17 : getParam<FileName>("prop_file_names_csv"));
130 17 : file_name_reader.setHeaderFlag(MooseUtils::DelimitedFileOfStringReader::HeaderFlag::OFF);
131 17 : file_name_reader.read();
132 17 : if (file_name_reader.getData().size())
133 : {
134 13 : const auto fdata = file_name_reader.getData(0);
135 39 : for (const auto & fname : fdata)
136 26 : prop_file_names.push_back(fname);
137 13 : }
138 17 : }
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 374 : if (prop_file_names.empty())
144 4 : paramError("prop_file_name", "A property file should have been specified.");
145 370 : return prop_file_names;
146 0 : }
147 :
148 : void
149 312 : PropertyReadFile::initialize()
150 : {
151 : // Since we read at construction, no need to re-read the file on first initialize
152 312 : if (!_initialize_called_once)
153 : {
154 288 : _initialize_called_once = true;
155 288 : return;
156 : }
157 :
158 : // Set then read new file, only if we have not reached the last file
159 24 : if (_current_file_index < _prop_file_names.size())
160 : {
161 24 : _reader.setFileName(_prop_file_names[_current_file_index]);
162 24 : 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 394 : PropertyReadFile::readData()
170 : {
171 394 : 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 394 : 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 394 : _reader.setFormatFlag(MooseUtils::DelimitedFileReader::FormatFlag::ROWS);
177 394 : _reader.read();
178 :
179 394 : unsigned int nobjects = 0;
180 :
181 394 : switch (_read_type)
182 : {
183 140 : case ReadTypeEnum::ELEMENT:
184 140 : nobjects = _mesh.nElem();
185 140 : break;
186 :
187 86 : case ReadTypeEnum::NODE:
188 86 : nobjects = _mesh.nNodes();
189 86 : break;
190 :
191 86 : case ReadTypeEnum::VORONOI:
192 86 : if (_nvoronoi <= 0)
193 4 : paramError("nvoronoi", "Provide non-zero number of voronoi tesselations/grains.");
194 82 : nobjects = _nvoronoi;
195 82 : 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 82 : case ReadTypeEnum::BLOCK:
205 82 : if (_nblock <= 0)
206 4 : paramError("nblock", "Provide non-zero number of blocks.");
207 78 : nobjects = _nblock;
208 78 : break;
209 : }
210 :
211 : // make sure the data from file has enough rows and columns
212 386 : 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 382 : 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 5172 : for (unsigned int i = 0; i < nobjects; i++)
225 4794 : 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 378 : if (_read_type == ReadTypeEnum::VORONOI || _read_type == ReadTypeEnum::GRAIN)
234 82 : initVoronoiCenterPoints();
235 :
236 378 : _current_file_index++;
237 378 : }
238 :
239 : void
240 82 : PropertyReadFile::initVoronoiCenterPoints()
241 : {
242 82 : _center.resize(_nvoronoi);
243 :
244 : // Generate a random tesselation
245 82 : if (_use_random_tesselation)
246 : {
247 11 : MooseRandom::seed(_rand_seed);
248 44 : for (const auto i : make_range(_nvoronoi))
249 132 : for (const auto j : make_range(Moose::dim))
250 198 : _center[i](j) = _bounding_box.min()(j) +
251 198 : MooseRandom::rand() * (_bounding_box.max() - _bounding_box.min())(j);
252 : }
253 : // Read tesselation from file
254 : else
255 : {
256 284 : for (const auto i : make_range(_nvoronoi))
257 852 : for (const auto j : make_range(Moose::dim))
258 639 : _center[i](j) = _reader.getData(i)[j];
259 : }
260 82 : }
261 :
262 : Real
263 1888 : PropertyReadFile::getData(const Elem * const elem, const unsigned int prop_num) const
264 : {
265 1888 : if (prop_num > _nprop)
266 0 : paramError(
267 0 : "nprop", "Property number ", prop_num, " greater than total number of properties ", _nprop);
268 :
269 1888 : Real data = 0.0;
270 1888 : switch (_read_type)
271 : {
272 1376 : case ReadTypeEnum::ELEMENT:
273 1376 : data = getElementData(elem, prop_num);
274 1376 : 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 512 : case ReadTypeEnum::BLOCK:
286 512 : data = getBlockData(elem, prop_num);
287 512 : 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 1888 : return data;
295 : }
296 :
297 : Real
298 1376 : PropertyReadFile::getElementData(const Elem * elem, unsigned int prop_num) const
299 : {
300 1376 : unsigned int jelem = elem->id();
301 1376 : 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 1376 : return _reader.getData(jelem)[prop_num];
308 : }
309 :
310 : Real
311 512 : PropertyReadFile::getNodeData(const Node * const node, const unsigned int prop_num) const
312 : {
313 512 : unsigned int jnode = node->id();
314 512 : 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 512 : return _reader.getData(jnode)[prop_num];
321 : }
322 :
323 : Real
324 512 : PropertyReadFile::getBlockData(const Elem * elem, unsigned int prop_num) const
325 : {
326 512 : unsigned int offset = 0;
327 512 : if (!_block_zero)
328 0 : offset = 1;
329 :
330 512 : unsigned int elem_subdomain_id = elem->subdomain_id();
331 512 : 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 512 : return _reader.getData(elem_subdomain_id - offset)[prop_num];
339 : }
340 :
341 : Real
342 1536 : PropertyReadFile::getVoronoiData(const Point & point, const unsigned int prop_num) const
343 : {
344 1536 : Real min_dist = std::numeric_limits<Real>::max();
345 1536 : unsigned int ivoronoi = 0;
346 :
347 6144 : for (const auto i : make_range(_nvoronoi))
348 : {
349 4608 : Real dist = 0.0;
350 4608 : switch (_rve_type)
351 : {
352 1536 : case 0:
353 : // Calculates minimum periodic distance when "periodic" is specified
354 : // for rve_type
355 1536 : dist = minPeriodicDistance(_center[i], point);
356 1536 : break;
357 :
358 3072 : default:
359 : // Calculates minimum distance when nothing is specified
360 : // for rve_type
361 3072 : Point dist_vec = _center[i] - point;
362 3072 : dist = dist_vec.norm();
363 : }
364 :
365 4608 : if (dist < min_dist)
366 : {
367 2760 : min_dist = dist;
368 2760 : ivoronoi = i;
369 : }
370 : }
371 1536 : return _reader.getData(ivoronoi)[prop_num];
372 : }
373 :
374 : // TODO: this should probably use the built-in min periodic distance!
375 : Real
376 1536 : PropertyReadFile::minPeriodicDistance(const Point & c, const Point & p) const
377 : {
378 1536 : Point dist_vec = c - p;
379 1536 : Real min_dist = dist_vec.norm();
380 :
381 1536 : Real fac[3] = {-1.0, 0.0, 1.0};
382 6144 : for (unsigned int i = 0; i < 3; i++)
383 18432 : for (unsigned int j = 0; j < 3; j++)
384 55296 : for (unsigned int k = 0; k < 3; k++)
385 : {
386 41472 : Point p1;
387 41472 : p1(0) = p(0) + fac[i] * (_bounding_box.max() - _bounding_box.min())(0);
388 41472 : p1(1) = p(1) + fac[j] * (_bounding_box.max() - _bounding_box.min())(1);
389 41472 : p1(2) = p(2) + fac[k] * (_bounding_box.max() - _bounding_box.min())(2);
390 :
391 41472 : dist_vec = c - p1;
392 41472 : Real dist = dist_vec.norm();
393 :
394 41472 : if (dist < min_dist)
395 1424 : min_dist = dist;
396 : }
397 :
398 1536 : return min_dist;
399 : }
|
