added benchmarks

# auxiliary build files

*.o

*.d

# final binaries

/moab-benchmark-measures

/moab-benchmark-surfaces-1

/moab-benchmark-surfaces-2

/moab-benchmark-surfaces-3

CXX := g++

CPPFLAGS := -MD -MP

CXXFLAGS := -std=c++11 -Wall -Wextra -pedantic -O3 -march=native -mtune=native -flto -fopenmp

LDFLAGS := -lm -lMOAB -lgomp -flto

SOURCES = $(wildcard *.cpp)

TARGETS = $(SOURCES:%.cpp=%)

all: $(TARGETS)

%.o: %.cpp

	$(CXX) $(CPPFLAGS) $(CXXFLAGS) -c -o $@ $<

$(TARGETS): %: %.o

	$(CXX) -o $@ $^ $(LDFLAGS)

clean:

	$(RM) *.[od] $(TARGETS)

-include $(SOURCES:%.cpp=%.d)

#include <iostream>

#include <stdexcept>

#include <cmath>

#include <vector>

#include <chrono>

#include "moab/Core.hpp"

#include "moab/Range.hpp"

#include "moab/CN.hpp"

#include <omp.h>

using namespace moab;

using namespace std;

using clock_type = std::chrono::steady_clock;

double getEntityMeasure(const Interface* mb, const EntityHandle handle)

{

    const EntityType entity_type = mb->type_from_handle(handle);

    if( entity_type == MBVERTEX )

    {

        // 0-dimensional measure of a vertex equals 1 by convention

        return 1;

    }

    if( entity_type == MBEDGE )

    {

        // extract subvertices of the triangle

        const EntityHandle* verts;

        int num_verts;

        ErrorCode rval = mb->get_connectivity( handle, verts, num_verts );

        if( rval != ErrorCode::MB_SUCCESS )

            throw std::runtime_error("mb->get_connectivity failed");

        if( num_verts != 2 )

            throw std::runtime_error(string("mb->get_connectivity returned unexpected number of vertices for an edge: ") + to_string(num_verts));

        // extract coordinates of the vertices

        double coords[2][3];

        mb->get_coords( verts, 3, reinterpret_cast<double*>(coords) );

        // subtract the 0-th coordinate from the other coordinates

        for (int i = 0; i < 3; i++)

            coords[1][i] -= coords[0][i];

        // calculate segment length

        return std::sqrt( coords[1][0]*coords[1][0] + coords[1][1]*coords[1][1] + coords[1][2]*coords[1][2] );

    }

    if( entity_type == MBTRI )

    {

        // extract subvertices of the triangle

        const EntityHandle* verts;

        int num_verts;

        ErrorCode rval = mb->get_connectivity( handle, verts, num_verts );

        if( rval != ErrorCode::MB_SUCCESS )

            throw std::runtime_error("mb->get_connectivity failed");

        if( num_verts != 3 )

            throw std::runtime_error(string("mb->get_connectivity returned unexpected number of vertices for a triangle: ") + to_string(num_verts));

        // extract coordinates of the vertices

        double coords[3][3];

        mb->get_coords( verts, 3, reinterpret_cast<double*>(coords) );

        // subtract the 0-th coordinate from the other coordinates

        for (int i = 0; i < 3; i++)

        {

            coords[1][i] -= coords[0][i];

            coords[2][i] -= coords[0][i];

        }

        // formula from http://math.stackexchange.com/a/128999

        const double c1 = coords[1][1] * coords[2][2] - coords[1][2] * coords[2][1];   // first component of the cross product

        const double c2 = coords[1][2] * coords[2][0] - coords[1][0] * coords[2][2];   // second component of the cross product

        const double c3 = coords[1][0] * coords[2][1] - coords[1][1] * coords[2][0];   // third component of the cross product

        return 0.5 * std::sqrt( c1 * c1 + c2 * c2 + c3 * c3 );

    }

    if( entity_type == MBTET )

    {

        // extract subvertices of the triangle

        const EntityHandle* verts;

        int num_verts;

        ErrorCode rval = mb->get_connectivity( handle, verts, num_verts );

        if( rval != ErrorCode::MB_SUCCESS )

            throw std::runtime_error("mb->get_connectivity failed");

        if( num_verts != 4 )

            throw std::runtime_error(string("mb->get_connectivity returned unexpected number of vertices for a tetrahedron: ") + to_string(num_verts));

        // extract coordinates of the vertices

        double coords[4][3];

        mb->get_coords( verts, 4, reinterpret_cast<double*>(coords) );

        // subtract the 0-th coordinate from the other coordinates

        for (int i = 0; i < 3; i++)

        {

            coords[1][i] -= coords[0][i];

            coords[2][i] -= coords[0][i];

            coords[3][i] -= coords[0][i];

        }

        // V = (1/6) * det(v1, v2, v3)

        const double det = coords[1][0] * coords[2][1] * coords[3][2] +

                           coords[1][1] * coords[2][2] * coords[3][0] +

                           coords[1][2] * coords[2][0] * coords[3][1] -

                         ( coords[1][2] * coords[2][1] * coords[3][0] +

                           coords[1][1] * coords[2][0] * coords[3][2] +

                           coords[1][0] * coords[2][2] * coords[3][1] );

        return 1.0 / 6.0 * std::abs( det );

    }

    else

    {

        throw std::runtime_error(string("getEntityMeasure is not implemented for entities of type '") + CN::EntityTypeName(entity_type) + "'.");

    }

}

void run_benchmark(Interface* mb, Range& ents)

{

    const int loops = 100;

    std::chrono::duration<double> duration;

    std::vector<double> results;

    for( int i = 0; i < loops; i++ )

    {

        // reset the results vector

        results.resize(ents.size(), 0);

        auto start = clock_type::now();

        if( omp_get_max_threads() == 1 )

        {

            for( Range::iterator it = ents.begin(); it != ents.end(); ++it )

            {

                const int idx = mb->id_from_handle(*it) - 1;

                results[idx] = getEntityMeasure(mb, *it);

            }

        }

        else

        {

            #pragma omp parallel for default(none) shared(ents, mb, results)

            for( Range::iterator it = ents.begin(); it != ents.end(); ++it )

            {

                const int idx = mb->id_from_handle(*it) - 1;

                results[idx] = getEntityMeasure(mb, *it);

            }

        }

        auto end = clock_type::now();

        duration += end - start;

    }

    cout << "benchmark time: " << duration.count() * 1e3 / loops << " ms" << endl;

    double total_measure = 0;

    for (auto v : results) total_measure += v;

    cout << "total measure = " << total_measure << endl;

}

int main( int argc, char** argv )

{

    string test_file_name;

    // Need option handling here for input filename

    if( argc == 2 )

        test_file_name = argv[1];

    else

    {

        cerr << "usage: " << argv[0] << " FILE.vtk" << endl;

        return 1;

    }

    // Get MOAB instance

    Interface* mb = new( std::nothrow ) Core;

    if( NULL == mb ) return 1;

    // Load the mesh from vtk file

    ErrorCode rval = mb->load_mesh( test_file_name.c_str() );MB_CHK_ERR( rval );

    // Get verts entities, by type

    Range verts;

    rval = mb->get_entities_by_type( 0, MBVERTEX, verts );MB_CHK_ERR( rval );

    // Get edge entities, by type

    Range edges;

    rval = mb->get_entities_by_type( 0, MBEDGE, edges );MB_CHK_ERR( rval );

    // Get faces, by dimension, so we stay generic to entity type

    Range faces;

    rval = mb->get_entities_by_dimension( 0, 2, faces );MB_CHK_ERR( rval );

    // Get cells, by dimension, so we stay generic to entity type

    Range cells;

    rval = mb->get_entities_by_dimension( 0, 3, cells );MB_CHK_ERR( rval );

    // Output the number of entities

    if( cells.size() > 0 )

    {

        cout << "Number of vertices:   " << verts.size() << endl;

        cout << "Number of edges:      " << edges.size() << endl;

        cout << "Number of faces:      " << faces.size() << endl;

        cout << "Number of cells (3D): " << cells.size() << endl;

    }

    else

    {

        cout << "Number of vertices:   " << verts.size() << endl;

        cout << "Number of edges:      " << edges.size() << endl;

        cout << "Number of cells (2D): " << faces.size() << endl;

    }

    if( cells.size() > 0 )

        run_benchmark(mb, cells);

    else

        run_benchmark(mb, faces);

    delete mb;

    return 0;

}

moab-benchmark-surfaces-2

 No newline at end of file

#include <iostream>

#include <stdexcept>

#include <cmath>

#include <vector>

#include <chrono>

#include "moab/Core.hpp"

#include "moab/Range.hpp"

#include "moab/CN.hpp"

#include <omp.h>

using namespace moab;

using namespace std;

using clock_type = std::chrono::steady_clock;

double getEntityMeasure(const Interface* mb, const EntityHandle handle)

{

    const EntityType entity_type = mb->type_from_handle(handle);

    if( entity_type == MBVERTEX )

    {

        // 0-dimensional measure of a vertex equals 1 by convention

        return 1;

    }

    if( entity_type == MBEDGE )

    {

        // extract subvertices of the triangle

        const EntityHandle* verts;

        int num_verts;

        ErrorCode rval = mb->get_connectivity( handle, verts, num_verts );

        if( rval != ErrorCode::MB_SUCCESS )

            throw std::runtime_error("mb->get_connectivity failed");

        if( num_verts != 2 )

            throw std::runtime_error(string("mb->get_connectivity returned unexpected number of vertices for an edge: ") + to_string(num_verts));

        // extract coordinates of the vertices

        double coords[2][3];

        mb->get_coords( verts, 3, reinterpret_cast<double*>(coords) );

        // subtract the 0-th coordinate from the other coordinates

        for (int i = 0; i < 3; i++)

            coords[1][i] -= coords[0][i];

        // calculate segment length

        return std::sqrt( coords[1][0]*coords[1][0] + coords[1][1]*coords[1][1] + coords[1][2]*coords[1][2] );

    }

    if( entity_type == MBTRI )

    {

        // extract subvertices of the triangle

        const EntityHandle* verts;

        int num_verts;

        ErrorCode rval = mb->get_connectivity( handle, verts, num_verts );

        if( rval != ErrorCode::MB_SUCCESS )

            throw std::runtime_error("mb->get_connectivity failed");

        if( num_verts != 3 )

            throw std::runtime_error(string("mb->get_connectivity returned unexpected number of vertices for a triangle: ") + to_string(num_verts));

        // extract coordinates of the vertices

        double coords[3][3];

        mb->get_coords( verts, 3, reinterpret_cast<double*>(coords) );

        // subtract the 0-th coordinate from the other coordinates

        for (int i = 0; i < 3; i++)

        {

            coords[1][i] -= coords[0][i];

            coords[2][i] -= coords[0][i];

        }

        // formula from http://math.stackexchange.com/a/128999

        const double c1 = coords[1][1] * coords[2][2] - coords[1][2] * coords[2][1];   // first component of the cross product

        const double c2 = coords[1][2] * coords[2][0] - coords[1][0] * coords[2][2];   // second component of the cross product

        const double c3 = coords[1][0] * coords[2][1] - coords[1][1] * coords[2][0];   // third component of the cross product

        return 0.5 * std::sqrt( c1 * c1 + c2 * c2 + c3 * c3 );

    }

    if( entity_type == MBTET )

    {

        // extract subvertices of the triangle

        const EntityHandle* verts;

        int num_verts;

        ErrorCode rval = mb->get_connectivity( handle, verts, num_verts );

        if( rval != ErrorCode::MB_SUCCESS )

            throw std::runtime_error("mb->get_connectivity failed");

        if( num_verts != 4 )

            throw std::runtime_error(string("mb->get_connectivity returned unexpected number of vertices for a tetrahedron: ") + to_string(num_verts));

        // extract coordinates of the vertices

        double coords[4][3];

        mb->get_coords( verts, 4, reinterpret_cast<double*>(coords) );

        // subtract the 0-th coordinate from the other coordinates

        for (int i = 0; i < 3; i++)

        {

            coords[1][i] -= coords[0][i];

            coords[2][i] -= coords[0][i];

            coords[3][i] -= coords[0][i];

        }

        // V = (1/6) * det(v1, v2, v3)

        const double det = coords[1][0] * coords[2][1] * coords[3][2] +

                           coords[1][1] * coords[2][2] * coords[3][0] +

                           coords[1][2] * coords[2][0] * coords[3][1] -

                         ( coords[1][2] * coords[2][1] * coords[3][0] +

                           coords[1][1] * coords[2][0] * coords[3][2] +

                           coords[1][0] * coords[2][2] * coords[3][1] );

        return 1.0 / 6.0 * std::abs( det );

    }

    else

    {

        throw std::runtime_error(string("getEntityMeasure is not implemented for entities of type '") + CN::EntityTypeName(entity_type) + "'.");

    }

}

bool hasSubvertex( const Interface* mb, const EntityHandle& entity, const EntityHandle& vertex )

{

    // extract subvertices of the entity

    const EntityHandle* subverts;

    int num_verts;

    ErrorCode rval = mb->get_connectivity( entity, subverts, num_verts );

    if( rval != ErrorCode::MB_SUCCESS )

        throw std::runtime_error("mb->get_connectivity failed");

    for( int v = 0; v < num_verts; v++ ) {

        if( subverts[v] == vertex )

            return true;

    }

    return false;

}

double getSphereMeasure(Interface* mb, const EntityHandle vertex, int mesh_dim)

{

    double s = 0.0;

    Range adj_cells;

    ErrorCode rval = mb->get_adjacencies( &vertex, 1, mesh_dim, false, adj_cells );

    if( rval != ErrorCode::MB_SUCCESS )

        throw std::runtime_error("mb->get_adjacencies failed");

    for( auto& cell : adj_cells )

    {

        Range adj_faces;

        ErrorCode rval = mb->get_adjacencies( &cell, 1, mesh_dim - 1, false, adj_faces );

        if( rval != ErrorCode::MB_SUCCESS )

            throw std::runtime_error("mb->get_adjacencies failed");

        // NOTE: This is a general version, which is slow because of the two nested loops.

        //       The TNL benchmark uses a shortcut based on local numbering on simplices

        //       (assuming that opposite vertex and face have the same local index) which

        //       involves only one nested loop to find the local index of the vertex on

        //       the cell.

        for( auto& face : adj_faces )

        {

            if( ! hasSubvertex( mb, face, vertex ) ) {

                s += getEntityMeasure( mb, face );

            }

        }

    }

    return s;

}

void run_benchmark(Interface* mb, const Range& verts, int mesh_dim)

{

    const int loops = 100;

    std::chrono::duration<double> duration;

    std::vector<double> results;

    for( int i = 0; i < loops; i++ )

    {

        // reset the results vector

        results.clear();

        results.resize(verts.size(), 0);

        auto start = clock_type::now();

        if( omp_get_max_threads() == 1 )

        {

            for( Range::iterator it = verts.begin(); it != verts.end(); ++it )

            {

                const int idx = mb->id_from_handle(*it) - 1;

                results[idx] = getSphereMeasure(mb, *it, mesh_dim);

            }

        }

        else

        {

            #pragma omp parallel for default(none) shared(verts, mb, results, mesh_dim)

            for( Range::iterator it = verts.begin(); it != verts.end(); ++it )

            {

                const int idx = mb->id_from_handle(*it) - 1;

                results[idx] = getSphereMeasure(mb, *it, mesh_dim);

            }

        }

        auto end = clock_type::now();

        duration += end - start;

    }

    cout << "benchmark time: " << duration.count() * 1e3 / loops << " ms" << endl;

    double total_measure = 0;

    for (auto v : results) total_measure += v;

    cout << "total measure = " << total_measure << endl;

}

int main( int argc, char** argv )

{

    string test_file_name;

    // Need option handling here for input filename

    if( argc == 2 )

        test_file_name = argv[1];

    else

    {

        cerr << "usage: " << argv[0] << " FILE.vtk" << endl;

        return 1;

    }

    // Get MOAB instance

    Interface* mb = new( std::nothrow ) Core;

    if( NULL == mb ) return 1;

    // Load the mesh from vtk file

    ErrorCode rval = mb->load_mesh( test_file_name.c_str() );MB_CHK_ERR( rval );

    // Get verts entities, by type

    Range verts;

    rval = mb->get_entities_by_type( 0, MBVERTEX, verts );MB_CHK_ERR( rval );

    // Get cells, by dimension, so we stay generic to entity type

    Range cells;

    rval = mb->get_entities_by_dimension( 0, 3, cells );MB_CHK_ERR( rval );

    // Create missing faces

    for( auto& cell : cells )

    {

        Range adj_faces;

        rval = mb->get_adjacencies( &cell, 1, 2, true, adj_faces );MB_CHK_ERR( rval );

    }

    // Get faces, by dimension, so we stay generic to entity type

    Range faces;

    rval = mb->get_entities_by_dimension( 0, 2, faces );MB_CHK_ERR( rval );

    // Create missing edges

    for( auto& face : faces )

    {

        Range adj_edges;

        rval = mb->get_adjacencies( &face, 1, 1, true, adj_edges );MB_CHK_ERR( rval );

    }

    // Get edge entities, by type

    Range edges;

    rval = mb->get_entities_by_type( 0, MBEDGE, edges );MB_CHK_ERR( rval );

    // Output the number of entities

    if( cells.size() > 0 )

    {

        cout << "Number of vertices:   " << verts.size() << endl;

        cout << "Number of edges:      " << edges.size() << endl;

        cout << "Number of faces:      " << faces.size() << endl;

        cout << "Number of cells (3D): " << cells.size() << endl;

    }

    else

    {

        cout << "Number of vertices:   " << verts.size() << endl;

        cout << "Number of edges:      " << edges.size() << endl;

        cout << "Number of cells (2D): " << faces.size() << endl;

    }

    run_benchmark(mb, verts, (cells.size() > 0) ? 3 : 2);

    delete mb;

    return 0;

}