Speed up irregular wave animation in VSG

Author: dav-og

include/hydroc/waves/irregular_wave.h

@@ -52,6 +52,21 @@ class IrregularWaves : public WaveBase {
     Eigen::Vector3d GetAcceleration(const Eigen::Vector3d& position, double time) override;
     Eigen::Vector2d GetElevationGradientXY(const Eigen::Vector3d& position, double time) const override;
 
+    /// @brief Compute elevation using only the first N frequency components (for visualization).
+    /// 
+    /// This method is optimized for real-time visualization where visual quality can be
+    /// traded for performance. Using fewer components (e.g., 50-100 instead of 1000)
+    /// provides a significant speedup while maintaining a visually acceptable wave surface.
+    /// 
+    /// @param position  World coordinates [m]. Only x-coordinate is used (waves propagate in +X).
+    /// @param time      Simulation time [s].
+    /// @param max_components  Maximum number of frequency components to use.
+    ///                        If <= 0 or >= num_components, all components are used.
+    /// @return Free-surface elevation [m] above mean water level.
+    /// 
+    /// @note Physics calculations should use GetElevation() which uses all components.
+    double GetElevationForVisualization(const Eigen::Vector3d& position, double time, int max_components) const;
+
   private:
     IrregularWaveParams params_;
     std::vector<double> spectrum_;

src/gui/vsg_gui_component.cpp

@@ -1,10 +1,22 @@
 // =============================================================================
 // HydroChrono VSG GUI Component - Implementation
 // =============================================================================
+
+// Suppress MSVC warning about IMGUI_API macro redefinition in vsgImGui headers.
+// This is a harmless conflict between vsgImGui/imgui.h and vsgImGui/Export.h.
+#ifdef _MSC_VER
+#pragma warning(push)
+#pragma warning(disable : 4005)  // macro redefinition
+#endif
+
 #include "vsg_gui_component.h"
 
 #include <vsgImGui/imgui.h>
 
+#ifdef _MSC_VER
+#pragma warning(pop)
+#endif
+
 namespace hydroc {
 namespace gui {
 

src/gui/vsg_water_surface.cpp

@@ -13,6 +13,7 @@
 
 #include <chrono/core/ChTypes.h>
 #include <chrono/physics/ChBodyEasy.h>
+#include <hydroc/waves/irregular_wave.h>
 
 #include <Eigen/Dense>
 
@@ -692,13 +693,27 @@ void AnimatedWaterSurface::Update(const std::shared_ptr<WaveBase>& wave, double
     // -------------------------------------------------------------------------
     std::vector<float> grid_eta(orig_verts.size(), 0.0f);
 
+    // For irregular waves, use limited frequency components for visualization.
+    // 500 components provides excellent visual quality while being faster than
+    // using all components (which can be 1000+). Physics always uses all components.
+    constexpr int kVizFrequencyLimit = 500;
+    
+    // Try to cast to IrregularWaves for the optimized visualization path.
+    IrregularWaves* irreg_wave = nullptr;
+    if (wave && wave->GetWaveMode() == WaveMode::irregular) {
+        irreg_wave = dynamic_cast<IrregularWaves*>(wave.get());
+    }
+
     for (size_t i = 0; i < orig_verts.size(); ++i) {
         const auto& grid_pos = orig_verts[i];
+        const Eigen::Vector3d pos(grid_pos.x(), grid_pos.y(), 0.0);
 
-        // Incident wave elevation from the wave model (regular or irregular).
+        // Incident wave elevation from the wave model.
+        // For irregular waves, use the visualization-optimized method with limited components.
         double eta_incident = 0.0;
-        if (wave) {
-            Eigen::Vector3d pos(grid_pos.x(), grid_pos.y(), 0.0);
+        if (irreg_wave) {
+            eta_incident = irreg_wave->GetElevationForVisualization(pos, t, kVizFrequencyLimit);
+        } else if (wave) {
             eta_incident = wave->GetElevation(pos, t);
         }
 

src/hydro/waves/irregular_wave.cpp

@@ -207,19 +207,18 @@ double IrregularWaves::GetElevation(const Eigen::Vector3d& position, double time
     //     x     = position along wave propagation direction [m]
     //     t     = time [s]
     //
-    // Performance: amplitudes_ and angular_freqs_ are pre-computed to avoid
-    // repeated sqrt() and multiplication calls (critical for visualization).
+    // Performance: Uses Eigen's vectorized operations (SIMD) for the summation.
+    // The cos() is applied element-wise and the dot product sums the result.
     // -------------------------------------------------------------------------
     const double x = position.x();
-    double eta = 0.0;
 
-    const Eigen::Index num_components = amplitudes_.size();
-    for (Eigen::Index i = 0; i < num_components; ++i) {
-        const double phase = wavenumbers_[i] * x - angular_freqs_[i] * time + wave_phases_[i];
-        eta += amplitudes_[i] * std::cos(phase);
-    }
+    // Vectorized phase computation: phase_i = k_i*x - ω_i*t + φ_i
+    const Eigen::ArrayXd phases = wavenumbers_.array() * x 
+                                - angular_freqs_.array() * time 
+                                + wave_phases_.array();
 
-    return eta;
+    // Vectorized elevation: η = Σ A_i * cos(phase_i)
+    return (amplitudes_.array() * phases.cos()).sum();
 }
 
 Eigen::Vector2d IrregularWaves::GetElevationGradientXY(const Eigen::Vector3d& position, double time) const {
@@ -240,19 +239,46 @@ Eigen::Vector2d IrregularWaves::GetElevationGradientXY(const Eigen::Vector3d& po
     //
     // Since waves propagate only in the +X direction, ∂η/∂y = 0.
     // The surface normal can be computed as: n = normalize(-∂η/∂x, -∂η/∂y, 1)
+    //
+    // Performance: Uses Eigen's vectorized operations (SIMD) for the summation.
     // -------------------------------------------------------------------------
     const double x = position.x();
-    double deta_dx = 0.0;
 
-    const Eigen::Index num_components = amplitudes_.size();
-    for (Eigen::Index i = 0; i < num_components; ++i) {
-        const double phase = wavenumbers_[i] * x - angular_freqs_[i] * time + wave_phases_[i];
-        deta_dx -= amplitudes_[i] * wavenumbers_[i] * std::sin(phase);
-    }
+    // Vectorized phase computation: phase_i = k_i*x - ω_i*t + φ_i
+    const Eigen::ArrayXd phases = wavenumbers_.array() * x 
+                                - angular_freqs_.array() * time 
+                                + wave_phases_.array();
+
+    // Vectorized gradient: ∂η/∂x = -Σ A_i * k_i * sin(phase_i)
+    const double deta_dx = -(amplitudes_.array() * wavenumbers_.array() * phases.sin()).sum();
 
     return Eigen::Vector2d(deta_dx, 0.0);
 }
 
+double IrregularWaves::GetElevationForVisualization(const Eigen::Vector3d& position, 
+                                                     double time, 
+                                                     int max_components) const {
+    // If no pre-computed amplitudes or max_components covers all, use full calculation.
+    const Eigen::Index num_total = amplitudes_.size();
+    if (num_total == 0) {
+        return GetEtaIrregular(position, time, spectrum_frequencies_, spectral_densities_,
+                               spectral_widths_, wave_phases_, wavenumbers_);
+    }
+    
+    // Determine how many components to use.
+    const Eigen::Index n = (max_components <= 0 || max_components >= num_total) 
+                         ? num_total 
+                         : static_cast<Eigen::Index>(max_components);
+    
+    // Use Eigen head() to get first n elements - still vectorized (SIMD).
+    const double x = position.x();
+    const Eigen::ArrayXd phases = wavenumbers_.head(n).array() * x 
+                                - angular_freqs_.head(n).array() * time 
+                                + wave_phases_.head(n).array();
+
+    return (amplitudes_.head(n).array() * phases.cos()).sum();
+}
+
 Eigen::VectorXd IrregularWaves::GetForceAtTime(double t) {
     unsigned int total_dofs = params_.num_bodies_ * 6;
     Eigen::VectorXd f(total_dofs);