[cig-commits] r18742 - in mc/3D/CitcomS/trunk: CitcomS/Components lib module
tan2 at geodynamics.org
tan2 at geodynamics.org
Tue Jul 12 09:27:58 PDT 2011
Author: tan2
Date: 2011-07-12 09:27:58 -0700 (Tue, 12 Jul 2011)
New Revision: 18742
Modified:
mc/3D/CitcomS/trunk/CitcomS/Components/Output.py
mc/3D/CitcomS/trunk/lib/Instructions.c
mc/3D/CitcomS/trunk/lib/Output.c
mc/3D/CitcomS/trunk/lib/Output_vtk.c
mc/3D/CitcomS/trunk/lib/global_defs.h
mc/3D/CitcomS/trunk/module/setProperties.c
Log:
enhancing VTK output, mostly contributed by Rene Gassmoeller at GFZ.
- A new input parameter 'vtk_format', which can be either 'ascii' (default) or
'binary'. When 'binary', gzip and base64 encoding are used.
- The node ordering is changed to CitcomS native ordering. The consequence is the
axis are rotated as z in CitcomS is mapped to X in vtk, and x to Y, y to Z.
- More data can be outputted in vtk format
- Using .pvts file for regional model; .vtm and .visit file for global model
(.vtm format provided by Tobias Hoeink at Rice, .visit format provided by Kat
at visitusers.org/forum)
Modified: mc/3D/CitcomS/trunk/CitcomS/Components/Output.py
===================================================================
--- mc/3D/CitcomS/trunk/CitcomS/Components/Output.py 2011-07-12 15:18:21 UTC (rev 18741)
+++ mc/3D/CitcomS/trunk/CitcomS/Components/Output.py 2011-07-12 16:27:58 UTC (rev 18742)
@@ -56,6 +56,9 @@
output_optional = inv.str("output_optional",
default="surf,botm,tracer")
+ vtk_format = inv.str("vtk_format", default="ascii",
+ validator=inv.choice(["binary", "ascii"]))
+
# experimental vtk output
gzdir_vtkio = inv.int("gzdir_vtkio", default=0)
# remove net rotation
Modified: mc/3D/CitcomS/trunk/lib/Instructions.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Instructions.c 2011-07-12 15:18:21 UTC (rev 18741)
+++ mc/3D/CitcomS/trunk/lib/Instructions.c 2011-07-12 16:27:58 UTC (rev 18742)
@@ -809,6 +809,16 @@
*/
}
+ if (strcmp(E->output.vtk_format, "binary") == 0) {
+#ifndef USE_GZDIR
+ /* zlib is required for vtk binary output */
+ if(E->parallel.me == 0) {
+ fputs("VTK binary output requires zlib, but couldn't find it at 'configure' time. Please either use VTK ascii output or re-run 'configure' with suitable flags.\n", stderr);
+ }
+ parallel_process_termination();
+#endif
+ }
+
return;
}
Modified: mc/3D/CitcomS/trunk/lib/Output.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Output.c 2011-07-12 15:18:21 UTC (rev 18741)
+++ mc/3D/CitcomS/trunk/lib/Output.c 2011-07-12 16:27:58 UTC (rev 18742)
@@ -93,6 +93,17 @@
//fprintf(stderr,"gzdir: vtkio: %i save basis vectors: %i\n",
// E->output.gzdir.vtk_io,E->output.gzdir.vtk_base_save);
}
+
+ if(strcmp(E->output.format, "vtk") == 0) {
+ input_string("vtk_format", E->output.vtk_format, "ascii",m);
+ if (strcmp(E->output.vtk_format, "binary") != 0 &&
+ strcmp(E->output.vtk_format, "ascii") != 0) {
+ if(E->parallel.me == 0) {
+ fprintf(stderr, "Unknown vtk_format: %s\n", E->output.vtk_format);
+ }
+ parallel_process_termination();
+ }
+ }
}
Modified: mc/3D/CitcomS/trunk/lib/Output_vtk.c
===================================================================
--- mc/3D/CitcomS/trunk/lib/Output_vtk.c 2011-07-12 15:18:21 UTC (rev 18741)
+++ mc/3D/CitcomS/trunk/lib/Output_vtk.c 2011-07-12 16:27:58 UTC (rev 18742)
@@ -25,31 +25,41 @@
*
*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
*/
+
/* Routine to process the output of the finite element cycles
- and to turn them into a coherent suite files */
+ and to turn them into a coherent suite of files */
-#include <stdlib.h>
#include <math.h>
#include "element_definitions.h"
#include "global_defs.h"
#include "output.h"
+#include "string.h"
+#ifdef USE_GZDIR
+#include "zlib.h"
+#endif
+
+#define CHUNK 16384
+
+static void write_binary_array(int nn, float* array, FILE * f);
+static void write_ascii_array(int nn, int perLine, float *array, FILE *fp);
+
static void vts_file_header(struct All_variables *E, FILE *fp)
{
const char format[] =
"<?xml version=\"1.0\"?>\n"
- "<VTKFile type=\"StructuredGrid\" version=\"0.1\">\n"
+ "<VTKFile type=\"StructuredGrid\" version=\"0.1\" compressor=\"vtkZLibDataCompressor\" byte_order=\"LittleEndian\">\n"
" <StructuredGrid WholeExtent=\"%s\">\n"
" <Piece Extent=\"%s\">\n";
char extent[64], header[1024];
snprintf(extent, 64, "%d %d %d %d %d %d",
+ E->lmesh.ezs, E->lmesh.ezs + E->lmesh.elz,
E->lmesh.exs, E->lmesh.exs + E->lmesh.elx,
- E->lmesh.eys, E->lmesh.eys + E->lmesh.ely,
- E->lmesh.ezs, E->lmesh.ezs + E->lmesh.elz);
+ E->lmesh.eys, E->lmesh.eys + E->lmesh.ely);
snprintf(header, 1024, format, extent, extent);
@@ -75,8 +85,6 @@
static void vtk_point_data_header(struct All_variables *E, FILE *fp)
{
fputs(" <PointData Scalars=\"temperature\" Vectors=\"velocity\">\n", fp);
-
-
return;
}
@@ -104,25 +112,22 @@
static void vtk_output_temp(struct All_variables *E, FILE *fp)
{
- int i, j;
- int m,k,l,nox,noy,noz,noxnoz;
- noy=E->lmesh.noy;
- nox=E->lmesh.nox;
- noz=E->lmesh.noz;
- noxnoz = nox * noz;
+ int i;
+ int nodes = E->sphere.caps_per_proc*E->lmesh.nno;
+ float* floattemp = malloc(nodes*sizeof(float));
- fputs(" <DataArray type=\"Float32\" Name=\"temperature\" format=\"ascii\">\n", fp);
+ fprintf(fp, " <DataArray type=\"Float32\" Name=\"temperature\" format=\"%s\">\n", E->output.vtk_format);
- for(j=1; j<=E->sphere.caps_per_proc; j++) {
- for(m=1;m <= noz;m++)
- for(k=1;k <= noy;k++)
- for(l=1;l <= nox;l++){
- i = m + (l-1) * noz + (k-1) * noxnoz;
- fprintf(fp, "%.6e\n", E->T[j][i]);
- }
+ for(i=0;i <= nodes;i++)
+ floattemp[i] = (float) *(E->T[1]+i+1);
+
+ if (strcmp(E->output.vtk_format,"binary") == 0) {
+ write_binary_array(nodes,floattemp,fp);
+ } else {
+ write_ascii_array(nodes,1,floattemp,fp);
}
-
fputs(" </DataArray>\n", fp);
+ free(floattemp);
return;
}
@@ -130,65 +135,54 @@
static void vtk_output_velo(struct All_variables *E, FILE *fp)
{
int i, j;
+ int nodes=E->sphere.caps_per_proc*E->lmesh.nno;
double sint, sinf, cost, cosf;
float *V[4];
const int lev = E->mesh.levmax;
- int m,k,l,nox,noy,noz,noxnoz;
- noy=E->lmesh.noy;
- nox=E->lmesh.nox;
- noz=E->lmesh.noz;
- noxnoz = nox * noz;
+ float* floatvel = malloc(nodes*3*sizeof(float));
- fputs(" <DataArray type=\"Float32\" Name=\"velocity\" NumberOfComponents=\"3\" format=\"ascii\">\n", fp);
+ fprintf(fp, " <DataArray type=\"Float32\" Name=\"velocity\" NumberOfComponents=\"3\" format=\"%s\">\n", E->output.vtk_format);
for(j=1; j<=E->sphere.caps_per_proc; j++) {
V[1] = E->sphere.cap[j].V[1];
V[2] = E->sphere.cap[j].V[2];
V[3] = E->sphere.cap[j].V[3];
- for(m=1;m <= noz;m++)
- for(k=1;k <= noy;k++)
- for(l=1;l <= nox;l++){
- i = m + (l-1) * noz + (k-1) * noxnoz;
-
- sint = E->SinCos[lev][j][0][i];
- sinf = E->SinCos[lev][j][1][i];
- cost = E->SinCos[lev][j][2][i];
- cosf = E->SinCos[lev][j][3][i];
-
- fprintf(fp, "%.6e %.6e %.6e\n",
- V[1][i]*cost*cosf - V[2][i]*sinf + V[3][i]*sint*cosf,
- V[1][i]*cost*sinf + V[2][i]*cosf + V[3][i]*sint*sinf,
- -V[1][i]*sint + V[3][i]*cost);
- }
+ for(i=1; i<=E->lmesh.nno; i++) {
+ sint = E->SinCos[lev][j][0][i];
+ sinf = E->SinCos[lev][j][1][i];
+ cost = E->SinCos[lev][j][2][i];
+ cosf = E->SinCos[lev][j][3][i];
+
+ floatvel[(((j-1)*E->sphere.caps_per_proc)+i-1)*3+0] = (float)(V[1][i]*cost*cosf - V[2][i]*sinf + V[3][i]*sint*cosf);
+ floatvel[(((j-1)*E->sphere.caps_per_proc)+i-1)*3+1] = (float)(V[1][i]*cost*sinf + V[2][i]*cosf + V[3][i]*sint*sinf);
+ floatvel[(((j-1)*E->sphere.caps_per_proc)+i-1)*3+2] = (float)(-V[1][i]*sint + V[3][i]*cost);
+ }
}
+ if (strcmp(E->output.vtk_format, "binary") == 0)
+ write_binary_array(nodes*3,floatvel,fp);
+ else
+ write_ascii_array(nodes*3,3,floatvel,fp);
fputs(" </DataArray>\n", fp);
+
+ free(floatvel);
return;
}
static void vtk_output_visc(struct All_variables *E, FILE *fp)
{
- int i, j;
+ int nodes = E->sphere.caps_per_proc*E->lmesh.nno;
int lev = E->mesh.levmax;
- int m,k,l,nox,noy,noz,noxnoz;
- noy=E->lmesh.noy;
- nox=E->lmesh.nox;
- noz=E->lmesh.noz;
- noxnoz = nox * noz;
- fputs(" <DataArray type=\"Float32\" Name=\"viscosity\" format=\"ascii\">\n", fp);
+ fprintf(fp, " <DataArray type=\"Float32\" Name=\"viscosity\" format=\"%s\">\n", E->output.vtk_format);
+ if (strcmp(E->output.vtk_format, "binary") == 0) {
+ write_binary_array(nodes,&E->VI[lev][1][1],fp);
+ } else {
+ write_ascii_array(nodes,1,&E->VI[lev][1][1],fp);
+ }
- for(j=1; j<=E->sphere.caps_per_proc; j++) {
- for(m=1;m <= noz;m++)
- for(k=1;k <= noy;k++)
- for(l=1;l <= nox;l++){
- i = m + (l-1) * noz + (k-1) * noxnoz;
- fprintf(fp, "%.4e\n", E->VI[lev][j][i]);
- }
- }
-
fputs(" </DataArray>\n", fp);
return;
}
@@ -199,72 +193,550 @@
/* Output Cartesian coordinates as most VTK visualization softwares
assume it. */
int i, j;
- int m,k,l,nox,noy,noz,noxnoz;
- noy=E->lmesh.noy;
- nox=E->lmesh.nox;
- noz=E->lmesh.noz;
- noxnoz = nox * noz;
+ int nodes = E->sphere.caps_per_proc*E->lmesh.nno;
+ float* floatpos = malloc(nodes*3*sizeof(float));
fputs(" <Points>\n", fp);
- fputs(" <DataArray type=\"Float32\" Name=\"coordinate\" NumberOfComponents=\"3\" format=\"ascii\">\n", fp);
+ fprintf(fp, " <DataArray type=\"Float32\" Name=\"coordinate\" NumberOfComponents=\"3\" format=\"%s\">\n", E->output.vtk_format);
for(j=1; j<=E->sphere.caps_per_proc; j++) {
- for(m=1;m <= noz;m++)
- for(k=1;k <= noy;k++)
- for(l=1;l <= nox;l++){
- i = m + (l-1) * noz + (k-1) * noxnoz;
-
- fprintf(fp,"%.6e %.6e %.6e\n",
- E->x[j][1][i],
- E->x[j][2][i],
- E->x[j][3][i]);
- }
+ for(i=1; i<=E->lmesh.nno; i++){
+ floatpos[((j-1)*E->lmesh.nno+i-1)*3] = (float)(E->x[j][1][i]);
+ floatpos[((j-1)*E->lmesh.nno+i-1)*3+1]=(float)(E->x[j][2][i]);
+ floatpos[((j-1)*E->lmesh.nno+i-1)*3+2]=(float)(E->x[j][3][i]);
+ }
}
+ if (strcmp(E->output.vtk_format, "binary") == 0)
+ write_binary_array(nodes*3,floatpos,fp);
+ else
+ write_ascii_array(nodes*3,3,floatpos,fp);
fputs(" </DataArray>\n", fp);
fputs(" </Points>\n", fp);
+ free(floatpos);
+ return;
+}
+static void vtk_output_stress(struct All_variables *E, FILE *fp)
+{
+ int nodes = E->sphere.caps_per_proc*E->lmesh.nno;
+ /* for stress computation */
+ void allocate_STD_mem();
+ void compute_nodal_stress();
+ void free_STD_mem();
+ float *SXX[NCS],*SYY[NCS],*SXY[NCS],*SXZ[NCS],*SZY[NCS],*SZZ[NCS];
+ float *divv[NCS],*vorv[NCS];
+
+ /* those are sorted like stt spp srr stp str srp */
+ allocate_STD_mem(E, SXX, SYY, SZZ, SXY, SXZ, SZY, divv, vorv);
+ compute_nodal_stress(E, SXX, SYY, SZZ, SXY, SXZ, SZY, divv, vorv);
+ free_STD_mem(E, SXX, SYY, SZZ, SXY, SXZ, SZY, divv, vorv);
+
+ fprintf(fp, " <DataArray type=\"Float32\" Name=\"stress\" NumberOfComponents=\"6\" format=\"%s\">\n", E->output.vtk_format);
+
+ if (strcmp(E->output.vtk_format, "binary") == 0) {
+ write_binary_array(nodes*6,&E->gstress[1][1],fp);
+ } else {
+ write_ascii_array(nodes*6,6,&E->gstress[1][1],fp);
+ }
+
+ fputs(" </DataArray>\n", fp);
return;
}
+static void vtk_output_comp_nd(struct All_variables *E, FILE *fp)
+{
+ int i, j, k;
+ char name[255];
+ int nodes = E->sphere.caps_per_proc*E->lmesh.nno;
+ float* floatcompo = malloc (nodes*sizeof(float));
+ for(k=0;k<E->composition.ncomp;k++) {
+ fprintf(fp, " <DataArray type=\"Float32\" Name=\"composition%d\" format=\"%s\">\n", k+1, E->output.vtk_format);
+
+ for(j=1; j<=E->sphere.caps_per_proc; j++) {
+ for(i=1; i<=E->lmesh.nno; i++) {
+ floatcompo[(j-1)*E->lmesh.nno+i-1] = (float) (E->composition.comp_node[j][k][i]);
+ }
+ }
+
+ if (strcmp(E->output.vtk_format, "binary") == 0)
+ write_binary_array(nodes,floatcompo,fp);
+ else
+ write_ascii_array(nodes,1,floatcompo,fp);
+ fputs(" </DataArray>\n", fp);
+ }
+ free(floatcompo);
+ return;
+}
+
+
+static void vtk_output_surf(struct All_variables *E, FILE *fp, int cycles)
+{
+ int i, j, k;
+ int nodes = E->sphere.caps_per_proc*E->lmesh.nno;
+ char output_file[255];
+ float* floattopo = malloc (nodes*sizeof(float));
+
+ if((E->output.write_q_files == 0) || (cycles == 0) ||
+ (cycles % E->output.write_q_files)!=0)
+ heat_flux(E);
+ /* else, the heat flux will have been computed already */
+
+ if(E->control.use_cbf_topo){
+ get_CBF_topo(E,E->slice.tpg,E->slice.tpgb);
+ }
+ else{
+ get_STD_topo(E,E->slice.tpg,E->slice.tpgb,E->slice.divg,E->slice.vort,cycles);
+ }
+
+ fprintf(fp," <DataArray type=\"Float32\" Name=\"surface\" format=\"%s\">\n", E->output.vtk_format);
+
+ for(j=1;j<=E->sphere.caps_per_proc;j++){
+ for(i=1;i<=E->lmesh.nsf;i++){
+ for(k=1;k<=E->lmesh.noz;k++){
+ floattopo[(j-1)*E->lmesh.nno + (i-1)*E->lmesh.noz + k-1] = 0.0;
+ }
+
+ if (E->parallel.me_loc[3]==E->parallel.nprocz-1) {
+
+ /* choose either STD topo or pseudo-free-surf topo */
+ if(E->control.pseudo_free_surf)
+ floattopo[(j-1)*E->lmesh.nno + i*E->lmesh.noz-1] = E->slice.freesurf[j][i];
+ else
+ floattopo[(j-1)*E->lmesh.nno + i*E->lmesh.noz-1] = E->slice.tpg[j][i];
+
+ }
+ }
+ }
+
+ if (strcmp(E->output.vtk_format, "binary") == 0)
+ write_binary_array(nodes,floattopo,fp);
+ else
+ write_ascii_array(nodes,1,floattopo,fp);
+
+ fputs(" </DataArray>\n", fp);
+ return;
+}
+
+
+static void write_vtm(struct All_variables *E, int cycles)
+{
+ FILE *fp;
+ char vtm_file[255];
+ int n;
+
+ const char header[] =
+ "<?xml version=\"1.0\"?>\n"
+ "<VTKFile type=\"vtkMultiBlockDataSet\" version=\"1.0\" compressor=\"vtkZLibDataCompressor\" byte_order=\"LittleEndian\">\n"
+ " <vtkMultiBlockDataSet>\n";
+
+ snprintf(vtm_file, 255, "%s.%d.vtm",
+ E->control.data_file, cycles);
+ fp = output_open(vtm_file, "w");
+ fputs(header, fp);
+
+ for(n=0; n<E->parallel.nproc; n++) {
+ fprintf(fp, " <DataSet index=\"%d\" file=\"%s.proc%d.%d.vts\"/>\n",
+ n, E->control.data_prefix, n, cycles);
+ }
+ fputs(" </vtkMultiBlockDataSet>\n",fp);
+ fputs("</VTKFile>",fp);
+
+ fclose(fp);
+}
+
+static void write_visit(struct All_variables *E, int cycles)
+{
+ FILE *fp;
+ char visit_file[255];
+ int n;
+
+ const char header[] = "!NBLOCKS %d\n";
+
+ snprintf(visit_file, 255, "%s.%d.visit",
+ E->control.data_file, cycles);
+ fp = output_open(visit_file, "w");
+ fprintf(fp, header, E->parallel.nproc);
+
+ for(n=0; n<E->parallel.nproc; n++) {
+ fprintf(fp, "%s.proc%d.%d.vts\n",
+ E->control.data_prefix, n, cycles);
+ }
+ fclose(fp);
+}
+
+static void write_pvts(struct All_variables *E, int cycles)
+{
+ FILE *fp;
+ char pvts_file[255];
+ int i,j,k;
+ snprintf(pvts_file, 255, "%s.%d.pvts",
+ E->control.data_file,cycles);
+ fp = output_open(pvts_file, "w");
+
+ const char format[] =
+ "<?xml version=\"1.0\"?>\n"
+ "<VTKFile type=\"PStructuredGrid\" version=\"0.1\" compressor=\"vtkZLibDataCompressor\" byte_order=\"LittleEndian\">\n"
+ " <PStructuredGrid WholeExtent=\"%s\" GhostLevel=\"#\">\n"
+ " <PPointData Scalars=\"temperature\" Vectors=\"velocity\">\n"
+ " <DataArray type=\"Float32\" Name=\"temperature\" format=\"%s\"/>\n"
+ " <DataArray type=\"Float32\" Name=\"velocity\" NumberOfComponents=\"3\" format=\"%s\"/>\n"
+ " <DataArray type=\"Float32\" Name=\"viscosity\" format=\"%s\"/>\n";
+
+ char extent[64], header[1024];
+
+ snprintf(extent, 64, "%d %d %d %d %d %d",
+ E->lmesh.ezs, E->lmesh.ezs + E->lmesh.elz*E->parallel.nprocz,
+ E->lmesh.exs, E->lmesh.exs + E->lmesh.elx*E->parallel.nprocx,
+ E->lmesh.eys, E->lmesh.eys + E->lmesh.ely*E->parallel.nprocy);
+
+ snprintf(header, 1024, format, extent, E->output.vtk_format,
+ E->output.vtk_format, E->output.vtk_format);
+ fputs(header, fp);
+
+ if (E->output.stress){
+ fprintf(fp," <DataArray type=\"Float32\" Name=\"stress\" NumberOfComponents=\"6\" format=\"%s\"/>\n", E->output.vtk_format);
+ }
+ if (E->output.comp_nd && E->composition.on){
+ fprintf(fp," <DataArray type=\"Float32\" Name=\"composition1\" format=\"%s\"/>\n", E->output.vtk_format);
+ }
+ if (E->output.surf){
+ fprintf(fp," <DataArray type=\"Float32\" Name=\"surface\" format=\"%s\"/>\n", E->output.vtk_format);
+ }
+
+ fputs(" </PPointData>\n \n"
+ " <PCellData>\n"
+ " </PCellData>\n \n"
+ " <PPoints>\n"
+ " <DataArray type=\"Float32\" Name=\"coordinate\" NumberOfComponents=\"3\" format=\"binary\" />\n"
+ " </PPoints>\n", fp);
+
+ for(i=0; i < E->parallel.nprocy;i++){
+ for(j=0; j < E->parallel.nprocx;j++){
+ for(k=0; k < E->parallel.nprocz;k++){
+ fprintf(fp, " <Piece Extent=\"%d %d %d %d %d %d\" Source=\"%s.proc%d.%d.vts\"/>\n",
+ (k%E->parallel.nprocz)*E->lmesh.elz,
+ (k%E->parallel.nprocz+1)*E->lmesh.elz,
+ (j%E->parallel.nprocx)*E->lmesh.elx, (j%E->parallel.nprocx+1)*E->lmesh.elx,
+ (i%E->parallel.nprocy)*E->lmesh.ely, (i%E->parallel.nprocy+1)*E->lmesh.ely,
+ E->control.data_prefix,
+ i*E->parallel.nprocx*E->parallel.nprocz+j*E->parallel.nprocz+k, cycles);
+ }
+ }
+ }
+
+ fputs(" </PStructuredGrid>\n",fp);
+ fputs("</VTKFile>",fp);
+
+ fclose(fp);
+}
+
+static void write_ascii_array(int nn, int perLine, float *array, FILE *fp)
+{
+ int i;
+
+ switch (perLine) {
+ case 1:
+ for(i=0; i<nn; i++)
+ fprintf(fp, "%.4e\n", array[i]);
+ break;
+ case 3:
+ for(i=0; i < nn/3; i++)
+ fprintf(fp,"%.4e %.4e %.4e\n",array[3*i],array[3*i+1],array[3*i+2]);
+ break;
+ case 6:
+ for(i=0; i < nn/6; i++)
+ fprintf(fp,"%.4e %.4e %.4e %.4e %.4e %.4e\n",
+ array[6*i],array[6*i+1],array[6*i+2],
+ array[6*i+3],array[6*i+4],array[6*i+5]);
+ break;
+ }
+ return;
+}
+
+static void FloatToUnsignedChar(float * floatarray, int nn, unsigned char * chararray)
+{
+ /* simple float to unsigned chararray routine via union
+ nn=length(intarray) chararray has to be BIG ENOUGH! */
+ int i;
+ union FloatToUnsignedChars
+ {
+ float input;
+ unsigned char output[4];
+ } floattransform;
+
+ for (i=0; i<nn; i++){
+ floattransform.input=floatarray[i];
+ chararray[4*i]=floattransform.output[0];
+ chararray[4*i+1]=floattransform.output[1];
+ chararray[4*i+2]=floattransform.output[2];
+ chararray[4*i+3]=floattransform.output[3];
+ }
+ return;
+}
+
+static void IntToUnsignedChar(int * intarray, int nn, unsigned char * chararray)
+{
+ /* simple int - to unsigned chararray routine via union
+ nn=length(intarray) chararray has to be BIG ENOUGH! */
+ int i;
+ union IntToUnsignedChars
+ {
+ int input;
+ unsigned char output[4];
+ } inttransform;
+
+ for (i=0; i<nn; i++){
+ inttransform.input=intarray[i];
+ chararray[4*i]=inttransform.output[0];
+ chararray[4*i+1]=inttransform.output[1];
+ chararray[4*i+2]=inttransform.output[2];
+ chararray[4*i+3]=inttransform.output[3];
+ }
+}
+
+
+static void zlibcompress(unsigned char* in, int nn, unsigned char** out, int *nn2)
+/* function to compress "in" to "out" reducing size from nn to nn2 */
+{
+#ifdef USE_GZDIR
+ int ntemp=0;
+
+ /* in and out of z-stream */
+ unsigned char inz[CHUNK];
+ unsigned char outz[CHUNK];
+
+ /* compression level */
+ int level = Z_DEFAULT_COMPRESSION;
+ int ret,flush;
+ int i,j,k;
+
+ /* zlib compression stream */
+ z_stream strm;
+
+ /* hope compressed data will be <= uncompressed */
+ *out = malloc(sizeof(unsigned char)*nn);
+
+ strm.zalloc = Z_NULL;
+ strm.zfree = Z_NULL;
+ strm.opaque = Z_NULL;
+
+ /* zlib init */
+ ret = deflateInit(&strm, level);
+ if (ret == Z_OK){
+ i=0; // position in "in" array
+ do{
+ j=0; // position in "inz"
+ do{
+ inz[j++]=in[i++];
+ } while((j<CHUNK) && (i<nn)); // stopps if "inz"-buffer is full or "in" array empty
+ strm.avail_in=j; // set number of input chars
+
+ flush = (i==nn) ? Z_FINISH : Z_NO_FLUSH; // done?
+ strm.next_in = inz; // set input buffer
+
+ do{
+ strm.avail_out = CHUNK; // set number of max output chars
+ strm.next_out = outz; // set output buffer
+
+ /* zlib compress */
+ ret = deflate(&strm, flush);
+ assert(ret != Z_STREAM_ERROR);
+
+ /* zlib changed strm.avail_out=CHUNK
+ to the number of chars we can NOT use
+ in outz */
+
+ for (k=0;k<CHUNK-strm.avail_out;k++){
+ (*out)[ntemp+k]=outz[k];
+ }
+
+ /* increase position in "out" */
+ ntemp+=(CHUNK-strm.avail_out);
+ }while(strm.avail_out==0);
+ assert(strm.avail_in == 0);
+
+ }while (flush != Z_FINISH);
+ }
+ else{fprintf(stderr,"Error during compression init\n");}
+
+ // now we know how short "out" should be!
+ *nn2=ntemp;
+ *out = realloc(*out,sizeof(unsigned char)*ntemp);
+
+ (void)deflateEnd(&strm);
+#endif
+
+ return;
+}
+
+static void base64(unsigned char * in, int nn, unsigned char* out)
+{
+ /*takes *in*-array and "in"-length-"nn" and fills "out"-array
+ with base64(in) "out" needs to be big enough!!!
+ length(out) >= 4* |^ nn/3.0 ^| */
+ char cb64[]="ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789+/";
+ int len;
+ int i;
+
+ for (i=0; i < nn; i+=3){
+
+ len = (3 < nn-i ? 3 : nn-i);
+ if (len >= 3){
+ /* normal base64 encoding */
+ out[i/3*4+0] = cb64[ in[i] >> 2 ];
+ out[i/3*4+1] = cb64[ ((in[i] & 0x03) << 4) | ((in[i+1] & 0xf0) >> 4) ];
+ out[i/3*4+2] = cb64[ ((in[i+1] & 0x0f) << 2) | ((in[i+2] & 0xc0) >> 6)];
+ out[i/3*4+3] = cb64[ in[i+2] & 0x3f ];
+ } else if (len == 2){
+ /* at the end of array fill up with '=' */
+ out[i/3*4+0] = cb64[ in[i] >> 2 ];
+ out[i/3*4+1] = cb64[ ((in[i] & 0x03) << 4) | ((in[i+1] & 0xf0) >> 4) ];
+ out[i/3*4+2] = cb64[((in[i+1] & 0x0f) << 2)];
+ out[i/3*4+3] = (unsigned char) '=';
+ } else if (len == 1){
+ /* at the end of array fill up with '=' */
+ out[i/3*4+0] = cb64[ in[i] >> 2 ];
+ out[i/3*4+1] = cb64[ ((in[i] & 0x03) << 4) ];
+ out[i/3*4+2] = (unsigned char) '=';
+ out[i/3*4+3] = (unsigned char) '=';
+ }
+ }
+}
+
+
+static void base64plushead(unsigned char * in, int nn, int orinn, unsigned char* out)
+{
+ /* writing vtk compatible zlib compressed base64 encoded data to "out" */
+ int i;
+ unsigned char * b64head;
+ int b64bodylength;
+ unsigned char * b64body;
+ /* header of data */
+ unsigned char * charhead = malloc(sizeof(unsigned char)*16);
+ /* - consists of "1" (number of pieces) */
+ /* - original datalength in byte */
+ /* - original datalength in byte */
+ /* - new datalength after z-lib compression */
+ int * headInts= malloc(sizeof(int)*4);
+ headInts[0]=1;
+ headInts[1]=orinn;
+ headInts[2]=orinn;
+ headInts[3]=nn;
+ // transform to unsigned char
+ IntToUnsignedChar(headInts,4,charhead);
+
+ // base64: 16byte -> 24byte
+ b64head = malloc(sizeof(unsigned char)*24);
+ // fills b64head
+ base64(charhead, 16, b64head);
+
+ // base64 data
+ b64bodylength = 4*ceil((float) nn/3.0);
+ b64body = malloc(sizeof(unsigned char)*b64bodylength);
+ // writes base64 data to b64body
+ base64(in,nn,b64body);
+
+ // combines header and body
+ for (i=0; i<24 ; i++){
+ out[i]=b64head[i];
+ }
+
+ for (i=0; i<b64bodylength ; i++){
+ out[24+i]=b64body[i];
+ }
+
+ if(b64body){free(b64body);}
+ if(b64head){free(b64head);}
+ if(headInts){free(headInts);}
+ if(charhead){free(charhead);}
+}
+
+static void write_binary_array(int nn, float* array, FILE * f)
+{
+ /* writes vtk-data array of floats and performs zip and base64 encoding */
+ int chararraylength=4*nn; /* nn floats -> 4*nn unsigned chars */
+ unsigned char * chararray = malloc (chararraylength * sizeof(unsigned char));
+ int compressedarraylength = 0;
+ unsigned char * compressedarray;
+ unsigned char ** pointertocompressedarray= &compressedarray;
+ int base64plusheadlength;
+ unsigned char * base64plusheadarray;
+
+ FloatToUnsignedChar(array,nn,chararray);
+
+ /* compression routine */
+ zlibcompress(chararray,chararraylength,pointertocompressedarray,&compressedarraylength);
+
+ /* special header for zip compressed and bas64 encoded data
+ header needs 4 int32 = 16 byte -> 24 byte due to base64 (4*16/3) */
+ base64plusheadlength = 24 + 4*ceil((float) compressedarraylength/3.0);
+ base64plusheadarray = malloc(sizeof(unsigned char)* base64plusheadlength);
+
+ /* fills base64plusheadarray with everything ready for simple writing */
+ base64plushead(compressedarray,compressedarraylength, chararraylength, base64plusheadarray);
+
+ fwrite(base64plusheadarray,sizeof(unsigned char),base64plusheadlength,f);
+ fprintf(f,"\n");
+ free(chararray);
+ free(base64plusheadarray);
+}
+
/**********************************************************************/
void vtk_output(struct All_variables *E, int cycles)
{
char output_file[255];
FILE *fp;
-
- snprintf(output_file, 255, "%s.%d.step%d.vts",
+ snprintf(output_file, 255, "%s.proc%d.%d.vts",
E->control.data_file, E->parallel.me, cycles);
fp = output_open(output_file, "w");
-
/* first, write volume data to vts file */
vts_file_header(E, fp);
/* write node-based field */
vtk_point_data_header(E, fp);
vtk_output_temp(E, fp);
+
vtk_output_velo(E, fp);
+
vtk_output_visc(E, fp);
+
+ if (E->output.stress)
+ vtk_output_stress(E, fp);
+
+ if (E->output.comp_nd && E->composition.on)
+ vtk_output_comp_nd(E, fp);
+
+ if (E->output.surf)
+ vtk_output_surf(E, fp, cycles);
+
vtk_point_data_trailer(E, fp);
/* write element-based field */
vtk_cell_data_header(E, fp);
- /**/
+ /* TODO: comp_el, heating */
vtk_cell_data_trailer(E, fp);
/* write coordinate */
vtk_output_coord(E, fp);
vts_file_trailer(E, fp);
+ fclose(fp);
/* then, write other type of data */
- //vtk_output_surf_botm(E, );
-
- fclose(fp);
-
+ if (E->parallel.me == 0) {
+ if (E->sphere.caps == 12) {
+ /* VTK multiblock file */
+ write_vtm(E, cycles);
+ /* LLNL VisIt multiblock file */
+ write_visit(E, cycles);
+ }
+ else
+ write_pvts(E, cycles);
+ }
return;
}
Modified: mc/3D/CitcomS/trunk/lib/global_defs.h
===================================================================
--- mc/3D/CitcomS/trunk/lib/global_defs.h 2011-07-12 15:18:21 UTC (rev 18741)
+++ mc/3D/CitcomS/trunk/lib/global_defs.h 2011-07-12 16:27:58 UTC (rev 18742)
@@ -624,6 +624,7 @@
struct Output {
char format[20]; /* ascii or hdf5 */
char optional[1000]; /* comma-delimited list of objects to output */
+ char vtk_format[10]; /*ascii or binary */
int llmax; /* max degree of spherical harmonics output */
Modified: mc/3D/CitcomS/trunk/module/setProperties.c
===================================================================
--- mc/3D/CitcomS/trunk/module/setProperties.c 2011-07-12 15:18:21 UTC (rev 18741)
+++ mc/3D/CitcomS/trunk/module/setProperties.c 2011-07-12 16:27:58 UTC (rev 18742)
@@ -311,6 +311,7 @@
getStringProperty(properties, "output_format", E->output.format, fp);
getStringProperty(properties, "output_optional", E->output.optional, fp);
+ getStringProperty(properties, "vtk_format", E->output.vtk_format, fp);
getIntProperty(properties, "gzdir_vtkio", E->output.gzdir.vtk_io, fp);
getIntProperty(properties, "gzdir_rnr", E->output.gzdir.rnr, fp);
More information about the CIG-COMMITS
mailing list