Merge pull request #5929 from dschwen/ebsd_5502

Add phase selection to ReconVarIC

modules/phase_field/include/ics/ReconVarIC.h

@@ -41,7 +41,6 @@ class ReconVarIC : public InitialCondition
   unsigned int _op_index;
 
   unsigned int _grain_num;
-  int _grain_index_offset;
 
   std::vector<Point> _centerpoints;
   std::vector<Real> _assigned_op;

modules/phase_field/include/userobjects/EBSDReader.h

@@ -18,6 +18,14 @@ InputParameters validParams<EBSDReader>();
 /**
  * A GeneralUserObject that reads an EBSD file and stores the centroid
  * data in a data structure which indexes on element centroids.
+ *
+ * Grains are indexed through multiple schemes:
+ *  * feature_id   The grain number in the EBSD data file
+ *  * global_id    The index into the global average data (this is feature_id shifted back by feature_id_origin)
+ *  * local_id     The index into the per-phase grain list. This is only unique when combined with phase number
+ *
+ * Phases are referred to using the numbers in the EBSD data file. In case the phase number in the data file
+ * starts at 1 the phase 0 will simply contain no grains.
  */
 class EBSDReader : public EulerAngleProvider, public EBSDAccessFunctors
 {
@@ -55,7 +63,7 @@ class EBSDReader : public EulerAngleProvider, public EBSDAccessFunctors
   /**
    * Get the requested type of average data for a given phase and grain.
    */
-  const EBSDAvgData &  getAvgData(unsigned int phase, unsigned int grain) const;
+  const EBSDAvgData &  getAvgData(unsigned int phase, unsigned int local_id) const;
 
   /**
    * EulerAngleProvider interface implementation to fetch a triplet of Euler angles
@@ -70,13 +78,23 @@ class EBSDReader : public EulerAngleProvider, public EBSDAccessFunctors
   /**
    * Return the total number of phases
    */
-  virtual unsigned int getPhaseNum() const { return _feature_id.size(); }
+  virtual unsigned int getPhaseNum() const { return _global_id.size(); }
 
   /**
    * Return the number of grains in a given phase
    */
   unsigned int getGrainNum(unsigned int phase) const;
 
+  /**
+   * Return the feature id (global grain number) for a given phase and phase grain number
+   */
+  unsigned int getFeatureID(unsigned int phase, unsigned int local_id) const { return _avg_data[_global_id[phase][local_id]].grain; }
+
+  /**
+   * Return the feature id (global grain number) for a given phase and phase grain number
+   */
+  unsigned int getGlobalID(unsigned int phase, unsigned int local_id) const { return _global_id[phase][local_id]; }
+
   /// Factory function to return a point functor specified by name
   EBSDPointDataFunctor * getPointDataAccessFunctor(const MooseEnum & field_name) const;
   /// Factory function to return a average functor specified by name
@@ -117,8 +135,11 @@ class EBSDReader : public EulerAngleProvider, public EBSDAccessFunctors
   /// Euler Angles by feature ID
   std::vector<EulerAngles> _avg_angles;
 
-  /// feature ID for given phases and grains
-  std::vector<std::vector<unsigned int> > _feature_id;
+  /// map from feature_id to global_id
+  std::map<unsigned int, unsigned int> _global_id_map;
+
+  /// global ID for given phases and grains
+  std::vector<std::vector<unsigned int> > _global_id;
 
   /// Map of grain weights per node
   std::map<dof_id_type, std::vector<Real> > _node_to_grain_weight_map;

modules/phase_field/include/utils/EBSDAccessFunctors.h

@@ -16,15 +16,15 @@ class EBSDAccessFunctors
   /// Per element EBSD data point
   struct EBSDPointData {
     Real phi1, phi, phi2, symmetry;
-    unsigned int grain, phase, op;
+    unsigned int grain, phase, op, global;
     Point p;
     std::vector<Real> custom;
   };
 
   /// Averaged EBSD data
   struct EBSDAvgData {
     EulerAngles * angles;
-    unsigned int phase, symmetry, grain, n;
+    unsigned int phase, local, symmetry, grain, n;
     Point p;
     std::vector<Real> custom;
   };
@@ -84,6 +84,12 @@ class EBSDAccessFunctors
   struct EBSDAvgDataPhase : EBSDAvgDataFunctor {
     virtual Real operator () (const EBSDAvgData & a) { return a.phase; };
   };
+  struct EBSDAvgDataLocalID : EBSDAvgDataFunctor {
+    virtual Real operator () (const EBSDAvgData & a) { return a.local; };
+  };
+  struct EBSDAvgDataGrain : EBSDAvgDataFunctor {
+    virtual Real operator () (const EBSDAvgData & a) { return a.grain; };
+  };
   struct EBSDAvgDataSymmetry : EBSDAvgDataFunctor {
     virtual Real operator () (const EBSDAvgData & a) { return a.symmetry; };
   };

modules/phase_field/src/ics/ReconVarIC.C

@@ -28,23 +28,18 @@ ReconVarIC::ReconVarIC(const InputParameters & parameters) :
     _op_index(getParam<unsigned int>("op_index")),
     _node_to_grain_weight_map(_ebsd_reader.getNodeToGrainWeightMap())
 {
-  if (_consider_phase)
-    mooseError("In ReconvarIC, setting ics from grain(phase) information is not currently supported");
 }
 
 void
 ReconVarIC::initialSetup()
 {
-  // grain index number offset
-  _grain_index_offset = _consider_phase ? 1 : 0;
-
   // number of grains this ICs deals with
   _grain_num = _consider_phase ?_ebsd_reader.getGrainNum(_phase) : _ebsd_reader.getGrainNum();
 
   // fetch all center points
   _centerpoints.resize(_grain_num);
-  for (unsigned int grain = 0; grain < _grain_num; ++grain)
-    _centerpoints[grain] = getCenterPoint(grain);
+  for (unsigned int index = 0; index < _grain_num; ++index)
+    _centerpoints[index] = getCenterPoint(index);
 
   // We do not want to have more order parameters than grains. That would leave some unused.
   if (_op_num > _grain_num)
@@ -71,26 +66,24 @@ ReconVarIC::value(const Point & /*p*/)
   //Get local information from EBSDReader object
   EBSDReader::EBSDPointData local_ebsd_data; // = _ebsd_reader.getData(p);
 
-  // Increment through all grains at node_index
-  for (unsigned int grain = 0; grain < _grain_num; ++grain)
+  // Increment through all grains at node_index (these are global IDs if consider_phase is false and local IDs otherwise)
+  for (unsigned int index = 0; index < _grain_num; ++index)
   {
     // If the current order parameter index (_op_index) is equal to the assinged index (_assigned_op),
-    // set the value from node_to_grain_weight_map
-    if (!_consider_phase)
-      if (_assigned_op[grain] == _op_index && (it->second)[grain] > 0.0)
-        return (it->second)[grain];
-
-    //TODO: Make this work when considering phase
+    // set the value from node_to_grain_weight_map=
+    Real value = (it->second)[_consider_phase ? _ebsd_reader.getGlobalID(_phase, index) : index];
+    if (_assigned_op[index] == _op_index && value > 0.0)
+      return value;
   }
 
   return 0.0;
 }
 
 Point
-ReconVarIC::getCenterPoint(unsigned int grain)
+ReconVarIC::getCenterPoint(unsigned int index)
 {
   if (_consider_phase)
-    return _ebsd_reader.getAvgData(_phase, grain + _grain_index_offset).p;
+    return _ebsd_reader.getAvgData(_phase, index).p;
   else
-    return _ebsd_reader.getAvgData(grain + _grain_index_offset).p;
+    return _ebsd_reader.getAvgData(index).p;
 }

modules/phase_field/src/userobjects/EBSDReader.C

@@ -101,21 +101,20 @@ EBSDReader::readFile()
       d.p = Point(x, y, z);
 
       // determine number of grains in the dataset
-      if (d.grain > _feature_num) _feature_num = d.grain;
+      if (_global_id_map.find(d.grain) == _global_id_map.end())
+        _global_id_map[d.grain] = _feature_num++;
+      d.global = _global_id_map[d.grain];
 
       // The Order parameter is not yet assigned.
       // We initialize it to zero in order not to have undefined values that break the testing.
       d.op = 0;
 
-      unsigned global_index = indexFromPoint(Point(x, y, z));
+      unsigned int global_index = indexFromPoint(Point(x, y, z));
       _data[global_index] = d;
     }
   }
   stream_in.close();
 
-  // total number of grains is one higher than the maximum grain id
-  _feature_num += 1;
-
   // Resize the variables
   _avg_data.resize(_feature_num);
   _avg_angles.resize(_feature_num);
@@ -135,8 +134,8 @@ EBSDReader::readFile()
   // Iterate through data points to get average variable values for each grain
   for (std::vector<EBSDPointData>::iterator j = _data.begin(); j != _data.end(); ++j)
   {
-    EBSDAvgData & a = _avg_data[j->grain];
-    EulerAngles & b = _avg_angles[j->grain];
+    EBSDAvgData & a = _avg_data[_global_id_map[j->grain]];
+    EulerAngles & b = _avg_angles[_global_id_map[j->grain]];
 
     //use Eigen::Quaternion<Real> here?
     b.phi1 += j->phi1;
@@ -158,6 +157,9 @@ EBSDReader::readFile()
     for (unsigned int i = 0; i < _custom_columns; ++i)
       a.custom[i] += j->custom[i];
 
+    // store the feature (or grain) ID
+    a.grain = j->grain;
+
     a.p += j->p;
     a.n++;
   }
@@ -169,18 +171,21 @@ EBSDReader::readFile()
 
     if (a.n == 0) continue;
 
+    // TODO: need betetr way to average angles
     b.phi1 /= Real(a.n);
     b.Phi  /= Real(a.n);
     b.phi2 /= Real(a.n);
 
     // link the EulerAngles into the EBSDAvgData for access via the functors
     a.angles = &b;
 
-    if (a.phase >= _feature_id.size())
-      _feature_id.resize(a.phase + 1);
+    if (a.phase >= _global_id.size())
+      _global_id.resize(a.phase + 1);
 
-    a.grain = _feature_id[a.phase].size();
-    _feature_id[a.phase].push_back(i);
+    // 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 = _global_id[a.phase].size();
+    _global_id[a.phase].push_back(i);
 
     a.p *= 1.0/Real(a.n);
 
@@ -216,9 +221,9 @@ EBSDReader::getEulerAngles(unsigned int var) const
 }
 
 const EBSDReader::EBSDAvgData &
-EBSDReader::getAvgData(unsigned int phase, unsigned int grain) const
+EBSDReader::getAvgData(unsigned int phase, unsigned int local_id) const
 {
-  return _avg_data[indexFromIndex(_feature_id[phase][grain])];
+  return _avg_data[indexFromIndex(_global_id[phase][local_id])];
 }
 
 unsigned int
@@ -230,7 +235,7 @@ EBSDReader::getGrainNum() const
 unsigned int
 EBSDReader::getGrainNum(unsigned int phase) const
 {
-  return _feature_id[phase].size();
+  return _global_id[phase].size();
 }
 
 unsigned int
@@ -320,7 +325,7 @@ EBSDReader::buildNodeWeightMaps()
       const EBSDReader::EBSDPointData & d = getData(elem->centroid());
 
       // Calculate eta value and add to map
-      _node_to_grain_weight_map[node_id][d.grain] += 1.0 / n_elems;
+      _node_to_grain_weight_map[node_id][d.global] += 1.0 / n_elems;
       _node_to_phase_weight_map[node_id][d.phase] += 1.0 / n_elems;
     }
   }

modules/phase_field/src/utils/EBSDAccessFunctors.C

@@ -9,5 +9,5 @@ EBSDAccessFunctors::getPointDataFieldType()
 MooseEnum
 EBSDAccessFunctors::getAvgDataFieldType()
 {
-  return MooseEnum("phi1 phi phi2 phase symmetry", "", true);
+  return MooseEnum("phi1 phi phi2 phase symmetry local grain", "", true);
 }

modules/phase_field/tests/reconstruction_test/1phase_reconstruction_test.i

@@ -1,3 +1,8 @@
+#
+# In this test we set the initial condition of a set of order parameters
+# by pulling out the grain data from given EBSD data file ignoring the phase completely.
+#
+
 [Problem]
   type = FEProblem
   solve = false

modules/phase_field/tests/reconstruction_test/2phase_reconstruction_test.i

@@ -1,3 +1,9 @@
+#
+# In this test we set the initial condition of two variables
+# based on solely the phase information in a given EBSD data file,
+# ignoring the feature IDs entirely
+#
+
 [Problem]
   type = FEProblem
   solve = false

modules/phase_field/tests/reconstruction_test/2phase_reconstruction_test2.i

@@ -0,0 +1,49 @@
+#
+# In this test we set the initial condition of a set of order parameters
+# by pulling out the only grains from given EBSD data file that belong to a specified phase
+#
+
+[Problem]
+  type = FEProblem
+  solve = false
+  kernel_coverage_check = false
+[]
+
+[Mesh]
+  type = EBSDMesh
+  filename = Ti_2Phase_28x28_Sqr_Marmot.txt
+[]
+
+[GlobalParams]
+  op_num = 3
+  var_name_base = gr
+[]
+
+[UserObjects]
+  [./ebsd]
+    type = EBSDReader
+  [../]
+[]
+
+[ICs]
+  [./PolycrystalICs]
+    [./ReconVarIC]
+      ebsd_reader = ebsd
+      phase = 1
+    [../]
+  [../]
+[]
+
+[Variables]
+  [./PolycrystalVariables]
+  [../]
+[]
+
+[Executioner]
+  type = Transient
+  num_steps = 0
+[]
+
+[Outputs]
+  exodus = true
+[]