Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
NCHelperDomain.cpp
Go to the documentation of this file.
00001 #include "NCHelperDomain.hpp"
00002 #include "moab/FileOptions.hpp"
00003 #include "moab/ReadUtilIface.hpp"
00004 #include "moab/IntxMesh/IntxUtils.hpp"
00005 #include "AEntityFactory.hpp"
00006 #ifdef MOAB_HAVE_MPI
00007 #include "moab/ParallelMergeMesh.hpp"
00008 #endif
00009 #include <cmath>
00010 #include <sstream>
00011 
00012 namespace moab
00013 {
00014 
00015 bool NCHelperDomain::can_read_file( ReadNC* readNC, int fileId )
00016 {
00017     std::vector< std::string >& dimNames = readNC->dimNames;
00018 
00019     // If dimension names "n" AND "ni" AND "nj" AND "nv" exist then it should be the Domain grid
00020     if( ( std::find( dimNames.begin(), dimNames.end(), std::string( "n" ) ) != dimNames.end() ) &&
00021         ( std::find( dimNames.begin(), dimNames.end(), std::string( "ni" ) ) != dimNames.end() ) &&
00022         ( std::find( dimNames.begin(), dimNames.end(), std::string( "nj" ) ) != dimNames.end() ) &&
00023         ( std::find( dimNames.begin(), dimNames.end(), std::string( "nv" ) ) != dimNames.end() ) )
00024     {
00025         // Make sure it is CAM grid
00026         std::map< std::string, ReadNC::AttData >::iterator attIt = readNC->globalAtts.find( "source" );
00027         if( attIt == readNC->globalAtts.end() ) return false;
00028         unsigned int sz = attIt->second.attLen;
00029         std::string att_data;
00030         att_data.resize( sz + 1 );
00031         att_data[sz] = '\000';
00032         int success =
00033             NCFUNC( get_att_text )( fileId, attIt->second.attVarId, attIt->second.attName.c_str(), &att_data[0] );
00034         if( success ) return false;
00035         /*if (att_data.find("CAM") == std::string::npos)
00036           return false;*/
00037 
00038         return true;
00039     }
00040 
00041     return false;
00042 }
00043 
00044 ErrorCode NCHelperDomain::init_mesh_vals()
00045 {
00046     Interface*& mbImpl                                = _readNC->mbImpl;
00047     std::vector< std::string >& dimNames              = _readNC->dimNames;
00048     std::vector< int >& dimLens                       = _readNC->dimLens;
00049     std::map< std::string, ReadNC::VarData >& varInfo = _readNC->varInfo;
00050     DebugOutput& dbgOut                               = _readNC->dbgOut;
00051 #ifdef MOAB_HAVE_MPI
00052     bool& isParallel = _readNC->isParallel;
00053 #endif
00054     int& partMethod     = _readNC->partMethod;
00055     ScdParData& parData = _readNC->parData;
00056 
00057     ErrorCode rval;
00058 
00059     // Look for names of i/j dimensions
00060     // First i
00061     std::vector< std::string >::iterator vit;
00062     unsigned int idx;
00063     if( ( vit = std::find( dimNames.begin(), dimNames.end(), "ni" ) ) != dimNames.end() )
00064         idx = vit - dimNames.begin();
00065     else
00066     {
00067         MB_SET_ERR( MB_FAILURE, "Couldn't find 'ni' variable" );
00068     }
00069     iDim     = idx;
00070     gDims[0] = 0;
00071     gDims[3] = dimLens[idx];
00072 
00073     // Then j
00074     if( ( vit = std::find( dimNames.begin(), dimNames.end(), "nj" ) ) != dimNames.end() )
00075         idx = vit - dimNames.begin();
00076     else
00077     {
00078         MB_SET_ERR( MB_FAILURE, "Couldn't find 'nj' variable" );
00079     }
00080     jDim     = idx;
00081     gDims[1] = 0;
00082     gDims[4] = dimLens[idx];  // Add 2 for the pole points ? not needed
00083 
00084     // do not use gcdims ? or use only gcdims?
00085 
00086     // Try a truly 2D mesh
00087     gDims[2] = -1;
00088     gDims[5] = -1;
00089 
00090     // Get number of vertices per cell
00091     if( ( vit = std::find( dimNames.begin(), dimNames.end(), "nv" ) ) != dimNames.end() )
00092         idx = vit - dimNames.begin();
00093     else
00094     {
00095         MB_SET_ERR( MB_FAILURE, "Couldn't find 'nv' dimension" );
00096     }
00097     nvDim = idx;
00098     nv    = dimLens[idx];
00099 
00100     // Parse options to get subset
00101     int rank = 0, procs = 1;
00102 #ifdef MOAB_HAVE_MPI
00103     if( isParallel )
00104     {
00105         ParallelComm*& myPcomm = _readNC->myPcomm;
00106         rank                   = myPcomm->proc_config().proc_rank();
00107         procs                  = myPcomm->proc_config().proc_size();
00108     }
00109 #endif
00110     if( procs > 1 )
00111     {
00112         for( int i = 0; i < 6; i++ )
00113             parData.gDims[i] = gDims[i];
00114         parData.partMethod = partMethod;
00115         int pdims[3];
00116 
00117         locallyPeriodic[0] = locallyPeriodic[1] = locallyPeriodic[2] = 0;
00118         rval = ScdInterface::compute_partition( procs, rank, parData, lDims, locallyPeriodic, pdims );MB_CHK_ERR( rval );
00119         for( int i = 0; i < 3; i++ )
00120             parData.pDims[i] = pdims[i];
00121 
00122         dbgOut.tprintf( 1, "Partition: %dx%d (out of %dx%d)\n", lDims[3] - lDims[0], lDims[4] - lDims[1],
00123                         gDims[3] - gDims[0], gDims[4] - gDims[1] );
00124         if( 0 == rank )
00125             dbgOut.tprintf( 1, "Contiguous chunks of size %d bytes.\n",
00126                             8 * ( lDims[3] - lDims[0] ) * ( lDims[4] - lDims[1] ) );
00127     }
00128     else
00129     {
00130         for( int i = 0; i < 6; i++ )
00131             lDims[i] = gDims[i];
00132         locallyPeriodic[0] = globallyPeriodic[0];
00133     }
00134 
00135     // Now get actual coordinate values for vertices and cell centers
00136     lCDims[0] = lDims[0];
00137 
00138     lCDims[3] = lDims[3];
00139 
00140     // For FV models, will always be non-periodic in j
00141     lCDims[1] = lDims[1];
00142     lCDims[4] = lDims[4];
00143 
00144 #if 0
00145   // Resize vectors to store values later
00146   if (-1 != lDims[0])
00147     ilVals.resize(lDims[3] - lDims[0] + 1);
00148   if (-1 != lCDims[0])
00149     ilCVals.resize(lCDims[3] - lCDims[0] + 1);
00150   if (-1 != lDims[1])
00151     jlVals.resize(lDims[4] - lDims[1] + 1);
00152   if (-1 != lCDims[1])
00153     jlCVals.resize(lCDims[4] - lCDims[1] + 1);
00154   if (nTimeSteps > 0)
00155     tVals.resize(nTimeSteps);
00156 #endif
00157 
00158     dbgOut.tprintf( 1, "I=%d-%d, J=%d-%d\n", lDims[0], lDims[3], lDims[1], lDims[4] );
00159     dbgOut.tprintf( 1, "%d elements, %d vertices\n", ( lDims[3] - lDims[0] ) * ( lDims[4] - lDims[1] ),
00160                     ( lDims[3] - lDims[0] ) * ( lDims[4] - lDims[1] ) * nv );
00161 
00162     // For each variable, determine the entity location type and number of levels
00163     std::map< std::string, ReadNC::VarData >::iterator mit;
00164     for( mit = varInfo.begin(); mit != varInfo.end(); ++mit )
00165     {
00166         ReadNC::VarData& vd = ( *mit ).second;
00167 
00168         // Default entLoc is ENTLOCSET
00169         if( std::find( vd.varDims.begin(), vd.varDims.end(), tDim ) != vd.varDims.end() )
00170         {
00171             if( ( std::find( vd.varDims.begin(), vd.varDims.end(), iCDim ) != vd.varDims.end() ) &&
00172                 ( std::find( vd.varDims.begin(), vd.varDims.end(), jCDim ) != vd.varDims.end() ) )
00173                 vd.entLoc = ReadNC::ENTLOCFACE;
00174             else if( ( std::find( vd.varDims.begin(), vd.varDims.end(), jDim ) != vd.varDims.end() ) &&
00175                      ( std::find( vd.varDims.begin(), vd.varDims.end(), iCDim ) != vd.varDims.end() ) )
00176                 vd.entLoc = ReadNC::ENTLOCNSEDGE;
00177             else if( ( std::find( vd.varDims.begin(), vd.varDims.end(), jCDim ) != vd.varDims.end() ) &&
00178                      ( std::find( vd.varDims.begin(), vd.varDims.end(), iDim ) != vd.varDims.end() ) )
00179                 vd.entLoc = ReadNC::ENTLOCEWEDGE;
00180         }
00181 
00182         // Default numLev is 0
00183         if( std::find( vd.varDims.begin(), vd.varDims.end(), levDim ) != vd.varDims.end() ) vd.numLev = nLevels;
00184     }
00185 
00186     std::vector< std::string > ijdimNames( 2 );
00187     ijdimNames[0] = "__ni";
00188     ijdimNames[1] = "__nj";
00189     std::string tag_name;
00190     Tag tagh;
00191 
00192     // __<dim_name>_LOC_MINMAX (for slon, slat, lon and lat)
00193     for( unsigned int i = 0; i != ijdimNames.size(); i++ )
00194     {
00195         std::vector< int > val( 2, 0 );
00196         if( ijdimNames[i] == "__ni" )
00197         {
00198             val[0] = lDims[0];
00199             val[1] = lDims[3];
00200         }
00201         else if( ijdimNames[i] == "__nj" )
00202         {
00203             val[0] = lDims[1];
00204             val[1] = lDims[4];
00205         }
00206 
00207         std::stringstream ss_tag_name;
00208         ss_tag_name << ijdimNames[i] << "_LOC_MINMAX";
00209         tag_name = ss_tag_name.str();
00210         rval     = mbImpl->tag_get_handle( tag_name.c_str(), 2, MB_TYPE_INTEGER, tagh, MB_TAG_SPARSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Trouble creating conventional tag " << tag_name );
00211         rval = mbImpl->tag_set_data( tagh, &_fileSet, 1, &val[0] );MB_CHK_SET_ERR( rval, "Trouble setting data to conventional tag " << tag_name );
00212         if( MB_SUCCESS == rval ) dbgOut.tprintf( 2, "Conventional tag %s is created.\n", tag_name.c_str() );
00213     }
00214 
00215     // __<dim_name>_LOC_VALS (for slon, slat, lon and lat)
00216     // Assume all have the same data type as lon (expected type is float or double)
00217     switch( varInfo["xc"].varDataType )
00218     {
00219         case NC_FLOAT:
00220         case NC_DOUBLE:
00221             break;
00222         default:
00223             MB_SET_ERR( MB_FAILURE, "Unexpected variable data type for 'lon'" );
00224     }
00225 
00226     // do not need conventional tags
00227     Tag convTagsCreated = 0;
00228     int def_val         = 0;
00229     rval                = mbImpl->tag_get_handle( "__CONV_TAGS_CREATED", 1, MB_TYPE_INTEGER, convTagsCreated,
00230                                                   MB_TAG_SPARSE | MB_TAG_CREAT, &def_val );MB_CHK_SET_ERR( rval, "Trouble getting _CONV_TAGS_CREATED tag" );
00231     int create_conv_tags_flag = 1;
00232     rval                      = mbImpl->tag_set_data( convTagsCreated, &_fileSet, 1, &create_conv_tags_flag );MB_CHK_SET_ERR( rval, "Trouble setting _CONV_TAGS_CREATED tag" );
00233 
00234     return MB_SUCCESS;
00235 }
00236 
00237 ErrorCode NCHelperDomain::create_mesh( Range& faces )
00238 {
00239     Interface*& mbImpl = _readNC->mbImpl;
00240     // std::string& fileName = _readNC->fileName;
00241     Tag& mGlobalIdTag = _readNC->mGlobalIdTag;
00242     // const Tag*& mpFileIdTag = _readNC->mpFileIdTag;
00243     DebugOutput& dbgOut = _readNC->dbgOut;
00244     /*int& gatherSetRank = _readNC->gatherSetRank;
00245     int& trivialPartitionShift = _readNC->trivialPartitionShift;*/
00246     /*
00247 
00248       int rank = 0;
00249       int procs = 1;
00250     #ifdef MOAB_HAVE_MPI
00251       bool& isParallel = _readNC->isParallel;
00252       if (isParallel) {
00253         ParallelComm*& myPcomm = _readNC->myPcomm;
00254         rank = myPcomm->proc_config().proc_rank();
00255         procs = myPcomm->proc_config().proc_size();
00256       }
00257     #endif
00258     */
00259 
00260     ErrorCode rval;
00261     int success = 0;
00262 
00263     /*
00264       bool create_gathers = false;
00265       if (rank == gatherSetRank)
00266         create_gathers = true;
00267 
00268       // Shift rank to obtain a rotated trivial partition
00269       int shifted_rank = rank;
00270       if (procs >= 2 && trivialPartitionShift > 0)
00271         shifted_rank = (rank + trivialPartitionShift) % procs;*/
00272 
00273     // how many will have mask 0 or 1
00274     // how many will have a fraction  ? we will not instantiate all elements; only those with mask 1
00275     // ? also, not all vertices, only those that belong to mask 1 elements ? we will not care about
00276     // duplicate vertices; maybe another time ? we will start reading masks, vertices
00277     int local_elems = ( lDims[4] - lDims[1] ) * ( lDims[3] - lDims[0] );
00278     dbgOut.tprintf( 1, "local cells: %d \n", local_elems );
00279 
00280     // count how many will be with mask 1 here
00281     // basically, read the mask variable on the local elements;
00282     std::string maskstr( "mask" );
00283     ReadNC::VarData& vmask = _readNC->varInfo[maskstr];
00284 
00285     // mask is (nj, ni)
00286     vmask.readStarts.push_back( lDims[1] );
00287     vmask.readStarts.push_back( lDims[0] );
00288     vmask.readCounts.push_back( lDims[4] - lDims[1] );
00289     vmask.readCounts.push_back( lDims[3] - lDims[0] );
00290     std::vector< int > mask( local_elems );
00291     success = NCFUNCAG( _vara_int )( _fileId, vmask.varId, &vmask.readStarts[0], &vmask.readCounts[0], &mask[0] );
00292     if( success ) MB_SET_ERR( MB_FAILURE, "Failed to read int data for mask variable " );
00293 
00294     int nb_with_mask1 = 0;
00295     for( int i = 0; i < local_elems; i++ )
00296         if( 1 == mask[i] ) nb_with_mask1++;
00297 
00298     dbgOut.tprintf( 1, "local cells with mask 1: %d \n", nb_with_mask1 );
00299     std::vector< NCDF_SIZE > startsv( 3 );
00300     startsv[0] = vmask.readStarts[0];
00301     startsv[1] = vmask.readStarts[1];
00302     startsv[2] = 0;
00303     std::vector< NCDF_SIZE > countsv( 3 );
00304     countsv[0] = vmask.readCounts[0];
00305     countsv[1] = vmask.readCounts[1];
00306     countsv[2] = nv;  // number of vertices per element
00307 
00308     // read xv and yv coords for vertices, and create elements;
00309     std::string xvstr( "xv" );
00310     ReadNC::VarData& var_xv = _readNC->varInfo[xvstr];
00311     std::vector< double > xv( local_elems * nv );
00312     success = NCFUNCAG( _vara_double )( _fileId, var_xv.varId, &startsv[0], &countsv[0], &xv[0] );
00313     if( success ) MB_SET_ERR( MB_FAILURE, "Failed to read double data for xv variable " );
00314 
00315     std::string yvstr( "yv" );
00316     ReadNC::VarData& var_yv = _readNC->varInfo[yvstr];
00317     std::vector< double > yv( local_elems * nv );
00318     success = NCFUNCAG( _vara_double )( _fileId, var_yv.varId, &startsv[0], &countsv[0], &yv[0] );
00319     if( success ) MB_SET_ERR( MB_FAILURE, "Failed to read double data for yv variable " );
00320 
00321     // read other variables, like xc, yc, frac, area
00322     std::string xcstr( "xc" );
00323     ReadNC::VarData& var_xc = _readNC->varInfo[xcstr];
00324     std::vector< double > xc( local_elems );
00325     success = NCFUNCAG( _vara_double )( _fileId, var_xc.varId, &vmask.readStarts[0], &vmask.readCounts[0], &xc[0] );
00326     if( success ) MB_SET_ERR( MB_FAILURE, "Failed to read double data for xc variable " );
00327 
00328     std::string ycstr( "yc" );
00329     ReadNC::VarData& var_yc = _readNC->varInfo[ycstr];
00330     std::vector< double > yc( local_elems );
00331     success = NCFUNCAG( _vara_double )( _fileId, var_yc.varId, &vmask.readStarts[0], &vmask.readCounts[0], &yc[0] );
00332     if( success ) MB_SET_ERR( MB_FAILURE, "Failed to read double data for yc variable " );
00333 
00334     std::string fracstr( "frac" );
00335     ReadNC::VarData& var_frac = _readNC->varInfo[fracstr];
00336     std::vector< double > frac( local_elems );
00337     success = NCFUNCAG( _vara_double )( _fileId, var_frac.varId, &vmask.readStarts[0], &vmask.readCounts[0], &frac[0] );
00338     if( success ) MB_SET_ERR( MB_FAILURE, "Failed to read double data for frac variable " );
00339     std::string areastr( "area" );
00340     ReadNC::VarData& var_area = _readNC->varInfo[areastr];
00341     std::vector< double > area( local_elems );
00342     success = NCFUNCAG( _vara_double )( _fileId, var_area.varId, &vmask.readStarts[0], &vmask.readCounts[0], &area[0] );
00343     if( success ) MB_SET_ERR( MB_FAILURE, "Failed to read double data for area variable " );
00344     // create tags for them
00345     Tag areaTag, fracTag, xcTag, ycTag;
00346     rval = mbImpl->tag_get_handle( "area", 1, MB_TYPE_DOUBLE, areaTag, MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Trouble creating area tag" );
00347     rval = mbImpl->tag_get_handle( "frac", 1, MB_TYPE_DOUBLE, fracTag, MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Trouble creating frac tag" );
00348     rval = mbImpl->tag_get_handle( "xc", 1, MB_TYPE_DOUBLE, xcTag, MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Trouble creating xc tag" );
00349     rval = mbImpl->tag_get_handle( "yc", 1, MB_TYPE_DOUBLE, ycTag, MB_TAG_DENSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Trouble creating yc tag" );
00350 
00351     //
00352     EntityHandle* conn_arr;
00353     EntityHandle vtx_handle;
00354     Range tmp_range;
00355 
00356     // set connectivity into that space
00357 
00358     EntityHandle start_cell;
00359     EntityType mdb_type = MBVERTEX;
00360     if( nv == 3 )
00361         mdb_type = MBTRI;
00362     else if( nv == 4 )
00363         mdb_type = MBQUAD;
00364     else if( nv > 4 )  // (nv > 4)
00365         mdb_type = MBPOLYGON;
00366     // for nv = 1 , type is vertex
00367 
00368     if( nv > 1 && nb_with_mask1 > 0 )
00369     {
00370         rval = _readNC->readMeshIface->get_element_connect( nb_with_mask1, nv, mdb_type, 0, start_cell, conn_arr );MB_CHK_SET_ERR( rval, "Failed to create local cells" );
00371         tmp_range.insert( start_cell, start_cell + nb_with_mask1 - 1 );
00372     }
00373 
00374     // Create vertices; first identify different ones, with a tolerance
00375     std::map< Node3D, EntityHandle > vertex_map;
00376 
00377     if (nb_with_mask1 > 0) {
00378         // Set vertex coordinates
00379         // will read all xv, yv, but use only those with correct mask on
00380 
00381         int elem_index     = 0;  // total index in netcdf arrays
00382         const double pideg = acos( -1.0 ) / 180.0;
00383 
00384         for( ; elem_index < local_elems; elem_index++ )
00385         {
00386             if( 0 == mask[elem_index] ) continue;  // nothing to do, do not advance elem_index in actual moab arrays
00387             // set area and fraction on those elements too
00388             for( int k = 0; k < nv; k++ )
00389             {
00390                 int index_v_arr = nv * elem_index + k;
00391                 double x, y;
00392                 if( nv > 1 )
00393                 {
00394                     x             = xv[index_v_arr];
00395                     y             = yv[index_v_arr];
00396                     double cosphi = cos( pideg * y );
00397                     double zmult  = sin( pideg * y );
00398                     double xmult  = cosphi * cos( x * pideg );
00399                     double ymult  = cosphi * sin( x * pideg );
00400                     Node3D pt( xmult, ymult, zmult );
00401                     vertex_map[pt] = 0;
00402                 }
00403                 else
00404                 {
00405                     x = xc[elem_index];
00406                     y = yc[elem_index];
00407                     Node3D pt( x, y, 0 );
00408                     vertex_map[pt] = 0;
00409                 }
00410             }
00411         }
00412         int nLocalVertices = (int)vertex_map.size();
00413         std::vector< double* > arrays;
00414         EntityHandle start_vertex;
00415         rval = _readNC->readMeshIface->get_node_coords( 3, nLocalVertices, 0, start_vertex, arrays );MB_CHK_SET_ERR( rval, "Failed to create local vertices" );
00416 
00417         vtx_handle = start_vertex;
00418         // Copy vertex coordinates into entity sequence coordinate arrays
00419         // and copy handle into vertex_map.
00420         double *x = arrays[0], *y = arrays[1], *z = arrays[2];
00421         for( auto i = vertex_map.begin(); i != vertex_map.end(); ++i )
00422         {
00423             i->second = vtx_handle;
00424             ++vtx_handle;
00425             *x = i->first.coords[0];
00426             ++x;
00427             *y = i->first.coords[1];
00428             ++y;
00429             *z = i->first.coords[2];
00430             ++z;
00431         }
00432 
00433         // int nj = gDims[4]-gDims[1]; // is it about 1 in irregular cases
00434 
00435         int local_row_size = lCDims[3] - lCDims[0];
00436         int global_row_size = gDims[3] - gDims[0]; // this is along
00437         elem_index         = -1;
00438         int index          = 0;  // consider the mask for advancing in moab arrays;
00439         //printf(" map size :%ld \n", vertex_map.size());
00440         // create now vertex arrays, size vertex_map.size()
00441         for (int j = lCDims[1]; j<lCDims[4]; j++)
00442             for (int i = lCDims[0]; i<lCDims[3]; i++)
00443             {
00444                 elem_index++;
00445                 if( 0 == mask[elem_index] ) continue;  // nothing to do, do not advance elem_index in actual moab arrays
00446                 // set area and fraction on those elements too
00447                 for( int k = 0; k < nv; k++ )
00448                 {
00449                     int index_v_arr = nv * elem_index + k;
00450                     if( nv > 1 )
00451                     {
00452                         double x      = xv[index_v_arr];
00453                         double y      = yv[index_v_arr];
00454                         double cosphi = cos( pideg * y );
00455                         double zmult  = sin( pideg * y );
00456                         double xmult  = cosphi * cos( x * pideg );
00457                         double ymult  = cosphi * sin( x * pideg );
00458                         Node3D pt( xmult, ymult, zmult );
00459                         conn_arr[index * nv + k] = vertex_map[pt];
00460                     }
00461                 }
00462                 EntityHandle cell = start_vertex + index;
00463                 if( nv > 1 ) cell = start_cell + index;
00464                 // set other tags, like xc, yc, frac, area
00465                 rval = mbImpl->tag_set_data( xcTag, &cell, 1, &xc[elem_index] );MB_CHK_SET_ERR( rval, "Failed to set xc tag" );
00466                 rval = mbImpl->tag_set_data( ycTag, &cell, 1, &yc[elem_index] );MB_CHK_SET_ERR( rval, "Failed to set yc tag" );
00467                 rval = mbImpl->tag_set_data( areaTag, &cell, 1, &area[elem_index] );MB_CHK_SET_ERR( rval, "Failed to set area tag" );
00468                 rval = mbImpl->tag_set_data( fracTag, &cell, 1, &frac[elem_index] );MB_CHK_SET_ERR( rval, "Failed to set frac tag" );
00469 
00470                 // set the global id too:
00471                 int globalId = j * global_row_size + i + 1;
00472                 rval = mbImpl->tag_set_data( mGlobalIdTag, &cell, 1, &globalId );MB_CHK_SET_ERR( rval, "Failed to set global id tag" );
00473                 index++;
00474             }
00475 
00476 
00477         rval = mbImpl->add_entities( _fileSet, tmp_range );MB_CHK_SET_ERR( rval, "Failed to add new cells to current file set" );
00478 
00479         // modify local file set, to merge coincident vertices, and to correct repeated vertices in elements
00480         std::vector< Tag > tagList;
00481         tagList.push_back( mGlobalIdTag );
00482         tagList.push_back( xcTag );
00483         tagList.push_back( ycTag );
00484         tagList.push_back( areaTag );
00485         tagList.push_back( fracTag );
00486         rval = IntxUtils::remove_padded_vertices( mbImpl, _fileSet, tagList );MB_CHK_SET_ERR( rval, "Failed to remove duplicate vertices" );
00487 
00488         rval = mbImpl->get_entities_by_dimension( _fileSet, 2, faces );MB_CHK_ERR( rval );
00489         Range all_verts;
00490         rval = mbImpl->get_connectivity( faces, all_verts );MB_CHK_ERR( rval );
00491         //printf(" range vert size :%ld \n", all_verts.size());
00492         rval = mbImpl->add_entities( _fileSet, all_verts );MB_CHK_ERR( rval );
00493 
00494         // need to add adjacencies; TODO: fix this for all nc readers
00495         // copy this logic from migrate mesh in par comm graph
00496         Core* mb                 = (Core*)mbImpl;
00497         AEntityFactory* adj_fact = mb->a_entity_factory();
00498         if( !adj_fact->vert_elem_adjacencies() )
00499             adj_fact->create_vert_elem_adjacencies();
00500         else
00501         {
00502             for( Range::iterator it = faces.begin(); it != faces.end(); ++it )
00503             {
00504                 EntityHandle eh          = *it;
00505                 const EntityHandle* conn = NULL;
00506                 int num_nodes            = 0;
00507                 rval                     = mb->get_connectivity( eh, conn, num_nodes );MB_CHK_ERR( rval );
00508                 adj_fact->notify_create_entity( eh, conn, num_nodes );
00509             }
00510         }
00511     }
00512 
00513 
00514 
00515 #ifdef MOAB_HAVE_MPI
00516     ParallelComm*& myPcomm = _readNC->myPcomm;
00517     if( myPcomm )
00518     {
00519         double tol = 1.e-12;  // this is the same as static tolerance in NCHelper
00520         ParallelMergeMesh pmm( myPcomm, tol );
00521         rval = pmm.merge( _fileSet,
00522                           /* do not do local merge*/ false,
00523                           /*  2d cells*/ 2 );MB_CHK_SET_ERR( rval, "Failed to merge vertices in parallel" );
00524 
00525         // assign global ids only for vertices, cells have them fine
00526         rval = myPcomm->assign_global_ids( _fileSet, /*dim*/ 0 );MB_CHK_ERR( rval );
00527     }
00528 #endif
00529 
00530     return MB_SUCCESS;
00531 }
00532 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines