diff --git a/include/fesa/analysis/analysis_model.h b/include/fesa/analysis/analysis_model.h
index d3fddde..1068a3a 100644
--- a/include/fesa/analysis/analysis_model.h
+++ b/include/fesa/analysis/analysis_model.h
@@ -20,8 +20,16 @@ class AnalysisModel {
   /// @brief Returns the non-owned Domain backing this view.
   const Domain& GetDomain() const noexcept;
 
-  /// @brief Returns the sole active static step.
-  const StepDefinition& Step() const noexcept;
+  /// @brief Returns the sole active structural step.
+  /// @throws std::logic_error when this is a thermal view.
+  /// @note Thermal callers must use ThermalStep().
+  const StepDefinition& Step() const;
+
+  /// @brief Reports whether this view selects the uncoupled thermal step.
+  bool IsThermal() const noexcept;
+
+  /// @brief Returns the selected thermal step, or null for structural views.
+  const ThermalStepDefinition* ThermalStep() const noexcept;
 
   /// @brief Returns active element-definition indices in stable Domain order.
   const std::vector<EntityIndex>& ActiveElements() const noexcept;
@@ -44,6 +52,18 @@ class AnalysisModel {
   /// @brief Returns concentrated-load indices in source order.
   const std::vector<EntityIndex>& ActiveLoads() const noexcept;
 
+  /// @brief Returns connected DS4 element indices in stable Domain order.
+  const std::vector<EntityIndex>& ActiveThermalElements() const noexcept;
+
+  /// @brief Returns connected thermal node indices in stable Domain order.
+  const std::vector<EntityIndex>& ActiveThermalNodes() const noexcept;
+
+  /// @brief Returns thermal-material collection positions in source order.
+  const std::vector<EntityIndex>& ActiveThermalMaterials() const noexcept;
+
+  /// @brief Returns thermal-section collection positions in source order.
+  const std::vector<EntityIndex>& ActiveThermalSections() const noexcept;
+
   /// @brief Returns active surface pressure records without copying Domain data.
   const std::vector<SurfacePressureLoadDefinition>& ActiveSurfacePressureLoads()
       const noexcept;
@@ -60,6 +80,11 @@ class AnalysisModel {
   std::vector<EntityIndex> active_sections_;
   std::vector<EntityIndex> active_boundary_conditions_;
   std::vector<EntityIndex> active_loads_;
+  std::vector<EntityIndex> active_thermal_elements_;
+  std::vector<EntityIndex> active_thermal_nodes_;
+  std::vector<EntityIndex> active_thermal_materials_;
+  std::vector<EntityIndex> active_thermal_sections_;
+  bool is_thermal_{false};
 };
 
 }  // namespace fesa
diff --git a/include/fesa/elements/ds4_thermal.h b/include/fesa/elements/ds4_thermal.h
new file mode 100644
index 0000000..1c6f0af
--- /dev/null
+++ b/include/fesa/elements/ds4_thermal.h
@@ -0,0 +1,99 @@
+#ifndef FESA_ELEMENTS_DS4_THERMAL_H_
+#define FESA_ELEMENTS_DS4_THERMAL_H_
+
+#include <array>
+#include <limits>
+#include <utility>
+#include <vector>
+
+#include "fesa/elements/element_definition.h"
+
+namespace fesa {
+
+/// @brief Defines one four-node uncoupled transient thermal DS4 element.
+/// @note This semantic definition contains no geometry kernel or equation
+/// state.
+class Ds4ThermalDefinition final : public ElementDefinition {
+ public:
+  /// @brief Marks that DS4 has no entry in Domain::Properties().
+  static constexpr EntityIndex kNoPropertyIndex =
+      (std::numeric_limits<EntityIndex>::max)();
+  Ds4ThermalDefinition(SourceEntityId source_id,
+                       std::array<EntityIndex, 4> node_indices,
+                       EntityIndex section_index, SourceLocation location);
+
+  ElementDefinitionKind Kind() const noexcept override;
+  const SourceEntityId& SourceId() const noexcept override;
+  std::string_view SourceElementType() const noexcept override;
+  const std::vector<EntityIndex>& NodeIndices() const noexcept override;
+  EntityIndex MaterialIndex() const noexcept override;
+  EntityIndex PropertyIndex() const noexcept override;
+
+  /// @brief Returns the stable thermal-section collection position.
+  EntityIndex ThermalSectionIndex() const noexcept;
+  /// @brief Returns the source location defining this element.
+  const SourceLocation& Location() const noexcept;
+
+ private:
+  friend class Domain;
+
+  void BindMaterialIndex(EntityIndex material_index) noexcept;
+
+  SourceEntityId source_id_;
+  std::array<EntityIndex, 4> node_indices_;
+  std::vector<EntityIndex> node_indices_view_;
+  EntityIndex section_index_;
+  EntityIndex material_index_{0U};
+  SourceLocation location_;
+};
+
+inline Ds4ThermalDefinition::Ds4ThermalDefinition(
+    SourceEntityId source_id, std::array<EntityIndex, 4> node_indices,
+    EntityIndex section_index, SourceLocation location)
+    : source_id_{std::move(source_id)},
+      node_indices_{node_indices},
+      node_indices_view_{node_indices.begin(), node_indices.end()},
+      section_index_{section_index},
+      location_{std::move(location)} {}
+
+inline ElementDefinitionKind Ds4ThermalDefinition::Kind() const noexcept {
+  return ElementDefinitionKind::kDs4Thermal;
+}
+
+inline const SourceEntityId& Ds4ThermalDefinition::SourceId() const noexcept {
+  return source_id_;
+}
+
+inline std::string_view Ds4ThermalDefinition::SourceElementType() const noexcept {
+  return "DS4";
+}
+
+inline const std::vector<EntityIndex>& Ds4ThermalDefinition::NodeIndices()
+    const noexcept {
+  return node_indices_view_;
+}
+
+inline EntityIndex Ds4ThermalDefinition::MaterialIndex() const noexcept {
+  return material_index_;
+}
+
+inline EntityIndex Ds4ThermalDefinition::PropertyIndex() const noexcept {
+  return kNoPropertyIndex;
+}
+
+inline EntityIndex Ds4ThermalDefinition::ThermalSectionIndex() const noexcept {
+  return section_index_;
+}
+
+inline const SourceLocation& Ds4ThermalDefinition::Location() const noexcept {
+  return location_;
+}
+
+inline void Ds4ThermalDefinition::BindMaterialIndex(
+    EntityIndex material_index) noexcept {
+  material_index_ = material_index;
+}
+
+}  // namespace fesa
+
+#endif  // FESA_ELEMENTS_DS4_THERMAL_H_
diff --git a/include/fesa/elements/element_definition.h b/include/fesa/elements/element_definition.h
index b8baca6..ed832a2 100644
--- a/include/fesa/elements/element_definition.h
+++ b/include/fesa/elements/element_definition.h
@@ -9,7 +9,7 @@
 namespace fesa {
 
 /// @brief Identifies the supported semantic element-definition kinds.
-enum class ElementDefinitionKind { kEulerBeam3D, kMitc4Shell };
+enum class ElementDefinitionKind { kEulerBeam3D, kMitc4Shell, kDs4Thermal };
 
 /// @brief Provides immutable source identity and topology for one element.
 /// @note Numerical stiffness, recovery, and equation ids are intentionally
diff --git a/include/fesa/materials/isotropic_thermal_material.h b/include/fesa/materials/isotropic_thermal_material.h
new file mode 100644
index 0000000..561af16
--- /dev/null
+++ b/include/fesa/materials/isotropic_thermal_material.h
@@ -0,0 +1,48 @@
+#ifndef FESA_MATERIALS_ISOTROPIC_THERMAL_MATERIAL_H_
+#define FESA_MATERIALS_ISOTROPIC_THERMAL_MATERIAL_H_
+
+#include <string>
+
+#include "fesa/core/status.h"
+#include "fesa/materials/material.h"
+
+namespace fesa {
+
+/// @brief Stores homogeneous isotropic heat-conduction material properties.
+class IsotropicThermalMaterial final : public Material {
+ public:
+  /// @brief Creates a material with finite positive conductivity, density, and
+  /// specific heat.
+  static Result<IsotropicThermalMaterial> Create(
+      SourceEntityId source_id, std::string name, double conductivity,
+      double density, double specific_heat, SourceLocation location);
+
+  MaterialKind Kind() const noexcept override;
+  const SourceEntityId& SourceId() const noexcept override;
+  const SourceLocation& Location() const noexcept override;
+
+  /// @brief Returns the source material name.
+  const std::string& Name() const noexcept;
+  /// @brief Returns conductivity in the active consistent unit system.
+  double Conductivity() const noexcept;
+  /// @brief Returns mass density in the active consistent unit system.
+  double Density() const noexcept;
+  /// @brief Returns specific heat in the active consistent unit system.
+  double SpecificHeat() const noexcept;
+
+ private:
+  IsotropicThermalMaterial(SourceEntityId source_id, std::string name,
+                           double conductivity, double density,
+                           double specific_heat, SourceLocation location);
+
+  SourceEntityId source_id_;
+  std::string name_;
+  double conductivity_;
+  double density_;
+  double specific_heat_;
+  SourceLocation location_;
+};
+
+}  // namespace fesa
+
+#endif  // FESA_MATERIALS_ISOTROPIC_THERMAL_MATERIAL_H_
diff --git a/include/fesa/materials/material.h b/include/fesa/materials/material.h
index 5643b3d..c93b46c 100644
--- a/include/fesa/materials/material.h
+++ b/include/fesa/materials/material.h
@@ -6,7 +6,7 @@
 namespace fesa {
 
 /// @brief Identifies the supported concrete material semantics.
-enum class MaterialKind { kIsotropicLinearElastic };
+enum class MaterialKind { kIsotropicLinearElastic, kIsotropicThermal };
 
 /// @brief Provides stable identity for a Domain-owned material definition.
 class Material {
diff --git a/include/fesa/model/domain.h b/include/fesa/model/domain.h
index c63f8cd..b7afd9b 100644
--- a/include/fesa/model/domain.h
+++ b/include/fesa/model/domain.h
@@ -183,6 +183,30 @@ class Domain {
   const DomainCollectionView<LinearElasticMaterial>& LinearElasticMaterials()
       const noexcept;
 
+  /// @brief Returns thermal materials in stable thermal declaration order.
+  const DomainCollectionView<IsotropicThermalMaterial>& ThermalMaterials()
+      const noexcept;
+
+  /// @brief Returns DS4 definitions in stable thermal declaration order.
+  const DomainCollectionView<Ds4ThermalDefinition>& ThermalElements()
+      const noexcept;
+
+  /// @brief Returns centered thermal shell assignments in declaration order.
+  const std::vector<ThermalShellSectionDefinition>& ThermalSections()
+      const noexcept;
+
+  /// @brief Returns resolved initial temperatures in source order.
+  const std::vector<ThermalInitialConditionDefinition>&
+  ThermalInitialConditions() const noexcept;
+
+  /// @brief Returns model-level thermal boundary definitions in source order.
+  const std::vector<ThermalBoundaryConditionDefinition>&
+  ThermalModelBoundaries() const noexcept;
+
+  /// @brief Returns the sole thermal step, or null when the model is
+  /// structural.
+  const ThermalStepDefinition* ThermalStep() const noexcept;
+
   /// @brief Returns all properties in stable Domain index order.
   const DomainCollectionView<ElementProperty>& Properties() const noexcept;
 
@@ -243,6 +267,8 @@ class Domain {
   DomainCollectionView<ShellSection> shell_sections_view_;
   DomainCollectionView<Material> materials_view_;
   DomainCollectionView<LinearElasticMaterial> linear_materials_view_;
+  DomainCollectionView<IsotropicThermalMaterial> thermal_materials_view_;
+  DomainCollectionView<Ds4ThermalDefinition> thermal_elements_view_;
   DomainCollectionView<StepDefinition> steps_view_;
 };
 
diff --git a/include/fesa/model/model_types.h b/include/fesa/model/model_types.h
index fdfd184..fd5845e 100644
--- a/include/fesa/model/model_types.h
+++ b/include/fesa/model/model_types.h
@@ -12,9 +12,11 @@
 
 #include "fesa/core/diagnostic.h"
 #include "fesa/core/source_identity.h"
+#include "fesa/elements/ds4_thermal.h"
 #include "fesa/elements/euler_beam_3d.h"
 #include "fesa/elements/mitc4_shell.h"
 #include "fesa/materials/isotropic_linear_elastic_material.h"
+#include "fesa/materials/isotropic_thermal_material.h"
 #include "fesa/properties/general_beam_section.h"
 #include "fesa/properties/shell_section.h"
 
@@ -226,6 +228,101 @@ struct DynamicStepDefinition {
   std::vector<NodalLoad> loads;
 };
 
+/// @brief Stores homogeneous isotropic thermal material source data.
+struct ThermalMaterialDefinition {
+  SourceEntityId source_id;
+  std::string name;
+  double conductivity;
+  double density;
+  double specific_heat;
+  SourceLocation location;
+};
+
+/// @brief Stores one centered homogeneous thermal shell assignment.
+/// @note material_index addresses ModelDefinition::thermal_materials and,
+/// after Domain construction, Domain::ThermalMaterials(), never Properties().
+struct ThermalShellSectionDefinition {
+  SourceEntityId source_id;
+  std::string name;
+  double thickness;
+  EntityIndex material_index;
+  SourceLocation location;
+};
+
+/// @brief Stores one DS4 definition before Domain ownership.
+struct Ds4ThermalElementDefinition {
+  SourceEntityId source_id;
+  std::array<EntityIndex, 4> node_indices;
+  EntityIndex section_index;
+  SourceLocation location;
+};
+
+/// @brief Stores one resolved initial temperature target at five stations.
+struct ThermalInitialConditionDefinition {
+  ThermalInitialConditionDefinition(EntityIndex target_node, double value,
+                                    SourceLocation source_location)
+      : node_index{target_node},
+        station_values{{value, value, value, value, value}},
+        location{std::move(source_location)} {}
+
+  EntityIndex node_index;
+  std::array<double, 5> station_values;
+  SourceLocation location;
+};
+
+/// @brief Identifies the fixed bottom-to-top thermal station source DOFs.
+/// @note This is separate from the six-component structural DofComponent.
+enum class ThermalDofComponent : int {
+  kBottom = 11,
+  kLower = 12,
+  kMiddle = 13,
+  kUpper = 14,
+  kTop = 15
+};
+
+/// @brief Stores one prescribed thermal source DOF value.
+struct ThermalBoundaryConditionDefinition {
+  EntityIndex node_index;
+  int dof;
+  double value;
+  SourceLocation location;
+};
+
+/// @brief Stores one signed inward CFlux source contribution.
+struct ThermalFluxLoadDefinition {
+  EntityIndex node_index;
+  int dof;
+  double power;
+  SourceLocation location;
+  std::size_t source_order{0U};
+};
+
+/// @brief Stores uncoupled transient heat-step controls without a time grid.
+struct ThermalStepDefinition {
+  std::string name;
+  double nominal_increment;
+  double time_period;
+  std::uint32_t maximum_increments{100U};
+  bool nlgeom{false};
+  std::vector<ThermalBoundaryConditionDefinition> boundaries;
+  std::vector<ThermalFluxLoadDefinition> loads;
+  SourceLocation location;
+
+  /// @brief Returns thermal boundaries in stable source order.
+  const std::vector<ThermalBoundaryConditionDefinition>& Boundaries() const
+      noexcept {
+    return boundaries;
+  }
+
+  /// @brief Returns additive thermal fluxes in stable source order.
+  const std::vector<ThermalFluxLoadDefinition>& Loads() const noexcept {
+    return loads;
+  }
+
+  static constexpr std::array<double, 5> kStationZeta{
+      {-1.0, -0.5, 0.0, 0.5, 1.0}};
+};
+
 /// @brief Stores a stable resolved node-set membership list.
 struct NodeSet {
   std::string name;
@@ -278,8 +375,11 @@ struct ModelDefinition {
   std::vector<EulerBeam3DDefinition> elements;
   std::vector<Mitc4ShellDefinition> shell_elements;
   std::vector<LinearElasticMaterial> materials;
+  std::vector<ThermalMaterialDefinition> thermal_materials;
   std::vector<GeneralBeamSection> sections;
   std::vector<ShellSection> shell_sections;
+  std::vector<ThermalShellSectionDefinition> thermal_sections;
+  std::vector<Ds4ThermalElementDefinition> thermal_elements;
   std::vector<ShellNodeInitialFrame> shell_node_initial_frames;
   std::vector<NodeSet> node_sets;
   std::vector<ElementSet> element_sets;
@@ -290,6 +390,9 @@ struct ModelDefinition {
   std::vector<MaterialDensityDefinition> densities;
   std::optional<FrequencyStepDefinition> frequency_step;
   std::optional<DynamicStepDefinition> dynamic_step;
+  std::vector<ThermalInitialConditionDefinition> thermal_initial_conditions;
+  std::vector<ThermalBoundaryConditionDefinition> thermal_model_boundaries;
+  std::optional<ThermalStepDefinition> thermal_step;
   std::vector<Diagnostic> warnings;
 };
 
diff --git a/src/fesa/analysis/analysis_model.cpp b/src/fesa/analysis/analysis_model.cpp
index e183d9b..759c19e 100644
--- a/src/fesa/analysis/analysis_model.cpp
+++ b/src/fesa/analysis/analysis_model.cpp
@@ -1,11 +1,24 @@
 #include "fesa/analysis/analysis_model.h"
 
+#include <algorithm>
+#include <cstddef>
+#include <stdexcept>
 #include <string>
 #include <vector>
 
 namespace fesa {
 
 Result<AnalysisModel> AnalysisModel::Create(const Domain& domain) {
+  if (domain.ThermalStep() != nullptr) {
+    if (!domain.Steps().Empty()) {
+      return Result<AnalysisModel>::Failure(Status::Failure(
+          FailureCategory::kModel,
+          {{Severity::kError, "mixed-thermal-procedure",
+            domain.ThermalStep()->location, "STEP", domain.ThermalStep()->name,
+            "Thermal and structural steps cannot coexist."}}));
+    }
+    return Result<AnalysisModel>::Success(AnalysisModel{domain});
+  }
   if (domain.Steps().Empty()) {
     return Result<AnalysisModel>::Failure(
         Status::Failure(FailureCategory::kInput,
@@ -29,10 +42,19 @@ Result<AnalysisModel> AnalysisModel::Create(const Domain& domain) {
 
 const Domain& AnalysisModel::GetDomain() const noexcept { return *domain_; }
 
-const StepDefinition& AnalysisModel::Step() const noexcept {
+const StepDefinition& AnalysisModel::Step() const {
+  if (is_thermal_) {
+    throw std::logic_error{"Thermal AnalysisModel requires ThermalStep()."};
+  }
   return domain_->Steps()[0U];
 }
 
+bool AnalysisModel::IsThermal() const noexcept { return is_thermal_; }
+
+const ThermalStepDefinition* AnalysisModel::ThermalStep() const noexcept {
+  return is_thermal_ ? domain_->ThermalStep() : nullptr;
+}
+
 const std::vector<EntityIndex>& AnalysisModel::ActiveElements() const noexcept {
   return active_elements_;
 }
@@ -65,12 +87,72 @@ const std::vector<EntityIndex>& AnalysisModel::ActiveLoads() const noexcept {
   return active_loads_;
 }
 
+const std::vector<EntityIndex>& AnalysisModel::ActiveThermalElements()
+    const noexcept {
+  return active_thermal_elements_;
+}
+
+const std::vector<EntityIndex>& AnalysisModel::ActiveThermalNodes()
+    const noexcept {
+  return active_thermal_nodes_;
+}
+
+const std::vector<EntityIndex>& AnalysisModel::ActiveThermalMaterials()
+    const noexcept {
+  return active_thermal_materials_;
+}
+
+const std::vector<EntityIndex>& AnalysisModel::ActiveThermalSections()
+    const noexcept {
+  return active_thermal_sections_;
+}
+
 const std::vector<SurfacePressureLoadDefinition>&
 AnalysisModel::ActiveSurfacePressureLoads() const noexcept {
+  static const std::vector<SurfacePressureLoadDefinition> kNoSurfacePressures;
+  if (is_thermal_) return kNoSurfacePressures;
   return Step().SurfacePressureLoads();
 }
 
 AnalysisModel::AnalysisModel(const Domain& domain) : domain_{&domain} {
+  if (domain.ThermalStep() != nullptr) {
+    is_thermal_ = true;
+    std::vector<bool> active_nodes(domain.Nodes().size(), false);
+    const std::size_t first_element = domain.Elements().Size() -
+                                      domain.ThermalElements().Size();
+    for (std::size_t index = 0U; index < domain.ThermalElements().Size();
+         ++index) {
+      active_thermal_elements_.push_back(
+          static_cast<EntityIndex>(first_element + index));
+      active_thermal_sections_.push_back(
+          domain.ThermalElements()[index].ThermalSectionIndex());
+      active_thermal_materials_.push_back(
+          domain.ThermalSections()[domain.ThermalElements()[index]
+                                       .ThermalSectionIndex()]
+              .material_index);
+      for (const EntityIndex node :
+           domain.ThermalElements()[index].NodeIndices()) {
+        active_nodes[node] = true;
+      }
+    }
+    for (std::size_t index = 0U; index < active_nodes.size(); ++index) {
+      if (active_nodes[index]) {
+        active_thermal_nodes_.push_back(static_cast<EntityIndex>(index));
+      }
+    }
+    std::sort(active_thermal_sections_.begin(), active_thermal_sections_.end());
+    active_thermal_sections_.erase(
+        std::unique(active_thermal_sections_.begin(),
+                    active_thermal_sections_.end()),
+        active_thermal_sections_.end());
+    std::sort(active_thermal_materials_.begin(),
+              active_thermal_materials_.end());
+    active_thermal_materials_.erase(
+        std::unique(active_thermal_materials_.begin(),
+                    active_thermal_materials_.end()),
+        active_thermal_materials_.end());
+    return;
+  }
   std::vector<bool> reachable_materials(domain.Materials().Size(), false);
   std::vector<bool> reachable_properties(domain.Properties().Size(), false);
   std::vector<bool> reachable_sections(domain.Sections().Size(), false);
diff --git a/src/fesa/materials/isotropic_thermal_material.cpp b/src/fesa/materials/isotropic_thermal_material.cpp
new file mode 100644
index 0000000..6b39c74
--- /dev/null
+++ b/src/fesa/materials/isotropic_thermal_material.cpp
@@ -0,0 +1,76 @@
+#include "fesa/materials/isotropic_thermal_material.h"
+
+#include <cmath>
+#include <string>
+#include <utility>
+
+namespace fesa {
+namespace {
+
+Status InvalidThermalMaterial(const SourceEntityId& source_id,
+                              const SourceLocation& location,
+                              const std::string& message) {
+  return Status::Failure(
+      FailureCategory::kModel,
+      {{Severity::kError, "invalid-thermal-material", location, "*MATERIAL",
+        source_id.source_label_text, message}});
+}
+
+}  // namespace
+
+Result<IsotropicThermalMaterial> IsotropicThermalMaterial::Create(
+    SourceEntityId source_id, std::string name, double conductivity,
+    double density, double specific_heat, SourceLocation location) {
+  if (name.empty() || source_id.source_label_text.empty()) {
+    return Result<IsotropicThermalMaterial>::Failure(InvalidThermalMaterial(
+        source_id, location, "Material source identity must not be empty."));
+  }
+  if (!std::isfinite(conductivity) || !(conductivity > 0.0) ||
+      !std::isfinite(density) || !(density > 0.0) ||
+      !std::isfinite(specific_heat) || !(specific_heat > 0.0)) {
+    return Result<IsotropicThermalMaterial>::Failure(InvalidThermalMaterial(
+        source_id, location,
+        "Conductivity, density, and specific heat must be finite and positive."));
+  }
+  return Result<IsotropicThermalMaterial>::Success(IsotropicThermalMaterial{
+      std::move(source_id), std::move(name), conductivity, density,
+      specific_heat, std::move(location)});
+}
+
+MaterialKind IsotropicThermalMaterial::Kind() const noexcept {
+  return MaterialKind::kIsotropicThermal;
+}
+
+const SourceEntityId& IsotropicThermalMaterial::SourceId() const noexcept {
+  return source_id_;
+}
+
+const SourceLocation& IsotropicThermalMaterial::Location() const noexcept {
+  return location_;
+}
+
+const std::string& IsotropicThermalMaterial::Name() const noexcept {
+  return name_;
+}
+
+double IsotropicThermalMaterial::Conductivity() const noexcept {
+  return conductivity_;
+}
+
+double IsotropicThermalMaterial::Density() const noexcept { return density_; }
+
+double IsotropicThermalMaterial::SpecificHeat() const noexcept {
+  return specific_heat_;
+}
+
+IsotropicThermalMaterial::IsotropicThermalMaterial(
+    SourceEntityId source_id, std::string name, double conductivity,
+    double density, double specific_heat, SourceLocation location)
+    : source_id_{std::move(source_id)},
+      name_{std::move(name)},
+      conductivity_{conductivity},
+      density_{density},
+      specific_heat_{specific_heat},
+      location_{std::move(location)} {}
+
+}  // namespace fesa
diff --git a/src/fesa/model/domain.cpp b/src/fesa/model/domain.cpp
index 9048948..1d30af1 100644
--- a/src/fesa/model/domain.cpp
+++ b/src/fesa/model/domain.cpp
@@ -1,11 +1,163 @@
 #include "fesa/model/domain.h"
 
 #include <cstddef>
+#include <cmath>
 #include <functional>
+#include <map>
 #include <memory>
+#include <string>
 #include <utility>
 
 namespace fesa {
+namespace {
+
+Status InvalidThermalDefinition(const SourceLocation& location,
+                                const std::string& identity,
+                                const std::string& message) {
+  return Status::Failure(
+      FailureCategory::kModel,
+      {{Severity::kError, "invalid-thermal-definition", location,
+        "THERMAL", identity, message}});
+}
+
+bool IsThermalDof(const int dof) { return dof == 11 || dof == 15; }
+
+Status ValidateThermalTargets(
+    const std::vector<ThermalBoundaryConditionDefinition>& boundaries,
+    const std::vector<bool>& active_nodes, const char* identity) {
+  std::map<std::pair<EntityIndex, int>, double> values;
+  for (const auto& boundary : boundaries) {
+    if (boundary.node_index >= active_nodes.size() ||
+        !active_nodes[boundary.node_index] || !IsThermalDof(boundary.dof) ||
+        !std::isfinite(boundary.value)) {
+      return InvalidThermalDefinition(boundary.location, identity,
+                                      "Thermal target is invalid or inactive.");
+    }
+    const auto inserted = values.emplace(
+        std::make_pair(boundary.node_index, boundary.dof), boundary.value);
+    if (!inserted.second && inserted.first->second != boundary.value) {
+      return InvalidThermalDefinition(
+          boundary.location, identity,
+          "Conflicting thermal boundary values target one node and DOF.");
+    }
+  }
+  return Status::Ok();
+}
+
+Status ValidateThermalDefinition(const ModelDefinition& definition) {
+  const bool has_structural = !definition.elements.empty() ||
+                              !definition.shell_elements.empty();
+  const bool has_thermal_data = !definition.thermal_materials.empty() ||
+                                !definition.thermal_sections.empty() ||
+                                definition.thermal_step.has_value() ||
+                                !definition.thermal_initial_conditions.empty() ||
+                                !definition.thermal_model_boundaries.empty();
+  if (has_structural && !definition.thermal_elements.empty()) {
+    return InvalidThermalDefinition({}, "DS4",
+                                    "Structural and DS4 elements cannot mix.");
+  }
+  if (definition.thermal_elements.empty()) {
+    if (has_thermal_data) {
+      return InvalidThermalDefinition({}, "DS4",
+                                      "Thermal definitions require a DS4 element.");
+    }
+    return Status::Ok();
+  }
+  if (!definition.thermal_step) {
+    return InvalidThermalDefinition({}, "DS4",
+                                    "A DS4 model requires one thermal step.");
+  }
+  std::map<std::string, bool> material_identities;
+  for (const auto& material : definition.thermal_materials) {
+    const std::string identity = material.source_id.instance_name + ":" +
+                                 material.source_id.source_label_text;
+    if (!material_identities.emplace(identity, true).second) {
+      return InvalidThermalDefinition(material.location, identity,
+                                      "Thermal material identity is duplicated.");
+    }
+    const auto candidate = IsotropicThermalMaterial::Create(
+        material.source_id, material.name, material.conductivity,
+        material.density, material.specific_heat, material.location);
+    if (!candidate.HasValue()) return candidate.GetStatus();
+  }
+  std::map<std::string, bool> section_identities;
+  for (const auto& section : definition.thermal_sections) {
+    const std::string identity = section.source_id.instance_name + ":" +
+                                 section.source_id.source_label_text;
+    if (section.source_id.source_label_text.empty() || section.name.empty() ||
+        !std::isfinite(section.thickness) || !(section.thickness > 0.0) ||
+        section.material_index >= definition.thermal_materials.size() ||
+        !section_identities.emplace(identity, true).second) {
+      return InvalidThermalDefinition(section.location, section.name,
+                                      "Thermal section is invalid.");
+    }
+  }
+
+  std::vector<bool> active_nodes(definition.nodes.size(), false);
+  std::map<std::string, bool> identities;
+  for (const auto& element : definition.thermal_elements) {
+    const std::string identity = element.source_id.instance_name + ":" +
+                                 element.source_id.source_label_text;
+    if (element.source_id.source_label_text.empty() ||
+        !identities.emplace(identity, true).second ||
+        element.section_index >= definition.thermal_sections.size()) {
+      return InvalidThermalDefinition(element.location, identity,
+                                      "DS4 section assignment is invalid.");
+    }
+    for (const EntityIndex node : element.node_indices) {
+      if (node >= active_nodes.size()) {
+        return InvalidThermalDefinition(element.location, identity,
+                                        "DS4 node target is invalid.");
+      }
+      active_nodes[node] = true;
+    }
+  }
+  std::map<EntityIndex, std::array<double, 5>> initial_values;
+  for (const auto& initial : definition.thermal_initial_conditions) {
+    if (initial.node_index >= active_nodes.size() ||
+        !active_nodes[initial.node_index]) {
+      return InvalidThermalDefinition(initial.location, "INITIAL CONDITIONS",
+                                      "Thermal initial target is inactive.");
+    }
+    for (const double value : initial.station_values) {
+      if (!std::isfinite(value) || value != initial.station_values[0U]) {
+        return InvalidThermalDefinition(
+            initial.location, "INITIAL CONDITIONS",
+            "Initial temperatures must be finite and equal at five stations.");
+      }
+    }
+    const auto inserted = initial_values.emplace(initial.node_index,
+                                                  initial.station_values);
+    if (!inserted.second && inserted.first->second != initial.station_values) {
+      return InvalidThermalDefinition(initial.location, "INITIAL CONDITIONS",
+                                      "Conflicting initial temperatures target "
+                                      "one node.");
+    }
+  }
+  Status model_status = ValidateThermalTargets(
+      definition.thermal_model_boundaries, active_nodes, "BOUNDARY");
+  if (!model_status.IsOk()) return model_status;
+  const auto& step = *definition.thermal_step;
+  if (!std::isfinite(step.nominal_increment) || !(step.nominal_increment > 0.0) ||
+      !std::isfinite(step.time_period) || !(step.time_period > 0.0) ||
+      step.maximum_increments == 0U || step.nlgeom) {
+    return InvalidThermalDefinition(step.location, step.name,
+                                    "Thermal step controls are invalid.");
+  }
+  Status step_status =
+      ValidateThermalTargets(step.boundaries, active_nodes, "STEP BOUNDARY");
+  if (!step_status.IsOk()) return step_status;
+  for (const auto& load : step.loads) {
+    if (load.node_index >= active_nodes.size() || !active_nodes[load.node_index] ||
+        !IsThermalDof(load.dof) || !std::isfinite(load.power)) {
+      return InvalidThermalDefinition(load.location, "CFLUX",
+                                      "Thermal flux target is invalid or inactive.");
+    }
+  }
+  return Status::Ok();
+}
+
+}  // namespace
 
 StepDefinition::StepDefinition(StaticStepDefinition definition)
     : name_{std::move(definition.name)},
@@ -174,6 +326,8 @@ const SourceLocation& StepDefinition::Location() const noexcept {
 }
 
 Result<Domain> Domain::Create(ModelDefinition definition) {
+  const Status thermal_status = ValidateThermalDefinition(definition);
+  if (!thermal_status.IsOk()) return Result<Domain>::Failure(thermal_status);
   return Result<Domain>::Success(Domain{std::move(definition)});
 }
 
@@ -205,6 +359,35 @@ Domain::LinearElasticMaterials() const noexcept {
   return linear_materials_view_;
 }
 
+const DomainCollectionView<IsotropicThermalMaterial>& Domain::ThermalMaterials()
+    const noexcept {
+  return thermal_materials_view_;
+}
+
+const DomainCollectionView<Ds4ThermalDefinition>& Domain::ThermalElements()
+    const noexcept {
+  return thermal_elements_view_;
+}
+
+const std::vector<ThermalShellSectionDefinition>& Domain::ThermalSections()
+    const noexcept {
+  return definition_.thermal_sections;
+}
+
+const std::vector<ThermalInitialConditionDefinition>&
+Domain::ThermalInitialConditions() const noexcept {
+  return definition_.thermal_initial_conditions;
+}
+
+const std::vector<ThermalBoundaryConditionDefinition>&
+Domain::ThermalModelBoundaries() const noexcept {
+  return definition_.thermal_model_boundaries;
+}
+
+const ThermalStepDefinition* Domain::ThermalStep() const noexcept {
+  return definition_.thermal_step ? &*definition_.thermal_step : nullptr;
+}
+
 const DomainCollectionView<ElementProperty>& Domain::Properties()
     const noexcept {
   return properties_view_;
@@ -278,6 +461,20 @@ Domain::Domain(ModelDefinition definition)
   }
   definition_.materials.clear();
 
+  const std::size_t structural_material_count = materials_.size();
+  materials_.reserve(materials_.size() + definition_.thermal_materials.size());
+  for (const auto& material : definition_.thermal_materials) {
+    const auto candidate = IsotropicThermalMaterial::Create(
+        material.source_id, material.name, material.conductivity,
+        material.density, material.specific_heat, material.location);
+    auto owned = std::make_unique<IsotropicThermalMaterial>(
+        std::move(candidate.Value()));
+    materials_view_.Add(*owned);
+    thermal_materials_view_.Add(*owned);
+    materials_.push_back(std::move(owned));
+  }
+  definition_.thermal_materials.clear();
+
   element_properties_.reserve(definition_.sections.size() +
                               definition_.shell_sections.size());
   for (auto& section : definition_.sections) {
@@ -316,6 +513,25 @@ Domain::Domain(ModelDefinition definition)
     element_definitions_.push_back(std::move(owned));
   }
   definition_.shell_elements.clear();
+  for (auto& element : definition_.thermal_elements) {
+    const EntityIndex material_index = static_cast<EntityIndex>(
+        structural_material_count +
+        definition_.thermal_sections[element.section_index].material_index);
+    auto owned = std::make_unique<Ds4ThermalDefinition>(
+        std::move(element.source_id), element.node_indices, element.section_index,
+        std::move(element.location));
+    owned->BindMaterialIndex(material_index);
+    elements_view_.Add(*owned);
+    thermal_elements_view_.Add(*owned);
+    element_definitions_.push_back(std::move(owned));
+  }
+  definition_.thermal_elements.clear();
+  if (definition_.thermal_step) {
+    for (std::size_t index = 0U;
+         index < definition_.thermal_step->loads.size(); ++index) {
+      definition_.thermal_step->loads[index].source_order = index;
+    }
+  }
 
   step_definitions_.reserve(definition_.steps.size() +
                             (definition_.frequency_step ? 1U : 0U) +
diff --git a/tests/unit/model/thermal_semantics_test.cpp b/tests/unit/model/thermal_semantics_test.cpp
new file mode 100644
index 0000000..cfa65ec
--- /dev/null
+++ b/tests/unit/model/thermal_semantics_test.cpp
@@ -0,0 +1,261 @@
+#include "fesa/model/domain.h"
+
+#include <gtest/gtest.h>
+
+#include <array>
+#include <limits>
+#include <stdexcept>
+#include <type_traits>
+#include <utility>
+
+#include "fesa/analysis/analysis_model.h"
+#include "fesa/elements/ds4_thermal.h"
+#include "fesa/materials/isotropic_thermal_material.h"
+
+namespace {
+
+fesa::SourceLocation Location() { return {"models/thermal.inp", 10U}; }
+
+fesa::ModelDefinition ValidDefinition() {
+  fesa::ModelDefinition definition{};
+  definition.source_path = "models/thermal.inp";
+  definition.nodes = {
+      {{"Plate", 1, "1"}, {0.0, 0.0, 0.0}, Location()},
+      {{"Plate", 2, "2"}, {1.0, 0.0, 0.0}, Location()},
+      {{"Plate", 3, "3"}, {1.0, 1.0, 0.0}, Location()},
+      {{"Plate", 4, "4"}, {0.0, 1.0, 0.0}, Location()},
+      {{"Plate", 5, "5"}, {2.0, 2.0, 0.0}, Location()}};
+  definition.thermal_materials = {
+      {{"", 1, "Thermal-Steel"}, "Thermal-Steel", 45.0, 7800.0, 500.0,
+       Location()}};
+  definition.thermal_sections = {
+      {{"", 1, "Section"}, "Section", 0.01, 0U, Location()}};
+  definition.thermal_elements = {{{"Plate", 8, "8"}, {0U, 1U, 2U, 3U},
+                                  0U, Location()}};
+  definition.thermal_initial_conditions = {{0U, 12.5, Location()}};
+  definition.thermal_model_boundaries = {{0U, 11, 10.0, Location()}};
+  definition.thermal_step = {"Heat", 0.2, 1.0, 100U, false,
+                             {{0U, 11, 20.0, Location()}},
+                             {{0U, 15, -5.0, Location()}}, Location()};
+  return definition;
+}
+
+}  // namespace
+
+TEST(ThermalDomain, StructuralDefinitionHasNoThermalStep) {
+  fesa::ModelDefinition definition{};
+  definition.steps = {{"Static", {}, {}, 0.1, 1.0, 0.01, 1.0, Location()}};
+  const auto result = fesa::Domain::Create(std::move(definition));
+  ASSERT_TRUE(result.HasValue());
+  EXPECT_EQ(result.Value().ThermalStep(), nullptr);
+}
+
+TEST(IsotropicThermalMaterial, RequiresFinitePositivePropertiesAndIdentity) {
+  const auto valid = fesa::IsotropicThermalMaterial::Create(
+      {"", 7, "Thermal-Steel"}, "Thermal-Steel", 45.0, 7800.0, 500.0,
+      Location());
+  ASSERT_TRUE(valid.HasValue());
+  EXPECT_EQ(valid.Value().Kind(), fesa::MaterialKind::kIsotropicThermal);
+  EXPECT_DOUBLE_EQ(valid.Value().Conductivity(), 45.0);
+  EXPECT_DOUBLE_EQ(valid.Value().Density(), 7800.0);
+  EXPECT_DOUBLE_EQ(valid.Value().SpecificHeat(), 500.0);
+
+  const auto invalid = fesa::IsotropicThermalMaterial::Create(
+      {"", 7, "Thermal-Steel"}, "Thermal-Steel",
+      std::numeric_limits<double>::infinity(), 7800.0, 500.0, Location());
+  ASSERT_FALSE(invalid.HasValue());
+  EXPECT_EQ(invalid.GetStatus().Category(), fesa::FailureCategory::kModel);
+
+  EXPECT_FALSE(fesa::IsotropicThermalMaterial::Create(
+                   {"", 7, "Thermal-Steel"}, "Thermal-Steel", 45.0, 0.0,
+                   500.0, Location())
+                   .HasValue());
+  EXPECT_FALSE(fesa::IsotropicThermalMaterial::Create(
+                   {"", 7, "Thermal-Steel"}, "Thermal-Steel", 45.0, 7800.0,
+                   std::numeric_limits<double>::quiet_NaN(), Location())
+                   .HasValue());
+}
+
+TEST(ThermalDomain, OwnsFiveStationDefinitionsWithoutEquationState) {
+  static_assert(!std::is_copy_constructible_v<fesa::Domain>);
+  static_assert(std::is_move_constructible_v<fesa::Domain>);
+  static_assert(std::has_virtual_destructor_v<fesa::ElementDefinition>);
+
+  const auto result = fesa::Domain::Create(ValidDefinition());
+  ASSERT_TRUE(result.HasValue());
+  const fesa::Domain& domain = result.Value();
+  ASSERT_NE(domain.ThermalStep(), nullptr);
+  EXPECT_EQ(domain.ThermalElements().Size(), 1U);
+  EXPECT_EQ(domain.ThermalElements()[0U].Kind(),
+            fesa::ElementDefinitionKind::kDs4Thermal);
+  EXPECT_EQ(domain.ThermalElements()[0U].PropertyIndex(),
+            fesa::Ds4ThermalDefinition::kNoPropertyIndex);
+  EXPECT_EQ(domain.ThermalElements()[0U].NodeIndices(),
+            (std::vector<fesa::EntityIndex>{0U, 1U, 2U, 3U}));
+  EXPECT_EQ(domain.ThermalInitialConditions()[0U].station_values,
+            (std::array<double, 5>{12.5, 12.5, 12.5, 12.5, 12.5}));
+  EXPECT_EQ(domain.ThermalStep()->Boundaries()[0U].value, 20.0);
+  EXPECT_EQ(domain.ThermalStep()->Loads()[0U].source_order, 0U);
+  EXPECT_EQ(static_cast<int>(fesa::ThermalDofComponent::kBottom), 11);
+  EXPECT_EQ(static_cast<int>(fesa::ThermalDofComponent::kTop), 15);
+}
+
+TEST(ThermalDomain, RejectsInactiveTargetsConflictsAndMixedElements) {
+  auto invalid = ValidDefinition();
+  invalid.thermal_initial_conditions.push_back({4U, 12.5, Location()});
+  const auto inactive = fesa::Domain::Create(std::move(invalid));
+  ASSERT_FALSE(inactive.HasValue());
+  EXPECT_EQ(inactive.GetStatus().Category(), fesa::FailureCategory::kModel);
+
+  invalid = ValidDefinition();
+  invalid.thermal_step->boundaries.push_back({0U, 11, 21.0, Location()});
+  const auto conflicting = fesa::Domain::Create(std::move(invalid));
+  ASSERT_FALSE(conflicting.HasValue());
+  EXPECT_EQ(conflicting.GetStatus().Category(), fesa::FailureCategory::kModel);
+
+  invalid = ValidDefinition();
+  invalid.elements.push_back(
+      {{"Plate", 9, "9"}, {0U, 1U}, 0U, 0U, Location()});
+  const auto mixed = fesa::Domain::Create(std::move(invalid));
+  ASSERT_FALSE(mixed.HasValue());
+  EXPECT_EQ(mixed.GetStatus().Category(), fesa::FailureCategory::kModel);
+
+  invalid = ValidDefinition();
+  invalid.thermal_sections[0U].thickness = 0.0;
+  const auto invalid_thickness = fesa::Domain::Create(std::move(invalid));
+  ASSERT_FALSE(invalid_thickness.HasValue());
+  EXPECT_EQ(invalid_thickness.GetStatus().Category(),
+            fesa::FailureCategory::kModel);
+}
+
+TEST(ThermalDomain, AllowsEqualInitialTargetsAndStepBoundaryOverride) {
+  auto definition = ValidDefinition();
+  definition.thermal_initial_conditions.push_back({0U, 12.5, Location()});
+  const auto equal_initial = fesa::Domain::Create(std::move(definition));
+  EXPECT_TRUE(equal_initial.HasValue());
+
+  definition = ValidDefinition();
+  definition.thermal_initial_conditions.push_back({0U, 13.0, Location()});
+  const auto conflicting_initial = fesa::Domain::Create(std::move(definition));
+  ASSERT_FALSE(conflicting_initial.HasValue());
+  EXPECT_EQ(conflicting_initial.GetStatus().Category(),
+            fesa::FailureCategory::kModel);
+
+  definition = ValidDefinition();
+  definition.thermal_step->boundaries.push_back({0U, 11, 20.0, Location()});
+  const auto overriding_step = fesa::Domain::Create(std::move(definition));
+  EXPECT_TRUE(overriding_step.HasValue());
+}
+
+TEST(ThermalDomain, RejectsInteriorThermalDofsAndInactiveFluxTargets) {
+  for (const int dof : {12, 13, 14}) {
+    auto definition = ValidDefinition();
+    definition.thermal_step->boundaries = {{0U, dof, 1.0, Location()}};
+    const auto boundary = fesa::Domain::Create(std::move(definition));
+    EXPECT_FALSE(boundary.HasValue());
+
+    definition = ValidDefinition();
+    definition.thermal_step->loads = {{0U, dof, 1.0, Location()}};
+    const auto flux = fesa::Domain::Create(std::move(definition));
+    EXPECT_FALSE(flux.HasValue());
+  }
+
+  auto definition = ValidDefinition();
+  definition.thermal_step->loads = {{4U, 11, 1.0, Location()}};
+  const auto inactive_flux = fesa::Domain::Create(std::move(definition));
+  EXPECT_FALSE(inactive_flux.HasValue());
+
+  definition = ValidDefinition();
+  definition.thermal_model_boundaries = {{4U, 11, 1.0, Location()}};
+  const auto inactive_boundary = fesa::Domain::Create(std::move(definition));
+  EXPECT_FALSE(inactive_boundary.HasValue());
+}
+
+TEST(ThermalDomain, PreservesAdditiveFluxSourceOrder) {
+  auto definition = ValidDefinition();
+  definition.thermal_step->loads = {{0U, 11, -5.0, Location()},
+                                    {1U, 11, 3.0, Location()}};
+  const auto result = fesa::Domain::Create(std::move(definition));
+  ASSERT_TRUE(result.HasValue());
+  const auto& loads = result.Value().ThermalStep()->Loads();
+  ASSERT_EQ(loads.size(), 2U);
+  EXPECT_EQ(loads[0U].source_order, 0U);
+  EXPECT_EQ(loads[1U].source_order, 1U);
+  EXPECT_DOUBLE_EQ(loads[0U].power, -5.0);
+  EXPECT_DOUBLE_EQ(loads[1U].power, 3.0);
+}
+
+TEST(ThermalAnalysisModel,
+     KeepsThermalSectionMaterialIndicesSeparateFromDomainMaterials) {
+  auto definition = ValidDefinition();
+  definition.materials = {{"Unused-Structural", 210.0e9, 0.3, Location()}};
+  definition.thermal_materials.push_back(
+      {{"", 2, "Thermal-Aluminum"}, "Thermal-Aluminum", 120.0, 2700.0,
+       900.0, Location()});
+  definition.thermal_sections = {
+      {{"", 1, "Section-A"}, "Section-A", 0.01, 1U, Location()},
+      {{"", 2, "Section-B"}, "Section-B", 0.02, 0U, Location()}};
+  definition.nodes.push_back(
+      {{"Plate", 6, "6"}, {3.0, 0.0, 0.0}, Location()});
+  definition.nodes.push_back(
+      {{"Plate", 7, "7"}, {4.0, 0.0, 0.0}, Location()});
+  definition.nodes.push_back(
+      {{"Plate", 8, "8"}, {4.0, 1.0, 0.0}, Location()});
+  definition.nodes.push_back(
+      {{"Plate", 9, "9"}, {3.0, 1.0, 0.0}, Location()});
+  definition.thermal_elements = {
+      {{"Plate", 8, "8"}, {0U, 1U, 2U, 3U}, 0U, Location()},
+      {{"Plate", 9, "9"}, {5U, 6U, 7U, 8U}, 1U, Location()}};
+
+  const auto domain_result = fesa::Domain::Create(std::move(definition));
+  ASSERT_TRUE(domain_result.HasValue());
+  const fesa::Domain& domain = domain_result.Value();
+  ASSERT_EQ(domain.Materials().Size(), 3U);
+  ASSERT_EQ(domain.ThermalMaterials().Size(), 2U);
+  EXPECT_EQ(domain.ThermalSections()[0U].material_index, 1U);
+  EXPECT_EQ(domain.ThermalSections()[1U].material_index, 0U);
+  EXPECT_EQ(domain.ThermalElements()[0U].MaterialIndex(), 2U);
+  EXPECT_EQ(domain.ThermalElements()[1U].MaterialIndex(), 1U);
+  EXPECT_EQ(domain.ThermalElements()[0U].PropertyIndex(),
+            fesa::Ds4ThermalDefinition::kNoPropertyIndex);
+
+  const auto model_result = fesa::AnalysisModel::Create(domain);
+  ASSERT_TRUE(model_result.HasValue());
+  EXPECT_EQ(model_result.Value().ActiveThermalMaterials(),
+            (std::vector<fesa::EntityIndex>{0U, 1U}));
+  EXPECT_EQ(model_result.Value().ActiveThermalSections(),
+            (std::vector<fesa::EntityIndex>{0U, 1U}));
+  EXPECT_TRUE(model_result.Value().ActiveProperties().empty());
+}
+
+TEST(ThermalAnalysisModel, ViewsConnectedDs4EntitiesWithoutOwningDomain) {
+  auto domain_result = fesa::Domain::Create(ValidDefinition());
+  ASSERT_TRUE(domain_result.HasValue());
+  const fesa::Domain& domain = domain_result.Value();
+
+  const auto model_result = fesa::AnalysisModel::Create(domain);
+  ASSERT_TRUE(model_result.HasValue());
+  const fesa::AnalysisModel& model = model_result.Value();
+  EXPECT_TRUE(model.IsThermal());
+  EXPECT_EQ(model.ThermalStep(), domain.ThermalStep());
+  EXPECT_THROW(model.Step(), std::logic_error);
+  EXPECT_EQ(model.ActiveThermalElements(),
+            (std::vector<fesa::EntityIndex>{0U}));
+  EXPECT_EQ(model.ActiveThermalNodes(),
+            (std::vector<fesa::EntityIndex>{0U, 1U, 2U, 3U}));
+  EXPECT_EQ(model.ActiveThermalMaterials(),
+            (std::vector<fesa::EntityIndex>{0U}));
+  EXPECT_EQ(model.ActiveThermalSections(),
+            (std::vector<fesa::EntityIndex>{0U}));
+  EXPECT_TRUE(model.ActiveProperties().empty());
+  EXPECT_EQ(domain.ThermalSections()[0U].material_index, 0U);
+  EXPECT_EQ(domain.ThermalElements()[0U].MaterialIndex(), 0U);
+  EXPECT_EQ(&model.GetDomain(), &domain);
+
+  fesa::Domain moved_domain = std::move(domain_result.Value());
+  const auto moved_model = fesa::AnalysisModel::Create(moved_domain);
+  ASSERT_TRUE(moved_model.HasValue());
+  EXPECT_EQ(&moved_model.Value().GetDomain(), &moved_domain);
+  EXPECT_DOUBLE_EQ(moved_domain.ThermalMaterials()[0U].Conductivity(), 45.0);
+  EXPECT_EQ(moved_domain.ThermalElements()[0U].SourceId().source_label, 8);
+}