MOAB: Mesh Oriented datABase  (version 5.2.1)
addncdata.cpp File Reference
#include "moab/ProgOptions.hpp"
#include "moab/Core.hpp"
#include "netcdf.h"
#include <math.h>
#include <sstream>
+ Include dependency graph for addncdata.cpp:

Go to the source code of this file.

Defines

#define INS_ID(stringvar, prefix, id)   sprintf( stringvar, prefix, id )
#define GET_DIM(ncdim, name, val)
#define GET_DIMB(ncdim, name, varname, id, val)
#define GET_VAR(name, id, dims)
#define GET_1D_INT_VAR(name, id, vals)
#define GET_1D_DBL_VAR(name, id, vals)
#define GET_1D_FLT_VAR(name, id, vals)

Functions

int main (int argc, char *argv[])

Variables

int ncFile

Define Documentation

#define GET_1D_DBL_VAR (   name,
  id,
  vals 
)
Value:
{                                                                                                               \
        std::vector< int > dum_dims;                                                                                \
        GET_VAR( name, id, dum_dims );                                                                              \
        if( -1 != id )                                                                                              \
        {                                                                                                           \
            size_t ntmp;                                                                                            \
            int dvfail = nc_inq_dimlen( ncFile, dum_dims[0], &ntmp );                                               \
            if( NC_NOERR != dvfail ) { MB_SET_ERR( MB_FAILURE, "addncdata:: Couldn't get dimension length" ); }     \
            vals.resize( ntmp );                                                                                    \
            size_t ntmp1 = 0;                                                                                       \
            dvfail       = nc_get_vara_double( ncFile, id, &ntmp1, &ntmp, &vals[0] );                               \
            if( NC_NOERR != dvfail ) { MB_SET_ERR( MB_FAILURE, "addncdata:: Problem getting variable " << name ); } \
        }                                                                                                           \
    }

Definition at line 89 of file addncdata.cpp.

Referenced by main().

#define GET_1D_FLT_VAR (   name,
  id,
  vals 
)
Value:
{                                                                                                               \
        std::vector< int > dum_dims;                                                                                \
        GET_VAR( name, id, dum_dims );                                                                              \
        if( -1 != id )                                                                                              \
        {                                                                                                           \
            size_t ntmp;                                                                                            \
            int dvfail = nc_inq_dimlen( ncFile, dum_dims[0], &ntmp );                                               \
            if( NC_NOERR != dvfail ) { MB_SET_ERR( MB_FAILURE, "addncdata:: Couldn't get dimension length" ); }     \
            vals.resize( ntmp );                                                                                    \
            size_t ntmp1 = 0;                                                                                       \
            dvfail       = nc_get_vara_float( ncFile, id, &ntmp1, &ntmp, &vals[0] );                                \
            if( NC_NOERR != dvfail ) { MB_SET_ERR( MB_FAILURE, "addncdata:: Problem getting variable " << name ); } \
        }                                                                                                           \
    }

Definition at line 105 of file addncdata.cpp.

Referenced by main().

#define GET_1D_INT_VAR (   name,
  id,
  vals 
)
Value:
{                                                                                                               \
        GET_VAR( name, id, vals );                                                                                  \
        if( -1 != id )                                                                                              \
        {                                                                                                           \
            size_t ntmp;                                                                                            \
            int ivfail = nc_inq_dimlen( ncFile, vals[0], &ntmp );                                                   \
            if( NC_NOERR != ivfail ) { MB_SET_ERR( MB_FAILURE, "addncdata:: Couldn't get dimension length" ); }     \
            vals.resize( ntmp );                                                                                    \
            size_t ntmp1 = 0;                                                                                       \
            ivfail       = nc_get_vara_int( ncFile, id, &ntmp1, &ntmp, &vals[0] );                                  \
            if( NC_NOERR != ivfail ) { MB_SET_ERR( MB_FAILURE, "addncdata:: Problem getting variable " << name ); } \
        }                                                                                                           \
    }

Definition at line 74 of file addncdata.cpp.

Referenced by main().

#define GET_DIM (   ncdim,
  name,
  val 
)
Value:
{                                                                                                           \
        int gdfail = nc_inq_dimid( ncFile, name, &ncdim );                                                      \
        if( NC_NOERR == gdfail )                                                                                \
        {                                                                                                       \
            size_t tmp_val;                                                                                     \
            gdfail = nc_inq_dimlen( ncFile, ncdim, &tmp_val );                                                  \
            if( NC_NOERR != gdfail ) { MB_SET_ERR( MB_FAILURE, "addncdata:: Couldn't get dimension length" ); } \
            else                                                                                                \
                val = tmp_val;                                                                                  \
        }                                                                                                       \
        else                                                                                                    \
            val = 0;                                                                                            \
    }

Definition at line 37 of file addncdata.cpp.

#define GET_DIMB (   ncdim,
  name,
  varname,
  id,
  val 
)
Value:
INS_ID( name, varname, id );                  \
    GET_DIM( ncdim, name, val );

Definition at line 52 of file addncdata.cpp.

#define GET_VAR (   name,
  id,
  dims 
)
Value:
{                                                                                            \
        id         = -1;                                                                         \
        int gvfail = nc_inq_varid( ncFile, name, &id );                                          \
        if( NC_NOERR == gvfail )                                                                 \
        {                                                                                        \
            int ndims;                                                                           \
            gvfail = nc_inq_varndims( ncFile, id, &ndims );                                      \
            if( NC_NOERR == gvfail )                                                             \
            {                                                                                    \
                dims.resize( ndims );                                                            \
                gvfail = nc_inq_vardimid( ncFile, id, &dims[0] );                                \
                if( NC_NOERR != gvfail )                                                         \
                { MB_SET_ERR( MB_FAILURE, "addncdata:: Couldn't get variable dimension IDs" ); } \
            }                                                                                    \
        }                                                                                        \
    }

Definition at line 56 of file addncdata.cpp.

Referenced by main().

#define INS_ID (   stringvar,
  prefix,
  id 
)    sprintf( stringvar, prefix, id )

Definition at line 35 of file addncdata.cpp.


Function Documentation

int main ( int  argc,
char *  argv[] 
)

Definition at line 121 of file addncdata.cpp.

References ProgOptions::addOpt(), moab::Range::begin(), edges, moab::Range::end(), moab::fail(), GET_1D_DBL_VAR, GET_1D_FLT_VAR, GET_1D_INT_VAR, moab::Core::get_entities_by_dimension(), GET_VAR, moab::Range::index(), moab::Core::load_file(), mb, MB_CHK_SET_ERR, MB_FAILURE, MB_FILE_DOES_NOT_EXIST, MB_SET_ERR, MB_TAG_CREAT, MB_TAG_DENSE, MB_TYPE_DOUBLE, MB_TYPE_INTEGER, ncFile, np, outfile, ProgOptions::parseCommandLine(), moab::Range::size(), moab::Core::tag_get_data(), moab::Core::tag_get_handle(), moab::Core::tag_get_length(), moab::Core::tag_set_data(), and moab::Core::write_file().

{

    ProgOptions opts;

    std::string inputfile, outfile( "out.h5m" ), netcdfFile, variable_name, sefile_name;

    opts.addOpt< std::string >( "input,i", "input mesh filename", &inputfile );
    opts.addOpt< std::string >( "netcdfFile,n", "netcdf file aligned with the mesh input file", &netcdfFile );
    opts.addOpt< std::string >( "output,o", "output mesh filename", &outfile );

    opts.addOpt< std::string >( "var,v", "variable to extract and add to output file", &variable_name );

    opts.addOpt< std::string >( "sefile,s", "spectral elements file (coarse SE mesh)", &sefile_name );
    opts.parseCommandLine( argc, argv );

    ErrorCode rval;
    Core* mb = new Core();

    rval = mb->load_file( inputfile.c_str() );MB_CHK_SET_ERR( rval, "can't load input file" );

    std::cout << " opened " << inputfile << " with initial h5m data.\n";
    // open the netcdf file, and see if it has that variable we are looking for

    Range nodes;
    rval = mb->get_entities_by_dimension( 0, 0, nodes );MB_CHK_SET_ERR( rval, "can't get nodes" );

    Range edges;
    rval = mb->get_entities_by_dimension( 0, 1, edges );MB_CHK_SET_ERR( rval, "can't get edges" );

    Range cells;
    rval = mb->get_entities_by_dimension( 0, 2, cells );MB_CHK_SET_ERR( rval, "can't get cells" );

    std::cout << " it has " << nodes.size() << " vertices " << edges.size() << " edges " << cells.size() << " cells\n";

    // construct maps between global id and handles
    std::map< int, EntityHandle > vGidHandle;
    std::map< int, EntityHandle > eGidHandle;
    std::map< int, EntityHandle > cGidHandle;
    std::vector< int > gids;
    Tag gid;
    rval = mb->tag_get_handle( "GLOBAL_ID", gid );MB_CHK_SET_ERR( rval, "can't get global id tag" );
    gids.resize( nodes.size() );
    rval = mb->tag_get_data( gid, nodes, &gids[0] );MB_CHK_SET_ERR( rval, "can't get global id on vertices" );
    int i = 0;
    for( Range::iterator vit = nodes.begin(); vit != nodes.end(); vit++ )
    {
        vGidHandle[gids[i++]] = *vit;
    }

    gids.resize( edges.size() );
    rval = mb->tag_get_data( gid, edges, &gids[0] );MB_CHK_SET_ERR( rval, "can't get global id on edges" );
    i = 0;
    for( Range::iterator vit = edges.begin(); vit != edges.end(); vit++ )
    {
        eGidHandle[gids[i++]] = *vit;
    }

    gids.resize( cells.size() );
    rval = mb->tag_get_data( gid, cells, &gids[0] );MB_CHK_SET_ERR( rval, "can't get global id on cells" );
    i = 0;
    for( Range::iterator vit = cells.begin(); vit != cells.end(); vit++ )
    {
        cGidHandle[gids[i++]] = *vit;
    }

    // Open netcdf/exodus file
    int fail = nc_open( netcdfFile.c_str(), 0, &ncFile );
    if( NC_NOWRITE != fail )
    { MB_SET_ERR( MB_FILE_DOES_NOT_EXIST, "ReadNCDF:: problem opening Netcdf II file " << netcdfFile ); }

    std::cout << " opened " << netcdfFile << " with new data \n";
    std::vector< int > dims;
    int nc_var;

    size_t recs;
    char recname[NC_MAX_NAME + 1];

    std::cout << " looking for variable " << variable_name << "\n";
    GET_VAR( variable_name.c_str(), nc_var, dims );
    std::cout << " it has " << dims.size() << " dimensions\n";

    int dimIndex     = -1;  // index of the dimension of interest
    bool vertex_data = false;
    bool cell_data   = false;
    for( size_t j = 0; j < dims.size(); j++ )
    {
        fail = nc_inq_dim( ncFile, dims[j], recname, &recs );
        std::string name_dim( recname );
        std::cout << " dimension index " << j << " in file: " << dims[j] << " name: " << name_dim << " recs:" << recs
                  << "\n";
        if( recs == nodes.size() )
        {
            dimIndex    = j;
            vertex_data = true;
        }
        if( recs == cells.size() )
        {
            dimIndex  = j;
            cell_data = true;
        }
    }
    int otherDim    = 1 - dimIndex;  // used only if 2 dimensions ; could be 0 or 1;
    size_t size_tag = 1;             // good for one dimension
    std::vector< int > evals;        // size of size_tag
    std::vector< double > dvals;     // size of size_tag
    nc_type dataType;
    // read the variable, and set it to the tag
    fail            = nc_inq_vartype( ncFile, nc_var, &dataType );
    DataType mbtype = MB_TYPE_DOUBLE;
    bool float_var  = false;
    if( NC_INT == dataType )
        mbtype = MB_TYPE_INTEGER;
    else if( NC_DOUBLE != dataType && NC_FLOAT != dataType )
        MB_CHK_SET_ERR( MB_FAILURE, "unknown type" );

    if( NC_FLOAT == dataType ) float_var = true;

    int time_id = -1;
    fail        = nc_inq_varid( ncFile, "time", &time_id );
    std::vector< float > times;
    if( NC_NOERR == fail )
    {
        int ii;
        GET_1D_FLT_VAR( "time", ii, times );
    }
    Tag newTag;

    if( ( dims.size() >= 1 && dims.size() <= 2 ) && ( vertex_data || cell_data ) )
    {

        if( dims.size() == 2 ) { fail = nc_inq_dim( ncFile, dims[otherDim], recname, &size_tag ); }

        int def_val = 0;
        rval = mb->tag_get_handle( variable_name.c_str(), (int)size_tag, mbtype, newTag, MB_TAG_CREAT | MB_TAG_DENSE,
                                   &def_val );MB_CHK_SET_ERR( rval, "can't define new tag" );

        if( NC_INT == dataType )
        {
            std::vector< int > vals;
            if( 1 == dims.size() )
            {
                GET_1D_INT_VAR( variable_name.c_str(), dims[0], vals );
                if( vertex_data )
                {
                    for( size_t k = 0; k < vals.size(); k++ )
                    {
                        EntityHandle vh = vGidHandle[k + 1];
                        rval            = mb->tag_set_data( newTag, &vh, 1, &vals[k] );MB_CHK_SET_ERR( rval, "can't set tag on vertex" );
                    }
                }
                else  // cell_data
                {
                    for( size_t k = 0; k < vals.size(); k++ )
                    {
                        EntityHandle ch = cGidHandle[k + 1];  // cell handle
                        rval            = mb->tag_set_data( newTag, &ch, 1, &vals[k] );MB_CHK_SET_ERR( rval, "can't set tag on cell" );
                    }
                }
            }
            else  // dims.size() == 2
            {
                // Single var for all coords
                size_t start[2] = { 0, 0 }, count[2] = { 1, 1 };
                count[dimIndex] = recs;
                count[otherDim] = size_tag;
                vals.resize( recs * size_tag );
                fail = nc_get_vara_int( ncFile, nc_var, start, count, &vals[0] );
                evals.resize( size_tag );

                if( vertex_data )
                {
                    for( size_t k = 0; k < recs; k++ )
                    {
                        EntityHandle vh      = vGidHandle[k + 1];
                        size_t start_in_vals = k * size_tag, stride = 1;
                        for( size_t j = 0; j < size_tag; j++ )
                            evals[j] = vals[start_in_vals + j * stride];
                        rval = mb->tag_set_data( newTag, &vh, 1, &evals[0] );MB_CHK_SET_ERR( rval, "can't set tag on vertex" );
                    }
                }
                else  // cell_data
                {
                    for( size_t k = 0; k < recs; k++ )
                    {
                        EntityHandle ch      = cGidHandle[k + 1];  // cell handle
                        size_t start_in_vals = k * size_tag, stride = 1;
                        for( size_t j = 0; j < size_tag; j++ )
                            evals[j] = vals[start_in_vals + j * stride];
                        rval = mb->tag_set_data( newTag, &ch, 1, &evals[0] );MB_CHK_SET_ERR( rval, "can't set tag on cell" );
                    }
                }
            }
        }
        else
        {
            std::vector< double > vals;
            if( 1 == dims.size() )
            {
                GET_1D_DBL_VAR( variable_name.c_str(), dims[0], vals );
                if( vertex_data )
                {
                    for( size_t k = 0; k < vals.size(); k++ )
                    {
                        EntityHandle vh = vGidHandle[k + 1];
                        rval            = mb->tag_set_data( newTag, &vh, 1, &vals[k] );MB_CHK_SET_ERR( rval, "can't set tag on vertex" );
                    }
                }
                else  // cell_data
                {
                    for( size_t k = 0; k < vals.size(); k++ )
                    {
                        EntityHandle ch = cGidHandle[k + 1];  // cell handle
                        rval            = mb->tag_set_data( newTag, &ch, 1, &vals[k] );MB_CHK_SET_ERR( rval, "can't set tag on vertex" );
                    }
                }
            }
            else  // dims.size() == 2
            {
                // Single var for all coords
                size_t start[2] = { 0, 0 }, count[2] = { 1, 1 };
                count[dimIndex] = recs;
                count[otherDim] = size_tag;
                vals.resize( recs * size_tag );
                fail = nc_get_vara_double( ncFile, nc_var, start, count, &vals[0] );
                dvals.resize( size_tag );

                if( vertex_data )
                {
                    for( size_t k = 0; k < recs; k++ )
                    {
                        EntityHandle vh      = vGidHandle[k + 1];
                        size_t start_in_vals = k * size_tag, stride = 1;
                        for( size_t j = 0; j < size_tag; j++ )
                            dvals[j] = vals[start_in_vals + j * stride];
                        rval = mb->tag_set_data( newTag, &vh, 1, &dvals[0] );MB_CHK_SET_ERR( rval, "can't set tag on vertex" );
                    }
                }
                else  // cell_data
                {
                    for( size_t k = 0; k < recs; k++ )
                    {
                        EntityHandle ch      = cGidHandle[k + 1];  // cell handle
                        size_t start_in_vals = k * size_tag, stride = 1;
                        for( size_t j = 0; j < size_tag; j++ )
                            dvals[j] = vals[start_in_vals + j * stride];
                        rval = mb->tag_set_data( newTag, &ch, 1, &dvals[0] );MB_CHK_SET_ERR( rval, "can't set tag on cell" );
                    }
                }
            }
        }
    }
    else if( ( dims.size() == 3 ) && vertex_data && dimIndex == 2 &&
             mbtype == MB_TYPE_DOUBLE )  // the last one is the vertex
    {
        // the case when the last dim is ncol (for homme type mesh))
        size_t dim0, dim1;  // dim 2 is ncol..
        char recname0[NC_MAX_NAME + 1], recname1[NC_MAX_NAME + 1];
        fail = nc_inq_dim( ncFile, dims[0], recname0, &dim0 );
        fail = nc_inq_dim( ncFile, dims[1], recname1, &dim1 );
        std::string name0( recname0 );
        std::string name1( recname1 );
        std::string timestr( "time" );
        size_t start[3] = { 0, 0, 0 };
        size_t count[3] = { 1, 0, 0 };
        count[1]        = dim1;
        count[2]        = nodes.size();
        // create a few tags, with name inserted
        if( name0.compare( "time" ) == 0 )
        {

            std::vector< double > dvalues( dim1 * count[2] );
            std::vector< float > fvalues( dim1 * count[2] );
            std::vector< double > transp( dim1 * count[2] );
            // get the variable time
            for( size_t k = 0; k < dim0; k++ )
            {
                // create a tag for each time, and
                std::stringstream tag_name;
                tag_name << variable_name << "_t" << times[k];
                std::vector< double > defvals( dim1, 0. );
                rval = mb->tag_get_handle( tag_name.str().c_str(), (int)dim1, mbtype, newTag,
                                           MB_TAG_CREAT | MB_TAG_DENSE, &defvals[0] );MB_CHK_SET_ERR( rval, "can't define new tag" );
                start[0] = k;

                if( float_var )
                {
                    fail = nc_get_vara_float( ncFile, nc_var, start, count, &fvalues[0] );
                    // now arrange them in a tag, transpose data
                    for( size_t ii = 0; ii < dim1; ii++ )
                        for( size_t j = 0; j < count[2]; j++ )
                        {
                            transp[j * dim1 + ii] = fvalues[ii * count[2] + j];
                        }
                }
                else  // double
                {
                    fail = nc_get_vara_double( ncFile, nc_var, start, count, &dvalues[0] );
                    // now arrange them in a tag, transpose data
                    for( size_t ii = 0; ii < dim1; ii++ )
                        for( size_t j = 0; j < count[2]; j++ )
                        {
                            transp[j * dim1 + ii] = dvalues[ii * count[2] + j];
                        }
                }
                for( size_t ii = 0; ii < nodes.size(); ii++ )
                {
                    EntityHandle vh = vGidHandle[ii + 1];
                    rval            = mb->tag_set_data( newTag, &vh, 1, &transp[ii * dim1] );MB_CHK_SET_ERR( rval, "can't set tag on nodes" );
                }
            }
        }
    }
    rval = mb->write_file( outfile.c_str() );MB_CHK_SET_ERR( rval, "can't write file" );
    std::cout << " wrote file " << outfile << "\n";

    // now, if s option, load the coarse mesh and put data on each element, according to a matrix
    if( !sefile_name.empty() )
    {
        // load the file, check for GLOBAL_DOFS tag, and create a new tag with the data associated
        Core* mb2 = new Core();
        rval      = mb2->load_file( sefile_name.c_str() );MB_CHK_SET_ERR( rval, "can't load spectral element file" );
        std::cout << " loaded spectral file " << sefile_name << "\n";
        // look for GLOBAL_DOFS tag
        Tag gdofeTag;
        rval = mb2->tag_get_handle( "GLOBAL_DOFS", gdofeTag );MB_CHK_SET_ERR( rval, "file does not have GLOBAL_DOFS file" );
        int sizeTag;
        rval = mb2->tag_get_length( gdofeTag, sizeTag );MB_CHK_SET_ERR( rval, "can't get size of tag" );
        int np = (int)sqrt( 1.0 * sizeTag );
        std::cout << " size of tag: " << sizeTag << " np = " << np << "\n";
        std::vector< int > gdofs;
        Range cells2;
        rval = mb2->get_entities_by_dimension( 0, 2, cells2 );MB_CHK_SET_ERR( rval, "can't get cells on spectral mesh" );
        gdofs.resize( cells2.size() * sizeTag );
        rval = mb2->tag_get_data( gdofeTag, cells2, &gdofs[0] );MB_CHK_SET_ERR( rval, "can't get global dofs tag" );
        // create a new tag for element data arranged as DFIELD

        std::vector< double > dfield;
        dfield.resize( sizeTag, 0.0 );
        Tag newTag2;
        rval = mb2->tag_get_handle( variable_name.c_str(), (int)sizeTag, mbtype, newTag2, MB_TAG_CREAT | MB_TAG_DENSE,
                                    &dfield[0] );MB_CHK_SET_ERR( rval, "can't define new tag" );

        int i1 = 0;  // index in the gdofs array, per element

        // get the tag values from the other moab core, for newTag
        int dataTagLen;
        rval = mb->tag_get_length( newTag, dataTagLen );MB_CHK_SET_ERR( rval, "can't get size of newTag" );
        //
        std::vector< double > oldData;
        oldData.resize( dataTagLen * nodes.size() );  //

        // get the "old" values
        rval = mb->tag_get_data( newTag, nodes, &oldData[0] );MB_CHK_SET_ERR( rval, "can't get old values" );
        for( Range::iterator it = cells2.begin(); it != cells2.end(); it++ )
        {
            EntityHandle cel = *it;
            // gdofs per element are gdofs[i:i+np*np];
            for( int k = 0; k < sizeTag; k++ )
            {
                int gdof          = gdofs[i1 + k];
                EntityHandle node = vGidHandle[gdof];
                int index         = nodes.index( node );
                // get last layer of data
                double val = oldData[index * dataTagLen + dataTagLen - 1];  //
                dfield[k]  = val;
            }
            i1   = i1 + sizeTag;
            rval = mb2->tag_set_data( newTag2, &cel, 1, &dfield[0] );MB_CHK_SET_ERR( rval, "can't set new tag" );
        }

        // write the appended file with the new field:
        rval = mb2->write_file( "atm2.h5m" );MB_CHK_SET_ERR( rval, "can't write new spectral file" );
        std::cout << " wrote file atm2.h5m \n";
    }
    return 0;
}

Variable Documentation

int ncFile

Definition at line 33 of file addncdata.cpp.

Referenced by main().

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines