DirectAccessNoHoles.cpp File Reference

Example demonstrating direct access to MOAB data for performance optimization. More...

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

Include dependency graph for DirectAccessNoHoles.cpp:

Go to the source code of this file.

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

Detailed Description

Example demonstrating direct access to MOAB data for performance optimization.

This example shows how to:

• Create a structured mesh with contiguous entity handles
• Use direct access iterators to avoid API overhead
• Access vertex coordinates, element connectivity, and tag data directly
• Work with adjacency lists for vertex-to-element relationships
• Perform efficient mesh operations using pointer arithmetic

The program creates a 1D row of quad elements with contiguous handles, then demonstrates direct access to MOAB's internal data structures for maximum performance when the mesh topology is static.

Author

MOAB Team

Date

2024

Parameters

[in]	argc	Number of command line arguments
[in]	argv	Command line arguments array
[in]	-nquads	Number of quads in the mesh (default: 1000)

Returns

0 on success, 1 on failure

Usage:

./DirectAccessNoHoles [-nquads <# quads>]

Example:

./DirectAccessNoHoles -nquads 500

See also

Core, Interface, ReadUtilIface, ProgOptions

Use direct access to MOAB data to avoid calling through API
This example creates a 1d row of quad elements, such that all quad and vertex handles are contiguous in the handle space and in the database. Then it shows how to get access to pointers to MOAB-native data for vertex coordinates, quad connectivity, tag storage, and vertex to quad adjacency lists. This allows applications to access this data directly without going through MOAB's API. In cases where the mesh is not changing (or only mesh vertices are moving), this can save significant execution time in applications.

*  ----------------------
*  |      |      |      |
*  |      |      |      | ...
*  |      |      |      |
*  ----------------------
* 

1. Initialize MOAB
   
2. Create a quad mesh, 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. Iterate through vertices, summing positions into tag2 on connected quads and incrementing vertex count
7. Iterate through quads, normalizing tag2 by vertex count and comparing values of tag1 and tag2

To compile:
make DirectAccessNoHoles
To run: ./DirectAccessNoHoles [-nquads <# quads>]

Definition in file DirectAccessNoHoles.cpp.

Function Documentation

create_mesh_no_holes()

ErrorCode create_mesh_no_holes	(	Interface *	mbImpl,
		int	nquads
	)

Definition at line 206 of file DirectAccessNoHoles.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 2(nquads+1) because they're in a 1d row; will initialize coords in loop
    // over quads later
    MB_CHK_SET_ERR( read_iface->get_node_coords( 3, 2 * ( nquads + 1 ), 0, start_vert, coords ),
                    "Error in get_node_arrays" );
    // Create quads
    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;
        connect[4 * i + 0]                      = quad_v;
        connect[4 * i + 1]                      = quad_v + 2;
        connect[4 * i + 2]                      = quad_v + 3;
        connect[4 * i + 3]                      = quad_v + 1;
    }
 
    // 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)
 
    // Call a vertex-quad adjacencies function to generate vertex-element adjacencies in MOAB
    Range dum_range;
    MB_CHK_SET_ERR( mbImpl->get_adjacencies( &start_vert, 1, 2, false, dum_range ), "Error in get_adjacencies" );
    assert( !dum_range.empty() );
 
    return MB_SUCCESS;
}

References moab::Range::empty(), moab::Interface::get_adjacencies(), 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 83 of file DirectAccessNoHoles.cpp.

{
    // Get MOAB instance
    Interface* mbImpl = new( std::nothrow ) Core;
    if( NULL == mbImpl ) return 1;
 
    int nquads = 1000;
 
    // Parse options
    ProgOptions opts;
    opts.addOpt< int >( string( "nquads,n" ), string( "Number of quads in the mesh (default = 1000" ), &nquads );
    opts.parseCommandLine( argc, argv );
 
    // Create simple structured mesh with hole, but using unstructured representation
    MB_CHK_SET_ERR( create_mesh_no_holes( mbImpl, nquads ), "Trouble creating mesh" );
 
    // Get all vertices and non-vertex entities
    Range verts, quads;
    MB_CHK_SET_ERR( mbImpl->get_entities_by_handle( 0, quads ), "Trouble getting all entities" );
    verts = quads.subset_by_dimension( 0 );
    quads -= 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 pointers to connectivity, coordinate, tag, and adjacency arrays; each of these returns a
    // count, which should be compared to the # entities you expect to verify there's only one chunk
    // (no holes)
    int count, vpere;
    EntityHandle* conn_ptr;
    MB_CHK_SET_ERR( mbImpl->connect_iterate( quads.begin(), quads.end(), conn_ptr, vpere, count ),
                    "Error in connect_iterate" );
    assert( count == (int)quads.size() );  // Should end up with just one contiguous chunk of quads
 
    double *x_ptr, *y_ptr, *z_ptr;
    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
 
    double *tag1_ptr, *tag2_ptr;
    int* tag3_ptr;
    MB_CHK_SET_ERR( mbImpl->tag_iterate( tag1, quads.begin(), quads.end(), count, (void*&)tag1_ptr ),
                    "Error in tag1_iterate" );
    assert( count == (int)quads.size() );  // Should end up with just one contiguous chunk of quads
    MB_CHK_SET_ERR( mbImpl->tag_iterate( tag2, quads.begin(), quads.end(), count, (void*&)tag2_ptr ),
                    "Error in tag2_iterate" );
    assert( count == (int)quads.size() );  // Should end up with just one contiguous chunk of quads
    MB_CHK_SET_ERR( mbImpl->tag_iterate( tag3, quads.begin(), quads.end(), count, (void*&)tag3_ptr ),
                    "Error in tag3_iterate" );
    assert( count == (int)quads.size() );  // Should end up with just one contiguous chunk of quads
 
    const vector< EntityHandle >** adjs_ptr;
    MB_CHK_SET_ERR( mbImpl->adjacencies_iterate( verts.begin(), verts.end(), adjs_ptr, count ),
                    "Error in adjacencies_iterate" );
    assert( count == (int)verts.size() );  // Should end up with just one contiguous chunk of vertices
    // Start_ handles used to compute indices into vertex/quad arrays; can use direct subtraction
    // because we know there aren't any holes in the handle spaces for verts or quads
    EntityHandle start_vert = *verts.begin(), start_quad = *quads.begin();
 
    // Iterate over elements, computing tag1 from coords positions
    for( int i = 0; i < nquads; i++ )
    {
        tag1_ptr[3 * i + 0] = tag1_ptr[3 * i + 1] = tag1_ptr[3 * i + 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[vpere * i + j] - start_vert;  // vert index is just the offset from start vertex
            tag1_ptr[3 * i + 0] += x_ptr[v_index];
            tag1_ptr[3 * i + 1] += y_ptr[v_index];  // Sum vertex positions into tag1...
            tag1_ptr[3 * i + 2] += z_ptr[v_index];
        }
        tag1_ptr[3 * i + 0] /= vpere;
        tag1_ptr[3 * i + 1] /= vpere;  // Then normalize
        tag1_ptr[3 * i + 2] /= vpere;
    }  // Loop over elements in chunk
 
    // Iterate through vertices, summing positions into tag2 on connected elements and incrementing
    // vertex count
    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 )
        {
            // *ait is the quad handle, its index is computed by subtracting the start quad handle
            int a_ind = *ait - start_quad;
            tag2_ptr[3 * a_ind + 0] += x_ptr[v];  // Tag on each element is 3 doubles, x/y/z
            tag2_ptr[3 * a_ind + 1] += y_ptr[v];
            tag2_ptr[3 * a_ind + 2] += z_ptr[v];
            tag3_ptr[a_ind]++;  // Increment the vertex count
        }
    }
 
    // Normalize tag2 by vertex count (tag3); loop over elements using same approach as before
    // At the same time, compare values of tag1 and tag2
    int n_dis = 0;
    for( Range::iterator q_it = quads.begin(); q_it != quads.end(); ++q_it )
    {
        int i = *q_it - start_quad;
        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 " << *q_it + i << endl;
            n_dis++;
        }
    }
    if( !n_dis ) cout << "All tags agree, success!" << endl;
 
    // Ok, we're done, shut down MOAB
    delete mbImpl;
 
    return 0;
}

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