#
libsff

Home | Modules | Examples

# readexample.cpp

A a self-contained example that reads the mesh and post data.

    #include <stdlib.h>     // For exit()
    #include <ccmio.h>
    #include <vector>
    #include <map>
    #include <iostream>

    typedef std::vector< unsigned int > UIntArray;
    typedef std::map< int, std::vector< int > > FaceMap;

    // If you are unfamiliar with the Standard Template Library (STL), please read.
    //   vector:    A dynamic array.  vector<int> is an array of ints.  Suppose
    //              we have vector<int> iarr;
    //                  iarr.push_back( ... ); adds an item to the end of the
    //                                         array, resizing if necessary.
    //                  iarr.size();           Returns the number of items.
    //                  iarr[10];              The 11th item in the array.
    //                                         No bounds checking is performed.
    //                  iarr.resize(10);       Makes the array 10 items long.
    //   map:       An associative container.   map< int, float > is an sorted
    //              based on an integer key and contains a value float.
    //              Suppose we have map<int, float> m;
    //                  m[18] = 3.141592;      Associates Pi with the key 18.
    //                                         Technically m[18] returns a reference
    //                                         to the value associated with 18,
    //                                         which is then assigned to.
    //                                         (A reference is like a pointer that
    //                                         cannot change where it is pointed to
    //                                         and is treated syntactically like a
    //                                         normal variable).
    //   iterators: The way to increment through the container.
    //                  container.begin()      Returns an iterator to the first item
    //                  container.end()        Returns an iterator just past the end
    //                  *iter                  The item being iterated to.

    enum DataType { kScalar, kVector, kVertex, kCell, kInternalFace, kBoundaryFace,
                    kBoundaryData, kBoundaryFaceData, kCellType };

    using namespace std;

    static const int kNValues = 10; // Number of values of each element to print
    static const char kDefaultState[] = "default";
    static const char kUnitsName[] = "Units";
    static const int kVertOffset = 2;
    static const int kCellInc = 4;

    static void ReadMesh( CCMIOError& err, CCMIOID& vertices, CCMIOID& topology );
    static void ReadVertices( CCMIOError& err, CCMIOID& vertices,
                              const char *counter, int offset = 0 );
    static void ReadSolverInfo( CCMIOError& err, CCMIOID& solution);
    static void ReadPost( CCMIOError& err, CCMIOID& solution );
    static void ReadScalar( CCMIOError& err, CCMIOID field, vector<int>& mapData,
                            vector<float>& data, bool readingVector = false );
    static void CheckError( const CCMIOError& err, const char *str );
    static void PrintData( int n, int id, DataType type, void *data,
                           void *data2 = NULL );

    int main( int argc, char *argv[])
    {
        int i = 0;
        CCMIOID root, state, processor, vertices, topology, solution, problem, next;
        CCMIOError err;

        // Open the file.  Because we did not initialize 'err' we need to pass
        // in NULL (which always means kCCMIONoErr) and then assign the return
        // value to 'err'.).
        err = CCMIOOpenFile(NULL, argv[1], kCCMIORead, &root);

        // We are going to assume that we have a state with a known name.
        // We could instead use CCMIONextEntity() to walk through all the states in
        // the file and present the list to the user for selection.
        CCMIOGetState(&err, root, kDefaultState, &problem, &state);

        unsigned int size;
        CCMIOEntityDescription(&err, state, &size, NULL);
        char *desc = new char[size + 1];
        CCMIOEntityDescription(&err, state, NULL, desc);
        cout << "Reading state '" << kDefaultState << "' (" << desc << ")" << endl;
        delete [] desc;

        // Find the first processor (i has previously been initialized to 0) and
        // read the mesh and solution information.
        CCMIONextEntity(&err, state, kCCMIOProcessor, &i, &processor);
        CCMIOReadProcessor(&err, processor, &vertices, &topology, NULL, &solution);

        if (err != kCCMIONoErr)
        {
            cout << "This example requires that the CCMIO file have both a mesh and post data" << endl;
            return(0);
        }

        // Print out the problem description (just cell types, and model constants)
        i = 0;
        if (CCMIOIsValidEntity(problem))   // if we have a problem description...
        {
            // ... walk through each cell type and print it...
            while (CCMIONextEntity(NULL, problem, kCCMIOCellType, &i, &next)
                                                                     == kCCMIONoErr)
            {
                char *name;
                int size, cellType;

                // ... if it has a material type.  (Note that we do not pass in
                // an array to get the name because we do not know how long the
                // string is yet.  Many parameters to CCMIO functions that return
                // data can be NULL if that data is not needed.)
                if (CCMIOReadOptstr(NULL, next, "MaterialType", &size, NULL)
                                                                       == kCCMIONoErr)
                {
                    name = new char[size + 1];
                    CCMIOReadOptstr(&err, next, "MaterialType", &size, name);
                    CCMIOGetEntityIndex(&err, next, &cellType);
                    PrintData(i, cellType, kCellType, name);
                    delete [] name;
                }
            }

            float value;
            CCMIOID constants;

            // Read some constants of particular interest, fall back to defaults
            // if they don't exist (or at least, that is what we would do if we
            // actually used them).  If you want to iterate over all the constants
            // instead, you will need to use CCMIOGetEntityNode() and use
            // CCMIOGetNextChild(), and CCMIOReadData() in the Core API.
            if (CCMIOGetEntity(NULL, problem, kCCMIOModelConstants, 0, &constants)
                                                                     == kCCMIONoErr)
            {
                if (CCMIOReadOptf(NULL, constants, "Gravity", &value) == kCCMIONoErr)
                    cout << "Gravity is " << value << " m/s^2" << endl;
                else
                    cout << "Gravity is unspecified" << endl;

                if (CCMIOReadOptf(NULL, constants, "Hubble constant", &value) == kCCMIONoErr)
                    cout << "Hubble constant is " << value << endl;
                else
                    cout << "Hubble constant is unspecified" << endl;
            }
        }

        ReadMesh(err, vertices, topology);
        ReadSolverInfo(err, solution);
        ReadPost(err, solution);

        // Read Lagrangian data
        vector<CCMIOID> ids;
        CCMIOID lagrangian;
        i = 0;
        while (CCMIONextEntity(&err, processor, kCCMIOLagrangianData, &i,
                               &lagrangian) == kCCMIONoErr)
        {
            CCMIOID solutionID, positionsID;
            CCMIOReadLagrangianData(&err, lagrangian, &positionsID, &solutionID);

            CCMIOEntityDescription(&err, lagrangian, &size, NULL);
            char *label = new char[size + 1];
            CCMIOEntityDescription(&err, lagrangian, NULL, label);
            cout << "Lagrangian data '" << label << "'" << endl;
            delete [] label;

            ReadVertices(err, positionsID, "positions");
            ReadPost(err, solutionID);
            ids.push_back(positionsID);
            ids.push_back(solutionID);
        }

        for (vector<CCMIOID>::iterator it = ids.begin();  it != ids.end();  ++it)
            CCMIOCloseFile(&err, *it);

        CCMIOCloseFile(&err, vertices);
        CCMIOCloseFile(&err, topology);
        CCMIOCloseFile(&err, solution);
        CCMIOCloseFile(&err, root);
    }

    void ReadMesh( CCMIOError& err, CCMIOID& vertices, CCMIOID& topology )
    {
        unsigned int i, j, nCells, nFaces, size;
        CCMIOID mapID, id;
        vector<int> mapData, faces, cells, faceCells;
        vector<float> verts;

        ReadVertices(err, vertices, "vertices", kVertOffset);

        // Read the cells.  Store cell IDs so that we know what cells are in
        // this post data.
        CCMIOGetEntity(&err, topology, kCCMIOCells, 0, &id);
        CCMIOEntitySize(&err, id, &nCells, NULL);
        cells.resize(nCells);
        mapData.resize(nCells);
        for (unsigned int k = 0;  k < nCells;  k += kCellInc)
        {
            CCMIOReadCells(&err, id, &mapID, &cells[k], k, k + kCellInc);
            CCMIOReadMap(&err, mapID, &mapData[k], k, k + kCellInc);
        }

        cout << "\t" << nCells << " cells" << endl;
        for (i = 0;  i < nCells && err == kCCMIONoErr;  ++i)
            PrintData(i, mapData[i], kCell, &cells[i]);
        CheckError(err, "Error reading cells");

        // Read the internal faces.
        CCMIOGetEntity(&err, topology, kCCMIOInternalFaces, 0, &id);
        CCMIOEntitySize(&err, id, &nFaces, NULL);
        mapData.resize(nFaces);
        faceCells.resize(2 * nFaces);
        CCMIOReadFaces(&err, id, kCCMIOInternalFaces, NULL, &size, NULL,
                       kCCMIOStart, kCCMIOEnd);
        faces.resize(size);
        CCMIOReadFaces(&err, id, kCCMIOInternalFaces, &mapID, NULL, &faces[0],
                       kCCMIOStart, kCCMIOEnd);
        CCMIOReadFaceCells(&err, id, kCCMIOInternalFaces, &faceCells[0],
                           kCCMIOStart, kCCMIOEnd);
        CCMIOReadMap(&err, mapID, &mapData[0], kCCMIOStart, kCCMIOEnd);

        cout << "\t" << nFaces << " faces" << endl;
        unsigned int pos = 0;
        i = 0;
        while (pos < faces.size())
        {
            PrintData(i, mapData[i], kInternalFace, &faces[pos], &faceCells[2 * i]);
            pos += faces[pos] + 1;
            i++;
        }
        CheckError(err, "Error reading internal faces");

        // Read the boundary faces.
        int index = 0;
        while (CCMIONextEntity(NULL, topology, kCCMIOBoundaryFaces, &index, &id)
                                                                     == kCCMIONoErr)
        {
            int boundaryVal;
            vector<int> prostarIDs;

            CCMIOEntitySize(&err, id, &nFaces, NULL);
            mapData.resize(nFaces);
            faceCells.resize(nFaces);
            CCMIOReadFaces(&err, id, kCCMIOBoundaryFaces, NULL, &size, NULL,
                           kCCMIOStart, kCCMIOEnd);
            faces.resize(size);
            CCMIOReadFaces(&err, id, kCCMIOBoundaryFaces, &mapID, NULL, &faces[0],
                           kCCMIOStart, kCCMIOEnd);
            CCMIOReadFaceCells(&err, id, kCCMIOBoundaryFaces, &faceCells[0],
                               kCCMIOStart, kCCMIOEnd);
            CCMIOReadMap(&err, mapID, &mapData[0], kCCMIOStart, kCCMIOEnd);

            // If the optional ProstarFaceId is there, read it
            if (CCMIOReadOpt1i(NULL, id, "ProstarFaceId", NULL,
                               kCCMIOStart, kCCMIOEnd) == kCCMIONoErr)
            {
                prostarIDs.resize(nFaces);
                CCMIOReadOpt1i(NULL, id, "ProstarFaceId", &prostarIDs[0],
                               kCCMIOStart, kCCMIOEnd);
            }

            CCMIOGetEntityIndex(&err, id, &boundaryVal);
            cout << "\t" << nFaces << " boundary faces (boundary " << boundaryVal
                 << ")" << endl;
            unsigned int pos = 0;
            j = 0;
            while (pos < faces.size())
            {
                PrintData(j, mapData[j], kBoundaryFace, &faces[pos],
                          &faceCells[j]);
                j++;
                pos += faces[pos] + 1;
            }
            CheckError(err, "Error reading boundary faces");
        }
    }

    void ReadVertices( CCMIOError& err, CCMIOID& vertices, const char *counter,
                       int offset )
    {
        int dims = 1;
        unsigned int i, nVertices, size;
        float scale;
        int mapData[kNValues] = { 0 };
        float verts[3 * kNValues] = { 0 };
        CCMIOID mapID;

        // Read the vertices.  This involves reading both the vertex data and
        // the map, which maps the index into the data array with the ID number.
        // As we process the vertices we need to be sure to scale them by the
        // appropriate scaling factor.  The offset is just to show you can read
        // any chunk.  Normally this would be in a for loop.
        CCMIOEntitySize(&err, vertices, &nVertices, NULL);
        CCMIOReadVerticesf(&err, vertices, &dims, &scale, &mapID, verts,
                           offset, offset + kNValues);
        CCMIOReadMap(&err, mapID, mapData, offset, offset + kNValues);

        CCMIOEntityDescription(&err, vertices, &size, NULL);
        char *label = new char[size + 1];
        CCMIOEntityDescription(&err, vertices, NULL, label);
        cout << "label: '" << label << "'" << endl;
        delete [] label;

        cout << "\t" << nVertices << " " << counter << endl;
        if (nVertices > (unsigned int)kNValues)
            nVertices = kNValues;
        for (i = 0;  i < nVertices && err == kCCMIONoErr;  ++i)
        {
            verts[dims * i    ] *= scale;
            verts[dims * i + 1] *= scale;
            verts[dims * i + 2] *= scale;   // This example assumes three
                                            // dimensional vertices.
            PrintData(i, mapData[i], kVertex, &verts[dims * i]);
        }
        cout << "\t\t..." << endl;
        CheckError(err, "Error reading vertices");
    }

    void ReadSolverInfo( CCMIOError& err, CCMIOID& solution)
    {
        char solver[kCCMIOMaxStringLength+1], timeUnits[kCCMIOMaxStringLength+1];
        int iterations;
        float time, angle;
        CCMIOID restart;

        if (CCMIOGetEntity(NULL, solution, kCCMIORestart, 0, &restart)!=kCCMIONoErr)
        {
            cout << "(No solver information)" << endl;
            return;
        }

        CCMIOReadRestartInfo(&err, restart, solver, &iterations, &time, timeUnits,
                             &angle);
        cout << "The following post data is the result of iteration "
             << iterations << " (time: " << time << " " << timeUnits
             << ")" << endl;
        cout << "of solver '" << solver << "':" << endl;

        // You might want to read more solver data here, in which case, do
        // CCMIOGetEntity(&err, restart, kCCMIORestartData, 0, &restartData);
        // and proceed to read the information.  Do likewise if you want to read
        // the reference data.
    }

    void ReadPost( CCMIOError& err, CCMIOID& solution )
    {
        bool oldFile = false;
        char name[kCCMIOMaxStringLength + 1];
        int h = 0, i = 0, dimSizes[] = { 0, 1, 3, 9 };
        CCMIOID field, phase;
        CCMIODimensionality dims;
        vector<int> mapData;
        vector<float> data;

        oldFile = (CCMIONextEntity(NULL, solution, kCCMIOFieldPhase, &h, &phase)
                                                                    != kCCMIONoErr);
        h = 0;
        while (oldFile ||
               CCMIONextEntity(NULL, solution, kCCMIOFieldPhase, &h, &phase)
                                                                    == kCCMIONoErr)
        {
            if (oldFile)
                phase = solution;
            else
            {
                int phaseNum = 0;
                CCMIOGetEntityIndex(NULL, phase, &phaseNum);
                cout << "Phase " << phaseNum << endl;
            }
            // Walk through each field in this field set (i.e. the post data)
            // and print out each one.
            while (CCMIONextEntity(NULL, phase, kCCMIOField, &i, &field)
                                                                    == kCCMIONoErr)
            {
                // Read the information about the field so that we know how to
                // process this field.
                char shortName[kCCMIOProstarShortNameLength+1], *units = NULL;
                int usize;
                CCMIODataType datatype;
                CCMIOReadField(&err, field, name, shortName, &dims, &datatype);
                if (CCMIOReadOptstr(NULL, field, kUnitsName, &usize, NULL)
                                                                     == kCCMIONoErr)
                {
                    units = new char[usize + 1];
                    CCMIOReadOptstr(&err, field, kUnitsName, NULL, units);
                }
                cout << "\tPost field '" << name << "'  (" << shortName << "):\t(";
                if (dims == 1)    // datatype is not meaningful for vectors/tensors
                {
                    if (datatype == kCCMIOFloat32)
                        cout << "Float32";
                    else if (datatype == kCCMIOFloat64)
                        cout << "Float64";
                    else if (datatype == kCCMIOInt32)
                        cout << "Int32";
                    else
                        cout << "Invalid datatype";
                }
                cout << ")\tUnits: " << ((units) ? units : "<none>") << endl;
                delete [] units;

                switch (dims)
                {
                    case kCCMIOScalar:
                        {
                        ReadScalar(err, field, mapData, data);
                        for (unsigned int k = 0;  k < data.size();  ++k)
                            PrintData(k, mapData[i], kScalar, &data[k]);
                        }
                        break;
                    case kCCMIOVector:
                        {
                        vector<float> u, v, w;
                        CCMIOID scalar;
                        CCMIOReadMultiDimensionalFieldData(&err, field,
                                                           kCCMIOVectorX, &scalar);
                        if (err == kCCMIOVersionErr)
                        {
                            // If we are reading an older version of the file,
                            // where vectors are stored as vectors, not components,
                            // we need to call CCMIOReadFieldData*(), which is
                            // all that ReadScalar() does.
                            err = kCCMIONoErr;
                            ReadScalar(err, field, mapData, data, true);
                            const unsigned int size = data.size() / dimSizes[dims];
                            for (unsigned int k = 0;  k < size;  ++k)
                                PrintData(k, mapData[i], kVector,
                                      &data[dimSizes[dims] * k]);
                        }
                        else
                        {
                            ReadScalar(err, scalar, mapData, u);
                            CCMIOReadMultiDimensionalFieldData(&err, field,
                                                               kCCMIOVectorY,
                                                               &scalar);
                            ReadScalar(err, scalar, mapData, v);
                            CCMIOReadMultiDimensionalFieldData(&err, field,
                                                               kCCMIOVectorZ,
                                                               &scalar);
                            ReadScalar(err, scalar, mapData, w);
                            data.resize(3 * u.size());
                            for (unsigned int k = 0;  k < u.size();  ++k)
                            {
                                data[3 * k    ] = u[k];
                                data[3 * k + 1] = v[k];
                                data[3 * k + 2] = w[k];
                            }
                            for (unsigned int k = 0;  k < u.size();  ++k)
                                PrintData(k, mapData[i], kVector,
                                          &data[dimSizes[dims] * k]);
                        }
                        }
                        break;
                    case kCCMIOTensor:
                        cout << "Tensor data not supported.  Ignoring field "
                             << name << "." << endl;
                        continue;
                }
                CheckError(err, "Error reading post data");
            }
            oldFile = false;
        }
    }

    void ReadScalar( CCMIOError& err, CCMIOID field,
                     vector<int>& mapData, vector<float>& data,
                     bool readingVector /* = false */)
    {
        unsigned int n, max;
        int j = 0;
        CCMIOID fieldData, mapID;
        CCMIODataLocation type;

        // Read each piece of field data
        while (CCMIONextEntity(NULL, field, kCCMIOFieldData, &j, &fieldData)
                                                                     == kCCMIONoErr)
        {
            // Figure out how big this data is so we can read it. If we were
            // storing this information permanently we might use a sparse
            // array, in which case we would need to find the maximum ID and
            // make the array that size.
            CCMIOEntitySize(&err, fieldData, &n, &max);
            mapData.resize(n);
            CCMIOReadFieldDataf(&err, fieldData, &mapID, &type, NULL,
                                kCCMIOStart, kCCMIOEnd);
            CCMIOReadMap(&err, mapID, &mapData[0], kCCMIOStart, kCCMIOEnd);

            // We are only going to process cell data.  Vertex data would
            // be processed similarly. If your appliation has only one value
            // for boundary data, you would separate the face data into
            // each boundary and combine it together.  If your application
            // stores boundary data on each face you could read it in using
            // similar procedures for the cell and vertex data.  Note that
            // the file may not contain all types of data.
            if (type == kCCMIOCell)
            {
                if (readingVector)
                    data.resize(3 * n);
                else
                    data.resize(n);
                // If we want double precision we should use
                // CCMIOReadFieldDatad().
                CCMIOReadFieldDataf(&err, fieldData, &mapID, NULL,
                                    &data[0], kCCMIOStart, kCCMIOEnd);
            }

            CheckError(err, "Error reading post data");
        }
    }
    //--------------------------- Helper functions --------------------------------
    void CheckError( const CCMIOError& err, const char *str )
    {
        if (err == kCCMIONoErr)
            return;

        cout << str << " (error " << err << ")" << endl;
        exit(1);
    }

    void PrintData( int n, int id, DataType type, void *data, void *data2 = NULL )
    {
        if (n > kNValues)
            return;
        if (n == kNValues)
        {
            cout << "\t\t..." << endl;
            return;
        }

        switch (type)
        {
            case kScalar:
                cout << "\t\t" << *(float *)data << endl;
                break;
            case kVector:
                cout << "\t\t(" << ((float *)data)[0] << ", "
                     << ((float *)data)[1] << ", " << ((float *)data)[2] << ")"
                     << endl;
                break;
            case kVertex:
                cout << "\t\t" << "Vertex " << id << ":  ("
                     << ((float *)data)[0] << ", "
                     << ((float *)data)[1] << ", "
                     << ((float *)data)[2] << ")" << endl;
                break;
            case kCell:
                cout << "\t\t" << "Cell " << id << ":  type " << *(int *)data
                     << endl;
                break;
            case kInternalFace:
                {
                    int nVerts = *(int *)data;
                    cout << "\t\t" << "Face " << id << ": \tCells: "
                         << ((int *)data2)[0] << ", " << ((int *)data2)[1]
                         << "\tVertices: ";
                    for (int i = 0;  i < nVerts;  ++i)
                        cout << *(((int *)data) + 1 + i) << " ";
                    cout << endl;
                }
                break;
            case kBoundaryFace:
                {
                    int nVerts = *(int *)data;
                    cout << "\t\t" << "Face " << id << ": \tCell: "
                         << *(int *)data2 << "\tVertices: ";
                    for (int i = 0;  i < nVerts;  ++i)
                        cout << *(((int *)data) + 1 + i) << " ";
                    cout << endl;
                }
                break;
            case kBoundaryData:
                cout << "\t\tBoundary " << *(int *)data << "  priority "
                     << *(int *)data2 << endl;
                break;
            case kBoundaryFaceData:
                {
                    cout << "\t\tBoundary data [ ";
                    for (int i = 0;  i < *(int *)data2; ++i)
                        cout << ((float *)data)[i] << " ";
                    cout << "]" << endl;
                }
                break;
            case kCellType:
                cout << "Cell type " << id << ":  " << (char *)data << endl;
                break;
        }
    }