From 01279f6ccb975ee4477b646dfd74defcec48d447 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Wed, 24 Jun 2026 16:31:06 -0500 Subject: [PATCH 01/19] Move fcType to separate header and source files --- Code/Source/solver/CMakeLists.txt | 3 +- Code/Source/solver/ComMod.h | 37 +----------------------- Code/Source/solver/fcType.cpp | 4 +++ Code/Source/solver/fcType.h | 47 +++++++++++++++++++++++++++++++ 4 files changed, 54 insertions(+), 37 deletions(-) create mode 100644 Code/Source/solver/fcType.cpp create mode 100644 Code/Source/solver/fcType.h diff --git a/Code/Source/solver/CMakeLists.txt b/Code/Source/solver/CMakeLists.txt index c546c2822..c0e299a0d 100644 --- a/Code/Source/solver/CMakeLists.txt +++ b/Code/Source/solver/CMakeLists.txt @@ -156,7 +156,8 @@ set(CSRCS PetscLinearAlgebra.h PetscLinearAlgebra.cpp TrilinosLinearAlgebra.h TrilinosLinearAlgebra.cpp Tensor4.h Tensor4.cpp - Vector.h Vector.cpp + Vector.h Vector.cpp + fcType.cpp lapack_defs.h diff --git a/Code/Source/solver/ComMod.h b/Code/Source/solver/ComMod.h index 86db290f4..4ca4114ad 100644 --- a/Code/Source/solver/ComMod.h +++ b/Code/Source/solver/ComMod.h @@ -21,6 +21,7 @@ #include "CoupledBoundaryCondition.h" #include "Timer.h" #include "Vector.h" +#include "fcType.h" #include "DebugMsg.h" @@ -43,42 +44,6 @@ class LinearAlgebra; -/// @brief Fourier coefficients that are used to specify unsteady BCs -// -class fcType -{ - public: - - bool defined() { return n != 0; }; - - // If this is a ramp function - bool lrmp = false; - - // Number of Fourier coefficient - int n = 0; - - // No. of dimensions (scalar or vector) - int d = 0; - - // Initial value - Vector qi; - - // Time derivative of linear part - Vector qs; - - // Period - double T = 0.0; - - // Initial time - double ti = 0.0; - - // Imaginary part of coefficint - Array i; - - // Real part of coefficint - Array r; -}; - /// @brief Moving boundary data structure (used for general BC) // class MBType diff --git a/Code/Source/solver/fcType.cpp b/Code/Source/solver/fcType.cpp new file mode 100644 index 000000000..29e2b7959 --- /dev/null +++ b/Code/Source/solver/fcType.cpp @@ -0,0 +1,4 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the +// University of California, and others. SPDX-License-Identifier: BSD-3-Clause + +#include "fcType.h" \ No newline at end of file diff --git a/Code/Source/solver/fcType.h b/Code/Source/solver/fcType.h new file mode 100644 index 000000000..342103b7f --- /dev/null +++ b/Code/Source/solver/fcType.h @@ -0,0 +1,47 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the +// University of California, and others. SPDX-License-Identifier: BSD-3-Clause + +#ifndef FC_TYPE_H +#define FC_TYPE_H + +#include "Array.h" +#include "Vector.h" + +/// @brief Fourier coefficients that are used to specify time-dependent +/// functions. +/// +/// These are used for boundary conditions and for time-dependent prescribed +/// active stress. +class fcType { +public: + bool defined() { return n != 0; }; + + /// Toggle whether this is a ramp function or not. + bool lrmp = false; + + /// Number of Fourier coefficients. + int n = 0; + + /// Numbero of dimensions (scalar or vector). + int d = 0; + + /// Initial value. + Vector qi; + + /// Time derivative of linear part. + Vector qs; + + /// Period. + double T = 0.0; + + /// Initial time. + double ti = 0.0; + + /// Imaginary part of coefficients. + Array i; + + /// Real part of coefficients. + Array r; +}; + +#endif \ No newline at end of file From 0fb860b54ce978ff3b529d387c65d1357fc42673 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Fri, 26 Jun 2026 12:49:07 -0500 Subject: [PATCH 02/19] Move evaluation of inverse Fourier transform to fcType methods --- Code/Source/solver/bf.cpp | 3 +- Code/Source/solver/fcType.cpp | 238 +++++++++++++++++- Code/Source/solver/fcType.h | 56 ++++- Code/Source/solver/fft.cpp | 66 ----- Code/Source/solver/fft.h | 2 - Code/Source/solver/mat_models.cpp | 6 +- Code/Source/solver/read_files.cpp | 404 +++--------------------------- Code/Source/solver/read_files.h | 6 - Code/Source/solver/set_bc.cpp | 15 +- 9 files changed, 341 insertions(+), 455 deletions(-) diff --git a/Code/Source/solver/bf.cpp b/Code/Source/solver/bf.cpp index 1574b5df0..2832db0ee 100644 --- a/Code/Source/solver/bf.cpp +++ b/Code/Source/solver/bf.cpp @@ -151,8 +151,7 @@ void set_bf_l(ComMod& com_mod, bfType& lBf, mshType& lM, const SolutionStates& s f = lBf.b; } else if (utils::btest(lBf.bType, enum_int(BodyForceType::bfType_ustd))) { - Vector rtmp(1); - ifft(com_mod, lBf.bt, f, rtmp); + f = lBf.bt.value(com_mod.time); } else if (utils::btest(lBf.bType, enum_int(BodyForceType::bfType_gen))) { bfl.resize(idof,nNo); diff --git a/Code/Source/solver/fcType.cpp b/Code/Source/solver/fcType.cpp index 29e2b7959..0197a8cfa 100644 --- a/Code/Source/solver/fcType.cpp +++ b/Code/Source/solver/fcType.cpp @@ -1,4 +1,240 @@ // SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the // University of California, and others. SPDX-License-Identifier: BSD-3-Clause -#include "fcType.h" \ No newline at end of file +#include "fcType.h" +#include "fft.h" + +#include +#include +#include + +namespace { +/// Clean a line by removing carriage returns and leading/trailing spaces, and +/// replacing multiple spaces with a single space. +std::string clean_line(std::string line) { + line.erase(std::remove(line.begin(), line.end(), '\r'), line.end()); + return std::regex_replace(line, std::regex("^ +| +$|( ) +"), "$1"); +} +} // namespace + +fcType fcType::from_time_series_file(const std::string &file_name, + const unsigned int n_dimensions, + const bool is_ramp) { + fcType result; + result.d = n_dimensions; + result.lrmp = is_ramp; + + std::ifstream file(file_name); + + // @todo[michelebucelli] This should actually thrown an exception, ideally of + // a dedicated type such as FileNotFoundException (to be defined). + if (!file.is_open()) { + throw std::runtime_error("Could not open file: " + file_name); + } + + // Read the header of the file. + int n_time_points; + file >> n_time_points >> result.n; + + // @todo[michelebucelli] This should also be an appropriate exception. + if (n_time_points < 2 || result.n == 0) { + throw std::runtime_error( + "Error reading the first line of the temporal values file '" + + file_name + "'."); + } + + if (is_ramp) + result.n = 1; + + // Read the time-value pairs. + std::vector> values; + double tmp; + + std::string line; + int line_number = 1; + + while (std::getline(file, line)) { + line = clean_line(line); + if (line.empty()) + continue; + + std::istringstream line_string_stream(line); + std::vector line_values; + + while (!line_string_stream.eof()) { + line_string_stream >> tmp; + + // @todo[michelebucelli] This should also be an appropriate exception. + if (line_string_stream.fail()) { + throw std::runtime_error( + "Error reading values for the temporal values file '" + file_name + + "' for line " + std::to_string(line_number) + ": '" + line + + "'; value number " + std::to_string(line_values.size() + 1) + + " is not a double."); + } + + line_values.push_back(tmp); + } + + // @todo[michelebucelli] This should also be an appropriate exception. + if (line_values.size() != 1 + n_dimensions) { + throw std::runtime_error( + "Error reading values for the temporal values file '" + file_name + + "' for line " + std::to_string(line_number) + ": '" + line + + "'; expected " + std::to_string(1 + n_dimensions) + + " values, but got " + std::to_string(line_values.size()) + "."); + } + + values.push_back(line_values); + ++line_number; + } + + result.qi.resize(n_dimensions); + result.qs.resize(n_dimensions); + result.r.resize(n_dimensions, result.n); + result.i.resize(n_dimensions, result.n); + + fft(n_time_points, values, result); + + return result; +} + +fcType fcType::from_fourier_coefficients_file(const std::string &file_name, + const unsigned int n_dimensions) { + fcType result; + result.d = n_dimensions; + result.lrmp = false; + + std::ifstream file(file_name); + + // @todo[michelebucelli] This should actually thrown an exception, ideally of + // a dedicated type such as FileNotFoundException (to be defined). + if (!file.is_open()) { + throw std::runtime_error("Could not open file: " + file_name); + } + + // Read the initial time and period from the header. + file >> result.ti >> result.T; + + // Read the linear trend part. + result.qi.resize(n_dimensions); + result.qs.resize(n_dimensions); + + std::string line; + double tmp; + unsigned int current_dimension = 0; + + while (current_dimension < n_dimensions && std::getline(file, line)) { + line = clean_line(line); + if (line.empty()) + continue; + + std::istringstream line_string_stream(line); + std::vector values; + + while (line_string_stream >> tmp) { + values.push_back(tmp); + } + + result.qi[current_dimension] = values[0]; + result.qs[current_dimension] = values[1]; + ++current_dimension; + } + + // Read the Fourier coefficients. + file >> result.n; + result.r.resize(n_dimensions, result.n); + result.i.resize(n_dimensions, result.n); + + unsigned int current_coefficient = 0; + while (current_coefficient < result.n && std::getline(file, line)) { + line = clean_line(line); + if (line.empty()) + continue; + + std::istringstream line_string_stream(line); + std::vector values; + + while (line_string_stream >> tmp) { + values.push_back(tmp); + } + + for (unsigned int i = 0; i < n_dimensions; ++i) { + result.r(i, current_coefficient) = values[i]; + result.i(i, current_coefficient) = values[i + n_dimensions]; + } + + ++current_coefficient; + } + + return result; +} + +Vector fcType::value(const double time) const { + Vector result(d); + static Vector dummy; + + evaluate_internal(time, /* evaluate_derivative = */ false, result, dummy); + + return result; +} + +std::pair, Vector> +fcType::value_and_derivative(const double time) const { + Vector value(d), derivative(d); + + evaluate_internal(time, /* evaluate_derivative = */ true, value, derivative); + + return std::make_pair(value, derivative); +} + +void fcType::evaluate_internal(const double time, + const bool evaluate_derivative, + Vector &value, + Vector &derivative) const { + // @todo[michelebucelli] This should be an appropriate exception. + if (!defined()) { + throw std::runtime_error( + "Cannot evaluate fcType instance that has not been defined."); + } + + // Shifted and rescaled time. + // The input time is shifted by ti. Then, if using the ramp function, it is + // clamped to the interval [0, T]. Otherwise, the time is wrapped to the + // interval [0, T], to enable periodicity. + const double t = + lrmp ? std::max(std::min(time - ti, T), 0.0) : std::fmod(time - ti, T); + + // Linear trend. + for (int i = 0; i < d; ++i) { + value[i] = qi[i] + t * qs[i]; + + if (evaluate_derivative) + derivative[i] = qs[i]; + } + + // Fourier series. + if (!lrmp) { + const double tmp = 2.0 * consts::pi / T; + + // Fourier series. + for (int i = 0; i < n; ++i) { + const double dk = tmp * i; + const double K = t * dk; + + for (int j = 0; j < d; ++j) { + // Using value[j] = value[j] + ... instead of value[j] += ..., because + // the latter changes the order of operations enough to break some of + // the tests. + // @todo[michelebucelli] This seems pretty fragile! + value[j] = + value[j] + r(j, i) * std::cos(K) - this->i(j, i) * std::sin(K); + + if (evaluate_derivative) { + derivative[j] -= + (r(j, i) * std::sin(K) + this->i(j, i) * std::cos(K)) * dk; + } + } + } + } +} \ No newline at end of file diff --git a/Code/Source/solver/fcType.h b/Code/Source/solver/fcType.h index 342103b7f..c1dabe39e 100644 --- a/Code/Source/solver/fcType.h +++ b/Code/Source/solver/fcType.h @@ -7,6 +7,10 @@ #include "Array.h" #include "Vector.h" +#include +#include +#include + /// @brief Fourier coefficients that are used to specify time-dependent /// functions. /// @@ -14,7 +18,31 @@ /// active stress. class fcType { public: - bool defined() { return n != 0; }; + bool defined() const { return n != 0; }; + + /// @brief Read Fourier coefficients from file and return the corresponding + /// fcType instance. + /// + /// @todo[michelebucelli] More detailed documentation, especially on the file + /// format. + static fcType from_fourier_coefficients_file(const std::string &file_name, + const unsigned int n_dimensions); + + /// @brief Read a time series from file and return the corresponding fcType + /// instance. + /// + /// @todo[michelebucelli] More detailed documentation, especially on the file + /// format. + static fcType from_time_series_file(const std::string &file_name, + const unsigned int n_dimensions, + const bool is_ramp); + + /// @brief Return the interpolated value. + Vector value(const double time) const; + + /// @brief Return the interpolated value and time derivative. + std::pair, Vector> + value_and_derivative(const double time) const; /// Toggle whether this is a ramp function or not. bool lrmp = false; @@ -42,6 +70,32 @@ class fcType { /// Real part of coefficients. Array r; + +private: + /** @brief Internal evaluation function. + * + * Uses the inverse Fourier transform to evaluate the value, and optionally + * the derivative, of the interpolated data. This function is not meant to be + * used directly, but only as a backend to @ref value and @ref + * value_and_derivative. + * + * The vectors value and derivative are assumed to be of size @ref d. This is + * not checked by this function. + * + * @throws std::runtime_error if this fcType instance has not been + * initialized (i.e. if @ref defined returns false). + * + * @param[in] time The time at which to evaluate the interpolation. + * @param[in] evaluate_derivative Whether to also evaluate the time + * derivative of the interpolation. + * @param[out] value The interpolated value at the given time. + * @param[out] derivative The time derivative of the interpolated value at + * the given time. If evaluated_derivative is false, this will not be + * modified or accessed. + */ + void evaluate_internal(const double time, const bool evaluate_derivative, + Vector &value, + Vector &derivative) const; }; #endif \ No newline at end of file diff --git a/Code/Source/solver/fft.cpp b/Code/Source/solver/fft.cpp index 5a7c79bf9..4af250798 100644 --- a/Code/Source/solver/fft.cpp +++ b/Code/Source/solver/fft.cpp @@ -93,72 +93,6 @@ void fft(const int np, const std::vector>& temporal_values, } } -/// @brief This is to calculate flow rate and flow acceleration (IFFT) -// -void ifft(const ComMod& com_mod, const fcType& gt, Vector& Y, Vector& dY) -{ - using namespace consts; - #define n_debug_ifft - #ifdef debug_ifft - DebugMsg dmsg(__func__, com_mod.cm.idcm()); - dmsg.banner(); - #endif - - double time = com_mod.time; - #ifdef debug_ifft - dmsg << "time: " << time; - dmsg << "gt.lrmp: " << gt.lrmp; - dmsg << "gt.ti: " << gt.ti; - dmsg << "gt.n: " << gt.n; - dmsg << "pi: " << pi; - #endif - - if (gt.lrmp) { - double t = time - gt.ti; - if (t <= 0.0) { - t = fmax(t, 0.0); - } else { - t = fmin(t, gt.T); - } - #ifdef debug_ifft - dmsg << "t: " << t; - #endif - - for (int i = 0; i < gt.d; i++) { - #ifdef debug_ifft - dmsg << " gt.qi(i): " << gt.qi(i); - dmsg << " gt.qs(i): " << gt.qs(i); - #endif - Y(i) = gt.qi(i) + t*gt.qs(i); - dY(i) = gt.qs(i); - } - - } else { - double t = fmod(time - gt.ti, gt.T); - double tmp = 2.0 * pi / gt.T; - #ifdef debug_ifft - dmsg << "t: " << t; - dmsg << "tmp: " << tmp; - dmsg << "gt.d: " << gt.d; - #endif - - for (int j = 0; j < gt.d; j++) { - Y(j) = gt.qi(j) + t*gt.qs(j); - dY(j) = gt.qs(j); - } - - for (int i = 0; i < gt.n; i++) { - double dk = tmp * static_cast(i); - double K = t*dk; - - for (int j = 0; j < gt.d; j++) { - Y(j) = Y(j) + gt.r(j,i)*cos(K) - gt.i(j,i)*sin(K); - dY(j) = dY(j) - (gt.r(j,i)*sin(K) + gt.i(j,i)*cos(K))*dk; - } - } - } -} - /// @brief This routine is for calculating values by the inverse of general BC // void igbc(const ComMod& com_mod, const MBType& gm, Array& Y, Array& dY) diff --git a/Code/Source/solver/fft.h b/Code/Source/solver/fft.h index 5ffeb0fa2..dbfddeb91 100644 --- a/Code/Source/solver/fft.h +++ b/Code/Source/solver/fft.h @@ -10,8 +10,6 @@ void fft(const int np, const std::vector>& temporal_values, fcType& gt); -void ifft(const ComMod& com_mod, const fcType& gt, Vector& Y, Vector& dY); - void igbc(const ComMod& com_mod, const MBType& gm, Array& Y, Array& dY); #endif diff --git a/Code/Source/solver/mat_models.cpp b/Code/Source/solver/mat_models.cpp index dd8a608ce..58658ac54 100644 --- a/Code/Source/solver/mat_models.cpp +++ b/Code/Source/solver/mat_models.cpp @@ -216,10 +216,8 @@ void compute_fib_stress(const ComMod& com_mod, const CepMod& cep_mod, const fibS if (utils::btest(Tfl.fType, iBC_std)) { g = Tfl.g; - } else if (utils::btest(Tfl.fType, iBC_ustd)) { - Vector gv(1), tv(1); - ifft(com_mod, Tfl.gt, gv, tv); - g = gv[0]; + } else if (utils::btest(Tfl.fType, iBC_ustd)) { + g = Tfl.gt.value(com_mod.time)[0]; } } diff --git a/Code/Source/solver/read_files.cpp b/Code/Source/solver/read_files.cpp index 6c6723bce..01dfc4f9e 100644 --- a/Code/Source/solver/read_files.cpp +++ b/Code/Source/solver/read_files.cpp @@ -222,27 +222,26 @@ void read_bc(Simulation* simulation, EquationParameters* eq_params, eqType& lEq, lBc.g = bc_params->value.value(); } - } else if (ctmp == "Unsteady") { - lBc.bType = utils::ibset(lBc.bType, enum_int(BoundaryConditionType::bType_ustd)); - - if (utils::btest(lBc.bType, enum_int(BoundaryConditionType::bType_trac))) { - lBc.gt.d = com_mod.nsd; - } else { - lBc.gt.d = 1; - } - - if (bc_params->temporal_values_file_path.defined()) { - lBc.gt.lrmp = bc_params->ramp_function.value(); - auto file_name = bc_params->temporal_values_file_path.value(); - read_temporal_values(file_name, lBc); - - } else { - if (!bc_params->fourier_coefficients_file_path.defined()) { - throw std::runtime_error("[read_bc] Undefined data for boundary condition type '" + bc_type + "'."); - } - lBc.gt.lrmp = false; - auto file_name = bc_params->fourier_coefficients_file_path.value(); - read_fourier_coeff_values_file(file_name, lBc); + } else if (ctmp == "Unsteady") { + lBc.bType = + utils::ibset(lBc.bType, enum_int(BoundaryConditionType::bType_ustd)); + + const unsigned int n_dimensions = + utils::btest(lBc.bType, enum_int(BoundaryConditionType::bType_trac)) + ? com_mod.nsd + : 1; + + if (bc_params->temporal_values_file_path.defined()) { + lBc.gt = fcType::from_time_series_file( + bc_params->temporal_values_file_path.value(), n_dimensions, + bc_params->ramp_function.value()); + } else if (bc_params->fourier_coefficients_file_path.defined()) { + lBc.gt = fcType::from_fourier_coefficients_file( + bc_params->fourier_coefficients_file_path.value(), n_dimensions); + } else { + throw std::runtime_error( + "[read_bc] Undefined data for boundary condition type '" + bc_type + + "'."); } // Coupling to a 0D model: @@ -291,7 +290,7 @@ void read_bc(Simulation* simulation, EquationParameters* eq_params, eqType& lEq, } } - } else if (ctmp == "Resistance") { + } else if (ctmp == "Resistance") { lBc.bType = utils::ibset(lBc.bType, enum_int(BoundaryConditionType::bType_res)); if (!utils::btest(lBc.bType, enum_int(BoundaryConditionType::bType_Neu))) { throw std::runtime_error("[read_bc] Resistance is only defined for Neu BC."); @@ -337,7 +336,7 @@ void read_bc(Simulation* simulation, EquationParameters* eq_params, eqType& lEq, // S p a t i a l //// //////////////////////// - } else if (ctmp == "Spatial") { + } else if (ctmp == "Spatial") { lBc.bType = utils::ibset(lBc.bType, enum_int(BoundaryConditionType::bType_gen)); if (!utils::btest(lBc.bType, enum_int(BoundaryConditionType::bType_Neu))) { throw std::runtime_error("[read_bc] Spatial BC is only defined for Neu BC."); @@ -361,7 +360,7 @@ void read_bc(Simulation* simulation, EquationParameters* eq_params, eqType& lEq, // G e n e r a l //// //////////////////////// - } else if (ctmp == "General") { + } else if (ctmp == "General") { lBc.bType = utils::ibset(lBc.bType, enum_int(BoundaryConditionType::bType_gen)); int iM = lBc.iM; int iFa = lBc.iFa; @@ -947,16 +946,18 @@ void read_bf(ComMod& com_mod, BodyForceParameters* bf_params, bfType& lBf) // Unsteady // } else if (time_dependence == "unsteady") { - lBf.bType = utils::ibset(lBf.bType, enum_int(BodyForceType::bfType_ustd)); - lBf.bt.d = lBf.dof; + lBf.bType = utils::ibset(lBf.bType, enum_int(BodyForceType::bfType_ustd)); + if (bf_params->temporal_values_file_path.defined()) { - lBf.bt.lrmp = bf_params->ramp_function.value(); - auto file_name = bf_params->temporal_values_file_path.value(); - read_temporal_values(file_name, lBf); + lBf.bt = fcType::from_time_series_file( + bf_params->temporal_values_file_path.value(), lBf.dof, + bf_params->ramp_function.value()); + } else if (bf_params->fourier_coefficients_file_path.defined()) { + lBf.bt = fcType::from_fourier_coefficients_file( + bf_params->fourier_coefficients_file_path.value(), lBf.dof); } else { - lBf.bt.lrmp = false; - auto fTmp = bf_params->fourier_coefficients_file_path.value(); - read_fourier_coeff_values_file(fTmp, lBf); + throw std::runtime_error("No temporal values or Fourier coefficients " + "file provided for unsteady body force."); } // Spatial // @@ -1580,62 +1581,6 @@ void read_eq(Simulation* simulation, EquationParameters* eq_params, eqType& lEq) } } -//--------------------------------- -// read_fiber_temporal_values_file -//--------------------------------- -// Read fiber reinforcement stress the values from a file. -// -// Note: There is no equivalent Fortran function. -// -void read_fiber_temporal_values_file(FiberReinforcementStressParameters& fiber_params, dmnType& lDmn) -{ - auto file_name = fiber_params.temporal_values_file_path.value(); - std::ifstream temporal_values_file; - temporal_values_file.open(file_name); - if (!temporal_values_file.is_open()) { - throw std::runtime_error("Failed to open the fiber reinforcement stress the values file '" + file_name + "'."); - } - - int i, j; - temporal_values_file >> i >> j; - if ((i == 0) || (j == 0) || (i < 2)) { - throw std::runtime_error("Error reading the first line of the fiber reinforcement stress the values file '" + file_name + "' has an incorrect format."); - } - - lDmn.stM.Tf.gt.d = 1; - lDmn.stM.Tf.gt.n = j; - - // Read time/value pairs. - // - std::vector> temporal_values; - double time, value; - std::string line; - - while (std::getline(temporal_values_file, line)) { - line.erase(std::remove(line.begin(), line.end(), '\r'), line.end()); - if (line == "") { - continue; - } - std::istringstream line_input(line); - std::vector values; - while (line_input >> value) { - values.push_back(value); - } - temporal_values.push_back(values); - } - - if (lDmn.stM.Tf.gt.lrmp) { - lDmn.stM.Tf.gt.n = 1; - } - - lDmn.stM.Tf.gt.qi.resize(lDmn.stM.Tf.gt.d); - lDmn.stM.Tf.gt.qs.resize(lDmn.stM.Tf.gt.d); - lDmn.stM.Tf.gt.r.resize(lDmn.stM.Tf.gt.d,lDmn.stM.Tf.gt.n); - lDmn.stM.Tf.gt.i.resize(lDmn.stM.Tf.gt.d,lDmn.stM.Tf.gt.n); - - fft(i, temporal_values, lDmn.stM.Tf.gt); -} - //------------ // read_files //------------ @@ -1890,151 +1835,6 @@ void read_files(Simulation* simulation, const std::string& file_name) #endif } -//-------------------------------- -// read_fourier_coeff_values_file -//-------------------------------- -// Set boundary condition Fourier coefficients read in from a file. -// -// Note: There is no equivalent Fortran function. -// -// [TODO:DaveP] this is not fully implemented. -// -void read_fourier_coeff_values_file(const std::string& file_name, bcType& lBc) -{ - std::ifstream temporal_values_file; - temporal_values_file.open(file_name); - if (!temporal_values_file.is_open()) { - throw std::runtime_error("Failed to open the Fourier coefficients values file '" + file_name + "'."); - } - - temporal_values_file >> lBc.gt.ti >> lBc.gt.T; - - // Read time/value pairs. - // - lBc.gt.qi.resize(lBc.gt.d); - lBc.gt.qs.resize(lBc.gt.d); - - double time, value; - std::string line; - int n = 0; - - while (std::getline(temporal_values_file, line)) { - line.erase(std::remove(line.begin(), line.end(), '\r'), line.end()); - if (line == "") { - continue; - } - std::istringstream line_input(line); - std::vector values; - while (line_input >> value) { - values.push_back(value); - } - lBc.gt.qi[n] = values[0]; - lBc.gt.qs[n] = values[1]; - n += 1; - - if (n == lBc.gt.d) { - break; - } - } - - temporal_values_file >> lBc.gt.n; - - lBc.gt.r.resize(lBc.gt.d, lBc.gt.n); - lBc.gt.i.resize(lBc.gt.d, lBc.gt.n); - - int j = 0; - - while (std::getline(temporal_values_file, line)) { - line.erase(std::remove(line.begin(), line.end(), '\r'), line.end()); - if (line == "") { - continue; - } - std::istringstream line_input(line); - std::vector values; - while (line_input >> value) { - values.push_back(value); - } - - for (int i = 0; i < lBc.gt.d; i++) { - lBc.gt.r(i,j) = values[i]; - lBc.gt.i(i,j) = values[i + lBc.gt.d]; - } - - j += 1; - if (j == lBc.gt.n) { - break; - } - } -} - -void read_fourier_coeff_values_file(const std::string& file_name, bfType& lBf) -{ - std::ifstream temporal_values_file; - temporal_values_file.open(file_name); - if (!temporal_values_file.is_open()) { - throw std::runtime_error("Failed to open the Fourier coefficients values file '" + file_name + "'."); - } - - temporal_values_file >> lBf.bt.ti >> lBf.bt.T; - - // Read time/value pairs. - // - lBf.bt.qi.resize(lBf.bt.d); - lBf.bt.qs.resize(lBf.bt.d); - - double time, value; - std::string line; - int n = 0; - - while (std::getline(temporal_values_file, line)) { - line.erase(std::remove(line.begin(), line.end(), '\r'), line.end()); - if (line == "") { - continue; - } - std::istringstream line_input(line); - std::vector values; - while (line_input >> value) { - values.push_back(value); - } - lBf.bt.qi[n] = values[0]; - lBf.bt.qs[n] = values[1]; - n += 1; - - if (n == lBf.bt.d) { - break; - } - } - - temporal_values_file >> lBf.bt.n; - - lBf.bt.r.resize(lBf.bt.d, lBf.bt.n); - lBf.bt.i.resize(lBf.bt.d, lBf.bt.n); - - int j = 0; - - while (std::getline(temporal_values_file, line)) { - line.erase(std::remove(line.begin(), line.end(), '\r'), line.end()); - if (line == "") { - continue; - } - std::istringstream line_input(line); - std::vector values; - while (line_input >> value) { - values.push_back(value); - } - - for (int i = 0; i < lBf.bt.d; i++) { - lBf.bt.r(i,j) = values[i]; - lBf.bt.i(i,j) = values[i + lBf.bt.d]; - } - - j += 1; - if (j == lBf.bt.n) { - break; - } - } -} - //--------- // read_ls //--------- @@ -2230,9 +2030,13 @@ void read_mat_model(Simulation* simulation, EquationParameters* eq_params, Domai lDmn.stM.Tf.g = fiber_params.value.value(); } else if (fiber_stress == "unsteady") { - lDmn.stM.Tf.fType = utils::ibset(lDmn.stM.Tf.fType, static_cast(BoundaryConditionType::bType_ustd)); - lDmn.stM.Tf.gt.lrmp = fiber_params.ramp_function.value(); - read_fiber_temporal_values_file(fiber_params, lDmn); + lDmn.stM.Tf.fType = + utils::ibset(lDmn.stM.Tf.fType, + static_cast(BoundaryConditionType::bType_ustd)); + + lDmn.stM.Tf.gt = fcType::from_time_series_file( + fiber_params.temporal_values_file_path.value(), + /* n_dimensions = */ 1, fiber_params.ramp_function.value()); } // Read directional stress distribution parameters @@ -2664,134 +2468,6 @@ void read_temp_spat_values(const ComMod& com_mod, const mshType& msh, const std: } } -//---------------------- -// read_temporal_values -//---------------------- -// Set boundary condition temporal values read in from a file. -// -// Data modified: -// -// Set lBc.gt Fourier coefficients data (fcType) -// lBc.gt.n -// lBc.gt.qi -// lBc.gt.qs -// lBc.gt.r -// lBc.gt.i -// -// Note: There is no equivalent Fortran function. -// -void read_temporal_values(const std::string& file_name, bcType& lBc) -{ - std::ifstream temporal_values_file; - temporal_values_file.open(file_name); - if (!temporal_values_file.is_open()) { - throw std::runtime_error("Failed to open the temporal values file '" + file_name + "'."); - } - - int i, j; - temporal_values_file >> i >> j; - if ((i == 0) || (j == 0) || (i < 2)) { - throw std::runtime_error("Error reading the first line of the temporal values file '" + file_name + "'."); - } - lBc.gt.n = j; - - // Read time/value pairs. - // - std::vector> temporal_values; - double time, value; - std::string line; - int line_number = 1; - int num_values_per_line = lBc.gt.d + 1; - - while (std::getline(temporal_values_file, line)) { - line.erase(std::remove(line.begin(), line.end(), '\r'), line.end()); - if (line == "") { - continue; - } - // Remove leading and trailing spaces. - auto cleaned_line = std::regex_replace(line, std::regex("^ +| +$|( ) +"), "$1"); - std::istringstream line_input(cleaned_line); - std::vector values; - - while (!line_input.eof()) { - line_input >> value; - - if (line_input.fail()) { - throw std::runtime_error("Error reading values for the temporal values file '" + file_name + "' for line " + - std::to_string(line_number) + ": '" + line + "'; value number " + std::to_string(values.size()+1) + " is not a double."); - } - values.push_back(value); - } - - if (values.size() != num_values_per_line) { - throw std::runtime_error("Error reading values for the temporal values file '" + file_name + "' for line " + - std::to_string(line_number) + ": '" + line + "'; expected " + std::to_string(num_values_per_line) + " values per line."); - } - - temporal_values.push_back(values); - line_number += 1; - } - - if (lBc.gt.lrmp) { - lBc.gt.n = 1; - } - - lBc.gt.qi.resize(lBc.gt.d); - lBc.gt.qs.resize(lBc.gt.d); - lBc.gt.r.resize(lBc.gt.d,lBc.gt.n); - lBc.gt.i.resize(lBc.gt.d,lBc.gt.n); - - fft(i, temporal_values, lBc.gt); -} - -// [NOTE] This is a hack, should really just have a single function for this.. -// -void read_temporal_values(const std::string& file_name, bfType& lBf) -{ - std::ifstream temporal_values_file; - temporal_values_file.open(file_name); - if (!temporal_values_file.is_open()) { - throw std::runtime_error("Failed to open the temporal values file '" + file_name + "'."); - } - - int i, j; - temporal_values_file >> i >> j; - if (i < 2) { - throw std::runtime_error("The temporal values file '" + file_name + "' has an incorrect format."); - } - lBf.bt.n = j; - - // Read time/value pairs. - // - std::vector> temporal_values; - double time, value; - std::string line; - - while (std::getline(temporal_values_file, line)) { - line.erase(std::remove(line.begin(), line.end(), '\r'), line.end()); - if (line == "") { - continue; - } - std::istringstream line_input(line); - std::vector values; - while (line_input >> value) { - values.push_back(value); - } - temporal_values.push_back(values); - } - - if (lBf.bt.lrmp) { - lBf.bt.n = 1; - } - - lBf.bt.qi.resize(lBf.bt.d); - lBf.bt.qs.resize(lBf.bt.d); - lBf.bt.r.resize(lBf.bt.d, lBf.bt.n); - lBf.bt.i.resize(lBf.bt.d, lBf.bt.n); - - fft(i, temporal_values, lBf.bt); -} - //---------------- // read_trac_bcff //---------------- diff --git a/Code/Source/solver/read_files.h b/Code/Source/solver/read_files.h index 3cccfc469..197f83fdb 100644 --- a/Code/Source/solver/read_files.h +++ b/Code/Source/solver/read_files.h @@ -35,9 +35,6 @@ namespace read_files_ns { void read_files(Simulation* simulation, const std::string& file_name); - void read_fourier_coeff_values_file(const std::string& file_name, bcType& lBc); - void read_fourier_coeff_values_file(const std::string& file_name, bfType& lBf); - void read_ls(Simulation* simulation, EquationParameters* eq_params, consts::SolverType solver_type, eqType& lEq); void read_mat_model(Simulation* simulation, EquationParameters* eq_params, DomainParameters* domain_params, dmnType& lDmn); @@ -48,9 +45,6 @@ namespace read_files_ns { void read_spatial_values(const ComMod& com_mod, const mshType& msh, const faceType& lFa, const std::string& file_name, bcType& lBc); - void read_temporal_values(const std::string& file_name, bcType& lBc); - void read_temporal_values(const std::string& file_name, bfType& lBf); - void read_temp_spat_values(const ComMod& com_mod, const mshType& msh, const faceType& lFa, const std::string& file_name, bcType& lBc); void read_temp_spat_values(const ComMod& com_mod, const mshType& msh, const std::string& file_name, bfType& lBf); diff --git a/Code/Source/solver/set_bc.cpp b/Code/Source/solver/set_bc.cpp index 01d361ceb..300a365be 100644 --- a/Code/Source/solver/set_bc.cpp +++ b/Code/Source/solver/set_bc.cpp @@ -1075,10 +1075,9 @@ void set_bc_dir_l(ComMod& com_mod, const bcType& lBc, const faceType& lFa, Array #ifdef debug_set_bc_dir_l dmsg << "bType_ustd"; #endif - Vector dirY_v(1), dirA_v(1); - ifft(com_mod, lBc.gt, dirY_v, dirA_v); - dirY = dirY_v(0); - dirA = dirA_v(0); + const auto value_and_derivative = lBc.gt.value_and_derivative(com_mod.time); + dirY = value_and_derivative.first[0]; + dirA = value_and_derivative.second[0]; } else { dirA = 0.0; @@ -1427,7 +1426,7 @@ void set_bc_neu_l(ComMod& com_mod, const CmMod& cm_mod, const bcType& lBc, const int nsd = com_mod.nsd; int nNo = lFa.nNo; - Vector h(1), rtmp(1); + Vector h(1); Vector tmpA(nNo); // Geting the contribution of Neu BC @@ -1466,7 +1465,7 @@ void set_bc_neu_l(ComMod& com_mod, const CmMod& cm_mod, const bcType& lBc, const h(0) = lBc.g; } else if (utils::btest(lBc.bType,iBC_ustd)) { - ifft(com_mod, lBc.gt, h, rtmp); + h = lBc.gt.value(com_mod.time); } else { throw std::runtime_error("[set_bc_neu_l] Correction in SETBCNEU is needed"); @@ -1830,9 +1829,7 @@ void set_bc_trac_l(ComMod& com_mod, const CmMod& cm_mod, const bcType& lBc, cons if (lBc.gt.d != nsd) { throw std::runtime_error("[set_bc_trac_l] Traction dof not initialized properly"); } - Vector h(nsd); - Vector tmp(nsd); - ifft(com_mod, lBc.gt, h, tmp); + const Vector h = lBc.gt.value(com_mod.time); for (int a = 0; a < nNo; a++) { int Ac = lFa.gN(a); hg.set_col(Ac, h); From ba5234c5ba691151672407f8a2decead39e2efef Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Fri, 26 Jun 2026 13:23:01 -0500 Subject: [PATCH 03/19] Distribution of fcType members encapsulated in fcType::distribute --- Code/Source/solver/distribute.cpp | 85 ++----------------------------- Code/Source/solver/fcType.cpp | 30 +++++++++++ Code/Source/solver/fcType.h | 7 +++ 3 files changed, 41 insertions(+), 81 deletions(-) diff --git a/Code/Source/solver/distribute.cpp b/Code/Source/solver/distribute.cpp index 64acce552..742f01d31 100644 --- a/Code/Source/solver/distribute.cpp +++ b/Code/Source/solver/distribute.cpp @@ -644,44 +644,10 @@ void dist_bc(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, bcType& lBc } // Communicating time-dependent BC data - // - // lBc.gt is declare ALLOCATABLE in MOD.f but we don't - // want to use pointers so use the define() method - // to check if it has data define. - // - bool flag = lBc.gt.defined(); - cm.bcast(cm_mod, &flag); - - if (flag) { - if (is_slave) { - // [NOTE] This is allocated in ComMod. - //lBc.gt = new fcType; - } - - cm.bcast(cm_mod, &lBc.gt.lrmp); - cm.bcast(cm_mod, &lBc.gt.d); - cm.bcast(cm_mod, &lBc.gt.n); - - int j = lBc.gt.d; - int i = lBc.gt.n; - - if (is_slave) { - lBc.gt.qi.resize(j); - lBc.gt.qs.resize(j); - lBc.gt.r.resize(j,i); - lBc.gt.i.resize(j,i); - } - - cm.bcast(cm_mod, &lBc.gt.ti); - cm.bcast(cm_mod, &lBc.gt.T); - cm.bcast(cm_mod, lBc.gt.qi); - cm.bcast(cm_mod, lBc.gt.qs); - cm.bcast(cm_mod, lBc.gt.r); - cm.bcast(cm_mod, lBc.gt.i); - } + lBc.gt.distribute(cm_mod, cm); // Communicating moving BC data - flag = lBc.gm.defined(); + bool flag = lBc.gm.defined(); cm.bcast(cm_mod, &flag); if (flag) { @@ -917,31 +883,7 @@ void dist_bf(ComMod& com_mod, const CmMod& cm_mod, const cmType& cm, bfType& lBf } // Communicating time-dependent BF data - flag = lBf.bt.defined(); - cm.bcast(cm_mod, &flag); - if (flag) { - if (is_slave) { - //ALLOCATE(lBf.bt) - } - cm.bcast(cm_mod, &lBf.bt.lrmp); - cm.bcast(cm_mod, &lBf.bt.d); - cm.bcast(cm_mod, &lBf.bt.n); - - if (is_slave) { - int j = lBf.bt.d; - int i = lBf.bt.n; - lBf.bt.qi.resize(j); - lBf.bt.qs.resize(j); - lBf.bt.r.resize(j,i); - lBf.bt.i.resize(j,i); - } - cm.bcast(cm_mod, &lBf.bt.ti); - cm.bcast(cm_mod, &lBf.bt.T); - cm.bcast(cm_mod, lBf.bt.qi); - cm.bcast(cm_mod, lBf.bt.qs); - cm.bcast(cm_mod, lBf.bt.r); - cm.bcast(cm_mod, lBf.bt.i); - } + lBf.bt.distribute(cm_mod, cm); // Communicating moving BF data // @@ -1731,26 +1673,7 @@ void dist_mat_consts(const ComMod& com_mod, const CmMod& cm_mod, const cmType& c cm.bcast(cm_mod, &lStM.Tf.g); } else if (utils::btest(lStM.Tf.fType, static_cast(BoundaryConditionType::bType_ustd))) { - cm.bcast(cm_mod, &lStM.Tf.gt.lrmp); - cm.bcast(cm_mod, &lStM.Tf.gt.d); - cm.bcast(cm_mod, &lStM.Tf.gt.n); - - if (cm.slv(cm_mod)) { - int j = lStM.Tf.gt.d; - int i = lStM.Tf.gt.n; - lStM.Tf.gt.qi.resize(j); - lStM.Tf.gt.qs.resize(j); - lStM.Tf.gt.r.resize(j,i); - lStM.Tf.gt.i.resize(j,i); - } - - cm.bcast(cm_mod, &lStM.Tf.gt.ti); - cm.bcast(cm_mod, &lStM.Tf.gt.T); - - cm.bcast(cm_mod, lStM.Tf.gt.qi, "lStM.Tf.gt.qi"); - cm.bcast(cm_mod, lStM.Tf.gt.qs, "lStM.Tf.gt.qs"); - cm.bcast(cm_mod, lStM.Tf.gt.r, "lStM.Tf.gt.r"); - cm.bcast(cm_mod, lStM.Tf.gt.i, "lStM.Tf.gt.i"); + lStM.Tf.gt.distribute(cm_mod, cm); } // Broadcast directional stress distribution parameters diff --git a/Code/Source/solver/fcType.cpp b/Code/Source/solver/fcType.cpp index 0197a8cfa..5a1cb8e0c 100644 --- a/Code/Source/solver/fcType.cpp +++ b/Code/Source/solver/fcType.cpp @@ -170,6 +170,36 @@ fcType fcType::from_fourier_coefficients_file(const std::string &file_name, return result; } +void fcType::distribute(const CmMod &cm_mod, const cmType &cm) { + // Only the master knows whether the object has been initialized. Therefore, + // we broadcast the initialization flag to all ranks, so that the following if + // statement can be run correctly by all. + bool initialized = defined(); + cm.bcast(cm_mod, &initialized); + + if (initialized) { + cm.bcast(cm_mod, &lrmp); + cm.bcast(cm_mod, &n); + cm.bcast(cm_mod, &d); + + // All ranks but the master need to allocate the arrays before receiving + // data. + if (cm.slv(cm_mod)) { + qi.resize(d); + qs.resize(d); + r.resize(d, n); + i.resize(d, n); + } + + cm.bcast(cm_mod, &ti); + cm.bcast(cm_mod, &T); + cm.bcast(cm_mod, qi); + cm.bcast(cm_mod, qs); + cm.bcast(cm_mod, r); + cm.bcast(cm_mod, i); + } +} + Vector fcType::value(const double time) const { Vector result(d); static Vector dummy; diff --git a/Code/Source/solver/fcType.h b/Code/Source/solver/fcType.h index c1dabe39e..a64aaa7be 100644 --- a/Code/Source/solver/fcType.h +++ b/Code/Source/solver/fcType.h @@ -5,6 +5,7 @@ #define FC_TYPE_H #include "Array.h" +#include "CmMod.h" #include "Vector.h" #include @@ -37,6 +38,12 @@ class fcType { const unsigned int n_dimensions, const bool is_ramp); + /// @brief Distribute the data to all parallel processes. + /// + /// Does nothing if the object has not been initialized (i.e. if @ref defined + /// returns false) on the master rank. + void distribute(const CmMod &cm_mod, const cmType &cm); + /// @brief Return the interpolated value. Vector value(const double time) const; From 6f460abd2a664cc02570c29939e2da89c363e07b Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Fri, 26 Jun 2026 13:32:30 -0500 Subject: [PATCH 04/19] Introduce getter method fcType::dimension --- Code/Source/solver/fcType.h | 8 ++++++-- Code/Source/solver/set_bc.cpp | 7 ++++--- 2 files changed, 10 insertions(+), 5 deletions(-) diff --git a/Code/Source/solver/fcType.h b/Code/Source/solver/fcType.h index a64aaa7be..e668f40bc 100644 --- a/Code/Source/solver/fcType.h +++ b/Code/Source/solver/fcType.h @@ -19,8 +19,6 @@ /// active stress. class fcType { public: - bool defined() const { return n != 0; }; - /// @brief Read Fourier coefficients from file and return the corresponding /// fcType instance. /// @@ -51,6 +49,12 @@ class fcType { std::pair, Vector> value_and_derivative(const double time) const; + /// @brief Return whether this object has been initialized. + bool defined() const { return n != 0; }; + + /// @brief Get the dimension of the data interpolated by this object. + unsigned int dimension() const { return d; } + /// Toggle whether this is a ramp function or not. bool lrmp = false; diff --git a/Code/Source/solver/set_bc.cpp b/Code/Source/solver/set_bc.cpp index 300a365be..89e691a60 100644 --- a/Code/Source/solver/set_bc.cpp +++ b/Code/Source/solver/set_bc.cpp @@ -1826,9 +1826,10 @@ void set_bc_trac_l(ComMod& com_mod, const CmMod& cm_mod, const bcType& lBc, cons } } else if (utils::btest(lBc.bType,iBC_ustd)) { - if (lBc.gt.d != nsd) { - throw std::runtime_error("[set_bc_trac_l] Traction dof not initialized properly"); - } + if (lBc.gt.dimension() != nsd) { + throw std::runtime_error( + "[set_bc_trac_l] Traction dof not initialized properly"); + } const Vector h = lBc.gt.value(com_mod.time); for (int a = 0; a < nNo; a++) { int Ac = lFa.gN(a); From bc6813c6aa3d4f29d18f8860d4a8e3f71da241fe Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Fri, 26 Jun 2026 14:32:37 -0500 Subject: [PATCH 05/19] Move fft into fcType, and cleaner separation of fcType input from file and construction --- Code/Source/solver/fcType.cpp | 173 +++++++++++++++++++----- Code/Source/solver/fcType.h | 29 +++- Code/Source/solver/fft.cpp | 93 +------------ Code/Source/solver/fft.h | 2 - tests/unitTests/math_tests/test_fft.cpp | 3 +- 5 files changed, 168 insertions(+), 132 deletions(-) diff --git a/Code/Source/solver/fcType.cpp b/Code/Source/solver/fcType.cpp index 5a1cb8e0c..c40bbb4b4 100644 --- a/Code/Source/solver/fcType.cpp +++ b/Code/Source/solver/fcType.cpp @@ -17,13 +17,111 @@ std::string clean_line(std::string line) { } } // namespace -fcType fcType::from_time_series_file(const std::string &file_name, - const unsigned int n_dimensions, - const bool is_ramp) { +fcType fcType::from_time_series( + const unsigned int n_fourier_coefficients, + const std::vector> &temporal_values, + const bool is_ramp) { + const unsigned int n_time_points = temporal_values.size(); + + // @todo[michelebucelli] This should be an appropriate exception. + if (n_time_points < 2) { + throw std::runtime_error( + "At least two time points are needed to construct an fcType object."); + } + + const unsigned int n_dimensions = temporal_values[0].size() - 1; + + // Check that all entries of temporal_values have the same number of elements. + // @todo[michelebucelli] This should probably be an Array, in which this is + // enforced a priori, rather than a vector of vectors. This change would + // also save us the extraction for loop below. + for (const auto &row : temporal_values) { + if (row.size() != n_dimensions + 1) { + throw std::runtime_error( + "All rows of temporal_values must have the same number of elements."); + } + } + + Vector times(n_time_points); + Array values(n_dimensions, n_time_points); + + for (unsigned int i = 0; i < n_time_points; ++i) { + times[i] = temporal_values[i][0]; + for (unsigned int j = 0; j < n_dimensions; ++j) { + values(j, i) = temporal_values[i][j + 1]; + } + } + fcType result; - result.d = n_dimensions; + result.lrmp = is_ramp; + result.d = n_dimensions; + result.n = is_ramp ? 1 : n_fourier_coefficients; + result.ti = times[0]; + result.T = times[n_time_points - 1] - times[0]; + result.qi.resize(n_dimensions); + result.qs.resize(n_dimensions); + result.r.resize(n_dimensions, result.n); + result.i.resize(n_dimensions, result.n); + + // Compute the linear trend part. + for (unsigned int j = 0; j < n_dimensions; ++j) { + result.qi[j] = values(j, 0); + result.qs[j] = (values(j, n_time_points - 1) - values(j, 0)) / result.T; + } + + // Subtract the linear trend part from the values. + for (unsigned int i = 0; i < n_time_points; ++i) { + times[i] -= result.ti; + for (unsigned int j = 0; j < n_dimensions; ++j) { + values(j, i) = values(j, i) - result.qi[j] - result.qs[j] * times[i]; + } + } + + // Compute the Fourier coefficients. + for (int n = 0; n < result.n; ++n) { + const double tmp = static_cast(n); + result.r.set_col(n, 0.0); + result.i.set_col(n, 0.0); + + for (int i = 0; i < n_time_points - 1; ++i) { + const double ko = 2.0 * consts::pi * tmp * times[i] / result.T; + const double kn = 2.0 * consts::pi * tmp * times[i + 1] / result.T; + + for (int j = 0; j < result.d; j++) { + const double s = + (values(j, i + 1) - values(j, i)) / (times[i + 1] - times[i]); + + if (n == 0) { + result.r(j, n) += 0.5 * (times[i + 1] - times[i]) * + (values(j, i + 1) + values(j, i)); + } else { + result.r(j, n) += s * (std::cos(kn) - std::cos(ko)); + result.i(j, n) -= s * (std::sin(kn) - std::sin(ko)); + } + } + } + + if (n == 0) { + for (int k = 0; k < result.d; k++) { + result.r(k, n) /= result.T; + } + } else { + for (int k = 0; k < result.d; k++) { + result.r(k, n) = 0.5 * result.r(k, n) * result.T / + (consts::pi * consts::pi * tmp * tmp); + result.i(k, n) = 0.5 * result.i(k, n) * result.T / + (consts::pi * consts::pi * tmp * tmp); + } + } + } + + return result; +} +fcType fcType::from_time_series_file(const std::string &file_name, + const unsigned int n_dimensions, + const bool is_ramp) { std::ifstream file(file_name); // @todo[michelebucelli] This should actually thrown an exception, ideally of @@ -33,19 +131,16 @@ fcType fcType::from_time_series_file(const std::string &file_name, } // Read the header of the file. - int n_time_points; - file >> n_time_points >> result.n; + int n_time_points, n_fourier_coefficients; + file >> n_time_points >> n_fourier_coefficients; // @todo[michelebucelli] This should also be an appropriate exception. - if (n_time_points < 2 || result.n == 0) { + if (n_time_points < 2 || n_fourier_coefficients == 0) { throw std::runtime_error( "Error reading the first line of the temporal values file '" + file_name + "'."); } - if (is_ramp) - result.n = 1; - // Read the time-value pairs. std::vector> values; double tmp; @@ -89,22 +184,33 @@ fcType fcType::from_time_series_file(const std::string &file_name, ++line_number; } - result.qi.resize(n_dimensions); - result.qs.resize(n_dimensions); - result.r.resize(n_dimensions, result.n); - result.i.resize(n_dimensions, result.n); + return fcType::from_time_series(n_fourier_coefficients, values, is_ramp); +} + +fcType fcType::from_fourier_coefficients(const Vector &qi, + const Vector &qs, + const Array &r, + const Array &i, + const double ti, const double T) { + // @todo[michelebucelli] Add some correctness checks on the input arguments. + + fcType result; - fft(n_time_points, values, result); + result.lrmp = false; + result.n = r.ncols(); + result.d = qi.size(); + result.qi = qi; + result.qs = qs; + result.r = r; + result.i = i; + result.ti = ti; + result.T = T; return result; } fcType fcType::from_fourier_coefficients_file(const std::string &file_name, const unsigned int n_dimensions) { - fcType result; - result.d = n_dimensions; - result.lrmp = false; - std::ifstream file(file_name); // @todo[michelebucelli] This should actually thrown an exception, ideally of @@ -114,11 +220,11 @@ fcType fcType::from_fourier_coefficients_file(const std::string &file_name, } // Read the initial time and period from the header. - file >> result.ti >> result.T; + double ti, T; + file >> ti >> T; // Read the linear trend part. - result.qi.resize(n_dimensions); - result.qs.resize(n_dimensions); + Vector qi(n_dimensions), qs(n_dimensions); std::string line; double tmp; @@ -136,18 +242,21 @@ fcType fcType::from_fourier_coefficients_file(const std::string &file_name, values.push_back(tmp); } - result.qi[current_dimension] = values[0]; - result.qs[current_dimension] = values[1]; + qi[current_dimension] = values[0]; + qs[current_dimension] = values[1]; ++current_dimension; } // Read the Fourier coefficients. - file >> result.n; - result.r.resize(n_dimensions, result.n); - result.i.resize(n_dimensions, result.n); + unsigned int n_fourier_coefficients; + file >> n_fourier_coefficients; + + Array r(n_dimensions, n_fourier_coefficients); + Array i(n_dimensions, n_fourier_coefficients); unsigned int current_coefficient = 0; - while (current_coefficient < result.n && std::getline(file, line)) { + while (current_coefficient < n_fourier_coefficients && + std::getline(file, line)) { line = clean_line(line); if (line.empty()) continue; @@ -159,15 +268,15 @@ fcType fcType::from_fourier_coefficients_file(const std::string &file_name, values.push_back(tmp); } - for (unsigned int i = 0; i < n_dimensions; ++i) { - result.r(i, current_coefficient) = values[i]; - result.i(i, current_coefficient) = values[i + n_dimensions]; + for (unsigned int j = 0; j < n_dimensions; ++j) { + r(j, current_coefficient) = values[j]; + i(j, current_coefficient) = values[j + n_dimensions]; } ++current_coefficient; } - return result; + return fcType::from_fourier_coefficients(qi, qs, r, i, ti, T); } void fcType::distribute(const CmMod &cm_mod, const cmType &cm) { diff --git a/Code/Source/solver/fcType.h b/Code/Source/solver/fcType.h index e668f40bc..0fe3e9ad1 100644 --- a/Code/Source/solver/fcType.h +++ b/Code/Source/solver/fcType.h @@ -17,15 +17,17 @@ /// /// These are used for boundary conditions and for time-dependent prescribed /// active stress. +/// +/// @todo[michelebucelli] More detailed documentation. class fcType { public: - /// @brief Read Fourier coefficients from file and return the corresponding - /// fcType instance. + /// @brief Construct an fcType from a time series. /// - /// @todo[michelebucelli] More detailed documentation, especially on the file - /// format. - static fcType from_fourier_coefficients_file(const std::string &file_name, - const unsigned int n_dimensions); + /// @todo[michelebucelli] More detailed documentation. + static fcType + from_time_series(const unsigned int n_fourier_coefficients, + const std::vector> &temporal_values, + const bool is_ramp); /// @brief Read a time series from file and return the corresponding fcType /// instance. @@ -36,6 +38,21 @@ class fcType { const unsigned int n_dimensions, const bool is_ramp); + /// @brief Construct an fcType instance from Fourier coefficients. + static fcType from_fourier_coefficients(const Vector &qi, + const Vector &qs, + const Array &r, + const Array &i, + const double ti, const double T); + + /// @brief Read Fourier coefficients from file and return the corresponding + /// fcType instance. + /// + /// @todo[michelebucelli] More detailed documentation, especially on the file + /// format. + static fcType from_fourier_coefficients_file(const std::string &file_name, + const unsigned int n_dimensions); + /// @brief Distribute the data to all parallel processes. /// /// Does nothing if the object has not been initialized (i.e. if @ref defined diff --git a/Code/Source/solver/fft.cpp b/Code/Source/solver/fft.cpp index 4af250798..5a362ebd0 100644 --- a/Code/Source/solver/fft.cpp +++ b/Code/Source/solver/fft.cpp @@ -6,93 +6,6 @@ #include "fft.h" #include -/// @brief Replicates Fortran 'SUBROUTINE FFT(fid, np, gt)'. -/// -/// temporal_values contains all of the values read in from -/// the temporal values file. -// -void fft(const int np, const std::vector>& temporal_values, fcType& gt) -{ - using namespace consts; - - #define n_debug_fft - #ifdef debug_fft - DebugMsg dmsg(__func__, 0); - dmsg.banner(); - dmsg << "np: " << np; - #endif - - Vector t(np); - Array q(gt.d, np); - - for (int i = 0; i < np; i++) { - t[i] = temporal_values[i][0]; - for (int j = 0; j < gt.d; j++) { - q(j,i) = temporal_values[i][j+1]; - } - } - - - gt.ti = t[0]; - gt.T = t[np-1] - t[0]; - #ifdef debug_fft - dmsg << "pi: " << pi; - dmsg << "gt.ti: " << gt.ti; - dmsg << "gt.T: " << gt.T; - dmsg << "gt.d: " << gt.d; - dmsg << "gt.n: " << gt.n; - dmsg << "np: " << np; - dmsg << "gt.r.nrows(): " << gt.r.nrows(); - #endif - - for (int j = 0; j < gt.d; j++) { - gt.qi[j] = q(j,0); - gt.qs[j] = (q(j,np-1) - q(j,0)) / gt.T; - } - - for (int i = 0; i < np; i++) { - t[i] = t[i] - gt.ti; - for (int j = 0; j < gt.d; j++) { - q(j,i) = q(j,i) - gt.qi[j] - gt.qs[j]*t[i]; - } - } - - int ns = 0; - Vector ns_array(512); - - for (int n = 0; n < gt.n; n++) { - double tmp = static_cast(n); - gt.r.set_col(n, 0.0); - gt.i.set_col(n, 0.0); - - for (int i = 0; i < np-1; i++) { - double ko = 2.0 * pi * tmp * t[i] / gt.T; - double kn = 2.0 * pi * tmp * t[i+1] / gt.T; - - for (int j = 0; j < gt.d; j++) { - double s = (q(j,i+1) - q(j,i)) / (t[i+1] - t[i]); - if (n == 0) { - gt.r(j,n) = gt.r(j,n) + 0.5*(t[i+1] - t[i]) * (q(j,i+1) + q(j,i)); - } else { - gt.r(j,n) = gt.r(j,n) + s*(cos(kn) - cos(ko)); - gt.i(j,n) = gt.i(j,n) - s*(sin(kn) - sin(ko)); - } - } - } - - if (n == 0) { - for (int k = 0; k < gt.d; k++) { - gt.r(k,n) = gt.r(k,n) / gt.T; - } - } else { - for (int k = 0; k < gt.d; k++) { - gt.r(k,n) = 0.5 * gt.r(k,n) * gt.T / (pi*pi*tmp*tmp); - gt.i(k,n) = 0.5 * gt.i(k,n) * gt.T / (pi*pi*tmp*tmp); - } - } - } -} - /// @brief This routine is for calculating values by the inverse of general BC // void igbc(const ComMod& com_mod, const MBType& gm, Array& Y, Array& dY) @@ -117,7 +30,5 @@ void igbc(const ComMod& com_mod, const MBType& gm, Array& Y, Array -void fft(const int np, const std::vector>& temporal_values, fcType& gt); - void igbc(const ComMod& com_mod, const MBType& gm, Array& Y, Array& dY); #endif diff --git a/tests/unitTests/math_tests/test_fft.cpp b/tests/unitTests/math_tests/test_fft.cpp index 044472f29..f4d3dd1db 100644 --- a/tests/unitTests/math_tests/test_fft.cpp +++ b/tests/unitTests/math_tests/test_fft.cpp @@ -78,7 +78,8 @@ TEST_F(FFTTest, SinCosLinearCombination) { InitializeFourierCoefficients(gt, 1, 16); // Compute the Fourier coefficients - fft(N, temporal_values, gt); + gt = fcType::from_time_series(/* n_fourier_coefficients = */ N, + temporal_values, /* is_ramp = */ false); // Check the slope (first Fourier coefficient) ASSERT_NEAR(gt.qs[0], -0.13830, 1e-2) << "Expected slope ~-0.13830"; From 474031619bba6ef842462b635d3e51387c670a51 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Fri, 26 Jun 2026 14:41:30 -0500 Subject: [PATCH 06/19] Make data members of fcType private --- Code/Source/solver/fcType.h | 84 +++++++++++++++++-------- tests/unitTests/math_tests/test_fft.cpp | 34 ++++------ 2 files changed, 69 insertions(+), 49 deletions(-) diff --git a/Code/Source/solver/fcType.h b/Code/Source/solver/fcType.h index 0fe3e9ad1..c83a56929 100644 --- a/Code/Source/solver/fcType.h +++ b/Code/Source/solver/fcType.h @@ -66,38 +66,43 @@ class fcType { std::pair, Vector> value_and_derivative(const double time) const; + /// @name Data member access. + /// @{ + /// @brief Return whether this object has been initialized. bool defined() const { return n != 0; }; /// @brief Get the dimension of the data interpolated by this object. unsigned int dimension() const { return d; } - /// Toggle whether this is a ramp function or not. - bool lrmp = false; - - /// Number of Fourier coefficients. - int n = 0; - - /// Numbero of dimensions (scalar or vector). - int d = 0; - - /// Initial value. - Vector qi; - - /// Time derivative of linear part. - Vector qs; - - /// Period. - double T = 0.0; - - /// Initial time. - double ti = 0.0; - - /// Imaginary part of coefficients. - Array i; - - /// Real part of coefficients. - Array r; + /// @brief Get the initial value of the linear trend part for one component. + const double get_qi(const unsigned int component) const { + // @todo[michelebucelli] Add check on the index. + return qi[component]; + } + + /// @brief Get the slope of the linear trend part for one component. + const double get_qs(const unsigned int component) const { + // @todo[michelebucelli] Add check on the index. + return qs[component]; + } + + /// @brief Get the real part of the Fourier coefficients for one component. + const double get_r(const unsigned int component, + const unsigned int frequency) const { + // @todo[michelebucelli] Add check on the indices. + return r(component, frequency); + } + + /// @brief Get the imaginary part of the Fourier coefficients for one + /// component. + const double get_i(const unsigned int component, + const unsigned int frequency) const { + // @todo[michelebucelli] Add check on the indices. + return i(component, frequency); + } + + /// @} private: /** @brief Internal evaluation function. @@ -124,6 +129,33 @@ class fcType { void evaluate_internal(const double time, const bool evaluate_derivative, Vector &value, Vector &derivative) const; + + /// Toggle whether this is a ramp function or not. + bool lrmp = false; + + /// Number of Fourier coefficients. + int n = 0; + + /// Numbero of dimensions (scalar or vector). + int d = 0; + + /// Initial value. + Vector qi; + + /// Time derivative of linear part. + Vector qs; + + /// Period. + double T = 0.0; + + /// Initial time. + double ti = 0.0; + + /// Imaginary part of coefficients. + Array i; + + /// Real part of coefficients. + Array r; }; #endif \ No newline at end of file diff --git a/tests/unitTests/math_tests/test_fft.cpp b/tests/unitTests/math_tests/test_fft.cpp index f4d3dd1db..554180577 100644 --- a/tests/unitTests/math_tests/test_fft.cpp +++ b/tests/unitTests/math_tests/test_fft.cpp @@ -52,15 +52,6 @@ class FFTTest : public ::testing::Test { temporal_values.push_back({t, y}); } } - - void InitializeFourierCoefficients(fcType& gt, int d = 1, int n = 16) { - gt.d = d; - gt.n = n; - gt.qi.resize(gt.d); - gt.qs.resize(gt.d); - gt.r.resize(gt.d, gt.n); - gt.i.resize(gt.d, gt.n); - } }; TEST_F(FFTTest, SinCosLinearCombination) { @@ -72,23 +63,20 @@ TEST_F(FFTTest, SinCosLinearCombination) { std::vector> temporal_values; CreateTemporalValues(N, x_start, x_end, temporal_values); - - fcType gt; - // Apply interpolation to 1 dimensional data with n = 16 Fourier modes - InitializeFourierCoefficients(gt, 1, 16); - + // Compute the Fourier coefficients - gt = fcType::from_time_series(/* n_fourier_coefficients = */ N, - temporal_values, /* is_ramp = */ false); + const fcType gt = + fcType::from_time_series(/* n_fourier_coefficients = */ 16, + temporal_values, /* is_ramp = */ false); // Check the slope (first Fourier coefficient) - ASSERT_NEAR(gt.qs[0], -0.13830, 1e-2) << "Expected slope ~-0.13830"; + ASSERT_NEAR(gt.get_qs(0), -0.13830, 1e-2) << "Expected slope ~-0.13830"; // Check the real and imaginary components of the first three Fourier coefficients - ASSERT_NEAR(gt.r(0, 0), 0.32094, 1e-2) << "Expected first real coefficient to be close to 0.32094"; - ASSERT_NEAR(gt.i(0, 0), 0.0, 1e-2) << "Expected first imaginary coefficient to be close to 0.0"; - ASSERT_NEAR(gt.r(0, 1), 0.42759, 1e-2) << "Expected second real coefficient to be close to 0.42759"; - ASSERT_NEAR(gt.i(0, 1), 1.25295, 1e-2) << "Expected second imaginary coefficient to be close to 1.25295"; - ASSERT_NEAR(gt.r(0, 2), -0.44685, 1e-2) << "Expected third real coefficient to be close to -0.44685"; - ASSERT_NEAR(gt.i(0, 2), -0.65403, 1e-2) << "Expected third imaginary coefficient to be close to -0.65403"; + ASSERT_NEAR(gt.get_r(0, 0), 0.32094, 1e-2) << "Expected first real coefficient to be close to 0.32094"; + ASSERT_NEAR(gt.get_i(0, 0), 0.0, 1e-2) << "Expected first imaginary coefficient to be close to 0.0"; + ASSERT_NEAR(gt.get_r(0, 1), 0.42759, 1e-2) << "Expected second real coefficient to be close to 0.42759"; + ASSERT_NEAR(gt.get_i(0, 1), 1.25295, 1e-2) << "Expected second imaginary coefficient to be close to 1.25295"; + ASSERT_NEAR(gt.get_r(0, 2), -0.44685, 1e-2) << "Expected third real coefficient to be close to -0.44685"; + ASSERT_NEAR(gt.get_i(0, 2), -0.65403, 1e-2) << "Expected third imaginary coefficient to be close to -0.65403"; } \ No newline at end of file From a371f0aad777130c12cce07f9f1dbdb4739a3a27 Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Fri, 26 Jun 2026 16:42:52 -0500 Subject: [PATCH 07/19] Add cmType::bcast overload for unsigned integers --- Code/Source/solver/CmMod.cpp | 5 +++++ Code/Source/solver/CmMod.h | 3 +++ 2 files changed, 8 insertions(+) diff --git a/Code/Source/solver/CmMod.cpp b/Code/Source/solver/CmMod.cpp index f3485a2bb..2fbfb1c86 100644 --- a/Code/Source/solver/CmMod.cpp +++ b/Code/Source/solver/CmMod.cpp @@ -138,6 +138,11 @@ void cmType::bcast(const CmMod& cm_mod, Vector& data) const MPI_Bcast(data.data(), data.size(), cm_mod::mpint, cm_mod.master, com()); } +/// @brief bcast unsigned int +void cmType::bcast(const CmMod &cm_mod, unsigned int *data) const { + MPI_Bcast(data, 1, cm_mod::mpuint, cm_mod.master, com()); +} + /// @brief gather int array void cmType::gather(const CmMod& cm_mod, const int* send_data, int send_count, int* recv_data, int recv_count, int root) const { diff --git a/Code/Source/solver/CmMod.h b/Code/Source/solver/CmMod.h index 36ac60e18..5d2587b67 100644 --- a/Code/Source/solver/CmMod.h +++ b/Code/Source/solver/CmMod.h @@ -24,6 +24,7 @@ using MpiCommWorldType = MPI_Comm; const decltype(MPI_CXX_BOOL) mplog = MPI_CXX_BOOL; //const decltype(MPI_LOGICAL) mplog = MPI_LOGICAL; const decltype(MPI_INTEGER) mpint = MPI_INTEGER; +const decltype(MPI_UNSIGNED) mpuint = MPI_UNSIGNED; const decltype(MPI_DOUBLE_PRECISION) mpreal = MPI_DOUBLE_PRECISION; const decltype(MPI_CHARACTER) mpchar = MPI_CHARACTER; }; @@ -83,6 +84,8 @@ class cmType { void bcast(const CmMod& cm_mod, int* data) const; void bcast(const CmMod& cm_mod, Vector& data) const; + void bcast(const CmMod &cm_mod, unsigned int *data) const; + void bcast(const CmMod& cm_mod, Array& data, const std::string& name="") const; // Gather operations From e3f4436d0762e3543090613a9b8fa598bd8361fb Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Fri, 26 Jun 2026 16:46:28 -0500 Subject: [PATCH 08/19] Rename fcType to FourierInterpolation, and rename its data members for clarity --- Code/Source/solver/CMakeLists.txt | 2 +- Code/Source/solver/ComMod.h | 14 +- Code/Source/solver/fcType.cpp | 379 ----------------- Code/Source/solver/fcType.h | 161 -------- Code/Source/solver/fourier_interpolation.cpp | 404 +++++++++++++++++++ Code/Source/solver/fourier_interpolation.h | 203 ++++++++++ Code/Source/solver/read_files.cpp | 10 +- Code/Source/solver/set_bc.cpp | 2 +- tests/unitTests/math_tests/test_fft.cpp | 23 +- tests/unitTests/math_tests/test_fft.h | 10 +- 10 files changed, 635 insertions(+), 573 deletions(-) delete mode 100644 Code/Source/solver/fcType.cpp delete mode 100644 Code/Source/solver/fcType.h create mode 100644 Code/Source/solver/fourier_interpolation.cpp create mode 100644 Code/Source/solver/fourier_interpolation.h diff --git a/Code/Source/solver/CMakeLists.txt b/Code/Source/solver/CMakeLists.txt index c0e299a0d..30294bf16 100644 --- a/Code/Source/solver/CMakeLists.txt +++ b/Code/Source/solver/CMakeLists.txt @@ -157,7 +157,7 @@ set(CSRCS TrilinosLinearAlgebra.h TrilinosLinearAlgebra.cpp Tensor4.h Tensor4.cpp Vector.h Vector.cpp - fcType.cpp + fourier_interpolation.cpp lapack_defs.h diff --git a/Code/Source/solver/ComMod.h b/Code/Source/solver/ComMod.h index 4ca4114ad..1d00fd230 100644 --- a/Code/Source/solver/ComMod.h +++ b/Code/Source/solver/ComMod.h @@ -8,20 +8,20 @@ // are defined here. #ifndef COMMOD_H -#define COMMOD_H +#define COMMOD_H #include "Array.h" #include "Array3.h" -#include "SolutionStates.h" #include "CepMod.h" #include "ChnlMod.h" #include "CmMod.h" +#include "CoupledBoundaryCondition.h" #include "Parameters.h" #include "RobinBoundaryCondition.h" -#include "CoupledBoundaryCondition.h" +#include "SolutionStates.h" #include "Timer.h" #include "Vector.h" -#include "fcType.h" +#include "fourier_interpolation.h" #include "DebugMsg.h" @@ -159,7 +159,7 @@ class bcType // // This is declare ALLOCATABLE in MOD.f. // - fcType gt; + FourierInterpolation gt; // Neu: RCR rcrType RCR; @@ -271,7 +271,7 @@ class bfType Array bx; // Time dependant (unsteady imposed value) - fcType bt; + FourierInterpolation bt; // General (unsteady and spatially dependent combination) MBType bm; @@ -299,7 +299,7 @@ class fibStrsType double eta_n = 0.0; // Unsteady time-dependent values - fcType gt; + FourierInterpolation gt; }; /// @brief Structural domain type diff --git a/Code/Source/solver/fcType.cpp b/Code/Source/solver/fcType.cpp deleted file mode 100644 index c40bbb4b4..000000000 --- a/Code/Source/solver/fcType.cpp +++ /dev/null @@ -1,379 +0,0 @@ -// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the -// University of California, and others. SPDX-License-Identifier: BSD-3-Clause - -#include "fcType.h" -#include "fft.h" - -#include -#include -#include - -namespace { -/// Clean a line by removing carriage returns and leading/trailing spaces, and -/// replacing multiple spaces with a single space. -std::string clean_line(std::string line) { - line.erase(std::remove(line.begin(), line.end(), '\r'), line.end()); - return std::regex_replace(line, std::regex("^ +| +$|( ) +"), "$1"); -} -} // namespace - -fcType fcType::from_time_series( - const unsigned int n_fourier_coefficients, - const std::vector> &temporal_values, - const bool is_ramp) { - const unsigned int n_time_points = temporal_values.size(); - - // @todo[michelebucelli] This should be an appropriate exception. - if (n_time_points < 2) { - throw std::runtime_error( - "At least two time points are needed to construct an fcType object."); - } - - const unsigned int n_dimensions = temporal_values[0].size() - 1; - - // Check that all entries of temporal_values have the same number of elements. - // @todo[michelebucelli] This should probably be an Array, in which this is - // enforced a priori, rather than a vector of vectors. This change would - // also save us the extraction for loop below. - for (const auto &row : temporal_values) { - if (row.size() != n_dimensions + 1) { - throw std::runtime_error( - "All rows of temporal_values must have the same number of elements."); - } - } - - Vector times(n_time_points); - Array values(n_dimensions, n_time_points); - - for (unsigned int i = 0; i < n_time_points; ++i) { - times[i] = temporal_values[i][0]; - for (unsigned int j = 0; j < n_dimensions; ++j) { - values(j, i) = temporal_values[i][j + 1]; - } - } - - fcType result; - - result.lrmp = is_ramp; - result.d = n_dimensions; - result.n = is_ramp ? 1 : n_fourier_coefficients; - result.ti = times[0]; - result.T = times[n_time_points - 1] - times[0]; - result.qi.resize(n_dimensions); - result.qs.resize(n_dimensions); - result.r.resize(n_dimensions, result.n); - result.i.resize(n_dimensions, result.n); - - // Compute the linear trend part. - for (unsigned int j = 0; j < n_dimensions; ++j) { - result.qi[j] = values(j, 0); - result.qs[j] = (values(j, n_time_points - 1) - values(j, 0)) / result.T; - } - - // Subtract the linear trend part from the values. - for (unsigned int i = 0; i < n_time_points; ++i) { - times[i] -= result.ti; - for (unsigned int j = 0; j < n_dimensions; ++j) { - values(j, i) = values(j, i) - result.qi[j] - result.qs[j] * times[i]; - } - } - - // Compute the Fourier coefficients. - for (int n = 0; n < result.n; ++n) { - const double tmp = static_cast(n); - result.r.set_col(n, 0.0); - result.i.set_col(n, 0.0); - - for (int i = 0; i < n_time_points - 1; ++i) { - const double ko = 2.0 * consts::pi * tmp * times[i] / result.T; - const double kn = 2.0 * consts::pi * tmp * times[i + 1] / result.T; - - for (int j = 0; j < result.d; j++) { - const double s = - (values(j, i + 1) - values(j, i)) / (times[i + 1] - times[i]); - - if (n == 0) { - result.r(j, n) += 0.5 * (times[i + 1] - times[i]) * - (values(j, i + 1) + values(j, i)); - } else { - result.r(j, n) += s * (std::cos(kn) - std::cos(ko)); - result.i(j, n) -= s * (std::sin(kn) - std::sin(ko)); - } - } - } - - if (n == 0) { - for (int k = 0; k < result.d; k++) { - result.r(k, n) /= result.T; - } - } else { - for (int k = 0; k < result.d; k++) { - result.r(k, n) = 0.5 * result.r(k, n) * result.T / - (consts::pi * consts::pi * tmp * tmp); - result.i(k, n) = 0.5 * result.i(k, n) * result.T / - (consts::pi * consts::pi * tmp * tmp); - } - } - } - - return result; -} - -fcType fcType::from_time_series_file(const std::string &file_name, - const unsigned int n_dimensions, - const bool is_ramp) { - std::ifstream file(file_name); - - // @todo[michelebucelli] This should actually thrown an exception, ideally of - // a dedicated type such as FileNotFoundException (to be defined). - if (!file.is_open()) { - throw std::runtime_error("Could not open file: " + file_name); - } - - // Read the header of the file. - int n_time_points, n_fourier_coefficients; - file >> n_time_points >> n_fourier_coefficients; - - // @todo[michelebucelli] This should also be an appropriate exception. - if (n_time_points < 2 || n_fourier_coefficients == 0) { - throw std::runtime_error( - "Error reading the first line of the temporal values file '" + - file_name + "'."); - } - - // Read the time-value pairs. - std::vector> values; - double tmp; - - std::string line; - int line_number = 1; - - while (std::getline(file, line)) { - line = clean_line(line); - if (line.empty()) - continue; - - std::istringstream line_string_stream(line); - std::vector line_values; - - while (!line_string_stream.eof()) { - line_string_stream >> tmp; - - // @todo[michelebucelli] This should also be an appropriate exception. - if (line_string_stream.fail()) { - throw std::runtime_error( - "Error reading values for the temporal values file '" + file_name + - "' for line " + std::to_string(line_number) + ": '" + line + - "'; value number " + std::to_string(line_values.size() + 1) + - " is not a double."); - } - - line_values.push_back(tmp); - } - - // @todo[michelebucelli] This should also be an appropriate exception. - if (line_values.size() != 1 + n_dimensions) { - throw std::runtime_error( - "Error reading values for the temporal values file '" + file_name + - "' for line " + std::to_string(line_number) + ": '" + line + - "'; expected " + std::to_string(1 + n_dimensions) + - " values, but got " + std::to_string(line_values.size()) + "."); - } - - values.push_back(line_values); - ++line_number; - } - - return fcType::from_time_series(n_fourier_coefficients, values, is_ramp); -} - -fcType fcType::from_fourier_coefficients(const Vector &qi, - const Vector &qs, - const Array &r, - const Array &i, - const double ti, const double T) { - // @todo[michelebucelli] Add some correctness checks on the input arguments. - - fcType result; - - result.lrmp = false; - result.n = r.ncols(); - result.d = qi.size(); - result.qi = qi; - result.qs = qs; - result.r = r; - result.i = i; - result.ti = ti; - result.T = T; - - return result; -} - -fcType fcType::from_fourier_coefficients_file(const std::string &file_name, - const unsigned int n_dimensions) { - std::ifstream file(file_name); - - // @todo[michelebucelli] This should actually thrown an exception, ideally of - // a dedicated type such as FileNotFoundException (to be defined). - if (!file.is_open()) { - throw std::runtime_error("Could not open file: " + file_name); - } - - // Read the initial time and period from the header. - double ti, T; - file >> ti >> T; - - // Read the linear trend part. - Vector qi(n_dimensions), qs(n_dimensions); - - std::string line; - double tmp; - unsigned int current_dimension = 0; - - while (current_dimension < n_dimensions && std::getline(file, line)) { - line = clean_line(line); - if (line.empty()) - continue; - - std::istringstream line_string_stream(line); - std::vector values; - - while (line_string_stream >> tmp) { - values.push_back(tmp); - } - - qi[current_dimension] = values[0]; - qs[current_dimension] = values[1]; - ++current_dimension; - } - - // Read the Fourier coefficients. - unsigned int n_fourier_coefficients; - file >> n_fourier_coefficients; - - Array r(n_dimensions, n_fourier_coefficients); - Array i(n_dimensions, n_fourier_coefficients); - - unsigned int current_coefficient = 0; - while (current_coefficient < n_fourier_coefficients && - std::getline(file, line)) { - line = clean_line(line); - if (line.empty()) - continue; - - std::istringstream line_string_stream(line); - std::vector values; - - while (line_string_stream >> tmp) { - values.push_back(tmp); - } - - for (unsigned int j = 0; j < n_dimensions; ++j) { - r(j, current_coefficient) = values[j]; - i(j, current_coefficient) = values[j + n_dimensions]; - } - - ++current_coefficient; - } - - return fcType::from_fourier_coefficients(qi, qs, r, i, ti, T); -} - -void fcType::distribute(const CmMod &cm_mod, const cmType &cm) { - // Only the master knows whether the object has been initialized. Therefore, - // we broadcast the initialization flag to all ranks, so that the following if - // statement can be run correctly by all. - bool initialized = defined(); - cm.bcast(cm_mod, &initialized); - - if (initialized) { - cm.bcast(cm_mod, &lrmp); - cm.bcast(cm_mod, &n); - cm.bcast(cm_mod, &d); - - // All ranks but the master need to allocate the arrays before receiving - // data. - if (cm.slv(cm_mod)) { - qi.resize(d); - qs.resize(d); - r.resize(d, n); - i.resize(d, n); - } - - cm.bcast(cm_mod, &ti); - cm.bcast(cm_mod, &T); - cm.bcast(cm_mod, qi); - cm.bcast(cm_mod, qs); - cm.bcast(cm_mod, r); - cm.bcast(cm_mod, i); - } -} - -Vector fcType::value(const double time) const { - Vector result(d); - static Vector dummy; - - evaluate_internal(time, /* evaluate_derivative = */ false, result, dummy); - - return result; -} - -std::pair, Vector> -fcType::value_and_derivative(const double time) const { - Vector value(d), derivative(d); - - evaluate_internal(time, /* evaluate_derivative = */ true, value, derivative); - - return std::make_pair(value, derivative); -} - -void fcType::evaluate_internal(const double time, - const bool evaluate_derivative, - Vector &value, - Vector &derivative) const { - // @todo[michelebucelli] This should be an appropriate exception. - if (!defined()) { - throw std::runtime_error( - "Cannot evaluate fcType instance that has not been defined."); - } - - // Shifted and rescaled time. - // The input time is shifted by ti. Then, if using the ramp function, it is - // clamped to the interval [0, T]. Otherwise, the time is wrapped to the - // interval [0, T], to enable periodicity. - const double t = - lrmp ? std::max(std::min(time - ti, T), 0.0) : std::fmod(time - ti, T); - - // Linear trend. - for (int i = 0; i < d; ++i) { - value[i] = qi[i] + t * qs[i]; - - if (evaluate_derivative) - derivative[i] = qs[i]; - } - - // Fourier series. - if (!lrmp) { - const double tmp = 2.0 * consts::pi / T; - - // Fourier series. - for (int i = 0; i < n; ++i) { - const double dk = tmp * i; - const double K = t * dk; - - for (int j = 0; j < d; ++j) { - // Using value[j] = value[j] + ... instead of value[j] += ..., because - // the latter changes the order of operations enough to break some of - // the tests. - // @todo[michelebucelli] This seems pretty fragile! - value[j] = - value[j] + r(j, i) * std::cos(K) - this->i(j, i) * std::sin(K); - - if (evaluate_derivative) { - derivative[j] -= - (r(j, i) * std::sin(K) + this->i(j, i) * std::cos(K)) * dk; - } - } - } - } -} \ No newline at end of file diff --git a/Code/Source/solver/fcType.h b/Code/Source/solver/fcType.h deleted file mode 100644 index c83a56929..000000000 --- a/Code/Source/solver/fcType.h +++ /dev/null @@ -1,161 +0,0 @@ -// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the -// University of California, and others. SPDX-License-Identifier: BSD-3-Clause - -#ifndef FC_TYPE_H -#define FC_TYPE_H - -#include "Array.h" -#include "CmMod.h" -#include "Vector.h" - -#include -#include -#include - -/// @brief Fourier coefficients that are used to specify time-dependent -/// functions. -/// -/// These are used for boundary conditions and for time-dependent prescribed -/// active stress. -/// -/// @todo[michelebucelli] More detailed documentation. -class fcType { -public: - /// @brief Construct an fcType from a time series. - /// - /// @todo[michelebucelli] More detailed documentation. - static fcType - from_time_series(const unsigned int n_fourier_coefficients, - const std::vector> &temporal_values, - const bool is_ramp); - - /// @brief Read a time series from file and return the corresponding fcType - /// instance. - /// - /// @todo[michelebucelli] More detailed documentation, especially on the file - /// format. - static fcType from_time_series_file(const std::string &file_name, - const unsigned int n_dimensions, - const bool is_ramp); - - /// @brief Construct an fcType instance from Fourier coefficients. - static fcType from_fourier_coefficients(const Vector &qi, - const Vector &qs, - const Array &r, - const Array &i, - const double ti, const double T); - - /// @brief Read Fourier coefficients from file and return the corresponding - /// fcType instance. - /// - /// @todo[michelebucelli] More detailed documentation, especially on the file - /// format. - static fcType from_fourier_coefficients_file(const std::string &file_name, - const unsigned int n_dimensions); - - /// @brief Distribute the data to all parallel processes. - /// - /// Does nothing if the object has not been initialized (i.e. if @ref defined - /// returns false) on the master rank. - void distribute(const CmMod &cm_mod, const cmType &cm); - - /// @brief Return the interpolated value. - Vector value(const double time) const; - - /// @brief Return the interpolated value and time derivative. - std::pair, Vector> - value_and_derivative(const double time) const; - - /// @name Data member access. - /// @{ - - /// @brief Return whether this object has been initialized. - bool defined() const { return n != 0; }; - - /// @brief Get the dimension of the data interpolated by this object. - unsigned int dimension() const { return d; } - - /// @brief Get the initial value of the linear trend part for one component. - const double get_qi(const unsigned int component) const { - // @todo[michelebucelli] Add check on the index. - return qi[component]; - } - - /// @brief Get the slope of the linear trend part for one component. - const double get_qs(const unsigned int component) const { - // @todo[michelebucelli] Add check on the index. - return qs[component]; - } - - /// @brief Get the real part of the Fourier coefficients for one component. - const double get_r(const unsigned int component, - const unsigned int frequency) const { - // @todo[michelebucelli] Add check on the indices. - return r(component, frequency); - } - - /// @brief Get the imaginary part of the Fourier coefficients for one - /// component. - const double get_i(const unsigned int component, - const unsigned int frequency) const { - // @todo[michelebucelli] Add check on the indices. - return i(component, frequency); - } - - /// @} - -private: - /** @brief Internal evaluation function. - * - * Uses the inverse Fourier transform to evaluate the value, and optionally - * the derivative, of the interpolated data. This function is not meant to be - * used directly, but only as a backend to @ref value and @ref - * value_and_derivative. - * - * The vectors value and derivative are assumed to be of size @ref d. This is - * not checked by this function. - * - * @throws std::runtime_error if this fcType instance has not been - * initialized (i.e. if @ref defined returns false). - * - * @param[in] time The time at which to evaluate the interpolation. - * @param[in] evaluate_derivative Whether to also evaluate the time - * derivative of the interpolation. - * @param[out] value The interpolated value at the given time. - * @param[out] derivative The time derivative of the interpolated value at - * the given time. If evaluated_derivative is false, this will not be - * modified or accessed. - */ - void evaluate_internal(const double time, const bool evaluate_derivative, - Vector &value, - Vector &derivative) const; - - /// Toggle whether this is a ramp function or not. - bool lrmp = false; - - /// Number of Fourier coefficients. - int n = 0; - - /// Numbero of dimensions (scalar or vector). - int d = 0; - - /// Initial value. - Vector qi; - - /// Time derivative of linear part. - Vector qs; - - /// Period. - double T = 0.0; - - /// Initial time. - double ti = 0.0; - - /// Imaginary part of coefficients. - Array i; - - /// Real part of coefficients. - Array r; -}; - -#endif \ No newline at end of file diff --git a/Code/Source/solver/fourier_interpolation.cpp b/Code/Source/solver/fourier_interpolation.cpp new file mode 100644 index 000000000..51bb207bf --- /dev/null +++ b/Code/Source/solver/fourier_interpolation.cpp @@ -0,0 +1,404 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the +// University of California, and others. SPDX-License-Identifier: BSD-3-Clause + +#include "fourier_interpolation.h" +#include "fft.h" + +#include +#include +#include + +namespace { +/// Clean a line by removing carriage returns and leading/trailing spaces, and +/// replacing multiple spaces with a single space. +std::string clean_line(std::string line) { + line.erase(std::remove(line.begin(), line.end(), '\r'), line.end()); + return std::regex_replace(line, std::regex("^ +| +$|( ) +"), "$1"); +} +} // namespace + +FourierInterpolation FourierInterpolation::from_time_series( + const unsigned int n_fourier_coefficients, + const std::vector> &temporal_values, + const bool use_ramp) { + const unsigned int n_time_points = temporal_values.size(); + + // @todo[michelebucelli] This should be an appropriate exception. + if (n_time_points < 2) { + throw std::runtime_error("At least two time points are needed to construct " + "a FourierInterpolation object."); + } + + const unsigned int n_components = temporal_values[0].size() - 1; + + // Check that all entries of temporal_values have the same number of elements. + // @todo[michelebucelli] This should probably be an Array, in which this is + // enforced a priori, rather than a vector of vectors. This change would + // also save us the extraction for loop below. + for (const auto &row : temporal_values) { + if (row.size() != n_components + 1) { + throw std::runtime_error( + "All rows of temporal_values must have the same number of elements."); + } + } + + Vector times(n_time_points); + Array values(n_components, n_time_points); + + for (unsigned int i = 0; i < n_time_points; ++i) { + times[i] = temporal_values[i][0]; + for (unsigned int j = 0; j < n_components; ++j) { + values(j, i) = temporal_values[i][j + 1]; + } + } + + FourierInterpolation result; + + result.use_ramp = use_ramp; + result.n_components = n_components; + result.n_fourier_coefficients = use_ramp ? 1 : n_fourier_coefficients; + result.initial_time = times[0]; + result.period = times[n_time_points - 1] - times[0]; + result.linear_trend_initial_values.resize(n_components); + result.linear_trend_slopes.resize(n_components); + result.fourier_coefficients_real.resize(n_components, + result.n_fourier_coefficients); + result.fourier_coefficients_imaginary.resize(n_components, + result.n_fourier_coefficients); + + // Compute the linear trend part. + for (unsigned int j = 0; j < n_components; ++j) { + result.linear_trend_initial_values[j] = values(j, 0); + result.linear_trend_slopes[j] = + (values(j, n_time_points - 1) - values(j, 0)) / result.period; + } + + // Subtract the linear trend part from the values. + for (unsigned int i = 0; i < n_time_points; ++i) { + times[i] -= result.initial_time; + for (unsigned int j = 0; j < n_components; ++j) { + values(j, i) = values(j, i) - result.linear_trend_initial_values[j] - + result.linear_trend_slopes[j] * times[i]; + } + } + + // Compute the Fourier coefficients. + for (int n = 0; n < result.n_fourier_coefficients; ++n) { + const double tmp = static_cast(n); + result.fourier_coefficients_real.set_col(n, 0.0); + result.fourier_coefficients_imaginary.set_col(n, 0.0); + + for (int i = 0; i < n_time_points - 1; ++i) { + const double ko = 2.0 * consts::pi * tmp * times[i] / result.period; + const double kn = 2.0 * consts::pi * tmp * times[i + 1] / result.period; + + for (int j = 0; j < result.n_components; j++) { + const double s = + (values(j, i + 1) - values(j, i)) / (times[i + 1] - times[i]); + + if (n == 0) { + result.fourier_coefficients_real(j, n) += + 0.5 * (times[i + 1] - times[i]) * + (values(j, i + 1) + values(j, i)); + } else { + result.fourier_coefficients_real(j, n) += + s * (std::cos(kn) - std::cos(ko)); + result.fourier_coefficients_imaginary(j, n) -= + s * (std::sin(kn) - std::sin(ko)); + } + } + } + + if (n == 0) { + for (int k = 0; k < result.n_components; k++) { + result.fourier_coefficients_real(k, n) /= result.period; + } + } else { + const double tmp_2 = (consts::pi * consts::pi * tmp * tmp); + + for (int k = 0; k < result.n_components; k++) { + // @todo[michelebucelli] Check if operator*= breaks something here. + result.fourier_coefficients_real(k, n) = + 0.5 * result.fourier_coefficients_real(k, n) * result.period / + tmp_2; + result.fourier_coefficients_imaginary(k, n) = + 0.5 * result.fourier_coefficients_imaginary(k, n) * result.period / + tmp_2; + } + } + } + + return result; +} + +FourierInterpolation +FourierInterpolation::from_time_series_file(const std::string &file_name, + const unsigned int n_components, + const bool use_ramp) { + std::ifstream file(file_name); + + // @todo[michelebucelli] This should actually thrown an exception, ideally of + // a dedicated type such as FileNotFoundException (to be defined). + if (!file.is_open()) { + throw std::runtime_error("Could not open file: " + file_name); + } + + // Read the header of the file. + int n_time_points, n_fourier_coefficients; + file >> n_time_points >> n_fourier_coefficients; + + // @todo[michelebucelli] This should also be an appropriate exception. + if (n_time_points < 2 || n_fourier_coefficients == 0) { + throw std::runtime_error( + "Error reading the first line of the temporal values file '" + + file_name + "'."); + } + + // Read the time-value pairs. + std::vector> values; + double tmp; + + std::string line; + int line_number = 1; + + while (std::getline(file, line)) { + line = clean_line(line); + if (line.empty()) + continue; + + std::istringstream line_string_stream(line); + std::vector line_values; + + while (!line_string_stream.eof()) { + line_string_stream >> tmp; + + // @todo[michelebucelli] This should also be an appropriate exception. + if (line_string_stream.fail()) { + throw std::runtime_error( + "Error reading values for the temporal values file '" + file_name + + "' for line " + std::to_string(line_number) + ": '" + line + + "'; value number " + std::to_string(line_values.size() + 1) + + " is not a double."); + } + + line_values.push_back(tmp); + } + + // @todo[michelebucelli] This should also be an appropriate exception. + if (line_values.size() != 1 + n_components) { + throw std::runtime_error( + "Error reading values for the temporal values file '" + file_name + + "' for line " + std::to_string(line_number) + ": '" + line + + "'; expected " + std::to_string(1 + n_components) + + " values, but got " + std::to_string(line_values.size()) + "."); + } + + values.push_back(line_values); + ++line_number; + } + + return FourierInterpolation::from_time_series(n_fourier_coefficients, values, + use_ramp); +} + +FourierInterpolation FourierInterpolation::from_fourier_coefficients( + const Vector &linear_trend_initial_values, + const Vector &linear_trend_slopes, + const Array &fourier_coefficients_real, + const Array &fourier_coefficients_imaginary, + const double initial_time, const double period) { + // @todo[michelebucelli] Add some correctness checks on the input arguments. + + FourierInterpolation result; + + result.use_ramp = false; + result.n_fourier_coefficients = fourier_coefficients_real.ncols(); + result.n_components = linear_trend_initial_values.size(); + result.linear_trend_initial_values = linear_trend_initial_values; + result.linear_trend_slopes = linear_trend_slopes; + result.fourier_coefficients_real = fourier_coefficients_real; + result.fourier_coefficients_imaginary = fourier_coefficients_imaginary; + result.initial_time = initial_time; + result.period = period; + + return result; +} + +FourierInterpolation FourierInterpolation::from_fourier_coefficients_file( + const std::string &file_name, const unsigned int n_components) { + std::ifstream file(file_name); + + // @todo[michelebucelli] This should actually thrown an exception, ideally of + // a dedicated type such as FileNotFoundException (to be defined). + if (!file.is_open()) { + throw std::runtime_error("Could not open file: " + file_name); + } + + double initial_time, period; + file >> initial_time >> period; + + // Read the linear trend part. + Vector linear_trend_initial_values(n_components); + Vector linear_trend_slopes(n_components); + + std::string line; + double tmp; + unsigned int current_component = 0; + + while (current_component < n_components && std::getline(file, line)) { + line = clean_line(line); + if (line.empty()) + continue; + + std::istringstream line_string_stream(line); + std::vector values; + + while (line_string_stream >> tmp) { + values.push_back(tmp); + } + + linear_trend_initial_values[current_component] = values[0]; + linear_trend_slopes[current_component] = values[1]; + ++current_component; + } + + // Read the Fourier coefficients. + unsigned int n_fourier_coefficients; + file >> n_fourier_coefficients; + + Array fourier_coefficients_real(n_components, n_fourier_coefficients); + Array fourier_coefficients_imaginary(n_components, + n_fourier_coefficients); + + unsigned int current_coefficient = 0; + while (current_coefficient < n_fourier_coefficients && + std::getline(file, line)) { + line = clean_line(line); + if (line.empty()) + continue; + + std::istringstream line_string_stream(line); + std::vector values; + + while (line_string_stream >> tmp) { + values.push_back(tmp); + } + + for (unsigned int j = 0; j < n_components; ++j) { + fourier_coefficients_real(j, current_coefficient) = values[j]; + fourier_coefficients_imaginary(j, current_coefficient) = + values[j + n_components]; + } + + ++current_coefficient; + } + + return FourierInterpolation::from_fourier_coefficients( + linear_trend_initial_values, linear_trend_slopes, + fourier_coefficients_real, fourier_coefficients_imaginary, initial_time, + period); +} + +void FourierInterpolation::distribute(const CmMod &cm_mod, const cmType &cm) { + // Only the master knows whether the object has been initialized. Therefore, + // we broadcast the initialization flag to all ranks, so that the following if + // statement can be run correctly by all. + bool initialized = defined(); + cm.bcast(cm_mod, &initialized); + + if (initialized) { + cm.bcast(cm_mod, &use_ramp); + cm.bcast(cm_mod, &n_fourier_coefficients); + cm.bcast(cm_mod, &n_components); + + // All ranks but the master need to allocate the arrays before receiving + // data. + if (cm.slv(cm_mod)) { + linear_trend_initial_values.resize(n_components); + linear_trend_slopes.resize(n_components); + fourier_coefficients_real.resize(n_components, n_fourier_coefficients); + fourier_coefficients_imaginary.resize(n_components, + n_fourier_coefficients); + } + + cm.bcast(cm_mod, &initial_time); + cm.bcast(cm_mod, &period); + cm.bcast(cm_mod, linear_trend_initial_values); + cm.bcast(cm_mod, linear_trend_slopes); + cm.bcast(cm_mod, fourier_coefficients_real); + cm.bcast(cm_mod, fourier_coefficients_imaginary); + } +} + +Vector FourierInterpolation::value(const double time) const { + Vector result(n_components); + static Vector dummy; + + evaluate_internal(time, /* evaluate_derivative = */ false, result, dummy); + + return result; +} + +std::pair, Vector> +FourierInterpolation::value_and_derivative(const double time) const { + Vector value(n_components); + Vector derivative(n_components); + + evaluate_internal(time, /* evaluate_derivative = */ true, value, derivative); + + return std::make_pair(value, derivative); +} + +void FourierInterpolation::evaluate_internal(const double time, + const bool evaluate_derivative, + Vector &value, + Vector &derivative) const { + // @todo[michelebucelli] This should be an appropriate exception. + if (!defined()) { + throw std::runtime_error("Cannot evaluate FourierInterpolation instance " + "that has not been defined."); + } + + // Shifted and rescaled time. + // The input time is shifted by ti. Then, if using the ramp function, it is + // clamped to the interval [0, T]. Otherwise, the time is wrapped to the + // interval [0, T], to enable periodicity. + const double t = use_ramp + ? std::max(std::min(time - initial_time, period), 0.0) + : std::fmod(time - initial_time, period); + + // Linear trend. + for (int i = 0; i < n_components; ++i) { + value[i] = linear_trend_initial_values[i] + t * linear_trend_slopes[i]; + + if (evaluate_derivative) + derivative[i] = linear_trend_slopes[i]; + } + + // Fourier series. + if (!use_ramp) { + const double tmp = 2.0 * consts::pi / period; + + // Fourier series. + for (int i = 0; i < n_fourier_coefficients; ++i) { + const double dk = tmp * i; + const double K = t * dk; + + for (int j = 0; j < n_components; ++j) { + // Using value[j] = value[j] + ... instead of value[j] += ..., because + // the latter changes the order of operations enough to break some of + // the tests. + // @todo[michelebucelli] This seems pretty fragile! + value[j] = value[j] + fourier_coefficients_real(j, i) * std::cos(K) - + fourier_coefficients_imaginary(j, i) * std::sin(K); + + if (evaluate_derivative) { + derivative[j] -= + (fourier_coefficients_real(j, i) * std::sin(K) + + fourier_coefficients_imaginary(j, i) * std::cos(K)) * + dk; + } + } + } + } +} \ No newline at end of file diff --git a/Code/Source/solver/fourier_interpolation.h b/Code/Source/solver/fourier_interpolation.h new file mode 100644 index 000000000..b6199e5e6 --- /dev/null +++ b/Code/Source/solver/fourier_interpolation.h @@ -0,0 +1,203 @@ +// SPDX-FileCopyrightText: Copyright (c) Stanford University, The Regents of the +// University of California, and others. SPDX-License-Identifier: BSD-3-Clause + +#ifndef FOURIER_INTERPOLATION_H +#define FOURIER_INTERPOLATION_H + +#include "Array.h" +#include "CmMod.h" +#include "Vector.h" + +#include +#include +#include + +/// @brief Fourier coefficients that are used to specify time-dependent +/// functions. +/// +/// These are used for boundary conditions and for time-dependent prescribed +/// active stress. +/// +/// @todo[michelebucelli] More detailed documentation. +class FourierInterpolation { +public: + /// @brief Construct an FourierInterpolation from a time series. + /// + /// @todo[michelebucelli] More detailed documentation. + static FourierInterpolation + from_time_series(const unsigned int n_fourier_coefficients, + const std::vector> &temporal_values, + const bool use_ramp); + + /// @brief Read a time series from file and return the corresponding + /// FourierInterpolation instance. + /// + /// @todo[michelebucelli] More detailed documentation, especially on the file + /// format. + static FourierInterpolation + from_time_series_file(const std::string &file_name, + const unsigned int n_components, const bool use_ramp); + + /// @brief Construct an FourierInterpolation instance from Fourier + /// coefficients. + static FourierInterpolation + from_fourier_coefficients(const Vector &linear_trend_initial_values, + const Vector &linear_trend_slopes, + const Array &fourier_coefficients_real, + const Array &fourier_coefficients_imaginary, + const double initial_time, const double period); + + /// @brief Read Fourier coefficients from file and return the corresponding + /// FourierInterpolation instance. + /// + /// @todo[michelebucelli] More detailed documentation, especially on the file + /// format. + static FourierInterpolation + from_fourier_coefficients_file(const std::string &file_name, + const unsigned int n_components); + + /// @brief Distribute the data to all parallel processes. + /// + /// Does nothing if the object has not been initialized (i.e. if @ref defined + /// returns false) on the master rank. + void distribute(const CmMod &cm_mod, const cmType &cm); + + /// @brief Return the interpolated value. + Vector value(const double time) const; + + /// @brief Return the interpolated value and time derivative. + std::pair, Vector> + value_and_derivative(const double time) const; + + /// @name Data member access. + /// @{ + + /// @brief Return whether this object has been initialized. + bool defined() const { return n_fourier_coefficients != 0; }; + + /// @brief Get the dimension of the data interpolated by this object. + unsigned int get_n_components() const { return n_components; } + + /// @brief Get the initial value of the linear trend part for one component. + const double + get_linear_trend_initial_value(const unsigned int component) const { + // @todo[michelebucelli] Add check on the index. + return linear_trend_initial_values[component]; + } + + /// @brief Get the slope of the linear trend part for one component. + const double get_linear_trend_slope(const unsigned int component) const { + // @todo[michelebucelli] Add check on the index. + return linear_trend_slopes[component]; + } + + /// @brief Get the real part of the Fourier coefficients for one component. + const double get_coefficient_real(const unsigned int component, + const unsigned int frequency) const { + // @todo[michelebucelli] Add check on the indices. + return fourier_coefficients_real(component, frequency); + } + + /// @brief Get the imaginary part of the Fourier coefficients for one + /// component. + const double get_coefficient_imaginary(const unsigned int component, + const unsigned int frequency) const { + // @todo[michelebucelli] Add check on the indices. + return fourier_coefficients_imaginary(component, frequency); + } + + /// @} + +private: + /** @brief Internal evaluation function. + * + * Uses the inverse Fourier transform to evaluate the value, and optionally + * the derivative, of the interpolated data. This function is not meant to be + * used directly, but only as a backend to @ref value and @ref + * value_and_derivative. + * + * The vectors value and derivative are assumed to be of size @ref d. This is + * not checked by this function. + * + * @throws std::runtime_error if this FourierInterpolation instance has not + * been initialized (i.e. if @ref defined returns false). + * + * @param[in] time The time at which to evaluate the interpolation. + * @param[in] evaluate_derivative Whether to also evaluate the time + * derivative of the interpolation. + * @param[out] value The interpolated value at the given time. + * @param[out] derivative The time derivative of the interpolated value at + * the given time. If evaluated_derivative is false, this will not be + * modified or accessed. + */ + void evaluate_internal(const double time, const bool evaluate_derivative, + Vector &value, + Vector &derivative) const; + + /** + * @brief Toggle whether this is a ramp function or not. + * + * See the general class documentation for details on the difference between + * ramp and periodic interpolation. + */ + bool use_ramp = false; + + /** + * @brief Number of Fourier coefficients. + */ + unsigned int n_fourier_coefficients = 0; + + /** + * @brief Number of components of the interpolated data. + */ + unsigned int n_components = 0; + + /** + * @brief Initial value for the linear trend. + * + * This is a vector with n_components entries, where each entry is the initial + * value (i.e. the value for time = ti) of the linear trend part of the + * interpolated data for that component. + */ + Vector linear_trend_initial_values; + + /** + * @brief Time derivative for the linear trend. + * + * This is a vector with n_components entries, where each entry is the slope + * of the linear trend part of the interpolated data for that component. + */ + Vector linear_trend_slopes; + + /** + * @brief Period of the interpolated data. + * + * This is disregarded if use_ramp is true. See the general class + * documentation for details on the difference between ramp and periodic + * interpolation. + */ + double period = 0.0; + + /** + * @brief Initial time. + */ + double initial_time = 0.0; + + /** + * @brief Real part of the Fourier series coefficients. + * + * This is a 2D array with n_components rows and n_fourier_coefficients + * columns. + */ + Array fourier_coefficients_real; + + /** + * @brief Imaginary part of the Fourier series coefficients. + * + * This is a 2D array with n_components rows and n_fourier_coefficients + * columns. + */ + Array fourier_coefficients_imaginary; +}; + +#endif \ No newline at end of file diff --git a/Code/Source/solver/read_files.cpp b/Code/Source/solver/read_files.cpp index 01dfc4f9e..4c70dbce3 100644 --- a/Code/Source/solver/read_files.cpp +++ b/Code/Source/solver/read_files.cpp @@ -232,11 +232,11 @@ void read_bc(Simulation* simulation, EquationParameters* eq_params, eqType& lEq, : 1; if (bc_params->temporal_values_file_path.defined()) { - lBc.gt = fcType::from_time_series_file( + lBc.gt = FourierInterpolation::from_time_series_file( bc_params->temporal_values_file_path.value(), n_dimensions, bc_params->ramp_function.value()); } else if (bc_params->fourier_coefficients_file_path.defined()) { - lBc.gt = fcType::from_fourier_coefficients_file( + lBc.gt = FourierInterpolation::from_fourier_coefficients_file( bc_params->fourier_coefficients_file_path.value(), n_dimensions); } else { throw std::runtime_error( @@ -949,11 +949,11 @@ void read_bf(ComMod& com_mod, BodyForceParameters* bf_params, bfType& lBf) lBf.bType = utils::ibset(lBf.bType, enum_int(BodyForceType::bfType_ustd)); if (bf_params->temporal_values_file_path.defined()) { - lBf.bt = fcType::from_time_series_file( + lBf.bt = FourierInterpolation::from_time_series_file( bf_params->temporal_values_file_path.value(), lBf.dof, bf_params->ramp_function.value()); } else if (bf_params->fourier_coefficients_file_path.defined()) { - lBf.bt = fcType::from_fourier_coefficients_file( + lBf.bt = FourierInterpolation::from_fourier_coefficients_file( bf_params->fourier_coefficients_file_path.value(), lBf.dof); } else { throw std::runtime_error("No temporal values or Fourier coefficients " @@ -2034,7 +2034,7 @@ void read_mat_model(Simulation* simulation, EquationParameters* eq_params, Domai utils::ibset(lDmn.stM.Tf.fType, static_cast(BoundaryConditionType::bType_ustd)); - lDmn.stM.Tf.gt = fcType::from_time_series_file( + lDmn.stM.Tf.gt = FourierInterpolation::from_time_series_file( fiber_params.temporal_values_file_path.value(), /* n_dimensions = */ 1, fiber_params.ramp_function.value()); } diff --git a/Code/Source/solver/set_bc.cpp b/Code/Source/solver/set_bc.cpp index 89e691a60..c9f92735f 100644 --- a/Code/Source/solver/set_bc.cpp +++ b/Code/Source/solver/set_bc.cpp @@ -1826,7 +1826,7 @@ void set_bc_trac_l(ComMod& com_mod, const CmMod& cm_mod, const bcType& lBc, cons } } else if (utils::btest(lBc.bType,iBC_ustd)) { - if (lBc.gt.dimension() != nsd) { + if (lBc.gt.get_n_components() != nsd) { throw std::runtime_error( "[set_bc_trac_l] Traction dof not initialized properly"); } diff --git a/tests/unitTests/math_tests/test_fft.cpp b/tests/unitTests/math_tests/test_fft.cpp index 554180577..ee28574c1 100644 --- a/tests/unitTests/math_tests/test_fft.cpp +++ b/tests/unitTests/math_tests/test_fft.cpp @@ -28,8 +28,7 @@ * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. */ -#include "fft.h" -#include "ComMod.h" +#include "fourier_interpolation.h" #include "../test_common.h" #include @@ -65,18 +64,18 @@ TEST_F(FFTTest, SinCosLinearCombination) { CreateTemporalValues(N, x_start, x_end, temporal_values); // Compute the Fourier coefficients - const fcType gt = - fcType::from_time_series(/* n_fourier_coefficients = */ 16, - temporal_values, /* is_ramp = */ false); + const FourierInterpolation gt = FourierInterpolation::from_time_series( + /* n_fourier_coefficients = */ 16, temporal_values, + /* is_ramp = */ false); // Check the slope (first Fourier coefficient) - ASSERT_NEAR(gt.get_qs(0), -0.13830, 1e-2) << "Expected slope ~-0.13830"; + ASSERT_NEAR(gt.get_linear_trend_slope(0), -0.13830, 1e-2) << "Expected slope ~-0.13830"; // Check the real and imaginary components of the first three Fourier coefficients - ASSERT_NEAR(gt.get_r(0, 0), 0.32094, 1e-2) << "Expected first real coefficient to be close to 0.32094"; - ASSERT_NEAR(gt.get_i(0, 0), 0.0, 1e-2) << "Expected first imaginary coefficient to be close to 0.0"; - ASSERT_NEAR(gt.get_r(0, 1), 0.42759, 1e-2) << "Expected second real coefficient to be close to 0.42759"; - ASSERT_NEAR(gt.get_i(0, 1), 1.25295, 1e-2) << "Expected second imaginary coefficient to be close to 1.25295"; - ASSERT_NEAR(gt.get_r(0, 2), -0.44685, 1e-2) << "Expected third real coefficient to be close to -0.44685"; - ASSERT_NEAR(gt.get_i(0, 2), -0.65403, 1e-2) << "Expected third imaginary coefficient to be close to -0.65403"; + ASSERT_NEAR(gt.get_coefficient_real(0, 0), 0.32094, 1e-2) << "Expected first real coefficient to be close to 0.32094"; + ASSERT_NEAR(gt.get_coefficient_imaginary(0, 0), 0.0, 1e-2) << "Expected first imaginary coefficient to be close to 0.0"; + ASSERT_NEAR(gt.get_coefficient_real(0, 1), 0.42759, 1e-2) << "Expected second real coefficient to be close to 0.42759"; + ASSERT_NEAR(gt.get_coefficient_imaginary(0, 1), 1.25295, 1e-2) << "Expected second imaginary coefficient to be close to 1.25295"; + ASSERT_NEAR(gt.get_coefficient_real(0, 2), -0.44685, 1e-2) << "Expected third real coefficient to be close to -0.44685"; + ASSERT_NEAR(gt.get_coefficient_imaginary(0, 2), -0.65403, 1e-2) << "Expected third imaginary coefficient to be close to -0.65403"; } \ No newline at end of file diff --git a/tests/unitTests/math_tests/test_fft.h b/tests/unitTests/math_tests/test_fft.h index 33c1c2611..f17c2082d 100644 --- a/tests/unitTests/math_tests/test_fft.h +++ b/tests/unitTests/math_tests/test_fft.h @@ -35,9 +35,6 @@ #include #include -// Forward declarations for types used in FFT testing -struct fcType; - // Test fixture class for FFT functionality class FFTTest : public ::testing::Test { protected: @@ -48,10 +45,9 @@ class FFTTest : public ::testing::Test { void TearDown() override; // Helper methods for common test operations - void CreateTemporalValues(int N, double x_start, double x_end, - std::vector>& temporal_values); - - void InitializeFourierCoefficients(fcType& gt, int d = 1, int n = 16); + void + CreateTemporalValues(int N, double x_start, double x_end, + std::vector> &temporal_values); }; #endif // FFT_TEST_H \ No newline at end of file From 03c224737948caac2ddf2509b8c9ddcfa58ab26b Mon Sep 17 00:00:00 2001 From: Michele Bucelli Date: Sun, 28 Jun 2026 14:29:52 -0500 Subject: [PATCH 09/19] Update and expand documentation of FourierInterpolation --- Code/Source/solver/fourier_interpolation.h | 257 ++++++++++++++++++--- 1 file changed, 229 insertions(+), 28 deletions(-) diff --git a/Code/Source/solver/fourier_interpolation.h b/Code/Source/solver/fourier_interpolation.h index b6199e5e6..336cd2bfd 100644 --- a/Code/Source/solver/fourier_interpolation.h +++ b/Code/Source/solver/fourier_interpolation.h @@ -12,34 +12,190 @@ #include #include -/// @brief Fourier coefficients that are used to specify time-dependent -/// functions. -/// -/// These are used for boundary conditions and for time-dependent prescribed -/// active stress. -/// -/// @todo[michelebucelli] More detailed documentation. +/** + * @brief Fourier interpolation of time dependent data. + * + * This class implements interpolation of time dependent data through either a + * Fourier series or a linear clamped ramp. Which of the two is used is + * determined by the boolean member variable @ref use_ramp. + * + * In what follows, let @f$\{t_i, \mathbf{v}_i\}_{i=0}^{N-1}@f$ be the + * interpolated time series, and @f$T = t_{N-1} - t_0@f$ be the period of the + * time series. The time series is assumed to be ordered, i.e. @f$t_i < + * t_{i+1}@f$ for all @f$i@f$. + * + * ## Using this class + * + * The FourierInterpolation class is not meant to be constructed directly, but + * rather through one of the static methods @ref from_time_series, @ref + * from_time_series_file, @ref from_fourier_coefficients, or @ref + * from_fourier_coefficients_file. + * + * If data is only read from the master rank in a parallel setting, the method + * @ref distribute can be used to broadcast the data to all ranks. + * + * Once the object has been constructed, the interpolated value at a given time + * can be obtained through the method @ref value, and the interpolated value and + * its time derivative can be obtained through the method @ref + * value_and_derivative. + * + * ## Fourier series interpolation + * + * This gives rise to a periodic interpolation, with period @f$T@f$. Let @f$M@f$ + * be the user-defined number of Fourier modes. The interpolated value at time + * @f$t@f$ is given by + * @f[ \begin{aligned} + * \tilde{\mathbf{v}}(t) &= + * \underbrace{\mathbf{q}_i + \mathbf{q}_s \tau(t)}_{\text{linear trend}} + + * \underbrace{ + * \sum_{k=0}^{M-1} \left( + * \mathbf{c}_k^{\text{real}} \cos\left(\frac{2 \pi k \tau(t)}{T}\right) + + * \mathbf{c}_k^{\text{imag}} \sin\left(\frac{2 \pi k \tau(t)}{T}\right) + * \right) + * }_{\text{Fourier series}} \\ + * \tau(t) &= (t - t_0)\;\text{mod}\;T + * \end{aligned} @f] + * The coefficients of the linear trend and Fourier series are determined from + * the input time series as follows: + * @f[ \begin{aligned} + * \mathbf{q}_i &= \mathbf{v}_0 \\ + * \mathbf{q}_s &= \frac{\mathbf{v}_{N-1} - \mathbf{v}_0}{T} \\ + * \mathbf{c}_0^{\text{real}} &= + * \frac{1}{2} \sum_{i=0}^{N-2} (\hat{t}_{i+1} - \hat{t}_i) + * (\hat{\mathbf{v}}_{i+1} + \hat{\mathbf{v}}_i) \\ + * \mathbf{c}_0^{\text{imag}} &= 0 \\ + * \mathbf{c}_k^{\text{real}} &= + * \sum_{i=0}^{N-2} \frac{\hat{\mathbf{v}}_{i+1} - \hat{\mathbf{v}}_i} + * {\hat{t}_{i+1} - \hat{t}_i} \left( + * \cos\left(\frac{2 \pi k \hat{t}_{i+1}}{T}\right) - + * \cos\left(\frac{2 \pi k \hat{t}_i}{T}\right) + * \right) \\ + * \mathbf{c}_k^{\text{imag}} &= + * -\sum_{i=0}^{N-2} \frac{\hat{\mathbf{v}}_{i+1} - \hat{\mathbf{v}}_i} + * {\hat{t}_{i+1} - \hat{t}_i} \left( + * \sin\left(\frac{2 \pi k \hat{t}_{i+1}}{T}\right) - + * \sin\left(\frac{2 \pi k \hat{t}_i}{T}\right) + * \right) + * \end{aligned} @f] + * The quantities @f$\hat{t}_i@f$ and @f$\hat{\mathbf{v}}_i@f$ are the time and + * value series after subtracting the linear trend: + * @f[ \begin{aligned} + * \hat{t}_i &= t_i - t_0 \\ + * \hat{\mathbf{v}}_i &= \mathbf{v}_i - \mathbf{q}_i - \mathbf{q}_s \hat{t}_i + * \end{aligned} @f] + * + * ## Ramp interpolation + * + * The interpolated value is equal to @f$\mathbf{v}_0@f$ until time @f$t_0@f$, + * then it follows a linear ramp from @f$\mathbf{v}_0@f$ to + * @f$\mathbf{v}_{N-1}@f$ at @f$t_{N-1}@f$, and then it remains constant at + * @f$\mathbf{v}_{N-1}@f$ for all times greater than @f$t_{N-1}@f$. Notice in + * particular that this interpolation is not periodic. + * + * The interpolated value is given by + * @f[ + * \tilde{\mathbf{v}}(t) = \begin{cases} + * \mathbf{v}_0 & t < t_0 \\ + * \mathbf{v}_0 + \frac{\mathbf{v}_{N-1} + * - \mathbf{v}_0}{t_{N-1} - t_0} (t - t_0) + * & t_0 \leq t < t_{N-1} \\ + * \mathbf{v}_{N-1} & t \geq t_{N-1} + * \end{cases} + * @f] + * + * This is equivalent to only taking the linear trend part of the Fourier series + * interpolation, and clamping the time to the interval @f$[t_0, t_{N-1}]@f$. + * + */ class FourierInterpolation { public: - /// @brief Construct an FourierInterpolation from a time series. - /// - /// @todo[michelebucelli] More detailed documentation. + /** + * @brief Default constructor. + * + * This constructor is not meant to be used directly, and is only here to + * facilitate storing objects of this type in STL containers. The object + * constructed this way is not initialized and will not be usable until it is + * assigned to a valid FourierInterpolation instance. Use the static methods + * @ref from_time_series, @ref from_time_series_file, @ref + * from_fourier_coefficients, or @ref from_fourier_coefficients_file to + * construct a valid instance. + */ + FourierInterpolation() = default; + + /** + * @brief Construct a FourierInterpolation from a time series. + * + * @param[in] n_fourier_coefficients The number @f$M@f$ of Fourier modes to + * use in the interpolation. + * @param[in] temporal_values The time series to interpolate. This is a vector + * of vectors, where each inner vector contains the time and the values of + * all components at that time. In other words, + * temporal_values[i][0] is the time @f$t_i@f$, and∏ + * temporal_values[i][j] is the value the @f$j@f$-th component of + * @f$\mathbf{v}_i@f$. + * @param[in] use_ramp Whether to use a ramp function for the interpolation. + * See the general class documentation for the precise meaning of this + * choice. + */ static FourierInterpolation from_time_series(const unsigned int n_fourier_coefficients, const std::vector> &temporal_values, const bool use_ramp); - /// @brief Read a time series from file and return the corresponding - /// FourierInterpolation instance. - /// - /// @todo[michelebucelli] More detailed documentation, especially on the file - /// format. + /** + * @brief Read a time series from file and return the corresponding instance + * of FourierInterpolation. + * + * The input file is expected to have the following format: + * ``` + * + *