ReadSms.cpp

Go to the documentation of this file.

/**
 * MOAB, a Mesh-Oriented datABase, is a software component for creating,
 * storing and accessing finite element mesh data.
 *
 * Copyright 2004 Sandia Corporation.  Under the terms of Contract
 * DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government
 * retains certain rights in this software.
 *
 * This library is free software; you can redistribute it and/or
 * modify it under the terms of the GNU Lesser General Public
 * License as published by the Free Software Foundation; either
 * version 2.1 of the License, or (at your option) any later version.
 *
 */

/**
 * \class ReadSms
 * \brief Sms (http://www.geuz.org/sms) file reader
 * \author Jason Kraftcheck
 */

#include "ReadSms.hpp"
#include "FileTokenizer.hpp"  // For file tokenizer
#include "Internals.hpp"
#include "moab/Interface.hpp"
#include "moab/ReadUtilIface.hpp"
#include "moab/Range.hpp"
#include "MBTagConventions.hpp"
#include "MBParallelConventions.h"
#include "moab/CN.hpp"

#include <cerrno>
#include <cstring>
#include <map>
#include <set>
#include <iostream>

#define CHECK( a )                       \
    if( MB_SUCCESS != result )           \
    {                                    \
        std::cerr << ( a ) << std::endl; \
        return result;                   \
    }

#define CHECKN( a ) \
    if( n != ( a ) ) return MB_FILE_WRITE_ERROR

namespace moab
{

ReaderIface* ReadSms::factory( Interface* iface )
{
    return new ReadSms( iface );
}

ReadSms::ReadSms( Interface* impl ) : mdbImpl( impl ), globalId( 0 ), paramCoords( 0 ), geomDimension( 0 ), setId( 0 )
{
    mdbImpl->query_interface( readMeshIface );
}

ReadSms::~ReadSms()
{
    if( readMeshIface )
    {
        mdbImpl->release_interface( readMeshIface );
        readMeshIface = 0;
    }
}

ErrorCode ReadSms::read_tag_values( const char* /* file_name */, const char* /* tag_name */,
                                    const FileOptions& /* opts */, std::vector< int >& /* tag_values_out */,
                                    const SubsetList* /* subset_list */ )
{
    return MB_NOT_IMPLEMENTED;
}

ErrorCode ReadSms::load_file( const char* filename, const EntityHandle* /* file_set */, const FileOptions& /* opts */,
                              const ReaderIface::SubsetList* subset_list, const Tag* file_id_tag )
{
    if( subset_list ) { MB_SET_ERR( MB_UNSUPPORTED_OPERATION, "Reading subset of files not supported for Sms" ); }

    setId = 1;

    // Open file
    FILE* file_ptr = fopen( filename, "r" );
    if( !file_ptr ) { MB_SET_ERR( MB_FILE_DOES_NOT_EXIST, filename << ": " << strerror( errno ) ); }

    const ErrorCode result = load_file_impl( file_ptr, file_id_tag );
    fclose( file_ptr );

    return result;
}

ErrorCode ReadSms::load_file_impl( FILE* file_ptr, const Tag* file_id_tag )
{
    bool warned         = false;
    double dum_params[] = { 0.0, 0.0, 0.0 };

    globalId = mdbImpl->globalId_tag();

    ErrorCode result =
        mdbImpl->tag_get_handle( "PARAMETER_COORDS", 3, MB_TYPE_DOUBLE, paramCoords, MB_TAG_DENSE | MB_TAG_CREAT );
    CHECK( "Failed to create param coords tag." );

    int negone = -1;
    result     = mdbImpl->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, geomDimension,
                                      MB_TAG_SPARSE | MB_TAG_CREAT, &negone );
    CHECK( "Failed to create geom dim tag." );

    int n;
    char line[256], all_line[1024];
    int file_type;

    if( fgets( all_line, sizeof( all_line ), file_ptr ) == NULL ) { return MB_FAILURE; }
    if( sscanf( all_line, "%s %d", line, &file_type ) != 2 ) { return MB_FAILURE; }

    if( 3 == file_type )
    {
        result = read_parallel_info( file_ptr );
        CHECK( "Failed to read parallel info." );
    }

    int nregions, nfaces, nedges, nvertices, npoints;
    n = fscanf( file_ptr, "%d %d %d %d %d", &nregions, &nfaces, &nedges, &nvertices, &npoints );
    CHECKN( 5 );
    if( nregions < 0 || nfaces < 0 || nedges < 0 || nvertices < 0 || npoints < 0 ) return MB_FILE_WRITE_ERROR;

    // Create the vertices
    std::vector< double* > coord_arrays;
    EntityHandle vstart = 0;
    result              = readMeshIface->get_node_coords( 3, nvertices, MB_START_ID, vstart, coord_arrays );
    CHECK( "Failed to get node arrays." );

    if( file_id_tag )
    {
        result = add_entities( vstart, nvertices, file_id_tag );MB_CHK_ERR( result );
    }

    EntityHandle this_gent, new_handle;
    std::vector< EntityHandle > gentities[4];
    int gent_id, dum_int;
    int gent_type, num_connections;

    for( int i = 0; i < nvertices; i++ )
    {
        n = fscanf( file_ptr, "%d", &gent_id );
        CHECKN( 1 );
        if( !gent_id ) continue;

        n = fscanf( file_ptr, "%d %d %lf %lf %lf", &gent_type, &num_connections, coord_arrays[0] + i,
                    coord_arrays[1] + i, coord_arrays[2] + i );
        CHECKN( 5 );

        result = get_set( gentities, gent_type, gent_id, geomDimension, this_gent, file_id_tag );MB_CHK_ERR( result );

        new_handle = vstart + i;
        result     = mdbImpl->add_entities( this_gent, &new_handle, 1 );
        CHECK( "Adding vertex to geom set failed." );

        switch( gent_type )
        {
            case 1:
                n = fscanf( file_ptr, "%le", dum_params );
                CHECKN( 1 );
                result = mdbImpl->tag_set_data( paramCoords, &new_handle, 1, dum_params );
                CHECK( "Failed to set param coords tag for vertex." );
                break;
            case 2:
                n = fscanf( file_ptr, "%le %le %d", dum_params, dum_params + 1, &dum_int );
                CHECKN( 3 );
                dum_params[2] = dum_int;
                result        = mdbImpl->tag_set_data( paramCoords, &new_handle, 1, dum_params );
                CHECK( "Failed to set param coords tag for vertex." );
                break;
            default:
                break;
        }
    }  // End of reading vertices

    // *******************************
    //  Read Edges
    // *******************************

    int vert1, vert2, num_pts;
    std::vector< EntityHandle > everts( 2 );
    EntityHandle estart, *connect;
    result = readMeshIface->get_element_connect( nedges, 2, MBEDGE, 1, estart, connect );
    CHECK( "Failed to create array of edges." );

    if( file_id_tag )
    {
        result = add_entities( estart, nedges, file_id_tag );
        if( MB_SUCCESS != result ) return result;
    }

    for( int i = 0; i < nedges; i++ )
    {
        n = fscanf( file_ptr, "%d", &gent_id );
        CHECKN( 1 );
        if( !gent_id ) continue;

        n = fscanf( file_ptr, "%d %d %d %d %d", &gent_type, &vert1, &vert2, &num_connections, &num_pts );
        CHECKN( 5 );
        if( vert1 < 1 || vert1 > nvertices ) return MB_FILE_WRITE_ERROR;
        if( vert2 < 1 || vert2 > nvertices ) return MB_FILE_WRITE_ERROR;

        connect[0] = vstart + vert1 - 1;
        connect[1] = vstart + vert2 - 1;
        if( num_pts > 1 && !warned )
        {
            std::cout << "Warning: num_points > 1 not supported; choosing last one." << std::endl;
            warned = true;
        }

        result = get_set( gentities, gent_type, gent_id, geomDimension, this_gent, file_id_tag );
        CHECK( "Problem getting geom set for edge." );

        new_handle = estart + i;
        result     = mdbImpl->add_entities( this_gent, &new_handle, 1 );
        CHECK( "Failed to add edge to geom set." );

        connect += 2;

        for( int j = 0; j < num_pts; j++ )
        {
            switch( gent_type )
            {
                case 1:
                    n = fscanf( file_ptr, "%le", dum_params );
                    CHECKN( 1 );
                    result = mdbImpl->tag_set_data( paramCoords, &new_handle, 1, dum_params );
                    CHECK( "Failed to set param coords tag for edge." );
                    break;
                case 2:
                    n = fscanf( file_ptr, "%le %le %d", dum_params, dum_params + 1, &dum_int );
                    CHECKN( 3 );
                    dum_params[2] = dum_int;
                    result        = mdbImpl->tag_set_data( paramCoords, &new_handle, 1, dum_params );
                    CHECK( "Failed to set param coords tag for edge." );
                    break;
                default:
                    break;
            }
        }
    }  // End of reading edges

    // *******************************
    //  Read Faces
    // *******************************
    std::vector< EntityHandle > bound_ents, bound_verts, new_faces;
    int bound_id;
    Range shverts;
    new_faces.resize( nfaces );
    int num_bounding;

    for( int i = 0; i < nfaces; i++ )
    {
        n = fscanf( file_ptr, "%d", &gent_id );
        CHECKN( 1 );
        if( !gent_id ) continue;

        n = fscanf( file_ptr, "%d %d", &gent_type, &num_bounding );
        CHECKN( 2 );

        result = get_set( gentities, gent_type, gent_id, geomDimension, this_gent, file_id_tag );
        CHECK( "Problem getting geom set for face." );

        bound_ents.resize( num_bounding + 1 );
        bound_verts.resize( num_bounding );
        for( int j = 0; j < num_bounding; j++ )
        {
            n = fscanf( file_ptr, "%d ", &bound_id );
            CHECKN( 1 );
            if( 0 > bound_id ) bound_id = abs( bound_id );
            assert( 0 < bound_id && bound_id <= nedges );
            if( bound_id < 1 || bound_id > nedges ) return MB_FILE_WRITE_ERROR;
            bound_ents[j] = estart + abs( bound_id ) - 1;
        }

        // Convert edge-based model to vertex-based one
        for( int j = 0; j < num_bounding; j++ )
        {
            if( j == num_bounding - 1 ) bound_ents[j + 1] = bound_ents[0];
            result = mdbImpl->get_adjacencies( &bound_ents[j], 2, 0, false, shverts );
            CHECK( "Failed to get vertices bounding edge." );
            assert( shverts.size() == 1 );
            bound_verts[j] = *shverts.begin();
            shverts.clear();
        }

        result = mdbImpl->create_element( ( EntityType )( MBTRI + num_bounding - 3 ), &bound_verts[0],
                                          bound_verts.size(), new_faces[i] );
        CHECK( "Failed to create edge." );

        result = mdbImpl->add_entities( this_gent, &new_faces[i], 1 );
        CHECK( "Failed to add edge to geom set." );

        int num_read = fscanf( file_ptr, "%d", &num_pts );
        if( !num_pts || !num_read ) continue;

        for( int j = 0; j < num_pts; j++ )
        {
            switch( gent_type )
            {
                case 1:
                    n = fscanf( file_ptr, "%le", dum_params );
                    CHECKN( 1 );
                    result = mdbImpl->tag_set_data( paramCoords, &new_faces[i], 1, dum_params );
                    CHECK( "Failed to set param coords tag for face." );
                    break;
                case 2:
                    n = fscanf( file_ptr, "%le %le %d", dum_params, dum_params + 1, &dum_int );
                    CHECKN( 3 );
                    dum_params[2] = dum_int;
                    result        = mdbImpl->tag_set_data( paramCoords, &new_faces[i], 1, dum_params );
                    CHECK( "Failed to set param coords tag for face." );
                    break;
                default:
                    break;
            }
        }
    }  // End of reading faces

    if( file_id_tag )
    {
        result = readMeshIface->assign_ids( *file_id_tag, &new_faces[0], new_faces.size(), 1 );
        if( MB_SUCCESS != result ) return result;
    }

    // *******************************
    //  Read Regions
    // *******************************
    int sense[MAX_SUB_ENTITIES];
    bound_verts.resize( MAX_SUB_ENTITIES );

    std::vector< EntityHandle > regions;
    if( file_id_tag ) regions.resize( nregions );
    for( int i = 0; i < nregions; i++ )
    {
        n = fscanf( file_ptr, "%d", &gent_id );
        CHECKN( 1 );
        if( !gent_id ) continue;
        result = get_set( gentities, 3, gent_id, geomDimension, this_gent, file_id_tag );
        CHECK( "Couldn't get geom set for region." );
        n = fscanf( file_ptr, "%d", &num_bounding );
        CHECKN( 1 );
        bound_ents.resize( num_bounding );
        for( int j = 0; j < num_bounding; j++ )
        {
            n = fscanf( file_ptr, "%d ", &bound_id );
            CHECKN( 1 );
            assert( abs( bound_id ) < (int)new_faces.size() + 1 && bound_id );
            if( !bound_id || abs( bound_id ) > nfaces ) return MB_FILE_WRITE_ERROR;
            sense[j]      = ( bound_id < 0 ) ? -1 : 1;
            bound_ents[j] = new_faces[abs( bound_id ) - 1];
        }

        EntityType etype;
        result = readMeshIface->get_ordered_vertices( &bound_ents[0], sense, num_bounding, 3, &bound_verts[0], etype );
        CHECK( "Failed in get_ordered_vertices." );

        // Make the element
        result = mdbImpl->create_element( etype, &bound_verts[0], CN::VerticesPerEntity( etype ), new_handle );
        CHECK( "Failed to create region." );

        result = mdbImpl->add_entities( this_gent, &new_handle, 1 );
        CHECK( "Failed to add region to geom set." );

        if( file_id_tag ) regions[i] = new_handle;

        n = fscanf( file_ptr, "%d ", &dum_int );
        CHECKN( 1 );
    }  // End of reading regions

    if( file_id_tag )
    {
        result = readMeshIface->assign_ids( *file_id_tag, &regions[0], regions.size(), 1 );
        if( MB_SUCCESS != result ) return result;
    }

    return MB_SUCCESS;
}

ErrorCode ReadSms::get_set( std::vector< EntityHandle >* sets, int set_dim, int set_id, Tag dim_tag,
                            EntityHandle& this_set, const Tag* file_id_tag )
{
    ErrorCode result = MB_SUCCESS;

    if( set_dim < 0 || set_dim > 3 ) return MB_FILE_WRITE_ERROR;

    if( (int)sets[set_dim].size() <= set_id || !sets[set_dim][set_id] )
    {
        if( (int)sets[set_dim].size() <= set_id ) sets[set_dim].resize( set_id + 1, 0 );

        if( !sets[set_dim][set_id] )
        {
            result = mdbImpl->create_meshset( MESHSET_SET, sets[set_dim][set_id] );
            if( MB_SUCCESS != result ) return result;
            result = mdbImpl->tag_set_data( globalId, &sets[set_dim][set_id], 1, &set_id );
            if( MB_SUCCESS != result ) return result;
            result = mdbImpl->tag_set_data( dim_tag, &sets[set_dim][set_id], 1, &set_dim );
            if( MB_SUCCESS != result ) return result;

            if( file_id_tag )
            {
                result = mdbImpl->tag_set_data( *file_id_tag, &sets[set_dim][set_id], 1, &setId );
                ++setId;
            }
        }
    }

    this_set = sets[set_dim][set_id];

    return result;
}

ErrorCode ReadSms::read_parallel_info( FILE* file_ptr )
{
    // ErrorCode result;

    // Read partition info
    int nparts, part_id, num_ifaces, num_corner_ents;
    int num_read = fscanf( file_ptr, "%d %d %d %d", &nparts, &part_id, &num_ifaces, &num_corner_ents );
    if( !num_read ) return MB_FAILURE;

    // Read interfaces
    int iface_id, iface_dim, iface_own, num_iface_corners;
    // EntityHandle this_iface;
    std::vector< int >* iface_corners = NULL;
    for( int i = 0; i < num_ifaces; i++ )
    {
        num_read = fscanf( file_ptr, "%d %d %d %d", &iface_id, &iface_dim, &iface_own, &num_iface_corners );
        if( !num_read ) return MB_FAILURE;

        // result = get_set(sets, iface_dim, iface_id, dim_tag, iface_own, this_iface);
        // CHECK("Failed to make iface set.");

        // Read the corner ids and store them on the set for now
        iface_corners = new std::vector< int >( num_iface_corners );
        for( int j = 0; j < num_iface_corners; j++ )
        {
            num_read = fscanf( file_ptr, "%d", &( *iface_corners )[j] );
            if( !num_read )
            {
                delete iface_corners;
                return MB_FAILURE;
            }
        }

        // result = tag_set_data(ifaceCornerTag, &this_iface, 1,
        //&iface_corners);
        // CHECK("Failed to set iface corner tag.");

        delete iface_corners;
        iface_corners = NULL;
    }

    // Interface data has been read
    return MB_SUCCESS;
}

ErrorCode ReadSms::add_entities( EntityHandle start, EntityHandle count, const Tag* file_id_tag )
{
    if( !count || !file_id_tag ) return MB_FAILURE;

    Range range;
    range.insert( start, start + count - 1 );
    return readMeshIface->assign_ids( *file_id_tag, range, 1 );
}

}  // namespace moab