diff --git a/Code/Source/svFSI/main.cpp b/Code/Source/svFSI/main.cpp index 09045754..6c81a569 100644 --- a/Code/Source/svFSI/main.cpp +++ b/Code/Source/svFSI/main.cpp @@ -89,6 +89,91 @@ void read_files(Simulation* simulation, const std::string& file_name) /// /// Reproduces the outer and inner loops in Fortan MAIN.f. // + +/// @brief Iterate the precomputed state-variables in time using linear interpolation to the current time step size +/// +void iterate_precomputed_time(Simulation* simulation) { + using namespace consts; + + auto& com_mod = simulation->com_mod; + auto& cm_mod = simulation->cm_mod; + auto& cm = com_mod.cm; + auto& cep_mod = simulation->get_cep_mod(); + + int nTS = com_mod.nTS; + int stopTS = nTS; + int tDof = com_mod.tDof; + int tnNo = com_mod.tnNo; + int nFacesLS = com_mod.nFacesLS; + int nsd = com_mod.nsd; + + auto& Ad = com_mod.Ad; // Time derivative of displacement + auto& Rd = com_mod.Rd; // Residual of the displacement equation + auto& Kd = com_mod.Kd; // LHS matrix for displacement equation + + auto& Ao = com_mod.Ao; // Old time derivative of variables (acceleration) + auto& Yo = com_mod.Yo; // Old variables (velocity) + auto& Do = com_mod.Do; // Old integrated variables (dissplacement) + + auto& An = com_mod.An; // New time derivative of variables + auto& Yn = com_mod.Yn; // New variables (velocity) + auto& Dn = com_mod.Dn; // New integrated variables + + int& cTS = com_mod.cTS; + int& nITs = com_mod.nITs; + double& dt = com_mod.dt; + + if (com_mod.usePrecomp) { +#ifdef debug_iterate_solution + dmsg << "Use precomputed values ..." << std::endl; +#endif + // This loop is used to interpolate between known time values of the precomputed + // state-variable solution + for (int l = 0; l < com_mod.nMsh; l++) { + auto lM = com_mod.msh[l]; + if (lM.Ys.nslices() > 1) { + // If there is only one temporal slice, then the solution is assumed constant + // in time and no interpolation is performed + // If there are multiple temporal slices, then the solution is linearly interpolated + // between the known time values and the current time. + double precompDt = com_mod.precompDt; + double preTT = precompDt * (lM.Ys.nslices() - 1); + double cT = cTS * dt; + double rT = std::fmod(cT, preTT); + int n1, n2; + double alpha; + if (precompDt == dt) { + alpha = 0.0; + if (cTS < lM.Ys.nslices()) { + n1 = cTS - 1; + } else { + n1 = cTS % lM.Ys.nslices() - 1; + } + } else { + n1 = static_cast(rT / precompDt) - 1; + alpha = std::fmod(rT, precompDt); + } + n2 = n1 + 1; + for (int i = 0; i < tnNo; i++) { + for (int j = 0; j < nsd; j++) { + if (alpha == 0.0) { + Yn(j, i) = lM.Ys(j, i, n2); + } else { + Yn(j, i) = (1.0 - alpha) * lM.Ys(j, i, n1) + alpha * lM.Ys(j, i, n2); + } + } + } + } else { + for (int i = 0; i < tnNo; i++) { + for (int j = 0; j < nsd; j++) { + Yn(j, i) = lM.Ys(j, i, 0); + } + } + } + } + } +} + void iterate_solution(Simulation* simulation) { using namespace consts; @@ -244,52 +329,8 @@ void iterate_solution(Simulation* simulation) #endif set_bc::set_bc_dir(com_mod, An, Yn, Dn); - - if (com_mod.usePrecomp) { - #ifdef debug_iterate_solution - dmsg << "Use precomputed values ..." << std::endl; - #endif - // This loop is used to interpolate between known time values of the precomputed - // state-variable solution - for (int l = 0; l < com_mod.nMsh; l++) { - auto lM = com_mod.msh[l]; - if (lM.Ys.nslices() == 1) { - // If there is only one temporal slice, then the solution is assumed constant - // in time and no interpolation is performed - continue; - } else { - // If there are multiple temporal slices, then the solution is linearly interpolated - // between the known time values and the current time. - double precompDt = com_mod.precompDt; - double preTT = precompDt * (lM.Ys.nslices() - 1); - double cT = cTS * dt; - double rT = std::fmod(cT, preTT); - int n1, n2; - double alpha; - if (precompDt == dt) { - alpha = 0.0; - if (cTS < lM.Ys.nslices()) { - n1 = cTS - 1; - } else { - n1 = cTS % lM.Ys.nslices() - 1; - } - } else { - n1 = static_cast(rT / precompDt) - 1; - alpha = std::fmod(rT, precompDt); - } - n2 = n1 + 1; - for (int i = 0; i < tnNo; i++) { - for (int j = 0; j < nsd; j++) { - if (alpha == 0.0) { - Yn(j, i) = lM.Ys(j, i, n2); - } else { - Yn(j, i) = (1.0 - alpha) * lM.Ys(j, i, n1) + alpha * lM.Ys(j, i, n2); - } - } - } - } - } - } + + iterate_precomputed_time(simulation); // Inner loop for Newton iteration // diff --git a/Code/Source/svFSI/vtk_xml.cpp b/Code/Source/svFSI/vtk_xml.cpp index cfaf0f76..2d5804e0 100644 --- a/Code/Source/svFSI/vtk_xml.cpp +++ b/Code/Source/svFSI/vtk_xml.cpp @@ -597,17 +597,12 @@ void read_vtu(const std::string& file_name, mshType& mesh) //---------- // read_precomputed_solution_vtu //---------- -// Read a mesh from a SimVascular .vtu or .vtp file. -// -// Mesh variables set -// mesh.gnNo - number of nodes -// mesh.gnEl - number of elements -// mesh.eNoN - number of noders per element -// mesh.x - node coordinates -// mesh.gIEN - element connectivity -// -// -// + +/// @brief Read a mesh file with state-variable solution fields from a .vtu or .vtp file. +/// +/// Mesh variables set +/// mesh.Ys = precomputed state-variable solutions (e.g. velocity, pressure, etc.) + void read_precomputed_solution_vtu(const std::string& file_name, const std::string& field_name, mshType& mesh) { using namespace vtk_xml_parser; diff --git a/Code/Source/svFSI/vtk_xml_parser.cpp b/Code/Source/svFSI/vtk_xml_parser.cpp index 4b84930d..3abd5eb6 100644 --- a/Code/Source/svFSI/vtk_xml_parser.cpp +++ b/Code/Source/svFSI/vtk_xml_parser.cpp @@ -755,20 +755,14 @@ void load_vtu(const std::string& file_name, mshType& mesh) } -/// $brief Read a time series field from a VTK .vtu file. +/// @brief Read a time series field from a VTK .vtu file. /// /// Mesh variables set -/// mesh.gnNo - number of nodes -/// mesh.x - node coordinates -/// mesh.gN - node IDs -/// mesh.gnEl - number of elements -/// mesh.eNoN - number of noders per element -/// mesh.gIEN - element connectivity (num_nodes_per_elem, num_elems) /// mesh.Ys - time series field data (num_components, num_nodes, num_time_steps) /// // -void load_time_varying_field_vtu(const std::string file_name, const std::string field_name, mshType& mesh) { - +void load_time_varying_field_vtu(const std::string file_name, const std::string field_name, mshType& mesh) +{ #define n_debug_load_vtu #ifdef debug_load_vtu std::cout << "[load_vtu] " << std::endl; @@ -785,49 +779,50 @@ void load_time_varying_field_vtu(const std::string file_name, const std::string std::vector> array_names; if (num_nodes == 0) { - throw std::runtime_error("Failed reading the VTK file '" + file_name + "'."); + throw std::runtime_error("Failed reading the VTK file '" + file_name + "'."); } // Store all array names for (int i = 0; i < vtk_ugrid->GetPointData()->GetNumberOfArrays(); i++) { - std::string array_name = vtk_ugrid->GetPointData()->GetArrayName(i); - size_t pos = array_name.find(field_name.c_str()); - if (pos != std::string::npos) { - auto not_digit = [](char c) { return !std::isdigit(c); }; - auto it = std::find_if(array_name.rbegin(), array_name.rend(), not_digit); - std::string time_step = std::string(it.base(), array_name.end()); - array_count++; - if (!time_step.empty()) { - array_names.push_back({array_name, std::stoi(time_step)}); - } else { - array_names.push_back({array_name, 0}); - } + std::string array_name = vtk_ugrid->GetPointData()->GetArrayName(i); + size_t pos = array_name.find(field_name.c_str()); + if (pos != std::string::npos) { + auto not_digit = [](char c) { return !std::isdigit(c); }; + auto it = std::find_if(array_name.rbegin(), array_name.rend(), not_digit); + std::string time_step = std::string(it.base(), array_name.end()); + array_count++; + if (!time_step.empty()) { + array_names.push_back({array_name, std::stoi(time_step)}); + } else { + array_names.push_back({array_name, 0}); } + } } - + // Check if there are any fields present in the VTK file if (array_count == 0) { - throw std::runtime_error("No '" + field_name + "' data found in the VTK file '" + file_name + "'."); + throw std::runtime_error("No '" + field_name + "' data found in the VTK file '" + file_name + "'."); } - // Order all array names + // Order all array names by time step std::sort(array_names.begin(), array_names.end(), [](const std::pair& a, const std::pair& b) { - return a.second < b.second; + return a.second < b.second; }); + // Get the expected number of state-variable components int num_components = vtk_ugrid->GetPointData()->GetArray(array_names[0].first.c_str())->GetNumberOfComponents(); mesh.Ys.resize(num_components, num_nodes, array_count); for (int i = 0; i < array_count; i++) { - auto array = vtk_ugrid->GetPointData()->GetArray(array_names[i].first.c_str()); - if (array == nullptr) { - throw std::runtime_error("No '" + array_names[i].first + "' data found in the VTK file '" + file_name + "'."); - } - if (array->GetNumberOfComponents() != num_components) { - throw std::runtime_error("The number of components in the field '" + array_names[i].first + "' is not equal to the number of components in the first field."); - } - for (int j = 0; j < num_nodes; j++) { - for (int k = 0; k < num_components; k++) { - mesh.Ys(k, j, i) = array->GetComponent(j, k); - } + auto array = vtk_ugrid->GetPointData()->GetArray(array_names[i].first.c_str()); + if (array == nullptr) { + throw std::runtime_error("No '" + array_names[i].first + "' data found in the VTK file '" + file_name + "'."); + } + if (array->GetNumberOfComponents() != num_components) { + throw std::runtime_error("The number of components in the field '" + array_names[i].first + "' is not equal to the number of components in the first field."); + } + for (int j = 0; j < num_nodes; j++) { + for (int k = 0; k < num_components; k++) { + mesh.Ys(k, j, i) = array->GetComponent(j, k); } + } } } } // namespace vtk_utils