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