MOAB: Mesh Oriented datABase  (version 5.2.1)
NCHelperHOMME.cpp
Go to the documentation of this file.
00001 #include "NCHelperHOMME.hpp"
00002 #include "moab/ReadUtilIface.hpp"
00003 #include "moab/FileOptions.hpp"
00004 #include "moab/SpectralMeshTool.hpp"
00005 
00006 #include <cmath>
00007 
00008 namespace moab
00009 {
00010 
00011 NCHelperHOMME::NCHelperHOMME( ReadNC* readNC, int fileId, const FileOptions& opts, EntityHandle fileSet )
00012     : UcdNCHelper( readNC, fileId, opts, fileSet ), _spectralOrder( -1 ), connectId( -1 ), isConnFile( false )
00013 {
00014     // Calculate spectral order
00015     std::map< std::string, ReadNC::AttData >::iterator attIt = readNC->globalAtts.find( "np" );
00016     if( attIt != readNC->globalAtts.end() )
00017     {
00018         int success = NCFUNC( get_att_int )( readNC->fileId, attIt->second.attVarId, attIt->second.attName.c_str(),
00019                                              &_spectralOrder );
00020         if( 0 == success ) _spectralOrder--;  // Spectral order is one less than np
00021     }
00022     else
00023     {
00024         // As can_read_file() returns true and there is no global attribute "np", it should be a
00025         // connectivity file
00026         isConnFile     = true;
00027         _spectralOrder = 3;  // Assume np is 4
00028     }
00029 }
00030 
00031 bool NCHelperHOMME::can_read_file( ReadNC* readNC, int fileId )
00032 {
00033     // If global attribute "np" exists then it should be the HOMME grid
00034     if( readNC->globalAtts.find( "np" ) != readNC->globalAtts.end() )
00035     {
00036         // Make sure it is CAM grid
00037         std::map< std::string, ReadNC::AttData >::iterator attIt = readNC->globalAtts.find( "source" );
00038         if( attIt == readNC->globalAtts.end() ) return false;
00039         unsigned int sz = attIt->second.attLen;
00040         std::string att_data;
00041         att_data.resize( sz + 1 );
00042         att_data[sz] = '\000';
00043         int success =
00044             NCFUNC( get_att_text )( fileId, attIt->second.attVarId, attIt->second.attName.c_str(), &att_data[0] );
00045         if( success ) return false;
00046         if( att_data.find( "CAM" ) == std::string::npos ) return false;
00047 
00048         return true;
00049     }
00050     else
00051     {
00052         // If dimension names "ncol" AND "ncorners" AND "ncells" exist, then it should be the HOMME
00053         // connectivity file In this case, the mesh can still be created
00054         std::vector< std::string >& dimNames = readNC->dimNames;
00055         if( ( std::find( dimNames.begin(), dimNames.end(), std::string( "ncol" ) ) != dimNames.end() ) &&
00056             ( std::find( dimNames.begin(), dimNames.end(), std::string( "ncorners" ) ) != dimNames.end() ) &&
00057             ( std::find( dimNames.begin(), dimNames.end(), std::string( "ncells" ) ) != dimNames.end() ) )
00058             return true;
00059     }
00060 
00061     return false;
00062 }
00063 
00064 ErrorCode NCHelperHOMME::init_mesh_vals()
00065 {
00066     std::vector< std::string >& dimNames              = _readNC->dimNames;
00067     std::vector< int >& dimLens                       = _readNC->dimLens;
00068     std::map< std::string, ReadNC::VarData >& varInfo = _readNC->varInfo;
00069 
00070     ErrorCode rval;
00071     unsigned int idx;
00072     std::vector< std::string >::iterator vit;
00073 
00074     // Look for time dimension
00075     if( isConnFile )
00076     {
00077         // Connectivity file might not have time dimension
00078     }
00079     else
00080     {
00081         if( ( vit = std::find( dimNames.begin(), dimNames.end(), "time" ) ) != dimNames.end() )
00082             idx = vit - dimNames.begin();
00083         else if( ( vit = std::find( dimNames.begin(), dimNames.end(), "t" ) ) != dimNames.end() )
00084             idx = vit - dimNames.begin();
00085         else
00086         {
00087             MB_SET_ERR( MB_FAILURE, "Couldn't find 'time' or 't' dimension" );
00088         }
00089         tDim       = idx;
00090         nTimeSteps = dimLens[idx];
00091     }
00092 
00093     // Get number of vertices (labeled as number of columns)
00094     if( ( vit = std::find( dimNames.begin(), dimNames.end(), "ncol" ) ) != dimNames.end() )
00095         idx = vit - dimNames.begin();
00096     else
00097     {
00098         MB_SET_ERR( MB_FAILURE, "Couldn't find 'ncol' dimension" );
00099     }
00100     vDim      = idx;
00101     nVertices = dimLens[idx];
00102 
00103     // Set number of cells
00104     nCells = nVertices - 2;
00105 
00106     // Get number of levels
00107     if( isConnFile )
00108     {
00109         // Connectivity file might not have level dimension
00110     }
00111     else
00112     {
00113         if( ( vit = std::find( dimNames.begin(), dimNames.end(), "lev" ) ) != dimNames.end() )
00114             idx = vit - dimNames.begin();
00115         else if( ( vit = std::find( dimNames.begin(), dimNames.end(), "ilev" ) ) != dimNames.end() )
00116             idx = vit - dimNames.begin();
00117         else
00118         {
00119             MB_SET_ERR( MB_FAILURE, "Couldn't find 'lev' or 'ilev' dimension" );
00120         }
00121         levDim  = idx;
00122         nLevels = dimLens[idx];
00123     }
00124 
00125     // Store lon values in xVertVals
00126     std::map< std::string, ReadNC::VarData >::iterator vmit;
00127     if( ( vmit = varInfo.find( "lon" ) ) != varInfo.end() && ( *vmit ).second.varDims.size() == 1 )
00128     {
00129         rval = read_coordinate( "lon", 0, nVertices - 1, xVertVals );MB_CHK_SET_ERR( rval, "Trouble reading 'lon' variable" );
00130     }
00131     else
00132     {
00133         MB_SET_ERR( MB_FAILURE, "Couldn't find 'lon' variable" );
00134     }
00135 
00136     // Store lat values in yVertVals
00137     if( ( vmit = varInfo.find( "lat" ) ) != varInfo.end() && ( *vmit ).second.varDims.size() == 1 )
00138     {
00139         rval = read_coordinate( "lat", 0, nVertices - 1, yVertVals );MB_CHK_SET_ERR( rval, "Trouble reading 'lat' variable" );
00140     }
00141     else
00142     {
00143         MB_SET_ERR( MB_FAILURE, "Couldn't find 'lat' variable" );
00144     }
00145 
00146     // Store lev values in levVals
00147     if( isConnFile )
00148     {
00149         // Connectivity file might not have level variable
00150     }
00151     else
00152     {
00153         if( ( vmit = varInfo.find( "lev" ) ) != varInfo.end() && ( *vmit ).second.varDims.size() == 1 )
00154         {
00155             rval = read_coordinate( "lev", 0, nLevels - 1, levVals );MB_CHK_SET_ERR( rval, "Trouble reading 'lev' variable" );
00156 
00157             // Decide whether down is positive
00158             char posval[10] = { 0 };
00159             int success     = NCFUNC( get_att_text )( _fileId, ( *vmit ).second.varId, "positive", posval );
00160             if( 0 == success && !strcmp( posval, "down" ) )
00161             {
00162                 for( std::vector< double >::iterator dvit = levVals.begin(); dvit != levVals.end(); ++dvit )
00163                     ( *dvit ) *= -1.0;
00164             }
00165         }
00166         else
00167         {
00168             MB_SET_ERR( MB_FAILURE, "Couldn't find 'lev' variable" );
00169         }
00170     }
00171 
00172     // Store time coordinate values in tVals
00173     if( isConnFile )
00174     {
00175         // Connectivity file might not have time variable
00176     }
00177     else
00178     {
00179         if( ( vmit = varInfo.find( "time" ) ) != varInfo.end() && ( *vmit ).second.varDims.size() == 1 )
00180         {
00181             rval = read_coordinate( "time", 0, nTimeSteps - 1, tVals );MB_CHK_SET_ERR( rval, "Trouble reading 'time' variable" );
00182         }
00183         else if( ( vmit = varInfo.find( "t" ) ) != varInfo.end() && ( *vmit ).second.varDims.size() == 1 )
00184         {
00185             rval = read_coordinate( "t", 0, nTimeSteps - 1, tVals );MB_CHK_SET_ERR( rval, "Trouble reading 't' variable" );
00186         }
00187         else
00188         {
00189             // If expected time variable does not exist, set dummy values to tVals
00190             for( int t = 0; t < nTimeSteps; t++ )
00191                 tVals.push_back( (double)t );
00192         }
00193     }
00194 
00195     // For each variable, determine the entity location type and number of levels
00196     std::map< std::string, ReadNC::VarData >::iterator mit;
00197     for( mit = varInfo.begin(); mit != varInfo.end(); ++mit )
00198     {
00199         ReadNC::VarData& vd = ( *mit ).second;
00200 
00201         // Default entLoc is ENTLOCSET
00202         if( std::find( vd.varDims.begin(), vd.varDims.end(), tDim ) != vd.varDims.end() )
00203         {
00204             if( std::find( vd.varDims.begin(), vd.varDims.end(), vDim ) != vd.varDims.end() )
00205                 vd.entLoc = ReadNC::ENTLOCVERT;
00206         }
00207 
00208         // Default numLev is 0
00209         if( std::find( vd.varDims.begin(), vd.varDims.end(), levDim ) != vd.varDims.end() ) vd.numLev = nLevels;
00210     }
00211 
00212     // Hack: create dummy variables for dimensions (like ncol) with no corresponding coordinate
00213     // variables
00214     rval = create_dummy_variables();MB_CHK_SET_ERR( rval, "Failed to create dummy variables" );
00215 
00216     return MB_SUCCESS;
00217 }
00218 
00219 // When noMesh option is used on this read, the old ReadNC class instance for last read can get out
00220 // of scope (and deleted). The old instance initialized localGidVerts properly when the mesh was
00221 // created, but it is now lost. The new instance (will not create the mesh with noMesh option) has
00222 // to restore it based on the existing mesh from last read
00223 ErrorCode NCHelperHOMME::check_existing_mesh()
00224 {
00225     Interface*& mbImpl = _readNC->mbImpl;
00226     Tag& mGlobalIdTag  = _readNC->mGlobalIdTag;
00227     bool& noMesh       = _readNC->noMesh;
00228 
00229     if( noMesh && localGidVerts.empty() )
00230     {
00231         // We need to populate localGidVerts range with the gids of vertices from current file set
00232         // localGidVerts is important in reading the variable data into the nodes
00233         // Also, for our purposes, localGidVerts is truly the GLOBAL_ID tag data, not other
00234         // file_id tags that could get passed around in other scenarios for parallel reading
00235 
00236         // Get all vertices from current file set (it is the input set in no_mesh scenario)
00237         Range local_verts;
00238         ErrorCode rval = mbImpl->get_entities_by_dimension( _fileSet, 0, local_verts );MB_CHK_SET_ERR( rval, "Trouble getting local vertices in current file set" );
00239 
00240         if( !local_verts.empty() )
00241         {
00242             std::vector< int > gids( local_verts.size() );
00243 
00244             // !IMPORTANT : this has to be the GLOBAL_ID tag
00245             rval = mbImpl->tag_get_data( mGlobalIdTag, local_verts, &gids[0] );MB_CHK_SET_ERR( rval, "Trouble getting local gid values of vertices" );
00246 
00247             // Restore localGidVerts
00248             std::copy( gids.rbegin(), gids.rend(), range_inserter( localGidVerts ) );
00249             nLocalVertices = localGidVerts.size();
00250         }
00251     }
00252 
00253     return MB_SUCCESS;
00254 }
00255 
00256 ErrorCode NCHelperHOMME::create_mesh( Range& faces )
00257 {
00258     Interface*& mbImpl         = _readNC->mbImpl;
00259     std::string& fileName      = _readNC->fileName;
00260     Tag& mGlobalIdTag          = _readNC->mGlobalIdTag;
00261     const Tag*& mpFileIdTag    = _readNC->mpFileIdTag;
00262     DebugOutput& dbgOut        = _readNC->dbgOut;
00263     bool& spectralMesh         = _readNC->spectralMesh;
00264     int& gatherSetRank         = _readNC->gatherSetRank;
00265     int& trivialPartitionShift = _readNC->trivialPartitionShift;
00266 
00267     int rank  = 0;
00268     int procs = 1;
00269 #ifdef MOAB_HAVE_MPI
00270     bool& isParallel = _readNC->isParallel;
00271     if( isParallel )
00272     {
00273         ParallelComm*& myPcomm = _readNC->myPcomm;
00274         rank                   = myPcomm->proc_config().proc_rank();
00275         procs                  = myPcomm->proc_config().proc_size();
00276     }
00277 #endif
00278 
00279     ErrorCode rval;
00280     int success = 0;
00281 
00282     // Need to get/read connectivity data before creating elements
00283     std::string conn_fname;
00284 
00285     if( isConnFile )
00286     {
00287         // Connectivity file has already been read
00288         connectId = _readNC->fileId;
00289     }
00290     else
00291     {
00292         // Try to open the connectivity file through CONN option, if used
00293         rval = _opts.get_str_option( "CONN", conn_fname );
00294         if( MB_SUCCESS != rval )
00295         {
00296             // Default convention for reading HOMME is a file HommeMapping.nc in same dir as data
00297             // file
00298             conn_fname = std::string( fileName );
00299             size_t idx = conn_fname.find_last_of( "/" );
00300             if( idx != std::string::npos )
00301                 conn_fname = conn_fname.substr( 0, idx ).append( "/HommeMapping.nc" );
00302             else
00303                 conn_fname = "HommeMapping.nc";
00304         }
00305 #ifdef MOAB_HAVE_PNETCDF
00306 #ifdef MOAB_HAVE_MPI
00307         if( isParallel )
00308         {
00309             ParallelComm*& myPcomm = _readNC->myPcomm;
00310             success =
00311                 NCFUNC( open )( myPcomm->proc_config().proc_comm(), conn_fname.c_str(), 0, MPI_INFO_NULL, &connectId );
00312         }
00313         else
00314             success = NCFUNC( open )( MPI_COMM_SELF, conn_fname.c_str(), 0, MPI_INFO_NULL, &connectId );
00315 #endif
00316 #else
00317         success = NCFUNC( open )( conn_fname.c_str(), 0, &connectId );
00318 #endif
00319         if( success ) MB_SET_ERR( MB_FAILURE, "Failed on open" );
00320     }
00321 
00322     std::vector< std::string > conn_names;
00323     std::vector< int > conn_vals;
00324     rval = _readNC->get_dimensions( connectId, conn_names, conn_vals );MB_CHK_SET_ERR( rval, "Failed to get dimensions for connectivity" );
00325 
00326     // Read connectivity into temporary variable
00327     int num_fine_quads   = 0;
00328     int num_coarse_quads = 0;
00329     int start_idx        = 0;
00330     std::vector< std::string >::iterator vit;
00331     int idx = 0;
00332     if( ( vit = std::find( conn_names.begin(), conn_names.end(), "ncells" ) ) != conn_names.end() )
00333         idx = vit - conn_names.begin();
00334     else if( ( vit = std::find( conn_names.begin(), conn_names.end(), "ncenters" ) ) != conn_names.end() )
00335         idx = vit - conn_names.begin();
00336     else
00337     {
00338         MB_SET_ERR( MB_FAILURE, "Failed to get number of quads" );
00339     }
00340     int num_quads = conn_vals[idx];
00341     if( !isConnFile && num_quads != nCells )
00342     {
00343         dbgOut.tprintf( 1,
00344                         "Warning: number of quads from %s and cells from %s are inconsistent; "
00345                         "num_quads = %d, nCells = %d.\n",
00346                         conn_fname.c_str(), fileName.c_str(), num_quads, nCells );
00347     }
00348 
00349     // Get the connectivity into tmp_conn2 and permute into tmp_conn
00350     int cornerVarId;
00351     success = NCFUNC( inq_varid )( connectId, "element_corners", &cornerVarId );
00352     if( success ) MB_SET_ERR( MB_FAILURE, "Failed to get variable id of 'element_corners'" );
00353     NCDF_SIZE tmp_starts[2] = { 0, 0 };
00354     NCDF_SIZE tmp_counts[2] = { 4, static_cast< NCDF_SIZE >( num_quads ) };
00355     std::vector< int > tmp_conn( 4 * num_quads ), tmp_conn2( 4 * num_quads );
00356     success = NCFUNCAG( _vara_int )( connectId, cornerVarId, tmp_starts, tmp_counts, &tmp_conn2[0] );
00357     if( success ) MB_SET_ERR( MB_FAILURE, "Failed to get temporary connectivity" );
00358     if( isConnFile )
00359     {
00360         // This data/connectivity file will be closed later in ReadNC::load_file()
00361     }
00362     else
00363     {
00364         success = NCFUNC( close )( connectId );
00365         if( success ) MB_SET_ERR( MB_FAILURE, "Failed on close" );
00366     }
00367     // Permute the connectivity
00368     for( int i = 0; i < num_quads; i++ )
00369     {
00370         tmp_conn[4 * i]     = tmp_conn2[i];
00371         tmp_conn[4 * i + 1] = tmp_conn2[i + 1 * num_quads];
00372         tmp_conn[4 * i + 2] = tmp_conn2[i + 2 * num_quads];
00373         tmp_conn[4 * i + 3] = tmp_conn2[i + 3 * num_quads];
00374     }
00375 
00376     // Need to know whether we'll be creating gather mesh later, to make sure
00377     // we allocate enough space in one shot
00378     bool create_gathers = false;
00379     if( rank == gatherSetRank ) create_gathers = true;
00380 
00381     // Shift rank to obtain a rotated trivial partition
00382     int shifted_rank = rank;
00383     if( procs >= 2 && trivialPartitionShift > 0 ) shifted_rank = ( rank + trivialPartitionShift ) % procs;
00384 
00385     // Compute the number of local quads, accounting for coarse or fine representation
00386     // spectral_unit is the # fine quads per coarse quad, or spectralOrder^2
00387     int spectral_unit = ( spectralMesh ? _spectralOrder * _spectralOrder : 1 );
00388     // num_coarse_quads is the number of quads instantiated in MOAB; if !spectralMesh,
00389     // num_coarse_quads = num_fine_quads
00390     num_coarse_quads = int( std::floor( 1.0 * num_quads / ( spectral_unit * procs ) ) );
00391     // start_idx is the starting index in the HommeMapping connectivity list for this proc, before
00392     // converting to coarse quad representation
00393     start_idx = 4 * shifted_rank * num_coarse_quads * spectral_unit;
00394     // iextra = # coarse quads extra after equal split over procs
00395     int iextra = num_quads % ( procs * spectral_unit );
00396     if( shifted_rank < iextra ) num_coarse_quads++;
00397     start_idx += 4 * spectral_unit * std::min( shifted_rank, iextra );
00398     // num_fine_quads is the number of quads in the connectivity list in HommeMapping file assigned
00399     // to this proc
00400     num_fine_quads = spectral_unit * num_coarse_quads;
00401 
00402     // Now create num_coarse_quads
00403     EntityHandle* conn_arr;
00404     EntityHandle start_vertex;
00405     Range tmp_range;
00406 
00407     // Read connectivity into that space
00408     EntityHandle* sv_ptr = NULL;
00409     EntityHandle start_quad;
00410     SpectralMeshTool smt( mbImpl, _spectralOrder );
00411     if( !spectralMesh )
00412     {
00413         rval = _readNC->readMeshIface->get_element_connect( num_coarse_quads, 4, MBQUAD, 0, start_quad, conn_arr,
00414                                                             // Might have to create gather mesh later
00415                                                             ( create_gathers ? num_coarse_quads + num_quads
00416                                                                              : num_coarse_quads ) );MB_CHK_SET_ERR( rval, "Failed to create local quads" );
00417         tmp_range.insert( start_quad, start_quad + num_coarse_quads - 1 );
00418         int* tmp_conn_end = ( &tmp_conn[start_idx + 4 * num_fine_quads - 1] ) + 1;
00419         std::copy( &tmp_conn[start_idx], tmp_conn_end, conn_arr );
00420         std::copy( conn_arr, conn_arr + 4 * num_fine_quads, range_inserter( localGidVerts ) );
00421     }
00422     else
00423     {
00424         rval = smt.create_spectral_elems( &tmp_conn[0], num_fine_quads, 2, tmp_range, start_idx, &localGidVerts );MB_CHK_SET_ERR( rval, "Failed to create spectral elements" );
00425         int count, v_per_e;
00426         rval = mbImpl->connect_iterate( tmp_range.begin(), tmp_range.end(), conn_arr, v_per_e, count );MB_CHK_SET_ERR( rval, "Failed to get connectivity of spectral elements" );
00427         rval = mbImpl->tag_iterate( smt.spectral_vertices_tag( true ), tmp_range.begin(), tmp_range.end(), count,
00428                                     (void*&)sv_ptr );MB_CHK_SET_ERR( rval, "Failed to get fine connectivity of spectral elements" );
00429     }
00430 
00431     // Create vertices
00432     nLocalVertices = localGidVerts.size();
00433     std::vector< double* > arrays;
00434     rval = _readNC->readMeshIface->get_node_coords( 3, nLocalVertices, 0, start_vertex, arrays,
00435                                                     // Might have to create gather mesh later
00436                                                     ( create_gathers ? nLocalVertices + nVertices : nLocalVertices ) );MB_CHK_SET_ERR( rval, "Failed to create local vertices" );
00437 
00438     // Set vertex coordinates
00439     Range::iterator rit;
00440     double* xptr = arrays[0];
00441     double* yptr = arrays[1];
00442     double* zptr = arrays[2];
00443     int i;
00444     for( i = 0, rit = localGidVerts.begin(); i < nLocalVertices; i++, ++rit )
00445     {
00446         assert( *rit < xVertVals.size() + 1 );
00447         xptr[i] = xVertVals[( *rit ) - 1];  // lon
00448         yptr[i] = yVertVals[( *rit ) - 1];  // lat
00449     }
00450 
00451     // Convert lon/lat/rad to x/y/z
00452     const double pideg = acos( -1.0 ) / 180.0;
00453     double rad         = ( isConnFile ) ? 8000.0 : 8000.0 + levVals[0];
00454     for( i = 0; i < nLocalVertices; i++ )
00455     {
00456         double cosphi = cos( pideg * yptr[i] );
00457         double zmult  = sin( pideg * yptr[i] );
00458         double xmult  = cosphi * cos( xptr[i] * pideg );
00459         double ymult  = cosphi * sin( xptr[i] * pideg );
00460         xptr[i]       = rad * xmult;
00461         yptr[i]       = rad * ymult;
00462         zptr[i]       = rad * zmult;
00463     }
00464 
00465     // Get ptr to gid memory for vertices
00466     Range vert_range( start_vertex, start_vertex + nLocalVertices - 1 );
00467     void* data;
00468     int count;
00469     rval = mbImpl->tag_iterate( mGlobalIdTag, vert_range.begin(), vert_range.end(), count, data );MB_CHK_SET_ERR( rval, "Failed to iterate global id tag on local vertices" );
00470     assert( count == nLocalVertices );
00471     int* gid_data = (int*)data;
00472     std::copy( localGidVerts.begin(), localGidVerts.end(), gid_data );
00473 
00474     // Duplicate global id data, which will be used to resolve sharing
00475     if( mpFileIdTag )
00476     {
00477         rval = mbImpl->tag_iterate( *mpFileIdTag, vert_range.begin(), vert_range.end(), count, data );MB_CHK_SET_ERR( rval, "Failed to iterate file id tag on local vertices" );
00478         assert( count == nLocalVertices );
00479         int bytes_per_tag = 4;
00480         rval              = mbImpl->tag_get_bytes( *mpFileIdTag, bytes_per_tag );MB_CHK_SET_ERR( rval, "Can't get number of bytes for file id tag" );
00481         if( 4 == bytes_per_tag )
00482         {
00483             gid_data = (int*)data;
00484             std::copy( localGidVerts.begin(), localGidVerts.end(), gid_data );
00485         }
00486         else if( 8 == bytes_per_tag )
00487         {  // Should be a handle tag on 64 bit machine?
00488             long* handle_tag_data = (long*)data;
00489             std::copy( localGidVerts.begin(), localGidVerts.end(), handle_tag_data );
00490         }
00491     }
00492 
00493     // Create map from file ids to vertex handles, used later to set connectivity
00494     std::map< EntityHandle, EntityHandle > vert_handles;
00495     for( rit = localGidVerts.begin(), i = 0; rit != localGidVerts.end(); ++rit, i++ )
00496         vert_handles[*rit] = start_vertex + i;
00497 
00498     // Compute proper handles in connectivity using offset
00499     for( int q = 0; q < 4 * num_coarse_quads; q++ )
00500     {
00501         conn_arr[q] = vert_handles[conn_arr[q]];
00502         assert( conn_arr[q] );
00503     }
00504     if( spectralMesh )
00505     {
00506         int verts_per_quad = ( _spectralOrder + 1 ) * ( _spectralOrder + 1 );
00507         for( int q = 0; q < verts_per_quad * num_coarse_quads; q++ )
00508         {
00509             sv_ptr[q] = vert_handles[sv_ptr[q]];
00510             assert( sv_ptr[q] );
00511         }
00512     }
00513 
00514     // Add new vertices and quads to current file set
00515     faces.merge( tmp_range );
00516     tmp_range.insert( start_vertex, start_vertex + nLocalVertices - 1 );
00517     rval = mbImpl->add_entities( _fileSet, tmp_range );MB_CHK_SET_ERR( rval, "Failed to add new vertices and quads to current file set" );
00518 
00519     // Mark the set with the spectral order
00520     Tag sporder;
00521     rval = mbImpl->tag_get_handle( "SPECTRAL_ORDER", 1, MB_TYPE_INTEGER, sporder, MB_TAG_SPARSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Trouble creating SPECTRAL_ORDER tag" );
00522     rval = mbImpl->tag_set_data( sporder, &_fileSet, 1, &_spectralOrder );MB_CHK_SET_ERR( rval, "Trouble setting data to SPECTRAL_ORDER tag" );
00523 
00524     if( create_gathers )
00525     {
00526         EntityHandle gather_set;
00527         rval = _readNC->readMeshIface->create_gather_set( gather_set );MB_CHK_SET_ERR( rval, "Failed to create gather set" );
00528 
00529         // Create vertices
00530         arrays.clear();
00531         // Don't need to specify allocation number here, because we know enough verts were created
00532         // before
00533         rval = _readNC->readMeshIface->get_node_coords( 3, nVertices, 0, start_vertex, arrays );MB_CHK_SET_ERR( rval, "Failed to create gather set vertices" );
00534 
00535         xptr = arrays[0];
00536         yptr = arrays[1];
00537         zptr = arrays[2];
00538         for( i = 0; i < nVertices; i++ )
00539         {
00540             double cosphi = cos( pideg * yVertVals[i] );
00541             double zmult  = sin( pideg * yVertVals[i] );
00542             double xmult  = cosphi * cos( xVertVals[i] * pideg );
00543             double ymult  = cosphi * sin( xVertVals[i] * pideg );
00544             xptr[i]       = rad * xmult;
00545             yptr[i]       = rad * ymult;
00546             zptr[i]       = rad * zmult;
00547         }
00548 
00549         // Get ptr to gid memory for vertices
00550         Range gather_set_verts_range( start_vertex, start_vertex + nVertices - 1 );
00551         rval = mbImpl->tag_iterate( mGlobalIdTag, gather_set_verts_range.begin(), gather_set_verts_range.end(), count,
00552                                     data );MB_CHK_SET_ERR( rval, "Failed to iterate global id tag on gather set vertices" );
00553         assert( count == nVertices );
00554         gid_data = (int*)data;
00555         for( int j = 1; j <= nVertices; j++ )
00556             gid_data[j - 1] = j;
00557         // Set the file id tag too, it should be bigger something not interfering with global id
00558         if( mpFileIdTag )
00559         {
00560             rval = mbImpl->tag_iterate( *mpFileIdTag, gather_set_verts_range.begin(), gather_set_verts_range.end(),
00561                                         count, data );MB_CHK_SET_ERR( rval, "Failed to iterate file id tag on gather set vertices" );
00562             assert( count == nVertices );
00563             int bytes_per_tag = 4;
00564             rval              = mbImpl->tag_get_bytes( *mpFileIdTag, bytes_per_tag );MB_CHK_SET_ERR( rval, "Can't get number of bytes for file id tag" );
00565             if( 4 == bytes_per_tag )
00566             {
00567                 gid_data = (int*)data;
00568                 for( int j = 1; j <= nVertices; j++ )
00569                     gid_data[j - 1] = nVertices + j;  // Bigger than global id tag
00570             }
00571             else if( 8 == bytes_per_tag )
00572             {  // Should be a handle tag on 64 bit machine?
00573                 long* handle_tag_data = (long*)data;
00574                 for( int j = 1; j <= nVertices; j++ )
00575                     handle_tag_data[j - 1] = nVertices + j;  // Bigger than global id tag
00576             }
00577         }
00578 
00579         rval = mbImpl->add_entities( gather_set, gather_set_verts_range );MB_CHK_SET_ERR( rval, "Failed to add vertices to the gather set" );
00580 
00581         // Create quads
00582         Range gather_set_quads_range;
00583         // Don't need to specify allocation number here, because we know enough quads were created
00584         // before
00585         rval = _readNC->readMeshIface->get_element_connect( num_quads, 4, MBQUAD, 0, start_quad, conn_arr );MB_CHK_SET_ERR( rval, "Failed to create gather set quads" );
00586         gather_set_quads_range.insert( start_quad, start_quad + num_quads - 1 );
00587         int* tmp_conn_end = ( &tmp_conn[4 * num_quads - 1] ) + 1;
00588         std::copy( &tmp_conn[0], tmp_conn_end, conn_arr );
00589         for( i = 0; i != 4 * num_quads; i++ )
00590             conn_arr[i] += start_vertex - 1;  // Connectivity array is shifted by where the gather verts start
00591         rval = mbImpl->add_entities( gather_set, gather_set_quads_range );MB_CHK_SET_ERR( rval, "Failed to add quads to the gather set" );
00592     }
00593 
00594     return MB_SUCCESS;
00595 }
00596 
00597 ErrorCode NCHelperHOMME::read_ucd_variables_to_nonset_allocate( std::vector< ReadNC::VarData >& vdatas,
00598                                                                 std::vector< int >& tstep_nums )
00599 {
00600     Interface*& mbImpl          = _readNC->mbImpl;
00601     std::vector< int >& dimLens = _readNC->dimLens;
00602     DebugOutput& dbgOut         = _readNC->dbgOut;
00603 
00604     Range* range = NULL;
00605 
00606     // Get vertices
00607     Range verts;
00608     ErrorCode rval = mbImpl->get_entities_by_dimension( _fileSet, 0, verts );MB_CHK_SET_ERR( rval, "Trouble getting vertices in current file set" );
00609     assert( "Should only have a single vertex subrange, since they were read in one shot" && verts.psize() == 1 );
00610 
00611     for( unsigned int i = 0; i < vdatas.size(); i++ )
00612     {
00613         // Support non-set variables with 3 dimensions like (time, lev, ncol)
00614         assert( 3 == vdatas[i].varDims.size() );
00615 
00616         // For a non-set variable, time should be the first dimension
00617         assert( tDim == vdatas[i].varDims[0] );
00618 
00619         // Set up readStarts and readCounts
00620         vdatas[i].readStarts.resize( 3 );
00621         vdatas[i].readCounts.resize( 3 );
00622 
00623         // First: time
00624         vdatas[i].readStarts[0] = 0;  // This value is timestep dependent, will be set later
00625         vdatas[i].readCounts[0] = 1;
00626 
00627         // Next: lev
00628         vdatas[i].readStarts[1] = 0;
00629         vdatas[i].readCounts[1] = vdatas[i].numLev;
00630 
00631         // Finally: ncol
00632         switch( vdatas[i].entLoc )
00633         {
00634             case ReadNC::ENTLOCVERT:
00635                 // Vertices
00636                 // Start from the first localGidVerts
00637                 // Actually, this will be reset later on in a loop
00638                 vdatas[i].readStarts[2] = localGidVerts[0] - 1;
00639                 vdatas[i].readCounts[2] = nLocalVertices;
00640                 range                   = &verts;
00641                 break;
00642             default:
00643                 MB_SET_ERR( MB_FAILURE, "Unexpected entity location type for variable " << vdatas[i].varName );
00644         }
00645 
00646         // Get variable size
00647         vdatas[i].sz = 1;
00648         for( std::size_t idx = 0; idx != 3; idx++ )
00649             vdatas[i].sz *= vdatas[i].readCounts[idx];
00650 
00651         for( unsigned int t = 0; t < tstep_nums.size(); t++ )
00652         {
00653             dbgOut.tprintf( 2, "Reading variable %s, time step %d\n", vdatas[i].varName.c_str(), tstep_nums[t] );
00654 
00655             if( tstep_nums[t] >= dimLens[tDim] )
00656             { MB_SET_ERR( MB_INDEX_OUT_OF_RANGE, "Wrong value for timestep number " << tstep_nums[t] ); }
00657 
00658             // Get the tag to read into
00659             if( !vdatas[i].varTags[t] )
00660             {
00661                 rval = get_tag_to_nonset( vdatas[i], tstep_nums[t], vdatas[i].varTags[t], vdatas[i].numLev );MB_CHK_SET_ERR( rval, "Trouble getting tag for variable " << vdatas[i].varName );
00662             }
00663 
00664             // Get ptr to tag space
00665             void* data;
00666             int count;
00667             rval = mbImpl->tag_iterate( vdatas[i].varTags[t], range->begin(), range->end(), count, data );MB_CHK_SET_ERR( rval, "Failed to iterate tag for variable " << vdatas[i].varName );
00668             assert( (unsigned)count == range->size() );
00669             vdatas[i].varDatas[t] = data;
00670         }
00671     }
00672 
00673     return rval;
00674 }
00675 
00676 #ifdef MOAB_HAVE_PNETCDF
00677 ErrorCode NCHelperHOMME::read_ucd_variables_to_nonset_async( std::vector< ReadNC::VarData >& vdatas,
00678                                                              std::vector< int >& tstep_nums )
00679 {
00680     DebugOutput& dbgOut = _readNC->dbgOut;
00681 
00682     ErrorCode rval = read_ucd_variables_to_nonset_allocate( vdatas, tstep_nums );MB_CHK_SET_ERR( rval, "Trouble allocating space to read non-set variables" );
00683 
00684     // Finally, read into that space
00685     int success;
00686 
00687     for( unsigned int i = 0; i < vdatas.size(); i++ )
00688     {
00689         std::size_t sz = vdatas[i].sz;
00690 
00691         // A typical supported variable: float T(time, lev, ncol)
00692         // For tag values, need transpose (lev, ncol) to (ncol, lev)
00693         size_t ni = vdatas[i].readCounts[2];  // ncol
00694         size_t nj = 1;                        // Here we should just set nj to 1
00695         size_t nk = vdatas[i].readCounts[1];  // lev
00696 
00697         for( unsigned int t = 0; t < tstep_nums.size(); t++ )
00698         {
00699             // We will synchronize all these reads with the other processors,
00700             // so the wait will be inside this double loop; is it too much?
00701             size_t nb_reads = localGidVerts.psize();
00702             std::vector< int > requests( nb_reads ), statuss( nb_reads );
00703             size_t idxReq = 0;
00704 
00705             // Tag data for this timestep
00706             void* data = vdatas[i].varDatas[t];
00707 
00708             // Set readStart for each timestep along time dimension
00709             vdatas[i].readStarts[0] = tstep_nums[t];
00710 
00711             switch( vdatas[i].varDataType )
00712             {
00713                 case NC_FLOAT:
00714                 case NC_DOUBLE: {
00715                     // Read float as double
00716                     std::vector< double > tmpdoubledata( sz );
00717 
00718                     // In the case of ucd mesh, and on multiple proc,
00719                     // we need to read as many times as subranges we have in the
00720                     // localGidVerts range;
00721                     // basically, we have to give a different point
00722                     // for data to start, for every subrange :(
00723                     size_t indexInDoubleArray = 0;
00724                     size_t ic                 = 0;
00725                     for( Range::pair_iterator pair_iter = localGidVerts.pair_begin();
00726                          pair_iter != localGidVerts.pair_end(); ++pair_iter, ic++ )
00727                     {
00728                         EntityHandle starth     = pair_iter->first;
00729                         EntityHandle endh       = pair_iter->second;  // Inclusive
00730                         vdatas[i].readStarts[2] = ( NCDF_SIZE )( starth - 1 );
00731                         vdatas[i].readCounts[2] = ( NCDF_SIZE )( endh - starth + 1 );
00732 
00733                         // Do a partial read, in each subrange
00734                         // Wait outside this loop
00735                         success =
00736                             NCFUNCREQG( _vara_double )( _fileId, vdatas[i].varId, &( vdatas[i].readStarts[0] ),
00737                                                         &( vdatas[i].readCounts[0] ),
00738                                                         &( tmpdoubledata[indexInDoubleArray] ), &requests[idxReq++] );
00739                         if( success )
00740                             MB_SET_ERR( MB_FAILURE,
00741                                         "Failed to read double data in a loop for variable " << vdatas[i].varName );
00742                         // We need to increment the index in double array for the
00743                         // next subrange
00744                         indexInDoubleArray += ( endh - starth + 1 ) * 1 * vdatas[i].numLev;
00745                     }
00746                     assert( ic == localGidVerts.psize() );
00747 
00748                     success = ncmpi_wait_all( _fileId, requests.size(), &requests[0], &statuss[0] );
00749                     if( success ) MB_SET_ERR( MB_FAILURE, "Failed on wait_all" );
00750 
00751                     if( vdatas[i].numLev > 1 )
00752                         // Transpose (lev, ncol) to (ncol, lev)
00753                         kji_to_jik_stride( ni, nj, nk, data, &tmpdoubledata[0], localGidVerts );
00754                     else
00755                     {
00756                         for( std::size_t idx = 0; idx != tmpdoubledata.size(); idx++ )
00757                             ( (double*)data )[idx] = tmpdoubledata[idx];
00758                     }
00759 
00760                     break;
00761                 }
00762                 default:
00763                     MB_SET_ERR( MB_FAILURE, "Unexpected data type for variable " << vdatas[i].varName );
00764             }
00765         }
00766     }
00767 
00768     // Debug output, if requested
00769     if( 1 == dbgOut.get_verbosity() )
00770     {
00771         dbgOut.printf( 1, "Read variables: %s", vdatas.begin()->varName.c_str() );
00772         for( unsigned int i = 1; i < vdatas.size(); i++ )
00773             dbgOut.printf( 1, ", %s ", vdatas[i].varName.c_str() );
00774         dbgOut.tprintf( 1, "\n" );
00775     }
00776 
00777     return rval;
00778 }
00779 #else
00780 ErrorCode NCHelperHOMME::read_ucd_variables_to_nonset( std::vector< ReadNC::VarData >& vdatas,
00781                                                        std::vector< int >& tstep_nums )
00782 {
00783     DebugOutput& dbgOut = _readNC->dbgOut;
00784 
00785     ErrorCode rval = read_ucd_variables_to_nonset_allocate( vdatas, tstep_nums );MB_CHK_SET_ERR( rval, "Trouble allocating space to read non-set variables" );
00786 
00787     // Finally, read into that space
00788     int success;
00789     for( unsigned int i = 0; i < vdatas.size(); i++ )
00790     {
00791         std::size_t sz = vdatas[i].sz;
00792 
00793         // A typical supported variable: float T(time, lev, ncol)
00794         // For tag values, need transpose (lev, ncol) to (ncol, lev)
00795         size_t ni = vdatas[i].readCounts[2];  // ncol
00796         size_t nj = 1;                        // Here we should just set nj to 1
00797         size_t nk = vdatas[i].readCounts[1];  // lev
00798 
00799         for( unsigned int t = 0; t < tstep_nums.size(); t++ )
00800         {
00801             // Tag data for this timestep
00802             void* data = vdatas[i].varDatas[t];
00803 
00804             // Set readStart for each timestep along time dimension
00805             vdatas[i].readStarts[0] = tstep_nums[t];
00806 
00807             switch( vdatas[i].varDataType )
00808             {
00809                 case NC_FLOAT:
00810                 case NC_DOUBLE: {
00811                     // Read float as double
00812                     std::vector< double > tmpdoubledata( sz );
00813 
00814                     // In the case of ucd mesh, and on multiple proc,
00815                     // we need to read as many times as subranges we have in the
00816                     // localGidVerts range;
00817                     // basically, we have to give a different point
00818                     // for data to start, for every subrange :(
00819                     size_t indexInDoubleArray = 0;
00820                     size_t ic = 0;
00821                     for( Range::pair_iterator pair_iter = localGidVerts.pair_begin();
00822                          pair_iter != localGidVerts.pair_end(); ++pair_iter, ic++ )
00823                     {
00824                         EntityHandle starth = pair_iter->first;
00825                         EntityHandle endh = pair_iter->second;  // Inclusive
00826                         vdatas[i].readStarts[2] = ( NCDF_SIZE )( starth - 1 );
00827                         vdatas[i].readCounts[2] = ( NCDF_SIZE )( endh - starth + 1 );
00828 
00829                         success = NCFUNCAG( _vara_double )( _fileId, vdatas[i].varId, &( vdatas[i].readStarts[0] ),
00830                                                             &( vdatas[i].readCounts[0] ),
00831                                                             &( tmpdoubledata[indexInDoubleArray] ) );
00832                         if( success )
00833                             MB_SET_ERR( MB_FAILURE,
00834                                         "Failed to read double data in a loop for variable " << vdatas[i].varName );
00835                         // We need to increment the index in double array for the
00836                         // next subrange
00837                         indexInDoubleArray += ( endh - starth + 1 ) * 1 * vdatas[i].numLev;
00838                     }
00839                     assert( ic == localGidVerts.psize() );
00840 
00841                     if( vdatas[i].numLev > 1 )
00842                         // Transpose (lev, ncol) to (ncol, lev)
00843                         kji_to_jik_stride( ni, nj, nk, data, &tmpdoubledata[0], localGidVerts );
00844                     else
00845                     {
00846                         for( std::size_t idx = 0; idx != tmpdoubledata.size(); idx++ )
00847                             ( (double*)data )[idx] = tmpdoubledata[idx];
00848                     }
00849 
00850                     break;
00851                 }
00852                 default:
00853                     MB_SET_ERR( MB_FAILURE, "Unexpected data type for variable " << vdatas[i].varName );
00854             }
00855         }
00856     }
00857 
00858     // Debug output, if requested
00859     if( 1 == dbgOut.get_verbosity() )
00860     {
00861         dbgOut.printf( 1, "Read variables: %s", vdatas.begin()->varName.c_str() );
00862         for( unsigned int i = 1; i < vdatas.size(); i++ )
00863             dbgOut.printf( 1, ", %s ", vdatas[i].varName.c_str() );
00864         dbgOut.tprintf( 1, "\n" );
00865     }
00866 
00867     return rval;
00868 }
00869 #endif
00870 
00871 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines