diff --git a/include/fesa/analysis/thermal_time_grid.h b/include/fesa/analysis/thermal_time_grid.h index 6a6ac90..d27a454 100644 --- a/include/fesa/analysis/thermal_time_grid.h +++ b/include/fesa/analysis/thermal_time_grid.h @@ -13,21 +13,25 @@ struct ThermalTimeFrame { double delta; }; -/// @brief Stores an O(1) fixed-increment grid with exact input-cell endpoint tests. +/// @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. + /// @brief Checks finite positive controls and the uint32 INC bound before + /// use. static Result Create(double period, double increment, - std::uint32_t maximum = 100U); + 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. + /// @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). + /// @brief Checks index, finite values and strict increase for one frame in + /// O(1). Result Frame(std::uint32_t index) const; private: diff --git a/include/fesa/elements/thermal_shell.h b/include/fesa/elements/thermal_shell.h index 3bb7a5f..8bd1197 100644 --- a/include/fesa/elements/thermal_shell.h +++ b/include/fesa/elements/thermal_shell.h @@ -2,26 +2,71 @@ #define FESA_ELEMENTS_THERMAL_SHELL_H_ #include -#include +#include "fesa/core/source_identity.h" +#include "fesa/core/status.h" #include "fesa/math/matrix.h" #include "fesa/math/vector3.h" namespace fesa { -/// @brief Local Q4 times two-segment P2 transient-heat shell kernel. +/// @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 global_flux; + std::array, 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: - ThermalShell(std::array node_indices, - std::array coordinates, double conductivity, - double density, double specific_heat, double thickness); - const std::array& NodeIndices() const noexcept; + /// @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 Create(SourceEntityId source, + std::array node_indices, + const std::array& 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& 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 Recover( + const std::array& temperature) const; + private: - std::array node_indices_; + /// @brief Initializes an unpublished candidate whose operators are filled by + /// Create. + ThermalShell(SourceEntityId source, std::array node_indices, + SourceLocation location, double conductivity); + SourceEntityId source_; + std::array 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, 24> gradients_{}; + std::array, 4> frames_{}; }; } // namespace fesa diff --git a/src/fesa/analysis/thermal_time_grid.cpp b/src/fesa/analysis/thermal_time_grid.cpp index a7499ab..16a3e9a 100644 --- a/src/fesa/analysis/thermal_time_grid.cpp +++ b/src/fesa/analysis/thermal_time_grid.cpp @@ -38,7 +38,8 @@ 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::uint64_t{a.digit[i]} - sub); + result.digit[i] = + static_cast(std::uint64_t{a.digit[i]} - sub); borrow = std::uint64_t{a.digit[i]} < sub ? 1U : 0U; } return borrow == 0; @@ -68,8 +69,8 @@ bool Multiply(const UInt128& a, std::uint64_t multiplier, UInt128& result) { 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; + const std::uint64_t product = + std::uint64_t{a.digit[i]} * factor + result.digit[i + j] + carry; result.digit[i + j] = static_cast(product); carry = product >> 32; } @@ -85,7 +86,8 @@ struct Cell { bool closed; }; -/// @brief Constructs ORIGINAL binary64 cells, including subnormal/maxfinite gaps. +/// @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)); @@ -108,21 +110,23 @@ bool MakeCell(double value, int period_exponent, Cell& cell) { Add(cell.center, upper_gap, cell.upper); } -/// @brief Returns source-neutral utility failure for parser/procedure decoration. +/// @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}}); + return Status::Failure( + FailureCategory::kModel, + {{Severity::kError, code, {}, "HEAT TRANSFER", "", message}}); } } // namespace -Result ThermalTimeGrid::Create( - double period, double increment, std::uint32_t maximum) { +Result 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::Failure(GridFailure( - "invalid-thermal-time-grid", "Time controls must be finite and positive.")); + return Result::Failure( + GridFailure("invalid-thermal-time-grid", + "Time controls must be finite and positive.")); } if (increment >= period) { return Result::Success( @@ -152,8 +156,7 @@ Result ThermalTimeGrid::Create( 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)) && + const bool close = nearest != 0 && (first > 0 || (first == 0 && closed)) && (second > 0 || (second == 0 && closed)); std::uint64_t count = nearest; if (!close && exact_ok) { @@ -180,11 +183,14 @@ Result ThermalTimeGrid::Create( return Result::Failure(GridFailure( "thermal-increment-limit", "Thermal increment count exceeds INC.")); } - const double last = close ? increment : - std::fma(-static_cast(count - 1), increment, period); - ThermalTimeGrid grid{period, increment, static_cast(count), last}; + const double last = + close ? increment + : std::fma(-static_cast(count - 1), increment, period); + ThermalTimeGrid grid{period, increment, static_cast(count), + last}; const auto frame = grid.Frame(grid.Count()); - if (!frame.HasValue()) return Result::Failure(frame.GetStatus()); + if (!frame.HasValue()) + return Result::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. @@ -192,8 +198,10 @@ Result ThermalTimeGrid::Create( } ThermalTimeGrid::ThermalTimeGrid(double period, double increment, - std::uint32_t count, double final_delta) - : period_{period}, increment_{increment}, count_{count}, + 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_; } @@ -205,16 +213,18 @@ double ThermalTimeGrid::Delta(std::uint32_t index) const noexcept { } Result ThermalTimeGrid::Frame(std::uint32_t index) const { if (index == 0 || index > count_) { - return Result::Failure(GridFailure( - "invalid-thermal-time-frame", "Completed frame index is outside the grid.")); + return Result::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::Failure(GridFailure( - "invalid-thermal-time-frame", "Stored times must be finite and strictly increasing.")); + if (!std::isfinite(time) || !std::isfinite(previous) || + !std::isfinite(delta) || !(time > previous) || !(delta > 0)) { + return Result::Failure( + GridFailure("invalid-thermal-time-frame", + "Stored times must be finite and strictly increasing.")); } return Result::Success({time, delta}); } diff --git a/src/fesa/elements/thermal_interval.h b/src/fesa/elements/thermal_interval.h new file mode 100644 index 0000000..9b7c836 --- /dev/null +++ b/src/fesa/elements/thermal_interval.h @@ -0,0 +1,108 @@ +#ifndef FESA_ELEMENTS_THERMAL_INTERVAL_H_ +#define FESA_ELEMENTS_THERMAL_INTERVAL_H_ + +#include +#include +#include +#include + +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::infinity()); +} +/// @brief Rounds an operation endpoint toward positive infinity. +inline double Up(double value) { + return std::nextafter(value, std::numeric_limits::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 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::infinity(), + std::numeric_limits::infinity()}; + const std::array 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; +/// @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_ diff --git a/src/fesa/elements/thermal_shell.cpp b/src/fesa/elements/thermal_shell.cpp index 417279f..35ccb7d 100644 --- a/src/fesa/elements/thermal_shell.cpp +++ b/src/fesa/elements/thermal_shell.cpp @@ -1,43 +1,236 @@ #include "fesa/elements/thermal_shell.h" +#include #include #include -#include +#include + +#include "thermal_shell_geometry.h" namespace fesa { namespace { -constexpr std::array kGauss{{-0.7745966692414834, 0.0, 0.7745966692414834}}; -constexpr std::array kWeight{{5.0 / 9.0, 8.0 / 9.0, 5.0 / 9.0}}; -constexpr std::array kXi{{-1., 1., 1., -1.}}; -constexpr std::array kEta{{-1., -1., 1., 1.}}; -void AddZ(Matrix& z0, Matrix& z1, std::size_t first) { - const double a[3][3]={{4,2,-1},{2,16,2},{-1,2,4}}; - const double b[3][3]={{7,-8,1},{-8,16,-8},{1,-8,7}}; - for (std::size_t i=0;i<3;++i) for (std::size_t j=0;j<3;++j) { z0(first+i,first+j)+=a[i][j]/30.; z1(first+i,first+j)+=b[i][j]/3.; } +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(std::array node_indices, - std::array coordinates, - double conductivity, double density, - double specific_heat, double thickness) - : node_indices_{node_indices} { - if (!(conductivity > 0.) || !(density > 0.) || !(specific_heat > 0.) || - !(thickness > 0.)) throw std::invalid_argument("invalid thermal shell material"); - Matrix m{4U,4U}; Matrix s{4U,4U}; Matrix z0{5U,5U}; Matrix z1{5U,5U}; - AddZ(z0,z1,0U); AddZ(z0,z1,2U); - for (std::size_t x=0;x<3;++x) for (std::size_t y=0;y<3;++y) { - const double xi=kGauss[x], eta=kGauss[y]; std::array n{}, dx{}, dy{}; - Vector3 a1{},a2{}; - for(std::size_t i=0;i<4;++i) { n[i]=.25*(1+kXi[i]*xi)*(1+kEta[i]*eta); dx[i]=.25*kXi[i]*(1+kEta[i]*eta); dy[i]=.25*kEta[i]*(1+kXi[i]*xi); a1=a1+dx[i]*coordinates[i]; a2=a2+dy[i]*coordinates[i]; } - const Vector3 w=a1.Cross(a2); const double jac=w.Norm(); const double g11=a1.Dot(a1),g12=a1.Dot(a2),g22=a2.Dot(a2), det=g11*g22-g12*g12; - if (!(jac>0.) || !(det>0.) || !std::isfinite(jac)) throw std::invalid_argument("invalid thermal shell geometry"); - const double wt=kWeight[x]*kWeight[y]*jac; - for(std::size_t i=0;i<4;++i) for(std::size_t j=0;j<4;++j) { m(i,j)+=n[i]*n[j]*wt; s(i,j)+=(dx[i]*(g22*dx[j]-g12*dy[j])+dy[i]*(g11*dy[j]-g12*dx[j]))/det*wt; } +ThermalShell::ThermalShell(SourceEntityId source, + std::array node_indices, + SourceLocation location, double conductivity) + : source_{std::move(source)}, + node_indices_{node_indices}, + location_{std::move(location)}, + conductivity_value_{conductivity} {} + +Result ThermalShell::Create( + SourceEntityId source, std::array node_indices, + const std::array& 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::Failure( + GeometryFailure(source, location, "invalid-thermal-shell-material", + "Thermal properties must be finite and positive.")); + } } - 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,q=5*j+b; capacity_(p,q)=density*specific_heat*thickness*.5*m(i,j)*z0(a,b); conductivity_(p,q)=conductivity*(thickness*.5*s(i,j)*z0(a,b)+2./thickness*m(i,j)*z1(a,b)); } + auto sorted = node_indices; + std::sort(sorted.begin(), sorted.end()); + if (std::adjacent_find(sorted.begin(), sorted.end()) != sorted.end()) { + return Result::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::Failure(geometry.GetStatus()); + ThermalShell kernel{source, node_indices, location, conductivity}; + Matrix mass{4, 4}; + Matrix surface{4, 4}; + std::array, 4> mass_interval{}; + std::array, 4> surface_interval{}; + const double gauss = std::sqrt(3. / 5.); + const Interval gauss_interval = Sqrt(Interval{3} / Interval{5}); + const std::array points{{-gauss, 0., gauss}}; + const std::array point_intervals{ + {Interval{0} - gauss_interval, Interval{0}, gauss_interval}}; + const std::array weights{{5. / 9., 8. / 9., 5. / 9.}}; + const std::array 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::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, 5> z0_interval{}; + std::array, 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::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 kXi{{-1, 1, 1, -1}}; + constexpr std::array kEta{{-1, -1, 1, 1}}; + constexpr std::array 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::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 shape{ + {.5 * s * (s - 1), 1 - s * s, .5 * s * (s + 1)}}; + const std::array 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::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::Success(std::move(kernel)); +} + +const SourceEntityId& ThermalShell::SourceId() const noexcept { + return source_; +} +const std::array& ThermalShell::NodeIndices() const noexcept { + return node_indices_; +} +const Matrix& ThermalShell::Conductivity() const noexcept { + return conductivity_; } -const std::array& 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 ThermalShell::Recover( + const std::array& temperature) const { + for (const double value : temperature) { + if (!std::isfinite(value)) { + return Result::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::Failure( + thermal_internal::GeometryFailure( + source_, location_, "nonfinite-thermal-shell-flux", + "Recovered heat flux is not finite.")); + } + } + return Result::Success(recovery); +} + } // namespace fesa diff --git a/src/fesa/elements/thermal_shell_geometry.cpp b/src/fesa/elements/thermal_shell_geometry.cpp index 053bbf3..8219714 100644 --- a/src/fesa/elements/thermal_shell_geometry.cpp +++ b/src/fesa/elements/thermal_shell_geometry.cpp @@ -1,3 +1,431 @@ #include "thermal_shell_geometry.h" + +#include +#include #include -namespace fesa { Status ValidateThermalShellGeometry(const std::array& c) { for(const auto& p:c) if(!p.IsFinite()) return Status::Failure(FailureCategory::kModel,{{Severity::kError,"invalid-thermal-shell-geometry",{},"DS4","","Coordinates are nonfinite."}}); const Vector3 b=.25*(-1.*c[0]+c[1]+c[2]-c[3]); const Vector3 d=.25*(c[0]-c[1]+c[2]-c[3]); const Vector3 e=.25*(-1.*c[0]-c[1]+c[2]+c[3]); Vector3 candidate=b.Cross(e); for(double xi:{-1.,1.}) for(double eta:{-1.,1.}) { const Vector3 w=(b+eta*d).Cross(e+xi*d); if(!(candidate.Dot(w)>0.)||!std::isfinite(candidate.Dot(w))) return Status::Failure(FailureCategory::kModel,{{Severity::kError,"invalid-thermal-shell-geometry",{},"DS4","","Q4 patch has no positive orientation certificate."}}); } return Status::Ok(); } } +#include +#include +#include +#include +#include +#include +#include + +#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::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 BuildGeometry(const std::array& 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::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::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 offsets{}; + std::array 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 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 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 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::Success(geometry); + } + } + return Result::Failure( + GeometryFailure(source, location, "uncertain-thermal-shell-geometry", + "No strictly positive whole-patch interval bound was " + "certified (lower <= 0).")); +} + +Result 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::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 kXi{{-1., 1., 1., -1.}}; + constexpr std::array 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 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::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::Failure(GeometryFailure( + source, location, "uncertain-thermal-shell-metric", + "Surface measure or frame is not finite and positive.")); + } + return Result::Success(point); +} +} // namespace thermal_internal + +namespace { +using Edge = std::pair; + +/// @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> edges; + std::map>> links; + std::set> patches; + std::map, 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 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 visited; + std::vector 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> 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 diff --git a/src/fesa/elements/thermal_shell_geometry.h b/src/fesa/elements/thermal_shell_geometry.h index b177b51..a31b339 100644 --- a/src/fesa/elements/thermal_shell_geometry.h +++ b/src/fesa/elements/thermal_shell_geometry.h @@ -1,10 +1,59 @@ #ifndef FESA_ELEMENTS_THERMAL_SHELL_GEOMETRY_H_ #define FESA_ELEMENTS_THERMAL_SHELL_GEOMETRY_H_ + #include + #include "fesa/core/status.h" #include "fesa/math/vector3.h" +#include "thermal_interval.h" + namespace fesa { -/// @brief Validates the finite, oriented Q4 geometry used by the thermal kernel. -Status ValidateThermalShellGeometry(const std::array& coordinates); -} -#endif +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 shape{}; + std::array shape_interval{}; + std::array gradient{}; + std::array gradient_interval{}; + std::array axes{}; + double jacobian = 0; + Interval jacobian_interval; +}; + +/// @brief Certifies whole-patch rank/injectivity from fixed candidate +/// directions. +Result BuildGeometry(const std::array& coordinates, + const SourceEntityId& source, + const SourceLocation& location); +/// @brief Checks finite metric inverse and gradients at one used surface point. +Result 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_ diff --git a/src/fesa/model/domain.cpp b/src/fesa/model/domain.cpp index 1d30af1..5989890 100644 --- a/src/fesa/model/domain.cpp +++ b/src/fesa/model/domain.cpp @@ -8,6 +8,8 @@ #include #include +#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> initial_values; for (const auto& initial : definition.thermal_initial_conditions) { if (initial.node_index >= active_nodes.size() || diff --git a/tests/unit/analysis/thermal_time_grid_test.cpp b/tests/unit/analysis/thermal_time_grid_test.cpp index 59448ee..3014445 100644 --- a/tests/unit/analysis/thermal_time_grid_test.cpp +++ b/tests/unit/analysis/thermal_time_grid_test.cpp @@ -13,8 +13,8 @@ TEST(ThermalTimeGrid, SingleStepUsesPeriodAsActualIncrement) { } TEST(ThermalTimeGrid, PreservesGenuineRemainderOutsideInputCells) { - const auto grid = fesa::ThermalTimeGrid::Create( - std::nextafter(.9, 1.0), .3, 4U); + 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); @@ -41,13 +41,14 @@ TEST(ThermalTimeGrid, DecimalCellsCounterexampleAndScale) { 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); + 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()); + EXPECT_FALSE( + fesa::ThermalTimeGrid::Create(5.351455173733616, .17262758624947144, 30U) + .HasValue()); } TEST(ThermalTimeGrid, BinadeLimitsAndMaximumFinitePreserveOriginalCells) { @@ -56,13 +57,16 @@ TEST(ThermalTimeGrid, BinadeLimitsAndMaximumFinitePreserveOriginalCells) { 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); + 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::denorm_min(), 4294967295U).HasValue()); + EXPECT_FALSE( + fesa::ThermalTimeGrid::Create(4294967296., 1., 4294967295U).HasValue()); + EXPECT_FALSE(fesa::ThermalTimeGrid::Create( + max, std::numeric_limits::denorm_min(), 4294967295U) + .HasValue()); } TEST(ThermalTimeGrid, ShortenedIncrementAndCheckedFrames) { @@ -83,7 +87,8 @@ TEST(ThermalTimeGrid, ShortenedIncrementAndCheckedFrames) { 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); + 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); diff --git a/tests/unit/elements/thermal_shell_test.cpp b/tests/unit/elements/thermal_shell_test.cpp new file mode 100644 index 0000000..fc68ad4 --- /dev/null +++ b/tests/unit/elements/thermal_shell_test.cpp @@ -0,0 +1,498 @@ +#include "fesa/elements/thermal_shell.h" + +#include + +#include +#include +#include +#include +#include + +namespace { +using fesa::Matrix; +using fesa::Vector3; +using Patch = std::array; +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 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 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 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 l{ + {s * (s - 1) / 2, 1 - s * s, s * (s + 1) / 2}}; + const std::array dl{{2 * s - 1, -4 * s, 2 * s + 1}}; + std::array shape{}; + std::array 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 value{}; + std::array scale{}; +}; +FluxOracle Flux(const Patch& patch, double h, const std::array& 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 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(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 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(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(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{{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((j + 1) * (j + 2))); + for (std::size_t i = 0; i <= j; ++i) basis(i, j) = 1 / norm; + basis(j + 1, j) = -static_cast(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 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 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(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 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::infinity(), + std::numeric_limits::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 temperature{}; + temperature[3] = std::numeric_limits::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::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 temperature{}; + for (std::size_t p = 0; p < 20; ++p) + temperature[p] = static_cast((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::denorm_min(), + std::numeric_limits::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::max(), 2, 1); + EXPECT_FALSE(overflow.HasValue()); +} diff --git a/tests/unit/model/thermal_shell_geometry_test.cpp b/tests/unit/model/thermal_shell_geometry_test.cpp new file mode 100644 index 0000000..4c7633d --- /dev/null +++ b/tests/unit/model/thermal_shell_geometry_test.cpp @@ -0,0 +1,212 @@ +#include + +#include +#include +#include +#include + +#include "fesa/elements/thermal_shell.h" +#include "fesa/model/domain.h" + +namespace { +using fesa::Vector3; +using Patch = std::array; +const Patch kSquare{{{0, 0, 0}, {1, 0, 0}, {1, 1, 0}, {0, 1, 0}}}; + +fesa::Result 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& nodes, + const std::vector>& elements) { + fesa::ModelDefinition model{}; + for (std::size_t i = 0; i < nodes.size(); ++i) { + model.nodes.push_back( + {{"part", static_cast(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(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 xi{{-1, 1, 1, -1}}; + const std::array 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::infinity(), + std::numeric_limits::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 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 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 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}}}))); +}