<p><b>cahrens@lanl.gov</b> 2010-04-15 18:07:00 -0600 (Thu, 15 Apr 2010)</p><p>Fixed transdual and projlatlon to handle quad grid (excluding land points -- Mat's grid) and put in command-line options to projlatlon to display either single layer or multilayer views and allow you to set layer thickness, vertical level and/or center longitude. Run without arguments to see usage.<br>
</p><hr noshade><pre><font color="gray">Modified: branches/ocean_projects/graphics/paraview/projection/projlatlon.cpp
===================================================================
--- branches/ocean_projects/graphics/paraview/projection/projlatlon.cpp        2010-04-15 22:21:11 UTC (rev 189)
+++ branches/ocean_projects/graphics/paraview/projection/projlatlon.cpp        2010-04-16 00:07:00 UTC (rev 190)
@@ -11,16 +11,16 @@
// 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
+// 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 Beta
+// Version 1.0
// Christine Ahrens
-// 3/8/2010
+// 4/13/2010
////////////////////////////////////////////////////////////////////////////////
@@ -38,10 +38,10 @@
using namespace std;
#define CHECK_MALLOC(ptr) \
-        if (ptr == NULL) { \
-                cerr << "malloc failed!</font>
<font color="red">"; \
-                exit(1); \
-        }
+ if (ptr == NULL) { \
+ cerr << "malloc failed!</font>
<font color="gray">"; \
+ exit(1); \
+ }
#define MAX_VARS 100
#define DEFAULT_LAYER_THICKNESS 10
@@ -86,7 +86,7 @@
// check for NaN
if (inputData != inputData) {
- cerr << "found NaN!" << endl;
+ //cerr << "found NaN!" << endl;
return -FLT_MAX;
}
@@ -114,253 +114,397 @@
return (float)inputData;
}
+void printUsage() {
+ cerr << "USAGE: projlatlon [OPTIONS] infile.nc outfile.vtk" << endl;
+ cerr << "Create lat/lon projection of MPAS-style data from input file infile.nc " << endl;
+ cerr << "and write to legacy VTK format for use in ParaView." << endl;
+ cerr << "A series of vtk files are created, one file for each time step." << endl;
+ cerr << "with time prepended to the file name (e.g. 0outfile.vtk)." << 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 0" << endl;
+ cerr << " -c floatval " << endl;
+ cerr << " center_longitude in degrees, default is 180" << 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, projlatlon Version 1.0]" << endl;
+ exit(1);
+}
+void getArgs(int argc, char* argv[], int &outputVertLevel,
+ float &layerThickness, double &centerLon, bool &isMultilayer,
+ char* &inputFileName, char* &outputFileName) {
+
+ int c;
+ bool selectedVar = false;
+ optind = 1;
+
+ bool lflag = false;
+ bool tflag = false;
+ bool cflag = false;
+ bool mflag = 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:m")) != -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;
+ default:
+ cerr << "Unrecognized option: -" << c << endl;
+ printUsage();
+ exit(1);
+ break;
+ }
+ //argc -= optind;
+ //argv += optind;
+ }
+
+ inputFileName = argv[optind++];
+ outputFileName = argv[optind];
+
+}
+
+
int main(int argc, char* argv[])
{
-        if ((argc < 3) || (argc > 4)) {
-                cerr << "Usage: projlatlon infile.nc infile.vtk [layer_thickness]" << endl;
-                cerr << "Note: layer_thickness defaults to 10." << endl;
-                exit(1);
-        }
-                
-        NcFile ncFile(argv[1]);
+ bool multilayer = false;
+ int outputVertLevel = 0;
+ char* inputFileName = NULL;
+ char* outputFileName = NULL;
+ float layerThickness = DEFAULT_LAYER_THICKNESS;
+ double centerLon = 180;
-        if (!ncFile.is_valid())
-        {
-                cerr << "Couldn't open file: " << argv[1] << endl;
-                exit(1);
-        }
+ if (argc < 3) {
+ printUsage();
+ 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");
+ getArgs(argc, argv, outputVertLevel, layerThickness, centerLon, multilayer,
+ inputFileName, outputFileName);
+
+ NcFile ncFile(inputFileName);
-        if (vertexDegree->size() != 3) {
-                cerr << "This code is only for hexagonal-primal/triangular-dual grid" << endl;
-                exit(1);
-        }
+ if (!ncFile.is_valid())
+ {
+ cerr << "Couldn't open file: " << inputFileName << endl;
+ exit(1);
+ }
-        // Can't check for this, b/c if not there it crashes program        
-        //NcDim* nVertLevelsP1 = ncFile.get_dim("nVertLevelsP1");
-        int maxNVertLevels = nVertLevels->size();
+ 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");
+ int maxNVertLevels = nVertLevels->size();
-        //cout << "maxNVertLevels: " << maxNVertLevels << endl;
-
-        int outputVertLevel = 0;
+ if ((vertexDegree->size() != 3) && (vertexDegree->size() != 4)) {
+ cerr << "This code is only for hexagonal or quad grid projection" << endl;
+ exit(1);
+ }
-        float layerThickness = DEFAULT_LAYER_THICKNESS;
-        if (argc == 4) {
-                layerThickness = atof(argv[3]);
-        }
+ // Can't check for this, b/c if not there it crashes program
+ //NcDim* nVertLevelsP1 = ncFile.get_dim("nVertLevelsP1");
-        // figure out what variables to visualize
-        NcVar* dualCellVars[MAX_VARS];
-        NcVar* dualPointVars[MAX_VARS];
-        int dualCellVarIndex = -1;
-        int dualPointVarIndex = -1;
-        int numDualCells = nVertices->size();
-        int numDualPoints = nCells->size()+1;
-        
-        int numVars = ncFile.num_vars();
+ bool singlelayer = !multilayer;
-        bool tracersExist = false;
+ // double-check we can do multilayer
+ if ((multilayer) && (maxNVertLevels == 1)) {
+ multilayer = false;
+ singlelayer = !multilayer;
+ }
-        for (int i = 0; i < numVars; i++) {
-                NcVar* aVar = ncFile.get_var(i);
+ if (singlelayer && (outputVertLevel > (maxNVertLevels-1))) {
+ cerr << "There is no data for level " << outputVertLevel << "." << endl;
+ exit(1);
+ }
-                // must have 3 dims
-                // (Time, nCells | nVertices, nVertLevels | nVertLevelsP1)
+ if (singlelayer && (outputVertLevel < 0)) {
+ cerr << "Level number must be 0 or greater." << endl;
+ exit(1);
+ }
-                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
+ if ((centerLon < 0) || (centerLon > 360)) {
+ cerr << "Center longitude must be between 0 and 360." << endl;
+ exit(1);
+ }
-                        // check for Time dim 0
-                        NcToken dim0Name = aVar->get_dim(0)->name();
-                        if (strcmp(dim0Name, "Time"))
-                                continue;
+ if (singlelayer) {
+ maxNVertLevels = 1;
+ }
-                        // 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;
+ cerr << "centerLon: " << centerLon << " singleLayer: " << singlelayer << " outputVertLevel: " << outputVertLevel << endl;
-                        // check if dim 2 is nVertLevels or nVertLevelsP1, too
-                        NcToken dim2Name = aVar->get_dim(2)->name();
-                        if ((strcmp(dim2Name, "nVertLevels"))
-                                        && (strcmp(dim2Name, "nVertLevelsP1"))) {
-                                continue;
-                        }
+ // figure out what variables to visualize
+ NcVar* dualCellVars[MAX_VARS];
+ NcVar* dualPointVars[MAX_VARS];
+ int dualCellVarIndex = -1;
+ int dualPointVarIndex = -1;
+ int numDualCells = nVertices->size();
+ int numDualPoints = nCells->size()+1;
+
+ int numVars = ncFile.num_vars();
-                        // 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
-                                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;
-                                        }
-                                }
-                        }
-                }
-        }
+ bool tracersExist = false;
-        // 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)
+ for (int i = 0; i < numVars; i++) {
+ NcVar* aVar = ncFile.get_var(i);
-        // TO DO check malloc return vals.
+ // 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
+ 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)
+
+ // TO DO check malloc return vals.
+
const float BLOATFACTOR = .5;
-        int myCellSize = floor(nCells->size()*(1.0 + BLOATFACTOR));
-        double *xCellData = (double*)malloc(myCellSize * sizeof(double));
-        CHECK_MALLOC(xCellData);
-        NcVar *xCellVar = ncFile.get_var("lonCell");
-        xCellVar->get(xCellData+1, nCells->size());
+ int myCellSize = (int)floor(nCells->size()*(1.0 + BLOATFACTOR));
+ 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;
-        double *yCellData = (double*)malloc(myCellSize * sizeof(double));
-        CHECK_MALLOC(yCellData);
-        NcVar *yCellVar = ncFile.get_var("latCell");
-        yCellVar->get(yCellData+1, nCells->size());
+ const double PI = 3.141592;
+
+ 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;
+ }
+ }
+ }
+ }
+
+ 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
+ // get dual-mesh cells
-        int *cellsOnVertex = (int *) malloc((nVertices->size()) * vertexDegree->size() *
-                sizeof(int));
-        CHECK_MALLOC(cellsOnVertex);
+ int *cellsOnVertex = (int *) malloc((nVertices->size()) * vertexDegree->size() *
+ sizeof(int));
+ CHECK_MALLOC(cellsOnVertex);
-        int myCOVSize = floor(nVertices->size()*(1.0 + BLOATFACTOR))+1;
-        int *myCellsOnVertex = (int *) malloc(myCOVSize * vertexDegree->size() * sizeof(int));
-        CHECK_MALLOC(myCellsOnVertex);
+ 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="red">";
-        cellsOnVertexVar->get(cellsOnVertex, nVertices->size(), vertexDegree->size());
+ 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(floor(nCells->size()*BLOATFACTOR) * sizeof(int));
-        CHECK_MALLOC(pointMap);
-        int *cellMap = (int*)malloc(floor(nVertices->size()*BLOATFACTOR) * sizeof(int));
-        CHECK_MALLOC(cellMap);
+ // 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;
- const double PI = 3.141592;
+
+ 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.
+ // 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());
-                int 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;
-                }
+ for (int j = 0; j < numDualCells; j++ ) {
+ int *dualCells = cellsOnVertex + (j * vertexDegree->size());
-                if (xWrap || yWrap) {
-                        //cerr << "Cell wrap: " << j << endl;
-                }
+ // 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;
+ }
+ }
-                if (xWrap) {
-                        //cerr << "It wrapped in x direction" << endl;
-                        double anchorX = xCellData[dualCells[0]];
-                        double anchorY = yCellData[dualCells[0]];
-                        *myptr = *(dualCells);
-                        myptr++;
+ int 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
+ // modify existing cell, so it doesn't wrap
+ // move points to one side
-                        // first point is anchor it doesn't move
+ // 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]];
+ 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 a new point, figure out east or west
+ if (abs(neighX - anchorX) > 5.5) {
+ double neighEastX;
+ double neighWestX;
-                                        // add on east
-                                        if (neighX < anchorX) {
+ // 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
+ 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;
-                                        }
-                        
+ neighWestX = neighX - (2*PI);
+ xCellData[currentExtraPoint] = neighWestX;
+ yCellData[currentExtraPoint] = neighY;
+ }
+
if (j == acell) {
//cout << "x: " << xCellData[currentExtraPoint] << " y: " << yCellData[currentExtraPoint] << endl;
}
@@ -369,11 +513,11 @@
//cout << "currentExtraPoint: " << currentExtraPoint << endl;
}
-                                        // add the new point to list of vertices
-                                        *myptr = currentExtraPoint;
+ // add the new point to list of vertices
+ *myptr = currentExtraPoint;
-                                        // record mapping
-                                        *(pointMap + (currentExtraPoint - numDualPoints)) = dualCells[k];
+ // record mapping
+ *(pointMap + (currentExtraPoint - numDualPoints)) = dualCells[k];
if (j == acell) {
mirrorpoint = currentExtraPoint;
@@ -381,28 +525,28 @@
//cout << "mirror point " << mirrorpoint << " has pointMap to: " << dualCells[k] << endl;
}
-                                        myptr++;
-                                        currentExtraPoint++;
-                                } else {
-                                        // use existing kth point
+ myptr++;
+ currentExtraPoint++;
+ } else {
+ // use existing kth point
if (j == acell) {
//cout << "use existing point" << endl;
}
-                                        *myptr = dualCells[k];
-                                        myptr++;
-                                }
-                        }
+ *myptr = dualCells[k];
+ myptr++;
+ }
+ }
-                        // add a mirror image cell on other side so there are no
-                        // gaps in the map
+ // add a mirror image cell on other side so there are no
+ // gaps in the map
-                        // move anchor to other side
-                        if (anchorX > PI) {
-                                anchorX = anchorX - (2*PI);
-                        } else {
-                                anchorX = anchorX + (2*PI);
-                        }
-        
+ // 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;
@@ -412,367 +556,467 @@
int* addedCellsPtr = myCellsOnVertex + (currentExtraCell * vertexDegree->size());
// add point coord and add to list of cells
-                        xCellData[currentExtraPoint] = anchorX;
-                        yCellData[currentExtraPoint] = anchorY;
-                        *addedCellsPtr = currentExtraPoint;
+ xCellData[currentExtraPoint] = anchorX;
+ yCellData[currentExtraPoint] = anchorY;
+ *addedCellsPtr = currentExtraPoint;
-                        // record mapping
-                        *(pointMap + (currentExtraPoint - numDualPoints)) = dualCells[0];
+ // record mapping
+ *(pointMap + (currentExtraPoint - numDualPoints)) = dualCells[0];
-                        addedCellsPtr++;
-                        currentExtraPoint++;
+ addedCellsPtr++;
+ currentExtraPoint++;
-                        for (int k = 1; k < vertexDegree->size(); k++) {
-                                double neighX = xCellData[dualCells[k]];
-                                double neighY = yCellData[dualCells[k]];
+ 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 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 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;
+ // add the new point to list of vertices
+ *addedCellsPtr = currentExtraPoint;
-                                        // record mapping
-                                        *(pointMap + (currentExtraPoint - numDualPoints)) = dualCells[k];
+ // record mapping
+ *(pointMap + (currentExtraPoint - numDualPoints)) = dualCells[k];
-                                        addedCellsPtr++;
-                                        currentExtraPoint++;
-                                } else {
-                                        // use existing kth point
-                                        *addedCellsPtr = dualCells[k];
-                                        addedCellsPtr++;
-                                }
-                        }
-                        *(cellMap + (currentExtraCell - numDualCells)) = j;
-                        currentExtraCell++;
-                }
+ 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 (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;
+ }
+ }
+
+ // decls for data storage
+ double* dualCellVarData;
+ double* dualPointVarData;
-                // 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;
-                }
-        }        
-                                                
-        // decls for data storage
-        double* dualCellVarData;
-        double* dualPointVarData;
+ // for each variable, allocate space for variables
-        // for each variable, allocate space for variables
+ //cout << "dualCellVarIndex: " << dualCellVarIndex << endl;
-        //cout << "dualCellVarIndex: " << dualCellVarIndex << endl;
+ int varVertLevels = 0;
+
+ // write a file with the geometry.
-        int varVertLevels = 0;
-        
-        // write a file with the geometry.
+ ostringstream geoFileName;
-        ostringstream geoFileName;
+ geoFileName << "geo_" << outputFileName;
-        geoFileName << "geo_" << argv[2];         
+ ofstream geoFile(geoFileName.str().c_str(), ios::out);
-        ofstream geoFile(geoFileName.str().c_str(), ios::out);
+ if (!geoFile)
+ {
+ cerr << "vtk output file could not be opened" <<endl;
+ exit(1);
+ }
-        if (!geoFile)
-        {
-                cerr << "vtk output file could not be opened" <<endl;
-                exit(1);
-        }
+ // write header
-        // write header
+ geoFile << "# vtk DataFile Version 2.0" << endl;
+ geoFile << "Project newpop geometry to lat/lon projection from netCDF by Christine Ahrens"
+ << endl;
+ geoFile << "ASCII" << endl;
+ geoFile << "DATASET UNSTRUCTURED_GRID" << endl;
-        geoFile << "# vtk DataFile Version 2.0" << endl;
-        geoFile << "Project newpop geometry to lat/lon projection from netCDF by Christine Ahrens"
-                << endl;
-        geoFile << "ASCII" << endl;
-        geoFile << "DATASET UNSTRUCTURED_GRID" << endl;
+ // write points (the points are the primal-grid cell centers)
-        // write points (the points are the primal-grid cell centers)
+ if (singlelayer) {
+ geoFile << "POINTS " << currentExtraPoint << " float" << endl;
+ } else {
+ geoFile << "POINTS " << currentExtraPoint*(maxNVertLevels+1) << " float" << endl;
+ }
-        geoFile << "POINTS " << currentExtraPoint*(maxNVertLevels+1) << " float" << endl;
+ // write the point at each vertical level, plus one level for last layer
-        // write the point at each vertical level, plus one level for last layer
+ //cout << "Writing points at each level" << endl;
-        //cout << "Writing points at each level" << endl;
+ for (int j = 0; j < currentExtraPoint; j++ )
+ {
+ geoFile.precision(16);
-        for (int j = 0; j < currentExtraPoint; j++ )
-        {
-                geoFile.precision(16);
-                const double PI = 3.141592;
-                xCellData[j] = xCellData[j] * 180.0 / PI;
-                yCellData[j] = yCellData[j] * 180.0 / PI;
-                if (abs(xCellData[j]) < 1e-126) xCellData[j] = 0;
-                if (abs(yCellData[j]) < 1e-126) yCellData[j] = 0;
+ if (abs(xCellData[j]) < 1e-126) xCellData[j] = 0;
+ if (abs(yCellData[j]) < 1e-126) yCellData[j] = 0;
-                for (int levelNum = 0; levelNum < maxNVertLevels+1; levelNum++) {
-                        geoFile << xCellData[j] << "\t" << yCellData[j] << "\t" << -(((float)levelNum)*layerThickness) << endl;
-                }
-        }        
+ float xLon = xCellData[j] * 180.0 / PI;
+ float yLon = yCellData[j] * 180.0 / PI;
-                geoFile << endl;
+ if (singlelayer) {
+ geoFile << xLon << "\t" << yLon << "\t" << 0 << endl;
+ } else {
+ for (int levelNum = 0; levelNum < maxNVertLevels+1; levelNum++) {
+ geoFile << xLon << "\t" << yLon << "\t" << -(((float)levelNum)*layerThickness) << endl;
+ }
+ }
+ }
+ geoFile << endl;
-                // 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.
-                int newDegree = vertexDegree->size()*2;
+ // 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 " << currentExtraCell*maxNVertLevels << " "
-                        << currentExtraCell * (newDegree + 1) * maxNVertLevels << endl;        
+ int degree = multilayer?vertexDegree->size()*2:vertexDegree->size();
-                // for each dual-mesh cell, write number of points for each
-                // and then list the points by number
+ if (singlelayer) {
+ geoFile << "CELLS " << currentExtraCell << " "
+ << currentExtraCell * (degree + 1) << endl;
+ } else {
+ geoFile << "CELLS " << currentExtraCell*maxNVertLevels << " "
+ << currentExtraCell * (degree + 1) * maxNVertLevels << endl;
+ }
-                //cout << "Writing Cells" << endl;
+ // for each dual-mesh cell, write number of points for each
+ // and then list the points by number
-                for (int j = 0; j < currentExtraCell ; j++) {
+ //cout << "Writing Cells" << endl;
-                        // 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 = myCellsOnVertex + (j * vertexDegree->size());
+ for (int j = 0; j < currentExtraCell ; j++) {
-                        // for each level, write the prism
-                        for (int levelNum = 0; levelNum < maxNVertLevels; levelNum++) {
-                                geoFile << newDegree << "\t" ;
-                                for (int k = 0; k < vertexDegree->size(); k++)
-                                {                
-                                        geoFile << (dualCells[k]*(maxNVertLevels+1)) + levelNum << "\t";
-                                }
-                                for (int k = 0; k < vertexDegree->size(); k++)
-                                {                
-                                        geoFile << (dualCells[k]*(maxNVertLevels+1)) + levelNum+1 << "\t";
-                                }
-                                geoFile << endl;
-                        }
-                }        
+ // 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 = myCellsOnVertex + (j * vertexDegree->size());
-                // write cell types
-                int cellType = VTK_WEDGE;
+ // for each level, write the cell
+ if (singlelayer) {
+ geoFile << degree << "\t" ;
+ for (int k = 0; k < vertexDegree->size(); k++)
+ geoFile << dualCells[k] << "\t";
+ geoFile << endl;
+ } else {
+ for (int levelNum = 0; levelNum < maxNVertLevels; levelNum++) {
-                geoFile << "CELL_TYPES " << currentExtraCell*maxNVertLevels << endl;
+ geoFile << degree << "\t" ;
+ for (int k = 0; k < vertexDegree->size(); k++)
+ {
+ geoFile << (dualCells[k]*(maxNVertLevels+1)) + levelNum << "\t";
+ }
+ for (int k = 0; k < vertexDegree->size(); k++)
+ {
+ geoFile << (dualCells[k]*(maxNVertLevels+1)) + levelNum+1 << "\t";
+ }
+ geoFile << endl;
+ }
+ }
+ }
+
+ // write cell types
+ int 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;
+ }
+
+ if (singlelayer) {
+ geoFile << "CELL_TYPES " << currentExtraCell << endl;
+ } else {
+ geoFile << "CELL_TYPES " << currentExtraCell*maxNVertLevels << endl;
+ }
+
//multiply by number of levels
-                for (int j = 0; j < currentExtraCell*maxNVertLevels; j++)
-                {
-                        geoFile << cellType << endl;
-                }
+ if (singlelayer) {
+ for (int j = 0; j < currentExtraCell; j++)
+ {
+ geoFile << cellType << endl;
+ }
+ } else {
+ for (int j = 0; j < currentExtraCell*maxNVertLevels; j++)
+ {
+ geoFile << cellType << endl;
+ }
+ }
-        // release resources
-        geoFile.close();
-        free(xCellData);
-        free(yCellData);
-        free(cellsOnVertex);
-        free(myCellsOnVertex);
+ // release resources
+ geoFile.close();
+ free(xCellData);
+ free(yCellData);
+ free(cellsOnVertex);
+ free(myCellsOnVertex);
-        // For each timestep, write data for each level
+ // For each timestep, write data for each level
-        dualCellVarData = (double*)malloc((sizeof(double)) * currentExtraCell * maxNVertLevels);
-        CHECK_MALLOC(dualCellVarData);
+ dualCellVarData = (double*)malloc((sizeof(double)) * currentExtraCell * maxNVertLevels);
+ CHECK_MALLOC(dualCellVarData);
-        dualPointVarData = (double*)malloc((sizeof(double)) * currentExtraPoint * maxNVertLevels);
-        CHECK_MALLOC(dualPointVarData);
+ dualPointVarData = (double*)malloc((sizeof(double)) * currentExtraPoint * maxNVertLevels);
+ CHECK_MALLOC(dualPointVarData);
-        // for each timestep, copy the geometry file, since we don't want to
-        // have to recompute the points
-        for (int i = 0; i < Time->size(); i++) {
-                ostringstream vtkFileName;
-                vtkFileName << i << argv[2];
+ // for each timestep, copy the geometry file, since we don't want to
+ // have to recompute the points
+ for (int i = 0; i < Time->size(); i++) {
+ ostringstream vtkFileName;
+ vtkFileName << i << outputFileName;
-                string geoFileNameString = geoFileName.str();
-                string vtkFileNameString = vtkFileName.str();
-                int copyRetVal = myCopyFile(&geoFileNameString, &vtkFileNameString);
-                //cout << "myCopyFile returned: " << copyRetVal << endl;
+ string geoFileNameString = geoFileName.str();
+ string vtkFileNameString = vtkFileName.str();
+ int copyRetVal = myCopyFile(&geoFileNameString, &vtkFileNameString);
+ //cout << "myCopyFile returned: " << copyRetVal << endl;
-                ofstream vtkFile(vtkFileName.str().c_str(), ios::out|ios::app);
-                if (!vtkFile)
-                {
-                        cerr << "vtk output file could not be opened" <<endl;
-                        exit(1);
-                }
+ ofstream vtkFile(vtkFileName.str().c_str(), ios::out|ios::app);
+ if (!vtkFile)
+ {
+ cerr << "vtk output file could not be opened" <<endl;
+ exit(1);
+ }
-                if (!vtkFile)
-                {
-                        cerr << "vtk output file could not be opened" <<endl;
-                        exit(1);
-                }
+ vtkFile.precision(16);
-                vtkFile.precision(16);
+ int vertLevelIndex = 0;
+ if (singlelayer) vertLevelIndex = outputVertLevel;
-                // If by point, write out point data
+ // If by point, write out point data
-                if (dualPointVarIndex >= 0) vtkFile << "POINT_DATA " << currentExtraPoint*(maxNVertLevels+1) << endl;
+ if (singlelayer) {
+ if (dualPointVarIndex >= 0) vtkFile << "POINT_DATA " << currentExtraPoint << endl;
+ } else {
+ if (dualPointVarIndex >= 0) vtkFile << "POINT_DATA " << currentExtraPoint*(maxNVertLevels+1) << endl;
+ }
-                int printstep = -1;
+ int printstep = -1;
-                if (i == printstep) cout << "TIME STEP: " << i << endl;
+ if (i == printstep) cout << "TIME STEP: " << i << endl;
-                int tracerNum = 0;
+ int tracerNum = 0;
-                for (int v = 0; v <= dualPointVarIndex; v++) {
+ for (int v = 0; v <= dualPointVarIndex; v++) {
-                        // Read variable number v data for that timestep
+ // Read variable number v data for that timestep
-                        varVertLevels = dualPointVars[v]->get_dim(2)->size();
+ varVertLevels = dualPointVars[v]->get_dim(2)->size();
-                        bool isTracer = false;
+ bool isTracer = false;
-                        if (!strcmp(dualPointVars[v]->name(), "tracers")) {
-                                isTracer = true;
-                        // Uncomment if want to exclude tracers.
-                        //        continue;
-                        }
+ if (!strcmp(dualPointVars[v]->name(), "tracers")) {
+ isTracer = true;
+ // Uncomment if want to exclude tracers.
+ // continue;
+ }
-                        // Write variable number v data for that timestep
-                        vtkFile << "SCALARS " << dualPointVars[v]->name();
-                        if (isTracer) vtkFile << tracerNum+1;
-                        vtkFile << " float 1" << endl;
-                        vtkFile << "LOOKUP_TABLE default" << endl;
+ // Write variable number v data for that timestep
+ vtkFile << "SCALARS " << dualPointVars[v]->name();
+ if (isTracer) vtkFile << tracerNum+1;
+ vtkFile << " float 1" << endl;
+ vtkFile << "LOOKUP_TABLE default" << endl;
+ 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);
+ }
+ }
-                        if (isTracer) {
-                                dualPointVars[v]->set_cur(i, 0, 0, tracerNum);
-                                dualPointVars[v]->get(dualPointVarData+maxNVertLevels, 1, nCells->size(), maxNVertLevels, 1);
-                        } else {
-                                dualPointVars[v]->set_cur(i, 0, 0);
-                                dualPointVars[v]->get(dualPointVarData+maxNVertLevels, 1, nCells->size(), maxNVertLevels);
-                        }
+ float defaultPointVal = 0.0;
-                        float defaultPointVal = 0.0;
+ double *var_target;
+ float validData;
-                        //write dummy
+ //write dummy
+ if (singlelayer) {
+ var_target = dualPointVarData + 1;
+ validData = convertDouble2ValidFloat (*var_target);
-                        double *var_target = dualPointVarData + maxNVertLevels;
-                        float validData;
+ // write dummy
+ vtkFile << validData << endl;
+ } else {
+ //write dummy
+ var_target = dualPointVarData + maxNVertLevels;
+ for (int levelNum = 0; levelNum < maxNVertLevels; levelNum++) {
-                        for (int levelNum = 0; levelNum < maxNVertLevels; levelNum++) {
+ validData = convertDouble2ValidFloat (*var_target);
-                                validData = convertDouble2ValidFloat (*var_target);
+ // write dummy
+ vtkFile << validData << endl;
+ var_target++;
+ }
-                                // write dummy
-                                vtkFile << validData << endl;
- var_target++;
-        
-                        }
+ // write highest level dummy point
+ if (multilayer) vtkFile << validData << endl;
-                        // write highest level dummy point
-                        vtkFile << validData << endl;
+ var_target = dualPointVarData + maxNVertLevels;
+ }
-                        var_target = dualPointVarData + maxNVertLevels;
+ for (int j = 1; j < numDualPoints; j++) {
-                        for (int j = 1; j < numDualPoints; j++) {
+ if (singlelayer) {
+ validData = convertDouble2ValidFloat (*var_target);
+ vtkFile << validData << endl;
+ var_target++;
+ } else {
-                                // write data for one point -- lowest level to highest
-                                for (int levelNum = 0; levelNum < maxNVertLevels; levelNum++) {
-                                        validData = convertDouble2ValidFloat (*var_target);
-                                        vtkFile << validData << endl;
-                                        var_target++;
-                                }
-                        
-                                // for last layer of dual points, repeat last level's values
-                                // Need Mark's input on this one
-                                vtkFile << validData << endl;
-                        }
+ // write data for one point -- lowest level to highest
+ for (int levelNum = 0; levelNum < maxNVertLevels; levelNum++) {
+ validData = convertDouble2ValidFloat (*var_target);
+ vtkFile << validData << endl;
+ var_target++;
+ }
+
+ // for last layer of dual points, repeat last level's values
+ // Need Mark's input on this one
+ vtkFile << validData << endl;
+ }
+ }
-                        // put out data for extra points
-                        for (int j = numDualPoints; j < currentExtraPoint; j++) {
+ // 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);
+ var_target = dualPointVarData + ((*(pointMap + j - numDualPoints))*maxNVertLevels);
-                                // write data for one point -- lowest level to highest
-                                for (int levelNum = 0; levelNum < maxNVertLevels; levelNum++) {
-                                        validData = convertDouble2ValidFloat (*var_target);
-                                        vtkFile << validData << endl;
-                                        var_target++;
-                                }
-                        
-                                // for last layer of dual points, repeat last level's values
-                                // Need Mark's input on this one
-                                vtkFile << validData << endl;
-                        }
-                        if (isTracer) tracerNum++;
-                }
+ if (singlelayer) {
+ validData = convertDouble2ValidFloat (*var_target);
+ vtkFile << validData << endl;
+ var_target++;
+ } else {
+ // write data for one point -- lowest level to highest
+ for (int levelNum = 0; levelNum < maxNVertLevels; levelNum++) {
+ validData = convertDouble2ValidFloat (*var_target);
+ vtkFile << validData << endl;
+ var_target++;
+ }
+
+ // for last layer of dual points, repeat last level's values
+ // Need Mark's input on this one
+ vtkFile << validData << endl;
+ }
+ }
+ if (isTracer) tracerNum++;
+ }
-                // if by cell, then write out cell data
+ // if by cell, then write out cell data
-                if (dualCellVarIndex >= 0) vtkFile << "CELL_DATA " << currentExtraCell*maxNVertLevels << endl;
+ if (singlelayer) {
+ if (dualCellVarIndex >= 0) vtkFile << "CELL_DATA " << currentExtraCell << endl;
+ } else {
+ if (dualCellVarIndex >= 0) vtkFile << "CELL_DATA " << currentExtraCell*maxNVertLevels << endl;
+ }
-                for (int v = 0; v <= dualCellVarIndex; v++) {
+ for (int v = 0; v <= dualCellVarIndex; v++) {
-                        // Write variable number v data for that timestep
-                        vtkFile << "SCALARS " << dualCellVars[v]->name() << " float 1" << endl;
-                        vtkFile << "LOOKUP_TABLE default" << endl;
+ // Write variable number v data for that timestep
+ vtkFile << "SCALARS " << dualCellVars[v]->name() << " float 1" << endl;
+ vtkFile << "LOOKUP_TABLE default" << endl;
-                        // Read variable number v data for that timestep and level
-                        dualCellVars[v]->set_cur(i, 0, 0);
-                        dualCellVars[v]->get(dualCellVarData, 1, numDualCells, maxNVertLevels);
+ // 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);
+ }
double *var_target = dualCellVarData;
-                        for (int j = 0; j < numDualCells; j++) {
+ for (int j = 0; j < numDualCells; j++) {
-                                for (int levelNum = 0; levelNum < maxNVertLevels; levelNum++)
-                                {
-                                        float validData = convertDouble2ValidFloat (*var_target);
-                                        vtkFile << validData << endl;
-                                        var_target++;
-                                }
-                        }
+ if (singlelayer) {
+ float validData = convertDouble2ValidFloat (*var_target);
+ vtkFile << validData << endl;
+ var_target++;
+ } else {
+ for (int levelNum = 0; levelNum < maxNVertLevels; levelNum++)
+ {
+ float validData = convertDouble2ValidFloat (*var_target);
+ vtkFile << validData << endl;
+ var_target++;
+ }
+ }
+ }
-                        for (int j = numDualCells; j < currentExtraCell; j++) {
+ 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;
- for (int levelNum = 0; levelNum < maxNVertLevels; levelNum++)
- {
+ if (singlelayer) {
float validData = convertDouble2ValidFloat (*var_target);
vtkFile << validData << endl;
var_target++;
- }
- }
+ } else {
+ for (int levelNum = 0; levelNum < maxNVertLevels; levelNum++)
+ {
+ float validData = convertDouble2ValidFloat (*var_target);
+ vtkFile << validData << endl;
+ var_target++;
+ }
+ }
+ }
}
}
}
Modified: branches/ocean_projects/graphics/paraview/translator/nc2vtk_dual.cpp
===================================================================
--- branches/ocean_projects/graphics/paraview/translator/nc2vtk_dual.cpp        2010-04-15 22:21:11 UTC (rev 189)
+++ branches/ocean_projects/graphics/paraview/translator/nc2vtk_dual.cpp        2010-04-16 00:07:00 UTC (rev 190)
@@ -17,9 +17,9 @@
// if outside of -FLT_MAX and FLT_MAX make -FLT_MAX or FLT_MAX
//
//
-// Version 1.3
+// Version 1.4
// Christine Ahrens
-// 3/8/2010
+// 4/15/2010
//
////////////////////////////////////////////////////////////////////////////////
@@ -90,7 +90,7 @@
                cerr << "If vertical level is not specified, default is 0." << endl;
                cerr << "A series of vtk files will be created, one file for each time step." << endl;
                cerr << "with time prepended to the file name (e.g. 0outfile.vtk)." << endl;
-                cerr << "vtk datafile Version 2.0, transdual Version 1.2" << endl;
+                cerr << "vtk datafile Version 2.0, transdual Version 1.4" << endl;
                exit(1);
        }
@@ -350,14 +350,29 @@
                        // 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());
+ int* dualCells = cellsOnVertex + (j * vertexDegree->size());
-                        for (int k = 0; k < vertexDegree->size(); k++)
-                        {                
-                                vtkFile << dualCells[k] << "\t";
-                        }
+ // go through and make sure none of the referenced points are
+ // out of range (<=0 or > nCells->size())
-                        vtkFile << endl;
+ 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++)
+ {
+ vtkFile << dualCells[k] << "\t";
+ }
+
+ vtkFile << endl;
                }        
                vtkFile << endl;
</font>
</pre>