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