Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
ReadCCMIO.cpp
Go to the documentation of this file.
00001 #include <cstdlib>  // For exit()
00002 #include <vector>
00003 #include <map>
00004 #include <iostream>
00005 #include <string>
00006 #include <algorithm>
00007 
00008 #include "moab/CN.hpp"
00009 #include "moab/Range.hpp"
00010 #include "moab/Interface.hpp"
00011 #include "MBTagConventions.hpp"
00012 #include "Internals.hpp"
00013 #include "moab/ReadUtilIface.hpp"
00014 #include "moab/FileOptions.hpp"
00015 #include "ReadCCMIO.hpp"
00016 #include "moab/MeshTopoUtil.hpp"
00017 
00018 #include "ccmio.h"
00019 
00020 /*
00021  * CCMIO file structure
00022  *
00023  * Root
00024  *   State(kCCMIOState)
00025  *     Processor*
00026  *       Vertices
00027  *         ->ReadVerticesx, ReadMap
00028  *       Topology
00029  *         Boundary faces*(kCCMIOBoundaryFaces)
00030  *            ->ReadFaces, ReadFaceCells, ReadMap
00031  *         Internal faces(kCCMIOInternalFaces)
00032  *         Cells (kCCMIOCells)
00033  *            ->ReadCells (mapID), ReadMap, ReadCells (cellTypes)
00034  *       Solution
00035  *         Phase
00036  *           Field
00037  *             FieldData
00038  *   Problem(kCCMIOProblemDescription)
00039  *     CellType* (kCCMIOCellType)
00040  *       Index (GetEntityIndex), MaterialId(ReadOpti), MaterialType(ReadOptstr),
00041  *         PorosityId(ReadOpti), SpinId(ReadOpti), GroupId(ReadOpti)
00042  *
00043  * MaterialType (CCMIOReadOptstr in readexample)
00044  * constants (see readexample)
00045  * lagrangian data (CCMIOReadLagrangianData)
00046  * vertices label (CCMIOEntityDescription)
00047  * restart info: char solver[], iteratoins, time, char timeUnits[], angle
00048  *      (CCMIOReadRestartInfo, kCCMIORestartData), reference data?
00049  * phase:
00050  *   field: char name[], dims, CCMIODataType datatype, char units[]
00051  *       dims = kCCMIOScalar (CCMIOReadFieldDataf),
00052  *              kCCMIOVector (CCMIOReadMultiDimensionalFieldData),
00053  *              kCCMIOTensor
00054  * MonitoringSets: num, name (CellSet, VertexSet, BoundarySet, BlockSet, SplineSet, CoupleSet)
00055  *      CCMIOGetProstarSet, CCMIOReadOpt1i,
00056  */
00057 
00058 enum DataType
00059 {
00060     kScalar,
00061     kVector,
00062     kVertex,
00063     kCell,
00064     kInternalFace,
00065     kBoundaryFace,
00066     kBoundaryData,
00067     kBoundaryFaceData,
00068     kCellType
00069 };
00070 
00071 namespace moab
00072 {
00073 
00074 static int const kNValues         = 10;  // Number of values of each element to print
00075 static char const kDefaultState[] = "default";
00076 static char const kUnitsName[]    = "Units";
00077 static int const kVertOffset      = 2;
00078 static int const kCellInc         = 4;
00079 
00080 #define CHK_SET_CCMERR( ccm_err_code, ccm_err_msg )                                   \
00081     {                                                                                 \
00082         if( kCCMIONoErr != ( ccm_err_code ) && kCCMIONoFileErr != ( ccm_err_code ) && \
00083             kCCMIONoNodeErr != ( ccm_err_code ) )                                     \
00084             MB_SET_ERR( MB_FAILURE, ccm_err_msg );                                    \
00085     }
00086 
00087 ReaderIface* ReadCCMIO::factory( Interface* iface )
00088 {
00089     return new ReadCCMIO( iface );
00090 }
00091 
00092 ReadCCMIO::ReadCCMIO( Interface* impl )
00093     : mMaterialIdTag( 0 ), mMaterialTypeTag( 0 ), mRadiationTag( 0 ), mPorosityIdTag( 0 ), mSpinIdTag( 0 ),
00094       mGroupIdTag( 0 ), mColorIdxTag( 0 ), mProcessorIdTag( 0 ), mLightMaterialTag( 0 ), mFreeSurfaceMaterialTag( 0 ),
00095       mThicknessTag( 0 ), mProstarRegionNumberTag( 0 ), mBoundaryTypeTag( 0 ), mCreatingProgramTag( 0 ), mbImpl( impl ),
00096       hasSolution( false )
00097 {
00098     assert( impl != NULL );
00099 
00100     impl->query_interface( readMeshIface );
00101 
00102     // Initialize in case tag_get_handle fails below
00103     mMaterialSetTag  = 0;
00104     mDirichletSetTag = 0;
00105     mNeumannSetTag   = 0;
00106     mHasMidNodesTag  = 0;
00107     mGlobalIdTag     = 0;
00108     mNameTag         = 0;
00109 
00110     //! Get and cache predefined tag handles
00111     const int negone = -1;
00112     ErrorCode result = impl->tag_get_handle( MATERIAL_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mMaterialSetTag,
00113                                              MB_TAG_CREAT | MB_TAG_SPARSE, &negone );MB_CHK_SET_ERR_RET( result, "Failed to get MATERIAL_SET tag" );
00114 
00115     result = impl->tag_get_handle( DIRICHLET_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mDirichletSetTag,
00116                                    MB_TAG_CREAT | MB_TAG_SPARSE, &negone );MB_CHK_SET_ERR_RET( result, "Failed to get DIRICHLET_SET tag" );
00117 
00118     result = impl->tag_get_handle( NEUMANN_SET_TAG_NAME, 1, MB_TYPE_INTEGER, mNeumannSetTag,
00119                                    MB_TAG_CREAT | MB_TAG_SPARSE, &negone );MB_CHK_SET_ERR_RET( result, "Failed to get NEUMANN_SET tag" );
00120 
00121     const int negonearr[] = { -1, -1, -1, -1 };
00122     result                = impl->tag_get_handle( HAS_MID_NODES_TAG_NAME, 4, MB_TYPE_INTEGER, mHasMidNodesTag,
00123                                                   MB_TAG_CREAT | MB_TAG_SPARSE, negonearr );MB_CHK_SET_ERR_RET( result, "Failed to get HAS_MID_NODES tag" );
00124 
00125     mGlobalIdTag = impl->globalId_tag();
00126 
00127     result =
00128         impl->tag_get_handle( NAME_TAG_NAME, NAME_TAG_SIZE, MB_TYPE_OPAQUE, mNameTag, MB_TAG_CREAT | MB_TAG_SPARSE );MB_CHK_SET_ERR_RET( result, "Failed to get NAME tag" );
00129 }
00130 
00131 ReadCCMIO::~ReadCCMIO()
00132 {
00133     mbImpl->release_interface( readMeshIface );
00134 }
00135 
00136 ErrorCode ReadCCMIO::load_file( const char* file_name,
00137                                 const EntityHandle* file_set,
00138                                 const FileOptions& /* opts */,
00139                                 const ReaderIface::SubsetList* subset_list,
00140                                 const Tag* /* file_id_tag */ )
00141 {
00142     CCMIOID rootID, problemID, stateID, processorID, verticesID, topologyID, solutionID;
00143     CCMIOError error = kCCMIONoErr;
00144 
00145     if( subset_list )
00146     {
00147         MB_SET_ERR( MB_UNSUPPORTED_OPERATION, "Reading subset of files not supported for CCMOI data" );
00148     }
00149 
00150     CCMIOOpenFile( &error, file_name, kCCMIORead, &rootID );
00151     CHK_SET_CCMERR( error, "Problem opening file" );
00152 
00153     // Get the file state
00154     ErrorCode rval = get_state( rootID, problemID, stateID );MB_CHK_SET_ERR( rval, "Failed to get state" );
00155 
00156     // Get processors
00157     std::vector< CCMIOSize_t > procs;
00158     bool has_solution = false;
00159     rval              = get_processors( stateID, processorID, verticesID, topologyID, solutionID, procs, has_solution );MB_CHK_SET_ERR( rval, "Failed to get processors" );
00160 
00161     std::vector< CCMIOSize_t >::iterator vit;
00162     Range new_ents, *new_ents_ptr = NULL;
00163     if( file_set ) new_ents_ptr = &new_ents;
00164 
00165     for( vit = procs.begin(); vit != procs.end(); ++vit )
00166     {
00167         rval = read_processor( stateID, problemID, processorID, verticesID, topologyID, *vit, new_ents_ptr );MB_CHK_SET_ERR( rval, "Failed to read processors" );
00168     }
00169 
00170     // Load some meta-data
00171     rval = load_metadata( rootID, problemID, stateID, processorID, file_set );MB_CHK_SET_ERR( rval, "Failed to load some meta-data" );
00172 
00173     // Now, put all this into the file set, if there is one
00174     if( file_set )
00175     {
00176         rval = mbImpl->add_entities( *file_set, new_ents );MB_CHK_SET_ERR( rval, "Failed to add new entities to file set" );
00177     }
00178 
00179     return rval;
00180 }
00181 
00182 ErrorCode ReadCCMIO::get_state( CCMIOID rootID, CCMIOID& problemID, CCMIOID& stateID )
00183 {
00184     CCMIOError error = kCCMIONoErr;
00185 
00186     // First try default
00187     CCMIOGetState( &error, rootID, "default", &problemID, &stateID );
00188     if( kCCMIONoErr != error )
00189     {
00190         CCMIOSize_t i        = CCMIOSIZEC( 0 );
00191         CCMIOError tmp_error = kCCMIONoErr;
00192         CCMIONextEntity( &tmp_error, rootID, kCCMIOState, &i, &stateID );
00193         if( kCCMIONoErr == tmp_error ) CCMIONextEntity( &error, rootID, kCCMIOProblemDescription, &i, &problemID );
00194     }
00195     CHK_SET_CCMERR( error, "Couldn't find state" );
00196 
00197     return MB_SUCCESS;
00198 }
00199 
00200 ErrorCode ReadCCMIO::load_metadata( CCMIOID rootID,
00201                                     CCMIOID problemID,
00202                                     CCMIOID /* stateID */,
00203                                     CCMIOID processorID,
00204                                     const EntityHandle* file_set )
00205 {
00206     // Read the simulation title.
00207     CCMIOError error = kCCMIONoErr;
00208     ErrorCode rval   = MB_SUCCESS;
00209     CCMIONode rootNode;
00210     if( kCCMIONoErr == CCMIOGetEntityNode( &error, rootID, &rootNode ) )
00211     {
00212         char* name = NULL;
00213         CCMIOGetTitle( &error, rootNode, &name );
00214 
00215         if( NULL != name && strlen( name ) != 0 )
00216         {
00217             // Make a tag for it and tag the read set
00218             Tag simname;
00219             rval = mbImpl->tag_get_handle( "Title", strlen( name ), MB_TYPE_OPAQUE, simname,
00220                                            MB_TAG_CREAT | MB_TAG_SPARSE );MB_CHK_SET_ERR( rval, "Simulation name tag not found or created" );
00221             EntityHandle set = file_set ? *file_set : 0;
00222             rval             = mbImpl->tag_set_data( simname, &set, 1, name );MB_CHK_SET_ERR( rval, "Problem setting simulation name tag" );
00223         }
00224         if( name ) free( name );
00225     }
00226 
00227     // Creating program
00228     EntityHandle dumh = ( file_set ? *file_set : 0 );
00229     rval              = get_str_option( "CreatingProgram", dumh, mCreatingProgramTag, processorID );MB_CHK_SET_ERR( rval, "Trouble getting CreatingProgram tag" );
00230 
00231     rval = load_matset_data( problemID );MB_CHK_SET_ERR( rval, "Failure loading matset data" );
00232 
00233     rval = load_neuset_data( problemID );MB_CHK_SET_ERR( rval, "Failure loading neuset data" );
00234 
00235     return rval;
00236 }
00237 
00238 ErrorCode ReadCCMIO::load_matset_data( CCMIOID problemID )
00239 {
00240     // Make sure there are matsets
00241     if( newMatsets.empty() ) return MB_SUCCESS;
00242 
00243     // ... walk through each cell type
00244     CCMIOSize_t i = CCMIOSIZEC( 0 );
00245     CCMIOID next;
00246     CCMIOError error = kCCMIONoErr;
00247 
00248     while( CCMIONextEntity( NULL, problemID, kCCMIOCellType, &i, &next ) == kCCMIONoErr )
00249     {
00250         // Get index, corresponding set, and label with material set tag
00251         int mindex;
00252         CCMIOGetEntityIndex( &error, next, &mindex );
00253         std::map< int, EntityHandle >::iterator mit = newMatsets.find( mindex );
00254         if( mit == newMatsets.end() )
00255             // No actual faces for this matset; continue to next
00256             continue;
00257 
00258         EntityHandle dum_ent = mit->second;
00259         ErrorCode rval       = mbImpl->tag_set_data( mMaterialSetTag, &dum_ent, 1, &mindex );MB_CHK_SET_ERR( rval, "Trouble setting material set tag" );
00260 
00261         // Set name
00262         CCMIOSize_t len;
00263         CCMIOEntityLabel( &error, next, &len, NULL );
00264         std::vector< char > opt_string2( GETINT32( len ) + 1, '\0' );
00265         CCMIOEntityLabel( &error, next, NULL, &opt_string2[0] );
00266         if( opt_string2.size() >= NAME_TAG_SIZE )
00267             opt_string2[NAME_TAG_SIZE - 1] = '\0';
00268         else
00269             ( opt_string2.resize( NAME_TAG_SIZE, '\0' ) );
00270         rval = mbImpl->tag_set_data( mNameTag, &dum_ent, 1, &opt_string2[0] );MB_CHK_SET_ERR( rval, "Trouble setting name tag for material set" );
00271 
00272         // Material id
00273         rval = get_int_option( "MaterialId", dum_ent, mMaterialIdTag, next );MB_CHK_SET_ERR( rval, "Trouble getting MaterialId tag" );
00274 
00275         rval = get_str_option( "MaterialType", dum_ent, mMaterialTypeTag, next );MB_CHK_SET_ERR( rval, "Trouble getting MaterialType tag" );
00276 
00277         rval = get_int_option( "Radiation", dum_ent, mRadiationTag, next );MB_CHK_SET_ERR( rval, "Trouble getting Radiation option" );
00278 
00279         rval = get_int_option( "PorosityId", dum_ent, mPorosityIdTag, next );MB_CHK_SET_ERR( rval, "Trouble getting PorosityId option" );
00280 
00281         rval = get_int_option( "SpinId", dum_ent, mSpinIdTag, next );MB_CHK_SET_ERR( rval, "Trouble getting SpinId option" );
00282 
00283         rval = get_int_option( "GroupId", dum_ent, mGroupIdTag, next );MB_CHK_SET_ERR( rval, "Trouble getting GroupId option" );
00284 
00285         rval = get_int_option( "ColorIdx", dum_ent, mColorIdxTag, next );MB_CHK_SET_ERR( rval, "Trouble getting ColorIdx option" );
00286 
00287         rval = get_int_option( "ProcessorId", dum_ent, mProcessorIdTag, next );MB_CHK_SET_ERR( rval, "Trouble getting ProcessorId option" );
00288 
00289         rval = get_int_option( "LightMaterial", dum_ent, mLightMaterialTag, next );MB_CHK_SET_ERR( rval, "Trouble getting LightMaterial option" );
00290 
00291         rval = get_int_option( "FreeSurfaceMaterial", dum_ent, mFreeSurfaceMaterialTag, next );MB_CHK_SET_ERR( rval, "Trouble getting FreeSurfaceMaterial option" );
00292 
00293         rval = get_dbl_option( "Thickness", dum_ent, mThicknessTag, next );MB_CHK_SET_ERR( rval, "Trouble getting Thickness option" );
00294     }
00295 
00296     return MB_SUCCESS;
00297 }
00298 
00299 ErrorCode ReadCCMIO::get_int_option( const char* opt_str, EntityHandle seth, Tag& tag, CCMIOID node )
00300 {
00301     int idum;
00302     ErrorCode rval;
00303     if( kCCMIONoErr == CCMIOReadOpti( NULL, node, opt_str, &idum ) )
00304     {
00305         if( !tag )
00306         {
00307             rval = mbImpl->tag_get_handle( opt_str, 1, MB_TYPE_INTEGER, tag, MB_TAG_SPARSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Failed to get tag handle" );
00308         }
00309 
00310         rval = mbImpl->tag_set_data( tag, &seth, 1, &idum );MB_CHK_SET_ERR( rval, "Failed to set tag data" );
00311     }
00312 
00313     return MB_SUCCESS;
00314 }
00315 
00316 ErrorCode ReadCCMIO::get_dbl_option( const char* opt_str, EntityHandle seth, Tag& tag, CCMIOID node )
00317 {
00318     float fdum;
00319     if( kCCMIONoErr == CCMIOReadOptf( NULL, node, opt_str, &fdum ) )
00320     {
00321         ErrorCode rval;
00322         if( !tag )
00323         {
00324             rval = mbImpl->tag_get_handle( opt_str, 1, MB_TYPE_DOUBLE, tag, MB_TAG_SPARSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Failed to get tag handle" );
00325         }
00326 
00327         double dum_dbl = fdum;
00328         rval           = mbImpl->tag_set_data( tag, &seth, 1, &dum_dbl );MB_CHK_SET_ERR( rval, "Failed to set tag data" );
00329     }
00330 
00331     return MB_SUCCESS;
00332 }
00333 
00334 ErrorCode ReadCCMIO::get_str_option( const char* opt_str,
00335                                      EntityHandle seth,
00336                                      Tag& tag,
00337                                      CCMIOID node,
00338                                      const char* other_tag_name )
00339 {
00340     int len;
00341     CCMIOError error = kCCMIONoErr;
00342     std::vector< char > opt_string;
00343     if( kCCMIONoErr != CCMIOReadOptstr( NULL, node, opt_str, &len, NULL ) ) return MB_SUCCESS;
00344 
00345     opt_string.resize( len );
00346     CCMIOReadOptstr( &error, node, opt_str, &len, &opt_string[0] );
00347     ErrorCode rval = MB_SUCCESS;
00348     if( !tag )
00349     {
00350         rval = mbImpl->tag_get_handle( other_tag_name ? other_tag_name : opt_str, NAME_TAG_SIZE, MB_TYPE_OPAQUE, tag,
00351                                        MB_TAG_SPARSE | MB_TAG_CREAT );MB_CHK_SET_ERR( rval, "Failed to get tag handle" );
00352     }
00353 
00354     if( opt_string.size() > NAME_TAG_SIZE )
00355         opt_string[NAME_TAG_SIZE - 1] = '\0';
00356     else
00357         ( opt_string.resize( NAME_TAG_SIZE, '\0' ) );
00358 
00359     rval = mbImpl->tag_set_data( tag, &seth, 1, &opt_string[0] );MB_CHK_SET_ERR( rval, "Failed to set tag data" );
00360 
00361     return MB_SUCCESS;
00362 }
00363 
00364 ErrorCode ReadCCMIO::load_neuset_data( CCMIOID problemID )
00365 {
00366     CCMIOSize_t i = CCMIOSIZEC( 0 );
00367     CCMIOID next;
00368 
00369     // Make sure there are matsets
00370     if( newNeusets.empty() ) return MB_SUCCESS;
00371 
00372     while( CCMIONextEntity( NULL, problemID, kCCMIOBoundaryRegion, &i, &next ) == kCCMIONoErr )
00373     {
00374         // Get index, corresponding set, and label with neumann set tag
00375         int mindex;
00376         CCMIOError error = kCCMIONoErr;
00377         CCMIOGetEntityIndex( &error, next, &mindex );
00378         std::map< int, EntityHandle >::iterator mit = newNeusets.find( mindex );
00379         if( mit == newNeusets.end() )
00380             // No actual faces for this neuset; continue to next
00381             continue;
00382 
00383         EntityHandle dum_ent = mit->second;
00384         ErrorCode rval       = mbImpl->tag_set_data( mNeumannSetTag, &dum_ent, 1, &mindex );MB_CHK_SET_ERR( rval, "Trouble setting neumann set tag" );
00385 
00386         // Set name
00387         rval = get_str_option( "BoundaryName", dum_ent, mNameTag, next, NAME_TAG_NAME );MB_CHK_SET_ERR( rval, "Trouble creating BoundaryName tag" );
00388 
00389         // BoundaryType
00390         rval = get_str_option( "BoundaryType", dum_ent, mBoundaryTypeTag, next );MB_CHK_SET_ERR( rval, "Trouble creating BoundaryType tag" );
00391 
00392         // ProstarRegionNumber
00393         rval = get_int_option( "ProstarRegionNumber", dum_ent, mProstarRegionNumberTag, next );MB_CHK_SET_ERR( rval, "Trouble creating ProstarRegionNumber tag" );
00394     }
00395 
00396     return MB_SUCCESS;
00397 }
00398 
00399 ErrorCode ReadCCMIO::read_processor( CCMIOID /* stateID */,
00400                                      CCMIOID problemID,
00401                                      CCMIOID processorID,
00402                                      CCMIOID verticesID,
00403                                      CCMIOID topologyID,
00404                                      CCMIOSize_t proc,
00405                                      Range* new_ents )
00406 {
00407     ErrorCode rval;
00408 
00409     // vert_map fields: s: none, i: gid, ul: vert handle, r: none
00410     // TupleList vert_map(0, 1, 1, 0, 0);
00411     TupleList vert_map;
00412     rval = read_vertices( proc, processorID, verticesID, topologyID, new_ents, vert_map );MB_CHK_SET_ERR( rval, "Failed to read vertices" );
00413 
00414     rval = read_cells( proc, problemID, verticesID, topologyID, vert_map, new_ents );MB_CHK_SET_ERR( rval, "Failed to read cells" );
00415 
00416     return rval;
00417 }
00418 
00419 ErrorCode ReadCCMIO::read_cells( CCMIOSize_t /* proc */,
00420                                  CCMIOID problemID,
00421                                  CCMIOID /* verticesID */,
00422                                  CCMIOID topologyID,
00423                                  TupleList& vert_map,
00424                                  Range* new_ents )
00425 {
00426     // Read the faces.
00427     // face_map fields: s:forward/reverse, i: cell id, ul: face handle, r: none
00428     ErrorCode rval;
00429 #ifdef TUPLE_LIST
00430     TupleList face_map( 1, 1, 1, 0, 0 );
00431 #else
00432     TupleList face_map;
00433     SenseList sense_map;
00434 #endif
00435     rval = read_all_faces( topologyID, vert_map, face_map,
00436 #ifndef TUPLE_LIST
00437                            sense_map,
00438 #endif
00439                            new_ents );MB_CHK_SET_ERR( rval, "Failed to read all cells" );
00440 
00441     // Read the cell topology types, if any exist in the file
00442     std::map< int, int > cell_topo_types;
00443     rval = read_topology_types( topologyID, cell_topo_types );MB_CHK_SET_ERR( rval, "Problem reading cell topo types" );
00444 
00445     // Now construct the cells; sort the face map by cell ids first
00446 #ifdef TUPLE_LIST
00447     rval = face_map.sort( 1 );MB_CHK_SET_ERR( rval, "Couldn't sort face map by cell id" );
00448 #endif
00449     std::vector< EntityHandle > new_cells;
00450     rval = construct_cells( face_map,
00451 #ifndef TUPLE_LIST
00452                             sense_map,
00453 #endif
00454                             vert_map, cell_topo_types, new_cells );MB_CHK_SET_ERR( rval, "Failed to construct cells" );
00455     if( new_ents )
00456     {
00457         Range::iterator rit = new_ents->end();
00458         std::vector< EntityHandle >::reverse_iterator vit;
00459         for( vit = new_cells.rbegin(); vit != new_cells.rend(); ++vit )
00460             rit = new_ents->insert( rit, *vit );
00461     }
00462 
00463     rval = read_gids_and_types( problemID, topologyID, new_cells );MB_CHK_SET_ERR( rval, "Failed to read gids and types" );
00464 
00465     return MB_SUCCESS;
00466 }
00467 
00468 ErrorCode ReadCCMIO::read_topology_types( CCMIOID& topologyID, std::map< int, int >& cell_topo_types )
00469 {
00470     CCMIOError error = kCCMIONoErr;
00471     CCMIOID cellID, mapID;
00472     CCMIOSize_t ncells;
00473     CCMIOGetEntity( &error, topologyID, kCCMIOCells, 0, &cellID );
00474     CCMIOEntitySize( &error, cellID, &ncells, NULL );
00475     int num_cells = GETINT32( ncells );
00476 
00477     // First, do a dummy read to see if we even have topo types in this mesh
00478     int dum_int;
00479     CCMIOReadOpt1i( &error, cellID, "CellTopologyType", &dum_int, CCMIOINDEXC( kCCMIOStart ),
00480                     CCMIOINDEXC( kCCMIOStart ) + 1 );
00481     if( kCCMIONoErr != error ) return MB_SUCCESS;
00482 
00483     // OK, we have topo types; first get the map node
00484     std::vector< int > dum_ints( num_cells );
00485     CCMIOReadCells( &error, cellID, &mapID, &dum_ints[0], CCMIOINDEXC( kCCMIOStart ), CCMIOINDEXC( kCCMIOStart ) + 1 );
00486     CHK_SET_CCMERR( error, "Failed to get the map node" );
00487 
00488     // Now read the map
00489     CCMIOReadMap( &error, mapID, &dum_ints[0], CCMIOINDEXC( kCCMIOStart ), CCMIOINDEXC( kCCMIOEnd ) );
00490     CHK_SET_CCMERR( error, "Failed to get cell ids" );
00491     int i;
00492     for( i = 0; i < num_cells; i++ )
00493         cell_topo_types[dum_ints[i]] = 0;
00494 
00495     // Now read the cell topo types for real, reusing cell_topo_types
00496     std::vector< int > topo_types( num_cells );
00497     CCMIOReadOpt1i( &error, cellID, "CellTopologyType", &topo_types[0], CCMIOINDEXC( kCCMIOStart ),
00498                     CCMIOINDEXC( kCCMIOEnd ) );
00499     CHK_SET_CCMERR( error, "Failed to get cell topo types" );
00500     for( i = 0; i < num_cells; i++ )
00501         cell_topo_types[dum_ints[i]] = topo_types[i];
00502 
00503     return MB_SUCCESS;
00504 }
00505 
00506 ErrorCode ReadCCMIO::read_gids_and_types( CCMIOID /* problemID */,
00507                                           CCMIOID topologyID,
00508                                           std::vector< EntityHandle >& cells )
00509 {
00510     // Get the cells entity and number of cells
00511     CCMIOSize_t dum_cells;
00512     int num_cells;
00513     CCMIOError error = kCCMIONoErr;
00514     CCMIOID cellsID, mapID;
00515     CCMIOGetEntity( &error, topologyID, kCCMIOCells, 0, &cellsID );
00516     CCMIOEntitySize( &error, cellsID, &dum_cells, NULL );
00517     num_cells = GETINT32( dum_cells );
00518 
00519     // Check the number of cells against how many are in the cell array
00520     if( num_cells != (int)cells.size() ) MB_SET_ERR( MB_FAILURE, "Number of cells doesn't agree" );
00521 
00522     // Read the gid map and set global ids
00523     std::vector< int > cell_gids( num_cells );
00524     CCMIOReadCells( &error, cellsID, &mapID, NULL, CCMIOINDEXC( kCCMIOStart ), CCMIOINDEXC( kCCMIOEnd ) );
00525     CHK_SET_CCMERR( error, "Couldn't read cells" );
00526     CCMIOReadMap( &error, mapID, &cell_gids[0], CCMIOINDEXC( kCCMIOStart ), CCMIOINDEXC( kCCMIOEnd ) );
00527     CHK_SET_CCMERR( error, "Couldn't read cell id map" );
00528 
00529     ErrorCode rval = mbImpl->tag_set_data( mGlobalIdTag, &cells[0], cells.size(), &cell_gids[0] );MB_CHK_SET_ERR( rval, "Couldn't set gids tag" );
00530 
00531     // Now read cell material types; reuse cell_gids
00532     CCMIOReadCells( &error, cellsID, NULL, &cell_gids[0], CCMIOINDEXC( kCCMIOStart ), CCMIOINDEXC( kCCMIOEnd ) );
00533     CHK_SET_CCMERR( error, "Trouble reading cell types" );
00534 
00535     // Create the matsets
00536     std::map< int, Range > matset_ents;
00537     for( int i = 0; i < num_cells; i++ )
00538         matset_ents[cell_gids[i]].insert( cells[i] );
00539 
00540     for( std::map< int, Range >::iterator mit = matset_ents.begin(); mit != matset_ents.end(); ++mit )
00541     {
00542         EntityHandle matset;
00543         rval = mbImpl->create_meshset( MESHSET_SET, matset );MB_CHK_SET_ERR( rval, "Couldn't create material set" );
00544         newMatsets[mit->first] = matset;
00545 
00546         rval = mbImpl->add_entities( matset, mit->second );MB_CHK_SET_ERR( rval, "Couldn't add entities to material set" );
00547     }
00548 
00549     return MB_SUCCESS;
00550 }
00551 
00552 ErrorCode ReadCCMIO::construct_cells( TupleList& face_map,
00553 #ifndef TUPLE_LIST
00554                                       SenseList& sense_map,
00555 #endif
00556                                       TupleList& /* vert_map */,
00557                                       std::map< int, int >& cell_topo_types,
00558                                       std::vector< EntityHandle >& new_cells )
00559 {
00560     std::vector< EntityHandle > facehs;
00561     std::vector< int > senses;
00562     EntityHandle cell;
00563     ErrorCode tmp_rval, rval = MB_SUCCESS;
00564     EntityType this_type = MBMAXTYPE;
00565     bool has_mid_nodes   = false;
00566 #ifdef TUPLE_LIST
00567     unsigned int i = 0;
00568     while( i < face_map.n )
00569     {
00570         // Pull out face handles bounding the same cell
00571         facehs.clear();
00572         int this_id        = face_map.get_int( i );
00573         unsigned int inext = i;
00574         while( face_map.get_int( inext ) == this_id && inext <= face_map.n )
00575         {
00576             inext++;
00577             EntityHandle face = face_map.get_ulong( inext );
00578             facehs.push_back( face );
00579             senses.push_back( face_map.get_short( inext ) );
00580         }
00581         this_type     = MBMAXTYPE;
00582         has_mid_nodes = false;
00583 #else
00584     std::map< int, std::vector< EntityHandle > >::iterator fmit;
00585     std::map< int, std::vector< int > >::iterator smit;
00586     std::map< int, int >::iterator typeit;
00587     for( fmit = face_map.begin(), smit = sense_map.begin(); fmit != face_map.end(); ++fmit, ++smit )
00588     {
00589         // Pull out face handles bounding the same cell
00590         facehs.clear();
00591         int this_id = ( *fmit ).first;
00592         facehs      = ( *fmit ).second;
00593         senses.clear();
00594         senses = ( *smit ).second;
00595         typeit = cell_topo_types.find( this_id );
00596         if( typeit != cell_topo_types.end() )
00597         {
00598             rval = ccmio_to_moab_type( typeit->second, this_type, has_mid_nodes );
00599         }
00600         else
00601         {
00602             this_type     = MBMAXTYPE;
00603             has_mid_nodes = false;
00604         }
00605 #endif
00606         tmp_rval = create_cell_from_faces( facehs, senses, this_type, has_mid_nodes, cell );
00607         if( MB_SUCCESS != tmp_rval )
00608             rval = tmp_rval;
00609         else
00610         {
00611             new_cells.push_back( cell );
00612             // Tag cell with global id
00613             tmp_rval = mbImpl->tag_set_data( mGlobalIdTag, &cell, 1, &this_id );
00614             if( MB_SUCCESS != tmp_rval ) rval = tmp_rval;
00615         }
00616     }
00617 
00618     return rval;
00619 }
00620 
00621 ErrorCode ReadCCMIO::ccmio_to_moab_type( int ccm_type, EntityType& moab_type, bool& has_mid_nodes )
00622 {
00623     switch( ccm_type )
00624     {
00625         case 1:
00626             moab_type = MBVERTEX;
00627             break;
00628         case 2:
00629         case 28:
00630             moab_type = MBEDGE;
00631             break;
00632         case 29:
00633             moab_type = MBMAXTYPE;
00634             break;
00635         case 3:
00636         case 4:
00637             moab_type = MBQUAD;
00638             break;
00639         case 11:
00640         case 21:
00641             moab_type = MBHEX;
00642             break;
00643         case 12:
00644         case 22:
00645             moab_type = MBPRISM;
00646             break;
00647         case 13:
00648         case 23:
00649             moab_type = MBTET;
00650             break;
00651         case 14:
00652         case 24:
00653             moab_type = MBPYRAMID;
00654             break;
00655         case 255:
00656             moab_type = MBPOLYHEDRON;
00657             break;
00658         default:
00659             moab_type = MBMAXTYPE;
00660     }
00661 
00662     switch( ccm_type )
00663     {
00664         case 28:
00665         case 4:
00666         case 21:
00667         case 22:
00668         case 23:
00669         case 24:
00670             has_mid_nodes = true;
00671             break;
00672         default:
00673             break;
00674     }
00675 
00676     return MB_SUCCESS;
00677 }
00678 
00679 ErrorCode ReadCCMIO::create_cell_from_faces( std::vector< EntityHandle >& facehs,
00680                                              std::vector< int >& senses,
00681                                              EntityType this_type,
00682                                              bool /* has_mid_nodes */,
00683                                              EntityHandle& cell )
00684 {
00685     ErrorCode rval;
00686 
00687     // Test up front to see if they're one type
00688     EntityType face_type = mbImpl->type_from_handle( facehs[0] );
00689     bool same_type       = true;
00690     for( std::vector< EntityHandle >::iterator vit = facehs.begin(); vit != facehs.end(); ++vit )
00691     {
00692         if( face_type != mbImpl->type_from_handle( *vit ) )
00693         {
00694             same_type = false;
00695             break;
00696         }
00697     }
00698 
00699     std::vector< EntityHandle > verts;
00700     EntityType input_type = this_type;
00701     std::vector< EntityHandle > storage;
00702     MeshTopoUtil mtu( mbImpl );
00703 
00704     // Preset this to maxtype, so we get an affirmative choice in loop
00705     this_type = MBMAXTYPE;
00706 
00707     if( ( MBTET == input_type || MBMAXTYPE == input_type ) && same_type && face_type == MBTRI && facehs.size() == 4 )
00708     {
00709         // Try to get proper connectivity for tet
00710 
00711         // Get connectivity of first face, and reverse it if sense is forward, since
00712         // base face always points into entity
00713         rval = mbImpl->get_connectivity( &facehs[0], 1, verts );MB_CHK_SET_ERR( rval, "Couldn't get connectivity" );
00714         if( senses[0] > 0 ) std::reverse( verts.begin(), verts.end() );
00715 
00716         // Get the 4th vertex through the next tri
00717         const EntityHandle* conn;
00718         int conn_size;
00719         rval = mbImpl->get_connectivity( facehs[1], conn, conn_size, true, &storage );MB_CHK_SET_ERR( rval, "Couldn't get connectivity" );
00720         int i = 0;
00721         while( std::find( verts.begin(), verts.end(), conn[i] ) != verts.end() && i < conn_size )
00722             i++;
00723 
00724         // If i is not at the end of the verts, found the apex; otherwise fall back to polyhedron
00725         if( conn_size != i )
00726         {
00727             this_type = MBTET;
00728             verts.push_back( conn[i] );
00729         }
00730     }
00731     else if( ( MBHEX == input_type || MBMAXTYPE == input_type ) && same_type && MBQUAD == face_type &&
00732              facehs.size() == 6 )
00733     {
00734         // Build hex from quads
00735         // Algorithm:
00736         // - verts = vertices from 1st quad
00737         // - Find quad q1 sharing verts[0] and verts[1]
00738         // - Find quad q2 sharing other 2 verts in q1
00739         // - Find v1 = opposite vert from verts[1] in q1 , v2 = opposite from verts[0]
00740         // - Get i = offset of v1 in verts2 of q2, rotate verts2 by i
00741         // - If verts2[(i + 1) % 4] != v2, flip verts2 by switching verts2[1] and verts2[3]
00742         // - append verts2 to verts
00743 
00744         // Get the other vertices for this hex; need to find the quad with no common vertices
00745         Range tmp_faces, tmp_verts;
00746         // Get connectivity of first face, and reverse it if sense is forward, since
00747         // base face always points into entity
00748         rval = mbImpl->get_connectivity( &facehs[0], 1, verts );MB_CHK_SET_ERR( rval, "Couldn't get connectivity" );
00749         if( senses[0] > 0 ) std::reverse( verts.begin(), verts.end() );
00750 
00751         // Get q1, which shares 2 vertices with q0
00752         std::copy( facehs.begin(), facehs.end(), range_inserter( tmp_faces ) );
00753         rval = mbImpl->get_adjacencies( &verts[0], 2, 2, false, tmp_faces );
00754         if( MB_SUCCESS != rval || tmp_faces.size() != 2 ) MB_SET_ERR( MB_FAILURE, "Couldn't get adj face" );
00755         tmp_faces.erase( facehs[0] );
00756         EntityHandle q1 = *tmp_faces.begin();
00757         // Get other 2 verts of q1
00758         rval = mbImpl->get_connectivity( &q1, 1, tmp_verts );MB_CHK_SET_ERR( rval, "Couldn't get adj verts" );
00759         tmp_verts.erase( verts[0] );
00760         tmp_verts.erase( verts[1] );
00761         // Get q2
00762         std::copy( facehs.begin(), facehs.end(), range_inserter( tmp_faces ) );
00763         rval = mbImpl->get_adjacencies( tmp_verts, 2, false, tmp_faces );
00764         if( MB_SUCCESS != rval || tmp_faces.size() != 2 ) MB_SET_ERR( MB_FAILURE, "Couldn't get adj face" );
00765         tmp_faces.erase( q1 );
00766         EntityHandle q2 = *tmp_faces.begin();
00767         // Get verts in q2
00768         rval = mbImpl->get_connectivity( &q2, 1, storage );MB_CHK_SET_ERR( rval, "Couldn't get adj vertices" );
00769 
00770         // Get verts in q1 opposite from v[1] and v[0] in q0
00771         EntityHandle v0 = 0, v1 = 0;
00772         rval = mtu.opposite_entity( q1, verts[1], v0 );MB_CHK_SET_ERR( rval, "Couldn't get the opposite side entity" );
00773         rval = mtu.opposite_entity( q1, verts[0], v1 );MB_CHK_SET_ERR( rval, "Couldn't get the opposite side entity" );
00774         if( v0 && v1 )
00775         {
00776             // Offset of v0 in q2, then rotate and flip
00777             unsigned int ioff = std::find( storage.begin(), storage.end(), v0 ) - storage.begin();
00778             if( 4 == ioff ) MB_SET_ERR( MB_FAILURE, "Trouble finding offset" );
00779 
00780             if( storage[( ioff + 1 ) % 4] != v1 )
00781             {
00782                 std::reverse( storage.begin(), storage.end() );
00783                 ioff = std::find( storage.begin(), storage.end(), v0 ) - storage.begin();
00784             }
00785             if( 0 != ioff ) std::rotate( storage.begin(), storage.begin() + ioff, storage.end() );
00786 
00787             // Copy into verts, and make hex
00788             std::copy( storage.begin(), storage.end(), std::back_inserter( verts ) );
00789             this_type = MBHEX;
00790         }
00791     }
00792 
00793     if( MBMAXTYPE == this_type && facehs.size() == 5 )
00794     {
00795         // Some preliminaries
00796         std::vector< EntityHandle > tris, quads;
00797         for( unsigned int i = 0; i < 5; i++ )
00798         {
00799             if( MBTRI == mbImpl->type_from_handle( facehs[i] ) )
00800                 tris.push_back( facehs[i] );
00801             else if( MBQUAD == mbImpl->type_from_handle( facehs[i] ) )
00802                 quads.push_back( facehs[i] );
00803         }
00804 
00805         // Check for prisms
00806         if( 2 == tris.size() && 3 == quads.size() )
00807         {
00808             // OK, we have the right number of tris and quads; try to find the proper verts
00809 
00810             // Get connectivity of first tri, and reverse if necessary
00811             int index = std::find( facehs.begin(), facehs.end(), tris[0] ) - facehs.begin();
00812             rval      = mbImpl->get_connectivity( &tris[0], 1, verts );MB_CHK_SET_ERR( rval, "Couldn't get connectivity" );
00813             if( senses[index] > 0 ) std::reverse( verts.begin(), verts.end() );
00814 
00815             // Now align vertices of other tri, through a quad, similar to how we did hexes
00816             // Get q1, which shares 2 vertices with t0
00817             Range tmp_faces, tmp_verts;
00818             std::copy( facehs.begin(), facehs.end(), range_inserter( tmp_faces ) );
00819             rval = mbImpl->get_adjacencies( &verts[0], 2, 2, false, tmp_faces );
00820             if( MB_SUCCESS != rval || tmp_faces.size() != 2 ) MB_SET_ERR( MB_FAILURE, "Couldn't get adj face" );
00821             tmp_faces.erase( tris[0] );
00822             EntityHandle q1 = *tmp_faces.begin();
00823             // Get verts in q1
00824             rval = mbImpl->get_connectivity( &q1, 1, storage );MB_CHK_SET_ERR( rval, "Couldn't get adj vertices" );
00825 
00826             // Get verts in q1 opposite from v[1] and v[0] in q0
00827             EntityHandle v0 = 0, v1 = 0;
00828             rval = mtu.opposite_entity( q1, verts[1], v0 );MB_CHK_SET_ERR( rval, "Couldn't get the opposite side entity" );
00829             rval = mtu.opposite_entity( q1, verts[0], v1 );MB_CHK_SET_ERR( rval, "Couldn't get the opposite side entity" );
00830             if( v0 && v1 )
00831             {
00832                 // Offset of v0 in t2, then rotate and flip
00833                 storage.clear();
00834                 rval = mbImpl->get_connectivity( &tris[1], 1, storage );MB_CHK_SET_ERR( rval, "Couldn't get connectivity" );
00835 
00836                 index = std::find( facehs.begin(), facehs.end(), tris[1] ) - facehs.begin();
00837                 if( senses[index] < 0 ) std::reverse( storage.begin(), storage.end() );
00838                 unsigned int ioff = std::find( storage.begin(), storage.end(), v0 ) - storage.begin();
00839                 if( 3 == ioff ) MB_SET_ERR( MB_FAILURE, "Trouble finding offset" );
00840                 for( unsigned int i = 0; i < 3; i++ )
00841                     verts.push_back( storage[( ioff + i ) % 3] );
00842 
00843                 this_type = MBPRISM;
00844             }
00845         }
00846         else if( tris.size() == 4 && quads.size() == 1 )
00847         {
00848             // Check for pyramid
00849             // Get connectivity of first tri, and reverse if necessary
00850             int index = std::find( facehs.begin(), facehs.end(), quads[0] ) - facehs.begin();
00851             rval      = mbImpl->get_connectivity( &quads[0], 1, verts );MB_CHK_SET_ERR( rval, "Couldn't get connectivity" );
00852             if( senses[index] > 0 ) std::reverse( verts.begin(), verts.end() );
00853 
00854             // Get apex node
00855             rval = mbImpl->get_connectivity( &tris[0], 1, storage );MB_CHK_SET_ERR( rval, "Couldn't get connectivity" );
00856             for( unsigned int i = 0; i < 3; i++ )
00857             {
00858                 if( std::find( verts.begin(), verts.end(), storage[i] ) == verts.end() )
00859                 {
00860                     verts.push_back( storage[i] );
00861                     break;
00862                 }
00863             }
00864 
00865             if( 5 == verts.size() ) this_type = MBPYRAMID;
00866         }
00867         else
00868         {
00869             // Dummy else clause to stop in debugger
00870             this_type = MBMAXTYPE;
00871         }
00872     }
00873 
00874     if( MBMAXTYPE != input_type && input_type != this_type && this_type != MBMAXTYPE )
00875         std::cerr << "Warning: types disagree (cell_topo_type = " << CN::EntityTypeName( input_type )
00876                   << ", faces indicate type " << CN::EntityTypeName( this_type ) << std::endl;
00877 
00878     if( MBMAXTYPE != input_type && this_type == MBMAXTYPE && input_type != MBPOLYHEDRON )
00879         std::cerr << "Warning: couldn't find proper connectivity for specified topo_type = "
00880                   << CN::EntityTypeName( input_type ) << std::endl;
00881 
00882     // Now make the element; if we fell back to polyhedron, use faces, otherwise use verts
00883     if( MBPOLYHEDRON == this_type || MBMAXTYPE == this_type )
00884     {
00885         rval = mbImpl->create_element( MBPOLYHEDRON, &facehs[0], facehs.size(), cell );MB_CHK_SET_ERR( rval, "create_element failed" );
00886     }
00887     else
00888     {
00889         rval = mbImpl->create_element( this_type, &verts[0], verts.size(), cell );MB_CHK_SET_ERR( rval, "create_element failed" );
00890     }
00891 
00892     return MB_SUCCESS;
00893 }
00894 
00895 ErrorCode ReadCCMIO::read_all_faces( CCMIOID topologyID,
00896                                      TupleList& vert_map,
00897                                      TupleList& face_map,
00898 #ifndef TUPLE_LIST
00899                                      SenseList& sense_map,
00900 #endif
00901                                      Range* new_faces )
00902 {
00903     CCMIOSize_t index;
00904     CCMIOID faceID;
00905     ErrorCode rval;
00906     CCMIOError error = kCCMIONoErr;
00907 
00908     // Get total # internal/bdy faces, size the face map accordingly
00909 #ifdef TUPLE_LIST
00910     index          = CCMIOSIZEC( 0 );
00911     int nbdy_faces = 0;
00912     CCMIOSize_t nf;
00913     error = kCCMIONoErr;
00914     while( kCCMIONoErr == CCMIONextEntity( NULL, topologyID, kCCMIOBoundaryFaces, &index, &faceID ) )
00915     {
00916         CCMIOEntitySize( &error, faceID, &nf, NULL );
00917         nbdy_faces += nf;
00918     }
00919     CCMIOGetEntity( &error, topologyID, kCCMIOInternalFaces, 0, &faceID );
00920     CCMIOEntitySize( &error, faceID, &nf, NULL );
00921 
00922     int nint_faces = nf;
00923     face_map.resize( 2 * nint_faces + nbdy_faces );
00924 #endif
00925 
00926     // Get multiple blocks of bdy faces
00927     index = CCMIOSIZEC( 0 );
00928     while( kCCMIONoErr == CCMIONextEntity( NULL, topologyID, kCCMIOBoundaryFaces, &index, &faceID ) )
00929     {
00930         rval = read_faces( faceID, kCCMIOBoundaryFaces, vert_map, face_map,
00931 #ifndef TUPLE_LIST
00932                            sense_map,
00933 #endif
00934                            new_faces );MB_CHK_SET_ERR( rval, "Trouble reading boundary faces" );
00935     }
00936 
00937     // Now get internal faces
00938     CCMIOGetEntity( &error, topologyID, kCCMIOInternalFaces, 0, &faceID );
00939     CHK_SET_CCMERR( error, "Couldn't get internal faces" );
00940 
00941     rval = read_faces( faceID, kCCMIOInternalFaces, vert_map, face_map,
00942 #ifndef TUPLE_LIST
00943                        sense_map,
00944 #endif
00945                        new_faces );MB_CHK_SET_ERR( rval, "Trouble reading internal faces" );
00946 
00947     return rval;
00948 }
00949 
00950 ErrorCode ReadCCMIO::read_faces( CCMIOID faceID,
00951                                  CCMIOEntity bdy_or_int,
00952                                  TupleList& vert_map,
00953                                  TupleList& face_map,
00954 #ifndef TUPLE_LIST
00955                                  SenseList& sense_map,
00956 #endif
00957                                  Range* new_faces )
00958 {
00959     if( kCCMIOInternalFaces != bdy_or_int && kCCMIOBoundaryFaces != bdy_or_int )
00960         MB_SET_ERR( MB_FAILURE, "Face type isn't boundary or internal" );
00961 
00962     CCMIOSize_t dum_faces;
00963     CCMIOError error = kCCMIONoErr;
00964     CCMIOEntitySize( &error, faceID, &dum_faces, NULL );
00965     int num_faces = GETINT32( dum_faces );
00966 
00967     // Get the size of the face connectivity array (not really a straight connect
00968     // array, has n, connect(n), ...)
00969     CCMIOSize_t farray_size = CCMIOSIZEC( 0 );
00970     CCMIOReadFaces( &error, faceID, bdy_or_int, NULL, &farray_size, NULL, CCMIOINDEXC( kCCMIOStart ),
00971                     CCMIOINDEXC( kCCMIOEnd ) );
00972     CHK_SET_CCMERR( error, "Trouble reading face connectivity length" );
00973 
00974     // Allocate vectors for holding farray and cells for each face; use new for finer
00975     // control of de-allocation
00976     int num_sides = ( kCCMIOInternalFaces == bdy_or_int ? 2 : 1 );
00977     int* farray   = new int[GETINT32( farray_size )];
00978 
00979     // Read farray and make the faces
00980     CCMIOID mapID;
00981     CCMIOReadFaces( &error, faceID, bdy_or_int, &mapID, NULL, farray, CCMIOINDEXC( kCCMIOStart ),
00982                     CCMIOINDEXC( kCCMIOEnd ) );
00983     CHK_SET_CCMERR( error, "Trouble reading face connectivity" );
00984 
00985     std::vector< EntityHandle > face_handles;
00986     ErrorCode rval = make_faces( farray, vert_map, face_handles, num_faces );MB_CHK_SET_ERR( rval, "Failed to make the faces" );
00987 
00988     // Read face cells and make tuples
00989     int* face_cells;
00990     if( num_sides * num_faces < farray_size )
00991         face_cells = new int[num_sides * num_faces];
00992     else
00993         face_cells = farray;
00994     CCMIOReadFaceCells( &error, faceID, bdy_or_int, face_cells, CCMIOINDEXC( kCCMIOStart ), CCMIOINDEXC( kCCMIOEnd ) );
00995     CHK_SET_CCMERR( error, "Trouble reading face cells" );
00996 
00997     int* tmp_ptr = face_cells;
00998     for( unsigned int i = 0; i < face_handles.size(); i++ )
00999     {
01000 #ifdef TUPLE_LIST
01001         short forward = 1, reverse = -1;
01002         face_map.push_back( &forward, tmp_ptr++, &face_handles[i], NULL );
01003         if( 2 == num_sides ) face_map.push_back( &reverse, tmp_ptr++, &face_handles[i], NULL );
01004 #else
01005         face_map[*tmp_ptr].push_back( face_handles[i] );
01006         sense_map[*tmp_ptr++].push_back( 1 );
01007         if( 2 == num_sides )
01008         {
01009             face_map[*tmp_ptr].push_back( face_handles[i] );
01010             sense_map[*tmp_ptr++].push_back( -1 );
01011         }
01012 #endif
01013     }
01014 
01015     // Now read & set face gids, reuse face_cells 'cuz we know it's big enough
01016     CCMIOReadMap( &error, mapID, face_cells, CCMIOINDEXC( kCCMIOStart ), CCMIOINDEXC( kCCMIOEnd ) );
01017     CHK_SET_CCMERR( error, "Trouble reading face gids" );
01018 
01019     rval = mbImpl->tag_set_data( mGlobalIdTag, &face_handles[0], face_handles.size(), face_cells );MB_CHK_SET_ERR( rval, "Couldn't set face global ids" );
01020 
01021     // Make a neumann set for these faces if they're all in a boundary face set
01022     if( kCCMIOBoundaryFaces == bdy_or_int )
01023     {
01024         EntityHandle neuset;
01025         rval = mbImpl->create_meshset( MESHSET_SET, neuset );MB_CHK_SET_ERR( rval, "Failed to create neumann set" );
01026 
01027         // Don't trust entity index passed in
01028         int index;
01029         CCMIOGetEntityIndex( &error, faceID, &index );
01030         newNeusets[index] = neuset;
01031 
01032         rval = mbImpl->add_entities( neuset, &face_handles[0], face_handles.size() );MB_CHK_SET_ERR( rval, "Failed to add faces to neumann set" );
01033 
01034         // Now tag as neumann set; will add id later
01035         int dum_val = 0;
01036         rval        = mbImpl->tag_set_data( mNeumannSetTag, &neuset, 1, &dum_val );MB_CHK_SET_ERR( rval, "Failed to tag neumann set" );
01037     }
01038 
01039     if( new_faces )
01040     {
01041         std::sort( face_handles.begin(), face_handles.end() );
01042         std::copy( face_handles.rbegin(), face_handles.rend(), range_inserter( *new_faces ) );
01043     }
01044 
01045     return MB_SUCCESS;
01046 }
01047 
01048 ErrorCode ReadCCMIO::make_faces( int* farray,
01049                                  TupleList& vert_map,
01050                                  std::vector< EntityHandle >& new_faces,
01051                                  int num_faces )
01052 {
01053     std::vector< EntityHandle > verts;
01054     ErrorCode tmp_rval = MB_SUCCESS, rval = MB_SUCCESS;
01055 
01056     for( int i = 0; i < num_faces; i++ )
01057     {
01058         int num_verts = *farray++;
01059         verts.resize( num_verts );
01060 
01061         // Fill in connectivity by looking up by gid in vert tuple_list
01062         for( int j = 0; j < num_verts; j++ )
01063         {
01064 #ifdef TUPLE_LIST
01065             int tindex = vert_map.find( 1, farray[j] );
01066             if( -1 == tindex )
01067             {
01068                 tmp_rval = MB_FAILURE;
01069                 break;
01070             }
01071             verts[j] = vert_map.get_ulong( tindex, 0 );
01072 #else
01073             verts[j] = ( vert_map[farray[j]] )[0];
01074 #endif
01075         }
01076         farray += num_verts;
01077 
01078         if( MB_SUCCESS == tmp_rval )
01079         {
01080             // Make face
01081             EntityType ftype = ( 3 == num_verts ? MBTRI : ( 4 == num_verts ? MBQUAD : MBPOLYGON ) );
01082             EntityHandle faceh;
01083             tmp_rval = mbImpl->create_element( ftype, &verts[0], num_verts, faceh );
01084             if( faceh ) new_faces.push_back( faceh );
01085         }
01086 
01087         if( MB_SUCCESS != tmp_rval ) rval = tmp_rval;
01088     }
01089 
01090     return rval;
01091 }
01092 
01093 ErrorCode ReadCCMIO::read_vertices( CCMIOSize_t /* proc */,
01094                                     CCMIOID /* processorID */,
01095                                     CCMIOID verticesID,
01096                                     CCMIOID /* topologyID */,
01097                                     Range* verts,
01098                                     TupleList& vert_map )
01099 {
01100     CCMIOError error = kCCMIONoErr;
01101 
01102     // Pre-read the number of vertices, so we can pre-allocate & read directly in
01103     CCMIOSize_t nverts = CCMIOSIZEC( 0 );
01104     CCMIOEntitySize( &error, verticesID, &nverts, NULL );
01105     CHK_SET_CCMERR( error, "Couldn't get number of vertices" );
01106 
01107     // Get # dimensions
01108     CCMIOSize_t dims;
01109     float scale;
01110     CCMIOReadVerticesf( &error, verticesID, &dims, NULL, NULL, NULL, CCMIOINDEXC( 0 ), CCMIOINDEXC( 1 ) );
01111     CHK_SET_CCMERR( error, "Couldn't get number of dimensions" );
01112 
01113     // Allocate vertex space
01114     EntityHandle node_handle = 0;
01115     std::vector< double* > arrays;
01116     readMeshIface->get_node_coords( 3, GETINT32( nverts ), MB_START_ID, node_handle, arrays );
01117 
01118     // Read vertex coords
01119     CCMIOID mapID;
01120     std::vector< double > tmp_coords( GETINT32( dims ) * GETINT32( nverts ) );
01121     CCMIOReadVerticesd( &error, verticesID, &dims, &scale, &mapID, &tmp_coords[0], CCMIOINDEXC( 0 ),
01122                         CCMIOINDEXC( 0 + nverts ) );
01123     CHK_SET_CCMERR( error, "Trouble reading vertex coordinates" );
01124 
01125     // Copy interleaved coords into moab blocked coordinate space
01126     int i = 0, threei = 0;
01127     for( ; i < nverts; i++ )
01128     {
01129         arrays[0][i] = tmp_coords[threei++];
01130         arrays[1][i] = tmp_coords[threei++];
01131         if( 3 == GETINT32( dims ) )
01132             arrays[2][i] = tmp_coords[threei++];
01133         else
01134             arrays[2][i] = 0.0;
01135     }
01136 
01137     // Scale, if necessary
01138     if( 1.0 != scale )
01139     {
01140         for( i = 0; i < nverts; i++ )
01141         {
01142             arrays[0][i] *= scale;
01143             arrays[1][i] *= scale;
01144             if( 3 == GETINT32( dims ) ) arrays[2][i] *= scale;
01145         }
01146     }
01147 
01148     // Read gids for vertices
01149     std::vector< int > gids( GETINT32( nverts ) );
01150     CCMIOReadMap( &error, mapID, &gids[0], CCMIOINDEXC( kCCMIOStart ), CCMIOINDEXC( kCCMIOEnd ) );
01151     CHK_SET_CCMERR( error, "Trouble reading vertex global ids" );
01152 
01153     // Put new vertex handles into range, and set gids for them
01154     Range new_verts( node_handle, node_handle + nverts - 1 );
01155     ErrorCode rval = mbImpl->tag_set_data( mGlobalIdTag, new_verts, &gids[0] );MB_CHK_SET_ERR( rval, "Couldn't set gids on vertices" );
01156 
01157     // Pack vert_map with global ids and handles for these vertices
01158 #ifdef TUPLE_LIST
01159     vert_map.resize( GETINT32( nverts ) );
01160     for( i = 0; i < GETINT32( nverts ); i++ )
01161     {
01162         vert_map.push_back( NULL, &gids[i], &node_handle, NULL );
01163 #else
01164     for( i = 0; i < GETINT32( nverts ); i++ )
01165     {
01166         ( vert_map[gids[i]] ).push_back( node_handle );
01167 #endif
01168         node_handle += 1;
01169     }
01170 
01171     if( verts ) verts->merge( new_verts );
01172 
01173     return MB_SUCCESS;
01174 }
01175 
01176 ErrorCode ReadCCMIO::get_processors( CCMIOID stateID,
01177                                      CCMIOID& processorID,
01178                                      CCMIOID& verticesID,
01179                                      CCMIOID& topologyID,
01180                                      CCMIOID& solutionID,
01181                                      std::vector< CCMIOSize_t >& procs,
01182                                      bool& /* has_solution */ )
01183 {
01184     CCMIOSize_t proc = CCMIOSIZEC( 0 );
01185     CCMIOError error = kCCMIONoErr;
01186 
01187     CCMIONextEntity( &error, stateID, kCCMIOProcessor, &proc, &processorID );
01188     CHK_SET_CCMERR( error, "CCMIONextEntity() failed" );
01189     if( CCMIOReadProcessor( NULL, processorID, &verticesID, &topologyID, NULL, &solutionID ) != kCCMIONoErr )
01190     {
01191         // Maybe no solution; try again
01192         CCMIOReadProcessor( &error, processorID, &verticesID, &topologyID, NULL, NULL );
01193         CHK_SET_CCMERR( error, "CCMIOReadProcessor() failed" );
01194         hasSolution = false;
01195     }
01196 
01197     procs.push_back( proc );
01198 
01199     return MB_SUCCESS;
01200 }
01201 
01202 ErrorCode ReadCCMIO::read_tag_values( const char* /* file_name */,
01203                                       const char* /* tag_name */,
01204                                       const FileOptions& /* opts */,
01205                                       std::vector< int >& /* tag_values_out */,
01206                                       const SubsetList* /* subset_list */ )
01207 {
01208     return MB_FAILURE;
01209 }
01210 
01211 }  // namespace moab
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines