MOAB: Mesh Oriented datABase  (version 5.3.0)
WriteNCDF.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 #ifdef _DEBUG
00018 // turn off warnings that say they debugging identifier has been truncated
00019 // this warning comes up when using some STL containers
00020 #pragma warning( disable : 4786 )
00021 #endif
00022 #endif
00023 
00024 #include "WriteNCDF.hpp"
00025 
00026 #include "netcdf.h"
00027 #include <utility>
00028 #include <algorithm>
00029 #include <ctime>
00030 #include <string>
00031 #include <vector>
00032 #include <cstdio>
00033 #include <cstring>
00034 #include <cassert>
00035 
00036 #include "moab/Interface.hpp"
00037 #include "moab/Range.hpp"
00038 #include "moab/CN.hpp"
00039 #include "moab/FileOptions.hpp"
00040 #include "MBTagConventions.hpp"
00041 #include "Internals.hpp"
00042 #include "ExoIIUtil.hpp"
00043 #include "moab/WriteUtilIface.hpp"
00044 #include "exodus_order.h"
00045 
00046 #ifndef MOAB_HAVE_NETCDF
00047 #error Attempt to compile WriteNCDF with NetCDF support disabled
00048 #endif
00049 
00050 namespace moab
00051 {
00052 
00053 const int TIME_STR_LEN = 11;
00054 
00055 #define INS_ID( stringvar, prefix, id ) sprintf( stringvar, prefix, id )
00056 
00057 #define GET_DIM( ncdim, name, val )                                                                             \
00058     {                                                                                                           \
00059         int gdfail = nc_inq_dimid( ncFile, name, &( ncdim ) );                                                  \
00060         if( NC_NOERR == gdfail )                                                                                \
00061         {                                                                                                       \
00062             size_t tmp_val;                                                                                     \
00063             gdfail = nc_inq_dimlen( ncFile, ncdim, &tmp_val );                                                  \
00064             if( NC_NOERR != gdfail ) { MB_SET_ERR( MB_FAILURE, "WriteNCDF:: couldn't get dimension length" ); } \
00065             else                                                                                                \
00066                 ( val ) = tmp_val;                                                                              \
00067         }                                                                                                       \
00068         else                                                                                                    \
00069             ( val ) = 0;                                                                                        \
00070     }
00071 
00072 #define GET_DIMB( ncdim, name, varname, id, val ) \
00073     INS_ID( name, varname, id );                  \
00074     GET_DIM( ncdim, name, val );
00075 
00076 #define GET_VAR( name, id, dims )                                     \
00077     {                                                                 \
00078         ( id )     = -1;                                              \
00079         int gvfail = nc_inq_varid( ncFile, name, &( id ) );           \
00080         if( NC_NOERR == gvfail )                                      \
00081         {                                                             \
00082             int ndims;                                                \
00083             gvfail = nc_inq_varndims( ncFile, id, &ndims );           \
00084             if( NC_NOERR == gvfail )                                  \
00085             {                                                         \
00086                 ( dims ).resize( ndims );                             \
00087                 gvfail = nc_inq_vardimid( ncFile, id, &( dims )[0] ); \
00088             }                                                         \
00089         }                                                             \
00090     }
00091 
00092 WriterIface* WriteNCDF::factory( Interface* iface )
00093 {
00094     return new WriteNCDF( iface );
00095 }
00096 
00097 WriteNCDF::WriteNCDF( Interface* impl )
00098     : mdbImpl( impl ), ncFile( 0 ), mCurrentMeshHandle( 0 ), mGeomDimensionTag( 0 ), repeat_face_blocks( 0 )
00099 {
00100     assert( impl != NULL );
00101 
00102     impl->query_interface( mWriteIface );
00103 
00104     // Initialize in case tag_get_handle fails below
00105     //! Get and cache predefined tag handles
00106     int negone = -1;
00107     impl->tag_get_handle( MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mMaterialSetTag, MB_TAG_SPARSE | MB_TAG_CREAT,
00108                           &negone );
00109 
00110     impl->tag_get_handle( DIRICHLET_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mDirichletSetTag, MB_TAG_SPARSE | MB_TAG_CREAT,
00111                           &negone );
00112 
00113     impl->tag_get_handle( NEUMANN_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mNeumannSetTag, MB_TAG_SPARSE | MB_TAG_CREAT,
00114                           &negone );
00115 
00116     mGlobalIdTag = impl->globalId_tag();
00117 
00118     int dum_val_array[] = { -1, -1, -1, -1 };
00119     impl->tag_get_handle( HAS_MID_NODES_TAG_NAME, 4, MB_TYPE_INTEGER, mHasMidNodesTag, MB_TAG_SPARSE | MB_TAG_CREAT,
00120                           dum_val_array );
00121 
00122     impl->tag_get_handle( "distFactor", 0, MB_TYPE_DOUBLE, mDistFactorTag,
00123                           MB_TAG_SPARSE | MB_TAG_VARLEN | MB_TAG_CREAT );
00124 
00125     impl->tag_get_handle( "qaRecord", 0, MB_TYPE_OPAQUE, mQaRecordTag, MB_TAG_SPARSE | MB_TAG_VARLEN | MB_TAG_CREAT );
00126 
00127     impl->tag_get_handle( "WriteNCDF element mark", 1, MB_TYPE_BIT, mEntityMark, MB_TAG_CREAT );
00128 }
00129 
00130 WriteNCDF::~WriteNCDF()
00131 {
00132     mdbImpl->release_interface( mWriteIface );
00133 
00134     mdbImpl->tag_delete( mEntityMark );
00135 
00136     if( 0 != ncFile ) ncFile = 0;
00137 }
00138 
00139 void WriteNCDF::reset_block( std::vector< MaterialSetData >& block_info )
00140 {
00141     std::vector< MaterialSetData >::iterator iter;
00142 
00143     for( iter = block_info.begin(); iter != block_info.end(); ++iter )
00144     {
00145         iter->elements.clear();
00146     }
00147 }
00148 
00149 void WriteNCDF::time_and_date( char* time_string, char* date_string )
00150 {
00151     struct tm* local_time;
00152     time_t calendar_time;
00153 
00154     calendar_time = time( NULL );
00155     local_time    = localtime( &calendar_time );
00156 
00157     assert( NULL != time_string && NULL != date_string );
00158 
00159     strftime( time_string, TIME_STR_LEN, "%H:%M:%S", local_time );
00160     strftime( date_string, TIME_STR_LEN, "%m/%d/%Y", local_time );
00161 
00162     // Terminate with NULL character
00163     time_string[10] = (char)NULL;
00164     date_string[10] = (char)NULL;
00165 }
00166 
00167 ErrorCode WriteNCDF::write_file( const char* exodus_file_name, const bool overwrite, const FileOptions& opts,
00168                                  const EntityHandle* ent_handles, const int num_sets,
00169                                  const std::vector< std::string >& qa_records, const Tag*, int, int user_dimension )
00170 {
00171     assert( 0 != mMaterialSetTag && 0 != mNeumannSetTag && 0 != mDirichletSetTag );
00172 
00173     if( user_dimension == 0 ) mdbImpl->get_dimension( user_dimension );
00174 
00175     if( opts.get_null_option( "REPEAT_FACE_BLOCKS" ) == MB_SUCCESS ) repeat_face_blocks = 1;
00176 
00177     std::vector< EntityHandle > blocks, nodesets, sidesets, entities;
00178 
00179     // Separate into blocks, nodesets, sidesets
00180 
00181     if( num_sets == 0 )
00182     {
00183         // Default to all defined block, nodeset and sideset-type sets
00184         Range this_range;
00185         mdbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &mMaterialSetTag, NULL, 1, this_range );
00186         std::copy( this_range.begin(), this_range.end(), std::back_inserter( blocks ) );
00187         this_range.clear();
00188         mdbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &mDirichletSetTag, NULL, 1, this_range );
00189         std::copy( this_range.begin(), this_range.end(), std::back_inserter( nodesets ) );
00190         this_range.clear();
00191         mdbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &mNeumannSetTag, NULL, 1, this_range );
00192         std::copy( this_range.begin(), this_range.end(), std::back_inserter( sidesets ) );
00193 
00194         // If there is nothing to write, write everything as one block.
00195         if( blocks.empty() && nodesets.empty() && sidesets.empty() )
00196         {
00197             this_range.clear();
00198             for( int d = user_dimension; d > 0 && this_range.empty(); --d )
00199                 mdbImpl->get_entities_by_dimension( 0, d, this_range, false );
00200             if( this_range.empty() ) return MB_FILE_WRITE_ERROR;
00201 
00202             EntityHandle block_handle;
00203             int block_id = 1;
00204             mdbImpl->create_meshset( MESHSET_SET, block_handle );
00205             mdbImpl->tag_set_data( mMaterialSetTag, &block_handle, 1, &block_id );
00206             mdbImpl->add_entities( block_handle, this_range );
00207             blocks.push_back( block_handle );
00208         }
00209     }
00210     else
00211     {
00212         int dummy;
00213         for( const EntityHandle* iter = ent_handles; iter < ent_handles + num_sets; ++iter )
00214         {
00215             if( MB_SUCCESS == mdbImpl->tag_get_data( mMaterialSetTag, &( *iter ), 1, &dummy ) && -1 != dummy )
00216                 blocks.push_back( *iter );
00217             else if( MB_SUCCESS == mdbImpl->tag_get_data( mDirichletSetTag, &( *iter ), 1, &dummy ) && -1 != dummy )
00218                 nodesets.push_back( *iter );
00219             else if( MB_SUCCESS == mdbImpl->tag_get_data( mNeumannSetTag, &( *iter ), 1, &dummy ) && -1 != dummy )
00220                 sidesets.push_back( *iter );
00221         }
00222     }
00223 
00224     // If there is nothing to write just return.
00225     if( blocks.empty() && nodesets.empty() && sidesets.empty() ) return MB_FILE_WRITE_ERROR;
00226 
00227     // Try to get mesh information
00228     ExodusMeshInfo mesh_info;
00229 
00230     std::vector< MaterialSetData > block_info;
00231     std::vector< NeumannSetData > sideset_info;
00232     std::vector< DirichletSetData > nodeset_info;
00233 
00234     mesh_info.num_dim = user_dimension;
00235 
00236     if( qa_records.empty() )
00237     {
00238         // qa records are empty - initialize some MB-standard ones
00239         mesh_info.qaRecords.push_back( "MB" );
00240         mesh_info.qaRecords.push_back( "0.99" );
00241         char string1[80], string2[80];
00242         time_and_date( string2, string1 );
00243         mesh_info.qaRecords.push_back( string2 );
00244         mesh_info.qaRecords.push_back( string1 );
00245     }
00246     else
00247     {
00248         // Constrained to multiples of 4 qa records
00249         assert( qa_records.size() % 4 == 0 );
00250 
00251         std::copy( qa_records.begin(), qa_records.end(), std::back_inserter( mesh_info.qaRecords ) );
00252     }
00253 
00254     block_info.clear();
00255     if( gather_mesh_information( mesh_info, block_info, sideset_info, nodeset_info, blocks, sidesets, nodesets ) !=
00256         MB_SUCCESS )
00257     {
00258         reset_block( block_info );
00259         return MB_FAILURE;
00260     }
00261 
00262     // Try to open the file after gather mesh info succeeds
00263     int fail = nc_create( exodus_file_name, overwrite ? NC_CLOBBER : NC_NOCLOBBER, &ncFile );
00264     if( NC_NOERR != fail )
00265     {
00266         reset_block( block_info );
00267         return MB_FAILURE;
00268     }
00269 
00270     if( write_header( mesh_info, block_info, sideset_info, nodeset_info, exodus_file_name ) != MB_SUCCESS )
00271     {
00272         reset_block( block_info );
00273         return MB_FAILURE;
00274     }
00275 
00276     {
00277         // write dummy time_whole
00278         double timev = 0.0;  // dummy, to make paraview happy
00279         size_t start = 0, count = 1;
00280         int nc_var;
00281         std::vector< int > dims;
00282         GET_VAR( "time_whole", nc_var, dims );
00283         fail = nc_put_vara_double( ncFile, nc_var, &start, &count, &timev );
00284         if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "Failed writing dist factor array" ); }
00285     }
00286 
00287     if( write_nodes( mesh_info.num_nodes, mesh_info.nodes, mesh_info.num_dim ) != MB_SUCCESS )
00288     {
00289         reset_block( block_info );
00290         return MB_FAILURE;
00291     }
00292 
00293     if( !mesh_info.polyhedronFaces.empty() )
00294     {
00295         if( write_poly_faces( mesh_info ) != MB_SUCCESS )
00296         {
00297             reset_block( block_info );
00298             return MB_FAILURE;
00299         }
00300     }
00301 
00302     if( write_elementblocks( mesh_info, block_info ) )
00303     {
00304         reset_block( block_info );
00305         return MB_FAILURE;
00306     }
00307 
00308     // Write the three maps
00309     if( write_global_node_order_map( mesh_info.num_nodes, mesh_info.nodes ) != MB_SUCCESS )
00310     {
00311         reset_block( block_info );
00312         return MB_FAILURE;
00313     }
00314 
00315     if( write_global_element_order_map( mesh_info.num_elements ) != MB_SUCCESS )
00316     {
00317         reset_block( block_info );
00318         return MB_FAILURE;
00319     }
00320 
00321     if( write_element_order_map( mesh_info.num_elements ) != MB_SUCCESS )
00322     {
00323         reset_block( block_info );
00324         return MB_FAILURE;
00325     }
00326 
00327     /*
00328      if (write_elementmap(mesh_info) != MB_SUCCESS)
00329        return MB_FAILURE;
00330     */
00331 
00332     if( write_BCs( sideset_info, nodeset_info ) != MB_SUCCESS )
00333     {
00334         reset_block( block_info );
00335         return MB_FAILURE;
00336     }
00337 
00338     if( write_qa_records( mesh_info.qaRecords ) != MB_SUCCESS ) return MB_FAILURE;
00339 
00340     // Copy the qa records into the argument
00341     // mesh_info.qaRecords.swap(qa_records);
00342     // Close the file
00343     fail = nc_close( ncFile );
00344     if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "Trouble closing file" ); }
00345 
00346     return MB_SUCCESS;
00347 }
00348 
00349 ErrorCode WriteNCDF::gather_mesh_information( ExodusMeshInfo& mesh_info, std::vector< MaterialSetData >& block_info,
00350                                               std::vector< NeumannSetData >& sideset_info,
00351                                               std::vector< DirichletSetData >& nodeset_info,
00352                                               std::vector< EntityHandle >& blocks,
00353                                               std::vector< EntityHandle >& sidesets,
00354                                               std::vector< EntityHandle >& nodesets )
00355 {
00356     ErrorCode rval;
00357     std::vector< EntityHandle >::iterator vector_iter, end_vector_iter;
00358 
00359     mesh_info.num_nodes            = 0;
00360     mesh_info.num_elements         = 0;
00361     mesh_info.num_elementblocks    = 0;
00362     mesh_info.num_polyhedra_blocks = 0;
00363 
00364     int id = 0;
00365 
00366     vector_iter     = blocks.begin();
00367     end_vector_iter = blocks.end();
00368 
00369     std::vector< EntityHandle > parent_meshsets;
00370 
00371     // Clean out the bits for the element mark
00372     rval = mdbImpl->tag_delete( mEntityMark );
00373     if( MB_SUCCESS != rval ) return rval;
00374     rval = mdbImpl->tag_get_handle( "WriteNCDF element mark", 1, MB_TYPE_BIT, mEntityMark, MB_TAG_CREAT );
00375     if( MB_SUCCESS != rval ) return rval;
00376 
00377     int highest_dimension_of_element_blocks = 0;
00378 
00379     for( vector_iter = blocks.begin(); vector_iter != blocks.end(); ++vector_iter )
00380     {
00381         MaterialSetData block_data;
00382 
00383         // For the purpose of qa records, get the parents of these blocks
00384         if( mdbImpl->get_parent_meshsets( *vector_iter, parent_meshsets ) != MB_SUCCESS ) return MB_FAILURE;
00385 
00386         // Get all Entity Handles in the mesh set
00387         Range dummy_range;
00388         rval = mdbImpl->get_entities_by_handle( *vector_iter, dummy_range, true );
00389         if( MB_SUCCESS != rval ) return rval;
00390 
00391         // Skip empty blocks
00392         if( dummy_range.empty() ) continue;
00393 
00394         // Get the block's id
00395         if( mdbImpl->tag_get_data( mMaterialSetTag, &( *vector_iter ), 1, &id ) != MB_SUCCESS )
00396         {
00397             MB_SET_ERR( MB_FAILURE, "Couldn't get block id from a tag for an element block" );
00398         }
00399 
00400         block_data.id                = id;
00401         block_data.number_attributes = 0;
00402 
00403         // Wait a minute, we are doing some filtering here that doesn't make sense at this level CJS
00404 
00405         // Find the dimension of the last entity in this range
00406         int this_dim = CN::Dimension( TYPE_FROM_HANDLE( dummy_range.back() ) );
00407         if( this_dim > 3 ) { MB_SET_ERR( MB_TYPE_OUT_OF_RANGE, "Block " << id << " contains entity sets" ); }
00408         block_data.elements = dummy_range.subset_by_dimension( this_dim );
00409 
00410         // End of -- wait a minute, we are doing some filtering here that doesn't make sense at this
00411         // level CJS
00412 
00413         // Get the entity type for this block, verifying that it's the same for all elements
00414         EntityType entity_type = TYPE_FROM_HANDLE( block_data.elements.front() );
00415         if( !block_data.elements.all_of_type( entity_type ) )
00416         {
00417             MB_SET_ERR( MB_FAILURE, "Entities in block " << id << " not of common type" );
00418         }
00419 
00420         int dimension = -1;
00421         if( entity_type == MBQUAD || entity_type == MBTRI )
00422             dimension = 2;  // Output shells by default
00423         else if( entity_type == MBEDGE )
00424             dimension = 1;
00425         else
00426             dimension = CN::Dimension( entity_type );
00427 
00428         if( dimension > highest_dimension_of_element_blocks ) highest_dimension_of_element_blocks = dimension;
00429 
00430         std::vector< EntityHandle > tmp_conn;
00431         rval = mdbImpl->get_connectivity( &( block_data.elements.front() ), 1, tmp_conn );
00432         if( MB_SUCCESS != rval ) return rval;
00433         block_data.element_type = ExoIIUtil::get_element_type_from_num_verts( tmp_conn.size(), entity_type, dimension );
00434 
00435         if( block_data.element_type == EXOII_MAX_ELEM_TYPE )
00436         {
00437             MB_SET_ERR( MB_FAILURE, "Element type in block " << id << " didn't get set correctly" );
00438         }
00439 
00440         if( block_data.element_type == EXOII_POLYGON )
00441         {
00442             // get all poly connectivity
00443             int numconn = 0;
00444             for( Range::iterator eit = block_data.elements.begin(); eit != block_data.elements.end(); eit++ )
00445             {
00446                 EntityHandle polg        = *eit;
00447                 int nnodes               = 0;
00448                 const EntityHandle* conn = NULL;
00449                 rval                     = mdbImpl->get_connectivity( polg, conn, nnodes );MB_CHK_ERR( rval );
00450                 numconn += nnodes;
00451             }
00452             block_data.number_nodes_per_element = numconn;
00453         }
00454         else
00455             block_data.number_nodes_per_element = ExoIIUtil::VerticesPerElement[block_data.element_type];
00456 
00457         // Number of nodes for this block
00458         block_data.number_elements = block_data.elements.size();
00459 
00460         // Total number of elements
00461         mesh_info.num_elements += block_data.number_elements;
00462 
00463         // Get the nodes for the elements
00464         rval = mWriteIface->gather_nodes_from_elements( block_data.elements, mEntityMark, mesh_info.nodes );
00465         if( MB_SUCCESS != rval ) return rval;
00466 
00467         // if polyhedra block
00468         if( EXOII_POLYHEDRON == block_data.element_type )
00469         {
00470             rval = mdbImpl->get_connectivity( block_data.elements, mesh_info.polyhedronFaces );MB_CHK_ERR( rval );
00471             mesh_info.num_polyhedra_blocks++;
00472         }
00473 
00474         if( !sidesets.empty() )
00475         {
00476             // If there are sidesets, keep track of which elements are being written out
00477             for( Range::iterator iter = block_data.elements.begin(); iter != block_data.elements.end(); ++iter )
00478             {
00479                 unsigned char bit = 0x1;
00480                 rval              = mdbImpl->tag_set_data( mEntityMark, &( *iter ), 1, &bit );
00481                 if( MB_SUCCESS != rval ) return rval;
00482             }
00483         }
00484 
00485         block_info.push_back( block_data );
00486 
00487         const void* data = NULL;
00488         int size         = 0;
00489         if( MB_SUCCESS == mdbImpl->tag_get_by_ptr( mQaRecordTag, &( *vector_iter ), 1, &data, &size ) && NULL != data )
00490         {
00491             // There are qa records on this block - copy over to mesh qa records
00492             const char* qa_rec = static_cast< const char* >( data );
00493             int start          = 0;
00494             int count          = 0;
00495             for( int i = 0; i < size; i++ )
00496             {
00497                 if( qa_rec[i] == '\0' )
00498                 {
00499                     std::string qa_string( &qa_rec[start], i - start );
00500                     mesh_info.qaRecords.push_back( qa_string );
00501                     start = i + 1;
00502                     count++;
00503                 }
00504             }
00505 
00506             // Constrained to multiples of 4 qa records
00507             if( count > 0 ) assert( count % 4 == 0 );
00508         }
00509     }
00510 
00511     mesh_info.num_elementblocks = block_info.size();
00512 
00513     for( std::vector< MaterialSetData >::iterator blit = block_info.begin(); blit != block_info.end(); blit++ )
00514     {
00515         MaterialSetData& block = *blit;
00516         if( block.element_type != EXOII_POLYHEDRON )
00517         {
00518             mesh_info.polyhedronFaces = subtract( mesh_info.polyhedronFaces, block.elements );
00519         }
00520     }
00521 
00522     // If user hasn't entered dimension, we figure it out
00523     if( mesh_info.num_dim == 0 )
00524     {
00525         // Never want 1 or zero dimensions
00526         if( highest_dimension_of_element_blocks < 2 )
00527             mesh_info.num_dim = 3;
00528         else
00529             mesh_info.num_dim = highest_dimension_of_element_blocks;
00530     }
00531 
00532     Range::iterator range_iter, end_range_iter;
00533     range_iter     = mesh_info.nodes.begin();
00534     end_range_iter = mesh_info.nodes.end();
00535 
00536     mesh_info.num_nodes = mesh_info.nodes.size();
00537 
00538     //------nodesets--------
00539 
00540     vector_iter     = nodesets.begin();
00541     end_vector_iter = nodesets.end();
00542 
00543     for( ; vector_iter != end_vector_iter; ++vector_iter )
00544     {
00545         DirichletSetData nodeset_data;
00546         nodeset_data.id           = 0;
00547         nodeset_data.number_nodes = 0;
00548 
00549         // Get the nodeset's id
00550         if( mdbImpl->tag_get_data( mDirichletSetTag, &( *vector_iter ), 1, &id ) != MB_SUCCESS )
00551         {
00552             MB_SET_ERR( MB_FAILURE, "Couldn't get id tag for nodeset " << id );
00553         }
00554 
00555         nodeset_data.id = id;
00556 
00557         std::vector< EntityHandle > node_vector;
00558         // Get the nodes of the nodeset that are in mesh_info.nodes
00559         if( mdbImpl->get_entities_by_handle( *vector_iter, node_vector, true ) != MB_SUCCESS )
00560         {
00561             MB_SET_ERR( MB_FAILURE, "Couldn't get nodes in nodeset " << id );
00562         }
00563 
00564         // Get the tag for distribution factors
00565         const double* dist_factor_vector;
00566         int dist_factor_size;
00567         const void* ptr = 0;
00568 
00569         int has_dist_factors = 0;
00570         if( mdbImpl->tag_get_by_ptr( mDistFactorTag, &( *vector_iter ), 1, &ptr, &dist_factor_size ) == MB_SUCCESS &&
00571             dist_factor_size )
00572             has_dist_factors = 1;
00573         dist_factor_size /= sizeof( double );
00574         dist_factor_vector = reinterpret_cast< const double* >( ptr );
00575         std::vector< EntityHandle >::iterator iter, end_iter;
00576         iter     = node_vector.begin();
00577         end_iter = node_vector.end();
00578 
00579         int j                     = 0;
00580         unsigned char node_marked = 0;
00581         ErrorCode result;
00582         for( ; iter != end_iter; ++iter )
00583         {
00584             if( TYPE_FROM_HANDLE( *iter ) != MBVERTEX ) continue;
00585             result = mdbImpl->tag_get_data( mEntityMark, &( *iter ), 1, &node_marked );MB_CHK_SET_ERR( result, "Couldn't get mark data" );
00586 
00587             if( 0x1 == node_marked )
00588             {
00589                 nodeset_data.nodes.push_back( *iter );
00590                 if( 0 != has_dist_factors )
00591                     nodeset_data.node_dist_factors.push_back( dist_factor_vector[j] );
00592                 else
00593                     nodeset_data.node_dist_factors.push_back( 1.0 );
00594             }
00595             j++;
00596         }
00597 
00598         nodeset_data.number_nodes = nodeset_data.nodes.size();
00599         nodeset_info.push_back( nodeset_data );
00600     }
00601 
00602     //------sidesets--------
00603     vector_iter     = sidesets.begin();
00604     end_vector_iter = sidesets.end();
00605 
00606     for( ; vector_iter != end_vector_iter; ++vector_iter )
00607     {
00608         NeumannSetData sideset_data;
00609 
00610         // Get the sideset's id
00611         if( mdbImpl->tag_get_data( mNeumannSetTag, &( *vector_iter ), 1, &id ) != MB_SUCCESS ) return MB_FAILURE;
00612 
00613         sideset_data.id              = id;
00614         sideset_data.mesh_set_handle = *vector_iter;
00615 
00616         // Get the sides in two lists, one forward the other reverse; starts with forward sense
00617         // by convention
00618         Range forward_elems, reverse_elems;
00619         if( get_sideset_elems( *vector_iter, 0, forward_elems, reverse_elems ) == MB_FAILURE ) return MB_FAILURE;
00620 
00621         ErrorCode result = get_valid_sides( forward_elems, mesh_info, 1, sideset_data );MB_CHK_SET_ERR( result, "Couldn't get valid sides data" );
00622         result = get_valid_sides( reverse_elems, mesh_info, -1, sideset_data );MB_CHK_SET_ERR( result, "Couldn't get valid sides data" );
00623 
00624         sideset_data.number_elements = sideset_data.elements.size();
00625         sideset_info.push_back( sideset_data );
00626     }
00627 
00628     return MB_SUCCESS;
00629 }
00630 
00631 ErrorCode WriteNCDF::get_valid_sides( Range& elems, ExodusMeshInfo& /*mesh_info*/, const int sense,
00632                                       NeumannSetData& sideset_data )
00633 {
00634     // This is where we see if underlying element of side set element is included in output
00635 
00636     // Get the sideset-based info for distribution factors
00637     const double* dist_factor_vector = 0;
00638     int dist_factor_size             = 0;
00639 
00640     // Initialize dist_fac_iter to get rid of compiler warning
00641     const double* dist_fac_iter = 0;
00642     const void* ptr             = 0;
00643     bool has_dist_factors       = false;
00644     if( mdbImpl->tag_get_by_ptr( mDistFactorTag, &( sideset_data.mesh_set_handle ), 1, &ptr, &dist_factor_size ) ==
00645             MB_SUCCESS &&
00646         dist_factor_size )
00647     {
00648         has_dist_factors   = true;
00649         dist_factor_vector = reinterpret_cast< const double* >( ptr );
00650         dist_fac_iter      = dist_factor_vector;
00651         dist_factor_size /= sizeof( double );
00652     }
00653 
00654     unsigned char element_marked = 0;
00655     ErrorCode result;
00656     for( Range::iterator iter = elems.begin(); iter != elems.end(); ++iter )
00657     {
00658         // Should insert here if "side" is a quad/tri on a quad/tri mesh
00659         result = mdbImpl->tag_get_data( mEntityMark, &( *iter ), 1, &element_marked );MB_CHK_SET_ERR( result, "Couldn't get mark data" );
00660 
00661         if( 0x1 == element_marked )
00662         {
00663             sideset_data.elements.push_back( *iter );
00664 
00665             // TJT TODO: the sense should really be # edges + 1or2
00666             sideset_data.side_numbers.push_back( ( sense == 1 ? 1 : 2 ) );
00667         }
00668         else
00669         {  // Then "side" is probably a quad/tri on a hex/tet mesh
00670             std::vector< EntityHandle > parents;
00671             int dimension = CN::Dimension( TYPE_FROM_HANDLE( *iter ) );
00672 
00673             // Get the adjacent parent element of "side"
00674             if( mdbImpl->get_adjacencies( &( *iter ), 1, dimension + 1, false, parents ) != MB_SUCCESS )
00675             {
00676 #if 0
00677         // This is not treated as an error, print warning messages for
00678         // debugging only
00679         fprintf(stderr, "[Warning]: Couldn't get adjacencies for sideset.\n");
00680 #endif
00681             }
00682 
00683             if( !parents.empty() )
00684             {
00685                 // Make sure the adjacent parent element will be output
00686                 for( unsigned int k = 0; k < parents.size(); k++ )
00687                 {
00688                     result = mdbImpl->tag_get_data( mEntityMark, &( parents[k] ), 1, &element_marked );MB_CHK_SET_ERR( result, "Couldn't get mark data" );
00689 
00690                     int side_no, this_sense, this_offset;
00691                     if( 0x1 == element_marked &&
00692                         mdbImpl->side_number( parents[k], *iter, side_no, this_sense, this_offset ) == MB_SUCCESS &&
00693                         this_sense == sense )
00694                     {
00695                         sideset_data.elements.push_back( parents[k] );
00696                         sideset_data.side_numbers.push_back( side_no + 1 );
00697                         break;
00698                     }
00699                 }
00700             }
00701             else
00702             {
00703 #if 0
00704         // This is not treated as an error, print warning messages for
00705         // debugging only
00706         fprintf(stderr, "[Warning]: No parent element exists for element in sideset %i\n", sideset_data.id);
00707 #endif
00708             }
00709         }
00710 
00711         if( sideset_data.elements.size() != 0 )
00712         {
00713             // Distribution factors
00714             int num_nodes = CN::VerticesPerEntity( TYPE_FROM_HANDLE( *iter ) );
00715             // put some dummy dist factors for polygons; why do we need them?
00716             if( TYPE_FROM_HANDLE( *iter ) == MBPOLYGON ) num_nodes = 1;  // dummy
00717             if( has_dist_factors )
00718             {
00719                 std::copy( dist_fac_iter, dist_fac_iter + num_nodes,
00720                            std::back_inserter( sideset_data.ss_dist_factors ) );
00721                 dist_fac_iter += num_nodes;
00722             }
00723             else
00724             {
00725                 for( int j = 0; j < num_nodes; j++ )
00726                     sideset_data.ss_dist_factors.push_back( 1.0 );
00727             }
00728         }
00729     }
00730 
00731     return MB_SUCCESS;
00732 }
00733 
00734 ErrorCode WriteNCDF::write_qa_records( std::vector< std::string >& qa_record_list )
00735 {
00736     int i = 0;
00737 
00738     for( std::vector< std::string >::iterator string_it = qa_record_list.begin(); string_it != qa_record_list.end(); )
00739     {
00740         for( int j = 0; j < 4; j++ )
00741             write_qa_string( ( *string_it++ ).c_str(), i, j );
00742         i++;
00743     }
00744 
00745     return MB_SUCCESS;
00746 }
00747 
00748 ErrorCode WriteNCDF::write_qa_string( const char* string, int record_number, int record_position )
00749 {
00750     // Get the variable id in the exodus file
00751 
00752     std::vector< int > dims;
00753     int temp_var = -1;
00754     GET_VAR( "qa_records", temp_var, dims );
00755     if( -1 == temp_var ) { MB_SET_ERR( MB_FAILURE, "WriteNCDF:: Problem getting qa record variable" ); }
00756     size_t count[3], start[3];
00757 
00758     // Write out the record
00759     start[0] = record_number;
00760     start[1] = record_position;
00761     start[2] = 0;
00762 
00763     count[0] = 1;
00764     count[1] = 1;
00765     count[2] = (long)strlen( string ) + 1;
00766     int fail = nc_put_vara_text( ncFile, temp_var, start, count, string );
00767     if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "Failed to position qa string variable" ); }
00768 
00769     return MB_SUCCESS;
00770 }
00771 
00772 ErrorCode WriteNCDF::write_nodes( int num_nodes, Range& nodes, int dimension )
00773 {
00774     // Write coordinates names
00775     int nc_var = -1;
00776     std::vector< int > dims;
00777     GET_VAR( "coor_names", nc_var, dims );
00778     if( -1 == nc_var ) { MB_SET_ERR( MB_FAILURE, "Trouble getting coordinate name variable" ); }
00779 
00780     size_t start[2] = { 0, 0 }, count[2] = { 1, ExoIIInterface::MAX_STR_LENGTH };
00781     char dum_str[ExoIIInterface::MAX_STR_LENGTH];
00782     strcpy( dum_str, "x" );
00783     int fail = nc_put_vara_text( ncFile, nc_var, start, count, dum_str );
00784     if( NC_NOERR != fail )
00785     {
00786         MB_SET_ERR( MB_FAILURE, "Trouble adding x coordinate name; netcdf message: " << nc_strerror( fail ) );
00787     }
00788 
00789     start[0] = 1;
00790     strcpy( dum_str, "y" );
00791     fail = nc_put_vara_text( ncFile, nc_var, start, count, dum_str );
00792     if( NC_NOERR != fail )
00793     {
00794         MB_SET_ERR( MB_FAILURE, "Trouble adding y coordinate name; netcdf message: " << nc_strerror( fail ) );
00795     }
00796 
00797     start[0] = 2;
00798     strcpy( dum_str, "z" );
00799     fail = nc_put_vara_text( ncFile, nc_var, start, count, dum_str );
00800     if( NC_NOERR != fail )
00801     {
00802         MB_SET_ERR( MB_FAILURE, "Trouble adding z coordinate name; netcdf message: " << nc_strerror( fail ) );
00803     }
00804 
00805     // See if should transform coordinates
00806     ErrorCode result;
00807     Tag trans_tag;
00808     result                = mdbImpl->tag_get_handle( MESH_TRANSFORM_TAG_NAME, 16, MB_TYPE_DOUBLE, trans_tag );
00809     bool transform_needed = true;
00810     if( result == MB_TAG_NOT_FOUND ) transform_needed = false;
00811 
00812     int num_coords_to_fill = transform_needed ? 3 : dimension;
00813 
00814     std::vector< double* > coord_arrays( 3 );
00815     coord_arrays[0] = new double[num_nodes];
00816     coord_arrays[1] = new double[num_nodes];
00817     coord_arrays[2] = NULL;
00818 
00819     if( num_coords_to_fill == 3 ) coord_arrays[2] = new double[num_nodes];
00820 
00821     result = mWriteIface->get_node_coords( dimension, num_nodes, nodes, mGlobalIdTag, 1, coord_arrays );
00822     if( result != MB_SUCCESS )
00823     {
00824         delete[] coord_arrays[0];
00825         delete[] coord_arrays[1];
00826         if( coord_arrays[2] ) delete[] coord_arrays[2];
00827         return result;
00828     }
00829 
00830     if( transform_needed )
00831     {
00832         double trans_matrix[16];
00833         const EntityHandle mesh = 0;
00834         result                  = mdbImpl->tag_get_data( trans_tag, &mesh, 0, trans_matrix );MB_CHK_SET_ERR( result, "Couldn't get transform data" );
00835 
00836         for( int i = 0; i < num_nodes; i++ )
00837         {
00838             double vec1[3];
00839             double vec2[3];
00840 
00841             vec2[0] = coord_arrays[0][i];
00842             vec2[1] = coord_arrays[1][i];
00843             vec2[2] = coord_arrays[2][i];
00844 
00845             for( int row = 0; row < 3; row++ )
00846             {
00847                 vec1[row] = 0.0;
00848                 for( int col = 0; col < 3; col++ )
00849                     vec1[row] += ( trans_matrix[( row * 4 ) + col] * vec2[col] );
00850             }
00851 
00852             coord_arrays[0][i] = vec1[0];
00853             coord_arrays[1][i] = vec1[1];
00854             coord_arrays[2][i] = vec1[2];
00855         }
00856     }
00857 
00858     // Write the nodes
00859     nc_var = -1;
00860     GET_VAR( "coord", nc_var, dims );
00861     if( -1 == nc_var ) { MB_SET_ERR( MB_FAILURE, "Trouble getting coordinate variable" ); }
00862     start[0] = 0;
00863     count[1] = num_nodes;
00864     fail     = nc_put_vara_double( ncFile, nc_var, start, count, &( coord_arrays[0][0] ) );
00865     if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "Trouble writing x coordinate" ); }
00866 
00867     start[0] = 1;
00868     fail     = nc_put_vara_double( ncFile, nc_var, start, count, &( coord_arrays[1][0] ) );
00869     if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "Trouble writing y coordinate" ); }
00870 
00871     start[0] = 2;
00872     fail     = nc_put_vara_double( ncFile, nc_var, start, count, &( coord_arrays[2][0] ) );
00873     if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "Trouble writing z coordinate" ); }
00874 
00875     delete[] coord_arrays[0];
00876     delete[] coord_arrays[1];
00877     if( coord_arrays[2] ) delete[] coord_arrays[2];
00878 
00879     return MB_SUCCESS;
00880 }
00881 
00882 ErrorCode WriteNCDF::write_poly_faces( ExodusMeshInfo& mesh_info )
00883 {
00884     // write all polygons that are not in another element block;
00885     // usually they are nowhere else, but be sure, write in this block only ones that are not in the
00886     // other blocks
00887     Range pfaces = mesh_info.polyhedronFaces;
00888 
00889     /*
00890      * int fbconn1(num_nod_per_fa1) ;
00891                 fbconn1:elem_type = "nsided" ;
00892         int fbepecnt1(num_fa_in_blk1) ;
00893                 fbepecnt1:entity_type1 = "NODE" ;
00894                 fbepecnt1:entity_type2 = "FACE" ;
00895      */
00896     if( pfaces.empty() ) return MB_SUCCESS;
00897     char wname[80];
00898     int nc_var = -1;
00899     std::vector< int > dims;
00900 
00901     // write one for each element block, to make paraview and visit happy
00902     int num_faces_in_block = (int)pfaces.size();
00903     for( unsigned int bl = 0; bl < mesh_info.num_polyhedra_blocks; bl++ )
00904     {
00905         INS_ID( wname, "fbconn%u", bl + 1 );  // it is the first block
00906         GET_VAR( wname, nc_var, dims );       // fbconn# variable, 1 dimensional
00907 
00908         INS_ID( wname, "num_nod_per_fa%u", bl + 1 );
00909         int ncdim, num_nod_per_face;
00910         GET_DIM( ncdim, wname, num_nod_per_face );
00911         int* connectivity = new int[num_nod_per_face];
00912         int ixcon = 0, j = 0;
00913         std::vector< int > fbepe( num_faces_in_block );  // fbepecnt1
00914         for( Range::iterator eit = pfaces.begin(); eit != pfaces.end(); eit++ )
00915         {
00916             EntityHandle polyg       = *eit;
00917             int nnodes               = 0;
00918             const EntityHandle* conn = NULL;
00919             ErrorCode rval           = mdbImpl->get_connectivity( polyg, conn, nnodes );MB_CHK_ERR( rval );
00920             for( int k = 0; k < nnodes; k++ )
00921                 connectivity[ixcon++] = conn[k];
00922             fbepe[j++] = nnodes;
00923         }
00924         size_t start[1] = { 0 }, count[1] = { 0 };
00925         count[0] = ixcon;
00926         int fail = nc_put_vara_int( ncFile, nc_var, start, count, connectivity );
00927         if( NC_NOERR != fail )
00928         {
00929             delete[] connectivity;
00930             MB_SET_ERR( MB_FAILURE, "Couldn't write fbconn variable" );
00931         }
00932 
00933         INS_ID( wname, "fbepecnt%u", bl + 1 );
00934         GET_VAR( wname, nc_var, dims );  // fbconn# variable, 1 dimensional
00935         count[0] = num_faces_in_block;
00936 
00937         fail = nc_put_vara_int( ncFile, nc_var, start, count, &fbepe[0] );
00938         if( NC_NOERR != fail )
00939         {
00940             delete[] connectivity;
00941             MB_SET_ERR( MB_FAILURE, "Couldn't write fbepecnt variable" );
00942         }
00943 
00944         int id = bl + 1;
00945         if( write_exodus_integer_variable( "fa_prop1", &id, bl, 1 ) != MB_SUCCESS )
00946         {
00947             MB_SET_ERR_CONT( "Problem writing element block id " << id );
00948         }
00949 
00950         int status = 1;
00951         if( write_exodus_integer_variable( "fa_status", &status, bl, 1 ) != MB_SUCCESS )
00952         {
00953             MB_SET_ERR( MB_FAILURE, "Problem writing face block status" );
00954         }
00955 
00956         delete[] connectivity;
00957         if( 0 == repeat_face_blocks ) break;  // do not repeat face blocks
00958     }
00959 
00960     return MB_SUCCESS;
00961 }
00962 ErrorCode WriteNCDF::write_header( ExodusMeshInfo& mesh_info, std::vector< MaterialSetData >& block_info,
00963                                    std::vector< NeumannSetData >& sideset_info,
00964                                    std::vector< DirichletSetData >& nodeset_info, const char* filename )
00965 {
00966     // Get the date and time
00967     char time[TIME_STR_LEN];
00968     char date[TIME_STR_LEN];
00969     time_and_date( time, date );
00970 
00971     std::string title_string = "MOAB";
00972     title_string.append( "(" );
00973     title_string.append( filename );
00974     title_string.append( "): " );
00975     title_string.append( date );
00976     title_string.append( ": " );
00977     title_string.append( "time " );
00978 
00979     if( title_string.length() > ExoIIInterface::MAX_LINE_LENGTH )
00980         title_string.resize( ExoIIInterface::MAX_LINE_LENGTH );
00981 
00982     // Initialize the exodus file
00983 
00984     int result = initialize_exodus_file( mesh_info, block_info, sideset_info, nodeset_info, title_string.c_str() );
00985 
00986     if( result == MB_FAILURE ) return MB_FAILURE;
00987 
00988     return MB_SUCCESS;
00989 }
00990 
00991 ErrorCode WriteNCDF::write_elementblocks( ExodusMeshInfo& mesh_info, std::vector< MaterialSetData >& block_data )
00992 {
00993     unsigned int i;
00994     int block_index = 0;  // Index into block list, may != 1 if there are inactive blocks
00995     int exodus_id   = 1;
00996 
00997     for( i = 0; i < block_data.size(); i++ )
00998     {
00999         MaterialSetData& block = block_data[i];
01000 
01001         unsigned int num_nodes_per_elem = block.number_nodes_per_element;
01002 
01003         // Write out the id for the block
01004 
01005         int id         = block.id;
01006         int num_values = 1;
01007 
01008         if( write_exodus_integer_variable( "eb_prop1", &id, block_index, num_values ) != MB_SUCCESS )
01009         {
01010             MB_SET_ERR_CONT( "Problem writing element block id " << id );
01011         }
01012 
01013         // Write out the block status
01014 
01015         int status = 1;
01016         if( 0 == block.number_elements ) { MB_SET_ERR( MB_FAILURE, "No elements in block " << id ); }
01017 
01018         if( write_exodus_integer_variable( "eb_status", &status, block_index, num_values ) != MB_SUCCESS )
01019         {
01020             MB_SET_ERR( MB_FAILURE, "Problem writing element block status" );
01021         }
01022 
01023         //
01024         // Map the connectivity to the new nodes
01025         const unsigned int num_elem = block.number_elements;
01026         unsigned int num_nodes      = num_nodes_per_elem * num_elem;
01027         if( EXOII_POLYGON == block.element_type || EXOII_POLYHEDRON == block.element_type )
01028         {
01029             num_nodes = num_nodes_per_elem;
01030         }
01031         int* connectivity = new int[num_nodes];
01032 
01033         ErrorCode result = MB_SUCCESS;
01034         if( block.element_type != EXOII_POLYHEDRON )
01035             mWriteIface->get_element_connect( num_elem, num_nodes_per_elem, mGlobalIdTag, block.elements, mGlobalIdTag,
01036                                               exodus_id, connectivity );
01037         if( result != MB_SUCCESS )
01038         {
01039             delete[] connectivity;
01040             MB_SET_ERR( result, "Couldn't get element array to write from" );
01041         }
01042 
01043         // If necessary, convert from EXODUS to CN node order
01044         const EntityType elem_type = ExoIIUtil::ExoIIElementMBEntity[block.element_type];
01045         assert( block.elements.all_of_type( elem_type ) );
01046         const int* reorder = 0;
01047         if( block.element_type != EXOII_POLYHEDRON && block.element_type != EXOII_POLYGON )
01048             reorder = exodus_elem_order_map[elem_type][block.number_nodes_per_element];
01049         if( reorder )
01050             WriteUtilIface::reorder( reorder, connectivity, block.number_elements, block.number_nodes_per_element );
01051 
01052         char wname[80];
01053         int nc_var = -1;
01054         std::vector< int > dims;
01055         if( block.element_type != EXOII_POLYHEDRON )
01056         {
01057             exodus_id += num_elem;
01058             INS_ID( wname, "connect%u", i + 1 );
01059 
01060             GET_VAR( wname, nc_var, dims );
01061             if( -1 == nc_var )
01062             {
01063                 delete[] connectivity;
01064                 MB_SET_ERR( MB_FAILURE, "Couldn't get connectivity variable" );
01065             }
01066         }
01067 
01068         if( EXOII_POLYGON == block.element_type )
01069         {
01070             size_t start[1] = { 0 }, count[1] = { num_nodes_per_elem };
01071             int fail = nc_put_vara_int( ncFile, nc_var, start, count, connectivity );
01072             if( NC_NOERR != fail )
01073             {
01074                 delete[] connectivity;
01075                 MB_SET_ERR( MB_FAILURE, "Couldn't write connectivity variable" );
01076             }
01077             // now put also number ebepecnt1
01078             INS_ID( wname, "ebepecnt%u", i + 1 );
01079             GET_VAR( wname, nc_var, dims );
01080             count[0] = block.number_elements;
01081             start[0] = 0;
01082             // reuse connectivity array, to not allocate another one
01083             int j = 0;
01084             for( Range::iterator eit = block.elements.begin(); eit != block.elements.end(); j++, eit++ )
01085             {
01086                 EntityHandle polg        = *eit;
01087                 int nnodes               = 0;
01088                 const EntityHandle* conn = NULL;
01089                 ErrorCode rval           = mdbImpl->get_connectivity( polg, conn, nnodes );MB_CHK_ERR( rval );
01090                 connectivity[j] = nnodes;
01091             }
01092             fail = nc_put_vara_int( ncFile, nc_var, start, count, connectivity );
01093             if( NC_NOERR != fail )
01094             {
01095                 delete[] connectivity;
01096                 MB_SET_ERR( MB_FAILURE, "Couldn't write ebepecnt variable" );
01097             }
01098         }
01099         else if( block.element_type != EXOII_POLYHEDRON )
01100         {
01101             size_t start[2] = { 0, 0 }, count[2] = { num_elem, num_nodes_per_elem };
01102             int fail = nc_put_vara_int( ncFile, nc_var, start, count, connectivity );
01103             if( NC_NOERR != fail )
01104             {
01105                 delete[] connectivity;
01106                 MB_SET_ERR( MB_FAILURE, "Couldn't write connectivity variable" );
01107             }
01108         }
01109         else  // if (block.element_type == EXOII_POLYHEDRON)
01110         {
01111             /* write a lot of stuff // faconn
01112             num_fa_in_blk1 = 15 ;
01113             num_nod_per_fa1 = 58 ;
01114             num_el_in_blk1 = 3 ;
01115             num_fac_per_el1 = 17 ;
01116             int fbconn1(num_nod_per_fa1) ;
01117                     fbconn1:elem_type = "nsided" ;
01118             int fbepecnt1(num_fa_in_blk1) ;
01119                     fbepecnt1:entity_type1 = "NODE" ;
01120                     fbepecnt1:entity_type2 = "FACE" ;
01121             int facconn1(num_fac_per_el1) ;
01122                     facconn1:elem_type = "NFACED" ;
01123             int ebepecnt1(num_el_in_blk1) ;
01124                     ebepecnt1:entity_type1 = "FACE" ;
01125                     ebepecnt1:entity_type2 = "ELEM" ;
01126 
01127              */
01128             Range& block_faces = mesh_info.polyhedronFaces;
01129             // ErrorCode rval = mdbImpl->get_connectivity(block.elements, block_faces);
01130             // MB_CHK_ERR(rval);
01131 
01132             // reuse now connectivity for facconn1
01133             INS_ID( wname, "facconn%u", i + 1 );
01134             GET_VAR( wname, nc_var, dims );  // fbconn# variable, 1 dimensional
01135 
01136             std::vector< int > ebepe( block.elements.size() );  // ebepecnt1
01137             int ixcon = 0, j = 0;
01138             size_t start[1] = { 0 }, count[1] = { 0 };
01139 
01140             for( Range::iterator eit = block.elements.begin(); eit != block.elements.end(); eit++ )
01141             {
01142                 EntityHandle polyh       = *eit;
01143                 int nfaces               = 0;
01144                 const EntityHandle* conn = NULL;
01145                 ErrorCode rval           = mdbImpl->get_connectivity( polyh, conn, nfaces );MB_CHK_ERR( rval );
01146                 for( int k = 0; k < nfaces; k++ )
01147                 {
01148                     int index = block_faces.index( conn[k] );
01149                     if( index == -1 ) MB_SET_ERR( MB_FAILURE, "Couldn't find face in polyhedron" );
01150                     connectivity[ixcon++] = index + 1;
01151                 }
01152                 ebepe[j++] = nfaces;
01153                 // num_faces+=nfaces;
01154             }
01155             count[0] = ixcon;  // facconn1
01156             int fail = nc_put_vara_int( ncFile, nc_var, start, count, connectivity );
01157             if( NC_NOERR != fail )
01158             {
01159                 delete[] connectivity;
01160                 MB_SET_ERR( MB_FAILURE, "Couldn't write fbconn variable" );
01161             }
01162 
01163             INS_ID( wname, "ebepecnt%u", i + 1 );
01164             GET_VAR( wname, nc_var, dims );  // ebepecnt# variable, 1 dimensional
01165             count[0] = block.elements.size();
01166 
01167             fail = nc_put_vara_int( ncFile, nc_var, start, count, &ebepe[0] );
01168             if( NC_NOERR != fail )
01169             {
01170                 delete[] connectivity;
01171                 MB_SET_ERR( MB_FAILURE, "Couldn't write fbepecnt variable" );
01172             }
01173         }
01174         block_index++;
01175         delete[] connectivity;
01176     }
01177 
01178     return MB_SUCCESS;
01179 }
01180 
01181 ErrorCode WriteNCDF::write_global_node_order_map( int num_nodes, Range& nodes )
01182 {
01183     // Note: this routine bypasses the standard exodusII interface for efficiency!
01184 
01185     // Node order map
01186     int* map = new int[num_nodes];
01187 
01188     // For now, output a dummy map!
01189 
01190     Range::iterator range_iter, end_iter;
01191     range_iter = nodes.begin();
01192     end_iter   = nodes.end();
01193 
01194     int i = 0;
01195 
01196     for( ; range_iter != end_iter; ++range_iter )
01197     {
01198         // TODO -- do we really want to cast this to an int?
01199         map[i++] = (int)ID_FROM_HANDLE( *range_iter );
01200     }
01201 
01202     // Output array and cleanup
01203 
01204     int error = write_exodus_integer_variable( "node_num_map", map, 0, num_nodes );
01205 
01206     if( map ) delete[] map;
01207 
01208     if( error < 0 ) { MB_SET_ERR( MB_FAILURE, "Failed writing global node order map" ); }
01209 
01210     return MB_SUCCESS;
01211 }
01212 
01213 ErrorCode WriteNCDF::write_global_element_order_map( int num_elements )
01214 {
01215     // Allocate map array
01216     int* map = new int[num_elements];
01217 
01218     // Many Sandia codes assume this map is unique, and CUBIT does not currently
01219     // have unique ids for all elements. Therefore, to make sure nothing crashes,
01220     // insert a dummy map...
01221 
01222     for( int i = 0; i < num_elements; i++ )
01223         map[i] = i + 1;
01224 
01225     // Output array and cleanup
01226 
01227     int error = write_exodus_integer_variable( "elem_num_map", map, 0, num_elements );
01228 
01229     if( map ) delete[] map;
01230 
01231     if( error < 0 ) { MB_SET_ERR( MB_FAILURE, "Failed writing global element order map" ); }
01232 
01233     return MB_SUCCESS;
01234 }
01235 
01236 ErrorCode WriteNCDF::write_element_order_map( int num_elements )
01237 {
01238     // Note: this routine bypasses the standard exodusII interface for efficiency!
01239 
01240     // Element order map
01241     int* map = new int[num_elements];
01242 
01243     // For now, output a dummy map!
01244 
01245     for( int i = 0; i < num_elements; i++ )
01246     {
01247         map[i] = i + 1;
01248     }
01249 
01250     // Output array and cleanup
01251 
01252     int error = write_exodus_integer_variable( "elem_map", map, 0, num_elements );
01253 
01254     if( map ) delete[] map;
01255 
01256     if( error < 0 ) { MB_SET_ERR( MB_FAILURE, "Failed writing element map" ); }
01257 
01258     return MB_SUCCESS;
01259 }
01260 
01261 ErrorCode WriteNCDF::write_exodus_integer_variable( const char* variable_name, int* variable_array, int start_position,
01262                                                     int number_values )
01263 {
01264     // Note: this routine bypasses the standard exodusII interface for efficiency!
01265 
01266     // Write directly to netcdf interface for efficiency
01267 
01268     // Get the variable id of the element map
01269     int nc_var = -1;
01270     std::vector< int > dims;
01271     GET_VAR( variable_name, nc_var, dims );
01272     if( -1 == nc_var )
01273     {
01274         MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to locate variable " << variable_name << " in file" );
01275     }
01276     // This contortion is necessary because netCDF is expecting nclongs;
01277     // fortunately it's necessary only when ints and nclongs aren't the same size
01278 
01279     size_t start[1], count[1];
01280     start[0] = start_position;
01281     count[0] = number_values;
01282 
01283     int fail = NC_NOERR;
01284     if( sizeof( int ) == sizeof( long ) ) { fail = nc_put_vara_int( ncFile, nc_var, start, count, variable_array ); }
01285     else
01286     {
01287         long* lptr = new long[number_values];
01288         for( int jj = 0; jj < number_values; jj++ )
01289             lptr[jj] = variable_array[jj];
01290         fail = nc_put_vara_long( ncFile, nc_var, start, count, lptr );
01291         delete[] lptr;
01292     }
01293 
01294     if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "Failed to store variable " << variable_name ); }
01295 
01296     return MB_SUCCESS;
01297 }
01298 
01299 ErrorCode WriteNCDF::write_BCs( std::vector< NeumannSetData >& sidesets, std::vector< DirichletSetData >& nodesets )
01300 {
01301     unsigned int i, j;
01302     int id;
01303     int ns_index = -1;
01304 
01305     for( std::vector< DirichletSetData >::iterator ns_it = nodesets.begin(); ns_it != nodesets.end(); ++ns_it )
01306     {
01307         // Get number of nodes in set
01308         int number_nodes = ( *ns_it ).number_nodes;
01309         if( 0 == number_nodes ) continue;
01310 
01311         // If we're here, we have a non-empty nodeset; increment the index
01312         ns_index++;
01313 
01314         // Get the node set id
01315         id = ( *ns_it ).id;
01316 
01317         // Build new array to old exodus ids
01318         int* exodus_id_array      = new int[number_nodes];
01319         double* dist_factor_array = new double[number_nodes];
01320 
01321         std::vector< EntityHandle >::iterator begin_iter, end_iter;
01322         std::vector< double >::iterator other_iter;
01323         begin_iter = ( *ns_it ).nodes.begin();
01324         end_iter   = ( *ns_it ).nodes.end();
01325         other_iter = ( *ns_it ).node_dist_factors.begin();
01326 
01327         j = 0;
01328         int exodus_id;
01329         ErrorCode result;
01330         // Fill up node array and dist. factor array at the same time
01331         for( ; begin_iter != end_iter; ++begin_iter )
01332         {
01333             result = mdbImpl->tag_get_data( mGlobalIdTag, &( *begin_iter ), 1, &exodus_id );MB_CHK_SET_ERR( result, "Problem getting id tag data" );
01334 
01335             exodus_id_array[j]   = exodus_id;
01336             dist_factor_array[j] = *( other_iter );
01337             ++other_iter;
01338             j++;
01339         }
01340 
01341         // Write out the id for the nodeset
01342 
01343         int num_values = 1;
01344 
01345         result = write_exodus_integer_variable( "ns_prop1", &id, ns_index, num_values );MB_CHK_SET_ERR_RET_VAL( result, "Problem writing node set id " << id, MB_FAILURE );
01346 
01347         // Write out the nodeset status
01348 
01349         int status = 1;
01350         if( !number_nodes ) status = 0;
01351 
01352         result = write_exodus_integer_variable( "ns_status", &status, ns_index, num_values );MB_CHK_SET_ERR_RET_VAL( result, "Problem writing node set status", MB_FAILURE );
01353 
01354         // Write it out
01355         char wname[80];
01356         int nc_var = -1;
01357         std::vector< int > dims;
01358         INS_ID( wname, "node_ns%d", ns_index + 1 );
01359         GET_VAR( wname, nc_var, dims );
01360         if( -1 == nc_var ) { MB_SET_ERR( MB_FAILURE, "Failed to get node_ns variable" ); }
01361 
01362         size_t start = 0, count = number_nodes;
01363         int fail = nc_put_vara_int( ncFile, nc_var, &start, &count, exodus_id_array );
01364         if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "Failed writing exodus id array" ); }
01365 
01366         // Write out nodeset distribution factors
01367         INS_ID( wname, "dist_fact_ns%d", ns_index + 1 );
01368         nc_var = -1;
01369         GET_VAR( wname, nc_var, dims );
01370         if( -1 == nc_var ) { MB_SET_ERR( MB_FAILURE, "Failed to get dist_fact variable" ); }
01371         fail = nc_put_vara_double( ncFile, nc_var, &start, &count, dist_factor_array );
01372         if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "Failed writing dist factor array" ); }
01373 
01374         delete[] dist_factor_array;
01375         delete[] exodus_id_array;
01376     }
01377 
01378     // Now do sidesets
01379     int ss_index = 0;  // Index of sideset - not the same as 'i' because
01380                        // only writing non-empty side sets
01381     for( i = 0; i < sidesets.size(); i++ )
01382     {
01383         NeumannSetData sideset_data = sidesets[i];
01384 
01385         // Get the side set id
01386         int side_set_id = sideset_data.id;
01387 
01388         // Get number of elements in set
01389         int number_elements = sideset_data.number_elements;
01390         if( 0 == number_elements ) continue;
01391 
01392         // Build new array to old exodus ids
01393         int* output_element_ids          = new int[number_elements];
01394         int* output_element_side_numbers = new int[number_elements];
01395 
01396         std::vector< EntityHandle >::iterator begin_iter, end_iter;
01397         begin_iter                             = sideset_data.elements.begin();
01398         end_iter                               = sideset_data.elements.end();
01399         std::vector< int >::iterator side_iter = sideset_data.side_numbers.begin();
01400 
01401         // Get the tag handle
01402         j = 0;
01403         int exodus_id;
01404 
01405         // For each "side"
01406         for( ; begin_iter != end_iter; ++begin_iter, ++side_iter )
01407         {
01408             ErrorCode result = mdbImpl->tag_get_data( mGlobalIdTag, &( *begin_iter ), 1, &exodus_id );MB_CHK_SET_ERR( result, "Problem getting exodus id for sideset element "
01409                                         << (long unsigned int)ID_FROM_HANDLE( *begin_iter ) );
01410 
01411             output_element_ids[j]            = exodus_id;
01412             output_element_side_numbers[j++] = *side_iter;
01413         }
01414 
01415         if( 0 != number_elements )
01416         {
01417             // Write out the id for the nodeset
01418 
01419             int num_values = 1;
01420 
01421             // ss_prop1[ss_index] = side_set_id
01422             ErrorCode result = write_exodus_integer_variable( "ss_prop1", &side_set_id, ss_index, num_values );MB_CHK_SET_ERR_RET_VAL( result, "Problem writing node set id " << id, MB_FAILURE );
01423 
01424             // FIXME : Something seems wrong here.  The we are within a block
01425             // started with if (0 != number_elements), so this condition is always
01426             // false.  This code seems to indicate that we want to write all
01427             // sidesets, not just empty ones.  But the code both here and in
01428             // initialize_exodus_file() skip empty side sets.
01429             //  - j.k. 2007-03-09
01430             int status = 1;
01431             if( 0 == number_elements ) status = 0;
01432 
01433             // ss_status[ss_index] = status
01434             result = write_exodus_integer_variable( "ss_status", &status, ss_index, num_values );MB_CHK_SET_ERR_RET_VAL( result, "Problem writing side set status", MB_FAILURE );
01435 
01436             // Increment ss_index now because we want a) we need to
01437             // increment it somewhere within the if (0 != number_elements)
01438             // block and b) the above calls need a zero-based index
01439             // while the following use a one-based index.
01440             ++ss_index;
01441 
01442             char wname[80];
01443             int nc_var;
01444             std::vector< int > dims;
01445             INS_ID( wname, "elem_ss%d", ss_index );
01446             GET_VAR( wname, nc_var, dims );
01447             if( -1 == nc_var ) { MB_SET_ERR( MB_FAILURE, "Failed to get elem_ss variable" ); }
01448             size_t start = 0, count = number_elements;
01449             int fail = nc_put_vara_int( ncFile, nc_var, &start, &count, output_element_ids );
01450             if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "Failed writing sideset element array" ); }
01451 
01452             INS_ID( wname, "side_ss%d", ss_index );
01453             nc_var = -1;
01454             GET_VAR( wname, nc_var, dims );
01455             if( -1 == nc_var ) { MB_SET_ERR( MB_FAILURE, "Failed to get side_ss variable" ); }
01456             fail = nc_put_vara_int( ncFile, nc_var, &start, &count, output_element_side_numbers );
01457             if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "Failed writing sideset side array" ); }
01458 
01459             INS_ID( wname, "dist_fact_ss%d", ss_index );
01460             nc_var = -1;
01461             GET_VAR( wname, nc_var, dims );
01462             if( -1 == nc_var ) { MB_SET_ERR( MB_FAILURE, "Failed to get sideset dist factors variable" ); }
01463             count = sideset_data.ss_dist_factors.size();
01464             fail  = nc_put_vara_double( ncFile, nc_var, &start, &count, &( sideset_data.ss_dist_factors[0] ) );
01465             if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "Failed writing sideset dist factors array" ); }
01466         }
01467 
01468         delete[] output_element_ids;
01469         delete[] output_element_side_numbers;
01470     }
01471 
01472     return MB_SUCCESS;
01473 }
01474 
01475 ErrorCode WriteNCDF::initialize_exodus_file( ExodusMeshInfo& mesh_info, std::vector< MaterialSetData >& block_data,
01476                                              std::vector< NeumannSetData >& sideset_data,
01477                                              std::vector< DirichletSetData >& nodeset_data, const char* title_string,
01478                                              bool write_maps, bool /* write_sideset_distribution_factors */ )
01479 {
01480     // This routine takes the place of the exodusii routine ex_put_init,
01481     // and additionally pre-defines variables such as qa, element blocks,
01482     // sidesets and nodesets in a single pass.
01483     //
01484     // This is necessary because of the way exodusII works.  Each time the
01485     // netcdf routine endef is called to take the file out of define mode,
01486     // the entire file is copied before the new information is added.
01487     //
01488     // With very large files, this is simply not workable.  This routine takes
01489     // the definition portions of all applicable exodus routines and puts them
01490     // in a single definition, preventing repeated copying of the file.
01491     //
01492     // Most of the code is copied directly from the applicable exodus routine,
01493     // and thus this routine may not seem very consistent in usage or variable
01494     // naming, etc.
01495 
01496     // Perform the initializations
01497 
01498     int element_block_index;
01499 
01500     // Inquire on defined string dimension and general dimension for qa
01501 
01502     int dim_str, dim_four, dim_line, dim_time;
01503     if( nc_def_dim( ncFile, "len_string", ExoIIInterface::MAX_STR_LENGTH, &dim_str ) != NC_NOERR )
01504     {
01505         MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to get string length in file" );
01506     }
01507 
01508     if( nc_def_dim( ncFile, "len_line", ExoIIInterface::MAX_STR_LENGTH, &dim_line ) != NC_NOERR )
01509     {
01510         MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to get line length in file" );
01511     }
01512 
01513     if( nc_def_dim( ncFile, "four", 4, &dim_four ) != NC_NOERR )
01514     {
01515         MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to locate four in file" );
01516     }
01517 
01518     if( nc_def_dim( ncFile, "time_step", 1, &dim_time ) != NC_NOERR )
01519     {
01520         MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to locate time step in file" );
01521     }
01522     // some whole_time dummy :(
01523     int dtime;
01524     if( NC_NOERR != nc_def_var( ncFile, "time_whole", NC_DOUBLE, 1, &dim_time, &dtime ) )
01525     {
01526         MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define time whole array" );
01527     }
01528 
01529     /* Put file into define mode */
01530 
01531     // It is possible that an NT filename using backslashes is in the title string
01532     // this can give fits to unix codes where the backslash is assumed to be an escape
01533     // sequence.  For the exodus file, just flip the backslash to a slash to prevent
01534     // this problem
01535 
01536     // Get a working copy of the title_string;
01537 
01538     char working_title[80];
01539     strncpy( working_title, title_string, 79 );
01540 
01541     int length = strlen( working_title );
01542     for( int pos = 0; pos < length; pos++ )
01543     {
01544         if( working_title[pos] == '\\' ) strncpy( &working_title[pos], "/", 1 );
01545     }
01546 
01547     if( NC_NOERR != nc_put_att_text( ncFile, NC_GLOBAL, "title", length, working_title ) )
01548     {
01549         MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define title attribute" );
01550     }
01551 
01552     // Add other attributes while we're at it
01553     float dum_vers = 6.28F;
01554     if( NC_NOERR != nc_put_att_float( ncFile, NC_GLOBAL, "api_version", NC_FLOAT, 1, &dum_vers ) )
01555     {
01556         MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define api_version attribute" );
01557     }
01558     dum_vers = 6.28F;
01559     if( NC_NOERR != nc_put_att_float( ncFile, NC_GLOBAL, "version", NC_FLOAT, 1, &dum_vers ) )
01560     {
01561         MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define version attribute" );
01562     }
01563     int dum_siz = sizeof( double );
01564     if( NC_NOERR != nc_put_att_int( ncFile, NC_GLOBAL, "floating_point_word_size", NC_INT, 1, &dum_siz ) )
01565     {
01566         MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define floating pt word size attribute" );
01567     }
01568 
01569     // Set up number of dimensions
01570 
01571     int num_el_blk, num_elem, num_nodes, num_dim, num_fa_blk, num_faces;
01572     if( nc_def_dim( ncFile, "num_dim", (size_t)mesh_info.num_dim, &num_dim ) != NC_NOERR )
01573     {
01574         MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of dimensions" );
01575     }
01576 
01577     if( nc_def_dim( ncFile, "num_nodes", mesh_info.num_nodes, &num_nodes ) != NC_NOERR )
01578     {
01579         MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of nodes" );
01580     }
01581 
01582     int num_nod_per_fa;  // it is needed for polyhedron only; need to compute it (connectivity of
01583                          // faces of polyhedra)
01584     if( mesh_info.polyhedronFaces.size() > 0 )
01585         if( nc_def_dim( ncFile, "num_faces", (int)mesh_info.polyhedronFaces.size(), &num_faces ) != NC_NOERR )
01586         {
01587             MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of nodes" );
01588         }
01589 
01590     if( nc_def_dim( ncFile, "num_elem", mesh_info.num_elements, &num_elem ) != NC_NOERR )
01591     {
01592         MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of elements" );
01593     }
01594 
01595     if( nc_def_dim( ncFile, "num_el_blk", mesh_info.num_elementblocks, &num_el_blk ) != NC_NOERR )
01596     {
01597         MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of element blocks" );
01598     }
01599 
01600     /* ...and some variables */
01601 
01602     /* Element block id status array */
01603     int idstat = -1;
01604     if( NC_NOERR != nc_def_var( ncFile, "eb_status", NC_LONG, 1, &num_el_blk, &idstat ) )
01605     {
01606         MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define element block status array" );
01607     }
01608 
01609     /* Element block id array */
01610 
01611     int idarr = -1;
01612     if( NC_NOERR != nc_def_var( ncFile, "eb_prop1", NC_LONG, 1, &num_el_blk, &idarr ) )
01613     {
01614         MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define element block id array" );
01615     }
01616 
01617     /*   store property name as attribute of property array variable */
01618     if( NC_NOERR != nc_put_att_text( ncFile, idarr, "name", strlen( "ID" ), "ID" ) )
01619     {
01620         MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store element block property name ID" );
01621     }
01622 
01623     // count how many are polyhedron blocks
01624     int num_fa_blocks = 0, num_polyh_blocks = 0;
01625     for( unsigned int i = 0; i < block_data.size(); i++ )
01626     {
01627         MaterialSetData& block = block_data[i];
01628         if( EXOII_POLYHEDRON == block.element_type )
01629         {
01630             num_fa_blocks++;
01631             num_polyh_blocks++;
01632         }
01633     }
01634     if( 0 == this->repeat_face_blocks && num_fa_blocks > 1 ) num_fa_blocks = 1;
01635     char wname[80];
01636 
01637     if( num_fa_blocks > 0 )
01638     {
01639         /* face block id status array */
01640         if( nc_def_dim( ncFile, "num_fa_blk", num_fa_blocks, &num_fa_blk ) != NC_NOERR )
01641         {
01642             MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of face blocks" );
01643         }
01644 
01645         int idstatf = -1;
01646         if( NC_NOERR != nc_def_var( ncFile, "fa_status", NC_LONG, 1, &num_fa_blk, &idstatf ) )
01647         {
01648             MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define face block status array" );
01649         }
01650 
01651         /* Element block id array */
01652 
01653         int idarrf = -1;
01654         if( NC_NOERR != nc_def_var( ncFile, "fa_prop1", NC_LONG, 1, &num_fa_blk, &idarrf ) )
01655         {
01656             MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define face block id array" );
01657         }
01658 
01659         /*   store property name as attribute of property array variable */
01660         if( NC_NOERR != nc_put_att_text( ncFile, idarrf, "name", strlen( "ID" ), "ID" ) )
01661         {
01662             MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store face block property name ID" );
01663         }
01664         // determine the number of num_nod_per_face
01665         /*
01666           num_fa_in_blk1 = 15 ;
01667           num_nod_per_fa1 = 58 ;
01668 
01669           int fbconn1(num_nod_per_fa1) ;
01670                   fbconn1:elem_type = "nsided" ;
01671           int fbepecnt1(num_fa_in_blk1) ;
01672                   fbepecnt1:entity_type1 = "NODE" ;
01673                   fbepecnt1:entity_type2 = "FACE" ;
01674            */
01675 
01676         int num_nodes_per_face = 0;
01677 
01678         int dims[1];  // maybe 1 is enough here
01679         for( Range::iterator eit = mesh_info.polyhedronFaces.begin(); eit != mesh_info.polyhedronFaces.end(); eit++ )
01680         {
01681             EntityHandle polyg       = *eit;
01682             int nnodes               = 0;
01683             const EntityHandle* conn = NULL;
01684             ErrorCode rval           = mdbImpl->get_connectivity( polyg, conn, nnodes );MB_CHK_ERR( rval );
01685             num_nodes_per_face += nnodes;
01686         }
01687 
01688         // duplicate if needed; default is not duplicate
01689         for( int j = 1; j <= num_fa_blocks; j++ )
01690         {
01691             INS_ID( wname, "num_nod_per_fa%d", j );
01692             if( nc_def_dim( ncFile, wname, (size_t)num_nodes_per_face, &num_nod_per_fa ) != NC_NOERR )
01693             {
01694                 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of nodes for face block " );
01695             }
01696             dims[0] = num_nod_per_fa;
01697             INS_ID( wname, "fbconn%d", j );  // first one, or more
01698             int fbconn;
01699             if( NC_NOERR != nc_def_var( ncFile, wname, NC_LONG, 1, dims, &fbconn ) )
01700             {
01701                 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create connectivity array for face block " << 1 );
01702             }
01703             std::string element_type_string( "nsided" );
01704             if( NC_NOERR != nc_put_att_text( ncFile, fbconn, "elem_type", element_type_string.length(),
01705                                              element_type_string.c_str() ) )
01706             {
01707                 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store element type nsided " );
01708             }
01709 
01710             INS_ID( wname, "num_fa_in_blk%d", j );  // first one, or more
01711             int num_fa_in_blk;
01712             if( nc_def_dim( ncFile, wname, (size_t)mesh_info.polyhedronFaces.size(), &num_fa_in_blk ) != NC_NOERR )
01713             {
01714                 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of nodes for face block " );
01715             }
01716 
01717             // fbepecnt
01718             INS_ID( wname, "fbepecnt%d", j );  // first one, or more
01719             int fbepecnt;
01720             dims[0] = num_fa_in_blk;
01721             if( NC_NOERR != nc_def_var( ncFile, wname, NC_LONG, 1, dims, &fbepecnt ) )
01722             {
01723                 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create fbepecnt array for block " << 1 );
01724             }
01725             std::string enttype1( "NODE" );
01726             if( NC_NOERR != nc_put_att_text( ncFile, fbepecnt, "entity_type1", enttype1.length(), enttype1.c_str() ) )
01727             {
01728                 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store entity type 1  " );
01729             }
01730             std::string enttype2( "FACE" );
01731             if( NC_NOERR != nc_put_att_text( ncFile, fbepecnt, "entity_type2", enttype2.length(), enttype2.c_str() ) )
01732             {
01733                 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store entity type 2  " );
01734             }
01735         }
01736     }
01737 
01738     // Define element blocks
01739 
01740     for( unsigned int i = 0; i < block_data.size(); i++ )
01741     {
01742         MaterialSetData& block = block_data[i];
01743 
01744         element_block_index = i + 1;
01745         int num_el_in_blk = -1, num_att_in_blk = -1;
01746         int blk_attrib, connect;
01747 
01748         /* Define number of elements in this block */
01749 
01750         INS_ID( wname, "num_el_in_blk%d", element_block_index );
01751         if( nc_def_dim( ncFile, wname, (size_t)block.number_elements, &num_el_in_blk ) != NC_NOERR )
01752         {
01753             MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of elements/block for block " << i + 1 );
01754         }
01755 
01756         /* Define number of nodes per element for this block */
01757         INS_ID( wname, "num_nod_per_el%d", element_block_index );
01758         int num_nod_per_el = -1;
01759         if( EXOII_POLYHEDRON != block.element_type )
01760             if( nc_def_dim( ncFile, wname, (size_t)block.number_nodes_per_element, &num_nod_per_el ) != NC_NOERR )
01761             {
01762                 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of nodes/element for block " << block.id );
01763             }
01764 
01765         /* Define element attribute array for this block */
01766         int dims[3];
01767         if( block.number_attributes > 0 )
01768         {
01769             INS_ID( wname, "num_att_in_blk%d", element_block_index );
01770             if( nc_def_dim( ncFile, wname, (size_t)block.number_attributes, &num_att_in_blk ) != NC_NOERR )
01771             {
01772                 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of attributes in block " << block.id );
01773             }
01774 
01775             INS_ID( wname, "attrib%d", element_block_index );
01776             dims[0] = num_el_in_blk;
01777             dims[1] = num_att_in_blk;
01778             if( NC_NOERR != nc_def_var( ncFile, wname, NC_DOUBLE, 2, dims, &blk_attrib ) )
01779             {
01780                 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define attributes for element block " << block.id );
01781             }
01782         }
01783 
01784         /* Define element connectivity array for this block */
01785 
01786         if( EXOII_POLYGON != block.element_type && EXOII_POLYHEDRON != block.element_type )
01787         {
01788             INS_ID( wname, "connect%d", element_block_index );
01789             dims[0] = num_el_in_blk;
01790             dims[1] = num_nod_per_el;
01791             if( NC_NOERR != nc_def_var( ncFile, wname, NC_LONG, 2, dims, &connect ) )
01792             {
01793                 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create connectivity array for block " << i + 1 );
01794             }
01795 
01796             /* Store element type as attribute of connectivity variable */
01797 
01798             std::string element_type_string( ExoIIUtil::ElementTypeNames[block.element_type] );
01799             if( NC_NOERR != nc_put_att_text( ncFile, connect, "elem_type", element_type_string.length(),
01800                                              element_type_string.c_str() ) )
01801             {
01802                 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store element type name " << (int)block.element_type );
01803             }
01804         }
01805         else if( EXOII_POLYGON == block.element_type )
01806         {
01807             INS_ID( wname, "connect%d", element_block_index );
01808             // need to define num_nod_per_el as total number of nodes
01809             // ebepecnt1 as number of nodes per polygon
01810             /*
01811              * int connect1(num_nod_per_el1) ;
01812                   connect1:elem_type = "nsided" ;
01813                int ebepecnt1(num_el_in_blk1) ;
01814                   ebepecnt1:entity_type1 = "NODE" ;
01815                   ebepecnt1:entity_type2 = "ELEM" ;*/
01816             dims[0] = num_nod_per_el;
01817             if( NC_NOERR != nc_def_var( ncFile, wname, NC_LONG, 1, dims, &connect ) )
01818             {
01819                 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create connectivity array for block " << i + 1 );
01820             }
01821             std::string element_type_string( "nsided" );
01822             if( NC_NOERR != nc_put_att_text( ncFile, connect, "elem_type", element_type_string.length(),
01823                                              element_type_string.c_str() ) )
01824             {
01825                 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store element type name " << (int)block.element_type );
01826             }
01827             INS_ID( wname, "ebepecnt%d", element_block_index );
01828             int ebepecnt;
01829             dims[0] = num_el_in_blk;
01830             if( NC_NOERR != nc_def_var( ncFile, wname, NC_LONG, 1, dims, &ebepecnt ) )
01831             {
01832                 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create ebepecnt array for block " << i + 1 );
01833             }
01834             std::string etype1( "NODE" );
01835             if( NC_NOERR != nc_put_att_text( ncFile, ebepecnt, "entity_type1", etype1.length(), etype1.c_str() ) )
01836             {
01837                 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store entity type1 " << (int)block.element_type );
01838             }
01839             std::string etype2( "ELEM" );
01840             if( NC_NOERR != nc_put_att_text( ncFile, ebepecnt, "entity_type2", etype2.length(), etype2.c_str() ) )
01841             {
01842                 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store entity type2 " << (int)block.element_type );
01843             }
01844         }
01845         else if( EXOII_POLYHEDRON == block.element_type )
01846         {
01847             // INS_ID(wname, "connect%d", element_block_index);
01848             /*
01849             testn  face   example: 3 polyh, 15 total faces, 2 shared; 15+2 = 17
01850             num_elem = 3 ;
01851             num_face = 15 ; // not needed?
01852             num_el_blk = 1 ;
01853 
01854             num_el_in_blk1 = 3 ;
01855             num_fac_per_el1 = 17 ;
01856 
01857              * num faces will be total face conn
01858              * num_faces_in_block will be number of faces (non repeated)
01859              * num_nodes_per_face will have total face connectivity
01860              */
01861             int num_faces2 = 0;
01862 
01863             for( Range::iterator eit = block.elements.begin(); eit != block.elements.end(); eit++ )
01864             {
01865                 EntityHandle polyh       = *eit;
01866                 int nfaces               = 0;
01867                 const EntityHandle* conn = NULL;
01868                 ErrorCode rval           = mdbImpl->get_connectivity( polyh, conn, nfaces );MB_CHK_ERR( rval );
01869                 num_faces2 += nfaces;
01870             }
01871 
01872             int num_fac_per_el;
01873 
01874             INS_ID( wname, "num_fac_per_el%d", element_block_index );
01875             if( nc_def_dim( ncFile, wname, (size_t)num_faces2, &num_fac_per_el ) != NC_NOERR )
01876             {
01877                 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of faces per block " << block.id );
01878             }
01879 
01880             /*
01881             int facconn1(num_fac_per_el1) ;
01882                     facconn1:elem_type = "NFACED" ;
01883             int ebepecnt1(num_el_in_blk1) ;
01884                     ebepecnt1:entity_type1 = "FACE" ;
01885                     ebepecnt1:entity_type2 = "ELEM" ;
01886              */
01887 
01888             // facconn
01889             INS_ID( wname, "facconn%d", element_block_index );
01890             int facconn;
01891             dims[0] = num_fac_per_el;
01892             if( NC_NOERR != nc_def_var( ncFile, wname, NC_LONG, 1, dims, &facconn ) )
01893             {
01894                 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create facconn array for block " << i + 1 );
01895             }
01896             std::string etype( "NFACED" );
01897             if( NC_NOERR != nc_put_att_text( ncFile, facconn, "elem_type", etype.length(), etype.c_str() ) )
01898             {
01899                 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store elem type " << (int)block.element_type );
01900             }
01901 
01902             // ebepecnt
01903             INS_ID( wname, "ebepecnt%d", element_block_index );
01904             int ebepecnt;
01905             dims[0] = num_el_in_blk;
01906             if( NC_NOERR != nc_def_var( ncFile, wname, NC_LONG, 1, dims, &ebepecnt ) )
01907             {
01908                 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create ebepecnt array for block " << i + 1 );
01909             }
01910             std::string etype1( "FACE" );
01911             if( NC_NOERR != nc_put_att_text( ncFile, ebepecnt, "entity_type1", etype1.length(), etype1.c_str() ) )
01912             {
01913                 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store entity type1 " << (int)block.element_type );
01914             }
01915             std::string etype2( "ELEM" );
01916             if( NC_NOERR != nc_put_att_text( ncFile, ebepecnt, "entity_type2", etype2.length(), etype2.c_str() ) )
01917             {
01918                 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store entity type2 " << (int)block.element_type );
01919             }
01920 
01921             block.number_nodes_per_element = num_faces2;  // connectivity for all polyhedra in block
01922         }
01923     }
01924 
01925     /* Node set id array: */
01926 
01927     int non_empty_nss = 0;
01928     // Need to go through nodesets to compute # nodesets, some might be empty
01929     std::vector< DirichletSetData >::iterator ns_it;
01930     for( ns_it = nodeset_data.begin(); ns_it != nodeset_data.end(); ++ns_it )
01931     {
01932         if( 0 != ( *ns_it ).number_nodes ) non_empty_nss++;
01933     }
01934 
01935     int num_ns    = -1;
01936     int ns_idstat = -1, ns_idarr = -1;
01937     if( non_empty_nss > 0 )
01938     {
01939         if( nc_def_dim( ncFile, "num_node_sets", ( size_t )( non_empty_nss ), &num_ns ) != NC_NOERR )
01940         {
01941             MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of node sets" );
01942         }
01943 
01944         /* Node set id status array: */
01945 
01946         if( NC_NOERR != nc_def_var( ncFile, "ns_status", NC_LONG, 1, &num_ns, &ns_idstat ) )
01947         {
01948             MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create node sets status array" );
01949         }
01950 
01951         /* Node set id array: */
01952         if( NC_NOERR != nc_def_var( ncFile, "ns_prop1", NC_LONG, 1, &num_ns, &ns_idarr ) )
01953         {
01954             MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create node sets property array" );
01955         }
01956 
01957         /* Store property name as attribute of property array variable */
01958         if( NC_NOERR != nc_put_att_text( ncFile, NC_GLOBAL, "name", strlen( "ID" ), "ID" ) )
01959         {
01960             MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store node set property name ID" );
01961         }
01962 
01963         // Now, define the arrays needed for each node set
01964 
01965         int index = 0;
01966 
01967         for( unsigned i = 0; i < nodeset_data.size(); i++ )
01968         {
01969             DirichletSetData node_set = nodeset_data[i];
01970 
01971             if( 0 == node_set.number_nodes )
01972             {
01973                 MB_SET_ERR_CONT( "WriteNCDF: empty nodeset " << node_set.id );
01974                 continue;
01975             }
01976             index++;
01977 
01978             int num_nod_ns = -1;
01979             INS_ID( wname, "num_nod_ns%d", index );
01980             if( nc_def_dim( ncFile, wname, (size_t)node_set.number_nodes, &num_nod_ns ) != NC_NOERR )
01981             {
01982                 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of nodes for set " << node_set.id );
01983             }
01984 
01985             /* Create variable array in which to store the node set node list */
01986             int node_ns = -1;
01987             INS_ID( wname, "node_ns%d", index );
01988             if( NC_NOERR != nc_def_var( ncFile, wname, NC_LONG, 1, &num_nod_ns, &node_ns ) )
01989             {
01990                 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create node set " << node_set.id << " node list" );
01991             }
01992 
01993             // Create distribution factor array
01994             int fact_ns = -1;
01995             INS_ID( wname, "dist_fact_ns%d", index );
01996             if( NC_NOERR != nc_def_var( ncFile, wname, NC_DOUBLE, 1, &num_nod_ns, &fact_ns ) )
01997             {
01998                 MB_SET_ERR( MB_FAILURE,
01999                             "WriteNCDF: failed to create node set " << node_set.id << " distribution factor list" );
02000             }
02001         }
02002     }
02003 
02004     /* Side set id array: */
02005 
02006     long non_empty_ss = 0;
02007     // Need to go through nodesets to compute # nodesets, some might be empty
02008     std::vector< NeumannSetData >::iterator ss_it;
02009     for( ss_it = sideset_data.begin(); ss_it != sideset_data.end(); ++ss_it )
02010     {
02011         if( 0 != ( *ss_it ).number_elements ) non_empty_ss++;
02012     }
02013 
02014     if( non_empty_ss > 0 )
02015     {
02016         int num_ss = -1;
02017         if( nc_def_dim( ncFile, "num_side_sets", non_empty_ss, &num_ss ) != NC_NOERR )
02018         {
02019             MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of side sets" );
02020         }
02021 
02022         /* Side set id status array: */
02023         int ss_idstat = -1, ss_idarr = -1;
02024         if( NC_NOERR != nc_def_var( ncFile, "ss_status", NC_LONG, 1, &num_ss, &ss_idstat ) )
02025         {
02026             MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define side set status" );
02027         }
02028 
02029         /* Side set id array: */
02030         if( NC_NOERR != nc_def_var( ncFile, "ss_prop1", NC_LONG, 1, &num_ss, &ss_idarr ) )
02031         {
02032             MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define side set property" );
02033         }
02034 
02035         /* Store property name as attribute of property array variable */
02036         if( NC_NOERR != nc_put_att_text( ncFile, ss_idarr, "name", strlen( "ID" ), "ID" ) )
02037         {
02038             MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to store side set property name ID" );
02039         }
02040 
02041         // Now, define the arrays needed for each side set
02042 
02043         int index = 0;
02044         for( unsigned int i = 0; i < sideset_data.size(); i++ )
02045         {
02046             NeumannSetData side_set = sideset_data[i];
02047 
02048             // Don't define an empty set
02049             if( 0 == side_set.number_elements ) continue;
02050 
02051             index++;
02052 
02053             int num_side_ss = -1;
02054             int elem_ss = -1, side_ss = -1;
02055             INS_ID( wname, "num_side_ss%d", index );
02056             if( nc_def_dim( ncFile, wname, (size_t)side_set.number_elements, &num_side_ss ) != NC_NOERR )
02057             {
02058                 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of sides in side set " << side_set.id );
02059             }
02060 
02061             INS_ID( wname, "elem_ss%d", index );
02062             if( NC_NOERR != nc_def_var( ncFile, wname, NC_LONG, 1, &num_side_ss, &elem_ss ) )
02063             {
02064                 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create element list for side set "
02065                                             << side_set.id ); /* Exit define mode and return */
02066             }
02067             INS_ID( wname, "side_ss%d", index );
02068             if( NC_NOERR != nc_def_var( ncFile, wname, NC_LONG, 1, &num_side_ss, &side_ss ) )
02069             {
02070                 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create side list for side set "
02071                                             << side_set.id ); /* Exit define mode and return */
02072             }
02073 
02074             // sideset distribution factors
02075             int num_df_ss = -1;
02076             INS_ID( wname, "num_df_ss%d", index );
02077             if( nc_def_dim( ncFile, wname, (size_t)side_set.ss_dist_factors.size(), &num_df_ss ) != NC_NOERR )
02078             {
02079                 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define number of dist factors in side set "
02080                                             << side_set.id ); /* Exit define mode and return */
02081             }
02082 
02083             /* Create variable array in which to store the side set distribution factors */
02084 
02085             int fact_ss = -1;
02086             INS_ID( wname, "dist_fact_ss%d", index );
02087             if( NC_NOERR != nc_def_var( ncFile, wname, NC_LONG, 1, &num_df_ss, &fact_ss ) )
02088             {
02089                 MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create dist factors list for side set "
02090                                             << side_set.id ); /* Exit define mode and return */
02091             }
02092         }
02093     }
02094 
02095     /* Node coordinate arrays: */
02096 
02097     int coord, name_coord, dims[3];
02098     dims[0] = num_dim;
02099     dims[1] = num_nodes;
02100     if( NC_NOERR != nc_def_var( ncFile, "coord", NC_DOUBLE, 2, dims, &coord ) )
02101     {
02102         MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define node coordinate array" );
02103     }
02104 
02105     /* Coordinate names array */
02106 
02107     dims[0] = num_dim;
02108     dims[1] = dim_str;
02109     if( NC_NOERR != nc_def_var( ncFile, "coor_names", NC_CHAR, 2, dims, &name_coord ) )
02110     {
02111         MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define coordinate name array" );
02112     }
02113 
02114     // Define genesis maps if required
02115 
02116     if( write_maps )
02117     {
02118         // Element map
02119         int elem_map = -1, elem_map2 = -1, node_map = -1;
02120         if( NC_NOERR != nc_def_var( ncFile, "elem_map", NC_LONG, 1, &num_elem, &elem_map ) )
02121         {
02122             MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create element map array" ); /* Exit define mode and return */
02123         }
02124 
02125         // Create the element number map
02126         if( NC_NOERR != nc_def_var( ncFile, "elem_num_map", NC_LONG, 1, &num_elem, &elem_map2 ) )
02127         {
02128             MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create element numbering map" ); /* Exit define mode
02129                                                                                               and return */
02130         }
02131 
02132         // Create node number map
02133         if( NC_NOERR != nc_def_var( ncFile, "node_num_map", NC_LONG, 1, &num_nodes, &node_map ) )
02134         {
02135             MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to create node numbering map array" ); /* Exit define mode and
02136                                                                                                  return */
02137         }
02138     }
02139 
02140     // Define qa records to be used
02141 
02142     int num_qa_rec = mesh_info.qaRecords.size() / 4;
02143     int num_qa     = -1;
02144 
02145     if( nc_def_dim( ncFile, "num_qa_rec", (long)num_qa_rec, &num_qa ) != NC_NOERR )
02146     {
02147         MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define qa record array size" );
02148     }
02149 
02150     // Define qa array
02151     int qa_title;
02152     dims[0] = num_qa;
02153     dims[1] = dim_four;
02154     dims[2] = dim_str;
02155     if( NC_NOERR != nc_def_var( ncFile, "qa_records", NC_CHAR, 3, dims, &qa_title ) )
02156     {
02157         MB_SET_ERR( MB_FAILURE, "WriteNCDF: failed to define qa record array" );
02158     }
02159 
02160     // Take it out of define mode
02161     if( NC_NOERR != nc_enddef( ncFile ) ) { MB_SET_ERR( MB_FAILURE, "WriteNCDF: Trouble leaving define mode" ); }
02162 
02163     return MB_SUCCESS;
02164 }
02165 
02166 ErrorCode WriteNCDF::open_file( const char* filename )
02167 {
02168     // Not a valid filename
02169     if( strlen( (const char*)filename ) == 0 ) { MB_SET_ERR( MB_FAILURE, "Output Exodus filename not specified" ); }
02170 
02171     int fail = nc_create( filename, NC_CLOBBER, &ncFile );
02172 
02173     // File couldn't be opened
02174     if( NC_NOERR != fail ) { MB_SET_ERR( MB_FAILURE, "Cannot open " << filename ); }
02175 
02176     return MB_SUCCESS;
02177 }
02178 
02179 ErrorCode WriteNCDF::get_sideset_elems( EntityHandle sideset, int current_sense, Range& forward_elems,
02180                                         Range& reverse_elems )
02181 {
02182     Range ss_elems, ss_meshsets;
02183 
02184     // Get the sense tag; don't need to check return, might be an error if the tag
02185     // hasn't been created yet
02186     Tag sense_tag = 0;
02187     mdbImpl->tag_get_handle( "SENSE", 1, MB_TYPE_INTEGER, sense_tag );
02188 
02189     // Get the entities in this set
02190     ErrorCode result = mdbImpl->get_entities_by_handle( sideset, ss_elems, true );
02191     if( MB_FAILURE == result ) return result;
02192 
02193     // Now remove the meshsets into the ss_meshsets; first find the first meshset,
02194     Range::iterator range_iter = ss_elems.begin();
02195     while( TYPE_FROM_HANDLE( *range_iter ) != MBENTITYSET && range_iter != ss_elems.end() )
02196         ++range_iter;
02197 
02198     // Then, if there are some, copy them into ss_meshsets and erase from ss_elems
02199     if( range_iter != ss_elems.end() )
02200     {
02201         std::copy( range_iter, ss_elems.end(), range_inserter( ss_meshsets ) );
02202         ss_elems.erase( range_iter, ss_elems.end() );
02203     }
02204 
02205     // OK, for the elements, check the sense of this set and copy into the right range
02206     // (if the sense is 0, copy into both ranges)
02207 
02208     // Need to step forward on list until we reach the right dimension
02209     Range::iterator dum_it = ss_elems.end();
02210     --dum_it;
02211     int target_dim = CN::Dimension( TYPE_FROM_HANDLE( *dum_it ) );
02212     dum_it         = ss_elems.begin();
02213     while( target_dim != CN::Dimension( TYPE_FROM_HANDLE( *dum_it ) ) && dum_it != ss_elems.end() )
02214         ++dum_it;
02215 
02216     if( current_sense == 1 || current_sense == 0 ) std::copy( dum_it, ss_elems.end(), range_inserter( forward_elems ) );
02217     if( current_sense == -1 || current_sense == 0 )
02218         std::copy( dum_it, ss_elems.end(), range_inserter( reverse_elems ) );
02219 
02220     // Now loop over the contained meshsets, getting the sense of those and calling this
02221     // function recursively
02222     for( range_iter = ss_meshsets.begin(); range_iter != ss_meshsets.end(); ++range_iter )
02223     {
02224         // First get the sense; if it's not there, by convention it's forward
02225         int this_sense;
02226         if( 0 == sense_tag || MB_FAILURE == mdbImpl->tag_get_data( sense_tag, &( *range_iter ), 1, &this_sense ) )
02227             this_sense = 1;
02228 
02229         // Now get all the entities on this meshset, with the proper (possibly reversed) sense
02230         get_sideset_elems( *range_iter, this_sense * current_sense, forward_elems, reverse_elems );
02231     }
02232 
02233     return result;
02234 }
02235 
02236 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines