Skip to content

Commit

Permalink
Implement read_bct() from svFSI, add some VTK funcions needed to proc…
Browse files Browse the repository at this point in the history 
…ess that data in the bc.vtp file. (#166)
  • Loading branch information
@ktbolt
ktbolt authored Dec 19, 2023
1 parent dbdb17c commit c4a1cc6
Show file tree
Hide file tree
Showing 3 changed files with 208 additions and 13 deletions.
69 changes: 69 additions & 0 deletions Code/Source/svFSI/VtkData.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -628,6 +628,26 @@ void VtkVtpData::copy_point_data(const std::string& data_name, Vector<double>& m
}
}

void VtkVtpData::copy_point_data(const std::string& data_name, Vector<int>& mesh_data)
{
auto vtk_data = vtkIntArray::SafeDownCast(impl->vtk_polydata->GetPointData()->GetArray(data_name.c_str()));
if (vtk_data == nullptr) {
return;
}

int num_data = vtk_data->GetNumberOfTuples();
if (num_data == 0) {
return;
}

int num_comp = vtk_data->GetNumberOfComponents();

// Set the data.
for (int i = 0; i < num_data; i++) {
mesh_data(i) = vtk_data->GetValue(i);
}
}

/// @brief Copy points into the given array.
//
void VtkVtpData::copy_points(Array<double>& points)
Expand Down Expand Up @@ -675,6 +695,20 @@ Array<double> VtkVtpData::get_point_data(const std::string& data_name)
return data;
}

/// @brief Get a list of point data names.
std::vector<std::string> VtkVtpData::get_point_data_names()
{
std::vector<std::string> data_names;
int num_arrays = impl->vtk_polydata->GetPointData()->GetNumberOfArrays();

for (int i = 0; i < num_arrays; i++) {
auto array_name = impl->vtk_polydata->GetPointData()->GetArrayName(i);
data_names.push_back(array_name);
}

return data_names;
}

/// @brief Get an array of point data from an unstructured grid.
//
Array<double> VtkVtpData::get_points()
Expand Down Expand Up @@ -812,6 +846,20 @@ Array<int> VtkVtuData::get_connectivity()
return conn;
}

/// @brief Get a list of point data names.
std::vector<std::string> VtkVtuData::get_point_data_names()
{
std::vector<std::string> data_names;
int num_arrays = impl->vtk_ugrid->GetPointData()->GetNumberOfArrays();

for (int i = 0; i < num_arrays; i++) {
auto array_name = impl->vtk_ugrid->GetPointData()->GetArrayName(i);
data_names.push_back(array_name);
}

return data_names;
}

/// @brief Copy an array of point data from an unstructured grid into the given Array.
//
void VtkVtuData::copy_point_data(const std::string& data_name, Array<double>& mesh_data)
Expand Down Expand Up @@ -858,6 +906,27 @@ void VtkVtuData::copy_point_data(const std::string& data_name, Vector<double>& m
}
}

void VtkVtuData::copy_point_data(const std::string& data_name, Vector<int>& mesh_data)
{
auto vtk_data = vtkIntArray::SafeDownCast(impl->vtk_ugrid->GetPointData()->GetArray(data_name.c_str()));
if (vtk_data == nullptr) {
return;
}

int num_data = vtk_data->GetNumberOfTuples();
if (num_data == 0) {
return;
}

int num_comp = vtk_data->GetNumberOfComponents();

// Set the data.
for (int i = 0; i < num_data; i++) {
mesh_data[i] = vtk_data->GetValue(i);
}
}


bool VtkVtuData::has_point_data(const std::string& data_name)
{
int num_arrays = impl->vtk_ugrid->GetPointData()->GetNumberOfArrays();
Expand Down
4 changes: 4 additions & 0 deletions Code/Source/svFSI/VtkData.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -86,7 +86,9 @@ class VtkVtpData : public VtkData {
void copy_points(Array<double>& points);
void copy_point_data(const std::string& data_name, Array<double>& mesh_data);
void copy_point_data(const std::string& data_name, Vector<double>& mesh_data);
void copy_point_data(const std::string& data_name, Vector<int>& mesh_data);
Array<double> get_point_data(const std::string& data_name);
std::vector<std::string> get_point_data_names();
bool has_point_data(const std::string& data_name);
virtual void set_connectivity(const int nsd, const Array<int>& conn, const int pid = 0);

Expand Down Expand Up @@ -119,8 +121,10 @@ class VtkVtuData : public VtkData {

void copy_point_data(const std::string& data_name, Array<double>& mesh_data);
void copy_point_data(const std::string& data_name, Vector<double>& mesh_data);
void copy_point_data(const std::string& data_name, Vector<int>& mesh_data);

Array<double> get_point_data(const std::string& data_name);
std::vector<std::string> get_point_data_names();
virtual Array<double> get_points();
bool has_point_data(const std::string& data_name);
virtual void set_connectivity(const int nsd, const Array<int>& conn, const int pid = 0);
Expand Down
148 changes: 135 additions & 13 deletions Code/Source/svFSI/read_files.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -340,10 +340,7 @@ void read_bc(Simulation* simulation, EquationParameters* eq_params, eqType& lEq,
// [NOTE] This is not implemented.
if (bc_params->bct_file_path.defined()) {
auto file_name = bc_params->bct_file_path.value();
throw std::runtime_error("[read_bc] read_bct is not imlemented.");
read_bct(com_mod, lBc.gm, com_mod.msh[iM].fa[iFa], file_name);
//ALLOCATE(lBc%gm)
//CALL READBCT(lBc%gm, msh(iM)%fa(iFa), ftmp%fname)

} else {
if (bc_params->temporal_and_spatial_values_file_path.defined()) {
Expand Down Expand Up @@ -638,33 +635,158 @@ void read_bc(Simulation* simulation, EquationParameters* eq_params, eqType& lEq,
//
// Reproduces 'SUBROUTINE READBCT(lMB, lFa, fName)' in READFILES.f.
//
// [TODO:DaveP] this is not implemented.
//
void read_bct(ComMod& com_mod, MBType& lMB, faceType& lFa, const std::string& fName)
{
// The name of the arrays to search for in the bct.vtp file.
static std::string shdr("velocity_");

static std::string file_desc("general velocity data file '" + fName + "'");

if (FILE *file = fopen(fName.c_str(), "r")) {
fclose(file);
} else {
throw std::runtime_error("The VTK VTP bct data file '" + fName + "' can't be read.");
throw std::runtime_error("The " + file_desc + " can't be read.");
}

// Read the vtp file.
//
VtkVtpData vtp_data(fName);
int num_points = vtp_data.num_points();
if (num_points == 0) {
throw std::runtime_error("The VTK VTP bct data file '" + fName + "' does not contain any points.");
int nNo = vtp_data.num_points();
if (nNo == 0) {
throw std::runtime_error("The " + file_desc + " does not contain any points.");
}

int num_elems = vtp_data.num_elems();
if (num_elems == 0) {
throw std::runtime_error("The VTK VTP bct data file '" + fName + "' does not contain any elements.");
throw std::runtime_error("The " + file_desc + " does not contain any elements.");
}

if (nNo != lFa.nNo) {
throw std::runtime_error("The number of points (" + std::to_string(nNo) + ") in the " + file_desc +
" does not match the number of points (" + std::to_string(nNo) + " for the face '" + lFa.name + "'.");
}

// Get all the point data starting with "velocity_"
//
auto namesL = vtp_data.get_point_data_names();
int n = namesL.size();

int ntime = 0;
int nj = 1;

for (int i = 0; i < n; i++) {
auto stmp = namesL[i];

for (int j = 0; j < 9; j++) {
if (shdr[j] != stmp[j]) {
break;
}
nj += 1;
}

if (nj < 9) {
namesL[i] = "";
continue;
}
ntime = ntime + 1;
}

// Initialize lMB data structure
//
const int nsd = com_mod.nsd;
lMB.dof = nsd;
lMB.nTP = ntime;
int iM = lFa.iM;
lMB.t.resize(ntime);
lMB.d.resize(nsd,nNo,ntime);
Vector<int> ptr(com_mod.msh[iM].gnNo);

ntime = 0;

for (int i = 0; i < n; i++) {
auto stmp = namesL[i];
if (stmp.size() == 0) {
continue;
}

stmp = stmp.substr(9);
double t = std::stod(stmp);
lMB.t(ntime) = t;

if (ntime == 0) {
if (!utils::is_zero(t)) {
throw std::runtime_error("The first time step in the " + file_desc + " should be zero.");
}
} else {
t = t - lMB.t(ntime-1);
if (utils::is_zero(t) || t < 0.0) {
throw std::runtime_error("There is a non-increaing series of times in the " + file_desc + ".");
}
}

ntime = ntime + 1;
}

if (num_points != lFa.nNo) {
throw std::runtime_error("The number of points (" + std::to_string(num_points) + " in the VTK VTP bct data file '" + fName +
"' does not match the number of points (" + std::to_string(num_points) + " for the face '" + lFa.name + "'.");
lMB.period = lMB.t(ntime-1);
auto& msh = com_mod.msh[iM];

// Prepare pointer array
//
for (int a = 0; a < lFa.nNo; a++) {
int Ac = lFa.gN(a);
Ac = msh.lN(Ac);

if (Ac == -1) {
throw std::runtime_error("Incorrect global node number detected in the " + file_desc +
". Mesh: '" + msh.name + "' Face: '" + lFa.name + "' Node: " + std::to_string(a+1) +
" gN: " + std::to_string(lFa.gN(a)));
}

ptr(Ac) = a;
}

// Get GlobalNodeID from vtp file and make sure it is consistent
// with mesh structure
//
Vector<int> gN(nNo);
vtp_data.copy_point_data("GlobalNodeID", gN);

for (int a = 0; a < nNo; a++) {
int Ac = gN(a) - 1;
if (Ac > msh.gnNo-1 || Ac < 0) {
throw std::runtime_error("Global node IDs in the '" + file_desc + " are not consistent with the volume mesh '" + msh.name + ".");
}
Ac = ptr(Ac);
if (Ac == -1) {
throw std::runtime_error("Global node IDs in the '" + file_desc + " are not consistent with the volume mesh '" + msh.name + ".");
}
}

// Load spatial data for each time point from vtp file
//
Array<double> tmpR(nsd,nNo);
ntime = 0;

for (int i = 0; i < n; i++) {
auto stmp = namesL[i];
if (stmp.size() == 0) {
continue;
}

tmpR = 0.0;
vtp_data.copy_point_data(stmp, tmpR);

for (int a = 0; a < nNo; a++) {
int Ac = gN(a) - 1;
Ac = ptr(Ac);
for (int j = 0; j < nsd; j++) {
lMB.d(j,Ac,ntime) = tmpR(j,a);
}
}

ntime = ntime + 1;
}

}

//---------
Expand Down

0 comments on commit c4a1cc6

Please sign in to comment.