1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72
|
#include <iostream>
#include <fstream>
#include <vector>
using namespace std;
//======================================================================
void writeVTK( ostream &out, // output file stream
const vector<double> &x, // x coordinates
const vector<double> &y, // y coordinates
const vector<double> &z, // z coordinates
const vector<double> &p, // pressure
const vector<double> &u, // velocity-x
const vector<double> &v, // velocity-y
const vector<double> &w, // velocity-z
int NX, int NY, int NZ ) // dimensions x, y, z
{
int N = NX * NY * NZ;
out << "# vtk DataFile Version 3.0\n";
out << "Test file\n";
out << "ASCII\n";
out << "DATASET STRUCTURED_GRID\n";
out << "DIMENSIONS " << ' ' << NX << ' ' << NY << ' ' << NZ << '\n';
out << "POINTS " << N << " double\n";
for ( int i = 0; i < N; i++ ) out << x[i] << " " << y[i] << " " << z[i] << '\n';
out << "POINT_DATA " << N << '\n';
out << "SCALARS pressure double 1\n";
out << "LOOKUP_TABLE default\n";
for ( int i = 0; i < N; i++ ) out << p[i] << '\n';
out << "VECTORS velocity double\n";
for ( int i = 0; i < N; i++ ) out << u[i] << " " << v[i] << " " << w[i] << '\n';
}
//======================================================================
int main()
{
int NX = 5, NY = 10, NZ = 3;
int N = NX * NY * NZ;
vector<double> x(N), y(N), z(N);
vector<double> u(N), v(N), w(N), p(N);
// Some data (plane Poiseuille flow)
double pressureGradient = 2.5;
double depth = 2.0;
double U0 = 1.0;
int n = 0;
for ( int k = 0; k < NZ; k++ )
{
for ( int j = 0; j < NY; j++ )
{
for ( int i = 0; i < NX; i++ )
{
x[n] = i;
y[n] = depth * j / ( NY - 1 );
z[n] = k;
u[n] = U0 * ( y[n] / depth ) * ( 1.0 - y[n] / depth );
v[n] = w[n] = 0.0;
p[n] = NX - 1 - pressureGradient * i;
n++;
}
}
}
ofstream out( "temp.vtk" );
writeVTK( out, x, y, z, p, u, v, w, NX, NY, NZ );
}
|