<p><b>cahrens@lanl.gov</b> 2010-08-16 16:07:43 -0600 (Mon, 16 Aug 2010)</p><p>Added xml/binary .vtu outputting versions: transdualbin and projlatlonbin<br>
</p><hr noshade><pre><font color="gray">Added: branches/ocean_projects/graphics/paraview/projection/EncodeBase64.cpp
===================================================================
--- branches/ocean_projects/graphics/paraview/projection/EncodeBase64.cpp (rev 0)
+++ branches/ocean_projects/graphics/paraview/projection/EncodeBase64.cpp 2010-08-16 22:07:43 UTC (rev 477)
@@ -0,0 +1,169 @@
+using namespace std;
+
+#include "EncodeBase64.h"
+
+
+EncodeBase64::EncodeBase64(ofstream* astream)
+{
+ outstream = astream;
+ this->BufferLength = 0;
+}
+
+//----------------------------------------------------------------------------
+EncodeBase64::~EncodeBase64()
+{
+}
+
+
+
+//----------------------------------------------------------------------------
+static const unsigned char Base64EncodeTable[65] =
+"ABCDEFGHIJKLMNOPQRSTUVWXYZ"
+"abcdefghijklmnopqrstuvwxyz"
+"0123456789+/";
+
+//----------------------------------------------------------------------------
+inline static unsigned char EncodeChar(unsigned char c)
+{
+ assert( c < 65 );
+ return Base64EncodeTable[c];
+}
+
+//----------------------------------------------------------------------------
+void EncodeBase64::EncodePair(unsigned char i0,
+ unsigned char i1,
+ unsigned char *o0,
+ unsigned char *o1,
+ unsigned char *o2,
+ unsigned char *o3)
+{
+ *o0 = EncodeChar((i0 >> 2) & 0x3F);
+ *o1 = EncodeChar(((i0 << 4) & 0x30)|((i1 >> 4) & 0x0F));
+ *o2 = EncodeChar(((i1 << 2) & 0x3C));
+ *o3 = '=';
+}
+
+//----------------------------------------------------------------------------
+void EncodeBase64::EncodeSingle(unsigned char i0,
+ unsigned char *o0,
+ unsigned char *o1,
+ unsigned char *o2,
+ unsigned char *o3)
+{
+ *o0 = EncodeChar((i0 >> 2) & 0x3F);
+ *o1 = EncodeChar(((i0 << 4) & 0x30));
+ *o2 = '=';
+ *o3 = '=';
+}
+
+//----------------------------------------------------------------------------
+void EncodeBase64::EncodeTriplet(unsigned char i0,
+ unsigned char i1,
+ unsigned char i2,
+ unsigned char *o0,
+ unsigned char *o1,
+ unsigned char *o2,
+ unsigned char *o3)
+{
+ *o0 = EncodeChar((i0 >> 2) & 0x3F);
+ *o1 = EncodeChar(((i0 << 4) & 0x30)|((i1 >> 4) & 0x0F));
+ *o2 = EncodeChar(((i1 << 2) & 0x3C)|((i2 >> 6) & 0x03));
+ *o3 = EncodeChar(i2 & 0x3F);
+}
+
+//----------------------------------------------------------------------------
+int EncodeBase64::EncodeTriplet(unsigned char c0,
+ unsigned char c1,
+ unsigned char c2)
+{
+ // Encodes 3 bytes into 4 bytes and writes them to the output stream.
+ unsigned char out[4];
+ EncodeTriplet(c0, c1, c2,
+ &out[0], &out[1], &out[2], &out[3]);
+ return (outstream->write(reinterpret_cast<char*>(out), 4)? 1:0);
+}
+
+//----------------------------------------------------------------------------
+int EncodeBase64::EncodeEnding(unsigned char c0,
+ unsigned char c1)
+{
+ // Encodes a 2-byte ending into 3 bytes and 1 pad byte and writes.
+ unsigned char out[4];
+ EncodePair(c0, c1,
+ &out[0], &out[1], &out[2], &out[3]);
+ return (outstream->write(reinterpret_cast<char*>(out), 4)? 1:0);
+}
+
+//----------------------------------------------------------------------------
+int EncodeBase64::EncodeEnding(unsigned char c0)
+{
+ // Encodes a 1-byte ending into 2 bytes and 2 pad bytes and writes.
+ unsigned char out[4];
+ EncodeSingle(c0, &out[0], &out[1], &out[2], &out[3]);
+ return (outstream->write(reinterpret_cast<char*>(out), 4)? 1:0);
+}
+
+//----------------------------------------------------------------------------
+int EncodeBase64::StartWriting()
+{
+ this->BufferLength = 0;
+ this->WriteCount = 0;
+ //cout << "StartWriting!" << endl;
+ return 1;
+}
+
+//----------------------------------------------------------------------------
+int EncodeBase64::EndWriting()
+{
+ if(this->BufferLength == 1)
+ {
+ if(!this->EncodeEnding(this->Buffer[0])) { return 0; }
+ this->BufferLength = 0;
+ }
+ else if(this->BufferLength == 2)
+ {
+ if(!this->EncodeEnding(this->Buffer[0], this->Buffer[1])) { return 0; }
+ this->BufferLength = 0;
+ }
+ //cout << "EndWriting: WriteCount: " << WriteCount << endl;
+ return 1;
+}
+
+//----------------------------------------------------------------------------
+int EncodeBase64::Write(const unsigned char* data,
+ unsigned long length)
+{
+ unsigned long totalLength = this->BufferLength + length;
+ const unsigned char* in = data;
+ const unsigned char* end = data+length;
+
+ if(totalLength >= 3)
+ {
+ if(this->BufferLength == 1)
+ {
+ if(!this->EncodeTriplet(this->Buffer[0], in[0], in[1])) { return 0; }
+ in += 2;
+ this->BufferLength = 0;
+ }
+ else if(this->BufferLength == 2)
+ {
+ if(!this->EncodeTriplet(this->Buffer[0], this->Buffer[1], in[0]))
+ { return 0; }
+ in += 1;
+ this->BufferLength = 0;
+ }
+ }
+
+ while((end - in) >= 3)
+ {
+ if(!this->EncodeTriplet(in[0], in[1], in[2])) { return 0; }
+ in += 3;
+ }
+
+ while(in != end)
+ {
+ this->Buffer[this->BufferLength++] = *in++;
+ }
+ WriteCount++;
+ return 1;
+}
Property changes on: branches/ocean_projects/graphics/paraview/projection/EncodeBase64.cpp
___________________________________________________________________
Name: svn:mime-type
+ text/x-c++src
Name: svn:keywords
+ Author Date ID Revision
Name: svn:eol-style
+ native
Added: branches/ocean_projects/graphics/paraview/projection/EncodeBase64.h
===================================================================
--- branches/ocean_projects/graphics/paraview/projection/EncodeBase64.h (rev 0)
+++ branches/ocean_projects/graphics/paraview/projection/EncodeBase64.h 2010-08-16 22:07:43 UTC (rev 477)
@@ -0,0 +1,72 @@
+using namespace std;
+
+#include <iostream>
+#include <fstream>
+#include "stdlib.h"
+#include <assert.h>
+
+class EncodeBase64
+{
+public:
+ EncodeBase64(ofstream *outstream);
+ ~EncodeBase64();
+
+ // Description:
+ // Called after the stream position has been set by the caller, but
+ // before any Write calls. The stream position should not be
+ // adjusted by the caller until after an EndWriting call.
+ int StartWriting();
+
+ // Description:
+ // Write output data of the given length.
+ int Write(const unsigned char* data, unsigned long length);
+
+ // Description:
+ // Called after all desired calls to Write have been made. After
+ // this call, the caller is free to change the position of the
+ // stream. Additional writes should not be done until after another
+ // call to StartWriting.
+ int EndWriting();
+
+protected:
+
+ // Methods to encode and write data.
+
+ // Number of un-encoded bytes left in Buffer from last call to Write.
+ unsigned int BufferLength;
+ unsigned int WriteCount;
+ unsigned char Buffer[2];
+ ofstream *outstream;
+
+ // Description:
+ // Encode 2 bytes into 4 bytes
+ static void EncodePair(unsigned char i0,
+ unsigned char i1,
+ unsigned char *o0,
+ unsigned char *o1,
+ unsigned char *o2,
+ unsigned char *o3);
+
+ // Description:
+ // Encode 1 byte into 4 bytes
+ static void EncodeSingle(unsigned char i0,
+ unsigned char *o0,
+ unsigned char *o1,
+ unsigned char *o2,
+ unsigned char *o3);
+ // Description:
+ // Encode 3 bytes into 4 bytes
+ static void EncodeTriplet(unsigned char i0,
+ unsigned char i1,
+ unsigned char i2,
+ unsigned char *o0,
+ unsigned char *o1,
+ unsigned char *o2,
+ unsigned char *o3);
+
+ int EncodeTriplet(unsigned char c0, unsigned char c1, unsigned char c2);
+ int EncodeEnding(unsigned char c0, unsigned char c1);
+ int EncodeEnding(unsigned char c0);
+
+};
+
Property changes on: branches/ocean_projects/graphics/paraview/projection/EncodeBase64.h
___________________________________________________________________
Name: svn:mime-type
+ text/x-c++hdr
Name: svn:keywords
+ Author Date ID Revision
Name: svn:eol-style
+ native
Modified: branches/ocean_projects/graphics/paraview/projection/Makefile
===================================================================
--- branches/ocean_projects/graphics/paraview/projection/Makefile 2010-08-12 22:02:03 UTC (rev 476)
+++ branches/ocean_projects/graphics/paraview/projection/Makefile 2010-08-16 22:07:43 UTC (rev 477)
@@ -1,12 +1,12 @@
-NETCDF = /usr/projects/climate/mhecht/netcdf-3.6.1
+#NETCDF = /usr/projects/climate/mhecht/netcdf-3.6.1
CC = gcc
INCLUDES = -I$(NETCDF)/include
LIBS = -L$(NETCDF)/lib -lnetcdf_c++ -lnetcdf \
-lstdc++ \
- -L/opt/Intel/cce/10.0.023/lib/ -lirc
+# -L/opt/Intel/cce/10.0.023/lib/ -lirc
+# uncomment above if on coyote/lobo
-# -L/usr/lib -lstdc++ \
#MYDEBUG = -g
MYDEBUG =
@@ -25,5 +25,14 @@
proj3dquadrect: proj3dquadrect.cpp
$(CC) $(INCLUDES) -o $@ $< $(LIBS)
+.cpp.o:
+ $(CC) $(MYDEBUG) $(INCLUDES) -c $<
+
+EncodeBase64.o: EncodeBase64.cpp
+ $(CC) $(INCLUDES) -c $<
+
+projlatlonbin: projlatlonbin.cpp EncodeBase64.o convert.o
+ $(CC) $(INCLUDES) -o $@ $< EncodeBase64.o convert.o $(LIBS)
+
clean:
- rm *.o proj3dsphere projlatlon proj3dquadrect
+ rm *.o proj3dsphere projlatlon proj3dquadrect projlatlonbin
Added: branches/ocean_projects/graphics/paraview/projection/convert.c
===================================================================
--- branches/ocean_projects/graphics/paraview/projection/convert.c (rev 0)
+++ branches/ocean_projects/graphics/paraview/projection/convert.c 2010-08-16 22:07:43 UTC (rev 477)
@@ -0,0 +1,16 @@
+#include <math.h>
+
+void CartesianToSpherical(float x, float y, float z, float* rho, float* phi, float* theta)
+{
+ *rho = sqrtf( (x*x) + (y*y) + (z*z));
+ *theta = atan2f(y, x);
+ *phi = acosf(z/(*rho));
+}
+
+void SphericalToCartesian(float rho, float phi, float theta, float* x, float* y, float* z)
+{
+ *x = rho* sinf(phi) * cosf(theta);
+ *y = rho* sinf(phi) * sinf(theta);
+ *z = rho* cosf(phi);
+}
+
Property changes on: branches/ocean_projects/graphics/paraview/projection/convert.c
___________________________________________________________________
Name: svn:mime-type
+ text/x-csrc
Name: svn:keywords
+ Author Date ID Revision
Name: svn:eol-style
+ native
Added: branches/ocean_projects/graphics/paraview/projection/projlatlonbin.cpp
===================================================================
--- branches/ocean_projects/graphics/paraview/projection/projlatlonbin.cpp (rev 0)
+++ branches/ocean_projects/graphics/paraview/projection/projlatlonbin.cpp 2010-08-16 22:07:43 UTC (rev 477)
@@ -0,0 +1,1600 @@
+////////////////////////////////////////////////////////////////////////////////
+//
+// This program translates a newpop netCDF data file to dual-grid lat/lon multilayer
+// projection in legacy, binary .vtu format. Assumes little-endian platform.
+//
+// Assume all variables are of interest.
+// Assume variable data type is double.
+// Assume variable dims are (Time, nCells|nVertices, nVertLevels|nVertLevelsP1)
+// Assume no more than 100 vars each for cell and point data
+// Deals with vectors.
+// Does not deal with edge data.
+// Only vis up to nVertLevels, not nVertLevelsP1.
+// Doubles converted to floats in .vtk files follow these rules:
+// if a NaN, then become -FLT_MAX.
+// if positive infinity, then become FLT_MAX
+// if negative infinity, then become -FLT_MAX
+// if smaller than a float, then become 0.
+// if not zero, but between -FLT_MIN and FLT_MIN make -FLT_MIN or FLT_MIN
+// if outside of -FLT_MAX and FLT_MAX make -FLT_MAX or FLT_MAX
+//
+// Version 1.3
+// Christine Ahrens
+// 7/26/2010
+////////////////////////////////////////////////////////////////////////////////
+
+
+#include <iostream>
+#include <fstream>
+#include <sstream>
+#include "stdlib.h"
+#include "vtkCellType.h"
+#include "netcdfcpp.h"
+#include <string>
+#include <cmath>
+//#include <math>
+#include <cfloat>
+#include "EncodeBase64.h"
+
+using namespace std;
+
+#define CHECK_MALLOC(ptr) \
+ if (ptr == NULL) { \
+ cerr << "malloc failed!</font>
<font color="blue">"; \
+ exit(1); \
+ }
+
+#define CHECK_DIM(ncFile, name) \
+ if (!isNcDim(ncFile, name)) { \
+ cerr << "Cannot find dimension: " << name << endl; \
+ exit(1); \
+ }
+
+#define CHECK_VAR(ncFile, name) \
+ if (!isNcVar(ncFile, name)) { \
+ cerr << "Cannot find variable: " << name << endl; \
+ exit(1); \
+ }
+
+#define MAX_VARS 100
+#define DEFAULT_LAYER_THICKNESS 10
+using namespace std;
+
+#define BUFF_SIZE 2048
+
+// append inFile to outFile
+int myAppendFile(string *inFile, string *outFile)
+{
+ char buff[BUFF_SIZE];
+ int readBytes = 1;
+
+ ifstream inFileStream(inFile->c_str(), ios::in|ios::binary);
+ if(!inFileStream)
+ {
+ return -1;
+ }
+
+ ofstream outFileStream(outFile->c_str(), ios::out|ios::app|ios::binary);
+ if(!outFileStream)
+ {
+ return -1;
+ }
+
+ while(readBytes != 0)
+ {
+ inFileStream.read((char*)buff, BUFF_SIZE);
+ readBytes = inFileStream.gcount();
+ outFileStream.write((char*)buff, readBytes);
+ }
+ return 0;
+}
+
+int my_isinf_f(float x){
+ if (isinf(x)) {
+ return (x < 0.0 ? -1 : 1);
+ } else return 0;
+}
+
+unsigned int padSize(unsigned int binarySize) {
+ switch ((binarySize*4)%3) {
+ case 1:
+ binarySize = binarySize + 2;
+ break;
+ case 2:
+ binarySize++;
+ break;
+ case 0:
+ break;
+ }
+ return binarySize;
+}
+
+float convertDouble2ValidFloat(double inputData) {
+
+ // check for NaN
+ if (inputData != inputData) {
+ //cerr << "found NaN!" << endl;
+ return -FLT_MAX;
+ }
+
+ // check for infinity
+ int retval = my_isinf_f((float)inputData);
+ if (retval < 0) {
+ return -FLT_MAX;
+ } else if (retval > 0) {
+ return FLT_MAX;
+ }
+
+ // check number too small for float
+ if (abs(inputData) < 1e-126) return 0.0;
+
+ if ((float)inputData == 0) return 0.0;
+
+ if ((abs(inputData) > 0) && (abs(inputData) < FLT_MIN)) {
+ if (inputData < 0) return -FLT_MIN; else return FLT_MIN;
+ }
+
+ if (abs(inputData) > FLT_MAX) {
+ if (inputData < 0) return -FLT_MAX; else return FLT_MAX;
+ }
+
+ return (float)inputData;
+}
+
+double convertXY2Lon(double x, double y, double lon) {
+ return ((-x * sin(lon)) + (y*cos(lon)));
+}
+
+double convertXYZ2Lat(double x, double y, double z, double lon, double lat) {
+ return ((((x*cos(lon)) + (y*sin(lon)))*sin(lat)) + (z*cos(lat)));
+}
+
+void writeFloat(EncodeBase64 &b64, ofstream &file, double data, bool writeBinary) {
+ float validData = convertDouble2ValidFloat (data);
+ if (writeBinary) {
+ b64.Write((const unsigned char*)&validData, sizeof(validData));
+ } else {
+ file << validData << " ";
+ }
+}
+
+void writeInt(EncodeBase64 &b64, ofstream &file, int data, bool writeBinary) {
+ if (writeBinary) {
+ b64.Write((const unsigned char*)&data, sizeof(data));
+ } else {
+ file << data << " ";
+ }
+}
+
+void writeChar(EncodeBase64 &b64, ofstream &file, unsigned char data, bool writeBinary) {
+ if (writeBinary) {
+ b64.Write((const unsigned char*)&data, sizeof(data));
+ } else {
+ file << data << " ";
+ }
+}
+
+void printUsage() {
+ cerr << "USAGE: projlatlonbin [OPTIONS] infile.nc outfile.vtu" << endl;
+ cerr << "Create lat/lon projection of MPAS-style data from input file infile.nc " << endl;
+ cerr << "and write to XML binary .vtu format for use in ParaView." << endl;
+ cerr << "A series of vtu files are created, one file for each time step." << endl;
+ cerr << "with time prepended to the file name (e.g. 0outfile.vtu)." << endl;
+ cerr << "OPTIONS:" << endl;
+ cerr << " -l intval " << endl;
+ cerr << " vertical level, in range 0 to nVertLevels-1, default is 0 " << endl;
+ cerr << " -m " << endl;
+ cerr << " multilayer view -- default is single layer view" << endl;
+ cerr << " -t floatval " << endl;
+ cerr << " layer_thickness -- use negative value for atmosphere, default is 10" << endl;
+ cerr << " -c floatval " << endl;
+ cerr << " display image with center_longitude in degrees, default is 180" << endl;
+ cerr << " -z " << endl;
+ cerr << " input grid has center longitude at 0, the default is 180 degree or PI radians" << endl;
+ cerr << " -a " << endl;
+ cerr << " atmosphere data (longitude data is in range -PI to PI) " << endl;
+ cerr << " The -a flag also converts the multilayer view to have the " << endl;
+ cerr << " vertical levels stack in the positive Z direction, rather " << endl;
+ cerr << " the negative Z direction, which is the default." << endl;
+ cerr << " -s " << endl;
+ cerr << " timestep to start on (defaults to 0) " << endl;
+ cerr << " -f " << endl;
+ cerr << " timestep to finish on (defaults to last timestep available) " << endl;
+ cerr << " -g " << endl;
+ cerr << " include topography (assumes you have maxLevelCell variable) " << endl;
+ cerr << " -b " << endl;
+ cerr << " do bugfix to ensure no connections to faraway points " << endl;
+ cerr << "Variables with time and vertical level are written out." << endl;
+ cerr << "Tracer vars are named tracer1, tracer2, etc." << endl;
+ cerr << "[vtk datafile Version 2.0, projlatlonbin Version 1.3]" << endl;
+ exit(1);
+}
+
+void getArgs(int argc, char* argv[], int &outputVertLevel,
+ float &layerThickness, double &centerLon, bool &isMultilayer,
+ bool &isAtmosphere, int &startTimeStep, int &finishTimeStep,
+ bool &isZeroCentered, bool &includeTopography, bool &doBugfix,
+ char* &inputFileName, char* &outputFileName) {
+
+ int c;
+ bool selectedVar = false;
+ optind = 1;
+
+ bool lflag = false;
+ bool tflag = false;
+ bool cflag = false;
+ bool mflag = false;
+ bool aflag = false;
+ bool sflag = false;
+ bool fflag = false;
+ bool zflag = false;
+ bool gflag = false;
+ bool bflag = false;
+
+ // -l vertical level (if not specified, assume multilayer)
+ // -t layer thickness (defaults to 1/10 * radius/#layers)
+ // -c center longitude in degrees for latlon projection
+ while ((c=getopt(argc, argv, "l:t:c:mazgbs:f:")) != -1) {
+ switch(c) {
+ case 'l':
+ lflag = true;
+ if (tflag) {
+ cerr << "Can't specify -l and -t at same time." << endl;
+ exit(1);
+ }
+ outputVertLevel = atoi(optarg);
+ isMultilayer = false;
+ break;
+ case 't':
+ tflag = true;
+ if (lflag) {
+ cerr << "Can't specify -l and -t at same time." << endl;
+ exit(1);
+ }
+ layerThickness = atof(optarg);
+ isMultilayer = true;
+ break;
+ case 'c':
+ cflag = true;
+ centerLon = (double)atof(optarg);
+ break;
+ case 'm':
+ mflag = true;
+ if (lflag) {
+ cerr << "Can't specify -l and -m at same time." << endl;
+ exit(1);
+ }
+ isMultilayer = true;
+ break;
+ case 'a':
+ aflag = true;
+ isAtmosphere = true;
+ break;
+ case 's':
+ sflag = true;
+ startTimeStep = atoi(optarg);
+ break;
+ case 'f':
+ sflag = true;
+ finishTimeStep = atoi(optarg);
+ break;
+ case 'z':
+ zflag = true;
+ isZeroCentered = true;
+ break;
+ case 'g':
+ gflag = true;
+ includeTopography = true;
+ break;
+ case 'b':
+ bflag = true;
+ doBugfix = true;
+ break;
+ default:
+ cerr << "Unrecognized option: -" << c << endl;
+ printUsage();
+ exit(1);
+ break;
+ }
+ //argc -= optind;
+ //argv += optind;
+ }
+
+ inputFileName = argv[optind++];
+ outputFileName = argv[optind];
+
+}
+
+// Check if there is a variable by that name
+
+bool isNcVar(NcFile &ncFile, NcToken name) {
+ int num_vars = ncFile.num_vars();
+ for (int i = 0; i < num_vars; i++) {
+ NcVar* ncVar = ncFile.get_var(i);
+ if ((strcmp(ncVar->name(), name)) == 0) {
+ // we have a match, so return
+ return true;
+ }
+ }
+ return false;
+}
+
+// Check if there is a dimension by that name
+
+bool isNcDim(NcFile &ncFile, NcToken name) {
+ int num_dims = ncFile.num_dims();
+ //cerr << "looking for: " << name << endl;
+ for (int i = 0; i < num_dims; i++) {
+ NcDim* ncDim = ncFile.get_dim(i);
+ //cerr << "checking " << ncDim->name() << endl;
+ if ((strcmp(ncDim->name(), name)) == 0) {
+ // we have a match, so return
+ return true;
+ }
+ }
+ return false;
+}
+
+
+int main(int argc, char* argv[])
+{
+
+ bool multilayer = false;
+ bool atmosphere = false;
+ bool includeTopography = false;
+ bool doBugfix = false;
+ int outputVertLevel = 0;
+ char* inputFileName = NULL;
+ char* outputFileName = NULL;
+ float layerThickness = DEFAULT_LAYER_THICKNESS;
+ double centerLon = 180;
+ int startTimeStep = 0;
+ int finishTimeStep = -1;
+ bool isZeroCentered = false;
+
+ if (argc < 3) {
+ printUsage();
+ exit(1);
+ }
+
+ getArgs(argc, argv, outputVertLevel, layerThickness, centerLon, multilayer,
+ atmosphere, startTimeStep, finishTimeStep, isZeroCentered,
+ includeTopography, doBugfix, inputFileName, outputFileName);
+
+ NcFile ncFile(inputFileName);
+
+ if (!ncFile.is_valid())
+ {
+ cerr << "Couldn't open file: " << inputFileName << endl;
+ exit(1);
+ }
+
+ CHECK_DIM(ncFile, "nCells");
+ NcDim* nCells = ncFile.get_dim("nCells");
+
+ CHECK_DIM(ncFile, "nVertices");
+ NcDim* nVertices = ncFile.get_dim("nVertices");
+
+ CHECK_DIM(ncFile, "vertexDegree");
+ NcDim* vertexDegree = ncFile.get_dim("vertexDegree");
+
+ CHECK_DIM(ncFile, "Time");
+ NcDim* Time = ncFile.get_dim("Time");
+
+ CHECK_DIM(ncFile, "nVertLevels");
+ NcDim* nVertLevels = ncFile.get_dim("nVertLevels");
+
+ int maxNVertLevels = nVertLevels->size();
+
+ if ((vertexDegree->size() != 3) && (vertexDegree->size() != 4)) {
+ cerr << "This code is only for hexagonal or quad grid projection" << endl;
+ exit(1);
+ }
+
+ // what to do about nVertLevelsP1?
+ //NcDim* nVertLevelsP1 = ncFile.get_dim("nVertLevelsP1");
+
+ bool singlelayer = !multilayer;
+
+ // double-check we can do multilayer
+ if ((multilayer) && (maxNVertLevels == 1)) {
+ multilayer = false;
+ singlelayer = !multilayer;
+ }
+
+ if (singlelayer && (outputVertLevel > (maxNVertLevels-1))) {
+ cerr << "There is no data for level " << outputVertLevel << "." << endl;
+ exit(1);
+ }
+
+ if (singlelayer && (outputVertLevel < 0)) {
+ cerr << "Level number must be 0 or greater." << endl;
+ exit(1);
+ }
+
+ if ((centerLon < 0) || (centerLon > 360)) {
+ cerr << "Center longitude must be between 0 and 360." << endl;
+ exit(1);
+ }
+
+ if (singlelayer) {
+ maxNVertLevels = 1;
+ }
+
+ if (atmosphere) {
+ layerThickness = -layerThickness;
+ }
+
+ int availTimeSteps = Time->size();
+
+ if (startTimeStep > (availTimeSteps-1)) {
+ cerr << "Start timestep is too large, last timestep is: " << availTimeSteps-1 << endl;
+ exit(1);
+ }
+
+ if (startTimeStep < 0) {
+ cerr << "Start timestep is too small, must be at least 0." << endl;
+ exit(1);
+ }
+
+ if (finishTimeStep > (availTimeSteps-1)) {
+ cerr << "Finish timestep is too large, last timestep is: " << availTimeSteps-1 << endl;
+ exit(1);
+ }
+
+ if (finishTimeStep < 0) {
+ finishTimeStep = availTimeSteps-1;
+ }
+
+ cerr << "centerLon: " << centerLon << " singleLayer: " << singlelayer << " outputVertLevel: " << outputVertLevel << " start: " << startTimeStep << " finish: " << finishTimeStep << " atmosphere: " << atmosphere << " includeTopography: " << includeTopography << " doBugfix: " << doBugfix << endl;
+
+ bool writeBinary = true;
+
+ // figure out what variables to visualize
+ NcVar* dualCellVars[MAX_VARS];
+ NcVar* dualPointVars[MAX_VARS];
+ NcVar* vectorVars[3];
+ int dualCellVarIndex = -1;
+ int dualPointVarIndex = -1;
+ int numDualCells = nVertices->size();
+ int numDualPoints = nCells->size()+1;
+
+ int numVars = ncFile.num_vars();
+
+ bool tracersExist = false;
+ bool isUVector = false;
+ bool isVVector = false;
+
+ for (int i = 0; i < numVars; i++) {
+ NcVar* aVar = ncFile.get_var(i);
+
+ // must have 3 dims
+ // (Time, nCells | nVertices, nVertLevels | nVertLevelsP1)
+
+ int numDims = aVar->num_dims();
+ //cout << "Num Dims of var: " << aVar->name() << " is " << numDims << endl;
+ if ((numDims != 3) && (strcmp(aVar->name(), "tracers"))) {
+ continue; // try the next var
+ } else {
+ // TODO, check if it is a double
+ // assume a double for now
+
+ // check for Time dim 0
+ NcToken dim0Name = aVar->get_dim(0)->name();
+ if (strcmp(dim0Name, "Time"))
+ continue;
+
+ // check for dim 1 being num vertices or cells
+ bool isVertexData = false;
+ bool isCellData = false;
+ NcToken dim1Name = aVar->get_dim(1)->name();
+ if (!strcmp(dim1Name, "nVertices"))
+ isVertexData = true;
+ else if (!strcmp(dim1Name, "nCells"))
+ isCellData = true;
+ else continue;
+
+ // check if dim 2 is nVertLevels or nVertLevelsP1, too
+ NcToken dim2Name = aVar->get_dim(2)->name();
+ if ((strcmp(dim2Name, "nVertLevels"))
+ && (strcmp(dim2Name, "nVertLevelsP1"))) {
+ continue;
+ }
+
+ // Add to cell or point var array
+ if (isVertexData) { // means it is dual cell data
+ dualCellVarIndex++;
+ if (dualCellVarIndex > MAX_VARS-1) {
+ cerr << "Exceeded number of cell vars." << endl;
+ exit(1);
+ }
+ dualCellVars[dualCellVarIndex] = aVar;
+ //cout << "Adding var " << aVar->name() << " to dualCellVars" << endl;
+ } else if (isCellData) { // means it is dual vertex data
+ // check for u/v/w vectors
+
+ if (!strncmp(aVar->name(), "uReconstruct", 12)) {
+ if (!strcmp(aVar->name(), "uReconstructZonal")) {
+ vectorVars[0] = aVar;
+ isUVector = true;
+ } else if (!strcmp(aVar->name(), "uReconstructMeridional")) {
+ vectorVars[1] = aVar;
+ isVVector = true;
+ }
+ } else if (strcmp(aVar->name(), "tracers")) {
+ dualPointVarIndex++;
+ if (dualPointVarIndex > MAX_VARS-1) {
+ cerr << "Exceeded number of point vars." << endl;
+ exit(1);
+ }
+ dualPointVars[dualPointVarIndex] = aVar;
+ //cout << "Adding var " << aVar->name() << " to dualPointVars" << endl;
+ } else { // case of tracers, add each as "tracer0", "tracer1", etc.
+ tracersExist = true;
+ int numTracers = aVar->get_dim(3)->size();
+ for (int t = 0; t < numTracers; t++) {
+ dualPointVarIndex++;
+ if (dualPointVarIndex > MAX_VARS-1) {
+ cerr << "Exceeded number of point vars." << endl;
+ exit(1);
+ }
+ dualPointVars[dualPointVarIndex] = aVar;
+ //cout << "Adding var " << aVar->name() << " to dualPointVars" << endl;
+ }
+ }
+ }
+ }
+ }
+
+ // TODO
+ // prompt the user to find out which fields are of interest?
+ // for now, assume all are of interest
+
+ // get points (centers of primal-mesh cells)
+
+ const float BLOATFACTOR = .5;
+ int myCellSize = (int)floor(nCells->size()*(1.0 + BLOATFACTOR));
+ CHECK_VAR(ncFile, "lonCell");
+ double *xCellData = (double*)malloc(myCellSize * sizeof(double));
+ CHECK_MALLOC(xCellData);
+ NcVar* xCellVar = ncFile.get_var("lonCell");
+ xCellVar->get(xCellData+1, nCells->size());
+ // point 0 is 0.0
+ *xCellData = 0.0;
+
+ const double PI = 3.141592;
+
+ // if atmospheric data, or zero centered, set center to 180 instead of 0
+ if (atmosphere || isZeroCentered) {
+ for (int j = 1; j <= nCells->size(); j++) {
+ // need to shift over the point so center is at PI
+ if (xCellData[j] < 0) {
+ xCellData[j] += 2*PI;
+ }
+ }
+ }
+
+ double centerRad = centerLon * PI / 180.0;
+
+ if (centerLon != 180) {
+ for (int j = 1; j <= nCells->size(); j++) {
+ // need to shift over the point if centerLon dictates
+ if (centerRad < PI) {
+ if (xCellData[j] > (centerRad + PI)) {
+ xCellData[j] = -((2*PI) - xCellData[j]);
+ }
+ } else if (centerRad > PI) {
+ if (xCellData[j] < (centerRad - PI)) {
+ xCellData[j] += 2*PI;
+ }
+ }
+ }
+ }
+
+ CHECK_VAR(ncFile, "latCell");
+ double *yCellData = (double*)malloc(myCellSize * sizeof(double));
+ CHECK_MALLOC(yCellData);
+ NcVar* yCellVar = ncFile.get_var("latCell");
+ yCellVar->get(yCellData+1, nCells->size());
+ // point 0 is 0.0
+ *yCellData = 0.0;
+
+ // get dual-mesh cells
+
+ CHECK_VAR(ncFile, "cellsOnVertex");
+ int *cellsOnVertex = (int *) malloc((nVertices->size()) * vertexDegree->size() *
+ sizeof(int));
+ CHECK_MALLOC(cellsOnVertex);
+ int myCOVSize = (int)floor(nVertices->size()*(1.0 + BLOATFACTOR))+1;
+ int *myCellsOnVertex = (int *) malloc(myCOVSize * vertexDegree->size() * sizeof(int));
+ CHECK_MALLOC(myCellsOnVertex);
+ NcVar *cellsOnVertexVar = ncFile.get_var("cellsOnVertex");
+ //cout << "getting cellsOnVertexVar</font>
<font color="gray">";
+ cellsOnVertexVar->get(cellsOnVertex, nVertices->size(), vertexDegree->size());
+
+ // allocate an array to map the extra points and cells to the original
+ // so that when obtaining data, we know where to get it
+ int *pointMap = (int*)malloc((int)floor(nCells->size()*BLOATFACTOR) * sizeof(int));
+ CHECK_MALLOC(pointMap);
+ int *cellMap = (int*)malloc((int)floor(nVertices->size()*BLOATFACTOR) * sizeof(int));
+ CHECK_MALLOC(cellMap);
+
+
+ int *myptr = myCellsOnVertex;
+ int *myExtraX;
+ int *myExtraY;
+ int currentExtraPoint = numDualPoints;
+ int currentExtraCell = numDualCells;
+
+ // For each cell, examine vertices
+ // Add new points and cells where needed to account for wraparound.
+
+ int acell = 775;
+ int mirrorcell;
+ int apoint;
+ int mirrorpoint;
+
+
+ for (int j = 0; j < numDualCells; j++ ) {
+ int *dualCells = cellsOnVertex + (j * vertexDegree->size());
+
+ // go through and make sure none of the referenced points are
+ // out of range (<=0 or >nCells->size())
+ // if so, set all to point 0
+ for (int k = 0; k < vertexDegree->size(); k++) {
+ if ((dualCells[k] <= 0) || (dualCells[k] > nCells->size())) {
+ for (int m = 0; m < vertexDegree->size(); m++) {
+ dualCells[m] = 0;
+ }
+ break;
+ }
+ }
+
+ int lastk;
+
+ if (doBugfix) {
+ /////BUG FIX for problem where cells are stretching to a faraway point
+ lastk = vertexDegree->size()-1;
+ //const double thresh = .139626340159546; // 8 degrees
+ const double thresh = .06981317007977; // 8 degrees
+ for (int k = 0; k < vertexDegree->size(); k++) {
+ double xdiff = abs(xCellData[dualCells[k]] - xCellData[dualCells[lastk]]);
+ double ydiff = abs(yCellData[dualCells[k]] - yCellData[dualCells[lastk]]);
+ // Don't look at cells at map border
+ //if ((xdiff <= 5.5) && (xdiff > (2*.06981317)))
+ //if (((xdiff <= 6) && (xdiff > 3)) || (ydiff > .06981317))
+ if (((xdiff > thresh) && (ydiff > thresh))
+ // || ((xdiff <= 6) && (xdiff > 3))
+ || (ydiff > thresh)
+ ) {
+ for (int m = 0; m < vertexDegree->size(); m++) {
+ dualCells[m] = 0;
+ }
+ break;
+ }
+ }
+ }
+
+ lastk = vertexDegree->size()-1;
+ bool xWrap = false;
+ bool yWrap = false;
+ for (int k = 0; k < vertexDegree->size(); k++) {
+ if (abs(xCellData[dualCells[k]] - xCellData[dualCells[lastk]]) > 5.5) xWrap = true;
+ //if (abs(yCellData[dualCells[k]] - yCellData[dualCells[lastk]]) > 2) yWrap = true;
+ lastk = k;
+ }
+
+ if (xWrap || yWrap) {
+ //cerr << "Cell wrap: " << j << endl;
+ }
+
+ if (xWrap) {
+ //cerr << "It wrapped in x direction" << endl;
+ double anchorX = xCellData[dualCells[0]];
+ double anchorY = yCellData[dualCells[0]];
+ *myptr = *(dualCells);
+ myptr++;
+
+ if (j == acell) {
+ //cout << "cell " << acell << " anchor x: " << anchorX << " y: " << anchorY << endl;
+ }
+
+ // modify existing cell, so it doesn't wrap
+ // move points to one side
+
+ // first point is anchor it doesn't move
+ for (int k = 1; k < vertexDegree->size(); k++) {
+ double neighX = xCellData[dualCells[k]];
+ double neighY = yCellData[dualCells[k]];
+
+ if (j == acell) {
+ //cout << "cell " << acell << " k: " << k << " x: " << neighX << " y: " << neighY << endl;
+ }
+
+ // add a new point, figure out east or west
+ if (abs(neighX - anchorX) > 5.5) {
+ double neighEastX;
+ double neighWestX;
+
+ // add on east
+ if (neighX < anchorX) {
+ if (j == acell) {
+ //cout << "add on east" << endl;
+ }
+ neighEastX = neighX + (2*PI);
+ xCellData[currentExtraPoint] = neighEastX;
+ yCellData[currentExtraPoint] = neighY;
+ } else {
+ // add on west
+ if (j == acell) {
+ //cout << "add on west" << endl;
+ }
+ neighWestX = neighX - (2*PI);
+ xCellData[currentExtraPoint] = neighWestX;
+ yCellData[currentExtraPoint] = neighY;
+ }
+
+ if (j == acell) {
+ //cout << "x: " << xCellData[currentExtraPoint] << " y: " << yCellData[currentExtraPoint] << endl;
+ }
+
+ if (j == acell) {
+ //cout << "currentExtraPoint: " << currentExtraPoint << endl;
+ }
+
+ // add the new point to list of vertices
+ *myptr = currentExtraPoint;
+
+ // record mapping
+ *(pointMap + (currentExtraPoint - numDualPoints)) = dualCells[k];
+
+ if (j == acell) {
+ mirrorpoint = currentExtraPoint;
+ apoint = dualCells[k];
+ //cout << "mirror point " << mirrorpoint << " has pointMap to: " << dualCells[k] << endl;
+ }
+
+ myptr++;
+ currentExtraPoint++;
+ } else {
+ // use existing kth point
+ if (j == acell) {
+ //cout << "use existing point" << endl;
+ }
+ *myptr = dualCells[k];
+ myptr++;
+ }
+ }
+
+ // add a mirror image cell on other side so there are no
+ // gaps in the map
+
+ // move anchor to other side
+ if (anchorX > centerRad) {
+ anchorX = anchorX - (2*PI);
+ } else {
+ anchorX = anchorX + (2*PI);
+ }
+
+ if (j == acell) {
+ //cout << "add new cell " << currentExtraCell << " to mirror " << acell << " anchorX: " << anchorX << " anchorY: " << anchorY << endl;
+ mirrorcell = currentExtraCell;
+ }
+
+ // move addedCellsPtr to myCellsOnVertex extra cells area
+ int* addedCellsPtr = myCellsOnVertex + (currentExtraCell * vertexDegree->size());
+
+ // add point coord and add to list of cells
+ xCellData[currentExtraPoint] = anchorX;
+ yCellData[currentExtraPoint] = anchorY;
+ *addedCellsPtr = currentExtraPoint;
+
+ // record mapping
+ *(pointMap + (currentExtraPoint - numDualPoints)) = dualCells[0];
+
+ addedCellsPtr++;
+ currentExtraPoint++;
+
+ for (int k = 1; k < vertexDegree->size(); k++) {
+ double neighX = xCellData[dualCells[k]];
+ double neighY = yCellData[dualCells[k]];
+
+ // add a new point, figure out east or west
+ if (abs(neighX - anchorX) > 5.5) {
+ double neighEastX;
+ double neighWestX;
+
+ // add on east
+ if (neighX < anchorX) {
+ neighEastX = neighX + (2*PI);
+ xCellData[currentExtraPoint] = neighEastX;
+ yCellData[currentExtraPoint] = neighY;
+ } else {
+ // add on west
+ neighWestX = neighX - (2*PI);
+ xCellData[currentExtraPoint] = neighWestX;
+ yCellData[currentExtraPoint] = neighY;
+ }
+
+ // add the new point to list of vertices
+ *addedCellsPtr = currentExtraPoint;
+
+ // record mapping
+ *(pointMap + (currentExtraPoint - numDualPoints)) = dualCells[k];
+
+ addedCellsPtr++;
+ currentExtraPoint++;
+ } else {
+ // use existing kth point
+ *addedCellsPtr = dualCells[k];
+ addedCellsPtr++;
+ }
+ }
+ *(cellMap + (currentExtraCell - numDualCells)) = j;
+ currentExtraCell++;
+ }
+ if (yWrap) {
+ //cerr << "It wrapped in y direction" << endl;
+ }
+
+ // if cell doesn't extend past lat/lon perimeter, then add it to myCellsOnVertex
+ if (!xWrap && !yWrap) {
+ for (int k=0; k< vertexDegree->size(); k++) {
+ *myptr = *(dualCells+k);
+ myptr++;
+ }
+ }
+ if (currentExtraCell > myCOVSize) {
+ cerr << "Exceeded storage for extra cells!" << endl;
+ return 1;
+ }
+ if (currentExtraPoint > myCellSize) {
+ cerr << "Exceeded storage for extra points!" << endl;
+ return 1;
+ }
+ }
+
+
+ int varVertLevels = 0;
+
+ // write a file with the geometry.
+
+ ostringstream geoFileName;
+ geoFileName << "geo_" << outputFileName;
+ ofstream geoFile(geoFileName.str().c_str(), ios::out);
+ if (!geoFile)
+ {
+ cerr << "vtk output file could not be opened" <<endl;
+ exit(1);
+ }
+
+ EncodeBase64 b64(&geoFile);
+
+ // write header
+
+ string formatString;
+ if (writeBinary) {
+ formatString = "binary";
+ } else {
+ formatString = "ascii";
+ }
+
+ unsigned int binarySize = 0;
+
+ // write points (the points are the primal-grid cell centers)
+
+
+ // write the point at each vertical level, plus one level for last layer
+
+ geoFile << " <Points>" << endl;
+
+ // see if we can change to double
+
+ geoFile << " <DataArray type=\"Float32\" Name=\"Points\" NumberOfComponents=\"3\" format=\"" << formatString << "\">" << endl;
+
+ geoFile.precision(16);
+
+ float tmp = 0;
+
+ geoFile << " ";
+
+ if (writeBinary) {
+ b64.StartWriting();
+ if (singlelayer) {
+ binarySize = currentExtraPoint*sizeof(float)*3;
+ } else {
+ binarySize = currentExtraPoint*(maxNVertLevels+1)*sizeof(float)*3;
+ }
+ binarySize = padSize(binarySize);
+ b64.Write((const unsigned char*)&binarySize, sizeof(int));
+ }
+
+ //cout << "Writing points at each level" << endl;
+
+ for (int j = 0; j < currentExtraPoint; j++ )
+ {
+ if (abs(xCellData[j]) < 1e-126) xCellData[j] = 0;
+ if (abs(yCellData[j]) < 1e-126) yCellData[j] = 0;
+
+ float xLon = xCellData[j] * 180.0 / PI;
+ float yLat = yCellData[j] * 180.0 / PI;
+
+ if (singlelayer) {
+ if (writeBinary) {
+ tmp = 0;
+ b64.Write((const unsigned char*)&xLon, sizeof(float));
+ b64.Write((const unsigned char*)&yLat, sizeof(float));
+ b64.Write((const unsigned char*)&tmp, sizeof(float));
+ } else {
+ geoFile << " " << xLon << " " << yLat << " " << 0 << endl;
+ }
+ } else {
+ for (int levelNum = 0; levelNum < maxNVertLevels+1; levelNum++) {
+ tmp = -(((float)levelNum)*layerThickness);
+ if (writeBinary) {
+ b64.Write((const unsigned char*)&xLon, sizeof(float));
+ b64.Write((const unsigned char*)&yLat, sizeof(float));
+ b64.Write((const unsigned char*)&tmp, sizeof(float));
+ } else {
+ geoFile << " " << xLon << " " << yLat << " " << tmp << endl;
+ }
+ }
+ }
+ }
+
+ if (writeBinary) {
+ b64.EndWriting();
+ geoFile << endl;
+ }
+
+ geoFile << " </DataArray>" << endl;
+ geoFile << " </Points>" << endl;
+
+ free (xCellData);
+ free (yCellData);
+
+ // Write dual-mesh cells
+ // Dual-mesh cells are triangles with primal-mesh cell
+ // centers as the vertices.
+ // The number of dual-mesh cells is the number of vertices in the
+ // primal mesh.
+
+ geoFile << " <Cells>" << endl;
+ geoFile << " <DataArray type=\"Int32\" Name=\"connectivity\" format=\"" << formatString << "\">" << endl;
+
+ // for each dual-mesh cell, write the points for each
+
+ if (writeBinary) {
+ geoFile << " ";
+ b64.StartWriting();
+ if (singlelayer) {
+ binarySize = currentExtraCell*vertexDegree->size()*sizeof(int);
+ } else {
+ binarySize = currentExtraCell*vertexDegree->size()*2*sizeof(int)
+ *maxNVertLevels;
+ }
+ binarySize = padSize(binarySize);
+ //cout << "BinarySize: " << binarySize << endl;
+ b64.Write((const unsigned char*)&binarySize, sizeof(int));
+ }
+
+ // for each dual-mesh cell, write number of points for each
+ // and then list the points by number
+
+ //cout << "Writing connectivity" << endl;
+
+ unsigned int val;
+
+ int* maxLevelCell;
+
+ if (includeTopography) {
+ CHECK_VAR(ncFile, "maxLevelCell");
+ maxLevelCell = (int*)malloc((nCells->size()+1) * sizeof(int));
+ CHECK_MALLOC(maxLevelCell);
+ NcVar *maxLevelCellVar = ncFile.get_var("maxLevelCell");
+ maxLevelCellVar->get(maxLevelCell+1, nCells->size());
+ }
+
+ for (int j = 0; j < currentExtraCell ; j++) {
+
+ // since primal vertex(pt) numbers == dual cell numbers
+ // we go through the primal vertices, find the cells around
+ // them, and since those primal cell numbers are dual
+ // point numbers,
+ // we can write the cell numbers for the cellsOnVertex
+ // and those will be the numbers of the dual vertices (pts).
+
+ //cout << "j: " << j << endl;
+
+ int* dualCells = myCellsOnVertex + (j * vertexDegree->size());
+
+ int minLevel;
+
+ if (includeTopography) {
+ int* cells;
+ //check if it is a mirror cell, if so, get original
+ if (j >= numDualCells) {
+ //cout << "setting origCell" << endl;
+ int origCell = *(cellMap + (j - numDualCells));
+ //cout << "mirror cell: " <<j<< " origCell: "<< origCell << endl;
+ cells = cellsOnVertex + (origCell*vertexDegree->size());
+ } else cells = cellsOnVertex + (j * vertexDegree->size());
+
+ //cout << "calc minLevel" << endl;
+ minLevel = maxLevelCell[cells[0]];
+ //cout << "initial minLevel:" << minLevel << endl;
+ // Take the min of the maxLevelCell of each primal cell (dual point).
+ for (int k = 1; k < vertexDegree->size(); k++) {
+ //cout << " " << maxLevelCell[cells[k]];
+ minLevel = min(minLevel, maxLevelCell[cells[k]]);
+ }
+ //cout << endl;
+ }
+
+ if (singlelayer) {
+ if (!writeBinary) geoFile << " ";
+
+ // If that min is greater than or equal to this output level,
+ // include the dual cell, otherwise set all points to zero.
+
+ if (includeTopography && ((minLevel-1) < outputVertLevel)) {
+ //cerr << "Setting all points to zero" << endl;
+ val = 0;
+ for (int k = 0; k < vertexDegree->size(); k++) {
+ writeInt(b64, geoFile, val, writeBinary);
+ }
+ } else {
+ for (int k = 0; k < vertexDegree->size(); k++) {
+ writeInt(b64, geoFile, dualCells[k], writeBinary);
+ }
+ }
+
+ if (!writeBinary) geoFile << endl;
+ } else {
+
+ // for each level, write the cell
+ for (int levelNum = 0; levelNum < maxNVertLevels; levelNum++) {
+ if(!writeBinary) geoFile << " ";
+
+ if (includeTopography && ((minLevel-1) < levelNum)) {
+ // setting all points to zero
+ val = 0;
+ for (int m = 0; m < vertexDegree->size()*2; m++) {
+ writeInt(b64, geoFile, val, writeBinary);
+ }
+
+ } else {
+
+ for (int k = 0; k < vertexDegree->size(); k++)
+ {
+ val = (dualCells[k]*(maxNVertLevels+1)) + levelNum;
+ writeInt(b64, geoFile, val, writeBinary);
+ }
+
+ for (int k = 0; k < vertexDegree->size(); k++)
+ {
+ val = (dualCells[k]*(maxNVertLevels+1)) + levelNum + 1;
+ writeInt(b64, geoFile, val, writeBinary);
+ }
+ }
+ if (!writeBinary) geoFile << endl;
+ }
+ }
+ }
+
+ if (writeBinary) {
+ b64.EndWriting();
+ geoFile << endl;
+ }
+
+ geoFile << " </DataArray>" << endl;
+
+ if (includeTopography) {
+ free (maxLevelCell);
+ }
+ free(cellsOnVertex);
+ free(myCellsOnVertex);
+
+
+ // Specify the offset into connectivity array for the end of each cell
+ geoFile << " <DataArray type=\"Int32\" Name=\"offsets\" format=\"" << formatString << "\">";
+
+ int maxCells;
+ int maxPoints;
+
+ if (singlelayer) {
+ maxCells = currentExtraCell;
+ } else {
+ maxCells = currentExtraCell*maxNVertLevels;
+ }
+
+ if (writeBinary) {
+ geoFile << endl << " ";
+ b64.StartWriting();
+ binarySize = maxCells*sizeof(int);
+ binarySize = padSize(binarySize);
+ //cout << "BinarySize: " << binarySize << endl;
+ b64.Write((const unsigned char*)&binarySize, sizeof(int));
+ }
+
+ int counter = 0;
+ unsigned int offset = 0;
+
+ for (int j = 0; j < maxCells ; j++) {
+ if ((!writeBinary) && (counter%6 == 0)) geoFile << endl << " ";
+ if (singlelayer) {
+ offset += vertexDegree->size();
+ } else {
+ offset += 2*vertexDegree->size();
+ }
+ writeInt(b64, geoFile, offset, writeBinary);
+ counter++;
+ }
+
+ if (writeBinary) {
+ b64.EndWriting();
+ }
+
+ geoFile << endl;
+ geoFile << " </DataArray>" << endl;
+
+ // write cell types
+ unsigned char cellType;
+ switch (vertexDegree->size()) {
+ case 3:
+ if (singlelayer) {
+ cellType = VTK_TRIANGLE;
+ } else {
+ cellType = VTK_WEDGE;
+ }
+ break;
+ case 4:
+ if (singlelayer) {
+ cellType = VTK_QUAD;
+ } else {
+ cellType = VTK_HEXAHEDRON;
+ }
+ break;
+ default:
+ break;
+ }
+
+ geoFile << " <DataArray type=\"UInt8\" Name=\"types\" format=\"" << formatString << "\">";
+
+ if (writeBinary) {
+ geoFile << endl << " ";
+ b64.StartWriting();
+ binarySize = maxCells*sizeof(char);
+ binarySize = padSize(binarySize);
+ //cout << "BinarySize: " << binarySize << endl;
+ b64.Write((const unsigned char*)&binarySize, sizeof(int));
+ }
+
+ for (int j = 0; j < maxCells ; j++) {
+ if ((!writeBinary) && (j%6 == 0)) {
+ geoFile << endl << " ";
+ }
+ writeChar(b64, geoFile, cellType, writeBinary);
+ }
+
+ if (writeBinary) b64.EndWriting();
+
+ geoFile << endl;
+ geoFile << " </DataArray>" << endl;
+ geoFile << " </Cells>" << endl;
+ geoFile << " </Piece>" << endl;
+ geoFile << " </UnstructuredGrid>" << endl;
+ geoFile << "</VTKFile>" << endl;
+
+ geoFile.close();
+
+ // For each timestep, write data for each level
+
+
+ // for each timestep, copy the geometry file, since we don't want to
+ // have to recompute the points
+ for (int i = startTimeStep; i <= finishTimeStep; i++) {
+
+ ostringstream vtkFileName;
+ vtkFileName << i << outputFileName;
+ ofstream vtkFile(vtkFileName.str().c_str(), ios::out);
+ if (!vtkFile)
+ {
+ cerr << "vtk output file could not be opened" <<endl;
+ exit(1);
+ }
+
+ vtkFile.precision(16);
+ vtkFile << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">" << endl;
+ vtkFile << " <UnstructuredGrid>" << endl;
+
+ if (singlelayer) {
+ maxPoints = currentExtraPoint;
+ } else {
+ maxPoints = currentExtraPoint*(maxNVertLevels+1);
+ }
+
+ vtkFile << " <Piece NumberOfPoints=\"" << maxPoints << "\" NumberOfCells=\"" << maxCells << "\">" << endl;
+
+ int vertLevelIndex = 0;
+ if (singlelayer) vertLevelIndex = outputVertLevel;
+
+ // If by point, write out point data
+ vtkFile << " <PointData>" << endl;
+
+ int printstep = -1;
+
+ if (i == printstep) cout << "TIME STEP: " << i << endl;
+
+ int tracerNum = 0;
+
+ double* dualPointVarData;
+ dualPointVarData = (double*)malloc((sizeof(double)) * currentExtraPoint * maxNVertLevels);
+ CHECK_MALLOC(dualPointVarData);
+
+ EncodeBase64 vtkb64(&vtkFile);
+
+ for (int v = 0; v <= dualPointVarIndex; v++) {
+
+ // Read variable number v data for that timestep
+
+ varVertLevels = dualPointVars[v]->get_dim(2)->size();
+
+ bool isTracer = false;
+
+ if (!strcmp(dualPointVars[v]->name(), "tracers")) {
+ isTracer = true;
+ // Uncomment if want to exclude tracers.
+ // continue;
+ }
+
+ // Write variable number v data for that timestep
+
+ if (isTracer) {
+ dualPointVars[v]->set_cur(i, 0, vertLevelIndex, tracerNum);
+ if (singlelayer) {
+ dualPointVars[v]->get(dualPointVarData+1, 1, nCells->size(), 1, 1);
+ } else {
+ dualPointVars[v]->get(dualPointVarData+maxNVertLevels, 1, nCells->size(), maxNVertLevels, 1);
+ }
+ } else {
+ dualPointVars[v]->set_cur(i, 0, vertLevelIndex);
+ if (singlelayer) {
+ dualPointVars[v]->get(dualPointVarData+1, 1, nCells->size(), 1);
+ } else {
+ dualPointVars[v]->get(dualPointVarData+maxNVertLevels, 1, nCells->size(), maxNVertLevels);
+ }
+ }
+
+ vtkFile << " <DataArray type=\"Float32\" Name=\"" << dualPointVars[v]->name();
+ if(isTracer) vtkFile << tracerNum+1;
+ vtkFile << "\" format=\"" << formatString << "\"> " << endl;
+
+ // Write variable number v data for that timestep
+
+ vtkFile << " ";
+
+ float defaultPointVal = 0.0;
+
+ double *var_target;
+ float validData;
+
+ if (writeBinary) {
+ vtkb64.StartWriting();
+ binarySize = maxPoints*sizeof(float);
+ binarySize = padSize(binarySize);
+ //cout << "BinarySize: " << binarySize << endl;
+ vtkb64.Write((const unsigned char*)&binarySize, sizeof(int));
+ }
+
+ counter = 0;
+
+ //write dummy
+ // 0th dualPointVarData value is a filler, zero, which is not necessarily in
+ // the range of the other data, so we look at the 1st value instead for our dummy
+ if (singlelayer) {
+ var_target = dualPointVarData + 1;
+ writeFloat(vtkb64, vtkFile, *var_target, writeBinary);
+ } else {
+ var_target = dualPointVarData + maxNVertLevels;
+ for (int levelNum = 0; levelNum < maxNVertLevels; levelNum++) {
+ writeFloat(vtkb64, vtkFile, *var_target, writeBinary);
+ var_target++;
+ }
+
+ // write highest level dummy point (duplicate of last level)
+ var_target--;
+ writeFloat(vtkb64, vtkFile, *var_target, writeBinary);
+
+ // Now set the var_target to the 1st data value
+ var_target = dualPointVarData + maxNVertLevels;
+ }
+
+ for (int j = 1; j < numDualPoints; j++) {
+
+ if (singlelayer) {
+ writeFloat(vtkb64, vtkFile, *var_target, writeBinary);
+ var_target++;
+ } else {
+
+ // write data for one point -- lowest level to highest
+ for (int levelNum = 0; levelNum < maxNVertLevels; levelNum++) {
+ writeFloat(vtkb64, vtkFile, *var_target, writeBinary);
+ var_target++;
+ }
+
+ // for last layer of dual points, repeat last level's values
+ // Need Mark's input on this one
+ var_target--;
+ writeFloat(vtkb64, vtkFile, *var_target, writeBinary);
+ var_target++;
+ }
+ }
+
+ // put out data for extra points
+ for (int j = numDualPoints; j < currentExtraPoint; j++) {
+ if (j == mirrorpoint) {
+ //cout << "data for mirror point " << mirrorpoint << " from point " << *(pointMap + j - numDualPoints) << endl;
+ }
+ // use map to find out what point data we are using
+ var_target = dualPointVarData + ((*(pointMap + j - numDualPoints))*maxNVertLevels);
+
+ if (singlelayer) {
+ writeFloat(vtkb64, vtkFile, *var_target, writeBinary);
+ validData = convertDouble2ValidFloat (*var_target);
+ var_target++;
+ } else {
+ // write data for one point -- lowest level to highest
+ for (int levelNum = 0; levelNum < maxNVertLevels; levelNum++) {
+ writeFloat(vtkb64, vtkFile, *var_target, writeBinary);
+ var_target++;
+ }
+
+ // for last layer of dual points, repeat last level's values
+ // Need Mark's input on this one
+ var_target--;
+ writeFloat(vtkb64, vtkFile, *var_target, writeBinary);
+ }
+ }
+ if (writeBinary) vtkb64.EndWriting();
+ vtkFile << endl << " </DataArray>" << endl;
+
+ if (isTracer) tracerNum++;
+ }
+
+ free (dualPointVarData);
+
+ float validData;
+ const int maxDim = 2; //z is always 0 for vectors in lat/lon projection
+
+ if (isUVector) {
+ vtkFile << " <DataArray type=\"Float32\" Name=\"uReconstruct\" NumberOfComponents=\"3\" format=\"" << formatString << "\"> " << endl;
+
+ double *vectorData[maxDim] = {NULL, NULL};
+ for (int d = 0; d < maxDim; d++) {
+ vectorVars[d]->set_cur(i, 0, vertLevelIndex);
+ if (singlelayer) {
+ vectorData[d] = (double*)malloc((sizeof(double)) * (nCells->size() + 1));
+ vectorVars[d]->get(vectorData[d]+1, 1, nCells->size(), 1);
+ } else {
+ vectorData[d] = (double*)malloc((sizeof(double)) * (nCells->size() + 1) * maxNVertLevels);
+ vectorVars[d]->get(vectorData[d]+maxNVertLevels, 1, nCells->size(), maxNVertLevels);
+ }
+ CHECK_MALLOC(vectorData[d]);
+ }
+
+ vtkFile << " ";
+ if (writeBinary) {
+ vtkb64.StartWriting();
+ binarySize = maxPoints*sizeof(float)*3;
+ binarySize = padSize(binarySize);
+ //cout << "BinarySize: " << binarySize << endl;
+ vtkb64.Write((const unsigned char*)&binarySize, sizeof(int));
+ }
+
+ double* var_target;
+
+ // dummy first
+ if (singlelayer) {
+ for (int d= 0; d < maxDim; d++) {
+ var_target = vectorData[d] + 1;
+ writeFloat(vtkb64, vtkFile, *var_target, writeBinary);
+ }
+ //third dim is 0
+ writeFloat(vtkb64, vtkFile, 0.0, writeBinary);
+ if (!writeBinary) vtkFile << endl;
+ } else {
+ for (int levelNum = 0; levelNum < maxNVertLevels; levelNum++) {
+ for (int d= 0; d < maxDim; d++) {
+ var_target = vectorData[d]+maxNVertLevels+levelNum;
+ writeFloat(vtkb64, vtkFile, *var_target, writeBinary);
+ }
+ //third dim is 0
+ writeFloat(vtkb64, vtkFile, 0.0, writeBinary);
+ if (!writeBinary) vtkFile << endl;
+ }
+ }
+ // write highest level dummy point
+ if (multilayer) {
+ for (int d= 0; d < maxDim; d++) {
+ var_target = vectorData[d]+maxNVertLevels+(maxNVertLevels-1);
+ writeFloat(vtkb64, vtkFile, *var_target, writeBinary);
+ }
+ //third dim is 0
+ writeFloat(vtkb64, vtkFile, 0.0, writeBinary);
+ if (!writeBinary) vtkFile << endl;
+ }
+
+ // rest of data
+ for (int j = 1; j < numDualPoints; j++)
+ {
+ if (singlelayer) {
+ for (int d= 0; d < maxDim; d++) {
+ writeFloat(vtkb64, vtkFile, *(vectorData[d] + j), writeBinary);
+ }
+ //third dim is 0
+ writeFloat(vtkb64, vtkFile, 0.0, writeBinary);
+ if (!writeBinary) vtkFile << endl;
+ } else {
+ for (int levelNum = 0; levelNum < maxNVertLevels; levelNum++) {
+ for (int d= 0; d < maxDim; d++) {
+ writeFloat(vtkb64, vtkFile, *(vectorData[d] + (j*maxNVertLevels) + levelNum), writeBinary);
+ }
+ //third dim is 0
+ writeFloat(vtkb64, vtkFile, 0.0, writeBinary);
+ if (!writeBinary) vtkFile << endl;
+ }
+ // for last layer of dual points, repeat last level's values
+ // Need Mark's input on this one
+ for (int d= 0; d < maxDim; d++) {
+ writeFloat(vtkb64, vtkFile, *(vectorData[d]+(j*maxNVertLevels)+maxNVertLevels-1), writeBinary);
+ }
+ //third dim is 0
+ writeFloat(vtkb64, vtkFile, 0.0, writeBinary);
+ if (!writeBinary) vtkFile << endl;
+ }
+ }
+
+ // put out data for extra points
+ for (int j = numDualPoints; j < currentExtraPoint; j++) {
+
+ // use map to find out what point data we are using
+ int var_target = ((*(pointMap + j - numDualPoints))*maxNVertLevels);
+
+ if (singlelayer) {
+ for (int d= 0; d < maxDim; d++) {
+ writeFloat(vtkb64, vtkFile, *(vectorData[d]+var_target), writeBinary);
+ }
+ //third dim is 0
+ validData = 0;
+ writeFloat(vtkb64, vtkFile, 0.0, writeBinary);
+ if (!writeBinary) vtkFile << endl;
+ } else {
+ // write data for one point -- lowest level to highest
+ for (int levelNum = 0; levelNum < maxNVertLevels; levelNum++) {
+ for (int d= 0; d < maxDim; d++) {
+ writeFloat(vtkb64, vtkFile, *(vectorData[d]+var_target+levelNum), writeBinary);
+ }
+ //third dim is 0
+ writeFloat(vtkb64, vtkFile, 0.0, writeBinary);
+ if (!writeBinary) vtkFile << endl;
+ }
+
+ // for last layer of dual points, repeat last level's values
+ for (int d= 0; d < maxDim; d++) {
+ writeFloat(vtkb64, vtkFile, *(vectorData[d]+var_target+maxNVertLevels-1), writeBinary);
+ }
+ //third dim is 0
+ writeFloat(vtkb64, vtkFile, 0.0, writeBinary);
+ if (!writeBinary) vtkFile << endl;
+ }
+ }
+ if (writeBinary) {
+ vtkb64.EndWriting();
+ vtkFile << endl;
+ }
+ vtkFile << " </DataArray>" << endl;
+ for (int d = 0; d < maxDim; d++) {
+ free(vectorData[d]);
+ }
+ }
+
+ vtkFile << " </PointData>" << endl;
+
+ // if by cell, then write out cell data
+ double* dualCellVarData;
+ dualCellVarData = (double*)malloc((sizeof(double)) * currentExtraCell * maxNVertLevels);
+ CHECK_MALLOC(dualCellVarData);
+
+ vtkFile << " <CellData>" << endl;
+
+
+ for (int v = 0; v <= dualCellVarIndex; v++) {
+
+ // Write variable number v data for that timestep
+ vtkFile << " <DataArray type=\"Float32\" Name=\"" << dualCellVars[v]->name() << "\" format=\"" << formatString << "\"> " ;
+
+ // Read variable number v data for that timestep and level
+ dualCellVars[v]->set_cur(i, 0, vertLevelIndex);
+ if (singlelayer) {
+ dualCellVars[v]->get(dualCellVarData, 1, numDualCells, 1);
+ } else {
+ dualCellVars[v]->get(dualCellVarData, 1, numDualCells, maxNVertLevels);
+ }
+
+ if (writeBinary) {
+ vtkFile << endl << " ";
+ vtkb64.StartWriting();
+ binarySize = maxCells*sizeof(float);
+ binarySize = padSize(binarySize);
+ //cout << "BinarySize: " << binarySize << endl;
+ vtkb64.Write((const unsigned char*)&binarySize, sizeof(int));
+ }
+
+ double *var_target = dualCellVarData;
+ counter = 0;
+
+ for (int j = 0; j < numDualCells; j++) {
+
+ if (singlelayer) {
+ if ((!writeBinary) && (j%6 == 0)) vtkFile << endl << " ";
+ writeFloat(vtkb64, vtkFile, *var_target, writeBinary);
+ var_target++;
+ } else {
+ for (int levelNum = 0; levelNum < maxNVertLevels; levelNum++)
+ {
+ if ((!writeBinary) && (counter%6 == 0)) vtkFile << endl << " ";
+ writeFloat(vtkb64, vtkFile, *var_target, writeBinary);
+ var_target++;
+ counter++;
+ }
+ }
+ }
+
+ counter = 0;
+ for (int j = numDualCells; j < currentExtraCell; j++) {
+ if (j == mirrorcell) {
+ //cout << "data for mirror cell " << mirrorcell << " from cell " << *(cellMap + j - numDualCells) << endl;
+ }
+
+
+ var_target = dualCellVarData + (*(cellMap + j - numDualCells))*maxNVertLevels;
+ if (singlelayer) {
+ if ((!writeBinary) && (j%6 == 0)) vtkFile << endl << " ";
+ writeFloat(vtkb64, vtkFile, *var_target, writeBinary);
+ } else {
+ for (int levelNum = 0; levelNum < maxNVertLevels; levelNum++)
+ {
+ if ((!writeBinary) && (counter%6 == 0)) vtkFile << endl << " ";
+ writeFloat(vtkb64, vtkFile, *var_target, writeBinary);
+ var_target++;
+ counter++;
+ }
+ }
+ }
+
+ if (writeBinary) {
+ vtkb64.EndWriting();
+ }
+ vtkFile << endl << " </DataArray>" << endl;
+
+ }
+
+ vtkFile << " </CellData>" << endl;
+
+ free(dualCellVarData);
+
+ vtkFile.close();
+
+ // append geoFile to vtkFile
+
+ string geoFileNameString = geoFileName.str();
+ string vtkFileNameString = vtkFileName.str();
+ int appRetVal = myAppendFile(&geoFileNameString, &vtkFileNameString);
+
+ }
+
+ return(0);
+
+}
+
Property changes on: branches/ocean_projects/graphics/paraview/projection/projlatlonbin.cpp
___________________________________________________________________
Name: svn:mime-type
+ text/x-c++src
Name: svn:keywords
+ Author Date ID Revision
Name: svn:eol-style
+ native
Added: branches/ocean_projects/graphics/paraview/translator/EncodeBase64.cpp
===================================================================
--- branches/ocean_projects/graphics/paraview/translator/EncodeBase64.cpp (rev 0)
+++ branches/ocean_projects/graphics/paraview/translator/EncodeBase64.cpp 2010-08-16 22:07:43 UTC (rev 477)
@@ -0,0 +1,169 @@
+using namespace std;
+
+#include "EncodeBase64.h"
+
+
+EncodeBase64::EncodeBase64(ofstream* astream)
+{
+ outstream = astream;
+ this->BufferLength = 0;
+}
+
+//----------------------------------------------------------------------------
+EncodeBase64::~EncodeBase64()
+{
+}
+
+
+
+//----------------------------------------------------------------------------
+static const unsigned char Base64EncodeTable[65] =
+"ABCDEFGHIJKLMNOPQRSTUVWXYZ"
+"abcdefghijklmnopqrstuvwxyz"
+"0123456789+/";
+
+//----------------------------------------------------------------------------
+inline static unsigned char EncodeChar(unsigned char c)
+{
+ assert( c < 65 );
+ return Base64EncodeTable[c];
+}
+
+//----------------------------------------------------------------------------
+void EncodeBase64::EncodePair(unsigned char i0,
+ unsigned char i1,
+ unsigned char *o0,
+ unsigned char *o1,
+ unsigned char *o2,
+ unsigned char *o3)
+{
+ *o0 = EncodeChar((i0 >> 2) & 0x3F);
+ *o1 = EncodeChar(((i0 << 4) & 0x30)|((i1 >> 4) & 0x0F));
+ *o2 = EncodeChar(((i1 << 2) & 0x3C));
+ *o3 = '=';
+}
+
+//----------------------------------------------------------------------------
+void EncodeBase64::EncodeSingle(unsigned char i0,
+ unsigned char *o0,
+ unsigned char *o1,
+ unsigned char *o2,
+ unsigned char *o3)
+{
+ *o0 = EncodeChar((i0 >> 2) & 0x3F);
+ *o1 = EncodeChar(((i0 << 4) & 0x30));
+ *o2 = '=';
+ *o3 = '=';
+}
+
+//----------------------------------------------------------------------------
+void EncodeBase64::EncodeTriplet(unsigned char i0,
+ unsigned char i1,
+ unsigned char i2,
+ unsigned char *o0,
+ unsigned char *o1,
+ unsigned char *o2,
+ unsigned char *o3)
+{
+ *o0 = EncodeChar((i0 >> 2) & 0x3F);
+ *o1 = EncodeChar(((i0 << 4) & 0x30)|((i1 >> 4) & 0x0F));
+ *o2 = EncodeChar(((i1 << 2) & 0x3C)|((i2 >> 6) & 0x03));
+ *o3 = EncodeChar(i2 & 0x3F);
+}
+
+//----------------------------------------------------------------------------
+int EncodeBase64::EncodeTriplet(unsigned char c0,
+ unsigned char c1,
+ unsigned char c2)
+{
+ // Encodes 3 bytes into 4 bytes and writes them to the output stream.
+ unsigned char out[4];
+ EncodeTriplet(c0, c1, c2,
+ &out[0], &out[1], &out[2], &out[3]);
+ return (outstream->write(reinterpret_cast<char*>(out), 4)? 1:0);
+}
+
+//----------------------------------------------------------------------------
+int EncodeBase64::EncodeEnding(unsigned char c0,
+ unsigned char c1)
+{
+ // Encodes a 2-byte ending into 3 bytes and 1 pad byte and writes.
+ unsigned char out[4];
+ EncodePair(c0, c1,
+ &out[0], &out[1], &out[2], &out[3]);
+ return (outstream->write(reinterpret_cast<char*>(out), 4)? 1:0);
+}
+
+//----------------------------------------------------------------------------
+int EncodeBase64::EncodeEnding(unsigned char c0)
+{
+ // Encodes a 1-byte ending into 2 bytes and 2 pad bytes and writes.
+ unsigned char out[4];
+ EncodeSingle(c0, &out[0], &out[1], &out[2], &out[3]);
+ return (outstream->write(reinterpret_cast<char*>(out), 4)? 1:0);
+}
+
+//----------------------------------------------------------------------------
+int EncodeBase64::StartWriting()
+{
+ this->BufferLength = 0;
+ this->WriteCount = 0;
+ //cout << "StartWriting!" << endl;
+ return 1;
+}
+
+//----------------------------------------------------------------------------
+int EncodeBase64::EndWriting()
+{
+ if(this->BufferLength == 1)
+ {
+ if(!this->EncodeEnding(this->Buffer[0])) { return 0; }
+ this->BufferLength = 0;
+ }
+ else if(this->BufferLength == 2)
+ {
+ if(!this->EncodeEnding(this->Buffer[0], this->Buffer[1])) { return 0; }
+ this->BufferLength = 0;
+ }
+ //cout << "EndWriting: WriteCount: " << WriteCount << endl;
+ return 1;
+}
+
+//----------------------------------------------------------------------------
+int EncodeBase64::Write(const unsigned char* data,
+ unsigned long length)
+{
+ unsigned long totalLength = this->BufferLength + length;
+ const unsigned char* in = data;
+ const unsigned char* end = data+length;
+
+ if(totalLength >= 3)
+ {
+ if(this->BufferLength == 1)
+ {
+ if(!this->EncodeTriplet(this->Buffer[0], in[0], in[1])) { return 0; }
+ in += 2;
+ this->BufferLength = 0;
+ }
+ else if(this->BufferLength == 2)
+ {
+ if(!this->EncodeTriplet(this->Buffer[0], this->Buffer[1], in[0]))
+ { return 0; }
+ in += 1;
+ this->BufferLength = 0;
+ }
+ }
+
+ while((end - in) >= 3)
+ {
+ if(!this->EncodeTriplet(in[0], in[1], in[2])) { return 0; }
+ in += 3;
+ }
+
+ while(in != end)
+ {
+ this->Buffer[this->BufferLength++] = *in++;
+ }
+ WriteCount++;
+ return 1;
+}
Property changes on: branches/ocean_projects/graphics/paraview/translator/EncodeBase64.cpp
___________________________________________________________________
Name: svn:mime-type
+ text/x-c++src
Name: svn:keywords
+ Author Date ID Revision
Name: svn:eol-style
+ native
Added: branches/ocean_projects/graphics/paraview/translator/EncodeBase64.h
===================================================================
--- branches/ocean_projects/graphics/paraview/translator/EncodeBase64.h (rev 0)
+++ branches/ocean_projects/graphics/paraview/translator/EncodeBase64.h 2010-08-16 22:07:43 UTC (rev 477)
@@ -0,0 +1,72 @@
+using namespace std;
+
+#include <iostream>
+#include <fstream>
+#include "stdlib.h"
+#include <assert.h>
+
+class EncodeBase64
+{
+public:
+ EncodeBase64(ofstream *outstream);
+ ~EncodeBase64();
+
+ // Description:
+ // Called after the stream position has been set by the caller, but
+ // before any Write calls. The stream position should not be
+ // adjusted by the caller until after an EndWriting call.
+ int StartWriting();
+
+ // Description:
+ // Write output data of the given length.
+ int Write(const unsigned char* data, unsigned long length);
+
+ // Description:
+ // Called after all desired calls to Write have been made. After
+ // this call, the caller is free to change the position of the
+ // stream. Additional writes should not be done until after another
+ // call to StartWriting.
+ int EndWriting();
+
+protected:
+
+ // Methods to encode and write data.
+
+ // Number of un-encoded bytes left in Buffer from last call to Write.
+ unsigned int BufferLength;
+ unsigned int WriteCount;
+ unsigned char Buffer[2];
+ ofstream *outstream;
+
+ // Description:
+ // Encode 2 bytes into 4 bytes
+ static void EncodePair(unsigned char i0,
+ unsigned char i1,
+ unsigned char *o0,
+ unsigned char *o1,
+ unsigned char *o2,
+ unsigned char *o3);
+
+ // Description:
+ // Encode 1 byte into 4 bytes
+ static void EncodeSingle(unsigned char i0,
+ unsigned char *o0,
+ unsigned char *o1,
+ unsigned char *o2,
+ unsigned char *o3);
+ // Description:
+ // Encode 3 bytes into 4 bytes
+ static void EncodeTriplet(unsigned char i0,
+ unsigned char i1,
+ unsigned char i2,
+ unsigned char *o0,
+ unsigned char *o1,
+ unsigned char *o2,
+ unsigned char *o3);
+
+ int EncodeTriplet(unsigned char c0, unsigned char c1, unsigned char c2);
+ int EncodeEnding(unsigned char c0, unsigned char c1);
+ int EncodeEnding(unsigned char c0);
+
+};
+
Property changes on: branches/ocean_projects/graphics/paraview/translator/EncodeBase64.h
___________________________________________________________________
Name: svn:mime-type
+ text/x-c++hdr
Name: svn:keywords
+ Author Date ID Revision
Name: svn:eol-style
+ native
Modified: branches/ocean_projects/graphics/paraview/translator/Makefile
===================================================================
--- branches/ocean_projects/graphics/paraview/translator/Makefile 2010-08-12 22:02:03 UTC (rev 476)
+++ branches/ocean_projects/graphics/paraview/translator/Makefile 2010-08-16 22:07:43 UTC (rev 477)
@@ -1,18 +1,26 @@
# NETCDF should point to the directory holding the netcdf lib and bin subdirs
#NETCDF = /usr/projects/climate/bzhao/netcdf-3.6.1
-NETCDF = /usr/projects/climate/mhecht/netcdf-3.6.1
+#NETCDF = /usr/projects/climate/mhecht/netcdf-3.6.1
#CC = gcc -m64
CC = gcc
INCLUDES = -I$(NETCDF)/include
-#LIBS = -L$(NETCDF)/lib -lnetcdf_c++ -lnetcdf -L/usr/lib -lstdc++
+LIBS = -L$(NETCDF)/lib -lnetcdf_c++ -lnetcdf -L/usr/lib -lstdc++
-LIBS = -L$(NETCDF)/lib -lnetcdf_c++ -lnetcdf \
- -lstdc++ \
- -L/opt/Intel/cce/10.0.023/lib/ -lirc
-# -L/opt/Intel/cce/9.1.039/lib/ -lirc
+#LIBS = -L$(NETCDF)/lib -lnetcdf_c++ -lnetcdf \
+# -lstdc++ \
+# -L/opt/Intel/cce/10.0.023/lib/ -lirc
transdual: nc2vtk_dual.cpp
$(CC) $(INCLUDES) -o $@ $< $(LIBS)
+.cpp.o:
+ $(CC) $(MYDEBUG) $(INCLUDES) -c $<
+
+EncodeBase64.o: EncodeBase64.cpp
+ $(CC) $(INCLUDES) -c $<
+
+transdualbin: nc2vtu_dualbin.cpp EncodeBase64.o
+ $(CC) $(INCLUDES) -o $@ $< EncodeBase64.o $(LIBS)
+
clean:
- rm transdual
+ rm transdual transdualbin *.o
Added: branches/ocean_projects/graphics/paraview/translator/nc2vtu_dualbin.cpp
===================================================================
--- branches/ocean_projects/graphics/paraview/translator/nc2vtu_dualbin.cpp (rev 0)
+++ branches/ocean_projects/graphics/paraview/translator/nc2vtu_dualbin.cpp 2010-08-16 22:07:43 UTC (rev 477)
@@ -0,0 +1,853 @@
+////////////////////////////////////////////////////////////////////////////////
+//
+// This program translates a newpop netCDF data file to dual-grid binary VTU format
+// The variables that have time dim are automatically written.
+//
+// Assume all variables are of interest.
+// Assume variable data type is double.
+// Assume variable dims are (Time, nCells|nVertices, nVertLevels|nVertLevelsP1, [nTracers])
+// Assume no more than 100 vars each for cell and point data
+// Does not deal with edge data.
+// Handles vectors.
+// Doubles converted to floats in .vtu files follow these rules:
+// if a NaN, then become -FLT_MAX.
+// if positive infinity, then become FLT_MAX
+// if negative infinity, then become -FLT_MAX
+// if smaller than a float, then become 0.
+// if not zero, but between -FLT_MIN and FLT_MIN make -FLT_MIN or FLT_MIN
+// if outside of -FLT_MAX and FLT_MAX make -FLT_MAX or FLT_MAX
+//
+//
+// Version 1.7
+// Christine Ahrens
+// 7/23/2010
+//
+////////////////////////////////////////////////////////////////////////////////
+
+
+#include <iostream>
+#include <fstream>
+#include <sstream>
+#include "stdlib.h"
+#include "vtkCellType.h"
+#include "netcdfcpp.h"
+#include <string>
+#include <cmath>
+#include <cfloat>
+#include "EncodeBase64.h"
+
+using namespace std;
+
+#define CHECK_MALLOC(ptr) \
+ if (ptr == NULL) { \
+ cerr << "malloc failed!</font>
<font color="blue">"; \
+ exit(1); \
+ }
+
+#define BUFF_SIZE 2048
+
+// append inFile to outFile
+int myAppendFile(string *inFile, string *outFile)
+{
+ char buff[BUFF_SIZE];
+ int readBytes = 1;
+
+ ifstream inFileStream(inFile->c_str(), ios::in|ios::binary);
+ if(!inFileStream)
+ {
+ return -1;
+ }
+
+ ofstream outFileStream(outFile->c_str(), ios::out|ios::app|ios::binary);
+ if(!outFileStream)
+ {
+ return -1;
+ }
+
+ while(readBytes != 0)
+ {
+ inFileStream.read((char*)buff, BUFF_SIZE);
+ readBytes = inFileStream.gcount();
+ outFileStream.write((char*)buff, readBytes);
+ }
+ return 0;
+}
+
+#define MAX_VARS 100
+
+int my_isinf_f(float x){
+ if (isinf(x)) {
+ return (x < 0.0 ? -1 : 1);
+ } else return 0;
+}
+
+unsigned int padSize(unsigned int binarySize) {
+ switch ((binarySize*4)%3) {
+ case 1:
+ binarySize = binarySize + 2;
+ break;
+ case 2:
+ binarySize++;
+ break;
+ case 0:
+ break;
+ }
+ return binarySize;
+}
+
+float convertDouble2ValidFloat(double inputData) {
+
+ // check for NaN
+ if (inputData != inputData) {
+ //cerr << "found NaN!" << endl;
+ return -FLT_MAX;
+ }
+
+ // check for infinity
+ int retval = my_isinf_f((float)inputData);
+ if (retval < 0) {
+ return -FLT_MAX;
+ } else if (retval > 0) {
+ return FLT_MAX;
+ }
+
+ // check number too small for float
+ if (abs(inputData) < 1e-126) return 0.0;
+
+ if ((float)inputData == 0) return 0.0;
+
+ if ((abs(inputData) > 0) && (abs(inputData) < FLT_MIN)) {
+ if (inputData < 0) return -FLT_MIN; else return FLT_MIN;
+ }
+
+ if (abs(inputData) > FLT_MAX) {
+ if (inputData < 0) return -FLT_MAX; else return FLT_MAX;
+ }
+
+ return (float)inputData;
+}
+
+int main(int argc, char* argv[])
+{
+ if ((argc < 3) || (argc > 4)) {
+ cerr << "Usage: transdualbin infile.nc outfile.vtu [verticalLevel]" << endl;
+ cerr << "Variables with time and vertical level are written out." << endl;
+ cerr << "Tracer vars are named tracer1, tracer2, etc." << endl;
+ cerr << "If vertical level is not specified, default is 0." << endl;
+ cerr << "A series of vtu files will be created, one file for each time step." << endl;
+ cerr << "with time prepended to the file name (e.g. 0outfile.vtu)." << endl;
+ cerr << "vtk datafile Version 2.0, transdualbin Version 1.7" << endl;
+ exit(1);
+ }
+
+
+ NcFile ncFile(argv[1]);
+
+ if (!ncFile.is_valid())
+ {
+ cerr << "Couldn't open file: " << argv[1] << endl;
+ exit(1);
+ }
+
+ NcDim* nCells = ncFile.get_dim("nCells");
+ NcDim* nVertices = ncFile.get_dim("nVertices");
+ NcDim* vertexDegree = ncFile.get_dim("vertexDegree");
+ NcDim* Time = ncFile.get_dim("Time");
+ NcDim* nVertLevels = ncFile.get_dim("nVertLevels");
+
+ // Can't check for this, b/c if not there it crashes program
+ //NcDim* nVertLevelsP1 = ncFile.get_dim("nVertLevelsP1");
+ int maxNVertLevels = nVertLevels->size() + 1;
+ /*if (nVertLevelsP1 != NULL) {
+ maxNVertLevels = nVertLevelsP1->size();
+ }
+*/
+
+ int outputVertLevel = 0;
+
+ if (argc == 4) {
+ outputVertLevel = atoi(argv[3]);
+ if (outputVertLevel > (maxNVertLevels-1)) {
+ cerr << "Specified vertical level " << outputVertLevel;
+ cerr << " doesn't exist. The highest level is level ";
+ cerr << maxNVertLevels-1 << "." << endl;
+ exit(1);
+ }
+ }
+
+
+ // figure out what variables to visualize
+ NcVar* dualCellVars[MAX_VARS];
+ NcVar* dualPointVars[MAX_VARS];
+ NcVar* vectorVars[3][3];
+ int dualCellVarIndex = -1;
+ int dualPointVarIndex = -1;
+ int numDualCells = nVertices->size();
+ int numDualPoints = nCells->size()+1;
+
+ int numVars = ncFile.num_vars();
+
+ bool tracersExist = false;
+ bool isUVector = false;
+ bool isVVector = false;
+ bool isWVector = false;
+
+ bool writeBinary = true;
+
+ for (int i = 0; i < numVars; i++) {
+ NcVar* aVar = ncFile.get_var(i);
+
+ // must have 3 dims
+ // (Time, nCells | nVertices, nVertLevels | nVertLevelsP1)
+
+ int numDims = aVar->num_dims();
+
+ if ((numDims != 3) && (strcmp(aVar->name(), "tracers"))) {
+ continue;
+ } else {
+ // TODO, check if it is a double
+ // assume a double for now
+
+ // check for Time dim 0
+ NcToken dim0Name = aVar->get_dim(0)->name();
+ if (strcmp(dim0Name, "Time"))
+ continue;
+
+ // check for dim 1 being num vertices or cells
+ bool isVertexData = false;
+ bool isCellData = false;
+ NcToken dim1Name = aVar->get_dim(1)->name();
+ if (!strcmp(dim1Name, "nVertices"))
+ isVertexData = true;
+ else if (!strcmp(dim1Name, "nCells"))
+ isCellData = true;
+ else continue;
+
+ // check if dim 2 is nVertLevels or nVertLevelsP1, too
+ NcToken dim2Name = aVar->get_dim(2)->name();
+ if ((strcmp(dim2Name, "nVertLevels"))
+ && (strcmp(dim2Name, "nVertLevelsP1"))) {
+ continue;
+ }
+
+ // check if we have data for the selected level
+ if (aVar->get_dim(2)->size()-1 < outputVertLevel) {
+ cout << "No data found for level ";
+ cout << outputVertLevel << " for variable ";
+ cout << aVar->name() << endl;
+ continue;
+ }
+
+ // Add to cell or point var array
+ if (isVertexData) { // means it is dual cell data
+ dualCellVarIndex++;
+ if (dualCellVarIndex > MAX_VARS-1) {
+ cerr << "Exceeded number of cell vars." << endl;
+ exit(1);
+ }
+ dualCellVars[dualCellVarIndex] = aVar;
+ } else if (isCellData) { // means it is dual point data
+ // check for u/v/w vectors
+
+ if (!strncmp(aVar->name(), "uReconstruct", 12)) {
+ isUVector = true;
+ switch (aVar->name()[12]) {
+ case 'X':
+ vectorVars[0][0] = aVar;
+ case 'Y':
+ vectorVars[0][1] = aVar;
+ case 'Z':
+ vectorVars[0][2] = aVar;
+ }
+ } else if (!strncmp(aVar->name(), "vReconstruct", 12)) {
+ isVVector = true;
+ switch (aVar->name()[12]) {
+ case 'X':
+ vectorVars[1][0] = aVar;
+ case 'Y':
+ vectorVars[1][1] = aVar;
+ case 'Z':
+ vectorVars[1][2] = aVar;
+ }
+ } else if (!strncmp(aVar->name(), "wReconstruct", 12)) {
+ isWVector = true;
+ switch (aVar->name()[12]) {
+ case 'X':
+ vectorVars[2][0] = aVar;
+ case 'Y':
+ vectorVars[2][1] = aVar;
+ case 'Z':
+ vectorVars[2][2] = aVar;
+ }
+ } else if (strcmp(aVar->name(), "tracers")) {
+ dualPointVarIndex++;
+ if (dualPointVarIndex > MAX_VARS-1) {
+ cerr << "Exceeded number of point vars." << endl;
+ exit(1);
+ }
+ dualPointVars[dualPointVarIndex] = aVar;
+ } else { // case of tracers, add each as "tracer0", "tracer1", etc.
+ tracersExist = true;
+ int numTracers = aVar->get_dim(3)->size();
+ for (int t = 0; t < numTracers; t++) {
+ dualPointVarIndex++;
+ if (dualPointVarIndex > MAX_VARS-1) {
+ cerr << "Exceeded number of point vars." << endl;
+ exit(1);
+ }
+ dualPointVars[dualPointVarIndex] = aVar;
+ }
+ }
+ }
+ }
+ }
+
+ // TODO
+ // prompt the user to find out which fields are of interest?
+ // for now, assume all are of interest
+
+ // Write geometry to a geo+outputfile.vtu
+ ostringstream geoFileName;
+ geoFileName << "geo_" << argv[2];
+ ofstream geoFile(geoFileName.str().c_str(), ios::out|ios::binary);
+ if (!geoFile)
+ {
+ cerr << "vtu geo output file could not be opened" <<endl;
+ exit(1);
+ }
+
+ // get points (centers of primal-mesh cells)
+
+ // TO DO check malloc return vals.
+
+ double *xCellData = (double*)malloc(nCells->size()
+ * sizeof(double));
+ CHECK_MALLOC(xCellData);
+ NcVar *xCellVar = ncFile.get_var("xCell");
+ xCellVar->get(xCellData, nCells->size());
+
+ double *yCellData = (double*)malloc(nCells->size()
+ * sizeof(double));
+ CHECK_MALLOC(yCellData);
+ NcVar *yCellVar = ncFile.get_var("yCell");
+ yCellVar->get(yCellData, nCells->size());
+
+ double *zCellData = (double*)malloc(nCells->size()
+ * sizeof(double));
+ CHECK_MALLOC(zCellData);
+ NcVar *zCellVar = ncFile.get_var("zCell");
+ zCellVar->get(zCellData, nCells->size());
+
+ // get num cells on each vertex
+ // how do I do that?
+ /*
+ int *cellsOnVertex = (int *) malloc((nVertices->size()) * vertexDegree->size() *
+ sizeof(int));
+ //cout << "ptr for cellsOnVertex" << cellsOnVertex << endl;
+ CHECK_MALLOC(cellsOnVertex);
+ NcVar *cellsOnVertexVar = ncFile.get_var("cellsOnVertex");
+ //cout << "getting cellsOnVertexVar</font>
<font color="blue">";
+ cellsOnVertexVar->get(cellsOnVertex, nVertices->size(), vertexDegree->size());
+*/
+
+ EncodeBase64 b64(&geoFile);
+
+ // write header
+
+ string formatString;
+ if (writeBinary) {
+ formatString = "binary";
+ } else {
+ formatString = "ascii";
+ }
+
+ unsigned int binarySize = 0;
+
+ // write points (the points are the primal-grid cell centers)
+
+ geoFile << " <Points>" << endl;
+
+ // see if we can change to double
+
+ geoFile << " <DataArray type=\"Float32\" Name=\"Points\" NumberOfComponents=\"3\" format=\"" << formatString << "\">" << endl;
+
+ // first write a dummy point, because the climate code
+ // starts their cell numbering at 1 and VTK starts it at
+ // 0
+ geoFile.precision(16);
+ float tmp = 0;
+
+ geoFile << " ";
+
+ if (writeBinary) {
+ b64.StartWriting();
+ binarySize = numDualPoints*sizeof(float)*3;
+ binarySize = padSize(binarySize);
+ //cout << "BinarySize: " << binarySize << endl;
+ b64.Write((const unsigned char*)&binarySize, sizeof(int));
+ b64.Write((const unsigned char*)&tmp, sizeof(float));
+ b64.Write((const unsigned char*)&tmp, sizeof(float));
+ b64.Write((const unsigned char*)&tmp, sizeof(float));
+ } else {
+ geoFile << (float)0.0 << " " << (float)0.0 << " " << (float)0.0 << endl;
+ }
+
+ for (int j = 0; j < nCells->size(); j++ )
+ {
+ geoFile.precision(16);
+ if (abs(xCellData[j]) < 1e-126) xCellData[j] = 0;
+ if (abs(yCellData[j]) < 1e-126) yCellData[j] = 0;
+ if (abs(zCellData[j]) < 1e-126) zCellData[j] = 0;
+ if (writeBinary) {
+ tmp = (float)xCellData[j];
+ b64.Write((const unsigned char*)&tmp, sizeof(float));
+ tmp = (float)yCellData[j];
+ b64.Write((const unsigned char*)&tmp, sizeof(float));
+ tmp = (float)zCellData[j];
+ b64.Write((const unsigned char*)&tmp, sizeof(float));
+ } else {
+ geoFile << " " << (float)xCellData[j] << " " << (float)yCellData[j] << " "
+ << (float)zCellData[j] << endl;
+ }
+ }
+ if (writeBinary) {
+ b64.EndWriting();
+ geoFile << endl;
+ }
+ geoFile << " </DataArray>" << endl;
+ geoFile << " </Points>" << endl;
+
+ free (xCellData);
+ free (yCellData);
+ free (zCellData);
+
+ // Write dual-mesh cells
+ // Dual-mesh cells are triangles with primal-mesh cell
+ // centers as the vertices.
+ // The number of dual-mesh cells is the number of vertices in the
+ // primal mesh.
+
+ // get dual-mesh cells
+
+ int *cellsOnVertex = (int *) malloc((nVertices->size()) * vertexDegree->size() *
+ sizeof(int));
+ CHECK_MALLOC(cellsOnVertex);
+ NcVar *cellsOnVertexVar = ncFile.get_var("cellsOnVertex");
+ cellsOnVertexVar->get(cellsOnVertex, nVertices->size(), vertexDegree->size());
+
+ geoFile << " <Cells>" << endl;
+ geoFile << " <DataArray type=\"Int32\" Name=\"connectivity\" format=\"" << formatString << "\">" << endl;
+
+ // for each dual-mesh cell, write the points for each
+
+ if (writeBinary) {
+ geoFile << " ";
+ b64.StartWriting();
+ binarySize = numDualCells*sizeof(int)*vertexDegree->size();
+ binarySize = padSize(binarySize);
+ //cout << "BinarySize: " << binarySize << endl;
+ b64.Write((const unsigned char*)&binarySize, sizeof(int));
+ }
+
+ for (int j = 0; j < numDualCells ; j++) {
+
+ if (!writeBinary) geoFile << " ";
+
+ // since primal vertex(pt) numbers == dual cell numbers
+ // we go through the primal vertices, find the cells around
+ // them, and since those primal cell numbers are dual
+ // point numbers,
+ // we can write the cell numbers for the cellsOnVertex
+ // and those will be the numbers of the dual vertices (pts).
+
+ int* dualCells = cellsOnVertex + (j * vertexDegree->size());
+
+ // go through and make sure none of the referenced points are
+ // out of range (<=0 or > nCells->size())
+
+ for (int k = 0; k < vertexDegree->size(); k++) {
+ if ((dualCells[k] <= 0) || (dualCells[k] > nCells->size())) {
+
+ // if any out of range exclude entire cell by setting
+ // all of its points to 0
+ for (int m = 0; m < vertexDegree->size(); m++) {
+ dualCells[m] = 0;
+ }
+ break;
+ }
+ }
+
+ for (int k = 0; k < vertexDegree->size(); k++)
+ {
+ if (writeBinary) {
+ b64.Write((const unsigned char*)&dualCells[k], sizeof(int));
+ } else {
+ geoFile << dualCells[k] << " ";
+ }
+ }
+ if (!writeBinary) geoFile << endl;
+ }
+
+ if (writeBinary) {
+ b64.EndWriting();
+ geoFile << endl;
+ }
+
+ geoFile << " </DataArray>" << endl;
+
+ free(cellsOnVertex);
+
+
+ // Specify the offset into connectivity array for the end of each cell
+ geoFile << " <DataArray type=\"Int32\" Name=\"offsets\" format=\"" << formatString << "\">";
+
+ if (writeBinary) {
+ geoFile << endl << " ";
+ b64.StartWriting();
+ binarySize = numDualCells*sizeof(int);
+ binarySize = padSize(binarySize);
+ //cout << "BinarySize: " << binarySize << endl;
+ b64.Write((const unsigned char*)&binarySize, sizeof(int));
+ }
+ int counter = 0;
+ for (int j = 0; j < numDualCells ; j++) {
+ if ((!writeBinary) && (counter%6 == 0)) {
+ geoFile << endl << " ";
+ }
+ counter += vertexDegree->size();
+ if (writeBinary) {
+ b64.Write((const unsigned char*)&counter, sizeof(int));
+ } else {
+ geoFile << counter << " ";
+ }
+ }
+ if (writeBinary) {
+ b64.EndWriting();
+ }
+ geoFile << endl;
+ geoFile << " </DataArray>" << endl;
+
+ // write cell types
+ int cellType;
+ if (vertexDegree->size() == 3)
+ cellType = VTK_TRIANGLE;
+ else if (vertexDegree->size() == 4)
+ cellType = VTK_QUAD;
+ else cellType = VTK_POLYGON;
+
+ geoFile << " <DataArray type=\"UInt8\" Name=\"types\" format=\"" << formatString << "\">";
+ if (writeBinary) {
+ geoFile << endl << " ";
+ b64.StartWriting();
+ binarySize = numDualCells*sizeof(char);
+ binarySize = padSize(binarySize);
+ //cout << "BinarySize: " << binarySize << endl;
+ b64.Write((const unsigned char*)&binarySize, sizeof(int));
+ }
+ for (int j = 0; j < numDualCells ; j++) {
+ if ((!writeBinary) && (j%6 == 0)) {
+ geoFile << endl << " ";
+ }
+ if (writeBinary) {
+ b64.Write((const unsigned char*)&cellType, sizeof(unsigned char));
+ } else {
+ geoFile << cellType << " ";
+ }
+ }
+ if (writeBinary) b64.EndWriting();
+ geoFile << endl;
+ geoFile << " </DataArray>" << endl;
+ geoFile << " </Cells>" << endl;
+ geoFile << " </Piece>" << endl;
+ geoFile << " </UnstructuredGrid>" << endl;
+ geoFile << "</VTKFile>" << endl;
+
+ geoFile.close();
+
+ // decls for data storage
+ double* dualCellVarData;
+ double* dualPointVarData;
+
+ // allocate space for variables
+
+ int varVertLevels = 0;
+
+
+
+ // for each time step, write a file with the time prepended to filename
+ for (int i = 0; i < Time->size(); i++)
+ {
+ ostringstream vtkFileName;
+ vtkFileName << i << argv[2];
+ ofstream vtkFile(vtkFileName.str().c_str(), ios::out|ios::binary);
+
+ vtkFile << "<VTKFile type=\"UnstructuredGrid\" version=\"0.1\" byte_order=\"LittleEndian\">" << endl;
+ vtkFile << " <UnstructuredGrid>" << endl;
+ vtkFile << " <Piece NumberOfPoints=\"" << numDualPoints << "\" NumberOfCells=\"" << numDualCells << "\">" << endl;
+
+ if (!vtkFile)
+ {
+ cerr << "vtk output file could not be opened" <<endl;
+ exit(1);
+ }
+
+
+
+ // Write attributes of dual-mesh cell data (attributes of primal-mesh
+ // vertex data)
+
+ vtkFile.precision(16);
+
+ // If by point, write out point data
+ vtkFile << " <PointData>" << endl;
+
+ dualPointVarData = (double*)malloc((sizeof(double)) * nCells->size());
+ CHECK_MALLOC(dualPointVarData);
+
+ int tracerNum = 0;
+ EncodeBase64 vtkb64(&vtkFile);
+
+ for (int v = 0; v <= dualPointVarIndex; v++) {
+
+ // Read variable number v data for that timestep
+
+ varVertLevels = dualPointVars[v]->get_dim(2)->size();
+
+ bool isTracer = false;
+
+ if (!strcmp(dualPointVars[v]->name(), "tracers")) {
+ isTracer = true;
+ dualPointVars[v]->set_cur(i, 0, outputVertLevel, tracerNum);
+ dualPointVars[v]->get(dualPointVarData, 1, nCells->size(), 1, 1);
+ } else {
+ dualPointVars[v]->set_cur(i, 0, outputVertLevel);
+ dualPointVars[v]->get(dualPointVarData, 1, nCells->size(), 1);
+ }
+
+ vtkFile << " <DataArray type=\"Float32\" Name=\"" << dualPointVars[v]->name();
+ if(isTracer) vtkFile << tracerNum+1;
+ vtkFile << "\" format=\"" << formatString << "\"> " << endl;
+
+ // Write variable number v data for that timestep
+
+ // get starting point
+ double *var_target = dualPointVarData;
+
+ float validData;
+ validData = convertDouble2ValidFloat (*var_target);
+
+ vtkFile << " ";
+
+ // write dummy
+ if (writeBinary) {
+ vtkb64.StartWriting();
+ binarySize = numDualPoints*sizeof(float);
+ binarySize = padSize(binarySize);
+ //cout << "BinarySize: " << binarySize << endl;
+ vtkb64.Write((const unsigned char*)&binarySize, sizeof(int));
+ vtkb64.Write((const unsigned char*)&validData, sizeof(float));
+ } else {
+ vtkFile << validData;
+ }
+
+ // write data
+ for (int j = 0; j < nCells->size(); j++)
+ {
+ if ((!writeBinary) && (j%6==0)) {
+ vtkFile << endl << " ";
+ }
+ validData = convertDouble2ValidFloat (*var_target);
+ if (writeBinary) {
+ vtkb64.Write((const unsigned char*)&validData, sizeof(float));
+ } else {
+ vtkFile << validData << " ";
+ }
+ var_target++;
+ }
+
+ if (writeBinary) vtkb64.EndWriting();
+ if (isTracer) tracerNum++;
+ vtkFile << endl << " </DataArray>" << endl;
+ }
+
+ free (dualPointVarData);
+
+ // write out vectors
+
+ float validData;
+ if (isUVector) {
+
+ vtkFile << " <DataArray type=\"Float32\" Name=\"uReconstruct\" NumberOfComponents=\"3\" format=\"" << formatString << "\"> " << endl;
+ // cout << "isUVector:" << endl;
+
+ double *vectorData[3] = {NULL, NULL, NULL};
+
+ for (int d = 0; d < 3; d++) {
+ vectorData[d] = (double*)malloc((sizeof(double)) * nCells->size());
+ CHECK_MALLOC(vectorData[d]);
+ }
+
+ for (int d = 0; d < 3; d++) {
+ // read data in
+ vectorVars[0][d]->set_cur(i, 0, outputVertLevel);
+ vectorVars[0][d]->get(vectorData[d], 1, nCells->size(), 1);
+ }
+
+ vtkFile << " ";
+ if (writeBinary) {
+ vtkb64.StartWriting();
+ binarySize = numDualPoints*sizeof(float)*3;
+ binarySize = padSize(binarySize);
+ //cout << "BinarySize: " << binarySize << endl;
+ vtkb64.Write((const unsigned char*)&binarySize, sizeof(int));
+ }
+
+ // dummy first
+
+ for (int d= 0; d < 3; d++) {
+ validData = convertDouble2ValidFloat (vectorData[d][0]);
+ if (writeBinary) {
+ vtkb64.Write((const unsigned char*)&validData, sizeof(validData));
+ } else {
+ vtkFile << validData << " ";
+ }
+ }
+ if (!writeBinary) vtkFile << endl;
+
+ // rest of data
+ for (int j = 0; j < nCells->size(); j++)
+ {
+ if (!writeBinary) vtkFile << " ";
+ for (int d= 0; d < 3; d++) {
+ validData = convertDouble2ValidFloat (vectorData[d][j]);
+ if (writeBinary) {
+ vtkb64.Write((const unsigned char*)&validData, sizeof(validData));
+ } else {
+ vtkFile << validData << " ";
+ }
+ }
+ if (!writeBinary) vtkFile << endl;
+ }
+ if (writeBinary) {
+ vtkb64.EndWriting();
+ vtkFile << endl;
+ }
+ vtkFile << " </DataArray>" << endl;
+ for (int d = 0; d < 3; d++) {
+ free(vectorData[d]);
+ }
+ }
+
+
+#ifdef NOT IMPLEMENTED YET
+ if (isVVector) {
+ for (int d = 0; d < 3; d++) {
+ // read data in
+ vectorVars[1][d]->set_cur(i, 0, outputVertLevel);
+ vectorVars[1][d]->get(vectorData[d], 1, nCells->size(), 1);
+ }
+ vtkFile << "VECTORS vReconstruct double" << endl;
+ for (int d= 0; d < 3; d++) {
+ validData = convertDouble2ValidFloat (vectorData[d][0]);
+ vtkFile << validData << " ";
+ }
+ vtkFile << endl;
+
+ for (int j = 0; j < nCells->size(); j++)
+ {
+ for (int d= 0; d < 3; d++) {
+ validData = convertDouble2ValidFloat (vectorData[d][j]);
+ vtkFile << validData << endl;
+ }
+ vtkFile << endl;
+ }
+ }
+
+ if (isWVector) {
+ for (int d = 0; d < 3; d++) {
+ // read data in
+ vectorVars[2][d]->set_cur(i, 0, outputVertLevel);
+ vectorVars[2][d]->get(vectorData[d], 1, nCells->size(), 1);
+ }
+ vtkFile << "VECTORS vReconstruct double" << endl;
+ for (int d= 0; d < 3; d++) {
+ validData = convertDouble2ValidFloat (vectorData[d][0]);
+ vtkFile << validData << " ";
+ }
+ vtkFile << endl;
+
+ for (int j = 0; j < nCells->size(); j++)
+ {
+ for (int d= 0; d < 3; d++) {
+ validData = convertDouble2ValidFloat (vectorData[d][j]);
+ vtkFile << validData << " ";
+ }
+ vtkFile << endl;
+ }
+ }
+#endif
+
+ vtkFile << " </PointData>" << endl;
+
+ // if by cell, then write out cell data
+
+ vtkFile << " <CellData>" << endl;
+
+ dualCellVarData = (double*)malloc((sizeof(double)) * numDualCells);
+ CHECK_MALLOC(dualCellVarData);
+
+ for (int v = 0; v <= dualCellVarIndex; v++) {
+
+ // Read variable number v data for that timestep
+ varVertLevels = dualCellVars[v]->get_dim(2)->size();
+ dualCellVars[v]->set_cur(i, 0, outputVertLevel);
+ dualCellVars[v]->get(dualCellVarData, 1, numDualCells, 1);
+
+ // Write variable number v data for that timestep
+ vtkFile << " <DataArray type=\"Float32\" Name=\"" << dualCellVars[v]->name() << "\" format=\"" << formatString << "\"> " ;
+
+ double *var_target = dualCellVarData;
+ float validData;
+ if (writeBinary) {
+ vtkFile << endl << " ";
+ vtkb64.StartWriting();
+ binarySize = numDualCells*sizeof(float);
+ binarySize = padSize(binarySize);
+ //cout << "BinarySize: " << binarySize << endl;
+ vtkb64.Write((const unsigned char*)&binarySize, sizeof(int));
+ }
+ for (int j = 0; j < numDualCells; j++)
+ {
+ if ((!writeBinary) && (j%6 == 0)) {
+ vtkFile << endl << " ";
+ }
+ validData = convertDouble2ValidFloat (*var_target);
+ if (writeBinary) {
+ vtkb64.Write((const unsigned char*)&validData, sizeof(float));
+ } else {
+ vtkFile << validData << " ";
+ }
+ var_target++;
+ }
+ if (writeBinary) {
+ vtkb64.EndWriting();
+ }
+ vtkFile << endl << " </DataArray>" << endl;
+ }
+ vtkFile << " </CellData>" << endl;
+
+ free(dualCellVarData);
+
+ vtkFile.close();
+
+ string geoFileNameString = geoFileName.str();
+ string vtkFileNameString = vtkFileName.str();
+ // append geoFile to vtkFile
+ int appRetVal = myAppendFile(&geoFileNameString, &vtkFileNameString);
+
+ }
+
+ return(0);
+
+}
+
Property changes on: branches/ocean_projects/graphics/paraview/translator/nc2vtu_dualbin.cpp
___________________________________________________________________
Name: svn:mime-type
+ text/x-c++src
Name: svn:keywords
+ Author Date ID Revision
Name: svn:eol-style
+ native
</font>
</pre>