Legacy VTK format is simple, easier and pdf documented if you want to save your time go for it, but if you need to run Paraview in Parallel you need to write in XML format. Though guys got such a skill probably never refers to my blog.
It's C coded although the main calculation bits is written in Fortran 90. I don't think there is particular issue but, always the rule, be aware with loops. Its assumed Staggard grid and that's why linear interpolation is embedded.
Here is the code in C
日本語 / Japanese
自作プログラム用にVTK 2D quadで結果をだす部分を書いてみた。需要ないしかなり汚いんだけどXMLサンプルとかほとんどwebになかったのでおいとく。スタガードグリッドをデータとして過程しているので2D quadに対応するために線形補間してます。あとmainで読んでる側をF90で書いてるので配列とか使うとき注意、特に問題ないと思うけど。
#include<stdio.h>
#include<stdlib.h>
#include<string.h>
void vtk_headder( FILE *wfp, int ni, int nj, int nk);
void vtk_cell_scalar_data( FILE *wfp, int ni, int nj, int nk, double *var, char name[32]);
void vtk_cell_scalar_data_int( FILE *wfp, int ni, int nj, int nk, int *var, char name[32]);
void vtk_nodal_vector_data( FILE *wfp, int ni, int nj, int nk, double *var1, double *var2, double *var3);
void vtk_coordinate( FILE *wfp, int ni, int nj, int nk, double dx, double dy, double dz);
void vtk_connectivity( FILE *wfp, int ni, int nj, int nk);
void vtk_offset( FILE *wfp, int ni, int nj, int nk);
void vtk_cell_type( FILE *wfp, int ni, int nj, int nk, int cell_type);
void vtk_footer( FILE *wfp );
void velocity_interpolation( int ni, int nj, int nk, double *var1, double *var2, double *var3);
void vtk_headder( FILE *wfp, int ni, int nj, int nk ){
//-----------------------------------------
// VTK header
//
// Varibales:
// number of nodes
// number of cells
//-----------------------------------------
fprintf(wfp, "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">\n");
fprintf(wfp, "<UnstructuredGrid>\n");
fprintf(wfp, "<Piece NumberOfPoints=\" %d \" NumberOfCells=\" %d \">\n", (ni-3)*(nj-3), (ni-4)*(nj-4));
}
void vtk_cell_scalar_data_int( FILE *wfp,int ni, int nj, int nk, int *var, char name[32]){
//-----------------------------------------
// Scalar cell data for int, material use
//
// ==Varibales==
// char name: (data's name)
// int *var: each cell
//-----------------------------------------
int i,j;
fprintf(wfp,"<CellData Scalars=\"scalar\">\n");
fprintf(wfp, "<DataArray type=\"Float32\" Name=\"%s\" format=\"ascii\">\n", name);
for( i=2; i<ni-2; i++){
for( j=2; j<nj-2; j++){
fprintf(wfp, "%d ", *(var+j+i*ni));
}
fprintf(wfp,"\n");
}
fprintf(wfp,"</DataArray>\n");
}
void vtk_cell_scalar_data( FILE *wfp,int ni, int nj, int nk, double *var, char name[32]){
//-----------------------------------------
// Scalar cell data, Pressure use
//
// ==Varibales==
// char name: (data's name)
// int *var: for each cell
//-----------------------------------------
int i,j;
fprintf(wfp, "<DataArray type=\"Float32\" Name=\"%s\" format=\"ascii\">\n", name);
for( i=2; i<ni-2; i++){
for( j=2; j<nj-2; j++){
fprintf(wfp, "%lf ", *(var+j+i*ni));
}
fprintf(wfp,"\n");
}
fprintf(wfp,"</DataArray>\n");
fprintf(wfp,"</CellData>\n");
}
void vtk_nodal_vector_data( FILE *wfp, int ni, int nj, int nk, double *var1, double *var2, double *var3){
int i,j;
fprintf(wfp, "<PointData Scalars=\"point_data\">\n");
fprintf(wfp, "<DataArray type=\"Float32\" Name=\"vars\" NumberOfComponents=\"3\" format=\"ascii\">\n");
//for( i=2; i<ni-1; i++){
// for( j=2 ; j<nj-1; j++){
for( i=1; i<ni-2; i++){
for( j=1 ; j<nj-2; j++){
fprintf( wfp, "%lf %lf %lf\n",
*(var1+j+i*ni),
*(var2+j+i*ni),
*(var3+j+i*ni)
);
}
}
fprintf(wfp,"</DataArray>\n");
fprintf(wfp,"</PointData>\n");
}
void vtk_coordinate( FILE *wfp, int ni, int nj, int nk, double dx, double dy, double dz){
int i,j;
double x;
double y;
double z;
fprintf(wfp,"<Points>\n");
fprintf(wfp,"<DataArray type=\"Float32\" Name=\"Points\" NumberOfComponents=\"3\" format=\"ascii\">\n");
for( i=2; i<ni-1; i++){
for( j=2; j<nj-1; j++){
x = (double)j*dx;
y = (double)i*dy;
z = 0.0;
fprintf(wfp, "%lf %lf %lf\n",
x,
y,
z
);
}
}
fprintf(wfp,"</DataArray>\n");
fprintf(wfp,"</Points>\n");
}
void vtk_connectivity( FILE *wfp, int ni, int nj, int nk){
int i,j;
int vtx1;
int vtx2;
int vtx3;
int vtx4;
fprintf(wfp,"<Cells>\n");
fprintf(wfp,"<DataArray type=\"Int64\" Name=\"connectivity\" format=\"ascii\">\n");
for( i=2; i<ni-2; i++){
for( j=2; j<nj-2; j++){
vtx1 = (j-2)+(ni-3)*(i-2);
vtx2 = (j-2)+1+(ni-3)*(i-2);
vtx3 = (j-2)+(ni-3)*(i-1);
vtx4 = (j-2)+1+(ni-3)*(i-1);
fprintf( wfp, "%d %d %d %d\n",
vtx3,
vtx4,
vtx2,
vtx1
);
}
}
fprintf(wfp,"</DataArray>\n");
}
void vtk_offset( FILE *wfp, int ni, int nj, int nk){
int i;
int offset;
fprintf(wfp,"<DataArray type=\"Int64\" Name=\"offsets\" format=\"ascii\" >\n");
offset = 1;
for( i=0; i<(ni-4)*(nj-4); i++){
offset = 4+4*i;
fprintf(wfp,"%d\n",offset);
}
fprintf(wfp,"</DataArray>\n");
}
void vtk_cell_type( FILE *wfp, int ni, int nj, int nk, int cell_type){
/*----------------------------------------------------------------
VTK Cell type use 9 for quad element
------------------------------------------------------------------*/
int i;
fprintf(wfp,"<DataArray type=\"UInt8\" Name=\"types\" format=\"ascii\">\n");
for( i=0; i<(ni-4)*(nj-4); i++){
fprintf(wfp,"%d\n",cell_type);
}
fprintf(wfp,"</DataArray>\n");
}
void vtk_footer( FILE *wfp){
/*----------------------------------------------------------------
VTK footter
------------------------------------------------------------------*/
fprintf(wfp,"</Cells>\n");
fprintf(wfp,"</Piece>\n");
fprintf(wfp,"</UnstructuredGrid>\n");
fprintf(wfp,"</VTKFile>\n");
}
void velocity_interpolation( int ni, int nj, int nk, double *var1, double *var2, double *var3){
int i,j;
double *tmp_var1;
double *tmp_var2;
double *tmp_var3;
tmp_var1 = (double*)malloc(sizeof(double)*ni*nj*nk);
tmp_var2 = (double*)malloc(sizeof(double)*ni*nj*nk);
tmp_var3 = (double*)malloc(sizeof(double)*ni*nj*nk);
for ( i=1; i<ni-2; i++){
for ( j=1; j<nj-2; j++){
*(tmp_var1+j+i*nj) = 0.5*( *(var1+j+i*nj) + *(var1+j+(i+1)*nj) );
*(tmp_var2+j+i*nj) = 0.5*( *(var2+j+i*nj) + *(var2+j+1+i*nj) );
//*(tmp_var3+i+j*ni);
}
}
for ( i=1; i<ni-2; i++){
for ( j=1; j<nj-2; j++){
*(var1+j+i*nj) = *(tmp_var1+j+i*nj);
*(var2+j+i*nj) = *(tmp_var2+j+i*nj);
}
}
free( tmp_var1 );
free( tmp_var2 );
free( tmp_var3 );
}
void yosuke_vtk_file_write_2d_quad_(
int *ni,
int *nj,
int *nk,
double *dx,
double *dy,
double *dz,
double *u,
double *v,
double *w,
double *p,
int *image_step,
int *fluid,
int *bound
){
int i,j,k;
char target[32];
char target_file_name[32]="result/vtk_quad_2d";
char scalar_var_name[32];
int cell_type;
/*----------------------------------------------------------------
Debug
------------------------------------------------------------------*/
//printf("%d %d %d\n", *ni, *nj, *nk);
//printf("%d\n",*image_step);
/*for( i=2; i<*ni-2; i++){
for( j=2; j<*nj-2; j++){
//printf("%d ",*(fluid+j+i**ni));
printf("%f ",*(u+j+i**ni));
}
printf("\n");
}*/
/*----------------------------------------------------------------
File name update
------------------------------------------------------------------*/
sprintf(target,"%s-%d.vtu",target_file_name,*image_step);
velocity_interpolation( *ni, *nj, *nk, u, v, w);
/*----------------------------------------------------------------
File write
------------------------------------------------------------------*/
FILE *wfp;
if( ( wfp=fopen( target, "w") ) == NULL ){
printf("ERROR: file write error can not open the target\n");
}
/*----------------------------------------------------------
VTK file for 2d quad element
This file format is supported by Paraview
-----------------------------------------------------------*/
vtk_headder(&*wfp, *ni,*nj,*nk);
/*----------------------------------------------------------
Material (Cell Data)
-----------------------------------------------------------*/
strcpy( scalar_var_name, "material" );
vtk_cell_scalar_data_int(&*wfp, *ni, *nj, *nk, fluid, scalar_var_name);
/*----------------------------------------------------------
Pressure (Cell Data)
-----------------------------------------------------------*/
strcpy( scalar_var_name, "pressure" );
vtk_cell_scalar_data(&*wfp, *ni, *nj, *nk, p, scalar_var_name);
/*----------------------------------------------------------
Velocity (Nodal Data)
-----------------------------------------------------------*/
vtk_nodal_vector_data(&*wfp, *ni, *nj, *nk, u, v, w);
/*----------------------------------------------------------
Nodal Coordinate
-----------------------------------------------------------*/
vtk_coordinate(&*wfp, *ni, *nj, *nk, *dx, *dy, *dz);
/*----------------------------------------------------------
Connectivity
-----------------------------------------------------------*/
vtk_connectivity(&*wfp, *ni, *nj, *nk);
/*----------------------------------------------------------
Offset (4 each)
-----------------------------------------------------------*/
vtk_offset(&*wfp, *ni, *nj, *nk);
/*----------------------------------------------------------
Cell Type ( Use 9 for quad element )
-----------------------------------------------------------*/
cell_type=9;
vtk_cell_type(&*wfp, *ni, *nj, *nk, cell_type);
/*----------------------------------------------------------
VTK footter
-----------------------------------------------------------*/
vtk_footer(&*wfp);
fclose(wfp);
}
0 件のコメント:
コメントを投稿