MOAB: Mesh Oriented datABase  (version 5.4.1)
WriteCCMIO.cpp
Go to the documentation of this file.
00001 /*
00002  * CCMIO file structure
00003  *
00004  * Root
00005  *   State(kCCMIOState)
00006  *     Processor*
00007  *       VerticesID
00008  *       TopologyID
00009  *       InitialID
00010  *       SolutionID
00011  *   Vertices*
00012  *     ->WriteVerticesx, WriteMap
00013  *   Topology*
00014  *     Boundary faces*(kCCMIOBoundaryFaces)
00015  *        ->WriteFaces, WriteFaceCells, WriteMap
00016  *     Internal faces(kCCMIOInternalFaces)
00017  *     Cells (kCCMIOCells)
00018  *        ->WriteCells (mapID), WriteMap, WriteCells
00019  *   Solution
00020  *     Phase
00021  *       Field
00022  *         FieldData
00023  *   Problem(kCCMIOProblemDescription)
00024  *     CellType* (kCCMIOCellType)
00025  *       Index (GetEntityIndex), MaterialId(WriteOpti), MaterialType(WriteOptstr),
00026  *         PorosityId(WriteOpti), SpinId(WriteOpti), GroupId(WriteOpti)
00027  *
00028  * MaterialType (CCMIOWriteOptstr in readexample)
00029  * constants (see readexample)
00030  * lagrangian data (CCMIOWriteLagrangianData)
00031  * vertices label (CCMIOEntityDescription)
00032  * restart info: char solver[], iterations, time, char timeUnits[], angle
00033  *      (CCMIOWriteRestartInfo, kCCMIORestartData), reference data?
00034  * phase:
00035  *   field: char name[], dims, CCMIODataType datatype, char units[]
00036  *       dims = kCCMIOScalar (CCMIOWriteFieldDataf),
00037  *              kCCMIOVector (CCMIOWriteMultiDimensionalFieldData),
00038  *              kCCMIOTensor
00039  * MonitoringSets: num, name (CellSet, VertexSet, BoundarySet, BlockSet, SplineSet, CoupleSet)
00040  *      CCMIOGetProstarSet, CCMIOWriteOpt1i,
00041  */
00042 
00043 #ifdef WIN32
00044 #ifdef _DEBUG
00045 // turn off warnings that say they debugging identifier has been truncated
00046 // this warning comes up when using some STL containers
00047 #pragma warning( disable : 4786 )
00048 #endif
00049 #endif
00050 
00051 #include "WriteCCMIO.hpp"
00052 #include "ccmio.h"
00053 #include "ccmioutility.h"
00054 #include "ccmiocore.h"
00055 #include <utility>
00056 #include <algorithm>
00057 #include <ctime>
00058 #include <string>
00059 #include <vector>
00060 #include <cstdio>
00061 #include <iostream>
00062 #include <algorithm>
00063 #include <sstream>
00064 
00065 #include "moab/Interface.hpp"
00066 #include "moab/Range.hpp"
00067 #include "moab/CN.hpp"
00068 #include "moab/Skinner.hpp"
00069 #include <cassert>
00070 #include "Internals.hpp"
00071 #include "ExoIIUtil.hpp"
00072 #include "MBTagConventions.hpp"
00073 #ifdef MOAB_HAVE_MPI
00074 #include "MBParallelConventions.h"
00075 #endif
00076 #include "moab/WriteUtilIface.hpp"
00077 
00078 namespace moab
00079 {
00080 
00081 static char const kStateName[] = "default";
00082 
00083 /*
00084   static const int ccm_types[] = {
00085     1,   // MBVERTEX
00086     2,   // MBEDGE
00087     -1,  // MBTRI
00088     -1,  // MBQUAD
00089     -1,  // MBPOLYGON
00090     13,  // MBTET
00091     14,  // MBPYRAMID
00092     12,  // MBPRISM
00093     -1,  // MBKNIFE
00094     11,  // MBHEX
00095     255  // MBPOLYHEDRON
00096   };
00097 */
00098 
00099 #define INS_ID( stringvar, prefix, id ) sprintf( stringvar, prefix, id )
00100 
00101 #define CHK_SET_CCMERR( ccm_err_code, ccm_err_msg )                                  \
00102     {                                                                                \
00103         if( kCCMIONoErr != ( ccm_err_code ) ) MB_SET_ERR( MB_FAILURE, ccm_err_msg ); \
00104     }
00105 
00106 WriterIface* WriteCCMIO::factory( Interface* iface )
00107 {
00108     return new WriteCCMIO( iface );
00109 }
00110 
00111 WriteCCMIO::WriteCCMIO( Interface* impl )
00112     : mbImpl( impl ), mCurrentMeshHandle( 0 ), mPartitionSetTag( 0 ), mNameTag( 0 ), mMaterialIdTag( 0 ),
00113       mMaterialTypeTag( 0 ), mRadiationTag( 0 ), mPorosityIdTag( 0 ), mSpinIdTag( 0 ), mGroupIdTag( 0 ),
00114       mColorIdxTag( 0 ), mProcessorIdTag( 0 ), mLightMaterialTag( 0 ), mFreeSurfaceMaterialTag( 0 ), mThicknessTag( 0 ),
00115       mProstarRegionNumberTag( 0 ), mBoundaryTypeTag( 0 ), mCreatingProgramTag( 0 ), mDimension( 0 ),
00116       mWholeMesh( false )
00117 {
00118     assert( impl != NULL );
00119 
00120     impl->query_interface( mWriteIface );
00121 
00122     // Initialize in case tag_get_handle fails below
00123     //! Get and cache predefined tag handles
00124     int negone = -1;
00125     impl->tag_get_handle( MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mMaterialSetTag, MB_TAG_SPARSE | MB_TAG_CREAT,
00126                           &negone );
00127 
00128     impl->tag_get_handle( DIRICHLET_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mDirichletSetTag, MB_TAG_SPARSE | MB_TAG_CREAT,
00129                           &negone );
00130 
00131     impl->tag_get_handle( NEUMANN_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mNeumannSetTag, MB_TAG_SPARSE | MB_TAG_CREAT,
00132                           &negone );
00133 
00134     mGlobalIdTag = impl->globalId_tag();
00135 
00136 #ifdef MOAB_HAVE_MPI
00137     impl->tag_get_handle( PARALLEL_PARTITION_TAG_NAME, 1, MB_TYPE_INTEGER, mPartitionSetTag, MB_TAG_SPARSE );
00138     // No need to check result, if it's not there, we don't create one
00139 #endif
00140 
00141     int dum_val_array[] = { -1, -1, -1, -1 };
00142     impl->tag_get_handle( HAS_MID_NODES_TAG_NAME, 4, MB_TYPE_INTEGER, mHasMidNodesTag, MB_TAG_SPARSE | MB_TAG_CREAT,
00143                           dum_val_array );
00144 
00145     impl->tag_get_handle( "__WriteCCMIO element mark", 1, MB_TYPE_BIT, mEntityMark, MB_TAG_CREAT );
00146 
00147     // Don't need to check return of following, since it doesn't matter if there isn't one
00148     mbImpl->tag_get_handle( NAME_TAG_NAME, NAME_TAG_SIZE, MB_TYPE_OPAQUE, mNameTag );
00149 }
00150 
00151 WriteCCMIO::~WriteCCMIO()
00152 {
00153     mbImpl->release_interface( mWriteIface );
00154     mbImpl->tag_delete( mEntityMark );
00155 }
00156 
00157 ErrorCode WriteCCMIO::write_file( const char* file_name,
00158                                   const bool overwrite,
00159                                   const FileOptions&,
00160                                   const EntityHandle* ent_handles,
00161                                   const int num_sets,
00162                                   const std::vector< std::string >& /* qa_list */,
00163                                   const Tag* /* tag_list */,
00164                                   int /* num_tags */,
00165                                   int /* export_dimension */ )
00166 {
00167     assert( 0 != mMaterialSetTag && 0 != mNeumannSetTag && 0 != mDirichletSetTag );
00168 
00169     ErrorCode result;
00170 
00171     // Check overwrite flag and file existence
00172     if( !overwrite )
00173     {
00174         FILE* file = fopen( file_name, "r" );
00175         if( file )
00176         {
00177             fclose( file );
00178             MB_SET_ERR( MB_FILE_WRITE_ERROR, "File exists but overwrite set to false" );
00179         }
00180     }
00181 
00182     mDimension = 3;
00183 
00184     std::vector< EntityHandle > matsets, dirsets, neusets, partsets;
00185 
00186     // Separate into material, dirichlet, neumann, partition sets
00187     result = get_sets( ent_handles, num_sets, matsets, dirsets, neusets, partsets );MB_CHK_SET_ERR( result, "Failed to get material/etc. sets" );
00188 
00189     // If entity handles were input but didn't contain matsets, return error
00190     if( ent_handles && matsets.empty() )
00191     {
00192         MB_SET_ERR( MB_FILE_WRITE_ERROR, "Sets input to write but no material sets found" );
00193     }
00194 
00195     // Otherwise, if no matsets, use root set
00196     if( matsets.empty() ) matsets.push_back( 0 );
00197 
00198     std::vector< MaterialSetData > matset_info;
00199     Range all_verts;
00200     result = gather_matset_info( matsets, matset_info, all_verts );MB_CHK_SET_ERR( result, "gathering matset info failed" );
00201 
00202     // Assign vertex gids
00203     result = mWriteIface->assign_ids( all_verts, mGlobalIdTag, 1 );MB_CHK_SET_ERR( result, "Failed to assign vertex global ids" );
00204 
00205     // Some CCMIO descriptors
00206     CCMIOID rootID, topologyID, stateID, problemID, verticesID, processorID;
00207 
00208     // Try to open the file and establish state
00209     result = open_file( file_name, overwrite, rootID );MB_CHK_SET_ERR( result, "Couldn't open file or create state" );
00210 
00211     result = create_ccmio_structure( rootID, stateID, processorID );MB_CHK_SET_ERR( result, "Problem creating CCMIO file structure" );
00212 
00213     result = write_nodes( rootID, all_verts, mDimension, verticesID );MB_CHK_SET_ERR( result, "write_nodes failed" );
00214 
00215     std::vector< NeumannSetData > neuset_info;
00216     result = gather_neuset_info( neusets, neuset_info );MB_CHK_SET_ERR( result, "Failed to get neumann set info" );
00217 
00218     result = write_cells_and_faces( rootID, matset_info, neuset_info, all_verts, topologyID );MB_CHK_SET_ERR( result, "write_cells_and_faces failed" );
00219 
00220     result = write_problem_description( rootID, stateID, problemID, processorID, matset_info, neuset_info );MB_CHK_SET_ERR( result, "write_problem_description failed" );
00221 
00222     result = write_solution_data();MB_CHK_SET_ERR( result, "Trouble writing solution data" );
00223 
00224     result = write_processor( processorID, verticesID, topologyID );MB_CHK_SET_ERR( result, "Trouble writing processor" );
00225 
00226     result = close_and_compress( file_name, rootID );MB_CHK_SET_ERR( result, "Close or compress failed" );
00227 
00228     return MB_SUCCESS;
00229 }
00230 
00231 ErrorCode WriteCCMIO::write_solution_data()
00232 {
00233     // For now, no solution (tag) data
00234     return MB_SUCCESS;
00235 }
00236 
00237 ErrorCode WriteCCMIO::write_processor( CCMIOID processorID, CCMIOID verticesID, CCMIOID topologyID )
00238 {
00239     CCMIOError error = kCCMIONoErr;
00240 
00241     // Now we have the mesh (vertices and topology) and the post data written.
00242     // Since we now have their IDs, we can write out the processor information.
00243     CCMIOWriteProcessor( &error, processorID, NULL, &verticesID, NULL, &topologyID, NULL, NULL, NULL, NULL );
00244     CHK_SET_CCMERR( error, "Problem writing CCMIO processor" );
00245 
00246     return MB_SUCCESS;
00247 }
00248 
00249 ErrorCode WriteCCMIO::create_ccmio_structure( CCMIOID rootID, CCMIOID& stateID, CCMIOID& processorID )
00250 {
00251     // Create problem state and other CCMIO nodes under it
00252     CCMIOError error = kCCMIONoErr;
00253 
00254     // Create a new state (or re-use an existing one).
00255     if( CCMIOGetState( NULL, rootID, kStateName, NULL, &stateID ) != kCCMIONoErr )
00256     {
00257         CCMIONewState( &error, rootID, kStateName, NULL, NULL, &stateID );
00258         CHK_SET_CCMERR( error, "Trouble creating state" );
00259     }
00260 
00261     // Create or get an old processor for this state
00262     CCMIOSize_t i = CCMIOSIZEC( 0 );
00263     if( CCMIONextEntity( NULL, stateID, kCCMIOProcessor, &i, &processorID ) != kCCMIONoErr )
00264     {
00265         CCMIONewEntity( &error, stateID, kCCMIOProcessor, NULL, &processorID );
00266         CHK_SET_CCMERR( error, "Trouble creating processor node" );
00267     }
00268     // Get rid of any data that may be in this processor (if the state was
00269     // not new).
00270     else
00271     {
00272         CCMIOClearProcessor( &error, stateID, processorID, TRUE, TRUE, TRUE, TRUE, TRUE );
00273         CHK_SET_CCMERR( error, "Trouble clearing processor data" );
00274     }
00275 
00276     /*
00277      //  for (; i < CCMIOSIZEC(partsets.size()); i++) {
00278      CCMIOSize_t id = CCMIOSIZEC(0);
00279      if (CCMIONextEntity(NULL, stateID, kCCMIOProcessor, &id, &processorID) != kCCMIONoErr)
00280      CCMIONewEntity(&error, stateID, kCCMIOProcessor, NULL, &processorID);
00281      CHKCCMERR(error, "Trouble creating processor node.");
00282      */
00283     return MB_SUCCESS;
00284 }
00285 
00286 ErrorCode WriteCCMIO::close_and_compress( const char*, CCMIOID rootID )
00287 {
00288     CCMIOError error = kCCMIONoErr;
00289     CCMIOCloseFile( &error, rootID );
00290     CHK_SET_CCMERR( error, "File close failed" );
00291 
00292     // The CCMIO library uses ADF to store the actual data.  Unfortunately,
00293     // ADF leaks disk space;  deleting a node does not recover all the disk
00294     // space.  Now that everything is successfully written it might be useful
00295     // to call CCMIOCompress() here to ensure that the file is as small as
00296     // possible.  Please see the Core API documentation for caveats on its
00297     // usage.
00298     // CCMIOCompress(&error, const_cast<char*>(filename));CHK_SET_CCMERR(error, "Error compressing
00299     // file");
00300 
00301     return MB_SUCCESS;
00302 }
00303 
00304 ErrorCode WriteCCMIO::open_file( const char* filename, bool, CCMIOID& rootID )
00305 {
00306     CCMIOError error = kCCMIONoErr;
00307     CCMIOOpenFile( &error, filename, kCCMIOWrite, &rootID );
00308     CHK_SET_CCMERR( error, "Cannot open file" );
00309 
00310     return MB_SUCCESS;
00311 }
00312 
00313 ErrorCode WriteCCMIO::get_sets( const EntityHandle* ent_handles,
00314                                 int num_sets,
00315                                 std::vector< EntityHandle >& matsets,
00316                                 std::vector< EntityHandle >& dirsets,
00317                                 std::vector< EntityHandle >& neusets,
00318                                 std::vector< EntityHandle >& partsets )
00319 {
00320     if( num_sets == 0 )
00321     {
00322         // Default to all defined sets
00323         Range this_range;
00324         mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &mMaterialSetTag, NULL, 1, this_range );
00325         std::copy( this_range.begin(), this_range.end(), std::back_inserter( matsets ) );
00326         this_range.clear();
00327         mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &mDirichletSetTag, NULL, 1, this_range );
00328         std::copy( this_range.begin(), this_range.end(), std::back_inserter( dirsets ) );
00329         this_range.clear();
00330         mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &mNeumannSetTag, NULL, 1, this_range );
00331         std::copy( this_range.begin(), this_range.end(), std::back_inserter( neusets ) );
00332         if( mPartitionSetTag )
00333         {
00334             this_range.clear();
00335             mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &mPartitionSetTag, NULL, 1, this_range );
00336             std::copy( this_range.begin(), this_range.end(), std::back_inserter( partsets ) );
00337         }
00338     }
00339     else
00340     {
00341         int dummy;
00342         for( const EntityHandle* iter = ent_handles; iter < ent_handles + num_sets; ++iter )
00343         {
00344             if( MB_SUCCESS == mbImpl->tag_get_data( mMaterialSetTag, &( *iter ), 1, &dummy ) )
00345                 matsets.push_back( *iter );
00346             else if( MB_SUCCESS == mbImpl->tag_get_data( mDirichletSetTag, &( *iter ), 1, &dummy ) )
00347                 dirsets.push_back( *iter );
00348             else if( MB_SUCCESS == mbImpl->tag_get_data( mNeumannSetTag, &( *iter ), 1, &dummy ) )
00349                 neusets.push_back( *iter );
00350             else if( mPartitionSetTag && MB_SUCCESS == mbImpl->tag_get_data( mPartitionSetTag, &( *iter ), 1, &dummy ) )
00351                 partsets.push_back( *iter );
00352         }
00353     }
00354 
00355     return MB_SUCCESS;
00356 }
00357 
00358 ErrorCode WriteCCMIO::write_problem_description( CCMIOID rootID,
00359                                                  CCMIOID stateID,
00360                                                  CCMIOID& problemID,
00361                                                  CCMIOID processorID,
00362                                                  std::vector< WriteCCMIO::MaterialSetData >& matset_data,
00363                                                  std::vector< WriteCCMIO::NeumannSetData >& neuset_data )
00364 {
00365     // Write out a dummy problem description.  If we happen to know that
00366     // there already is a problem description previously recorded that
00367     // is valid we could skip this step.
00368     CCMIOID id;
00369     CCMIOError error = kCCMIONoErr;
00370     ErrorCode rval;
00371     const EntityHandle mesh = 0;
00372 
00373     bool root_tagged = false, other_set_tagged = false;
00374     Tag simname;
00375     Range dum_sets;
00376     rval = mbImpl->tag_get_handle( "Title", 0, MB_TYPE_OPAQUE, simname, MB_TAG_ANY );
00377     if( MB_SUCCESS == rval )
00378     {
00379         int tag_size;
00380         rval = mbImpl->tag_get_bytes( simname, tag_size );
00381         if( MB_SUCCESS == rval )
00382         {
00383             std::vector< char > title_tag( tag_size + 1 );
00384             rval = mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &simname, NULL, 1, dum_sets );
00385             if( MB_SUCCESS == rval && !dum_sets.empty() )
00386             {
00387                 rval = mbImpl->tag_get_data( simname, &( *dum_sets.begin() ), 1, &title_tag[0] );MB_CHK_SET_ERR( rval, "Problem getting simulation name tag" );
00388                 other_set_tagged = true;
00389             }
00390             else if( MB_SUCCESS == rval )
00391             {
00392                 // Check to see if interface was tagged
00393                 rval = mbImpl->tag_get_data( simname, &mesh, 1, &title_tag[0] );
00394                 if( MB_SUCCESS == rval )
00395                     root_tagged = true;
00396                 else
00397                     rval = MB_SUCCESS;
00398             }
00399             *title_tag.rbegin() = '\0';
00400             if( root_tagged || other_set_tagged )
00401             {
00402                 CCMIONode rootNode;
00403                 if( kCCMIONoErr == CCMIOGetEntityNode( &error, rootID, &rootNode ) )
00404                 {
00405                     CCMIOSetTitle( &error, rootNode, &title_tag[0] );
00406                     CHK_SET_CCMERR( error, "Trouble setting title" );
00407                 }
00408             }
00409         }
00410     }
00411 
00412     rval = mbImpl->tag_get_handle( "CreatingProgram", 0, MB_TYPE_OPAQUE, mCreatingProgramTag, MB_TAG_ANY );
00413     if( MB_SUCCESS == rval )
00414     {
00415         int tag_size;
00416         rval = mbImpl->tag_get_bytes( mCreatingProgramTag, tag_size );
00417         if( MB_SUCCESS == rval )
00418         {
00419             std::vector< char > cp_tag( tag_size + 1 );
00420             rval = mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET, &mCreatingProgramTag, NULL, 1, dum_sets );
00421             if( MB_SUCCESS == rval && !dum_sets.empty() )
00422             {
00423                 rval = mbImpl->tag_get_data( mCreatingProgramTag, &( *dum_sets.begin() ), 1, &cp_tag[0] );MB_CHK_SET_ERR( rval, "Problem getting creating program tag" );
00424                 other_set_tagged = true;
00425             }
00426             else if( MB_SUCCESS == rval )
00427             {
00428                 // Check to see if interface was tagged
00429                 rval = mbImpl->tag_get_data( mCreatingProgramTag, &mesh, 1, &cp_tag[0] );
00430                 if( MB_SUCCESS == rval )
00431                     root_tagged = true;
00432                 else
00433                     rval = MB_SUCCESS;
00434             }
00435             *cp_tag.rbegin() = '\0';
00436             if( root_tagged || other_set_tagged )
00437             {
00438                 CCMIONode rootNode;
00439                 if( kCCMIONoErr == CCMIOGetEntityNode( &error, rootID, &rootNode ) )
00440                 {
00441                     CCMIOWriteOptstr( &error, processorID, "CreatingProgram", &cp_tag[0] );
00442                     CHK_SET_CCMERR( error, "Trouble setting creating program" );
00443                 }
00444             }
00445         }
00446     }
00447 
00448     CCMIONewEntity( &error, rootID, kCCMIOProblemDescription, NULL, &problemID );
00449     CHK_SET_CCMERR( error, "Trouble creating problem node" );
00450 
00451     // Write material types and other info
00452     for( unsigned int i = 0; i < matset_data.size(); i++ )
00453     {
00454         if( !matset_data[i].setName.empty() )
00455         {
00456             CCMIONewIndexedEntity( &error, problemID, kCCMIOCellType, matset_data[i].matsetId,
00457                                    matset_data[i].setName.c_str(), &id );
00458             CHK_SET_CCMERR( error, "Failure creating celltype node" );
00459 
00460             CCMIOWriteOptstr( &error, id, "MaterialType", matset_data[i].setName.c_str() );
00461             CHK_SET_CCMERR( error, "Error assigning material name" );
00462         }
00463         else
00464         {
00465             char dum_name[NAME_TAG_SIZE];
00466             std::ostringstream os;
00467             std::string mat_name = "Material", temp_str;
00468             os << mat_name << ( i + 1 );
00469             temp_str = os.str();
00470             strcpy( dum_name, temp_str.c_str() );
00471             CCMIONewIndexedEntity( &error, problemID, kCCMIOCellType, matset_data[i].matsetId, dum_name, &id );
00472             CHK_SET_CCMERR( error, "Failure creating celltype node" );
00473 
00474             CCMIOWriteOptstr( &error, id, "MaterialType", dum_name );
00475             CHK_SET_CCMERR( error, "Error assigning material name" );
00476 
00477             os.str( "" );
00478         }
00479         rval = write_int_option( "MaterialId", matset_data[i].setHandle, mMaterialIdTag, id );MB_CHK_SET_ERR( rval, "Trouble writing MaterialId option" );
00480 
00481         rval = write_int_option( "Radiation", matset_data[i].setHandle, mRadiationTag, id );MB_CHK_SET_ERR( rval, "Trouble writing Radiation option" );
00482 
00483         rval = write_int_option( "PorosityId", matset_data[i].setHandle, mPorosityIdTag, id );MB_CHK_SET_ERR( rval, "Trouble writing PorosityId option" );
00484 
00485         rval = write_int_option( "SpinId", matset_data[i].setHandle, mSpinIdTag, id );MB_CHK_SET_ERR( rval, "Trouble writing SpinId option" );
00486 
00487         rval = write_int_option( "GroupId", matset_data[i].setHandle, mGroupIdTag, id );MB_CHK_SET_ERR( rval, "Trouble writing GroupId option" );
00488 
00489         rval = write_int_option( "ColorIdx", matset_data[i].setHandle, mColorIdxTag, id );MB_CHK_SET_ERR( rval, "Trouble writing ColorIdx option" );
00490 
00491         rval = write_int_option( "ProcessorId", matset_data[i].setHandle, mProcessorIdTag, id );MB_CHK_SET_ERR( rval, "Trouble writing ProcessorId option" );
00492 
00493         rval = write_int_option( "LightMaterial", matset_data[i].setHandle, mLightMaterialTag, id );MB_CHK_SET_ERR( rval, "Trouble writing LightMaterial option." );
00494 
00495         rval = write_int_option( "FreeSurfaceMaterial", matset_data[i].setHandle, mFreeSurfaceMaterialTag, id );MB_CHK_SET_ERR( rval, "Trouble writing FreeSurfaceMaterial option" );
00496 
00497         rval = write_dbl_option( "Thickness", matset_data[i].setHandle, mThicknessTag, id );MB_CHK_SET_ERR( rval, "Trouble writing Thickness option" );
00498 
00499         rval = write_str_option( "MaterialType", matset_data[i].setHandle, mMaterialTypeTag, id );MB_CHK_SET_ERR( rval, "Trouble writing MaterialType option" );
00500     }
00501 
00502     // Write neumann set info
00503     for( unsigned int i = 0; i < neuset_data.size(); i++ )
00504     {
00505         // Use the label to encode the id
00506         std::ostringstream dum_id;
00507         dum_id << neuset_data[i].neusetId;
00508         CCMIONewIndexedEntity( &error, problemID, kCCMIOBoundaryRegion, neuset_data[i].neusetId, dum_id.str().c_str(),
00509                                &id );
00510         CHK_SET_CCMERR( error, "Failure creating BoundaryRegion node" );
00511 
00512         rval = write_str_option( "BoundaryName", neuset_data[i].setHandle, mNameTag, id );MB_CHK_SET_ERR( rval, "Trouble writing boundary type number" );
00513 
00514         rval = write_str_option( "BoundaryType", neuset_data[i].setHandle, mBoundaryTypeTag, id );MB_CHK_SET_ERR( rval, "Trouble writing boundary type number" );
00515 
00516         rval = write_int_option( "ProstarRegionNumber", neuset_data[i].setHandle, mProstarRegionNumberTag, id );MB_CHK_SET_ERR( rval, "Trouble writing prostar region number" );
00517     }
00518 
00519     CCMIOWriteState( &error, stateID, problemID, "Example state" );
00520     CHK_SET_CCMERR( error, "Failure writing problem state" );
00521 
00522     // Get cell types; reuse cell ids array
00523     //  for (i = 0, rit = all_elems.begin(); i < num_elems; i++, ++rit) {
00524     //    egids[i] = ccm_types[mbImpl->type_from_handle(*rit)];
00525     //    assert(-1 != egids[i]);
00526     //  }
00527 
00528     return MB_SUCCESS;
00529 }
00530 
00531 ErrorCode WriteCCMIO::write_int_option( const char* opt_name, EntityHandle seth, Tag& tag, CCMIOID& node )
00532 {
00533     ErrorCode rval;
00534 
00535     if( !tag )
00536     {
00537         rval = mbImpl->tag_get_handle( opt_name, 1, MB_TYPE_INTEGER, tag );
00538         // Return success since that just means we don't have to write this option
00539         if( MB_SUCCESS != rval ) return MB_SUCCESS;
00540     }
00541 
00542     int dum_val;
00543     rval = mbImpl->tag_get_data( tag, &seth, 1, &dum_val );
00544     // Return success since that just means we don't have to write this option
00545     if( MB_SUCCESS != rval ) return MB_SUCCESS;
00546 
00547     CCMIOError error = kCCMIONoErr;
00548     CCMIOWriteOpti( &error, node, opt_name, dum_val );
00549     CHK_SET_CCMERR( error, "Trouble writing int option" );
00550 
00551     return MB_SUCCESS;
00552 }
00553 
00554 ErrorCode WriteCCMIO::write_dbl_option( const char* opt_name, EntityHandle seth, Tag& tag, CCMIOID& node )
00555 {
00556     ErrorCode rval;
00557 
00558     if( !tag )
00559     {
00560         rval = mbImpl->tag_get_handle( opt_name, 1, MB_TYPE_DOUBLE, tag );
00561         // Return success since that just means we don't have to write this option
00562         if( MB_SUCCESS != rval ) return MB_SUCCESS;
00563     }
00564 
00565     double dum_val;
00566     rval = mbImpl->tag_get_data( tag, &seth, 1, &dum_val );
00567     // Return success since that just means we don't have to write this option
00568     if( MB_SUCCESS != rval ) return MB_SUCCESS;
00569 
00570     CCMIOError error = kCCMIONoErr;
00571     CCMIOWriteOptf( &error, node, opt_name, dum_val );
00572     CHK_SET_CCMERR( error, "Trouble writing int option" );
00573 
00574     return MB_SUCCESS;
00575 }
00576 
00577 ErrorCode WriteCCMIO::write_str_option( const char* opt_name,
00578                                         EntityHandle seth,
00579                                         Tag& tag,
00580                                         CCMIOID& node,
00581                                         const char* other_name )
00582 {
00583     int tag_size;
00584     ErrorCode rval;
00585 
00586     if( !tag )
00587     {
00588         rval = mbImpl->tag_get_handle( opt_name, 0, MB_TYPE_OPAQUE, tag, MB_TAG_ANY );
00589         // Return success since that just means we don't have to write this option
00590         if( MB_SUCCESS != rval ) return MB_SUCCESS;
00591     }
00592 
00593     rval = mbImpl->tag_get_bytes( tag, tag_size );
00594     if( MB_SUCCESS != rval ) return MB_SUCCESS;
00595     std::vector< char > opt_val( tag_size + 1 );
00596 
00597     rval = mbImpl->tag_get_data( tag, &seth, 1, &opt_val[0] );
00598     if( MB_SUCCESS != rval ) return MB_SUCCESS;
00599 
00600     // Null-terminate if necessary
00601     if( std::find( opt_val.begin(), opt_val.end(), '\0' ) == opt_val.end() ) *opt_val.rbegin() = '\0';
00602 
00603     CCMIOError error = kCCMIONoErr;
00604     if( other_name )
00605     {
00606         CCMIOWriteOptstr( &error, node, other_name, &opt_val[0] );
00607         CHK_SET_CCMERR( error, "Failure writing an option string MaterialType" );
00608     }
00609     else
00610     {
00611         CCMIOWriteOptstr( &error, node, opt_name, &opt_val[0] );
00612         CHK_SET_CCMERR( error, "Failure writing an option string MaterialType" );
00613     }
00614 
00615     return MB_SUCCESS;
00616 }
00617 
00618 ErrorCode WriteCCMIO::gather_matset_info( std::vector< EntityHandle >& matsets,
00619                                           std::vector< MaterialSetData >& matset_data,
00620                                           Range& all_verts )
00621 {
00622     ErrorCode result;
00623     matset_data.resize( matsets.size() );
00624     if( 1 == matsets.size() && 0 == matsets[0] )
00625     {
00626         // Whole mesh
00627         mWholeMesh = true;
00628 
00629         result = mbImpl->get_entities_by_dimension( 0, mDimension, matset_data[0].elems );MB_CHK_SET_ERR( result, "Trouble getting all elements in mesh" );
00630         result = mWriteIface->gather_nodes_from_elements( matset_data[0].elems, mEntityMark, all_verts );MB_CHK_SET_ERR( result, "Trouble gathering nodes from elements" );
00631 
00632         return result;
00633     }
00634 
00635     std::vector< unsigned char > marks;
00636     for( unsigned int i = 0; i < matsets.size(); i++ )
00637     {
00638         EntityHandle this_set = matset_data[i].setHandle = matsets[i];
00639 
00640         // Get all Entity Handles in the set
00641         result = mbImpl->get_entities_by_dimension( this_set, mDimension, matset_data[i].elems, true );MB_CHK_SET_ERR( result, "Trouble getting m-dimensional ents" );
00642 
00643         // Get all connected vertices
00644         result = mWriteIface->gather_nodes_from_elements( matset_data[i].elems, mEntityMark, all_verts );MB_CHK_SET_ERR( result, "Trouble getting vertices for a matset" );
00645 
00646         // Check for consistent entity type
00647         EntityType start_type = mbImpl->type_from_handle( *matset_data[i].elems.begin() );
00648         if( start_type == mbImpl->type_from_handle( *matset_data[i].elems.rbegin() ) )
00649             matset_data[i].entityType = start_type;
00650 
00651         // Mark elements in this matset
00652         marks.resize( matset_data[i].elems.size(), 0x1 );
00653         result = mbImpl->tag_set_data( mEntityMark, matset_data[i].elems, &marks[0] );MB_CHK_SET_ERR( result, "Couln't mark entities being output" );
00654 
00655         // Get id for this matset
00656         result = mbImpl->tag_get_data( mMaterialSetTag, &this_set, 1, &matset_data[i].matsetId );MB_CHK_SET_ERR( result, "Couln't get global id for material set" );
00657 
00658         // Get name for this matset
00659         if( mNameTag )
00660         {
00661             char dum_name[NAME_TAG_SIZE];
00662             result = mbImpl->tag_get_data( mNameTag, &this_set, 1, dum_name );
00663             if( MB_SUCCESS == result ) matset_data[i].setName = dum_name;
00664 
00665             // Reset success, so later checks don't fail
00666             result = MB_SUCCESS;
00667         }
00668     }
00669 
00670     if( all_verts.empty() )
00671     {
00672         MB_SET_ERR( MB_FILE_WRITE_ERROR, "No vertices from elements" );
00673     }
00674 
00675     return MB_SUCCESS;
00676 }
00677 
00678 ErrorCode WriteCCMIO::gather_neuset_info( std::vector< EntityHandle >& neusets,
00679                                           std::vector< NeumannSetData >& neuset_info )
00680 {
00681     ErrorCode result;
00682 
00683     neuset_info.resize( neusets.size() );
00684     for( unsigned int i = 0; i < neusets.size(); i++ )
00685     {
00686         EntityHandle this_set = neuset_info[i].setHandle = neusets[i];
00687 
00688         // Get all Entity Handles of one less dimension than that being output
00689         result = mbImpl->get_entities_by_dimension( this_set, mDimension - 1, neuset_info[i].elems, true );MB_CHK_SET_ERR( result, "Trouble getting (m-1)-dimensional ents for neuset" );
00690 
00691         result = mbImpl->tag_get_data( mGlobalIdTag, &this_set, 1, &neuset_info[i].neusetId );
00692         if( MB_TAG_NOT_FOUND == result )
00693         {
00694             result = mbImpl->tag_get_data( mNeumannSetTag, &this_set, 1, &neuset_info[i].neusetId );
00695             if( MB_SUCCESS != result )
00696                 // Need some id; use the loop iteration number
00697                 neuset_info[i].neusetId = i;
00698         }
00699 
00700         // Get name for this neuset
00701         if( mNameTag )
00702         {
00703             char dum_name[NAME_TAG_SIZE];
00704             result = mbImpl->tag_get_data( mNameTag, &this_set, 1, dum_name );
00705             if( MB_SUCCESS == result ) neuset_info[i].setName = dum_name;
00706 
00707             // Reset success, so later checks don't fail
00708             result = MB_SUCCESS;
00709         }
00710     }
00711 
00712     return MB_SUCCESS;
00713 }
00714 
00715 ErrorCode WriteCCMIO::get_gids( const Range& ents, int*& gids, int& minid, int& maxid )
00716 {
00717     int num_ents     = ents.size();
00718     gids             = new int[num_ents];
00719     ErrorCode result = mbImpl->tag_get_data( mGlobalIdTag, ents, &gids[0] );MB_CHK_SET_ERR( result, "Couldn't get global id data" );
00720     minid = *std::min_element( gids, gids + num_ents );
00721     maxid = *std::max_element( gids, gids + num_ents );
00722     if( 0 == minid )
00723     {
00724         // gids need to be assigned
00725         for( int i = 1; i <= num_ents; i++ )
00726             gids[i] = i;
00727         result = mbImpl->tag_set_data( mGlobalIdTag, ents, &gids[0] );MB_CHK_SET_ERR( result, "Couldn't set global id data" );
00728         maxid = num_ents;
00729     }
00730 
00731     return MB_SUCCESS;
00732 }
00733 
00734 ErrorCode WriteCCMIO::write_nodes( CCMIOID rootID, const Range& verts, const int dimension, CCMIOID& verticesID )
00735 {
00736     // Get/write map (global ids) first (gids already assigned)
00737     unsigned int num_verts = verts.size();
00738     std::vector< int > vgids( num_verts );
00739     ErrorCode result = mbImpl->tag_get_data( mGlobalIdTag, verts, &vgids[0] );MB_CHK_SET_ERR( result, "Failed to get global ids for vertices" );
00740 
00741     // Create the map node for vertex ids, and write them to that node
00742     CCMIOID mapID;
00743     CCMIOError error = kCCMIONoErr;
00744     CCMIONewEntity( &error, rootID, kCCMIOMap, "Vertex map", &mapID );
00745     CHK_SET_CCMERR( error, "Failure creating Vertex map node" );
00746 
00747     int maxid = *std::max_element( vgids.begin(), vgids.end() );
00748 
00749     CCMIOWriteMap( &error, mapID, CCMIOSIZEC( num_verts ), CCMIOSIZEC( maxid ), &vgids[0], CCMIOINDEXC( kCCMIOStart ),
00750                    CCMIOINDEXC( kCCMIOEnd ) );
00751     CHK_SET_CCMERR( error, "Problem writing node map" );
00752 
00753     // Create the vertex coordinate node, and write it
00754     CCMIONewEntity( &error, rootID, kCCMIOVertices, "Vertices", &verticesID );
00755     CHK_SET_CCMERR( error, "Trouble creating vertices node" );
00756 
00757     // Get the vertex locations
00758     double* coords = new double[3 * num_verts];
00759     std::vector< double* > coord_arrays( 3 );
00760     // Cppcheck warning (false positive): variable coord_arrays is assigned a value that is never
00761     // used
00762     coord_arrays[0] = coords;
00763     coord_arrays[1] = coords + num_verts;
00764     coord_arrays[2] = ( dimension == 3 ? coords + 2 * num_verts : NULL );
00765     result          = mWriteIface->get_node_coords( -1, verts.begin(), verts.end(), 3 * num_verts, coords );
00766     if( result != MB_SUCCESS )
00767     {
00768         delete[] coords;
00769         return result;
00770     }
00771 
00772     // Transform coordinates, if necessary
00773     result = transform_coords( dimension, num_verts, coords );
00774     if( result != MB_SUCCESS )
00775     {
00776         delete[] coords;
00777         MB_SET_ERR( result, "Trouble transforming vertex coordinates" );
00778     }
00779 
00780     // Write the vertices
00781     CCMIOWriteVerticesd( &error, verticesID, CCMIOSIZEC( dimension ), 1.0, mapID, coords, CCMIOINDEXC( kCCMIOStart ),
00782                          CCMIOINDEXC( kCCMIOEnd ) );
00783     CHK_SET_CCMERR( error, "CCMIOWriteVertices failed" );
00784 
00785     // Clean up
00786     delete[] coords;
00787 
00788     return MB_SUCCESS;
00789 }
00790 
00791 ErrorCode WriteCCMIO::transform_coords( const int dimension, const int num_nodes, double* coords )
00792 {
00793     Tag trans_tag;
00794     ErrorCode result = mbImpl->tag_get_handle( MESH_TRANSFORM_TAG_NAME, 16, MB_TYPE_DOUBLE, trans_tag );
00795     if( result == MB_TAG_NOT_FOUND )
00796         return MB_SUCCESS;
00797     else if( MB_SUCCESS != result )
00798         return result;
00799     double trans_matrix[16];
00800     const EntityHandle mesh = 0;
00801     result                  = mbImpl->tag_get_data( trans_tag, &mesh, 1, trans_matrix );MB_CHK_SET_ERR( result, "Couldn't get transform data" );
00802 
00803     double* tmp_coords = coords;
00804     for( int i = 0; i < num_nodes; i++, tmp_coords += 1 )
00805     {
00806         double vec1[3] = { 0.0, 0.0, 0.0 };
00807         for( int row = 0; row < 3; row++ )
00808         {
00809             vec1[row] += ( trans_matrix[( row * 4 ) + 0] * coords[0] );
00810             vec1[row] += ( trans_matrix[( row * 4 ) + 1] * coords[num_nodes] );
00811             if( 3 == dimension ) vec1[row] += ( trans_matrix[( row * 4 ) + 2] * coords[2 * num_nodes] );
00812         }
00813 
00814         coords[0]             = vec1[0];
00815         coords[num_nodes]     = vec1[1];
00816         coords[2 * num_nodes] = vec1[2];
00817     }
00818 
00819     return MB_SUCCESS;
00820 }
00821 
00822 ErrorCode WriteCCMIO::write_cells_and_faces( CCMIOID rootID,
00823                                              std::vector< MaterialSetData >& matset_data,
00824                                              std::vector< NeumannSetData >& neuset_data,
00825                                              Range& /* verts */,
00826                                              CCMIOID& topologyID )
00827 {
00828     std::vector< int > connect;
00829     ErrorCode result;
00830     CCMIOID cellMapID, cells;
00831     CCMIOError error = kCCMIONoErr;
00832 
00833     // Don't usually have anywhere near 31 nodes per element
00834     connect.reserve( 31 );
00835     Range::const_iterator rit;
00836 
00837     // Create the topology node, and the cell and cell map nodes
00838     CCMIONewEntity( &error, rootID, kCCMIOTopology, "Topology", &topologyID );
00839     CHK_SET_CCMERR( error, "Trouble creating topology node" );
00840 
00841     CCMIONewEntity( &error, rootID, kCCMIOMap, "Cell map", &cellMapID );
00842     CHK_SET_CCMERR( error, "Failure creating Cell Map node" );
00843 
00844     CCMIONewEntity( &error, topologyID, kCCMIOCells, "Cells", &cells );
00845     CHK_SET_CCMERR( error, "Trouble creating Cell node under Topology node" );
00846 
00847     //================================================
00848     // Loop over material sets, doing each one at a time
00849     //================================================
00850     Range all_elems;
00851     unsigned int i, num_elems = 0;
00852     int max_id = 1;
00853     std::vector< int > egids;
00854     int tot_elems = 0;
00855 
00856     for( unsigned int m = 0; m < matset_data.size(); m++ )
00857         tot_elems += matset_data[m].elems.size();
00858 
00859     for( unsigned int m = 0; m < matset_data.size(); m++ )
00860     {
00861         unsigned int this_num = matset_data[m].elems.size();
00862 
00863         //================================================
00864         // Save all elements being output
00865         //================================================
00866         all_elems.merge( matset_data[m].elems );
00867 
00868         //================================================
00869         // Assign global ids for elements being written
00870         //================================================
00871         egids.resize( matset_data[m].elems.size() );
00872         for( i = 0; i < this_num; i++ )
00873             egids[i] = max_id++;
00874         result = mbImpl->tag_set_data( mGlobalIdTag, matset_data[m].elems, &egids[0] );MB_CHK_SET_ERR( result, "Failed to assign global ids for all elements being written" );
00875 
00876         //================================================
00877         // Write cell ids and material types for this matset; reuse egids for cell mat type
00878         //================================================
00879         CCMIOWriteMap( &error, cellMapID, CCMIOSIZEC( tot_elems ), CCMIOSIZEC( tot_elems ), &egids[0],
00880                        CCMIOINDEXC( 0 == m ? kCCMIOStart : num_elems ),
00881                        CCMIOINDEXC( matset_data.size() == m ? kCCMIOEnd : num_elems + this_num ) );
00882         CHK_SET_CCMERR( error, "Trouble writing cell map" );
00883 
00884         if( -1 == matset_data[m].matsetId )
00885         {
00886             for( i = 0; i < this_num; i++ )
00887                 egids[i] = m;
00888         }
00889         else
00890         {
00891             for( i = 0; i < this_num; i++ )
00892                 egids[i] = matset_data[m].matsetId;
00893         }
00894 
00895         CCMIOWriteCells( &error, cells, cellMapID, &egids[0], CCMIOINDEXC( 0 == m ? kCCMIOStart : num_elems ),
00896                          CCMIOINDEXC( matset_data.size() == m ? kCCMIOEnd : num_elems + this_num ) );
00897         CHK_SET_CCMERR( error, "Trouble writing Cell node" );
00898 
00899         //================================================
00900         // Write cell entity types
00901         //================================================
00902         const EntityHandle* conn;
00903         int num_conn;
00904         int has_mid_nodes[4];
00905         std::vector< EntityHandle > storage;
00906         for( i = 0, rit = matset_data[m].elems.begin(); i < this_num; i++, ++rit )
00907         {
00908             result = mbImpl->get_connectivity( *rit, conn, num_conn, false, &storage );MB_CHK_SET_ERR( result, "Trouble getting connectivity for entity type check" );
00909             CN::HasMidNodes( mbImpl->type_from_handle( *rit ), num_conn, has_mid_nodes );
00910             egids[i] = moab_to_ccmio_type( mbImpl->type_from_handle( *rit ), has_mid_nodes );
00911         }
00912 
00913         CCMIOWriteOpt1i( &error, cells, "CellTopologyType", CCMIOSIZEC( tot_elems ), &egids[0],
00914                          CCMIOINDEXC( 0 == m ? kCCMIOStart : num_elems ),
00915                          CCMIOINDEXC( matset_data.size() == m ? kCCMIOEnd : num_elems + this_num ) );
00916         CHK_SET_CCMERR( error, "Failed to write cell topo types" );
00917 
00918         num_elems += this_num;
00919     }
00920 
00921     //================================================
00922     // Get skin and neumann set faces
00923     //================================================
00924     Range neuset_facets, skin_facets;
00925     Skinner skinner( mbImpl );
00926     result = skinner.find_skin( 0, all_elems, mDimension - 1, skin_facets );MB_CHK_SET_ERR( result, "Failed to get skin facets" );
00927 
00928     // Remove neumann set facets from skin facets, we have to output these
00929     // separately
00930     for( i = 0; i < neuset_data.size(); i++ )
00931         neuset_facets.merge( neuset_data[i].elems );
00932 
00933     skin_facets -= neuset_facets;
00934     // Make neuset_facets the union, and get ids for them
00935     neuset_facets.merge( skin_facets );
00936     result = mWriteIface->assign_ids( neuset_facets, mGlobalIdTag, 1 );
00937 
00938     int fmaxid = neuset_facets.size();
00939 
00940     //================================================
00941     // Write external faces
00942     //================================================
00943     for( i = 0; i < neuset_data.size(); i++ )
00944     {
00945         Range::reverse_iterator rrit;
00946         unsigned char cmarks[2];
00947         Range ext_faces;
00948         std::vector< EntityHandle > mcells;
00949         // Removing the faces connected to two regions
00950         for( rrit = neuset_data[i].elems.rbegin(); rrit != neuset_data[i].elems.rend(); ++rrit )
00951         {
00952             mcells.clear();
00953             result = mbImpl->get_adjacencies( &( *rrit ), 1, mDimension, false, mcells );MB_CHK_SET_ERR( result, "Trouble getting bounding cells" );
00954 
00955             result = mbImpl->tag_get_data( mEntityMark, &mcells[0], mcells.size(), cmarks );MB_CHK_SET_ERR( result, "Trouble getting mark tags on cells bounding facets" );
00956 
00957             if( mcells.size() == 2 && ( mWholeMesh || ( cmarks[0] && cmarks[1] ) ) )
00958             {
00959             }
00960             else
00961             {
00962                 // External face
00963                 ext_faces.insert( *rrit );
00964             }
00965         }
00966         if( ext_faces.size() != 0 && neuset_data[i].neusetId != 0 )
00967         {
00968             result = write_external_faces( rootID, topologyID, neuset_data[i].neusetId, ext_faces );MB_CHK_SET_ERR( result, "Trouble writing Neumann set facets" );
00969         }
00970         ext_faces.clear();
00971     }
00972 
00973     if( !skin_facets.empty() )
00974     {
00975         result = write_external_faces( rootID, topologyID, 0, skin_facets );MB_CHK_SET_ERR( result, "Trouble writing skin facets" );
00976     }
00977 
00978     //================================================
00979     // Now internal faces; loop over elements, do each face on the element
00980     //================================================
00981     // Mark tag, for face marking on each non-polyhedral element
00982 
00983     if( num_elems > 1 )
00984     {  // No internal faces for just one element
00985         Tag fmark_tag;
00986         unsigned char mval = 0x0, omval;
00987         result = mbImpl->tag_get_handle( "__fmark", 1, MB_TYPE_OPAQUE, fmark_tag, MB_TAG_DENSE | MB_TAG_CREAT, &mval );MB_CHK_SET_ERR( result, "Couldn't create mark tag" );
00988 
00989         std::vector< EntityHandle > tmp_face_cells, storage;
00990         std::vector< int > iface_connect, iface_cells;
00991         EntityHandle tmp_connect[CN::MAX_NODES_PER_ELEMENT];  // tmp connect vector
00992         const EntityHandle *connectc, *oconnectc;
00993         int num_connectc;  // Cell connectivity
00994         const EntityHandle* connectf;
00995         int num_connectf;  // Face connectivity
00996 
00997         for( i = 0, rit = all_elems.begin(); i < num_elems; i++, ++rit )
00998         {
00999             EntityType etype = TYPE_FROM_HANDLE( *rit );
01000 
01001             //-----------------------
01002             // If not polyh, get mark
01003             //-----------------------
01004             if( MBPOLYHEDRON != etype && MBPOLYGON != etype )
01005             {
01006                 result = mbImpl->tag_get_data( fmark_tag, &( *rit ), 1, &mval );MB_CHK_SET_ERR( result, "Couldn't get mark data" );
01007             }
01008 
01009             //-----------------------
01010             // Get cell connectivity, and whether it's a polyhedron
01011             //-----------------------
01012             result = mbImpl->get_connectivity( *rit, connectc, num_connectc, false, &storage );MB_CHK_SET_ERR( result, "Couldn't get entity connectivity" );
01013 
01014             // If polyh, write faces directly
01015             bool is_polyh = ( MBPOLYHEDRON == etype );
01016 
01017             int num_facets = ( is_polyh ? num_connectc : CN::NumSubEntities( etype, mDimension - 1 ) );
01018 
01019             //----------------------------------------------------------
01020             // Loop over each facet of element, outputting it if not marked
01021             //----------------------------------------------------------
01022             for( int f = 0; f < num_facets; f++ )
01023             {
01024                 //.............................................
01025                 // If this face marked, skip
01026                 //.............................................
01027                 if( !is_polyh && ( ( mval >> f ) & 0x1 ) ) continue;
01028 
01029                 //.................
01030                 // Get face connect and adj cells
01031                 //.................
01032                 if( !is_polyh )
01033                 {
01034                     // (from CN)
01035                     CN::SubEntityConn( connectc, etype, mDimension - 1, f, tmp_connect, num_connectf );
01036                     connectf = tmp_connect;
01037                 }
01038                 else
01039                 {
01040                     // Directly
01041                     result = mbImpl->get_connectivity( connectc[f], connectf, num_connectf, false );MB_CHK_SET_ERR( result, "Couldn't get polyhedron connectivity" );
01042                 }
01043 
01044                 //............................
01045                 // Get adj cells from face connect (same for poly's and not, since both usually
01046                 // go through vertices anyway)
01047                 //............................
01048                 tmp_face_cells.clear();
01049                 result = mbImpl->get_adjacencies( connectf, num_connectf, mDimension, false, tmp_face_cells );MB_CHK_SET_ERR( result, "Error getting adj hexes" );
01050 
01051                 //...............................
01052                 // If this face only bounds one cell, skip, since we exported external faces
01053                 // before this loop
01054                 //...............................
01055                 if( tmp_face_cells.size() != 2 ) continue;
01056 
01057                 //.................
01058                 // Switch cells so that *rit is always 1st (face connectivity is always written such
01059                 // that that one is with forward sense)
01060                 //.................
01061                 int side_num = 0, sense = 0, offset = 0;
01062                 if( !is_polyh && tmp_face_cells[0] != *rit )
01063                 {
01064                     EntityHandle tmph = tmp_face_cells[0];
01065                     tmp_face_cells[0] = tmp_face_cells[1];
01066                     tmp_face_cells[1] = tmph;
01067                 }
01068 
01069                 //.................
01070                 // Save ids of cells
01071                 //.................
01072                 assert( tmp_face_cells[0] != tmp_face_cells[1] );
01073                 iface_cells.resize( iface_cells.size() + 2 );
01074                 result = mbImpl->tag_get_data( mGlobalIdTag, &tmp_face_cells[0], tmp_face_cells.size(),
01075                                                &iface_cells[iface_cells.size() - 2] );MB_CHK_SET_ERR( result, "Trouble getting global ids for bounded cells" );
01076                 iface_connect.push_back( num_connectf );
01077 
01078                 //.................
01079                 // Save indices of face vertices
01080                 //.................
01081                 unsigned int tmp_size = iface_connect.size();
01082                 iface_connect.resize( tmp_size + num_connectf );
01083                 result = mbImpl->tag_get_data( mGlobalIdTag, connectf, num_connectf, &iface_connect[tmp_size] );MB_CHK_SET_ERR( result, "Trouble getting global id for internal face" );
01084 
01085                 //.................
01086                 // Mark other cell with the right side #
01087                 //.................
01088                 if( !is_polyh )
01089                 {
01090                     // Mark other cell for this face, if there is another cell
01091 
01092                     result = mbImpl->get_connectivity( tmp_face_cells[1], oconnectc, num_connectc, false, &storage );MB_CHK_SET_ERR( result, "Couldn't get other entity connectivity" );
01093 
01094                     // Get side number in other cell
01095                     CN::SideNumber( TYPE_FROM_HANDLE( tmp_face_cells[1] ), oconnectc, connectf, num_connectf,
01096                                     mDimension - 1, side_num, sense, offset );
01097                     // Set mark for that face on the other cell
01098                     result = mbImpl->tag_get_data( fmark_tag, &tmp_face_cells[1], 1, &omval );MB_CHK_SET_ERR( result, "Couldn't get mark data for other cell" );
01099                 }
01100 
01101                 omval |= ( 0x1 << (unsigned int)side_num );
01102                 result = mbImpl->tag_set_data( fmark_tag, &tmp_face_cells[1], 1, &omval );MB_CHK_SET_ERR( result, "Couldn't set mark data for other cell" );
01103             }  // Loop over faces in elem
01104         }      // Loop over elems
01105 
01106         //================================================
01107         // Write internal faces
01108         //================================================
01109 
01110         CCMIOID mapID;
01111         CCMIONewEntity( &error, rootID, kCCMIOMap, NULL, &mapID );
01112         CHK_SET_CCMERR( error, "Trouble creating Internal Face map node" );
01113 
01114         unsigned int num_ifaces = iface_cells.size() / 2;
01115 
01116         // Set gids for internal faces; reuse egids
01117         egids.resize( num_ifaces );
01118         for( i = 1; i <= num_ifaces; i++ )
01119             egids[i - 1] = fmaxid + i;
01120         CCMIOWriteMap( &error, mapID, CCMIOSIZEC( num_ifaces ), CCMIOSIZEC( fmaxid + num_ifaces ), &egids[0],
01121                        CCMIOINDEXC( kCCMIOStart ), CCMIOINDEXC( kCCMIOEnd ) );
01122         CHK_SET_CCMERR( error, "Trouble writing Internal Face map node" );
01123 
01124         CCMIOID id;
01125         CCMIONewEntity( &error, topologyID, kCCMIOInternalFaces, "Internal faces", &id );
01126         CHK_SET_CCMERR( error, "Failed to create Internal face node under Topology node" );
01127         CCMIOWriteFaces( &error, id, kCCMIOInternalFaces, mapID, CCMIOSIZEC( iface_connect.size() ), &iface_connect[0],
01128                          CCMIOINDEXC( kCCMIOStart ), CCMIOINDEXC( kCCMIOEnd ) );
01129         CHK_SET_CCMERR( error, "Failure writing Internal face connectivity" );
01130         CCMIOWriteFaceCells( &error, id, kCCMIOInternalFaces, mapID, &iface_cells[0], CCMIOINDEXC( kCCMIOStart ),
01131                              CCMIOINDEXC( kCCMIOEnd ) );
01132         CHK_SET_CCMERR( error, "Failure writing Internal face cells" );
01133     }
01134 
01135     return MB_SUCCESS;
01136 }
01137 
01138 int WriteCCMIO::moab_to_ccmio_type( EntityType etype, int has_mid_nodes[] )
01139 {
01140     int ctype = -1;
01141     if( has_mid_nodes[0] || has_mid_nodes[2] || has_mid_nodes[3] ) return ctype;
01142 
01143     switch( etype )
01144     {
01145         case MBVERTEX:
01146             ctype = 1;
01147             break;
01148         case MBEDGE:
01149             if( !has_mid_nodes[1] )
01150                 ctype = 2;
01151             else
01152                 ctype = 28;
01153             break;
01154         case MBQUAD:
01155             if( has_mid_nodes[1] )
01156                 ctype = 4;
01157             else
01158                 ctype = 3;
01159             break;
01160         case MBTET:
01161             if( has_mid_nodes[1] )
01162                 ctype = 23;
01163             else
01164                 ctype = 13;
01165             break;
01166         case MBPRISM:
01167             if( has_mid_nodes[1] )
01168                 ctype = 22;
01169             else
01170                 ctype = 12;
01171             break;
01172         case MBPYRAMID:
01173             if( has_mid_nodes[1] )
01174                 ctype = 24;
01175             else
01176                 ctype = 14;
01177             break;
01178         case MBHEX:
01179             if( has_mid_nodes[1] )
01180                 ctype = 21;
01181             else
01182                 ctype = 11;
01183             break;
01184         case MBPOLYHEDRON:
01185             ctype = 255;
01186             break;
01187         default:
01188             break;
01189     }
01190 
01191     return ctype;
01192 }
01193 
01194 ErrorCode WriteCCMIO::write_external_faces( CCMIOID rootID, CCMIOID topologyID, int set_num, Range& facets )
01195 {
01196     CCMIOError error = kCCMIONoErr;
01197     CCMIOID mapID, id;
01198 
01199     // Get gids for these faces
01200     int *gids        = NULL, minid, maxid;
01201     ErrorCode result = get_gids( facets, gids, minid, maxid );MB_CHK_SET_ERR( result, "Trouble getting global ids for facets" );
01202 
01203     // Write the face id map
01204     CCMIONewEntity( &error, rootID, kCCMIOMap, NULL, &mapID );
01205     CHK_SET_CCMERR( error, "Problem creating face id map" );
01206 
01207     CCMIOWriteMap( &error, mapID, CCMIOSIZEC( facets.size() ), CCMIOSIZEC( maxid ), gids, CCMIOINDEXC( kCCMIOStart ),
01208                    CCMIOINDEXC( kCCMIOEnd ) );
01209     CHK_SET_CCMERR( error, "Problem writing face id map" );
01210 
01211     // Get the connectivity of the faces; set size by how many verts in last facet
01212     const EntityHandle* connect;
01213     int num_connect;
01214     result = mbImpl->get_connectivity( *facets.rbegin(), connect, num_connect );MB_CHK_SET_ERR( result, "Failed to get connectivity of last facet" );
01215     std::vector< int > fconnect( facets.size() * ( num_connect + 1 ) );
01216 
01217     result = mWriteIface->get_element_connect( facets.begin(), facets.end(), num_connect, mGlobalIdTag, fconnect.size(),
01218                                                &fconnect[0], true );MB_CHK_SET_ERR( result, "Failed to get facet connectivity" );
01219 
01220     // Get and write a new external face entity
01221     CCMIONewIndexedEntity( &error, topologyID, kCCMIOBoundaryFaces, set_num, "Boundary faces", &id );
01222     CHK_SET_CCMERR( error, "Problem creating boundary face entity" );
01223 
01224     CCMIOWriteFaces( &error, id, kCCMIOBoundaryFaces, mapID, CCMIOSIZEC( fconnect.size() ), &fconnect[0],
01225                      CCMIOINDEXC( kCCMIOStart ), CCMIOINDEXC( kCCMIOEnd ) );
01226     CHK_SET_CCMERR( error, "Problem writing boundary faces" );
01227 
01228     // Get info on bounding cells; reuse fconnect
01229     std::vector< EntityHandle > cells;
01230     unsigned char cmarks[2];
01231     int i, j = 0;
01232     Range dead_facets;
01233     Range::iterator rit;
01234 
01235     // About error checking in this loop: if any facets have no bounding cells,
01236     // this is an error, since global ids for facets are computed outside this loop
01237     for( rit = facets.begin(), i = 0; rit != facets.end(); ++rit, i++ )
01238     {
01239         cells.clear();
01240 
01241         // Get cell then gid of cell
01242         result = mbImpl->get_adjacencies( &( *rit ), 1, mDimension, false, cells );MB_CHK_SET_ERR( result, "Trouble getting bounding cells" );
01243         if( cells.empty() )
01244         {
01245             MB_SET_ERR( MB_FILE_WRITE_ERROR, "External facet with no output bounding cell" );
01246         }
01247 
01248         // Check we don't bound more than one cell being output
01249         result = mbImpl->tag_get_data( mEntityMark, &cells[0], cells.size(), cmarks );MB_CHK_SET_ERR( result, "Trouble getting mark tags on cells bounding facets" );
01250         if( cells.size() == 2 && ( mWholeMesh || ( cmarks[0] && cmarks[1] ) ) )
01251         {
01252             MB_SET_ERR( MB_FILE_WRITE_ERROR, "External facet with two output bounding cells" );
01253         }
01254         else if( 1 == cells.size() && !mWholeMesh && !cmarks[0] )
01255         {
01256             MB_SET_ERR( MB_FILE_WRITE_ERROR, "External facet with no output bounding cells" );
01257         }
01258 
01259         // Make sure 1st cell is the one being output
01260         if( 2 == cells.size() && !( cmarks[0] | 0x0 ) && ( cmarks[1] & 0x1 ) ) cells[0] = cells[1];
01261 
01262         // Get gid for bounded cell
01263         result = mbImpl->tag_get_data( mGlobalIdTag, &cells[0], 1, &fconnect[j] );MB_CHK_SET_ERR( result, "Couldn't get global id tag for bounded cell" );
01264 
01265         j++;
01266     }
01267 
01268     // Write the bounding cell data
01269     CCMIOWriteFaceCells( &error, id, kCCMIOBoundaryFaces, mapID, &fconnect[0], CCMIOINDEXC( kCCMIOStart ),
01270                          CCMIOINDEXC( kCCMIOEnd ) );
01271     CHK_SET_CCMERR( error, "Problem writing boundary cell data" );
01272 
01273     return MB_SUCCESS;
01274 }
01275 
01276 ErrorCode WriteCCMIO::get_neuset_elems( EntityHandle neuset,
01277                                         int current_sense,
01278                                         Range& forward_elems,
01279                                         Range& reverse_elems )
01280 {
01281     Range neuset_elems, neuset_meshsets;
01282 
01283     // Get the sense tag; don't need to check return, might be an error if the tag
01284     // hasn't been created yet
01285     Tag sense_tag = 0;
01286     mbImpl->tag_get_handle( "SENSE", 1, MB_TYPE_INTEGER, sense_tag );
01287 
01288     // Get the entities in this set, non-recursive
01289     ErrorCode result = mbImpl->get_entities_by_handle( neuset, neuset_elems );
01290     if( MB_FAILURE == result ) return result;
01291 
01292     // Now remove the meshsets into the neuset_meshsets; first find the first meshset,
01293     Range::iterator range_iter = neuset_elems.begin();
01294     while( TYPE_FROM_HANDLE( *range_iter ) != MBENTITYSET && range_iter != neuset_elems.end() )
01295         ++range_iter;
01296 
01297     // Then, if there are some, copy them into neuset_meshsets and erase from neuset_elems
01298     if( range_iter != neuset_elems.end() )
01299     {
01300         std::copy( range_iter, neuset_elems.end(), range_inserter( neuset_meshsets ) );
01301         neuset_elems.erase( range_iter, neuset_elems.end() );
01302     }
01303 
01304     // OK, for the elements, check the sense of this set and copy into the right range
01305     // (if the sense is 0, copy into both ranges)
01306 
01307     // Need to step forward on list until we reach the right dimension
01308     Range::iterator dum_it = neuset_elems.end();
01309     --dum_it;
01310     int target_dim = CN::Dimension( TYPE_FROM_HANDLE( *dum_it ) );
01311     dum_it         = neuset_elems.begin();
01312     while( target_dim != CN::Dimension( TYPE_FROM_HANDLE( *dum_it ) ) && dum_it != neuset_elems.end() )
01313         ++dum_it;
01314 
01315     if( current_sense == 1 || current_sense == 0 )
01316         std::copy( dum_it, neuset_elems.end(), range_inserter( forward_elems ) );
01317     if( current_sense == -1 || current_sense == 0 )
01318         std::copy( dum_it, neuset_elems.end(), range_inserter( reverse_elems ) );
01319 
01320     // Now loop over the contained meshsets, getting the sense of those and calling this
01321     // function recursively
01322     for( range_iter = neuset_meshsets.begin(); range_iter != neuset_meshsets.end(); ++range_iter )
01323     {
01324         // First get the sense; if it's not there, by convention it's forward
01325         int this_sense;
01326         if( 0 == sense_tag || MB_FAILURE == mbImpl->tag_get_data( sense_tag, &( *range_iter ), 1, &this_sense ) )
01327             this_sense = 1;
01328 
01329         // Now get all the entities on this meshset, with the proper (possibly reversed) sense
01330         get_neuset_elems( *range_iter, this_sense * current_sense, forward_elems, reverse_elems );
01331     }
01332 
01333     return result;
01334 }
01335 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines