Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Implement read_bct() from svFSI, add some VTK funcions needed to proc… #166

Merged
merged 1 commit into from
Dec 19, 2023
Merged

Implement read_bct() from svFSI, add some VTK funcions needed to proc… #166

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
69 changes: 69 additions & 0 deletions Code/Source/svFSI/VtkData.cpp
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -629,6 +629,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 @@ -676,6 +696,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 @@ -813,6 +847,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 @@ -859,6 +907,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
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -87,7 +87,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 @@ -120,8 +122,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
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -341,10 +341,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 @@ -639,33 +636,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