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