EBSDReader.C

Go to the documentation of this file.

//* This file is part of the MOOSE framework
//* https://www.mooseframework.org
//*
//* All rights reserved, see COPYRIGHT for full restrictions
//* https://github.com/idaholab/moose/blob/master/COPYRIGHT
//*
//* Licensed under LGPL 2.1, please see LICENSE for details
//* https://www.gnu.org/licenses/lgpl-2.1.html
 
#include "EBSDReader.h"
#include "EBSDMesh.h"
#include "MooseMesh.h"
#include "Conversion.h"
#include "NonlinearSystem.h"
#include <Eigen/Dense>
 
#include <fstream>
 
registerMooseObject("PhaseFieldApp", EBSDReader);
 
template <>
InputParameters
validParams<EBSDReader>()
{
  InputParameters params = validParams<EulerAngleProvider>();
  params.addClassDescription("Load and manage DREAM.3D EBSD data files for running simulations on "
                             "reconstructed microstructures.");
  params.addParam<unsigned int>(
      "custom_columns", 0, "Number of additional custom data columns to read from the EBSD file");
  params.addParam<unsigned int>("bins", 20, "Number of bins to segregate quaternions");
  params.addParam<Real>("L_norm", 1, "Specifies the type of average the user intends to perform");
  return params;
}
 
EBSDReader::EBSDReader(const InputParameters & params)
  : EulerAngleProvider(params),
    _mesh(_fe_problem.mesh()),
    _nl(_fe_problem.getNonlinearSystemBase()),
    _grain_num(0),
    _custom_columns(getParam<unsigned int>("custom_columns")),
    _time_step(_fe_problem.timeStep()),
    _mesh_dimension(_mesh.dimension()),
    _bins(getParam<unsigned int>("bins")),
    _L_norm(getParam<Real>("L_norm")),
    _nx(0),
    _ny(0),
    _nz(0),
    _dx(0.),
    _dy(0.),
    _dz(0.)
{
  readFile();
  // throws an error for zero bins
  if (_bins == 0)
    mooseError("One cannot have zero bins");
}
 
void
EBSDReader::readFile()
{
  // No need to re-read data upon recovery
  if (_app.isRecovering())
    return;
 
  // Fetch and check mesh
  EBSDMesh * mesh = dynamic_cast<EBSDMesh *>(&_mesh);
  if (mesh == NULL)
    mooseError("Please use an EBSDMesh in your simulation.");
  std::ifstream stream_in(mesh->getEBSDFilename().c_str());
  if (!stream_in)
    mooseError("Can't open EBSD file: ", mesh->getEBSDFilename());
 
  const EBSDMesh::EBSDMeshGeometry & g = mesh->getEBSDGeometry();
  // Copy file header data from the EBSDMesh
  _dx = g.d[0];
  _nx = g.n[0];
  _minx = g.min[0];
  _maxx = _minx + _dx * _nx;
 
  _dy = g.d[1];
  _ny = g.n[1];
  _miny = g.min[1];
  _maxy = _miny + _dy * _ny;
 
  _dz = g.d[2];
  _nz = g.n[2];
  _minz = g.min[2];
  _maxz = _minz + _dz * _nz;
 
  // Resize the _data array
  unsigned total_size = g.dim < 3 ? _nx * _ny : _nx * _ny * _nz;
  _data.resize(total_size);
 
  std::string line;
  while (std::getline(stream_in, line))
  {
    if (line.find("#") != 0)
    {
      // Temporary variables to read in on each line
      EBSDPointData d;
      Real x, y, z;
 
      std::istringstream iss(line);
      iss >> d._phi1 >> d._Phi >> d._phi2 >> x >> y >> z >> d._feature_id >> d._phase >>
          d._symmetry;
 
      // Transform angles to degrees
      d._phi1 *= 180.0 / libMesh::pi;
      d._Phi *= 180.0 / libMesh::pi;
      d._phi2 *= 180.0 / libMesh::pi;
 
      // Custom columns
      d._custom.resize(_custom_columns);
      for (unsigned int i = 0; i < _custom_columns; ++i)
        if (!(iss >> d._custom[i]))
          mooseError("Unable to read in EBSD custom data column #", i);
 
      if (x < _minx || y < _miny || x > _maxx || y > _maxy ||
          (g.dim == 3 && (z < _minz || z > _maxz)))
        mooseError("EBSD Data ouside of the domain declared in the header ([",
                   _minx,
                   ':',
                   _maxx,
                   "], [",
                   _miny,
                   ':',
                   _maxy,
                   "], [",
                   _minz,
                   ':',
                   _maxz,
                   "]) dim=",
                   g.dim,
                   "\n",
                   line);
 
      d._p = Point(x, y, z);
 
      // determine number of grains in the dataset
      if (_global_id_map.find(d._feature_id) == _global_id_map.end())
        _global_id_map[d._feature_id] = _grain_num++;
 
      unsigned int global_index = indexFromPoint(Point(x, y, z));
      _data[global_index] = d;
    }
  }
  stream_in.close();
 
  // Resize the variables
  _avg_data.resize(_grain_num);
  _avg_angles.resize(_grain_num);
 
  // clear the averages
  for (unsigned int i = 0; i < _grain_num; ++i)
  {
    EBSDAvgData & a = _avg_data[i];
    a._symmetry = a._phase = a._n = 0;
    a._p = 0.0;
    a._custom.assign(_custom_columns, 0.0);
 
    EulerAngles & b = _avg_angles[i];
    b.phi1 = b.Phi = b.phi2 = 0.0;
  }
 
  // Array of vectors to store quaternions of each grain
  std::vector<std::vector<Eigen::Quaternion<Real>>> quat(_grain_num);
 
  // Iterate through data points to store orientation information for each grain
  for (auto & j : _data)
  {
    EBSDAvgData & a = _avg_data[_global_id_map[j._feature_id]];
    EulerAngles angles;
 
    angles.phi1 = j._phi1;
    angles.Phi = j._Phi;
    angles.phi2 = j._phi2;
 
    // convert Euler angles to quaternions
    Eigen::Quaternion<Real> q = angles.toQuaternion();
 
    quat[_global_id_map[j._feature_id]].push_back(q);
 
    if (a._n == 0)
      a._phase = j._phase;
    else if (a._phase != j._phase)
      mooseError("An EBSD feature needs to have a uniform phase.");
 
    if (a._n == 0)
      a._symmetry = j._symmetry;
    else if (a._symmetry != j._symmetry)
      mooseError("An EBSD feature needs to have a uniform symmetry parameter.");
 
    for (unsigned int i = 0; i < _custom_columns; ++i)
      a._custom[i] += j._custom[i];
 
    // store the feature (or grain) ID
    a._feature_id = j._feature_id;
 
    a._p += j._p;
    a._n++;
  }
  // creating a vector of maps to store the quaternion count for each bin index
  std::vector<
      std::map<std::tuple<unsigned int, unsigned int, unsigned int, unsigned int>, unsigned int>>
      feature_weights(_grain_num);
  // creating a vector of vectors to store the bin index corresponding to a quaternion
  std::vector<std::vector<std::tuple<unsigned int, unsigned int, unsigned int, unsigned int>>>
      quat_to_bin(_grain_num);
 
  for (unsigned int i = 0; i < _grain_num; i++)
  {
    EBSDAvgData & a = _avg_data[i];
    EulerAngles & b = _avg_angles[i];
 
    if (a._n == 0)
      continue;
 
    Real w_index, x_index, y_index, z_index;
    Real w, x, y, z;
    std::vector<std::tuple<double, Eigen::VectorXd>> eigen_vectors_and_values;
 
    // creating an empty map
    feature_weights.push_back(
        std::map<std::tuple<unsigned int, unsigned int, unsigned int, unsigned int>,
                 unsigned int>());
 
    // looping through quaternions of each grain
    for (unsigned int k = 0; k < quat[i].size(); k++)
    {
      std::tuple<unsigned int, unsigned int, unsigned int, unsigned int> bin;
      w_index = int((quat[i][k].w() + 1) * 0.5 * _bins);
      x_index = int((quat[i][k].x() + 1) * 0.5 * _bins);
      y_index = int((quat[i][k].y() + 1) * 0.5 * _bins);
      z_index = int((quat[i][k].z() + 1) * 0.5 * _bins);
 
      bin = std::make_tuple(w_index, x_index, y_index, z_index);
 
      // storing the bincorresponding to the quaternion quat[i][k]
      quat_to_bin[i].push_back(bin);
 
      // storing the bin name and its corresponding size value
      auto j = feature_weights[i].find(bin);
      if (j == feature_weights[i].end())
        feature_weights[i].emplace_hint(j, bin, 1);
      else
        j->second++;
    }
 
    // creating a quaternion matrix
    Eigen::MatrixXd quat_mat(4, quat[i].size());
    // creating the weight matrix
    Eigen::MatrixXd weight = Eigen::MatrixXd::Identity(quat[i].size(), quat[i].size());
 
    unsigned int bin_size;
    bool data_quality_ok = false;
    Real total_weight = 0.0;
 
    for (unsigned int j = 0; j < quat[i].size(); j++)
    {
      bin_size = feature_weights[i][quat_to_bin[i][j]];
      w = quat[i][j].w();
      x = quat[i][j].x();
      y = quat[i][j].y();
      z = quat[i][j].z();
 
      // instantiating columns of the matrix
      quat_mat.col(j) << w, x, y, z;
      // assigning weights to each quaternion
      weight(j, j) = std::pow(bin_size, _L_norm);
      total_weight += weight(j, j);
      if (bin_size / quat[i].size() > 0.5)
        data_quality_ok = true;
    }
 
    weight = weight / total_weight;
 
    // throws a warning if EBSD data is not reliable
    if (!data_quality_ok)
      _console << COLOR_YELLOW << "EBSD data may not be reliable"
               << "\n"
               << COLOR_DEFAULT;
 
    // compute eigen values and eigen vectors
    Eigen::EigenSolver<Eigen::MatrixXd> EigenSolver(quat_mat * weight * quat_mat.transpose());
    Eigen::VectorXd eigen_values = EigenSolver.eigenvalues().real();
    Eigen::MatrixXd eigen_vectors = EigenSolver.eigenvectors().real();
 
    // creating a tuple of eigen values and eigen vectors
    for (unsigned int i = 0; i < eigen_values.size(); i++)
    {
      std::tuple<double, Eigen::VectorXd> vec_and_val(eigen_values[i], eigen_vectors.col(i));
      eigen_vectors_and_values.push_back(vec_and_val);
    }
    // Sorting the eigen values and vectors in ascending order of eigen values
    auto j = std::max_element(eigen_vectors_and_values.begin(),
                              eigen_vectors_and_values.end(),
                              [&](const std::tuple<double, Eigen::VectorXd> & a,
                                  const std::tuple<double, Eigen::VectorXd> & b) -> bool {
                                return std::get<0>(a) < std::get<0>(b);
                              });
    // Selecting eigen vector corresponding to max eigen value to compute average euler angle
    Eigen::Quaternion<Real> q(
        std::get<1>(*j)(0), std::get<1>(*j)(1), std::get<1>(*j)(2), std::get<1>(*j)(3));
 
    b = EulerAngles(q);
 
    // link the EulerAngles into the EBSDAvgData for access via the functors
    a._angles = &b;
 
    if (a._phase >= _global_id.size())
      _global_id.resize(a._phase + 1);
 
    // The averaged per grain data locally contains the phase id, local id, and
    // original feature id. It is stored contiguously indexed by global id.
    a._local_id = _global_id[a._phase].size();
    _global_id[a._phase].push_back(i);
 
    a._p *= 1.0 / Real(a._n);
 
    for (unsigned int i = 0; i < _custom_columns; ++i)
      a._custom[i] /= Real(a._n);
  }
  // Build maps to indicate the weights with which grain and phase data
  // from the surrounding elements contributes to a node fo IC purposes
  buildNodeWeightMaps();
}
 
EBSDReader::~EBSDReader() {}
 
const EBSDReader::EBSDPointData &
EBSDReader::getData(const Point & p) const
{
  return _data[indexFromPoint(p)];
}
 
const EBSDReader::EBSDAvgData &
EBSDReader::getAvgData(unsigned int var) const
{
  return _avg_data[indexFromIndex(var)];
}
 
const EulerAngles &
EBSDReader::getEulerAngles(unsigned int var) const
{
  return _avg_angles[indexFromIndex(var)];
}
 
const EBSDReader::EBSDAvgData &
EBSDReader::getAvgData(unsigned int phase, unsigned int local_id) const
{
  return _avg_data[indexFromIndex(_global_id[phase][local_id])];
}
 
unsigned int
EBSDReader::getGrainNum() const
{
  return _grain_num;
}
 
unsigned int
EBSDReader::getGrainNum(unsigned int phase) const
{
  return _global_id[phase].size();
}
 
unsigned int
EBSDReader::indexFromPoint(const Point & p) const
{
  // Don't assume an ordering on the input data, use the (x, y,
  // z) values of this centroid to determine the index.
  unsigned int x_index, y_index, z_index, global_index;
 
  x_index = (unsigned int)((p(0) - _minx) / _dx);
  y_index = (unsigned int)((p(1) - _miny) / _dy);
  if (p(0) <= _minx || p(0) >= _maxx || p(1) <= _miny || p(1) >= _maxy)
    mooseError("Data points must be on the interior of the mesh elements. In EBSDReader ", name());
 
  if (_mesh_dimension == 3)
  {
    z_index = (unsigned int)((p(2) - _minz) / _dz);
    global_index = z_index * _ny;
    if (p(2) <= _minz || p(2) >= _maxz)
      mooseError("Data points must be on the interior of the mesh elements. In EBSDReader ",
                 name());
  }
  else
    global_index = 0;
 
  // Compute the global index into the _data array.  This stores points
  // in a [z][y][x] ordering.
  global_index = (global_index + y_index) * _nx + x_index;
 
  // Don't access out of range!
  mooseAssert(global_index < _data.size(),
              "global_index " << global_index << " points out of _data range: " << _data.size());
 
  return global_index;
}
 
unsigned int
EBSDReader::indexFromIndex(unsigned int var) const
{
 
  // Transfer the index into the _avg_data array.
  unsigned avg_index = var;
 
  // Don't access out of range!
  if (avg_index >= _avg_data.size())
    mooseError("Error! Index out of range in EBSDReader::indexFromIndex(), index: ",
               avg_index,
               " size: ",
               _avg_data.size());
 
  return avg_index;
}
 
const std::map<dof_id_type, std::vector<Real>> &
EBSDReader::getNodeToGrainWeightMap() const
{
  return _node_to_grain_weight_map;
}
 
const std::map<dof_id_type, std::vector<Real>> &
EBSDReader::getNodeToPhaseWeightMap() const
{
  return _node_to_phase_weight_map;
}
 
unsigned int
EBSDReader::getGlobalID(unsigned int feature_id) const
{
  auto it = _global_id_map.find(feature_id);
  if (it == _global_id_map.end())
    mooseError("Invalid Feature ID");
  return it->second;
}
 
void
EBSDReader::meshChanged()
{
  // maps are only rebuild for use in initial conditions, which happens in time step zero
  if (_time_step == 0)
    buildNodeWeightMaps();
}
 
void
EBSDReader::buildNodeWeightMaps()
{
  // Import nodeToElemMap from MooseMesh for current node
  // This map consists of the node index followed by a vector of element indices that are associated
  // with that node
  const std::map<dof_id_type, std::vector<dof_id_type>> & node_to_elem_map =
      _mesh.nodeToActiveSemilocalElemMap();
  libMesh::MeshBase & mesh = _mesh.getMesh();
 
  // Loop through each node in mesh and calculate eta values for each grain associated with the node
  for (const auto & node : as_range(mesh.active_nodes_begin(), mesh.active_nodes_end()))
  {
    // Get node_id
    const dof_id_type node_id = node->id();
 
    // Initialize map entries for current node
    _node_to_grain_weight_map[node_id].assign(getGrainNum(), 0.0);
    _node_to_phase_weight_map[node_id].assign(getPhaseNum(), 0.0);
 
    // Loop through element indices associated with the current node and record weighted eta value
    // in new map
    const auto & node_to_elem_pair = node_to_elem_map.find(node_id);
    if (node_to_elem_pair != node_to_elem_map.end())
    {
      unsigned int n_elems =
          node_to_elem_pair->second
              .size(); // n_elems can range from 1 to 4 for 2D and 1 to 8 for 3D problems
 
      for (unsigned int ne = 0; ne < n_elems; ++ne)
      {
        // Current element index
        unsigned int elem_id = (node_to_elem_pair->second)[ne];
 
        // Retrieve EBSD grain number for the current element index
        const Elem * elem = mesh.elem_ptr(elem_id);
        const EBSDReader::EBSDPointData & d = getData(elem->centroid());
 
        // get the (global) grain ID for the EBSD feature ID
        const unsigned int global_id = getGlobalID(d._feature_id);
 
        // Calculate eta value and add to map
        _node_to_grain_weight_map[node_id][global_id] += 1.0 / n_elems;
        _node_to_phase_weight_map[node_id][d._phase] += 1.0 / n_elems;
      }
    }
  }
}
 
MooseSharedPointer<EBSDAccessFunctors::EBSDPointDataFunctor>
EBSDReader::getPointDataAccessFunctor(const MooseEnum & field_name) const
{
  EBSDPointDataFunctor * ret_val = NULL;
 
  switch (field_name)
  {
    case 0: // phi1
      ret_val = new EBSDPointDataPhi1();
      break;
    case 1: // phi
      ret_val = new EBSDPointDataPhi();
      break;
    case 2: // phi2
      ret_val = new EBSDPointDataPhi2();
      break;
    case 3: // grain
      ret_val = new EBSDPointDataFeatureID();
      break;
    case 4: // phase
      ret_val = new EBSDPointDataPhase();
      break;
    case 5: // symmetry
      ret_val = new EBSDPointDataSymmetry();
      break;
    default:
    {
      // check for custom columns
      for (unsigned int i = 0; i < _custom_columns; ++i)
        if (field_name == "CUSTOM" + Moose::stringify(i))
        {
          ret_val = new EBSDPointDataCustom(i);
          break;
        }
    }
  }
 
  // If ret_val was not set by any of the above cases, throw an error.
  if (!ret_val)
    mooseError("Error:  Please input supported EBSD_param");
 
  // If we made it here, wrap up the the ret_val in a
  // MooseSharedPointer and ship it off.
  return MooseSharedPointer<EBSDPointDataFunctor>(ret_val);
}
 
MooseSharedPointer<EBSDAccessFunctors::EBSDAvgDataFunctor>
EBSDReader::getAvgDataAccessFunctor(const MooseEnum & field_name) const
{
  EBSDAvgDataFunctor * ret_val = NULL;
 
  switch (field_name)
  {
    case 0: // phi1
      ret_val = new EBSDAvgDataPhi1();
      break;
    case 1: // phi
      ret_val = new EBSDAvgDataPhi();
      break;
    case 2: // phi2
      ret_val = new EBSDAvgDataPhi2();
      break;
    case 3: // phase
      ret_val = new EBSDAvgDataPhase();
      break;
    case 4: // symmetry
      ret_val = new EBSDAvgDataSymmetry();
      break;
    case 5: // local_id
      ret_val = new EBSDAvgDataLocalID();
      break;
    case 6: // feature_id
      ret_val = new EBSDAvgDataFeatureID();
      break;
    default:
    {
      // check for custom columns
      for (unsigned int i = 0; i < _custom_columns; ++i)
        if (field_name == "CUSTOM" + Moose::stringify(i))
        {
          ret_val = new EBSDAvgDataCustom(i);
          break;
        }
    }
  }
 
  // If ret_val was not set by any of the above cases, throw an error.
  if (!ret_val)
    mooseError("Error:  Please input supported EBSD_param");
 
  // If we made it here, wrap up the the ret_val in a
  // MooseSharedPointer and ship it off.
  return MooseSharedPointer<EBSDAvgDataFunctor>(ret_val);
}