1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
/**
 * 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 );<--- %lu in format string (no. 1) requires 'unsigned long' but the argument type is 'signed long'.
    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 );<--- %lu in format string (no. 1) requires 'unsigned long' but the argument type is 'signed long'.
    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;
}