moab
moab::NCHelperHOMME Class Reference

Child helper class for HOMME grid (CAM_SE) More...

#include <NCHelperHOMME.hpp>

Inheritance diagram for moab::NCHelperHOMME:
moab::UcdNCHelper moab::NCHelper

List of all members.

Public Member Functions

 NCHelperHOMME (ReadNC *readNC, int fileId, const FileOptions &opts, EntityHandle fileSet)

Static Public Member Functions

static bool can_read_file (ReadNC *readNC, int fileId)

Private Member Functions

virtual ErrorCode init_mesh_vals ()
 Implementation of NCHelper::init_mesh_vals()
virtual ErrorCode check_existing_mesh ()
 Implementation of NCHelper::check_existing_mesh()
virtual ErrorCode create_mesh (Range &faces)
 Implementation of NCHelper::create_mesh()
virtual std::string get_mesh_type_name ()
 Implementation of NCHelper::get_mesh_type_name()
virtual ErrorCode read_ucd_variable_to_nonset_allocate (std::vector< ReadNC::VarData > &vdatas, std::vector< int > &tstep_nums)
 Implementation of UcdNCHelper::read_ucd_variable_to_nonset_allocate()
virtual ErrorCode read_ucd_variable_to_nonset (std::vector< ReadNC::VarData > &vdatas, std::vector< int > &tstep_nums)
 Implementation of UcdNCHelper::read_ucd_variable_to_nonset()

Private Attributes

int _spectralOrder
int connectId
bool isConnFile

Detailed Description

Child helper class for HOMME grid (CAM_SE)

Definition at line 17 of file NCHelperHOMME.hpp.


Constructor & Destructor Documentation

moab::NCHelperHOMME::NCHelperHOMME ( ReadNC readNC,
int  fileId,
const FileOptions opts,
EntityHandle  fileSet 
)

Definition at line 16 of file NCHelperHOMME.cpp.

: UcdNCHelper(readNC, fileId, opts, fileSet),
_spectralOrder(-1), connectId(-1), isConnFile(false)
{
  // Calculate spectral order
  std::map<std::string, ReadNC::AttData>::iterator attIt = readNC->globalAtts.find("np");
  if (attIt != readNC->globalAtts.end()) {
    int success = NCFUNC(get_att_int)(readNC->fileId, attIt->second.attVarId, attIt->second.attName.c_str(), &_spectralOrder);
    if (success != 0)
      readNC->readMeshIface->report_error("%s", "Failed to read np global attribute int data.");
    else
      _spectralOrder--; // Spectral order is one less than np
  }
  else {
    // As can_read_file() returns true and there is no global attribute "np", it should be a connectivity file
    isConnFile = true;
    _spectralOrder = 3; // Assume np is 4
  }
}

Member Function Documentation

bool moab::NCHelperHOMME::can_read_file ( ReadNC readNC,
int  fileId 
) [static]

Definition at line 36 of file NCHelperHOMME.cpp.

{
  // If global attribute "np" exists then it should be the HOMME grid
  if (readNC->globalAtts.find("np") != readNC->globalAtts.end()) {
    // Make sure it is CAM grid
    std::map<std::string, ReadNC::AttData>::iterator attIt = readNC->globalAtts.find("source");
    if (attIt == readNC->globalAtts.end()) {
      readNC->readMeshIface->report_error("%s", "File does not have source global attribute.");
      return false;
    }
    unsigned int sz = attIt->second.attLen;
    std::string att_data;
    att_data.resize(sz + 1);
    att_data[sz] = '\000';
    int success = NCFUNC(get_att_text)(fileId, attIt->second.attVarId, attIt->second.attName.c_str(), &att_data[0]);
    if (success != 0) {
      readNC->readMeshIface->report_error("%s", "Failed to read source global attribute char data.");
      return false;
    }
    if (att_data.find("CAM") == std::string::npos)
      return false;

    return true;
  }
  else {
    // If dimension names "ncol" AND "ncorners" AND "ncells" exist, then it should be the HOMME connectivity file
    // In this case, the mesh can still be created
    std::vector<std::string>& dimNames = readNC->dimNames;
    if ((std::find(dimNames.begin(), dimNames.end(), std::string("ncol")) != dimNames.end()) &&
        (std::find(dimNames.begin(), dimNames.end(), std::string("ncorners")) != dimNames.end()) &&
        (std::find(dimNames.begin(), dimNames.end(), std::string("ncells")) != dimNames.end()))
      return true;
  }

  return false;
}

Implementation of NCHelper::check_existing_mesh()

Implements moab::NCHelper.

Definition at line 215 of file NCHelperHOMME.cpp.

{
  Interface*& mbImpl = _readNC->mbImpl;
  Tag& mGlobalIdTag = _readNC->mGlobalIdTag;
  bool& noMesh = _readNC->noMesh;

  if (noMesh && localGidVerts.empty()) {
    // We need to populate localGidVerts range with the gids of vertices from the tmp_set
    // localGidVerts is important in reading the variable data into the nodes
    // also, for our purposes, localGidVerts is truly the GLOBAL_ID tag data, not other
    // file_id tags that could get passed around in other scenarios for parallel reading
    // for nodal_partition, this local gid is easier, should be initialized with only
    // the owned nodes

    // We need to get all vertices from tmp_set (it is the input set in no_mesh scenario)
    Range local_verts;
    ErrorCode rval = mbImpl->get_entities_by_dimension(_fileSet, 0, local_verts);
    if (MB_FAILURE == rval)
      return rval;

    if (!local_verts.empty()) {
      std::vector<int> gids(local_verts.size());

      // !IMPORTANT : this has to be the GLOBAL_ID tag
      rval = mbImpl->tag_get_data(mGlobalIdTag, local_verts, &gids[0]);
      if (MB_FAILURE == rval)
        return rval;

      // Restore localGidVerts
      std::copy(gids.rbegin(), gids.rend(), range_inserter(localGidVerts));
      nLocalVertices = localGidVerts.size();
    }
  }

  return MB_SUCCESS;
}
ErrorCode moab::NCHelperHOMME::create_mesh ( Range faces) [private, virtual]

Implementation of NCHelper::create_mesh()

Implements moab::NCHelper.

Definition at line 252 of file NCHelperHOMME.cpp.

{
  Interface*& mbImpl = _readNC->mbImpl;
  std::string& fileName = _readNC->fileName;
  Tag& mGlobalIdTag = _readNC->mGlobalIdTag;
  const Tag*& mpFileIdTag = _readNC->mpFileIdTag;
  DebugOutput& dbgOut = _readNC->dbgOut;
  bool& spectralMesh = _readNC->spectralMesh;
  int& gatherSetRank = _readNC->gatherSetRank;

  int rank = 0;
  int procs = 1;
#ifdef USE_MPI
  bool& isParallel = _readNC->isParallel;
  if (isParallel) {
    ParallelComm*& myPcomm = _readNC->myPcomm;
    rank = myPcomm->proc_config().proc_rank();
    procs = myPcomm->proc_config().proc_size();
  }
#endif

  ErrorCode rval;
  int success = 0;

  // Need to get/read connectivity data before creating elements
  std::string conn_fname;

  if (isConnFile) {
    // Connectivity file has already been read
    connectId = _readNC->fileId;
  }
  else {
    // Try to open the connectivity file through CONN option, if used
    rval = _opts.get_str_option("CONN", conn_fname);
    if (MB_SUCCESS != rval) {
      // Default convention for reading HOMME is a file HommeMapping.nc in same dir as data file
      conn_fname = std::string(fileName);
      size_t idx = conn_fname.find_last_of("/");
      if (idx != std::string::npos)
        conn_fname = conn_fname.substr(0, idx).append("/HommeMapping.nc");
      else
        conn_fname = "HommeMapping.nc";
    }
#ifdef PNETCDF_FILE
#ifdef USE_MPI
    if (isParallel) {
      ParallelComm*& myPcomm = _readNC->myPcomm;
      success = NCFUNC(open)(myPcomm->proc_config().proc_comm(), conn_fname.c_str(), 0, MPI_INFO_NULL, &connectId);
    }
    else
      success = NCFUNC(open)(MPI_COMM_SELF, conn_fname.c_str(), 0, MPI_INFO_NULL, &connectId);
#endif
#else
    success = NCFUNC(open)(conn_fname.c_str(), 0, &connectId);
#endif
    ERRORS(success, "Failed on open.");
  }

  std::vector<std::string> conn_names;
  std::vector<int> conn_vals;
  rval = _readNC->get_dimensions(connectId, conn_names, conn_vals);
  ERRORR(rval, "Failed to get dimensions for connectivity.");

  // Read connectivity into temporary variable
  int num_fine_quads = 0;
  int num_coarse_quads = 0;
  int start_idx = 0;
  std::vector<std::string>::iterator vit;
  int idx = 0;
  if ((vit = std::find(conn_names.begin(), conn_names.end(), "ncells")) != conn_names.end())
    idx = vit - conn_names.begin();
  else if ((vit = std::find(conn_names.begin(), conn_names.end(), "ncenters")) != conn_names.end())
    idx = vit - conn_names.begin();
  else {
    ERRORR(MB_FAILURE, "Failed to get number of quads.");
  }
  int num_quads = conn_vals[idx];
  if (!isConnFile && num_quads != nCells) {
    dbgOut.tprintf(1, "Warning: number of quads from %s and cells from %s are inconsistent; num_quads = %d, nCells = %d.\n",
        conn_fname.c_str(), fileName.c_str(), num_quads, nCells);
  }

  // Get the connectivity into tmp_conn2 and permute into tmp_conn
  int cornerVarId;
  success = NCFUNC(inq_varid)(connectId, "element_corners", &cornerVarId);
  ERRORS(success, "Failed to get variable id.");
  NCDF_SIZE tmp_starts[2] = {0, 0};
  NCDF_SIZE tmp_counts[2] = {4, static_cast<NCDF_SIZE>(num_quads)};
  std::vector<int> tmp_conn(4 * num_quads), tmp_conn2(4 * num_quads);
  success = NCFUNCAG(_vara_int)(connectId, cornerVarId, tmp_starts, tmp_counts, &tmp_conn2[0]);
  ERRORS(success, "Failed to get temporary connectivity.");
  if (isConnFile) {
    // This data/connectivity file will be closed later in ReadNC::load_file()
  }
  else {
    success = NCFUNC(close)(connectId);
    ERRORS(success, "Failed on close.");
  }
  // Permute the connectivity
  for (int i = 0; i < num_quads; i++) {
    tmp_conn[4 * i] = tmp_conn2[i];
    tmp_conn[4 * i + 1] = tmp_conn2[i + 1 * num_quads];
    tmp_conn[4 * i + 2] = tmp_conn2[i + 2 * num_quads];
    tmp_conn[4 * i + 3] = tmp_conn2[i + 3 * num_quads];
  }

  // Need to know whether we'll be creating gather mesh later, to make sure we allocate enough space
  // in one shot
  bool create_gathers = false;
  if (rank == gatherSetRank)
    create_gathers = true;

  // Compute the number of local quads, accounting for coarse or fine representation
  // spectral_unit is the # fine quads per coarse quad, or spectralOrder^2
  int spectral_unit = (spectralMesh ? _spectralOrder * _spectralOrder : 1);
  // num_coarse_quads is the number of quads instantiated in MOAB; if !spectralMesh, num_coarse_quads = num_fine_quads
  num_coarse_quads = int(std::floor(1.0 * num_quads / (spectral_unit * procs)));
  // start_idx is the starting index in the HommeMapping connectivity list for this proc, before converting to coarse quad representation
  start_idx = 4 * rank * num_coarse_quads * spectral_unit;
  // iextra = # coarse quads extra after equal split over procs
  int iextra = num_quads % (procs * spectral_unit);
  if (rank < iextra)
    num_coarse_quads++;
  start_idx += 4 * spectral_unit * std::min(rank, iextra);
  // num_fine_quads is the number of quads in the connectivity list in HommeMapping file assigned to this proc
  num_fine_quads = spectral_unit * num_coarse_quads;

  // Now create num_coarse_quads
  EntityHandle* conn_arr;
  EntityHandle start_vertex;
  Range tmp_range;

  // Read connectivity into that space
  EntityHandle* sv_ptr = NULL;
  EntityHandle start_quad;
  SpectralMeshTool smt(mbImpl, _spectralOrder);
  if (!spectralMesh) {
    rval = _readNC->readMeshIface->get_element_connect(num_coarse_quads, 4,
                                              MBQUAD, 0, start_quad, conn_arr,
                                              // Might have to create gather mesh later
                                              (create_gathers ? num_coarse_quads + num_quads : num_coarse_quads));
    ERRORR(rval, "Failed to create quads.");
    tmp_range.insert(start_quad, start_quad + num_coarse_quads - 1);
    std::copy(&tmp_conn[start_idx], &tmp_conn[start_idx + 4 * num_fine_quads], conn_arr);
    std::copy(conn_arr, conn_arr + 4 * num_fine_quads, range_inserter(localGidVerts));
  }
  else {
    rval = smt.create_spectral_elems(&tmp_conn[0], num_fine_quads, 2, tmp_range, start_idx, &localGidVerts);
    ERRORR(rval, "Failed to create spectral elements.");
    int count, v_per_e;
    rval = mbImpl->connect_iterate(tmp_range.begin(), tmp_range.end(), conn_arr, v_per_e, count);
    ERRORR(rval, "Failed to get connectivity of spectral elements.");
    rval = mbImpl->tag_iterate(smt.spectral_vertices_tag(true), tmp_range.begin(), tmp_range.end(),
                               count, (void*&)sv_ptr);
    ERRORR(rval, "Failed to get fine connectivity of spectral elements.");
  }

  // Create vertices
  nLocalVertices = localGidVerts.size();
  std::vector<double*> arrays;
  rval = _readNC->readMeshIface->get_node_coords(3, nLocalVertices, 0, start_vertex, arrays,
                                        // Might have to create gather mesh later
                                        (create_gathers ? nLocalVertices + nVertices : nLocalVertices));
  ERRORR(rval, "Couldn't create vertices in HOMME mesh.");

  // Set vertex coordinates
  Range::iterator rit;
  double* xptr = arrays[0];
  double* yptr = arrays[1];
  double* zptr = arrays[2];
  int i;
  for (i = 0, rit = localGidVerts.begin(); i < nLocalVertices; i++, ++rit) {
    assert(*rit < xVertVals.size() + 1);
    xptr[i] = xVertVals[(*rit) - 1]; // lon
    yptr[i] = yVertVals[(*rit) - 1]; // lat
  }

  // Convert lon/lat/rad to x/y/z
  const double pideg = acos(-1.0) / 180.0;
  double rad = (isConnFile) ? 8000.0 : 8000.0 + levVals[0];
  for (i = 0; i < nLocalVertices; i++) {
    double cosphi = cos(pideg * yptr[i]);
    double zmult = sin(pideg * yptr[i]);
    double xmult = cosphi * cos(xptr[i] * pideg);
    double ymult = cosphi * sin(xptr[i] * pideg);
    xptr[i] = rad * xmult;
    yptr[i] = rad * ymult;
    zptr[i] = rad * zmult;
  }

  // Get ptr to gid memory for vertices
  Range vert_range(start_vertex, start_vertex + nLocalVertices - 1);
  void* data;
  int count;
  rval = mbImpl->tag_iterate(mGlobalIdTag, vert_range.begin(), vert_range.end(), count, data);
  ERRORR(rval, "Failed to get tag iterator.");
  assert(count == nLocalVertices);
  int* gid_data = (int*) data;
  std::copy(localGidVerts.begin(), localGidVerts.end(), gid_data);

  // Duplicate global id data, which will be used to resolve sharing
  if (mpFileIdTag) {
    rval = mbImpl->tag_iterate(*mpFileIdTag, vert_range.begin(), vert_range.end(), count, data);
    ERRORR(rval, "Failed to get tag iterator on file id tag.");
    assert(count == nLocalVertices);
    gid_data = (int*) data;
    std::copy(localGidVerts.begin(), localGidVerts.end(), gid_data);
  }

  // Create map from file ids to vertex handles, used later to set connectivity
  std::map<EntityHandle, EntityHandle> vert_handles;
  for (rit = localGidVerts.begin(), i = 0; rit != localGidVerts.end(); ++rit, i++)
    vert_handles[*rit] = start_vertex + i;

  // Compute proper handles in connectivity using offset
  for (int q = 0; q < 4 * num_coarse_quads; q++) {
    conn_arr[q] = vert_handles[conn_arr[q]];
    assert(conn_arr[q]);
  }
  if (spectralMesh) {
    int verts_per_quad = (_spectralOrder + 1) * (_spectralOrder + 1);
    for (int q = 0; q < verts_per_quad * num_coarse_quads; q++) {
      sv_ptr[q] = vert_handles[sv_ptr[q]];
      assert(sv_ptr[q]);
    }
  }

  // Add new vertices and elements to the set
  faces.merge(tmp_range);
  tmp_range.insert(start_vertex, start_vertex + nLocalVertices - 1);
  rval = mbImpl->add_entities(_fileSet, tmp_range);
  ERRORR(rval, "Couldn't add new vertices and quads to file set.");

  // Mark the set with the spectral order
  Tag sporder;
  rval = mbImpl->tag_get_handle("SPECTRAL_ORDER", 1, MB_TYPE_INTEGER, sporder, MB_TAG_CREAT | MB_TAG_SPARSE);
  ERRORR(rval, "Couldn't create spectral order tag.");
  rval = mbImpl->tag_set_data(sporder, &_fileSet, 1, &_spectralOrder);
  ERRORR(rval, "Couldn't set value for spectral order tag.");

  if (create_gathers) {
    EntityHandle gather_set;
    rval = _readNC->readMeshIface->create_gather_set(gather_set);
    ERRORR(rval, "Trouble creating gather set.");

    // Create vertices
    arrays.clear();
    // Don't need to specify allocation number here, because we know enough verts were created before
    rval = _readNC->readMeshIface->get_node_coords(3, nVertices, 0, start_vertex, arrays);
    ERRORR(rval, "Couldn't create vertices in HOMME mesh for gather set.");

    xptr = arrays[0];
    yptr = arrays[1];
    zptr = arrays[2];
    for (i = 0; i < nVertices; i++) {
      double cosphi = cos(pideg * yVertVals[i]);
      double zmult = sin(pideg * yVertVals[i]);
      double xmult = cosphi * cos(xVertVals[i] * pideg);
      double ymult = cosphi * sin(xVertVals[i] * pideg);
      xptr[i] = rad * xmult;
      yptr[i] = rad * ymult;
      zptr[i] = rad * zmult;
    }

    // Get ptr to gid memory for vertices
    Range gather_verts(start_vertex, start_vertex + nVertices - 1);
    rval = mbImpl->tag_iterate(mGlobalIdTag, gather_verts.begin(), gather_verts.end(), count, data);
    ERRORR(rval, "Failed to get tag iterator.");
    assert(count == nVertices);
    gid_data = (int*) data;
    for (int j = 1; j <= nVertices; j++)
      gid_data[j - 1] = j;
    // Set the file id tag too, it should be bigger something not interfering with global id
    if (mpFileIdTag) {
      rval = mbImpl->tag_iterate(*mpFileIdTag, gather_verts.begin(), gather_verts.end(), count, data);
      ERRORR(rval, "Failed to get tag iterator in file id tag.");
      assert(count == nVertices);
      gid_data = (int*) data;
      for (int j = 1; j <= nVertices; j++)
        gid_data[j - 1] = nVertices + j; // bigger than global id tag
    }

    rval = mbImpl->add_entities(gather_set, gather_verts);
    ERRORR(rval, "Couldn't add vertices to gather set.");

    // Create quads
    Range gather_quads;
    // Don't need to specify allocation number here, because we know enough quads were created before
    rval = _readNC->readMeshIface->get_element_connect(num_quads, 4,
                                              MBQUAD, 0, start_quad, conn_arr);
    ERRORR(rval, "Failed to create quads.");
    gather_quads.insert(start_quad, start_quad + num_quads - 1);
    std::copy(&tmp_conn[0], &tmp_conn[4 * num_quads], conn_arr);
    for (i = 0; i != 4 * num_quads; i++)
      conn_arr[i] += start_vertex - 1; // Connectivity array is shifted by where the gather verts start
    rval = mbImpl->add_entities(gather_set, gather_quads);
    ERRORR(rval, "Couldn't add quads to gather set.");
  }

  return MB_SUCCESS;
}
virtual std::string moab::NCHelperHOMME::get_mesh_type_name ( ) [inline, private, virtual]

Implementation of NCHelper::get_mesh_type_name()

Implements moab::NCHelper.

Definition at line 31 of file NCHelperHOMME.hpp.

{ return "CAM_SE"; }

Implementation of NCHelper::init_mesh_vals()

Implements moab::NCHelper.

Definition at line 73 of file NCHelperHOMME.cpp.

{
  std::vector<std::string>& dimNames = _readNC->dimNames;
  std::vector<int>& dimLens = _readNC->dimLens;
  std::map<std::string, ReadNC::VarData>& varInfo = _readNC->varInfo;

  ErrorCode rval;
  unsigned int idx;
  std::vector<std::string>::iterator vit;

  // Look for time dimension
  if (isConnFile) {
    // Connectivity file might not have time dimension
  }
  else {
    if ((vit = std::find(dimNames.begin(), dimNames.end(), "time")) != dimNames.end())
      idx = vit - dimNames.begin();
    else if ((vit = std::find(dimNames.begin(), dimNames.end(), "t")) != dimNames.end())
      idx = vit - dimNames.begin();
    else {
      ERRORR(MB_FAILURE, "Couldn't find 'time' or 't' dimension.");
    }
    tDim = idx;
    nTimeSteps = dimLens[idx];
  }

  // Get number of vertices (labeled as number of columns)
  if ((vit = std::find(dimNames.begin(), dimNames.end(), "ncol")) != dimNames.end())
    idx = vit - dimNames.begin();
  else {
    ERRORR(MB_FAILURE, "Couldn't find 'ncol' dimension.");
  }
  vDim = idx;
  nVertices = dimLens[idx];

  // Set number of cells
  nCells = nVertices - 2;

  // Get number of levels
  if (isConnFile) {
    // Connectivity file might not have level dimension
  }
  else {
    if ((vit = std::find(dimNames.begin(), dimNames.end(), "lev")) != dimNames.end())
      idx = vit - dimNames.begin();
    else if ((vit = std::find(dimNames.begin(), dimNames.end(), "ilev")) != dimNames.end())
      idx = vit - dimNames.begin();
    else {
      ERRORR(MB_FAILURE, "Couldn't find 'lev' or 'ilev' dimension.");
    }
    levDim = idx;
    nLevels = dimLens[idx];
  }

  // Store lon values in xVertVals
  std::map<std::string, ReadNC::VarData>::iterator vmit;
  if ((vmit = varInfo.find("lon")) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
    rval = read_coordinate("lon", 0, nVertices - 1, xVertVals);
    ERRORR(rval, "Trouble reading 'lon' variable.");
  }
  else {
    ERRORR(MB_FAILURE, "Couldn't find 'lon' variable.");
  }

  // Store lat values in yVertVals
  if ((vmit = varInfo.find("lat")) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
    rval = read_coordinate("lat", 0, nVertices - 1, yVertVals);
    ERRORR(rval, "Trouble reading 'lat' variable.");
  }
  else {
    ERRORR(MB_FAILURE, "Couldn't find 'lat' variable.");
  }

  // Store lev values in levVals
  if (isConnFile) {
    // Connectivity file might not have level variable
  }
  else {
    if ((vmit = varInfo.find("lev")) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
      rval = read_coordinate("lev", 0, nLevels - 1, levVals);
      ERRORR(rval, "Trouble reading 'lev' variable.");

      // Decide whether down is positive
      char posval[10];
      int success = NCFUNC(get_att_text)(_fileId, (*vmit).second.varId, "positive", posval);
      if (0 == success && !strcmp(posval, "down")) {
        for (std::vector<double>::iterator dvit = levVals.begin(); dvit != levVals.end(); ++dvit)
          (*dvit) *= -1.0;
      }
    }
    else {
      ERRORR(MB_FAILURE, "Couldn't find 'lev' variable.");
    }
  }

  // Store time coordinate values in tVals
  if (isConnFile) {
    // Connectivity file might not have time variable
  }
  else {
    if ((vmit = varInfo.find("time")) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
      rval = read_coordinate("time", 0, nTimeSteps - 1, tVals);
      ERRORR(rval, "Trouble reading 'time' variable.");
    }
    else if ((vmit = varInfo.find("t")) != varInfo.end() && (*vmit).second.varDims.size() == 1) {
      rval = read_coordinate("t", 0, nTimeSteps - 1, tVals);
      ERRORR(rval, "Trouble reading 't' variable.");
    }
    else {
      // If expected time variable is not available, set dummy time coordinate values to tVals
      for (int t = 0; t < nTimeSteps; t++)
        tVals.push_back((double)t);
    }
  }

  // For each variable, determine the entity location type and number of levels
  std::map<std::string, ReadNC::VarData>::iterator mit;
  for (mit = varInfo.begin(); mit != varInfo.end(); ++mit) {
    ReadNC::VarData& vd = (*mit).second;

    vd.entLoc = ReadNC::ENTLOCSET;
    if (std::find(vd.varDims.begin(), vd.varDims.end(), tDim) != vd.varDims.end()) {
      if (std::find(vd.varDims.begin(), vd.varDims.end(), vDim) != vd.varDims.end())
        vd.entLoc = ReadNC::ENTLOCVERT;
    }

    vd.numLev = 1;
    if (std::find(vd.varDims.begin(), vd.varDims.end(), levDim) != vd.varDims.end())
      vd.numLev = nLevels;
  }

  // Hack: create dummy variables for dimensions (like ncol) with no corresponding coordinate variables
  rval = create_dummy_variables();
  ERRORR(rval, "Failed to create dummy variables.");

  return MB_SUCCESS;
}
ErrorCode moab::NCHelperHOMME::read_ucd_variable_to_nonset ( std::vector< ReadNC::VarData > &  vdatas,
std::vector< int > &  tstep_nums 
) [private, virtual]

Implementation of UcdNCHelper::read_ucd_variable_to_nonset()

Implements moab::UcdNCHelper.

Definition at line 799 of file NCHelperHOMME.cpp.

{
  DebugOutput& dbgOut = _readNC->dbgOut;

  ErrorCode rval = read_ucd_variable_to_nonset_allocate(vdatas, tstep_nums);
  ERRORR(rval, "Trouble allocating read variables.");

  // Finally, read into that space
  int success;
  for (unsigned int i = 0; i < vdatas.size(); i++) {
    std::size_t sz = vdatas[i].sz;

    // A typical supported variable: float T(time, lev, ncol)
    // For tag values, need transpose (lev, ncol) to (ncol, lev)
    size_t ni = vdatas[i].readCounts[2]; // ncol
    size_t nj = 1; // Here we should just set nj to 1
    size_t nk = vdatas[i].readCounts[1]; // lev

    for (unsigned int t = 0; t < tstep_nums.size(); t++) {
      // Tag data for this timestep
      void* data = vdatas[i].varDatas[t];

      // Set readStart for each timestep along time dimension
      vdatas[i].readStarts[0] = tstep_nums[t];

      switch (vdatas[i].varDataType) {
        case NC_BYTE:
        case NC_CHAR: {
          ERRORR(MB_FAILURE, "not implemented");
          break;
        }
        case NC_DOUBLE: {
          // Copied from float case
          std::vector<double> tmpdoubledata(sz);

          // In the case of ucd mesh, and on multiple proc,
          // we need to read as many times as subranges we have in the
          // localGidVerts range;
          // basically, we have to give a different point
          // for data to start, for every subrange :(
          size_t indexInDoubleArray = 0;
          size_t ic = 0;
          for (Range::pair_iterator pair_iter = localGidVerts.pair_begin();
              pair_iter != localGidVerts.pair_end();
              pair_iter++, ic++) {
            EntityHandle starth = pair_iter->first;
            EntityHandle endh = pair_iter->second; // Inclusive
            vdatas[i].readStarts[2] = (NCDF_SIZE) (starth - 1);
            vdatas[i].readCounts[2] = (NCDF_SIZE) (endh - starth + 1);

            success = NCFUNCAG(_vara_double)(_fileId, vdatas[i].varId,
                            &(vdatas[i].readStarts[0]), &(vdatas[i].readCounts[0]),
                            &(tmpdoubledata[indexInDoubleArray]));
            ERRORS(success, "Failed to read float data in loop");
            // We need to increment the index in double array for the
            // next subrange
            indexInDoubleArray += (endh - starth + 1) * 1 * vdatas[i].numLev;
          }
          assert(ic == localGidVerts.psize());

          if (vdatas[i].numLev != 1)
            // Transpose (lev, ncol) to (ncol, lev)
            success = kji_to_jik_stride(ni, nj, nk, data, &tmpdoubledata[0], localGidVerts);
          else {
            for (std::size_t idx = 0; idx != tmpdoubledata.size(); idx++)
              ((double*) data)[idx] = tmpdoubledata[idx];
          }
          ERRORS(success, "Failed to read double data.");
          break;
        }
        case NC_FLOAT: {
          std::vector<float> tmpfloatdata(sz);

          // In the case of ucd mesh, and on multiple proc,
          // we need to read as many times as subranges we have in the
          // localGidVerts range;
          // basically, we have to give a different point
          // for data to start, for every subrange :(
          size_t indexInFloatArray = 0;
          size_t ic = 0;
          for (Range::pair_iterator pair_iter = localGidVerts.pair_begin();
              pair_iter != localGidVerts.pair_end();
              pair_iter++, ic++) {
            EntityHandle starth = pair_iter->first;
            EntityHandle endh = pair_iter->second; // Inclusive
            vdatas[i].readStarts[2] = (NCDF_SIZE) (starth - 1);
            vdatas[i].readCounts[2] = (NCDF_SIZE) (endh - starth + 1);

            success = NCFUNCAG(_vara_float)(_fileId, vdatas[i].varId,
                            &(vdatas[i].readStarts[0]), &(vdatas[i].readCounts[0]),
                            &(tmpfloatdata[indexInFloatArray]));
            ERRORS(success, "Failed to read float data in loop");
            // We need to increment the index in float array for the
            // next subrange
            indexInFloatArray += (endh - starth + 1) * 1 * vdatas[i].numLev;
          }
          assert(ic == localGidVerts.psize());

          if (vdatas[i].numLev != 1)
            // Transpose (lev, ncol) to (ncol, lev)
            success = kji_to_jik_stride(ni, nj, nk, data, &tmpfloatdata[0], localGidVerts);
          else {
            for (std::size_t idx = 0; idx != tmpfloatdata.size(); idx++)
              ((float*) data)[idx] = tmpfloatdata[idx];
          }
          ERRORS(success, "Failed to read float data.");
          break;
        }
        case NC_INT: {
          ERRORR(MB_FAILURE, "not implemented");
          break;
        }
        case NC_SHORT: {
          ERRORR(MB_FAILURE, "not implemented");
          break;
        }
        default:
          success = 1;
      }

      if (success)
        ERRORR(MB_FAILURE, "Trouble reading variable.");
    }
  }

  for (unsigned int i = 0; i < vdatas.size(); i++) {
    for (unsigned int t = 0; t < tstep_nums.size(); t++) {
      dbgOut.tprintf(2, "Converting variable %s, time step %d\n", vdatas[i].varName.c_str(), tstep_nums[t]);
      ErrorCode tmp_rval = convert_variable(vdatas[i], t);
      if (MB_SUCCESS != tmp_rval)
        rval = tmp_rval;
    }
  }
  // Debug output, if requested
  if (1 == dbgOut.get_verbosity()) {
    dbgOut.printf(1, "Read variables: %s", vdatas.begin()->varName.c_str());
    for (unsigned int i = 1; i < vdatas.size(); i++)
      dbgOut.printf(1, ", %s ", vdatas[i].varName.c_str());
    dbgOut.tprintf(1, "\n");
  }

  return rval;
}
ErrorCode moab::NCHelperHOMME::read_ucd_variable_to_nonset_allocate ( std::vector< ReadNC::VarData > &  vdatas,
std::vector< int > &  tstep_nums 
) [private, virtual]

Implementation of UcdNCHelper::read_ucd_variable_to_nonset_allocate()

Implements moab::UcdNCHelper.

Definition at line 554 of file NCHelperHOMME.cpp.

{
  Interface*& mbImpl = _readNC->mbImpl;
  std::vector<int>& dimLens = _readNC->dimLens;
  DebugOutput& dbgOut = _readNC->dbgOut;

  ErrorCode rval = MB_SUCCESS;

  Range* range = NULL;

  // Get vertices in set
  Range verts;
  rval = mbImpl->get_entities_by_dimension(_fileSet, 0, verts);
  ERRORR(rval, "Trouble getting vertices in set.");
  assert("Should only have a single vertex subrange, since they were read in one shot" &&
      verts.psize() == 1);

  for (unsigned int i = 0; i < vdatas.size(); i++) {
    // Support non-set variables with 3 dimensions like (time, lev, ncol)
    assert(3 == vdatas[i].varDims.size());

    // For a non-set variable, time should be the first dimension
    assert(tDim == vdatas[i].varDims[0]);

    // Set up readStarts and readCounts
    vdatas[i].readStarts.resize(3);
    vdatas[i].readCounts.resize(3);

    // First: time
    vdatas[i].readStarts[0] = 0; // This value is timestep dependent, will be set later
    vdatas[i].readCounts[0] = 1;

    // Next: lev
    vdatas[i].readStarts[1] = 0;
    vdatas[i].readCounts[1] = vdatas[i].numLev;

    // Finally: ncol
    switch (vdatas[i].entLoc) {
      case ReadNC::ENTLOCVERT:
        // Vertices
        // Start from the first localGidVerts
        // Actually, this will be reset later on in a loop
        vdatas[i].readStarts[2] = localGidVerts[0] - 1;
        vdatas[i].readCounts[2] = nLocalVertices;
        range = &verts;
        break;
      default:
        ERRORR(MB_FAILURE, "Unexpected entity location type for HOMME non-set variable.");
        break;
    }

    // Get variable size
    vdatas[i].sz = 1;
    for (std::size_t idx = 0; idx != 3; idx++)
      vdatas[i].sz *= vdatas[i].readCounts[idx];

    for (unsigned int t = 0; t < tstep_nums.size(); t++) {
      dbgOut.tprintf(2, "Reading variable %s, time step %d\n", vdatas[i].varName.c_str(), tstep_nums[t]);

      if (tstep_nums[t] >= dimLens[tDim]) {
        ERRORR(MB_INDEX_OUT_OF_RANGE, "Wrong value for a timestep number.");
      }

      // Get the tag to read into
      if (!vdatas[i].varTags[t]) {
        rval = get_tag_to_nonset(vdatas[i], tstep_nums[t], vdatas[i].varTags[t], vdatas[i].numLev);
        ERRORR(rval, "Trouble getting tag.");
      }

      // Get ptr to tag space
      void* data;
      int count;
      rval = mbImpl->tag_iterate(vdatas[i].varTags[t], range->begin(), range->end(), count, data);
      ERRORR(rval, "Failed to get tag iterator.");
      assert((unsigned)count == range->size());
      vdatas[i].varDatas[t] = data;
    }
  }

  return rval;
}

Member Data Documentation

Definition at line 47 of file NCHelperHOMME.hpp.

Definition at line 48 of file NCHelperHOMME.hpp.

Definition at line 49 of file NCHelperHOMME.hpp.


The documentation for this class was generated from the following files:
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines