improve TestHydro docs and indexing clarity

Author: dav-og

include/hydroc/hydro_forces.h

@@ -90,6 +90,12 @@ class ChronoHydroCoupler;
 class ForceFunc6d;
 class TestHydro;
 
+// ─────────────────────────────────────────────────────────────────────────────
+// Degrees of Freedom Constant
+// ─────────────────────────────────────────────────────────────────────────────
+/// Number of degrees of freedom per rigid body: surge, sway, heave, roll, pitch, yaw.
+constexpr int kDofPerBody = 6;
+
 /**
  * @brief Per-DOF force function for Chrono's ChForce system.
  *
@@ -193,9 +199,13 @@ class ForceFunc6d {
     void ApplyForceAndTorqueToBody();
 
     std::shared_ptr<chrono::ChBody> body_;          ///< Pointer to the body this force is applied to.
-    int b_num_;                                     ///< Body's index in the system. Currently 1-indexed.
-    ComponentFunc forces_[6];                       ///< Forces for each degree of freedom.
-    std::shared_ptr<ComponentFunc> force_ptrs_[6];  ///< Pointers to the forces.
+    
+    /// Body index in the system (1-indexed, parsed from Chrono body name "bodyN").
+    /// This matches Chrono's body naming convention where bodies are named body1, body2, etc.
+    int b_num_;
+    
+    ComponentFunc forces_[kDofPerBody];                       ///< Forces for each degree of freedom.
+    std::shared_ptr<ComponentFunc> force_ptrs_[kDofPerBody];  ///< Pointers to the forces.
     std::shared_ptr<chrono::ChForce> chrono_force_;  ///< Chrono force for the body.
     std::shared_ptr<chrono::ChForce> chrono_torque_; ///< Chrono torque for the body.
     TestHydro* all_hydro_forces_;                   ///< Pointer to TestHydro for calculations.
@@ -214,8 +224,65 @@ struct HydroProfileStats {
     int waves_calls             = 0;
 };
 
-// TODO: Rename TestHydro for clarity, perhaps to HydroForces?
-// TODO: Split TestHydro class from its helper classes for clearer code structure.
+// ═══════════════════════════════════════════════════════════════════════════════
+//
+//  TestHydro — Primary Hydrodynamics Façade / Adapter
+//
+// ═══════════════════════════════════════════════════════════════════════════════
+/**
+ * @brief Primary hydrodynamics façade for HydroChrono applications.
+ *
+ * TestHydro is the main user-facing object for attaching hydrodynamic forces to
+ * Chrono bodies. Despite its name (a historical artifact), it is NOT a test-only
+ * class — it is the production hydrodynamics adapter used by all workflows.
+ *
+ * @note The name "TestHydro" is retained for backward compatibility. New code
+ *       may use the alias `HydroForces` for clarity (see below).
+ *
+ * ## What It Does
+ *
+ * TestHydro bridges three layers:
+ *   1. **Chrono bodies** — The physical bodies in the Chrono simulation.
+ *   2. **HydroSystem (core)** — Chrono-free hydrodynamic force computation.
+ *   3. **Chrono force callbacks** — ChForce functions that Chrono calls each timestep.
+ *
+ * ## Typical Usage
+ *
+ * @code{.cpp}
+ *   // Create Chrono bodies
+ *   auto body1 = chrono_types::make_shared<ChBody>();
+ *   body1->SetName("body1");
+ *   system.Add(body1);
+ *
+ *   // Attach hydrodynamic forces (with or without waves)
+ *   std::vector<std::shared_ptr<ChBody>> bodies = {body1};
+ *   TestHydro hydro(bodies, "path/to/hydro_data.h5");
+ *   hydro.AddWaves(std::make_shared<RegularWave>(...));
+ *
+ *   // Run simulation — Chrono automatically queries TestHydro for forces
+ * @endcode
+ *
+ * ## Used By
+ *
+ *   - **C++ tests and demos** — Directly instantiate TestHydro with bodies and H5 file.
+ *   - **YAML runner (hydrochrono.exe)** — SetupHydroFromYAML creates TestHydro internally.
+ *
+ * ## Architecture Notes
+ *
+ *   - Internally delegates to HydroSystem + ChronoHydroCoupler for force computation.
+ *   - Maintains velocity history for radiation damping convolution.
+ *   - Caches forces per timestep to avoid redundant computation.
+ *
+ * ## Body Indexing Convention
+ *
+ *   Chrono bodies must be named "body1", "body2", etc. The numeric suffix is parsed
+ *   to determine body order. Internally, body indices are **1-based** in the callback
+ *   interface (CoordinateFuncForBody) to match this naming convention.
+ *
+ * @see HydroForces (alias below)
+ * @see HydroSystem (Chrono-free hydrodynamic core)
+ * @see ChronoHydroCoupler (extracts state from Chrono bodies)
+ */
 class TestHydro {
   public:
     TestHydro() = delete;
@@ -328,13 +395,20 @@ class TestHydro {
     /**
      * @brief Calculates or retrieves the total force on a specific body in a particular degree of freedom.
      *
-     * If the total force for the body and DOF was computed for the current timestep, it's retrieved.
-     * Otherwise, the function calculates it. Note: Body index is 1-based here due to its origin from ForceFunc6d.
+     * This is the main entry point called by Chrono's ChForce callbacks during simulation.
+     * Forces are computed once per timestep and cached; subsequent calls within the same
+     * timestep return the cached value.
+     *
+     * @note **1-Based Body Indexing**: The body index `b` is 1-based (body1 → b=1, body2 → b=2, etc.)
+     *       to match Chrono's body naming convention. This is intentional and matches how
+     *       ForceFunc6d parses body names like "body1", "body2", etc.
      *
-     * @param b Body index (1-based due to source from ForceFunc6d).
-     * @param i Degree of Freedom (DOF) index, ranging from [0,...5].
+     * @param b Body index (1-based: valid range is [1, num_bodies]).
+     * @param i Degree of Freedom (DOF) index (0-based: 0=surge, 1=sway, 2=heave, 3=roll, 4=pitch, 5=yaw).
      *
-     * @return Component of the force vector for body 'b' and DOF 'i'.
+     * @return Force component for body 'b' in DOF 'i' (Newtons for forces, Newton-meters for torques).
+     *
+     * @throws std::out_of_range if b or i are out of valid range.
      */
     double CoordinateFuncForBody(int b, int i);
 
@@ -459,4 +533,21 @@ class TestHydro {
     void EnsureHydroSystemAndCoupler();
 };
 
+// ─────────────────────────────────────────────────────────────────────────────
+// Public Alias for Clarity
+// ─────────────────────────────────────────────────────────────────────────────
+/**
+ * @brief Alias for TestHydro — the primary hydrodynamics façade.
+ *
+ * HydroForces is the recommended name for new code. The underlying class
+ * is TestHydro (historical name retained for backward compatibility).
+ *
+ * Usage:
+ * @code{.cpp}
+ *   HydroForces hydro(bodies, "hydro_data.h5");
+ *   hydro.AddWaves(my_waves);
+ * @endcode
+ */
+using HydroForces = TestHydro;
+
 #endif

src/hydro/chrono/chrono_hydro_coupler.cpp

@@ -6,7 +6,7 @@
 #include <hydroc/coupling/chrono_coupler.h>
 #include "chrono_state_utils.h"
 
-#include <cassert>
+#include <stdexcept>
 
 namespace hydrochrono::hydro {
 
@@ -15,8 +15,17 @@ ChronoHydroCoupler::ChronoHydroCoupler(
     std::vector<std::shared_ptr<chrono::ChBody>> bodies)
     : hydro_system_(hydro_system),
       bodies_(std::move(bodies)) {
-    assert(hydro_system_ != nullptr);
-    assert(!bodies_.empty());
+    // HydroSystem is required for force evaluation.
+    if (hydro_system_ == nullptr) {
+        throw std::invalid_argument(
+            "ChronoHydroCoupler: hydro_system pointer must not be null");
+    }
+
+    // At least one body is required for hydrodynamic coupling.
+    if (bodies_.empty()) {
+        throw std::invalid_argument(
+            "ChronoHydroCoupler: bodies vector must not be empty");
+    }
 }
 
 BodyForces ChronoHydroCoupler::Evaluate(double time) {

src/hydro/core/hydro_system.cpp

@@ -5,8 +5,8 @@
 
 #include <hydroc/core/hydro_system.h>
 
-#include <cassert>
 #include <chrono>
+#include <stdexcept>
 #include <Eigen/Dense>
 
 namespace hydrochrono::hydro {
@@ -17,7 +17,11 @@ HydroSystem::HydroSystem(int num_bodies,
       components_(std::move(components)),
       profile_stats_(),
       profiling_enabled_(false) {
-    assert(num_bodies_ > 0);
+    // Require at least one body; otherwise hydrodynamic forces are meaningless.
+    if (num_bodies_ <= 0) {
+        throw std::invalid_argument(
+            "HydroSystem: num_bodies must be > 0 (got " + std::to_string(num_bodies_) + ")");
+    }
 }
 
 BodyForces HydroSystem::Evaluate(const SystemState& state, double time) {

src/hydro/force_components/excitation_component.cpp

@@ -5,21 +5,35 @@
 
 #include "excitation_component.h"
 
-#include <cassert>
+#include <stdexcept>
 #include <Eigen/Dense>
 
 namespace hydrochrono::hydro {
 
 ExcitationComponent::ExcitationComponent(std::shared_ptr<WaveBase> waves, int num_bodies)
     : waves_(waves), num_bodies_(num_bodies) {
-    assert(waves_ != nullptr);
-    assert(num_bodies_ > 0);
+    // Wave model is required for excitation force computation.
+    if (waves_ == nullptr) {
+        throw std::invalid_argument(
+            "ExcitationComponent: waves pointer must not be null");
+    }
+
+    // Require at least one body for wave excitation computation.
+    if (num_bodies_ <= 0) {
+        throw std::invalid_argument(
+            "ExcitationComponent: num_bodies must be > 0 (got " + std::to_string(num_bodies_) + ")");
+    }
 }
 
 void ExcitationComponent::Compute(const SystemState& state,
                               double time,
                               BodyForces& inout_forces) {
-    assert(static_cast<int>(inout_forces.size()) == num_bodies_);
+    // Internal consistency check: forces must match expected body count.
+    if (static_cast<int>(inout_forces.size()) != num_bodies_) {
+        throw std::runtime_error(
+            "ExcitationComponent::Compute: inout_forces.size() mismatch (expected " +
+            std::to_string(num_bodies_) + ", got " + std::to_string(inout_forces.size()) + ")");
+    }
 
     // Get wave excitation forces from WaveBase (returns flat 6N vector)
     Eigen::VectorXd force_flat = waves_->GetForceAtTime(time);

src/hydro/force_components/hydrostatics_component.cpp

@@ -7,7 +7,7 @@
 #include <hydroc/io/h5_reader.h>
 
 #include <Eigen/Dense>
-#include <cassert>
+#include <stdexcept>
 
 namespace hydrochrono::hydro {
 
@@ -22,14 +22,46 @@ HydrostaticsComponent::HydrostaticsComponent(
       equilibrium_(equilibrium),
       cb_minus_cg_(cb_minus_cg),
       gravitational_acceleration_(gravitational_acceleration) {
-    assert(num_bodies_ > 0);
+    // Require at least one body for hydrostatic force computation.
+    if (num_bodies_ <= 0) {
+        throw std::invalid_argument(
+            "HydrostaticsComponent: num_bodies must be > 0 (got " + std::to_string(num_bodies_) + ")");
+    }
+
+    // Validate equilibrium array size: must have 6 DOF per body.
+    const int expected_equilibrium_size = kDofPerBody * num_bodies_;
+    if (static_cast<int>(equilibrium_.size()) != expected_equilibrium_size) {
+        throw std::invalid_argument(
+            "HydrostaticsComponent: equilibrium array size mismatch (expected " +
+            std::to_string(expected_equilibrium_size) + ", got " +
+            std::to_string(equilibrium_.size()) + ")");
+    }
+
+    // Validate cb_minus_cg array size: must have 3 components per body.
+    const int expected_cb_cg_size = kDofLinOrRot * num_bodies_;
+    if (static_cast<int>(cb_minus_cg_.size()) != expected_cb_cg_size) {
+        throw std::invalid_argument(
+            "HydrostaticsComponent: cb_minus_cg array size mismatch (expected " +
+            std::to_string(expected_cb_cg_size) + ", got " +
+            std::to_string(cb_minus_cg_.size()) + ")");
+    }
 }
 
 void HydrostaticsComponent::Compute(const SystemState& state,
                                 double time,
                                 BodyForces& inout_forces) {
-    assert(static_cast<int>(state.bodies.size()) == num_bodies_);
-    assert(static_cast<int>(inout_forces.size()) == num_bodies_);
+    // Internal consistency check: state and forces must match expected body count.
+    // This can fail if HydroSystem is misused (integration bug, not user config error).
+    if (static_cast<int>(state.bodies.size()) != num_bodies_) {
+        throw std::runtime_error(
+            "HydrostaticsComponent::Compute: state.bodies.size() mismatch (expected " +
+            std::to_string(num_bodies_) + ", got " + std::to_string(state.bodies.size()) + ")");
+    }
+    if (static_cast<int>(inout_forces.size()) != num_bodies_) {
+        throw std::runtime_error(
+            "HydrostaticsComponent::Compute: inout_forces.size() mismatch (expected " +
+            std::to_string(num_bodies_) + ", got " + std::to_string(inout_forces.size()) + ")");
+    }
 
     const double rho = file_info_.GetRhoVal();
     const double rho_times_g = rho * gravitational_acceleration_.norm();

src/hydro/force_components/radiation_component.cpp

@@ -8,7 +8,6 @@
 #include <hydroc/io/h5_reader.h>
 
 #include <Eigen/Dense>
-#include <cassert>
 #include <algorithm>
 #include <stdexcept>
 
@@ -33,7 +32,30 @@ RadiationComponent::RadiationComponent(
       rirf_time_vector_(rirf_time_vector),
       rirf_width_vector_(rirf_width_vector),
       rirf_processed_ready_(false) {
-    assert(num_bodies_ > 0);
+    // Require at least one body for radiation damping computation.
+    if (num_bodies_ <= 0) {
+        throw std::invalid_argument(
+            "RadiationComponent: num_bodies must be > 0 (got " + std::to_string(num_bodies_) + ")");
+    }
+
+    // Validate RIRF time series data from HDF5.
+    if (rirf_steps <= 0) {
+        throw std::invalid_argument(
+            "RadiationComponent: rirf_steps must be > 0 (got " + std::to_string(rirf_steps) +
+            "). Check that HDF5 file contains valid RIRF data.");
+    }
+
+    // Validate RIRF time and width vectors match the expected size.
+    if (static_cast<int>(rirf_time_vector_.size()) != rirf_steps) {
+        throw std::invalid_argument(
+            "RadiationComponent: rirf_time_vector size mismatch (expected " +
+            std::to_string(rirf_steps) + ", got " + std::to_string(rirf_time_vector_.size()) + ")");
+    }
+    if (static_cast<int>(rirf_width_vector_.size()) != rirf_steps) {
+        throw std::invalid_argument(
+            "RadiationComponent: rirf_width_vector size mismatch (expected " +
+            std::to_string(rirf_steps) + ", got " + std::to_string(rirf_width_vector_.size()) + ")");
+    }
 }
 
 void RadiationComponent::EnsureProcessedRIRF() {
@@ -78,13 +100,17 @@ double RadiationComponent::GetRIRFval(int row, int col, int st) {
 void RadiationComponent::Compute(const SystemState& state,
                              double time,
                              BodyForces& inout_forces) {
-    assert(static_cast<int>(state.bodies.size()) == num_bodies_);
-    assert(static_cast<int>(inout_forces.size()) == num_bodies_);
-
-    const int rirf_steps = file_info_.GetRIRFDims(2);
-    const int total_dofs = kDofPerBody * num_bodies_;
-
-    assert(total_dofs > 0 && rirf_steps > 0);
+    // Internal consistency check: state and forces must match expected body count.
+    if (static_cast<int>(state.bodies.size()) != num_bodies_) {
+        throw std::runtime_error(
+            "RadiationComponent::Compute: state.bodies.size() mismatch (expected " +
+            std::to_string(num_bodies_) + ", got " + std::to_string(state.bodies.size()) + ")");
+    }
+    if (static_cast<int>(inout_forces.size()) != num_bodies_) {
+        throw std::runtime_error(
+            "RadiationComponent::Compute: inout_forces.size() mismatch (expected " +
+            std::to_string(num_bodies_) + ", got " + std::to_string(inout_forces.size()) + ")");
+    }
 
     // If using TaperedDirect, ensure processed kernel is ready
     if (convolution_mode_ == RadiationConvolutionMode::TaperedDirect) {