Skip to content

Commit

Permalink
adding Doxygen compatible comments for new functions in vtk_xml.cpp a…
Browse files Browse the repository at this point in the history 
…nd vtk_xml_parser.cpp, encapsulated precomputed linear interpolation time advancement in a function to remove from main.cpp (SimVascular#202)
  • Loading branch information
@zasexton
zasexton committed Mar 29, 2024
1 parent 27fe339 commit f819abd
Show file tree
Hide file tree
Showing 3 changed files with 125 additions and 94 deletions.
133 changes: 87 additions & 46 deletions Code/Source/svFSI/main.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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<int>(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;
Expand Down Expand Up @@ -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<int>(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
//
Expand Down
17 changes: 6 additions & 11 deletions Code/Source/svFSI/vtk_xml.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand Down
69 changes: 32 additions & 37 deletions Code/Source/svFSI/vtk_xml_parser.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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;
Expand All @@ -785,49 +779,50 @@ void load_time_varying_field_vtu(const std::string file_name, const std::string
std::vector<std::pair<std::string, int>> 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<std::string, int>& a, const std::pair<std::string, int>& 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
Expand Down

0 comments on commit f819abd

Please sign in to comment.