WriteCGNS.cpp

Go to the documentation of this file.

#include "WriteCGNS.hpp"
#include "moab/CN.hpp"
#include "MBTagConventions.hpp"
#include "MBParallelConventions.h"
#include "moab/Interface.hpp"
#include "moab/Range.hpp"
#include "moab/WriteUtilIface.hpp"
#include "moab/FileOptions.hpp"
#include "GmshUtil.hpp"
 
#include <fstream>
#include <map>
#include <set>
 
#include <iostream>
 
namespace moab
{
 
WriterIface* WriteCGNS::factory( Interface* iface )
{
    return new WriteCGNS( iface );
}
 
WriteCGNS::WriteCGNS( Interface* impl )
    : mbImpl( impl ), fileName( NULL ), IndexFile( 0 ), BaseName( NULL ), IndexBase( 0 ), ZoneName( NULL ),
      IndexZone( 0 ), IndexSection( 0 ), celldim( 0 ), physdim( 0 ), VrtSize( 0 ), EdgeSize( 0 ), FaceSize( 0 ),
      CellSize( 0 )
{
    impl->query_interface( mWriteIface );
    IndexCoord[0] = IndexCoord[1] = IndexCoord[2] = 0;
    isize[0] = isize[1] = isize[2] = 0;
}
 
WriteCGNS::~WriteCGNS()
{
    mbImpl->release_interface( mWriteIface );
}
 
//! writes out a file
ErrorCode WriteCGNS::write_file( const char* file_name,
                                 const bool overwrite,
                                 const FileOptions& /*options*/,
                                 const EntityHandle* /*output_list*/,
                                 const int /*num_sets*/,
                                 const std::vector< std::string >&,
                                 const Tag*,
                                 int,
                                 int )
{
    ErrorCode rval;
 
    if( !overwrite )
    {
        rval = mWriteIface->check_doesnt_exist( file_name );
        if( MB_SUCCESS != rval ) return rval;
    }
    std::cout << "THE CGNS CONVERSION ONLY WORKS FOR ENTITIES CITED BELOW:";
    std::cout << "\n   -MBVERTEX\n   -MBEDGE\n   -MBTRI\n   -MBQUAD";
    std::cout << "\n   -MBTET\n   -MBPYRAMID\n   -MBHEX\n";
 
    // Get entities to write
    // Get and count vertex entities
    rval = get_vertex_entities( VrtSize, Nodes );
    if( rval != MB_SUCCESS )
    {
        return rval;
    }
    // Get and count edge entities
    rval = get_edge_entities( EdgeSize, Edges );
    if( rval != MB_SUCCESS )
    {
        return rval;
    }
    // Get and count face entities
    rval = get_face_entities( FaceSize, Faces );
    if( rval != MB_SUCCESS )
    {
        return rval;
    }
    // Get and count cell entities
    rval = get_cell_entities( CellSize, Cells );
    if( rval != MB_SUCCESS )
    {
        return rval;
    }
    std::cout << "\nThe Number of Vertex is " << VrtSize << ".\n";
    std::cout << "The Number of Edges is " << EdgeSize << ".\n";
    std::cout << "The Number of Faces is " << FaceSize << ".\n";
    std::cout << "The Number of Cells is " << CellSize << ".\n\n";
 
    // save filename to member variable so we don't need to pass as an argument
    // to called functions
    fileName = file_name;
    std::cout << fileName << " file is a " << physdim << "-D mesh.\n";
 
    // Open file
    IndexFile = 0;
 
    // open the cgns file
    // filename:      (input) Name of the CGNS file, including path name if necessary. There is no
    // limit on the
    //                length of this character variable.
    // CG_MODE_WRITE: (input) Mode used for opening the file. The modes currently supported are
    // CG_MODE_READ,
    //                CG_MODE_WRITE, and CG_MODE_MODIFY.
    // filePtr:       (output) CGNS file index number.
    if( cg_open( fileName, CG_MODE_WRITE, &IndexFile ) )
    {
        std::cout << "Error opening file\n";
        cg_error_exit();
    }
    // Give a base name
    BaseName = "Cgns Base";
    if( cg_base_write( IndexFile, BaseName, celldim, physdim, &IndexBase ) )
    {
        std::cout << "Error creating CGNS base";
    }
    // Give a zone name
    ZoneName = "Cgns Zone";
    // isize array contains the total vertex size, cell size, and boundary
    // vertex size for the zone
    // Note that for unstructured zones, the index dimension is always 1
    isize[0] = VrtSize;   // isize[0] contains the total vertex size
    isize[1] = CellSize;  // isize[1] contains the total cell size
    isize[2] = 0;         // isize[2] = 0 for unsorted elements
    // Create zone */
    // ZoneType_t: Unstructured
    if( cg_zone_write( IndexFile, IndexBase, ZoneName, isize, Unstructured, &IndexZone ) )
    {
        std::cout << "Error creating CGNS zone\n";
        cg_error_exit();
    }
    // Write the vertex coordinates
    rval = write_coord_cgns( Nodes );
    if( rval != MB_SUCCESS )
    {
        return rval;
    }
 
    // Create a vector to hold the Tags
    std::vector< moab::Tag > TagHandles;
    // Get Tags
    rval = mbImpl->tag_get_tags( TagHandles );
    if( rval != MB_SUCCESS )
    {
        return rval;
    }
    // Get the number of Tags in the mesh
    int NbTags = TagHandles.size();
    std::cout << "\nThe mesh has " << NbTags << " Tags.\n";
 
    // Create a vector of size NbTags
    // Sets have informations about the entity set
    std::vector< SetStruct > Sets;
    Sets.reserve( NbTags );
    // Fill Sets with all information needed
    rval = set_tag_values( TagHandles, Edges, Faces, Cells, Sets );
    if( rval != MB_SUCCESS )
    {
        std::cout << "Problem to set tag values\n";
        return rval;
    }
 
    // Create a matrix to hold connectivity
    std::vector< std::vector< cgsize_t > > ConnTable;
    ConnTable.resize( NbTags );
 
    std::vector< int > Begin( NbTags, 0 );
    std::vector< int > End( NbTags, 0 );
 
    // Take the connectivity of higher dimension entities
    cgsize_t BeginSetsIndex = 1;
    cgsize_t EndSetsIndex;
    switch( physdim )
    {
        case 1:
            rval = get_conn_table( Edges, Begin, End, TagHandles, Sets, ConnTable );
            if( rval != MB_SUCCESS )
            {
                std::cout << "Problem to fill the connectivity table for 1-D entities\n";
                return rval;
            }
            for( int i = 0; i < NbTags; ++i )
            {
                if( Sets[i].IdSet != -1 )
                {
                    const char* SectionName = Sets[i].TagName.c_str();
                    EndSetsIndex            = BeginSetsIndex + Sets[i].NbEdges - 1;
                    // Write the section in CGNS file
                    if( cg_section_write( IndexFile, IndexBase, IndexBase, SectionName, Sets[i].CGNSType,
                                          BeginSetsIndex, EndSetsIndex, 0, &ConnTable[i][0], &IndexSection ) )
                    {
                        std::cout << "Issue on writing connectivity - 3-D\n";
                        cg_error_exit();
                    }
                    BeginSetsIndex = EndSetsIndex + 1;
                }
            }
            break;
        case 2:
            rval = get_conn_table( Edges, Begin, End, TagHandles, Sets, ConnTable );
            if( rval != MB_SUCCESS )
            {
                std::cout << "Problem to fill the connectivity table for 1-D entities\n";
                return rval;
            }
            rval = get_conn_table( Faces, Begin, End, TagHandles, Sets, ConnTable );
            if( rval != MB_SUCCESS )
            {
                std::cout << "Problem to fill the connectivity table for 2-D entities\n";
                return rval;
            }
            for( int i = 0; i < NbTags; ++i )
            {
                if( Sets[i].IdSet != -1 )
                {
                    const char* SectionName = Sets[i].TagName.c_str();
                    EndSetsIndex            = BeginSetsIndex + Sets[i].NbEdges + Sets[i].NbFaces - 1;
                    // Write the section in CGNS file
                    if( cg_section_write( IndexFile, IndexBase, IndexBase, SectionName, Sets[i].CGNSType,
                                          BeginSetsIndex, EndSetsIndex, 0, &ConnTable[i][0], &IndexSection ) )
                    {
                        std::cout << "Issue on writing connectivity -- 2-D\n";
                        cg_error_exit();
                    }
                    BeginSetsIndex = EndSetsIndex + 1;
                }
            }
            break;
        case 3:
            rval = get_conn_table( Faces, Begin, End, TagHandles, Sets, ConnTable );
            if( rval != MB_SUCCESS )
            {
                std::cout << "Problem to fill the connectivity table for 2-D entities\n";
                return rval;
            }
            rval = get_conn_table( Cells, Begin, End, TagHandles, Sets, ConnTable );
            if( rval != MB_SUCCESS )
            {
                std::cout << "Problem to fill the connectivity table for 3-D entities\n";
                return rval;
            }
            for( int i = 0; i < NbTags; ++i )
            {
                if( Sets[i].IdSet != -1 )
                {
                    const char* SectionName = Sets[i].TagName.c_str();
                    EndSetsIndex            = BeginSetsIndex + Sets[i].NbFaces + Sets[i].NbCells - 1;
                    std::cout << "BeginSetsIndex = " << BeginSetsIndex << "\tEndSetsIndex = " << EndSetsIndex << "\n";
                    // Write the section in CGNS file
                    if( cg_section_write( IndexFile, IndexBase, IndexBase, SectionName, Sets[i].CGNSType,
                                          BeginSetsIndex, EndSetsIndex, 0, &ConnTable[i][0], &IndexSection ) )
                    {
                        std::cout << "Issue on writing connectivity -- 3-D\n";
                        cg_error_exit();
                    }
                    BeginSetsIndex = EndSetsIndex + 1;
                }
            }
            break;
        default:
            std::cout << "Issue on Physical dimension\n";
            return MB_FAILURE;
    }
 
    // Close the CGNS mesh file
    if( cg_close( IndexFile ) )
    {
        std::cout << "Error closing file\n";
        cg_error_exit();
    }
    // done
    return MB_SUCCESS;
}
 
// Get and count vertex entities
ErrorCode WriteCGNS::get_vertex_entities( cgsize_t& VrtSize_, std::vector< moab::EntityHandle >& Nodes_ )
{
    ErrorCode rval;
    // Get vertex entities
    // The first input is 0 because one queries the entire mesh.
    // Retrieves all entities of dimension = 0 in "Nodes"
    rval = mbImpl->get_entities_by_dimension( 0, 0, Nodes_, false );
    if( Nodes.size() > 0 )
    {
        celldim = 0;
        physdim = 0;
        // get the amout of vertex
        VrtSize_ = Nodes_.size();
    }
    else
    {
        std::cout << "The mesh has not node points.\n";
    }
    // done
    return rval;
}
 
// Get and count edge entities
ErrorCode WriteCGNS::get_edge_entities( cgsize_t& EdgeSize_, std::vector< moab::EntityHandle >& Edges_ )
{
    ErrorCode rval;
    // The first input is 0 because one queries the entire mesh.
    // Get all entities of dimension = 1 in Edges
    rval = mbImpl->get_entities_by_dimension( 0, 1, Edges_, false );
    if( Edges_.size() > 0 )
    {
        celldim = 1;
        physdim = 1;
        // get the amout of edges
        EdgeSize_ = Edges_.size();
    }
    // done
    return rval;
}
 
// Get and count face entities
ErrorCode WriteCGNS::get_face_entities( cgsize_t& FaceSize_, std::vector< moab::EntityHandle >& Faces_ )
{
    ErrorCode rval;
    // The first input is 0 because one queries the entire mesh.
    // Get all entities of dimension = 2 in Faces
    rval = mbImpl->get_entities_by_dimension( 0, 2, Faces_, false );
    if( Faces_.size() )
    {
        celldim = 2;
        physdim = 2;
        // get the amout of faces
        FaceSize_ = Faces_.size();
    }
    // done
    return rval;
}
 
// Get and count cell entities
ErrorCode WriteCGNS::get_cell_entities( cgsize_t& CellSize_, std::vector< moab::EntityHandle >& Cells_ )
{
    ErrorCode rval;
    // The first input is 0 because one queries the entire mesh.
    // Get all entities of dimension = 3 in Cell
    rval = mbImpl->get_entities_by_dimension( 0, 3, Cells_, false );
    if( Cells_.size() )
    {
        celldim = 3;
        physdim = 3;
        // get the amout of volumes
        CellSize_ = Cells_.size();
    }
    // done
    return rval;
}
 
ErrorCode WriteCGNS::write_coord_cgns( std::vector< moab::EntityHandle >& Nodes_ )
{
    ErrorCode rval;
 
    const int num_entities = (int)Nodes_.size();
 
    // Moab works with one vector for the threee coordinates
    std::vector< double > Coords( 3 * num_entities );
    std::vector< double >::iterator c = Coords.begin();
 
    // CGNS uses one vector for each coordinate
    std::vector< double > CoordX;
    std::vector< double > CoordY;
    std::vector< double > CoordZ;
 
    // Summ the values of all coordinates to be sure if it is not zero
    double SumX = 0;
    double SumY = 0;
    double SumZ = 0;
 
    // Get the moab coordinates - Coords is the output
    rval = mbImpl->get_coords( &Nodes_[0], num_entities, &Coords[0] );
    if( MB_SUCCESS != rval )
    {
        std::cout << "Error getting coordinates from nodes.\n";
        return rval;
    }
 
    // Reserve the size of nodes
    CoordX.reserve( Nodes_.size() );
    CoordY.reserve( Nodes_.size() );
    CoordZ.reserve( Nodes_.size() );
 
    for( std::vector< moab::EntityHandle >::iterator i = Nodes_.begin(); i != Nodes_.end(); ++i )
    {
        CoordX.push_back( *c );  // Put the X coordinate in CoordX vector
        SumX += abs( *c );       // Sum all X coordinates
        ++c;                     // Move to Y coordinate
        CoordY.push_back( *c );  // Put the Y coordinate in CoordY vector
        SumY += abs( *c );       // Sum all Y coordinates
        ++c;                     // Move to Z coordinate
        CoordZ.push_back( *c );  // Put the Z coordinate in CoordZ vector
        SumZ += abs( *c );       // Sum all Z coordinates
        ++c;                     // Move to X coordinate
    }
 
    // If X coordinate is not empty then write CoordX (user must use SIDS-standard names here)
    if( SumX != 0 )
    {
        if( cg_coord_write( IndexFile, IndexBase, IndexZone, RealDouble, "CoordinateX", &CoordX[0], &IndexCoord[0] ) )
        {
            std::cout << " Error writing X coordinates.\n";
            cg_error_exit();
        }
    }
    // If Y coordinate is not empty then write CoordY (user must use SIDS-standard names here)
    if( SumY != 0 )
    {
        if( cg_coord_write( IndexFile, IndexBase, IndexZone, RealDouble, "CoordinateY", &CoordY[0], &IndexCoord[1] ) )
        {
            std::cout << " Error writing Y coordinates.\n";
            cg_error_exit();
        }
    }
    // If Z coordinate is not empty then write CoordZ (user must use SIDS-standard names here)
    if( SumZ != 0 )
    {
        if( cg_coord_write( IndexFile, IndexBase, IndexZone, RealDouble, "CoordinateZ", &CoordZ[0], &IndexCoord[2] ) )
        {
            std::cout << " Error writing Z coordinates.\n";
            cg_error_exit();
        }
    }
 
    // Clear vectors
    Coords.clear();
    CoordX.clear();
    CoordY.clear();
    CoordZ.clear();
 
    // done
    return MB_SUCCESS;
}
 
ErrorCode WriteCGNS::set_tag_values( std::vector< Tag >& TagHandles,
                                     std::vector< moab::EntityHandle >& Edges_,
                                     std::vector< moab::EntityHandle >& Faces_,
                                     std::vector< moab::EntityHandle >& Cells_,
                                     std::vector< WriteCGNS::SetStruct >& Sets )
{
    ErrorCode rval;
 
    // Get the number of Tags in the mesh
    int NbTags = TagHandles.size();
 
    // Loop over all Tags
    for( int i = 0; i < NbTags; ++i )
    {
 
        // Allocate another position in the vector of SetStruct using a default constructor
        Sets.push_back( SetStruct() );
 
        // Get the Tag name
        rval = mbImpl->tag_get_name( TagHandles[i], Sets[i].TagName );
        if( rval != MB_SUCCESS )
        {
            std::cout << "Problem to get Tag Name\n";
            return rval;
        }
        std::cout << "Tag name= " << Sets[i].TagName << "\n";
 
        // Count all entities by type and put in Sets[i].NbEntities vector
        rval = get_set_entities( i, TagHandles, Sets );
        if( rval != MB_SUCCESS )
        {
            std::cout << "Problem to get Set entities\n";
            return rval;
        }
 
        // Get the CGNSTYpe of the Set
        rval = get_cgns_type( i, Sets );
        if( rval != MB_SUCCESS )
        {
            std::cout << "Problem to get CGNSType\n";
            return rval;
        }
        std::cout << "\tSets[" << i << "].CGNSType= " << Sets[i].CGNSType << "\n";
 
        // int Number;
 
        // Set a data index for Edges and TagHandles[i]
        if( Sets[i].CGNSType == BAR_2 || Sets[i].CGNSType == TRI_3 || Sets[i].CGNSType == QUAD_4 ||
            Sets[i].CGNSType == TETRA_4 || Sets[i].CGNSType == PYRA_5 || Sets[i].CGNSType == PENTA_6 ||
            Sets[i].CGNSType == HEXA_8 || Sets[i].CGNSType == MIXED )
        {
 
            if( Sets[i].NbEdges > 0 && physdim < 3 )
            {
                // Set a data index for Edges and TagHandles[i]
                const std::vector< int > tag_values( Edges_.size(), i );
                rval = mbImpl->tag_set_data( TagHandles[i], &Edges_[0], Edges_.size(), &tag_values[0] );
                if( rval != MB_SUCCESS )
                {
                    std::cout << "Problem to set data for 1-D entities\n";
                    return rval;
                }
            }
            if( Sets[i].NbFaces > 0 && physdim > 1 )
            {
                // Set a data index for Faces and TagHandles[i]
                const std::vector< int > tag_values( Faces.size(), i );
                rval = mbImpl->tag_set_data( TagHandles[i], &Faces_[0], Faces_.size(), &tag_values[0] );
                if( rval != MB_SUCCESS )
                {
                    std::cout << "Problem to set data for 2-D entities\n";
                    return rval;
                }
            }
            if( Sets[i].NbCells > 0 && physdim > 2 )
            {
                // Set a data index for Cells and TagHandles[i]
                const std::vector< int > tag_values( Cells.size(), i );
                rval = mbImpl->tag_set_data( TagHandles[i], &Cells_[0], Cells_.size(), &tag_values[0] );
                if( rval != MB_SUCCESS )
                {
                    std::cout << "Problem to set data for 3-D entities\n";
                    return rval;
                }
            }
            // IdSet gets the Set Index indicating that it is not empty
            Sets[i].IdSet = i;
        }
        std::cout << "\tSets[" << i << "].IdSet = " << Sets[i].IdSet << "\n\n";
    }
    return MB_SUCCESS;
}
 
// Get Entities in the set
ErrorCode WriteCGNS::get_set_entities( int i,
                                       std::vector< Tag >& TagHandles,
                                       std::vector< WriteCGNS::SetStruct >& Sets )
{
    ErrorCode rval;
 
    // Get the number of MBEDGE entities
    // NbEntities[0] holds the number of MBEDGE in the "Sets"
    int Number = 0;
    rval       = mbImpl->get_number_entities_by_type_and_tag( 0, MBEDGE, &TagHandles[i], 0, 1, Number );
    if( rval != MB_SUCCESS )
    {
        std::cout << "Problem to get the number of entities by type and tag\n";
        return rval;
    }
    Sets[i].NbEntities.push_back( Number );  // MBEDGE == Sets[i].NbEntities[0]
    Sets[i].NbEdges += Number;
    std::cout << "\tNumber of MBEDGE = " << Number << "\n";
 
    // Get the number of MBTRI entities
    // NbEntities[1] holds the number of MBTRI in the "Sets"
    Number = 0;
    rval   = mbImpl->get_number_entities_by_type_and_tag( 0, MBTRI, &TagHandles[i], 0, 1, Number );
    if( rval != MB_SUCCESS )
    {
        std::cout << "Problem to get the number of entities by type and tag\n";
        return rval;
    }
    Sets[i].NbEntities.push_back( Number );  // MBTRI == Sets[i].NbEntities[1]
    Sets[i].NbFaces += Number;
    std::cout << "\tNumber of MBTRI = " << Number << "\n";
 
    // Get the number of MBQUAD entities
    // NbEntities[2] holds the number of MBQUAD in the "Sets"
    Number = 0;
    rval   = mbImpl->get_number_entities_by_type_and_tag( 0, MBQUAD, &TagHandles[i], 0, 1, Number );
    if( rval != MB_SUCCESS )
    {
        std::cout << "Problem to get the number of entities by type and tag\n";
        return rval;
    }
    Sets[i].NbEntities.push_back( Number );  // MBQUAD == Sets[i].NbEntities[2]
    Sets[i].NbFaces += Number;
    std::cout << "\tNumber of MBQUAD = " << Number << "\n";
 
    // Get the number of MBTET entities
    // NbEntities[3] holds the number of MBTET in the "Sets"
    Number = 0;
    rval   = mbImpl->get_number_entities_by_type_and_tag( 0, MBTET, &TagHandles[i], 0, 1, Number );
    if( rval != MB_SUCCESS )
    {
        std::cout << "Problem to get the number of entities by type and tag\n";
        return rval;
    }
    Sets[i].NbEntities.push_back( Number );  // MBTET == Sets[i].NbEntities[3]
    Sets[i].NbCells += Number;
    std::cout << "\tNumber of MBTET = " << Number << "\n";
 
    // Get the number of MBPYRAMID entities
    // NbEntities[4] holds the number of MBPYRAMID in the "Sets"
    Number = 0;
    rval   = mbImpl->get_number_entities_by_type_and_tag( 0, MBPYRAMID, &TagHandles[i], 0, 1, Number );
    if( rval != MB_SUCCESS )
    {
        std::cout << "Problem to get the number of entities by type and tag\n";
        return rval;
    }
    Sets[i].NbEntities.push_back( Number );  // MBPYRAMID == Sets[i].NbEntities[4]
    Sets[i].NbCells += Number;
    std::cout << "\tNumber of MBPYRAMID = " << Number << "\n";
 
    // Get the number of MBPRISM entities - MBPRISM == PENTA_6
    // NbEntities[5] holds the number of MBPRISM in the "Sets"
    Number = 0;
    rval   = mbImpl->get_number_entities_by_type_and_tag( 0, MBPRISM, &TagHandles[i], 0, 1, Number );
    if( rval != MB_SUCCESS )
    {
        std::cout << "Problem to get the number of entities by type and tag\n";
        return rval;
    }
    Sets[i].NbEntities.push_back( Number );  // MBPRISM == Sets[i].NbEntities[5]
    Sets[i].NbCells += Number;
    std::cout << "\tNumber of MBPRISM = " << Number << "\n";
 
    // Get the number of MBHEX entities
    // NbEntities[6] holds the number of MBHEX in the "Sets"
    Number = 0;
    rval   = mbImpl->get_number_entities_by_type_and_tag( 0, MBHEX, &TagHandles[i], 0, 1, Number );
    if( rval != MB_SUCCESS )
    {
        std::cout << "Problem to get the number of entities by type and tag\n";
        return rval;
    }
    Sets[i].NbEntities.push_back( Number );  // MBHEX == Sets[i].NbEntities[6]
    Sets[i].NbCells += Number;
    std::cout << "\tNumber of MBHEX = " << Number << "\n";
 
    std::cout << "\tTotal number of Edges = " << Sets[i].NbEdges << "\n";
    std::cout << "\tTotal number of Faces = " << Sets[i].NbFaces << "\n";
    std::cout << "\tTotal number of Cells = " << Sets[i].NbCells << "\n";
 
    return MB_SUCCESS;
}
 
// Get CGNSType
ErrorCode WriteCGNS::get_cgns_type( int i, std::vector< WriteCGNS::SetStruct >& Sets )
{
    std::vector< int > Test;
    int Sum = 0;
 
    // NbEntities is a vector which has the number of entities of each type
    // 0-MBEDGE | 1-MBTRI | 2-MBQUAD | 3-MBTET | 4-MBPYRAMID | 5-MBPRISM | 6-MBHEX
    // if NbEntities[i]>0 then Test[i]=1
    // else then Test[i]=0
    for( int j = 0; j < (int)Sets[i].NbEntities.size(); ++j )
    {
        if( Sets[i].NbEntities[j] > 0 )
        {
            Test.push_back( 1 );
            Sum++;
        }
        else
        {
            Test.push_back( 0 );
        }
    }
 
    // Test the Sum
    // if Sum > 1 then the Set is MIXED
    // if Sum = 0 then the Set is Homogeneous
    // else then the Set is empty
    if( Sum > 1 )
    {
        Sets[i].CGNSType = MIXED;
    }
    else if( Sum == 1 )
    {
        int j = 0;
        std::cout << "Homogeneous Type\n";
        while( j < (int)Sets[i].NbEntities.size() && Sets[i].NbEntities[j] != 1 )
        {
            ++j;
        }
        switch( j )
        {
            case 0:
                Sets[i].CGNSType = BAR_2;
                break;
            case 1:
                Sets[i].CGNSType = TRI_3;
                break;
            case 2:
                Sets[i].CGNSType = QUAD_4;
                break;
            case 3:
                Sets[i].CGNSType = TETRA_4;
                break;
            case 4:
                Sets[i].CGNSType = PYRA_5;
                break;
            case 5:
                Sets[i].CGNSType = PENTA_6;
                break;
            case 6:
                Sets[i].CGNSType = HEXA_8;
                break;
            default:
                std::cout << "It was not possible to identify the CGNSType\n";
                return MB_FAILURE;
        }
    }
    else
    {
        Sets[i].CGNSType = ElementTypeNull;
    }  // NOT SURE IF THAT'S THE RIGHT WAY.......
 
    // Clear the test vector
    Test.clear();
 
    return MB_SUCCESS;
}
 
// Fill the connectivity table
ErrorCode WriteCGNS::get_conn_table( std::vector< moab::EntityHandle >& Elements,
                                     std::vector< int >& Begin,
                                     std::vector< int >& End,
                                     std::vector< moab::Tag >& TagHandles,
                                     std::vector< WriteCGNS::SetStruct >& Sets,
                                     std::vector< std::vector< cgsize_t > >& ConnTable )
{
    ErrorCode rval;
 
    //   int Begin = 0; // GOT TO WORK ON THIS
    //   int End;
 
    // Get the number of Tags in the mesh
    int NbTags = TagHandles.size();
 
    // Test all Elements, get their ids and connectivity
    // to fill ConnTable
    for( std::vector< moab::EntityHandle >::iterator i = Elements.begin(); i != Elements.end(); ++i )
    {
        int id;
        // Test all Tags
        for( int j = 0; j < NbTags; ++j )
        {
            // Test if the Tag has data
            if( Sets[j].IdSet != -1 )
            {
                // Try to get data from entity
                rval = mbImpl->tag_get_data( TagHandles[j], &*i, 1, &id );
                if( MB_SUCCESS != rval )
                {
                    return rval;
                }
                // If successful id==j
                if( id == j )
                {
                    // Get the entity type of the EntityHandle wich points to Cells
                    int num_vtx;               // Number of MeshVertices in array connectivity.
                    const EntityHandle* conn;  // Array in which connectivity of entity_handle is returned.
                    // Gets a pointer to constant connectivity data of entity_handle
                    rval = mbImpl->get_connectivity( *i, conn, num_vtx );
                    if( MB_SUCCESS != rval )
                    {
                        return rval;
                    }
                    // If the Set is MIXED type
                    // push CGNS ENUM type of the entity
                    // before the connectivity
                    if( Sets[j].CGNSType == MIXED )
                    {
                        ConnTable[j].push_back( moab_cgns_conv( *i ) );  // moab_cgns_conv return an int which
                                                                         // represents the CGNS type
                        Begin[j]++;
                    }
                    End[j] = Begin[j] + num_vtx;
                    // Push conn in ConnTable in which "j" is the Set Index
                    for( int k = Begin[j]; k < End[j]; ++k )
                    {
                        ConnTable[j].push_back( (cgsize_t)conn[k - Begin[j]] );
                    }
                    Begin[j] = End[j];
                }
            }
        }
    }
    return MB_SUCCESS;
}
 
// Read the Moab type and return CGNS type
int WriteCGNS::moab_cgns_conv( const EntityHandle handle )
{
    EntityType MoabType = mbImpl->type_from_handle( handle );
    switch( MoabType )
    {
        case MBEDGE: /**< Mesh Edge */
            return BAR_2;
        case MBTRI: /**< Triangular element (including shells) */
            return TRI_3;
        case MBQUAD: /**< Quadrilateral element (including shells) */
            return QUAD_4;
        case MBTET: /**< Tetrahedral element */
            return TETRA_4;
        case MBPYRAMID: /**< Pyramid element */
            return PYRA_5;
        case MBPRISM: /**< Wedge element */
            return PENTA_6;
        case MBHEX: /**< Hexahedral element */
            return HEXA_8;
        default:
            std::cout << "It was not possible to identify the CGNSType\n";
            return 0;
    }
}
 
}  // namespace moab