33 *
44 * @brief implementation file for Wavebase and classes inheriting from WaveBase.
55 *********************************************************************/
6+ #include < hydroc/helper.h>
67#include < hydroc/wave_types.h>
78#include < unsupported/Eigen/Splines>
8- #include < hydroc/helper.h>
99
1010Eigen::VectorXd NoWave::GetForceAtTime (double t) {
1111 unsigned int dof = num_bodies_ * 6 ;
@@ -89,15 +89,15 @@ double RegularWave::GetExcitationPhaseInterp(int b, int i, int j, double freq_in
8989// ///////////////////////////////////////////////////////////////////////////////////////////////////////
9090
9191// // Irregular wave class definitions:
92- // IrregularWave::IrregularWave() {
92+ // IrregularWave::IrregularWave() {
9393// num_bodies = 1;
9494// }
9595//
96- // IrregularWave::IrregularWave(unsigned int num_b) {
96+ // IrregularWave::IrregularWave(unsigned int num_b) {
9797// num_bodies = num_b;
9898// }
9999//
100- // void IrregularWave::Initialize() {
100+ // void IrregularWave::Initialize() {
101101// std::vector<Eigen::MatrixXd> ex_irf_old(num_bodies);
102102// std::vector<Eigen::VectorXd> ex_irf_time_old(num_bodies);
103103// for (int b = 0; b < num_bodies; b++) {
@@ -120,7 +120,7 @@ double RegularWave::GetExcitationPhaseInterp(int b, int i, int j, double freq_in
120120// CreateFreeSurfaceElevation(); // eta initialized in here, and output to eta.txt
121121// }
122122//
123- // Eigen::MatrixXd IrregularWave::ResampleVals(const Eigen::VectorXd& t_old,
123+ // Eigen::MatrixXd IrregularWave::ResampleVals(const Eigen::VectorXd& t_old,
124124// Eigen::MatrixXd& vals_old,
125125// const Eigen::VectorXd& t_new) {
126126// assert(vals_old.rows() == 6);
@@ -140,7 +140,7 @@ double RegularWave::GetExcitationPhaseInterp(int b, int i, int j, double freq_in
140140// return vals_new;
141141// }
142142//
143- // Eigen::VectorXd IrregularWave::ResampleTime(const Eigen::VectorXd& t_old, const double dt_new) {
143+ // Eigen::VectorXd IrregularWave::ResampleTime(const Eigen::VectorXd& t_old, const double dt_new) {
144144// double dt_old = t_old[1] - t_old[0];
145145// int size_new = static_cast<int>(ceil(t_old.size() * dt_old / dt_new));
146146// double t_initial = t_old[0];
@@ -161,13 +161,13 @@ double RegularWave::GetExcitationPhaseInterp(int b, int i, int j, double freq_in
161161// return t_new;
162162// }
163163//
164- // void IrregularWave::AddH5Data(std::vector<HydroData::IrregularWaveInfo>& irreg_h5_data,
164+ // void IrregularWave::AddH5Data(std::vector<HydroData::IrregularWaveInfo>& irreg_h5_data,
165165// HydroData::SimulationParameters& sim_data) {
166166// wave_info = irreg_h5_data;
167167// sim_data = sim_data;
168168// }
169169//
170- // Eigen::VectorXd IrregularWave::GetForceAtTime(double t) {
170+ // Eigen::VectorXd IrregularWave::GetForceAtTime(double t) {
171171// unsigned int total_dofs = num_bodies * 6;
172172// Eigen::VectorXd f(total_dofs);
173173// // initialize the force here:
@@ -199,7 +199,7 @@ double RegularWave::GetExcitationPhaseInterp(int b, int i, int j, double freq_in
199199// // return wave_info[b].excitation_irf_matrix;
200200// //}
201201//
202- // double IrregularWave::ExcitationConvolution(int body, int dof, double time) {
202+ // double IrregularWave::ExcitationConvolution(int body, int dof, double time) {
203203// double f_ex = 0.0;
204204// double width = ex_irf_time_resampled[body][1] - ex_irf_time_resampled[body][0];
205205// //std::cout << "width=" << width << std::endl;
@@ -218,7 +218,7 @@ double RegularWave::GetExcitationPhaseInterp(int b, int i, int j, double freq_in
218218// return f_ex;
219219// }
220220//
221- // Eigen::VectorXd IrregularWave::SetSpectrumFrequencies(double start, double end, int num_points) {
221+ // Eigen::VectorXd IrregularWave::SetSpectrumFrequencies(double start, double end, int num_points) {
222222// Eigen::VectorXd result(num_points);
223223// double step = (end - start) / (num_points - 1);
224224//
@@ -231,7 +231,7 @@ double RegularWave::GetExcitationPhaseInterp(int b, int i, int j, double freq_in
231231// return result;
232232// }
233233//
234- // void IrregularWave::CreateSpectrum() {
234+ // void IrregularWave::CreateSpectrum() {
235235// // Define the frequency vector
236236// spectrum_frequencies = Eigen::VectorXd::LinSpaced(1000, 0.001, 1.0);
237237//
@@ -255,7 +255,7 @@ double RegularWave::GetExcitationPhaseInterp(int b, int i, int j, double freq_in
255255// }
256256// }
257257//
258- // void IrregularWave::CreateFreeSurfaceElevation() {
258+ // void IrregularWave::CreateFreeSurfaceElevation() {
259259// // Create a time index vector
260260// // UpdateNumTimesteps();
261261// int num_timesteps = static_cast<int>(simulation_duration / simulation_dt) + 1;
@@ -296,7 +296,6 @@ double RegularWave::GetExcitationPhaseInterp(int b, int i, int j, double freq_in
296296// // WriteFreeSurfaceMeshObj(free_surface_3d_pts, free_surface_triangles, "fse_mesh.obj");
297297// }
298298
299-
300299std::vector<double > ComputeWaveNumbers (const std::vector<double >& omegas,
301300 double water_depth,
302301 double g = 9.81 ,
@@ -330,10 +329,10 @@ std::vector<double> ComputeWaveNumbers(const std::vector<double>& omegas,
330329}
331330
332331std::vector<double > FreeSurfaceElevation (const Eigen::VectorXd& freqs_hz,
333- const Eigen::VectorXd& spectral_densities,
334- const Eigen::VectorXd& time_index,
335- double water_depth,
336- int seed) {
332+ const Eigen::VectorXd& spectral_densities,
333+ const Eigen::VectorXd& time_index,
334+ double water_depth,
335+ int seed) {
337336 double delta_f = freqs_hz (Eigen::last) / freqs_hz.size ();
338337 std::vector<double > omegas (freqs_hz.size ());
339338
@@ -377,7 +376,7 @@ std::vector<double> FreeSurfaceElevation(const Eigen::VectorXd& freqs_hz,
377376 return eta;
378377}
379378
380- // void IrregularWave::SetUpWaveMesh(std::string filename) {
379+ // void IrregularWave::SetUpWaveMesh(std::string filename) {
381380// mesh_file_name = filename;
382381// int num_timesteps = static_cast<int>(simulation_duration / simulation_dt) + 1;
383382// Eigen::VectorXd time_index = Eigen::VectorXd::LinSpaced(num_timesteps, 0, simulation_duration);
@@ -452,16 +451,15 @@ void WriteFreeSurfaceMeshObj(const std::vector<std::array<double, 3>>& points,
452451 out.close ();
453452}
454453
455- // std::string IrregularWave::GetMeshFile() {
454+ // std::string IrregularWave::GetMeshFile() {
456455// return mesh_file_name;
457456// }
458457//
459- // Eigen::Vector3<double> IrregularWave::GetWaveMeshVelocity() {
458+ // Eigen::Vector3<double> IrregularWave::GetWaveMeshVelocity() {
460459// return Eigen::Vector3d(1.0, 0, 0);
461460// }
462461
463-
464- // IrregularWaves::IrregularWaves(const std::string& eta_file_path, double wave_height, double wave_period)
462+ // IrregularWaves::IrregularWaves(const std::string& eta_file_path, double wave_height, double wave_period)
465463// : eta_file_path(eta_file_path), wave_height(wave_height), wave_period(wave_period) {
466464// if (!eta_file_path.empty()) {
467465// ReadEtaFromFile();
@@ -501,16 +499,16 @@ IrregularWaves::IrregularWaves(const IrregularWaveParams& params)
501499 std::vector<Eigen::MatrixXd> ex_irf_old (num_bodies_);
502500 std::vector<Eigen::VectorXd> ex_irf_time_old (num_bodies_);
503501
504- // for (unsigned int b = 0; b < num_bodies; b++) {
502+ // for (unsigned int b = 0; b < num_bodies; b++) {
505503 // std::cout << "get excitation IRFs..." << std::endl;
506504 // ex_irf_old[b] = GetExcitationIRF(b);
507505 // std::cout << "get excitation IRF time..." << std::endl;
508506 // ex_irf_time_old[b] = wave_info[b].excitation_irf_time;
509507 // }
510508
511509 // // Resample excitation IRF time series
512- // std::cout << "resample IRFs..." << std::endl;
513- // for (unsigned int b = 0; b < num_bodies; b++) {
510+ // std::cout << "resample IRFs..." << std::endl;
511+ // for (unsigned int b = 0; b < num_bodies; b++) {
514512 // ex_irf_time_resampled[b] = ResampleTime(ex_irf_time_old[b], simulation_dt);
515513 // ex_irf_resampled[b] = ResampleVals(ex_irf_time_old[b], ex_irf_old[b], ex_irf_time_resampled[b]);
516514 // }
@@ -582,10 +580,10 @@ Eigen::MatrixXd IrregularWaves::GetExcitationIRF(int b) const {
582580void IrregularWaves::AddH5Data (std::vector<HydroData::IrregularWaveInfo>& irreg_h5_data,
583581 HydroData::SimulationParameters& sim_data) {
584582 wave_info_ = irreg_h5_data;
585- sim_data = sim_data; // Note: This line seems to be redundant as mentioned before.
583+ sim_data = sim_data; // Note: This line seems to be redundant as mentioned before.
586584
587585 // // Add this loop to print the values of wave_info[b].excitation_irf_matrix for each body
588- // for (unsigned int b = 0; b < num_bodies; ++b) {
586+ // for (unsigned int b = 0; b < num_bodies; ++b) {
589587 // std::cout << "wave_info[" << b << "].excitation_irf_matrix:\n";
590588 // std::cout << wave_info[b].excitation_irf_matrix << std::endl;
591589 // }
@@ -614,13 +612,13 @@ Eigen::VectorXd IrregularWaves::GetForceAtTime(double t) {
614612 f[b_offset + dof] = f_dof;
615613 }
616614 }
617-
615+
618616 return f;
619617}
620618
621619Eigen::MatrixXd IrregularWaves::ResampleVals (const Eigen::VectorXd& t_old,
622- Eigen::MatrixXd& vals_old,
623- const Eigen::VectorXd& t_new) {
620+ Eigen::MatrixXd& vals_old,
621+ const Eigen::VectorXd& t_new) {
624622 assert (vals_old.rows () == 6 );
625623
626624 Eigen::MatrixXd vals_new (6 , t_new.size ());
@@ -635,11 +633,10 @@ Eigen::MatrixXd IrregularWaves::ResampleVals(const Eigen::VectorXd& t_old,
635633 vals_new.col (i) = spline (t_new_scaled[i]);
636634 }
637635
638- // std::cout << "Old time vector:\n" << t_old << std::endl;
639- // std::cout << "Old values:\n" << vals_old << std::endl;
640- // std::cout << "New time vector:\n" << t_new << std::endl;
641- // std::cout << "New values:\n" << vals_new << std::endl; // assuming new_vals is the output of the function
642-
636+ // std::cout << "Old time vector:\n" << t_old << std::endl;
637+ // std::cout << "Old values:\n" << vals_old << std::endl;
638+ // std::cout << "New time vector:\n" << t_new << std::endl;
639+ // std::cout << "New values:\n" << vals_new << std::endl; // assuming new_vals is the output of the function
643640
644641 return vals_new;
645642}
@@ -677,7 +674,8 @@ void IrregularWaves::CreateSpectrum() {
677674 spectrum_frequencies_ = Eigen::VectorXd::LinSpaced (1000 , 0.001 , 1.0 );
678675
679676 // Calculate the Pierson-Moskowitz Spectrum
680- spectral_densities_ = JONSWAPSpectrumHz (spectrum_frequencies_, wave_height_, wave_period_, peak_enhancement_factor_);
677+ spectral_densities_ =
678+ JONSWAPSpectrumHz (spectrum_frequencies_, wave_height_, wave_period_, peak_enhancement_factor_);
681679
682680 // Open a file stream for writing
683681 std::ofstream outputFile (" spectral_densities.txt" );
@@ -777,18 +775,18 @@ double IrregularWaves::ExcitationConvolution(int body, int dof, double time) {
777775 double width = ex_irf_time_resampled_[body][1 ] - ex_irf_time_resampled_[body][0 ];
778776 // std::cout << "width=" << width << std::endl;
779777 for (size_t j = 0 ; j < ex_irf_time_resampled_[0 ].size (); ++j) {
780- double tau = ex_irf_time_resampled_[0 ][j];
781- double t_tau = time - tau;
782-
778+ double tau = ex_irf_time_resampled_[0 ][j];
779+ double t_tau = time - tau;
780+
783781 double ex_irf_val = ex_irf_resampled_[body](dof, j);
784782
785783 if (0.0 < t_tau && t_tau < free_surface_elevation_.size () * width) {
786784 size_t eta_index = static_cast <size_t >(t_tau / width);
787785 double eta_val = free_surface_elevation_[eta_index - 1 ];
788786
789- // std::cout << "ex_irf_val = " << ex_irf_val << std::endl;
790- // std::cout << "eta_val = " << eta_val << std::endl;
791- // std::cout << "width = " << width << std::endl;
787+ // std::cout << "ex_irf_val = " << ex_irf_val << std::endl;
788+ // std::cout << "eta_val = " << eta_val << std::endl;
789+ // std::cout << "width = " << width << std::endl;
792790
793791 f_ex += ex_irf_val * eta_val * width; // eta is wave elevation
794792 }
@@ -797,11 +795,12 @@ double IrregularWaves::ExcitationConvolution(int body, int dof, double time) {
797795 return f_ex;
798796}
799797
800- void IrregularWaves::SetUpWaveMesh (std::string filename) {
801- mesh_file_name_ = filename;
798+ void IrregularWaves::SetUpWaveMesh (std::string filename) {
799+ mesh_file_name_ = filename;
802800 int num_timesteps = static_cast <int >(simulation_duration_ / simulation_dt_) + 1 ;
803801 Eigen::VectorXd time_index = Eigen::VectorXd::LinSpaced (num_timesteps, 0 , simulation_duration_);
804- std::vector<std::array<double , 3 >> free_surface_3d_pts = CreateFreeSurface3DPts (free_surface_elevation_, time_index);
802+ std::vector<std::array<double , 3 >> free_surface_3d_pts =
803+ CreateFreeSurface3DPts (free_surface_elevation_, time_index);
805804 std::vector<std::array<size_t , 3 >> free_surface_triangles = CreateFreeSurfaceTriangles (time_index.size ());
806805
807806 WriteFreeSurfaceMeshObj (free_surface_3d_pts, free_surface_triangles, mesh_file_name_);
0 commit comments