feat: map thermal DS4 Abaqus input

src/fesa/io/abaqus/domain_mapping_context.cpp

@@ -17,6 +17,7 @@
 #include <utility>
 #include <vector>
 
+#include "fesa/analysis/thermal_time_grid.h"
 #include "fesa/core/ascii.h"
 #include "fesa/math/vector3.h"
 #include "fesa/model/shell_geometry.h"
@@ -49,7 +50,7 @@ struct RawNode {
 };
 
 struct RawElement {
-  enum class Type { kB33, kS4, kS4r };
+  enum class Type { kB33, kS4, kS4r, kDs4 };
 
   std::int64_t label;
   std::string label_text;
@@ -58,7 +59,7 @@ struct RawElement {
   SourceLocation location;
 };
 
-enum class ElementFamily { kBeam, kShell };
+enum class ElementFamily { kBeam, kShell, kThermalShell };
 
 struct RawSet {
   std::string name;
@@ -137,11 +138,33 @@ struct RawMaterial {
   double poisson_ratio{0.0};
   bool has_elastic{false};
   std::optional<double> density;
+  std::optional<double> conductivity;
+  std::optional<double> specific_heat;
   SourceLocation location;
   SourceLocation elastic_location;
   SourceLocation density_location;
 };
 
+struct RawThermalBoundary {
+  std::string target;
+  int dof;
+  double value;
+  SourceLocation location;
+};
+
+struct RawThermalInitial {
+  std::string target;
+  double value;
+  SourceLocation location;
+};
+
+struct RawThermalFlux {
+  std::string target;
+  int dof;
+  double power;
+  SourceLocation location;
+};
+
 struct RawStep {
   SourceLocation location;
   std::string name_raw;
@@ -156,6 +179,11 @@ struct RawStep {
   std::uint32_t maximum_increments{0U};
   std::array<double, 4> static_values{};
   std::optional<FrequencyStepDefinition> frequency;
+  bool has_heat_transfer{false};
+  double thermal_increment{0.0};
+  double thermal_period{0.0};
+  std::vector<RawThermalBoundary> thermal_boundaries;
+  std::vector<RawThermalFlux> thermal_fluxes;
   std::vector<PrescribedDisplacementDefinition> boundaries;
   std::vector<NodalLoad> loads;
   std::vector<RawSurfacePressureLoad> surface_pressure_loads;
@@ -630,7 +658,9 @@ class DomainMappingContext::Impl {
     }
 
     if (block.canonical_name != "ELASTIC" &&
-        block.canonical_name != "DENSITY") {
+        block.canonical_name != "DENSITY" &&
+        block.canonical_name != "CONDUCTIVITY" &&
+        block.canonical_name != "SPECIFIC HEAT") {
       material_eligible_.reset();
     }
     pending_section_.reset();
@@ -667,6 +697,18 @@ class DomainMappingContext::Impl {
       AddIgnoredWarning(block);
       return;
     }
+    if (block.canonical_name == "PHYSICAL CONSTANTS") {
+      if (!ValidateParameters(block, {"ABSOLUTE ZERO"}) ||
+          !RequireNoData(block))
+        return;
+      const auto* absolute_zero =
+          RequiredParameterValue(block, "ABSOLUTE ZERO");
+      if (absolute_zero == nullptr) return;
+      double value = 0.0;
+      ParseDouble(*absolute_zero, value, block.location, block.canonical_name);
+      physical_constants_seen_ = !failure_;
+      return;
+    }
     if (block.canonical_name == "PART") {
       ParsePartStart(block);
       return;
@@ -691,6 +733,18 @@ class DomainMappingContext::Impl {
       ParseDensity(block);
       return;
     }
+    if (block.canonical_name == "CONDUCTIVITY") {
+      ParseThermalProperty(block, true);
+      return;
+    }
+    if (block.canonical_name == "SPECIFIC HEAT") {
+      ParseThermalProperty(block, false);
+      return;
+    }
+    if (block.canonical_name == "INITIAL CONDITIONS") {
+      ParseInitialConditions(block);
+      return;
+    }
     if (block.canonical_name == "BOUNDARY") {
       if (!assembly_seen_ || materials_.empty()) {
         InvalidKeywordLocation(
@@ -698,7 +752,11 @@ class DomainMappingContext::Impl {
         return;
       }
       model_boundary_seen_ = true;
-      ParseBoundary(block, model_boundaries_);
+      if (model_element_family_ == ElementFamily::kThermalShell) {
+        ParseThermalBoundary(block, thermal_model_boundaries_);
+      } else {
+        ParseBoundary(block, model_boundaries_);
+      }
       return;
     }
     if (block.canonical_name == "STEP") {
@@ -750,7 +808,9 @@ class DomainMappingContext::Impl {
         return;
       }
       const RawPart& part = parts_[*current_part_];
-      const bool shell_part = model_element_family_ == ElementFamily::kShell;
+      const bool shell_part =
+          model_element_family_ == ElementFamily::kShell ||
+          model_element_family_ == ElementFamily::kThermalShell;
       if (part.nodes.empty() || part.elements.empty() ||
           (!shell_part && part.sections.empty())) {
         InvalidKeywordLocation(block,
@@ -939,12 +999,16 @@ class DomainMappingContext::Impl {
       element_type = RawElement::Type::kS4r;
       element_family = ElementFamily::kShell;
       expected_field_count = 5U;
+    } else if (AsciiCaseInsensitiveEquals(*type, "DS4")) {
+      element_type = RawElement::Type::kDs4;
+      element_family = ElementFamily::kThermalShell;
+      expected_field_count = 5U;
     } else {
-      InputFailure(
+      InputFailure("unsupported-element-formulation", block.location,
-          "unsupported-element-formulation", block.location,
+                   block.canonical_name, *type,
-          block.canonical_name, *type,
+                   "Only TYPE=B33, TYPE=S4, TYPE=S4R, and TYPE=DS4 belong to "
-          "Only TYPE=B33, TYPE=S4, and TYPE=S4R belong to the approved "
+                   "the approved "
-          "element subsets.");
+                   "element subsets.");
       return;
     }
     if (model_element_family_ && *model_element_family_ != element_family) {
@@ -992,7 +1056,7 @@ class DomainMappingContext::Impl {
         }
         element.node_labels.push_back(node_label);
       }
-      if (element_family == ElementFamily::kShell) {
+      if (element_family != ElementFamily::kBeam) {
         const std::set<std::int64_t> distinct_nodes{element.node_labels.begin(),
                                                     element.node_labels.end()};
         if (distinct_nodes.size() != element.node_labels.size()) {
@@ -1056,13 +1120,23 @@ class DomainMappingContext::Impl {
       return;
     }
     if (block.data[0].fields.size() == 2U) {
-      if (!ParsePositiveSourceLabel(block.data[0].fields[1]).HasValue()) {
+      const auto integration_points =
+          ParsePositiveSourceLabel(block.data[0].fields[1]);
+      if (!integration_points.HasValue()) {
         InputFailure("unsupported-shell-section-option", block.data[0].location,
                      block.canonical_name, block.data[0].fields[1],
                      "The optional shell integration-point field must be a "
                      "positive integer.");
         return;
       }
+      if (model_element_family_ == ElementFamily::kThermalShell &&
+          integration_points.Value() != 5) {
+        InputFailure(
+            "unsupported-shell-section-option", block.data[0].location,
+            block.canonical_name, block.data[0].fields[1],
+            "Thermal DS4 sections require exactly five integration points.");
+        return;
+      }
     }
     part.shell_sections.push_back(
         {*element_set, *material, thickness, block.location});
@@ -1551,8 +1625,16 @@ class DomainMappingContext::Impl {
                    *name, "Material names are globally unique.");
       return;
     }
-    materials_.push_back(
+    materials_.push_back({*name,
-        {*name, 0.0, 0.0, false, std::nullopt, block.location, {}, {}});
+                          0.0,
+                          0.0,
+                          false,
+                          std::nullopt,
+                          std::nullopt,
+                          std::nullopt,
+                          block.location,
+                          {},
+                          {}});
     material_eligible_ = materials_.size() - 1U;
   }
 
@@ -1701,15 +1783,123 @@ class DomainMappingContext::Impl {
     material.density_location = block.data[0].location;
   }
 
+  void ParseThermalProperty(const KeywordBlock& block, bool conductivity) {
+    if (model_element_family_ != ElementFamily::kThermalShell ||
+        !material_eligible_ || !ValidateParameters(block, {})) {
+      if (!failure_) {
+        InputFailure("invalid-keyword-location", block.location,
+                     block.canonical_name, "",
+                     "Thermal properties belong to a DS4 material.");
+      }
+      return;
+    }
+    RawMaterial& material = materials_[*material_eligible_];
+    std::optional<double>& destination =
+        conductivity ? material.conductivity : material.specific_heat;
+    if (destination || block.data.size() != 1U ||
+        (block.data[0].fields.size() != 1U &&
+         (block.data[0].fields.size() != 2U ||
+          !block.data[0].fields[1].empty()))) {
+      InputFailure("invalid-data-arity", block.location, block.canonical_name,
+                   material.name,
+                   "Thermal property requires one constant row.");
+      return;
+    }
+    double value = 0.0;
+    if (!ParseModelDouble(
+            block.data[0].fields[0], value, block.data[0].location,
+            block.canonical_name, "invalid-thermal-material",
+            "Thermal material values must be finite and positive.") ||
+        !(value > 0.0)) {
+      if (!failure_) {
+        ModelFailure("invalid-thermal-material", block.data[0].location,
+                     block.canonical_name, material.name,
+                     "Thermal material values must be finite and positive.");
+      }
+      return;
+    }
+    destination = value;
+  }
+
+  void ParseInitialConditions(const KeywordBlock& block) {
+    if (model_element_family_ != ElementFamily::kThermalShell ||
+        !ValidateParameters(block, {"TYPE"}) || block.data.empty()) {
+      if (!failure_) {
+        InputFailure("invalid-data-arity", block.location, block.canonical_name,
+                     "", "Thermal initial conditions require data rows.");
+      }
+      return;
+    }
+    const auto* type = RequiredParameterValue(block, "TYPE");
+    if (type == nullptr || !AsciiCaseInsensitiveEquals(*type, "TEMPERATURE")) {
+      if (!failure_)
+        InputFailure("unsupported-keyword", block.location,
+                     block.canonical_name, "TYPE",
+                     "Only TYPE=TEMPERATURE is supported.");
+      return;
+    }
+    for (const auto& row : block.data) {
+      if (row.fields.size() != 2U || row.fields[0].empty()) {
+        InputFailure("invalid-data-arity", row.location, block.canonical_name,
+                     "", "Initial temperatures require target, value.");
+        return;
+      }
+      double value = 0.0;
+      if (!ParseDouble(row.fields[1], value, row.location,
+                       block.canonical_name))
+        return;
+      thermal_initial_conditions_.push_back(
+          {row.fields[0], value, row.location});
+    }
+  }
+
+  void ParseThermalBoundary(const KeywordBlock& block,
+                            std::vector<RawThermalBoundary>& destination) {
+    if (!ValidateParameters(block, {}) || block.data.empty()) {
+      if (!failure_)
+        InputFailure("invalid-data-arity", block.location, block.canonical_name,
+                     "", "Thermal BOUNDARY requires data rows.");
+      return;
+    }
+    for (const auto& row : block.data) {
+      if ((row.fields.size() != 3U && row.fields.size() != 4U) ||
+          row.fields[0].empty()) {
+        InputFailure("invalid-data-arity", row.location, block.canonical_name,
+                     "", "Thermal boundaries require target, DOF, DOF, value.");
+        return;
+      }
+      std::int64_t first = 0;
+      std::int64_t last = 0;
+      if (!ParseInteger(row.fields[1], first, row.location,
+                        block.canonical_name, true) ||
+          !ParseInteger(row.fields[2], last, row.location, block.canonical_name,
+                        true) ||
+          first != last || (first != 11 && first != 15)) {
+        InputFailure("invalid-dof", row.location, block.canonical_name,
+                     row.fields[0], "Thermal boundaries support DOF 11 or 15.");
+        return;
+      }
+      double value = 0.0;
+      if (row.fields.size() == 4U &&
+          !ParseDouble(row.fields[3], value, row.location,
+                       block.canonical_name))
+        return;
+      destination.push_back(
+          {row.fields[0], static_cast<int>(first), value, row.location});
+    }
+  }
+
   void ParseBoundary(
       const KeywordBlock& block,
       std::vector<PrescribedDisplacementDefinition>& destination) {
     const bool dynamic_input = std::any_of(
-        input_.blocks.begin(), input_.blocks.end(),
+        input_.blocks.begin(), input_.blocks.end(), [](const auto& candidate) {
-        [](const auto& candidate) { return candidate.canonical_name == "DYNAMIC"; });
+          return candidate.canonical_name == "DYNAMIC";
-    if (!ValidateParameters(block, {}, dynamic_input
+        });
-                                ? "unsupported-dynamic-boundary"
+    if (!ValidateParameters(block, {},
-                                : "invalid-keyword-parameter") || block.data.empty()) {
+                            dynamic_input ? "unsupported-dynamic-boundary"
+                                          : "invalid-keyword-parameter") ||
+        block.data.empty()) {
       if (!failure_) {
         InputFailure("invalid-data-arity", block.location, block.canonical_name,
                      "", "BOUNDARY requires one or more data rows.");
@@ -1790,7 +1980,8 @@ class DomainMappingContext::Impl {
       if (step_.dynamic && dof > 3) {
         InputFailure("unsupported-dynamic-load", row.location,
                      block.canonical_name, row.fields[1],
-                     "Dynamic CLOAD supports global translational DOF 1 through 3 only.");
+                     "Dynamic CLOAD supports global translational DOF 1 "
+                     "through 3 only.");
         return;
       }
       loads.push_back(
@@ -1927,8 +2118,8 @@ class DomainMappingContext::Impl {
       std::int64_t parsed_limit = 0;
       if (!ParseInteger(*increment_limit->value, parsed_limit, block.location,
                         block.canonical_name, true) ||
-          parsed_limit >
+          parsed_limit > static_cast<std::int64_t>(
-              static_cast<std::int64_t>(std::numeric_limits<std::uint32_t>::max())) {
+                             std::numeric_limits<std::uint32_t>::max())) {
         if (!failure_) {
           InputFailure("invalid-dynamic-data", block.location,
                        block.canonical_name, "INC",
@@ -1959,7 +2150,8 @@ class DomainMappingContext::Impl {
       if (!ValidateParameters(block, {}) || !RequireNoData(block)) {
         return;
       }
-      if (!step_.has_static && !step_.frequency && !step_.dynamic) {
+      if (!step_.has_static && !step_.frequency && !step_.dynamic &&
+          !step_.has_heat_transfer) {
         InvalidKeywordLocation(
             block, "The sole step requires one supported procedure keyword.");
         return;
@@ -1973,7 +2165,17 @@ class DomainMappingContext::Impl {
                    "A second or nested step is unsupported.");
       return;
     }
-    if (block.canonical_name == "STATIC") {
+    if (block.canonical_name == "HEAT TRANSFER") {
+      active_output_ = false;
+      if (step_.has_static || step_.frequency || step_.dynamic ||
+          step_.has_heat_transfer || !step_.boundaries.empty() ||
+          !step_.loads.empty() || step_no_op_seen_) {
+        InvalidKeywordLocation(block,
+                               "HEAT TRANSFER must be the first procedure.");
+        return;
+      }
+      ParseHeatTransfer(block);
+    } else if (block.canonical_name == "STATIC") {
       active_output_ = false;
       if (!step_.boundaries.empty() || !step_.loads.empty() ||
           !step_.surface_pressure_loads.empty() || step_no_op_seen_) {
@@ -2004,8 +2206,9 @@ class DomainMappingContext::Impl {
       ParseDynamic(block);
     } else if (block.canonical_name == "BOUNDARY") {
       active_output_ = false;
-      if ((!step_.has_static && !step_.frequency && !step_.dynamic) || step_load_seen_ ||
+      if ((!step_.has_static && !step_.frequency && !step_.dynamic &&
-          step_no_op_seen_) {
+           !step_.has_heat_transfer) ||
+          step_load_seen_ || step_no_op_seen_) {
         InvalidKeywordLocation(
             block,
             "Step boundaries follow STATIC and precede loads "
@@ -2018,7 +2221,11 @@ class DomainMappingContext::Impl {
                      "Dynamic BOUNDARY does not support keyword parameters.");
         return;
       }
-      ParseBoundary(block, step_.boundaries);
+      if (step_.has_heat_transfer) {
+        ParseThermalBoundary(block, step_.thermal_boundaries);
+      } else {
+        ParseBoundary(block, step_.boundaries);
+      }
       if (!failure_ && step_.frequency) {
         ValidateEigenBoundaries(step_.boundaries);
       }
@@ -2038,6 +2245,15 @@ class DomainMappingContext::Impl {
         step_.next_load_ordinal += block.data.size();
       }
       step_load_seen_ = !failure_;
+    } else if (block.canonical_name == "CFLUX") {
+      active_output_ = false;
+      if (!step_.has_heat_transfer || step_no_op_seen_) {
+        InvalidKeywordLocation(block,
+                               "CFLUX follows HEAT TRANSFER and boundaries.");
+        return;
+      }
+      ParseCflux(block);
+      step_load_seen_ = !failure_;
     } else if (block.canonical_name == "DSLOAD") {
       active_output_ = false;
       if (!step_.has_static || step_.dynamic || step_no_op_seen_) {
@@ -2050,7 +2266,8 @@ class DomainMappingContext::Impl {
       step_load_seen_ = !failure_;
     } else if (block.canonical_name == "RESTART") {
       active_output_ = false;
-      if (!step_.has_static && !step_.frequency && !step_.dynamic) {
+      if (!step_.has_static && !step_.frequency && !step_.dynamic &&
+          !step_.has_heat_transfer) {
         InvalidKeywordLocation(
             block, "Step no-op requests follow STATIC, boundaries, and loads.");
         return;
@@ -2065,7 +2282,8 @@ class DomainMappingContext::Impl {
           if (!failure_) {
             InputFailure("unsupported-keyword", block.location,
                          block.canonical_name, "",
-                         "Dynamic and frequency RESTART accept only WRITE, FREQUENCY=0.");
+                         "Dynamic and frequency RESTART accept only WRITE, "
+                         "FREQUENCY=0.");
           }
           return;
         }
@@ -2074,7 +2292,8 @@ class DomainMappingContext::Impl {
         AddIgnoredWarning(block);
       }
     } else if (block.canonical_name == "OUTPUT") {
-      if (!step_.has_static && !step_.frequency && !step_.dynamic) {
+      if (!step_.has_static && !step_.frequency && !step_.dynamic &&
+          !step_.has_heat_transfer) {
         InvalidKeywordLocation(
             block, "Step no-op requests follow STATIC, boundaries, and loads.");
         return;
@@ -2256,11 +2475,13 @@ class DomainMappingContext::Impl {
   void ParseDynamic(const KeywordBlock& block) {
     if (step_.perturbation_explicit || !step_.nlgeom ||
         model_element_family_ != ElementFamily::kShell ||
-        !ValidateParameters(block, {"DIRECT"}, "unsupported-dynamic-parameter")) {
+        !ValidateParameters(block, {"DIRECT"},
+                            "unsupported-dynamic-parameter")) {
       if (!failure_) {
-        InputFailure("unsupported-analysis-procedure", block.location,
+        InputFailure(
-                     block.canonical_name, "",
+            "unsupported-analysis-procedure", block.location,
-                     "DYNAMIC requires the approved nonlinear shell configuration.");
+            block.canonical_name, "",
+            "DYNAMIC requires the approved nonlinear shell configuration.");
       }
       return;
     }
@@ -2270,8 +2491,9 @@ class DomainMappingContext::Impl {
         std::any_of(block.data[0].fields.begin() + 2,
                     block.data[0].fields.end(),
                     [](const auto& field) { return !field.empty(); })) {
-      InputFailure("invalid-dynamic-data", block.location, block.canonical_name,
+      InputFailure(
-                   "", "DYNAMIC,DIRECT requires one row with increment and time period.");
+          "invalid-dynamic-data", block.location, block.canonical_name, "",
+          "DYNAMIC,DIRECT requires one row with increment and time period.");
       return;
     }
     DynamicStepDefinition dynamic{step_.name_raw,
@@ -2294,9 +2516,10 @@ class DomainMappingContext::Impl {
                             block.data[0].location, block.canonical_name) ||
         !(dynamic.initial_increment > 0.0) || !(dynamic.time_period > 0.0)) {
       if (!failure_) {
-        InputFailure("invalid-dynamic-data", block.data[0].location,
+        InputFailure(
-                     block.canonical_name, "",
+            "invalid-dynamic-data", block.data[0].location,
-                     "Dynamic increment and time period must be finite and positive.");
+            block.canonical_name, "",
+            "Dynamic increment and time period must be finite and positive.");
       }
       return;
     }
@@ -2308,7 +2531,8 @@ class DomainMappingContext::Impl {
         dynamic.initial_increment * rounded_count != dynamic.time_period) {
       InputFailure("invalid-dynamic-schedule", block.data[0].location,
                    block.canonical_name, "",
-                   "DYNAMIC,DIRECT requires an exact positive integer number of increments.");
+                   "DYNAMIC,DIRECT requires an exact positive integer number "
+                   "of increments.");
       return;
     }
     dynamic.increment_count = static_cast<std::uint32_t>(rounded_count);
@@ -2321,6 +2545,85 @@ class DomainMappingContext::Impl {
     step_.dynamic = std::move(dynamic);
   }
 
+  void ParseHeatTransfer(const KeywordBlock& block) {
+    if (model_element_family_ != ElementFamily::kThermalShell || step_.nlgeom ||
+        step_.perturbation_explicit ||
+        !ValidateParameters(block, {"END"}, "unsupported-thermal-option") ||
+        block.data.size() != 1U) {
+      if (!failure_)
+        InputFailure("unsupported-analysis-procedure", block.location,
+                     block.canonical_name, "",
+                     "HEAT TRANSFER requires the DS4 thermal subset.");
+      return;
+    }
+    const auto* end = Parameter(block, "END");
+    if (end != nullptr &&
+        (!end->value || !AsciiCaseInsensitiveEquals(*end->value, "PERIOD"))) {
+      InputFailure("unsupported-thermal-option", block.location,
+                   block.canonical_name, "END",
+                   "HEAT TRANSFER accepts only END=PERIOD.");
+      return;
+    }
+    const auto fields = WithoutTrailingEmpty(block.data[0].fields);
+    if (fields.size() != 2U) {
+      InputFailure("invalid-thermal-time-grid", block.data[0].location,
+                   block.canonical_name, "",
+                   "HEAT TRANSFER requires increment and time period.");
+      return;
+    }
+    if (!ParseDouble(fields[0], step_.thermal_increment, block.data[0].location,
+                     block.canonical_name) ||
+        !ParseDouble(fields[1], step_.thermal_period, block.data[0].location,
+                     block.canonical_name) ||
+        !(step_.thermal_increment > 0.0) || !(step_.thermal_period > 0.0)) {
+      if (!failure_)
+        InputFailure(
+            "invalid-thermal-time-grid", block.data[0].location,
+            block.canonical_name, "",
+            "Thermal increment and period must be finite and positive.");
+      return;
+    }
+    const auto grid =
+        ThermalTimeGrid::Create(step_.thermal_period, step_.thermal_increment,
+                                step_.maximum_increments);
+    if (!grid.HasValue()) {
+      failure_ = MappingFailure{
+          grid.GetStatus().Category().value_or(FailureCategory::kInput),
+          grid.GetStatus().Diagnostics().front()};
+      return;
+    }
+    step_.has_heat_transfer = true;
+  }
+
+  void ParseCflux(const KeywordBlock& block) {
+    if (!ValidateParameters(block, {}) || block.data.empty()) {
+      if (!failure_)
+        InputFailure("invalid-data-arity", block.location, block.canonical_name,
+                     "", "CFLUX requires target, DOF, power rows.");
+      return;
+    }
+    for (const auto& row : block.data) {
+      if (row.fields.size() != 3U || row.fields[0].empty()) {
+        InputFailure("invalid-data-arity", row.location, block.canonical_name,
+                     "", "CFLUX rows require target, DOF, power.");
+        return;
+      }
+      std::int64_t dof = 0;
+      double power = 0.0;
+      if (!ParseInteger(row.fields[1], dof, row.location, block.canonical_name,
+                        true) ||
+          !ParseDouble(row.fields[2], power, row.location,
+                       block.canonical_name) ||
+          (dof != 11 && dof != 15)) {
+        InputFailure("invalid-dof", row.location, block.canonical_name,
+                     row.fields[0], "CFLUX supports DOF 11 or 15.");
+        return;
+      }
+      step_.thermal_fluxes.push_back(
+          {row.fields[0], static_cast<int>(dof), power, row.location});
+    }
+  }
+
   void ParseStatic(const KeywordBlock& block) {
     if (step_.perturbation_explicit || step_.frequency) {
       InputFailure("unsupported-analysis-procedure", block.location,
@@ -2376,9 +2679,11 @@ class DomainMappingContext::Impl {
 
   void ParseOutputRoot(const KeywordBlock& block) {
     if (!ValidateParameters(
-            block, step_.dynamic
+            block, (step_.dynamic || step_.has_heat_transfer)
-                       ? std::vector<std::string_view>{"FIELD", "HISTORY", "VARIABLE", "FREQUENCY"}
+                       ? std::vector<std::string_view>{"FIELD", "HISTORY",
-                       : std::vector<std::string_view>{"FIELD", "HISTORY", "VARIABLE"})) {
+                                                       "VARIABLE", "FREQUENCY"}
+                       : std::vector<std::string_view>{"FIELD", "HISTORY",
+                                                       "VARIABLE"})) {
       return;
     }
     const auto* field = Parameter(block, "FIELD");
@@ -2403,24 +2708,25 @@ class DomainMappingContext::Impl {
     }
     if (step_.dynamic) {
       const auto* frequency = Parameter(block, "FREQUENCY");
-      if (!ValidateParameters(block, {"FIELD", "HISTORY", "VARIABLE", "FREQUENCY"}) ||
+      if (!ValidateParameters(block,
+                              {"FIELD", "HISTORY", "VARIABLE", "FREQUENCY"}) ||
           variable == nullptr || !variable->value ||
           !AsciiCaseInsensitiveEquals(*variable->value, "PRESELECT") ||
           (frequency != nullptr &&
            (!frequency->value || *frequency->value != "1")) ||
           !RequireNoData(block)) {
         if (!failure_) {
-          InputFailure("unsupported-keyword", block.location,
+          InputFailure(
-                       block.canonical_name, "",
+              "unsupported-keyword", block.location, block.canonical_name, "",
-                       "Dynamic OUTPUT accepts only VARIABLE=PRESELECT, FREQUENCY=1.");
+              "Dynamic OUTPUT accepts only VARIABLE=PRESELECT, FREQUENCY=1.");
         }
         return;
       }
     } else if (step_.frequency &&
-        ((variable != nullptr &&
+               ((variable != nullptr &&
-          (!variable->value ||
+                 (!variable->value || !AsciiCaseInsensitiveEquals(
-           !AsciiCaseInsensitiveEquals(*variable->value, "PRESELECT"))) ||
+                                          *variable->value, "PRESELECT"))) ||
-         !RequireNoData(block))) {
+                !RequireNoData(block))) {
       if (!failure_) {
         InputFailure("unsupported-keyword", block.location,
                      block.canonical_name, "",
@@ -2467,12 +2773,12 @@ class DomainMappingContext::Impl {
   void AddIgnoredWarning(const KeywordBlock& block) {
     // One warning per allowlisted keyword keeps no-op provenance stable;
     // subordinate variable rows remain attached to that keyword record.
-    definition_.warnings.push_back({Severity::kWarning,
+    definition_.warnings.push_back(
-                                    step_.dynamic ? "ignored-output-request"
+        {Severity::kWarning,
-                                                  : "ignored-input-keyword",
+         step_.dynamic ? "ignored-output-request" : "ignored-input-keyword",
-                                    block.location, block.canonical_name, "",
+         block.location, block.canonical_name, "",
-                                    "The allowlisted Abaqus keyword is "
+         "The allowlisted Abaqus keyword is "
-                                    "ignored without semantic effect."});
+         "ignored without semantic effect."});
   }
 
   void RejectUnknown(const KeywordBlock& block) {
@@ -2540,7 +2846,8 @@ class DomainMappingContext::Impl {
   void ValidateCardinality() {
     if (parts_.empty() || !assembly_seen_ || instances_.empty() ||
         materials_.empty() || !step_seen_ ||
-        (!step_.has_static && !step_.frequency && !step_.dynamic)) {
+        (!step_.has_static && !step_.frequency && !step_.dynamic &&
+         !step_.has_heat_transfer)) {
       const auto location = input_.blocks.empty()
                                 ? SourceLocation{input_.source_path, 0U}
                                 : input_.blocks.back().location;
@@ -2552,6 +2859,23 @@ class DomainMappingContext::Impl {
 
   void FinalizeMaterials() {
     for (const auto& material : materials_) {
+      if (model_element_family_ == ElementFamily::kThermalShell) {
+        if (!material.density || !material.conductivity ||
+            !material.specific_heat) {
+          ModelFailure("missing-thermal-material-property", material.location,
+                       "MATERIAL", material.name,
+                       "DS4 materials require DENSITY, CONDUCTIVITY, and "
+                       "SPECIFIC HEAT.");
+          return;
+        }
+        definition_.thermal_materials.push_back({{"", 0, material.name},
+                                                 material.name,
+                                                 *material.conductivity,
+                                                 *material.density,
+                                                 *material.specific_heat,
+                                                 material.location});
+        continue;
+      }
       if (!material.has_elastic) {
         InputFailure("unresolved-reference", material.location, "MATERIAL",
                      material.name,
@@ -2591,6 +2915,31 @@ class DomainMappingContext::Impl {
   }
 
   void FinalizeDensities() {
+    if (model_element_family_ == ElementFamily::kThermalShell) {
+      for (const auto& block : input_.blocks) {
+        if (block.canonical_name != "PREPRINT") continue;
+        if (!ValidateParameters(block,
+                                {"ECHO", "MODEL", "HISTORY", "CONTACT"}) ||
+            block.parameters.size() != 4U)
+          return;
+        for (const auto& parameter : block.parameters) {
+          if (!parameter.value ||
+              !AsciiCaseInsensitiveEquals(*parameter.value, "NO")) {
+            InputFailure("unsupported-keyword", block.location,
+                         block.canonical_name, parameter.name,
+                         "Thermal PREPRINT parameters must equal NO.");
+            return;
+          }
+        }
+      }
+      return;
+    }
+    if (physical_constants_seen_) {
+      InputFailure("unsupported-keyword", step_.location, "PHYSICAL CONSTANTS",
+                   "",
+                   "PHYSICAL CONSTANTS is supported only by thermal input.");
+      return;
+    }
     if (!step_.frequency && !step_.dynamic) {
       if (std::any_of(materials_.begin(), materials_.end(),
                       [](const RawMaterial& material) {
@@ -2641,9 +2990,10 @@ class DomainMappingContext::Impl {
           }
         }
         if (active) {
-          ModelFailure("missing-material-density", material.location,
+          ModelFailure(
-                       "MATERIAL", material.name,
+              "missing-material-density", material.location, "MATERIAL",
-                       "Each active dynamic or frequency material requires DENSITY.");
+              material.name,
+              "Each active dynamic or frequency material requires DENSITY.");
           return;
         }
         continue;
@@ -2696,7 +3046,8 @@ class DomainMappingContext::Impl {
         if (!failure_) {
           InputFailure("unsupported-keyword", block.location,
                        block.canonical_name, "",
-                       "DYNAMIC requires PREPRINT with ECHO, MODEL, HISTORY, and CONTACT parameters.");
+                       "DYNAMIC requires PREPRINT with ECHO, MODEL, HISTORY, "
+                       "and CONTACT parameters.");
         }
         return;
       }
@@ -2715,7 +3066,9 @@ class DomainMappingContext::Impl {
   void FinalizePartsAndSections() {
     part_section_assignments_.resize(parts_.size());
     part_shell_section_assignments_.resize(parts_.size());
-    const bool shell_model = model_element_family_ == ElementFamily::kShell;
+    const bool shell_model =
+        model_element_family_ == ElementFamily::kShell ||
+        model_element_family_ == ElementFamily::kThermalShell;
     for (std::size_t part_index = 0U; part_index < parts_.size();
          ++part_index) {
       const RawPart& part = parts_[part_index];
@@ -2783,10 +3136,23 @@ class DomainMappingContext::Impl {
             return;
           }
           const EntityIndex section_index =
-              static_cast<EntityIndex>(definition_.shell_sections.size());
+              model_element_family_ == ElementFamily::kThermalShell
-          definition_.shell_sections.push_back(
+                  ? static_cast<EntityIndex>(
-              {section.element_set_name, section.thickness, *material_index,
+                        definition_.thermal_sections.size())
-               section.location});
+                  : static_cast<EntityIndex>(definition_.shell_sections.size());
+          if (model_element_family_ == ElementFamily::kThermalShell) {
+            definition_.thermal_sections.push_back(
+                {{part.name, static_cast<std::int64_t>(section_index),
+                  section.element_set_name},
+                 section.element_set_name,
+                 section.thickness,
+                 *material_index,
+                 section.location});
+          } else {
+            definition_.shell_sections.push_back(
+                {section.element_set_name, section.thickness, *material_index,
+                 section.location});
+          }
           for (const auto element_label : element_set->members) {
             if (!assignments
                      .emplace(element_label,
@@ -3008,6 +3374,26 @@ class DomainMappingContext::Impl {
                 "Expanded shell section assignment must resolve exactly once.");
             return;
           }
+          if (raw_element.type == RawElement::Type::kDs4) {
+            if (definition_.thermal_elements.size() >
+                std::numeric_limits<EntityIndex>::max()) {
+              InputFailure("entity-index-overflow", raw_element.location,
+                           "ELEMENT", raw_element.label_text,
+                           "Expanded thermal element count exceeds EntityIndex "
+                           "capacity.");
+              return;
+            }
+            index =
+                static_cast<EntityIndex>(definition_.thermal_elements.size());
+            definition_.thermal_elements.push_back(
+                {{raw_instance.name, raw_element.label, raw_element.label_text},
+                 connected_nodes,
+                 assignment->second.first,
+                 raw_element.location});
+            element_indices.emplace(raw_element.label, index);
+            instance.element_mappings.push_back({raw_element.label, index});
+            continue;
+          }
           if (definition_.shell_elements.size() >
               std::numeric_limits<EntityIndex>::max()) {
             InputFailure(
@@ -3309,6 +3695,9 @@ class DomainMappingContext::Impl {
   }
 
   void ResolveElementSurfaces() {
+    if (model_element_family_ == ElementFamily::kThermalShell) {
+      return;
+    }
     std::vector<ElementSurfaceDefinition> resolved_surfaces;
     resolved_surfaces.reserve(element_surfaces_.size());
     for (const auto& surface : element_surfaces_) {
@@ -3423,6 +3812,72 @@ class DomainMappingContext::Impl {
   }
 
   void FinalizeStep() {
+    if (step_.has_heat_transfer) {
+      const SourceTargetIndex target_index = BuildSourceTargetIndex();
+      const SourceTargetResolver resolver{target_index};
+      const auto append_boundaries =
+          [this, &resolver](
+              const std::vector<RawThermalBoundary>& rows,
+              std::map<std::pair<EntityIndex, int>, double>* values,
+              std::vector<ThermalBoundaryConditionDefinition>* out) {
+            for (const auto& row : rows) {
+              const auto nodes = ResolveNodeTarget(resolver, row.target,
+                                                   row.location, "BOUNDARY");
+              if (!nodes) return false;
+              for (const EntityIndex node : *nodes) {
+                const auto key = std::make_pair(node, row.dof);
+                const auto existing = values->find(key);
+                if (existing != values->end() &&
+                    existing->second != row.value) {
+                  return InputFailure("conflicting-boundary-condition",
+                                      row.location, "BOUNDARY", row.target,
+                                      "Expanded thermal boundaries prescribe "
+                                      "different values.");
+                }
+                (*values)[key] = row.value;
+                out->push_back({node, row.dof, row.value, row.location});
+              }
+            }
+            return true;
+          };
+      std::map<std::pair<EntityIndex, int>, double> model_values;
+      if (!append_boundaries(thermal_model_boundaries_, &model_values,
+                             &definition_.thermal_model_boundaries))
+        return;
+      std::vector<ThermalBoundaryConditionDefinition> step_boundaries;
+      std::map<std::pair<EntityIndex, int>, double> step_values;
+      if (!append_boundaries(step_.thermal_boundaries, &step_values,
+                             &step_boundaries))
+        return;
+      for (const auto& row : thermal_initial_conditions_) {
+        const auto nodes = ResolveNodeTarget(resolver, row.target, row.location,
+                                             "INITIAL CONDITIONS");
+        if (!nodes) return;
+        for (const EntityIndex node : *nodes) {
+          definition_.thermal_initial_conditions.push_back(
+              {node, row.value, row.location});
+        }
+      }
+      std::vector<ThermalFluxLoadDefinition> loads;
+      for (const auto& row : step_.thermal_fluxes) {
+        const auto nodes =
+            ResolveNodeTarget(resolver, row.target, row.location, "CFLUX");
+        if (!nodes) return;
+        for (const EntityIndex node : *nodes) {
+          loads.push_back(
+              {node, row.dof, row.power, row.location, loads.size()});
+        }
+      }
+      definition_.thermal_step = {step_.name_raw,
+                                  step_.thermal_increment,
+                                  step_.thermal_period,
+                                  step_.maximum_increments,
+                                  false,
+                                  std::move(step_boundaries),
+                                  std::move(loads),
+                                  step_.location};
+      return;
+    }
     std::vector<PrescribedDisplacementDefinition> boundaries =
         model_boundaries_;
     boundaries.insert(boundaries.end(), step_.boundaries.begin(),
@@ -3525,6 +3980,8 @@ class DomainMappingContext::Impl {
   std::vector<RawElementSurface> element_surfaces_;
   std::vector<RawMaterial> materials_;
   std::vector<PrescribedDisplacementDefinition> model_boundaries_;
+  std::vector<RawThermalBoundary> thermal_model_boundaries_;
+  std::vector<RawThermalInitial> thermal_initial_conditions_;
   RawStep step_{};
   std::vector<std::map<std::int64_t, std::pair<EntityIndex, EntityIndex>>>
       part_section_assignments_;
@@ -3551,6 +4008,7 @@ class DomainMappingContext::Impl {
   bool step_no_op_seen_{false};
   bool active_output_{false};
   bool active_output_is_field_{false};
+  bool physical_constants_seen_{false};
 };
 
 DomainMappingContext::DomainMappingContext(const ParsedInput& input)

tests/unit/io/abaqus/thermal_domain_test.cpp

@@ -0,0 +1,136 @@
+#include <gtest/gtest.h>
+
+#include <filesystem>
+#include <fstream>
+#include <stdexcept>
+#include <string>
+
+#include "fesa/io/abaqus/domain_mapper.h"
+#include "fesa/io/abaqus/input_reader.h"
+
+namespace {
+
+class TemporaryInput final {
+ public:
+  explicit TemporaryInput(const std::string& content)
+      : path_{std::filesystem::temp_directory_path() /
+              "fesa-thermal-domain.inp"} {
+    std::ofstream output{path_, std::ios::binary | std::ios::trunc};
+    output << content;
+    if (!output) throw std::runtime_error{"Unable to write thermal fixture."};
+  }
+  ~TemporaryInput() {
+    std::error_code error;
+    std::filesystem::remove(path_, error);
+  }
+  const std::filesystem::path& Path() const noexcept { return path_; }
+
+ private:
+  std::filesystem::path path_;
+};
+
+fesa::Result<fesa::Domain> MapText(const std::string& text) {
+  const TemporaryInput input{text};
+  const auto parsed = fesa::AbaqusInputReader{}.Read(input.Path());
+  if (!parsed.HasValue())
+    return fesa::Result<fesa::Domain>::Failure(parsed.GetStatus());
+  return fesa::AbaqusDomainMapper{}.Map(parsed.Value());
+}
+
+std::string ThermalDeck() {
+  return R"inp(*Part, name=P
+*Node
+1, 0., 0., 0.
+2, 1., 0., 0.
+3, 1., 1., 0.
+4, 0., 1., 0.
+*Element, type=DS4
+1, 1, 2, 3, 4
+*Elset, elset=E
+1
+*Shell Section, elset=E, material=M
+.5, 5
+*End Part
+*Assembly, name=A
+*Instance, name=I, part=P
+*End Instance
+*Nset, nset=Bottom, instance=I
+1
+2
+*End Assembly
+*Material, name=M
+*Density
+7850.,
+*Conductivity
+54.,
+*Specific Heat
+460.,
+*Initial Conditions, type=TEMPERATURE
+Bottom, 20.
+*Boundary
+Bottom, 11, 11
+*Step, name=Heat, inc=100
+*Heat Transfer
+1., 100.
+*Boundary
+Bottom, 15, 15, 30.
+*Cflux
+Bottom, 11, -2.
+*Cflux
+Bottom, 11, -3.
+*End Step
+)inp";
+}
+
+TEST(ThermalDomainMapping, MapsDs4MaterialsTargetsAndAdditiveFluxes) {
+  const auto result = MapText(ThermalDeck());
+  ASSERT_TRUE(result.HasValue())
+      << (result.HasValue() ? ""
+                            : result.GetStatus().Diagnostics().front().message);
+  const auto& domain = result.Value();
+  ASSERT_EQ(domain.ThermalElements().Size(), 1U);
+  ASSERT_EQ(domain.ThermalMaterials().Size(), 1U);
+  EXPECT_DOUBLE_EQ(domain.ThermalMaterials()[0U].Conductivity(), 54.0);
+  ASSERT_NE(domain.ThermalStep(), nullptr);
+  EXPECT_EQ(domain.ThermalStep()->name, "Heat");
+  EXPECT_EQ(domain.ThermalStep()->maximum_increments, 100U);
+  EXPECT_EQ(domain.ThermalModelBoundaries().size(), 2U);
+  EXPECT_EQ(domain.ThermalStep()->Boundaries().size(), 2U);
+  ASSERT_EQ(domain.ThermalStep()->Loads().size(), 4U);
+  EXPECT_DOUBLE_EQ(domain.ThermalStep()->Loads()[0U].power, -2.0);
+  EXPECT_DOUBLE_EQ(domain.ThermalStep()->Loads()[3U].power, -3.0);
+}
+
+TEST(ThermalDomainMapping, RejectsInteriorAndConflictingThermalBoundaries) {
+  auto interior = ThermalDeck();
+  interior.replace(interior.find("Bottom, 11, 11\n"), 14U, "Bottom, 12, 12\n");
+  const auto invalid_dof = MapText(interior);
+  ASSERT_FALSE(invalid_dof.HasValue());
+  EXPECT_EQ(invalid_dof.GetStatus().Diagnostics()[0U].code, "invalid-dof");
+
+  auto conflicting = ThermalDeck();
+  conflicting.replace(conflicting.find("*Cflux"), 0U,
+                      "*Boundary\nBottom, 15, 15, 99.\n");
+  const auto invalid_conflict = MapText(conflicting);
+  ASSERT_FALSE(invalid_conflict.HasValue());
+  EXPECT_EQ(invalid_conflict.GetStatus().Diagnostics()[0U].code,
+            "conflicting-boundary-condition");
+}
+
+TEST(ThermalDomainMapping, MapsApprovedShellThermalReference) {
+  auto path = std::filesystem::path{__FILE__}.parent_path();
+  for (int index = 0; index < 4; ++index) path = path.parent_path();
+  path /= "reference/shell-thermal/shell-thermal.inp";
+  const auto parsed = fesa::AbaqusInputReader{}.Read(path);
+  ASSERT_TRUE(parsed.HasValue());
+  const auto result = fesa::AbaqusDomainMapper{}.Map(parsed.Value());
+  ASSERT_TRUE(result.HasValue())
+      << (result.HasValue() ? ""
+                            : result.GetStatus().Diagnostics().front().message);
+  EXPECT_EQ(result.Value().Nodes().size(), 729U);
+  EXPECT_EQ(result.Value().ThermalElements().Size(), 676U);
+  ASSERT_NE(result.Value().ThermalStep(), nullptr);
+  EXPECT_EQ(result.Value().ThermalStep()->maximum_increments, 100000U);
+}
+
+}  // namespace