diff --git a/include/fesa/analysis/thermal_time_grid.h b/include/fesa/analysis/thermal_time_grid.h
new file mode 100644
index 0000000..8a16f32
--- /dev/null
+++ b/include/fesa/analysis/thermal_time_grid.h
@@ -0,0 +1,11 @@
+#ifndef FESA_ANALYSIS_THERMAL_TIME_GRID_H_
+#define FESA_ANALYSIS_THERMAL_TIME_GRID_H_
+#include <cstdint>
+#include "fesa/core/status.h"
+namespace fesa {
+class ThermalTimeGrid final { public:
+ static Result<ThermalTimeGrid> Create(double period,double increment,std::uint32_t maximum);
+ std::uint32_t Count() const noexcept; double Time(std::uint32_t index) const noexcept; double Delta(std::uint32_t index) const noexcept;
+ private: ThermalTimeGrid(double p,double d,std::uint32_t n,bool close):period_(p),increment_(d),count_(n),close_(close){} double period_,increment_; std::uint32_t count_; bool close_; };
+}
+#endif
diff --git a/include/fesa/elements/thermal_shell.h b/include/fesa/elements/thermal_shell.h
new file mode 100644
index 0000000..3bb7a5f
--- /dev/null
+++ b/include/fesa/elements/thermal_shell.h
@@ -0,0 +1,28 @@
+#ifndef FESA_ELEMENTS_THERMAL_SHELL_H_
+#define FESA_ELEMENTS_THERMAL_SHELL_H_
+
+#include <array>
+#include <cstddef>
+
+#include "fesa/math/matrix.h"
+#include "fesa/math/vector3.h"
+
+namespace fesa {
+
+/// @brief Local Q4 times two-segment P2 transient-heat shell kernel.
+class ThermalShell final {
+ public:
+  ThermalShell(std::array<std::size_t, 4> node_indices,
+               std::array<Vector3, 4> coordinates, double conductivity,
+               double density, double specific_heat, double thickness);
+  const std::array<std::size_t, 4>& NodeIndices() const noexcept;
+  const Matrix& Conductivity() const noexcept;
+  const Matrix& Capacity() const noexcept;
+ private:
+  std::array<std::size_t, 4> node_indices_;
+  Matrix conductivity_{20U, 20U};
+  Matrix capacity_{20U, 20U};
+};
+
+}  // namespace fesa
+#endif  // FESA_ELEMENTS_THERMAL_SHELL_H_
diff --git a/src/fesa/analysis/thermal_time_grid.cpp b/src/fesa/analysis/thermal_time_grid.cpp
new file mode 100644
index 0000000..2f14576
--- /dev/null
+++ b/src/fesa/analysis/thermal_time_grid.cpp
@@ -0,0 +1,4 @@
+#include "fesa/analysis/thermal_time_grid.h"
+#include <cmath>
+#include <limits>
+namespace fesa { Result<ThermalTimeGrid> ThermalTimeGrid::Create(double p,double d,std::uint32_t max) { if(!std::isfinite(p)||!std::isfinite(d)||p<=0||d<=0||max==0) return Result<ThermalTimeGrid>::Failure(Status::Failure(FailureCategory::kModel,{{Severity::kError,"invalid-thermal-time-grid",{},"STEP","","Invalid thermal time grid."}})); if(d>=p) return Result<ThermalTimeGrid>::Success(ThermalTimeGrid{p,d,1,true}); const double q=p/d; const auto m=static_cast<std::uint64_t>(std::nearbyint(q)); const bool close=m>0 && std::abs(std::fma(-static_cast<double>(m),d,p)) <= .5*(std::nextafter(p,std::numeric_limits<double>::infinity())-p)+.5*static_cast<double>(m)*(std::nextafter(d,std::numeric_limits<double>::infinity())-d); const std::uint64_t n=close?m:static_cast<std::uint64_t>(std::ceil(q)); if(n>max||n>std::numeric_limits<std::uint32_t>::max()) return Result<ThermalTimeGrid>::Failure(Status::Failure(FailureCategory::kModel,{{Severity::kError,"thermal-increment-limit",{},"STEP","","Thermal increment count exceeds INC."}})); return Result<ThermalTimeGrid>::Success(ThermalTimeGrid{p,d,static_cast<std::uint32_t>(n),close}); } std::uint32_t ThermalTimeGrid::Count()const noexcept{return count_;} double ThermalTimeGrid::Time(std::uint32_t i)const noexcept{return i==count_?period_:static_cast<double>(i)*increment_;} double ThermalTimeGrid::Delta(std::uint32_t i)const noexcept{return (i==count_&&!close_)?std::fma(-static_cast<double>(count_-1),increment_,period_):increment_;} }
diff --git a/src/fesa/elements/thermal_shell.cpp b/src/fesa/elements/thermal_shell.cpp
new file mode 100644
index 0000000..417279f
--- /dev/null
+++ b/src/fesa/elements/thermal_shell.cpp
@@ -0,0 +1,43 @@
+#include "fesa/elements/thermal_shell.h"
+
+#include <array>
+#include <cmath>
+#include <stdexcept>
+
+namespace fesa {
+namespace {
+constexpr std::array<double, 3> kGauss{{-0.7745966692414834, 0.0, 0.7745966692414834}};
+constexpr std::array<double, 3> kWeight{{5.0 / 9.0, 8.0 / 9.0, 5.0 / 9.0}};
+constexpr std::array<double, 4> kXi{{-1., 1., 1., -1.}};
+constexpr std::array<double, 4> 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.; }
+}
+}  // namespace
+
+ThermalShell::ThermalShell(std::array<std::size_t, 4> node_indices,
+                           std::array<Vector3, 4> 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<double,4> 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; }
+  }
+  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)); }
+}
+const std::array<std::size_t,4>& ThermalShell::NodeIndices() const noexcept { return node_indices_; }
+const Matrix& ThermalShell::Conductivity() const noexcept { return conductivity_; }
+const Matrix& ThermalShell::Capacity() const noexcept { return capacity_; }
+}  // namespace fesa
diff --git a/src/fesa/elements/thermal_shell_geometry.cpp b/src/fesa/elements/thermal_shell_geometry.cpp
new file mode 100644
index 0000000..053bbf3
--- /dev/null
+++ b/src/fesa/elements/thermal_shell_geometry.cpp
@@ -0,0 +1,3 @@
+#include "thermal_shell_geometry.h"
+#include <cmath>
+namespace fesa { Status ValidateThermalShellGeometry(const std::array<Vector3,4>& 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(); } }
diff --git a/src/fesa/elements/thermal_shell_geometry.h b/src/fesa/elements/thermal_shell_geometry.h
new file mode 100644
index 0000000..b177b51
--- /dev/null
+++ b/src/fesa/elements/thermal_shell_geometry.h
@@ -0,0 +1,10 @@
+#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"
+namespace fesa {
+/// @brief Validates the finite, oriented Q4 geometry used by the thermal kernel.
+Status ValidateThermalShellGeometry(const std::array<Vector3,4>& coordinates);
+}
+#endif