feat: integrate checked thermal shell operators

include/fesa/analysis/thermal_time_grid.h

@@ -0,0 +1,48 @@
+#ifndef FESA_ANALYSIS_THERMAL_TIME_GRID_H_
+#define FESA_ANALYSIS_THERMAL_TIME_GRID_H_
+
+#include <cstdint>
+
+#include "fesa/core/status.h"
+
+namespace fesa {
+
+/// @brief Carries one stored time and its actual integration coefficient.
+struct ThermalTimeFrame {
+  double time;
+  double delta;
+};
+
+/// @brief Stores an O(1) fixed-increment grid with exact input-cell endpoint
+/// tests.
+class ThermalTimeGrid final {
+ public:
+  /// @brief Checks finite positive controls and the uint32 INC bound before
+  /// use.
+  static Result<ThermalTimeGrid> Create(double period, double increment,
+                                        std::uint32_t maximum = 100U);
+  /// @brief Returns the number of completed increments, excluding initial time.
+  std::uint32_t Count() const noexcept;
+  /// @brief Returns fl(index * nominal increment), or the exact final input
+  /// time.
+  /// @pre index is in [0, Count()].
+  double Time(std::uint32_t index) const noexcept;
+  /// @brief Returns the actual coefficient, never a difference of stored times.
+  /// @pre index is in [1, Count()].
+  double Delta(std::uint32_t index) const noexcept;
+  /// @brief Checks index, finite values and strict increase for one frame in
+  /// O(1).
+  Result<ThermalTimeFrame> Frame(std::uint32_t index) const;
+
+ private:
+  /// @brief Stores controls only after count and terminal interval validation.
+  ThermalTimeGrid(double period, double increment, std::uint32_t count,
+                  double final_delta);
+  double period_;
+  double increment_;
+  std::uint32_t count_;
+  double final_delta_;
+};
+
+}  // namespace fesa
+#endif  // FESA_ANALYSIS_THERMAL_TIME_GRID_H_

include/fesa/elements/thermal_shell.h

@@ -0,0 +1,73 @@
+#ifndef FESA_ELEMENTS_THERMAL_SHELL_H_
+#define FESA_ELEMENTS_THERMAL_SHELL_H_
+
+#include <array>
+
+#include "fesa/core/source_identity.h"
+#include "fesa/core/status.h"
+#include "fesa/math/matrix.h"
+#include "fesa/math/vector3.h"
+
+namespace fesa {
+
+/// @brief Holds global heat flux at four CCW surface points times six traces.
+/// @note Trace order is -1+, -.5, 0-, 0+, .5, 1-. Frame columns are e1,e2,n.
+struct ThermalShellRecovery {
+  std::array<Vector3, 24> global_flux;
+  std::array<std::array<Vector3, 3>, 4> frames;
+};
+
+/// @brief Owns checked local Q4 times two-segment P2 transient-heat operators.
+/// @note Twenty DOFs are node-major, five stations bottom to top per node.
+class ThermalShell final {
+ public:
+  /// @brief Builds finite consistent capacity/conduction and flux operators.
+  /// @param coordinates Original midsurface Cartesian coordinates, without
+  /// projection.
+  /// @return A complete kernel or a source-identified geometry/material
+  /// failure.
+  static Result<ThermalShell> Create(SourceEntityId source,
+                                     std::array<EntityIndex, 4> node_indices,
+                                     const std::array<Vector3, 4>& coordinates,
+                                     double conductivity, double density,
+                                     double specific_heat, double thickness,
+                                     SourceLocation location = {});
+  /// @brief Returns the original element identity.
+  const SourceEntityId& SourceId() const noexcept;
+  /// @brief Returns Domain node indices, not equation numbers.
+  const std::array<EntityIndex, 4>& NodeIndices() const noexcept;
+  /// @brief Returns the node-major conductivity matrix [power/temperature].
+  const Matrix& Conductivity() const noexcept;
+  /// @brief Returns the consistent capacity matrix [energy/temperature].
+  const Matrix& Capacity() const noexcept;
+  /// @brief Returns componentwise outward roundoff bounds, excluding quadrature
+  /// error.
+  const Matrix& ConductivityError() const noexcept;
+  /// @brief Returns componentwise outward roundoff bounds, excluding quadrature
+  /// error.
+  const Matrix& CapacityError() const noexcept;
+  /// @brief Evaluates q=-k*grad(T) without averaging the two midsurface traces.
+  /// @return Finite global fluxes and right-handed frames or source-identified
+  /// failure.
+  Result<ThermalShellRecovery> Recover(
+      const std::array<double, 20>& temperature) const;
+
+ private:
+  /// @brief Initializes an unpublished candidate whose operators are filled by
+  /// Create.
+  ThermalShell(SourceEntityId source, std::array<EntityIndex, 4> node_indices,
+               SourceLocation location, double conductivity);
+  SourceEntityId source_;
+  std::array<EntityIndex, 4> node_indices_;
+  SourceLocation location_;
+  double conductivity_value_;
+  Matrix conductivity_{20U, 20U};
+  Matrix capacity_{20U, 20U};
+  Matrix conductivity_error_{20U, 20U};
+  Matrix capacity_error_{20U, 20U};
+  std::array<std::array<Vector3, 20>, 24> gradients_{};
+  std::array<std::array<Vector3, 3>, 4> frames_{};
+};
+
+}  // namespace fesa
+#endif  // FESA_ELEMENTS_THERMAL_SHELL_H_

src/fesa/analysis/thermal_time_grid.cpp

@@ -0,0 +1,232 @@
+#include "fesa/analysis/thermal_time_grid.h"
+
+#include <array>
+#include <cmath>
+#include <cstdint>
+#include <cstring>
+#include <limits>
+
+namespace fesa {
+namespace {
+
+// Fixed four base-2^32 digits. Contract bounds all exact values to 122 bits.
+struct UInt128 {
+  std::array<std::uint32_t, 4> digit{};
+};
+
+/// @brief Orders unsigned integers from their most significant digit.
+int Compare(const UInt128& a, const UInt128& b) {
+  for (int i = 3; i >= 0; --i) {
+    if (a.digit[i] != b.digit[i]) return a.digit[i] < b.digit[i] ? -1 : 1;
+  }
+  return 0;
+}
+
+/// @brief Adds with explicit carry; false means the fixed storage overflowed.
+bool Add(const UInt128& a, const UInt128& b, UInt128& result) {
+  std::uint64_t carry = 0;
+  for (std::size_t i = 0; i < 4; ++i) {
+    const std::uint64_t sum = std::uint64_t{a.digit[i]} + b.digit[i] + carry;
+    result.digit[i] = static_cast<std::uint32_t>(sum);
+    carry = sum >> 32;
+  }
+  return carry == 0;
+}
+
+/// @brief Subtracts nonnegative magnitudes with explicit borrow.
+bool Subtract(const UInt128& a, const UInt128& b, UInt128& result) {
+  std::uint64_t borrow = 0;
+  for (std::size_t i = 0; i < 4; ++i) {
+    const std::uint64_t sub = std::uint64_t{b.digit[i]} + borrow;
+    result.digit[i] =
+        static_cast<std::uint32_t>(std::uint64_t{a.digit[i]} - sub);
+    borrow = std::uint64_t{a.digit[i]} < sub ? 1U : 0U;
+  }
+  return borrow == 0;
+}
+
+/// @brief Places a significand on the exact 2^-88 scaled lattice.
+bool Shift(std::uint64_t value, int shift, UInt128& result) {
+  result = {};
+  if (shift < 0) return false;
+  for (int bit = 0; bit < 64; ++bit) {
+    if ((value & (std::uint64_t{1} << bit)) == 0) continue;
+    const int target = bit + shift;
+    if (target >= 128) return false;
+    result.digit[target / 32] |= std::uint32_t{1} << (target % 32);
+  }
+  return true;
+}
+
+/// @brief Multiplies by at most 2^32+1 using bounded schoolbook products.
+bool Multiply(const UInt128& a, std::uint64_t multiplier, UInt128& result) {
+  result = {};
+  for (std::size_t j = 0; j < 2; ++j) {
+    const auto factor = static_cast<std::uint32_t>(multiplier >> (32 * j));
+    std::uint64_t carry = 0;
+    for (std::size_t i = 0; i < 4; ++i) {
+      if (i + j >= 4) {
+        if (a.digit[i] != 0 && factor != 0) return false;
+        continue;
+      }
+      const std::uint64_t product =
+          std::uint64_t{a.digit[i]} * factor + result.digit[i + j] + carry;
+      result.digit[i + j] = static_cast<std::uint32_t>(product);
+      carry = product >> 32;
+    }
+    if (carry != 0) return false;
+  }
+  return true;
+}
+
+struct Cell {
+  UInt128 center;
+  UInt128 lower;
+  UInt128 upper;
+  bool closed;
+};
+
+/// @brief Constructs ORIGINAL binary64 cells, including subnormal/maxfinite
+/// gaps.
+bool MakeCell(double value, int period_exponent, Cell& cell) {
+  std::uint64_t bits = 0;
+  std::memcpy(&bits, &value, sizeof(bits));
+  const int exponent_bits = static_cast<int>((bits >> 52) & 0x7ffU);
+  const std::uint64_t fraction = bits & ((std::uint64_t{1} << 52) - 1);
+  const std::uint64_t significand =
+      exponent_bits == 0 ? fraction : fraction + (std::uint64_t{1} << 52);
+  const int exponent = exponent_bits == 0 ? -1074 : exponent_bits - 1075;
+  const int lattice_shift = exponent - period_exponent + 88;
+  UInt128 lower_gap;
+  UInt128 upper_gap;
+  // At a normal binade boundary only the predecessor gap is halved. The
+  // smallest normal has subnormal neighbors with the same original spacing.
+  const int boundary = fraction == 0 && exponent_bits > 1 ? 1 : 0;
+  cell.closed = (bits & 1) == 0;
+  return Shift(significand, lattice_shift, cell.center) &&
+         Shift(1, lattice_shift - 1 - boundary, lower_gap) &&
+         Shift(1, lattice_shift - 1, upper_gap) &&
+         Subtract(cell.center, lower_gap, cell.lower) &&
+         Add(cell.center, upper_gap, cell.upper);
+}
+
+/// @brief Returns source-neutral utility failure for parser/procedure
+/// decoration.
+Status GridFailure(const char* code, const char* message) {
+  return Status::Failure(
+      FailureCategory::kModel,
+      {{Severity::kError, code, {}, "HEAT TRANSFER", "", message}});
+}
+
+}  // namespace
+
+Result<ThermalTimeGrid> ThermalTimeGrid::Create(double period, double increment,
+                                                std::uint32_t maximum) {
+  if (!std::isfinite(period) || !std::isfinite(increment) || period <= 0 ||
+      increment <= 0 || maximum == 0) {
+    return Result<ThermalTimeGrid>::Failure(
+        GridFailure("invalid-thermal-time-grid",
+                    "Time controls must be finite and positive."));
+  }
+  if (increment >= period) {
+    return Result<ThermalTimeGrid>::Success(
+        ThermalTimeGrid{period, increment, 1U, period});
+  }
+  const double quotient = period / increment;
+  const std::uint64_t limit = std::uint64_t{maximum} + 1U;
+  if (!std::isfinite(quotient) || quotient > static_cast<double>(limit)) {
+    return Result<ThermalTimeGrid>::Failure(GridFailure(
+        "thermal-increment-limit", "Thermal increment count exceeds INC."));
+  }
+  // No floating-to-integer cast precedes the finite, positive count bound.
+  std::uint64_t floor = static_cast<std::uint64_t>(std::floor(quotient));
+  const double fraction = quotient - static_cast<double>(floor);
+  const std::uint64_t nearest =
+      floor + ((fraction > .5 || (fraction == .5 && (floor & 1) != 0)) ? 1 : 0);
+  int exponent = 0;
+  std::frexp(period, &exponent);
+  Cell period_cell{};
+  Cell increment_cell{};
+  UInt128 lower;
+  UInt128 upper;
+  bool exact_ok = MakeCell(period, exponent, period_cell) &&
+                  MakeCell(increment, exponent, increment_cell) &&
+                  Multiply(increment_cell.lower, nearest, lower) &&
+                  Multiply(increment_cell.upper, nearest, upper);
+  const int first = Compare(period_cell.upper, lower);
+  const int second = Compare(upper, period_cell.lower);
+  const bool closed = period_cell.closed && increment_cell.closed;
+  const bool close = nearest != 0 && (first > 0 || (first == 0 && closed)) &&
+                     (second > 0 || (second == 0 && closed));
+  std::uint64_t count = nearest;
+  if (!close && exact_ok) {
+    UInt128 product;
+    UInt128 next;
+    exact_ok = Multiply(increment_cell.center, floor, product) &&
+               Multiply(increment_cell.center, floor + 1, next);
+    if (Compare(product, period_cell.center) > 0) {
+      --floor;
+    } else if (Compare(next, period_cell.center) <= 0) {
+      ++floor;
+    }
+    exact_ok = exact_ok && Multiply(increment_cell.center, floor, product) &&
+               Multiply(increment_cell.center, floor + 1, next) &&
+               Compare(product, period_cell.center) <= 0 &&
+               Compare(next, period_cell.center) > 0;
+    count = floor + (Compare(product, period_cell.center) == 0 ? 0U : 1U);
+  }
+  if (!exact_ok) {
+    return Result<ThermalTimeGrid>::Failure(GridFailure(
+        "invalid-thermal-time-grid", "Exact time lattice bounds failed."));
+  }
+  if (count == 0 || count > maximum) {
+    return Result<ThermalTimeGrid>::Failure(GridFailure(
+        "thermal-increment-limit", "Thermal increment count exceeds INC."));
+  }
+  const double last =
+      close ? increment
+            : std::fma(-static_cast<double>(count - 1), increment, period);
+  ThermalTimeGrid grid{period, increment, static_cast<std::uint32_t>(count),
+                       last};
+  const auto frame = grid.Frame(grid.Count());
+  if (!frame.HasValue())
+    return Result<ThermalTimeGrid>::Failure(frame.GetStatus());
+  // For j <= 2^32, adjacent exact j*d are separated by more than 2^19 ULPs
+  // (subnormal products are exact). Thus the terminal check also suffices for
+  // O(1) preflight; every requested frame still checks its own stored interval.
+  return Result<ThermalTimeGrid>::Success(grid);
+}
+
+ThermalTimeGrid::ThermalTimeGrid(double period, double increment,
+                                 std::uint32_t count, double final_delta)
+    : period_{period},
+      increment_{increment},
+      count_{count},
+      final_delta_{final_delta} {}
+
+std::uint32_t ThermalTimeGrid::Count() const noexcept { return count_; }
+double ThermalTimeGrid::Time(std::uint32_t index) const noexcept {
+  return index == count_ ? period_ : static_cast<double>(index) * increment_;
+}
+double ThermalTimeGrid::Delta(std::uint32_t index) const noexcept {
+  return index == count_ ? final_delta_ : increment_;
+}
+Result<ThermalTimeFrame> ThermalTimeGrid::Frame(std::uint32_t index) const {
+  if (index == 0 || index > count_) {
+    return Result<ThermalTimeFrame>::Failure(
+        GridFailure("invalid-thermal-time-frame",
+                    "Completed frame index is outside the grid."));
+  }
+  const double time = Time(index);
+  const double delta = Delta(index);
+  const double previous = Time(index - 1);
+  if (!std::isfinite(time) || !std::isfinite(previous) ||
+      !std::isfinite(delta) || !(time > previous) || !(delta > 0)) {
+    return Result<ThermalTimeFrame>::Failure(
+        GridFailure("invalid-thermal-time-frame",
+                    "Stored times must be finite and strictly increasing."));
+  }
+  return Result<ThermalTimeFrame>::Success({time, delta});
+}
+
+}  // namespace fesa

src/fesa/elements/thermal_interval.h

@@ -0,0 +1,108 @@
+#ifndef FESA_ELEMENTS_THERMAL_INTERVAL_H_
+#define FESA_ELEMENTS_THERMAL_INTERVAL_H_
+
+#include <algorithm>
+#include <array>
+#include <cmath>
+#include <limits>
+
+namespace fesa::thermal_internal {
+
+/// @brief Encloses binary64 operations with one outward nextafter per endpoint.
+struct Interval {
+  double lower = 0;
+  double upper = 0;
+  Interval() = default;
+  explicit Interval(double value) : lower{value}, upper{value} {}
+  Interval(double low, double high) : lower{low}, upper{high} {}
+  bool IsFinite() const {
+    return std::isfinite(lower) && std::isfinite(upper) && lower <= upper;
+  }
+  bool ContainsZero() const { return lower <= 0 && upper >= 0; }
+};
+
+/// @brief Rounds an operation endpoint toward negative infinity.
+inline double Down(double value) {
+  return std::nextafter(value, -std::numeric_limits<double>::infinity());
+}
+/// @brief Rounds an operation endpoint toward positive infinity.
+inline double Up(double value) {
+  return std::nextafter(value, std::numeric_limits<double>::infinity());
+}
+inline Interval operator+(Interval a, Interval b) {
+  return {Down(a.lower + b.lower), Up(a.upper + b.upper)};
+}
+inline Interval operator-(Interval a, Interval b) {
+  return {Down(a.lower - b.upper), Up(a.upper - b.lower)};
+}
+inline Interval operator*(Interval a, Interval b) {
+  const std::array<double, 4> products{{a.lower * b.lower, a.lower * b.upper,
+                                        a.upper * b.lower, a.upper * b.upper}};
+  return {Down(*std::min_element(products.begin(), products.end())),
+          Up(*std::max_element(products.begin(), products.end()))};
+}
+inline Interval operator/(Interval a, Interval b) {
+  if (b.ContainsZero())
+    return {-std::numeric_limits<double>::infinity(),
+            std::numeric_limits<double>::infinity()};
+  const std::array<double, 4> quotients{{a.lower / b.lower, a.lower / b.upper,
+                                         a.upper / b.lower, a.upper / b.upper}};
+  return {Down(*std::min_element(quotients.begin(), quotients.end())),
+          Up(*std::max_element(quotients.begin(), quotients.end()))};
+}
+/// @brief Encloses the square without introducing negative cancellation bounds.
+inline Interval Square(Interval value) {
+  const double a = value.lower * value.lower;
+  const double b = value.upper * value.upper;
+  return {value.ContainsZero() ? 0 : Down(std::min(a, b)), Up(std::max(a, b))};
+}
+/// @brief Encloses a real square root; negative upper endpoints remain invalid.
+inline Interval Sqrt(Interval value) {
+  return {value.lower <= 0 ? 0 : Down(std::sqrt(value.lower)),
+          Up(std::sqrt(value.upper))};
+}
+/// @brief Scales exactly by powers of two, enclosing subnormal loss if present.
+inline Interval Scale(Interval value, int exponent) {
+  const double low = std::ldexp(value.lower, exponent);
+  const double high = std::ldexp(value.upper, exponent);
+  const bool low_exact =
+      std::isfinite(low) && std::ldexp(low, -exponent) == value.lower;
+  const bool high_exact =
+      std::isfinite(high) && std::ldexp(high, -exponent) == value.upper;
+  return {low_exact ? low : Down(low), high_exact ? high : Up(high)};
+}
+using IntervalVector = std::array<Interval, 3>;
+/// @brief Performs componentwise vector addition in fixed Cartesian order.
+inline IntervalVector operator+(const IntervalVector& a,
+                                const IntervalVector& b) {
+  return {{a[0] + b[0], a[1] + b[1], a[2] + b[2]}};
+}
+/// @brief Performs componentwise vector subtraction in fixed Cartesian order.
+inline IntervalVector operator-(const IntervalVector& a,
+                                const IntervalVector& b) {
+  return {{a[0] - b[0], a[1] - b[1], a[2] - b[2]}};
+}
+/// @brief Scales a vector by an interval scalar.
+inline IntervalVector operator*(Interval s, const IntervalVector& a) {
+  return {{s * a[0], s * a[1], s * a[2]}};
+}
+/// @brief Accumulates a three-component inner product in fixed order.
+inline Interval Dot(const IntervalVector& a, const IntervalVector& b) {
+  return (a[0] * b[0] + a[1] * b[1]) + a[2] * b[2];
+}
+/// @brief Encloses a right-handed Cartesian cross product.
+inline IntervalVector Cross(const IntervalVector& a, const IntervalVector& b) {
+  return {{a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2],
+           a[0] * b[1] - a[1] * b[0]}};
+}
+/// @brief Encloses a squared norm as a sum of nonnegative squares.
+inline Interval NormSquared(const IntervalVector& a) {
+  return (Square(a[0]) + Square(a[1])) + Square(a[2]);
+}
+/// @brief Checks all Cartesian components for finite ordered endpoints.
+inline bool IsFinite(const IntervalVector& a) {
+  return a[0].IsFinite() && a[1].IsFinite() && a[2].IsFinite();
+}
+
+}  // namespace fesa::thermal_internal
+#endif  // FESA_ELEMENTS_THERMAL_INTERVAL_H_

src/fesa/elements/thermal_shell.cpp

@@ -0,0 +1,236 @@
+#include "fesa/elements/thermal_shell.h"
+
+#include <algorithm>
+#include <array>
+#include <cmath>
+#include <utility>
+
+#include "thermal_shell_geometry.h"
+
+namespace fesa {
+namespace {
+using thermal_internal::Interval;
+
+/// @brief Measures the full endpoint deviation from a computed coefficient.
+double ErrorBound(Interval interval, double value) {
+  return thermal_internal::Up(std::max(std::abs(interval.lower - value),
+                                       std::abs(interval.upper - value)));
+}
+
+}  // namespace
+
+ThermalShell::ThermalShell(SourceEntityId source,
+                           std::array<EntityIndex, 4> node_indices,
+                           SourceLocation location, double conductivity)
+    : source_{std::move(source)},
+      node_indices_{node_indices},
+      location_{std::move(location)},
+      conductivity_value_{conductivity} {}
+
+Result<ThermalShell> ThermalShell::Create(
+    SourceEntityId source, std::array<EntityIndex, 4> node_indices,
+    const std::array<Vector3, 4>& coordinates, double conductivity,
+    double density, double specific_heat, double thickness,
+    SourceLocation location) {
+  using namespace thermal_internal;
+  for (const double property :
+       {conductivity, density, specific_heat, thickness}) {
+    if (!(property > 0) || !std::isfinite(property)) {
+      return Result<ThermalShell>::Failure(
+          GeometryFailure(source, location, "invalid-thermal-shell-material",
+                          "Thermal properties must be finite and positive."));
+    }
+  }
+  auto sorted = node_indices;
+  std::sort(sorted.begin(), sorted.end());
+  if (std::adjacent_find(sorted.begin(), sorted.end()) != sorted.end()) {
+    return Result<ThermalShell>::Failure(
+        GeometryFailure(source, location, "invalid-thermal-shell-geometry",
+                        "Patch contains repeated source nodes."));
+  }
+  const auto geometry = BuildGeometry(coordinates, source, location);
+  if (!geometry.HasValue())
+    return Result<ThermalShell>::Failure(geometry.GetStatus());
+  ThermalShell kernel{source, node_indices, location, conductivity};
+  Matrix mass{4, 4};
+  Matrix surface{4, 4};
+  std::array<std::array<Interval, 4>, 4> mass_interval{};
+  std::array<std::array<Interval, 4>, 4> surface_interval{};
+  const double gauss = std::sqrt(3. / 5.);
+  const Interval gauss_interval = Sqrt(Interval{3} / Interval{5});
+  const std::array<double, 3> points{{-gauss, 0., gauss}};
+  const std::array<Interval, 3> point_intervals{
+      {Interval{0} - gauss_interval, Interval{0}, gauss_interval}};
+  const std::array<double, 3> weights{{5. / 9., 8. / 9., 5. / 9.}};
+  const std::array<Interval, 3> weight_intervals{{Interval{5} / Interval{9},
+                                                  Interval{8} / Interval{9},
+                                                  Interval{5} / Interval{9}}};
+  for (std::size_t x = 0; x < 3; ++x) {
+    for (std::size_t y = 0; y < 3; ++y) {
+      const auto result = EvaluateSurface(geometry.Value(), points[x],
+                                          points[y], point_intervals[x],
+                                          point_intervals[y], source, location);
+      if (!result.HasValue())
+        return Result<ThermalShell>::Failure(result.GetStatus());
+      const auto& point = result.Value();
+      const double weight = weights[x] * weights[y] * point.jacobian;
+      const Interval weight_interval =
+          weight_intervals[x] * weight_intervals[y] * point.jacobian_interval;
+      for (std::size_t i = 0; i < 4; ++i) {
+        for (std::size_t j = 0; j < 4; ++j) {
+          mass(i, j) += point.shape[i] * point.shape[j] * weight;
+          surface(i, j) += point.gradient[i].Dot(point.gradient[j]) * weight;
+          mass_interval[i][j] = mass_interval[i][j] +
+                                point.shape_interval[i] *
+                                    point.shape_interval[j] * weight_interval;
+          surface_interval[i][j] =
+              surface_interval[i][j] +
+              Dot(point.gradient_interval[i], point.gradient_interval[j]) *
+                  weight_interval;
+        }
+      }
+    }
+  }
+  Matrix z0{5, 5};
+  Matrix z1{5, 5};
+  std::array<std::array<Interval, 5>, 5> z0_interval{};
+  std::array<std::array<Interval, 5>, 5> z1_interval{};
+  constexpr double kZ0[3][3] = {{4, 2, -1}, {2, 16, 2}, {-1, 2, 4}};
+  constexpr double kZ1[3][3] = {{7, -8, 1}, {-8, 16, -8}, {1, -8, 7}};
+  for (const std::size_t first : {0U, 2U}) {
+    for (std::size_t a = 0; a < 3; ++a) {
+      for (std::size_t b = 0; b < 3; ++b) {
+        z0(first + a, first + b) += kZ0[a][b] / 30.;
+        z1(first + a, first + b) += kZ1[a][b] / 3.;
+        z0_interval[first + a][first + b] = z0_interval[first + a][first + b] +
+                                            Interval{kZ0[a][b]} / Interval{30};
+        z1_interval[first + a][first + b] = z1_interval[first + a][first + b] +
+                                            Interval{kZ1[a][b]} / Interval{3};
+      }
+    }
+  }
+  for (std::size_t i = 0; i < 4; ++i) {
+    for (std::size_t j = 0; j < 4; ++j) {
+      for (std::size_t a = 0; a < 5; ++a) {
+        for (std::size_t b = 0; b < 5; ++b) {
+          const std::size_t p = 5 * i + a;
+          const std::size_t q = 5 * j + b;
+          kernel.capacity_(p, q) =
+              density * specific_heat * (thickness / 2) * mass(i, j) * z0(a, b);
+          kernel.conductivity_(p, q) =
+              conductivity * ((thickness / 2) * surface(i, j) * z0(a, b) +
+                              (2 / thickness) * mass(i, j) * z1(a, b));
+          const Interval capacity = Interval{density} *
+                                    Interval{specific_heat} *
+                                    (Interval{thickness} / Interval{2}) *
+                                    mass_interval[i][j] * z0_interval[a][b];
+          const Interval conduction =
+              Interval{conductivity} *
+              ((Interval{thickness} / Interval{2}) * surface_interval[i][j] *
+                   z0_interval[a][b] +
+               (Interval{2} / Interval{thickness}) * mass_interval[i][j] *
+                   z1_interval[a][b]);
+          kernel.capacity_error_(p, q) =
+              ErrorBound(capacity, kernel.capacity_(p, q));
+          kernel.conductivity_error_(p, q) =
+              ErrorBound(conduction, kernel.conductivity_(p, q));
+          if (!capacity.IsFinite() || !conduction.IsFinite() ||
+              !std::isfinite(kernel.capacity_(p, q)) ||
+              !std::isfinite(kernel.conductivity_(p, q)) ||
+              !std::isfinite(kernel.capacity_error_(p, q)) ||
+              !std::isfinite(kernel.conductivity_error_(p, q)) ||
+              (p == q && (!(kernel.capacity_(p, q) > 0) ||
+                          !(kernel.conductivity_(p, q) > 0)))) {
+            return Result<ThermalShell>::Failure(GeometryFailure(
+                source, location, "nonfinite-thermal-shell-operator",
+                "Local matrix/enclosure is not finite or a diagonal lost "
+                "positivity."));
+          }
+        }
+      }
+    }
+  }
+  const double recovery = 1 / std::sqrt(3.);
+  const Interval recovery_interval = Interval{1} / Sqrt(Interval{3});
+  constexpr std::array<double, 4> kXi{{-1, 1, 1, -1}};
+  constexpr std::array<double, 4> kEta{{-1, -1, 1, 1}};
+  constexpr std::array<double, 6> kLocal{{-1, 0, 1, -1, 0, 1}};
+  for (std::size_t ip = 0; ip < 4; ++ip) {
+    const auto result = EvaluateSurface(
+        geometry.Value(), kXi[ip] * recovery, kEta[ip] * recovery,
+        Interval{kXi[ip]} * recovery_interval,
+        Interval{kEta[ip]} * recovery_interval, source, location);
+    if (!result.HasValue())
+      return Result<ThermalShell>::Failure(result.GetStatus());
+    const auto& point = result.Value();
+    kernel.frames_[ip] = point.axes;
+    for (std::size_t trace = 0; trace < 6; ++trace) {
+      const double s = kLocal[trace];
+      const std::array<double, 3> shape{
+          {.5 * s * (s - 1), 1 - s * s, .5 * s * (s + 1)}};
+      const std::array<double, 3> derivative{{2 * s - 1, -4 * s, 2 * s + 1}};
+      const std::size_t first = trace < 3 ? 0 : 2;
+      for (std::size_t i = 0; i < 4; ++i) {
+        for (std::size_t a = 0; a < 3; ++a) {
+          const Vector3 gradient =
+              shape[a] * point.gradient[i] +
+              (point.shape[i] * (2 / thickness) * derivative[a]) *
+                  point.axes[2];
+          if (!gradient.IsFinite()) {
+            return Result<ThermalShell>::Failure(GeometryFailure(
+                source, location, "nonfinite-thermal-shell-recovery",
+                "Flux gradient is not finite."));
+          }
+          kernel.gradients_[6 * ip + trace][5 * i + first + a] = gradient;
+        }
+      }
+    }
+  }
+  return Result<ThermalShell>::Success(std::move(kernel));
+}
+
+const SourceEntityId& ThermalShell::SourceId() const noexcept {
+  return source_;
+}
+const std::array<EntityIndex, 4>& ThermalShell::NodeIndices() const noexcept {
+  return node_indices_;
+}
+const Matrix& ThermalShell::Conductivity() const noexcept {
+  return conductivity_;
+}
+const Matrix& ThermalShell::Capacity() const noexcept { return capacity_; }
+const Matrix& ThermalShell::ConductivityError() const noexcept {
+  return conductivity_error_;
+}
+const Matrix& ThermalShell::CapacityError() const noexcept {
+  return capacity_error_;
+}
+
+Result<ThermalShellRecovery> ThermalShell::Recover(
+    const std::array<double, 20>& temperature) const {
+  for (const double value : temperature) {
+    if (!std::isfinite(value)) {
+      return Result<ThermalShellRecovery>::Failure(
+          thermal_internal::GeometryFailure(
+              source_, location_, "nonfinite-thermal-shell-temperature",
+              "Recovery temperature is not finite."));
+    }
+  }
+  ThermalShellRecovery recovery{};
+  recovery.frames = frames_;
+  for (std::size_t row = 0; row < 24; ++row) {
+    Vector3 gradient;
+    for (std::size_t p = 0; p < 20; ++p)
+      gradient = gradient + temperature[p] * gradients_[row][p];
+    recovery.global_flux[row] = -conductivity_value_ * gradient;
+    if (!recovery.global_flux[row].IsFinite()) {
+      return Result<ThermalShellRecovery>::Failure(
+          thermal_internal::GeometryFailure(
+              source_, location_, "nonfinite-thermal-shell-flux",
+              "Recovered heat flux is not finite."));
+    }
+  }
+  return Result<ThermalShellRecovery>::Success(recovery);
+}
+
+}  // namespace fesa

src/fesa/elements/thermal_shell_geometry.cpp

@@ -0,0 +1,431 @@
+#include "thermal_shell_geometry.h"
+
+#include <algorithm>
+#include <array>
+#include <cmath>
+#include <cstdint>
+#include <limits>
+#include <map>
+#include <set>
+#include <string>
+#include <utility>
+#include <vector>
+
+#include "fesa/model/model_types.h"
+
+namespace fesa {
+namespace thermal_internal {
+namespace {
+
+/// @brief Represents an exact input vector by three singleton intervals.
+IntervalVector Singleton(const Vector3& point) {
+  return {{Interval{point.X()}, Interval{point.Y()}, Interval{point.Z()}}};
+}
+
+/// @brief Checks a candidate direction against the ORIGINAL interval geometry.
+bool Certify(const Geometry& geometry, Vector3 direction, double& lower) {
+  const double magnitude = direction.Norm();
+  if (!(magnitude > 0) || !std::isfinite(magnitude)) return false;
+  direction = direction / magnitude;
+  const auto s = Singleton(direction);
+  lower = std::numeric_limits<double>::infinity();
+  for (double xi : {-1., 1.}) {
+    for (double eta : {-1., 1.}) {
+      const auto w = Cross(geometry.bi + Interval{eta} * geometry.di,
+                           geometry.ci + Interval{xi} * geometry.di);
+      const Interval projection = Dot(s, w);
+      if (!projection.IsFinite() || !(projection.lower > 0)) return false;
+      lower = std::min(lower, projection.lower);
+    }
+  }
+  const Interval bound = Interval{lower} / Sqrt(NormSquared(s));
+  lower = bound.lower;
+  return bound.IsFinite() && lower > 0;
+}
+
+}  // namespace
+
+Status GeometryFailure(const SourceEntityId& source,
+                       const SourceLocation& location, const char* code,
+                       const char* message) {
+  return Status::Failure(
+      FailureCategory::kModel,
+      {{Severity::kError, code, location, "DS4",
+        source.instance_name + ":" + source.source_label_text, message}});
+}
+
+Result<Geometry> BuildGeometry(const std::array<Vector3, 4>& coordinates,
+                               const SourceEntityId& source,
+                               const SourceLocation& location) {
+  double largest = 0;
+  for (std::size_t i = 0; i < 4; ++i) {
+    if (!coordinates[i].IsFinite()) {
+      return Result<Geometry>::Failure(
+          GeometryFailure(source, location, "invalid-thermal-shell-geometry",
+                          "Coordinates must be finite."));
+    }
+    for (std::size_t j = 0; j < i; ++j) {
+      if (coordinates[i] == coordinates[j]) {
+        return Result<Geometry>::Failure(
+            GeometryFailure(source, location, "invalid-thermal-shell-geometry",
+                            "Patch has repeated coordinates."));
+      }
+    }
+    for (double component : coordinates[i].Components()) {
+      largest = std::max(largest, std::abs(component));
+    }
+  }
+  int coordinate_exponent = 0;
+  std::frexp(largest, &coordinate_exponent);
+  std::array<Vector3, 4> offsets{};
+  std::array<IntervalVector, 4> intervals{};
+  double length = 0;
+  // Scaling before subtraction prevents overflow even for opposite huge inputs.
+  for (std::size_t i = 1; i < 4; ++i) {
+    std::array<double, 3> offset{};
+    for (std::size_t j = 0; j < 3; ++j) {
+      offset[j] = std::ldexp(coordinates[i][j], -coordinate_exponent) -
+                  std::ldexp(coordinates[0][j], -coordinate_exponent);
+      intervals[i][j] =
+          Scale(Interval{coordinates[i][j]}, -coordinate_exponent) -
+          Scale(Interval{coordinates[0][j]}, -coordinate_exponent);
+      length = std::max(length, std::abs(offset[j]));
+    }
+    offsets[i] = Vector3{offset};
+  }
+  int offset_exponent = 0;
+  std::frexp(length, &offset_exponent);
+  for (std::size_t i = 1; i < 4; ++i) {
+    std::array<double, 3> offset{};
+    for (std::size_t j = 0; j < 3; ++j) {
+      offset[j] = std::ldexp(offsets[i][j], -offset_exponent);
+      intervals[i][j] = Scale(intervals[i][j], -offset_exponent);
+    }
+    offsets[i] = Vector3{offset};
+  }
+  Geometry geometry{};
+  geometry.length_exponent = coordinate_exponent + offset_exponent;
+  geometry.b = .25 * (offsets[1] + offsets[2] - offsets[3]);
+  geometry.c = .25 * (offsets[2] + offsets[3] - offsets[1]);
+  geometry.d = .25 * (offsets[2] - offsets[1] - offsets[3]);
+  geometry.bi = Interval{.25} * (intervals[1] + intervals[2] - intervals[3]);
+  geometry.ci = Interval{.25} * (intervals[2] + intervals[3] - intervals[1]);
+  geometry.di = Interval{.25} * (intervals[2] - intervals[1] - intervals[3]);
+  const Vector3 w0 = geometry.b.Cross(geometry.c);
+  const Vector3 u = geometry.b.Cross(geometry.d);
+  const Vector3 v = geometry.d.Cross(geometry.c);
+  std::vector<Vector3> candidates;
+  candidates.reserve(10);
+  for (double xi : {-1., 1.}) {
+    for (double eta : {-1., 1.}) candidates.push_back(w0 + xi * u + eta * v);
+  }
+  // Projection here is the exact constrained-edge minimizer, not a geometry
+  // clamp.
+  for (double edge : {-1., 1.}) {
+    if (v.Dot(v) > 0) {
+      const double eta = std::clamp(-v.Dot(w0 + edge * u) / v.Dot(v), -1., 1.);
+      candidates.push_back(w0 + edge * u + eta * v);
+    }
+    if (u.Dot(u) > 0) {
+      const double xi = std::clamp(-u.Dot(w0 + edge * v) / u.Dot(u), -1., 1.);
+      candidates.push_back(w0 + xi * u + edge * v);
+    }
+  }
+  const double uu = u.Dot(u);
+  const double vv = v.Dot(v);
+  const double uv = u.Dot(v);
+  const double determinant = uu * vv - uv * uv;
+  if (determinant > 0 && std::isfinite(determinant)) {
+    const double xi = (uv * v.Dot(w0) - vv * u.Dot(w0)) / determinant;
+    const double eta = (uv * u.Dot(w0) - uu * v.Dot(w0)) / determinant;
+    if (std::abs(xi) <= 1 && std::abs(eta) <= 1) {
+      candidates.push_back(w0 + xi * u + eta * v);
+    }
+  }
+  // The center is also useful for affine or rank-one normal images.
+  candidates.push_back(w0);
+  for (const auto& candidate : candidates) {
+    if (Certify(geometry, candidate, geometry.scaled_jacobian_lower)) {
+      return Result<Geometry>::Success(geometry);
+    }
+  }
+  return Result<Geometry>::Failure(
+      GeometryFailure(source, location, "uncertain-thermal-shell-geometry",
+                      "No strictly positive whole-patch interval bound was "
+                      "certified (lower <= 0)."));
+}
+
+Result<SurfacePoint> EvaluateSurface(const Geometry& geometry, double xi,
+                                     double eta, Interval xi_interval,
+                                     Interval eta_interval,
+                                     const SourceEntityId& source,
+                                     const SourceLocation& location) {
+  const Vector3 a1 = geometry.b + eta * geometry.d;
+  const Vector3 a2 = geometry.c + xi * geometry.d;
+  const auto a1i = geometry.bi + eta_interval * geometry.di;
+  const auto a2i = geometry.ci + xi_interval * geometry.di;
+  const Vector3 w = a1.Cross(a2);
+  const auto wi = Cross(a1i, a2i);
+  const double jacobian = w.Norm();
+  const Interval determinant_interval = NormSquared(wi);
+  const double determinant = w.Dot(w);
+  // |a1 x a2|^2 is the Gram determinant identity, avoiding subtraction of
+  // nearly equal g11*g22 and g12*g12. This is still the actual-point metric.
+  const Interval inverse11 = Dot(a2i, a2i) / determinant_interval;
+  const Interval inverse12 =
+      (Interval{0} - Dot(a1i, a2i)) / determinant_interval;
+  const Interval inverse22 = Dot(a1i, a1i) / determinant_interval;
+  if (!(determinant_interval.lower > 0) || !inverse11.IsFinite() ||
+      !inverse12.IsFinite() || !inverse22.IsFinite() || !(jacobian > 0) ||
+      !std::isfinite(determinant) || determinant <= 0 ||
+      !Scale(inverse11, -2 * geometry.length_exponent).IsFinite() ||
+      !Scale(inverse12, -2 * geometry.length_exponent).IsFinite() ||
+      !Scale(inverse22, -2 * geometry.length_exponent).IsFinite() ||
+      !Scale(Dot(a1i, a1i), 2 * geometry.length_exponent).IsFinite() ||
+      !Scale(Dot(a2i, a2i), 2 * geometry.length_exponent).IsFinite()) {
+    return Result<SurfacePoint>::Failure(GeometryFailure(
+        source, location, "uncertain-thermal-shell-metric",
+        "Metric inverse enclosure is not finite or excludes no zero divisor."));
+  }
+  const Vector3 reciprocal1 = a2.Cross(w) / determinant;
+  const Vector3 reciprocal2 = w.Cross(a1) / determinant;
+  const auto reciprocal1i =
+      (Interval{1} / determinant_interval) * Cross(a2i, wi);
+  const auto reciprocal2i =
+      (Interval{1} / determinant_interval) * Cross(wi, a1i);
+  SurfacePoint point{};
+  point.jacobian = std::ldexp(jacobian, 2 * geometry.length_exponent);
+  point.jacobian_interval =
+      Scale(Sqrt(determinant_interval), 2 * geometry.length_exponent);
+  point.axes[0] = a1 / a1.Norm();
+  point.axes[2] = w / jacobian;
+  point.axes[1] = point.axes[2].Cross(point.axes[0]);
+  constexpr std::array<double, 4> kXi{{-1., 1., 1., -1.}};
+  constexpr std::array<double, 4> kEta{{-1., -1., 1., 1.}};
+  for (std::size_t i = 0; i < 4; ++i) {
+    point.shape[i] = .25 * (1 + kXi[i] * xi) * (1 + kEta[i] * eta);
+    point.shape_interval[i] = Interval{.25} *
+                              (Interval{1} + Interval{kXi[i]} * xi_interval) *
+                              (Interval{1} + Interval{kEta[i]} * eta_interval);
+    const double dx = .25 * kXi[i] * (1 + kEta[i] * eta);
+    const double dy = .25 * kEta[i] * (1 + kXi[i] * xi);
+    const Interval dxi = Interval{.25 * kXi[i]} *
+                         (Interval{1} + Interval{kEta[i]} * eta_interval);
+    const Interval dyi = Interval{.25 * kEta[i]} *
+                         (Interval{1} + Interval{kXi[i]} * xi_interval);
+    const Vector3 gradient = dx * reciprocal1 + dy * reciprocal2;
+    const auto gradient_interval = dxi * reciprocal1i + dyi * reciprocal2i;
+    std::array<double, 3> components{};
+    for (std::size_t j = 0; j < 3; ++j) {
+      components[j] = std::ldexp(gradient[j], -geometry.length_exponent);
+      point.gradient_interval[i][j] =
+          Scale(gradient_interval[j], -geometry.length_exponent);
+    }
+    point.gradient[i] = Vector3{components};
+    if (!point.gradient[i].IsFinite() ||
+        !IsFinite(point.gradient_interval[i])) {
+      return Result<SurfacePoint>::Failure(
+          GeometryFailure(source, location, "uncertain-thermal-shell-metric",
+                          "Surface gradient enclosure is not finite."));
+    }
+  }
+  if (!(point.jacobian > 0) || !std::isfinite(point.jacobian) ||
+      !point.jacobian_interval.IsFinite() ||
+      !(point.jacobian_interval.lower > 0) || !point.axes[0].IsFinite() ||
+      !point.axes[1].IsFinite() || !point.axes[2].IsFinite()) {
+    return Result<SurfacePoint>::Failure(GeometryFailure(
+        source, location, "uncertain-thermal-shell-metric",
+        "Surface measure or frame is not finite and positive."));
+  }
+  return Result<SurfacePoint>::Success(point);
+}
+}  // namespace thermal_internal
+
+namespace {
+using Edge = std::pair<EntityIndex, EntityIndex>;
+
+/// @brief Forms a deterministic undirected source edge key.
+Edge EdgeKey(EntityIndex first, EntityIndex second) {
+  return std::minmax(first, second);
+}
+
+/// @brief Rejects possible open-edge contact only for an actual source
+/// candidate.
+Status CheckPartialEdge(const ModelDefinition& definition, EntityIndex p,
+                        EntityIndex q, EntityIndex v,
+                        const Ds4ThermalElementDefinition& owner) {
+  using namespace thermal_internal;
+  const Vector3 x{definition.nodes[p].coordinates};
+  const Vector3 y{definition.nodes[q].coordinates};
+  const Vector3 z{definition.nodes[v].coordinates};
+  double largest = 0;
+  for (std::size_t j = 0; j < 3; ++j) {
+    largest =
+        std::max({largest, std::abs(x[j]), std::abs(y[j]), std::abs(z[j])});
+  }
+  int exponent = 0;
+  std::frexp(largest, &exponent);
+  IntervalVector a;
+  IntervalVector b;
+  std::size_t axis = 0;
+  double best = 0;
+  for (std::size_t j = 0; j < 3; ++j) {
+    a[j] = Scale(Interval{y[j]}, -exponent) - Scale(Interval{x[j]}, -exponent);
+    b[j] = Scale(Interval{z[j]}, -exponent) - Scale(Interval{x[j]}, -exponent);
+    const double lower = a[j].ContainsZero() ? 0
+                                             : std::min(std::abs(a[j].lower),
+                                                        std::abs(a[j].upper));
+    if (lower > best) {
+      best = lower;
+      axis = j;
+    }
+  }
+  if (best > 0) {
+    const Interval lambda = b[axis] / a[axis];
+    if (lambda.upper < 0 || lambda.lower > 1) return Status::Ok();
+    for (std::size_t j = 0; j < 3; ++j) {
+      if (!(b[j] - lambda * a[j]).ContainsZero()) return Status::Ok();
+    }
+  }
+  return GeometryFailure(owner.source_id, owner.location,
+                         "uncertain-thermal-shell-connection",
+                         "Source-sharing coarse/partial edges have unresolved "
+                         "or nonconforming open-edge contact.");
+}
+}  // namespace
+
+Status ValidateThermalShellMesh(const ModelDefinition& definition) {
+  using namespace thermal_internal;
+  struct Incidence {
+    EntityIndex first;
+    EntityIndex second;
+    std::size_t element;
+  };
+  std::map<Edge, std::vector<Incidence>> edges;
+  std::map<EntityIndex, std::map<EntityIndex, std::vector<EntityIndex>>> links;
+  std::set<std::array<EntityIndex, 4>> patches;
+  std::map<std::pair<std::string, std::int64_t>, EntityIndex> source_nodes;
+  for (std::size_t index = 0; index < definition.thermal_elements.size();
+       ++index) {
+    const auto& element = definition.thermal_elements[index];
+    for (const EntityIndex node : element.node_indices) {
+      const auto& source = definition.nodes[node].source_id;
+      const auto inserted = source_nodes.emplace(
+          std::make_pair(source.instance_name, source.source_label), node);
+      if (!inserted.second && inserted.first->second != node) {
+        return GeometryFailure(
+            element.source_id, element.location,
+            "invalid-thermal-shell-topology",
+            "One source node identity has multiple indices.");
+      }
+    }
+    auto sorted = element.node_indices;
+    std::sort(sorted.begin(), sorted.end());
+    if (std::adjacent_find(sorted.begin(), sorted.end()) != sorted.end() ||
+        !patches.insert(sorted).second) {
+      return GeometryFailure(element.source_id, element.location,
+                             "invalid-thermal-shell-topology",
+                             "Repeated node or duplicate element patch.");
+    }
+    std::array<Vector3, 4> coordinates{};
+    for (std::size_t i = 0; i < 4; ++i)
+      coordinates[i] =
+          Vector3{definition.nodes[element.node_indices[i]].coordinates};
+    const auto geometry =
+        BuildGeometry(coordinates, element.source_id, element.location);
+    if (!geometry.HasValue()) return geometry.GetStatus();
+    // Domain never publishes geometry that fails at any kernel/recovery point.
+    const double gauss = std::sqrt(3. / 5.);
+    for (double xi : {-gauss, 0., gauss}) {
+      for (double eta : {-gauss, 0., gauss}) {
+        const auto point = EvaluateSurface(
+            geometry.Value(), xi, eta, Interval{Down(xi), Up(xi)},
+            Interval{Down(eta), Up(eta)}, element.source_id, element.location);
+        if (!point.HasValue()) return point.GetStatus();
+      }
+    }
+    const double recovery = 1 / std::sqrt(3.);
+    for (double xi : {-recovery, recovery}) {
+      for (double eta : {-recovery, recovery}) {
+        const auto point = EvaluateSurface(
+            geometry.Value(), xi, eta, Interval{Down(xi), Up(xi)},
+            Interval{Down(eta), Up(eta)}, element.source_id, element.location);
+        if (!point.HasValue()) return point.GetStatus();
+      }
+    }
+    for (std::size_t i = 0; i < 4; ++i) {
+      const EntityIndex node = element.node_indices[i];
+      const EntityIndex next = element.node_indices[(i + 1) % 4];
+      const EntityIndex previous = element.node_indices[(i + 3) % 4];
+      edges[EdgeKey(node, next)].push_back({node, next, index});
+      links[node][previous].push_back(next);
+      links[node][next].push_back(previous);
+    }
+  }
+  for (const auto& item : edges) {
+    const auto& incidence = item.second;
+    if (incidence.size() > 2 ||
+        (incidence.size() == 2 && incidence[0].first == incidence[1].first)) {
+      const auto& element =
+          definition.thermal_elements[incidence.back().element];
+      return GeometryFailure(
+          element.source_id, element.location, "invalid-thermal-shell-topology",
+          "Shared edge must have at most two opposite incidences.");
+    }
+  }
+  for (const auto& vertex : links) {
+    const auto& graph = vertex.second;
+    std::set<EntityIndex> visited;
+    std::vector<EntityIndex> pending{graph.begin()->first};
+    std::size_t ends = 0;
+    bool invalid = false;
+    for (const auto& item : graph) {
+      const auto& adjacent = item.second;
+      if (adjacent.size() == 1) ++ends;
+      if (adjacent.size() > 2 ||
+          (adjacent.size() == 2 && adjacent[0] == adjacent[1]))
+        invalid = true;
+    }
+    while (!pending.empty()) {
+      const auto node = pending.back();
+      pending.pop_back();
+      if (!visited.insert(node).second) continue;
+      for (const auto next : graph.at(node)) pending.push_back(next);
+    }
+    if (invalid || visited.size() != graph.size() || (ends != 0 && ends != 2)) {
+      const auto& incidence =
+          edges.at(EdgeKey(vertex.first, graph.begin()->first))[0];
+      const auto& element = definition.thermal_elements[incidence.element];
+      return GeometryFailure(element.source_id, element.location,
+                             "invalid-thermal-shell-topology",
+                             "Vertex link is not one path or cycle.");
+    }
+  }
+  // Enumerate only incident SOURCE edges, never all mesh coordinates vs edges.
+  std::map<EntityIndex, std::set<EntityIndex>> neighbors;
+  for (const auto& item : edges) {
+    neighbors[item.first.first].insert(item.first.second);
+    neighbors[item.first.second].insert(item.first.first);
+  }
+  for (const auto& item : edges) {
+    for (const auto endpoint : {item.first.first, item.first.second}) {
+      const EntityIndex other =
+          endpoint == item.first.first ? item.first.second : item.first.first;
+      const auto& owner = definition.thermal_elements[item.second[0].element];
+      for (const auto partial : neighbors[endpoint]) {
+        if (partial == other) continue;
+        const auto& instance =
+            definition.nodes[endpoint].source_id.instance_name;
+        if (definition.nodes[other].source_id.instance_name != instance ||
+            definition.nodes[partial].source_id.instance_name != instance)
+          continue;
+        const Status status =
+            CheckPartialEdge(definition, endpoint, other, partial, owner);
+        if (!status.IsOk()) return status;
+      }
+    }
+  }
+  return Status::Ok();
+}
+
+}  // namespace fesa

src/fesa/elements/thermal_shell_geometry.h

@@ -0,0 +1,59 @@
+#ifndef FESA_ELEMENTS_THERMAL_SHELL_GEOMETRY_H_
+#define FESA_ELEMENTS_THERMAL_SHELL_GEOMETRY_H_
+
+#include <array>
+
+#include "fesa/core/status.h"
+#include "fesa/math/vector3.h"
+#include "thermal_interval.h"
+
+namespace fesa {
+struct ModelDefinition;
+namespace thermal_internal {
+
+/// @brief Holds translated, power-of-two scaled tangents and their enclosures.
+struct Geometry {
+  Vector3 b;
+  Vector3 c;
+  Vector3 d;
+  IntervalVector bi;
+  IntervalVector ci;
+  IntervalVector di;
+  int length_exponent = 0;
+  double scaled_jacobian_lower = 0;
+};
+
+/// @brief Holds actual-point product metric, gradients and orthonormal columns.
+struct SurfacePoint {
+  std::array<double, 4> shape{};
+  std::array<Interval, 4> shape_interval{};
+  std::array<Vector3, 4> gradient{};
+  std::array<IntervalVector, 4> gradient_interval{};
+  std::array<Vector3, 3> axes{};
+  double jacobian = 0;
+  Interval jacobian_interval;
+};
+
+/// @brief Certifies whole-patch rank/injectivity from fixed candidate
+/// directions.
+Result<Geometry> BuildGeometry(const std::array<Vector3, 4>& coordinates,
+                               const SourceEntityId& source,
+                               const SourceLocation& location);
+/// @brief Checks finite metric inverse and gradients at one used surface point.
+Result<SurfacePoint> EvaluateSurface(const Geometry& geometry, double xi,
+                                     double eta, Interval xi_interval,
+                                     Interval eta_interval,
+                                     const SourceEntityId& source,
+                                     const SourceLocation& location);
+/// @brief Builds a geometry failure retaining original source element identity.
+Status GeometryFailure(const SourceEntityId& source,
+                       const SourceLocation& location, const char* code,
+                       const char* message);
+}  // namespace thermal_internal
+
+/// @brief Validates DS4 whole-patch geometry and source-only manifold topology.
+/// @pre Node indices and section assignments have been checked by Domain.
+Status ValidateThermalShellMesh(const ModelDefinition& definition);
+
+}  // namespace fesa
+#endif  // FESA_ELEMENTS_THERMAL_SHELL_GEOMETRY_H_

src/fesa/model/domain.cpp

@@ -8,6 +8,8 @@
 #include <string>
 #include <utility>
 
+#include "../elements/thermal_shell_geometry.h"
+
 namespace fesa {
 namespace {
 
@@ -112,6 +114,8 @@ Status ValidateThermalDefinition(const ModelDefinition& definition) {
       active_nodes[node] = true;
     }
   }
+  const Status geometry_status = ValidateThermalShellMesh(definition);
+  if (!geometry_status.IsOk()) return geometry_status;
   std::map<EntityIndex, std::array<double, 5>> initial_values;
   for (const auto& initial : definition.thermal_initial_conditions) {
     if (initial.node_index >= active_nodes.size() ||

tests/unit/analysis/thermal_time_grid_test.cpp

@@ -0,0 +1,108 @@
+#include "fesa/analysis/thermal_time_grid.h"
+
+#include <gtest/gtest.h>
+
+#include <cmath>
+#include <limits>
+
+TEST(ThermalTimeGrid, SingleStepUsesPeriodAsActualIncrement) {
+  const auto grid = fesa::ThermalTimeGrid::Create(2.0, 3.0, 1U);
+  ASSERT_TRUE(grid.HasValue());
+  EXPECT_EQ(grid.Value().Count(), 1U);
+  EXPECT_EQ(grid.Value().Delta(1U), 2.0);
+}
+
+TEST(ThermalTimeGrid, PreservesGenuineRemainderOutsideInputCells) {
+  const auto grid =
+      fesa::ThermalTimeGrid::Create(std::nextafter(.9, 1.0), .3, 4U);
+  ASSERT_TRUE(grid.HasValue());
+  EXPECT_EQ(grid.Value().Count(), 4U);
+  EXPECT_EQ(grid.Value().Delta(4U), 1.6653345369377348e-16);
+}
+
+TEST(ThermalTimeGrid, OriginalSubnormalCellsRespectOpenAndClosedContact) {
+  const double tiny = std::numeric_limits<double>::denorm_min();
+  const auto even = fesa::ThermalTimeGrid::Create(20 * tiny, 6 * tiny, 4U);
+  const auto odd = fesa::ThermalTimeGrid::Create(17 * tiny, 5 * tiny, 4U);
+  ASSERT_TRUE(even.HasValue());
+  ASSERT_TRUE(odd.HasValue());
+  EXPECT_EQ(even.Value().Count(), 3U);
+  EXPECT_EQ(odd.Value().Count(), 4U);
+  EXPECT_EQ(odd.Value().Delta(4U), 2 * tiny);
+}
+
+TEST(ThermalTimeGrid, DecimalCellsCounterexampleAndScale) {
+  for (const int exponent : {-1000, 0, 1000}) {
+    const double period = std::ldexp(.9, exponent);
+    const double delta = std::ldexp(.3, exponent);
+    const auto grid = fesa::ThermalTimeGrid::Create(period, delta, 3U);
+    ASSERT_TRUE(grid.HasValue());
+    EXPECT_EQ(grid.Value().Count(), 3U);
+    EXPECT_EQ(grid.Value().Delta(3U), delta);
+    EXPECT_EQ(grid.Value().Time(3U), period);
+  }
+  const auto old_counterexample =
+      fesa::ThermalTimeGrid::Create(5.351455173733616, .17262758624947144, 31U);
+  ASSERT_TRUE(old_counterexample.HasValue());
+  EXPECT_EQ(old_counterexample.Value().Count(), 31U);
+  EXPECT_EQ(old_counterexample.Value().Delta(31U), .17262758624947144);
+  EXPECT_FALSE(
+      fesa::ThermalTimeGrid::Create(5.351455173733616, .17262758624947144, 30U)
+          .HasValue());
+}
+
+TEST(ThermalTimeGrid, BinadeLimitsAndMaximumFinitePreserveOriginalCells) {
+  const double max = std::numeric_limits<double>::max();
+  const auto grid = fesa::ThermalTimeGrid::Create(max, max / 3, 3U);
+  ASSERT_TRUE(grid.HasValue());
+  EXPECT_EQ(grid.Value().Count(), 3U);
+  EXPECT_EQ(grid.Value().Delta(3U), max / 3);
+  const auto upper =
+      fesa::ThermalTimeGrid::Create(4294967295., 1., 4294967295U);
+  ASSERT_TRUE(upper.HasValue());
+  EXPECT_EQ(upper.Value().Count(), 4294967295U);
+  EXPECT_TRUE(upper.Value().Frame(4294967295U).HasValue());
+  EXPECT_FALSE(
+      fesa::ThermalTimeGrid::Create(4294967296., 1., 4294967295U).HasValue());
+  EXPECT_FALSE(fesa::ThermalTimeGrid::Create(
+                   max, std::numeric_limits<double>::denorm_min(), 4294967295U)
+                   .HasValue());
+}
+
+TEST(ThermalTimeGrid, ShortenedIncrementAndCheckedFrames) {
+  const auto grid = fesa::ThermalTimeGrid::Create(10., 3., 4U);
+  ASSERT_TRUE(grid.HasValue());
+  EXPECT_EQ(grid.Value().Count(), 4U);
+  EXPECT_EQ(grid.Value().Time(0U), 0.);
+  for (std::uint32_t i = 1; i <= 4; ++i) {
+    ASSERT_TRUE(grid.Value().Frame(i).HasValue());
+    EXPECT_EQ(grid.Value().Time(i), i == 4 ? 10. : 3. * i);
+    EXPECT_EQ(grid.Value().Delta(i), i == 4 ? 1. : 3.);
+  }
+  EXPECT_FALSE(grid.Value().Frame(0U).HasValue());
+  EXPECT_FALSE(grid.Value().Frame(5U).HasValue());
+  const auto decimals = fesa::ThermalTimeGrid::Create(1., .1);
+  ASSERT_TRUE(decimals.HasValue());
+  EXPECT_EQ(decimals.Value().Time(3U), .30000000000000004);
+  EXPECT_EQ(decimals.Value().Delta(10U), .1);
+  const double below = std::nextafter(.9, 0.);
+  const auto close = fesa::ThermalTimeGrid::Create(below, .3, 3U);
+  const auto short_last =
+      fesa::ThermalTimeGrid::Create(std::nextafter(below, 0.), .3, 3U);
+  ASSERT_TRUE(close.HasValue());
+  ASSERT_TRUE(short_last.HasValue());
+  EXPECT_EQ(close.Value().Delta(3U), .3);
+  EXPECT_EQ(short_last.Value().Delta(3U), .2999999999999998);
+}
+
+TEST(ThermalTimeGrid, RejectsInvalidControlsBeforeIntegerConversion) {
+  for (double bad : {0., -1., std::numeric_limits<double>::infinity(),
+                     std::numeric_limits<double>::quiet_NaN()}) {
+    EXPECT_FALSE(fesa::ThermalTimeGrid::Create(bad, 1.).HasValue());
+    EXPECT_FALSE(fesa::ThermalTimeGrid::Create(1., bad).HasValue());
+  }
+  EXPECT_FALSE(fesa::ThermalTimeGrid::Create(1., 1., 0U).HasValue());
+  const auto equal = fesa::ThermalTimeGrid::Create(2., 2., 1U);
+  ASSERT_TRUE(equal.HasValue());
+  EXPECT_EQ(equal.Value().Delta(1U), 2.);
+}

tests/unit/elements/thermal_shell_test.cpp

@@ -0,0 +1,498 @@
+#include "fesa/elements/thermal_shell.h"
+
+#include <gtest/gtest.h>
+
+#include <algorithm>
+#include <array>
+#include <cmath>
+#include <limits>
+#include <numeric>
+
+namespace {
+using fesa::Matrix;
+using fesa::Vector3;
+using Patch = std::array<Vector3, 4>;
+const Patch kSquare{{{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}}};
+
+double Gamma(int operations) {
+  const double product = operations * std::ldexp(1., -53);
+  return product / (1 - product);
+}
+double Norm(const Matrix& matrix) {
+  double norm = 0;
+  for (std::size_t i = 0; i < matrix.Rows(); ++i)
+    for (std::size_t j = 0; j < matrix.Columns(); ++j)
+      norm = std::hypot(norm, matrix(i, j));
+  return norm;
+}
+double Difference(const Matrix& a, const Matrix& b) {
+  double norm = 0;
+  for (std::size_t i = 0; i < a.Rows(); ++i)
+    for (std::size_t j = 0; j < a.Columns(); ++j)
+      norm = std::hypot(norm, a(i, j) - b(i, j));
+  return norm;
+}
+fesa::Result<fesa::ThermalShell> Kernel(const Patch& patch = kSquare,
+                                        double h = 2, double k = 1) {
+  return fesa::ThermalShell::Create({"patch", 19, "0019"}, {{7, 2, 9, 4}},
+                                    patch, k, 1, 1, h, {"kernel.inp", 42});
+}
+struct Oracle {
+  Matrix c{20, 20};
+  Matrix k{20, 20};
+  Matrix pc{20, 20};
+  Matrix pk{20, 20};
+};
+
+// Independently integrated exact rational unit-square matrices, h variable.
+// The signed products are kept separate from their absolute contribution sums.
+Oracle Rational(double h) {
+  constexpr int mass[4][4] = {
+      {4, 2, 1, 2}, {2, 4, 2, 1}, {1, 2, 4, 2}, {2, 1, 2, 4}};
+  constexpr int stiffness[4][4] = {
+      {4, -1, -2, -1}, {-1, 4, -1, -2}, {-2, -1, 4, -1}, {-1, -2, -1, 4}};
+  constexpr int z0[5][5] = {{4, 2, -1, 0, 0},
+                            {2, 16, 2, 0, 0},
+                            {-1, 2, 8, 2, -1},
+                            {0, 0, 2, 16, 2},
+                            {0, 0, -1, 2, 4}};
+  constexpr int z1[5][5] = {{7, -8, 1, 0, 0},
+                            {-8, 16, -8, 0, 0},
+                            {1, -8, 14, -8, 1},
+                            {0, 0, -8, 16, -8},
+                            {0, 0, 1, -8, 7}};
+  Oracle result;
+  for (std::size_t p = 0; p < 20; ++p)
+    for (std::size_t q = 0; q < 20; ++q) {
+      const auto i = p / 5, j = q / 5, a = p % 5, b = q % 5;
+      result.c(p, q) = h * mass[i][j] * z0[a][b] / 2160.;
+      const double tangential = h * stiffness[i][j] * z0[a][b] / 360.;
+      const double normal = mass[i][j] * z1[a][b] / (54. * h);
+      result.k(p, q) = tangential + normal;
+      result.pc(p, q) = std::abs(result.c(p, q));
+      const double sx[4] = {-1, 1, 1, -1}, sy[4] = {-1, -1, 1, 1};
+      const double ps = (sx[i] == sx[j] ? 1. / 3 : 1. / 6) +
+                        (sy[i] == sy[j] ? 1. / 3 : 1. / 6);
+      result.pk(p, q) = h * ps * std::abs(z0[a][b]) / 60 + std::abs(normal);
+    }
+  return result;
+}
+
+// Independent 54-point B^T B/H^T H oracle: no exact-Z or production helper.
+Oracle Direct54(const Patch& coordinates, double h) {
+  Oracle result;
+  const double g = std::sqrt(.6);
+  const std::array<double, 3> gauss{{-g, 0, g}},
+      weight{{5. / 9, 8. / 9, 5. / 9}};
+  constexpr double sx[4] = {-1, 1, 1, -1}, sy[4] = {-1, -1, 1, 1};
+  for (std::size_t ix = 0; ix < 3; ++ix)
+    for (std::size_t iy = 0; iy < 3; ++iy) {
+      std::array<double, 4> n{}, dx{}, dy{};
+      Vector3 a1, a2;
+      for (std::size_t i = 0; i < 4; ++i) {
+        n[i] = (1 + sx[i] * gauss[ix]) * (1 + sy[i] * gauss[iy]) / 4;
+        dx[i] = sx[i] * (1 + sy[i] * gauss[iy]) / 4;
+        dy[i] = sy[i] * (1 + sx[i] * gauss[ix]) / 4;
+        a1 = a1 + dx[i] * coordinates[i];
+        a2 = a2 + dy[i] * coordinates[i];
+      }
+      const double aa = a1.Dot(a1), bb = a2.Dot(a2), ab = a1.Dot(a2);
+      const double det = aa * bb - ab * ab;
+      const Vector3 normal = a1.Cross(a2) / std::sqrt(det);
+      const Vector3 reciprocal1 = (bb * a1 - ab * a2) / det;
+      const Vector3 reciprocal2 = (aa * a2 - ab * a1) / det;
+      for (std::size_t segment = 0; segment < 2; ++segment)
+        for (std::size_t iz = 0; iz < 3; ++iz) {
+          const double s = gauss[iz];
+          const std::array<double, 3> l{
+              {s * (s - 1) / 2, 1 - s * s, s * (s + 1) / 2}};
+          const std::array<double, 3> dl{{2 * s - 1, -4 * s, 2 * s + 1}};
+          std::array<double, 20> shape{};
+          std::array<Vector3, 20> gradient{};
+          const double w =
+              weight[ix] * weight[iy] * weight[iz] * std::sqrt(det) * h / 4;
+          for (std::size_t i = 0; i < 4; ++i)
+            for (std::size_t a = 0; a < 3; ++a) {
+              const auto p = 5 * i + 2 * segment + a;
+              shape[p] = n[i] * l[a];
+              gradient[p] = l[a] * (dx[i] * reciprocal1 + dy[i] * reciprocal2) +
+                            (2 * n[i] * dl[a] / h) * normal;
+            }
+          for (std::size_t p = 0; p < 20; ++p)
+            for (std::size_t q = 0; q < 20; ++q) {
+              const double c = shape[p] * shape[q] * w;
+              const double k = gradient[p].Dot(gradient[q]) * w;
+              result.c(p, q) += c;
+              result.k(p, q) += k;
+              result.pc(p, q) += std::abs(c);
+              // Component products, not absolute signed dot products.
+              for (std::size_t axis = 0; axis < 3; ++axis)
+                result.pk(p, q) +=
+                    std::abs(gradient[p][axis] * gradient[q][axis] * w);
+            }
+        }
+    }
+  return result;
+}
+
+// Direct recovery operator and absolute propagation scale, independent of the
+// production cache. Rows preserve the two one-sided center derivatives.
+struct FluxOracle {
+  std::array<Vector3, 24> value{};
+  std::array<double, 24> scale{};
+};
+FluxOracle Flux(const Patch& patch, double h, const std::array<double, 20>& t) {
+  FluxOracle result;
+  const double g = 1 / std::sqrt(3.);
+  const double sx[4] = {-1, 1, 1, -1}, sy[4] = {-1, -1, 1, 1};
+  for (std::size_t ip = 0; ip < 4; ++ip) {
+    std::array<double, 4> n{}, dx{}, dy{};
+    Vector3 a1, a2;
+    for (std::size_t i = 0; i < 4; ++i) {
+      n[i] = (1 + sx[i] * sx[ip] * g) * (1 + sy[i] * sy[ip] * g) / 4;
+      dx[i] = sx[i] * (1 + sy[i] * sy[ip] * g) / 4;
+      dy[i] = sy[i] * (1 + sx[i] * sx[ip] * g) / 4;
+      a1 = a1 + dx[i] * patch[i];
+      a2 = a2 + dy[i] * patch[i];
+    }
+    const double aa = a1.Dot(a1), bb = a2.Dot(a2), ab = a1.Dot(a2),
+                 det = aa * bb - ab * ab;
+    const Vector3 normal = a1.Cross(a2) / std::sqrt(det);
+    const Vector3 r1 = (bb * a1 - ab * a2) / det,
+                  r2 = (aa * a2 - ab * a1) / det;
+    for (std::size_t trace = 0; trace < 6; ++trace) {
+      const double x = static_cast<double>(trace % 3) - 1;
+      const double l[3] = {x * (x - 1) / 2, 1 - x * x, x * (x + 1) / 2};
+      const double dl[3] = {2 * x - 1, -4 * x, 2 * x + 1};
+      const auto row = 6 * ip + trace;
+      for (std::size_t i = 0; i < 4; ++i)
+        for (std::size_t a = 0; a < 3; ++a) {
+          const auto p = 5 * i + (trace < 3 ? 0 : 2) + a;
+          const Vector3 first = (l[a] * dx[i]) * r1;
+          const Vector3 second = (l[a] * dy[i]) * r2;
+          const Vector3 third = (2 * n[i] * dl[a] / h) * normal;
+          result.value[row] =
+              result.value[row] - t[p] * (first + second + third);
+          result.scale[row] +=
+              std::abs(t[p]) * (first.Norm() + second.Norm() + third.Norm());
+        }
+    }
+  }
+  return result;
+}
+
+// Jacobi orthogonal diagonalization plus a posteriori residual and
+// orthogonality budget. No guessed eigenvalue epsilon is used for the rank
+// decision.
+std::pair<double, double> EigenExtrema(const Matrix& original,
+                                       double& backward) {
+  Matrix a = original;
+  const auto n = a.Rows();
+  Matrix v{n, n};
+  for (std::size_t i = 0; i < n; ++i) v(i, i) = 1;
+  for (int sweep = 0; sweep < 100; ++sweep) {
+    for (std::size_t p = 0; p < n; ++p)
+      for (std::size_t q = p + 1; q < n; ++q) {
+        if (a(p, q) == 0) continue;
+        const double tau = (a(q, q) - a(p, p)) / (2 * a(p, q));
+        const double t =
+            std::copysign(1., tau) / (std::abs(tau) + std::hypot(1., tau));
+        const double c = 1 / std::hypot(1., t), s = t * c;
+        const double off = a(p, q);
+        a(p, p) -= t * off;
+        a(q, q) += t * off;
+        a(p, q) = a(q, p) = 0;
+        for (std::size_t j = 0; j < n; ++j) {
+          if (j != p && j != q) {
+            const double x = a(j, p), y = a(j, q);
+            a(j, p) = a(p, j) = c * x - s * y;
+            a(j, q) = a(q, j) = s * x + c * y;
+          }
+          const double x = v(j, p), y = v(j, q);
+          v(j, p) = c * x - s * y;
+          v(j, q) = s * x + c * y;
+        }
+      }
+  }
+  double residual = 0, orthogonality = 0;
+  for (std::size_t i = 0; i < n; ++i)
+    for (std::size_t j = 0; j < n; ++j) {
+      double av = 0, vv = 0, absolute = 0;
+      for (std::size_t k = 0; k < n; ++k) {
+        av += original(i, k) * v(k, j);
+        absolute += std::abs(original(i, k) * v(k, j));
+        vv += v(k, i) * v(k, j);
+      }
+      residual = std::hypot(residual,
+                            std::abs(av - v(i, j) * a(j, j)) +
+                                Gamma(2 * static_cast<int>(n) + 2) *
+                                    (absolute + std::abs(v(i, j) * a(j, j))));
+      orthogonality =
+          std::hypot(orthogonality, std::abs(vv - (i == j ? 1. : 0.)) +
+                                        Gamma(2 * static_cast<int>(n) + 2));
+    }
+  EXPECT_LT(orthogonality, 1.);
+  backward =
+      (residual + 2 * Norm(original) * orthogonality) / (1 - orthogonality);
+  double low = a(0, 0), high = a(0, 0);
+  for (std::size_t i = 1; i < n; ++i) {
+    low = std::min(low, a(i, i));
+    high = std::max(high, a(i, i));
+  }
+  return {low, high};
+}
+}  // namespace
+
+TEST(ThermalShell, RationalMatricesSymmetryAndConstantModeUseFixedBudgets) {
+  const auto created = Kernel();
+  ASSERT_TRUE(created.HasValue());
+  const auto& kernel = created.Value();
+  const auto exact = Rational(2);
+  const double ec = Gamma(4096) * Norm(exact.pc),
+               ek = Gamma(4096) * Norm(exact.pk);
+  EXPECT_LE(Difference(kernel.Capacity(), exact.c), ec);
+  EXPECT_LE(Difference(kernel.Conductivity(), exact.k), ek);
+  double cs = 0, ks = 0, null = 0, total = 0;
+  for (std::size_t i = 0; i < 20; ++i) {
+    double row = 0;
+    for (std::size_t j = 0; j < 20; ++j) {
+      cs = std::hypot(cs, kernel.Capacity()(i, j) - kernel.Capacity()(j, i));
+      ks = std::hypot(
+          ks, kernel.Conductivity()(i, j) - kernel.Conductivity()(j, i));
+      row += kernel.Conductivity()(i, j);
+      total += kernel.Capacity()(i, j);
+    }
+    null = std::hypot(null, row);
+  }
+  EXPECT_LE(cs, 2 * ec);
+  EXPECT_LE(ks, 2 * ek);
+  EXPECT_LE(null, std::sqrt(20.) * ek);
+  EXPECT_LE(std::abs(total - 2), 20 * ec);
+  EXPECT_EQ(kernel.SourceId().source_label_text, "0019");
+  EXPECT_EQ(kernel.NodeIndices(),
+            (std::array<fesa::EntityIndex, 4>{{7, 2, 9, 4}}));
+}
+
+TEST(ThermalShell, CapacityIsPositiveAndConductionHasExactlyOneNullMode) {
+  const auto created = Kernel();
+  ASSERT_TRUE(created.HasValue());
+  const auto exact = Rational(2);
+  double backward = 0;
+  const auto capacity = EigenExtrema(created.Value().Capacity(), backward);
+  EXPECT_GT(capacity.first, Gamma(4096) * Norm(exact.pc) + backward);
+  const auto conduction =
+      EigenExtrema(created.Value().Conductivity(), backward);
+  EXPECT_LE(std::abs(conduction.first),
+            Gamma(4096) * Norm(exact.pk) + backward);
+  // Helmert columns form an independent orthonormal basis of 1-perp.
+  Matrix basis{20, 19};
+  for (std::size_t j = 0; j < 19; ++j) {
+    const double norm = std::sqrt(static_cast<double>((j + 1) * (j + 2)));
+    for (std::size_t i = 0; i <= j; ++i) basis(i, j) = 1 / norm;
+    basis(j + 1, j) = -static_cast<double>(j + 1) / norm;
+  }
+  Matrix restricted{19, 19};
+  for (std::size_t i = 0; i < 19; ++i)
+    for (std::size_t j = 0; j < 19; ++j)
+      for (std::size_t p = 0; p < 20; ++p)
+        for (std::size_t q = 0; q < 20; ++q)
+          restricted(i, j) +=
+              basis(p, i) * created.Value().Conductivity()(p, q) * basis(q, j);
+  const auto positive = EigenExtrema(restricted, backward);
+  const double projection_budget = Gamma(8192) * Norm(exact.pk) * 20;
+  EXPECT_GT(positive.first,
+            Gamma(4096) * Norm(exact.pk) + backward + projection_budget);
+}
+
+TEST(ThermalShell,
+     Direct54MatchesSeparatedOperatorsOnActualThreeDimensionalPatches) {
+  const std::array<Patch, 5> patches{
+      {kSquare,
+       {{{0, 0, 0}, {2, 0, 0}, {1.5, 1, 0}, {0, 1, 0}}},
+       {{{-1, -1, .1}, {1, -1, -.1}, {1, 1, .1}, {-1, 1, -.1}}},
+       {{{-1, -1, .5}, {1, -1, -.5}, {1, 1, .5}, {-1, 1, -.5}}},
+       {{{-1, -1, 2}, {1, -1, -2}, {1, 1, 2}, {-1, 1, -2}}}}};
+  for (const auto& patch : patches) {
+    const auto created = Kernel(patch);
+    ASSERT_TRUE(created.HasValue());
+    const auto direct = Direct54(patch, 2);
+    const auto& kernel = created.Value();
+    // General-geometry separation error is the propagated interval enclosure;
+    // direct 54's independent arithmetic has the prescribed gamma8192 budget.
+    EXPECT_LE(Difference(kernel.Capacity(), direct.c),
+              Norm(kernel.CapacityError()) + Gamma(8192) * Norm(direct.pc));
+    EXPECT_LE(Difference(kernel.Conductivity(), direct.k),
+              Norm(kernel.ConductivityError()) + Gamma(8192) * Norm(direct.pk));
+  }
+}
+
+TEST(ThermalShell, ThicknessScalingPreservesBothConductionTermsAndCapacity) {
+  for (double h : {.5, 1., 2.}) {
+    const auto created = Kernel(kSquare, h);
+    ASSERT_TRUE(created.HasValue());
+    const auto exact = Rational(h);
+    EXPECT_LE(Difference(created.Value().Capacity(), exact.c),
+              Gamma(4096) * Norm(exact.pc));
+    EXPECT_LE(Difference(created.Value().Conductivity(), exact.k),
+              Gamma(4096) * Norm(exact.pk));
+  }
+}
+
+TEST(ThermalShell, AffineLinearFluxFramesAndSeparateMidsurfaceTraces) {
+  const auto created = Kernel();
+  ASSERT_TRUE(created.HasValue());
+  std::array<double, 20> temperature{};
+  for (std::size_t i = 0; i < 4; ++i)
+    for (std::size_t a = 0; a < 5; ++a)
+      temperature[5 * i + a] =
+          2 + 3 * kSquare[i].X() - 2 * kSquare[i].Y() + 4 * (-1 + .5 * a);
+  const auto recovery = created.Value().Recover(temperature);
+  ASSERT_TRUE(recovery.HasValue());
+  const auto linear_scale = Flux(kSquare, 2, temperature);
+  for (std::size_t row = 0; row < 24; ++row)
+    EXPECT_LE(
+        (recovery.Value().global_flux[row] - Vector3{-3, 2, -4}).Norm(),
+        Gamma(2048) * (linear_scale.scale[row] + Vector3{3, -2, 4}.Norm()));
+  for (const auto& frame : recovery.Value().frames) {
+    for (std::size_t i = 0; i < 3; ++i)
+      for (std::size_t j = 0; j < 3; ++j)
+        EXPECT_NEAR(frame[i].Dot(frame[j]), i == j ? 1. : 0., Gamma(2048));
+    EXPECT_NEAR(frame[0].Cross(frame[1]).Dot(frame[2]), 1., Gamma(2048));
+  }
+  for (std::size_t i = 0; i < 4; ++i)
+    for (std::size_t a = 0; a < 5; ++a)
+      temperature[5 * i + a] = a < 3 ? 0. : static_cast<double>(a - 2);
+  const auto traces = created.Value().Recover(temperature);
+  ASSERT_TRUE(traces.HasValue());
+  for (std::size_t ip = 0; ip < 4; ++ip) {
+    EXPECT_EQ(traces.Value().global_flux[6 * ip + 2].Z(), 0.);
+    EXPECT_NEAR(traces.Value().global_flux[6 * ip + 3].Z(), -2.,
+                Gamma(2048) * 2);
+  }
+  temperature.fill(2);
+  const auto constant = created.Value().Recover(temperature);
+  ASSERT_TRUE(constant.HasValue());
+  const auto constant_scale = Flux(kSquare, 2, temperature);
+  for (std::size_t row = 0; row < 24; ++row)
+    EXPECT_LE(constant.Value().global_flux[row].Norm(),
+              Gamma(2048) * constant_scale.scale[row]);
+}
+
+TEST(ThermalShell, TranslationRotationAndReversedNormalPreservePhysicalFlux) {
+  const auto base = Kernel();
+  ASSERT_TRUE(base.HasValue());
+  Patch translated = kSquare, rotated = kSquare, reversed{};
+  const auto rotate = [](Vector3 p) { return Vector3{p.Z(), p.X(), p.Y()}; };
+  for (std::size_t i = 0; i < 4; ++i) {
+    translated[i] = kSquare[i] + Vector3{1048576, -1048576, 1048576};
+    rotated[i] = rotate(kSquare[i]);
+    reversed[i] = kSquare[(4 - i) % 4];
+  }
+  const auto moved = Kernel(translated), turned = Kernel(rotated),
+             flipped = Kernel(reversed);
+  ASSERT_TRUE(moved.HasValue());
+  ASSERT_TRUE(turned.HasValue());
+  ASSERT_TRUE(flipped.HasValue());
+  EXPECT_LE(
+      Difference(base.Value().Capacity(), moved.Value().Capacity()),
+      Norm(base.Value().CapacityError()) + Norm(moved.Value().CapacityError()));
+  EXPECT_LE(
+      Difference(base.Value().Conductivity(), turned.Value().Conductivity()),
+      Norm(base.Value().ConductivityError()) +
+          Norm(turned.Value().ConductivityError()));
+  std::array<double, 20> temperature{}, reverse_temperature{};
+  for (std::size_t i = 0; i < 4; ++i)
+    for (std::size_t a = 0; a < 5; ++a) {
+      temperature[5 * i + a] =
+          2 + 3 * kSquare[i].X() - 2 * kSquare[i].Y() + 4 * (-1 + .5 * a);
+      reverse_temperature[5 * i + a] =
+          2 + 3 * reversed[i].X() - 2 * reversed[i].Y() - 4 * (-1 + .5 * a);
+    }
+  const auto flux = turned.Value().Recover(temperature),
+             reverse = flipped.Value().Recover(reverse_temperature);
+  ASSERT_TRUE(flux.HasValue());
+  ASSERT_TRUE(reverse.HasValue());
+  const auto turn_scale = Flux(rotated, 2, temperature);
+  const auto reverse_scale = Flux(reversed, 2, reverse_temperature);
+  for (std::size_t row = 0; row < 24; ++row) {
+    EXPECT_LE((flux.Value().global_flux[row] - rotate({-3, 2, -4})).Norm(),
+              Gamma(2048) * (turn_scale.scale[row] + Vector3{3, -2, 4}.Norm()));
+    EXPECT_LE(
+        (reverse.Value().global_flux[row] - Vector3{-3, 2, -4}).Norm(),
+        Gamma(2048) * (reverse_scale.scale[row] + Vector3{3, -2, 4}.Norm()));
+  }
+}
+
+TEST(ThermalShell, CheckedCreationAndRecoveryRetainSourceAndRejectNonfinite) {
+  for (double bad : {0., -1., std::numeric_limits<double>::infinity(),
+                     std::numeric_limits<double>::quiet_NaN()}) {
+    EXPECT_FALSE(Kernel(kSquare, bad).HasValue());
+    EXPECT_FALSE(Kernel(kSquare, 2, bad).HasValue());
+  }
+  const auto created = Kernel();
+  ASSERT_TRUE(created.HasValue());
+  std::array<double, 20> temperature{};
+  temperature[3] = std::numeric_limits<double>::quiet_NaN();
+  const auto bad = created.Value().Recover(temperature);
+  ASSERT_FALSE(bad.HasValue());
+  EXPECT_EQ(bad.GetStatus().Diagnostics()[0].entity_identity, "patch:0019");
+  EXPECT_EQ(bad.GetStatus().Diagnostics()[0].location.line, 42U);
+  temperature.fill(std::numeric_limits<double>::max());
+  temperature[0] = -temperature[0];
+  EXPECT_FALSE(created.Value().Recover(temperature).HasValue());
+}
+
+TEST(ThermalShell, GeneralWarpedFluxUsesActualReciprocalTangentsAndTraceOrder) {
+  const Patch saddle{{{-1, -1, 2}, {1, -1, -2}, {1, 1, 2}, {-1, 1, -2}}};
+  const auto created = Kernel(saddle, 1);
+  ASSERT_TRUE(created.HasValue());
+  std::array<double, 20> temperature{};
+  for (std::size_t p = 0; p < 20; ++p)
+    temperature[p] = static_cast<double>((p * 7) % 11) - 4;
+  const auto actual = created.Value().Recover(temperature);
+  ASSERT_TRUE(actual.HasValue());
+  const auto expected = Flux(saddle, 1, temperature);
+  for (std::size_t row = 0; row < 24; ++row)
+    EXPECT_LE((actual.Value().global_flux[row] - expected.value[row]).Norm(),
+              2 * Gamma(2048) * expected.scale[row]);
+  // Local normals vary around this saddle; recovery must retain those frames.
+  EXPECT_GT((actual.Value().frames[0][2] - actual.Value().frames[1][2]).Norm(),
+            1.);
+}
+
+TEST(ThermalShell, PowerOfTwoCoordinateScalingPreservesDimensionedOperators) {
+  for (int exponent : {-200, 200}) {
+    Patch patch = kSquare;
+    const double scale = std::ldexp(1., exponent);
+    for (auto& point : patch) point = scale * point;
+    const auto created = Kernel(patch, 2);
+    ASSERT_TRUE(created.HasValue());
+    const auto unit = Rational(2);
+    const auto very_thick = Rational(2 / scale);
+    // Area multiplies C by s^2; surface conduction is invariant, while the
+    // normal term multiplies by s^2. Equivalent expression uses h/s and k*s.
+    Matrix c{20, 20}, k{20, 20};
+    for (std::size_t p = 0; p < 20; ++p)
+      for (std::size_t q = 0; q < 20; ++q) {
+        c(p, q) = std::ldexp(unit.c(p, q), 2 * exponent);
+        k(p, q) = scale * very_thick.k(p, q);
+      }
+    EXPECT_LE(Difference(created.Value().Capacity(), c),
+              Norm(created.Value().CapacityError()) + Gamma(4096) * Norm(c));
+    EXPECT_LE(
+        Difference(created.Value().Conductivity(), k),
+        Norm(created.Value().ConductivityError()) + Gamma(4096) * Norm(k));
+  }
+}
+
+TEST(ThermalShell, RejectsUnrepresentablePositivePropertiesWithoutClamping) {
+  const auto underflow =
+      fesa::ThermalShell::Create({"patch", 19, "0019"}, {{0, 1, 2, 3}}, kSquare,
+                                 1, std::numeric_limits<double>::denorm_min(),
+                                 std::numeric_limits<double>::denorm_min(), 1);
+  EXPECT_FALSE(underflow.HasValue());
+  const auto overflow =
+      fesa::ThermalShell::Create({"patch", 19, "0019"}, {{0, 1, 2, 3}}, kSquare,
+                                 1, std::numeric_limits<double>::max(), 2, 1);
+  EXPECT_FALSE(overflow.HasValue());
+}

tests/unit/model/thermal_shell_geometry_test.cpp

@@ -0,0 +1,212 @@
+#include <gtest/gtest.h>
+
+#include <array>
+#include <cmath>
+#include <limits>
+#include <string>
+
+#include "fesa/elements/thermal_shell.h"
+#include "fesa/model/domain.h"
+
+namespace {
+using fesa::Vector3;
+using Patch = std::array<Vector3, 4>;
+const Patch kSquare{{{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}}};
+
+fesa::Result<fesa::ThermalShell> Kernel(const Patch& patch) {
+  return fesa::ThermalShell::Create({"part", 27, "0027"}, {{0, 1, 2, 3}}, patch,
+                                    1., 1., 1., 1., {"mesh.inp", 17});
+}
+
+fesa::ModelDefinition Mesh(
+    const std::vector<Vector3>& nodes,
+    const std::vector<std::array<fesa::EntityIndex, 4>>& elements) {
+  fesa::ModelDefinition model{};
+  for (std::size_t i = 0; i < nodes.size(); ++i) {
+    model.nodes.push_back(
+        {{"part", static_cast<std::int64_t>(i + 1), std::to_string(i + 1)},
+         nodes[i].Components(),
+         {"mesh.inp", i + 1}});
+  }
+  model.thermal_materials.push_back({{"part", 1, "1"}, "heat", 1, 1, 1, {}});
+  model.thermal_sections.push_back({{"part", 1, "1"}, "section", 1, 0, {}});
+  for (std::size_t i = 0; i < elements.size(); ++i) {
+    model.thermal_elements.push_back(
+        {{"part", static_cast<std::int64_t>(i + 1), std::to_string(i + 1)},
+         elements[i],
+         0,
+         {"mesh.inp", i + 20}});
+  }
+  model.thermal_step = fesa::ThermalStepDefinition{"heat", .1, 1., 10};
+  return model;
+}
+
+bool Accepted(fesa::ModelDefinition model) {
+  return fesa::Domain::Create(std::move(model)).HasValue();
+}
+}  // namespace
+
+TEST(ThermalShellGeometry, CertifiesAffineDistortedSaddleAndSlenderPatches) {
+  EXPECT_TRUE(Kernel(kSquare).HasValue());
+  EXPECT_TRUE(
+      Kernel({{{0, 0, 0}, {2, 0, 0}, {3, 1, 0}, {1, 1, 0}}}).HasValue());
+  EXPECT_TRUE(
+      Kernel({{{0, 0, 0}, {2, 0, 0}, {1.5, 1, 0}, {0, 1, 0}}}).HasValue());
+  EXPECT_TRUE(Kernel({{{0, 0, 0},
+                       {1, 0, 0},
+                       {1, std::ldexp(1., -40), 0},
+                       {0, std::ldexp(1., -40), 0}}})
+                  .HasValue());
+  for (double a : {.1, .5, 2.}) {
+    EXPECT_TRUE(
+        Kernel({{{-1, -1, a}, {1, -1, -a}, {1, 1, a}, {-1, 1, -a}}}).HasValue())
+        << a;
+  }
+}
+
+TEST(ThermalShellGeometry, PowerOfTwoScalingTranslationAndRotation) {
+  for (int exponent : {-200, 0, 200}) {
+    Patch patch = kSquare;
+    for (auto& point : patch) point = std::ldexp(1., exponent) * point;
+    EXPECT_TRUE(Kernel(patch).HasValue()) << exponent;
+  }
+  Patch translated = kSquare;
+  Patch rotated = kSquare;
+  for (std::size_t i = 0; i < 4; ++i) {
+    translated[i] = translated[i] + Vector3{1048576., -1048576., 1048576.};
+    rotated[i] = Vector3{.6 * kSquare[i].X() - .8 * kSquare[i].Y(),
+                         .8 * kSquare[i].X() + .6 * kSquare[i].Y(), 0};
+  }
+  EXPECT_TRUE(Kernel(translated).HasValue());
+  EXPECT_TRUE(Kernel(rotated).HasValue());
+}
+
+TEST(ThermalShellGeometry,
+     RejectsWholePatchSingularityDespiteNonzeroGaussPoints) {
+  Patch patch{};
+  const std::array<double, 4> xi{{-1, 1, 1, -1}};
+  const std::array<double, 4> eta{{-1, -1, 1, 1}};
+  for (std::size_t i = 0; i < 4; ++i)
+    patch[i] = {xi[i], (xi[i] - .25) * eta[i], 0};
+  const auto result = Kernel(patch);
+  ASSERT_FALSE(result.HasValue());
+  EXPECT_EQ(result.GetStatus().Diagnostics()[0].entity_identity, "part:0027");
+  EXPECT_EQ(result.GetStatus().Diagnostics()[0].location.line, 17U);
+  for (double q : {-std::sqrt(.6), 0., std::sqrt(.6)})
+    EXPECT_GT(std::abs(q - .25), 0.);
+}
+
+TEST(ThermalShellGeometry, RejectsInvalidAndDoesNotRepairUncertainGeometry) {
+  EXPECT_FALSE(
+      Kernel({{{0, 0, 0}, {1, 1, 0}, {1, 0, 0}, {0, 1, 0}}}).HasValue());
+  EXPECT_FALSE(
+      Kernel({{{0, 0, 0}, {1, 0, 0}, {2, 0, 0}, {3, 0, 0}}}).HasValue());
+  Patch repeated = kSquare;
+  repeated[2] = repeated[1];
+  EXPECT_FALSE(Kernel(repeated).HasValue());
+  for (double bad : {std::numeric_limits<double>::infinity(),
+                     std::numeric_limits<double>::quiet_NaN()}) {
+    Patch patch = kSquare;
+    patch[0] = {bad, 0, 0};
+    EXPECT_FALSE(Kernel(patch).HasValue());
+  }
+  const auto duplicate = fesa::ThermalShell::Create(
+      {"part", 1, "1"}, {{0, 1, 1, 3}}, kSquare, 1, 1, 1, 1);
+  EXPECT_FALSE(duplicate.HasValue());
+  const double e = std::ldexp(1., -52);
+  const auto uncertain =
+      Kernel({{{-2, -2 - e, 0}, {0, -e, 0}, {2, 2 + e, 0}, {0, e, 0}}});
+  if (!uncertain.HasValue())
+    EXPECT_EQ(uncertain.GetStatus().Diagnostics()[0].code,
+              "uncertain-thermal-shell-geometry");
+}
+
+TEST(ThermalShellGeometry, DomainEnforcesOppositeEdgesAndOneVertexLink) {
+  const std::vector<Vector3> nodes{{0, 0, 0}, {1, 0, 0}, {1, 1, 0},
+                                   {0, 1, 0}, {2, 0, 0}, {2, 1, 0}};
+  EXPECT_TRUE(Accepted(Mesh(nodes, {{{0, 1, 2, 3}}, {{1, 4, 5, 2}}})));
+  EXPECT_FALSE(Accepted(Mesh(nodes, {{{0, 1, 2, 3}}, {{2, 5, 4, 1}}})));
+  EXPECT_FALSE(Accepted(Mesh(nodes, {{{0, 1, 2, 3}}, {{0, 1, 2, 3}}})));
+  EXPECT_FALSE(Accepted(Mesh({{0, 0, 0},
+                              {1, 0, 0},
+                              {1, 1, 0},
+                              {0, 1, 0},
+                              {-1, -1, 0},
+                              {0, -1, 0},
+                              {-1, 0, 0}},
+                             {{{0, 1, 2, 3}}, {{0, 6, 4, 5}}})));
+  EXPECT_FALSE(
+      Accepted(Mesh({{0, 0, 0},
+                     {1, 0, 0},
+                     {1, 1, 0},
+                     {0, 1, 0},
+                     {2, 0, 0},
+                     {2, 1, 0},
+                     {1, 0, 1},
+                     {1, 1, 1}},
+                    {{{0, 1, 2, 3}}, {{1, 4, 5, 2}}, {{1, 6, 7, 2}}})));
+}
+
+TEST(ThermalShellGeometry, DomainAllowsFoldAndDisconnectedCoordinateContact) {
+  auto folded =
+      Mesh({{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}, {1, 0, 1}, {1, 1, 1}},
+           {{{0, 1, 2, 3}}, {{1, 4, 5, 2}}});
+  folded.thermal_materials.push_back({{"part", 2, "2"}, "heat2", 10, 1, 1, {}});
+  folded.thermal_sections.push_back({{"part", 2, "2"}, "section2", .5, 1, {}});
+  folded.thermal_elements[1].section_index = 1;
+  EXPECT_TRUE(Accepted(std::move(folded)));
+  const std::vector<Vector3> nodes{{0, 0, 0},  {2, 0, 0}, {2, 1, 0},
+                                   {0, 1, 0},  {0, 0, 0}, {1, 0, 0},
+                                   {1, -1, 0}, {0, -1, 0}};
+  auto disconnected = Mesh(nodes, {{{0, 1, 2, 3}}, {{4, 5, 6, 7}}});
+  EXPECT_TRUE(Accepted(disconnected));
+  for (std::size_t i = 4; i < 8; ++i)
+    disconnected.nodes[i].source_id.instance_name = "other";
+  disconnected.thermal_elements[1].source_id.instance_name = "other";
+  EXPECT_TRUE(Accepted(disconnected));
+}
+
+TEST(ThermalShellGeometry,
+     DomainRejectsSourceSharingPartialEdgeWithoutGlobalScan) {
+  // The patches share the full 3->0 / 0->3 edge: connected vertex links remain
+  // paths, while the additional 0->4 edge lies inside coarse 0->1.
+  const std::vector<Vector3> nodes{{0, 0, 0}, {2, 0, 0}, {2, 1, 0},
+                                   {0, 1, 0}, {1, 0, 0}, {1, 1, 1}};
+  const auto model = Mesh(nodes, {{{0, 1, 2, 3}}, {{0, 3, 5, 4}}});
+  const auto result = fesa::Domain::Create(model);
+  ASSERT_FALSE(result.HasValue());
+  EXPECT_EQ(result.GetStatus().Diagnostics()[0].code,
+            "uncertain-thermal-shell-connection");
+  auto invalid =
+      Mesh({{0, 0, 0}, {1, 1, 0}, {1, 0, 0}, {0, 1, 0}}, {{{0, 1, 2, 3}}});
+  EXPECT_FALSE(Accepted(invalid));
+}
+
+TEST(ThermalShellGeometry, SearchesSeparatingDirectionsBeyondCenterAndCorners) {
+  // Exact raw certificate s=(0,-3,6) has min(s.w)=9; center has zero margin.
+  EXPECT_TRUE(
+      Kernel({{{3, -9, -3}, {1, 3, -3}, {-1, 3, 3}, {-3, 3, 3}}}).HasValue());
+  // Exact s=(1.5,2.5,1) gives min(s.w)=3.5, whereas center/corner candidates
+  // cannot separate the normal image. Edge/interior minimization is needed.
+  EXPECT_TRUE(
+      Kernel({{{-1, 4, -4}, {-1, -2, 2}, {3, 2, 4}, {-1, -4, -2}}}).HasValue());
+}
+
+TEST(ThermalShellGeometry,
+     RejectsAliasedSourceIdentityAndAcceptsClosedVertexLink) {
+  auto aliases =
+      Mesh({{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}}, {{{0, 1, 2, 3}}});
+  aliases.nodes[2].source_id = aliases.nodes[1].source_id;
+  EXPECT_FALSE(Accepted(aliases));
+  EXPECT_TRUE(Accepted(
+      Mesh({{-1, -1, 0},
+            {0, -1, 0},
+            {1, -1, 0},
+            {-1, 0, 0},
+            {0, 0, 0},
+            {1, 0, 0},
+            {-1, 1, 0},
+            {0, 1, 0},
+            {1, 1, 0}},
+           {{{0, 1, 4, 3}}, {{1, 2, 5, 4}}, {{3, 4, 7, 6}}, {{4, 5, 8, 7}}})));
+}