fix: validate thermal dof layouts

src/fesa/fem/dof_manager.cpp

@@ -53,7 +53,9 @@ std::vector<DofComponent> FullNodeComponents() {
 bool IsThermalSurfaceDof(const int dof) { return dof == 11 || dof == 15; }
 
 std::size_t ThermalComponentIndex(const ThermalDofComponent component) {
-  return static_cast<std::size_t>(static_cast<int>(component) - 11);
+  const int source_dof = static_cast<int>(component);
+  if (source_dof < 11 || source_dof > 15) return kThermalDofsPerNode;
+  return static_cast<std::size_t>(source_dof - 11);
 }
 
 void AppendScatter(std::vector<std::vector<std::size_t>>& columns_by_row,
@@ -157,7 +159,8 @@ Result<DofManager> DofManager::CreateThermal(const AnalysisModel& model) {
 Status DofManager::Build(const AnalysisModel& analysis_model,
                          const ElementView& elements) {
   if (analysis_model.IsThermal()) {
-    return DofFailure("invalid-dof-profile", {analysis_model.GetDomain().SourcePath(), 0U},
+    return DofFailure(
+        "invalid-dof-profile", {analysis_model.GetDomain().SourcePath(), 0U},
         "structural", "Structural element layouts cannot build a thermal mapping.");
   }
   return Build(analysis_model, elements,
@@ -169,7 +172,8 @@ Status DofManager::Build(const AnalysisModel& analysis_model,
                          const BoundaryConditionView& boundaries) {
   std::vector<ElementDofLayout> layouts;
   if (analysis_model.IsThermal()) {
-    return DofFailure("invalid-dof-profile", {analysis_model.GetDomain().SourcePath(), 0U},
+    return DofFailure(
+        "invalid-dof-profile", {analysis_model.GetDomain().SourcePath(), 0U},
         "structural", "Structural element layouts cannot build a thermal mapping.");
   }
   layouts.reserve(elements.size());
@@ -183,13 +187,16 @@ Status DofManager::BuildThermal(
     const AnalysisModel& analysis_model,
     const std::vector<ThermalElementDofLayout>& layouts) {
   if (!analysis_model.IsThermal() || analysis_model.ThermalStep() == nullptr) {
-    return DofFailure("invalid-dof-profile", {analysis_model.GetDomain().SourcePath(), 0U},
+    return DofFailure(
+        "invalid-dof-profile", {analysis_model.GetDomain().SourcePath(), 0U},
         "thermal", "Thermal layouts require a thermal AnalysisModel.");
   }
   const Domain& domain = analysis_model.GetDomain();
   if (layouts.size() != analysis_model.ActiveThermalElements().size()) {
-    return DofFailure("invalid-element-layout-order", {domain.SourcePath(), 0U},
+    return DofFailure(
-                      std::to_string(layouts.size()), "Thermal layouts must match the active DS4 inventory.");
+        "invalid-element-layout-order", {domain.SourcePath(), 0U},
+        std::to_string(layouts.size()),
+        "Thermal layouts must match the active DS4 inventory.");
   }
   DofManager candidate;
   candidate.profile_ = DofProfile::kThermal;
@@ -198,7 +205,8 @@ Status DofManager::BuildThermal(
        ++row) {
     const EntityIndex node = analysis_model.ActiveThermalNodes()[row];
     if (node >= candidate.thermal_node_rows_.size()) {
-      return DofFailure("invalid-thermal-node", {domain.SourcePath(), 0U},
+      return DofFailure(
+          "invalid-thermal-node", {domain.SourcePath(), 0U},
           std::to_string(node), "Active thermal node is outside Domain nodes.");
     }
     candidate.thermal_node_rows_[node] = row;
@@ -209,7 +217,8 @@ Status DofManager::BuildThermal(
   for (std::size_t index = 0U; index < layouts.size(); ++index) {
     const auto& definition = domain.ThermalElements()[index];
     if (!SameSourceIdentity(layouts[index].source_id, definition.SourceId())) {
-      return DofFailure("invalid-element-layout-order", {domain.SourcePath(), 0U},
+      return DofFailure(
+          "invalid-element-layout-order", {domain.SourcePath(), 0U},
           std::to_string(index), "Thermal layout identity changed source order.");
     }
     for (std::size_t node = 0U; node < layouts[index].node_indices.size();
@@ -225,7 +234,8 @@ Status DofManager::BuildThermal(
     candidate.element_scatters_.push_back(scatter.Value());
   }
   std::map<std::size_t, double> effective;
-  const auto apply = [&candidate, &effective](
+  const auto apply =
+      [&candidate, &effective](
           const std::vector<ThermalBoundaryConditionDefinition>& boundaries,
           const bool overrides) -> Status {
     for (const auto& boundary : boundaries) {
@@ -521,6 +531,15 @@ Result<std::vector<std::size_t>> DofManager::ThermalElementScatter(
         "Thermal layouts require a thermal mapping."));
   }
   std::array<bool, kThermalDofsPerNode> components{};
+  const std::array<ThermalDofComponent, kThermalDofsPerNode> canonical{
+      {ThermalDofComponent::kBottom, ThermalDofComponent::kLower,
+       ThermalDofComponent::kMiddle, ThermalDofComponent::kUpper,
+       ThermalDofComponent::kTop}};
+  if (layout.components != canonical) {
+    return Result<std::vector<std::size_t>>::Failure(DofFailure(
+        "invalid-thermal-element-layout", {}, layout.source_id.source_label_text,
+        "Thermal layouts must preserve canonical bottom-to-top stations."));
+  }
   for (const auto component : layout.components) {
     const std::size_t index = ThermalComponentIndex(component);
     if (index >= components.size() || components[index]) {
@@ -543,12 +562,17 @@ Result<std::vector<std::size_t>> DofManager::ThermalElementScatter(
       }
     }
   }
-  if (std::adjacent_find(layout.node_indices.begin(), layout.node_indices.end()) !=
+  for (std::size_t index = 0U; index < layout.node_indices.size(); ++index) {
-      layout.node_indices.end()) {
+    if (std::find(layout.node_indices.begin(),
+                  layout.node_indices.begin() +
+                      static_cast<std::ptrdiff_t>(index),
+                  layout.node_indices[index]) !=
+        layout.node_indices.begin() + static_cast<std::ptrdiff_t>(index)) {
       return Result<std::vector<std::size_t>>::Failure(DofFailure(
           "duplicate-element-dof", {}, layout.source_id.source_label_text,
           "A thermal layout must not repeat a node."));
     }
+  }
   return Result<std::vector<std::size_t>>::Success(std::move(scatter));
 }
 

tests/unit/fem/dof_manager_test.cpp

@@ -149,18 +149,20 @@ fesa::ModelDefinition MakeThermalDefinition() {
   const std::filesystem::path source{"models/thermal-dof-manager.inp"};
   fesa::ModelDefinition definition{};
   definition.source_path = source;
-  for (std::int64_t label = 1; label <= 6; ++label) {
+  definition.nodes = {{{"Plate", 1, "1"}, {0.0, 0.0, 0.0}, {source, 1U}},
-    definition.nodes.push_back({{"Plate", label, std::to_string(label)},
+                      {{"Plate", 2, "2"}, {1.0, 0.0, 0.0}, {source, 2U}},
-                                {static_cast<double>(label), 0.0, 0.0},
+                      {{"Plate", 3, "3"}, {1.0, 1.0, 0.0}, {source, 3U}},
-                                {source, static_cast<std::size_t>(label)}});
+                      {{"Plate", 7, "7"}, {3.0, 3.0, 0.0}, {source, 7U}},
-  }
+                      {{"Plate", 4, "4"}, {0.0, 1.0, 0.0}, {source, 4U}},
+                      {{"Plate", 5, "5"}, {2.0, 0.0, 0.0}, {source, 5U}},
+                      {{"Plate", 6, "6"}, {2.0, 1.0, 0.0}, {source, 6U}}};
   definition.thermal_materials = {
       {{"", 1, "Thermal"}, "Thermal", 40.0, 7800.0, 500.0, {source, 20U}}};
   definition.thermal_sections = {
       {{"", 1, "Section"}, "Section", 0.01, 0U, {source, 21U}}};
   definition.thermal_elements = {
-      {{"Plate", 10, "10"}, {0U, 1U, 2U, 3U}, 0U, {source, 30U}},
+      {{"Plate", 10, "10"}, {0U, 1U, 2U, 4U}, 0U, {source, 30U}},
-      {{"Plate", 20, "20"}, {1U, 2U, 3U, 4U}, 0U, {source, 31U}}};
+      {{"Plate", 20, "20"}, {1U, 5U, 6U, 2U}, 0U, {source, 31U}}};
   definition.thermal_model_boundaries = {{0U, 11, 1.0, {source, 40U}},
                                          {1U, 15, 2.0, {source, 41U}}};
   definition.thermal_step = {"Heat", 0.1, 1.0, 100U, false,
@@ -445,39 +447,166 @@ TEST(ThermalDofManager, MapsOnlyConnectedNodesAndOverridesModelBoundary) {
   auto dofs = std::move(result.Value());
 
   EXPECT_EQ(dofs.Profile(), fesa::DofProfile::kThermal);
-  EXPECT_EQ(dofs.FullDofCount(), 25U);
+  EXPECT_EQ(dofs.FullDofCount(), 30U);
+  EXPECT_THROW(dofs.FullDof(0U, fesa::DofComponent::kUx), std::logic_error);
   EXPECT_EQ(dofs.ConstrainedDofs(), (std::vector<std::size_t>{0U, 9U}));
   ASSERT_EQ(dofs.PrescribedValues().Size(), 2U);
   EXPECT_DOUBLE_EQ(dofs.PrescribedValues()[0U], 3.0);
   EXPECT_DOUBLE_EQ(dofs.PrescribedValues()[1U], 2.0);
-  EXPECT_EQ(dofs.FullThermalDof(4U, fesa::ThermalDofComponent::kTop), 24U);
+  EXPECT_EQ(dofs.FullThermalDof(6U, fesa::ThermalDofComponent::kTop), 29U);
-  EXPECT_THROW(dofs.FullThermalDof(5U, fesa::ThermalDofComponent::kBottom),
+  EXPECT_THROW(dofs.FullThermalDof(3U, fesa::ThermalDofComponent::kBottom),
+               std::out_of_range);
+  EXPECT_THROW(
+      dofs.FullThermalDof(
+          0U, static_cast<fesa::ThermalDofComponent>(
+                  (std::numeric_limits<int>::min)())),
       std::out_of_range);
 
   const auto first = dofs.ThermalElementScatter(
-      {{"Plate", 10, "10"}, {0U, 1U, 2U, 3U}});
+      {{"Plate", 10, "10"}, {0U, 1U, 2U, 4U}});
   ASSERT_TRUE(first.HasValue());
   EXPECT_EQ(first.Value(),
             (std::vector<std::size_t>{0U, 1U, 2U, 3U, 4U, 5U, 6U, 7U, 8U,
                                       9U, 10U, 11U, 12U, 13U, 14U, 15U, 16U,
                                       17U, 18U, 19U}));
+  const auto second = dofs.ThermalElementScatter(
+      {{"Plate", 20, "20"}, {1U, 5U, 6U, 2U}});
+  ASSERT_TRUE(second.HasValue());
+  EXPECT_EQ(std::vector<std::size_t>(second.Value().begin(),
+                                     second.Value().begin() + 5U),
+            (std::vector<std::size_t>{5U, 6U, 7U, 8U, 9U}));
+  EXPECT_EQ(std::vector<std::size_t>(second.Value().end() - 5U,
+                                     second.Value().end()),
+            (std::vector<std::size_t>{10U, 11U, 12U, 13U, 14U}));
+  const fesa::ElementDofLayout structural_layout{
+      {"Plate", 10, "10"}, {0U}, {fesa::DofComponent::kUx}};
+  EXPECT_FALSE(dofs.ElementScatter(structural_layout).HasValue());
+  const auto repeated = dofs.ThermalElementScatter(
+      {{"Plate", 10, "10"}, {0U, 1U, 0U, 4U}});
+  ASSERT_FALSE(repeated.HasValue());
+  EXPECT_EQ(repeated.GetStatus().Diagnostics()[0U].code, "duplicate-element-dof");
+  const auto reordered = dofs.ThermalElementScatter(
+      {{"Plate", 10, "10"},
+       {0U, 1U, 2U, 4U},
+       {fesa::ThermalDofComponent::kTop, fesa::ThermalDofComponent::kLower,
+        fesa::ThermalDofComponent::kMiddle, fesa::ThermalDofComponent::kUpper,
+        fesa::ThermalDofComponent::kBottom}});
+  ASSERT_FALSE(reordered.HasValue());
+  const auto invalid_component = dofs.ThermalElementScatter(
+      {{"Plate", 10, "10"},
+       {0U, 1U, 2U, 4U},
+       {static_cast<fesa::ThermalDofComponent>((std::numeric_limits<int>::min)()),
+        fesa::ThermalDofComponent::kLower, fesa::ThermalDofComponent::kMiddle,
+        fesa::ThermalDofComponent::kUpper, fesa::ThermalDofComponent::kTop}});
+  ASSERT_FALSE(invalid_component.HasValue());
   EXPECT_EQ(RowColumns(dofs.GetSparsePattern(), 0U),
             (std::vector<std::size_t>{0U, 1U, 2U, 3U, 4U, 5U, 6U, 7U, 8U,
                                       9U, 10U, 11U, 12U, 13U, 14U, 15U, 16U,
                                       17U, 18U, 19U}));
 
+  const auto profile_before = dofs.Profile();
+  const auto full_count_before = dofs.FullDofCount();
+  const auto free_count_before = dofs.FreeDofCount();
+  const auto constrained_count_before = dofs.ConstrainedDofCount();
   const auto free_before = dofs.FreeDofs();
   const auto constrained_before = dofs.ConstrainedDofs();
+  const auto prescribed_before = dofs.PrescribedValues();
+  const auto pattern_before = dofs.GetSparsePattern();
+  const auto scatter_before = dofs.ThermalElementScatter(
+      {{"Plate", 10, "10"}, {0U, 1U, 2U, 4U}});
+  ASSERT_TRUE(scatter_before.HasValue());
   const fesa::Status invalid = dofs.BuildThermal(
-      model.Value(), {{{"Plate", 10, "10"}, {0U, 1U, 2U, 3U},
+      model.Value(), {{{"Plate", 10, "10"}, {0U, 1U, 2U, 4U},
                        {fesa::ThermalDofComponent::kBottom,
                         fesa::ThermalDofComponent::kBottom,
                         fesa::ThermalDofComponent::kMiddle,
                         fesa::ThermalDofComponent::kUpper,
                         fesa::ThermalDofComponent::kTop}},
-                      {{"Plate", 20, "20"}, {1U, 2U, 3U, 4U}}});
+                      {{"Plate", 20, "20"}, {1U, 5U, 6U, 2U}}});
   ASSERT_FALSE(invalid.IsOk());
   EXPECT_EQ(invalid.Diagnostics()[0U].code, "invalid-thermal-element-layout");
+  EXPECT_EQ(dofs.Profile(), profile_before);
+  EXPECT_EQ(dofs.FullDofCount(), full_count_before);
+  EXPECT_EQ(dofs.FreeDofCount(), free_count_before);
+  EXPECT_EQ(dofs.ConstrainedDofCount(), constrained_count_before);
   EXPECT_EQ(dofs.FreeDofs(), free_before);
   EXPECT_EQ(dofs.ConstrainedDofs(), constrained_before);
+  ASSERT_EQ(dofs.PrescribedValues().Size(), prescribed_before.Size());
+  for (std::size_t index = 0U; index < prescribed_before.Size(); ++index) {
+    EXPECT_DOUBLE_EQ(dofs.PrescribedValues()[index], prescribed_before[index]);
+  }
+  EXPECT_EQ(dofs.GetSparsePattern().row_offsets, pattern_before.row_offsets);
+  EXPECT_EQ(dofs.GetSparsePattern().column_indices,
+            pattern_before.column_indices);
+  const auto scatter_after = dofs.ThermalElementScatter(
+      {{"Plate", 10, "10"}, {0U, 1U, 2U, 4U}});
+  ASSERT_TRUE(scatter_after.HasValue());
+  EXPECT_EQ(scatter_after.Value(), scatter_before.Value());
+
+  const fesa::ElementView no_structural_elements;
+  const fesa::Status structural_misuse =
+      dofs.Build(model.Value(), no_structural_elements);
+  ASSERT_FALSE(structural_misuse.IsOk());
+  EXPECT_EQ(structural_misuse.Diagnostics()[0U].code, "invalid-dof-profile");
+  EXPECT_EQ(dofs.GetSparsePattern().row_offsets, pattern_before.row_offsets);
+}
+
+TEST(ThermalDofManager, RejectsChangedThermalLayoutInventory) {
+  auto domain = fesa::Domain::Create(MakeThermalDefinition());
+  ASSERT_TRUE(domain.HasValue());
+  auto model = fesa::AnalysisModel::Create(domain.Value());
+  ASSERT_TRUE(model.HasValue());
+  auto result = fesa::DofManager::CreateThermal(model.Value());
+  ASSERT_TRUE(result.HasValue());
+  auto dofs = std::move(result.Value());
+
+  const std::vector<fesa::ThermalElementDofLayout> layouts{
+      {{"Plate", 10, "10"}, {0U, 1U, 2U, 4U}},
+      {{"Plate", 20, "20"}, {1U, 5U, 6U, 2U}}};
+  auto changed_identity = layouts;
+  changed_identity[0U].source_id = {"Plate", 99, "99"};
+  EXPECT_FALSE(dofs.BuildThermal(model.Value(), changed_identity).IsOk());
+
+  auto changed_order = layouts;
+  std::swap(changed_order[0U], changed_order[1U]);
+  EXPECT_FALSE(dofs.BuildThermal(model.Value(), changed_order).IsOk());
+
+  auto changed_nodes = layouts;
+  changed_nodes[1U].node_indices = {1U, 6U, 5U, 2U};
+  EXPECT_FALSE(dofs.BuildThermal(model.Value(), changed_nodes).IsOk());
+}
+
+TEST(ThermalDofManager, KeepsDisconnectedThermalComponentsInSeparateCsrBlocks) {
+  auto definition = MakeThermalDefinition();
+  definition.nodes.push_back({{"Plate", 8, "8"}, {3.0, 0.0, 0.0},
+                              {definition.source_path, 8U}});
+  definition.nodes.push_back({{"Plate", 9, "9"}, {3.0, 1.0, 0.0},
+                              {definition.source_path, 9U}});
+  definition.thermal_elements[1U].node_indices = {5U, 7U, 8U, 6U};
+  definition.thermal_model_boundaries.clear();
+  definition.thermal_step->boundaries.clear();
+  auto domain = fesa::Domain::Create(std::move(definition));
+  ASSERT_TRUE(domain.HasValue());
+  auto model = fesa::AnalysisModel::Create(domain.Value());
+  ASSERT_TRUE(model.HasValue());
+  auto result = fesa::DofManager::Create(model.Value());
+  ASSERT_TRUE(result.HasValue());
+  EXPECT_EQ(result.Value().FullDofCount(), 40U);
+  EXPECT_EQ(RowColumns(result.Value().GetSparsePattern(), 0U).size(), 20U);
+  EXPECT_EQ(RowColumns(result.Value().GetSparsePattern(), 20U).front(), 20U);
+}
+
+TEST(ThermalDofManager, RejectsThermalLayoutsForStructuralModels) {
+  auto domain = fesa::Domain::Create(MakeDefinition());
+  ASSERT_TRUE(domain.HasValue());
+  auto model = fesa::AnalysisModel::Create(domain.Value());
+  ASSERT_TRUE(model.HasValue());
+  fesa::DofManager dofs;
+  const fesa::Status status = dofs.BuildThermal(model.Value(), {});
+  ASSERT_FALSE(status.IsOk());
+  EXPECT_EQ(status.Diagnostics()[0U].code, "invalid-dof-profile");
+  const auto thermal_create = fesa::DofManager::CreateThermal(model.Value());
+  ASSERT_FALSE(thermal_create.HasValue());
+  EXPECT_EQ(thermal_create.GetStatus().Diagnostics()[0U].code,
+            "invalid-dof-profile");
 }