DirectAccessWithHoles.cpp File Reference

Example demonstrating direct access to MOAB data with holes in entity handle space. More...

#include "moab/Core.hpp"
#include "moab/ProgOptions.hpp"
#include "moab/ReadUtilIface.hpp"
#include <map>
#include <iostream>
#include <cassert>

Include dependency graph for DirectAccessWithHoles.cpp:

Go to the source code of this file.

Classes
struct	tag_struct

Functions
ErrorCode	create_mesh_with_holes (Interface *mbImpl, int nquads, int nholes)
int	main (int argc, char **argv)

Detailed Description

Example demonstrating direct access to MOAB data with holes in entity handle space.

This example shows how to:

• Create a structured mesh with holes in the entity handle space
• Use direct access iterators to avoid API overhead
• Handle non-contiguous entity handles efficiently
• Access vertex coordinates, element connectivity, and tag data directly
• Work with adjacency lists for vertex-to-element relationships
• Use maps to track data chunks for non-contiguous handles

The program creates a 1D row of quad elements with specified holes, then demonstrates direct access to MOAB's internal data structures even when entity handles are not contiguous.

Author

MOAB Development Team

Date

2024

Use direct access to MOAB data to avoid calling through API
This example creates a 1d row of quad elements, with a user-specified number of "holes" (missing quads) in the row:

*  ----------------------      ----------------------      --------
*  |      |      |      |      |      |      |      |      |      |
*  |      |      |      |(hole)|      |      |      |(hole)|      | ...
*  |      |      |      |      |      |      |      |      |      |
*  ----------------------      ----------------------      --------
* 

This makes (nholes+1) contiguous runs of quad handles in the handle space This example shows how to use the xxx_iterate functions in MOAB (xxx = coords, connect, tag, adjacencies) to get direct pointer access to MOAB internal storage, which can be used without calling through the MOAB API.

1. Initialize MOAB
   
2. Create a quad mesh with holes, as depicted above
3. Create 2 dense tags (tag1, tag2) for avg position to assign to quads, and # verts per quad (tag3)
4. Get connectivity, coordinate, tag1 iterators
5. Iterate through quads, computing midpoint based on vertex positions, set on quad-based tag1
6. Set up map from starting quad handle for a chunk to struct of (tag1_ptr, tag2_ptr, tag3_ptr), pointers to the dense tag storage for those tags for the chunk
7. Iterate through vertices, summing positions into tag2 on connected quads and incrementing vertex count
8. Iterate through quads, normalizing tag2 by vertex count and comparing values of tag1 and tag2

To compile:
make DirectAccessWithHoles
To run: ./DirectAccess [-nquads <# quads>] [-holes <# holes>]

Parameters

argc	Number of command line arguments
argv	Command line arguments array

Returns

0 on success, 1 on failure

Definition in file DirectAccessWithHoles.cpp.

Function Documentation

create_mesh_with_holes()

ErrorCode create_mesh_with_holes	(	Interface *	mbImpl,
		int	nquads,
		int	nholes
	)

Definition at line 266 of file DirectAccessWithHoles.cpp.

{
    // First make the mesh, a 1d array of quads with left hand side x = elem_num; vertices are
    // numbered in layers
    ReadUtilIface* read_iface;
    MB_CHK_SET_ERR( mbImpl->query_interface( read_iface ), "Error in query_interface" );
    vector< double* > coords;
    EntityHandle start_vert, start_elem, *connect;
    // Create verts, num is 4(nquads+1) because they're in a 1d row; will initialize coords in loop
    // over elems later
    MB_CHK_SET_ERR( read_iface->get_node_coords( 3, 2 * ( nquads + 1 ), 0, start_vert, coords ),
                    "Error in get_node_arrays" );
    // Create elems
    MB_CHK_SET_ERR( read_iface->get_element_connect( nquads, 4, MBQUAD, 0, start_elem, connect ),
                    "Error in get_element_connect" );
    for( int i = 0; i < nquads; i++ )
    {
        coords[0][2 * i] = coords[0][2 * i + 1] = (double)i;  // x values are all i
        coords[1][2 * i]                        = 0.0;
        coords[1][2 * i + 1]                    = 1.0;          // y coords
        coords[2][2 * i] = coords[2][2 * i + 1] = (double)0.0;  // z values, all zero (2d mesh)
        EntityHandle quad_v                     = start_vert + 2 * i;
        for( int j = 0; j < 4; j++ )
            connect[4 * i + j] = quad_v + j;  // Connectivity of each quad is a sequence starting from quad_v
    }
    // Last two vertices
    // Cppcheck warning (false positive): variable coords is assigned a value that is never used
    coords[0][2 * nquads] = coords[0][2 * nquads + 1] = (double)nquads;
    coords[1][2 * nquads]                             = 0.0;
    coords[1][2 * nquads + 1]                         = 1.0;          // y coords
    coords[2][2 * nquads] = coords[2][2 * nquads + 1] = (double)0.0;  // z values, all zero (2d mesh)
 
    // Now delete nholes elements, spaced approximately equally through mesh, so contiguous size is
    // about nquads/(nholes + 1) reinterpret start_elem as the next element to be deleted
    int de = nquads / ( nholes + 1 );
    for( int i = 0; i < nholes; i++ )
    {
        start_elem += de;
        MB_CHK_SET_ERR( mbImpl->delete_entities( &start_elem, 1 ), "Error in delete_entities" );
    }
 
    return MB_SUCCESS;
}

References moab::Interface::delete_entities(), moab::ReadUtilIface::get_element_connect(), moab::ReadUtilIface::get_node_coords(), MB_CHK_SET_ERR, MB_SUCCESS, MBQUAD, and moab::Interface::query_interface().

Referenced by main().

main()

int main	(	int	argc,
		char **	argv
	)

Definition at line 77 of file DirectAccessWithHoles.cpp.

{
    // Get MOAB instance
    Interface* mbImpl = new( std::nothrow ) Core;
    if( NULL == mbImpl ) return 1;
 
    int nquads = 1000, nholes = 1;
 
    // Parse options
    ProgOptions opts;
    opts.addOpt< int >( string( "nquads,n" ), string( "Number of quads in the mesh (default = 1000" ), &nquads );
    opts.addOpt< int >( string( "holes,H" ), string( "Number of holes in the element handle space (default = 1" ),
                        &nholes );
    opts.parseCommandLine( argc, argv );
    if( nholes >= nquads )
    {
        cerr << "Number of holes needs to be < number of elements." << endl;
        return 1;
    }
 
    // Create simple structured mesh with hole, but using unstructured representation
    MB_CHK_SET_ERR( create_mesh_with_holes( mbImpl, nquads, nholes ), "Trouble creating mesh" );
 
    // Get all vertices and non-vertex entities
    Range verts, elems;
    MB_CHK_SET_ERR( mbImpl->get_entities_by_handle( 0, elems ), "Trouble getting all entities" );
    verts = elems.subset_by_dimension( 0 );
    elems -= verts;
 
    // Create tag1 (element-based avg), tag2 (vertex-based avg), tag3 (# connected verts)
    Tag tag1, tag2, tag3;
    MB_CHK_SET_ERR( mbImpl->tag_get_handle( "tag1", 3, MB_TYPE_DOUBLE, tag1, MB_TAG_DENSE | MB_TAG_CREAT ),
                    "Trouble creating tag1" );
    double def_val[3] = { 0.0, 0.0, 0.0 };  // Need a default value for tag2 because we sum into it
    MB_CHK_SET_ERR( mbImpl->tag_get_handle( "tag2", 3, MB_TYPE_DOUBLE, tag2, MB_TAG_DENSE | MB_TAG_CREAT, def_val ),
                    "Trouble creating tag2" );
    int def_val_int = 0;  // Need a default value for tag3 because we increment it
    MB_CHK_SET_ERR( mbImpl->tag_get_handle( "tag3", 1, MB_TYPE_INTEGER, tag3, MB_TAG_DENSE | MB_TAG_CREAT,
                                            &def_val_int ),
                    "Trouble creating tag3" );
 
    // Get connectivity, coordinate, tag, and adjacency iterators
    EntityHandle* conn_ptr;
    double *x_ptr, *y_ptr, *z_ptr, *tag1_ptr, *tag2_ptr;
    int* tag3_ptr;
 
    // First vertex is at start of range (ranges are sorted), and is offset for vertex index
    // calculation
    EntityHandle first_vert = *verts.begin();
 
    // When iterating over elements, each chunk can have a different # vertices; also, count tells
    // you how many elements are in the current chunk
    int vpere, count;
 
    // Get coordinates iterator, just need this once because we know verts handle space doesn't have
    // holes
    MB_CHK_SET_ERR( mbImpl->coords_iterate( verts.begin(), verts.end(), x_ptr, y_ptr, z_ptr, count ),
                    "Error in coords_iterate" );
    assert( count == (int)verts.size() );  // Should end up with just one contiguous chunk of vertices
 
    // Iterate through elements, computing midpoint based on vertex positions, set on element-based
    // tag1 Control loop by iterator over elem range
    Range::iterator e_it = elems.begin();
 
    while( e_it != elems.end() )
    {
        // Get conn_ptr, tag1_ptr for next contiguous chunk of element handles, and coords pointers
        // for all verts
        MB_CHK_SET_ERR( mbImpl->connect_iterate( e_it, elems.end(), conn_ptr, vpere, count ),
                        "Error in connect_iterate" );
        MB_CHK_SET_ERR( mbImpl->tag_iterate( tag1, e_it, elems.end(), count, (void*&)tag1_ptr ),
                        "Error in tag1_iterate" );
 
        // Iterate over elements in this chunk
        for( int i = 0; i < count; i++ )
        {
            tag1_ptr[0] = tag1_ptr[1] = tag1_ptr[2] = 0.0;  // Initialize position for this element
            for( int j = 0; j < vpere; j++ )
            {                                            // Loop over vertices in this element
                int v_index = conn_ptr[j] - first_vert;  // vert index is just the offset from first vertex
                tag1_ptr[0] += x_ptr[v_index];
                tag1_ptr[1] += y_ptr[v_index];  // Sum vertex positions into tag1...
                tag1_ptr[2] += z_ptr[v_index];
            }
            tag1_ptr[0] /= vpere;
            tag1_ptr[1] /= vpere;  // Then normalize
            tag1_ptr[2] /= vpere;
 
            // Done with this element; advance connect_ptr and tag1_ptr to next element
            conn_ptr += vpere;
            tag1_ptr += 3;
        }  // Loop over elements in chunk
 
        // Done with chunk; advance range iterator by count; will skip over gaps in range
        e_it += count;
    }  // While loop over all elements
 
    // Iterate through vertices, summing positions into tag2 on connected elements and incrementing
    // vertex count Iterate over chunks the same as elements, even though we know we have only one
    // chunk here, just to show how it's done
 
    // Create a std::map from EntityHandle (first entity handle in chunk) to
    // tag_struct (ptrs to start of avg/#verts tags for that chunk); then for a given entity handle,
    // we can quickly find the chunk it's in using map::lower_bound; could have set up this map in
    // earlier loop over elements, but do it here for clarity
 
    map< EntityHandle, tag_struct > elem_map;
    e_it = elems.begin();
    while( e_it != elems.end() )
    {
        tag_struct ts = { NULL, NULL };
        MB_CHK_SET_ERR( mbImpl->tag_iterate( tag2, e_it, elems.end(), count, (void*&)ts.avg_ptr ),
                        "Error in tag2_iterate" );
        MB_CHK_SET_ERR( mbImpl->tag_iterate( tag3, e_it, elems.end(), count, (void*&)ts.nv_ptr ),
                        "Error in tag3_iterate" );
        elem_map[*e_it] = ts;
        e_it += count;
    }
 
    // Call a vertex-quad adjacencies function to generate vertex-element adjacencies in MOAB
    Range::iterator v_it = verts.begin();
    Range dum_range;
    MB_CHK_SET_ERR( mbImpl->get_adjacencies( &( *v_it ), 1, 2, false, dum_range ), "Error in get_adjacencies" );
    const vector< EntityHandle >** adjs_ptr;
    while( v_it != verts.end() )
    {
        // Get coords ptrs, adjs_ptr; can't set tag2_ptr by direct access, because of hole in
        // element handle space
        MB_CHK_SET_ERR( mbImpl->coords_iterate( v_it, verts.end(), x_ptr, y_ptr, z_ptr, count ),
                        "Error in coords_iterate" );
        MB_CHK_SET_ERR( mbImpl->adjacencies_iterate( v_it, verts.end(), adjs_ptr, count ),
                        "Error in adjacencies_iterate" );
 
        for( int v = 0; v < count; v++ )
        {
            const vector< EntityHandle >* avec = *( adjs_ptr + v );
            for( vector< EntityHandle >::const_iterator ait = avec->begin(); ait != avec->end(); ++ait )
            {
                // Get chunk that this element resides in; upper_bound points to the first element
                // strictly > key, so get that and decrement (would work the same as lower_bound
                // with an if-test and decrement)
                map< EntityHandle, tag_struct >::iterator mit = elem_map.upper_bound( *ait );
                --mit;
                // Index of *ait in that chunk
                int a_ind     = *ait - ( *mit ).first;
                tag_struct ts = ( *mit ).second;
                ts.avg_ptr[3 * a_ind + 0] += x_ptr[v];  // Tag on each element is 3 doubles, x/y/z
                ts.avg_ptr[3 * a_ind + 1] += y_ptr[v];
                ts.avg_ptr[3 * a_ind + 2] += z_ptr[v];
                ts.nv_ptr[a_ind]++;  // Increment the vertex count
            }
        }
 
        v_it += count;
    }
 
    // Normalize tag2 by vertex count; loop over elements using same approach as before
    // At the same time, compare values of tag1 and tag2
    e_it = elems.begin();
    while( e_it != elems.end() )
    {
        // Get conn_ptr, tag1_ptr for next contiguous chunk of element handles, and coords pointers
        // for all verts
        MB_CHK_SET_ERR( mbImpl->tag_iterate( tag1, e_it, elems.end(), count, (void*&)tag1_ptr ),
                        "Error in tag1_iterate" );
        MB_CHK_SET_ERR( mbImpl->tag_iterate( tag2, e_it, elems.end(), count, (void*&)tag2_ptr ),
                        "Error in tag2_iterate" );
        MB_CHK_SET_ERR( mbImpl->tag_iterate( tag3, e_it, elems.end(), count, (void*&)tag3_ptr ),
                        "Error in tag3_iterate" );
 
        // Iterate over elements in this chunk
        for( int i = 0; i < count; i++ )
        {
            for( int j = 0; j < 3; j++ )
                tag2_ptr[3 * i + j] /= (double)tag3_ptr[i];  // Normalize by # verts
            if( tag1_ptr[3 * i] != tag2_ptr[3 * i] || tag1_ptr[3 * i + 1] != tag2_ptr[3 * i + 1] ||
                tag1_ptr[3 * i + 2] != tag2_ptr[3 * i + 2] )
                cout << "Tag1, tag2 disagree for element " << *e_it + i << endl;
        }
 
        e_it += count;
    }
 
    // Ok, we're done, shut down MOAB
    delete mbImpl;
 
    return 0;
}

References ProgOptions::addOpt(), moab::Interface::adjacencies_iterate(), tag_struct::avg_ptr, moab::Range::begin(), moab::Interface::connect_iterate(), moab::Interface::coords_iterate(), create_mesh_with_holes(), moab::Range::end(), moab::Interface::get_adjacencies(), moab::Interface::get_entities_by_handle(), MB_CHK_SET_ERR, MB_TAG_CREAT, MB_TAG_DENSE, MB_TYPE_DOUBLE, MB_TYPE_INTEGER, tag_struct::nv_ptr, ProgOptions::parseCommandLine(), moab::Range::size(), moab::Range::subset_by_dimension(), moab::Interface::tag_get_handle(), and moab::Interface::tag_iterate().