feat(solid): add Hex8 input and three-DOF semantics

src/fesa/model/domain.cpp

@@ -1,10 +1,13 @@
 #include "fesa/model/domain.h"
 
-#include <cstddef>
+#include <algorithm>
+#include <array>
 #include <cmath>
+#include <cstddef>
 #include <functional>
 #include <map>
 #include <memory>
+#include <set>
 #include <string>
 #include <utility>
 
@@ -48,12 +51,14 @@ Status ValidateThermalTargets(
 
 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();
+                              !definition.shell_elements.empty() ||
+                              !definition.solid_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.");
@@ -161,6 +166,180 @@ Status ValidateThermalDefinition(const ModelDefinition& definition) {
   return Status::Ok();
 }
 
+/// @brief Checks Solid semantic candidates before immutable ownership starts.
+Status ValidateSolidDefinition(ModelDefinition& definition) {
+  const auto fail = [](const SourceLocation& location,
+                       const std::string& identity,
+                       const std::string& message) {
+    return Status::Failure(FailureCategory::kModel,
+                           {{Severity::kError, "invalid-solid-definition",
+                             location, "SOLID", identity, message}});
+  };
+  if (definition.solid_elements.empty()) {
+    if (!definition.solid_sections.empty() || !definition.amplitudes.empty())
+      return fail({}, "C3D8", "Solid definitions require a C3D8 model.");
+    return Status::Ok();
+  }
+  if (!definition.elements.empty() || !definition.shell_elements.empty() ||
+      !definition.thermal_elements.empty() || definition.thermal_step ||
+      !definition.sections.empty() || !definition.shell_sections.empty() ||
+      !definition.element_surfaces.empty()) {
+    return fail({}, "C3D8",
+                "Solid models cannot mix element or section families.");
+  }
+  if (definition.steps.size() + (definition.frequency_step ? 1U : 0U) +
+          (definition.dynamic_step ? 1U : 0U) !=
+      1U) {
+    return fail({}, "STEP", "Solid models require exactly one procedure.");
+  }
+  std::vector<bool> active(definition.nodes.size(), false);
+  std::set<std::pair<std::string, std::int64_t>> node_ids;
+  for (const auto& node : definition.nodes) {
+    if (!node_ids
+             .emplace(node.source_id.instance_name, node.source_id.source_label)
+             .second)
+      return fail(node.location, node.source_id.source_label_text,
+                  "Node identity is duplicated.");
+    for (const double coordinate : node.coordinates)
+      if (!std::isfinite(coordinate))
+        return fail(node.location, node.source_id.source_label_text,
+                    "Solid coordinates must be finite.");
+  }
+  std::set<std::pair<std::string, std::int64_t>> element_ids;
+  std::set<EntityIndex> active_materials;
+  for (const auto& element : definition.solid_elements) {
+    if (!element_ids
+             .emplace(element.source_id.instance_name,
+                      element.source_id.source_label)
+             .second ||
+        element.section_index >= definition.solid_sections.size() ||
+        element.material_index >= definition.materials.size()) {
+      return fail(
+          element.location, element.source_id.source_label_text,
+          "Solid identity, section, and material must resolve uniquely.");
+    }
+    if (definition.solid_sections[element.section_index].MaterialIndex() !=
+        element.material_index)
+      return fail(element.location, element.source_id.source_label_text,
+                  "Solid material assignment conflicts.");
+    active_materials.insert(element.material_index);
+    std::set<EntityIndex> distinct;
+    for (const auto node : element.node_indices) {
+      if (node >= active.size() || !distinct.insert(node).second)
+        return fail(element.location, element.source_id.source_label_text,
+                    "Solid connectivity requires eight distinct active nodes.");
+      active[node] = true;
+    }
+  }
+  for (std::size_t node = 0U; node < active.size(); ++node)
+    if (!active[node])
+      return fail(definition.nodes[node].location,
+                  definition.nodes[node].source_id.source_label_text,
+                  "Solid nodes must be connected.");
+  for (const auto& section : definition.solid_sections)
+    if (section.MaterialIndex() >= definition.materials.size() ||
+        section.Name().empty())
+      return fail(section.Location(), section.Name(),
+                  "Solid material assignment is invalid.");
+  for (const auto index : active_materials) {
+    const auto& material = definition.materials[index];
+    const double young = material.YoungsModulus();
+    const double nu = material.PoissonsRatio();
+    const double mu = young / (2.0 * (1.0 + nu));
+    const double lambda = young * nu / ((1.0 + nu) * (1.0 - 2.0 * nu));
+    if (!std::isfinite(young) || !(young > 0.0) || !std::isfinite(nu) ||
+        !(nu > -1.0 && nu < 0.5) || !std::isfinite(mu) ||
+        !std::isfinite(lambda))
+      return fail(material.Location(), material.Name(),
+                  "Solid isotropic constants are invalid.");
+  }
+  std::set<EntityIndex> density_materials;
+  for (const auto& density : definition.densities) {
+    if (density.material_index >= definition.materials.size() ||
+        !std::isfinite(density.rho) || !(density.rho > 0.0) ||
+        !density_materials.insert(density.material_index).second)
+      return fail(density.location, "DENSITY",
+                  "Solid density must be finite positive and unique.");
+  }
+  if (definition.dynamic_step || definition.frequency_step)
+    for (const auto material : active_materials)
+      if (density_materials.count(material) == 0U)
+        return fail({}, "DENSITY",
+                    "Solid frequency and dynamic materials require density.");
+  std::set<std::string> amplitude_names;
+  for (const auto& amplitude : definition.amplitudes) {
+    if (!definition.dynamic_step || amplitude.name_raw.empty() ||
+        amplitude.name_canonical.empty() ||
+        !amplitude_names.insert(amplitude.name_canonical).second ||
+        amplitude.points.empty())
+      return fail(amplitude.location, amplitude.name_raw,
+                  "Solid amplitude identity is invalid.");
+    double previous = -1.0;
+    for (const auto& point : amplitude.points) {
+      if (!std::isfinite(point[0]) || !std::isfinite(point[1]) ||
+          point[0] < 0.0 || !(point[0] > previous))
+        return fail(amplitude.location, amplitude.name_raw,
+                    "Solid amplitude pairs are invalid.");
+      previous = point[0];
+    }
+  }
+  const auto validate_boundaries = [&fail](const auto& boundaries,
+                                           bool homogeneous) {
+    for (const auto& boundary : boundaries) {
+      if (boundary.first_dof < 1 || boundary.last_dof < boundary.first_dof ||
+          boundary.last_dof > 3 || !std::isfinite(boundary.value) ||
+          (homogeneous && boundary.value != 0.0))
+        return fail(boundary.location, boundary.target,
+                    "Solid physical boundaries require DOF1..3.");
+    }
+    return Status::Ok();
+  };
+  const auto validate_loads = [&fail, &definition](const auto& loads,
+                                                   bool dynamic) {
+    for (const auto& load : loads) {
+      if (load.dof < 1 || load.dof > 3 || !std::isfinite(load.magnitude) ||
+          (load.amplitude_index &&
+           (!dynamic || *load.amplitude_index >= definition.amplitudes.size())))
+        return fail(load.location, load.target,
+                    "Solid load component or amplitude is invalid.");
+    }
+    return Status::Ok();
+  };
+  if (definition.frequency_step)
+    return validate_boundaries(definition.frequency_step->boundaries, true);
+  if (definition.dynamic_step) {
+    const auto& dynamic = *definition.dynamic_step;
+    if (!std::isfinite(dynamic.initial_increment) ||
+        !(dynamic.initial_increment > 0.0) ||
+        !std::isfinite(dynamic.time_period) || !(dynamic.time_period > 0.0) ||
+        dynamic.increment_count == 0U ||
+        dynamic.increment_count > dynamic.maximum_increments ||
+        dynamic.initial_increment *
+                static_cast<double>(dynamic.increment_count) !=
+            dynamic.time_period ||
+        dynamic.alpha != -0.05 || dynamic.beta != 0.275625 ||
+        dynamic.gamma != 0.55)
+      return fail(dynamic.data_location, dynamic.source_step_name_raw,
+                  "Solid dynamic controls are invalid.");
+    const auto status = validate_boundaries(dynamic.boundaries, true);
+    if (!status.IsOk()) return status;
+    return validate_loads(dynamic.loads, true);
+  }
+  const auto& step = definition.steps[0];
+  if (step.direct &&
+      (!std::isfinite(step.initial_increment) ||
+       !(step.initial_increment > 0.0) || !std::isfinite(step.time_period) ||
+       !(step.time_period > 0.0) || step.increment_count == 0U ||
+       step.increment_count > step.maximum_increments ||
+       step.initial_increment * static_cast<double>(step.increment_count) !=
+           step.time_period))
+    return fail(step.location, step.name,
+                "Solid static fixed schedule is invalid.");
+  const auto status = validate_boundaries(step.boundaries, false);
+  if (!status.IsOk()) return status;
+  return validate_loads(step.loads, false);
+}
+
 }  // namespace
 
 StepDefinition::StepDefinition(StaticStepDefinition definition)
@@ -192,12 +371,17 @@ StepDefinition::StepDefinition(StaticStepDefinition definition)
     auto& load = definition.loads[source_order];
     auto owned = std::make_unique<ConcentratedNodalLoad>(
         SourceTargetQuery{SourceEntityKind::kNode, "", load.target}, load.dof,
-        load.magnitude, source_order, std::move(load.location));
+        load.magnitude, source_order, std::move(load.location),
+        load.amplitude_index, std::move(load.amplitude_name_raw),
+        std::move(load.amplitude_location));
     loads_view_.push_back(std::cref(*owned));
     concentrated_loads_view_.Add(*owned);
     loads_.push_back(std::move(owned));
   }
   surface_pressure_loads_ = std::move(definition.surface_pressure_loads);
+  direct_static_ = definition.direct;
+  increment_count_ = definition.increment_count;
+  maximum_increments_ = definition.maximum_increments;
 }
 
 StepDefinition::StepDefinition(FrequencyStepDefinition definition)
@@ -227,8 +411,10 @@ StepDefinition::StepDefinition(FrequencyStepDefinition definition)
 
 StepDefinition::StepDefinition(DynamicStepDefinition definition)
     : name_{definition.source_step_name_raw},
-      procedure_kind_{StaticProcedureKind::kNonlinearFullNewton},
-      nlgeom_{true},
+      procedure_kind_{definition.nlgeom
+                          ? StaticProcedureKind::kNonlinearFullNewton
+                          : StaticProcedureKind::kLinearSmallDisplacement},
+      nlgeom_{definition.nlgeom},
       initial_increment_{definition.initial_increment},
       time_period_{definition.time_period},
       minimum_increment_{definition.initial_increment},
@@ -255,7 +441,9 @@ StepDefinition::StepDefinition(DynamicStepDefinition definition)
     auto& load = dynamic_->loads[source_order];
     auto owned = std::make_unique<ConcentratedNodalLoad>(
         SourceTargetQuery{SourceEntityKind::kNode, "", load.target}, load.dof,
-        load.magnitude, source_order, std::move(load.location));
+        load.magnitude, source_order, std::move(load.location),
+        load.amplitude_index, std::move(load.amplitude_name_raw),
+        std::move(load.amplitude_location));
     loads_view_.push_back(std::cref(*owned));
     concentrated_loads_view_.Add(*owned);
     loads_.push_back(std::move(owned));
@@ -325,11 +513,21 @@ double StepDefinition::MaximumIncrement() const noexcept {
   return maximum_increment_;
 }
 
+bool StepDefinition::IsDirectStatic() const noexcept { return direct_static_; }
+std::uint32_t StepDefinition::IncrementCount() const noexcept {
+  return increment_count_;
+}
+std::uint32_t StepDefinition::MaximumIncrements() const noexcept {
+  return maximum_increments_;
+}
+
 const SourceLocation& StepDefinition::Location() const noexcept {
   return location_;
 }
 
 Result<Domain> Domain::Create(ModelDefinition definition) {
+  const Status solid_status = ValidateSolidDefinition(definition);
+  if (!solid_status.IsOk()) return Result<Domain>::Failure(solid_status);
   const Status thermal_status = ValidateThermalDefinition(definition);
   if (!thermal_status.IsOk()) return Result<Domain>::Failure(thermal_status);
   return Result<Domain>::Success(Domain{std::move(definition)});
@@ -354,6 +552,23 @@ const DomainCollectionView<Mitc4ShellDefinition>& Domain::ShellElements()
   return shell_elements_view_;
 }
 
+const DomainCollectionView<Hex8Definition>& Domain::SolidElements()
+    const noexcept {
+  return solid_elements_view_;
+}
+const DomainCollectionView<SolidSection>& Domain::SolidSections()
+    const noexcept {
+  return solid_sections_view_;
+}
+bool Domain::IsSolid() const noexcept {
+  return !solid_elements_view_.Empty() &&
+         solid_elements_view_.Size() == elements_view_.Size();
+}
+const std::vector<TabularAmplitudeDefinition>& Domain::Amplitudes()
+    const noexcept {
+  return definition_.amplitudes;
+}
+
 const DomainCollectionView<Material>& Domain::Materials() const noexcept {
   return materials_view_;
 }
@@ -496,6 +711,14 @@ Domain::Domain(ModelDefinition definition)
     element_properties_.push_back(std::move(owned));
   }
   definition_.shell_sections.clear();
+  const std::size_t solid_property_offset = element_properties_.size();
+  for (auto& section : definition_.solid_sections) {
+    auto owned = std::make_unique<SolidSection>(std::move(section));
+    properties_view_.Add(*owned);
+    solid_sections_view_.Add(*owned);
+    element_properties_.push_back(std::move(owned));
+  }
+  definition_.solid_sections.clear();
 
   element_definitions_.reserve(definition_.elements.size() +
                                definition_.shell_elements.size());
@@ -517,6 +740,16 @@ Domain::Domain(ModelDefinition definition)
     element_definitions_.push_back(std::move(owned));
   }
   definition_.shell_elements.clear();
+  for (auto& element : definition_.solid_elements) {
+    element.SynchronizeNodeIndices();
+    element.SetPropertyIndex(static_cast<EntityIndex>(solid_property_offset +
+                                                      element.section_index));
+    auto owned = std::make_unique<Hex8Definition>(std::move(element));
+    elements_view_.Add(*owned);
+    solid_elements_view_.Add(*owned);
+    element_definitions_.push_back(std::move(owned));
+  }
+  definition_.solid_elements.clear();
   for (auto& element : definition_.thermal_elements) {
     const EntityIndex material_index = static_cast<EntityIndex>(
         structural_material_count +