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