Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
WriteSLAC.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 "WriteSLAC.hpp"
00025 
00026 #include <utility>
00027 #include <algorithm>
00028 #include <ctime>
00029 #include <string>
00030 #include <vector>
00031 #include <cstdio>
00032 #include <cstring>
00033 #include <iostream>
00034 #include <cassert>
00035 
00036 #include "netcdf.h"
00037 #include "moab/Interface.hpp"
00038 #include "moab/Range.hpp"
00039 #include "moab/CN.hpp"
00040 #include "Internals.hpp"
00041 #include "ExoIIUtil.hpp"
00042 #include "MBTagConventions.hpp"
00043 #include "moab/WriteUtilIface.hpp"
00044 
00045 #ifndef MOAB_HAVE_NETCDF
00046 #error Attempt to compile WriteSLAC with NetCDF disabled.
00047 #endif
00048 
00049 namespace moab
00050 {
00051 
00052 #define INS_ID( stringvar, prefix, id ) sprintf( stringvar, prefix, id )
00053 
00054 #define GET_VAR( name, id, dims )                                     \
00055     {                                                                 \
00056         ( id )     = -1;                                              \
00057         int gvfail = nc_inq_varid( ncFile, name, &( id ) );           \
00058         if( NC_NOERR == gvfail )                                      \
00059         {                                                             \
00060             int ndims;                                                \
00061             gvfail = nc_inq_varndims( ncFile, id, &ndims );           \
00062             if( NC_NOERR == gvfail )                                  \
00063             {                                                         \
00064                 ( dims ).resize( ndims );                             \
00065                 gvfail = nc_inq_vardimid( ncFile, id, &( dims )[0] ); \
00066             }                                                         \
00067         }                                                             \
00068     }
00069 
00070 WriterIface* WriteSLAC::factory( Interface* iface )
00071 {
00072     return new WriteSLAC( iface );
00073 }
00074 
00075 WriteSLAC::WriteSLAC( Interface* impl ) : mbImpl( impl ), ncFile( 0 )
00076 {
00077     assert( impl != NULL );
00078 
00079     impl->query_interface( mWriteIface );
00080 
00081     // Initialize in case tag_get_handle fails below
00082     //! get and cache predefined tag handles
00083     int negone = -1;
00084     impl->tag_get_handle( MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mMaterialSetTag, MB_TAG_SPARSE | MB_TAG_CREAT,
00085                           &negone );
00086 
00087     impl->tag_get_handle( DIRICHLET_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mDirichletSetTag, MB_TAG_SPARSE | MB_TAG_CREAT,
00088                           &negone );
00089 
00090     impl->tag_get_handle( NEUMANN_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mNeumannSetTag, MB_TAG_SPARSE | MB_TAG_CREAT,
00091                           &negone );
00092 
00093     mGlobalIdTag = impl->globalId_tag();
00094 
00095     int dum_val = -1;
00096     impl->tag_get_handle( "__matSetIdTag", 1, MB_TYPE_INTEGER, mMatSetIdTag, MB_TAG_DENSE | MB_TAG_CREAT, &dum_val );
00097 
00098     impl->tag_get_handle( "WriteSLAC element mark", 1, MB_TYPE_BIT, mEntityMark, MB_TAG_CREAT );
00099 }
00100 
00101 WriteSLAC::~WriteSLAC()
00102 {
00103     mbImpl->release_interface( mWriteIface );
00104     mbImpl->tag_delete( mEntityMark );
00105 }
00106 
00107 void WriteSLAC::reset_matset( std::vector< WriteSLAC::MaterialSetData >& matset_info )
00108 {
00109     std::vector< WriteSLAC::MaterialSetData >::iterator iter;
00110 
00111     for( iter = matset_info.begin(); iter != matset_info.end(); ++iter )
00112         delete( *iter ).elements;
00113 }
00114 
00115 ErrorCode WriteSLAC::write_file( const char* file_name,
00116                                  const bool overwrite,
00117                                  const FileOptions&,
00118                                  const EntityHandle* ent_handles,
00119                                  const int num_sets,
00120                                  const std::vector< std::string >& /* qa_list */,
00121                                  const Tag* /* tag_list */,
00122                                  int /* num_tags */,
00123                                  int /* export_dimension */ )
00124 {
00125     assert( 0 != mMaterialSetTag && 0 != mNeumannSetTag && 0 != mDirichletSetTag );
00126 
00127     // Check the file name
00128     if( NULL == strstr( file_name, ".ncdf" ) ) return MB_FAILURE;
00129 
00130     std::vector< EntityHandle > matsets, dirsets, neusets, entities;
00131 
00132     fileName = file_name;
00133 
00134     // Separate into material sets, dirichlet sets, neumann sets
00135 
00136     if( num_sets == 0 )
00137     {
00138         // Default to all defined sets
00139         Range this_range;
00140         mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &mMaterialSetTag, NULL, 1, this_range );
00141         std::copy( this_range.begin(), this_range.end(), std::back_inserter( matsets ) );
00142         this_range.clear();
00143         mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &mDirichletSetTag, NULL, 1, this_range );
00144         std::copy( this_range.begin(), this_range.end(), std::back_inserter( dirsets ) );
00145         this_range.clear();
00146         mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &mNeumannSetTag, NULL, 1, this_range );
00147         std::copy( this_range.begin(), this_range.end(), std::back_inserter( neusets ) );
00148     }
00149     else
00150     {
00151         int dummy;
00152         for( const EntityHandle* iter = ent_handles; iter < ent_handles + num_sets; ++iter )
00153         {
00154             if( MB_SUCCESS == mbImpl->tag_get_data( mMaterialSetTag, &( *iter ), 1, &dummy ) )
00155                 matsets.push_back( *iter );
00156             else if( MB_SUCCESS == mbImpl->tag_get_data( mDirichletSetTag, &( *iter ), 1, &dummy ) )
00157                 dirsets.push_back( *iter );
00158             else if( MB_SUCCESS == mbImpl->tag_get_data( mNeumannSetTag, &( *iter ), 1, &dummy ) )
00159                 neusets.push_back( *iter );
00160         }
00161     }
00162 
00163     // If there is nothing to write just return.
00164     if( matsets.empty() && dirsets.empty() && neusets.empty() ) return MB_FILE_WRITE_ERROR;
00165 
00166     std::vector< WriteSLAC::MaterialSetData > matset_info;
00167     std::vector< WriteSLAC::DirichletSetData > dirset_info;
00168     std::vector< WriteSLAC::NeumannSetData > neuset_info;
00169 
00170     MeshInfo mesh_info;
00171 
00172     matset_info.clear();
00173     if( gather_mesh_information( mesh_info, matset_info, neuset_info, dirset_info, matsets, neusets, dirsets ) !=
00174         MB_SUCCESS )
00175     {
00176         reset_matset( matset_info );
00177         return MB_FAILURE;
00178     }
00179 
00180     // Try to open the file after gather mesh info succeeds
00181     int fail = nc_create( file_name, overwrite ? NC_CLOBBER : NC_NOCLOBBER, &ncFile );
00182     if( NC_NOERR != fail )
00183     {
00184         reset_matset( matset_info );
00185         return MB_FAILURE;
00186     }
00187 
00188     if( initialize_file( mesh_info ) != MB_SUCCESS )
00189     {
00190         reset_matset( matset_info );
00191         return MB_FAILURE;
00192     }
00193 
00194     if( write_nodes( mesh_info.num_nodes, mesh_info.nodes, mesh_info.num_dim ) != MB_SUCCESS )
00195     {
00196         reset_matset( matset_info );
00197         return MB_FAILURE;
00198     }
00199 
00200     if( write_matsets( mesh_info, matset_info, neuset_info ) )
00201     {
00202         reset_matset( matset_info );
00203         return MB_FAILURE;
00204     }
00205 
00206     fail = nc_close( ncFile );
00207     if( NC_NOERR != fail ) return MB_FAILURE;
00208 
00209     return MB_SUCCESS;
00210 }
00211 
00212 ErrorCode WriteSLAC::gather_mesh_information( MeshInfo& mesh_info,
00213                                               std::vector< WriteSLAC::MaterialSetData >& matset_info,
00214                                               std::vector< WriteSLAC::NeumannSetData >& neuset_info,
00215                                               std::vector< WriteSLAC::DirichletSetData >& dirset_info,
00216                                               std::vector< EntityHandle >& matsets,
00217                                               std::vector< EntityHandle >& neusets,
00218                                               std::vector< EntityHandle >& dirsets )
00219 {
00220     std::vector< EntityHandle >::iterator vector_iter, end_vector_iter;
00221 
00222     mesh_info.num_nodes    = 0;
00223     mesh_info.num_elements = 0;
00224     mesh_info.num_matsets  = 0;
00225 
00226     int id = 0;
00227 
00228     vector_iter     = matsets.begin();
00229     end_vector_iter = matsets.end();
00230 
00231     mesh_info.num_matsets = matsets.size();
00232 
00233     std::vector< EntityHandle > parent_meshsets;
00234 
00235     // Clean out the bits for the element mark
00236     mbImpl->tag_delete( mEntityMark );
00237     mbImpl->tag_get_handle( "WriteSLAC element mark", 1, MB_TYPE_BIT, mEntityMark, MB_TAG_CREAT );
00238 
00239     int highest_dimension_of_element_matsets = 0;
00240 
00241     for( vector_iter = matsets.begin(); vector_iter != matsets.end(); ++vector_iter )
00242     {
00243         WriteSLAC::MaterialSetData matset_data;
00244         matset_data.elements = new Range;
00245 
00246         // For the purpose of qa records, get the parents of these matsets
00247         if( mbImpl->get_parent_meshsets( *vector_iter, parent_meshsets ) != MB_SUCCESS ) return MB_FAILURE;
00248 
00249         // Get all Entity Handles in the mesh set
00250         Range dummy_range;
00251         mbImpl->get_entities_by_handle( *vector_iter, dummy_range, true );
00252 
00253         // Wait a minute, we are doing some filtering here that doesn't make sense at this level CJS
00254 
00255         // Find the dimension of the last entity in this range
00256         Range::iterator entity_iter = dummy_range.end();
00257         entity_iter                 = dummy_range.end();
00258         --entity_iter;
00259         int this_dim = CN::Dimension( TYPE_FROM_HANDLE( *entity_iter ) );
00260         entity_iter  = dummy_range.begin();
00261         while( entity_iter != dummy_range.end() && CN::Dimension( TYPE_FROM_HANDLE( *entity_iter ) ) != this_dim )
00262             ++entity_iter;
00263 
00264         if( entity_iter != dummy_range.end() )
00265             std::copy( entity_iter, dummy_range.end(), range_inserter( *( matset_data.elements ) ) );
00266 
00267         assert( matset_data.elements->begin() == matset_data.elements->end() ||
00268                 CN::Dimension( TYPE_FROM_HANDLE( *( matset_data.elements->begin() ) ) ) == this_dim );
00269 
00270         // Get the matset's id
00271         if( mbImpl->tag_get_data( mMaterialSetTag, &( *vector_iter ), 1, &id ) != MB_SUCCESS )
00272         {
00273             MB_SET_ERR( MB_FAILURE, "Couldn't get matset id from a tag for an element matset" );
00274         }
00275 
00276         matset_data.id                = id;
00277         matset_data.number_attributes = 0;
00278 
00279         // Iterate through all the elements in the meshset
00280         Range::iterator elem_range_iter, end_elem_range_iter;
00281         elem_range_iter     = matset_data.elements->begin();
00282         end_elem_range_iter = matset_data.elements->end();
00283 
00284         // Get the entity type for this matset, verifying that it's the same for all elements
00285         // THIS ASSUMES HANDLES SORT BY TYPE!!!
00286         EntityType entity_type = TYPE_FROM_HANDLE( *elem_range_iter );
00287         --end_elem_range_iter;
00288         if( entity_type != TYPE_FROM_HANDLE( *( end_elem_range_iter++ ) ) )
00289         {
00290             MB_SET_ERR( MB_FAILURE, "Entities in matset " << id << " not of common type" );
00291         }
00292 
00293         int dimension = -1;
00294         if( entity_type == MBQUAD || entity_type == MBTRI )
00295             dimension = 3;  // Output shells by default
00296         else if( entity_type == MBEDGE )
00297             dimension = 2;
00298         else
00299             dimension = CN::Dimension( entity_type );
00300 
00301         if( dimension > highest_dimension_of_element_matsets ) highest_dimension_of_element_matsets = dimension;
00302 
00303         matset_data.moab_type = mbImpl->type_from_handle( *( matset_data.elements->begin() ) );
00304         if( MBMAXTYPE == matset_data.moab_type ) return MB_FAILURE;
00305 
00306         std::vector< EntityHandle > tmp_conn;
00307         mbImpl->get_connectivity( &( *( matset_data.elements->begin() ) ), 1, tmp_conn );
00308         matset_data.element_type =
00309             ExoIIUtil::get_element_type_from_num_verts( tmp_conn.size(), entity_type, dimension );
00310 
00311         if( matset_data.element_type == EXOII_MAX_ELEM_TYPE )
00312         {
00313             MB_SET_ERR( MB_FAILURE, "Element type in matset " << id << " didn't get set correctly" );
00314         }
00315 
00316         matset_data.number_nodes_per_element = ExoIIUtil::VerticesPerElement[matset_data.element_type];
00317 
00318         // Number of nodes for this matset
00319         matset_data.number_elements = matset_data.elements->size();
00320 
00321         // Total number of elements
00322         mesh_info.num_elements += matset_data.number_elements;
00323 
00324         // Get the nodes for the elements
00325         mWriteIface->gather_nodes_from_elements( *matset_data.elements, mEntityMark, mesh_info.nodes );
00326 
00327         if( !neusets.empty() )
00328         {
00329             // If there are neusets, keep track of which elements are being written out
00330             for( Range::iterator iter = matset_data.elements->begin(); iter != matset_data.elements->end(); ++iter )
00331             {
00332                 unsigned char bit = 0x1;
00333                 mbImpl->tag_set_data( mEntityMark, &( *iter ), 1, &bit );
00334             }
00335         }
00336 
00337         matset_info.push_back( matset_data );
00338     }
00339 
00340     // If user hasn't entered dimension, we figure it out
00341     if( mesh_info.num_dim == 0 )
00342     {
00343         // Never want 1 or zero dimensions
00344         if( highest_dimension_of_element_matsets < 2 )
00345             mesh_info.num_dim = 3;
00346         else
00347             mesh_info.num_dim = highest_dimension_of_element_matsets;
00348     }
00349 
00350     Range::iterator range_iter, end_range_iter;
00351     range_iter     = mesh_info.nodes.begin();
00352     end_range_iter = mesh_info.nodes.end();
00353 
00354     mesh_info.num_nodes = mesh_info.nodes.size();
00355 
00356     //------dirsets--------
00357 
00358     vector_iter     = dirsets.begin();
00359     end_vector_iter = dirsets.end();
00360 
00361     for( ; vector_iter != end_vector_iter; ++vector_iter )
00362     {
00363         WriteSLAC::DirichletSetData dirset_data;
00364         dirset_data.id           = 0;
00365         dirset_data.number_nodes = 0;
00366 
00367         // Get the dirset's id
00368         if( mbImpl->tag_get_data( mDirichletSetTag, &( *vector_iter ), 1, &id ) != MB_SUCCESS )
00369         {
00370             MB_SET_ERR( MB_FAILURE, "Couldn't get id tag for dirset " << id );
00371         }
00372 
00373         dirset_data.id = id;
00374 
00375         std::vector< EntityHandle > node_vector;
00376         // Get the nodes of the dirset that are in mesh_info.nodes
00377         if( mbImpl->get_entities_by_handle( *vector_iter, node_vector, true ) != MB_SUCCESS )
00378         {
00379             MB_SET_ERR( MB_FAILURE, "Couldn't get nodes in dirset " << id );
00380         }
00381 
00382         std::vector< EntityHandle >::iterator iter, end_iter;
00383         iter     = node_vector.begin();
00384         end_iter = node_vector.end();
00385 
00386         int j                     = 0;
00387         unsigned char node_marked = 0;
00388         ErrorCode result;
00389         for( ; iter != end_iter; ++iter )
00390         {
00391             if( TYPE_FROM_HANDLE( *iter ) != MBVERTEX ) continue;
00392             result = mbImpl->tag_get_data( mEntityMark, &( *iter ), 1, &node_marked );MB_CHK_SET_ERR( result, "Couldn't get mark data" );
00393 
00394             if( 0x1 == node_marked ) dirset_data.nodes.push_back( *iter );
00395             j++;
00396         }
00397 
00398         dirset_data.number_nodes = dirset_data.nodes.size();
00399         dirset_info.push_back( dirset_data );
00400     }
00401 
00402     //------neusets--------
00403     vector_iter     = neusets.begin();
00404     end_vector_iter = neusets.end();
00405 
00406     for( ; vector_iter != end_vector_iter; ++vector_iter )
00407     {
00408         WriteSLAC::NeumannSetData neuset_data;
00409 
00410         // Get the neuset's id
00411         if( mbImpl->tag_get_data( mNeumannSetTag, &( *vector_iter ), 1, &id ) != MB_SUCCESS ) return MB_FAILURE;
00412 
00413         neuset_data.id              = id;
00414         neuset_data.mesh_set_handle = *vector_iter;
00415 
00416         // Get the sides in two lists, one forward the other reverse; starts with forward sense
00417         // by convention
00418         Range forward_elems, reverse_elems;
00419         if( get_neuset_elems( *vector_iter, 0, forward_elems, reverse_elems ) == MB_FAILURE ) return MB_FAILURE;
00420 
00421         ErrorCode result = get_valid_sides( forward_elems, 1, neuset_data );MB_CHK_SET_ERR( result, "Couldn't get valid sides data" );
00422         result = get_valid_sides( reverse_elems, -1, neuset_data );MB_CHK_SET_ERR( result, "Couldn't get valid sides data" );
00423 
00424         neuset_data.number_elements = neuset_data.elements.size();
00425         neuset_info.push_back( neuset_data );
00426     }
00427 
00428     // Get information about interior/exterior tets/hexes, and mark matset ids
00429     return gather_interior_exterior( mesh_info, matset_info, neuset_info );
00430 }
00431 
00432 ErrorCode WriteSLAC::get_valid_sides( Range& elems, const int sense, WriteSLAC::NeumannSetData& neuset_data )
00433 {
00434     // This is where we see if underlying element of side set element is included in output
00435 
00436     unsigned char element_marked = 0;
00437     ErrorCode result;
00438     for( Range::iterator iter = elems.begin(); iter != elems.end(); ++iter )
00439     {
00440         // Should insert here if "side" is a quad/tri on a quad/tri mesh
00441         result = mbImpl->tag_get_data( mEntityMark, &( *iter ), 1, &element_marked );MB_CHK_SET_ERR( result, "Couldn't get mark data" );
00442 
00443         if( 0x1 == element_marked )
00444         {
00445             neuset_data.elements.push_back( *iter );
00446 
00447             // TJT TODO: the sense should really be # edges + 1or2
00448             neuset_data.side_numbers.push_back( ( sense == 1 ? 1 : 2 ) );
00449         }
00450         else
00451         {  // Then "side" is probably a quad/tri on a hex/tet mesh
00452             std::vector< EntityHandle > parents;
00453             int dimension = CN::Dimension( TYPE_FROM_HANDLE( *iter ) );
00454 
00455             // Get the adjacent parent element of "side"
00456             if( mbImpl->get_adjacencies( &( *iter ), 1, dimension + 1, false, parents ) != MB_SUCCESS )
00457             {
00458                 MB_SET_ERR( MB_FAILURE, "Couldn't get adjacencies for neuset" );
00459             }
00460 
00461             if( !parents.empty() )
00462             {
00463                 // Make sure the adjacent parent element will be output
00464                 for( unsigned int k = 0; k < parents.size(); k++ )
00465                 {
00466                     result = mbImpl->tag_get_data( mEntityMark, &( parents[k] ), 1, &element_marked );MB_CHK_SET_ERR( result, "Couldn't get mark data" );
00467 
00468                     int side_no, this_sense, this_offset;
00469                     if( 0x1 == element_marked &&
00470                         mbImpl->side_number( parents[k], *iter, side_no, this_sense, this_offset ) == MB_SUCCESS &&
00471                         this_sense == sense )
00472                     {
00473                         neuset_data.elements.push_back( parents[k] );
00474                         neuset_data.side_numbers.push_back( side_no + 1 );
00475                         break;
00476                     }
00477                 }
00478             }
00479             else
00480             {
00481                 MB_SET_ERR( MB_FAILURE, "No parent element exists for element in neuset " << neuset_data.id );
00482             }
00483         }
00484     }
00485 
00486     return MB_SUCCESS;
00487 }
00488 
00489 ErrorCode WriteSLAC::write_nodes( const int num_nodes, const Range& nodes, const int dimension )
00490 {
00491     // See if should transform coordinates
00492     ErrorCode result;
00493     Tag trans_tag;
00494     result                = mbImpl->tag_get_handle( MESH_TRANSFORM_TAG_NAME, 16, MB_TYPE_DOUBLE, trans_tag );
00495     bool transform_needed = true;
00496     if( result == MB_TAG_NOT_FOUND ) transform_needed = false;
00497 
00498     int num_coords_to_fill = transform_needed ? 3 : dimension;
00499 
00500     std::vector< double* > coord_arrays( 3 );
00501     coord_arrays[0] = new double[num_nodes];
00502     coord_arrays[1] = new double[num_nodes];
00503     coord_arrays[2] = NULL;
00504 
00505     if( num_coords_to_fill == 3 ) coord_arrays[2] = new double[num_nodes];
00506 
00507     result = mWriteIface->get_node_coords( dimension, num_nodes, nodes, mGlobalIdTag, 0, coord_arrays );
00508     if( result != MB_SUCCESS )
00509     {
00510         delete[] coord_arrays[0];
00511         delete[] coord_arrays[1];
00512         if( coord_arrays[2] ) delete[] coord_arrays[2];
00513         return result;
00514     }
00515 
00516     if( transform_needed )
00517     {
00518         double trans_matrix[16];
00519         const EntityHandle mesh = 0;
00520         result                  = mbImpl->tag_get_data( trans_tag, &mesh, 1, trans_matrix );MB_CHK_SET_ERR( result, "Couldn't get transform data" );
00521 
00522         for( int i = 0; i < num_nodes; i++ )
00523         {
00524             double vec1[3];
00525             double vec2[3];
00526 
00527             vec2[0] = coord_arrays[0][i];
00528             vec2[1] = coord_arrays[1][i];
00529             vec2[2] = coord_arrays[2][i];
00530 
00531             for( int row = 0; row < 3; row++ )
00532             {
00533                 vec1[row] = 0.0;
00534                 for( int col = 0; col < 3; col++ )
00535                     vec1[row] += ( trans_matrix[( row * 4 ) + col] * vec2[col] );
00536             }
00537 
00538             coord_arrays[0][i] = vec1[0];
00539             coord_arrays[1][i] = vec1[1];
00540             coord_arrays[2][i] = vec1[2];
00541         }
00542     }
00543 
00544     // Write the nodes
00545     int nc_var = -1;
00546     std::vector< int > dims;
00547     GET_VAR( "coords", nc_var, dims );
00548     if( -1 == nc_var ) return MB_FAILURE;
00549     size_t start[2] = { 0, 0 }, count[2] = { static_cast< size_t >( num_nodes ), 1 };
00550     int fail = nc_put_vara_double( ncFile, nc_var, start, count, coord_arrays[0] );
00551     if( NC_NOERR != fail ) return MB_FAILURE;
00552     start[1] = 1;
00553     fail     = nc_put_vara_double( ncFile, nc_var, start, count, coord_arrays[1] );
00554     if( NC_NOERR != fail ) return MB_FAILURE;
00555     start[1] = 2;
00556     fail     = nc_put_vara_double( ncFile, nc_var, start, count, coord_arrays[2] );
00557     if( NC_NOERR != fail ) return MB_FAILURE;
00558 
00559     delete[] coord_arrays[0];
00560     delete[] coord_arrays[1];
00561     if( coord_arrays[2] ) delete[] coord_arrays[2];
00562 
00563     return MB_SUCCESS;
00564 }
00565 
00566 ErrorCode WriteSLAC::gather_interior_exterior( MeshInfo& mesh_info,
00567                                                std::vector< WriteSLAC::MaterialSetData >& matset_data,
00568                                                std::vector< WriteSLAC::NeumannSetData >& neuset_data )
00569 {
00570     // Need to assign a tag with the matset id
00571     Tag matset_id_tag;
00572     unsigned int i;
00573     int dum = -1;
00574     ErrorCode result =
00575         mbImpl->tag_get_handle( "__matset_id", 4, MB_TYPE_INTEGER, matset_id_tag, MB_TAG_DENSE | MB_TAG_CREAT, &dum );
00576     if( MB_SUCCESS != result ) return result;
00577 
00578     Range::iterator rit;
00579     mesh_info.num_int_hexes = mesh_info.num_int_tets = 0;
00580 
00581     for( i = 0; i < matset_data.size(); i++ )
00582     {
00583         WriteSLAC::MaterialSetData matset = matset_data[i];
00584         if( matset.moab_type == MBHEX )
00585             mesh_info.num_int_hexes += matset.elements->size();
00586         else if( matset.moab_type == MBTET )
00587             mesh_info.num_int_tets += matset.elements->size();
00588         else
00589         {
00590             std::cout << "WriteSLAC doesn't support elements of type " << CN::EntityTypeName( matset.moab_type )
00591                       << std::endl;
00592             continue;
00593         }
00594 
00595         for( rit = matset.elements->begin(); rit != matset.elements->end(); ++rit )
00596         {
00597             result = mbImpl->tag_set_data( mMatSetIdTag, &( *rit ), 1, &( matset.id ) );
00598             if( MB_SUCCESS != result ) return result;
00599         }
00600     }
00601 
00602     // Now go through the neumann sets, pulling out the hexes with faces on the
00603     // boundary
00604     std::vector< EntityHandle >::iterator vit;
00605     for( i = 0; i < neuset_data.size(); i++ )
00606     {
00607         WriteSLAC::NeumannSetData neuset = neuset_data[i];
00608         for( vit = neuset.elements.begin(); vit != neuset.elements.end(); ++vit )
00609         {
00610             if( TYPE_FROM_HANDLE( *vit ) == MBHEX )
00611                 mesh_info.bdy_hexes.insert( *vit );
00612             else if( TYPE_FROM_HANDLE( *vit ) == MBTET )
00613                 mesh_info.bdy_tets.insert( *vit );
00614         }
00615     }
00616 
00617     // Now we have the number of bdy hexes and tets, we know how many interior ones
00618     // there are too
00619     mesh_info.num_int_hexes -= mesh_info.bdy_hexes.size();
00620     mesh_info.num_int_tets -= mesh_info.bdy_tets.size();
00621 
00622     return MB_SUCCESS;
00623 }
00624 
00625 ErrorCode WriteSLAC::write_matsets( MeshInfo& mesh_info,
00626                                     std::vector< WriteSLAC::MaterialSetData >& matset_data,
00627                                     std::vector< WriteSLAC::NeumannSetData >& neuset_data )
00628 {
00629     unsigned int i;
00630     std::vector< int > connect;
00631     const EntityHandle* connecth;
00632     int num_connecth;
00633     ErrorCode result;
00634 
00635     // First write the interior hexes
00636     int hex_conn = -1;
00637     std::vector< int > dims;
00638     if( mesh_info.bdy_hexes.size() != 0 || mesh_info.num_int_hexes != 0 )
00639     {
00640         GET_VAR( "hexahedron_interior", hex_conn, dims );
00641         if( -1 == hex_conn ) return MB_FAILURE;
00642     }
00643     connect.reserve( 13 );
00644     Range::iterator rit;
00645 
00646     int elem_num = 0;
00647     WriteSLAC::MaterialSetData matset;
00648     size_t start[2] = { 0, 0 }, count[2] = { 1, 1 };
00649     int fail;
00650     for( i = 0; i < matset_data.size(); i++ )
00651     {
00652         matset = matset_data[i];
00653         if( matset.moab_type != MBHEX ) continue;
00654 
00655         int id     = matset.id;
00656         connect[0] = id;
00657 
00658         for( rit = matset.elements->begin(); rit != matset.elements->end(); ++rit )
00659         {
00660             // Skip if it's on the bdy
00661             if( mesh_info.bdy_hexes.find( *rit ) != mesh_info.bdy_hexes.end() ) continue;
00662 
00663             // Get the connectivity of this element
00664             result = mbImpl->get_connectivity( *rit, connecth, num_connecth );
00665             if( MB_SUCCESS != result ) return result;
00666 
00667             // Get the vertex ids
00668             result = mbImpl->tag_get_data( mGlobalIdTag, connecth, num_connecth, &connect[1] );
00669             if( MB_SUCCESS != result ) return result;
00670 
00671             // Put the variable at the right position
00672             start[0] = elem_num++;
00673             count[1] = 9;
00674 
00675             // Write the data
00676             fail = nc_put_vara_int( ncFile, hex_conn, start, count, &connect[0] );
00677             if( NC_NOERR != fail ) return MB_FAILURE;
00678         }
00679     }
00680 
00681     int tet_conn = -1;
00682     if( mesh_info.bdy_tets.size() != 0 || mesh_info.num_int_tets != 0 )
00683     {
00684         GET_VAR( "tetrahedron_interior", tet_conn, dims );
00685         if( -1 == tet_conn ) return MB_FAILURE;
00686     }
00687 
00688     // Now the interior tets
00689     elem_num = 0;
00690     for( i = 0; i < matset_data.size(); i++ )
00691     {
00692         matset = matset_data[i];
00693         if( matset.moab_type != MBTET ) continue;
00694 
00695         int id     = matset.id;
00696         connect[0] = id;
00697         elem_num   = 0;
00698         for( rit = matset.elements->begin(); rit != matset.elements->end(); ++rit )
00699         {
00700             // Skip if it's on the bdy
00701             if( mesh_info.bdy_tets.find( *rit ) != mesh_info.bdy_tets.end() ) continue;
00702 
00703             // Get the connectivity of this element
00704             result = mbImpl->get_connectivity( *rit, connecth, num_connecth );
00705             if( MB_SUCCESS != result ) return result;
00706 
00707             // Get the vertex ids
00708             result = mbImpl->tag_get_data( mGlobalIdTag, connecth, num_connecth, &connect[1] );
00709             if( MB_SUCCESS != result ) return result;
00710 
00711             // Put the variable at the right position
00712             start[0] = elem_num++;
00713             count[1] = 5;
00714             fail     = nc_put_vara_int( ncFile, tet_conn, start, count, &connect[0] );
00715             // Write the data
00716             if( NC_NOERR != fail ) return MB_FAILURE;
00717         }
00718     }
00719 
00720     // Now the exterior hexes
00721     if( mesh_info.bdy_hexes.size() != 0 )
00722     {
00723         hex_conn = -1;
00724         GET_VAR( "hexahedron_exterior", hex_conn, dims );
00725         if( -1 == hex_conn ) return MB_FAILURE;
00726 
00727         connect.reserve( 15 );
00728         elem_num = 0;
00729 
00730         // Write the elements
00731         for( rit = mesh_info.bdy_hexes.begin(); rit != mesh_info.bdy_hexes.end(); ++rit )
00732         {
00733             // Get the material set for this hex
00734             result = mbImpl->tag_get_data( mMatSetIdTag, &( *rit ), 1, &connect[0] );
00735             if( MB_SUCCESS != result ) return result;
00736 
00737             // Get the connectivity of this element
00738             result = mbImpl->get_connectivity( *rit, connecth, num_connecth );
00739             if( MB_SUCCESS != result ) return result;
00740 
00741             // Get the vertex ids
00742             result = mbImpl->tag_get_data( mGlobalIdTag, connecth, num_connecth, &connect[1] );
00743             if( MB_SUCCESS != result ) return result;
00744 
00745             // Preset side numbers
00746             for( i = 9; i < 15; i++ )
00747                 connect[i] = -1;
00748 
00749             // Now write the side numbers
00750             for( i = 0; i < neuset_data.size(); i++ )
00751             {
00752                 std::vector< EntityHandle >::iterator vit =
00753                     std::find( neuset_data[i].elements.begin(), neuset_data[i].elements.end(), *rit );
00754                 while( vit != neuset_data[i].elements.end() )
00755                 {
00756                     // Have a side - get the side # and put in connect array
00757                     int side_no          = neuset_data[i].side_numbers[vit - neuset_data[i].elements.begin()];
00758                     connect[9 + side_no] = neuset_data[i].id;
00759                     ++vit;
00760                     vit = std::find( vit, neuset_data[i].elements.end(), *rit );
00761                 }
00762             }
00763 
00764             // Put the variable at the right position
00765             start[0] = elem_num++;
00766             count[1] = 15;
00767             fail     = nc_put_vara_int( ncFile, hex_conn, start, count, &connect[0] );
00768             // Write the data
00769             if( NC_NOERR != fail ) return MB_FAILURE;
00770         }
00771     }
00772 
00773     // Now the exterior tets
00774     if( mesh_info.bdy_tets.size() != 0 )
00775     {
00776         tet_conn = -1;
00777         GET_VAR( "tetrahedron_exterior", tet_conn, dims );
00778         if( -1 == tet_conn ) return MB_FAILURE;
00779 
00780         connect.reserve( 9 );
00781         elem_num = 0;
00782 
00783         // Write the elements
00784         for( rit = mesh_info.bdy_tets.begin(); rit != mesh_info.bdy_tets.end(); ++rit )
00785         {
00786             // Get the material set for this tet
00787             result = mbImpl->tag_get_data( mMatSetIdTag, &( *rit ), 1, &connect[0] );
00788             if( MB_SUCCESS != result ) return result;
00789 
00790             // Get the connectivity of this element
00791             result = mbImpl->get_connectivity( *rit, connecth, num_connecth );
00792             if( MB_SUCCESS != result ) return result;
00793 
00794             // Get the vertex ids
00795             result = mbImpl->tag_get_data( mGlobalIdTag, connecth, num_connecth, &connect[1] );
00796             if( MB_SUCCESS != result ) return result;
00797 
00798             // Preset side numbers
00799             for( i = 5; i < 9; i++ )
00800                 connect[i] = -1;
00801 
00802             // Now write the side numbers
00803             for( i = 0; i < neuset_data.size(); i++ )
00804             {
00805                 std::vector< EntityHandle >::iterator vit =
00806                     std::find( neuset_data[i].elements.begin(), neuset_data[i].elements.end(), *rit );
00807                 while( vit != neuset_data[i].elements.end() )
00808                 {
00809                     // Have a side - get the side # and put in connect array
00810                     int side_no          = neuset_data[i].side_numbers[vit - neuset_data[i].elements.begin()];
00811                     connect[5 + side_no] = neuset_data[i].id;
00812                     ++vit;
00813                     vit = std::find( vit, neuset_data[i].elements.end(), *rit );
00814                 }
00815             }
00816 
00817             // Put the variable at the right position
00818             start[0] = elem_num++;
00819             count[1] = 9;
00820             fail     = nc_put_vara_int( ncFile, tet_conn, start, count, &connect[0] );
00821             // Write the data
00822             if( NC_NOERR != fail ) return MB_FAILURE;
00823         }
00824     }
00825 
00826     return MB_SUCCESS;
00827 }
00828 
00829 ErrorCode WriteSLAC::initialize_file( MeshInfo& mesh_info )
00830 {
00831     // Perform the initializations
00832 
00833     int coord_size = -1, ncoords = -1;
00834     // Initialization to avoid warnings on Linux
00835     int hexinterior = -1, hexinteriorsize, hexexterior = -1, hexexteriorsize = -1;
00836     int tetinterior = -1, tetinteriorsize, tetexterior = -1, tetexteriorsize = -1;
00837 
00838     if( nc_def_dim( ncFile, "coord_size", (size_t)mesh_info.num_dim, &coord_size ) != NC_NOERR )
00839     {
00840         MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to define number of dimensions" );
00841     }
00842 
00843     if( nc_def_dim( ncFile, "ncoords", (size_t)mesh_info.num_nodes, &ncoords ) != NC_NOERR )
00844     {
00845         MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to define number of nodes" );
00846     }
00847 
00848     if( 0 != mesh_info.num_int_hexes &&
00849         nc_def_dim( ncFile, "hexinterior", (size_t)mesh_info.num_int_hexes, &hexinterior ) != NC_NOERR )
00850     {
00851         MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to define number of interior hex elements" );
00852     }
00853 
00854     if( nc_def_dim( ncFile, "hexinteriorsize", (size_t)9, &hexinteriorsize ) != NC_NOERR )
00855     {
00856         MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to define interior hex element size" );
00857     }
00858 
00859     if( 0 != mesh_info.bdy_hexes.size() &&
00860         nc_def_dim( ncFile, "hexexterior", (size_t)mesh_info.bdy_hexes.size(), &hexexterior ) != NC_NOERR )
00861     {
00862         MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to define number of exterior hex elements" );
00863     }
00864 
00865     if( nc_def_dim( ncFile, "hexexteriorsize", (size_t)15, &hexexteriorsize ) != NC_NOERR )
00866     {
00867         MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to define exterior hex element size" );
00868     }
00869 
00870     if( 0 != mesh_info.num_int_tets &&
00871         nc_def_dim( ncFile, "tetinterior", (size_t)mesh_info.num_int_tets, &tetinterior ) != NC_NOERR )
00872     {
00873         MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to define number of interior tet elements" );
00874     }
00875 
00876     if( nc_def_dim( ncFile, "tetinteriorsize", (size_t)5, &tetinteriorsize ) != NC_NOERR )
00877     {
00878         MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to define interior tet element size" );
00879     }
00880 
00881     if( 0 != mesh_info.bdy_tets.size() &&
00882         nc_def_dim( ncFile, "tetexterior", (size_t)mesh_info.bdy_tets.size(), &tetexterior ) != NC_NOERR )
00883     {
00884         MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to define number of exterior tet elements" );
00885     }
00886 
00887     if( nc_def_dim( ncFile, "tetexteriorsize", (size_t)9, &tetexteriorsize ) != NC_NOERR )
00888     {
00889         MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to define exterior tet element size" );
00890     }
00891 
00892     /* ...and some variables */
00893 
00894     int dims[2];
00895     dims[0] = hexinterior;
00896     dims[1] = hexinteriorsize;
00897     int dum_var;
00898     if( 0 != mesh_info.num_int_hexes &&
00899         NC_NOERR != nc_def_var( ncFile, "hexahedron_interior", NC_LONG, 2, dims, &dum_var ) )
00900     {
00901         MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to create connectivity array for interior hexes" );
00902     }
00903 
00904     dims[0] = hexexterior;
00905     dims[1] = hexexteriorsize;
00906     if( 0 != mesh_info.bdy_hexes.size() &&
00907         NC_NOERR != nc_def_var( ncFile, "hexahedron_exterior", NC_LONG, 2, dims, &dum_var ) )
00908     {
00909         MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to create connectivity array for exterior hexes" );
00910     }
00911 
00912     dims[0] = tetinterior;
00913     dims[1] = tetinteriorsize;
00914     if( 0 != mesh_info.num_int_tets &&
00915         NC_NOERR != nc_def_var( ncFile, "tetrahedron_exterior", NC_LONG, 2, dims, &dum_var ) )
00916     {
00917         MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to create connectivity array for interior tets" );
00918     }
00919 
00920     dims[0] = tetexterior;
00921     dims[1] = tetexteriorsize;
00922     if( 0 != mesh_info.bdy_tets.size() &&
00923         NC_NOERR != nc_def_var( ncFile, "tetrahedron_exterior", NC_LONG, 2, dims, &dum_var ) )
00924     {
00925         MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to create connectivity array for exterior tets" );
00926     }
00927 
00928     /* Node coordinate arrays: */
00929 
00930     dims[0] = ncoords;
00931     dims[1] = coord_size;
00932     if( NC_NOERR != nc_def_var( ncFile, "coords", NC_DOUBLE, 2, dims, &dum_var ) )
00933     {
00934         MB_SET_ERR( MB_FAILURE, "WriteSLAC: failed to define node coordinate array" );
00935     }
00936 
00937     return MB_SUCCESS;
00938 }
00939 
00940 ErrorCode WriteSLAC::open_file( const char* filename )
00941 {
00942     // Not a valid filname
00943     if( strlen( (const char*)filename ) == 0 )
00944     {
00945         MB_SET_ERR( MB_FAILURE, "Output filename not specified" );
00946     }
00947 
00948     int fail = nc_create( filename, NC_CLOBBER, &ncFile );
00949     // File couldn't be opened
00950     if( NC_NOERR != fail )
00951     {
00952         MB_SET_ERR( MB_FAILURE, "Cannot open " << filename );
00953     }
00954 
00955     return MB_SUCCESS;
00956 }
00957 
00958 ErrorCode WriteSLAC::get_neuset_elems( EntityHandle neuset,
00959                                        int current_sense,
00960                                        Range& forward_elems,
00961                                        Range& reverse_elems )
00962 {
00963     Range ss_elems, ss_meshsets;
00964 
00965     // Get the sense tag; don't need to check return, might be an error if the tag
00966     // hasn't been created yet
00967     Tag sense_tag = 0;
00968     mbImpl->tag_get_handle( "SENSE", 1, MB_TYPE_INTEGER, sense_tag );
00969 
00970     // Get the entities in this set
00971     ErrorCode result = mbImpl->get_entities_by_handle( neuset, ss_elems, true );
00972     if( MB_FAILURE == result ) return result;
00973 
00974     // Now remove the meshsets into the ss_meshsets; first find the first meshset,
00975     Range::iterator range_iter = ss_elems.begin();
00976     while( TYPE_FROM_HANDLE( *range_iter ) != MBENTITYSET && range_iter != ss_elems.end() )
00977         ++range_iter;
00978 
00979     // Then, if there are some, copy them into ss_meshsets and erase from ss_elems
00980     if( range_iter != ss_elems.end() )
00981     {
00982         std::copy( range_iter, ss_elems.end(), range_inserter( ss_meshsets ) );
00983         ss_elems.erase( range_iter, ss_elems.end() );
00984     }
00985 
00986     // OK, for the elements, check the sense of this set and copy into the right range
00987     // (if the sense is 0, copy into both ranges)
00988 
00989     // Need to step forward on list until we reach the right dimension
00990     Range::iterator dum_it = ss_elems.end();
00991     --dum_it;
00992     int target_dim = CN::Dimension( TYPE_FROM_HANDLE( *dum_it ) );
00993     dum_it         = ss_elems.begin();
00994     while( target_dim != CN::Dimension( TYPE_FROM_HANDLE( *dum_it ) ) && dum_it != ss_elems.end() )
00995         ++dum_it;
00996 
00997     if( current_sense == 1 || current_sense == 0 ) std::copy( dum_it, ss_elems.end(), range_inserter( forward_elems ) );
00998     if( current_sense == -1 || current_sense == 0 )
00999         std::copy( dum_it, ss_elems.end(), range_inserter( reverse_elems ) );
01000 
01001     // Now loop over the contained meshsets, getting the sense of those and calling this
01002     // function recursively
01003     for( range_iter = ss_meshsets.begin(); range_iter != ss_meshsets.end(); ++range_iter )
01004     {
01005         // First get the sense; if it's not there, by convention it's forward
01006         int this_sense;
01007         if( 0 == sense_tag || MB_FAILURE == mbImpl->tag_get_data( sense_tag, &( *range_iter ), 1, &this_sense ) )
01008             this_sense = 1;
01009 
01010         // Now get all the entities on this meshset, with the proper (possibly reversed) sense
01011         get_neuset_elems( *range_iter, this_sense * current_sense, forward_elems, reverse_elems );
01012     }
01013 
01014     return result;
01015 }
01016 
01017 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines