CN.cpp

Go to the documentation of this file.

/**
 * MOAB, a Mesh-Oriented datABase, is a software component for creating,
 * storing and accessing finite element mesh data.
 *
 * Copyright 2004 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.
 *
 */

#include "moab/CN.hpp"
#include "MBCNArrays.hpp"
#include "MBCN.h"
#include <cassert>
#include <cstring>
#include <iterator>

namespace moab
{

const char* CN::entityTypeNames[] = { "Vertex", "Edge",  "Tri", "Quad",       "Polygon",   "Tet",    "Pyramid",
                                      "Prism",  "Knife", "Hex", "Polyhedron", "EntitySet", "MaxType" };

short int CN::numberBasis = 0;

short int CN::permuteVec[MBMAXTYPE][3][MAX_SUB_ENTITIES + 1];
short int CN::revPermuteVec[MBMAXTYPE][3][MAX_SUB_ENTITIES + 1];

const DimensionPair CN::TypeDimensionMap[] = {
    DimensionPair( MBVERTEX, MBVERTEX ),       DimensionPair( MBEDGE, MBEDGE ),
    DimensionPair( MBTRI, MBPOLYGON ),         DimensionPair( MBTET, MBPOLYHEDRON ),
    DimensionPair( MBENTITYSET, MBENTITYSET ), DimensionPair( MBMAXTYPE, MBMAXTYPE ) };

short CN::increasingInts[] = { 0,  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 };

//! set the basis of the numbering system; may or may not do things besides setting the
//! member variable
void CN::SetBasis( const int in_basis )
{
    numberBasis = in_basis;
}

/// Get the dimension pair corresponding to a dimension
DimensionPair CN::getDimPair( int entity_type )
{
    return TypeDimensionMap[entity_type];
}

//! return a type for the given name
EntityType CN::EntityTypeFromName( const char* name )
{
    for( EntityType i = MBVERTEX; i < MBMAXTYPE; i++ )
    {
        if( 0 == strcmp( name, entityTypeNames[i] ) ) return i;
    }

    return MBMAXTYPE;
}

void CN::SubEntityNodeIndices( const EntityType this_topo,
                               const int num_nodes,
                               const int sub_dimension,
                               const int sub_index,
                               EntityType& subentity_topo,
                               int& num_sub_entity_nodes,
                               int sub_entity_conn[] )
{
    // If asked for a node, the special case...
    if( sub_dimension == 0 )
    {
        assert( sub_index < num_nodes );
        subentity_topo       = MBVERTEX;
        num_sub_entity_nodes = 1;
        sub_entity_conn[0]   = sub_index;
        return;
    }

    const int ho_bits    = HasMidNodes( this_topo, num_nodes );
    subentity_topo       = SubEntityType( this_topo, sub_dimension, sub_index );
    num_sub_entity_nodes = VerticesPerEntity( subentity_topo );
    const short* corners = mConnectivityMap[this_topo][sub_dimension - 1].conn[sub_index];
    std::copy( corners, corners + num_sub_entity_nodes, sub_entity_conn );

    int sub_sub_corners[MAX_SUB_ENTITY_VERTICES];
    int side, sense, offset;
    for( int dim = 1; dim <= sub_dimension; ++dim )
    {
        if( !( ho_bits & ( 1 << dim ) ) ) continue;

        const short num_mid = NumSubEntities( subentity_topo, dim );
        for( int i = 0; i < num_mid; ++i )
        {
            const EntityType sub_sub_topo = SubEntityType( subentity_topo, dim, i );
            const int sub_sub_num_vert    = VerticesPerEntity( sub_sub_topo );
            SubEntityVertexIndices( subentity_topo, dim, i, sub_sub_corners );

            for( int j = 0; j < sub_sub_num_vert; ++j )
                sub_sub_corners[j] = corners[sub_sub_corners[j]];
            SideNumber( this_topo, sub_sub_corners, sub_sub_num_vert, dim, side, sense, offset );
            sub_entity_conn[num_sub_entity_nodes++] = HONodeIndex( this_topo, num_nodes, dim, side );
        }
    }
}

//! return the vertices of the specified sub entity
//! \param parent_conn Connectivity of parent entity
//! \param parent_type Entity type of parent entity
//! \param sub_dimension Dimension of sub-entity being queried
//! \param sub_index Index of sub-entity being queried
//! \param sub_entity_conn Connectivity of sub-entity, based on parent_conn and canonical
//!           ordering for parent_type
//! \param num_sub_vertices Number of vertices in sub-entity
void CN::SubEntityConn( const void* parent_conn,
                        const EntityType parent_type,
                        const int sub_dimension,
                        const int sub_index,
                        void* sub_entity_conn,
                        int& num_sub_vertices )
{
    static int sub_indices[MAX_SUB_ENTITY_VERTICES];

    SubEntityVertexIndices( parent_type, sub_dimension, sub_index, sub_indices );

    num_sub_vertices       = VerticesPerEntity( SubEntityType( parent_type, sub_dimension, sub_index ) );
    void** parent_conn_ptr = static_cast< void** >( const_cast< void* >( parent_conn ) );
    void** sub_conn_ptr    = static_cast< void** >( sub_entity_conn );
    for( int i = 0; i < num_sub_vertices; i++ )
        sub_conn_ptr[i] = parent_conn_ptr[sub_indices[i]];
}

//! given an entity and a target dimension & side number, get that entity
short int CN::AdjacentSubEntities( const EntityType this_type,
                                   const int* source_indices,
                                   const int num_source_indices,
                                   const int source_dim,
                                   const int target_dim,
                                   std::vector< int >& index_list,
                                   const int operation_type )
{
    // first get all the vertex indices
    std::vector< int > tmp_indices;
    const int* it1 = source_indices;

    assert( source_dim <= 3 && target_dim >= 0 && target_dim <= 3 &&
            // make sure we're not stepping off the end of the array;
            ( ( source_dim > 0 && *it1 < mConnectivityMap[this_type][source_dim - 1].num_sub_elements ) ||
              ( source_dim == 0 &&
                *it1 < mConnectivityMap[this_type][Dimension( this_type ) - 1].num_corners_per_sub_element[0] ) ) &&
            *it1 >= 0 );

#define MUC CN::mUpConnMap[this_type][source_dim][target_dim]

    // if we're looking for the vertices of a single side, return them in
    // the canonical ordering; otherwise, return them in sorted order
    if( num_source_indices == 1 && 0 == target_dim && source_dim != target_dim )
    {

        // element of mConnectivityMap should be for this type and for one
        // less than source_dim, which should give the connectivity of that sub element
        const ConnMap& cm = mConnectivityMap[this_type][source_dim - 1];
        std::copy( cm.conn[source_indices[0]],
                   cm.conn[source_indices[0]] + cm.num_corners_per_sub_element[source_indices[0]],
                   std::back_inserter( index_list ) );
        return 0;
    }

    // now go through source indices, folding adjacencies into target list
    for( it1 = source_indices; it1 != source_indices + num_source_indices; it1++ )
    {
        // *it1 is the side index
        // at start of iteration, index_list has the target list

        // if a union, or first iteration and index list was empty, copy the list
        if( operation_type == CN::UNION || ( it1 == source_indices && index_list.empty() ) )
        {
            std::copy( MUC.targets_per_source_element[*it1],
                       MUC.targets_per_source_element[*it1] + MUC.num_targets_per_source_element[*it1],
                       std::back_inserter( index_list ) );
        }
        else
        {
            // else we're intersecting, and have a non-empty list; intersect with this target list
            tmp_indices.clear();
            for( int i = MUC.num_targets_per_source_element[*it1] - 1; i >= 0; i-- )
                if( std::find( index_list.begin(), index_list.end(), MUC.targets_per_source_element[*it1][i] ) !=
                    index_list.end() )
                    tmp_indices.push_back( MUC.targets_per_source_element[*it1][i] );
            //      std::set_intersection(MUC.targets_per_source_element[*it1],
            //                            MUC.targets_per_source_element[*it1]+
            //                            MUC.num_targets_per_source_element[*it1],
            //                            index_list.begin(), index_list.end(),
            //                            std::back_inserter(tmp_indices));
            index_list.swap( tmp_indices );

            // if we're at this point and the list is empty, the intersection will be NULL;
            // return if so
            if( index_list.empty() ) return 0;
        }
    }

    if( operation_type == CN::UNION && num_source_indices != 1 )
    {
        // need to sort then unique the list
        std::sort( index_list.begin(), index_list.end() );
        index_list.erase( std::unique( index_list.begin(), index_list.end() ), index_list.end() );
    }

    return 0;
}

template < typename T >
static short int side_number( const T* parent_conn,
                              const EntityType parent_type,
                              const T* child_conn,
                              const int child_num_verts,
                              const int child_dim,
                              int& side_no,
                              int& sense,
                              int& offset )
{
    int parent_num_verts = CN::VerticesPerEntity( parent_type );
    int side_indices[8];
    assert( sizeof( side_indices ) / sizeof( side_indices[0] ) >= (size_t)child_num_verts );

    for( int i = 0; i < child_num_verts; i++ )
    {
        side_indices[i] = std::find( parent_conn, parent_conn + parent_num_verts, child_conn[i] ) - parent_conn;
        if( side_indices[i] == parent_num_verts ) return -1;
    }

    return CN::SideNumber( parent_type, &side_indices[0], child_num_verts, child_dim, side_no, sense, offset );
}

short int CN::SideNumber( const EntityType parent_type,
                          const int* parent_conn,
                          const int* child_conn,
                          const int child_num_verts,
                          const int child_dim,
                          int& side_no,
                          int& sense,
                          int& offset )
{
    return side_number( parent_conn, parent_type, child_conn, child_num_verts, child_dim, side_no, sense, offset );
}

short int CN::SideNumber( const EntityType parent_type,
                          const unsigned int* parent_conn,
                          const unsigned int* child_conn,
                          const int child_num_verts,
                          const int child_dim,
                          int& side_no,
                          int& sense,
                          int& offset )
{
    return side_number( parent_conn, parent_type, child_conn, child_num_verts, child_dim, side_no, sense, offset );
}
short int CN::SideNumber( const EntityType parent_type,
                          const long* parent_conn,
                          const long* child_conn,
                          const int child_num_verts,
                          const int child_dim,
                          int& side_no,
                          int& sense,
                          int& offset )
{
    return side_number( parent_conn, parent_type, child_conn, child_num_verts, child_dim, side_no, sense, offset );
}
short int CN::SideNumber( const EntityType parent_type,
                          const unsigned long* parent_conn,
                          const unsigned long* child_conn,
                          const int child_num_verts,
                          const int child_dim,
                          int& side_no,
                          int& sense,
                          int& offset )
{
    return side_number( parent_conn, parent_type, child_conn, child_num_verts, child_dim, side_no, sense, offset );
}

short int CN::SideNumber( const EntityType parent_type,
                          const unsigned long long* parent_conn,
                          const unsigned long long* child_conn,
                          const int child_num_verts,
                          const int child_dim,
                          int& side_no,
                          int& sense,
                          int& offset )

{
    return side_number( parent_conn, parent_type, child_conn, child_num_verts, child_dim, side_no, sense, offset );
}

short int CN::SideNumber( const EntityType parent_type,
                          void* const* parent_conn,
                          void* const* child_conn,
                          const int child_num_verts,
                          const int child_dim,
                          int& side_no,
                          int& sense,
                          int& offset )
{
    return side_number( parent_conn, parent_type, child_conn, child_num_verts, child_dim, side_no, sense, offset );
}

short int CN::SideNumber( const EntityType parent_type,
                          const int* child_conn_indices,
                          const int child_num_verts,
                          const int child_dim,
                          int& side_no,
                          int& sense,
                          int& offset )
{
    int parent_dim       = Dimension( parent_type );
    int parent_num_verts = VerticesPerEntity( parent_type );

    // degenerate case (vertex), output == input
    if( child_dim == 0 )
    {
        if( child_num_verts != 1 ) return -1;
        side_no = *child_conn_indices;
        sense = offset = 0;
    }

    // given a parent and child element, find the corresponding side number

    // dim_diff should be -1, 0 or 1 (same dimension, one less dimension, two less, resp.)
    if( child_dim > parent_dim || child_dim < 0 ) return -1;

    // different types of same dimension won't be the same
    if( parent_dim == child_dim && parent_num_verts != child_num_verts )
    {
        side_no = -1;
        sense   = 0;
        return 0;
    }

    // loop over the sub-elements, comparing to child connectivity
    int sub_conn_indices[10];
    for( int i = 0; i < NumSubEntities( parent_type, child_dim ); i++ )
    {
        int sub_size = VerticesPerEntity( SubEntityType( parent_type, child_dim, i ) );
        if( sub_size != child_num_verts ) continue;

        SubEntityVertexIndices( parent_type, child_dim, i, sub_conn_indices );
        bool they_match = ConnectivityMatch( child_conn_indices, sub_conn_indices, sub_size, sense, offset );
        if( they_match )
        {
            side_no = i;
            return 0;
        }
    }

    // if we've gotten here, we don't match
    side_no = -1;

    // return value is no success
    return 1;
}

//! return the dimension and index of the opposite side, given parent entity type and child
//! dimension and index.  This function is only defined for certain types of parent/child types:
//! (Parent, Child dim->Opposite dim):
//!  (Tri, 1->0), (Tri, 0->1), (Quad, 1->1), (Quad, 0->0),
//!  (Tet, 2->0), (Tet, 1->1), (Tet, 0->2),
//!  (Hex, 2->2), (Hex, 1->1)(diagonally across element), (Hex, 0->0) (diagonally across element)
//! All other parent types and child dimensions return an error.
//!
//! \param parent_type The type of parent element
//! \param child_type The type of child element
//! \param child_index The index of the child element
//! \param opposite_index The index of the opposite element
//! \return status Returns 0 if successful, -1 if not
short int CN::OppositeSide( const EntityType parent_type,
                            const int child_index,
                            const int child_dim,
                            int& opposite_index,
                            int& opposite_dim )
{
    switch( parent_type )
    {
        case MBEDGE:
            if( 0 != child_dim )
                return -1;
            else
                opposite_index = 1 - child_index;
            opposite_dim = 0;
            break;

        case MBTRI:
            switch( child_dim )
            {
                case 0:
                    opposite_dim   = 1;
                    opposite_index = ( child_index + 1 ) % 3;
                    break;
                case 1:
                    opposite_dim   = 0;
                    opposite_index = ( child_index + 2 ) % 3;
                    break;
                default:
                    return -1;
            }
            break;

        case MBQUAD:
            switch( child_dim )
            {
                case 0:
                case 1:
                    opposite_dim   = child_dim;
                    opposite_index = ( child_index + 2 ) % 4;
                    break;
                default:
                    return -1;
            }
            break;

        case MBTET:
            switch( child_dim )
            {
                case 0:
                    opposite_dim   = 2;
                    opposite_index = ( child_index + 1 ) % 3 + 2 * ( child_index / 3 );
                    break;
                case 1:
                    opposite_dim   = 1;
                    opposite_index = child_index < 3 ? 3 + ( child_index + 2 ) % 3 : ( child_index + 1 ) % 3;
                    break;
                case 2:
                    opposite_dim   = 0;
                    opposite_index = ( child_index + 2 ) % 3 + child_index / 3;
                    break;
                default:
                    return -1;
            }
            break;
        case MBHEX:
            opposite_dim = child_dim;
            switch( child_dim )
            {
                case 0:
                    opposite_index = child_index < 4 ? 4 + ( child_index + 2 ) % 4 : ( child_index - 2 ) % 4;
                    break;
                case 1:
                    opposite_index = 4 * ( 2 - child_index / 4 ) + ( child_index + 2 ) % 4;
                    break;
                case 2:
                    opposite_index = child_index < 4 ? ( child_index + 2 ) % 4 : 9 - child_index;
                    break;
                default:
                    return -1;
            }
            break;

        default:
            return -1;
    }

    return 0;
}

template < typename T >
inline bool connectivity_match( const T* conn1_i, const T* conn2_i, const int num_vertices, int& direct, int& offset )
{

    bool they_match;

    // special test for 2 handles, since we don't want to wrap the list in this
    // case
    if( num_vertices == 2 )
    {
        they_match = false;
        if( conn1_i[0] == conn2_i[0] && conn1_i[1] == conn2_i[1] )
        {
            direct     = 1;
            they_match = true;
            offset     = 0;
        }
        else if( conn1_i[0] == conn2_i[1] && conn1_i[1] == conn2_i[0] )
        {
            they_match = true;
            direct     = -1;
            offset     = 1;
        }
    }

    else
    {
        const T* iter;
        iter = std::find( &conn2_i[0], &conn2_i[num_vertices], conn1_i[0] );
        if( iter == &conn2_i[num_vertices] ) return false;

        they_match = true;

        offset = iter - conn2_i;
        int i;

        // first compare forward
        for( i = 1; i < num_vertices; ++i )
        {
            if( conn1_i[i] != conn2_i[( offset + i ) % num_vertices] )
            {
                they_match = false;
                break;
            }
        }

        if( they_match == true )
        {
            direct = 1;
            return they_match;
        }

        they_match = true;

        // then compare reverse
        for( i = 1; i < num_vertices; i++ )
        {
            if( conn1_i[i] != conn2_i[( offset + num_vertices - i ) % num_vertices] )
            {
                they_match = false;
                break;
            }
        }
        if( they_match )
        {
            direct = -1;
        }
    }

    return they_match;
}

bool CN::ConnectivityMatch( const int* conn1_i, const int* conn2_i, const int num_vertices, int& direct, int& offset )
{
    return connectivity_match< int >( conn1_i, conn2_i, num_vertices, direct, offset );
}

bool CN::ConnectivityMatch( const unsigned int* conn1_i,
                            const unsigned int* conn2_i,
                            const int num_vertices,
                            int& direct,
                            int& offset )
{
    return connectivity_match< unsigned int >( conn1_i, conn2_i, num_vertices, direct, offset );
}

bool CN::ConnectivityMatch( const long* conn1_i, const long* conn2_i, const int num_vertices, int& direct, int& offset )
{
    return connectivity_match< long >( conn1_i, conn2_i, num_vertices, direct, offset );
}

bool CN::ConnectivityMatch( const unsigned long* conn1_i,
                            const unsigned long* conn2_i,
                            const int num_vertices,
                            int& direct,
                            int& offset )
{
    return connectivity_match< unsigned long >( conn1_i, conn2_i, num_vertices, direct, offset );
}

bool CN::ConnectivityMatch( const unsigned long long* conn1_i,
                            const unsigned long long* conn2_i,
                            const int num_vertices,
                            int& direct,
                            int& offset )
{
    return connectivity_match< unsigned long long >( conn1_i, conn2_i, num_vertices, direct, offset );
}

bool CN::ConnectivityMatch( void* const* conn1_i,
                            void* const* conn2_i,
                            const int num_vertices,
                            int& direct,
                            int& offset )
{
    return connectivity_match< void* >( conn1_i, conn2_i, num_vertices, direct, offset );
}

//! for an entity of this type and a specified subfacet (dimension and index), return
//! the index of the higher order node for that entity in this entity's connectivity array
short int CN::HONodeIndex( const EntityType this_type,
                           const int num_verts,
                           const int subfacet_dim,
                           const int subfacet_index )
{
    int i;
    int has_mids[4];
    HasMidNodes( this_type, num_verts, has_mids );

    // if we have no mid nodes on the subfacet_dim, we have no index
    if( subfacet_index != -1 && !has_mids[subfacet_dim] ) return -1;

    // put start index at last index (one less than the number of vertices
    // plus the index basis)
    int index = VerticesPerEntity( this_type ) - 1 + numberBasis;

    // for each subfacet dimension less than the target subfacet dim which has mid nodes,
    // add the number of subfacets of that dimension to the index
    for( i = 1; i < subfacet_dim; i++ )
        if( has_mids[i] ) index += NumSubEntities( this_type, i );

    // now add the index of this subfacet, or one if we're asking about the entity as a whole
    if( subfacet_index == -1 && has_mids[subfacet_dim] )
        // want the index of the last ho node on this subfacet
        index += NumSubEntities( this_type, subfacet_dim );

    else if( subfacet_index != -1 && has_mids[subfacet_dim] )
        index += subfacet_index + 1 - numberBasis;

    // that's it
    return index;
}

//! given data about an element and a vertex in that element, return the dimension
//! and index of the sub-entity that the vertex resolves.  If it does not resolve a
//! sub-entity, either because it's a corner node or it's not in the element, -1 is
//! returned in both return values
void CN::HONodeParent( EntityType elem_type, int num_verts, int ho_index, int& parent_dim, int& parent_index )
{
    // begin with error values
    parent_dim = parent_index = -1;

    // given the number of verts and the element type, get the hasmidnodes solution
    int has_mids[4];
    HasMidNodes( elem_type, num_verts, has_mids );

    int index     = VerticesPerEntity( elem_type ) - 1;
    const int dim = Dimension( elem_type );

    // keep a running sum of the ho node indices for this type of element, and stop
    // when you get to the dimension which has the ho node
    for( int i = 1; i < dim; i++ )
    {
        if( has_mids[i] )
        {
            if( ho_index <= index + NumSubEntities( elem_type, i ) )
            {
                // the ho_index resolves an entity of dimension i, so set the return values
                // and break out of the loop
                parent_dim   = i;
                parent_index = ho_index - index - 1;
                return;
            }
            else
            {
                index += NumSubEntities( elem_type, i );
            }
        }
    }

    // mid region node case
    if( has_mids[dim] && ho_index == index + 1 )
    {
        parent_dim   = dim;
        parent_index = 0;
    }
}

const char* CN::EntityTypeName( const EntityType this_type )
{
    return entityTypeNames[this_type];
}

}  // namespace moab

using moab::CN;
using moab::EntityType;

//! get the basis of the numbering system
void MBCN_GetBasis( int* rval )
{
    *rval = CN::GetBasis();
}

//! set the basis of the numbering system
void MBCN_SetBasis( const int in_basis )
{
    CN::SetBasis( in_basis );
}

//! return the string type name for this type
void MBCN_EntityTypeName( const int this_type, char* rval, int rval_len )
{
    const char* rval_tmp = CN::EntityTypeName( (EntityType)this_type );
    int rval_len_tmp     = strlen( rval_tmp );
    rval_len_tmp         = ( rval_len_tmp < rval_len ? rval_len_tmp : rval_len );
    strncpy( rval, rval_tmp, rval_len_tmp );
}

//! given a name, find the corresponding entity type
void MBCN_EntityTypeFromName( const char* name, int* rval )
{
    *rval = CN::EntityTypeFromName( name );
}

//! return the topological entity dimension
void MBCN_Dimension( const int t, int* rval )
{
    *rval = CN::Dimension( (EntityType)t );
}

//! return the number of (corner) vertices contained in the specified type.
void MBCN_VerticesPerEntity( const int t, int* rval )
{
    *rval = CN::VerticesPerEntity( (EntityType)t );
}

//! return the number of sub-entities bounding the entity.
void MBCN_NumSubEntities( const int t, const int d, int* rval )
{
    *rval = CN::NumSubEntities( (EntityType)t, d );
}

//! return the type of a particular sub-entity.
//! \param this_type Type of entity for which sub-entity type is being queried
//! \param sub_dimension Topological dimension of sub-entity whose type is being queried
//! \param index Index of sub-entity whose type is being queried
//! \return type Entity type of sub-entity with specified dimension and index
void MBCN_SubEntityType( const int this_type, const int sub_dimension, const int index, int* rval )

{

    *rval = CN::SubEntityType( (EntityType)this_type, sub_dimension, index );
}

//! return the vertex indices of the specified sub-entity.
//! \param this_type Type of entity for which sub-entity connectivity is being queried
//! \param sub_dimension Dimension of sub-entity
//! \param sub_index Index of sub-entity
//! \param sub_entity_conn Connectivity of sub-entity (returned to calling function)
void MBCN_SubEntityVertexIndices( const int this_type,
                                  const int sub_dimension,
                                  const int sub_index,
                                  int sub_entity_conn[] )
{
    CN::SubEntityVertexIndices( (EntityType)this_type, sub_dimension, sub_index, sub_entity_conn );
}

//! return the vertices of the specified sub entity
//! \param parent_conn Connectivity of parent entity
//! \param parent_type Entity type of parent entity
//! \param sub_dimension Dimension of sub-entity being queried
//! \param sub_index Index of sub-entity being queried
//! \param sub_entity_conn Connectivity of sub-entity, based on parent_conn and canonical
//!           ordering for parent_type
//! \param num_sub_vertices Number of vertices in sub-entity
//  void MBCN_SubEntityConn(const void *parent_conn, const int parent_type,
//                            const int sub_dimension,
//                            const int sub_index,
//                            void *sub_entity_conn, int &num_sub_vertices) {return
//                            CN::SubEntityConn();}

//! For a specified set of sides of given dimension, return the intersection
//! or union of all sides of specified target dimension adjacent to those sides.
//! \param this_type Type of entity for which sub-entity connectivity is being queried
//! \param source_indices Indices of sides being queried
//! \param num_source_indices Number of entries in <em>source_indices</em>
//! \param source_dim Dimension of source entity
//! \param target_dim Dimension of target entity
//! \param index_list Indices of target entities (returned)
//! \param num_indices Number of indices of target entities (returned)
//! \param operation_type Specify either CN::INTERSECT (0) or CN::UNION (1) to get intersection
//!        or union of target entity lists over source entities
//! \param rval Error code indicating success or failure (returned)
void MBCN_AdjacentSubEntities( const int this_type,
                               const int* source_indices,
                               const int num_source_indices,
                               const int source_dim,
                               const int target_dim,
                               int* index_list,
                               int* num_indices,
                               const int operation_type,
                               int* rval )
{
    std::vector< int > tmp_index_list;
    *rval = CN::AdjacentSubEntities( (EntityType)this_type, source_indices, num_source_indices, source_dim, target_dim,
                                     tmp_index_list, operation_type );
    std::copy( tmp_index_list.begin(), tmp_index_list.end(), index_list );
    *num_indices = tmp_index_list.size();
}

//! return the side index represented in the input sub-entity connectivity
//! \param parent_type Entity type of parent entity
//! \param child_conn_indices Child connectivity to query, specified as indices
//!                           into the connectivity list of the parent.
//! \param child_num_verts Number of values in <em>child_conn_indices</em>
//! \param child_dim Dimension of child entity being queried
//! \param side_no Side number of child entity (returned)
//! \param sense Sense of child entity with respect to order in <em>child_conn</em> (returned)
//! \param offset Offset of <em>child_conn</em> with respect to canonical ordering data (returned)
//! \return status Returns zero if successful, -1 if not
void MBCN_SideNumber( const int parent_type,
                      const int* child_conn_indices,
                      const int child_num_verts,
                      const int child_dim,
                      int* side_no,
                      int* sense,
                      int* offset )
{
    CN::SideNumber( (EntityType)parent_type, child_conn_indices, child_num_verts, child_dim, *side_no, *sense,
                    *offset );
}

void MBCN_SideNumberInt( const int* parent_conn,
                         const EntityType parent_type,
                         const int* child_conn,
                         const int child_num_verts,
                         const int child_dim,
                         int* side_no,
                         int* sense,
                         int* offset )
{
    moab::side_number( parent_conn, parent_type, child_conn, child_num_verts, child_dim, *side_no, *sense, *offset );
}

void MBCN_SideNumberUint( const unsigned int* parent_conn,
                          const EntityType parent_type,
                          const unsigned int* child_conn,
                          const int child_num_verts,
                          const int child_dim,
                          int* side_no,
                          int* sense,
                          int* offset )
{
    moab::side_number( parent_conn, parent_type, child_conn, child_num_verts, child_dim, *side_no, *sense, *offset );
}

void MBCN_SideNumberLong( const long* parent_conn,
                          const EntityType parent_type,
                          const long* child_conn,
                          const int child_num_verts,
                          const int child_dim,
                          int* side_no,
                          int* sense,
                          int* offset )
{
    moab::side_number( parent_conn, parent_type, child_conn, child_num_verts, child_dim, *side_no, *sense, *offset );
}

void MBCN_SideNumberUlong( const unsigned long* parent_conn,
                           const EntityType parent_type,
                           const unsigned long* child_conn,
                           const int child_num_verts,
                           const int child_dim,
                           int* side_no,
                           int* sense,
                           int* offset )
{
    moab::side_number( parent_conn, parent_type, child_conn, child_num_verts, child_dim, *side_no, *sense, *offset );
}

void MBCN_SideNumberVoid( void* const* parent_conn,
                          const EntityType parent_type,
                          void* const* child_conn,
                          const int child_num_verts,
                          const int child_dim,
                          int* side_no,
                          int* sense,
                          int* offset )
{
    moab::side_number( parent_conn, parent_type, child_conn, child_num_verts, child_dim, *side_no, *sense, *offset );
}

//! return the dimension and index of the opposite side, given parent entity type and child
//! dimension and index.  This function is only defined for certain types of parent/child types:
//! (Parent, Child dim->Opposite dim):
//!  (Tri, 1->0), (Tri, 0->1), (Quad, 1->1), (Quad, 0->0),
//!  (Tet, 2->0), (Tet, 1->1), (Tet, 0->2),
//!  (Hex, 2->2), (Hex, 1->1)(diagonally across element), (Hex, 0->0) (diagonally across element)
//! All other parent types and child dimensions return an error.
//!
//! \param parent_type The type of parent element
//! \param child_type The type of child element
//! \param child_index The index of the child element
//! \param opposite_index The index of the opposite element
//! \return status Returns 0 if successful, -1 if not
void MBCN_OppositeSide( const int parent_type,
                        const int child_index,
                        const int child_dim,
                        int* opposite_index,
                        int* opposite_dim,
                        int* rval )
{
    *rval = CN::OppositeSide( (EntityType)parent_type, child_index, child_dim, *opposite_index, *opposite_dim );
}

//! given two connectivity arrays, determine whether or not they represent the same entity.
//! \param conn1 Connectivity array of first entity
//! \param conn2 Connectivity array of second entity
//! \param num_vertices Number of entries in <em>conn1</em> and <em>conn2</em>
//! \param direct If positive, entities have the same sense (returned)
//! \param offset Offset of <em>conn2</em>'s first vertex in <em>conn1</em>
//! \return rval Returns true if <em>conn1</em> and <em>conn2</em> match
void MBCN_ConnectivityMatchInt( const int* conn1,
                                const int* conn2,
                                const int num_vertices,
                                int* direct,
                                int* offset,
                                int* rval )
{
    *rval = CN::ConnectivityMatch( conn1, conn2, num_vertices, *direct, *offset );
}

void MBCN_ConnectivityMatchUint( const unsigned int* conn1,
                                 const unsigned int* conn2,
                                 const int num_vertices,
                                 int* direct,
                                 int* offset,
                                 int* rval )
{
    *rval = CN::ConnectivityMatch( conn1, conn2, num_vertices, *direct, *offset );
}

void MBCN_ConnectivityMatchLong( const long* conn1,
                                 const long* conn2,
                                 const int num_vertices,
                                 int* direct,
                                 int* offset,
                                 int* rval )
{
    *rval = CN::ConnectivityMatch( conn1, conn2, num_vertices, *direct, *offset );
}

void MBCN_ConnectivityMatchUlong( const unsigned long* conn1,
                                  const unsigned long* conn2,
                                  const int num_vertices,
                                  int* direct,
                                  int* offset,
                                  int* rval )
{
    *rval = CN::ConnectivityMatch( conn1, conn2, num_vertices, *direct, *offset );
}

void MBCN_ConnectivityMatchVoid( void* const* conn1,
                                 void* const* conn2,
                                 const int num_vertices,
                                 int* direct,
                                 int* offset,
                                 int* rval )
{
    *rval = CN::ConnectivityMatch( conn1, conn2, num_vertices, *direct, *offset );
}

//! true if entities of a given type and number of nodes indicates mid edge nodes are present.
//! \param this_type Type of entity for which sub-entity connectivity is being queried
//! \param num_verts Number of nodes defining entity
//! \return int Returns true if <em>this_type</em> combined with <em>num_nodes</em> indicates
//!  mid-edge nodes are likely
void MBCN_HasMidEdgeNodes( const int this_type, const int num_verts, int* rval )
{
    *rval = CN::HasMidEdgeNodes( (EntityType)this_type, num_verts );
}

//! true if entities of a given type and number of nodes indicates mid face nodes are present.
//! \param this_type Type of entity for which sub-entity connectivity is being queried
//! \param num_verts Number of nodes defining entity
//! \return int Returns true if <em>this_type</em> combined with <em>num_nodes</em> indicates
//!  mid-face nodes are likely
void MBCN_HasMidFaceNodes( const int this_type, const int num_verts, int* rval )
{
    *rval = CN::HasMidFaceNodes( (EntityType)this_type, num_verts );
}

//! true if entities of a given type and number of nodes indicates mid region nodes are present.
//! \param this_type Type of entity for which sub-entity connectivity is being queried
//! \param num_verts Number of nodes defining entity
//! \return int Returns true if <em>this_type</em> combined with <em>num_nodes</em> indicates
//!  mid-region nodes are likely
void MBCN_HasMidRegionNodes( const int this_type, const int num_verts, int* rval )
{
    *rval = CN::HasMidRegionNodes( (EntityType)this_type, num_verts );
}

//! true if entities of a given type and number of nodes indicates mid edge/face/region nodes
//! are present.
//! \param this_type Type of entity for which sub-entity connectivity is being queried
//! \param num_verts Number of nodes defining entity
//! \param mid_nodes If <em>mid_nodes[i], i=1..3</em> is true, indicates that mid-edge
//!    (i=1), mid-face (i=2), and/or mid-region (i=3) nodes are likely
void MBCN_HasMidNodes( const int this_type, const int num_verts, int mid_nodes[4] )
{
    return CN::HasMidNodes( (EntityType)this_type, num_verts, mid_nodes );
}

//! given data about an element and a vertex in that element, return the dimension
//! and index of the sub-entity that the vertex resolves.  If it does not resolve a
//! sub-entity, either because it's a corner node or it's not in the element, -1 is
//! returned in both return values.
//! \param elem_type Type of entity being queried
//! \param num_nodes The number of nodes in the element connectivity
//! \param ho_node_index The position of the HO node in the connectivity list (zero based)
//! \param parent_dim Dimension of sub-entity high-order node resolves (returned)
//! \param parent_index Index of sub-entity high-order node resolves (returned)
void MBCN_HONodeParent( int elem_type, int num_nodes, int ho_node_index, int* parent_dim, int* parent_index )
{
    return CN::HONodeParent( (EntityType)elem_type, num_nodes, ho_node_index, *parent_dim, *parent_index );
}

//! for an entity of this type with num_verts vertices, and a specified subfacet
//! (dimension and index), return the index of the higher order node for that entity
//! in this entity's connectivity array
//! \param this_type Type of entity being queried
//! \param num_verts Number of vertices for the entity being queried
//! \param subfacet_dim Dimension of sub-entity being queried
//! \param subfacet_index Index of sub-entity being queried
//! \return index Index of sub-entity's higher-order node
void MBCN_HONodeIndex( const int this_type,
                       const int num_verts,
                       const int subfacet_dim,
                       const int subfacet_index,
                       int* rval )

{

    *rval = CN::HONodeIndex( (EntityType)this_type, num_verts, subfacet_dim, subfacet_index );
}

namespace moab
{

template < typename T >
inline int permute_this( EntityType t, const int dim, T* conn, const int indices_per_ent, const int num_entries )
{
    T tmp_conn[MAX_SUB_ENTITIES];
    assert( indices_per_ent <= CN::permuteVec[t][dim][MAX_SUB_ENTITIES] );
    if( indices_per_ent > CN::permuteVec[t][dim][MAX_SUB_ENTITIES] ) return 1;
    short int* tvec = CN::permuteVec[t][dim];
    T* pvec         = conn;
    for( int j = 0; j < num_entries; j++ )
    {
        for( int i = 0; i < indices_per_ent; i++ )
            tmp_conn[tvec[i]] = pvec[i];
        memcpy( pvec, tmp_conn, indices_per_ent * sizeof( T ) );
        pvec += indices_per_ent;
    }

    return 0;
}

template < typename T >
inline int rev_permute_this( EntityType t, const int dim, T* conn, const int indices_per_ent, const int num_entries )
{
    T tmp_conn[MAX_SUB_ENTITIES];
    assert( indices_per_ent <= CN::revPermuteVec[t][dim][MAX_SUB_ENTITIES] );
    if( indices_per_ent > CN::revPermuteVec[t][dim][MAX_SUB_ENTITIES] ) return 1;
    short int* tvec = CN::revPermuteVec[t][dim];
    T* pvec         = conn;
    for( int j = 0; j < num_entries; j++ )
    {
        for( int i = 0; i < indices_per_ent; i++ )
            tmp_conn[i] = pvec[tvec[i]];
        memcpy( pvec, tmp_conn, indices_per_ent * sizeof( T ) );
        pvec += indices_per_ent;
    }

    return 0;
}

short int CN::Dimension( const EntityType t )
{
    return mConnectivityMap[t][0].topo_dimension;
}

short int CN::VerticesPerEntity( const EntityType t )
{
    return ( MBVERTEX == t
                 ? (short int)1
                 : mConnectivityMap[t][mConnectivityMap[t][0].topo_dimension - 1].num_corners_per_sub_element[0] );
}

short int CN::NumSubEntities( const EntityType t, const int d )
{
    return ( t != MBVERTEX && d > 0 ? mConnectivityMap[t][d - 1].num_sub_elements
                                    : ( d ? (short int)-1 : VerticesPerEntity( t ) ) );
}

//! return the type of a particular sub-entity.
EntityType CN::SubEntityType( const EntityType this_type, const int sub_dimension, const int index )
{
    return ( !sub_dimension ? MBVERTEX
                            : ( Dimension( this_type ) == sub_dimension && 0 == index
                                    ? this_type
                                    : mConnectivityMap[this_type][sub_dimension - 1].target_type[index] ) );
}

const short* CN::SubEntityVertexIndices( const EntityType this_type,
                                         const int sub_dimension,
                                         const int index,
                                         EntityType& sub_type,
                                         int& n )
{
    if( sub_dimension == 0 )
    {
        n        = 1;
        sub_type = MBVERTEX;
        return increasingInts + index;
    }
    else
    {
        const CN::ConnMap& map = mConnectivityMap[this_type][sub_dimension - 1];
        sub_type               = map.target_type[index];
        n                      = map.num_corners_per_sub_element[index];
        return map.conn[index];
    }
}

//! Permute this vector
inline int CN::permuteThis( const EntityType t, const int dim, int* pvec, const int num_indices, const int num_entries )
{
    return permute_this( t, dim, pvec, num_indices, num_entries );
}
inline int CN::permuteThis( const EntityType t,
                            const int dim,
                            unsigned int* pvec,
                            const int num_indices,
                            const int num_entries )
{
    return permute_this( t, dim, pvec, num_indices, num_entries );
}
inline int CN::permuteThis( const EntityType t,
                            const int dim,
                            long* pvec,
                            const int num_indices,
                            const int num_entries )
{
    return permute_this( t, dim, pvec, num_indices, num_entries );
}
inline int CN::permuteThis( const EntityType t,
                            const int dim,
                            void** pvec,
                            const int num_indices,
                            const int num_entries )
{
    return permute_this( t, dim, pvec, num_indices, num_entries );
}

//! Reverse permute this vector
inline int CN::revPermuteThis( const EntityType t,
                               const int dim,
                               int* pvec,
                               const int num_indices,
                               const int num_entries )
{
    return rev_permute_this( t, dim, pvec, num_indices, num_entries );
}
inline int CN::revPermuteThis( const EntityType t,
                               const int dim,
                               unsigned int* pvec,
                               const int num_indices,
                               const int num_entries )
{
    return rev_permute_this( t, dim, pvec, num_indices, num_entries );
}
inline int CN::revPermuteThis( const EntityType t,
                               const int dim,
                               long* pvec,
                               const int num_indices,
                               const int num_entries )
{
    return rev_permute_this( t, dim, pvec, num_indices, num_entries );
}
inline int CN::revPermuteThis( const EntityType t,
                               const int dim,
                               void** pvec,
                               const int num_indices,
                               const int num_entries )
{
    return rev_permute_this( t, dim, pvec, num_indices, num_entries );
}

}  // namespace moab