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