<p><b>cahrens@lanl.gov</b> 2010-04-08 13:38:18 -0600 (Thu, 08 Apr 2010)</p><p>Vis tools. paraview dir has code for tools that convert MPAS-style .nc files to ParaView .vtk file formats. Included is a translator that does the conversion for a single layer sphere and also tools to do three projections: a multilayer sphere, a multilayered lat/lon projection of a sphere and a multilayered quad rectangle.<br>
</p><hr noshade><pre><font color="gray">Added: branches/ocean_projects/graphics/paraview/projection/Makefile
===================================================================
--- branches/ocean_projects/graphics/paraview/projection/Makefile         (rev 0)
+++ branches/ocean_projects/graphics/paraview/projection/Makefile        2010-04-08 19:38:18 UTC (rev 186)
@@ -0,0 +1,29 @@
+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/usr/lib -lstdc++ \
+#MYDEBUG = -g
+MYDEBUG =
+
+proj3dsphere: proj3dsphere.o convert.o
+        $(CC) $(MYDEBUG) $(INCLUDES) -o $@ $< convert.o $(LIBS)
+
+proj3dsphere.o: proj3dsphere.cpp
+        $(CC) $(MYDEBUG) $(INCLUDES) -c $<
+
+convert.o: convert.c
+        $(CC) $(MYDEBUG) -c $<
+
+projlatlon: projlatlon.cpp
+        $(CC) $(INCLUDES) -o $@ $< $(LIBS)
+
+proj3dquadrect: proj3dquadrect.cpp
+        $(CC) $(INCLUDES) -o $@ $< $(LIBS)
+
+clean:
+        rm *.o proj3dsphere projlatlon proj3dquadrect
Added: branches/ocean_projects/graphics/paraview/projection/proj3dquadrect.cpp
===================================================================
--- branches/ocean_projects/graphics/paraview/projection/proj3dquadrect.cpp         (rev 0)
+++ branches/ocean_projects/graphics/paraview/projection/proj3dquadrect.cpp        2010-04-08 19:38:18 UTC (rev 186)
@@ -0,0 +1,611 @@
+////////////////////////////////////////////////////////////////////////////////
+//
+// This program translates a newpop netCDF data file to dual-grid lat/lon multilayer
+// projection in legacy, ascii VTK format.
+//
+// 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
+// Does not deal with tracers.
+// 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 Beta
+// Christine Ahrens
+// 3/8/2010
+////////////////////////////////////////////////////////////////////////////////
+
+
+#include <iostream>
+#include <fstream>
+#include <sstream>
+#include "stdlib.h"
+#include "vtkCellType.h"
+#include "netcdfcpp.h"
+#include <string>
+#include <cmath>
+//#include <math>
+#include <cfloat>
+
+using namespace std;
+
+#define CHECK_MALLOC(ptr) \
+        if (ptr == NULL) { \
+                cerr << "malloc failed!</font>
<font color="blue">"; \
+                exit(1); \
+        }
+
+#define MAX_VARS 100
+#define DEFAULT_LAYER_THICKNESS 60000
+#define MY_PRECISION 12
+using namespace std;
+
+#define BUFF_SIZE 2048
+
+// overwrite outFile if it already exists
+int myCopyFile(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::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;
+}
+
+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: projlatlon infile.nc infile.vtk [layer_thickness]" << endl;
+                cerr << "Note: layer_thickness defaults to 10." << 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");
+
+        if ((vertexDegree->size() < 3) || (vertexDegree->size() > 4)) {
+                        cerr << "This code is only for hexagonal-primal/triangular-dual grid or quad-rect grid" << 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();
+
+        //cout << "maxNVertLevels: " << maxNVertLevels << endl;
+
+        float layerThickness = DEFAULT_LAYER_THICKNESS;
+        if (argc == 4) {
+                layerThickness = atof(argv[3]);
+        }
+
+        // figure out what variables to visualize
+        NcVar* dualCellVars[MAX_VARS];
+        NcVar* dualPointVars[MAX_VARS];
+        int dualCellVarIndex = -1;
+        int dualPointVarIndex = -1;
+        int numDualCells = nVertices->size();
+        //cout << "numDualCells: " << numDualCells << endl;
+        int numDualPoints = nCells->size()+1;
+        //cout << "numDualPoints: " << numDualPoints << endl;
+        
+        int numVars = ncFile.num_vars();
+
+        bool tracersExist = 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
+                                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.
+
+        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());
+
+        // get dual-mesh cells
+
+        int *cellsOnVertex = (int *) malloc((nVertices->size()) * vertexDegree->size() *
+                sizeof(int));
+        CHECK_MALLOC(cellsOnVertex);
+        int *myCellsOnVertex = (int *) malloc((nVertices->size()) * 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());
+
+
+        // Go through all cellsOnVertex and "zero out" the cell if the triangle
+        // formed by the cell centers goes outside the lat/lon perimeter.
+        // "zero out" means make all the points in the cell (0,0,0).
+        // This is useful, because when we fetch the data, we don't have to skip
+        // over unused data, it will just be invisible when displayed.
+        // Store new list in myCellsOnVertex
+        
+        int *myptr = myCellsOnVertex;
+
+        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]-1] - xCellData[dualCells[lastk]-1]) > 2000000) xWrap = true;
+                        if (abs(yCellData[dualCells[k]-1] - yCellData[dualCells[lastk]]-1) > 4000000) yWrap = true;
+                        lastk = k;
+                }
+
+                if (xWrap || yWrap) {
+                        //cerr << "Cell wrap: " << j << endl;
+                }
+
+                if (xWrap) {
+                        //cerr << "It wrapped in x 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++;
+                        }
+                } else { // add (0,0,0) point
+                        for (int k=0; k< vertexDegree->size(); k++) {
+                                *myptr = 0;
+                                myptr++;
+                        }
+                }
+        }        
+                                                
+        // decls for data storage
+        double* dualCellVarData;
+        double* dualPointVarData;
+
+        // for each variable, allocate space for variables
+
+        //cout << "dualCellVarIndex: " << dualCellVarIndex << endl;
+
+        int varVertLevels = 0;
+/*
+        dualCellVarData = (double*)malloc((sizeof(double))
+                        * numDualCells * maxNVertLevels);
+        CHECK_MALLOC(dualCellVarData);
+        
+        dualPointVarData = (double*)malloc((sizeof(double))
+                        * nCells->size() * maxNVertLevels);
+        CHECK_MALLOC(dualPointVarData);
+*/
+        
+// write a file with the geometry.
+
+                ostringstream geoFileName;
+
+                geoFileName << "geo_" << argv[2];         
+
+                ofstream geoFile(geoFileName.str().c_str(), ios::out);
+
+                if (!geoFile)
+                {
+                        cerr << "vtk output file could not be opened" <<endl;
+                        exit(1);
+                }
+
+
+                // 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;
+
+
+                // write points (the points are the primal-grid cell centers)
+
+                geoFile << "POINTS " << numDualPoints*(maxNVertLevels+1) << " float" << endl;
+
+        // write the point at each vertical level, plus one level for last layer
+
+        //cout << "Writing dummy points" << endl;
+        for (int levelNum = 0; levelNum < maxNVertLevels+1; levelNum++) {
+                
+                // first write a dummy point, because the climate code
+                // starts their cell numbering at 1 and VTK starts it at
+                // 0
+                geoFile.precision(MY_PRECISION);
+                geoFile << (float)0.0 << "\t" << (float)0.0 << "\t" << (float)0.0
+                        << endl;
+        }
+
+        //cout << "Writing points at each level" << endl;
+
+        for (int j = 0; j < nCells->size(); j++ )
+        {
+                geoFile.precision(MY_PRECISION);
+        /*        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;
+
+                for (int levelNum = 0; levelNum < maxNVertLevels+1; levelNum++) {
+                        geoFile << xCellData[j] << "\t" << yCellData[j] << "\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 = 2*vertexDegree->size();
+
+                geoFile << "CELLS " << numDualCells*maxNVertLevels << " "
+                        << numDualCells * (newDegree + 1) * maxNVertLevels << endl;        
+
+                // for each dual-mesh cell, write number of points for each
+                // and then list the points by number
+
+                //cout << "Writing Cells" << endl;
+
+                for (int j = 0; j < numDualCells ; 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).
+                        
+                        int* dualCells = myCellsOnVertex + (j * vertexDegree->size());
+
+                        // 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;
+                        }
+                }        
+
+                // write cell types
+                int cellType;
+                if (vertexDegree->size()==3) {
+                        cellType = VTK_WEDGE;
+                } else {
+                        cellType = VTK_HEXAHEDRON;
+                }
+
+                geoFile << "CELL_TYPES " << numDualCells*maxNVertLevels << endl;
+
+//multiply by number of levels
+                for (int j = 0; j < numDualCells*maxNVertLevels; j++)
+                {
+                        geoFile << cellType << endl;
+                }
+
+        // release resources
+        geoFile.close();
+        free(xCellData);
+        free(yCellData);
+        free(cellsOnVertex);
+        free(myCellsOnVertex);
+
+        // For each timestep, write data for each level
+
+        dualCellVarData = (double*)malloc((sizeof(double)) * numDualCells * maxNVertLevels);
+        CHECK_MALLOC(dualCellVarData);
+
+        dualPointVarData = (double*)malloc((sizeof(double)) * nCells->size() * 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];
+
+                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);
+                }
+
+                if (!vtkFile)
+                {
+                        cerr << "vtk output file could not be opened" <<endl;
+                        exit(1);
+                }
+
+                vtkFile.precision(MY_PRECISION);
+
+                // If by point, write out point data
+
+                if (dualPointVarIndex >= 0) vtkFile << "POINT_DATA " << numDualPoints*(maxNVertLevels+1) << endl;
+
+                int printstep = -1;
+
+                if (i == printstep) cout << "TIME STEP: " << i << endl;
+
+                int tracerNum = 0;
+
+                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
+                        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, 0, tracerNum);
+                                dualPointVars[v]->get(dualPointVarData, 1, nCells->size(), maxNVertLevels, 1);
+                        } else {
+                                dualPointVars[v]->set_cur(i, 0, 0);
+                                dualPointVars[v]->get(dualPointVarData, 1, nCells->size(), maxNVertLevels);
+                        }
+
+
+                        float defaultPointVal = 0.0;
+
+                        vtkFile.precision(MY_PRECISION);
+                        //write dummy
+
+                        double *var_target = dualPointVarData;
+                        float validData;
+
+                        for (int levelNum = 0; levelNum < maxNVertLevels; levelNum++) {
+
+                                validData = convertDouble2ValidFloat (*var_target);
+
+                                // write dummy
+                                vtkFile << validData << endl;
+        
+                                var_target++;
+                        }
+
+                        // write highest level dummy point
+                        vtkFile << validData << endl;
+
+                        var_target = dualPointVarData;
+
+                        for (int j = 0; j < nCells->size(); j++) {
+
+                                // 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 (dualCellVarIndex >= 0) vtkFile << "CELL_DATA " << numDualCells*maxNVertLevels << endl;
+
+                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;
+
+                        // Read variable number v data for that timestep and level
+                        dualCellVars[v]->set_cur(i, 0, 0);
+                        dualCellVars[v]->get(dualCellVarData, 1, numDualCells, maxNVertLevels);
+
+                        vtkFile.precision(MY_PRECISION);
+                        for (int j = 0; j < numDualCells; j++) {
+
+                                double *var_target = dualCellVarData;
+
+                                for (int levelNum = 0; levelNum < maxNVertLevels; levelNum++)
+                                {
+                                        float validData = convertDouble2ValidFloat (*var_target);
+                                        vtkFile << validData << endl;
+
+                                        var_target++;
+                                }
+                        }
+                }
+        }
+}
+
Property changes on: branches/ocean_projects/graphics/paraview/projection/proj3dquadrect.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/proj3dsphere.cpp
===================================================================
--- branches/ocean_projects/graphics/paraview/projection/proj3dsphere.cpp         (rev 0)
+++ branches/ocean_projects/graphics/paraview/projection/proj3dsphere.cpp        2010-04-08 19:38:18 UTC (rev 186)
@@ -0,0 +1,569 @@
+////////////////////////////////////////////////////////////////////////////////
+//
+// This program translates a newpop netCDF data file to dual-grid legacy, ascii VTK format
+// The variables that have time dim are automatically written.
+// Version Beta
+//
+// 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.
+// 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
+//
+// Christine Ahrens
+// 3/8/2010
+//
+////////////////////////////////////////////////////////////////////////////////
+
+
+#include <iostream>
+#include <fstream>
+#include <sstream>
+#include "stdlib.h"
+#include "vtkCellType.h"
+#include "netcdfcpp.h"
+#include <string>
+#include <cmath>
+#include <cfloat>
+
+using namespace std;
+
+#define CHECK_MALLOC(ptr) \
+        if (ptr == NULL) { \
+                cerr << "malloc failed!</font>
<font color="blue">"; \
+                exit(1); \
+        }
+
+#define MAX_VARS 100
+#define DEFAULT_LAYER_THICKNESS 100000
+
+extern "C" {
+        void CartesianToSpherical(float x, float y, float z, float* rho, float* phi, float* theta);
+        void SphericalToCartesian(float rho, float phi, float theta, float* x, float* y, float* z);
+}        
+
+// START
+#include <fstream>
+#include <iostream>
+#include <string>
+
+using namespace std;
+
+#define BUFF_SIZE 2048
+
+// overwrite outFile if it already exists
+int myCopyFile(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::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;
+}
+
+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: proj3dsphere infile.nc infile.vtk [layer_thickness]" << endl;
+                cerr << "Note: layer_thickness defaults to 100000." << endl;
+                exit(1);
+        }
+
+        //cerr << "FLT_MAX: " << FLT_MAX << " FLT_MIN: " << FLT_MIN << endl;
+
+        NcFile ncFile(argv[1]);
+
+        if (!ncFile.is_valid())
+        {
+                cerr << "Couldn't open file: " << argv[1] << endl;
+                exit(1);
+        }
+
+        int layerThickness = DEFAULT_LAYER_THICKNESS;
+        if (argc == 4) {
+                layerThickness = atoi(argv[3]);
+        }
+
+        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();
+
+        //cout << "maxNVertLevels: " << maxNVertLevels << endl;
+
+        // 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 tracersExist = 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
+                                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.
+
+        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));
+        //cout << "ptr for zCellData" << zCellData << endl;
+        CHECK_MALLOC(zCellData);
+        NcVar *zCellVar = ncFile.get_var("zCell");
+        zCellVar->get(zCellData, nCells->size());
+
+        // get dual-mesh cells
+
+        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="gray">";
+        cellsOnVertexVar->get(cellsOnVertex, nVertices->size(), vertexDegree->size());
+
+        // decls for data storage
+        double* dualCellVarData;
+        double* dualPointVarData;
+
+        // for each variable, allocate space for variables
+
+        //cout << "dualCellVarIndex: " << dualCellVarIndex << endl;
+
+        int varVertLevels = 0;
+
+
+        // write a file with the geometry.
+
+        ostringstream geoFileName;
+
+        geoFileName << "geo_" << argv[2];         
+
+        ofstream geoFile(geoFileName.str().c_str(), ios::out);
+
+        if (!geoFile)
+        {
+                cerr << "vtk output file could not be opened" <<endl;
+                exit(1);
+        }
+
+
+        // write header
+
+        geoFile << "# vtk DataFile Version 2.0" << endl;
+        geoFile << "Translated newpop geometry to dual grid on 3D sphere from netCDF by Christine Ahrens"
+                << endl;
+        geoFile << "ASCII" << endl;
+        geoFile << "DATASET UNSTRUCTURED_GRID" << endl;
+
+
+        // write points (the points are the primal-grid cell centers)
+
+        geoFile << "POINTS " << numDualPoints*(maxNVertLevels+1) << " float" << endl;
+
+        // write the point at each vertical level, plus one level for last layer
+
+        //cout << "Writing dummy points" << endl;
+        for (int levelNum = 0; levelNum < maxNVertLevels+1; levelNum++) {
+
+                // first write a dummy point, because the climate code
+                // starts their cell numbering at 1 and VTK starts it at
+                // 0
+                geoFile.precision(16);
+                geoFile << (float)0.0 << "\t" << (float)0.0 << "\t" << (float)0.0
+                        << endl;
+        }
+
+        //cout << "Writing points at each level" << 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;
+
+                float rho, rholevel, theta, phi, x, y, z;
+
+                CartesianToSpherical((float)xCellData[j], (float)yCellData[j], (float)zCellData[j],
+                                &rho, &phi, &theta);
+
+                for (int levelNum = 0; levelNum < maxNVertLevels+1; levelNum++) {
+
+                        // decrease rho, so we go down (inwards towards the center of the sphere)
+                        rholevel = rho - (layerThickness * levelNum);
+
+                        SphericalToCartesian(rholevel, phi, theta, &x, &y, &z);
+
+                        geoFile << x << "\t" << y << "\t" << z << 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 = 6;
+
+        geoFile << "CELLS " << numDualCells*maxNVertLevels << " "
+                << numDualCells * (newDegree + 1) * maxNVertLevels << endl;        
+
+        // for each dual-mesh cell, write number of points for each
+        // and then list the points by number
+
+        //cout << "Writing Cells" << endl;
+
+        for (int j = 0; j < numDualCells ; 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).
+
+                int* dualCells = cellsOnVertex + (j * vertexDegree->size());
+
+                // 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;
+                }
+        }        
+        geoFile << endl;
+
+
+        // write cell types
+        int cellType;
+
+        if (vertexDegree->size() == 3)
+                cellType = VTK_WEDGE;
+
+        geoFile << "CELL_TYPES " << numDualCells*maxNVertLevels << endl;
+
+        //multiply by number of levels
+        for (int j = 0; j < numDualCells*maxNVertLevels; j++)
+        {
+                geoFile << cellType << endl;
+        }
+        geoFile << endl;
+
+
+        // release resources
+        geoFile.close();
+        free(xCellData);
+        free(yCellData);
+        free(zCellData);
+        free(cellsOnVertex);
+
+        // For each timestep, write data for each level
+
+        dualCellVarData = (double*)malloc((sizeof(double)) * numDualCells * maxNVertLevels);
+        CHECK_MALLOC(dualCellVarData);
+
+        dualPointVarData = (double*)malloc((sizeof(double)) * nCells->size() * 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];
+
+                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);
+                }
+
+                if (!vtkFile)
+                {
+                        cerr << "vtk output file could not be opened" <<endl;
+                        exit(1);
+                }
+
+                vtkFile.precision(16);
+
+                // If by point, write out point data
+
+                if (dualPointVarIndex >= 0) vtkFile << "POINT_DATA " << numDualPoints*(maxNVertLevels+1) << endl;
+
+                int printstep = -1;
+
+                if (i == printstep) cout << "TIME STEP: " << i << endl;
+
+                int tracerNum = 0;
+
+                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
+                        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, 0, tracerNum);
+                                dualPointVars[v]->get(dualPointVarData, 1, nCells->size(), maxNVertLevels, 1);
+                        } else {
+                                dualPointVars[v]->set_cur(i, 0, 0);
+                                dualPointVars[v]->get(dualPointVarData, 1, nCells->size(), maxNVertLevels);
+                        }
+
+
+                        float defaultPointVal = 0.0;
+
+                        //write dummy
+
+                        double *var_target = dualPointVarData;
+                        float validData;
+
+                        for (int levelNum = 0; levelNum < maxNVertLevels; levelNum++) {
+
+                                validData = convertDouble2ValidFloat (*var_target);
+
+                                // write dummy
+                                vtkFile << validData << endl;
+        
+                                var_target++;
+                        }
+
+                        // write highest level dummy point
+                        vtkFile << validData << endl;
+
+                        var_target = dualPointVarData;
+
+                        for (int j = 0; j < nCells->size(); j++) {
+
+                                // 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 (dualCellVarIndex >= 0) vtkFile << "CELL_DATA " << numDualCells*maxNVertLevels << endl;
+
+                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;
+
+                        // Read variable number v data for that timestep and level
+                        dualCellVars[v]->set_cur(i, 0, 0);
+                        dualCellVars[v]->get(dualCellVarData, 1, numDualCells, maxNVertLevels);
+
+                        for (int j = 0; j < numDualCells; j++) {
+
+                                double *var_target = dualCellVarData;
+
+                                for (int levelNum = 0; levelNum < maxNVertLevels; levelNum++)
+                                {
+                                        float validData = convertDouble2ValidFloat (*var_target);
+                                        vtkFile << validData << endl;
+
+                                        var_target++;
+                                }
+                        }
+                }
+        }
+}
Property changes on: branches/ocean_projects/graphics/paraview/projection/proj3dsphere.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/projlatlon.cpp
===================================================================
--- branches/ocean_projects/graphics/paraview/projection/projlatlon.cpp         (rev 0)
+++ branches/ocean_projects/graphics/paraview/projection/projlatlon.cpp        2010-04-08 19:38:18 UTC (rev 186)
@@ -0,0 +1,779 @@
+////////////////////////////////////////////////////////////////////////////////
+//
+// This program translates a newpop netCDF data file to dual-grid lat/lon multilayer
+// projection in legacy, ascii VTK format.
+//
+// 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
+// Does not deal with tracers.
+// 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 Beta
+// Christine Ahrens
+// 3/8/2010
+////////////////////////////////////////////////////////////////////////////////
+
+
+#include <iostream>
+#include <fstream>
+#include <sstream>
+#include "stdlib.h"
+#include "vtkCellType.h"
+#include "netcdfcpp.h"
+#include <string>
+#include <cmath>
+//#include <math>
+#include <cfloat>
+
+using namespace std;
+
+#define CHECK_MALLOC(ptr) \
+        if (ptr == NULL) { \
+                cerr << "malloc failed!</font>
<font color="blue">"; \
+                exit(1); \
+        }
+
+#define MAX_VARS 100
+#define DEFAULT_LAYER_THICKNESS 10
+using namespace std;
+
+#define BUFF_SIZE 2048
+
+// overwrite outFile if it already exists
+int myCopyFile(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::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;
+}
+
+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: projlatlon infile.nc infile.vtk [layer_thickness]" << endl;
+                cerr << "Note: layer_thickness defaults to 10." << 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");
+
+        if (vertexDegree->size() != 3) {
+                cerr << "This code is only for hexagonal-primal/triangular-dual grid" << 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();
+
+        //cout << "maxNVertLevels: " << maxNVertLevels << endl;
+
+        int outputVertLevel = 0;
+
+        float layerThickness = DEFAULT_LAYER_THICKNESS;
+        if (argc == 4) {
+                layerThickness = atof(argv[3]);
+        }
+
+        // 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 tracersExist = 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
+                                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());
+ // 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());
+ // point 0 is 0.0
+ *yCellData = 0.0;
+
+        // get dual-mesh cells
+
+        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);
+
+        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);
+
+        
+        int *myptr = myCellsOnVertex;
+        int *myExtraX;
+        int *myExtraY;
+        int currentExtraPoint = numDualPoints;
+        int currentExtraCell = numDualCells;
+ const double PI = 3.141592;
+
+ // 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());
+                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
+
+                        // 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 > PI) {
+                                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;
+                }
+        }        
+                                                
+        // decls for data storage
+        double* dualCellVarData;
+        double* dualPointVarData;
+
+        // for each variable, allocate space for variables
+
+        //cout << "dualCellVarIndex: " << dualCellVarIndex << endl;
+
+        int varVertLevels = 0;
+        
+        // write a file with the geometry.
+
+        ostringstream geoFileName;
+
+        geoFileName << "geo_" << argv[2];         
+
+        ofstream geoFile(geoFileName.str().c_str(), ios::out);
+
+        if (!geoFile)
+        {
+                cerr << "vtk output file could not be opened" <<endl;
+                exit(1);
+        }
+
+
+        // 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;
+
+
+        // write points (the points are the primal-grid cell centers)
+
+        geoFile << "POINTS " << currentExtraPoint*(maxNVertLevels+1) << " float" << endl;
+
+        // write the point at each vertical level, plus one level for last layer
+
+        //cout << "Writing points at each level" << endl;
+
+        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;
+
+                for (int levelNum = 0; levelNum < maxNVertLevels+1; levelNum++) {
+                        geoFile << xCellData[j] << "\t" << yCellData[j] << "\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;
+
+                geoFile << "CELLS " << currentExtraCell*maxNVertLevels << " "
+                        << currentExtraCell * (newDegree + 1) * maxNVertLevels << endl;        
+
+                // for each dual-mesh cell, write number of points for each
+                // and then list the points by number
+
+                //cout << "Writing Cells" << endl;
+
+                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).
+                        
+                        int* dualCells = myCellsOnVertex + (j * vertexDegree->size());
+
+                        // 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;
+                        }
+                }        
+
+                // write cell types
+                int cellType = VTK_WEDGE;
+
+                geoFile << "CELL_TYPES " << currentExtraCell*maxNVertLevels << endl;
+
+//multiply by number of levels
+                for (int j = 0; j < currentExtraCell*maxNVertLevels; j++)
+                {
+                        geoFile << cellType << endl;
+                }
+
+        // release resources
+        geoFile.close();
+        free(xCellData);
+        free(yCellData);
+        free(cellsOnVertex);
+        free(myCellsOnVertex);
+
+        // For each timestep, write data for each level
+
+        dualCellVarData = (double*)malloc((sizeof(double)) * currentExtraCell * maxNVertLevels);
+        CHECK_MALLOC(dualCellVarData);
+
+        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];
+
+                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);
+                }
+
+                if (!vtkFile)
+                {
+                        cerr << "vtk output file could not be opened" <<endl;
+                        exit(1);
+                }
+
+                vtkFile.precision(16);
+
+                // If by point, write out point data
+
+                if (dualPointVarIndex >= 0) vtkFile << "POINT_DATA " << currentExtraPoint*(maxNVertLevels+1) << endl;
+
+                int printstep = -1;
+
+                if (i == printstep) cout << "TIME STEP: " << i << endl;
+
+                int tracerNum = 0;
+
+                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
+                        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, 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;
+
+                        //write dummy
+
+                        double *var_target = dualPointVarData + maxNVertLevels;
+                        float validData;
+
+                        for (int levelNum = 0; levelNum < maxNVertLevels; levelNum++) {
+
+                                validData = convertDouble2ValidFloat (*var_target);
+
+                                // write dummy
+                                vtkFile << validData << endl;
+ var_target++;
+        
+                        }
+
+                        // write highest level dummy point
+                        vtkFile << validData << endl;
+
+                        var_target = dualPointVarData + maxNVertLevels;
+
+                        for (int j = 1; j < numDualPoints; j++) {
+
+                                // 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++) {
+ 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);
+
+                                // 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 (dualCellVarIndex >= 0) vtkFile << "CELL_DATA " << currentExtraCell*maxNVertLevels << endl;
+
+                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;
+
+                        // Read variable number v data for that timestep and level
+                        dualCellVars[v]->set_cur(i, 0, 0);
+                        dualCellVars[v]->get(dualCellVarData, 1, numDualCells, maxNVertLevels);
+
+ double *var_target = dualCellVarData;
+
+                        for (int j = 0; j < numDualCells; j++) {
+
+                                for (int levelNum = 0; levelNum < maxNVertLevels; levelNum++)
+                                {
+                                        float validData = convertDouble2ValidFloat (*var_target);
+                                        vtkFile << validData << endl;
+                                        var_target++;
+                                }
+                        }
+
+                        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++)
+ {
+ float validData = convertDouble2ValidFloat (*var_target);
+ vtkFile << validData << endl;
+ var_target++;
+ }
+ }
+ }
+ }
+}
+
Property changes on: branches/ocean_projects/graphics/paraview/projection/projlatlon.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/vtkCellType.h
===================================================================
--- branches/ocean_projects/graphics/paraview/projection/vtkCellType.h         (rev 0)
+++ branches/ocean_projects/graphics/paraview/projection/vtkCellType.h        2010-04-08 19:38:18 UTC (rev 186)
@@ -0,0 +1,98 @@
+/*=========================================================================
+
+ Program: Visualization Toolkit
+ Module: $RCSfile: vtkCellType.h,v $
+
+ Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
+ All rights reserved.
+ See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
+
+ This software is distributed WITHOUT ANY WARRANTY; without even
+ the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ PURPOSE. See the above copyright notice for more information.
+
+=========================================================================*/
+// .NAME vtkCellType - define types of cells
+// .SECTION Description
+// vtkCellType defines the allowable cell types in the visualization
+// library (vtk). In vtk, datasets consist of collections of cells.
+// Different datasets consist of different cell types. The cells may be
+// explicitly represented (as in vtkPolyData), or may be implicit to the
+// data type (as in vtkStructuredPoints).
+
+#ifndef __vtkCellType_h
+#define __vtkCellType_h
+
+// To add a new cell type, define a new integer type flag here, then
+// create a subclass of vtkCell to implement the proper behavior. You
+// may have to modify the following methods: vtkDataSet (and subclasses)
+// GetCell() and vtkGenericCell::SetCellType(). Also, to do the job right,
+// you'll also have to modify some filters (vtkGeometryFilter...) and
+// regression tests (example scripts) to reflect the new cell addition.
+// Also, make sure to update vtkCellTypesStrings in vtkCellTypes.cxx.
+
+// .SECTION Caveats
+// An unstructured grid stores the types of its cells as a
+// unsigned char array. Therefore, the maximum encoding number for a cell type
+// is 255.
+
+typedef enum {
+ // Linear cells
+ VTK_EMPTY_CELL = 0,
+ VTK_VERTEX = 1,
+ VTK_POLY_VERTEX = 2,
+ VTK_LINE = 3,
+ VTK_POLY_LINE = 4,
+ VTK_TRIANGLE = 5,
+ VTK_TRIANGLE_STRIP = 6,
+ VTK_POLYGON = 7,
+ VTK_PIXEL = 8,
+ VTK_QUAD = 9,
+ VTK_TETRA = 10,
+ VTK_VOXEL = 11,
+ VTK_HEXAHEDRON = 12,
+ VTK_WEDGE = 13,
+ VTK_PYRAMID = 14,
+ VTK_PENTAGONAL_PRISM = 15,
+ VTK_HEXAGONAL_PRISM = 16,
+
+ // Quadratic, isoparametric cells
+ VTK_QUADRATIC_EDGE = 21,
+ VTK_QUADRATIC_TRIANGLE = 22,
+ VTK_QUADRATIC_QUAD = 23,
+ VTK_QUADRATIC_TETRA = 24,
+ VTK_QUADRATIC_HEXAHEDRON = 25,
+ VTK_QUADRATIC_WEDGE = 26,
+ VTK_QUADRATIC_PYRAMID = 27,
+ VTK_BIQUADRATIC_QUAD = 28,
+ VTK_TRIQUADRATIC_HEXAHEDRON = 29,
+ VTK_QUADRATIC_LINEAR_QUAD = 30,
+ VTK_QUADRATIC_LINEAR_WEDGE = 31,
+ VTK_BIQUADRATIC_QUADRATIC_WEDGE = 32,
+ VTK_BIQUADRATIC_QUADRATIC_HEXAHEDRON = 33,
+
+ // Special class of cells formed by convex group of points
+ VTK_CONVEX_POINT_SET = 41,
+
+ // Higher order cells in parametric form
+ VTK_PARAMETRIC_CURVE = 51,
+ VTK_PARAMETRIC_SURFACE = 52,
+ VTK_PARAMETRIC_TRI_SURFACE = 53,
+ VTK_PARAMETRIC_QUAD_SURFACE = 54,
+ VTK_PARAMETRIC_TETRA_REGION = 55,
+ VTK_PARAMETRIC_HEX_REGION = 56,
+
+ // Higher order cells
+ VTK_HIGHER_ORDER_EDGE = 60,
+ VTK_HIGHER_ORDER_TRIANGLE = 61,
+ VTK_HIGHER_ORDER_QUAD = 62,
+ VTK_HIGHER_ORDER_POLYGON = 63,
+ VTK_HIGHER_ORDER_TETRAHEDRON = 64,
+ VTK_HIGHER_ORDER_WEDGE = 65,
+ VTK_HIGHER_ORDER_PYRAMID = 66,
+ VTK_HIGHER_ORDER_HEXAHEDRON = 67,
+
+ VTK_NUMBER_OF_CELL_TYPES
+} VTKCellType;
+
+#endif
Property changes on: branches/ocean_projects/graphics/paraview/projection/vtkCellType.h
___________________________________________________________________
Name: svn:mime-type
+ text/x-c++hdr
Name: svn:keywords
+ Author Date ID Revision
Name: svn:eol-style
+ native
Added: branches/ocean_projects/graphics/paraview/translator/Makefile
===================================================================
--- branches/ocean_projects/graphics/paraview/translator/Makefile         (rev 0)
+++ branches/ocean_projects/graphics/paraview/translator/Makefile        2010-04-08 19:38:18 UTC (rev 186)
@@ -0,0 +1,18 @@
+# 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
+#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 \
+        -lstdc++ \
+        -L/opt/Intel/cce/10.0.023/lib/ -lirc
+#         -L/opt/Intel/cce/9.1.039/lib/ -lirc
+
+transdual: nc2vtk_dual.cpp
+        $(CC) $(INCLUDES) -o $@ $< $(LIBS)
+
+clean:
+        rm transdual
Added: branches/ocean_projects/graphics/paraview/translator/nc2vtk_dual.cpp
===================================================================
--- branches/ocean_projects/graphics/paraview/translator/nc2vtk_dual.cpp         (rev 0)
+++ branches/ocean_projects/graphics/paraview/translator/nc2vtk_dual.cpp        2010-04-08 19:38:18 UTC (rev 186)
@@ -0,0 +1,490 @@
+////////////////////////////////////////////////////////////////////////////////
+//
+// This program translates a newpop netCDF data file to dual-grid legacy, ascii VTK 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.
+// 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
+// 3/8/2010
+//
+////////////////////////////////////////////////////////////////////////////////
+
+
+#include <iostream>
+#include <fstream>
+#include <sstream>
+#include "stdlib.h"
+#include "vtkCellType.h"
+#include "netcdfcpp.h"
+#include <string>
+#include <cmath>
+#include <cfloat>
+
+using namespace std;
+
+#define CHECK_MALLOC(ptr) \
+        if (ptr == NULL) { \
+                cerr << "malloc failed!</font>
<font color="blue">"; \
+                exit(1); \
+        }
+#define MAX_VARS 100
+
+int my_isinf_f(float x){
+         if (isinf(x)) {
+                return (x < 0.0 ? -1 : 1);
+        } else return 0;
+}
+
+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: transdual infile.nc outfile.vtk [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 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;
+                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");
+        NcDim* nTracers = ncFile.get_dim("nTracers");
+
+        // 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]);
+         // cout << "outputVertLevel: " << outputVertLevel << endl;
+                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];
+        int dualCellVarIndex = -1;
+        int dualPointVarIndex = -1;
+        int numDualCells = nVertices->size();
+        int numDualPoints = nCells->size()+1;
+        
+        int numVars = ncFile.num_vars();
+
+        bool tracersExist = 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;
+                } 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;
+                                //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.
+
+        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));
+        //cout << "ptr for zCellData" << zCellData << endl;
+        CHECK_MALLOC(zCellData);
+        NcVar *zCellVar = ncFile.get_var("zCell");
+        zCellVar->get(zCellData, nCells->size());
+
+        // get dual-mesh cells
+
+        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="gray">";
+        cellsOnVertexVar->get(cellsOnVertex, nVertices->size(), vertexDegree->size());
+
+        // decls for data storage
+        double* dualCellVarData;
+        double* dualPointVarData;
+
+        // for each variable, allocate space for variables
+
+        //cout << "dualCellVarIndex: " << dualCellVarIndex << endl;
+
+        int varVertLevels = 0;
+
+        dualCellVarData = (double*)malloc((sizeof(double)) * numDualCells);
+        CHECK_MALLOC(dualCellVarData);
+
+        dualPointVarData = (double*)malloc((sizeof(double)) * nCells->size());
+        CHECK_MALLOC(dualPointVarData);
+        
+
+        // 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);
+
+                if (!vtkFile)
+                {
+                        cerr << "vtk output file could not be opened" <<endl;
+                        exit(1);
+                }
+
+
+                // write header
+
+                vtkFile << "# vtk DataFile Version 2.0, transdual Version 1.2" << endl;
+                vtkFile << "Translated newpop data to dual grid for timestep " << i << " from netCDF by Christine Ahrens"
+                        << endl;
+                vtkFile << "ASCII" << endl;
+                vtkFile << "DATASET UNSTRUCTURED_GRID" << endl;
+
+
+                // write points (the points are the primal-grid cell centers)
+
+                vtkFile << "POINTS " << numDualPoints << " float" << endl;
+
+                // first write a dummy point, because the climate code
+                // starts their cell numbering at 1 and VTK starts it at
+                // 0
+                vtkFile.precision(16);
+                vtkFile << (float)0.0 << "\t" << (float)0.0 << "\t" << (float)0.0
+                        << endl;
+
+                for (int j = 0; j < nCells->size(); j++ )
+                {
+                        vtkFile.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;
+                        vtkFile << (float)xCellData[j] << "\t" << (float)yCellData[j] << "\t"
+                                << (float)zCellData[j] << endl;
+                }        
+                vtkFile << 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.
+
+                vtkFile << "CELLS " << numDualCells << " "
+                        << numDualCells * (vertexDegree->size() + 1) << endl;        
+
+                // for each dual-mesh cell, write number of points for each
+                // and then list the points by number
+
+                for (int j = 0; j < numDualCells ; j++) {
+
+                        vtkFile << vertexDegree->size() << "\t" ;
+
+                        // 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());
+
+                        for (int k = 0; k < vertexDegree->size(); k++)
+                        {                
+                                vtkFile << dualCells[k] << "\t";
+                        }
+
+                        vtkFile << endl;
+                }        
+                vtkFile << 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;
+
+                vtkFile << "CELL_TYPES " << numDualCells << endl;
+
+                for (int j = 0; j < numDualCells; j++)
+                {
+                        vtkFile << cellType << endl;
+                }
+                vtkFile << endl;
+
+                // Write attributes of dual-mesh cell data (attributes of primal-mesh
+                // vertex data)
+
+                // for each var, figure out if it is by cell or by point
+
+
+                vtkFile.precision(16);
+
+                // If by point, write out point data
+
+                if (dualPointVarIndex >= 0) vtkFile << "POINT_DATA " << numDualPoints << endl;
+
+                int printstep = -1;
+
+                if (i == printstep) cout << "TIME STEP: " << i << endl;
+
+                int tracerNum = 0;
+
+                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);
+                        }
+
+                        // 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;
+
+                        //debugging
+                        if (i == printstep) cout << "SCALARS " << dualPointVars[v]->name() << " float 1" << endl;
+                        if (i==printstep) {
+                                for (int z = 0; z < varVertLevels; z++) {
+                                        cout << z << ":" << *(dualPointVarData + z) << endl;
+                                }
+                        }
+
+                        // get starting point
+                        double *var_target = dualPointVarData;
+
+                        float validData;
+
+                        validData = convertDouble2ValidFloat (*var_target);
+                        
+                        // write dummy
+                        vtkFile << validData << endl;
+                
+                        // write data        
+                        for (int j = 0; j < nCells->size(); j++)
+                        {
+                                validData = convertDouble2ValidFloat (*var_target);
+                                vtkFile << validData << endl;
+                                var_target++;
+                        }
+
+                        if (isTracer) tracerNum++;
+                }
+
+                // if by cell, then write out cell data
+
+                if (dualCellVarIndex >= 0) vtkFile << "CELL_DATA " << numDualCells << endl;
+
+                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 << "SCALARS " << dualCellVars[v]->name() << " float 1" << endl;
+                        vtkFile << "LOOKUP_TABLE default" << endl;
+
+                        // debugging
+                        if (i==printstep) cout << "SCALARS " << dualCellVars[v]->name() << " float 1" << endl;
+                        if (i==printstep) {
+                                for (int z = 0; z < varVertLevels; z++) {
+                                        cout << z << ":" << *(dualCellVarData + z) << endl;
+                                }
+                        }
+
+                        double *var_target = dualCellVarData;
+                        float validData;
+
+                        for (int j = 0; j < numDualCells; j++)
+                        {
+                                validData = convertDouble2ValidFloat (*var_target);
+                                vtkFile << validData << endl;
+                                var_target++;
+                        }
+                }
+
+        }
+
+        return(0);
+
+}        
+
Property changes on: branches/ocean_projects/graphics/paraview/translator/nc2vtk_dual.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/vtkCellType.h
===================================================================
--- branches/ocean_projects/graphics/paraview/translator/vtkCellType.h         (rev 0)
+++ branches/ocean_projects/graphics/paraview/translator/vtkCellType.h        2010-04-08 19:38:18 UTC (rev 186)
@@ -0,0 +1,98 @@
+/*=========================================================================
+
+ Program: Visualization Toolkit
+ Module: $RCSfile: vtkCellType.h,v $
+
+ Copyright (c) Ken Martin, Will Schroeder, Bill Lorensen
+ All rights reserved.
+ See Copyright.txt or http://www.kitware.com/Copyright.htm for details.
+
+ This software is distributed WITHOUT ANY WARRANTY; without even
+ the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
+ PURPOSE. See the above copyright notice for more information.
+
+=========================================================================*/
+// .NAME vtkCellType - define types of cells
+// .SECTION Description
+// vtkCellType defines the allowable cell types in the visualization
+// library (vtk). In vtk, datasets consist of collections of cells.
+// Different datasets consist of different cell types. The cells may be
+// explicitly represented (as in vtkPolyData), or may be implicit to the
+// data type (as in vtkStructuredPoints).
+
+#ifndef __vtkCellType_h
+#define __vtkCellType_h
+
+// To add a new cell type, define a new integer type flag here, then
+// create a subclass of vtkCell to implement the proper behavior. You
+// may have to modify the following methods: vtkDataSet (and subclasses)
+// GetCell() and vtkGenericCell::SetCellType(). Also, to do the job right,
+// you'll also have to modify some filters (vtkGeometryFilter...) and
+// regression tests (example scripts) to reflect the new cell addition.
+// Also, make sure to update vtkCellTypesStrings in vtkCellTypes.cxx.
+
+// .SECTION Caveats
+// An unstructured grid stores the types of its cells as a
+// unsigned char array. Therefore, the maximum encoding number for a cell type
+// is 255.
+
+typedef enum {
+ // Linear cells
+ VTK_EMPTY_CELL = 0,
+ VTK_VERTEX = 1,
+ VTK_POLY_VERTEX = 2,
+ VTK_LINE = 3,
+ VTK_POLY_LINE = 4,
+ VTK_TRIANGLE = 5,
+ VTK_TRIANGLE_STRIP = 6,
+ VTK_POLYGON = 7,
+ VTK_PIXEL = 8,
+ VTK_QUAD = 9,
+ VTK_TETRA = 10,
+ VTK_VOXEL = 11,
+ VTK_HEXAHEDRON = 12,
+ VTK_WEDGE = 13,
+ VTK_PYRAMID = 14,
+ VTK_PENTAGONAL_PRISM = 15,
+ VTK_HEXAGONAL_PRISM = 16,
+
+ // Quadratic, isoparametric cells
+ VTK_QUADRATIC_EDGE = 21,
+ VTK_QUADRATIC_TRIANGLE = 22,
+ VTK_QUADRATIC_QUAD = 23,
+ VTK_QUADRATIC_TETRA = 24,
+ VTK_QUADRATIC_HEXAHEDRON = 25,
+ VTK_QUADRATIC_WEDGE = 26,
+ VTK_QUADRATIC_PYRAMID = 27,
+ VTK_BIQUADRATIC_QUAD = 28,
+ VTK_TRIQUADRATIC_HEXAHEDRON = 29,
+ VTK_QUADRATIC_LINEAR_QUAD = 30,
+ VTK_QUADRATIC_LINEAR_WEDGE = 31,
+ VTK_BIQUADRATIC_QUADRATIC_WEDGE = 32,
+ VTK_BIQUADRATIC_QUADRATIC_HEXAHEDRON = 33,
+
+ // Special class of cells formed by convex group of points
+ VTK_CONVEX_POINT_SET = 41,
+
+ // Higher order cells in parametric form
+ VTK_PARAMETRIC_CURVE = 51,
+ VTK_PARAMETRIC_SURFACE = 52,
+ VTK_PARAMETRIC_TRI_SURFACE = 53,
+ VTK_PARAMETRIC_QUAD_SURFACE = 54,
+ VTK_PARAMETRIC_TETRA_REGION = 55,
+ VTK_PARAMETRIC_HEX_REGION = 56,
+
+ // Higher order cells
+ VTK_HIGHER_ORDER_EDGE = 60,
+ VTK_HIGHER_ORDER_TRIANGLE = 61,
+ VTK_HIGHER_ORDER_QUAD = 62,
+ VTK_HIGHER_ORDER_POLYGON = 63,
+ VTK_HIGHER_ORDER_TETRAHEDRON = 64,
+ VTK_HIGHER_ORDER_WEDGE = 65,
+ VTK_HIGHER_ORDER_PYRAMID = 66,
+ VTK_HIGHER_ORDER_HEXAHEDRON = 67,
+
+ VTK_NUMBER_OF_CELL_TYPES
+} VTKCellType;
+
+#endif
Property changes on: branches/ocean_projects/graphics/paraview/translator/vtkCellType.h
___________________________________________________________________
Name: svn:mime-type
+ text/x-c++hdr
Name: svn:keywords
+ Author Date ID Revision
Name: svn:eol-style
+ native
</font>
</pre>