VV: Compute Schmids Fully V&V'ed

docs/vv_templates/report_template.md

@@ -1,5 +1,7 @@
 # V&V Report: <FilterName>
 
+<!-- Write the report in ASD-STE100 Simplified Technical English. Use short sentences, active voice, one meaning per sentence, and consistent technical terms. -->
+
 |        |              |
 |--------|--------------|
 | Plugin | <PluginName> |
@@ -46,6 +48,17 @@ For worked instances see `src/Plugins/OrientationAnalysis/vv/BadDataNeighborOrie
 
 *Second-engineer review:* *<name, date>* OR *skipped — reason.*
 
+## Bugs found and fixed
+
+<!-- Include only confirmed defects that are fixed in the verified branch. Use the stable `<FilterName>-D<N>` deviation IDs. Do not use `SC-*` labels. State all affected released versions. Use `docs/dream3d_nx_release_dates.md` to identify the DREAM3D-NX releases. If the next release number is not known, state that the fix will be in the release after the latest affected release. If the V&V found no defects, replace this table with `None.` -->
+
+*This branch fixes all defects in this table. The fixes will be in the DREAM3D-NX release after v<latest affected release>.*
+
+| Deviation | Defect | Affected released versions | Resolution in this branch |
+|---|---|---|---|
+| `<FilterName>-D1` | *State the defect and its effect.* | *DREAM.3D 6.5.171; DREAM3D-NX v<first> through v<last>.* | *State how the verified branch fixes the defect.* |
+| `<FilterName>-D2` | *State the defect and its effect.* | *DREAM.3D 6.5.171 only. DREAM3D-NX was not affected.* | *State how the verified branch fixes or prevents the defect.* |
+
 ## Code path coverage
 
 *N of M paths exercised. If N < M: which paths are NOT covered, and why each gap is acceptable (or what would close it).*
@@ -78,13 +91,17 @@ When a path is intentionally not covered, write `*Not directly tested. <one-line
 - **SHA512:** *`<copy from test/CMakeLists.txt>`*
 - **Provenance:** *`src/Plugins/<P>/vv/provenance/<name>.md`*
 
-## Deviations from DREAM3D 6.5.171
+## Deviations from DREAM.3D 6.5.171
 
 *Either:*
 
 - No deviations observed. Comparison run on *<fixture>*.
 
 *Or:*
 
-- `<FilterName>-D1` — *one-line symptom* — see `vv/deviations/<FilterName>.md`
-- `<FilterName>-D2` — *one-line symptom* — see `vv/deviations/<FilterName>.md`
+See `vv/deviations/<FilterName>.md` for the root cause, affected users, and recommendation for each deviation.
+
+| Deviation | Observed difference |
+|---|---|
+| `<FilterName>-D1` | *State the user-visible difference.* |
+| `<FilterName>-D2` | *State the user-visible difference.* |

src/Plugins/OrientationAnalysis/docs/ComputeSchmidsFilter.md

@@ -6,47 +6,68 @@ Statistics (Crystallography)
 
 ## Description
 
-This **Filter** calculates the Schmid factor for each **Feature** (grain), which measures how favorably oriented that grain is for plastic deformation under a given loading direction. A higher Schmid factor means the grain is more likely to deform (slip) under the applied load.
+This **Filter** calculates the Schmid factor of each **Feature** given its average orientation and a user-defined loading axis. The Schmid factor combines the component of the axial force *F* that lies parallel to the slip direction with the component perpendicular to the slip plane. The equation is:
 
-### What is the Schmid Factor?
+Schmid Factor = (cos φ cos λ)
 
-When a force is applied to a polycrystalline material, each grain responds differently depending on how its internal crystal planes are oriented relative to the loading direction. Deformation occurs by *slip* -- atoms sliding along specific crystal planes in specific directions. Each combination of a slip plane and a slip direction is called a *slip system*.
+*The angle φ is the angle between the tensile axis and the slip plane normal, and λ is the angle between the tensile axis and the slip direction in the slip plane.*
 
-The Schmid factor quantifies how much of the applied force is resolved onto a given slip system:
+The **Filter** evaluates this equation for every enumerated slip system in each **Feature's** crystal structure. It stores the largest Schmid factor and the corresponding slip-system index. Critical resolved shear stress is not considered.
 
-**Schmid Factor = cos(φ) × cos(λ)**
+The user-supplied *Loading Direction* is normalized before use, so only its direction matters — `[1, 2, 3]` and `[3, 6, 9]` give identical results. It is a **sample-frame** direction; the **Filter** rotates it into each **Feature's** crystal frame using that **Feature's** average orientation before evaluating the slip systems. The loading direction must be non-zero. When *Override Default Slip System* is enabled, the supplied slip-plane normal and slip direction must also be non-zero and perpendicular.
 
-where:
+**Feature 0** is the conventional "unassigned" **Feature** and is never computed; all of its output values are zero. A **Feature** whose crystal-structure value is outside EbsdLib's supported Laue-group range is skipped, and all of its output values are also zero. For a valid Laue class that has no enumerated slip systems, *Schmids*, *Slip Systems*, *Phis*, and *Lambdas* are zero, while *Poles* still records the transformed loading direction.
 
-- **φ** is the angle between the loading axis and the *slip plane normal* (the direction perpendicular to the slip plane)
-- **λ** is the angle between the loading axis and the *slip direction* (the direction atoms slide along within the slip plane)
+### Ties Between Slip Systems
 
-The Schmid factor ranges from 0 to 0.5. A value of 0.5 means the slip system is optimally aligned with the applied load. A value near 0 means the slip system is poorly oriented for deformation.
+Candidate slip systems are compared with a strict greater-than, in the fixed order in which they are enumerated for the crystal structure. When two or more slip systems share the maximum Schmid factor — which happens for high-symmetry loading directions such as `[0, 0, 1]` or `[1, 1, 1]` in a cubic crystal — the **lowest-numbered** of the tied systems is reported. This is a property of the enumeration order, not a physical preference: the tied systems are equally favoured.
 
-![Fig. 1: The geometric relationship between the tensile axis, slip plane, and slip direction that defines the Schmid factor.](Images/ComputeSchmids_SchmidFactor.png)
+### Slip Systems Numbering Depends on Override Default Slip System
 
-### How This Filter Works
+The meaning of the *Slip Systems* output changes with the *Override Default Slip System* toggle:
 
-1. The user specifies a **Loading Direction** as a unit vector in the sample reference frame (e.g., [0, 0, 1] for loading along the Z-axis). This is a dimensionless direction, not a magnitude.
-2. For each **Feature**, the filter rotates the loading direction from the sample frame into the grain's crystal frame using the grain's average orientation.
-3. The filter evaluates the Schmid factor for every slip system available for that grain's crystal structure.
-4. By default, the slip system with the **largest** Schmid factor is reported.  Alternatively, the user can specify a particular slip plane and slip direction to evaluate.
+| Override Default Slip System | *Slip Systems* value |
+|---|---|
+| Off (default) | Index into the crystal structure's built-in slip system list — numbered from **0** (0–11) for Cubic-High and Cubic-Low, and from **1** (1–6) for Hexagonal-High and Hexagonal-Low |
+| On | Index of the **symmetry operator** (0–23 for Cubic-High) that maps the user-supplied slip plane and slip direction onto the winning variant |
 
-### Note
+The two numbering schemes are not comparable, and neither is the base of the *off* numbering comparable between the cubic and hexagonal classes.
 
-Only the geometric Schmid factor is considered. The critical resolved shear stress (CRSS), which varies between slip systems in real materials, is not taken into account. This means the reported "most favorable" slip system is based purely on geometric alignment, not on the actual stress required to activate slip.
+A *Slip Systems* value of `0` can mean **"no slip system found"** in two situations:
+
+- With *Override Default Slip System* **on** and no symmetry-operator variant producing a non-zero Schmid factor — for example a loading direction parallel to the slip plane normal, which puts the loading axis perpendicular to every slip direction.
+- With the toggle **off** on a **Hexagonal-High or Hexagonal-Low** phase, where the built-in systems are numbered 1–6 and no candidate exceeds the initial Schmid factor of `0`. Because that numbering starts at 1, `0` falls outside the valid range and is unambiguous. (For the cubic classes, whose numbering starts at 0, the reported `0` in this situation is indistinguishable from a genuine win on system 0; check the reported Schmid factor, which is `0` in the degenerate case.)
+
+In every one of these cases the reported Schmid factor is `0`.
+
+With *Override Default Slip System* on, the reported index is **relative to the symmetry-operator table** of the orientation library, and several operators typically tie at the maximum Schmid factor — for a `(001)[100]` system under `[1,2,3]` loading, six of the twenty-four do. The reported index is simply the first of those in table order. It identifies *a* maximizing variant, not a uniquely determined one, and the same input run through DREAM3D 6.5.x reports a different index for the same physical answer because that version's table is ordered differently. Treat this value as a handle into the current library's table, not as a portable label.
+
+### Phis and Lambdas Units Depend on Override Default Slip System
+
+The *Phis* and *Lambdas* arrays change **units** with the same toggle, which is easy to miss because the array names do not change:
+
+| Override Default Slip System | *Phis* / *Lambdas* contents |
+|---|---|
+| Off (default) | The direction **cosines** cos φ and cos λ (dimensionless, in [0, 1]) |
+| On | The **angles** φ and λ themselves, in **radians** |
+
+With the toggle off, the reported Schmid factor is the product of the two stored values. With the toggle on, it is the product of their cosines. If you need consistent units across both modes, convert explicitly rather than assuming.
+
+### Poles Is Not a Miller Index
+
+Despite its name, the *Poles* array does **not** contain a crystallographic index. It is the unit loading direction expressed in the **Feature's** crystal frame, multiplied by 100 and **truncated toward zero** to an integer — a compact fixed-point encoding of that unit vector, retaining two decimal places. Components can be negative, and the sum of squares is approximately 10 000 rather than 1. Because the conversion truncates rather than rounds, a component whose scaled value falls near an integer can differ by one from the value you would get by rounding. Reducing the triplet to a Miller index requires dividing by the greatest common divisor yourself, and the truncation means the result is only approximate.
 
 ### Required Input Sources
 
-- **Average Quaternions** -- the per-feature average orientation produced by [Compute Average Orientations](ComputeAvgOrientationsFilter.md).
+- **Average Quaternions** -- produced by [Compute Average Orientations](ComputeAvgOrientationsFilter.md).
 - **Feature Phases** -- produced by [Compute Feature Phases](../SimplnxCore/ComputeFeaturePhasesFilter.md).
-- **Crystal Structures** -- ensemble-level array read from EBSD data or created by [Create Ensemble Info](CreateEnsembleInfoFilter.md).
+- **Crystal Structures** -- read from EBSD data or produced by [Create Ensemble Info](CreateEnsembleInfoFilter.md).
 
 % Auto generated parameter table will be inserted here
 
 ## Example Pipelines
 
-+ `(05) SmallIN100 Crystallographic Statistics`
++ (04) Small IN100 Crystallographic Statistics
 
 ## License & Copyright
 

src/Plugins/OrientationAnalysis/src/OrientationAnalysis/Filters/Algorithms/ComputeSchmids.cpp

@@ -9,19 +9,19 @@
 using namespace nx::core;
 
 // -----------------------------------------------------------------------------
-ComputeSchmids::ComputeSchmids(DataStructure& dataStructure, const IFilter::MessageHandler& mesgHandler, const std::atomic_bool& shouldCancel, ComputeSchmidsInputValues* inputValues)
+ComputeSchmids::ComputeSchmids(DataStructure& dataStructure, const IFilter::MessageHandler& messageHandler, const std::atomic_bool& shouldCancel, const ComputeSchmidsInputValues* inputValues)
 : m_DataStructure(dataStructure)
 , m_InputValues(inputValues)
 , m_ShouldCancel(shouldCancel)
-, m_MessageHandler(mesgHandler)
+, m_MessageHandler(messageHandler)
 {
 }
 
 // -----------------------------------------------------------------------------
 ComputeSchmids::~ComputeSchmids() noexcept = default;
 
 // -----------------------------------------------------------------------------
-const std::atomic_bool& ComputeSchmids::getCancel()
+const std::atomic_bool& ComputeSchmids::getCancel() const
 {
   return m_ShouldCancel;
 }
@@ -31,36 +31,36 @@ Result<> ComputeSchmids::operator()()
 {
   std::vector<ebsdlib::LaueOps::Pointer> orientationOps = ebsdlib::LaueOps::GetAllOrientationOps();
 
-  const auto& avgQuatPtr = m_DataStructure.getDataRefAs<Float32Array>(m_InputValues->AvgQuatsArrayPath);
-  const auto& featurePhases = m_DataStructure.getDataRefAs<Int32Array>(m_InputValues->FeaturePhasesArrayPath);
-  const auto& crystalStructures = m_DataStructure.getDataRefAs<UInt32Array>(m_InputValues->CrystalStructuresArrayPath);
-  auto& schmidArray = m_DataStructure.getDataRefAs<Float32Array>(m_InputValues->SchmidsArrayName);
-  auto& slipSystems = m_DataStructure.getDataRefAs<Int32Array>(m_InputValues->SlipSystemsArrayName);
-  auto& poleArrays = m_DataStructure.getDataRefAs<Int32Array>(m_InputValues->PolesArrayName);
-
-  auto* phiArray = m_DataStructure.getDataAs<Float32Array>(m_InputValues->PhisArrayName);
-  auto* lambdaArray = m_DataStructure.getDataAs<Float32Array>(m_InputValues->LambdasArrayName);
+  const auto& avgQuatsRef = m_DataStructure.getDataRefAs<Float32Array>(m_InputValues->AvgQuatsArrayPath);
+  const auto& featurePhasesRef = m_DataStructure.getDataRefAs<Int32Array>(m_InputValues->FeaturePhasesArrayPath);
+  const auto& crystalStructuresRef = m_DataStructure.getDataRefAs<UInt32Array>(m_InputValues->CrystalStructuresArrayPath);
+  auto& schmidsRef = m_DataStructure.getDataRefAs<Float32Array>(m_InputValues->SchmidsArrayName);
+  auto& slipSystemsRef = m_DataStructure.getDataRefAs<Int32Array>(m_InputValues->SlipSystemsArrayName);
+  auto& polesRef = m_DataStructure.getDataRefAs<Int32Array>(m_InputValues->PolesArrayName);
+
+  auto* phisPtr = m_DataStructure.getDataAs<Float32Array>(m_InputValues->PhisArrayName);
+  auto* lambdasPtr = m_DataStructure.getDataAs<Float32Array>(m_InputValues->LambdasArrayName);
+  // Feature 0 is the conventional unassigned-Feature sentinel and is never computed; the loop below
+  // starts at 1. Preflight now creates all five arrays with a "0" fill value so every tuple is
+  // defined regardless, but these explicit writes are kept because they are what documents the
+  // sentinel row's meaning at the point of use.
   if(m_InputValues->StoreAngleComponents)
   {
-    (*phiArray)[0] = 0.0F;
-    (*lambdaArray)[0] = 0.0F;
+    (*phisPtr)[0] = 0.0F;
+    (*lambdasPtr)[0] = 0.0F;
   }
-  schmidArray[0] = 0.0F;
-  poleArrays[0] = 0;
-  poleArrays[1] = 0;
-  poleArrays[2] = 0;
-  slipSystems[0] = 0;
-
-  size_t totalFeatures = avgQuatPtr.getNumberOfTuples();
+  schmidsRef[0] = 0.0F;
+  polesRef[0] = 0;
+  polesRef[1] = 0;
+  polesRef[2] = 0;
+  slipSystemsRef[0] = 0;
 
-  int32_t slipSystem = 0;
+  const usize featureCount = avgQuatsRef.getNumberOfTuples();
+  const usize ensembleCount = crystalStructuresRef.getNumberOfTuples();
 
   Eigen::Vector3d sampleLoading = {m_InputValues->LoadingDirection[0], m_InputValues->LoadingDirection[1], m_InputValues->LoadingDirection[2]};
   sampleLoading.normalize();
 
-  double angleComps[2] = {0.0f, 0.0f};
-  double schmid = 0.0f;
-
   Eigen::Vector3d plane;
   Eigen::Vector3d direction;
   if(m_InputValues->OverrideSystem)
@@ -72,41 +72,67 @@ Result<> ComputeSchmids::operator()()
     direction.normalize();
   }
 
-  for(size_t i = 1; i < totalFeatures; i++)
+  for(usize featureIdx = 1; featureIdx < featureCount; featureIdx++)
   {
     if(m_ShouldCancel)
     {
       return {};
     }
 
-    uint32_t laueClass = crystalStructures[featurePhases[i]];
+    // The phase id indexes the Crystal Structures array directly, so it has to be validated before
+    // it is used. An id at or beyond the ensemble count, or a negative one, would read outside the
+    // array and dispatch on whatever Laue value that read produced.
+    const int32 phaseId = featurePhasesRef[featureIdx];
+    if(phaseId < 0)
+    {
+      return MakeErrorResult(-13502, fmt::format("Feature {} has a negative phase value of {}. Phase values index the Crystal Structures array '{}' and must be zero or greater.", featureIdx, phaseId,
+                                                 m_InputValues->CrystalStructuresArrayPath.toString()));
+    }
+    if(static_cast<usize>(phaseId) >= ensembleCount)
+    {
+      return MakeErrorResult(-13501, fmt::format("Feature {} has a phase value of {} but the Crystal Structures array '{}' only has {} entries (valid phase values are 0 through {}).", featureIdx,
+                                                 phaseId, m_InputValues->CrystalStructuresArrayPath.toString(), ensembleCount, ensembleCount - 1));
+    }
+
+    const uint32_t laueClass = crystalStructuresRef[phaseId];
     if(laueClass >= ebsdlib::CrystalStructure::LaueGroupEnd)
     {
       continue;
     }
-    auto om = ebsdlib::QuaternionDType(avgQuatPtr[i * 4 + 0], avgQuatPtr[i * 4 + 1], avgQuatPtr[i * 4 + 2], avgQuatPtr[i * 4 + 3]).toOrientationMatrix();
-    Eigen::Vector3d crystalLoading = om * sampleLoading;
+
+    // Re-initialized every iteration on purpose. EbsdLib Laue ops that enumerate no slip systems
+    // for their class return without writing every output, so a buffer hoisted out of the loop
+    // would carry the previous Feature's values into the current Feature's row. Fixed in EbsdLib
+    // 3.1.1 as well; keeping the locals loop-scoped means this filter is correct against
+    // any EbsdLib.
+    double schmid = 0.0;
+    double angleComponents[2] = {0.0, 0.0};
+    int32_t slipSystem = 0;
+
+    const auto orientationMatrix =
+        ebsdlib::QuaternionDType(avgQuatsRef[featureIdx * 4 + 0], avgQuatsRef[featureIdx * 4 + 1], avgQuatsRef[featureIdx * 4 + 2], avgQuatsRef[featureIdx * 4 + 3]).toOrientationMatrix();
+    Eigen::Vector3d crystalLoading = orientationMatrix * sampleLoading;
 
     if(!m_InputValues->OverrideSystem)
     {
-      orientationOps[laueClass]->getSchmidFactorAndSS(crystalLoading.data(), schmid, angleComps, slipSystem);
+      orientationOps[laueClass]->getSchmidFactorAndSS(crystalLoading.data(), schmid, angleComponents, slipSystem);
     }
     else
     {
-      orientationOps[laueClass]->getSchmidFactorAndSS(crystalLoading.data(), plane.data(), direction.data(), schmid, angleComps, slipSystem);
+      orientationOps[laueClass]->getSchmidFactorAndSS(crystalLoading.data(), plane.data(), direction.data(), schmid, angleComponents, slipSystem);
     }
 
-    schmidArray[i] = static_cast<float>(schmid);
+    schmidsRef[featureIdx] = static_cast<float>(schmid);
     if(m_InputValues->StoreAngleComponents)
     {
-      (*phiArray)[i] = angleComps[0];
-      (*lambdaArray)[i] = angleComps[1];
+      (*phisPtr)[featureIdx] = angleComponents[0];
+      (*lambdasPtr)[featureIdx] = angleComponents[1];
     }
 
-    poleArrays[3 * i] = static_cast<int32>(crystalLoading[0] * 100.0);
-    poleArrays[3 * i + 1] = static_cast<int32>(crystalLoading[1] * 100.0);
-    poleArrays[3 * i + 2] = static_cast<int32>(crystalLoading[2] * 100.0);
-    slipSystems[i] = slipSystem;
+    polesRef[3 * featureIdx] = static_cast<int32>(crystalLoading[0] * 100.0);
+    polesRef[3 * featureIdx + 1] = static_cast<int32>(crystalLoading[1] * 100.0);
+    polesRef[3 * featureIdx + 2] = static_cast<int32>(crystalLoading[2] * 100.0);
+    slipSystemsRef[featureIdx] = slipSystem;
   }
 
   return {};