MOAB: Mesh Oriented datABase  (version 5.3.0)
ReadNCDF.cpp
Go to the documentation of this file.
00001 /**
00002  * MOAB, a Mesh-Oriented datABase, is a software component for creating,
00003  * storing and accessing finite element mesh data.
00004  *
00005  * Copyright 2004 Sandia Corporation.  Under the terms of Contract
00006  * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
00007  * retains certain rights in this software.
00008  *
00009  * This library is free software; you can redistribute it and/or
00010  * modify it under the terms of the GNU Lesser General Public
00011  * License as published by the Free Software Foundation; either
00012  * version 2.1 of the License, or (at your option) any later version.
00013  *
00014  */
00015 
00016 #ifdef WIN32
00017 #pragma warning( disable : 4786 )
00018 #endif
00019 
00020 #include "ReadNCDF.hpp"
00021 #include "netcdf.h"
00022 
00023 #include <algorithm>
00024 #include <string>
00025 #include <cassert>
00026 #include <cstdio>
00027 #include <cstring>
00028 #include <cmath>
00029 #include <sstream>
00030 #include <iostream>
00031 #include <map>
00032 
00033 #include "moab/CN.hpp"
00034 #include "moab/Range.hpp"
00035 #include "moab/Interface.hpp"
00036 #include "ExoIIUtil.hpp"
00037 #include "MBTagConventions.hpp"
00038 #include "Internals.hpp"
00039 #include "moab/ReadUtilIface.hpp"
00040 #include "exodus_order.h"
00041 #include "moab/FileOptions.hpp"
00042 #include "moab/AdaptiveKDTree.hpp"
00043 #include "moab/CartVect.hpp"
00044 
00045 namespace moab
00046 {
00047 
00048 #define INS_ID( stringvar, prefix, id ) sprintf( stringvar, prefix, id )
00049 
00050 #define GET_DIM( ncdim, name, val )                                                                            \
00051     {                                                                                                          \
00052         int gdfail = nc_inq_dimid( ncFile, name, &( ncdim ) );                                                 \
00053         if( NC_NOERR == gdfail )                                                                               \
00054         {                                                                                                      \
00055             size_t tmp_val;                                                                                    \
00056             gdfail = nc_inq_dimlen( ncFile, ncdim, &tmp_val );                                                 \
00057             if( NC_NOERR != gdfail ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Couldn't get dimension length" ); } \
00058             else                                                                                               \
00059                 ( val ) = tmp_val;                                                                             \
00060         }                                                                                                      \
00061         else                                                                                                   \
00062             ( val ) = 0;                                                                                       \
00063     }
00064 
00065 #define GET_DIMB( ncdim, name, varname, id, val ) \
00066     INS_ID( name, varname, id );                  \
00067     GET_DIM( ncdim, name, val );
00068 
00069 #define GET_VAR( name, id, dims )                                                               \
00070     {                                                                                           \
00071         ( id )     = -1;                                                                        \
00072         int gvfail = nc_inq_varid( ncFile, name, &( id ) );                                     \
00073         if( NC_NOERR == gvfail )                                                                \
00074         {                                                                                       \
00075             int ndims;                                                                          \
00076             gvfail = nc_inq_varndims( ncFile, id, &ndims );                                     \
00077             if( NC_NOERR == gvfail )                                                            \
00078             {                                                                                   \
00079                 ( dims ).resize( ndims );                                                       \
00080                 gvfail = nc_inq_vardimid( ncFile, id, &( dims )[0] );                           \
00081                 if( NC_NOERR != gvfail )                                                        \
00082                 {                                                                               \
00083                     MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Couldn't get variable dimension IDs" ); \
00084                 }                                                                               \
00085             }                                                                                   \
00086         }                                                                                       \
00087     }
00088 
00089 #define GET_1D_INT_VAR( name, id, vals )                                                                               \
00090     {                                                                                                                  \
00091         GET_VAR( name, id, vals );                                                                                     \
00092         if( -1 != ( id ) )                                                                                             \
00093         {                                                                                                              \
00094             size_t ntmp;                                                                                               \
00095             int ivfail = nc_inq_dimlen( ncFile, ( vals )[0], &ntmp );                                                  \
00096             if( NC_NOERR != ivfail ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Couldn't get dimension length" ); }         \
00097             ( vals ).resize( ntmp );                                                                                   \
00098             size_t ntmp1 = 0;                                                                                          \
00099             ivfail       = nc_get_vara_int( ncFile, id, &ntmp1, &ntmp, &( vals )[0] );                                 \
00100             if( NC_NOERR != ivfail ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting variable " << ( name ) ); } \
00101         }                                                                                                              \
00102     }
00103 
00104 #define GET_1D_DBL_VAR( name, id, vals )                                                                               \
00105     {                                                                                                                  \
00106         std::vector< int > dum_dims;                                                                                   \
00107         GET_VAR( name, id, dum_dims );                                                                                 \
00108         if( -1 != ( id ) )                                                                                             \
00109         {                                                                                                              \
00110             size_t ntmp;                                                                                               \
00111             int dvfail = nc_inq_dimlen( ncFile, dum_dims[0], &ntmp );                                                  \
00112             if( NC_NOERR != dvfail ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Couldn't get dimension length" ); }         \
00113             ( vals ).resize( ntmp );                                                                                   \
00114             size_t ntmp1 = 0;                                                                                          \
00115             dvfail       = nc_get_vara_double( ncFile, id, &ntmp1, &ntmp, &( vals )[0] );                              \
00116             if( NC_NOERR != dvfail ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting variable " << ( name ) ); } \
00117         }                                                                                                              \
00118     }
00119 
00120 ReaderIface* ReadNCDF::factory( Interface* iface )
00121 {
00122     return new ReadNCDF( iface );
00123 }
00124 
00125 ReadNCDF::ReadNCDF( Interface* impl ) : mdbImpl( impl ), max_line_length( -1 ), max_str_length( -1 )
00126 {
00127     assert( impl != NULL );
00128     reset();
00129 
00130     impl->query_interface( readMeshIface );
00131 
00132     // Initialize in case tag_get_handle fails below
00133     mMaterialSetTag  = 0;
00134     mDirichletSetTag = 0;
00135     mNeumannSetTag   = 0;
00136     mHasMidNodesTag  = 0;
00137     mDistFactorTag   = 0;
00138     mQaRecordTag     = 0;
00139     mGlobalIdTag     = 0;
00140 
00141     //! Get and cache predefined tag handles
00142     ErrorCode result;
00143     const int negone = -1;
00144     result           = impl->tag_get_handle( MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mMaterialSetTag,
00145                                    MB_TAG_SPARSE | MB_TAG_CREAT, &negone );
00146     assert( MB_SUCCESS == result );
00147     result = impl->tag_get_handle( DIRICHLET_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mDirichletSetTag,
00148                                    MB_TAG_SPARSE | MB_TAG_CREAT, &negone );
00149     assert( MB_SUCCESS == result );
00150     result = impl->tag_get_handle( NEUMANN_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mNeumannSetTag,
00151                                    MB_TAG_SPARSE | MB_TAG_CREAT, &negone );
00152     assert( MB_SUCCESS == result );
00153     const int mids[] = { -1, -1, -1, -1 };
00154     result           = impl->tag_get_handle( HAS_MID_NODES_TAG_NAME, 4, MB_TYPE_INTEGER, mHasMidNodesTag,
00155                                    MB_TAG_SPARSE | MB_TAG_CREAT, mids );
00156     assert( MB_SUCCESS == result );
00157     result = impl->tag_get_handle( "distFactor", 0, MB_TYPE_DOUBLE, mDistFactorTag,
00158                                    MB_TAG_SPARSE | MB_TAG_VARLEN | MB_TAG_CREAT );
00159     assert( MB_SUCCESS == result );
00160     result = impl->tag_get_handle( "qaRecord", 0, MB_TYPE_OPAQUE, mQaRecordTag,
00161                                    MB_TAG_SPARSE | MB_TAG_VARLEN | MB_TAG_CREAT );
00162     assert( MB_SUCCESS == result );
00163     mGlobalIdTag = impl->globalId_tag();
00164 
00165     assert( MB_SUCCESS == result );
00166 #ifdef NDEBUG
00167     if( MB_SUCCESS == result ) {};  // Line to avoid compiler warning about unused variable
00168 #endif
00169     ncFile = 0;
00170 }
00171 
00172 void ReadNCDF::reset()
00173 {
00174     mCurrentMeshHandle = 0;
00175     vertexOffset       = 0;
00176 
00177     numberNodes_loading         = 0;
00178     numberElements_loading      = 0;
00179     numberElementBlocks_loading = 0;
00180     numberFaceBlocks_loading    = 0;
00181     numberNodeSets_loading      = 0;
00182     numberSideSets_loading      = 0;
00183     numberDimensions_loading    = 0;
00184 
00185     if( !blocksLoading.empty() ) blocksLoading.clear();
00186 
00187     if( !nodesInLoadedBlocks.empty() ) nodesInLoadedBlocks.clear();
00188 
00189     if( !faceBlocksLoading.empty() ) faceBlocksLoading.clear();
00190 }
00191 
00192 ReadNCDF::~ReadNCDF()
00193 {
00194     mdbImpl->release_interface( readMeshIface );
00195 }
00196 
00197 ErrorCode ReadNCDF::read_tag_values( const char* file_name, const char* tag_name, const FileOptions& /* opts */,
00198                                      std::vector< int >& id_array, const SubsetList* subset_list )
00199 {
00200     if( subset_list )
00201     {
00202         MB_SET_ERR( MB_UNSUPPORTED_OPERATION, "ExodusII reader supports subset read only by material ID" );
00203     }
00204 
00205     // Open netcdf/exodus file
00206     int fail = nc_open( file_name, 0, &ncFile );
00207     if( NC_NOWRITE != fail )
00208     {
00209         MB_SET_ERR( MB_FILE_DOES_NOT_EXIST, "ReadNCDF:: problem opening Netcdf/Exodus II file " << file_name );
00210     }
00211 
00212     // 1. Read the header
00213     ErrorCode rval = read_exodus_header();
00214     if( MB_FAILURE == rval ) return rval;
00215 
00216     int count = 0;
00217     const char* prop;
00218     const char* blocks   = "eb_prop1";
00219     const char* nodesets = "ns_prop1";
00220     const char* sidesets = "ss_prop1";
00221 
00222     if( !strcmp( tag_name, MATERIAL_SET_TAG_NAME ) )
00223     {
00224         count = numberElementBlocks_loading;
00225         prop  = blocks;
00226     }
00227     else if( !strcmp( tag_name, DIRICHLET_SET_TAG_NAME ) )
00228     {
00229         count = numberNodeSets_loading;
00230         prop  = nodesets;
00231     }
00232     else if( !strcmp( tag_name, NEUMANN_SET_TAG_NAME ) )
00233     {
00234         count = numberSideSets_loading;
00235         prop  = sidesets;
00236     }
00237     else
00238     {
00239         ncFile = 0;
00240         return MB_TAG_NOT_FOUND;
00241     }
00242 
00243     if( count )
00244     {
00245         int nc_var = -1;
00246         GET_1D_INT_VAR( prop, nc_var, id_array );
00247         if( !nc_var ) { MB_SET_ERR( MB_FAILURE, "Problem getting prop variable" ); }
00248     }
00249 
00250     // Close the file
00251     fail = nc_close( ncFile );
00252     if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "Trouble closing file" ); }
00253 
00254     ncFile = 0;
00255     return MB_SUCCESS;
00256 }
00257 
00258 ErrorCode ReadNCDF::load_file( const char* exodus_file_name, const EntityHandle* file_set, const FileOptions& opts,
00259                                const ReaderIface::SubsetList* subset_list, const Tag* file_id_tag )
00260 {
00261     ErrorCode status;
00262     int fail;
00263 
00264     int num_blocks            = 0;
00265     const int* blocks_to_load = 0;
00266     if( subset_list )
00267     {
00268         if( subset_list->tag_list_length > 1 || !strcmp( subset_list->tag_list[0].tag_name, MATERIAL_SET_TAG_NAME ) )
00269         {
00270             MB_SET_ERR( MB_UNSUPPORTED_OPERATION, "ExodusII reader supports subset read only by material ID" );
00271         }
00272         if( subset_list->num_parts )
00273         {
00274             MB_SET_ERR( MB_UNSUPPORTED_OPERATION, "ExodusII reader does not support mesh partitioning" );
00275         }
00276         blocks_to_load = subset_list->tag_list[0].tag_values;
00277         num_blocks     = subset_list->tag_list[0].num_tag_values;
00278     }
00279 
00280     // This function directs the reading of an exoii file, but doesn't do any of
00281     // the actual work
00282 
00283     // See if opts has tdata.
00284     ErrorCode rval;
00285     std::string s;
00286     rval = opts.get_str_option( "tdata", s );
00287     if( MB_SUCCESS == rval && !s.empty() )
00288         return update( exodus_file_name, opts, num_blocks, blocks_to_load, *file_set );
00289 
00290     reset();
00291 
00292     // 0. Open the file.
00293 
00294     // open netcdf/exodus file
00295     fail = nc_open( exodus_file_name, 0, &ncFile );
00296     if( NC_NOERR != fail )
00297     {
00298         MB_SET_ERR( MB_FILE_DOES_NOT_EXIST, "ReadNCDF:: problem opening Netcdf/Exodus II file " << exodus_file_name );
00299     }
00300 
00301     // 1. Read the header
00302     status = read_exodus_header();
00303     if( MB_FAILURE == status ) return status;
00304 
00305     status = mdbImpl->get_entities_by_handle( 0, initRange );
00306     if( MB_FAILURE == status ) return status;
00307 
00308     // 2. Read the nodes unless they've already been read before
00309     status = read_nodes( file_id_tag );
00310     if( MB_FAILURE == status ) return status;
00311 
00312     // 3.
00313     // extra for polyhedra blocks
00314     if( numberFaceBlocks_loading > 0 )
00315     {
00316         status = read_face_blocks_headers();
00317         if( MB_FAILURE == status ) return status;
00318     }
00319 
00320     status = read_block_headers( blocks_to_load, num_blocks );
00321     if( MB_FAILURE == status ) return status;
00322 
00323     // 4. Read elements (might not read them, depending on active blocks)
00324     if( numberFaceBlocks_loading > 0 )
00325     {
00326         status = read_polyhedra_faces();
00327         if( MB_FAILURE == status ) return status;
00328     }
00329     status = read_elements( file_id_tag );
00330     if( MB_FAILURE == status ) return status;
00331 
00332     // 5. Read global ids
00333     status = read_global_ids();
00334     if( MB_FAILURE == status ) return status;
00335 
00336     // 6. Read nodesets
00337     status = read_nodesets();
00338     if( MB_FAILURE == status ) return status;
00339 
00340     // 7. Read sidesets
00341     status = read_sidesets();
00342     if( MB_FAILURE == status ) return status;
00343 
00344     // 8. Read qa records
00345     if( file_set )
00346     {
00347         status = read_qa_records( *file_set );
00348         if( MB_FAILURE == status ) return status;
00349     }
00350 
00351     // What about properties???
00352 
00353     // Close the file
00354     fail = nc_close( ncFile );
00355     if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "Trouble closing file" ); }
00356 
00357     ncFile = 0;
00358     return MB_SUCCESS;
00359 }
00360 
00361 ErrorCode ReadNCDF::read_exodus_header()
00362 {
00363     CPU_WORD_SIZE = sizeof( double );  // With ExodusII version 2, all floats
00364     IO_WORD_SIZE  = sizeof( double );  // should be changed to doubles
00365 
00366     // NetCDF doesn't check its own limits on file read, so check
00367     // them here so it doesn't corrupt memory any more than absolutely
00368     // necessary.
00369     int ndims;
00370     int fail = nc_inq_ndims( ncFile, &ndims );
00371     if( NC_NOERR != fail || ndims > NC_MAX_DIMS )
00372     {
00373         MB_SET_ERR( MB_FAILURE, "ReadNCDF: File contains " << ndims << " dims but NetCDF library supports only "
00374                                                            << (int)NC_MAX_DIMS );
00375     }
00376     int nvars;
00377     fail = nc_inq_nvars( ncFile, &nvars );
00378     if( nvars > NC_MAX_VARS )
00379     {
00380         MB_SET_ERR( MB_FAILURE, "ReadNCDF: File contains " << nvars << " vars but NetCDF library supports only "
00381                                                            << (int)NC_MAX_VARS );
00382     }
00383 
00384     // Get the attributes
00385 
00386     // Get the word size, scalar value
00387     nc_type att_type;
00388     size_t att_len;
00389     fail = nc_inq_att( ncFile, NC_GLOBAL, "floating_point_word_size", &att_type, &att_len );
00390     if( NC_NOERR != fail )
00391     {
00392         MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting floating_point_word_size attribute" );
00393     }
00394     if( att_type != NC_INT || att_len != 1 )
00395     {
00396         MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Word size didn't have type int or size 1" );
00397     }
00398     fail = nc_get_att_int( ncFile, NC_GLOBAL, "floating_point_word_size", &IO_WORD_SIZE );
00399     if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Trouble getting word size" ); }
00400 
00401     // Exodus version
00402     fail = nc_inq_att( ncFile, NC_GLOBAL, "version", &att_type, &att_len );
00403     if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting version attribute" ); }
00404     if( att_type != NC_FLOAT || att_len != 1 )
00405     {
00406         MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Version didn't have type float or size 1" );
00407     }
00408     // float version = temp_att->as_float(0);
00409 
00410     // Read in initial variables
00411     int temp_dim;
00412     GET_DIM( temp_dim, "num_dim", numberDimensions_loading );
00413     GET_DIM( temp_dim, "num_nodes", numberNodes_loading );
00414     GET_DIM( temp_dim, "num_elem", numberElements_loading );
00415     GET_DIM( temp_dim, "num_el_blk", numberElementBlocks_loading );
00416     GET_DIM( temp_dim, "num_fa_blk", numberFaceBlocks_loading );  // for polyhedra blocks
00417     GET_DIM( temp_dim, "num_elem", numberElements_loading );
00418     GET_DIM( temp_dim, "num_node_sets", numberNodeSets_loading );
00419     GET_DIM( temp_dim, "num_side_sets", numberSideSets_loading );
00420     GET_DIM( temp_dim, "len_string", max_str_length );
00421     GET_DIM( temp_dim, "len_line", max_line_length );
00422 
00423     // Title; why are we even bothering if we're not going to keep it???
00424     fail = nc_inq_att( ncFile, NC_GLOBAL, "title", &att_type, &att_len );
00425     if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting title attribute" ); }
00426     if( att_type != NC_CHAR ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: title didn't have type char" ); }
00427     char* title = new char[att_len + 1];
00428     fail        = nc_get_att_text( ncFile, NC_GLOBAL, "title", title );
00429     if( NC_NOERR != fail )
00430     {
00431         delete[] title;
00432         MB_SET_ERR( MB_FAILURE, "ReadNCDF:: trouble getting title" );
00433     }
00434     delete[] title;
00435 
00436     return MB_SUCCESS;
00437 }
00438 
00439 ErrorCode ReadNCDF::read_nodes( const Tag* file_id_tag )
00440 {
00441     // Read the nodes into memory
00442 
00443     assert( 0 != ncFile );
00444 
00445     // Create a sequence to hold the node coordinates
00446     // Get the current number of entities and start at the next slot
00447 
00448     EntityHandle node_handle = 0;
00449     std::vector< double* > arrays;
00450     readMeshIface->get_node_coords( 3, numberNodes_loading, MB_START_ID, node_handle, arrays );
00451 
00452     vertexOffset = ID_FROM_HANDLE( node_handle ) - MB_START_ID;
00453 
00454     // Read in the coordinates
00455     int fail;
00456     int coord = 0;
00457     nc_inq_varid( ncFile, "coord", &coord );
00458 
00459     // Single var for all coords
00460     size_t start[2] = { 0, 0 }, count[2] = { 1, static_cast< size_t >( numberNodes_loading ) };
00461     if( coord )
00462     {
00463         for( int d = 0; d < numberDimensions_loading; ++d )
00464         {
00465             start[0] = d;
00466             fail     = nc_get_vara_double( ncFile, coord, start, count, arrays[d] );
00467             if( NC_NOERR != fail )
00468             {
00469                 MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting " << (char)( 'x' + d ) << " coord array" );
00470             }
00471         }
00472     }
00473     // Var for each coord
00474     else
00475     {
00476         char varname[] = "coord ";
00477         for( int d = 0; d < numberDimensions_loading; ++d )
00478         {
00479             varname[5] = 'x' + (char)d;
00480             fail       = nc_inq_varid( ncFile, varname, &coord );
00481             if( NC_NOERR != fail )
00482             {
00483                 MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting " << (char)( 'x' + d ) << " coord variable" );
00484             }
00485 
00486             fail = nc_get_vara_double( ncFile, coord, start, &count[1], arrays[d] );
00487             if( NC_NOERR != fail )
00488             {
00489                 MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting " << (char)( 'x' + d ) << " coord array" );
00490             }
00491         }
00492     }
00493 
00494     // Zero out any coord values that are in database but not in file
00495     // (e.g. if MOAB has 3D coords but file is 2D then set Z coords to zero.)
00496     for( unsigned d = numberDimensions_loading; d < arrays.size(); ++d )
00497         std::fill( arrays[d], arrays[d] + numberNodes_loading, 0.0 );
00498 
00499     if( file_id_tag )
00500     {
00501         Range nodes;
00502         nodes.insert( node_handle, node_handle + numberNodes_loading - 1 );
00503         readMeshIface->assign_ids( *file_id_tag, nodes, vertexOffset );
00504     }
00505 
00506     return MB_SUCCESS;
00507 }
00508 
00509 ErrorCode ReadNCDF::read_block_headers( const int* blocks_to_load, const int num_blocks )
00510 {
00511     // Get the element block ids; keep this in a separate list,
00512     // which is not offset by blockIdOffset; this list used later for
00513     // reading block connectivity
00514 
00515     // Get the ids of all the blocks of this file we're reading in
00516     std::vector< int > block_ids( numberElementBlocks_loading );
00517     int nc_block_ids = -1;
00518     GET_1D_INT_VAR( "eb_prop1", nc_block_ids, block_ids );
00519     if( -1 == nc_block_ids ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting eb_prop1 variable" ); }
00520 
00521     int exodus_id = 1;
00522 
00523     // If the active_block_id_list is NULL all blocks are active.
00524     if( NULL == blocks_to_load || 0 == num_blocks ) { blocks_to_load = &block_ids[0]; }
00525 
00526     std::vector< int > new_blocks( blocks_to_load, blocks_to_load + numberElementBlocks_loading );
00527 
00528     std::vector< int >::iterator iter, end_iter;
00529     iter     = block_ids.begin();
00530     end_iter = block_ids.end();
00531 
00532     // Read header information and initialize header-type block information
00533     int temp_dim;
00534     std::vector< char > temp_string_storage( max_str_length + 1 );
00535     char* temp_string = &temp_string_storage[0];
00536     int block_seq_id  = 1;
00537 
00538     for( ; iter != end_iter; ++iter, block_seq_id++ )
00539     {
00540         int num_elements;
00541 
00542         GET_DIMB( temp_dim, temp_string, "num_el_in_blk%d", block_seq_id, num_elements );
00543 
00544         // Don't read element type string for now, since it's an attrib
00545         // on the connectivity
00546         // Get the entity type corresponding to this ExoII element type
00547         // ExoIIElementType elem_type =
00548         // ExoIIUtil::static_element_name_to_type(element_type_string);
00549 
00550         // Tag each element block(mesh set) with enum for ElementType (SHELL, QUAD4, ....etc)
00551         ReadBlockData block_data;
00552         block_data.elemType    = EXOII_MAX_ELEM_TYPE;
00553         block_data.blockId     = *iter;
00554         block_data.startExoId  = exodus_id;
00555         block_data.numElements = num_elements;
00556 
00557         // If block is in 'blocks_to_load'----load it!
00558         if( std::find( new_blocks.begin(), new_blocks.end(), *iter ) != new_blocks.end() )
00559             block_data.reading_in = true;
00560         else
00561             block_data.reading_in = false;
00562 
00563         blocksLoading.push_back( block_data );
00564         exodus_id += num_elements;
00565     }
00566 
00567     return MB_SUCCESS;
00568 }
00569 
00570 ErrorCode ReadNCDF::read_face_blocks_headers()
00571 {
00572     // Get the ids of all the blocks of this file we're reading in
00573     std::vector< int > fblock_ids( numberFaceBlocks_loading );
00574     int nc_fblock_ids = -1;
00575     GET_1D_INT_VAR( "fa_prop1", nc_fblock_ids, fblock_ids );
00576     if( -1 == nc_fblock_ids ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting fa_prop1 variable" ); }
00577 
00578     int temp_dim;
00579     std::vector< char > temp_string_storage( max_str_length + 1 );
00580     char* temp_string = &temp_string_storage[0];
00581     // int block_seq_id = 1;
00582     int exodus_id = 1;
00583 
00584     for( int i = 1; i <= numberFaceBlocks_loading; i++ )
00585     {
00586         int num_elements;
00587 
00588         GET_DIMB( temp_dim, temp_string, "num_fa_in_blk%d", i, num_elements );
00589 
00590         // Don't read element type string for now, since it's an attrib
00591         // on the connectivity
00592         // Get the entity type corresponding to this ExoII element type
00593         // ExoIIElementType elem_type =
00594         // ExoIIUtil::static_element_name_to_type(element_type_string);
00595 
00596         // Tag each element block(mesh set) with enum for ElementType (SHELL, QUAD4, ....etc)
00597         ReadFaceBlockData fblock_data;
00598         fblock_data.faceBlockId = i;
00599         fblock_data.startExoId  = exodus_id;
00600         fblock_data.numElements = num_elements;
00601 
00602         faceBlocksLoading.push_back( fblock_data );
00603         exodus_id += num_elements;
00604     }
00605     return MB_SUCCESS;
00606 }
00607 ErrorCode ReadNCDF::read_polyhedra_faces()
00608 {
00609     int temp_dim;
00610     std::vector< char > temp_string_storage( max_str_length + 1 );
00611     char* temp_string = &temp_string_storage[0];
00612     nodesInLoadedBlocks.resize( numberNodes_loading + 1 );
00613     size_t start[1] = { 0 }, count[1];  // one dim arrays only here!
00614     std::vector< ReadFaceBlockData >::iterator this_it;
00615     this_it = faceBlocksLoading.begin();
00616 
00617     int fblock_seq_id = 1;
00618     std::vector< int > dims;
00619     int nc_var, fail;
00620 
00621     for( ; this_it != faceBlocksLoading.end(); ++this_it, fblock_seq_id++ )
00622     {
00623         int num_fa_in_blk;
00624         GET_DIMB( temp_dim, temp_string, "num_fa_in_blk%d", fblock_seq_id, num_fa_in_blk );
00625         int num_nod_per_fa;
00626         GET_DIMB( temp_dim, temp_string, "num_nod_per_fa%d", fblock_seq_id, num_nod_per_fa );
00627         // Get the ncdf connect variable and the element type
00628         INS_ID( temp_string, "fbconn%d", fblock_seq_id );
00629         GET_VAR( temp_string, nc_var, dims );
00630         std::vector< int > fbconn;
00631         fbconn.resize( num_nod_per_fa );
00632         count[0] = num_nod_per_fa;
00633 
00634         fail = nc_get_vara_int( ncFile, nc_var, start, count, &fbconn[0] );
00635         if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting face polyhedra connectivity " ); }
00636         std::vector< int > fbepecnt;
00637         fbepecnt.resize( num_fa_in_blk );
00638         INS_ID( temp_string, "fbepecnt%d", fblock_seq_id );
00639         GET_VAR( temp_string, nc_var, dims );
00640         count[0] = num_fa_in_blk;
00641         fail     = nc_get_vara_int( ncFile, nc_var, start, count, &fbepecnt[0] );
00642         if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting face polyhedra connectivity " ); }
00643         // now we can create some polygons
00644         std::vector< EntityHandle > polyConn;
00645         int ix = 0;
00646         for( int i = 0; i < num_fa_in_blk; i++ )
00647         {
00648             polyConn.resize( fbepecnt[i] );
00649             for( int j = 0; j < fbepecnt[i]; j++ )
00650             {
00651                 nodesInLoadedBlocks[fbconn[ix]] = 1;
00652                 polyConn[j]                     = vertexOffset + fbconn[ix++];
00653             }
00654             EntityHandle newp;
00655             /*
00656              *  ErrorCode create_element(const EntityType type,
00657                                          const EntityHandle *connectivity,
00658                                          const int num_vertices,
00659                                          EntityHandle &element_handle)
00660              */
00661             if( mdbImpl->create_element( MBPOLYGON, &polyConn[0], fbepecnt[i], newp ) != MB_SUCCESS ) return MB_FAILURE;
00662 
00663             // add the element in order
00664             polyfaces.push_back( newp );
00665         }
00666     }
00667     return MB_SUCCESS;
00668 }
00669 
00670 ErrorCode ReadNCDF::read_elements( const Tag* file_id_tag )
00671 {
00672     // Read in elements
00673 
00674     int result = 0;
00675 
00676     // Initialize the nodeInLoadedBlocks vector
00677     if( nodesInLoadedBlocks.size() < ( size_t )( numberNodes_loading + 1 ) )
00678         nodesInLoadedBlocks.resize( numberNodes_loading + 1 );  // this could be repeated?
00679     memset( &nodesInLoadedBlocks[0], 0, ( numberNodes_loading + 1 ) * sizeof( char ) );
00680 
00681     std::vector< ReadBlockData >::iterator this_it;
00682     this_it = blocksLoading.begin();
00683 
00684     int temp_dim;
00685     std::vector< char > temp_string_storage( max_str_length + 1 );
00686     char* temp_string = &temp_string_storage[0];
00687 
00688     int nc_var;
00689     int block_seq_id = 1;
00690     std::vector< int > dims;
00691     size_t start[2] = { 0, 0 }, count[2];
00692 
00693     for( ; this_it != blocksLoading.end(); ++this_it, block_seq_id++ )
00694     {
00695         // If this block isn't to be read in --- continue
00696         if( !( *this_it ).reading_in ) continue;
00697 
00698         // Get some information about this block
00699         int block_id       = ( *this_it ).blockId;
00700         EntityHandle* conn = 0;
00701 
00702         // Get the ncdf connect variable and the element type
00703         INS_ID( temp_string, "connect%d", block_seq_id );
00704         GET_VAR( temp_string, nc_var, dims );
00705         if( -1 == nc_var || 0 == nc_var )
00706         {  // try other var, for polyhedra blocks
00707             // it could be polyhedra block, defined by fbconn and NFACED attribute
00708             INS_ID( temp_string, "facconn%d", block_seq_id );
00709             GET_VAR( temp_string, nc_var, dims );
00710 
00711             if( -1 == nc_var || 0 == nc_var )
00712             {
00713                 MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting connec or faccon variable" );
00714             }
00715         }
00716         nc_type att_type;
00717         size_t att_len;
00718         int fail = nc_inq_att( ncFile, nc_var, "elem_type", &att_type, &att_len );
00719         if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting elem type attribute" ); }
00720         std::vector< char > dum_str( att_len + 1 );
00721         fail = nc_get_att_text( ncFile, nc_var, "elem_type", &dum_str[0] );
00722         if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting elem type" ); }
00723         ExoIIElementType elem_type = ExoIIUtil::static_element_name_to_type( &dum_str[0] );
00724         ( *this_it ).elemType      = elem_type;
00725 
00726         int verts_per_element    = ExoIIUtil::VerticesPerElement[( *this_it ).elemType];
00727         int number_nodes         = ( *this_it ).numElements * verts_per_element;
00728         const EntityType mb_type = ExoIIUtil::ExoIIElementMBEntity[( *this_it ).elemType];
00729 
00730         if( mb_type == MBPOLYGON )
00731         {
00732             // need to read another variable, num_nod_per_el, and
00733             int num_nodes_poly_conn;
00734             GET_DIMB( temp_dim, temp_string, "num_nod_per_el%d", block_seq_id, num_nodes_poly_conn );
00735             // read the connec1 array, which is of dimensions num_nodes_poly_conn (only one )
00736             std::vector< int > connec;
00737             connec.resize( num_nodes_poly_conn );
00738             count[0] = num_nodes_poly_conn;
00739 
00740             fail = nc_get_vara_int( ncFile, nc_var, start, count, &connec[0] );
00741             if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting polygon connectivity " ); }
00742             // ebepecnt1
00743             std::vector< int > ebec;
00744             ebec.resize( this_it->numElements );
00745             // Get the ncdf connect variable and the element type
00746             INS_ID( temp_string, "ebepecnt%d", block_seq_id );
00747             GET_VAR( temp_string, nc_var, dims );
00748             count[0] = this_it->numElements;
00749             fail     = nc_get_vara_int( ncFile, nc_var, start, count, &ebec[0] );
00750             if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting polygon nodes per elem " ); }
00751             EntityHandle ms_handle;
00752             if( mdbImpl->create_meshset( MESHSET_SET | MESHSET_TRACK_OWNER, ms_handle ) != MB_SUCCESS )
00753                 return MB_FAILURE;
00754             // create polygons one by one, and put them in the list
00755             // also need to read the number of nodes for each polygon, then create
00756             std::vector< EntityHandle > polyConn;
00757             int ix = 0;
00758             for( int i = 0; i < this_it->numElements; i++ )
00759             {
00760                 polyConn.resize( ebec[i] );
00761                 for( int j = 0; j < ebec[i]; j++ )
00762                 {
00763                     nodesInLoadedBlocks[connec[ix]] = 1;
00764                     polyConn[j]                     = vertexOffset + connec[ix++];
00765                 }
00766                 EntityHandle newp;
00767                 /*
00768                  *  ErrorCode create_element(const EntityType type,
00769                                              const EntityHandle *connectivity,
00770                                              const int num_vertices,
00771                                              EntityHandle &element_handle)
00772                  */
00773                 if( mdbImpl->create_element( MBPOLYGON, &polyConn[0], ebec[i], newp ) != MB_SUCCESS ) return MB_FAILURE;
00774 
00775                 // add the element in order
00776                 this_it->polys.push_back( newp );
00777                 if( mdbImpl->add_entities( ms_handle, &newp, 1 ) != MB_SUCCESS ) return MB_FAILURE;
00778             }
00779             // Set the block id with an offset
00780             if( mdbImpl->tag_set_data( mMaterialSetTag, &ms_handle, 1, &block_id ) != MB_SUCCESS ) return MB_FAILURE;
00781             if( mdbImpl->tag_set_data( mGlobalIdTag, &ms_handle, 1, &block_id ) != MB_SUCCESS ) return MB_FAILURE;
00782         }
00783         else if( mb_type == MBPOLYHEDRON )
00784         {
00785             // need to read another variable, num_fac_per_el
00786             int num_fac_per_el;
00787             GET_DIMB( temp_dim, temp_string, "num_fac_per_el%d", block_seq_id, num_fac_per_el );
00788             // read the fbconn array, which is of dimensions num_nod_per_fa (only one dimension)
00789             std::vector< int > facconn;
00790             facconn.resize( num_fac_per_el );
00791             count[0] = num_fac_per_el;
00792 
00793             fail = nc_get_vara_int( ncFile, nc_var, start, count, &facconn[0] );
00794             if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting polygon connectivity " ); }
00795             // ebepecnt1
00796             std::vector< int > ebepecnt;
00797             ebepecnt.resize( this_it->numElements );
00798             // Get the ncdf connect variable and the element type
00799             INS_ID( temp_string, "ebepecnt%d", block_seq_id );
00800             GET_VAR( temp_string, nc_var, dims );
00801             count[0] = this_it->numElements;
00802             fail     = nc_get_vara_int( ncFile, nc_var, start, count, &ebepecnt[0] );
00803             if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting polygon nodes per elem " ); }
00804             EntityHandle ms_handle;
00805             if( mdbImpl->create_meshset( MESHSET_SET | MESHSET_TRACK_OWNER, ms_handle ) != MB_SUCCESS )
00806                 return MB_FAILURE;
00807             // create polygons one by one, and put them in the list
00808             // also need to read the number of nodes for each polygon, then create
00809             std::vector< EntityHandle > polyConn;
00810             int ix = 0;
00811             for( int i = 0; i < this_it->numElements; i++ )
00812             {
00813                 polyConn.resize( ebepecnt[i] );
00814                 for( int j = 0; j < ebepecnt[i]; j++ )
00815                 {
00816                     polyConn[j] = polyfaces[facconn[ix++] - 1];
00817                 }
00818                 EntityHandle newp;
00819                 /*
00820                  *  ErrorCode create_element(const EntityType type,
00821                                              const EntityHandle *connectivity,
00822                                              const int num_vertices,
00823                                              EntityHandle &element_handle)
00824                  */
00825                 if( mdbImpl->create_element( MBPOLYHEDRON, &polyConn[0], ebepecnt[i], newp ) != MB_SUCCESS )
00826                     return MB_FAILURE;
00827 
00828                 // add the element in order
00829                 this_it->polys.push_back( newp );
00830                 if( mdbImpl->add_entities( ms_handle, &newp, 1 ) != MB_SUCCESS ) return MB_FAILURE;
00831             }
00832             // Set the block id with an offset
00833             if( mdbImpl->tag_set_data( mMaterialSetTag, &ms_handle, 1, &block_id ) != MB_SUCCESS ) return MB_FAILURE;
00834             if( mdbImpl->tag_set_data( mGlobalIdTag, &ms_handle, 1, &block_id ) != MB_SUCCESS ) return MB_FAILURE;
00835         }
00836         else  // this is regular
00837         {
00838             // Allocate an array to read in connectivity data
00839             readMeshIface->get_element_connect( this_it->numElements, verts_per_element, mb_type, this_it->startExoId,
00840                                                 this_it->startMBId, conn );
00841 
00842             // Create a range for this sequence of elements
00843             EntityHandle start_range, end_range;
00844             start_range = ( *this_it ).startMBId;
00845             end_range   = start_range + ( *this_it ).numElements - 1;
00846 
00847             Range new_range( start_range, end_range );
00848             // Range<EntityHandle> new_range((*this_it).startMBId,
00849             //                              (*this_it).startMBId + (*this_it).numElements - 1);
00850 
00851             // Create a MBSet for this block and set the material tag
00852 
00853             EntityHandle ms_handle;
00854             if( mdbImpl->create_meshset( MESHSET_SET | MESHSET_TRACK_OWNER, ms_handle ) != MB_SUCCESS )
00855                 return MB_FAILURE;
00856 
00857             if( mdbImpl->add_entities( ms_handle, new_range ) != MB_SUCCESS ) return MB_FAILURE;
00858 
00859             int mid_nodes[4];
00860             CN::HasMidNodes( mb_type, verts_per_element, mid_nodes );
00861             if( mdbImpl->tag_set_data( mHasMidNodesTag, &ms_handle, 1, mid_nodes ) != MB_SUCCESS ) return MB_FAILURE;
00862 
00863             // Just a check because the following code won't work if this case fails
00864             assert( sizeof( EntityHandle ) >= sizeof( int ) );
00865 
00866             // tmp_ptr is of type int* and points at the same place as conn
00867             int* tmp_ptr = reinterpret_cast< int* >( conn );
00868 
00869             // Read the connectivity into that memory,  this will take only part of the array
00870             // 1/2 if sizeof(EntityHandle) == 64 bits.
00871             count[0] = this_it->numElements;
00872             count[1] = verts_per_element;
00873             fail     = nc_get_vara_int( ncFile, nc_var, start, count, tmp_ptr );
00874             if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting connectivity" ); }
00875 
00876             // Convert from exodus indices to vertex handles.
00877             // Iterate backwards in case handles are larger than ints.
00878             for( int i = number_nodes - 1; i >= 0; --i )
00879             {
00880                 if( (unsigned)tmp_ptr[i] >= nodesInLoadedBlocks.size() )
00881                 {
00882                     MB_SET_ERR( MB_FAILURE, "Invalid node ID in block connectivity" );
00883                 }
00884                 nodesInLoadedBlocks[tmp_ptr[i]] = 1;
00885                 conn[i]                         = static_cast< EntityHandle >( tmp_ptr[i] ) + vertexOffset;
00886             }
00887 
00888             // Adjust connectivity order if necessary
00889             const int* reorder = exodus_elem_order_map[mb_type][verts_per_element];
00890             if( reorder ) ReadUtilIface::reorder( reorder, conn, this_it->numElements, verts_per_element );
00891 
00892             readMeshIface->update_adjacencies( ( *this_it ).startMBId, ( *this_it ).numElements,
00893                                                ExoIIUtil::VerticesPerElement[( *this_it ).elemType], conn );
00894 
00895             if( result == -1 )
00896             {
00897                 MB_SET_ERR( MB_FAILURE, "ReadNCDF:: error getting element connectivity for block " << block_id );
00898             }
00899 
00900             // Set the block id with an offset
00901             if( mdbImpl->tag_set_data( mMaterialSetTag, &ms_handle, 1, &block_id ) != MB_SUCCESS ) return MB_FAILURE;
00902             if( mdbImpl->tag_set_data( mGlobalIdTag, &ms_handle, 1, &block_id ) != MB_SUCCESS ) return MB_FAILURE;
00903 
00904             if( file_id_tag )
00905             {
00906                 Range range;
00907                 range.insert( this_it->startMBId, this_it->startMBId + this_it->numElements - 1 );
00908                 readMeshIface->assign_ids( *file_id_tag, range, this_it->startExoId );
00909             }
00910         }  // end regular block (not polygon)
00911     }
00912 
00913     return MB_SUCCESS;
00914 }
00915 
00916 ErrorCode ReadNCDF::read_global_ids()
00917 {
00918     // Read in the map from the exodus file
00919     std::vector< int > ids( std::max( numberElements_loading, numberNodes_loading ) );
00920 
00921     int varid = -1;
00922     GET_1D_INT_VAR( "elem_map", varid, ids );
00923     if( -1 != varid )
00924     {
00925         std::vector< ReadBlockData >::iterator iter;
00926         int id_pos = 0;
00927         for( iter = blocksLoading.begin(); iter != blocksLoading.end(); ++iter )
00928         {
00929             if( iter->reading_in )
00930             {
00931                 if( iter->elemType != EXOII_POLYGON && iter->elemType != EXOII_POLYHEDRON )
00932                 {
00933                     if( iter->startMBId != 0 )
00934                     {
00935                         Range range( iter->startMBId, iter->startMBId + iter->numElements - 1 );
00936                         ErrorCode error = mdbImpl->tag_set_data( mGlobalIdTag, range, &ids[id_pos] );
00937                         if( error != MB_SUCCESS ) return error;
00938                         id_pos += iter->numElements;
00939                     }
00940                     else
00941                         return MB_FAILURE;
00942                 }
00943                 else  // polygons or polyhedrons; use id from elements
00944                 {
00945                     for( std::vector< EntityHandle >::iterator eit = iter->polys.begin(); eit != iter->polys.end();
00946                          eit++, id_pos++ )
00947                     {
00948                         EntityHandle peh = *eit;
00949                         if( mdbImpl->tag_set_data( mGlobalIdTag, &peh, 1, &ids[id_pos] ) != MB_SUCCESS )
00950                             return MB_FAILURE;
00951                     }
00952                 }
00953             }
00954         }
00955     }
00956 
00957     // Read in node map next
00958     varid = -1;
00959     GET_1D_INT_VAR( "node_num_map", varid, ids );
00960     if( -1 != varid )
00961     {
00962         Range range( MB_START_ID + vertexOffset, MB_START_ID + vertexOffset + numberNodes_loading - 1 );
00963         ErrorCode error = mdbImpl->tag_set_data( mGlobalIdTag, range, &ids[0] );
00964         if( MB_SUCCESS != error ) MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem setting node global ids" );
00965     }
00966 
00967     return MB_SUCCESS;
00968 }
00969 
00970 ErrorCode ReadNCDF::read_nodesets()
00971 {
00972     // Read in the nodesets for the model
00973 
00974     if( 0 == numberNodeSets_loading ) return MB_SUCCESS;
00975     std::vector< int > id_array( numberNodeSets_loading );
00976 
00977     // Read in the nodeset ids
00978     int nc_var;
00979     GET_1D_INT_VAR( "ns_prop1", nc_var, id_array );
00980     if( -1 == nc_var ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting ns_prop1 variable" ); }
00981 
00982     // Use a vector of ints to read node handles
00983     std::vector< int > node_handles;
00984 
00985     int i;
00986     std::vector< char > temp_string_storage( max_str_length + 1 );
00987     char* temp_string = &temp_string_storage[0];
00988     int temp_dim;
00989     for( i = 0; i < numberNodeSets_loading; i++ )
00990     {
00991         // Get nodeset parameters
00992         int number_nodes_in_set;
00993         int number_dist_factors_in_set;
00994 
00995         GET_DIMB( temp_dim, temp_string, "num_nod_ns%d", i + 1, number_nodes_in_set );
00996         GET_DIMB( temp_dim, temp_string, "num_df_ns%d", i + 1, number_dist_factors_in_set );
00997 
00998         // Need to new a vector to store dist. factors
00999         // This vector * gets stored as a tag on the sideset meshset
01000         std::vector< double > temp_dist_factor_vector( number_nodes_in_set );
01001         if( number_dist_factors_in_set != 0 )
01002         {
01003             INS_ID( temp_string, "dist_fact_ns%d", i + 1 );
01004             GET_1D_DBL_VAR( temp_string, temp_dim, temp_dist_factor_vector );
01005             if( -1 == temp_dim ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting dist fact variable" ); }
01006         }
01007 
01008         // Size new arrays and get ids and distribution factors
01009         if( node_handles.size() < (unsigned int)number_nodes_in_set )
01010         {
01011             node_handles.reserve( number_nodes_in_set );
01012             node_handles.resize( number_nodes_in_set );
01013         }
01014 
01015         INS_ID( temp_string, "node_ns%d", i + 1 );
01016         int temp_var = -1;
01017         GET_1D_INT_VAR( temp_string, temp_var, node_handles );
01018         if( -1 == temp_var ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting nodeset node variable" ); }
01019 
01020         // Maybe there is already a nodesets meshset here we can append to
01021         Range child_meshsets;
01022         if( mdbImpl->get_entities_by_handle( 0, child_meshsets ) != MB_SUCCESS ) return MB_FAILURE;
01023 
01024         child_meshsets = subtract( child_meshsets, initRange );
01025 
01026         Range::iterator iter, end_iter;
01027         iter     = child_meshsets.begin();
01028         end_iter = child_meshsets.end();
01029 
01030         EntityHandle ns_handle = 0;
01031         for( ; iter != end_iter; ++iter )
01032         {
01033             int nodeset_id;
01034             if( mdbImpl->tag_get_data( mDirichletSetTag, &( *iter ), 1, &nodeset_id ) != MB_SUCCESS ) continue;
01035 
01036             if( id_array[i] == nodeset_id )
01037             {
01038                 // Found the meshset
01039                 ns_handle = *iter;
01040                 break;
01041             }
01042         }
01043 
01044         std::vector< EntityHandle > nodes_of_nodeset;
01045         if( ns_handle )
01046             if( mdbImpl->get_entities_by_handle( ns_handle, nodes_of_nodeset, true ) != MB_SUCCESS ) return MB_FAILURE;
01047 
01048         // Make these into entity handles
01049         // TODO: could we have read it into EntityHandle sized array in the first place?
01050         int j, temp;
01051         std::vector< EntityHandle > nodes;
01052         std::vector< double > dist_factor_vector;
01053         for( j = 0; j < number_nodes_in_set; j++ )
01054         {
01055             // See if this node is one we're currently reading in
01056             if( nodesInLoadedBlocks[node_handles[j]] == 1 )
01057             {
01058                 // Make sure that it already isn't in a nodeset
01059                 unsigned int node_id = CREATE_HANDLE( MBVERTEX, node_handles[j] + vertexOffset, temp );
01060                 if( !ns_handle ||
01061                     std::find( nodes_of_nodeset.begin(), nodes_of_nodeset.end(), node_id ) == nodes_of_nodeset.end() )
01062                 {
01063                     nodes.push_back( node_id );
01064 
01065                     if( number_dist_factors_in_set != 0 ) dist_factor_vector.push_back( temp_dist_factor_vector[j] );
01066                 }
01067             }
01068         }
01069 
01070         // No nodes to add
01071         if( nodes.empty() ) continue;
01072 
01073         // If there was no meshset found --> create one
01074         if( ns_handle == 0 )
01075         {
01076             if( mdbImpl->create_meshset( MESHSET_ORDERED | MESHSET_TRACK_OWNER, ns_handle ) != MB_SUCCESS )
01077                 return MB_FAILURE;
01078 
01079             // Set a tag signifying dirichlet bc
01080             // TODO: create this tag another way
01081 
01082             int nodeset_id = id_array[i];
01083             if( mdbImpl->tag_set_data( mDirichletSetTag, &ns_handle, 1, &nodeset_id ) != MB_SUCCESS ) return MB_FAILURE;
01084             if( mdbImpl->tag_set_data( mGlobalIdTag, &ns_handle, 1, &nodeset_id ) != MB_SUCCESS ) return MB_FAILURE;
01085 
01086             if( !dist_factor_vector.empty() )
01087             {
01088                 int size         = dist_factor_vector.size();
01089                 const void* data = &dist_factor_vector[0];
01090                 if( mdbImpl->tag_set_by_ptr( mDistFactorTag, &ns_handle, 1, &data, &size ) != MB_SUCCESS )
01091                     return MB_FAILURE;
01092             }
01093         }
01094         else if( !dist_factor_vector.empty() )
01095         {
01096             // Append dist factors to vector
01097             const void* ptr = 0;
01098             int size        = 0;
01099             if( mdbImpl->tag_get_by_ptr( mDistFactorTag, &ns_handle, 1, &ptr, &size ) != MB_SUCCESS ) return MB_FAILURE;
01100 
01101             const double* data = reinterpret_cast< const double* >( ptr );
01102             dist_factor_vector.reserve( dist_factor_vector.size() + size );
01103             std::copy( data, data + size, std::back_inserter( dist_factor_vector ) );
01104             size = dist_factor_vector.size();
01105             ptr  = &dist_factor_vector[0];
01106             if( mdbImpl->tag_set_by_ptr( mDistFactorTag, &ns_handle, 1, &ptr, &size ) != MB_SUCCESS ) return MB_FAILURE;
01107         }
01108 
01109         // Add the nodes to the meshset
01110         if( mdbImpl->add_entities( ns_handle, &nodes[0], nodes.size() ) != MB_SUCCESS ) return MB_FAILURE;
01111     }
01112 
01113     return MB_SUCCESS;
01114 }
01115 
01116 ErrorCode ReadNCDF::read_sidesets()
01117 {
01118     // Uhhh if you read the same file (previously_read==true) then blow away the
01119     // sidesets pertaining to this file?  How do you do that? If you read a
01120     // new file make sure all the offsets are correct.
01121 
01122     // If not loading any sidesets -- exit
01123     if( 0 == numberSideSets_loading ) return MB_SUCCESS;
01124 
01125     // Read in the sideset ids
01126     std::vector< int > id_array( numberSideSets_loading );
01127     int temp_var = -1;
01128     GET_1D_INT_VAR( "ss_prop1", temp_var, id_array );
01129     if( -1 == temp_var ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting ss_prop1 variable" ); }
01130 
01131     // Create a side set for each one
01132     int number_sides_in_set;
01133     int number_dist_factors_in_set;
01134 
01135     // Maybe there is already a sidesets meshset here we can append to
01136     Range child_meshsets;
01137     if( mdbImpl->get_entities_by_type( 0, MBENTITYSET, child_meshsets ) != MB_SUCCESS ) return MB_FAILURE;
01138 
01139     child_meshsets = subtract( child_meshsets, initRange );
01140 
01141     Range::iterator iter, end_iter;
01142 
01143     int i;
01144     std::vector< char > temp_string_storage( max_str_length + 1 );
01145     char* temp_string = &temp_string_storage[0];
01146     int temp_dim;
01147     for( i = 0; i < numberSideSets_loading; i++ )
01148     {
01149         // Get sideset parameters
01150         GET_DIMB( temp_dim, temp_string, "num_side_ss%d", i + 1, number_sides_in_set );
01151         GET_DIMB( temp_dim, temp_string, "num_df_ss%d", i + 1, number_dist_factors_in_set );
01152 
01153         // Size new arrays and get element and side lists
01154         std::vector< int > side_list( number_sides_in_set );
01155         std::vector< int > element_list( number_sides_in_set );
01156         INS_ID( temp_string, "side_ss%d", i + 1 );
01157         temp_var = -1;
01158         GET_1D_INT_VAR( temp_string, temp_var, side_list );
01159         if( -1 == temp_var ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting sideset side variable" ); }
01160 
01161         INS_ID( temp_string, "elem_ss%d", i + 1 );
01162         temp_var = -1;
01163         GET_1D_INT_VAR( temp_string, temp_var, element_list );
01164         if( -1 == temp_var ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting sideset elem variable" ); }
01165 
01166         std::vector< double > temp_dist_factor_vector;
01167         std::vector< EntityHandle > entities_to_add, reverse_entities;
01168         // Create the sideset entities
01169         if( create_ss_elements( &element_list[0], &side_list[0], number_sides_in_set, number_dist_factors_in_set,
01170                                 entities_to_add, reverse_entities, temp_dist_factor_vector, i + 1 ) != MB_SUCCESS )
01171             return MB_FAILURE;
01172 
01173         // If there are elements to add
01174         if( !entities_to_add.empty() || !reverse_entities.empty() )
01175         {
01176             iter     = child_meshsets.begin();
01177             end_iter = child_meshsets.end();
01178 
01179             EntityHandle ss_handle = 0;
01180             for( ; iter != end_iter; ++iter )
01181             {
01182                 int sideset_id;
01183                 if( mdbImpl->tag_get_data( mNeumannSetTag, &( *iter ), 1, &sideset_id ) != MB_SUCCESS ) continue;
01184 
01185                 if( id_array[i] == sideset_id )
01186                 {
01187                     // Found the meshset
01188                     ss_handle = *iter;
01189                     break;
01190                 }
01191             }
01192 
01193             // If we didn't find a sideset already
01194             if( ss_handle == 0 )
01195             {
01196                 if( mdbImpl->create_meshset( MESHSET_ORDERED | MESHSET_TRACK_OWNER, ss_handle ) != MB_SUCCESS )
01197                     return MB_FAILURE;
01198 
01199                 if( ss_handle == 0 ) return MB_FAILURE;
01200 
01201                 int sideset_id = id_array[i];
01202                 if( mdbImpl->tag_set_data( mNeumannSetTag, &ss_handle, 1, &sideset_id ) != MB_SUCCESS )
01203                     return MB_FAILURE;
01204                 if( mdbImpl->tag_set_data( mGlobalIdTag, &ss_handle, 1, &sideset_id ) != MB_SUCCESS ) return MB_FAILURE;
01205 
01206                 if( !reverse_entities.empty() )
01207                 {
01208                     // Also make a reverse set to put in this set
01209                     EntityHandle reverse_set;
01210                     if( mdbImpl->create_meshset( MESHSET_SET | MESHSET_TRACK_OWNER, reverse_set ) != MB_SUCCESS )
01211                         return MB_FAILURE;
01212 
01213                     // Add the reverse set to the sideset set and the entities to the reverse set
01214                     ErrorCode result = mdbImpl->add_entities( ss_handle, &reverse_set, 1 );
01215                     if( MB_SUCCESS != result ) return result;
01216 
01217                     result = mdbImpl->add_entities( reverse_set, &reverse_entities[0], reverse_entities.size() );
01218                     if( MB_SUCCESS != result ) return result;
01219 
01220                     // Set the reverse tag
01221                     Tag sense_tag;
01222                     int dum_sense = 0;
01223                     result        = mdbImpl->tag_get_handle( "NEUSET_SENSE", 1, MB_TYPE_INTEGER, sense_tag,
01224                                                       MB_TAG_SPARSE | MB_TAG_CREAT, &dum_sense );
01225                     if( result != MB_SUCCESS ) return result;
01226                     dum_sense = -1;
01227                     result    = mdbImpl->tag_set_data( sense_tag, &reverse_set, 1, &dum_sense );
01228                     if( result != MB_SUCCESS ) return result;
01229                 }
01230             }
01231 
01232             if( mdbImpl->add_entities( ss_handle, ( entities_to_add.empty() ) ? NULL : &entities_to_add[0],
01233                                        entities_to_add.size() ) != MB_SUCCESS )
01234                 return MB_FAILURE;
01235 
01236             // Distribution factor stuff
01237             if( number_dist_factors_in_set )
01238             {
01239                 // If this sideset does not already has a distribution factor array...set one
01240                 const void* ptr = 0;
01241                 int size        = 0;
01242                 if( MB_SUCCESS == mdbImpl->tag_get_by_ptr( mDistFactorTag, &ss_handle, 1, &ptr, &size ) )
01243                 {
01244                     const double* data = reinterpret_cast< const double* >( ptr );
01245                     std::copy( data, data + size, std::back_inserter( temp_dist_factor_vector ) );
01246                 }
01247 
01248                 ptr  = &temp_dist_factor_vector[0];
01249                 size = temp_dist_factor_vector.size();
01250                 if( mdbImpl->tag_set_by_ptr( mDistFactorTag, &ss_handle, 1, &ptr, &size ) != MB_SUCCESS )
01251                     return MB_FAILURE;
01252             }
01253         }
01254     }
01255 
01256     return MB_SUCCESS;
01257 }
01258 
01259 ErrorCode ReadNCDF::create_ss_elements( int* element_ids, int* side_list, int num_sides, int num_dist_factors,
01260                                         std::vector< EntityHandle >& entities_to_add,
01261                                         std::vector< EntityHandle >& reverse_entities,
01262                                         std::vector< double >& dist_factor_vector, int ss_seq_id )
01263 {
01264     // Determine entity type from element id
01265 
01266     // If there are dist. factors, create a vector to hold the array
01267     // and place this array as a tag onto the sideset meshset
01268 
01269     std::vector< double > temp_dist_factor_vector( num_dist_factors );
01270     std::vector< char > temp_string_storage( max_str_length + 1 );
01271     char* temp_string = &temp_string_storage[0];
01272     int temp_var;
01273     if( num_dist_factors )
01274     {
01275         INS_ID( temp_string, "dist_fact_ss%d", ss_seq_id );
01276         temp_var = -1;
01277         GET_1D_DBL_VAR( temp_string, temp_var, temp_dist_factor_vector );
01278         if( -1 == temp_var ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting dist fact variable" ); }
01279     }
01280 
01281     EntityType subtype;
01282     int num_side_nodes, num_elem_nodes;
01283     const EntityHandle* nodes;
01284     std::vector< EntityHandle > connectivity;
01285     int side_node_idx[32];
01286 
01287     int df_index = 0;
01288     int sense    = 0;
01289     for( int i = 0; i < num_sides; i++ )
01290     {
01291         ExoIIElementType exoii_type;
01292         ReadBlockData block_data;
01293         block_data.elemType = EXOII_MAX_ELEM_TYPE;
01294 
01295         if( find_side_element_type( element_ids[i], exoii_type, block_data, df_index, side_list[i] ) != MB_SUCCESS )
01296             continue;  // Isn't being read in this time
01297 
01298         EntityType type = ExoIIUtil::ExoIIElementMBEntity[exoii_type];
01299 
01300         EntityHandle ent_handle = element_ids[i] - block_data.startExoId + block_data.startMBId;
01301 
01302         const int side_num = side_list[i] - 1;
01303         if( type == MBHEX )
01304         {
01305             // Get the nodes of the element
01306             if( mdbImpl->get_connectivity( ent_handle, nodes, num_elem_nodes ) != MB_SUCCESS ) return MB_FAILURE;
01307 
01308             CN::SubEntityNodeIndices( type, num_elem_nodes, 2, side_num, subtype, num_side_nodes, side_node_idx );
01309             if( num_side_nodes <= 0 ) return MB_FAILURE;
01310 
01311             connectivity.resize( num_side_nodes );
01312             for( int k = 0; k < num_side_nodes; ++k )
01313                 connectivity[k] = nodes[side_node_idx[k]];
01314 
01315             if( MB_SUCCESS != create_sideset_element( connectivity, subtype, ent_handle, sense ) ) return MB_FAILURE;
01316             if( 1 == sense )
01317                 entities_to_add.push_back( ent_handle );
01318             else
01319                 reverse_entities.push_back( ent_handle );
01320 
01321             // Read in distribution factor array
01322             if( num_dist_factors )
01323             {
01324                 for( int k = 0; k < 4; k++ )
01325                     dist_factor_vector.push_back( temp_dist_factor_vector[df_index++] );
01326             }
01327         }
01328 
01329         // If it is a Tet
01330         else if( type == MBTET )
01331         {
01332             // Get the nodes of the element
01333             if( mdbImpl->get_connectivity( ent_handle, nodes, num_elem_nodes ) != MB_SUCCESS ) return MB_FAILURE;
01334 
01335             CN::SubEntityNodeIndices( type, num_elem_nodes, 2, side_num, subtype, num_side_nodes, side_node_idx );
01336             if( num_side_nodes <= 0 ) return MB_FAILURE;
01337 
01338             connectivity.resize( num_side_nodes );
01339             for( int k = 0; k < num_side_nodes; ++k )
01340                 connectivity[k] = nodes[side_node_idx[k]];
01341 
01342             if( MB_SUCCESS != create_sideset_element( connectivity, subtype, ent_handle, sense ) ) return MB_FAILURE;
01343             if( 1 == sense )
01344                 entities_to_add.push_back( ent_handle );
01345             else
01346                 reverse_entities.push_back( ent_handle );
01347 
01348             // Read in distribution factor array
01349             if( num_dist_factors )
01350             {
01351                 for( int k = 0; k < 3; k++ )
01352                     dist_factor_vector.push_back( temp_dist_factor_vector[df_index++] );
01353             }
01354         }
01355         else if( type == MBQUAD && exoii_type >= EXOII_SHELL && exoii_type <= EXOII_SHELL9 )
01356         {
01357             // ent_handle = CREATE_HANDLE(MBQUAD, base_id, error);
01358 
01359             // Just use this quad
01360             if( side_list[i] == 1 )
01361             {
01362                 if( 1 == sense )
01363                     entities_to_add.push_back( ent_handle );
01364                 else
01365                     reverse_entities.push_back( ent_handle );
01366 
01367                 if( num_dist_factors )
01368                 {
01369                     for( int k = 0; k < 4; k++ )
01370                         dist_factor_vector.push_back( temp_dist_factor_vector[df_index++] );
01371                 }
01372 
01373                 continue;
01374             }
01375             else if( side_list[i] == 2 )
01376             {
01377                 reverse_entities.push_back( ent_handle );
01378 
01379                 if( num_dist_factors )
01380                 {
01381                     for( int k = 0; k < 4; k++ )
01382                         dist_factor_vector.push_back( temp_dist_factor_vector[df_index++] );
01383                 }
01384 
01385                 continue;
01386             }
01387             else
01388             {
01389                 // Get the nodes of the element
01390                 if( mdbImpl->get_connectivity( ent_handle, nodes, num_elem_nodes ) != MB_SUCCESS ) return MB_FAILURE;
01391 
01392                 CN::SubEntityNodeIndices( type, num_elem_nodes, 1, side_num - 2, subtype, num_side_nodes,
01393                                           side_node_idx );
01394                 if( num_side_nodes <= 0 ) return MB_FAILURE;
01395 
01396                 connectivity.resize( num_side_nodes );
01397                 for( int k = 0; k < num_side_nodes; ++k )
01398                     connectivity[k] = nodes[side_node_idx[k]];
01399 
01400                 if( MB_SUCCESS != create_sideset_element( connectivity, subtype, ent_handle, sense ) )
01401                     return MB_FAILURE;
01402                 if( 1 == sense )
01403                     entities_to_add.push_back( ent_handle );
01404                 else
01405                     reverse_entities.push_back( ent_handle );
01406 
01407                 if( num_dist_factors )
01408                 {
01409                     for( int k = 0; k < 2; k++ )
01410                         dist_factor_vector.push_back( temp_dist_factor_vector[df_index++] );
01411                 }
01412             }
01413         }
01414         // If it is a Quad
01415         else if( type == MBQUAD )
01416         {
01417             // Get the nodes of the element
01418             if( mdbImpl->get_connectivity( ent_handle, nodes, num_elem_nodes ) != MB_SUCCESS ) return MB_FAILURE;
01419 
01420             CN::SubEntityNodeIndices( type, num_elem_nodes, 1, side_num, subtype, num_side_nodes, side_node_idx );
01421             if( num_side_nodes <= 0 ) return MB_FAILURE;
01422 
01423             connectivity.resize( num_side_nodes );
01424             for( int k = 0; k < num_side_nodes; ++k )
01425                 connectivity[k] = nodes[side_node_idx[k]];
01426 
01427             if( MB_SUCCESS != create_sideset_element( connectivity, subtype, ent_handle, sense ) ) return MB_FAILURE;
01428             if( 1 == sense )
01429                 entities_to_add.push_back( ent_handle );
01430             else
01431                 reverse_entities.push_back( ent_handle );
01432 
01433             // Read in distribution factor array
01434             if( num_dist_factors )
01435             {
01436                 for( int k = 0; k < 2; k++ )
01437                     dist_factor_vector.push_back( temp_dist_factor_vector[df_index++] );
01438             }
01439         }
01440         else if( type == MBTRI )
01441         {
01442             int side_offset = 0;
01443             if( number_dimensions() == 3 && side_list[i] <= 2 )
01444             {
01445                 if( 1 == sense )
01446                     entities_to_add.push_back( ent_handle );
01447                 else
01448                     reverse_entities.push_back( ent_handle );
01449                 if( num_dist_factors )
01450                 {
01451                     for( int k = 0; k < 3; k++ )
01452                         dist_factor_vector.push_back( temp_dist_factor_vector[df_index++] );
01453                 }
01454             }
01455             else
01456             {
01457                 if( number_dimensions() == 3 )
01458                 {
01459                     if( side_list[i] > 2 ) side_offset = 2;
01460                 }
01461 
01462                 // Get the nodes of the element
01463                 if( mdbImpl->get_connectivity( ent_handle, nodes, num_elem_nodes ) != MB_SUCCESS ) return MB_FAILURE;
01464 
01465                 CN::SubEntityNodeIndices( type, num_elem_nodes, 1, side_num - side_offset, subtype, num_side_nodes,
01466                                           side_node_idx );
01467                 if( num_side_nodes <= 0 ) return MB_FAILURE;
01468 
01469                 connectivity.resize( num_side_nodes );
01470                 for( int k = 0; k < num_side_nodes; ++k )
01471                     connectivity[k] = nodes[side_node_idx[k]];
01472 
01473                 if( MB_SUCCESS != create_sideset_element( connectivity, subtype, ent_handle, sense ) )
01474                     return MB_FAILURE;
01475                 if( 1 == sense )
01476                     entities_to_add.push_back( ent_handle );
01477                 else
01478                     reverse_entities.push_back( ent_handle );
01479 
01480                 if( num_dist_factors )
01481                 {
01482                     for( int k = 0; k < 2; k++ )
01483                         dist_factor_vector.push_back( temp_dist_factor_vector[df_index++] );
01484                 }
01485             }
01486         }
01487     }
01488 
01489     return MB_SUCCESS;
01490 }
01491 
01492 ErrorCode ReadNCDF::create_sideset_element( const std::vector< EntityHandle >& connectivity, EntityType type,
01493                                             EntityHandle& handle, int& sense )
01494 {
01495     // Get adjacent entities
01496     ErrorCode error = MB_SUCCESS;
01497     int to_dim      = CN::Dimension( type );
01498     // for matching purposes , consider only corner nodes
01499     // basically, assume everything is matching if the corner nodes are matching, and
01500     // the number of total nodes is the same
01501     int nb_corner_nodes = CN::VerticesPerEntity( type );
01502     std::vector< EntityHandle > adj_ent;
01503     mdbImpl->get_adjacencies( &( connectivity[0] ), 1, to_dim, false, adj_ent );
01504 
01505     // For each entity, see if we can find a match
01506     // If we find a match, return it
01507     bool match_found = false;
01508     std::vector< EntityHandle > match_conn;
01509     // By default, assume positive sense
01510     sense = 1;
01511     for( unsigned int i = 0; i < adj_ent.size() && match_found == false; i++ )
01512     {
01513         // Get the connectivity
01514         error = mdbImpl->get_connectivity( &( adj_ent[i] ), 1, match_conn );
01515         if( error != MB_SUCCESS ) continue;
01516 
01517         // Make sure they have the same number of vertices (higher order elements ?)
01518         if( match_conn.size() != connectivity.size() ) continue;
01519 
01520         // Find a matching node
01521         std::vector< EntityHandle >::iterator iter;
01522         std::vector< EntityHandle >::iterator end_corner = match_conn.begin() + nb_corner_nodes;
01523         iter                                             = std::find( match_conn.begin(), end_corner, connectivity[0] );
01524         if( iter == end_corner ) continue;
01525 
01526         // Rotate to match connectivity
01527         // rotate only corner nodes, ignore the rest from now on
01528         std::rotate( match_conn.begin(), iter, end_corner );
01529 
01530         bool they_match = true;
01531         for( int j = 1; j < nb_corner_nodes; j++ )
01532         {
01533             if( connectivity[j] != match_conn[j] )
01534             {
01535                 they_match = false;
01536                 break;
01537             }
01538         }
01539 
01540         // If we didn't get a match
01541         if( !they_match )
01542         {
01543             // Try the opposite sense
01544             they_match = true;
01545 
01546             int k = nb_corner_nodes - 1;
01547             for( int j = 1; j < nb_corner_nodes; )
01548             {
01549                 if( connectivity[j] != match_conn[k] )
01550                 {
01551                     they_match = false;
01552                     break;
01553                 }
01554                 ++j;
01555                 --k;
01556             }
01557             // If they matched here, sense is reversed
01558             if( they_match ) sense = -1;
01559         }
01560         match_found = they_match;
01561         if( match_found ) handle = adj_ent[i];
01562     }
01563 
01564     // If we didn't find a match, create an element
01565     if( !match_found ) error = mdbImpl->create_element( type, &connectivity[0], connectivity.size(), handle );
01566 
01567     return error;
01568 }
01569 
01570 ErrorCode ReadNCDF::find_side_element_type( const int exodus_id, ExoIIElementType& elem_type, ReadBlockData& block_data,
01571                                             int& df_index, int side_id )
01572 {
01573     std::vector< ReadBlockData >::iterator iter, end_iter;
01574     iter      = blocksLoading.begin();
01575     end_iter  = blocksLoading.end();
01576     elem_type = EXOII_MAX_ELEM_TYPE;
01577 
01578     for( ; iter != end_iter; ++iter )
01579     {
01580         if( exodus_id >= ( *iter ).startExoId && exodus_id < ( *iter ).startExoId + ( *iter ).numElements )
01581         {
01582             elem_type = ( *iter ).elemType;
01583 
01584             // If we're not reading this block in
01585             if( !( *iter ).reading_in )
01586             {
01587                 // Offset df_indes according to type
01588                 if( elem_type >= EXOII_HEX && elem_type <= EXOII_HEX27 )
01589                     df_index += 4;
01590                 else if( elem_type >= EXOII_TETRA && elem_type <= EXOII_TETRA14 )
01591                     df_index += 3;
01592                 else if( elem_type >= EXOII_QUAD && elem_type <= EXOII_QUAD9 )
01593                     df_index += 2;
01594                 else if( elem_type >= EXOII_SHELL && elem_type <= EXOII_SHELL9 )
01595                 {
01596                     if( side_id == 1 || side_id == 2 )
01597                         df_index += 4;
01598                     else
01599                         df_index += 2;
01600                 }
01601                 else if( elem_type >= EXOII_TRI && elem_type <= EXOII_TRI7 )
01602                     df_index += 3;
01603 
01604                 return MB_FAILURE;
01605             }
01606 
01607             block_data = *iter;
01608 
01609             return MB_SUCCESS;
01610         }
01611     }
01612     return MB_FAILURE;
01613 }
01614 
01615 ErrorCode ReadNCDF::read_qa_records( EntityHandle file_set )
01616 {
01617     std::vector< std::string > qa_records;
01618     read_qa_information( qa_records );
01619 
01620     std::vector< char > tag_data;
01621     for( std::vector< std::string >::iterator i = qa_records.begin(); i != qa_records.end(); ++i )
01622     {
01623         std::copy( i->begin(), i->end(), std::back_inserter( tag_data ) );
01624         tag_data.push_back( '\0' );
01625     }
01626 
01627     // If there were qa_records...tag them to the mCurrentMeshHandle
01628     if( !tag_data.empty() )
01629     {
01630         const void* ptr = &tag_data[0];
01631         int size        = tag_data.size();
01632         if( mdbImpl->tag_set_by_ptr( mQaRecordTag, &file_set, 1, &ptr, &size ) != MB_SUCCESS ) { return MB_FAILURE; }
01633     }
01634 
01635     return MB_SUCCESS;
01636 }
01637 
01638 ErrorCode ReadNCDF::read_qa_information( std::vector< std::string >& qa_record_list )
01639 {
01640     // Inquire on the genesis file to find the number of qa records
01641 
01642     int number_records = 0;
01643 
01644     int temp_dim;
01645     GET_DIM( temp_dim, "num_qa_rec", number_records );
01646     std::vector< char > data( max_str_length + 1 );
01647 
01648     for( int i = 0; i < number_records; i++ )
01649     {
01650         for( int j = 0; j < 4; j++ )
01651         {
01652             data[max_str_length] = '\0';
01653             if( read_qa_string( &data[0], i, j ) != MB_SUCCESS ) return MB_FAILURE;
01654 
01655             qa_record_list.push_back( &data[0] );
01656         }
01657     }
01658 
01659     return MB_SUCCESS;
01660 }
01661 
01662 ErrorCode ReadNCDF::read_qa_string( char* temp_string, int record_number, int record_position )
01663 {
01664     std::vector< int > dims;
01665     int temp_var = -1;
01666     GET_VAR( "qa_records", temp_var, dims );
01667     if( -1 == temp_var )
01668     {
01669         MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting qa record variable" );
01670         return MB_FAILURE;
01671     }
01672     size_t count[3], start[3];
01673     start[0] = record_number;
01674     start[1] = record_position;
01675     start[2] = 0;
01676     count[0] = 1;
01677     count[1] = 1;
01678     count[2] = max_str_length;
01679     int fail = nc_get_vara_text( ncFile, temp_var, start, count, temp_string );
01680     if( NC_NOERR != fail )
01681     {
01682         MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem setting current record number variable position" );
01683     }
01684     // Get the variable id in the exodus file
01685 
01686     return MB_SUCCESS;
01687 }
01688 
01689 // The cub_file_set contains the mesh to be updated. There could exist other
01690 // file sets that should be kept separate, such as the geometry file set from
01691 // ReadCGM.
01692 ErrorCode ReadNCDF::update( const char* exodus_file_name, const FileOptions& opts, const int num_blocks,
01693                             const int* blocks_to_load, const EntityHandle cub_file_set )
01694 {
01695     // Function : updating current database from new exodus_file.
01696     // Creator:   Jane Hu
01697     // opts is currently designed as following
01698     // tdata = <var_name>[, time][,op][,destination]
01699     // where var_name show the tag name to be updated, this version just takes
01700     // coord.
01701     // time is the optional, and it gives time step of each of the mesh
01702     // info in exodus file. It start from 1.
01703     // op is the operation that is going to be performed on the var_name info.
01704     // currently support 'set'
01705     // destination shows where to store the updated info, currently assume it is
01706     // stored in the same database by replacing the old info if there's no input
01707     // for destination, or the destination data is given in exodus format and just
01708     // need to update the coordinates.
01709     // Assumptions:
01710     // 1. Assume the num_el_blk's in both DB and update exodus file are the same.
01711     // 2. Assume num_el_in_blk1...num_el_in_blk(num_el_blk) numbers are matching, may in
01712     // different order. example: num_el_in_blk11 = num_el_in_blk22 && num_el_in_blk12 =
01713     // num_el_in_blk21.
01714     // 3. In exodus file, get node_num_map
01715     // 4. loop through the node_num_map, use it to find the node in the cub file.
01716     // 5. Replace coord[0][n] with coordx[m] + vals_nod_var1(time_step, m) for all directions for
01717     // matching nodes.
01718     //  Test: try to get hexes
01719 
01720     // *******************************************************************
01721     // Move nodes to their deformed locations.
01722     // *******************************************************************
01723     ErrorCode rval;
01724     std::string s;
01725     rval = opts.get_str_option( "tdata", s );
01726     if( MB_SUCCESS != rval ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem reading file options" ); }
01727     std::vector< std::string > tokens;
01728     tokenize( s, tokens, "," );
01729 
01730     // 1. Check for time step to find the match time
01731     int time_step = 1;
01732     if( tokens.size() > 1 && !tokens[1].empty() )
01733     {
01734         const char* time_s = tokens[1].c_str();
01735         char* endptr;
01736         long int pval   = strtol( time_s, &endptr, 0 );
01737         std::string end = endptr;
01738         if( !end.empty() )  // Syntax error
01739             return MB_TYPE_OUT_OF_RANGE;
01740 
01741         // Check for overflow (parsing long int, returning int)
01742         time_step = pval;
01743         if( pval != (long int)time_step ) return MB_TYPE_OUT_OF_RANGE;
01744         if( time_step <= 0 ) return MB_TYPE_OUT_OF_RANGE;
01745     }
01746 
01747     // 2. Check for the operations, currently support set.
01748     const char* op;
01749     if( tokens.size() < 3 || ( tokens[2] != "set" && tokens[2] != "add" ) )
01750     {
01751         MB_SET_ERR( MB_TYPE_OUT_OF_RANGE, "ReadNCDF: invalid operation specified for update" );
01752     }
01753     op = tokens[2].c_str();
01754 
01755     // Open netcdf/exodus file
01756     ncFile   = 0;
01757     int fail = nc_open( exodus_file_name, 0, &ncFile );
01758     if( !ncFile )
01759     {
01760         MB_SET_ERR( MB_FILE_DOES_NOT_EXIST, "ReadNCDF:: problem opening Netcdf/Exodus II file " << exodus_file_name );
01761     }
01762 
01763     rval = read_exodus_header();MB_CHK_ERR( rval );
01764 
01765     // Check to make sure that the requested time step exists
01766     int ncdim = -1;
01767     int max_time_steps;
01768     GET_DIM( ncdim, "time_step", max_time_steps );
01769     if( -1 == ncdim )
01770     {
01771         std::cout << "ReadNCDF: could not get number of time steps" << std::endl;
01772         return MB_FAILURE;
01773     }
01774     std::cout << "  Maximum time step=" << max_time_steps << std::endl;
01775     if( max_time_steps < time_step )
01776     {
01777         std::cout << "ReadNCDF: time step is greater than max_time_steps" << std::endl;
01778         return MB_FAILURE;
01779     }
01780 
01781     // Get the time
01782     std::vector< double > times( max_time_steps );
01783     int nc_var = -1;
01784     GET_1D_DBL_VAR( "time_whole", nc_var, times );
01785     if( -1 == nc_var ) { std::cout << "ReadNCDF: unable to get time variable" << std::endl; }
01786     else
01787     {
01788         std::cout << "  Step " << time_step << " is at " << times[time_step - 1] << " seconds" << std::endl;
01789     }
01790 
01791     // Read in the node_num_map.
01792     std::vector< int > ptr( numberNodes_loading );
01793 
01794     int varid = -1;
01795     GET_1D_INT_VAR( "node_num_map", varid, ptr );
01796     if( -1 == varid ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting node number map data" ); }
01797 
01798     // Read in the deformations.
01799     std::vector< std::vector< double > > deformed_arrays( 3 );
01800     std::vector< std::vector< double > > orig_coords( 3 );
01801     deformed_arrays[0].reserve( numberNodes_loading );
01802     deformed_arrays[1].reserve( numberNodes_loading );
01803     deformed_arrays[2].reserve( numberNodes_loading );
01804     orig_coords[0].reserve( numberNodes_loading );
01805     orig_coords[1].reserve( numberNodes_loading );
01806     orig_coords[2].reserve( numberNodes_loading );
01807     size_t start[2] = { static_cast< size_t >( time_step - 1 ), 0 },
01808            count[2] = { 1, static_cast< size_t >( numberNodes_loading ) };
01809     std::vector< int > dims;
01810     int coordx = -1, coordy = -1, coordz = -1;
01811     GET_VAR( "vals_nod_var1", coordx, dims );
01812     GET_VAR( "vals_nod_var2", coordy, dims );
01813     if( numberDimensions_loading == 3 ) GET_VAR( "vals_nod_var3", coordz, dims );
01814     if( -1 == coordx || -1 == coordy || ( numberDimensions_loading == 3 && -1 == coordz ) )
01815     {
01816         MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting coords variable" );
01817     }
01818 
01819     fail = nc_get_vara_double( ncFile, coordx, start, count, &deformed_arrays[0][0] );
01820     if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting x deformation array" ); }
01821 
01822     fail = nc_get_vara_double( ncFile, coordy, start, count, &deformed_arrays[1][0] );
01823     if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting y deformation array" ); }
01824     if( numberDimensions_loading == 3 )
01825     {
01826         fail = nc_get_vara_double( ncFile, coordz, start, count, &deformed_arrays[2][0] );
01827         if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting z deformation array" ); }
01828     }
01829 
01830     int coord1 = -1, coord2 = -1, coord3 = -1;
01831     GET_1D_DBL_VAR( "coordx", coord1, orig_coords[0] );
01832     if( -1 == coord1 ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting x coord array" ); }
01833     GET_1D_DBL_VAR( "coordy", coord2, orig_coords[1] );
01834     if( -1 == coord2 ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting y coord array" ); }
01835     if( numberDimensions_loading == 3 )
01836     {
01837         GET_1D_DBL_VAR( "coordz", coord3, orig_coords[2] );
01838         if( -1 == coord3 ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting z coord array" ); }
01839     }
01840 
01841     // b. Deal with DB file : get node info. according to node_num_map.
01842     if( tokens[0] != "coord" && tokens[0] != "COORD" ) return MB_NOT_IMPLEMENTED;
01843 
01844     if( strcmp( op, "set" ) != 0 && strcmp( op, " set" ) != 0 ) return MB_NOT_IMPLEMENTED;
01845 
01846     // Two methods of matching nodes (id vs. proximity)
01847     const bool match_node_ids = true;
01848 
01849     // Get nodes in cubit file
01850     Range cub_verts;
01851     rval = mdbImpl->get_entities_by_type( cub_file_set, MBVERTEX, cub_verts );MB_CHK_ERR( rval );
01852     std::cout << "  cub_file_set contains " << cub_verts.size() << " nodes." << std::endl;
01853 
01854     // Some accounting
01855     std::cout << "  exodus file contains " << numberNodes_loading << " nodes." << std::endl;
01856     double max_magnitude     = 0;
01857     double average_magnitude = 0;
01858     int found                = 0;
01859     int lost                 = 0;
01860     std::map< int, EntityHandle > cub_verts_id_map;
01861     AdaptiveKDTree kdtree( mdbImpl );
01862     EntityHandle root;
01863 
01864     // Should not use cub verts unless they have been matched. Place in a map
01865     // for fast handle_by_id lookup.
01866     std::map< int, EntityHandle > matched_cub_vert_id_map;
01867 
01868     // Place cub verts in a map for searching by id
01869     if( match_node_ids )
01870     {
01871         std::vector< int > cub_ids( cub_verts.size() );
01872         rval = mdbImpl->tag_get_data( mGlobalIdTag, cub_verts, &cub_ids[0] );MB_CHK_ERR( rval );
01873         for( unsigned i = 0; i != cub_verts.size(); ++i )
01874         {
01875             cub_verts_id_map.insert( std::pair< int, EntityHandle >( cub_ids[i], cub_verts[i] ) );
01876         }
01877 
01878         // Place cub verts in a kdtree for searching by proximity
01879     }
01880     else
01881     {
01882         FileOptions tree_opts( "MAX_PER_LEAF=1;SPLITS_PER_DIR=1;CANDIDATE_PLANE_SET=0" );
01883         rval = kdtree.build_tree( cub_verts, &root, &tree_opts );MB_CHK_ERR( rval );
01884         AdaptiveKDTreeIter tree_iter;
01885         rval = kdtree.get_tree_iterator( root, tree_iter );MB_CHK_ERR( rval );
01886     }
01887 
01888     // For each exo vert, find the matching cub vert
01889     for( int i = 0; i < numberNodes_loading; ++i )
01890     {
01891         int exo_id = ptr[i];
01892         CartVect exo_coords( orig_coords[0][i], orig_coords[1][i], orig_coords[2][i] );
01893         EntityHandle cub_vert = 0;
01894         bool found_match      = false;
01895 
01896         // By id
01897         if( match_node_ids )
01898         {
01899             std::map< int, EntityHandle >::iterator i_iter;
01900             i_iter = cub_verts_id_map.find( exo_id );
01901             if( i_iter != cub_verts_id_map.end() )
01902             {
01903                 found_match = true;
01904                 cub_vert    = i_iter->second;
01905             }
01906 
01907             // By proximity
01908         }
01909         else
01910         {
01911             // The MAX_NODE_DIST is the farthest distance to  search for a node.
01912             // For the 1/12th symmetry 85 pin model, the max node dist could not be less
01913             // than 1e-1 (March 26, 2010).
01914             const double MAX_NODE_DIST = 1e-1;
01915 
01916             std::vector< EntityHandle > leaves;
01917             double min_dist = MAX_NODE_DIST;
01918             rval            = kdtree.distance_search( exo_coords.array(), MAX_NODE_DIST, leaves );MB_CHK_ERR( rval );
01919             for( std::vector< EntityHandle >::const_iterator j = leaves.begin(); j != leaves.end(); ++j )
01920             {
01921                 std::vector< EntityHandle > leaf_verts;
01922                 rval = mdbImpl->get_entities_by_type( *j, MBVERTEX, leaf_verts );MB_CHK_ERR( rval );
01923                 for( std::vector< EntityHandle >::const_iterator k = leaf_verts.begin(); k != leaf_verts.end(); ++k )
01924                 {
01925                     CartVect orig_cub_coords, difference;
01926                     rval = mdbImpl->get_coords( &( *k ), 1, orig_cub_coords.array() );MB_CHK_ERR( rval );
01927                     difference  = orig_cub_coords - exo_coords;
01928                     double dist = difference.length();
01929                     if( dist < min_dist )
01930                     {
01931                         min_dist = dist;
01932                         cub_vert = *k;
01933                     }
01934                 }
01935             }
01936             if( 0 != cub_vert ) found_match = true;
01937         }
01938 
01939         // If a match is found, update it with the deformed coords from the exo file.
01940         if( found_match )
01941         {
01942             CartVect updated_exo_coords;
01943             matched_cub_vert_id_map.insert( std::pair< int, EntityHandle >( exo_id, cub_vert ) );
01944             updated_exo_coords[0] = orig_coords[0][i] + deformed_arrays[0][i];
01945             updated_exo_coords[1] = orig_coords[1][i] + deformed_arrays[1][i];
01946 
01947             if( numberDimensions_loading == 3 ) updated_exo_coords[2] = orig_coords[2][i] + deformed_arrays[2][i];
01948             rval = mdbImpl->set_coords( &cub_vert, 1, updated_exo_coords.array() );MB_CHK_ERR( rval );
01949             ++found;
01950 
01951             double magnitude =
01952                 sqrt( deformed_arrays[0][i] * deformed_arrays[0][i] + deformed_arrays[1][i] * deformed_arrays[1][i] +
01953                       deformed_arrays[2][i] * deformed_arrays[2][i] );
01954             if( magnitude > max_magnitude ) max_magnitude = magnitude;
01955             average_magnitude += magnitude;
01956         }
01957         else
01958         {
01959             ++lost;
01960             std::cout << "cannot match exo node " << exo_id << " " << exo_coords << std::endl;
01961         }
01962     }
01963 
01964     // Exo verts that could not be matched have already been printed. Now print the
01965     // cub verts that could not be matched.
01966     if( matched_cub_vert_id_map.size() < cub_verts.size() )
01967     {
01968         Range unmatched_cub_verts = cub_verts;
01969         for( std::map< int, EntityHandle >::const_iterator i = matched_cub_vert_id_map.begin();
01970              i != matched_cub_vert_id_map.end(); ++i )
01971         {
01972             unmatched_cub_verts.erase( i->second );
01973         }
01974 
01975         for( Range::const_iterator i = unmatched_cub_verts.begin(); i != unmatched_cub_verts.end(); ++i )
01976         {
01977             int cub_id;
01978             rval = mdbImpl->tag_get_data( mGlobalIdTag, &( *i ), 1, &cub_id );MB_CHK_ERR( rval );
01979 
01980             CartVect cub_coords;
01981             rval = mdbImpl->get_coords( &( *i ), 1, cub_coords.array() );MB_CHK_ERR( rval );
01982             std::cout << "cannot match cub node " << cub_id << " " << cub_coords << std::endl;
01983         }
01984 
01985         std::cout << "  " << unmatched_cub_verts.size() << " nodes from the cub file could not be matched."
01986                   << std::endl;
01987     }
01988 
01989     // Summarize statistics
01990     std::cout << "  " << found << " nodes from the exodus file were matched in the cub_file_set ";
01991     if( match_node_ids ) { std::cout << "by id." << std::endl; }
01992     else
01993     {
01994         std::cout << "by proximity." << std::endl;
01995     }
01996 
01997     // Fail if all of the nodes could not be matched.
01998     if( 0 != lost )
01999     {
02000         std::cout << "Error:  " << lost << " nodes from the exodus file could not be matched." << std::endl;
02001         // return MB_FAILURE;
02002     }
02003     std::cout << "  maximum node displacement magnitude: " << max_magnitude << " cm" << std::endl;
02004     std::cout << "  average node displacement magnitude: " << average_magnitude / found << " cm" << std::endl;
02005 
02006     // *******************************************************************
02007     // Remove dead elements from the MOAB instance.
02008     // *******************************************************************
02009 
02010     // How many element variables are in the file?
02011     int n_elem_var;
02012     GET_DIM( ncdim, "num_elem_var", n_elem_var );
02013 
02014     // Get element variable names
02015     varid       = -1;
02016     int cstatus = nc_inq_varid( ncFile, "name_elem_var", &varid );
02017     std::vector< char > names_memory( n_elem_var * max_str_length );
02018     std::vector< char* > names( n_elem_var );
02019     for( int i = 0; i < n_elem_var; ++i )
02020         names[i] = &names_memory[i * max_str_length];
02021     if( cstatus != NC_NOERR || varid == -1 )
02022     {
02023         std::cout << "ReadNCDF: name_elem_var does not exist" << std::endl;
02024         return MB_FAILURE;
02025     }
02026     int status = nc_get_var_text( ncFile, varid, &names_memory[0] );
02027     if( NC_NOERR != status ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF: Problem getting element variable names" ); }
02028 
02029     // Is one of the element variable names "death_status"? If so, get its index
02030     // in the element variable array.
02031     int death_index;
02032     bool found_death_index = false;
02033     for( int i = 0; i < n_elem_var; ++i )
02034     {
02035         std::string temp( names[i] );
02036         std::string::size_type pos0 = temp.find( "death_status" );
02037         std::string::size_type pos1 = temp.find( "Death_Status" );
02038         std::string::size_type pos2 = temp.find( "DEATH_STATUS" );
02039         if( std::string::npos != pos0 || std::string::npos != pos1 || std::string::npos != pos2 )
02040         {
02041             found_death_index = true;
02042             death_index       = i + 1;  // NetCDF variables start with 1
02043             break;
02044         }
02045     }
02046     if( !found_death_index ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF: Problem getting index of death_status variable" ); }
02047 
02048     // The exodus header has already been read. This contains the number of element
02049     // blocks.
02050 
02051     // Dead elements are listed by block. Read the block headers to determine how
02052     // many elements are in each block.
02053     rval = read_block_headers( blocks_to_load, num_blocks );
02054     if( MB_FAILURE == rval ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF: Problem reading block headers" ); }
02055 
02056     // Dead elements from the Exodus file can be located in the cub_file_set by id
02057     // or by connectivity. Currently, finding elements by id requires careful book
02058     // keeping when constructing the model in Cubit. To avoid this, one can match
02059     // elements by connectivity instead.
02060     const bool match_elems_by_connectivity = true;
02061 
02062     // Get the element id map. The ids in the map are from the elements in the blocks.
02063     // elem_num_map(blk1 elem ids, blk2 elem ids, blk3 elem ids, ...)
02064     std::vector< int > elem_ids( numberNodes_loading );
02065     if( !match_elems_by_connectivity )
02066     {
02067         GET_1D_INT_VAR( "elem_num_map", varid, elem_ids );
02068         if( -1 == varid ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF: Problem getting element number map data" ); }
02069     }
02070 
02071     // For each block
02072     int first_elem_id_in_block = 0;
02073     int block_count            = 1;  // NetCDF variables start with 1
02074     int total_elems            = 0;
02075     int total_dead_elems       = 0;
02076     for( std::vector< ReadBlockData >::iterator i = blocksLoading.begin(); i != blocksLoading.end(); ++i )
02077     {
02078         // Get the ncdf connect variable
02079         std::string temp_string( "connect" );
02080         std::stringstream temp_ss;
02081         temp_ss << block_count;
02082         temp_string += temp_ss.str();
02083         temp_string += "\0";
02084         std::vector< int > ldims;
02085         GET_VAR( temp_string.c_str(), nc_var, ldims );
02086         if( !nc_var ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting connect variable" ); }
02087         // The element type is an attribute of the connectivity variable
02088         nc_type att_type;
02089         size_t att_len;
02090         fail = nc_inq_att( ncFile, nc_var, "elem_type", &att_type, &att_len );
02091         if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting elem type attribute" ); }
02092         std::vector< char > dum_str( att_len + 1 );
02093         fail = nc_get_att_text( ncFile, nc_var, "elem_type", &dum_str[0] );
02094         if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting elem type" ); }
02095         ExoIIElementType elem_type = ExoIIUtil::static_element_name_to_type( &dum_str[0] );
02096         ( *i ).elemType            = elem_type;
02097         const EntityType mb_type   = ExoIIUtil::ExoIIElementMBEntity[( *i ).elemType];
02098 
02099         // Get the number of nodes per element
02100         unsigned int nodes_per_element = ExoIIUtil::VerticesPerElement[( *i ).elemType];
02101 
02102         // Read the connectivity into that memory.
02103         // int exo_conn[i->numElements][nodes_per_element];
02104         int* exo_conn = new int[i->numElements * nodes_per_element];
02105         // NcBool status = temp_var->get(&exo_conn[0][0], i->numElements, nodes_per_element);
02106         size_t lstart[2] = { 0, 0 }, lcount[2];
02107         lcount[0]        = i->numElements;
02108         lcount[1]        = nodes_per_element;
02109         fail             = nc_get_vara_int( ncFile, nc_var, lstart, lcount, exo_conn );
02110         if( NC_NOERR != fail )
02111         {
02112             delete[] exo_conn;
02113             MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting connectivity" );
02114         }
02115 
02116         // Get the death_status at the correct time step.
02117         std::vector< double > death_status( i->numElements );  // It seems wrong, but it uses doubles
02118         std::string array_name( "vals_elem_var" );
02119         temp_ss.str( "" );  // stringstream won't clear by temp.clear()
02120         temp_ss << death_index;
02121         array_name += temp_ss.str();
02122         array_name += "eb";
02123         temp_ss.str( "" );  // stringstream won't clear by temp.clear()
02124         temp_ss << block_count;
02125         array_name += temp_ss.str();
02126         array_name += "\0";
02127         GET_VAR( array_name.c_str(), nc_var, ldims );
02128         if( !nc_var ) { MB_SET_ERR( MB_FAILURE, "ReadNCDF:: Problem getting death status variable" ); }
02129         lstart[0] = time_step - 1;
02130         lstart[1] = 0;
02131         lcount[0] = 1;
02132         lcount[1] = i->numElements;
02133         status    = nc_get_vara_double( ncFile, nc_var, lstart, lcount, &death_status[0] );
02134         if( NC_NOERR != status )
02135         {
02136             delete[] exo_conn;
02137             MB_SET_ERR( MB_FAILURE, "ReadNCDF: Problem setting time step for death_status" );
02138         }
02139 
02140         // Look for dead elements. If there is too many dead elements and this starts
02141         // to take too long, I should put the elems in a kd-tree for more efficient
02142         // searching. Alternatively I could get the exo connectivity and match nodes.
02143         int dead_elem_counter = 0, missing_elem_counter = 0;
02144         for( int j = 0; j < i->numElements; ++j )
02145         {
02146             if( 1 != death_status[j] )
02147             {
02148                 Range cub_elem, cub_nodes;
02149                 if( match_elems_by_connectivity )
02150                 {
02151                     // Get exodus nodes for the element
02152                     std::vector< int > elem_conn( nodes_per_element );
02153                     for( unsigned int k = 0; k < nodes_per_element; ++k )
02154                     {
02155                         // elem_conn[k] = exo_conn[j][k];
02156                         elem_conn[k] = exo_conn[j * nodes_per_element + k];
02157                     }
02158                     // Get the ids of the nodes (assume we are matching by id)
02159                     // Remember that the exodus array locations start with 1 (not 0).
02160                     std::vector< int > elem_conn_node_ids( nodes_per_element );
02161                     for( unsigned int k = 0; k < nodes_per_element; ++k )
02162                     {
02163                         elem_conn_node_ids[k] = ptr[elem_conn[k] - 1];
02164                     }
02165                     // Get the cub_file_set nodes by id
02166                     // The map is a log search and takes almost no time.
02167                     // MOAB's linear tag search takes 5-10 minutes.
02168                     for( unsigned int k = 0; k < nodes_per_element; ++k )
02169                     {
02170                         std::map< int, EntityHandle >::iterator k_iter;
02171                         k_iter = matched_cub_vert_id_map.find( elem_conn_node_ids[k] );
02172 
02173                         if( k_iter == matched_cub_vert_id_map.end() )
02174                         {
02175                             std::cout << "ReadNCDF: Found no cub node with id=" << elem_conn_node_ids[k]
02176                                       << ", but expected to find only 1." << std::endl;
02177                             break;
02178                         }
02179                         cub_nodes.insert( k_iter->second );
02180                     }
02181 
02182                     if( nodes_per_element != cub_nodes.size() )
02183                     {
02184                         std::cout << "ReadNCDF: nodes_per_elemenet != cub_nodes.size()" << std::endl;
02185                         delete[] exo_conn;
02186                         return MB_INVALID_SIZE;
02187                     }
02188 
02189                     // Get the cub_file_set element with the same nodes
02190                     int to_dim = CN::Dimension( mb_type );
02191                     rval       = mdbImpl->get_adjacencies( cub_nodes, to_dim, false, cub_elem );
02192                     if( MB_SUCCESS != rval )
02193                     {
02194                         std::cout << "ReadNCDF: could not get dead cub element" << std::endl;
02195                         delete[] exo_conn;
02196                         return rval;
02197                     }
02198 
02199                     // Pronto/Presto renumbers elements, so matching cub and exo elements by
02200                     // id is not possible at this time.
02201                 }
02202                 else
02203                 {
02204                     // Get dead element's id
02205                     int elem_id = elem_ids[first_elem_id_in_block + j];
02206                     void* id[]  = { &elem_id };
02207                     // Get the element by id
02208                     rval = mdbImpl->get_entities_by_type_and_tag( cub_file_set, mb_type, &mGlobalIdTag, id, 1, cub_elem,
02209                                                                   Interface::INTERSECT );
02210                     if( MB_SUCCESS != rval )
02211                     {
02212                         delete[] exo_conn;
02213                         return rval;
02214                     }
02215                 }
02216 
02217                 if( 1 == cub_elem.size() )
02218                 {
02219                     // Delete the dead element from the cub file. It will be removed from sets
02220                     // ONLY if they are tracking meshsets.
02221                     rval = mdbImpl->remove_entities( cub_file_set, cub_elem );
02222                     if( MB_SUCCESS != rval )
02223                     {
02224                         delete[] exo_conn;
02225                         return rval;
02226                     }
02227                     rval = mdbImpl->delete_entities( cub_elem );
02228                     if( MB_SUCCESS != rval )
02229                     {
02230                         delete[] exo_conn;
02231                         return rval;
02232                     }
02233                 }
02234                 else
02235                 {
02236                     std::cout << "ReadNCDF: Should have found 1 element with  type=" << mb_type
02237                               << " in cub_file_set, but instead found " << cub_elem.size() << std::endl;
02238                     rval = mdbImpl->list_entities( cub_nodes );
02239                     ++missing_elem_counter;
02240                     delete[] exo_conn;
02241                     return MB_FAILURE;
02242                 }
02243                 ++dead_elem_counter;
02244             }
02245         }
02246         // Print some statistics
02247         temp_ss.str( "" );
02248         temp_ss << i->blockId;
02249         total_dead_elems += dead_elem_counter;
02250         total_elems += i->numElements;
02251         std::cout << "  Block " << temp_ss.str() << " has " << dead_elem_counter << "/" << i->numElements
02252                   << " dead elements." << std::endl;
02253         if( 0 != missing_elem_counter )
02254         {
02255             std::cout << "    " << missing_elem_counter
02256                       << " dead elements in this block were not found in the cub_file_set." << std::endl;
02257         }
02258 
02259         // Advance the pointers into element ids and block_count. Memory cleanup.
02260         first_elem_id_in_block += i->numElements;
02261         ++block_count;
02262         delete[] exo_conn;
02263     }
02264 
02265     std::cout << " Total: " << total_dead_elems << "/" << total_elems << " dead elements." << std::endl;
02266 
02267     // Close the file
02268     fail = nc_close( ncFile );
02269     if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "Trouble closing file" ); }
02270 
02271     ncFile = 0;
02272     return MB_SUCCESS;
02273 }
02274 
02275 void ReadNCDF::tokenize( const std::string& str, std::vector< std::string >& tokens, const char* delimiters )
02276 {
02277     std::string::size_type last = str.find_first_not_of( delimiters, 0 );
02278     std::string::size_type pos  = str.find_first_of( delimiters, last );
02279     while( std::string::npos != pos && std::string::npos != last )
02280     {
02281         tokens.push_back( str.substr( last, pos - last ) );
02282         last = str.find_first_not_of( delimiters, pos );
02283         pos  = str.find_first_of( delimiters, last );
02284         if( std::string::npos == pos ) pos = str.size();
02285     }
02286 }
02287 
02288 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines