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
/**
 * \class ReadTemplate
 * \brief Template for writing a new reader in MOAB
 *
 */

#include "ReadTemplate.hpp"
#include "moab/Interface.hpp"
#include "moab/ReadUtilIface.hpp"
#include "moab/Range.hpp"
#include "moab/FileOptions.hpp"

#include <cstdio>
#include <cassert>

namespace moab
{

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

ReadTemplate::ReadTemplate( Interface* impl ) : mbImpl( impl ), fileName( NULL )
{
    mbImpl->query_interface( readMeshIface );
}

ReadTemplate::~ReadTemplate()
{
    if( readMeshIface )
    {
        mbImpl->release_interface( readMeshIface );
        readMeshIface = 0;
    }
}

ErrorCode ReadTemplate::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 ReadTemplate::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 )
    {
        // See src/moab/ReaderIface.hpp, definition of SubsetList struct; this basically specifies
        // an integer tag and tag values for sets to read on this proc, or a part number and total #
        // parts for reading a trivial partition of entities
    }

    // Save filename to member variable so we don't need to pass as an argument
    // to called functions
    fileName = filename;

    // Process options; see src/FileOptions.hpp for API for FileOptions class, and
    // doc/metadata_info.doc for a description of various options used by some of the readers in
    // MOAB
    ErrorCode result = process_options( opts );MB_CHK_SET_ERR( result, fileName << ": problem reading options" );

    // Open file; filePtr is member of ReadTemplate, change to whatever mechanism is used to
    // identify file
    FILE* filePtr = fopen( fileName, "r" );
    if( !filePtr )
    {
        MB_SET_ERR( MB_FILE_DOES_NOT_EXIST, fileName << ": fopen returned error" );
    }

    // Read number of verts, elements, sets
    long num_verts = 0, num_elems = 0, num_sets = 0;

    // read_ents keeps a running set of entities read from this file, including vertices, elements,
    // and sets; these will get added to file_set (if input) at the end of the read
    Range read_ents;

    // start_vertex is passed back so we know how to convert indices from the file into vertex
    // handles; most of the time this is done by adding start_vertex to the (0-based) index; if the
    // index is 1-based, you also need to subtract one; see read_elements for details
    EntityHandle start_vertex;
    result = read_vertices( num_verts, start_vertex, read_ents );
    if( MB_SUCCESS != result )
    {
        fclose( filePtr );
        return result;
    }

    // Create/read elements; this template assumes that all elements are the same type, so can be
    // read in a single call to read_elements, and kept track of with a single start_elem handle. If
    // there are more entity types, might have to keep these start handles in an array/vector.
    // start_elem is only really needed if you're reading sets later, and need to convert some
    // file-based index to an entity handle
    EntityHandle start_elem;
    result = read_elements( num_elems, start_vertex, start_elem, read_ents );
    if( MB_SUCCESS != result )
    {
        fclose( filePtr );
        return result;
    }

    // Read/create entity sets; typically these sets have some tag identifying what they're for, see
    // doc/metadata_info.doc for examples of different kinds of sets and how they're marked
    result = create_sets( num_sets, start_vertex, num_verts, start_elem, num_elems, read_ents );
    if( MB_SUCCESS != result )
    {
        fclose( filePtr );
        return result;
    }

    // Finally, add all read_ents into the file set, if one was input
    if( file_set && *file_set )
    {
        result = mbImpl->add_entities( *file_set, read_ents );
        if( MB_SUCCESS != result )
        {
            fclose( filePtr );
            return result;
        }
    }

    fclose( filePtr );

    return result;
}

ErrorCode ReadTemplate::read_vertices( int num_verts, EntityHandle& start_vertex, Range& read_ents )
{
    // Allocate nodes; these are allocated in one shot, get contiguous handles starting with
    // start_handle, and the reader is passed back double*'s pointing to MOAB's native storage for
    // vertex coordinates for those verts
    std::vector< double* > coord_arrays;
    ErrorCode result = readMeshIface->get_node_coords( 3, num_verts, 1, start_vertex, coord_arrays );MB_CHK_SET_ERR( result, fileName << ": Trouble reading vertices" );

    // Fill in vertex coordinate arrays
    double *x = coord_arrays[0], *y = coord_arrays[1], *z = coord_arrays[2];
    for( long i = 0; i < num_verts; ++i )
    {
        // Read x/y/z; do something with them for now to avoid warning
        if( x || y || z )
        {
        }
    }

    if( num_verts ) read_ents.insert( start_vertex, start_vertex + num_verts - 1 );

    return result;
}

//! Read/create elements
ErrorCode ReadTemplate::read_elements( int num_elems,
                                       EntityHandle start_vertex,
                                       EntityHandle& start_elem,
                                       Range& read_ents )
{
    // Get the entity type being read
    EntityType ent_type = MBHEX;

    // Get the number of vertices per entity
    int verts_per_elem = 8;

    // Create the element sequence; passes back a pointer to the internal storage for connectivity
    // and the starting entity handle
    EntityHandle* conn_array;
    ErrorCode result =
        readMeshIface->get_element_connect( num_elems, verts_per_elem, ent_type, 1, start_elem, conn_array );MB_CHK_SET_ERR( result, fileName << ": Trouble reading elements" );

    // Read connectivity into conn_array directly
    for( long i = 0; i < num_elems; i++ )
    {
        // Read connectivity
    }

    // Convert file-based connectivity indices to vertex handles in-place; be careful, if indices
    // are smaller than handles, need to do from the end of the list so we don't overwrite data
    //
    // Here, we assume indices are smaller than handles, just for demonstration; create an
    // integer-type pointer to connectivity initialized to same start of connectivity array
    int* ind_array = reinterpret_cast< int* >( conn_array );
    // OFFSET is value of first vertex index in file; most files are 1-based, but some might be
    // 0-based
    int OFFSET = 1;
    for( long i = num_elems * verts_per_elem - 1; i >= 0; i-- )
    {
        conn_array[i] = ind_array[i] + start_vertex + OFFSET;

        // This assert assumes last handle in read_ents is highest vertex handle in this file
        assert( conn_array[i] >= start_vertex && conn_array[i] <= *read_ents.rbegin() );
    }

    // Notify MOAB of the new elements
    result = readMeshIface->update_adjacencies( start_elem, num_elems, verts_per_elem, conn_array );
    if( MB_SUCCESS != result ) return result;

    // Add elements to read_ents
    if( num_elems ) read_ents.insert( start_elem, start_elem + num_elems - 1 );

    return MB_SUCCESS;
}

//! Read/create sets
ErrorCode ReadTemplate::create_sets( int num_sets,
                                     EntityHandle /*start_vertex*/,
                                     int /*num_verts*/,
                                     EntityHandle /*start_elem*/,
                                     int /*num_elems*/,
                                     Range& read_ents )
{
    ErrorCode result = MB_SUCCESS;<--- Variable 'result' is assigned a value that is never used.
    EntityHandle this_set;

    for( int i = 0; i < num_sets; i++ )
    {
        // Create set
        result = mbImpl->create_meshset( MESHSET_SET, this_set );MB_CHK_SET_ERR( result, fileName << ": Trouble creating set" );

        Range set_ents;
        // Read/compute what's in this set; REMEMBER TO CONVERT THESE TO MOAB HANDLES

        // Add them to the set
        result = mbImpl->add_entities( this_set, set_ents );MB_CHK_SET_ERR( result, fileName << ": Trouble putting entities in set" );

        // Add the new set to read_ents
        read_ents.insert( this_set );
    }

    return MB_SUCCESS;
}

ErrorCode ReadTemplate::process_options( const FileOptions& opts )
{
    // Mark all options seen, to avoid compile warning on unused variable
    opts.mark_all_seen();
    return MB_SUCCESS;
}

}  // namespace moab