|
moab
|
Child helper class for HOMME grid (CAM_SE) More...
#include <NCHelperHOMME.hpp>
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 |
Child helper class for HOMME grid (CAM_SE)
Definition at line 17 of file NCHelperHOMME.hpp.
| 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 } }
| 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;
}
| ErrorCode moab::NCHelperHOMME::check_existing_mesh | ( | ) | [private, virtual] |
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"; }
| ErrorCode moab::NCHelperHOMME::init_mesh_vals | ( | ) | [private, virtual] |
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;
}
int moab::NCHelperHOMME::_spectralOrder [private] |
Definition at line 47 of file NCHelperHOMME.hpp.
int moab::NCHelperHOMME::connectId [private] |
Definition at line 48 of file NCHelperHOMME.hpp.
bool moab::NCHelperHOMME::isConnFile [private] |
Definition at line 49 of file NCHelperHOMME.hpp.
1.7.6.1