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