Refactor Wave classes.

- add elevation argument to virtual WaveBase functions for extracting wave velocity and acceleration (for calls where elevation was already querried).
- add non-virtual WaveBase functions for extracting wave velocity and acceleration which conditionally call GetElevation.
- add NoWaveParams and RegularWaveParams structs and corresponding constructors (for consistency with IrregularWaves).
- move utility function `get_lower_index` to irregular_wave.cpp (so that external projects need not include the HydroChrono configuration header which is included through helper.h).

Author: rserban

include/hydroc/helper.h

@@ -3,25 +3,20 @@
 
 #pragma once
 
-#include <Eigen/Dense>  // Need for the container function
 #include <fstream>
 #include <iostream>
-#include <hydroc/logging.h>
 #include <string>
 #include <filesystem>  // C++17
 
+#include <Eigen/Dense>  // Need for the container function
+
+#include <hydroc/logging.h>
+
+//// RADU - why not use constants from Chrono?
 #ifndef M_PI
     #define M_PI 3.14159265358979323846
 #endif
 
-/**@brief Returns last index of vector element below value.
- *
- * @param value Input value
- * @param ticks Array of ticks from which to find lower-bound index (assuming ascending order)
- *
- */
-size_t get_lower_index(double value, const std::vector<double>& ticks);
-
 /**@brief Base namespace for HydroChrono library
  *
  */

include/hydroc/waves/irregular_wave.h

@@ -51,8 +51,8 @@ class IrregularWaves : public WaveBase {
     void AddH5Data(std::vector<HydroData::IrregularWaveInfo>& irreg_h5_data, HydroData::SimulationParameters& sim_data);
 
     double GetElevation(const Eigen::Vector3d& position, double time) override;
-    Eigen::Vector3d GetVelocity(const Eigen::Vector3d& position, double time) override;
-    Eigen::Vector3d GetAcceleration(const Eigen::Vector3d& position, double time) override;
+    Eigen::Vector3d GetVelocity(const Eigen::Vector3d& position, double time, double elevation) override;
+    Eigen::Vector3d GetAcceleration(const Eigen::Vector3d& position, double time, double elevation) override;
 
   private:
     IrregularWaveParams params_;

include/hydroc/waves/regular_wave.h

@@ -14,25 +14,34 @@
 // Forward declaration for internal utilities (not exposed in public API)
 // Wave utilities are internal implementation details
 
+struct RegularWaveParams {
+    unsigned int num_bodies_;
+    double regular_wave_amplitude_;
+    double regular_wave_omega_;
+    double regular_wave_phase_ = 0.0;
+    bool wave_stretching_      = true;
+};
+
 class RegularWave : public WaveBase {
   public:
     RegularWave();
     explicit RegularWave(unsigned int num_b);
+    explicit RegularWave(const RegularWaveParams& params);
 
     void Initialize() override;
     Eigen::VectorXd GetForceAtTime(double t) override;
     WaveMode GetWaveMode() override { return mode_; }
 
+    //// RADU - eliminate these and use a RegularWaveParams struct (consistent with IrregularWaves)
     double regular_wave_amplitude_;
     double regular_wave_omega_;
     double regular_wave_phase_ = 0.0;
-    bool wave_stretching_      = true;
 
     void AddH5Data(std::vector<HydroData::RegularWaveInfo>& reg_h5_data, HydroData::SimulationParameters& sim_data);
 
     double GetElevation(const Eigen::Vector3d& position, double time) override;
-    Eigen::Vector3d GetVelocity(const Eigen::Vector3d& position, double time) override;
-    Eigen::Vector3d GetAcceleration(const Eigen::Vector3d& position, double time) override;
+    Eigen::Vector3d GetVelocity(const Eigen::Vector3d& position, double time, double elevation) override;
+    Eigen::Vector3d GetAcceleration(const Eigen::Vector3d& position, double time, double elevation) override;
 
   private:
     unsigned int num_bodies_;

include/hydroc/waves/wave_base.h

@@ -8,6 +8,10 @@
 
 #include <Eigen/Dense>
 
+//// RADU - why is this not in some namespace?
+////      - should these classes be DLL-exported (e.g., to allow AddWave(NoWave))?
+////        if not, why is this a public headers?
+
 enum class WaveMode {
     noWaveCIC = 0,
     regular   = 1,
@@ -18,29 +22,38 @@ class WaveBase {
   public:
     virtual ~WaveBase() = default;
 
-    virtual void Initialize()                                   = 0;
-    virtual Eigen::VectorXd GetForceAtTime(double t)            = 0;
-    virtual WaveMode         GetWaveMode()                      = 0;
-    virtual double           GetElevation(const Eigen::Vector3d& position, double time)    = 0;
-    virtual Eigen::Vector3d  GetVelocity(const Eigen::Vector3d& position, double time)     = 0;
-    virtual Eigen::Vector3d  GetAcceleration(const Eigen::Vector3d& position, double time) = 0;
+    virtual void Initialize()                                                                               = 0;
+    virtual Eigen::VectorXd GetForceAtTime(double t)                                                        = 0;
+    virtual WaveMode GetWaveMode()                                                                          = 0;
+    virtual double GetElevation(const Eigen::Vector3d& position, double time)                               = 0;
+    virtual Eigen::Vector3d GetVelocity(const Eigen::Vector3d& position, double time, double elevation)     = 0;
+    virtual Eigen::Vector3d GetAcceleration(const Eigen::Vector3d& position, double time, double elevation) = 0;
+
+    Eigen::Vector3d GetVelocity(const Eigen::Vector3d& position, double time);
+    Eigen::Vector3d GetAcceleration(const Eigen::Vector3d& position, double time);
 
-    double mwl_         = 0.0;
-    double g_           = 9.81;
-    double water_depth_ = 0.0;
+    double mwl_           = 0.0;
+    double g_             = 9.81;
+    double water_depth_   = 0.0;
+    bool wave_stretching_ = true;
+};
+
+struct NoWaveParams {
+    unsigned int num_bodies_;
 };
 
 class NoWave : public WaveBase {
   public:
     NoWave();
     explicit NoWave(unsigned int num_b);
+    explicit NoWave(const NoWaveParams& params);
 
     void Initialize() override {}
     Eigen::VectorXd GetForceAtTime(double t) override;
     WaveMode GetWaveMode() override { return mode_; }
     double GetElevation(const Eigen::Vector3d&, double) override { return 0.0; }
-    Eigen::Vector3d GetVelocity(const Eigen::Vector3d&, double) override { return Eigen::Vector3d(0.0, 0.0, 0.0); }
-    Eigen::Vector3d GetAcceleration(const Eigen::Vector3d&, double) override { return Eigen::Vector3d(0.0, 0.0, 0.0); }
+    Eigen::Vector3d GetVelocity(const Eigen::Vector3d&, double, double) override { return Eigen::Vector3d(0.0, 0.0, 0.0); }
+    Eigen::Vector3d GetAcceleration(const Eigen::Vector3d&, double, double) override { return Eigen::Vector3d(0.0, 0.0, 0.0); }
 
   private:
     unsigned int num_bodies_;

src/hydro/utils/helper.cpp

@@ -77,22 +77,6 @@ bool hydroc::GetCLIArguments(int argc,
     return true;
 }
 
-size_t get_lower_index(double value, const std::vector<double>& ticks) {
-    auto it = std::upper_bound(ticks.begin(), ticks.end(), value);
-    // get nearest-below index
-    size_t idx = it - ticks.begin() - 1;
-    // remove one if equal to value
-    if (ticks[idx] == value) {
-        idx -= 1;
-    }
-    if (idx <= 0 || idx >= ticks.size() - 1) {
-        throw std::runtime_error("Could not find index for value " + std::to_string(value) + " in array with bounds (" +
-                                 std::to_string(ticks.front()) + ", " + std::to_string(ticks.back()) + ").");
-    }
-    // return index
-    return idx;
-}
-
 using std::filesystem::path;
 
 static path DATADIR{};

src/hydro/waves/irregular_wave.cpp

@@ -18,10 +18,11 @@
 #include <sstream>
 #include <unsupported/Eigen/Splines>
 
-#include <hydroc/helper.h>
 #include <hydroc/logging.h>
 
-IrregularWaves::IrregularWaves(const IrregularWaveParams& params) : params_(params) {}
+IrregularWaves::IrregularWaves(const IrregularWaveParams& params) : params_(params) {
+    wave_stretching_ = params.wave_stretching_;
+}
 
 void IrregularWaves::InitializeIRFVectors() {
     ex_irf_sampled_.resize(params_.num_bodies_);
@@ -151,38 +152,18 @@ void IrregularWaves::AddH5Data(std::vector<HydroData::IrregularWaveInfo>& irreg_
     InitializeIRFVectors();
 }
 
-Eigen::Vector3d IrregularWaves::GetVelocity(const Eigen::Vector3d& position, double time) {
-    auto position_stretched = position;
-    if (params_.wave_stretching_) {
-        position_stretched = GetWheelerStretchedPosition(position, GetElevation(position, time), water_depth_, mwl_);
-    }
-
-    return GetWaterVelocityIrregular(position_stretched,
-                                     time,
-                                     spectrum_frequencies_,
-                                     spectral_densities_,
-                                     spectral_widths_,
-                                     wave_phases_,
-                                     wavenumbers_,
-                                     water_depth_,
-                                     mwl_);
+Eigen::Vector3d IrregularWaves::GetVelocity(const Eigen::Vector3d& position, double time, double elevation) {
+    auto position_stretched =
+        params_.wave_stretching_ ? GetWheelerStretchedPosition(position, elevation, water_depth_, mwl_) : position;
+    return GetWaterVelocityIrregular(position_stretched, time, spectrum_frequencies_, spectral_densities_,
+                                     spectral_widths_, wave_phases_, wavenumbers_, water_depth_, mwl_);
 }
 
-Eigen::Vector3d IrregularWaves::GetAcceleration(const Eigen::Vector3d& position, double time) {
-    auto position_stretched = position;
-    if (params_.wave_stretching_) {
-        position_stretched = GetWheelerStretchedPosition(position, GetElevation(position, time), water_depth_, mwl_);
-    }
-
-    return GetWaterAccelerationIrregular(position_stretched,
-                                         time,
-                                         spectrum_frequencies_,
-                                         spectral_densities_,
-                                         spectral_widths_,
-                                         wave_phases_,
-                                         wavenumbers_,
-                                         water_depth_,
-                                         mwl_);
+Eigen::Vector3d IrregularWaves::GetAcceleration(const Eigen::Vector3d& position, double time, double elevation) {
+    auto position_stretched =
+        params_.wave_stretching_ ? GetWheelerStretchedPosition(position, elevation, water_depth_, mwl_) : position;
+    return GetWaterAccelerationIrregular(position_stretched, time, spectrum_frequencies_, spectral_densities_,
+                                         spectral_widths_, wave_phases_, wavenumbers_, water_depth_, mwl_);
 }
 
 double IrregularWaves::GetElevation(const Eigen::Vector3d& position, double time) {
@@ -356,6 +337,25 @@ void IrregularWaves::CalculateWidthIRF() {
     }
 }
 
+// Return last index of vector element below value.
+// - value: Input value
+// - ticks: Array of ticks from which to find lower-bound index (assuming ascending order)
+static size_t get_lower_index(double value, const std::vector<double>& ticks) {
+    auto it = std::upper_bound(ticks.begin(), ticks.end(), value);
+    // get nearest-below index
+    size_t idx = it - ticks.begin() - 1;
+    // remove one if equal to value
+    if (ticks[idx] == value) {
+        idx -= 1;
+    }
+    if (idx <= 0 || idx >= ticks.size() - 1) {
+        throw std::runtime_error("Could not find index for value " + std::to_string(value) + " in array with bounds (" +
+                                 std::to_string(ticks.front()) + ", " + std::to_string(ticks.back()) + ").");
+    }
+    // return index
+    return idx;
+}
+
 double IrregularWaves::ExcitationConvolution(int body, int dof, double time) {
     double f_ex          = 0.0;
     auto&  irf_time_array  = ex_irf_time_sampled_[body];

src/hydro/waves/regular_wave.cpp

@@ -16,6 +16,14 @@ RegularWave::RegularWave(unsigned int num_b) {
     num_bodies_ = num_b;
 }
 
+RegularWave::RegularWave(const RegularWaveParams& params) {
+    num_bodies_ = params.num_bodies_;
+    regular_wave_amplitude_ = params.regular_wave_amplitude_;
+    regular_wave_omega_     = params.regular_wave_omega_;
+    regular_wave_phase_     = params.regular_wave_phase_;
+    wave_stretching_        = params.wave_stretching_;
+}
+
 void RegularWave::Initialize() {
     wavenumber_ = ComputeWaveNumber(regular_wave_omega_, water_depth_, g_);
 }
@@ -41,36 +49,20 @@ void RegularWave::AddH5Data(std::vector<HydroData::RegularWaveInfo>& reg_h5_data
     }
 }
 
-Eigen::Vector3d RegularWave::GetVelocity(const Eigen::Vector3d& position, double time) {
-    auto position_stretched = position;
-    if (wave_stretching_) {
-        position_stretched = GetWheelerStretchedPosition(position, GetElevation(position, time), water_depth_, mwl_);
-    }
-    return GetWaterVelocity(position_stretched,
-                            time,
-                            regular_wave_omega_,
-                            regular_wave_amplitude_,
-                            regular_wave_phase_,
-                            wavenumber_,
-                            water_depth_,
-                            mwl_);
+Eigen::Vector3d RegularWave::GetVelocity(const Eigen::Vector3d& position, double time, double elevation) {
+    auto position_stretched =
+        wave_stretching_ ? GetWheelerStretchedPosition(position, elevation, water_depth_, mwl_) : position;
+    return GetWaterVelocity(position_stretched, time, regular_wave_omega_, regular_wave_amplitude_, regular_wave_phase_,
+                            wavenumber_, water_depth_, mwl_);
 }
 
-Eigen::Vector3d RegularWave::GetAcceleration(const Eigen::Vector3d& position, double time) {
-    auto position_stretched = position;
-    if (wave_stretching_) {
-        position_stretched = GetWheelerStretchedPosition(position, GetElevation(position, time), water_depth_, mwl_);
-    }
-    return GetWaterAcceleration(position_stretched,
-                                time,
-                                regular_wave_omega_,
-                                regular_wave_amplitude_,
-                                regular_wave_phase_,
-                                wavenumber_,
-                                water_depth_,
-                                mwl_);
+Eigen::Vector3d RegularWave::GetAcceleration(const Eigen::Vector3d& position, double time, double elevation) {
+    auto position_stretched =
+        wave_stretching_ ? GetWheelerStretchedPosition(position, elevation, water_depth_, mwl_) : position;
+    return GetWaterAcceleration(position_stretched, time, regular_wave_omega_, regular_wave_amplitude_,
+                                regular_wave_phase_, wavenumber_, water_depth_, mwl_);
 }
-
+  
 double RegularWave::GetElevation(const Eigen::Vector3d& position, double time) {
     return GetEta(position, time, regular_wave_omega_, regular_wave_amplitude_, regular_wave_phase_, wavenumber_);
 }

src/hydro/waves/wave_base.cpp

@@ -5,10 +5,24 @@
 
 #include <hydroc/waves/wave_base.h>
 
+Eigen::Vector3d WaveBase::GetVelocity(const Eigen::Vector3d& position, double time) {
+    // Elevation needed only if stretching=true
+    double elevation = wave_stretching_ ? GetElevation(position, time) : 0.0;
+    return GetVelocity(position, time, elevation);
+}
+
+Eigen::Vector3d WaveBase::GetAcceleration(const Eigen::Vector3d& position, double time) {
+    // Elevation needed only if stretching=true
+    double elevation = wave_stretching_ ? GetElevation(position, time) : 0.0;
+    return GetAcceleration(position, time, elevation);
+}
+
 NoWave::NoWave() : num_bodies_(1) {}
 
 NoWave::NoWave(unsigned int num_b) : num_bodies_(num_b) {}
 
+NoWave::NoWave(const NoWaveParams& params) : num_bodies_(params.num_bodies_) {}
+
 Eigen::VectorXd NoWave::GetForceAtTime(double) {
     unsigned int dof = num_bodies_ * 6;
     Eigen::VectorXd f(dof);
@@ -17,4 +31,3 @@ Eigen::VectorXd NoWave::GetForceAtTime(double) {
     }
     return f;
 }
-