@@ -86,216 +86,6 @@ double RegularWave::GetExcitationPhaseInterp(int b, int i, int j, double freq_in
8686 return excitationPhase;
8787}
8888
89- // ///////////////////////////////////////////////////////////////////////////////////////////////////////
90-
91- // // Irregular wave class definitions:
92- // IrregularWave::IrregularWave() {
93- // num_bodies = 1;
94- // }
95- //
96- // IrregularWave::IrregularWave(unsigned int num_b) {
97- // num_bodies = num_b;
98- // }
99- //
100- // void IrregularWave::Initialize() {
101- // std::vector<Eigen::MatrixXd> ex_irf_old(num_bodies);
102- // std::vector<Eigen::VectorXd> ex_irf_time_old(num_bodies);
103- // for (int b = 0; b < num_bodies; b++) {
104- // ex_irf_old[b] = GetExcitationIRF(b);
105- // ex_irf_time_old[b] = wave_info[b].excitation_irf_time;
106- // }
107- //
108- // // resample excitation IRF time series
109- // // h5 file irf has different timestep, want to resample with interpolation (cubic spline?) once at start no
110- // // interpolation in convolution integral part
111- // // different one for each body it's 6x1x1000 so maybe switch to 2d reading
112- // ex_irf_resampled.resize(num_bodies);
113- // ex_irf_time_resampled.resize(num_bodies);
114- // for (int b = 0; b < num_bodies; b++) { // function call (time_in, vals_in, &time_out, &vals_out)
115- // ex_irf_time_resampled[b] = ResampleTime(ex_irf_time_old[b], simulation_dt);
116- // ex_irf_resampled[b] = ResampleVals(ex_irf_time_old[b], ex_irf_old[b], ex_irf_time_resampled[b]);
117- // }
118- //
119- // CreateSpectrum(); // output in spectral_densities.txt
120- // CreateFreeSurfaceElevation(); // eta initialized in here, and output to eta.txt
121- // }
122- //
123- // Eigen::MatrixXd IrregularWave::ResampleVals(const Eigen::VectorXd& t_old,
124- // Eigen::MatrixXd& vals_old,
125- // const Eigen::VectorXd& t_new) {
126- // assert(vals_old.rows() == 6);
127- //
128- // Eigen::MatrixXd vals_new(6, t_new.size());
129- // // we need to scale t to be [0,1] for spline use
130- // // TODO: change this to accomodate variable dt instead of const dt
131- // Eigen::VectorXd t_old_scaled = Eigen::VectorXd::LinSpaced(t_old.size(), 0, 1);
132- // Eigen::VectorXd t_new_scaled = Eigen::VectorXd::LinSpaced(t_new.size(), 0, 1);
133- //
134- // Eigen::Spline<double, 6> spline =
135- // Eigen::SplineFitting<Eigen::Spline<double, 6>>::Interpolate(vals_old, 3, t_old_scaled);
136- // for (int i = 0; i < t_new.rows(); i++) {
137- // vals_new.col(i) = spline(t_new_scaled[i]);
138- // }
139- //
140- // return vals_new;
141- // }
142- //
143- // Eigen::VectorXd IrregularWave::ResampleTime(const Eigen::VectorXd& t_old, const double dt_new) {
144- // double dt_old = t_old[1] - t_old[0];
145- // int size_new = static_cast<int>(ceil(t_old.size() * dt_old / dt_new));
146- // double t_initial = t_old[0];
147- // double t_final = t_old[t_old.size() - 1];
148- // // t_0 and t_final should be the same in old/new versions, use new time step, which means a different number of
149- // // timesteps
150- // Eigen::VectorXd t_new = Eigen::VectorXd::LinSpaced(size_new, t_initial, t_final);
151- //
152- // // print for testing
153- // int N_old = t_old.size() - 1;
154- // int N_new = t_new.size() - 1;
155- // //std::cout << "T_old = {t_i | i = 0 .. " << N_old << ", t_0 = " << t_old[0] << ", t_" << N_old << " = "
156- // // << t_old[N_old] << "}\n"
157- // // << "dt_old = " << dt_old << "\nT_new = {t_i | i = 0 .. " << N_new << ", t_0 = " << t_new[0] << ", t_"
158- // // << N_new << " = " << t_new[N_new] << "}\n"
159- // // << "dt_new = " << dt_new << std::endl;
160- //
161- // return t_new;
162- // }
163- //
164- // void IrregularWave::AddH5Data(std::vector<HydroData::IrregularWaveInfo>& irreg_h5_data,
165- // HydroData::SimulationParameters& sim_data) {
166- // wave_info = irreg_h5_data;
167- // sim_data = sim_data;
168- // }
169- //
170- // Eigen::VectorXd IrregularWave::GetForceAtTime(double t) {
171- // unsigned int total_dofs = num_bodies * 6;
172- // Eigen::VectorXd f(total_dofs);
173- // // initialize the force here:
174- // // for now set f to all zeros
175- // for (int i = 0; i < total_dofs; i++) {
176- // f[i] = 0.0;
177- // }
178- // // see ComputeForceExcitation and convolution functions
179- //
180- // // force_excitation.resize(total_dofs, 0.0);
181- //
182- // for (int body = 0; body < num_bodies; body++) {
183- // // Loop through the DOFs
184- // for (int dof = 0; dof < 6; ++dof) {
185- // // Compute the convolution for the current DOF
186- // double f_dof = ExcitationConvolution(body, dof, t);
187- // unsigned int b_offset = body * 6;
188- // f[b_offset + dof] = f_dof;
189- // }
190- // }
191- // return f;
192- // }
193- //
194- // ///*******************************************************************************
195- // // * IrregularWave::GetExcitationIRF()
196- // // * returns the std::vector of excitation_irf_matrix from h5 file
197- // // *******************************************************************************/
198- // //Eigen::MatrixXd IrregularWave::GetExcitationIRF(int b) const {
199- // // return wave_info[b].excitation_irf_matrix;
200- // //}
201- //
202- // double IrregularWave::ExcitationConvolution(int body, int dof, double time) {
203- // double f_ex = 0.0;
204- // double width = ex_irf_time_resampled[body][1] - ex_irf_time_resampled[body][0];
205- // //std::cout << "width=" << width << std::endl;
206- //
207- // for (size_t j = 0; j < ex_irf_time_resampled[0].size(); ++j) {
208- // double tau = ex_irf_time_resampled[0][j];
209- // double t_tau = time - tau;
210- // double ex_irf_val = ex_irf_resampled[body](dof, j);
211- // if (0.0 < t_tau && t_tau < eta.size() * width) {
212- // size_t eta_index = static_cast<size_t>(t_tau / width);
213- // double eta_val = eta[eta_index - 1];
214- // f_ex += ex_irf_val * eta_val * width; // eta is wave elevation
215- // }
216- // }
217- //
218- // return f_ex;
219- // }
220- //
221- // Eigen::VectorXd IrregularWave::SetSpectrumFrequencies(double start, double end, int num_points) {
222- // Eigen::VectorXd result(num_points);
223- // double step = (end - start) / (num_points - 1);
224- //
225- // for (int i = 0; i < num_points; ++i) {
226- // result[i] = start + i * step;
227- // }
228- //
229- // spectrum_frequencies = result;
230- //
231- // return result;
232- // }
233- //
234- // void IrregularWave::CreateSpectrum() {
235- // // Define the frequency vector
236- // spectrum_frequencies = Eigen::VectorXd::LinSpaced(1000, 0.001, 1.0);
237- //
238- // // Calculate the Pierson-Moskowitz Spectrum
239- // spectral_densities = PiersonMoskowitzSpectrumHz(spectrum_frequencies, wave_height, wave_period);
240- //
241- // // Open a file stream for writing
242- // std::ofstream outputFile("spectral_densities.txt");
243- //
244- // // Check if the file stream is open
245- // if (outputFile.is_open()) {
246- // // Write the spectral densities and their corresponding frequencies to the file
247- // for (size_t i = 0; i < spectral_densities.size(); ++i) {
248- // outputFile << spectrum_frequencies[i] << " : " << spectral_densities[i] << std::endl;
249- // }
250- //
251- // // Close the file stream
252- // outputFile.close();
253- // } else {
254- // std::cerr << "Unable to open file for writing." << std::endl;
255- // }
256- // }
257- //
258- // void IrregularWave::CreateFreeSurfaceElevation() {
259- // // Create a time index vector
260- // // UpdateNumTimesteps();
261- // int num_timesteps = static_cast<int>(simulation_duration / simulation_dt) + 1;
262- //
263- // Eigen::VectorXd time_index = Eigen::VectorXd::LinSpaced(num_timesteps, 0, simulation_duration);
264- //
265- // // Calculate the free surface elevation
266- // eta = FreeSurfaceElevation(spectrum_frequencies, spectral_densities, time_index, sim_data.water_depth);
267- //
268- // // Apply ramp if ramp_duration is greater than 0
269- // if (ramp_duration > 0.0) {
270- // // UpdateRampTimesteps();
271- // int ramp_timesteps = static_cast<int>(ramp_duration / simulation_dt) + 1;
272- // Eigen::VectorXd ramp = Eigen::VectorXd::LinSpaced(ramp_timesteps, 0.0, 1.0);
273- //
274- // for (size_t i = 0; i < ramp.size(); ++i) {
275- // eta[i] *= ramp[i];
276- // }
277- // }
278- //
279- // // Open a file stream for writing
280- // std::ofstream eta_output("eta.txt");
281- // // Check if the file stream is open
282- // if (eta_output.is_open()) {
283- // // Write the spectral densities and their corresponding frequencies to the file
284- // for (size_t i = 0; i < eta.size(); ++i) {
285- // eta_output << time_index[i] << " : " << eta[i] << std::endl;
286- // }
287- // // Close the file stream
288- // eta_output.close();
289- // } else {
290- // std::cerr << "Unable to open file for writing." << std::endl;
291- // }
292- //
293- // // std::vector<std::array<double, 3>> free_surface_3d_pts = CreateFreeSurface3DPts(eta, time_index);
294- // // std::vector<std::array<size_t, 3>> free_surface_triangles = CreateFreeSurfaceTriangles(time_index.size());
295- //
296- // // WriteFreeSurfaceMeshObj(free_surface_3d_pts, free_surface_triangles, "fse_mesh.obj");
297- // }
298-
29989std::vector<double > ComputeWaveNumbers (const std::vector<double >& omegas,
30090 double water_depth,
30191 double g = 9.81 ,
@@ -376,16 +166,6 @@ std::vector<double> FreeSurfaceElevation(const Eigen::VectorXd& freqs_hz,
376166 return eta;
377167}
378168
379- // void IrregularWave::SetUpWaveMesh(std::string filename) {
380- // mesh_file_name = filename;
381- // int num_timesteps = static_cast<int>(simulation_duration / simulation_dt) + 1;
382- // Eigen::VectorXd time_index = Eigen::VectorXd::LinSpaced(num_timesteps, 0, simulation_duration);
383- // std::vector<std::array<double, 3>> free_surface_3d_pts = CreateFreeSurface3DPts(eta, time_index);
384- // std::vector<std::array<size_t, 3>> free_surface_triangles = CreateFreeSurfaceTriangles(time_index.size());
385- //
386- // WriteFreeSurfaceMeshObj(free_surface_3d_pts, free_surface_triangles, mesh_file_name);
387- // }
388-
389169std::vector<std::array<double , 3 >> CreateFreeSurface3DPts (const std::vector<double >& eta,
390170 const Eigen::VectorXd& t_vec) {
391171 std::vector<std::array<double , 3 >> surface (t_vec.size () * 2 );
@@ -451,26 +231,6 @@ void WriteFreeSurfaceMeshObj(const std::vector<std::array<double, 3>>& points,
451231 out.close ();
452232}
453233
454- // std::string IrregularWave::GetMeshFile() {
455- // return mesh_file_name;
456- // }
457- //
458- // Eigen::Vector3<double> IrregularWave::GetWaveMeshVelocity() {
459- // return Eigen::Vector3d(1.0, 0, 0);
460- // }
461-
462- // IrregularWaves::IrregularWaves(const std::string& eta_file_path, double wave_height, double wave_period)
463- // : eta_file_path(eta_file_path), wave_height(wave_height), wave_period(wave_period) {
464- // if (!eta_file_path.empty()) {
465- // ReadEtaFromFile();
466- // spectrumCreated = false;
467- // } else {
468- // CreateSpectrum();
469- // CreateFreeSurfaceElevation();
470- // spectrumCreated = true;
471- // }
472- // }
473-
474234IrregularWaves::IrregularWaves (const IrregularWaveParams& params)
475235 : num_bodies_(params.num_bodies_),
476236 eta_file_path_(params.eta_file_path_),
@@ -500,20 +260,6 @@ IrregularWaves::IrregularWaves(const IrregularWaveParams& params)
500260 std::vector<Eigen::MatrixXd> ex_irf_old (num_bodies_);
501261 std::vector<Eigen::VectorXd> ex_irf_time_old (num_bodies_);
502262
503- // for (unsigned int b = 0; b < num_bodies; b++) {
504- // std::cout << "get excitation IRFs..." << std::endl;
505- // ex_irf_old[b] = GetExcitationIRF(b);
506- // std::cout << "get excitation IRF time..." << std::endl;
507- // ex_irf_time_old[b] = wave_info[b].excitation_irf_time;
508- // }
509-
510- // // Resample excitation IRF time series
511- // std::cout << "resample IRFs..." << std::endl;
512- // for (unsigned int b = 0; b < num_bodies; b++) {
513- // ex_irf_time_resampled[b] = ResampleTime(ex_irf_time_old[b], simulation_dt);
514- // ex_irf_resampled[b] = ResampleVals(ex_irf_time_old[b], ex_irf_old[b], ex_irf_time_resampled[b]);
515- // }
516-
517263 std::cout << " IrregularWaves instantiated..." << std::endl;
518264}
519265
@@ -570,39 +316,23 @@ void IrregularWaves::ReadEtaFromFile() {
570316 }
571317}
572318
573- /* ******************************************************************************
574- * IrregularWave::GetExcitationIRF()
575- * returns the std::vector of excitation_irf_matrix from h5 file
576- *******************************************************************************/
577319Eigen::MatrixXd IrregularWaves::GetExcitationIRF (int b) const {
578320 return wave_info_[b].excitation_irf_matrix ;
579321}
580322
581323void IrregularWaves::AddH5Data (std::vector<HydroData::IrregularWaveInfo>& irreg_h5_data,
582324 HydroData::SimulationParameters& sim_data) {
583325 wave_info_ = irreg_h5_data;
584- sim_data = sim_data; // Note: This line seems to be redundant as mentioned before.
585-
586- // // Add this loop to print the values of wave_info[b].excitation_irf_matrix for each body
587- // for (unsigned int b = 0; b < num_bodies; ++b) {
588- // std::cout << "wave_info[" << b << "].excitation_irf_matrix:\n";
589- // std::cout << wave_info[b].excitation_irf_matrix << std::endl;
590- // }
591326
592327 InitializeIRFVectors ();
593328}
594329
595330Eigen::VectorXd IrregularWaves::GetForceAtTime (double t) {
596331 unsigned int total_dofs = num_bodies_ * 6 ;
597332 Eigen::VectorXd f (total_dofs);
598- // initialize the force here:
599- // for now set f to all zeros
600333 for (int i = 0 ; i < total_dofs; i++) {
601334 f[i] = 0.0 ;
602335 }
603- // see ComputeForceExcitation and convolution functions
604-
605- // force_excitation.resize(total_dofs, 0.0);
606336
607337 for (int body = 0 ; body < num_bodies_; body++) {
608338 // Loop through the DOFs
@@ -619,33 +349,36 @@ Eigen::VectorXd IrregularWaves::GetForceAtTime(double t) {
619349
620350void IrregularWaves::ResampleIRF (double dt) {
621351 for (unsigned int b = 0 ; b < num_bodies_; b++) {
622- // 0) access arrays to modify at body index
623352 auto & time_array = ex_irf_time_sampled_[b];
624353 auto & width_array = ex_irf_width_sampled_[b];
625354 auto & val_array = ex_irf_sampled_[b];
626- // copy time array
355+
356+ // Copy time array
627357 auto time_array_old = time_array;
628358
629- // 1) resample time
359+ // 1) Resample time
630360 auto t0 = time_array_old[0 ];
631361 auto t1 = time_array_old[time_array_old.size () - 1 ];
632362 time_array = Eigen::VectorXd::LinSpaced (static_cast <int >(ceil ((t1 - t0) / dt)), t0, t1);
633363
634- // 2) resample width
364+ // 2) Resample width
635365 CalculateWidthIRF ();
636366
637- // 3) resample values
367+ // 3) Resample values
638368 assert (val_array.rows () == 6 );
639369 Eigen::MatrixXd vals_new (6 , time_array.size ());
640- // we need to scale t to be [0,1] for spline use
370+
371+ // We need to scale t to be [0,1] for spline use
641372 Eigen::VectorXd t_old_scaled = Eigen::VectorXd::LinSpaced (time_array_old.size (), 0 , 1 );
642373 Eigen::VectorXd t_new_scaled = Eigen::VectorXd::LinSpaced (time_array.size (), 0 , 1 );
643- // use spline to get new values
374+
375+ // Use spline to get new values
644376 Eigen::Spline<double , 6 > spline =
645377 Eigen::SplineFitting<Eigen::Spline<double , 6 >>::Interpolate (val_array, 3 , t_old_scaled);
646378 for (int i = 0 ; i < time_array.rows (); i++) {
647379 vals_new.col (i) = spline (t_new_scaled[i]);
648380 }
381+
649382 val_array = vals_new;
650383 }
651384}
@@ -655,7 +388,6 @@ void IrregularWaves::CalculateWidthIRF() {
655388 auto & time_array = ex_irf_time_sampled_[b];
656389 auto & width_array = ex_irf_width_sampled_[b];
657390
658- // 2) calculate width
659391 width_array.resize (time_array.size ());
660392 for (int ii = 0 ; ii < width_array.size (); ii++) {
661393 width_array[ii] = 0.0 ;
@@ -747,7 +479,8 @@ void IrregularWaves::CreateFreeSurfaceElevation() {
747479 int num_timesteps = static_cast <int >(simulation_duration_ / simulation_dt_) + 1 ;
748480
749481 auto time_array = Eigen::VectorXd::LinSpaced (num_timesteps, 0 , simulation_duration_);
750- // save time array as a std::vector
482+
483+ // Save time array as a std::vector
751484 free_surface_time_sampled_.resize (time_array.size ());
752485 Eigen::VectorXd::Map (&free_surface_time_sampled_[0 ], time_array.size ()) = time_array;
753486
@@ -768,6 +501,7 @@ void IrregularWaves::CreateFreeSurfaceElevation() {
768501
769502 // Open a file stream for writing
770503 std::ofstream eta_output (" eta.txt" );
504+
771505 // Check if the file stream is open
772506 if (eta_output.is_open ()) {
773507 // Write the spectral densities and their corresponding frequencies to the file
@@ -782,7 +516,6 @@ void IrregularWaves::CreateFreeSurfaceElevation() {
782516
783517 // std::vector<std::array<double, 3>> free_surface_3d_pts = CreateFreeSurface3DPts(eta, time_index);
784518 // std::vector<std::array<size_t, 3>> free_surface_triangles = CreateFreeSurfaceTriangles(time_index.size());
785-
786519 // WriteFreeSurfaceMeshObj(free_surface_3d_pts, free_surface_triangles, "fse_mesh.obj");
787520}
788521
0 commit comments