feat: add thermal dof mapping

src/fesa/fem/dof_manager.cpp

@@ -5,6 +5,7 @@
 #include <cmath>
 #include <cstddef>
 #include <limits>
+#include <map>
 #include <stdexcept>
 #include <string>
 #include <utility>
@@ -16,6 +17,7 @@ namespace fesa {
 namespace {
 
 constexpr std::size_t kDofsPerNode = 6U;
+constexpr std::size_t kThermalDofsPerNode = 5U;
 
 Status DofFailure(const std::string& code, const SourceLocation& location,
                   const std::string& identity, const std::string& message) {
@@ -48,6 +50,12 @@ std::vector<DofComponent> FullNodeComponents() {
           DofComponent::kUrx, DofComponent::kUry, DofComponent::kUrz};
 }
 
+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);
+}
+
 void AppendScatter(std::vector<std::vector<std::size_t>>& columns_by_row,
                    const std::vector<std::size_t>& scatter) {
   for (const std::size_t row : scatter) {
@@ -94,6 +102,7 @@ std::array<std::size_t, kSize> FixedScatter(
 }  // namespace
 
 Result<DofManager> DofManager::Create(const AnalysisModel& model) {
+  if (model.IsThermal()) return CreateThermal(model);
   const Domain& domain = model.GetDomain();
   std::vector<ElementDofLayout> layouts;
   layouts.reserve(model.ActiveElements().size());
@@ -118,8 +127,39 @@ Result<DofManager> DofManager::Create(const AnalysisModel& model) {
   return Result<DofManager>::Success(std::move(dofs));
 }
 
+Result<DofManager> DofManager::CreateThermal(const AnalysisModel& model) {
+  if (!model.IsThermal()) {
+    return Result<DofManager>::Failure(DofFailure(
+        "invalid-dof-profile", {model.GetDomain().SourcePath(), 0U}, "thermal",
+        "Thermal mapping requires a thermal AnalysisModel."));
+  }
+  const Domain& domain = model.GetDomain();
+  std::vector<ThermalElementDofLayout> layouts;
+  layouts.reserve(model.ActiveThermalElements().size());
+  for (std::size_t index = 0U; index < model.ActiveThermalElements().size();
+       ++index) {
+    if (index >= domain.ThermalElements().Size()) {
+      return Result<DofManager>::Failure(DofFailure(
+          "invalid-element-layout-order", {domain.SourcePath(), 0U},
+          std::to_string(index), "Thermal element inventory is invalid."));
+    }
+    const auto& element = domain.ThermalElements()[index];
+    layouts.push_back({element.SourceId(),
+                       {element.NodeIndices()[0U], element.NodeIndices()[1U],
+                        element.NodeIndices()[2U], element.NodeIndices()[3U]}});
+  }
+  DofManager dofs;
+  const Status status = dofs.BuildThermal(model, layouts);
+  if (!status.IsOk()) return Result<DofManager>::Failure(status);
+  return Result<DofManager>::Success(std::move(dofs));
+}
+
 Status DofManager::Build(const AnalysisModel& analysis_model,
                          const ElementView& elements) {
+  if (analysis_model.IsThermal()) {
+    return DofFailure("invalid-dof-profile", {analysis_model.GetDomain().SourcePath(), 0U},
+                      "structural", "Structural element layouts cannot build a thermal mapping.");
+  }
   return Build(analysis_model, elements,
                analysis_model.Step().BoundaryConditions());
 }
@@ -128,6 +168,10 @@ Status DofManager::Build(const AnalysisModel& analysis_model,
                          const ElementView& elements,
                          const BoundaryConditionView& boundaries) {
   std::vector<ElementDofLayout> layouts;
+  if (analysis_model.IsThermal()) {
+    return DofFailure("invalid-dof-profile", {analysis_model.GetDomain().SourcePath(), 0U},
+                      "structural", "Structural element layouts cannot build a thermal mapping.");
+  }
   layouts.reserve(elements.size());
   for (const auto& element : elements) {
     layouts.push_back(element.get().DofLayout());
@@ -135,6 +179,100 @@ Status DofManager::Build(const AnalysisModel& analysis_model,
   return BuildLayouts(analysis_model, layouts, boundaries);
 }
 
+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},
+                      "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},
+                      std::to_string(layouts.size()), "Thermal layouts must match the active DS4 inventory.");
+  }
+  DofManager candidate;
+  candidate.profile_ = DofProfile::kThermal;
+  candidate.thermal_node_rows_.resize(domain.Nodes().size());
+  for (std::size_t row = 0U; row < analysis_model.ActiveThermalNodes().size();
+       ++row) {
+    const EntityIndex node = analysis_model.ActiveThermalNodes()[row];
+    if (node >= candidate.thermal_node_rows_.size()) {
+      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;
+  }
+  candidate.full_dof_count_ =
+      analysis_model.ActiveThermalNodes().size() * kThermalDofsPerNode;
+  candidate.element_scatters_.reserve(layouts.size());
+  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},
+                        std::to_string(index), "Thermal layout identity changed source order.");
+    }
+    for (std::size_t node = 0U; node < layouts[index].node_indices.size();
+         ++node) {
+      if (layouts[index].node_indices[node] != definition.NodeIndices()[node]) {
+        return DofFailure("invalid-element-layout-order",
+                          {domain.SourcePath(), 0U}, std::to_string(index),
+                          "Thermal layout nodes must preserve DS4 order.");
+      }
+    }
+    const auto scatter = candidate.ThermalElementScatter(layouts[index]);
+    if (!scatter.HasValue()) return scatter.GetStatus();
+    candidate.element_scatters_.push_back(scatter.Value());
+  }
+  std::map<std::size_t, double> effective;
+  const auto apply = [&candidate, &effective](
+                         const std::vector<ThermalBoundaryConditionDefinition>& boundaries,
+                         const bool overrides) -> Status {
+    for (const auto& boundary : boundaries) {
+      if (!IsThermalSurfaceDof(boundary.dof) || !std::isfinite(boundary.value)) {
+        return DofFailure("invalid-thermal-boundary", boundary.location, "BOUNDARY",
+                          "Thermal boundaries require finite DOF 11 or 15 values.");
+      }
+      std::size_t full_dof = 0U;
+      try {
+        full_dof = candidate.FullThermalDof(
+            boundary.node_index, static_cast<ThermalDofComponent>(boundary.dof));
+      } catch (const std::out_of_range&) {
+        return DofFailure("invalid-thermal-boundary", boundary.location, "BOUNDARY",
+                          "Thermal boundary target is inactive.");
+      }
+      if (overrides || effective.find(full_dof) == effective.end()) {
+        effective[full_dof] = boundary.value;
+      }
+    }
+    return Status::Ok();
+  };
+  Status status = apply(domain.ThermalModelBoundaries(), false);
+  if (!status.IsOk()) return status;
+  status = apply(analysis_model.ThermalStep()->Boundaries(), true);
+  if (!status.IsOk()) return status;
+  candidate.free_equations_.resize(candidate.full_dof_count_);
+  for (std::size_t full = 0U; full < candidate.full_dof_count_; ++full) {
+    const auto prescribed = effective.find(full);
+    if (prescribed != effective.end()) {
+      candidate.constrained_dofs_.push_back(full);
+    } else {
+      candidate.free_equations_[full] = candidate.free_dofs_.size();
+      candidate.free_dofs_.push_back(full);
+    }
+  }
+  candidate.prescribed_values_ = Vector{candidate.constrained_dofs_.size()};
+  for (std::size_t index = 0U; index < candidate.constrained_dofs_.size(); ++index) {
+    candidate.prescribed_values_[index] = effective[candidate.constrained_dofs_[index]];
+  }
+  candidate.sparse_pattern_ =
+      BuildSparsePattern(candidate.full_dof_count_, candidate.element_scatters_);
+  status = candidate.ValidateInvariants();
+  if (!status.IsOk()) return status;
+  *this = std::move(candidate);
+  return Status::Ok();
+}
+
 Status DofManager::BuildLayouts(const AnalysisModel& analysis_model,
                                 const std::vector<ElementDofLayout>& layouts,
                                 const BoundaryConditionView& boundaries) {
@@ -292,6 +430,8 @@ std::size_t DofManager::FullDofCount() const noexcept {
   return full_dof_count_;
 }
 
+DofProfile DofManager::Profile() const noexcept { return profile_; }
+
 std::size_t DofManager::FreeDofCount() const noexcept {
   return free_dofs_.size();
 }
@@ -302,6 +442,9 @@ std::size_t DofManager::ConstrainedDofCount() const noexcept {
 
 std::size_t DofManager::FullDof(const EntityIndex node,
                                 const DofComponent component) const {
+  if (profile_ != DofProfile::kStructural) {
+    throw std::logic_error{"Structural components are invalid for thermal DOFs."};
+  }
   const std::size_t component_index = static_cast<std::size_t>(component);
   if (node >= full_dof_count_ / kDofsPerNode ||
       component_index >= kDofsPerNode) {
@@ -310,6 +453,17 @@ std::size_t DofManager::FullDof(const EntityIndex node,
   return static_cast<std::size_t>(node) * kDofsPerNode + component_index;
 }
 
+std::size_t DofManager::FullThermalDof(
+    const EntityIndex node, const ThermalDofComponent component) const {
+  if (profile_ != DofProfile::kThermal || node >= thermal_node_rows_.size() ||
+      !thermal_node_rows_[node].has_value() ||
+      ThermalComponentIndex(component) >= kThermalDofsPerNode) {
+    throw std::out_of_range{"Thermal node or component is out of range."};
+  }
+  return *thermal_node_rows_[node] * kThermalDofsPerNode +
+         ThermalComponentIndex(component);
+}
+
 std::optional<std::size_t> DofManager::FreeEquation(
     const std::size_t full_dof) const {
   return free_equations_.at(full_dof);
@@ -317,6 +471,11 @@ std::optional<std::size_t> DofManager::FreeEquation(
 
 Result<std::vector<std::size_t>> DofManager::ElementScatter(
     const ElementDofLayout& layout) const {
+  if (profile_ != DofProfile::kStructural) {
+    return Result<std::vector<std::size_t>>::Failure(DofFailure(
+        "invalid-dof-profile", {}, layout.source_id.source_label_text,
+        "Structural layouts are invalid for a thermal mapping."));
+  }
   if (layout.node_indices.empty() || layout.components_per_node.empty() ||
       layout.node_indices.size() > (std::numeric_limits<std::size_t>::max)() /
                                        layout.components_per_node.size()) {
@@ -354,6 +513,45 @@ Result<std::vector<std::size_t>> DofManager::ElementScatter(
   return Result<std::vector<std::size_t>>::Success(std::move(scatter));
 }
 
+Result<std::vector<std::size_t>> DofManager::ThermalElementScatter(
+    const ThermalElementDofLayout& layout) const {
+  if (profile_ != DofProfile::kThermal) {
+    return Result<std::vector<std::size_t>>::Failure(DofFailure(
+        "invalid-dof-profile", {}, layout.source_id.source_label_text,
+        "Thermal layouts require a thermal mapping."));
+  }
+  std::array<bool, kThermalDofsPerNode> components{};
+  for (const auto component : layout.components) {
+    const std::size_t index = ThermalComponentIndex(component);
+    if (index >= components.size() || components[index]) {
+      return Result<std::vector<std::size_t>>::Failure(DofFailure(
+          "invalid-thermal-element-layout", {}, layout.source_id.source_label_text,
+          "Thermal layouts require each fixed station exactly once."));
+    }
+    components[index] = true;
+  }
+  std::vector<std::size_t> scatter;
+  scatter.reserve(layout.node_indices.size() * layout.components.size());
+  for (const EntityIndex node : layout.node_indices) {
+    for (const auto component : layout.components) {
+      try {
+        scatter.push_back(FullThermalDof(node, component));
+      } catch (const std::out_of_range&) {
+        return Result<std::vector<std::size_t>>::Failure(DofFailure(
+            "invalid-thermal-element-layout", {}, layout.source_id.source_label_text,
+            "Thermal layouts must reference active thermal nodes."));
+      }
+    }
+  }
+  if (std::adjacent_find(layout.node_indices.begin(), layout.node_indices.end()) !=
+      layout.node_indices.end()) {
+    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));
+}
+
 std::array<std::size_t, 12> DofManager::ElementScatter(
     const EntityIndex element) const {
   return FixedScatter<12U>(element_scatters_.at(element));