test: add immutable shell Jet baseline oracle

FESA.sln

@@ -14,6 +14,8 @@ Project("{BC8A1FFA-BEE3-4634-8014-F334798102B3}") = "fesa_reference_tests", "msv
 EndProject
 Project("{BC8A1FFA-BEE3-4634-8014-F334798102B3}") = "gtest", "msvc\gtest.vcxproj", "{03194BA4-3F3B-5399-BCAF-FB4FA47300CE}"
 EndProject
+Project("{BC8A1FFA-BEE3-4634-8014-F334798102B3}") = "fesa_jet_baseline_tests", "msvc\fesa_jet_baseline_tests.vcxproj", "{5FD398DA-6058-4CE3-9440-C9E0FDD180C3}"
+EndProject
 Global
 	GlobalSection(SolutionConfigurationPlatforms) = preSolution
 		Debug|x64 = Debug|x64
@@ -44,5 +46,9 @@ Global
 		{03194BA4-3F3B-5399-BCAF-FB4FA47300CE}.Debug|x64.Build.0 = Debug|x64
 		{03194BA4-3F3B-5399-BCAF-FB4FA47300CE}.Release|x64.ActiveCfg = Release|x64
 		{03194BA4-3F3B-5399-BCAF-FB4FA47300CE}.Release|x64.Build.0 = Release|x64
+		{5FD398DA-6058-4CE3-9440-C9E0FDD180C3}.Debug|x64.ActiveCfg = Debug|x64
+		{5FD398DA-6058-4CE3-9440-C9E0FDD180C3}.Debug|x64.Build.0 = Debug|x64
+		{5FD398DA-6058-4CE3-9440-C9E0FDD180C3}.Release|x64.ActiveCfg = Release|x64
+		{5FD398DA-6058-4CE3-9440-C9E0FDD180C3}.Release|x64.Build.0 = Release|x64
 	EndGlobalSection
 EndGlobal

msvc/fesa_jet_baseline_tests.vcxproj

@@ -0,0 +1,21 @@
+<?xml version="1.0" encoding="utf-8"?>
+<Project DefaultTargets="Build" xmlns="http://schemas.microsoft.com/developer/msbuild/2003">
+  <ItemGroup Label="ProjectConfigurations"><ProjectConfiguration Include="Debug|x64"><Configuration>Debug</Configuration><Platform>x64</Platform></ProjectConfiguration><ProjectConfiguration Include="Release|x64"><Configuration>Release</Configuration><Platform>x64</Platform></ProjectConfiguration></ItemGroup>
+  <PropertyGroup Label="Globals"><ProjectGuid>{5FD398DA-6058-4CE3-9440-C9E0FDD180C3}</ProjectGuid><Keyword>Win32Proj</Keyword><WindowsTargetPlatformVersion>10.0</WindowsTargetPlatformVersion></PropertyGroup>
+  <Import Project="$(VCTargetsPath)\Microsoft.Cpp.Default.props" />
+  <PropertyGroup Label="Configuration"><ConfigurationType>Application</ConfigurationType><PlatformToolset>v145</PlatformToolset><UseDebugLibraries Condition="'$(Configuration)'=='Debug'">true</UseDebugLibraries></PropertyGroup>
+  <Import Project="$(VCTargetsPath)\Microsoft.Cpp.props" />
+  <Import Project="Fesa.Common.props" />
+  <ItemDefinitionGroup>
+    <ClCompile>
+      <AdditionalIncludeDirectories>$(FesaRoot)src\fesa;$(FESA_MKL_ROOT)\include;$(FESA_TBB_ROOT)\include;$(FESA_HDF5_ROOT)\include;%(AdditionalIncludeDirectories)</AdditionalIncludeDirectories>
+      <PreprocessorDefinitions>MKL_ILP64;__TBB_NO_IMPLICIT_LINKAGE=1;H5_BUILT_AS_DYNAMIC_LIB;%(PreprocessorDefinitions)</PreprocessorDefinitions>
+    </ClCompile>
+  </ItemDefinitionGroup>
+  <ItemGroup><ClCompile Include="..\tests\unit\elements\nonlinear_mitc4_shell_jet_baseline.cpp" /></ItemGroup>
+  <ItemGroup>
+    <ProjectReference Include="fesa_solver.vcxproj"><Project>{3AEC39D9-4641-5161-9107-265EB50CA882}</Project><LinkLibraryDependencies>true</LinkLibraryDependencies></ProjectReference>
+  </ItemGroup>
+  <Import Project="$(VCTargetsPath)\Microsoft.Cpp.targets" />
+  <Import Project="Fesa.Runtime.targets" />
+</Project>

tests/CMakeLists.txt

@@ -179,6 +179,17 @@ target_compile_definitions(
         FESA_TEST_BINARY_DIR="${PROJECT_BINARY_DIR}"
 )
 
+add_executable(fesa_jet_baseline_tests
+    unit/elements/nonlinear_mitc4_shell_jet_baseline.cpp)
+target_include_directories(fesa_jet_baseline_tests PRIVATE
+    "${PROJECT_SOURCE_DIR}/src/fesa")
+target_compile_definitions(fesa_jet_baseline_tests PRIVATE
+    MKL_ILP64 __TBB_NO_IMPLICIT_LINKAGE=1 H5_BUILT_AS_DYNAMIC_LIB)
+target_compile_options(fesa_jet_baseline_tests PRIVATE
+    $<$<CXX_COMPILER_ID:MSVC>:/fp:precise /EHsc>)
+target_link_libraries(fesa_jet_baseline_tests PRIVATE
+    fesa_solver Fesa::MKL Fesa::TBB Fesa::HDF5)
+
 if(WIN32)
     if(NOT TARGET ${FESA_MKL_THREAD_TARGET} OR NOT TARGET MKL::mkl_core)
         message(FATAL_ERROR "The oneMKL runtime files required by tests were not resolved")
@@ -226,6 +237,7 @@ if(WIN32)
     fesa_stage_test_runtime(fesa_unit_tests)
     fesa_stage_test_runtime(fesa_integration_tests)
     fesa_stage_test_runtime(fesa_reference_tests)
+    fesa_stage_test_runtime(fesa_jet_baseline_tests)
 
     unset(_fesa_mkl_dispatch_dll)
     unset(_fesa_mkl_dispatch_dlls)
@@ -258,4 +270,11 @@ add_custom_target(
         fesa_unit_tests
         fesa_integration_tests
         fesa_reference_tests
+        fesa_jet_baseline_tests
 )
+
+# Verification never creates or replaces an immutable baseline snapshot.
+add_test(NAME NonlinearMitc4Jet.IndependentOracle
+    COMMAND fesa_jet_baseline_tests --verify)
+set_tests_properties(NonlinearMitc4Jet.IndependentOracle PROPERTIES
+    LABELS "unit;jet-baseline")

tests/unit/elements/nonlinear_mitc4_shell_jet_baseline.cpp

@@ -0,0 +1,561 @@
+// Standalone private baseline fixture: compile this translation unit instead of
+// nonlinear_mitc4_shell.cpp, then link the ordinary solver archive. No
+// production source, public interface, or official CLI instrumentation is
+// needed. Build target: fesa_jet_baseline_tests (MSBuild/CMake). --verify runs
+// independent checks without writing evidence. Compare immutable evidence with
+// --compare PATH, Debug to Debug or Release to Release; never regenerate a
+// snapshot during regression. See nonlinear_mitc4_shell_jet_baseline.md.
+#include <float.h>
+#include <xmmintrin.h>
+
+#include <cfenv>
+#include <chrono>
+#include <cstdint>
+#include <cstring>
+#include <fstream>
+#include <future>
+#include <iomanip>
+#include <iostream>
+#include <sstream>
+#include <stdexcept>
+#include <vector>
+
+#include "../../../src/fesa/elements/nonlinear_mitc4_shell.cpp"
+
+namespace fesa {
+namespace {
+using Coefficients = std::array<double, 24>;
+std::fenv_t initial_environment;
+unsigned initial_mxcsr;
+
+std::string FpState() {
+  unsigned control = 0;
+  _controlfp_s(&control, 0, 0);
+  std::ostringstream out;
+  out << std::fetestexcept(FE_ALL_EXCEPT) << ':' << std::fegetround() << ':'
+      << _statusfp() << ':' << control << ':' << _mm_getcsr();
+  return out.str();
+}
+
+std::string ResetFp() {
+  std::fesetenv(&initial_environment);
+  std::feclearexcept(FE_ALL_EXCEPT);
+  _clearfp();
+  _mm_setcsr(initial_mxcsr & ~0x3fu);
+  return FpState();
+}
+
+struct Monomial {
+  double coefficient;
+  std::array<unsigned, 24> powers{};
+};
+
+// Independent oracle: differentiate monomial exponent lists, never Jet rules.
+double Differentiate(const std::vector<Monomial>& polynomial,
+                     const Coefficients& point,
+                     const std::vector<std::size_t>& derivative) {
+  long double sum = 0;
+  for (auto term : polynomial) {
+    long double value = term.coefficient;
+    for (auto index : derivative) {
+      value *= term.powers[index];
+      if (term.powers[index] == 0) {
+        value = 0;
+        break;
+      }
+      --term.powers[index];
+    }
+    for (std::size_t i = 0; i < 24; ++i)
+      for (unsigned p = 0; p < term.powers[i]; ++p) value *= point[i];
+    sum += value;
+  }
+  return static_cast<double>(sum);
+}
+
+std::uint64_t Bits(double value) {
+  std::uint64_t result;
+  std::memcpy(&result, &value, sizeof(result));
+  return result;
+}
+void Check(double actual, double expected, const char* channel) {
+  if (!std::isfinite(actual) ||
+      std::abs(actual - expected) > 2e-12 * std::max(1.0, std::abs(expected)))
+    throw std::runtime_error(std::string("independent oracle: ") + channel);
+}
+Jet Polynomial(const std::vector<Monomial>& terms,
+               const std::array<Jet, 24>& variables) {
+  Jet result = Constant(0);
+  for (const auto& term : terms) {
+    Jet value = Constant(term.coefficient);
+    for (std::size_t i = 0; i < 24; ++i)
+      for (unsigned p = 0; p < term.powers[i]; ++p)
+        value = value * variables[i];
+    result = result + value;
+  }
+  return result;
+}
+std::string Snapshot(const char* name, const Jet& jet, int exceptions,
+                     unsigned status, unsigned mxcsr) {
+  std::ostringstream out;
+  out << name << ' ' << jet.constant_derivatives << ' '
+      << jet.second.IsImplicitZero() << ' ' << IsFinite(jet) << ' '
+      << exceptions << ' ' << status << ' ' << mxcsr << std::hex;
+  auto put = [&](double x) {
+    if (std::isnan(x))
+      out << " nan";
+    else
+      out << ' ' << Bits(x);
+  };
+  put(jet.value);
+  for (auto x : jet.first) put(x);
+  for (auto x : jet.second) put(x);
+  put(jet.directional_first);
+  put(jet.directional_second);
+  for (auto x : jet.mixed_first) put(x);
+  for (auto x : jet.mixed_second) put(x);
+  return out.str() + "\n";
+}
+template <class F>
+std::string Observe(const char* name, F expression, Jet* returned = nullptr) {
+  const auto before = ResetFp();
+  Jet value = expression();
+  const int exceptions = std::fetestexcept(FE_ALL_EXCEPT);
+  const unsigned status = _statusfp();
+  const unsigned mxcsr = _mm_getcsr();
+  const auto after = FpState();
+  if (returned != nullptr) *returned = value;
+  return std::string("initial-fp=") + before + " final-fp=" + after + " " +
+         Snapshot(name, value, exceptions, status, mxcsr);
+}
+
+std::string Verify(bool corrupt_oracle) {
+  Coefficients x{}, direction{};
+  std::array<Jet, 24> variables{};
+  std::vector<Monomial> terms;
+  for (std::size_t i = 0; i < 24; ++i) {
+    x[i] = (static_cast<double>(i % 7) - 3) / 8;
+    direction[i] = (static_cast<double>(i % 5) - 2) / 16;
+    variables[i] = Variable(x[i], i, direction[i]);
+    Monomial term{(i % 2 ? -1.0 : 1.0) / 8};
+    term.powers[i] = 3;
+    term.powers[(i + 1) % 24] += 2;
+    terms.push_back(term);
+    Monomial cross{0.5};
+    cross.powers[i] = 1;
+    cross.powers[(i + 11) % 24] += 1;
+    terms.push_back(cross);
+  }
+  Jet polynomial = Polynomial(terms, variables);
+  if (corrupt_oracle) polynomial.value += 1;  // Deliberate fixture-only mutant.
+  Check(polynomial.value, Differentiate(terms, x, {}), "value");
+  double d1 = 0, d2 = 0;
+  for (std::size_t i = 0; i < 24; ++i) {
+    const double gradient = Differentiate(terms, x, {i});
+    Check(polynomial.first[i], gradient, "gradient");
+    d1 += gradient * direction[i];
+    double mixed1 = 0, mixed2 = 0;
+    for (std::size_t j = 0; j < 24; ++j) {
+      const double hessian = Differentiate(terms, x, {i, j});
+      Check(polynomial.second.ReadData()[24 * i + j], hessian, "Hessian");
+      d2 += hessian * direction[i] * direction[j];
+      mixed1 += hessian * direction[j];
+      for (std::size_t k = 0; k < 24; ++k)
+        mixed2 +=
+            Differentiate(terms, x, {i, j, k}) * direction[j] * direction[k];
+    }
+    Check(polynomial.mixed_first[i], mixed1, "mixed_first");
+    Check(polynomial.mixed_second[i], mixed2, "mixed_second");
+  }
+  Check(polynomial.directional_first, d1, "directional_first");
+  Check(polynomial.directional_second, d2, "directional_second");
+  // Affine composition has a rank-one Hessian; its third derivative contracts
+  // independently to the mixed-second channel.
+  Jet affine = Constant(2);
+  Coefficients weights{};
+  double slope = 0;
+  for (std::size_t i = 0; i < 24; ++i) {
+    weights[i] = (static_cast<double>(i % 3) - 1) / 64;
+    affine = affine + weights[i] * variables[i];
+    slope += weights[i] * direction[i];
+  }
+  auto unary = [&](const Jet& jet, double value, double a, double b, double c) {
+    Check(jet.value, value, "unary value");
+    Check(jet.directional_first, a * slope, "unary direction1");
+    Check(jet.directional_second, b * slope * slope, "unary direction2");
+    for (std::size_t i = 0; i < 24; ++i) {
+      Check(jet.first[i], a * weights[i], "unary gradient");
+      Check(jet.mixed_first[i], b * weights[i] * slope, "unary mixed1");
+      Check(jet.mixed_second[i], c * weights[i] * slope * slope,
+            "unary mixed2");
+      for (std::size_t j = 0; j < 24; ++j)
+        Check(jet.second.ReadData()[24 * i + j], b * weights[i] * weights[j],
+              "unary Hessian");
+    }
+  };
+  const double t = affine.value, r = std::sqrt(t);
+  unary(Sin(affine), std::sin(t), std::cos(t), -std::sin(t), -std::cos(t));
+  unary(Cos(affine), std::cos(t), -std::sin(t), -std::cos(t), std::sin(t));
+  unary(Reciprocal(affine), 1 / t, -1 / (t * t), 2 / (t * t * t),
+        -6 / (t * t * t * t));
+  unary(Sqrt(affine), r, 0.5 / r, -0.25 / (r * r * r),
+        0.375 / (r * r * r * r * r));
+  const double den = 1 + t * t;
+  unary(Atan2(affine, Constant(1)), std::atan(t), 1 / den, -2 * t / (den * den),
+        (6 * t * t - 2) / (den * den * den));
+  const Jet original = polynomial;
+  const Jet alias = original;
+  const auto before = Snapshot("alias", alias, 0, 0, 0);
+  const Jet transformed = alias * alias + original;
+  if (before != Snapshot("alias", original, 0, 0, 0))
+    throw std::runtime_error("alias modified");
+  std::string out;
+  std::vector<std::future<std::pair<Jet, std::string>>> futures;
+  for (int i = 0; i < 4; ++i)
+    futures.push_back(std::async(std::launch::async, [&] {
+      Jet value;
+      const auto evidence = Observe(
+          "shared-full", [&] { return original * original + original; },
+          &value);
+      return std::make_pair(std::move(value), evidence);
+    }));
+  for (std::size_t i = 0; i < futures.size(); ++i) {
+    const auto worker = futures[i].get();
+    if (Snapshot("shared", worker.first, 0, 0, 0) !=
+        Snapshot("shared", transformed, 0, 0, 0))
+      throw std::runtime_error("shared complete-channel snapshot");
+    out += "worker-" + std::to_string(i) + " " + worker.second;
+  }
+  auto add = [&](const char* name, auto f) { out += Observe(name, f); };
+  add("constant", [] { return Constant(-0.0); });
+  add("variable-zero-direction", [] { return Variable(-0.0, 23); });
+  add("variable-direction", [] { return Variable(0.25, 0, -0.5); });
+  add("polynomial", [&] { return Polynomial(terms, variables); });
+  add("polynomial-zero-direction", [&] {
+    auto zero_direction = variables;
+    for (auto& variable : zero_direction) variable.directional_first = 0;
+    return Polynomial(terms, zero_direction);
+  });
+  add("actual-zero-hessian", [&] { return original - original; });
+  add("implicit-zero-nonconstant",
+      [] { return Variable(1, 0) - Variable(1, 0); });
+  add("implicit-add", [] { return Variable(1, 0) + Constant(-1); });
+  add("implicit-subtract", [] { return Variable(1, 0) - Constant(1); });
+  for (int kind = 0; kind < 4; ++kind) {
+    const auto name = "surviving-copy-" + std::to_string(kind);
+    add(name.c_str(), [&] {
+      Jet survivor;
+      std::string expected;
+      {
+        Jet source = kind == 0   ? Constant(-0.0)
+                     : kind == 1 ? Variable(0.5, 7, -0.25)
+                                 : Polynomial(terms, variables);
+        if (kind == 3) source = source - source;
+        const Jet another_alias = source;
+        expected = Snapshot("survivor", source, 0, 0, 0);
+        survivor = another_alias;
+      }
+      if (Snapshot("survivor", survivor, 0, 0, 0) != expected)
+        throw std::runtime_error("surviving copy complete-channel snapshot");
+      return survivor;
+    });
+  }
+  add("constant-multiply", [] { return Constant(-0.0) * Constant(3); });
+  add("constant-compose", [] { return Sin(Constant(0)); });
+  add("scale", [&] { return -0.5 * original; });
+  add("multiply", [&] { return original * original; });
+  add("compose-sin", [&] { return Sin(original); });
+  add("compose-cos", [&] { return Cos(original); });
+  add("divide", [&] { return original / affine; });
+  add("sqrt", [&] { return Sqrt(affine); });
+  add("atan2", [&] { return Atan2(original, affine); });
+  add("dot", [&] {
+    return Dot(JetVector3{original, affine, variables[0]},
+               JetVector3{affine, variables[1], original});
+  });
+  add("cancellation",
+      [] { return (Variable(1e16, 0, 1) + Constant(1)) - Constant(1e16); });
+  add("overflow", [] {
+    return Variable(std::numeric_limits<double>::max(), 0, 1) * Constant(2);
+  });
+  add("subnormal", [] {
+    return std::numeric_limits<double>::denorm_min() * Variable(0.5, 0, 0.5);
+  });
+  add("nan", [] {
+    return Constant(std::numeric_limits<double>::quiet_NaN()) * Variable(1, 0);
+  });
+  add("inf", [] {
+    return Constant(std::numeric_limits<double>::infinity()) * Constant(0);
+  });
+  add("singular", [] { return Reciprocal(Variable(0, 0, 1)); });
+  add("singular-zero-direction", [] { return Reciprocal(Variable(0, 0)); });
+  add("constant-singular-reciprocal", [] { return Reciprocal(Constant(0)); });
+  add("constant-singular-sqrt", [] { return Sqrt(Constant(0)); });
+  add("negative-sqrt", [] { return Sqrt(Variable(-1, 0, 1)); });
+  add("after-failure", [&] { return Polynomial(terms, variables); });
+  // Exercise the production MKL k=1 outer path directly, including entries
+  // whose rounded bits agree while their IEEE exception state can differ.
+  auto outer = [](double a, double b) {
+    Jet jet{};
+    Coefficients left{}, right{};
+    for (std::size_t i = 0; i < 24; ++i) {
+      left[i] = i % 2 ? -a : a;
+      right[i] = i % 3 ? -b : b;
+    }
+    Outer(left, right, jet.second.OverwriteData());
+    return jet;
+  };
+  add("outer-signed-zero", [&] { return outer(-0.0, 1); });
+  add("outer-underflow",
+      [&] { return outer(std::numeric_limits<double>::denorm_min(), 0.5); });
+  add("outer-overflow",
+      [&] { return outer(std::numeric_limits<double>::max(), 2); });
+  add("outer-invalid",
+      [&] { return outer(std::numeric_limits<double>::infinity(), 0); });
+  if (!Constant(1).constant_derivatives ||
+      !Constant(1).second.IsImplicitZero() ||
+      Variable(1, 0).constant_derivatives ||
+      !Variable(1, 0).second.IsImplicitZero() ||
+      polynomial.second.IsImplicitZero())
+    throw std::runtime_error("provenance");
+  return out;
+}
+
+// Public request-mode evidence for later direction-channel separation. Existing
+// tangent tests check nonzero-direction derivatives but do not freeze the full
+// omitted-direction response inventory and invalid-direction diagnostics.
+std::string PublicSnapshots() {
+  std::array<Node, 4> nodes;
+  const std::array<std::array<double, 3>, 4> coordinates{
+      {{{-1, -1, 0}}, {{1, -1, 0}}, {{1, 1, 0}}, {{-1, 1, 0}}}};
+  std::array<const Node*, 4> pointers{};
+  std::array<ShellNodeInitialFrame, 4> frames{};
+  NonlinearMitc4TrialDofs trial{};
+  Vector zero(24), direction(24), wrong(23), nan(24), inf(24);
+  for (std::size_t i = 0; i < 4; ++i) {
+    nodes[i] = {{"Jet-Baseline", static_cast<std::int64_t>(i + 1),
+                 std::to_string(i + 1)},
+                coordinates[i],
+                {"jet-baseline.inp", i + 1}};
+    pointers[i] = &nodes[i];
+    frames[i] = {static_cast<EntityIndex>(i), {0, 0, 1}, {1, 0, 0}, {0, 1, 0}};
+    trial.displacement[i] = {0.01 * i, -0.02 * i, 0.015 * i};
+    trial.rotation_vector[i] = {0.02 * i, -0.01 * i, 0.005 * i};
+  }
+  for (std::size_t i = 0; i < 24; ++i)
+    direction[i] = (static_cast<double>(i % 5) - 2) / 8;
+  nan[7] = std::numeric_limits<double>::quiet_NaN();
+  inf[7] = std::numeric_limits<double>::infinity();
+  const auto created = NonlinearMitc4Shell::Create(
+      pointers, frames, {"Section", 0.2, 0, {"jet-baseline.inp", 20}},
+      {"Material", 120, 0.25, {"jet-baseline.inp", 21}});
+  if (!created.HasValue()) throw std::runtime_error("public fixture creation");
+  std::ostringstream out;
+  for (const auto mode : {Mitc4DrillingMode::kDirectSpring,
+                          Mitc4DrillingMode::kLegacyCorotatedTwist}) {
+    const auto shell =
+        mitc4_drilling_internal::TestAccess::WithMode(created.Value(), mode);
+    auto observe = [&](const char* name, const NonlinearMitc4TrialDofs& input,
+                       const Vector* request,
+                       std::optional<bool> expected_success) {
+      const auto before = ResetFp();
+      const auto result =
+          EvaluateInternalResidualAndTangent(shell, input, request);
+      const int exceptions = std::fetestexcept(FE_ALL_EXCEPT);
+      const unsigned status = _statusfp(), mxcsr = _mm_getcsr();
+      if (expected_success && result.HasValue() != *expected_success)
+        throw std::runtime_error(std::string("public outcome: ") + name);
+      out << "public " << static_cast<int>(mode) << ' ' << name << ' '
+          << result.HasValue() << ' ' << exceptions << ' ' << status << ' '
+          << mxcsr << " initial-fp=" << before << " final-fp=" << FpState();
+      auto put = [&](double value) {
+        if (!std::isfinite(value))
+          throw std::runtime_error("public nonfinite field");
+        out << ' ' << std::hex << Bits(value) << std::dec;
+      };
+      auto vector = [&](const Vector& value) {
+        out << " v" << value.Size();
+        for (std::size_t i = 0; i < value.Size(); ++i) put(value[i]);
+      };
+      auto matrix = [&](const Matrix& value) {
+        out << " m" << value.Rows() << 'x' << value.Columns();
+        for (std::size_t i = 0; i < value.Rows() * value.Columns(); ++i)
+          put(value.Data()[i]);
+      };
+      if (result.HasValue()) {
+        const auto& r = result.Value();
+        vector(r.internal_residual);
+        matrix(r.internal_tangent);
+        vector(r.internal_force);
+        put(r.physical_strain_energy);
+        vector(r.physical_force);
+        vector(r.drilling_force);
+        matrix(r.drilling_tangent);
+        matrix(r.drilling_metric);
+        put(r.drilling_energy);
+        vector(r.force_directional_second);
+        vector(r.drilling_force_directional_second);
+        vector(r.drilling_constraint);
+        matrix(r.drilling_constraint_gradient);
+        for (const auto& hessian : r.drilling_constraint_hessian)
+          matrix(hessian);
+        for (const auto& third :
+             r.drilling_constraint_third_directional_gradient)
+          vector(third);
+        put(r.recovery.physical_strain_energy);
+        for (const auto& point : r.recovery.gauss_points) {
+          out << " location=" << static_cast<int>(point.location);
+          for (double value : point.natural_coordinates) put(value);
+          for (const auto& axis : {point.local_frame.e1, point.local_frame.e2,
+                                   point.local_frame.e3})
+            for (double value : axis) put(value);
+          for (double value : point.generalized_strain) put(value);
+          for (double value : point.section_resultant) put(value);
+          for (const auto& stress : point.stress) {
+            out << " section=" << static_cast<int>(stress.position);
+            put(stress.zeta);
+            for (double value : stress.components) put(value);
+          }
+        }
+      } else {
+        const auto& failure = result.GetStatus();
+        out << " category="
+            << (failure.Category() ? static_cast<int>(*failure.Category())
+                                   : -1);
+        for (const auto& diagnostic : failure.Diagnostics())
+          out << ' ' << static_cast<int>(diagnostic.severity) << ' '
+              << std::quoted(diagnostic.code) << ' '
+              << std::quoted(diagnostic.location.file.generic_string()) << ' '
+              << diagnostic.location.line << ' '
+              << std::quoted(diagnostic.keyword) << ' '
+              << std::quoted(diagnostic.entity_identity) << ' '
+              << std::quoted(diagnostic.message);
+      }
+      out << '\n';
+    };
+    observe("none", trial, nullptr, true);
+    observe("zero", trial, &zero, true);
+    observe("nonzero", trial, &direction, true);
+    observe("invalid-dimension", trial, &wrong, false);
+    observe("nan-direction", trial, &nan, false);
+    observe("inf-direction", trial, &inf, false);
+    auto invalid = trial;
+    invalid.displacement[1] = {std::numeric_limits<double>::quiet_NaN(), 0, 0};
+    observe("nan-trial-none", invalid, nullptr, false);
+    observe("none-after-failure", trial, nullptr, true);
+    auto requests = [&](const std::string& name,
+                        const NonlinearMitc4TrialDofs& input) {
+      observe((name + "-none").c_str(), input, nullptr, std::nullopt);
+      observe((name + "-zero").c_str(), input, &zero, std::nullopt);
+      observe((name + "-nonzero").c_str(), input, &direction, std::nullopt);
+    };
+    requests("zero-rotation", NonlinearMitc4TrialDofs{});
+    for (const double threshold : {1e-8, 0.25}) {
+      const double angle = std::sqrt(threshold);
+      const std::array<double, 3> angles{
+          std::nextafter(angle, 0.0), angle,
+          std::nextafter(angle, std::numeric_limits<double>::infinity())};
+      for (std::size_t i = 0; i < angles.size(); ++i) {
+        NonlinearMitc4TrialDofs boundary{};
+        for (auto& rotation : boundary.rotation_vector)
+          rotation = {angles[i], 0, 0};
+        requests((threshold == 0.25 ? "rotation-squared-quarter-"
+                                    : "rotation-squared-small-") +
+                     std::to_string(i),
+                 boundary);
+      }
+    }
+    auto translated = trial;
+    for (auto& displacement : translated.displacement)
+      displacement = displacement + Vector3{1e9, -2e9, 3e9};
+    requests("large-translation", translated);
+  }
+  return out.str();
+}
+
+__declspec(noinline) double Consume(const Jet& value) {
+  double result =
+      value.value + value.directional_first + value.directional_second;
+  for (double entry : value.first) result += entry;
+  for (double entry : value.second) result += entry;
+  for (double entry : value.mixed_first) result += entry;
+  for (double entry : value.mixed_second) result += entry;
+  return result + value.constant_derivatives + value.second.IsImplicitZero();
+}
+
+void Bench() {
+  volatile double consumed = 0;
+  for (int kind = 0; kind < 8; ++kind) {
+    Jet left =
+        kind == 0 ? Constant(0.5) : Variable(0.5, 0, kind % 2 ? 0.0 : 0.25);
+    Jet right = Variable(0.75, 1, kind == 0 || kind % 2 ? 0.0 : -0.5);
+    if (kind >= 3) {
+      left = left * left;
+      right = right * right;
+    }
+    if (kind == 5) left = left - left;
+    if (kind == 7) {
+      left = Constant(0.5);
+      right = Constant(0.75);
+    }
+    const auto start = std::chrono::steady_clock::now();
+    for (int i = 0; i < 20000; ++i) {
+      left.value += (i % 2 ? -1.0 : 1.0) / 1048576;
+      const Jet value = kind == 6 ? Dot(JetVector3{left, right, left},
+                                        JetVector3{right, left, right})
+                                  : left * right;
+      consumed = consumed + Consume(value);
+    }
+    std::cout << "bench," << kind << ",20000,"
+              << std::chrono::duration<double>(
+                     std::chrono::steady_clock::now() - start)
+                     .count()
+              << ",allocation-included," << consumed << '\n';
+  }
+}
+}  // namespace
+}  // namespace fesa
+
+int main(int argc, char** argv) {
+  try {
+    std::fegetenv(&fesa::initial_environment);
+    fesa::initial_mxcsr = _mm_getcsr();
+    if (argc == 2 && std::string(argv[1]) == "--bench") {
+      fesa::Bench();
+      return 0;
+    }
+    const bool verify_only = argc == 2 && std::string(argv[1]) == "--verify";
+    if (argc < 3 && !verify_only)
+      throw std::runtime_error(
+          "--verify OR --write|--compare PATH [--mutate|--corrupt-oracle]");
+    std::string actual =
+        fesa::Verify(argc > 3 && std::string(argv[3]) == "--corrupt-oracle");
+    actual += fesa::PublicSnapshots();
+    if (argc > 3 && std::string(argv[3]) == "--mutate")
+      actual[actual.find('\n') - 1] ^= 1;
+    if (verify_only) {
+      // CTest and ordinary regression do not create baseline files.
+    } else if (std::string(argv[1]) == "--write") {
+      if (std::ifstream(argv[2], std::ios::binary))
+        throw std::runtime_error("snapshot already exists");
+      std::ofstream output(argv[2], std::ios::binary);
+      output << actual;
+      output.close();
+      if (!output) throw std::runtime_error("snapshot write failure");
+    } else if (std::string(argv[1]) == "--compare") {
+      std::ifstream input(argv[2], std::ios::binary);
+      const std::string expected((std::istreambuf_iterator<char>(input)), {});
+      if (expected != actual)
+        throw std::runtime_error("snapshot mismatch (bits/provenance/fenv)");
+    } else {
+      throw std::runtime_error("unknown verification mode");
+    }
+    std::cout
+        << argv[1]
+        << " passed; independent oracle and 42 private cases "
+           "(651 floating channels each), 64 public request/error cases\n";
+    return 0;
+  } catch (const std::exception& error) {
+    std::cerr << error.what() << '\n';
+    return 1;
+  }
+}

tests/unit/elements/nonlinear_mitc4_shell_jet_baseline.md

@@ -0,0 +1,40 @@
+# Nonlinear MITC4 Jet 기준 검사
+
+`fesa_jet_baseline_tests`는 제품의 private 구현을 별도 번역 단위에 포함하고
+일반 solver archive를 연결한다. 제품 source와 공개 API를 계측하지 않는다.
+독립 다항식 지수 미분 및 미분 항등식과 원래 경로의 비트 snapshot은 별도 근거다.
+
+VS x64 환경과 저장소의 승인된 dependency 환경변수에서 실행한다.
+
+```bat
+MSBuild msvc\fesa_jet_baseline_tests.vcxproj /m:2 /p:Configuration=Debug /p:Platform=x64
+.build\native\Debug\tests\fesa_jet_baseline_tests.exe --verify
+cmake --preset windows-msvc-ninja-debug
+cmake --build --preset windows-msvc-ninja-debug --target fesa_jet_baseline_tests --parallel 2
+ctest --preset windows-msvc-ninja-debug -R "^NonlinearMitc4Jet\.IndependentOracle$" --output-on-failure
+```
+
+Release는 Native의 `Configuration=Release` 또는 해당 CMake Release preset을 쓴다.
+기존 unit/integration/reference 세 GTest 실행물의 발견 identity와 이 전용 검사의
+선택 수는 구분한다. CTest의 `--verify`는 baseline 파일을 생성하지 않는다.
+
+승인된 기준 revision에서만 새 경로에 `--write <snapshot>`을 한 번 실행한다.
+기존 파일은 덮어쓰지 않는다. 이후에는 원래 파일과 hash를 보존하고 같은 빌드 경로와
+구성끼리 비교한다. 일반 회귀나 후보 코드에서 `--write`로 기대값을 만들지 않는다.
+
+```bat
+.build\native\Debug\tests\fesa_jet_baseline_tests.exe --compare "<immutable-Debug-snapshot>"
+.build\native\Release\tests\fesa_jet_baseline_tests.exe --compare "<immutable-Release-snapshot>"
+.build\cmake-debug\tests\fesa_jet_baseline_tests.exe --compare "<immutable-CMake-Debug-snapshot>"
+```
+
+`--compare <snapshot> --mutate`와 `--corrupt-oracle`는 각각 snapshot 비교와 독립
+oracle의 fixture 전용 오류 주입이며 exit 1이 기대값이다. `--bench`는 별도 micro
+측정이다. 연산의 할당·해제와 교대 입력값 갱신을 포함하고 fixture 생성은 제외한다.
+전체 채널을 읽는 noinline checksum과 volatile consume을 포함해 결과 제거를 막는다.
+공식 모델 CLI의 시간과 혼합하지 않는다.
+
+snapshot은 finite 값·signed zero·Inf 비트, NaN 분류, provenance와 implicit Hessian,
+복원한 초기 FP 환경 및 연산 뒤 FE/rounding/control word/full MXCSR를 기록한다.
+공개 반환의 방향 요청 유무·실패 진단·회복량 identity도 포함한다. 수식의 정확성은
+독립 oracle로 검사하고, snapshot은 기준의 계산 순서와 관측 의미 보존을 검사한다.