[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