www.mooseframework.org
EBSDReader.C
Go to the documentation of this file.
1 /****************************************************************/
2 /* MOOSE - Multiphysics Object Oriented Simulation Environment */
3 /* */
4 /* All contents are licensed under LGPL V2.1 */
5 /* See LICENSE for full restrictions */
6 /****************************************************************/
7 
8 #include "EBSDReader.h"
9 #include "EBSDMesh.h"
10 #include "MooseMesh.h"
11 #include "Conversion.h"
12 #include "NonlinearSystem.h"
13 
14 #include <fstream>
15 
16 template <>
17 InputParameters
19 {
20  InputParameters params = validParams<EulerAngleProvider>();
21  params.addClassDescription("Load and manage DREAM.3D EBSD data files for running simulations on "
22  "reconstructed microstructures.");
23  params.addParam<unsigned int>(
24  "custom_columns", 0, "Number of additional custom data columns to read from the EBSD file");
25  return params;
26 }
27 
28 EBSDReader::EBSDReader(const InputParameters & params)
29  : EulerAngleProvider(params),
30  _mesh(_fe_problem.mesh()),
31  _nl(_fe_problem.getNonlinearSystemBase()),
32  _grain_num(0),
33  _custom_columns(getParam<unsigned int>("custom_columns")),
34  _time_step(_fe_problem.timeStep()),
35  _mesh_dimension(_mesh.dimension()),
36  _nx(0),
37  _ny(0),
38  _nz(0),
39  _dx(0.),
40  _dy(0.),
41  _dz(0.)
42 {
43  readFile();
44 }
45 
46 void
48 {
49  // No need to re-read data upon recovery
50  if (_app.isRecovering())
51  return;
52 
53  // Fetch and check mesh
54  EBSDMesh * mesh = dynamic_cast<EBSDMesh *>(&_mesh);
55  if (mesh == NULL)
56  mooseError("Please use an EBSDMesh in your simulation.");
57 
58  std::ifstream stream_in(mesh->getEBSDFilename().c_str());
59  if (!stream_in)
60  mooseError("Can't open EBSD file: ", mesh->getEBSDFilename());
61 
62  const EBSDMesh::EBSDMeshGeometry & g = mesh->getEBSDGeometry();
63 
64  // Copy file header data from the EBSDMesh
65  _dx = g.d[0];
66  _nx = g.n[0];
67  _minx = g.min[0];
68  _maxx = _minx + _dx * _nx;
69 
70  _dy = g.d[1];
71  _ny = g.n[1];
72  _miny = g.min[1];
73  _maxy = _miny + _dy * _ny;
74 
75  _dz = g.d[2];
76  _nz = g.n[2];
77  _minz = g.min[2];
78  _maxz = _minz + _dz * _nz;
79 
80  // Resize the _data array
81  unsigned total_size = g.dim < 3 ? _nx * _ny : _nx * _ny * _nz;
82  _data.resize(total_size);
83 
84  std::string line;
85  while (std::getline(stream_in, line))
86  {
87  if (line.find("#") != 0)
88  {
89  // Temporary variables to read in on each line
90  EBSDPointData d;
91  Real x, y, z;
92 
93  std::istringstream iss(line);
94  iss >> d._phi1 >> d._Phi >> d._phi2 >> x >> y >> z >> d._feature_id >> d._phase >>
95  d._symmetry;
96 
97  // Transform angles to degrees
98  d._phi1 *= 180.0 / libMesh::pi;
99  d._Phi *= 180.0 / libMesh::pi;
100  d._phi2 *= 180.0 / libMesh::pi;
101 
102  // Custom columns
103  d._custom.resize(_custom_columns);
104  for (unsigned int i = 0; i < _custom_columns; ++i)
105  if (!(iss >> d._custom[i]))
106  mooseError("Unable to read in EBSD custom data column #", i);
107 
108  if (x < _minx || y < _miny || x > _maxx || y > _maxy ||
109  (g.dim == 3 && (z < _minz || z > _maxz)))
110  mooseError("EBSD Data ouside of the domain declared in the header ([",
111  _minx,
112  ':',
113  _maxx,
114  "], [",
115  _miny,
116  ':',
117  _maxy,
118  "], [",
119  _minz,
120  ':',
121  _maxz,
122  "]) dim=",
123  g.dim,
124  "\n",
125  line);
126 
127  d._p = Point(x, y, z);
128 
129  // determine number of grains in the dataset
130  if (_global_id_map.find(d._feature_id) == _global_id_map.end())
132 
133  unsigned int global_index = indexFromPoint(Point(x, y, z));
134  _data[global_index] = d;
135  }
136  }
137  stream_in.close();
138 
139  // Resize the variables
140  _avg_data.resize(_grain_num);
141  _avg_angles.resize(_grain_num);
142 
143  // clear the averages
144  for (unsigned int i = 0; i < _grain_num; ++i)
145  {
146  EBSDAvgData & a = _avg_data[i];
147  a._symmetry = a._phase = a._n = 0;
148  a._p = 0.0;
149  a._custom.assign(_custom_columns, 0.0);
150 
151  EulerAngles & b = _avg_angles[i];
152  b.phi1 = b.Phi = b.phi2 = 0.0;
153  }
154 
155  // Iterate through data points to get average variable values for each grain
156  for (auto & j : _data)
157  {
158  EBSDAvgData & a = _avg_data[_global_id_map[j._feature_id]];
159  EulerAngles & b = _avg_angles[_global_id_map[j._feature_id]];
160 
161  // use Eigen::Quaternion<Real> here?
162  b.phi1 += j._phi1;
163  b.Phi += j._Phi;
164  b.phi2 += j._phi2;
165 
166  if (a._n == 0)
167  a._phase = j._phase;
168  else if (a._phase != j._phase)
169  mooseError("An EBSD feature needs to have a uniform phase.");
170 
171  if (a._n == 0)
172  a._symmetry = j._symmetry;
173  else if (a._symmetry != j._symmetry)
174  mooseError("An EBSD feature needs to have a uniform symmetry parameter.");
175 
176  for (unsigned int i = 0; i < _custom_columns; ++i)
177  a._custom[i] += j._custom[i];
178 
179  // store the feature (or grain) ID
180  a._feature_id = j._feature_id;
181 
182  a._p += j._p;
183  a._n++;
184  }
185 
186  for (unsigned int i = 0; i < _grain_num; ++i)
187  {
188  EBSDAvgData & a = _avg_data[i];
189  EulerAngles & b = _avg_angles[i];
190 
191  if (a._n == 0)
192  continue;
193 
194  // TODO: need better way to average angles
195  b.phi1 /= Real(a._n);
196  b.Phi /= Real(a._n);
197  b.phi2 /= Real(a._n);
198 
199  // link the EulerAngles into the EBSDAvgData for access via the functors
200  a._angles = &b;
201 
202  if (a._phase >= _global_id.size())
203  _global_id.resize(a._phase + 1);
204 
205  // The averaged per grain data locally contains the phase id, local id, and
206  // original feature id. It is stored contiguously indexed by global id.
207  a._local_id = _global_id[a._phase].size();
208  _global_id[a._phase].push_back(i);
209 
210  a._p *= 1.0 / Real(a._n);
211 
212  for (unsigned int i = 0; i < _custom_columns; ++i)
213  a._custom[i] /= Real(a._n);
214  }
215 
216  // Build maps to indicate the weights with which grain and phase data
217  // from the surrounding elements contributes to a node fo IC purposes
219 }
220 
222 
224 EBSDReader::getData(const Point & p) const
225 {
226  return _data[indexFromPoint(p)];
227 }
228 
230 EBSDReader::getAvgData(unsigned int var) const
231 {
232  return _avg_data[indexFromIndex(var)];
233 }
234 
235 const EulerAngles &
236 EBSDReader::getEulerAngles(unsigned int var) const
237 {
238  return _avg_angles[indexFromIndex(var)];
239 }
240 
242 EBSDReader::getAvgData(unsigned int phase, unsigned int local_id) const
243 {
244  return _avg_data[indexFromIndex(_global_id[phase][local_id])];
245 }
246 
247 unsigned int
249 {
250  return _grain_num;
251 }
252 
253 unsigned int
254 EBSDReader::getGrainNum(unsigned int phase) const
255 {
256  return _global_id[phase].size();
257 }
258 
259 unsigned int
260 EBSDReader::indexFromPoint(const Point & p) const
261 {
262  // Don't assume an ordering on the input data, use the (x, y,
263  // z) values of this centroid to determine the index.
264  unsigned int x_index, y_index, z_index, global_index;
265 
266  x_index = (unsigned int)((p(0) - _minx) / _dx);
267  y_index = (unsigned int)((p(1) - _miny) / _dy);
268 
269  if (_mesh_dimension == 3)
270  {
271  z_index = (unsigned int)((p(2) - _minz) / _dz);
272  global_index = z_index * _ny;
273  }
274  else
275  global_index = 0;
276 
277  // Compute the global index into the _data array. This stores points
278  // in a [z][y][x] ordering.
279  global_index = (global_index + y_index) * _nx + x_index;
280 
281  // Don't access out of range!
282  mooseAssert(global_index < _data.size(),
283  "global_index " << global_index << " points out of _data range: " << _data.size());
284 
285  return global_index;
286 }
287 
288 unsigned int
289 EBSDReader::indexFromIndex(unsigned int var) const
290 {
291 
292  // Transfer the index into the _avg_data array.
293  unsigned avg_index = var;
294 
295  // Don't access out of range!
296  if (avg_index >= _avg_data.size())
297  mooseError("Error! Index out of range in EBSDReader::indexFromIndex(), index: ",
298  avg_index,
299  " size: ",
300  _avg_data.size());
301 
302  return avg_index;
303 }
304 
305 const std::map<dof_id_type, std::vector<Real>> &
307 {
309 }
310 
311 const std::map<dof_id_type, std::vector<Real>> &
313 {
315 }
316 
317 unsigned int
318 EBSDReader::getGlobalID(unsigned int feature_id) const
319 {
320  auto it = _global_id_map.find(feature_id);
321  if (it == _global_id_map.end())
322  mooseError("Invalid Feature ID");
323  return it->second;
324 }
325 
326 void
328 {
329  // maps are only rebuild for use in initial conditions, which happens in time step zero
330  if (_time_step == 0)
332 }
333 
334 void
336 {
337  // Import nodeToElemMap from MooseMesh for current node
338  // This map consists of the node index followed by a vector of element indices that are associated
339  // with that node
340  const std::map<dof_id_type, std::vector<dof_id_type>> & node_to_elem_map =
341  _mesh.nodeToActiveSemilocalElemMap();
342  libMesh::MeshBase & mesh = _mesh.getMesh();
343 
344  // Loop through each node in mesh and calculate eta values for each grain associated with the node
345  MeshBase::const_node_iterator ni = mesh.active_nodes_begin();
346  const MeshBase::const_node_iterator nend = mesh.active_nodes_end();
347  for (; ni != nend; ++ni)
348  {
349  // Get node_id
350  const dof_id_type node_id = (*ni)->id();
351 
352  // Initialize map entries for current node
353  _node_to_grain_weight_map[node_id].assign(getGrainNum(), 0.0);
354  _node_to_phase_weight_map[node_id].assign(getPhaseNum(), 0.0);
355 
356  // Loop through element indices associated with the current node and record weighted eta value
357  // in new map
358  const auto & node_to_elem_pair = node_to_elem_map.find(node_id);
359  if (node_to_elem_pair != node_to_elem_map.end())
360  {
361  unsigned int n_elems =
362  node_to_elem_pair->second
363  .size(); // n_elems can range from 1 to 4 for 2D and 1 to 8 for 3D problems
364 
365  for (unsigned int ne = 0; ne < n_elems; ++ne)
366  {
367  // Current element index
368  unsigned int elem_id = (node_to_elem_pair->second)[ne];
369 
370  // Retrieve EBSD grain number for the current element index
371  const Elem * elem = mesh.elem(elem_id);
372  const EBSDReader::EBSDPointData & d = getData(elem->centroid());
373 
374  // get the (global) grain ID for the EBSD feature ID
375  const unsigned int global_id = getGlobalID(d._feature_id);
376 
377  // Calculate eta value and add to map
378  _node_to_grain_weight_map[node_id][global_id] += 1.0 / n_elems;
379  _node_to_phase_weight_map[node_id][d._phase] += 1.0 / n_elems;
380  }
381  }
382  }
383 }
384 
385 MooseSharedPointer<EBSDAccessFunctors::EBSDPointDataFunctor>
386 EBSDReader::getPointDataAccessFunctor(const MooseEnum & field_name) const
387 {
388  EBSDPointDataFunctor * ret_val = NULL;
389 
390  switch (field_name)
391  {
392  case 0: // phi1
393  ret_val = new EBSDPointDataPhi1();
394  break;
395  case 1: // phi
396  ret_val = new EBSDPointDataPhi();
397  break;
398  case 2: // phi2
399  ret_val = new EBSDPointDataPhi2();
400  break;
401  case 3: // grain
402  ret_val = new EBSDPointDataFeatureID();
403  break;
404  case 4: // phase
405  ret_val = new EBSDPointDataPhase();
406  break;
407  case 5: // symmetry
408  ret_val = new EBSDPointDataSymmetry();
409  break;
410  default:
411  {
412  // check for custom columns
413  for (unsigned int i = 0; i < _custom_columns; ++i)
414  if (field_name == "CUSTOM" + Moose::stringify(i))
415  {
416  ret_val = new EBSDPointDataCustom(i);
417  break;
418  }
419  }
420  }
421 
422  // If ret_val was not set by any of the above cases, throw an error.
423  if (!ret_val)
424  mooseError("Error: Please input supported EBSD_param");
425 
426  // If we made it here, wrap up the the ret_val in a
427  // MooseSharedPointer and ship it off.
428  return MooseSharedPointer<EBSDPointDataFunctor>(ret_val);
429 }
430 
431 MooseSharedPointer<EBSDAccessFunctors::EBSDAvgDataFunctor>
432 EBSDReader::getAvgDataAccessFunctor(const MooseEnum & field_name) const
433 {
434  EBSDAvgDataFunctor * ret_val = NULL;
435 
436  switch (field_name)
437  {
438  case 0: // phi1
439  ret_val = new EBSDAvgDataPhi1();
440  break;
441  case 1: // phi
442  ret_val = new EBSDAvgDataPhi();
443  break;
444  case 2: // phi2
445  ret_val = new EBSDAvgDataPhi2();
446  break;
447  case 3: // phase
448  ret_val = new EBSDAvgDataPhase();
449  break;
450  case 4: // symmetry
451  ret_val = new EBSDAvgDataSymmetry();
452  break;
453  case 5: // local_id
454  ret_val = new EBSDAvgDataLocalID();
455  break;
456  case 6: // feature_id
457  ret_val = new EBSDAvgDataFeatureID();
458  break;
459  default:
460  {
461  // check for custom columns
462  for (unsigned int i = 0; i < _custom_columns; ++i)
463  if (field_name == "CUSTOM" + Moose::stringify(i))
464  {
465  ret_val = new EBSDAvgDataCustom(i);
466  break;
467  }
468  }
469  }
470 
471  // If ret_val was not set by any of the above cases, throw an error.
472  if (!ret_val)
473  mooseError("Error: Please input supported EBSD_param");
474 
475  // If we made it here, wrap up the the ret_val in a
476  // MooseSharedPointer and ship it off.
477  return MooseSharedPointer<EBSDAvgDataFunctor>(ret_val);
478 }
std::map< dof_id_type, std::vector< Real > > _node_to_phase_weight_map
Map of phase weights per node.
Definition: EBSDReader.h:157
unsigned int _custom_columns
number of additional custom data columns
Definition: EBSDReader.h:136
Real _minx
Grid origin.
Definition: EBSDReader.h:172
std::vector< Real > _custom
Custom data columns.
unsigned int getGlobalID(unsigned int phase, unsigned int local_id) const
Return the (global) grain id for a given phase and (local) grain number.
Definition: EBSDReader.h:93
Real _maxx
Maximum grid extent.
Definition: EBSDReader.h:175
unsigned indexFromPoint(const Point &p) const
Computes a global index in the _data array given an input centroid point.
Definition: EBSDReader.C:260
unsigned _nz
Definition: EBSDReader.h:166
virtual const EulerAngles & getEulerAngles(unsigned int) const
EulerAngleProvider interface implementation to fetch a triplet of Euler angles.
Definition: EBSDReader.C:236
Real _miny
Definition: EBSDReader.h:172
unsigned int _grain_num
Variables needed to determine reduced order parameter values.
Definition: EBSDReader.h:129
unsigned _nx
The number of values in the x, y and z directions.
Definition: EBSDReader.h:166
MooseMesh & _mesh
MooseMesh Variables.
Definition: EBSDReader.h:124
Real _minz
Definition: EBSDReader.h:172
const std::map< dof_id_type, std::vector< Real > > & getNodeToGrainWeightMap() const
Returns a map consisting of the node index followd by a vector of all grain weights for that node...
Definition: EBSDReader.C:306
unsigned int _feature_id
EBSD feature id, (gklobal) grain number, symmetry, and phase data.
Point _p
Element centroid position.
InputParameters validParams< EulerAngleProvider >()
Per element EBSD data point.
std::map< unsigned int, unsigned int > _global_id_map
map from feature_id to global_id
Definition: EBSDReader.h:148
const EBSDPointData & getData(const Point &p) const
Get the requested type of data at the point p.
Definition: EBSDReader.C:224
virtual unsigned int getPhaseNum() const
Return the total number of phases.
Definition: EBSDReader.h:74
Access functor base class for EBSDAvgData.
const int & _time_step
current timestep. Maps are only rebuild on mesh change during time step zero
Definition: EBSDReader.h:160
unsigned _ny
Definition: EBSDReader.h:166
Real _dx
The spacing of the values in x, y and z directions.
Definition: EBSDReader.h:169
unsigned int _mesh_dimension
Dimension of the problem domain.
Definition: EBSDReader.h:163
EBSDReader(const InputParameters &params)
Definition: EBSDReader.C:28
virtual unsigned int getGrainNum() const
Return the total number of grains.
Definition: EBSDReader.C:248
MooseSharedPointer< EBSDAvgDataFunctor > getAvgDataAccessFunctor(const MooseEnum &field_name) const
Factory function to return a average functor specified by name.
Definition: EBSDReader.C:432
std::vector< EBSDAvgData > _avg_data
Averages by (global) grain ID.
Definition: EBSDReader.h:142
const EBSDAvgData & getAvgData(unsigned int i) const
Get the requested type of average data for (global) grain number i.
Definition: EBSDReader.C:230
Real _maxz
Definition: EBSDReader.h:175
InputParameters validParams< EBSDReader >()
Definition: EBSDReader.C:18
std::map< dof_id_type, std::vector< Real > > _node_to_grain_weight_map
Map of grain weights per node.
Definition: EBSDReader.h:154
MooseSharedPointer< EBSDPointDataFunctor > getPointDataAccessFunctor(const MooseEnum &field_name) const
Factory function to return a point functor specified by name.
Definition: EBSDReader.C:386
Euler angle triplet.
Definition: EulerAngles.h:20
virtual ~EBSDReader()
Definition: EBSDReader.C:221
Real _maxy
Definition: EBSDReader.h:175
Abstract base class for user objects that implement the Euler Angle provider interface.
std::vector< EulerAngles > _avg_angles
Euler Angles by (global) grain ID.
Definition: EBSDReader.h:145
void buildNodeWeightMaps()
Build grain and phase weight maps.
Definition: EBSDReader.C:335
const EBSDMeshGeometry & getEBSDGeometry() const
Definition: EBSDMesh.h:43
const std::map< dof_id_type, std::vector< Real > > & getNodeToPhaseWeightMap() const
Returns a map consisting of the node index followd by a vector of all phase weights for that node...
Definition: EBSDReader.C:312
Mesh generated from parameters.
Definition: EBSDMesh.h:22
const std::string & getEBSDFilename() const
Definition: EBSDMesh.h:44
void meshChanged()
Maps need to be updated when the mesh changes.
Definition: EBSDReader.C:327
virtual void readFile()
Definition: EBSDReader.C:47
std::vector< std::vector< unsigned int > > _global_id
global ID for given phases and grains
Definition: EBSDReader.h:151
Access functor base class for EBSDPointData.
std::vector< EBSDPointData > _data
Logically three-dimensional data indexed by geometric points in a 1D vector.
Definition: EBSDReader.h:139
unsigned indexFromIndex(unsigned int var) const
Transfer the index into the _avg_data array from given index.
Definition: EBSDReader.C:289