WriteSLAC.cpp

Go to the documentation of this file.

/**
 * MOAB, a Mesh-Oriented datABase, is a software component for creating,
 * storing and accessing finite element mesh data.
 *
 * Copyright 2004 Sandia Corporation. Under the terms of Contract
 * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
 * retains certain rights in this software.
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 *
 */
 
#ifdef WIN32
#ifdef _DEBUG
// turn off warnings that say they debugging identifier has been truncated
// this warning comes up when using some STL containers
#pragma warning( disable : 4786 )
#endif
#endif
 
#include "WriteSLAC.hpp"
 
#include <utility>
#include <algorithm>
#include <ctime>
#include <string>
#include <vector>
#include <cstdio>
#include <cstring>
#include <iostream>
#include <cassert>
 
#include "netcdf.h"
#include "moab/Interface.hpp"
#include "moab/Range.hpp"
#include "moab/CN.hpp"
#include "Internals.hpp"
#include "ExoIIUtil.hpp"
#include "MBTagConventions.hpp"
#include "moab/WriteUtilIface.hpp"
 
#ifndef MOAB_HAVE_NETCDF
#error Attempt to compile WriteSLAC with NetCDF disabled.
#endif
 
namespace moab
{
 
#define INS_ID( stringvar, prefix, id ) sprintf( stringvar, prefix, id )
 
#define GET_VAR( name, id, dims ) \
 { \
 ( id ) = -1; \
 int gvfail = nc_inq_varid( ncFile, name, &( id ) ); \
 if( NC_NOERR == gvfail ) \
 { \
 int ndims; \
 gvfail = nc_inq_varndims( ncFile, id, &ndims ); \
 if( NC_NOERR == gvfail ) \
 { \
 ( dims ).resize( ndims ); \
 gvfail = nc_inq_vardimid( ncFile, id, &( dims )[0] ); \
 } \
 } \
 }
 
WriterIface* WriteSLAC::factory( Interface* iface )
{
 return new WriteSLAC( iface );
}
 
WriteSLAC::WriteSLAC( Interface* impl ) : mbImpl( impl ), ncFile( 0 )
{
 assert( impl != NULL );
 
 impl->query_interface( mWriteIface );
 
 // Initialize in case tag_get_handle fails below
 //! get and cache predefined tag handles
 int negone = -1;
 impl->tag_get_handle( MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mMaterialSetTag, MB_TAG_SPARSE | MB_TAG_CREAT,
 &negone );
 
 impl->tag_get_handle( DIRICHLET_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mDirichletSetTag, MB_TAG_SPARSE | MB_TAG_CREAT,
 &negone );
 
 impl->tag_get_handle( NEUMANN_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mNeumannSetTag, MB_TAG_SPARSE | MB_TAG_CREAT,
 &negone );
 
 mGlobalIdTag = impl->globalId_tag();
 
 int dum_val = -1;
 impl->tag_get_handle( "__matSetIdTag", 1, MB_TYPE_INTEGER, mMatSetIdTag, MB_TAG_DENSE | MB_TAG_CREAT, &dum_val );
 
 impl->tag_get_handle( "WriteSLAC element mark", 1, MB_TYPE_BIT, mEntityMark, MB_TAG_CREAT );
}
 
WriteSLAC::~WriteSLAC()
{
 mbImpl->release_interface( mWriteIface );
 mbImpl->tag_delete( mEntityMark );
}
 
void WriteSLAC::reset_matset( std::vector< WriteSLAC::MaterialSetData >& matset_info )
{
 std::vector< WriteSLAC::MaterialSetData >::iterator iter;
 
 for( iter = matset_info.begin(); iter != matset_info.end(); ++iter )
 delete( *iter ).elements;
}
 
ErrorCode WriteSLAC::write_file( const char* file_name,
 const bool overwrite,
 const FileOptions&,
 const EntityHandle* ent_handles,
 const int num_sets,
 const std::vector< std::string >& /* qa_list */,
 const Tag* /* tag_list */,
 int /* num_tags */,
 int /* export_dimension */ )
{
 assert( 0 != mMaterialSetTag && 0 != mNeumannSetTag && 0 != mDirichletSetTag );
 
 // Check the file name
 if( NULL == strstr( file_name, ".ncdf" ) ) return MB_FAILURE;
 
 std::vector< EntityHandle > matsets, dirsets, neusets, entities;
 
 fileName = file_name;
 
 // Separate into material sets, dirichlet sets, neumann sets
 
 if( num_sets == 0 )
 {
 // Default to all defined sets
 Range this_range;
 mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &mMaterialSetTag, NULL, 1, this_range );
 std::copy( this_range.begin(), this_range.end(), std::back_inserter( matsets ) );
 this_range.clear();
 mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &mDirichletSetTag, NULL, 1, this_range );
 std::copy( this_range.begin(), this_range.end(), std::back_inserter( dirsets ) );
 this_range.clear();
 mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &mNeumannSetTag, NULL, 1, this_range );
 std::copy( this_range.begin(), this_range.end(), std::back_inserter( neusets ) );
 }
 else
 {
 int dummy;
 for( const EntityHandle* iter = ent_handles; iter < ent_handles + num_sets; ++iter )
 {
 if( MB_SUCCESS == mbImpl->tag_get_data( mMaterialSetTag, &( *iter ), 1, &dummy ) )
 matsets.push_back( *iter );
 else if( MB_SUCCESS == mbImpl->tag_get_data( mDirichletSetTag, &( *iter ), 1, &dummy ) )
 dirsets.push_back( *iter );
 else if( MB_SUCCESS == mbImpl->tag_get_data( mNeumannSetTag, &( *iter ), 1, &dummy ) )
 neusets.push_back( *iter );
 }
 }
 
 // If there is nothing to write just return.
 if( matsets.empty() && dirsets.empty() && neusets.empty() ) return MB_FILE_WRITE_ERROR;
 
 std::vector< WriteSLAC::MaterialSetData > matset_info;
 std::vector< WriteSLAC::DirichletSetData > dirset_info;
 std::vector< WriteSLAC::NeumannSetData > neuset_info;
 
 MeshInfo mesh_info;
 
 matset_info.clear();
 if( gather_mesh_information( mesh_info, matset_info, neuset_info, dirset_info, matsets, neusets, dirsets ) !=
 MB_SUCCESS )
 {
 reset_matset( matset_info );
 return MB_FAILURE;
 }
 
 // Try to open the file after gather mesh info succeeds
 int fail = nc_create( file_name, overwrite ? NC_CLOBBER : NC_NOCLOBBER, &ncFile );
 if( NC_NOERR != fail )
 {
 reset_matset( matset_info );
 return MB_FAILURE;
 }
 
 if( initialize_file( mesh_info ) != MB_SUCCESS )
 {
 reset_matset( matset_info );
 return MB_FAILURE;
 }
 
 if( write_nodes( mesh_info.num_nodes, mesh_info.nodes, mesh_info.num_dim ) != MB_SUCCESS )
 {
 reset_matset( matset_info );
 return MB_FAILURE;
 }
 
 if( write_matsets( mesh_info, matset_info, neuset_info ) )
 {
 reset_matset( matset_info );
 return MB_FAILURE;
 }
 
 fail = nc_close( ncFile );
 if( NC_NOERR != fail ) return MB_FAILURE;
 
 return MB_SUCCESS;
}
 
ErrorCode WriteSLAC::gather_mesh_information( MeshInfo& mesh_info,
 std::vector< WriteSLAC::MaterialSetData >& matset_info,
 std::vector< WriteSLAC::NeumannSetData >& neuset_info,
 std::vector< WriteSLAC::DirichletSetData >& dirset_info,
 std::vector< EntityHandle >& matsets,
 std::vector< EntityHandle >& neusets,
 std::vector< EntityHandle >& dirsets )
{
 std::vector< EntityHandle >::iterator vector_iter, end_vector_iter;
 
 mesh_info.num_nodes = 0;
 mesh_info.num_elements = 0;
 mesh_info.num_matsets = 0;
 
 int id = 0;
 
 vector_iter = matsets.begin();
 end_vector_iter = matsets.end();
 
 mesh_info.num_matsets = matsets.size();
 
 std::vector< EntityHandle > parent_meshsets;
 
 // Clean out the bits for the element mark
 mbImpl->tag_delete( mEntityMark );
 mbImpl->tag_get_handle( "WriteSLAC element mark", 1, MB_TYPE_BIT, mEntityMark, MB_TAG_CREAT );
 
 int highest_dimension_of_element_matsets = 0;
 
 for( vector_iter = matsets.begin(); vector_iter != matsets.end(); ++vector_iter )
 {
 WriteSLAC::MaterialSetData matset_data;
 matset_data.elements = new Range;
 
 // For the purpose of qa records, get the parents of these matsets
 if( mbImpl->get_parent_meshsets( *vector_iter, parent_meshsets ) != MB_SUCCESS ) return MB_FAILURE;
 
 // Get all Entity Handles in the mesh set
 Range dummy_range;
 mbImpl->get_entities_by_handle( *vector_iter, dummy_range, true );
 
 // Wait a minute, we are doing some filtering here that doesn't make sense at this level CJS
 
 // Find the dimension of the last entity in this range
 Range::iterator entity_iter = dummy_range.end();
 entity_iter = dummy_range.end();
 --entity_iter;
 int this_dim = CN::Dimension( TYPE_FROM_HANDLE( *entity_iter ) );
 entity_iter = dummy_range.begin();
 while( entity_iter != dummy_range.end() && CN::Dimension( TYPE_FROM_HANDLE( *entity_iter ) ) != this_dim )
 ++entity_iter;
 
 if( entity_iter != dummy_range.end() )
 std::copy( entity_iter, dummy_range.end(), range_inserter( *( matset_data.elements ) ) );
 
 assert( matset_data.elements->begin() == matset_data.elements->end() ||
 CN::Dimension( TYPE_FROM_HANDLE( *( matset_data.elements->begin() ) ) ) == this_dim );
 
 // Get the matset's id
 if( mbImpl->tag_get_data( mMaterialSetTag, &( *vector_iter ), 1, &id ) != MB_SUCCESS )
 {
 MB_SET_ERR( MB_FAILURE, "Couldn't get matset id from a tag for an element matset" );
 }
 
 matset_data.id = id;
 matset_data.number_attributes = 0;
 
 // Iterate through all the elements in the meshset
 Range::iterator elem_range_iter, end_elem_range_iter;
 elem_range_iter = matset_data.elements->begin();
 end_elem_range_iter = matset_data.elements->end();
 
 // Get the entity type for this matset, verifying that it's the same for all elements
 // THIS ASSUMES HANDLES SORT BY TYPE!!!
 EntityType entity_type = TYPE_FROM_HANDLE( *elem_range_iter );
 --end_elem_range_iter;
 if( entity_type != TYPE_FROM_HANDLE( *( end_elem_range_iter++ ) ) )
 {
 MB_SET_ERR( MB_FAILURE, "Entities in matset " << id << " not of common type" );
 }
 
 int dimension = -1;
 if( entity_type == MBQUAD || entity_type == MBTRI )
 dimension = 3; // Output shells by default
 else if( entity_type == MBEDGE )
 dimension = 2;
 else
 dimension = CN::Dimension( entity_type );
 
 if( dimension > highest_dimension_of_element_matsets ) highest_dimension_of_element_matsets = dimension;
 
 matset_data.moab_type = mbImpl->type_from_handle( *( matset_data.elements->begin() ) );
 if( MBMAXTYPE == matset_data.moab_type ) return MB_FAILURE;
 
 std::vector< EntityHandle > tmp_conn;
 mbImpl->get_connectivity( &( *( matset_data.elements->begin() ) ), 1, tmp_conn );
 matset_data.element_type =
 ExoIIUtil::get_element_type_from_num_verts( tmp_conn.size(), entity_type, dimension );
 
 if( matset_data.element_type == EXOII_MAX_ELEM_TYPE )
 {
 MB_SET_ERR( MB_FAILURE, "Element type in matset " << id << " didn't get set correctly" );
 }
 
 matset_data.number_nodes_per_element = ExoIIUtil::VerticesPerElement[matset_data.element_type];
 
 // Number of nodes for this matset
 matset_data.number_elements = matset_data.elements->size();
 
 // Total number of elements
 mesh_info.num_elements += matset_data.number_elements;
 
 // Get the nodes for the elements
 mWriteIface->gather_nodes_from_elements( *matset_data.elements, mEntityMark, mesh_info.nodes );
 
 if( !neusets.empty() )
 {
 // If there are neusets, keep track of which elements are being written out
 for( Range::iterator iter = matset_data.elements->begin(); iter != matset_data.elements->end(); ++iter )
 {
 unsigned char bit = 0x1;
 mbImpl->tag_set_data( mEntityMark, &( *iter ), 1, &bit );
 }
 }
 
 matset_info.push_back( matset_data );
 }
 
 // If user hasn't entered dimension, we figure it out
 if( mesh_info.num_dim == 0 )
 {
 // Never want 1 or zero dimensions
 if( highest_dimension_of_element_matsets < 2 )
 mesh_info.num_dim = 3;
 else
 mesh_info.num_dim = highest_dimension_of_element_matsets;
 }
 
 Range::iterator range_iter, end_range_iter;
 range_iter = mesh_info.nodes.begin();
 end_range_iter = mesh_info.nodes.end();
 
 mesh_info.num_nodes = mesh_info.nodes.size();
 
 //------dirsets--------
 
 vector_iter = dirsets.begin();
 end_vector_iter = dirsets.end();
 
 for( ; vector_iter != end_vector_iter; ++vector_iter )
 {
 WriteSLAC::DirichletSetData dirset_data;
 dirset_data.id = 0;
 dirset_data.number_nodes = 0;
 
 // Get the dirset's id
 if( mbImpl->tag_get_data( mDirichletSetTag, &( *vector_iter ), 1, &id ) != MB_SUCCESS )
 {
 MB_SET_ERR( MB_FAILURE, "Couldn't get id tag for dirset " << id );
 }
 
 dirset_data.id = id;
 
 std::vector< EntityHandle > node_vector;
 // Get the nodes of the dirset that are in mesh_info.nodes
 if( mbImpl->get_entities_by_handle( *vector_iter, node_vector, true ) != MB_SUCCESS )
 {
 MB_SET_ERR( MB_FAILURE, "Couldn't get nodes in dirset " << id );
 }
 
 std::vector< EntityHandle >::iterator iter, end_iter;
 iter = node_vector.begin();
 end_iter = node_vector.end();
 
 int j = 0;
 unsigned char node_marked = 0;
 ErrorCode result;
 for( ; iter != end_iter; ++iter )
 {
 if( TYPE_FROM_HANDLE( *iter ) != MBVERTEX ) continue;
 result = mbImpl->tag_get_data( mEntityMark, &( *iter ), 1, &node_marked );MB_CHK_SET_ERR( result, "Couldn't get mark data" );
 
 if( 0x1 == node_marked ) dirset_data.nodes.push_back( *iter );
 j++;
 }
 
 dirset_data.number_nodes = dirset_data.nodes.size();
 dirset_info.push_back( dirset_data );
 }
 
 //------neusets--------
 vector_iter = neusets.begin();
 end_vector_iter = neusets.end();
 
 for( ; vector_iter != end_vector_iter; ++vector_iter )
 {
 WriteSLAC::NeumannSetData neuset_data;
 
 // Get the neuset's id
 if( mbImpl->tag_get_data( mNeumannSetTag, &( *vector_iter ), 1, &id ) != MB_SUCCESS ) return MB_FAILURE;
 
 neuset_data.id = id;
 neuset_data.mesh_set_handle = *vector_iter;
 
 // Get the sides in two lists, one forward the other reverse; starts with forward sense
 // by convention
 Range forward_elems, reverse_elems;
 if( get_neuset_elems( *vector_iter, 0, forward_elems, reverse_elems ) == MB_FAILURE ) return MB_FAILURE;
 
 ErrorCode result = get_valid_sides( forward_elems, 1, neuset_data );MB_CHK_SET_ERR( result, "Couldn't get valid sides data" );
 result = get_valid_sides( reverse_elems, -1, neuset_data );MB_CHK_SET_ERR( result, "Couldn't get valid sides data" );
 
 neuset_data.number_elements = neuset_data.elements.size();
 neuset_info.push_back( neuset_data );
 }
 
 // Get information about interior/exterior tets/hexes, and mark matset ids
 return gather_interior_exterior( mesh_info, matset_info, neuset_info );
}
 
ErrorCode WriteSLAC::get_valid_sides( Range& elems, const int sense, WriteSLAC::NeumannSetData& neuset_data )
{
 // This is where we see if underlying element of side set element is included in output
 
 unsigned char element_marked = 0;
 ErrorCode result;
 for( Range::iterator iter = elems.begin(); iter != elems.end(); ++iter )
 {
 // Should insert here if "side" is a quad/tri on a quad/tri mesh
 result = mbImpl->tag_get_data( mEntityMark, &( *iter ), 1, &element_marked );MB_CHK_SET_ERR( result, "Couldn't get mark data" );
 
 if( 0x1 == element_marked )
 {
 neuset_data.elements.push_back( *iter );
 
 // TJT TODO: the sense should really be # edges + 1or2
 neuset_data.side_numbers.push_back( ( sense == 1 ? 1 : 2 ) );
 }
 else
 { // Then "side" is probably a quad/tri on a hex/tet mesh
 std::vector< EntityHandle > parents;
 int dimension = CN::Dimension( TYPE_FROM_HANDLE( *iter ) );
 
 // Get the adjacent parent element of "side"
 if( mbImpl->get_adjacencies( &( *iter ), 1, dimension + 1, false, parents ) != MB_SUCCESS )
 {
 MB_SET_ERR( MB_FAILURE, "Couldn't get adjacencies for neuset" );
 }
 
 if( !parents.empty() )
 {
 // Make sure the adjacent parent element will be output
 for( unsigned int k = 0; k < parents.size(); k++ )
 {
 result = mbImpl->tag_get_data( mEntityMark, &( parents[k] ), 1, &element_marked );MB_CHK_SET_ERR( result, "Couldn't get mark data" );
 
 int side_no, this_sense, this_offset;
 if( 0x1 == element_marked &&
 mbImpl->side_number( parents[k], *iter, side_no, this_sense, this_offset ) == MB_SUCCESS &&
 this_sense == sense )
 {
 neuset_data.elements.push_back( parents[k] );
 neuset_data.side_numbers.push_back( side_no + 1 );
 break;
 }
 }
 }
 else
 {
 MB_SET_ERR( MB_FAILURE, "No parent element exists for element in neuset " << neuset_data.id );
 }
 }
 }
 
 return MB_SUCCESS;
}
 
ErrorCode WriteSLAC::write_nodes( const int num_nodes, const Range& nodes, const int dimension )
{
 // See if should transform coordinates
 ErrorCode result;
 Tag trans_tag;
 result = mbImpl->tag_get_handle( MESH_TRANSFORM_TAG_NAME, 16, MB_TYPE_DOUBLE, trans_tag );
 bool transform_needed = true;
 if( result == MB_TAG_NOT_FOUND ) transform_needed = false;
 
 int num_coords_to_fill = transform_needed ? 3 : dimension;
 
 std::vector< double* > coord_arrays( 3 );
 coord_arrays[0] = new double[num_nodes];
 coord_arrays[1] = new double[num_nodes];
 coord_arrays[2] = NULL;
 
 if( num_coords_to_fill == 3 ) coord_arrays[2] = new double[num_nodes];
 
 result = mWriteIface->get_node_coords( dimension, num_nodes, nodes, mGlobalIdTag, 0, coord_arrays );
 if( result != MB_SUCCESS )
 {
 delete[] coord_arrays[0];
 delete[] coord_arrays[1];
 if( coord_arrays[2] ) delete[] coord_arrays[2];
 return result;
 }
 
 if( transform_needed )
 {
 double trans_matrix[16];
 const EntityHandle mesh = 0;
 result = mbImpl->tag_get_data( trans_tag, &mesh, 1, trans_matrix );MB_CHK_SET_ERR( result, "Couldn't get transform data" );
 
 for( int i = 0; i < num_nodes; i++ )
 {
 double vec1[3];
 double vec2[3];
 
 vec2[0] = coord_arrays[0][i];
 vec2[1] = coord_arrays[1][i];
 vec2[2] = coord_arrays[2][i];
 
 for( int row = 0; row < 3; row++ )
 {
 vec1[row] = 0.0;
 for( int col = 0; col < 3; col++ )
 vec1[row] += ( trans_matrix[( row * 4 ) + col] * vec2[col] );
 }
 
 coord_arrays[0][i] = vec1[0];
 coord_arrays[1][i] = vec1[1];
 coord_arrays[2][i] = vec1[2];
 }
 }
 
 // Write the nodes
 int nc_var = -1;
 std::vector< int > dims;
 GET_VAR( "coords", nc_var, dims );
 if( -1 == nc_var ) return MB_FAILURE;
 size_t start[2] = { 0, 0 }, count[2] = { static_cast< size_t >( num_nodes ), 1 };
 int fail = nc_put_vara_double( ncFile, nc_var, start, count, coord_arrays[0] );
 if( NC_NOERR != fail ) return MB_FAILURE;
 start[1] = 1;
 fail = nc_put_vara_double( ncFile, nc_var, start, count, coord_arrays[1] );
 if( NC_NOERR != fail ) return MB_FAILURE;
 start[1] = 2;
 fail = nc_put_vara_double( ncFile, nc_var, start, count, coord_arrays[2] );
 if( NC_NOERR != fail ) return MB_FAILURE;
 
 delete[] coord_arrays[0];
 delete[] coord_arrays[1];
 if( coord_arrays[2] ) delete[] coord_arrays[2];
 
 return MB_SUCCESS;
}
 
ErrorCode WriteSLAC::gather_interior_exterior( MeshInfo& mesh_info,
 std::vector< WriteSLAC::MaterialSetData >& matset_data,
 std::vector< WriteSLAC::NeumannSetData >& neuset_data )
{
 // Need to assign a tag with the matset id
 Tag matset_id_tag;
 unsigned int i;
 int dum = -1;
 ErrorCode result =
 mbImpl->tag_get_handle( "__matset_id", 4, MB_TYPE_INTEGER, matset_id_tag, MB_TAG_DENSE | MB_TAG_CREAT, &dum );
 if( MB_SUCCESS != result ) return result;
 
 Range::iterator rit;
 mesh_info.num_int_hexes = mesh_info.num_int_tets = 0;
 
 for( i = 0; i < matset_data.size(); i++ )
 {
 WriteSLAC::MaterialSetData matset = matset_data[i];
 if( matset.moab_type == MBHEX )
 mesh_info.num_int_hexes += matset.elements->size();
 else if( matset.moab_type == MBTET )
 mesh_info.num_int_tets += matset.elements->size();
 else
 {
 std::cout << "WriteSLAC doesn't support elements of type " << CN::EntityTypeName( matset.moab_type )
 << std::endl;
 continue;
 }
 
 for( rit = matset.elements->begin(); rit != matset.elements->end(); ++rit )
 {
 result = mbImpl->tag_set_data( mMatSetIdTag, &( *rit ), 1, &( matset.id ) );
 if( MB_SUCCESS != result ) return result;
 }
 }
 
 // Now go through the neumann sets, pulling out the hexes with faces on the
 // boundary
 std::vector< EntityHandle >::iterator vit;
 for( i = 0; i < neuset_data.size(); i++ )
 {
 WriteSLAC::NeumannSetData neuset = neuset_data[i];
 for( vit = neuset.elements.begin(); vit != neuset.elements.end(); ++vit )
 {
 if( TYPE_FROM_HANDLE( *vit ) == MBHEX )
 mesh_info.bdy_hexes.insert( *vit );
 else if( TYPE_FROM_HANDLE( *vit ) == MBTET )
 mesh_info.bdy_tets.insert( *vit );
 }
 }
 
 // Now we have the number of bdy hexes and tets, we know how many interior ones
 // there are too
 mesh_info.num_int_hexes -= mesh_info.bdy_hexes.size();
 mesh_info.num_int_tets -= mesh_info.bdy_tets.size();
 
 return MB_SUCCESS;
}
 
ErrorCode WriteSLAC::write_matsets( MeshInfo& mesh_info,
 std::vector< WriteSLAC::MaterialSetData >& matset_data,
 std::vector< WriteSLAC::NeumannSetData >& neuset_data )
{
 unsigned int i;
 std::vector< int > connect;
 const EntityHandle* connecth;
 int num_connecth;
 ErrorCode result;
 
 // First write the interior hexes
 int hex_conn = -1;
 std::vector< int > dims;
 if( mesh_info.bdy_hexes.size() != 0 || mesh_info.num_int_hexes != 0 )
 {
 GET_VAR( "hexahedron_interior", hex_conn, dims );
 if( -1 == hex_conn ) return MB_FAILURE;
 }
 connect.reserve( 13 );
 Range::iterator rit;
 
 int elem_num = 0;
 WriteSLAC::MaterialSetData matset;
 size_t start[2] = { 0, 0 }, count[2] = { 1, 1 };
 int fail;
 for( i = 0; i < matset_data.size(); i++ )
 {
 matset = matset_data[i];
 if( matset.moab_type != MBHEX ) continue;
 
 int id = matset.id;
 connect[0] = id;
 
 for( rit = matset.elements->begin(); rit != matset.elements->end(); ++rit )
 {
 // Skip if it's on the bdy
 if( mesh_info.bdy_hexes.find( *rit ) != mesh_info.bdy_hexes.end() ) continue;
 
 // Get the connectivity of this element
 result = mbImpl->get_connectivity( *rit, connecth, num_connecth );
 if( MB_SUCCESS != result ) return result;
 
 // Get the vertex ids
 result = mbImpl->tag_get_data( mGlobalIdTag, connecth, num_connecth, &connect[1] );
 if( MB_SUCCESS != result ) return result;
 
 // Put the variable at the right position
 start[0] = elem_num++;
 count[1] = 9;
 
 // Write the data
 fail = nc_put_vara_int( ncFile, hex_conn, start, count, &connect[0] );
 if( NC_NOERR != fail ) return MB_FAILURE;
 }
 }
 
 int tet_conn = -1;
 if( mesh_info.bdy_tets.size() != 0 || mesh_info.num_int_tets != 0 )
 {
 GET_VAR( "tetrahedron_interior", tet_conn, dims );
 if( -1 == tet_conn ) return MB_FAILURE;
 }
 
 // Now the interior tets
 elem_num = 0;
 for( i = 0; i < matset_data.size(); i++ )
 {
 matset = matset_data[i];
 if( matset.moab_type != MBTET ) continue;
 
 int id = matset.id;
 connect[0] = id;
 elem_num = 0;
 for( rit = matset.elements->begin(); rit != matset.elements->end(); ++rit )
 {
 // Skip if it's on the bdy
 if( mesh_info.bdy_tets.find( *rit ) != mesh_info.bdy_tets.end() ) continue;
 
 // Get the connectivity of this element
 result = mbImpl->get_connectivity( *rit, connecth, num_connecth );
 if( MB_SUCCESS != result ) return result;
 
 // Get the vertex ids
 result = mbImpl->tag_get_data( mGlobalIdTag, connecth, num_connecth, &connect[1] );
 if( MB_SUCCESS != result ) return result;
 
 // Put the variable at the right position
 start[0] = elem_num++;
 count[1] = 5;
 fail = nc_put_vara_int( ncFile, tet_conn, start, count, &connect[0] );
 // Write the data
 if( NC_NOERR != fail ) return MB_FAILURE;
 }
 }
 
 // Now the exterior hexes
 if( mesh_info.bdy_hexes.size() != 0 )
 {
 hex_conn = -1;
 GET_VAR( "hexahedron_exterior", hex_conn, dims );
 if( -1 == hex_conn ) return MB_FAILURE;
 
 connect.reserve( 15 );
 elem_num = 0;
 
 // Write the elements
 for( rit = mesh_info.bdy_hexes.begin(); rit != mesh_info.bdy_hexes.end(); ++rit )
 {
 // Get the material set for this hex
 result = mbImpl->tag_get_data( mMatSetIdTag, &( *rit ), 1, &connect[0] );
 if( MB_SUCCESS != result ) return result;
 
 // Get the connectivity of this element
 result = mbImpl->get_connectivity( *rit, connecth, num_connecth );
 if( MB_SUCCESS != result ) return result;
 
 // Get the vertex ids
 result = mbImpl->tag_get_data( mGlobalIdTag, connecth, num_connecth, &connect[1] );
 if( MB_SUCCESS != result ) return result;
 
 // Preset side numbers
 for( i = 9; i < 15; i++ )
 connect[i] = -1;
 
 // Now write the side numbers
 for( i = 0; i < neuset_data.size(); i++ )
 {
 std::vector< EntityHandle >::iterator vit =
 std::find( neuset_data[i].elements.begin(), neuset_data[i].elements.end(), *rit );
 while( vit != neuset_data[i].elements.end() )
 {
 // Have a side - get the side # and put in connect array
 int side_no = neuset_data[i].side_numbers[vit - neuset_data[i].elements.begin()];
 connect[9 + side_no] = neuset_data[i].id;
 ++vit;
 vit = std::find( vit, neuset_data[i].elements.end(), *rit );
 }
 }
 
 // Put the variable at the right position
 start[0] = elem_num++;
 count[1] = 15;
 fail = nc_put_vara_int( ncFile, hex_conn, start, count, &connect[0] );
 // Write the data
 if( NC_NOERR != fail ) return MB_FAILURE;
 }
 }
 
 // Now the exterior tets
 if( mesh_info.bdy_tets.size() != 0 )
 {
 tet_conn = -1;
 GET_VAR( "tetrahedron_exterior", tet_conn, dims );
 if( -1 == tet_conn ) return MB_FAILURE;
 
 connect.reserve( 9 );
 elem_num = 0;
 
 // Write the elements
 for( rit = mesh_info.bdy_tets.begin(); rit != mesh_info.bdy_tets.end(); ++rit )
 {
 // Get the material set for this tet
 result = mbImpl->tag_get_data( mMatSetIdTag, &( *rit ), 1, &connect[0] );
 if( MB_SUCCESS != result ) return result;
 
 // Get the connectivity of this element
 result = mbImpl->get_connectivity( *rit, connecth, num_connecth );
 if( MB_SUCCESS != result ) return result;
 
 // Get the vertex ids
 result = mbImpl->tag_get_data( mGlobalIdTag, connecth, num_connecth, &connect[1] );
 if( MB_SUCCESS != result ) return result;
 
 // Preset side numbers
 for( i = 5; i < 9; i++ )
 connect[i] = -1;
 
 // Now write the side numbers
 for( i = 0; i < neuset_data.size(); i++ )
 {
 std::vector< EntityHandle >::iterator vit =
 std::find( neuset_data[i].elements.begin(), neuset_data[i].elements.end(), *rit );
 while( vit != neuset_data[i].elements.end() )
 {
 // Have a side - get the side # and put in connect array
 int side_no = neuset_data[i].side_numbers[vit - neuset_data[i].elements.begin()];
 connect[5 + side_no] = neuset_data[i].id;
 ++vit;
 vit = std::find( vit, neuset_data[i].elements.end(), *rit );
 }
 }
 
 // Put the variable at the right position
 start[0] = elem_num++;
 count[1] = 9;
 fail = nc_put_vara_int( ncFile, tet_conn, start, count, &connect[0] );
 // Write the data
 if( NC_NOERR != fail ) return MB_FAILURE;
 }
 }
 
 return MB_SUCCESS;
}
 
ErrorCode WriteSLAC::initialize_file( MeshInfo& mesh_info )
{
 // Perform the initializations
 
 int coord_size = -1, ncoords = -1;
 // Initialization to avoid warnings on Linux
 int hexinterior = -1, hexinteriorsize, hexexterior = -1, hexexteriorsize = -1;
 int tetinterior = -1, tetinteriorsize, tetexterior = -1, tetexteriorsize = -1;
 
 if( nc_def_dim( ncFile, "coord_size", (size_t)mesh_info.num_dim, &coord_size ) != NC_NOERR )
 {
 MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to define number of dimensions" );
 }
 
 if( nc_def_dim( ncFile, "ncoords", (size_t)mesh_info.num_nodes, &ncoords ) != NC_NOERR )
 {
 MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to define number of nodes" );
 }
 
 if( 0 != mesh_info.num_int_hexes &&
 nc_def_dim( ncFile, "hexinterior", (size_t)mesh_info.num_int_hexes, &hexinterior ) != NC_NOERR )
 {
 MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to define number of interior hex elements" );
 }
 
 if( nc_def_dim( ncFile, "hexinteriorsize", (size_t)9, &hexinteriorsize ) != NC_NOERR )
 {
 MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to define interior hex element size" );
 }
 
 if( 0 != mesh_info.bdy_hexes.size() &&
 nc_def_dim( ncFile, "hexexterior", (size_t)mesh_info.bdy_hexes.size(), &hexexterior ) != NC_NOERR )
 {
 MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to define number of exterior hex elements" );
 }
 
 if( nc_def_dim( ncFile, "hexexteriorsize", (size_t)15, &hexexteriorsize ) != NC_NOERR )
 {
 MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to define exterior hex element size" );
 }
 
 if( 0 != mesh_info.num_int_tets &&
 nc_def_dim( ncFile, "tetinterior", (size_t)mesh_info.num_int_tets, &tetinterior ) != NC_NOERR )
 {
 MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to define number of interior tet elements" );
 }
 
 if( nc_def_dim( ncFile, "tetinteriorsize", (size_t)5, &tetinteriorsize ) != NC_NOERR )
 {
 MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to define interior tet element size" );
 }
 
 if( 0 != mesh_info.bdy_tets.size() &&
 nc_def_dim( ncFile, "tetexterior", (size_t)mesh_info.bdy_tets.size(), &tetexterior ) != NC_NOERR )
 {
 MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to define number of exterior tet elements" );
 }
 
 if( nc_def_dim( ncFile, "tetexteriorsize", (size_t)9, &tetexteriorsize ) != NC_NOERR )
 {
 MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to define exterior tet element size" );
 }
 
 /* ...and some variables */
 
 int dims[2];
 dims[0] = hexinterior;
 dims[1] = hexinteriorsize;
 int dum_var;
 if( 0 != mesh_info.num_int_hexes &&
 NC_NOERR != nc_def_var( ncFile, "hexahedron_interior", NC_LONG, 2, dims, &dum_var ) )
 {
 MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to create connectivity array for interior hexes" );
 }
 
 dims[0] = hexexterior;
 dims[1] = hexexteriorsize;
 if( 0 != mesh_info.bdy_hexes.size() &&
 NC_NOERR != nc_def_var( ncFile, "hexahedron_exterior", NC_LONG, 2, dims, &dum_var ) )
 {
 MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to create connectivity array for exterior hexes" );
 }
 
 dims[0] = tetinterior;
 dims[1] = tetinteriorsize;
 if( 0 != mesh_info.num_int_tets &&
 NC_NOERR != nc_def_var( ncFile, "tetrahedron_exterior", NC_LONG, 2, dims, &dum_var ) )
 {
 MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to create connectivity array for interior tets" );
 }
 
 dims[0] = tetexterior;
 dims[1] = tetexteriorsize;
 if( 0 != mesh_info.bdy_tets.size() &&
 NC_NOERR != nc_def_var( ncFile, "tetrahedron_exterior", NC_LONG, 2, dims, &dum_var ) )
 {
 MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to create connectivity array for exterior tets" );
 }
 
 /* Node coordinate arrays: */
 
 dims[0] = ncoords;
 dims[1] = coord_size;
 if( NC_NOERR != nc_def_var( ncFile, "coords", NC_DOUBLE, 2, dims, &dum_var ) )
 {
 MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to define node coordinate array" );
 }
 
 return MB_SUCCESS;
}
 
ErrorCode WriteSLAC::open_file( const char* filename )
{
 // Not a valid filname
 if( strlen( (const char*)filename ) == 0 )
 {
 MB_SET_ERR( MB_FAILURE, "Output filename not specified" );
 }
 
 int fail = nc_create( filename, NC_CLOBBER, &ncFile );
 // File couldn't be opened
 if( NC_NOERR != fail )
 {
 MB_SET_ERR( MB_FAILURE, "Cannot open " << filename );
 }
 
 return MB_SUCCESS;
}
 
ErrorCode WriteSLAC::get_neuset_elems( EntityHandle neuset,
 int current_sense,
 Range& forward_elems,
 Range& reverse_elems )
{
 Range ss_elems, ss_meshsets;
 
 // Get the sense tag; don't need to check return, might be an error if the tag
 // hasn't been created yet
 Tag sense_tag = 0;
 mbImpl->tag_get_handle( "SENSE", 1, MB_TYPE_INTEGER, sense_tag );
 
 // Get the entities in this set
 ErrorCode result = mbImpl->get_entities_by_handle( neuset, ss_elems, true );
 if( MB_FAILURE == result ) return result;
 
 // Now remove the meshsets into the ss_meshsets; first find the first meshset,
 Range::iterator range_iter = ss_elems.begin();
 while( TYPE_FROM_HANDLE( *range_iter ) != MBENTITYSET && range_iter != ss_elems.end() )
 ++range_iter;
 
 // Then, if there are some, copy them into ss_meshsets and erase from ss_elems
 if( range_iter != ss_elems.end() )
 {
 std::copy( range_iter, ss_elems.end(), range_inserter( ss_meshsets ) );
 ss_elems.erase( range_iter, ss_elems.end() );
 }
 
 // OK, for the elements, check the sense of this set and copy into the right range
 // (if the sense is 0, copy into both ranges)
 
 // Need to step forward on list until we reach the right dimension
 Range::iterator dum_it = ss_elems.end();
 --dum_it;
 int target_dim = CN::Dimension( TYPE_FROM_HANDLE( *dum_it ) );
 dum_it = ss_elems.begin();
 while( target_dim != CN::Dimension( TYPE_FROM_HANDLE( *dum_it ) ) && dum_it != ss_elems.end() )
 ++dum_it;
 
 if( current_sense == 1 || current_sense == 0 ) std::copy( dum_it, ss_elems.end(), range_inserter( forward_elems ) );
 if( current_sense == -1 || current_sense == 0 )
 std::copy( dum_it, ss_elems.end(), range_inserter( reverse_elems ) );
 
 // Now loop over the contained meshsets, getting the sense of those and calling this
 // function recursively
 for( range_iter = ss_meshsets.begin(); range_iter != ss_meshsets.end(); ++range_iter )
 {
 // First get the sense; if it's not there, by convention it's forward
 int this_sense;
 if( 0 == sense_tag || MB_FAILURE == mbImpl->tag_get_data( sense_tag, &( *range_iter ), 1, &this_sense ) )
 this_sense = 1;
 
 // Now get all the entities on this meshset, with the proper (possibly reversed) sense
 get_neuset_elems( *range_iter, this_sense * current_sense, forward_elems, reverse_elems );
 }
 
 return result;
}
 
} // namespace moab