hexes_to_gmsh.c

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 2006 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.
 *
 */
 
/* Example file for using mhdf interface for reading MOAB native file foramt
 *
 * Requires libmhdf from MOAB source and HDF5 library
 *
 * Reads:
 *  - all node coordinates
 *  - a single group of hexahedral elements
 *  - GLOBAL_ID tag if present
 * and writes as Gmsh file version 1.0.
 *
 * Does not contain examples for reading meshsets or
 * reading polygon or polyhedron elements.
 *
 * A simple exorcise:
 *    Clean up handling of variable number of dimensions by reading
 *    node coordinate data one dimension at a time using
 *    mhdf_readNodeCoord
 *
 * A more complex exorcise:
 *    Read meshsets with MATERIAL_SET tag, which MOAB uses to
 *    represent element blocks, and assign appropriate element
 *    block for each hex in Gmsh file.
 *    Hint: look for and read MATERIAL_SET tag first to get list of entities
 *          then read set contents for sets with that tag
 *          then match hex file ids to lists of set contents
 *    Be careful of the mhdf_SET_RANGE_BIT for each set when reading set
 *    contents.
 */
 
#include <stdlib.h>
#include <stdio.h>
#include <assert.h>
#include <string.h>
 
/* mhdf API */
#include "mhdf.h"
/* need constants for native types (e.g. H5T_NATIVE_UINT)*/
#include <H5Tpublic.h>
 
/* Macro to check return status from mhdf calls.  Exit on error. */
#define CHK_ERR( A )                                             \
    do                                                           \
    {                                                            \
        if( mhdf_isError( A ) )                                  \
        {                                                        \
            fprintf( stderr, "Error: %s\n", mhdf_message( A ) ); \
            exit( 2 );                                           \
        }                                                        \
    } while( 0 )
 
int main( int argc, char* argv[] )
{
    /* input file */
    const char* filename;
    mhdf_FileHandle file;
    mhdf_Status status;
    mhdf_Status* const sptr = &status;
    hid_t handle; /* generic handle used to refer to any data block in file */
 
    /* output file */
    const char* gmsh_filename;
    FILE* gmsh;
    unsigned gmsh_type;            /* hexahedral element type number */
    double x, y, z;                /* temp storage of node coordinates */
    unsigned node_offset, node_id; /* temporary values */
    unsigned* connectivity;        /* temporary value */
 
    /* node data */
    long numnode;       /* total number of nodes */
    long nodestart;     /* file id of first node in list */
    int dimension;      /* coordinate values per node */
    double* nodecoords; /* interleaved node coordinates */
    unsigned* nodeids;  /* GLOBAL_ID value for nodes */
    int have_nodeids = 0;
 
    /* hex data */
    char* hexgroup = NULL;     /* name of element group containing hexes */
    long numhex;               /* total number of hexahedral elements */
    long hexstart;             /* file id of first hex in group */
    int nodes_per_hex;         /* length of connectivity list for a hex */
    unsigned* hexconnectivity; /* hex connectivity data */
    unsigned* hexids;          /* GLOBAL_ID value for hexes */
    int have_hexids = 0;
 
    /* list of element groups in file */
    char** elem_groups;
    unsigned num_elem_groups;
    char namebuffer[64];
 
    /* tag data for accessing GLOBAL_ID */
    int tagsize;                   /* number of values for each entity */
    int ts, td, tg;                /* unused tag properties */
    int havesparse, havedense;     /* Boolean values */
    enum mhdf_TagDataType tagtype; /* base data type of tag */
    hid_t sparse_handle[2];        /* handle pair for sparse tag data */
    unsigned* sparse_entities;     /* temp storage of sparse tag file ids */
    unsigned* sparse_ids;          /* temp storage of GLOBAL_ID values in spasre tag */
    long junk, numtag;             /* number of entities for which tag data is available */
    long fileid, globalid;         /* temporary values */
    long ncount = 0, hcount = 0;   /* temporary count of number of tag values */
 
    /* iteration */
    long i;
    int j;
    unsigned k;
 
    /* process CL args (expect input .h5m file and output .gmsh file name) */
    if( argc != 3 )
    {
        fprintf( stderr, "Usage: %s <input_file> <output_file>\n", argv[0] );
        return 1;
    }
    filename      = argv[1];
    gmsh_filename = argv[2];
 
    /* Open the file */
    file = mhdf_openFile( filename, 0, 0, sptr );CHK_ERR( sptr );
 
    /* Read node coordinates. */
    handle = mhdf_openNodeCoords( file, &numnode, &dimension, &nodestart, sptr );CHK_ERR( sptr );
    nodecoords = (double*)malloc( dimension * numnode * sizeof( double ) );
    mhdf_readNodeCoords( handle, 0, numnode, nodecoords, sptr );CHK_ERR( sptr );
    mhdf_closeData( file, handle, sptr );CHK_ERR( sptr );
 
    /* Find first element group containing hexahedra */
    elem_groups = mhdf_getElemHandles( file, &num_elem_groups, sptr );CHK_ERR( sptr );
    for( k = 0; k < num_elem_groups; ++k )
    {
        mhdf_getElemTypeName( file, elem_groups[k], namebuffer, sizeof( namebuffer ), sptr );CHK_ERR( sptr );
        if( !hexgroup && !strcmp( mdhf_HEX_TYPE_NAME, namebuffer ) )
            hexgroup = strdup( elem_groups[k] );
        else
            printf( "Skipping element group '%s' containing element of type '%s'\n", elem_groups[k], namebuffer );
    }
    free( elem_groups );
 
    if( !hexgroup )
    {
        fprintf( stderr, "No Hexahedra defined in file\n" );
        return 4;
    }
 
    /* Read Hexahedron connectivity */
    handle = mhdf_openConnectivity( file, hexgroup, &nodes_per_hex, &numhex, &hexstart, sptr );CHK_ERR( sptr );
    hexconnectivity = (unsigned*)malloc( numhex * nodes_per_hex * sizeof( unsigned ) );
    mhdf_readConnectivity( handle, 0, numhex, H5T_NATIVE_UINT, hexconnectivity, sptr );CHK_ERR( sptr );
    mhdf_closeData( file, handle, sptr );CHK_ERR( sptr );
    /* Note: hex connectivity list contains file-space node IDs, which are
             the nodes in the sequence they are read from the file, with
             the first node having an ID of 'nodestart' */
 
    /* Check for "GLOBAL_ID" tag */
    nodeids = (unsigned*)malloc( numnode * sizeof( unsigned ) );
    hexids  = (unsigned*)malloc( numhex * sizeof( unsigned ) );
    mhdf_getTagInfo( file, "GLOBAL_ID", &tagtype, &tagsize, &ts, &td, &tg, &havesparse, sptr );
 
    /* If have GLOBAL_ID tag, try to read values for nodes and hexes */
    if( !mhdf_isError( sptr ) )
    {
        /* Check that the tag contains what we expect */
        if( tagtype != mhdf_INTEGER || tagsize != 1 )
        {
            fprintf( stderr, "ERROR: Invalid data type for 'GLOBAL_ID' tag.\n" );
            exit( 3 );
        }
 
        /* Check for and read dense-format tag data for nodes */
        havedense = mhdf_haveDenseTag( file, "GLOBAL_ID", mhdf_node_type_handle(), sptr );CHK_ERR( sptr );
        if( havedense )
        {
            handle = mhdf_openDenseTagData( file, "GLOBAL_ID", mhdf_node_type_handle(), &numtag, sptr );CHK_ERR( sptr );
            assert( numtag == numnode );
            mhdf_readDenseTag( handle, 0, numtag, H5T_NATIVE_UINT, nodeids, sptr );CHK_ERR( sptr );
            mhdf_closeData( file, handle, sptr );CHK_ERR( sptr );
            have_nodeids = 1;
        }
        /* Check for and read dense-format tag data for hexes */
        havedense = mhdf_haveDenseTag( file, "GLOBAL_ID", hexgroup, sptr );CHK_ERR( sptr );
        if( havedense )
        {
            handle = mhdf_openDenseTagData( file, "GLOBAL_ID", hexgroup, &numtag, sptr );CHK_ERR( sptr );
            assert( numtag == numhex );
            mhdf_readDenseTag( handle, 0, numtag, H5T_NATIVE_UINT, hexids, sptr );CHK_ERR( sptr );
            mhdf_closeData( file, handle, sptr );CHK_ERR( sptr );
            have_hexids = 1;
        }
        /* Check for and read sparse-format tag data */
        if( havesparse )
        {
            mhdf_openSparseTagData( file, "GLOBAL_ID", &numtag, &junk, sparse_handle, sptr );CHK_ERR( sptr );
            sparse_entities = (unsigned*)malloc( numtag * sizeof( unsigned ) );
            mhdf_readSparseTagEntities( sparse_handle[0], 0, numtag, H5T_NATIVE_UINT, sparse_entities, sptr );CHK_ERR( sptr );
            sparse_ids = (unsigned*)malloc( numtag * sizeof( unsigned ) );
            mhdf_readSparseTagValues( sparse_handle[1], 0, numtag, H5T_NATIVE_UINT, sparse_ids, sptr );CHK_ERR( sptr );
            mhdf_closeData( file, sparse_handle[0], sptr );CHK_ERR( sptr );
            mhdf_closeData( file, sparse_handle[1], sptr );CHK_ERR( sptr );
 
            /* Set hex and node ids from sparse tag data */
            for( i = 0; i < numtag; ++i )
            {
                fileid   = sparse_entities[i];
                globalid = sparse_ids[i];
                if( fileid >= nodestart && fileid - nodestart < numnode )
                {
                    nodeids[fileid - nodestart] = globalid;
                    ++ncount;
                }
                else if( fileid >= hexstart && fileid - hexstart < numhex )
                {
                    hexids[fileid - hexstart] = globalid;
                    ++hcount;
                }
            }
            free( sparse_ids );
            free( sparse_entities );
 
            /* make sure there was an ID for each node and each hex */
            if( ncount == numnode ) have_nodeids = 1;
            if( hcount == numhex ) have_hexids = 1;
 
        } /* end have sparse tag for GLOBAL_ID */
    }     /* end have GLOBAL_ID tag */
 
    /* done with input file */
    free( hexgroup );
    mhdf_closeFile( file, sptr );CHK_ERR( sptr );
 
    /* if no GLOBAL_ID, just use incrementing values */
    if( !have_nodeids )
        for( i = 0; i < numnode; ++i )
            nodeids[i] = i + 1;
    if( !have_hexids )
        for( i = 0; i < numhex; ++i )
            hexids[i] = i + 1;
 
    /* write out as gmesh file version 1.0 */
 
    /* get gmsh type for hexahedrons */
    if( nodes_per_hex == 8 )
        gmsh_type = 5;
    else if( nodes_per_hex == 27 )
        gmsh_type = 12;
    else
    {
        fprintf( stderr, "Cannot store %d node hex in gmsh file.\n", nodes_per_hex );
        exit( 4 );
    }
 
    /* open file */
    gmsh = fopen( gmsh_filename, "w" );
 
    /* Write node data.  If dimension is less than 3,
       write zero for other coordinate values.  In the
       (highly unlikely) case that dimension is greater
       than three, disregard higher-dimension coordinate
       values. */
    fprintf( gmsh, "$NOD\n" );
    fprintf( gmsh, "%lu\n", numnode );
    for( i = 0; i < numnode; ++i )
    {
        x = nodecoords[dimension * i];
        y = z = 0.0;
        if( dimension > 1 )
        {
            y = nodecoords[dimension * i + 1];
            if( dimension > 2 )
            {
                z = nodecoords[dimension * i + 2];
            }
        }
        fprintf( gmsh, "%u %f %f %f\n", nodeids[i], x, y, z );
    }
 
    /* Write element connectivity data */
    fprintf( gmsh, "$ENDNOD\n$ELM\n" );
    fprintf( gmsh, "%lu\n", numhex );
    for( i = 0; i < numhex; ++i )
    {
        fprintf( gmsh, "%u %u 1 1 %d", hexids[i], gmsh_type, nodes_per_hex );
        /* connectivity list for this hex */
        connectivity = hexconnectivity + i * nodes_per_hex;
        for( j = 0; j < nodes_per_hex; ++j )
        {
            /* get offset in node list from file id */
            node_offset = connectivity[j] - nodestart;
            /* get node id from ID list */
            node_id = nodeids[node_offset];
            fprintf( gmsh, " %u", node_id );
        }
        fprintf( gmsh, "\n" );
    }
    fprintf( gmsh, "$ENDELM\n" );
    fclose( gmsh );
    return 0;
}