<p><b>cahrens@lanl.gov</b> 2010-09-27 17:20:18 -0600 (Mon, 27 Sep 2010)</p><p>New version of mpasnc2vtu.cpp whic also can visualize primal sphere view (not latlon though).<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/graphics/paraview/translator/mpasnc2vtu.cpp
===================================================================
--- branches/ocean_projects/graphics/paraview/translator/mpasnc2vtu.cpp 2010-09-24 22:57:36 UTC (rev 511)
+++ branches/ocean_projects/graphics/paraview/translator/mpasnc2vtu.cpp 2010-09-27 23:20:18 UTC (rev 512)
@@ -2,8 +2,8 @@
// mpasnc2vtu.cpp
//
// This program translates a MPAS netCDF data file to a dual grid in ParaView
-// .vtu format. Can do either sphere view or lat/lon view. Can show single
-// or multi-layers. Can translate to xml with either ascii or binary data.
+// .vtu format for sphere and latlon views or a primal grid for sphere view.
+// Can show single or multi-layers. Can translate to xml with either ascii or binary data.
//
// Assumes little-endian platform.
// Assume all variables are of interest.
@@ -19,7 +19,7 @@
//
// Version 1.5
// Christine Ahrens
-// 9/21/2010
+// 9/27/2010
////////////////////////////////////////////////////////////////////////////////
@@ -41,6 +41,7 @@
bool projectLatLon;
bool writeBinary;
bool singleLayer;
+ bool showPrimalGrid;
bool isAtmosphere;
bool includeTopography;
bool doBugFix;
@@ -57,26 +58,38 @@
// points and cells refer to DUAL points and cells
struct geometry {
+ // from .nc file
+ int nCells;
+ int nVertices;
+ int maxEdges;
int nVertLevels;
int maxNVertLevels;
int Time;
+ int vertexDegree;
+
+ // computed values
int numCells;
int numPoints;
+ int modNumCells;
+ int modNumPoints;
int cellOffset;
int pointOffset;
int pointsPerCell;
int currentExtraPoint; // current extra point
int currentExtraCell; // current extra cell
+
+ // from .nc file
double* pointX; // x coord of point
double* pointY; // y coord of point
double* pointZ; // z coord of point
- int modNumPoints;
- int modNumCells;
+ int* nEdgesOnCell;
int* origConnections; // original connections
+
+ //computed values
int* modConnections; // modified connections
int* cellMap; // maps from added cell to original cell #
int* pointMap; // maps from added point to original point #
- int* maxLevelPoint; //
+ int* maxLevel; //
int maxCells; // max cells
int maxPoints; // max points
int verticalIndex; // for singleLayer, which vertical level
@@ -149,38 +162,6 @@
return binarySize;
}
-float ConvertDouble2ValidFloat(double inputData) {
-
- // check for NaN
- if (inputData != inputData) {
- //cerr << "found NaN!" << endl;
- return -FLT_MAX;
- }
-
- // check for infinity
- int retval = 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 ConvertDouble2ValidDouble(double inputData) {
// check for NaN
@@ -245,15 +226,6 @@
}
}
-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));
@@ -283,8 +255,12 @@
cerr << " multilayer view -- default is single layer view" << endl;
cerr << " -p " << endl;
cerr << " project to lat/lon -- default is sphere view" << endl;
+ cerr << " -r " << endl;
+ cerr << " show primal grid (works only for sphere view) -- default is dual grid" << endl;
cerr << " -t floatval " << endl;
- cerr << " layer_thickness -- use negative value for atmosphere, default is 10" << endl;
+ cerr << " layer_thickness -- default is 10" << endl;
+ cerr << " For sphere view, 100000 is a good value to start with." << endl;
+ cerr << " For latlon view, 10 is a good value to start with." << endl;
cerr << " -c floatval " << endl;
cerr << " display image with center_longitude in degrees, default is 180" << endl;
cerr << " -z " << endl;
@@ -299,7 +275,7 @@
cerr << " -f " << endl;
cerr << " timestep to finish on (defaults to last timestep available) " << endl;
cerr << " -g " << endl;
- cerr << " include topography (assumes you have maxLevelPoint variable) " << 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 << " -x " << endl;
@@ -331,6 +307,7 @@
p.startTimeStep = 0;
p.finishTimeStep = -1;
p.projectLatLon = false;
+ p.showPrimalGrid = false;
bool lflag = false;
bool tflag = false;
@@ -339,7 +316,7 @@
// -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:mazpxgbs:f:")) != -1) {
+ while ((c=getopt(argc, argv, "l:t:c:mazrpxgbs:f:")) != -1) {
switch(c) {
case 'l':
lflag = true;
@@ -389,6 +366,9 @@
case 'p':
p.projectLatLon = true;
break;
+ case 'r':
+ p.showPrimalGrid = true;
+ break;
case 'g':
p.includeTopography = true;
break;
@@ -450,22 +430,20 @@
}
-// put dims into dual grid geometry
+// put dims into grid geometry
void GetNcDims(NcFile &ncFile, geometry &g)
{
CHECK_DIM(ncFile, "nCells");
NcDim* nCells = ncFile.get_dim("nCells");
- g.numPoints = nCells->size();
- g.pointOffset = 1;
+ g.nCells = nCells->size();
CHECK_DIM(ncFile, "nVertices");
NcDim* nVertices = ncFile.get_dim("nVertices");
- g.numCells = nVertices->size();
- g.cellOffset = 0;
+ g.nVertices = nVertices->size();
CHECK_DIM(ncFile, "vertexDegree");
- NcDim* vertexDegree = ncFile.get_dim("vertexDegree");
- g.pointsPerCell = vertexDegree->size();
+ NcDim* vertexDegreeDim = ncFile.get_dim("vertexDegree");
+ g.vertexDegree = vertexDegreeDim->size();
CHECK_DIM(ncFile, "Time");
NcDim* Time = ncFile.get_dim("Time");
@@ -476,16 +454,26 @@
g.nVertLevels = nVertLevels->size();
g.maxNVertLevels = g.nVertLevels;
+
+ CHECK_DIM(ncFile, "maxEdges");
+ NcDim* maxEdgesDim = ncFile.get_dim("maxEdges");
+ g.maxEdges = maxEdgesDim->size();
+
}
void CheckArgs(params &p, geometry &g)
{
- if ((g.pointsPerCell != 3) && (g.pointsPerCell != 4)) {
- cerr << "This code is only for hexagonal or quad grid projection" << endl;
+ if ((g.vertexDegree != 3) && (g.vertexDegree != 4)) {
+ cerr << "This code is only for grid with vertexDegree 3 or 4." << endl;
exit(1);
}
+ if (p.projectLatLon && p.showPrimalGrid) {
+ cerr << "Cannot currently view latlon projection with primal grid." << endl;
+ exit(1);
+ }
+
// double-check we can do multilayer
if ((!p.singleLayer) && (g.maxNVertLevels == 1)) {
p.singleLayer = true;
@@ -610,9 +598,79 @@
}
+// allocate into sphere view of primal geometry
+void AllocPrimalSphereGeometry(NcFile &ncFile, params &p, geometry &g)
+{
+ g.numPoints = g.nVertices;
+ g.numCells = g.nCells;
+ g.pointsPerCell = g.maxEdges;
+ g.pointOffset = 1;
+
+ CHECK_VAR(ncFile, "xVertex");
+ g.pointX = (double*)malloc((g.numPoints+g.pointOffset) * sizeof(double));
+ CHECK_MALLOC(g.pointX);
+ NcVar* xVar = ncFile.get_var("xVertex");
+ xVar->get(g.pointX + g.pointOffset, g.numPoints);
+ // point 0 is 0.0
+ g.pointX[0] = 0.0;
+
+ CHECK_VAR(ncFile, "yVertex");
+ g.pointY = (double*)malloc((g.numPoints+g.pointOffset) * sizeof(double));
+ CHECK_MALLOC(g.pointY);
+ NcVar* yVar = ncFile.get_var("yVertex");
+ yVar->get(g.pointY + g.pointOffset, g.numPoints);
+ // point 0 is 0.0
+ g.pointY[0] = 0.0;
+
+ CHECK_VAR(ncFile, "zVertex");
+ g.pointZ = (double*)malloc((g.numPoints+g.pointOffset) * sizeof(double));
+ CHECK_MALLOC(g.pointZ);
+ NcVar* zVar = ncFile.get_var("zVertex");
+ zVar->get(g.pointZ + g.pointOffset, g.numPoints);
+ // point 0 is 0.0
+ g.pointZ[0] = 0.0;
+
+ CHECK_VAR(ncFile, "nEdgesOnCell");
+ g.nEdgesOnCell = (int *) malloc(g.numCells * sizeof(int));
+ CHECK_MALLOC(g.nEdgesOnCell);
+ NcVar *eVar = ncFile.get_var("nEdgesOnCell");
+ eVar->get(g.nEdgesOnCell, g.numCells);
+
+ CHECK_VAR(ncFile, "verticesOnCell");
+ g.origConnections = (int*)malloc(sizeof(int) * g.pointsPerCell * g.numCells);
+ CHECK_MALLOC(g.origConnections);
+ NcVar *verticesOnCellVar = ncFile.get_var("verticesOnCell");
+ // We need to account for the 0th cell connections.
+ verticesOnCellVar->get(g.origConnections, g.numCells, g.pointsPerCell);
+
+ if (p.includeTopography) {
+ CHECK_VAR(ncFile, "maxLevelCell");
+ g.maxLevel = (int*)malloc((g.numCells) * sizeof(int));
+ CHECK_MALLOC(g.maxLevel);
+ NcVar *maxLevelVar = ncFile.get_var("maxLevelCell");
+ maxLevelVar->get(g.maxLevel, g.numCells);
+ }
+
+ g.currentExtraPoint = g.numPoints + g.pointOffset;
+ g.currentExtraCell = g.numCells;
+
+ if (p.singleLayer) {
+ g.maxCells = g.currentExtraCell;
+ g.maxPoints = g.currentExtraPoint;
+ } else {
+ g.maxCells = g.currentExtraCell*g.maxNVertLevels;
+ g.maxPoints = g.currentExtraPoint*(g.maxNVertLevels+1);
+ }
+}
+
// allocate into sphere view of dual geometry
-void AllocSphereGeometry(NcFile &ncFile, params &p, geometry &g)
+void AllocDualSphereGeometry(NcFile &ncFile, params &p, geometry &g)
{
+ g.numPoints = g.nCells;
+ g.pointOffset = 1;
+ g.numCells = g.nVertices;
+ g.pointsPerCell = g.vertexDegree;
+
CHECK_VAR(ncFile, "xCell");
g.pointX = (double*)malloc((g.numPoints + g.pointOffset) * sizeof(double));
CHECK_MALLOC(g.pointX);
@@ -646,14 +704,14 @@
if (p.includeTopography) {
CHECK_VAR(ncFile, "maxLevelCell");
- g.maxLevelPoint = (int*)malloc((g.numPoints + g.numPoints) * sizeof(int));
- CHECK_MALLOC(g.maxLevelPoint);
- NcVar *maxLevelPointVar = ncFile.get_var("maxLevelCell");
- maxLevelPointVar->get(g.maxLevelPoint + g.pointOffset, g.numPoints);
+ g.maxLevel = (int*)malloc((g.numPoints + g.pointOffset) * sizeof(int));
+ CHECK_MALLOC(g.maxLevel);
+ NcVar *maxLevelVar = ncFile.get_var("maxLevelCell");
+ maxLevelVar->get(g.maxLevel + g.pointOffset, g.numPoints);
}
g.currentExtraPoint = g.numPoints + g.pointOffset;
- g.currentExtraCell = g.numCells + g.cellOffset;
+ g.currentExtraCell = g.numCells;
if (p.singleLayer) {
g.maxCells = g.currentExtraCell;
@@ -667,6 +725,11 @@
// allocate into lat/lon projection of dual geometry
void AllocLatLonGeometry(NcFile &ncFile, params &p, geometry &g)
{
+ g.numPoints = g.nCells;
+ g.pointOffset = 1;
+ g.numCells = g.nVertices;
+ g.pointsPerCell = g.vertexDegree;
+
const float BLOATFACTOR = .5;
g.modNumPoints = (int)floor(g.numPoints*(1.0 + BLOATFACTOR));
g.modNumCells = (int)floor(g.numCells*(1.0 + BLOATFACTOR))+1;
@@ -713,14 +776,14 @@
CHECK_MALLOC(g.cellMap);
g.currentExtraPoint = g.numPoints + g.pointOffset;
- g.currentExtraCell = g.numCells + g.cellOffset;
+ g.currentExtraCell = g.numCells;
if (p.includeTopography) {
CHECK_VAR(ncFile, "maxLevelCell");
- g.maxLevelPoint = (int*)malloc((g.numPoints + g.numPoints) * sizeof(int));
- CHECK_MALLOC(g.maxLevelPoint);
- NcVar *maxLevelPointVar = ncFile.get_var("maxLevelCell");
- maxLevelPointVar->get(g.maxLevelPoint + g.pointOffset, g.numPoints);
+ g.maxLevel = (int*)malloc((g.numPoints + g.numPoints) * sizeof(int));
+ CHECK_MALLOC(g.maxLevel);
+ NcVar *maxLevelVar = ncFile.get_var("maxLevelCell");
+ maxLevelVar->get(g.maxLevel + g.pointOffset, g.numPoints);
}
}
@@ -782,15 +845,22 @@
void FixPoints(params &p, geometry &g)
{
- for (int j = g.cellOffset; j < g.numCells + g.cellOffset; j++ ) {
+ for (int j = 0; j < g.numCells; j++ ) {
int *conns = g.origConnections + (j * g.pointsPerCell);
// go through and make sure none of the referenced points are
// out of range
// if so, set all to point 0
- for (int k = 0; k < g.pointsPerCell; k++) {
+ int pointsPer2DCell;
+ if (p.showPrimalGrid) {
+ pointsPer2DCell = g.nEdgesOnCell[j];
+ } else {
+ pointsPer2DCell = g.pointsPerCell;
+ }
+
+ for (int k = 0; k < pointsPer2DCell; k++) {
if ((conns[k] <= 0) || (conns[k] > g.numPoints)) {
- for (int m = 0; m < g.pointsPerCell; m++) {
+ for (int m = 0; m < pointsPer2DCell; m++) {
conns[m] = 0;
}
break;
@@ -799,7 +869,7 @@
int lastk;
- if (p.doBugFix) {
+ if (p.doBugFix && p.projectLatLon) {
//BUG FIX for problem where cells are stretching to a faraway point
lastk = g.pointsPerCell-1;
const double thresh = .06981317007977; // 4 degrees
@@ -833,7 +903,7 @@
int mirrorpoint;
- for (int j = g.cellOffset; j < g.numCells + g.cellOffset; j++ ) {
+ for (int j = 0; j < g.numCells; j++ ) {
int *conns = g.origConnections + (j * g.pointsPerCell);
int *modConns = g.modConnections + (j * g.pointsPerCell);
int lastk;
@@ -891,7 +961,7 @@
addedConns[k] = neigh;
}
}
- *(g.cellMap + (g.currentExtraCell - g.numCells - g.cellOffset)) = j;
+ *(g.cellMap + (g.currentExtraCell - g.numCells)) = j;
g.currentExtraCell++;
} else {
@@ -937,7 +1007,6 @@
// see if we can change to double
- //geoFile << " <DataArray type=\"Float32\" Name=\"Points\" NumberOfComponents=\"3\" format=\"" << formatString << "\">" << endl;
geoFile << " <DataArray type=\"Float64\" Name=\"Points\" NumberOfComponents=\"3\" format=\"" << formatString << "\">" << endl;
geoFile.precision(16);
@@ -957,7 +1026,6 @@
for (int j = 0; j < g.currentExtraPoint; j++ )
{
- //float x, y, z;
double x, y, z;
if (p.projectLatLon) {
@@ -1046,23 +1114,43 @@
if (p.writeBinary) {
geoFile << " ";
b64.StartWriting();
- int binarySize;
- if (p.singleLayer) {
- binarySize = g.maxCells*g.pointsPerCell*sizeof(int);
+
+ // figure out the binary size
+ int binarySize = 0;
+
+ // for primal grid, add up number of edges of each cell
+ if (p.showPrimalGrid) {
+
+ for (int j=0; j < g.currentExtraCell; j++) {
+ binarySize += g.nEdgesOnCell[j];
+ }
+
+ binarySize *= sizeof(int);
+
+ if (!p.singleLayer) {
+ binarySize *= g.maxNVertLevels;
+ }
+
+ // otherwise compute based on set cell size
} else {
- binarySize = g.maxCells*g.pointsPerCell*2*sizeof(int);
+ binarySize = g.maxCells * g.pointsPerCell * sizeof(int);
}
+
+ // if multilayer, mult by 2 since cells are 3D
+ if (!p.singleLayer) {
+ binarySize *= 2;
+ }
+
binarySize = PadSize(binarySize);
//cout << "BinarySize: " << binarySize << endl;
b64.Write((const unsigned char*)&binarySize, sizeof(int));
}
- // for each cell, write number of points for each
- // and then list the points by number
+ // For each cell, list the points by number
//cout << "Writing connectivity" << endl;
- unsigned int val;
+ int val;
for (int j = 0; j < g.currentExtraCell ; j++) {
@@ -1079,23 +1167,36 @@
int* connections;
//check if it is a mirror cell, if so, get original
- if (j >= g.numCells + g.cellOffset) {
+ if (j >= g.numCells) {
//cout << "setting origCell" << endl;
- int origCellNum = *(g.cellMap + (j - g.numCells - g.cellOffset));
+ int origCellNum = *(g.cellMap + (j - g.numCells));
//cout << "mirror cell: " <<j<< " origCell: "<< origCell << endl;
conns = g.origConnections + (origCellNum*g.pointsPerCell);
} else connections = g.origConnections + (j * g.pointsPerCell);
- //cout << "calc minLevel" << endl;
- minLevel = g.maxLevelPoint[conns[0]];
- //cout << "initial minLevel:" << minLevel << endl;
- // Take the min of the maxLevelPoint of each point
- for (int k = 1; k < g.pointsPerCell; k++) {
- minLevel = min(minLevel, g.maxLevelPoint[connections[k]]);
+ // if primal, just use maxLevel to determine the level of the cell
+ if (p.showPrimalGrid) {
+ minLevel = g.maxLevel[j];
+
+ // otherwise, take the min of the levels of the points of the cell
+ } else {
+ //cout << "calc minLevel" << endl;
+ minLevel = g.maxLevel[conns[0]];
+ //cout << "initial minLevel:" << minLevel << endl;
+ // Take the min of the maxLevelof each point
+ for (int k = 1; k < g.pointsPerCell; k++) {
+ minLevel = min(minLevel, g.maxLevel[connections[k]]);
+ }
}
- //cout << endl;
}
+ int pointsPer2DCell;
+ if (p.showPrimalGrid) {
+ pointsPer2DCell = g.nEdgesOnCell[j];
+ } else {
+ pointsPer2DCell = g.pointsPerCell;
+ }
+
if (p.singleLayer) {
if (!p.writeBinary) geoFile << " ";
@@ -1105,11 +1206,11 @@
if (p.includeTopography && ((minLevel-1) < p.outputVertLevel)) {
//cerr << "Setting all points to zero" << endl;
val = 0;
- for (int k = 0; k < g.pointsPerCell; k++) {
+ for (int k = 0; k < pointsPer2DCell; k++) {
WriteInt(b64, geoFile, val, p.writeBinary);
}
} else {
- for (int k = 0; k < g.pointsPerCell; k++) {
+ for (int k = 0; k < pointsPer2DCell; k++) {
WriteInt(b64, geoFile, conns[k], p.writeBinary);
}
}
@@ -1124,21 +1225,27 @@
if (p.includeTopography && ((minLevel-1) < levelNum)) {
// setting all points to zero
val = 0;
- for (int m = 0; m < g.pointsPerCell*2; m++) {
+ for (int m = 0; m < pointsPer2DCell*2; m++) {
WriteInt(b64, geoFile, val, p.writeBinary);
}
} else {
- for (int k = 0; k < g.pointsPerCell; k++)
+ for (int k = 0; k < pointsPer2DCell; k++)
{
val = (conns[k]*(g.maxNVertLevels+1)) + levelNum;
+ if (false) {
+ cout << "conn pt: " << (signed int)val << " j: " << j << " k: " << k << " levelNum: " << levelNum << " conns[k] = " << conns[k] << endl;
+ }
WriteInt(b64, geoFile, val, p.writeBinary);
}
- for (int k = 0; k < g.pointsPerCell; k++)
+ for (int k = 0; k < pointsPer2DCell; k++)
{
val = (conns[k]*(g.maxNVertLevels+1)) + levelNum + 1;
+ if (false) {
+ cout << "conn pt: " << (signed int)val << " j: " << j << " k: " << k << " levelNum: " << levelNum << " conns[k] = " << conns[k] << endl;
+ }
WriteInt(b64, geoFile, val, p.writeBinary);
}
}
@@ -1181,15 +1288,28 @@
int counter = 0;
unsigned int offset = 0;
- for (int j = 0; j < g.maxCells ; j++) {
- if ((!p.writeBinary) && (counter%6 == 0)) geoFile << endl << " ";
+ for (int j = 0; j < g.currentExtraCell ; j++) {
+
+ int pointsPer2DCell;
+ if (p.showPrimalGrid) {
+ pointsPer2DCell = g.nEdgesOnCell[j];
+ } else {
+ pointsPer2DCell = g.pointsPerCell;
+ }
+
if (p.singleLayer) {
- offset += g.pointsPerCell;
+ if ((!p.writeBinary) && (counter%6 == 0)) geoFile << endl << " ";
+ offset += pointsPer2DCell;
+ WriteInt(b64, geoFile, offset, p.writeBinary);
+ counter++;
} else {
- offset += 2*g.pointsPerCell;
+ for (int levelNum = 0; levelNum < g.maxNVertLevels; levelNum++) {
+ if ((!p.writeBinary) && (counter%6 == 0)) geoFile << endl << " ";
+ offset += 2*pointsPer2DCell;
+ WriteInt(b64, geoFile, offset, p.writeBinary);
+ counter++;
+ }
}
- WriteInt(b64, geoFile, offset, p.writeBinary);
- counter++;
}
if (p.writeBinary) {
@@ -1200,30 +1320,54 @@
geoFile << " </DataArray>" << endl;
}
+unsigned char GetCellType(int numPointsPer2DCell, bool make3D) {
-void WriteCellTypes(ofstream &geoFile, EncodeBase64 &b64,params &p, geometry &g)
-{
- // write cell types
unsigned char cellType;
- switch (g.pointsPerCell) {
+ switch (numPointsPer2DCell) {
case 3:
- if (p.singleLayer) {
+ if (!make3D) {
cellType = VTK_TRIANGLE;
} else {
cellType = VTK_WEDGE;
}
break;
case 4:
- if (p.singleLayer) {
+ if (!make3D) {
cellType = VTK_QUAD;
} else {
cellType = VTK_HEXAHEDRON;
}
break;
+ case 5:
+ if (!make3D) {
+ cellType = VTK_POLYGON;
+ } else {
+ cellType = VTK_PENTAGONAL_PRISM;
+ }
+ break;
+ case 6:
+ if (!make3D) {
+ cellType = VTK_POLYGON;
+ } else {
+ cellType = VTK_HEXAGONAL_PRISM;
+ }
+ break;
default:
break;
}
+ return cellType;
+}
+
+void WriteCellTypes(ofstream &geoFile, EncodeBase64 &b64,params &p, geometry &g)
+{
+ // write cell types
+ unsigned char cellType;
+
+ if (!p.showPrimalGrid) {
+ cellType = GetCellType(g.pointsPerCell, !p.singleLayer);
+ }
+
string formatString;
if (p.writeBinary) {
formatString = "binary";
@@ -1243,11 +1387,25 @@
b64.Write((const unsigned char*)&binarySize, sizeof(int));
}
- for (int j = 0; j < g.maxCells ; j++) {
- if ((!p.writeBinary) && (j%6 == 0)) {
- geoFile << endl << " ";
+ int counter = 0;
+
+ for (int j = 0; j < g.currentExtraCell ; j++) {
+
+ if (p.showPrimalGrid) {
+ cellType = GetCellType(g.nEdgesOnCell[j], !p.singleLayer);
}
- WriteChar(b64, geoFile, cellType, p.writeBinary);
+
+ if (p.singleLayer) {
+ if ((!p.writeBinary) && (counter%6 == 0)) geoFile << endl << " ";
+ WriteChar(b64, geoFile, cellType, p.writeBinary);
+ counter++;
+ } else {
+ for (int levelNum = 0; levelNum < g.maxNVertLevels; levelNum++) {
+ if ((!p.writeBinary) && (counter%6 == 0)) geoFile << endl << " ";
+ WriteChar(b64, geoFile, cellType, p.writeBinary);
+ counter++;
+ }
+ }
}
if (p.writeBinary) b64.EndWriting();
@@ -1299,8 +1457,12 @@
}
if (p.includeTopography) {
- free(g.maxLevelPoint);
+ free(g.maxLevel);
}
+
+ if (p.showPrimalGrid) {
+ free(g.nEdgesOnCell);
+ }
}
@@ -1341,7 +1503,6 @@
// Write variable number i data for that timestep
- //vtkFile << " <DataArray type=\"Float32\" Name=\"" <<
vtkFile << " <DataArray type=\"Float64\" Name=\"" <<
v.pointVars[i]->name();
vtkFile << "\" format=\"" << formatString << "\"> " << endl;
@@ -1350,7 +1511,6 @@
vtkFile << " ";
double *var_target;
- //float validData;
double validData;
if (p.writeBinary) {
@@ -1368,20 +1528,17 @@
// the range of the other data, so we look at the 1st value instead for
// our dummy
if (p.singleLayer) {
- var_target = pointVarData + 1;
- //WriteFloat(vtkb64, vtkFile, *var_target, p.writeBinary);
+ var_target = pointVarData + g.pointOffset;
WriteDouble(vtkb64, vtkFile, *var_target, p.writeBinary);
} else {
var_target = pointVarData + g.maxNVertLevels;
for (int levelNum = 0; levelNum < g.maxNVertLevels; levelNum++) {
- //WriteFloat(vtkb64, vtkFile, *var_target, p.writeBinary);
WriteDouble(vtkb64, vtkFile, *var_target, p.writeBinary);
var_target++;
}
// write highest level dummy point (duplicate of last level)
var_target--;
- //WriteFloat(vtkb64, vtkFile, *var_target, p.writeBinary);
WriteDouble(vtkb64, vtkFile, *var_target, p.writeBinary);
// Now set the var_target to the 1st data value
@@ -1391,14 +1548,12 @@
for (int j = g.pointOffset; j < g.numPoints + g.pointOffset; j++) {
if (p.singleLayer) {
- //WriteFloat(vtkb64, vtkFile, *var_target, p.writeBinary);
WriteDouble(vtkb64, vtkFile, *var_target, p.writeBinary);
var_target++;
} else {
// write data for one point -- lowest level to highest
for (int levelNum = 0; levelNum < g.maxNVertLevels; levelNum++) {
- //WriteFloat(vtkb64, vtkFile, *var_target, p.writeBinary);
WriteDouble(vtkb64, vtkFile, *var_target, p.writeBinary);
var_target++;
}
@@ -1406,7 +1561,6 @@
// for last layer of points, repeat last level's values
// Need Mark's input on this one
var_target--;
- //WriteFloat(vtkb64, vtkFile, *var_target, p.writeBinary);
WriteDouble(vtkb64, vtkFile, *var_target, p.writeBinary);
var_target++;
}
@@ -1418,15 +1572,11 @@
var_target = pointVarData + ((*(g.pointMap + j - g.numPoints - g.pointOffset))*g.maxNVertLevels);
if (p.singleLayer) {
- //WriteFloat(vtkb64, vtkFile, *var_target, p.writeBinary);
WriteDouble(vtkb64, vtkFile, *var_target, p.writeBinary);
- //validData = ConvertDouble2ValidFloat (*var_target);
- validData = ConvertDouble2ValidDouble (*var_target);
var_target++;
} else {
// write data for one point -- lowest level to highest
for (int levelNum = 0; levelNum < g.maxNVertLevels; levelNum++) {
- //WriteFloat(vtkb64, vtkFile, *var_target, p.writeBinary);
WriteDouble(vtkb64, vtkFile, *var_target, p.writeBinary);
var_target++;
}
@@ -1434,7 +1584,6 @@
// for last layer of points, repeat last level's values
// Need Mark's input on this one
var_target--;
- //WriteFloat(vtkb64, vtkFile, *var_target, p.writeBinary);
WriteDouble(vtkb64, vtkFile, *var_target, p.writeBinary);
}
}
@@ -1476,17 +1625,15 @@
for (int i = 0; i < v.numCellVars; i++) {
// Write variable number v data for that timestep
- //vtkFile << " <DataArray type=\"Float32\" Name=\"" <<
vtkFile << " <DataArray type=\"Float64\" Name=\"" <<
v.cellVars[i]->name() << "\" format=\"" << formatString << "\"> " ;
// Read variable number v data for that timestep and level
v.cellVars[i]->set_cur(t, 0, vertLevelIndex);
if (p.singleLayer) {
- v.cellVars[i]->get(cellVarData + g.cellOffset, 1, g.numCells,1);
+ v.cellVars[i]->get(cellVarData, 1, g.numCells,1);
} else {
- v.cellVars[i]->get(cellVarData + (g.cellOffset * g.maxNVertLevels),
- 1, g.numCells, g.maxNVertLevels);
+ v.cellVars[i]->get(cellVarData, 1, g.numCells, g.maxNVertLevels);
}
if (p.writeBinary) {
@@ -1502,12 +1649,11 @@
double *var_target = cellVarData;
int counter = 0;
- for (int j = g.cellOffset; j < g.numCells + g.cellOffset; j++) {
+ for (int j = 0; j < g.numCells; j++) {
if (p.singleLayer) {
if ((!p.writeBinary) && (j%6 == 0)) vtkFile << endl
<< " ";
- //WriteFloat(vtkb64, vtkFile, *var_target, p.writeBinary);
WriteDouble(vtkb64, vtkFile, *var_target, p.writeBinary);
var_target++;
} else {
@@ -1515,7 +1661,6 @@
{
if ((!p.writeBinary) && (counter%6 == 0)) vtkFile << endl
<< " ";
- //WriteFloat(vtkb64, vtkFile, *var_target, p.writeBinary);
WriteDouble(vtkb64, vtkFile, *var_target, p.writeBinary);
var_target++;
counter++;
@@ -1525,22 +1670,20 @@
// put out for extra cells
counter = 0;
- for (int j = g.numCells + g.cellOffset; j < g.currentExtraCell; j++) {
+ for (int j = g.numCells; j < g.currentExtraCell; j++) {
var_target = cellVarData
- + (*(g.cellMap + j - g.numCells - g.cellOffset))
+ + (*(g.cellMap + j - g.numCells))
*g.maxNVertLevels;
if (p.singleLayer) {
if ((!p.writeBinary) && (j%6 == 0))
vtkFile << endl << " ";
- //WriteFloat(vtkb64, vtkFile, *var_target, p.writeBinary);
WriteDouble(vtkb64, vtkFile, *var_target, p.writeBinary);
} else {
for (int levelNum = 0; levelNum < g.maxNVertLevels; levelNum++)
{
if ((!p.writeBinary) && (counter%6 == 0))
vtkFile << endl << " ";
- //WriteFloat(vtkb64, vtkFile, *var_target, p.writeBinary);
WriteDouble(vtkb64, vtkFile, *var_target, p.writeBinary);
var_target++;
counter++;
@@ -1608,6 +1751,7 @@
// figure out geometry of a dual-grid lat/lon projection
geometry g;
GetNcDims(ncFile, g);
+
CheckArgs(p, g);
if (p.projectLatLon) {
AllocLatLonGeometry(ncFile, p, g);
@@ -1615,7 +1759,12 @@
FixPoints(p, g);
EliminateXWrap(p, g);
} else {
- AllocSphereGeometry(ncFile, p, g);
+ if (p.showPrimalGrid) {
+ AllocPrimalSphereGeometry(ncFile, p, g);
+ } else {
+ AllocDualSphereGeometry(ncFile, p, g);
+ }
+ FixPoints(p, g);
}
// write a file with the geometry.
@@ -1625,7 +1774,11 @@
// figure out what variables to visualize
vars v;
- GetNcVars(ncFile, v, "nVertices", "nCells");
+ if (p.showPrimalGrid) {
+ GetNcVars(ncFile, v, "nCells", "nVertices");
+ } else {
+ GetNcVars(ncFile, v, "nVertices", "nCells");
+ }
// For each timestep, write data for each variable
for (int t = p.startTimeStep; t <= p.finishTimeStep; t++) {
</font>
</pre>