Skip to content

Commit

Permalink
Improved vtk file reading using rectilinear mesh (#286)
Browse files Browse the repository at this point in the history 
* Add files via upload

* Add files via upload

* Add files via upload

* Add files via upload

* Add files via upload

* Add files via upload

* Add files via upload

* Update parameters.prm

added "set vtk file type" to select what type of vtk file to read

* This is the 64X64 2D rectilinear mesh

* Updated Body.hh for switch between RL and US mesh

* Instructions on how to create a rectilinear mesh using dream3d and convert it

* autoformatting

* fixing warnings

---------

Co-authored-by: Jason Landini <[email protected]>
  • Loading branch information
@supriyoumich @landinjm
supriyoumich and landinjm authored Nov 4, 2024
1 parent 4a5742f commit 2b66dfe
Show file tree
Hide file tree
Showing 10 changed files with 888 additions and 136 deletions.
492 changes: 492 additions & 0 deletions applications/grainGrowth_dream3d/20_grain_2D_RL.vtk
View file

Large diffs are not rendered by default.

2 changes: 1 addition & 1 deletion applications/grainGrowth_dream3d/parameters.prm
View file
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,7 @@ set Number of dimensions = 2
# Each axes spans from zero to the specified length
set Domain size X = 64
set Domain size Y = 64
set Domain size Z = 100
set Domain size Z = 1

# =================================================================================
# Set the element parameters
Expand Down
Binary file added applications/grainGrowth_dream3d/tex_files/rectilinear_grid_instructions.pdf
View file
Binary file not shown.
192 changes: 169 additions & 23 deletions include/IntegrationTools/pfield/Body.hh
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,3 @@

#ifndef Body_HH
#define Body_HH

Expand Down Expand Up @@ -39,7 +38,7 @@ namespace PRISMS
void
read_vtk(const std::string &vtkfile)
{
std::cout << "Begin reading vtk file" << std::endl;
std::cout << "Begin reading unstructured vtk file" << std::endl;

// read in vtk file here
std::ifstream infile_mesh(vtkfile.c_str());
Expand All @@ -53,7 +52,7 @@ namespace PRISMS
std::istringstream ss;
std::string str, name, type, line;
int numcomp;
unsigned long int Npoints;
unsigned long int Npoints, u, p;

while (!infile.eof())
{
Expand All @@ -69,7 +68,8 @@ namespace PRISMS
ss >> str >> Npoints;
}
}
else if (line[0] == 'S')

if (line[0] == 'S')
{
if (line.size() > 6 && line.substr(0, 7) == "SCALARS")
{
Expand All @@ -82,10 +82,11 @@ namespace PRISMS

// read data
std::cout << "begin reading data" << std::endl;
std::vector<double> data(Npoints);

std::vector<double> gid(Npoints);
for (unsigned int i = 0; i < Npoints; i++)
{
infile >> data[i];
infile >> gid[i];
// std::cout << data[i] << std::endl;
}
std::cout << " done" << std::endl;
Expand All @@ -94,14 +95,14 @@ namespace PRISMS
std::vector<std::string> var_name(DIM);
std::vector<std::string> var_description(DIM);

if (DIM >= 2)
if (DIM == 2)
{
var_name[0] = "x";
var_description[0] = "x coordinate";
var_name[1] = "y";
var_description[1] = "y coordinate";
}
if (DIM >= 3)
if (DIM > 2)
{
var_name[2] = "z";
var_description[2] = "z coordinate";
Expand All @@ -112,49 +113,191 @@ namespace PRISMS
var_name,
var_description,
mesh,
data,
gid,
0.0));
std::cout << " done" << std::endl;

gid.clear();
//
}
}
}

infile.close();
}

void
read_RL_vtk(const std::string &vtkfile)
{
std::cout << "Begin reading vtk file" << std::endl;

// read in vtk file here
std::ifstream infile_mesh(vtkfile.c_str());

// read mesh info
mesh.read_RL_vtk(infile_mesh);

std::ifstream infile(vtkfile.c_str());

// read point data
std::istringstream ss;
std::string str, name, type, line;
int numcomp;
unsigned long int N_points, Npoints_x, Npoints_y, Npoints_z, Npoints, u, p;

while (!infile.eof())
{
std::getline(infile, line);

if (line[0] == 'X')
{
if (line.size() > 12 && line.substr(0, 13) == "X_COORDINATES")
{
// read header line
// std::cout << line << "\n";
ss.clear();
ss.str(line);
ss >> str >> Npoints_x >> type;

std::cout << "Read X_COORDINATES: " << Npoints_x << std::endl;

// std::cout << " reserve OK" << std::endl;
}
}
// Alternative field descriptor used by ParaView (holds the same information as
// the "SCALAR" line above)
else if (line[0] == 'F')
if (line[0] == 'Y')
{
if (line.size() > 14 && line.substr(0, 15) == "FIELD FieldData")
if (line.size() > 12 && line.substr(0, 13) == "Y_COORDINATES")
{
// read header line
// std::cout << line << "\n";
ss.clear();
ss.str(line);
ss >> str >> numcomp;
ss >> str >> Npoints_y >> type;

// read LOOKUP_TABLE line
std::getline(infile, line);
// read points

std::cout << "Read Y_COORDINATES: " << Npoints_y << std::endl;

// std::cout << " reserve OK" << std::endl;
}
}
if (line[0] == 'Z')
{
if (line.size() > 12 && line.substr(0, 13) == "Z_COORDINATES")
{
// read header line
// std::cout << line << "\n";
ss.clear();
ss.str(line);
ss >> str >> Npoints_z >> type;

// read points

std::cout << "Read Z_COORDINATES: " << Npoints_z << std::endl;

// std::cout << " reserve OK" << std::endl;
}
}

if (line[0] == 'S')
{
if (line.size() > 6 && line.substr(0, 7) == "SCALARS")
{
ss.clear();
ss.str(line);
ss >> name >> numcomp >> Npoints >> type;
ss >> str >> name >> type >> numcomp;

// read LOOKUP_TABLE line
std::getline(infile, line);

// read data
std::cout << "begin reading data" << std::endl;
std::vector<double> data(Npoints);
for (unsigned int i = 0; i < Npoints; i++)

N_points = (Npoints_x) * (Npoints_y) * (Npoints_z);

if (DIM > 2)
{
Npoints = 8 * (Npoints_x - 1) * (Npoints_y - 1) * (Npoints_z - 1);
}
if (DIM == 2)
{
Npoints = 4 * (Npoints_x - 1) * (Npoints_y - 1);
}

std::vector<float> data(
N_points); // NB: data and gid are of different size
std::vector<double> gid(Npoints); // gid is the grainID for each node
// unsigned int gid[Npoints];

for (unsigned int i = 0; i < N_points; i++)
{
infile >> data[i];
// std::cout << data[i] << std::endl;
}
std::cout << " done" << std::endl;

std::cout << "beginning grain_id stencil" << std::endl;

u = 0;
if (DIM > 2)
{
for (unsigned int i = 0; i < (Npoints_z - 1); i++)
{
for (unsigned int j = 0; j < (Npoints_y - 1); j++)
{
for (unsigned int k = 0; k < (Npoints_x - 1); k++)
{
p = k + j * Npoints_x + i * Npoints_x * Npoints_y;

gid[u] = data[p];
gid[u + 1] = data[p + 1];
gid[u + 2] = data[p + Npoints_x];
gid[u + 3] = data[p + Npoints_x + 1];
gid[u + 4] = data[p + Npoints_x * Npoints_y];
gid[u + 5] = data[p + Npoints_x * Npoints_y + 1];
gid[u + 6] =
data[p + Npoints_x + Npoints_x * Npoints_y];
gid[u + 7] =
data[p + Npoints_x + Npoints_x * Npoints_y + 1];

u += 8;
}
}
}
}

if (DIM == 2)
{
for (unsigned int j = 0; j < (Npoints_y - 1); j++)
{
for (unsigned int k = 0; k < (Npoints_x - 1); k++)
{
p = k + j * Npoints_x;

gid[u] = data[p];
gid[u + 1] = data[p + 1];
gid[u + 2] = data[p + Npoints_x];
gid[u + 3] = data[p + Npoints_x + 1];

u += 4;
}
}
}

data.clear();
// MPI_Barrier(MPI_COMM_WORLD);

// construct field
std::vector<std::string> var_name(DIM);
std::vector<std::string> var_description(DIM);

if (DIM >= 2)
if (DIM == 2)
{
var_name[0] = "x";
var_description[0] = "x coordinate";
var_name[1] = "y";
var_description[1] = "y coordinate";
}
if (DIM >= 3)
if (DIM > 2)
{
var_name[2] = "z";
var_description[2] = "z coordinate";
Expand All @@ -165,9 +308,12 @@ namespace PRISMS
var_name,
var_description,
mesh,
data,
gid,
0.0));
std::cout << " done" << std::endl;

gid.clear();
//
}
}
}
Expand Down
6 changes: 3 additions & 3 deletions include/IntegrationTools/pfield/Coordinate.hh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,7 +13,7 @@ namespace PRISMS
template <int DIM>
class Coordinate
{
double _coord[DIM];
float _coord[DIM];

public:
int
Expand All @@ -22,13 +22,13 @@ namespace PRISMS
return DIM;
}

double &
float &
operator[](int i)
{
return _coord[i];
}

const double &
const float &
operator[](int i) const
{
return _coord[i];
Expand Down
Loading

0 comments on commit 2b66dfe

Please sign in to comment.