Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
cub2h5m.cpp File Reference
#include "moab/GeomQueryTool.hpp"
#include "moab/GeomTopoTool.hpp"
#include "InitCGMA.hpp"
#include "CGMApp.hpp"
#include "moab/CN.hpp"
#include "moab/Core.hpp"
#include "moab/CartVect.hpp"
#include "moab/FileOptions.hpp"
#include "moab/Skinner.hpp"
#include "quads_to_tris.hpp"
#include <limits>
#include <cstdlib>
#include <sstream>
#include <ctime>
+ Include dependency graph for cub2h5m.cpp:

Go to the source code of this file.

Defines

#define GF_CUBIT_FILE_TYPE   "CUBIT"
#define GF_STEP_FILE_TYPE   "STEP"
#define GF_IGES_FILE_TYPE   "IGES"
#define GF_ACIS_TXT_FILE_TYPE   "ACIS_SAT"
#define GF_ACIS_BIN_FILE_TYPE   "ACIS_SAB"
#define GF_OCC_BREP_FILE_TYPE   "OCC"

Functions

void tokenize (const std::string &str, std::vector< std::string > &tokens, const char *delimiters)
ErrorCode get_group_names (Interface *MBI, const EntityHandle group_set, const Tag nameTag, std::vector< std::string > &grp_names)
ErrorCode summarize_cell_volume_change (Interface *MBI, const EntityHandle cgm_file_set, const Tag categoryTag, const Tag dimTag, const Tag sizeTag, const Tag nameTag, const Tag idTag, const bool conserve_mass, const bool debug)
ErrorCode remove_empty_cgm_surfs_and_vols (Interface *MBI, const EntityHandle cgm_file_set, const Tag idTag, const Tag dimTag, const bool)
ErrorCode build_new_surface (Interface *MBI, EntityHandle &new_surf, const EntityHandle forward_parent_vol, const EntityHandle reverse_parent_vol, const int new_surf_id, const Tag dimTag, const Tag idTag, const Tag categoryTag, const Tag senseTag)
ErrorCode orient_faces_outward (Interface *MBI, const Range &faces, const bool)
void plot_histogram (const std::string &title, const std::string &x_axis_label, const std::string &y_axis_label, const int n_bins, const double data[], const int n_data)
void generate_plots (const double orig[], const double defo[], const int n_elems, const std::string &time_step)
static double tet_volume (const CartVect &v0, const CartVect &v1, const CartVect &v2, const CartVect &v3)
double measure (Interface *MBI, const EntityHandle element)
ErrorCode get_signed_volume (Interface *MBI, const EntityHandle surf_set, const CartVect &offset, double &signed_volume)
ErrorCode fix_surface_senses (Interface *MBI, const EntityHandle cgm_file_set, const EntityHandle cub_file_set, const Tag idTag, const Tag dimTag, const Tag senseTag, const bool debug)
ErrorCode replace_faceted_cgm_surfs (Interface *MBI, const EntityHandle cgm_file_set, const EntityHandle cub_file_set, const Tag idTag, const Tag dimTag, const bool debug)
ErrorCode add_dead_elems_to_impl_compl (Interface *MBI, const EntityHandle cgm_file_set, const EntityHandle cub_file_set, const Tag idTag, const Tag dimTag, const Tag categoryTag, const Tag senseTag, const bool debug)
const char * get_geom_file_type (const char *filename)
const char * get_geom_fptr_type (FILE *file)
int is_cubit_file (FILE *file)
int is_step_file (FILE *file)
int is_iges_file (FILE *file)
int is_acis_txt_file (FILE *file)
int is_acis_bin_file (FILE *file)
int is_occ_brep_file (FILE *file)
int main (int argc, char *argv[])

Variables

double DEFAULT_DISTANCE = 0.001
double DEFAULT_LEN = 0.0
int DEFAULT_NORM = 5

Define Documentation

#define GF_ACIS_BIN_FILE_TYPE   "ACIS_SAB"

Definition at line 20 of file cub2h5m.cpp.

Referenced by get_geom_fptr_type().

#define GF_ACIS_TXT_FILE_TYPE   "ACIS_SAT"

Definition at line 19 of file cub2h5m.cpp.

Referenced by get_geom_fptr_type().

#define GF_CUBIT_FILE_TYPE   "CUBIT"

Definition at line 16 of file cub2h5m.cpp.

Referenced by get_geom_fptr_type().

#define GF_IGES_FILE_TYPE   "IGES"

Definition at line 18 of file cub2h5m.cpp.

Referenced by get_geom_fptr_type().

#define GF_OCC_BREP_FILE_TYPE   "OCC"

Definition at line 21 of file cub2h5m.cpp.

Referenced by get_geom_fptr_type().

#define GF_STEP_FILE_TYPE   "STEP"

Definition at line 17 of file cub2h5m.cpp.

Referenced by get_geom_fptr_type().


Function Documentation

ErrorCode add_dead_elems_to_impl_compl ( Interface MBI,
const EntityHandle  cgm_file_set,
const EntityHandle  cub_file_set,
const Tag  idTag,
const Tag  dimTag,
const Tag  categoryTag,
const Tag  senseTag,
const bool  debug 
)

Definition at line 860 of file cub2h5m.cpp.

References moab::Interface::add_entities(), moab::Range::begin(), build_new_surface(), moab::Range::empty(), moab::Range::end(), ErrorCode, moab::Skinner::find_skin(), moab::Range::front(), moab::Interface::get_child_meshsets(), moab::Interface::get_entities_by_type(), moab::Interface::get_entities_by_type_and_tag(), idTag, moab::intersect(), MB_INVALID_SIZE, MB_SUCCESS, MBENTITYSET, MBHEX, MBQUAD, orient_faces_outward(), moab::Interface::remove_entities(), moab::Range::size(), moab::subtract(), and moab::Interface::tag_get_data().

Referenced by main().

{

    // Get the cgm surfaces
    ErrorCode result;
    const int two               = 2;
    const void* const two_val[] = { &two };
    Range cgm_surfs;
    result = MBI->get_entities_by_type_and_tag( cgm_file_set, MBENTITYSET, &dimTag, two_val, 1, cgm_surfs );
    if( MB_SUCCESS != result ) return result;

    // Get the maximum surface id. This is so that new surfaces do not have
    // duplicate ids.
    int max_surf_id = -1;
    for( Range::const_iterator i = cgm_surfs.begin(); i != cgm_surfs.end(); ++i )
    {
        int surf_id;
        result = MBI->tag_get_data( idTag, &( *i ), 1, &surf_id );
        if( MB_SUCCESS != result ) return result;
        if( max_surf_id < surf_id ) max_surf_id = surf_id;
    }
    std::cout << "  Maximum surface id=" << max_surf_id << std::endl;

    // For each cgm volume, does a cub volume with the same id exist?
    const int three               = 3;
    const void* const three_val[] = { &three };
    Range cgm_vols;
    result = MBI->get_entities_by_type_and_tag( cgm_file_set, MBENTITYSET, &dimTag, three_val, 1, cgm_vols );
    if( MB_SUCCESS != result ) return result;

    // get the corresponding cub volume
    for( Range::iterator i = cgm_vols.begin(); i != cgm_vols.end(); ++i )
    {
        int vol_id;
        result = MBI->tag_get_data( idTag, &( *i ), 1, &vol_id );
        assert( MB_SUCCESS == result );
        if( MB_SUCCESS != result ) return result;
        std::cout << "  Volume " << vol_id;

        // Find the meshed vol set with the same id
        Range cub_vol;
        const Tag tags[]             = { idTag, dimTag };
        const void* const tag_vals[] = { &vol_id, &three };
        result = MBI->get_entities_by_type_and_tag( cub_file_set, MBENTITYSET, tags, tag_vals, 2, cub_vol );
        assert( MB_SUCCESS == result );
        if( MB_SUCCESS != result ) return result;
        if( 1 != cub_vol.size() )
        {
            std::cout << ": no meshed representation found" << std::endl;
            continue;
        }
        else
        {
            std::cout << std::endl;
        }

        // get the mesh elements of the volume.
        Range elems;
        result = MBI->get_entities_by_type( cub_vol.front(), MBHEX, elems );
        assert( MB_SUCCESS == result );
        if( MB_SUCCESS != result ) return result;
        if( debug ) std::cout << "    found " << elems.size() << " hex elems" << std::endl;

        // skin the volumes
        Skinner tool( MBI );
        Range skin_faces;
        result = tool.find_skin( 0, elems, 2, skin_faces, true );
        assert( MB_SUCCESS == result );
        if( MB_SUCCESS != result ) return result;

        // Reconcile the difference between faces of the cub file surfaces and skin
        // faces of the 3D mesh. The difference exists because dead elements have been
        // removed. Faces are divided into:
        //   cub_faces    (in) - the faces in the cub file surface
        //   skin_faces   (in) - the faces of the 3D mesh elements
        //   common_faces (out)- the faces common to both the cub file surface and 3D mesh
        //                       they are still adjacent to this volume
        //   old_faces    (out)- the faces of the cub surface not on the 3D mesh skin
        //                       they are no longer adjacent to this vol
        //   new_faces    (out)- the faces of the 3D mesh skin not on the cub surface
        //                       they are now adjacent to this volume

        // get cub child surfaces.
        Range cub_surfs;
        result = MBI->get_child_meshsets( cub_vol.front(), cub_surfs );
        assert( MB_SUCCESS == result );
        if( MB_SUCCESS != result ) return result;
        for( Range::iterator j = cub_surfs.begin(); j != cub_surfs.end(); ++j )
        {
            int surf_id;
            result = MBI->tag_get_data( idTag, &( *j ), 1, &surf_id );
            assert( MB_SUCCESS == result );
            if( MB_SUCCESS != result ) return result;
            // get the quads on each surface
            Range cub_faces;
            result = MBI->get_entities_by_type( *j, MBQUAD, cub_faces );
            assert( MB_SUCCESS == result );
            if( MB_SUCCESS != result ) return result;
            // Meshed volumes must have meshed surfaces
            if( cub_faces.empty() )
            {
                std::cout << "    Surface " << surf_id << ": contains no meshed faces" << std::endl;
                // return MB_ENTITY_NOT_FOUND;
            }
            // get the faces common to both the skin and this surface
            Range common_faces = intersect( cub_faces, skin_faces );
            // find the surface faces not on the skin - these are old and need removed
            Range old_faces = subtract( cub_faces, common_faces );
            result          = MBI->remove_entities( *j, old_faces );
            assert( MB_SUCCESS == result );
            if( MB_SUCCESS != result ) return result;

            // remove the common faces from the skin faces
            skin_faces = subtract( skin_faces, common_faces );
            // If no old faces exist we are done
            if( old_faces.empty() ) continue;
            std::cout << "    Surface " << surf_id << ": " << old_faces.size() << " old faces removed" << std::endl;
            // Place the old faces in a new surface, because they may still be adjacent
            // to 3D mesh in another volume. Get the parent vols of the surface.
            Range cgm_surf;
            const Tag tags2[]             = { idTag, dimTag };
            const void* const tag_vals2[] = { &surf_id, &two };
            result = MBI->get_entities_by_type_and_tag( cgm_file_set, MBENTITYSET, tags2, tag_vals2, 2, cgm_surf );
            assert( MB_SUCCESS == result );
            if( MB_SUCCESS != result ) return result;
            if( 1 != cgm_surf.size() )
            {
                std::cout << "invalid size" << std::endl;
                return MB_INVALID_SIZE;
            }
            EntityHandle cgm_vols2[2], cub_vols[2];
            result = MBI->tag_get_data( senseTag, cgm_surf, &cgm_vols2 );
            assert( MB_SUCCESS == result );
            if( MB_SUCCESS != result ) return result;
            result = MBI->tag_get_data( senseTag, &( *j ), 1, &cub_vols );
            assert( MB_SUCCESS == result );
            if( MB_SUCCESS != result ) return result;
            // for the new surf, replace the current volume with the impl compl vol.
            // This is because the faces that no longer exist will become adjacent to
            // the impl compl
            if( *i == cgm_vols2[0] )
            {
                cgm_vols2[0] = 0;
                cub_vols[0]  = 0;
            }
            if( *i == cgm_vols2[1] )
            {
                cgm_vols2[1] = 0;
                cub_vols[1]  = 0;
            }
            // If both sides of the surface are the impl comp, do not create the surface.
            if( 0 == cgm_vols2[0] && 0 == cgm_vols2[1] )
            {
                std::cout << "    New surface was not created for old faces because both parents "
                             "are impl_compl volume "
                          << std::endl;
                continue;
            }
            // build the new surface.
            EntityHandle new_cgm_surf, new_cub_surf;
            ++max_surf_id;
            result = build_new_surface( MBI, new_cgm_surf, cgm_vols2[0], cgm_vols2[1], max_surf_id, dimTag, idTag,
                                        categoryTag, senseTag );
            assert( MB_SUCCESS == result );
            if( MB_SUCCESS != result ) return result;
            result = build_new_surface( MBI, new_cub_surf, cub_vols[0], cub_vols[1], max_surf_id, dimTag, idTag,
                                        categoryTag, senseTag );
            assert( MB_SUCCESS == result );
            if( MB_SUCCESS != result ) return result;
            // add the new surface to the file set and populate it with the old faces
            result = MBI->add_entities( cgm_file_set, &new_cgm_surf, 1 );
            assert( MB_SUCCESS == result );
            if( MB_SUCCESS != result ) return result;
            assert( MB_SUCCESS == result );
            result = MBI->add_entities( cub_file_set, &new_cub_surf, 1 );
            if( MB_SUCCESS != result ) return result;
            assert( MB_SUCCESS == result );
            result = MBI->add_entities( new_cub_surf, old_faces );
            assert( MB_SUCCESS == result );
            if( MB_SUCCESS != result ) return result;

            std::cout << "    Surface " << max_surf_id << ": created for " << old_faces.size() << " old faces"
                      << std::endl;
        }

        // the remaining skin faces are newly exposed faces
        Range new_faces = skin_faces;

        // new skin faces must be assigned to a surface
        if( new_faces.empty() ) continue;
        std::cout << "    Surface " << max_surf_id + 1 << ": created for " << new_faces.size() << " new faces"
                  << std::endl;

        // Ensure that faces are oriented outwards
        result = orient_faces_outward( MBI, new_faces, debug );
        assert( MB_SUCCESS == result );
        if( MB_SUCCESS != result ) return result;

        // Create the new surface.
        EntityHandle new_cgm_surf, new_cub_surf;
        ++max_surf_id;
        result = build_new_surface( MBI, new_cgm_surf, *i, 0, max_surf_id, dimTag, idTag, categoryTag, senseTag );
        assert( MB_SUCCESS == result );
        if( MB_SUCCESS != result ) return result;
        result = build_new_surface( MBI, new_cub_surf, cub_vol.front(), 0, max_surf_id, dimTag, idTag, categoryTag,
                                    senseTag );
        assert( MB_SUCCESS == result );
        if( MB_SUCCESS != result ) return result;
        // Insert the new surf into file sets and populate it with faces.
        result = MBI->add_entities( cgm_file_set, &new_cgm_surf, 1 );
        assert( MB_SUCCESS == result );
        if( MB_SUCCESS != result ) return result;
        result = MBI->add_entities( cub_file_set, &new_cub_surf, 1 );
        assert( MB_SUCCESS == result );
        if( MB_SUCCESS != result ) return result;
        result = MBI->add_entities( new_cub_surf, new_faces );
        assert( MB_SUCCESS == result );
        if( MB_SUCCESS != result ) return result;
    }

    return MB_SUCCESS;
}
ErrorCode build_new_surface ( Interface MBI,
EntityHandle new_surf,
const EntityHandle  forward_parent_vol,
const EntityHandle  reverse_parent_vol,
const int  new_surf_id,
const Tag  dimTag,
const Tag  idTag,
const Tag  categoryTag,
const Tag  senseTag 
)

Definition at line 300 of file cub2h5m.cpp.

References moab::Interface::add_parent_child(), CATEGORY_TAG_SIZE, moab::Interface::create_meshset(), ErrorCode, moab::geom_category, MB_SUCCESS, and moab::Interface::tag_set_data().

Referenced by add_dead_elems_to_impl_compl().

{

    ErrorCode result;
    result = MBI->create_meshset( 0, new_surf );
    if( MB_SUCCESS != result ) return result;
    if( 0 != forward_parent_vol )
    {
        result = MBI->add_parent_child( forward_parent_vol, new_surf );
        if( MB_SUCCESS != result ) return result;
    }
    if( 0 != reverse_parent_vol )
    {
        result = MBI->add_parent_child( reverse_parent_vol, new_surf );
        if( MB_SUCCESS != result ) return result;
    }
    const int two = 2;
    result        = MBI->tag_set_data( dimTag, &new_surf, 1, &two );
    if( MB_SUCCESS != result ) return result;
    result = MBI->tag_set_data( idTag, &new_surf, 1, &new_surf_id );
    if( MB_SUCCESS != result ) return result;
    const char geom_category[CATEGORY_TAG_SIZE] = { "Surface\0" };
    result                                      = MBI->tag_set_data( categoryTag, &new_surf, 1, &geom_category );
    if( MB_SUCCESS != result ) return result;
    EntityHandle vols[2] = { forward_parent_vol, reverse_parent_vol };
    result               = MBI->tag_set_data( senseTag, &new_surf, 1, vols );
    if( MB_SUCCESS != result ) return result;

    return MB_SUCCESS;
}
ErrorCode fix_surface_senses ( Interface MBI,
const EntityHandle  cgm_file_set,
const EntityHandle  cub_file_set,
const Tag  idTag,
const Tag  dimTag,
const Tag  senseTag,
const bool  debug 
)

Definition at line 694 of file cub2h5m.cpp.

References moab::Interface::add_entities(), moab::Range::begin(), moab::Range::end(), ErrorCode, moab::Range::front(), moab::Interface::get_entities_by_type(), moab::Interface::get_entities_by_type_and_tag(), get_signed_volume(), idTag, make_tris_from_quads(), MB_SUCCESS, MBENTITYSET, MBQUAD, moab::Range::size(), moab::Interface::tag_get_data(), and moab::Interface::tag_set_data().

Referenced by main().

{
    ErrorCode result;
    const int two               = 2;
    const void* const two_val[] = { &two };
    Range cgm_surfs;
    result = MBI->get_entities_by_type_and_tag( cgm_file_set, MBENTITYSET, &dimTag, two_val, 1, cgm_surfs );
    if( MB_SUCCESS != result ) return result;
    for( Range::iterator i = cgm_surfs.begin(); i != cgm_surfs.end(); ++i )
    {
        int surf_id;
        result = MBI->tag_get_data( idTag, &( *i ), 1, &surf_id );
        if( MB_SUCCESS != result ) return result;

        // Find the meshed surface set with the same id
        Range cub_surf;
        const Tag tags[]             = { idTag, dimTag };
        const void* const tag_vals[] = { &surf_id, &two };
        result = MBI->get_entities_by_type_and_tag( cub_file_set, MBENTITYSET, tags, tag_vals, 2, cub_surf );
        if( MB_SUCCESS != result ) return result;
        if( 1 != cub_surf.size() )
        {
            std::cout << "  Surface " << surf_id << ": no meshed representation found, using CAD representation instead"
                      << std::endl;
            continue;
        }

        // Get tris that represent the quads of the cub surf
        Range quads;
        result = MBI->get_entities_by_type( cub_surf.front(), MBQUAD, quads );
        if( MB_SUCCESS != result ) return result;
        Range cub_tris;
        result = make_tris_from_quads( MBI, quads, cub_tris );
        if( MB_SUCCESS != result ) return result;

        // Add the tris to the same surface meshset as the quads are inside.
        result = MBI->add_entities( cub_surf.front(), cub_tris );
        if( MB_SUCCESS != result ) return result;

        // get the signed volume for each surface representation. Keep trying until
        // The signed volumes are much greater than numerical precision. Planar
        // surfaces will have a signed volume of zero if the plane goes through the
        // origin, unless we apply an offset.
        const int n_attempts        = 100;
        const int max_random        = 10;
        const double min_signed_vol = 0.1;
        double cgm_signed_vol, cub_signed_vol;
        for( int j = 0; j < n_attempts; ++j )
        {
            cgm_signed_vol = 0;
            cub_signed_vol = 0;
            CartVect offset( std::rand() % max_random, std::rand() % max_random, std::rand() % max_random );
            result = get_signed_volume( MBI, *i, offset, cgm_signed_vol );
            if( MB_SUCCESS != result ) return result;
            result = get_signed_volume( MBI, cub_surf.front(), offset, cub_signed_vol );
            if( MB_SUCCESS != result ) return result;
            if( debug )
                std::cout << "  surf_id=" << surf_id << " cgm_signed_vol=" << cgm_signed_vol
                          << " cub_signed_vol=" << cub_signed_vol << std::endl;
            if( fabs( cgm_signed_vol ) > min_signed_vol && fabs( cub_signed_vol ) > min_signed_vol ) break;
            if( n_attempts == j + 1 )
            {
                std::cout << "error: signed volume could not be calculated unambiguously" << std::endl;
                return MB_FAILURE;
            }
        }

        // If the sign is different, reverse the cgm senses so that both
        // representations have the same signed volume.
        if( ( cgm_signed_vol < 0 && cub_signed_vol > 0 ) || ( cgm_signed_vol > 0 && cub_signed_vol < 0 ) )
        {
            EntityHandle cgm_surf_volumes[2], reversed_cgm_surf_volumes[2];
            result = MBI->tag_get_data( senseTag, &( *i ), 1, cgm_surf_volumes );
            if( MB_SUCCESS != result ) return result;
            if( MB_SUCCESS != result ) return result;

            reversed_cgm_surf_volumes[0] = cgm_surf_volumes[1];
            reversed_cgm_surf_volumes[1] = cgm_surf_volumes[0];

            result = MBI->tag_set_data( senseTag, &( *i ), 1, reversed_cgm_surf_volumes );
            if( MB_SUCCESS != result ) return result;
        }
    }

    return MB_SUCCESS;
}
void generate_plots ( const double  orig[],
const double  defo[],
const int  n_elems,
const std::string &  time_step 
)

Definition at line 570 of file cub2h5m.cpp.

References plot_histogram().

Referenced by main().

{

    // find volume ratio then max and min
    double* ratio = new double[n_elems];
    for( int i = 0; i < n_elems; ++i )
        ratio[i] = ( defo[i] - orig[i] ) / orig[i];

    plot_histogram( "Predeformed Element Volume", "Num_Elems", "Volume [cc]", 10, orig, n_elems );
    std::string title = "Element Volume Change Ratio at Time Step " + time_step;
    plot_histogram( title, "Num_Elems", "Volume Ratio", 10, ratio, n_elems );
    delete[] ratio;
}
const char * get_geom_file_type ( const char *  filename)

Definition at line 1440 of file cub2h5m.cpp.

References get_geom_fptr_type().

Referenced by main().

{
    FILE* file;
    const char* result = 0;

    file = fopen( name, "r" );
    if( file )
    {
        result = get_geom_fptr_type( file );
        fclose( file );
    }

    return result;
}
const char * get_geom_fptr_type ( FILE *  file)

Definition at line 1455 of file cub2h5m.cpp.

References GF_ACIS_BIN_FILE_TYPE, GF_ACIS_TXT_FILE_TYPE, GF_CUBIT_FILE_TYPE, GF_IGES_FILE_TYPE, GF_OCC_BREP_FILE_TYPE, GF_STEP_FILE_TYPE, is_acis_bin_file(), is_acis_txt_file(), is_cubit_file(), is_iges_file(), is_occ_brep_file(), and is_step_file().

Referenced by get_geom_file_type().

{
    static const char* CUBIT_NAME = GF_CUBIT_FILE_TYPE;
    static const char* STEP_NAME  = GF_STEP_FILE_TYPE;
    static const char* IGES_NAME  = GF_IGES_FILE_TYPE;
    static const char* SAT_NAME   = GF_ACIS_TXT_FILE_TYPE;
    static const char* SAB_NAME   = GF_ACIS_BIN_FILE_TYPE;
    static const char* BREP_NAME  = GF_OCC_BREP_FILE_TYPE;

    if( is_cubit_file( file ) )
        return CUBIT_NAME;
    else if( is_step_file( file ) )
        return STEP_NAME;
    else if( is_iges_file( file ) )
        return IGES_NAME;
    else if( is_acis_bin_file( file ) )
        return SAB_NAME;
    else if( is_acis_txt_file( file ) )
        return SAT_NAME;
    else if( is_occ_brep_file( file ) )
        return BREP_NAME;
    else
        return 0;
}
ErrorCode get_group_names ( Interface MBI,
const EntityHandle  group_set,
const Tag  nameTag,
std::vector< std::string > &  grp_names 
)

Definition at line 41 of file cub2h5m.cpp.

References ErrorCode, MB_SUCCESS, MB_TAG_NOT_FOUND, NAME_TAG_SIZE, and moab::Interface::tag_get_data().

Referenced by summarize_cell_volume_change().

{
    // get names
    char name0[NAME_TAG_SIZE];
    std::fill( name0, name0 + NAME_TAG_SIZE, '\0' );
    ErrorCode result = MBI->tag_get_data( nameTag, &group_set, 1, &name0 );
    if( MB_SUCCESS != result && MB_TAG_NOT_FOUND != result ) return MB_FAILURE;

    if( MB_TAG_NOT_FOUND != result ) grp_names.push_back( std::string( name0 ) );

    return MB_SUCCESS;
}
ErrorCode get_signed_volume ( Interface MBI,
const EntityHandle  surf_set,
const CartVect offset,
double &  signed_volume 
)

Definition at line 653 of file cub2h5m.cpp.

References moab::Range::begin(), moab::Range::end(), ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), moab::Interface::get_entities_by_type(), MB_INVALID_SIZE, MB_SUCCESS, and MBTRI.

Referenced by fix_surface_senses().

{
    ErrorCode rval;
    Range tris;
    rval = MBI->get_entities_by_type( surf_set, MBTRI, tris );
    if( MB_SUCCESS != rval ) return rval;
    signed_volume = 0.0;
    const EntityHandle* conn;
    int len;
    CartVect coords[3];
    for( Range::iterator j = tris.begin(); j != tris.end(); ++j )
    {
        rval = MBI->get_connectivity( *j, conn, len, true );
        if( MB_SUCCESS != rval ) return rval;
        if( 3 != len ) return MB_INVALID_SIZE;
        rval = MBI->get_coords( conn, 3, coords[0].array() );
        if( MB_SUCCESS != rval ) return rval;

        // Apply offset to avoid calculating 0 for cases when the origin is in the
        // plane of the surface.
        for( unsigned int k = 0; k < 3; ++k )
        {
            coords[k][0] += offset[0];
            coords[k][1] += offset[1];
            coords[k][2] += offset[2];
        }

        coords[1] -= coords[0];
        coords[2] -= coords[0];
        signed_volume += ( coords[0] % ( coords[1] * coords[2] ) );
    }
    return MB_SUCCESS;
}
int is_acis_bin_file ( FILE *  file)

Definition at line 1498 of file cub2h5m.cpp.

References buffer.

Referenced by get_geom_fptr_type().

{
    char buffer[15];
    return !fseek( file, 0, SEEK_SET ) && fread( buffer, 15, 1, file ) && !memcmp( buffer, "ACIS BinaryFile", 9 );
}
int is_acis_txt_file ( FILE *  file)

Definition at line 1504 of file cub2h5m.cpp.

References buffer, and length().

Referenced by get_geom_fptr_type().

{
    char buffer[5];
    int version, length;

    if( fseek( file, 0, SEEK_SET ) || 2 != fscanf( file, "%d %*d %*d %*d %d ", &version, &length ) ) return 0;

    if( version < 1 || version > 0xFFFF ) return 0;

    // Skip application name
    if( fseek( file, length, SEEK_CUR ) ) return 0;

    // Read length of version string followed by first 5 characters
    if( 2 != fscanf( file, "%d %4s", &length, buffer ) ) return 0;

    return !strcmp( buffer, "ACIS" );
}
int is_cubit_file ( FILE *  file)

Definition at line 1480 of file cub2h5m.cpp.

References buffer.

Referenced by get_geom_fptr_type().

{
    unsigned char buffer[4];
    return !fseek( file, 0, SEEK_SET ) && fread( buffer, 4, 1, file ) && !memcmp( buffer, "CUBE", 4 );
}
int is_iges_file ( FILE *  file)

Definition at line 1492 of file cub2h5m.cpp.

References buffer.

Referenced by get_geom_fptr_type().

{
    unsigned char buffer[10];
    return !fseek( file, 72, SEEK_SET ) && fread( buffer, 10, 1, file ) && !memcmp( buffer, "S      1\r\n", 10 );
}
int is_occ_brep_file ( FILE *  file)

Definition at line 1522 of file cub2h5m.cpp.

References buffer.

Referenced by get_geom_fptr_type().

{
    unsigned char buffer[6];
    return !fseek( file, 0, SEEK_SET ) && fread( buffer, 6, 1, file ) && !memcmp( buffer, "DBRep_", 6 );
}
int is_step_file ( FILE *  file)

Definition at line 1486 of file cub2h5m.cpp.

References buffer.

Referenced by get_geom_fptr_type().

{
    unsigned char buffer[9];
    return !fseek( file, 0, SEEK_SET ) && fread( buffer, 9, 1, file ) && !memcmp( buffer, "ISO-10303", 9 );
}
int main ( int  argc,
char *  argv[] 
)

Definition at line 1129 of file cub2h5m.cpp.

References add_dead_elems_to_impl_compl(), moab::Range::begin(), CATEGORY_TAG_NAME, CATEGORY_TAG_SIZE, moab::Core::create_meshset(), moab::debug, moab::Range::end(), ErrorCode, fix_surface_senses(), generate_plots(), GEOM_DIMENSION_TAG_NAME, moab::Core::get_entities_by_dimension(), moab::Core::get_entities_by_type_and_tag(), get_geom_file_type(), moab::Core::get_last_error(), moab::Core::globalId_tag(), idTag, moab::Core::load_file(), MB_SUCCESS, MB_TAG_CREAT, MB_TAG_DENSE, MB_TAG_SPARSE, MB_TYPE_DOUBLE, MB_TYPE_HANDLE, MB_TYPE_INTEGER, MB_TYPE_OPAQUE, MBENTITYSET, moab::measure(), moab::GeomQueryTool::measure_volume(), MESHSET_TRACK_OWNER, NAME_TAG_NAME, NAME_TAG_SIZE, nameTag, remove_empty_cgm_surfs_and_vols(), replace_faceted_cgm_surfs(), moab::Core::set_meshset_options(), size, moab::Range::size(), start_time, summarize_cell_volume_change(), moab::Core::tag_get_handle(), moab::Core::tag_set_data(), and moab::Core::write_mesh().

{
    clock_t start_time = clock();

    const bool debug      = false;
    const char* file_type = NULL;

    const char* cub_name  = 0;
    const char* exo_name  = 0;
    const char* out_name  = 0;
    const char* time_step = 0;
    const char* sat_name  = 0;
    double dist_tol = 0.001, len_tol = 0.0;
    int norm_tol = 5;

    if( 6 != argc && 9 != argc )
    {
        std::cerr << "To read meshed geometry for MOAB:" << std::endl;
        std::cerr << "$> <cub_file.cub> <acis_file.sat> <facet_tol> <output_file.h5m> conserve_mass<bool>" << std::endl;
        std::cerr << "To read meshed geometry for MOAB and update node coordinates:" << std::endl;
        std::cerr << "$> <cub_file.cub> <acis_file.sat> <facet_tol> <output_file.h5m> "
                     "<deformed_exo_file.e> time_step<int> check_vol_change<bool> conserve_mass<bool>"
                  << std::endl;
        exit( 4 );
    }

    // check filenames for proper suffix
    std::string temp;
    cub_name = argv[1];
    temp.assign( cub_name );
    if( std::string::npos == temp.find( ".cub" ) )
    {
        std::cerr << "cub_file does not have correct suffix" << std::endl;
        return 1;
    }
    sat_name = argv[2];  // needed because the cub file's embedded sat file does not have groups
    temp.assign( sat_name );
    if( std::string::npos == temp.find( ".sat" ) )
    {
        std::cerr << "sat_file does not have correct suffix" << std::endl;
        return 1;
    }
    out_name = argv[4];
    temp.assign( out_name );
    if( std::string::npos == temp.find( ".h5m" ) )
    {
        std::cerr << "out_file does not have correct suffix" << std::endl;
        return 1;
    }

    // get the facet tolerance
    dist_tol = atof( argv[3] );
    if( 0 > dist_tol || 1 < dist_tol )
    {
        std::cout << "error: facet_tolerance=" << dist_tol << std::endl;
        return 1;
    }

    // Should the nodes be updated?
    bool update_coords = false;
    if( 9 == argc )
    {
        exo_name = argv[5];
        temp.assign( exo_name );
        if( std::string::npos == temp.find( ".e" ) )
        {
            std::cerr << "e_file does not have correct suffix" << std::endl;
            return 1;
        }
        time_step     = argv[6];
        update_coords = true;
    }

    // Should the volume change be determined?
    bool determine_volume_change = false;
    if( 9 == argc )
    {
        temp.assign( argv[7] );
        if( std::string::npos != temp.find( "true" ) ) determine_volume_change = true;
    }

    // Should densities be changed to conserve mass?
    bool conserve_mass = false;
    if( 9 == argc )
    {
        temp.assign( argv[8] );
        if( std::string::npos != temp.find( "true" ) ) conserve_mass = true;
    }
    else if( 6 == argc )
    {
        temp.assign( argv[5] );
        if( std::string::npos != temp.find( "true" ) ) conserve_mass = true;
    }

    // Get CGM file type
    if( !file_type )
    {
        file_type = get_geom_file_type( cub_name );
        if( !file_type )
        {
            std::cerr << cub_name << " : unknown file type, try '-t'" << std::endl;
            exit( 1 );
        }
    }

    // Read the mesh from the cub file with Tqcdfr
    Core* MBI = new Core();
    ErrorCode result;
    EntityHandle cub_file_set;
    result = MBI->create_meshset( 0, cub_file_set );
    if( MB_SUCCESS != result ) return result;
    // Do not ignore the Cubit file version. In testing, a cub file from Cubit12
    // did not work.
    // char cub_options[256] = "120";
    // char cub_options[256] = "IGNORE_VERSION";
    // result = MBI->load_file(cub_name, &cub_file_set, cub_options, NULL, 0, 0);
    result = MBI->load_file( cub_name, &cub_file_set, 0, NULL, 0, 0 );
    if( MB_SUCCESS != result )
    {
        std::cout << "error: problem reading cub file" << std::endl;
        return result;
    }
    std::cout << "Mesh file read." << std::endl;

    // Read the ACIS file with ReadCGM
    char cgm_options[256];
    std::cout << "  facet tolerance=" << dist_tol << std::endl;
    sprintf( cgm_options,
             "CGM_ATTRIBS=yes;FACET_DISTANCE_TOLERANCE=%g;FACET_NORMAL_TOLERANCE=%d;MAX_FACET_EDGE_"
             "LENGTH=%g;",
             dist_tol, norm_tol, len_tol );
    EntityHandle cgm_file_set;
    result = MBI->create_meshset( 0, cgm_file_set );
    if( MB_SUCCESS != result ) return result;
    result = MBI->load_file( sat_name, &cgm_file_set, cgm_options, NULL, 0, 0 );
    if( MB_SUCCESS != result )
    {
        std::cout << "error: problem reading sat file" << std::endl;
        return result;
    }
    std::cout << "CAD file read." << std::endl;

    // Create tags
    Tag dimTag, idTag, categoryTag, senseTag, sizeTag, nameTag;
    result = MBI->tag_get_handle( GEOM_DIMENSION_TAG_NAME, 1, MB_TYPE_INTEGER, dimTag, MB_TAG_SPARSE | MB_TAG_CREAT );
    if( MB_SUCCESS != result ) return result;
    idTag  = MBI->globalId_tag();
    result = MBI->tag_get_handle( CATEGORY_TAG_NAME, CATEGORY_TAG_SIZE, MB_TYPE_OPAQUE, categoryTag,
                                  MB_TAG_SPARSE | MB_TAG_CREAT );
    if( MB_SUCCESS != result ) return result;
    result = MBI->tag_get_handle( "GEOM_SENSE_2", 2, MB_TYPE_HANDLE, senseTag, MB_TAG_SPARSE | MB_TAG_CREAT );
    if( MB_SUCCESS != result ) return result;
    result = MBI->tag_get_handle( "GEOM_SIZE", 1, MB_TYPE_DOUBLE, sizeTag, MB_TAG_DENSE | MB_TAG_CREAT );
    if( MB_SUCCESS != result ) return result;
    result = MBI->tag_get_handle( NAME_TAG_NAME, NAME_TAG_SIZE, MB_TYPE_OPAQUE, nameTag, MB_TAG_SPARSE | MB_TAG_CREAT );
    if( MB_SUCCESS != result ) return result;

    // Create triangles from the quads of the cub surface sets and add them to the
    // cub surface sets. Get the signed volume of each surface for both cgm and
    // cub representations. Change the sense of the cgm representation to match
    // the cub representation.
    result = fix_surface_senses( MBI, cgm_file_set, cub_file_set, idTag, dimTag, senseTag, debug );
    if( MB_SUCCESS != result ) return result;
    std::cout << "Fixed CAD surface senses to match meshed surface senses." << std::endl;

    // Get the 3D elements in the cub file and measure their volume.
    Range orig_elems;
    std::vector< double > orig_size;
    if( determine_volume_change )
    {
        result = MBI->get_entities_by_dimension( 0, 3, orig_elems );
        if( MB_SUCCESS != result ) return result;
        orig_size.resize( orig_elems.size() );
        for( unsigned int i = 0; i < orig_elems.size(); ++i )
        {
            orig_size[i] = measure( MBI, orig_elems[i] );
        }
    }

    // Before updating the nodes and removing dead elements, force the cub volume
    // sets to track ownership so that dead elements will be deleted from the sets.
    const int three               = 3;
    const void* const three_val[] = { &three };
    Range cub_vols;
    result = MBI->get_entities_by_type_and_tag( cub_file_set, MBENTITYSET, &dimTag, three_val, 1, cub_vols );
    if( MB_SUCCESS != result ) return result;
    for( Range::const_iterator i = cub_vols.begin(); i != cub_vols.end(); ++i )
    {
        result = MBI->set_meshset_options( *i, MESHSET_TRACK_OWNER );
        if( MB_SUCCESS != result ) return result;
    }

    // Tag volume sets with the undeformed size of each volume.
    Range vol_sets;
    result = MBI->get_entities_by_type_and_tag( cgm_file_set, MBENTITYSET, &dimTag, three_val, 1, vol_sets );
    if( MB_SUCCESS != result ) return result;

    moab::GeomTopoTool gtt  = moab::GeomTopoTool( MBI, false );
    moab::GeomQueryTool gqt = moab::GeomQueryTool( &gtt );

    for( Range::const_iterator i = vol_sets.begin(); i != vol_sets.end(); ++i )
    {
        double size;
        result = gqt.measure_volume( *i, size );
        if( MB_SUCCESS != result ) return result;
        result = MBI->tag_set_data( sizeTag, &( *i ), 1, &size );
        if( MB_SUCCESS != result ) return result;
    }

    // Update the coordinates if needed. Do not do this before checking surface
    // sense, because the coordinate update could deform the surfaces too much
    // to make an accurate comparison.
    // The cub node ids are unique because cgm vertex ids are tagged on the vertex
    // meshset and not the vertex itself.
    // result = MBI->delete_entities( &cub_file_set, 1 );
    // if(MB_SUCCESS != result) return result;
    // Assume dead elements exist until I think of something better.
    bool dead_elements_exist = true;
    if( update_coords )
    {
        // ReadNCDF my_ex_reader(MBI);
        char exo_options[120] = "tdata=coord,";
        strcat( exo_options, time_step );
        strcat( exo_options, ",set" );
        // FileOptions exo_opts(exo_options)  ;
        // opts = "tdata=coord, 100, sum, temp.exo";
        // result =  my_ex_reader.load_file(exo_name, cgm_file_set, exo_opts, NULL, 0 , 0);
        // result =  my_ex_reader.load_file(exo_name, cub_file_set, exo_opts, NULL, 0 , 0);
        // result = my_ex_reader.load_file(exo_name, &cub_file_set, exo_opts, NULL, 0 , 0);
        MBI->load_file( exo_name, &cub_file_set, exo_options );
        if( MB_SUCCESS != result )
        {
            std::string last_error;
            MBI->get_last_error( last_error );
            std::cout << "coordinate update failed, " << last_error << std::endl;
            return result;
        }
        std::cout << "Updated mesh nodes with deformed coordinates from exodus file." << std::endl;
    }

    if( determine_volume_change )
    {
        // Dead elements have been removed by the deformation. Get the elements that
        // still exist.
        Range defo_elems;
        result = MBI->get_entities_by_dimension( 0, 3, defo_elems );
        if( MB_SUCCESS != result ) return result;

        // Determine the volume of the elements now that a deformation has been
        // applied. Condense the original array by removing dead elements.
        double* orig_size_condensed = new double[defo_elems.size()];
        double* defo_size_condensed = new double[defo_elems.size()];
        int j                       = 0;
        for( unsigned int i = 0; i < orig_elems.size(); ++i )
        {
            if( orig_elems[i] == defo_elems[j] )
            {
                orig_size_condensed[j] = orig_size[i];
                defo_size_condensed[j] = measure( MBI, defo_elems[j] );
                ++j;
            }
        }
        generate_plots( orig_size_condensed, defo_size_condensed, defo_elems.size(), std::string( time_step ) );
        delete[] orig_size_condensed;  // can't use the same delete[] for both
        delete[] defo_size_condensed;
    }

    // Deal with dead elements. For now, add them to the impl_compl volume.
    // Extra surfaces are created to do this.
    if( update_coords && dead_elements_exist )
    {
        result = add_dead_elems_to_impl_compl( MBI, cgm_file_set, cub_file_set, idTag, dimTag, categoryTag, senseTag,
                                               debug );
        if( MB_SUCCESS != result ) return result;
        std::cout << "Placed dead elements in implicit complement volume and added required surfaces." << std::endl;
    }

    // The quads in the cub_file_set have been updated for dead elements. For each
    // cgm_surf, if there exists a cub_surf with the same id, replace the cgm tris
    // with cub_tris (created from the quads). Note the a surface that is not
    // meshed (in cub file) will not be effected.
    result = replace_faceted_cgm_surfs( MBI, cgm_file_set, cub_file_set, idTag, dimTag, debug );
    if( MB_SUCCESS != result ) return result;
    std::cout << "Replaced faceted CAD surfaces with meshed surfaces of triangles." << std::endl;

    result = remove_empty_cgm_surfs_and_vols( MBI, cgm_file_set, idTag, dimTag, debug );
    if( MB_SUCCESS != result ) return result;
    std::cout << "Removed surfaces and volumes that no longer have any triangles." << std::endl;

    // For each material, sum the volume. If the coordinates were updated for
    // deformation, summarize the volume change.
    result = summarize_cell_volume_change( MBI, cgm_file_set, categoryTag, dimTag, sizeTag, nameTag, idTag,
                                           conserve_mass, debug );
    if( MB_SUCCESS != result ) return result;
    std::cout << "Summarized the volume change of each material, with respect to the solid model." << std::endl;

    result = MBI->write_mesh( out_name, &cgm_file_set, 1 );
    if( MB_SUCCESS != result )
    {
        std::cout << "write mesh failed" << std::endl;
        return result;
    }
    std::cout << "Saved output file for mesh-based analysis." << std::endl;

    clock_t end_time = clock();
    std::cout << "  " << (double)( end_time - start_time ) / CLOCKS_PER_SEC << " seconds" << std::endl;
    std::cout << std::endl;

    return 0;
}
double measure ( Interface MBI,
const EntityHandle  element 
)

Definition at line 591 of file cub2h5m.cpp.

References ErrorCode, moab::Interface::get_connectivity(), moab::Interface::get_coords(), length(), MB_SUCCESS, MBEDGE, MBHEX, MBPOLYGON, MBPRISM, MBPYRAMID, MBQUAD, MBTET, MBTRI, moab::sum(), tet_volume(), and moab::Interface::type_from_handle().

{
    EntityType type = MBI->type_from_handle( element );

    const EntityHandle* conn;
    int num_vertices;
    ErrorCode result = MBI->get_connectivity( element, conn, num_vertices );
    if( MB_SUCCESS != result ) return result;

    std::vector< CartVect > coords( num_vertices );
    result = MBI->get_coords( conn, num_vertices, coords[0].array() );
    if( MB_SUCCESS != result ) return result;

    switch( type )
    {
        case MBEDGE:
            return ( coords[0] - coords[1] ).length();
        case MBTRI:
            return 0.5 * ( ( coords[1] - coords[0] ) * ( coords[2] - coords[0] ) ).length();
        case MBQUAD:
            num_vertices = 4;
        case MBPOLYGON: {
            CartVect mid( 0, 0, 0 );
            for( int i = 0; i < num_vertices; ++i )
                mid += coords[i];
            mid /= num_vertices;

            double sum = 0.0;
            for( int i = 0; i < num_vertices; ++i )
            {
                int j = ( i + 1 ) % num_vertices;
                sum += ( ( mid - coords[i] ) * ( mid - coords[j] ) ).length();
            }
            return 0.5 * sum;
        }
        case MBTET:
            return tet_volume( coords[0], coords[1], coords[2], coords[3] );
        case MBPYRAMID:
            return tet_volume( coords[0], coords[1], coords[2], coords[4] ) +
                   tet_volume( coords[0], coords[2], coords[3], coords[4] );
        case MBPRISM:
            return tet_volume( coords[0], coords[1], coords[2], coords[5] ) +
                   tet_volume( coords[3], coords[5], coords[4], coords[0] ) +
                   tet_volume( coords[1], coords[4], coords[5], coords[0] );
        case MBHEX:
            return tet_volume( coords[0], coords[1], coords[3], coords[4] ) +
                   tet_volume( coords[7], coords[3], coords[6], coords[4] ) +
                   tet_volume( coords[4], coords[5], coords[1], coords[6] ) +
                   tet_volume( coords[1], coords[6], coords[3], coords[4] ) +
                   tet_volume( coords[2], coords[6], coords[3], coords[1] );
        default:
            return 0.0;
    }
}
ErrorCode orient_faces_outward ( Interface MBI,
const Range faces,
const bool   
)

Definition at line 341 of file cub2h5m.cpp.

References moab::Range::begin(), moab::CN::Dimension(), moab::Range::end(), ErrorCode, moab::Range::front(), moab::Interface::get_adjacencies(), moab::Interface::get_connectivity(), MB_INVALID_SIZE, MB_SUCCESS, moab::Interface::set_connectivity(), moab::CN::SideNumber(), moab::Range::size(), and moab::Interface::type_from_handle().

Referenced by add_dead_elems_to_impl_compl().

{

    ErrorCode result;
    for( Range::const_iterator i = faces.begin(); i != faces.end(); ++i )
    {
        Range adj_elem;
        result = MBI->get_adjacencies( &( *i ), 1, 3, false, adj_elem );
        if( MB_SUCCESS != result ) return result;
        if( 1 != adj_elem.size() ) return MB_INVALID_SIZE;

        // get connectivity for element and face
        const EntityHandle* elem_conn;
        int elem_n_nodes;
        result = MBI->get_connectivity( adj_elem.front(), elem_conn, elem_n_nodes );
        if( MB_SUCCESS != result ) return result;
        const EntityHandle* face_conn;
        int face_n_nodes;
        result = MBI->get_connectivity( *i, face_conn, face_n_nodes );
        if( MB_SUCCESS != result ) return result;

        // Get the sense of the face wrt the element
        EntityType elem_type = MBI->type_from_handle( adj_elem.front() );
        EntityType face_type = MBI->type_from_handle( *i );
        int face_num, offset;
        int sense          = 0;
        const int face_dim = CN::Dimension( face_type );
        int rval = CN::SideNumber( elem_type, elem_conn, face_conn, face_n_nodes, face_dim, face_num, sense, offset );
        if( 0 != rval ) return MB_FAILURE;

        // If the face is not oriented outward wrt the element, reverse it
        if( -1 == sense )
        {
            EntityHandle new_face_conn[4] = { face_conn[3], face_conn[2], face_conn[1], face_conn[0] };
            result                        = MBI->set_connectivity( *i, new_face_conn, 4 );
            if( MB_SUCCESS != result ) return result;
        }
    }
    return MB_SUCCESS;
}
void plot_histogram ( const std::string &  title,
const std::string &  x_axis_label,
const std::string &  y_axis_label,
const int  n_bins,
const double  data[],
const int  n_data 
)

Definition at line 500 of file cub2h5m.cpp.

Referenced by generate_plots().

{
    // find max and min
    double min = std::numeric_limits< double >::max();
    double max = -std::numeric_limits< double >::max();
    for( int i = 0; i < n_data; ++i )
    {
        if( min > data[i] ) min = data[i];
        if( max < data[i] ) max = data[i];
    }

    // make bins for histogram
    double bin_width = ( max - min ) / n_bins;
    std::vector< int > bins( n_bins );
    for( int i = 0; i < n_bins; ++i )
        bins[i] = 0;

    // fill the bins
    for( int i = 0; i < n_data; ++i )
    {
        double diff = data[i] - min;
        int bin     = diff / bin_width;
        if( 9 < bin ) bin = 9;  // cheap fix for numerical precision error
        if( 0 > bin ) bin = 0;  // cheap fix for numerical precision error
        ++bins[bin];
    }

    // create bars
    int max_bin = 0;
    for( int i = 0; i < n_bins; ++i )
        if( max_bin < bins[i] ) max_bin = bins[i];
    int bar_height;
    int max_bar_chars = 72;
    std::vector< std::string > bars( n_bins );
    for( int i = 0; i < n_bins; ++i )
    {
        bar_height = ( max_bar_chars * bins[i] ) / max_bin;
        for( int j = 0; j < bar_height; ++j )
            bars[i] += "*";
    }

    // print histogram header
    std::cout << std::endl;
    std::cout << "                                 " << title << std::endl;

    // print results
    std::cout.width( 15 );
    std::cout << min << std::endl;
    for( int i = 0; i < n_bins; ++i )
    {
        std::cout.width( 15 );
        std::cout << min + ( ( i + 1 ) * bin_width );
        std::cout.width( max_bar_chars );
        std::cout << bars[i] << bins[i] << std::endl;
    }

    // print histogram footer
    std::cout.width( 15 );
    std::cout << y_axis_label;
    std::cout.width( max_bar_chars );
    std::cout << " " << x_axis_label << std::endl;
    std::cout << std::endl;
}
ErrorCode remove_empty_cgm_surfs_and_vols ( Interface MBI,
const EntityHandle  cgm_file_set,
const Tag  idTag,
const Tag  dimTag,
const bool   
)

Definition at line 188 of file cub2h5m.cpp.

References moab::Range::begin(), moab::Interface::contains_entities(), moab::Interface::delete_entities(), moab::Range::end(), ErrorCode, moab::Interface::get_child_meshsets(), moab::Interface::get_entities_by_type(), moab::Interface::get_entities_by_type_and_tag(), moab::Interface::get_number_entities_by_type(), moab::Interface::get_parent_meshsets(), MB_SUCCESS, MBENTITYSET, MBTRI, moab::Interface::num_child_meshsets(), moab::Interface::remove_entities(), moab::Interface::remove_parent_child(), and moab::Interface::tag_get_data().

Referenced by main().

{

    ErrorCode result;
    const int two               = 2;
    const void* const two_val[] = { &two };
    Range cgm_surfs;
    result = MBI->get_entities_by_type_and_tag( cgm_file_set, MBENTITYSET, &dimTag, two_val, 1, cgm_surfs );
    if( MB_SUCCESS != result ) return result;

    for( Range::iterator i = cgm_surfs.begin(); i != cgm_surfs.end(); ++i )
    {
        int n_tris;
        result = MBI->get_number_entities_by_type( *i, MBTRI, n_tris );
        if( MB_SUCCESS != result ) return result;

        if( 0 == n_tris )
        {
            int surf_id;
            result = MBI->tag_get_data( idTag, &( *i ), 1, &surf_id );
            assert( MB_SUCCESS == result );

            Range parent_vols;
            result = MBI->get_parent_meshsets( *i, parent_vols );
            assert( MB_SUCCESS == result );
            for( Range::iterator j = parent_vols.begin(); j != parent_vols.end(); ++j )
            {
                result = MBI->remove_parent_child( *j, *i );
                assert( MB_SUCCESS == result );
            }
            Range child_curves;
            result = MBI->get_child_meshsets( *i, child_curves );
            assert( MB_SUCCESS == result );
            for( Range::iterator j = child_curves.begin(); j != child_curves.end(); ++j )
            {
                result = MBI->remove_parent_child( *i, *j );
                assert( MB_SUCCESS == result );
            }
            result = MBI->remove_entities( cgm_file_set, &( *i ), 1 );
            assert( MB_SUCCESS == result );

            // Is the set contained anywhere else? If the surface is in a CUBIT group,
            // such as "unmerged_surfs" it will cause write_mesh to fail. This should
            // be a MOAB bug.
            Range all_sets;
            result = MBI->get_entities_by_type( 0, MBENTITYSET, all_sets );
            assert( MB_SUCCESS == result );
            for( Range::iterator j = all_sets.begin(); j != all_sets.end(); ++j )
            {
                if( MBI->contains_entities( *j, &( *i ), 1 ) )
                {
                    result = MBI->remove_entities( *j, &( *i ), 1 );
                    assert( MB_SUCCESS == result );
                }
            }

            result = MBI->delete_entities( &( *i ), 1 );
            assert( MB_SUCCESS == result );
            std::cout << "  Surface " << surf_id << ": removed because all of its mesh faces have been removed"
                      << std::endl;
        }
    }

    // get all volumes
    const int three               = 3;
    const void* const three_val[] = { &three };
    Range cgm_vols;
    result = MBI->get_entities_by_type_and_tag( cgm_file_set, MBENTITYSET, &dimTag, three_val, 1, cgm_vols );
    if( MB_SUCCESS != result ) return result;

    for( Range::iterator i = cgm_vols.begin(); i != cgm_vols.end(); ++i )
    {
        // get the volume's number of surfaces
        int n_surfs;
        result = MBI->num_child_meshsets( *i, &n_surfs );
        assert( MB_SUCCESS == result );

        // if no surfaces remain, remove the volume
        if( 0 == n_surfs )
        {
            int vol_id;
            result = MBI->tag_get_data( idTag, &( *i ), 1, &vol_id );
            assert( MB_SUCCESS == result );
            // Is the set contained anywhere else? If the surface is in a CUBIT group,
            // such as "unmerged_surfs" it will cause write_mesh to fail. This should
            // be a MOAB bug.
            Range all_sets;
            result = MBI->get_entities_by_type( 0, MBENTITYSET, all_sets );
            assert( MB_SUCCESS == result );
            for( Range::iterator j = all_sets.begin(); j != all_sets.end(); ++j )
            {
                if( MBI->contains_entities( *j, &( *i ), 1 ) )
                {
                    result = MBI->remove_entities( *j, &( *i ), 1 );
                    assert( MB_SUCCESS == result );
                }
            }
            result = MBI->delete_entities( &( *i ), 1 );
            assert( MB_SUCCESS == result );
            std::cout << "  Volume " << vol_id << ": removed because all of its surfaces have been removed"
                      << std::endl;
        }
    }
    return MB_SUCCESS;
}
ErrorCode replace_faceted_cgm_surfs ( Interface MBI,
const EntityHandle  cgm_file_set,
const EntityHandle  cub_file_set,
const Tag  idTag,
const Tag  dimTag,
const bool  debug 
)

Definition at line 791 of file cub2h5m.cpp.

References moab::Interface::add_entities(), moab::Range::begin(), moab::Interface::delete_entities(), moab::Range::end(), ErrorCode, moab::Range::front(), moab::Interface::get_entities_by_type(), moab::Interface::get_entities_by_type_and_tag(), idTag, make_tris_from_quads(), MB_SUCCESS, MBENTITYSET, MBQUAD, MBTRI, moab::Interface::remove_entities(), moab::Range::size(), and moab::Interface::tag_get_data().

Referenced by main().

{
    ErrorCode result;
    const int two               = 2;
    const void* const two_val[] = { &two };
    Range cgm_surfs;
    result = MBI->get_entities_by_type_and_tag( cgm_file_set, MBENTITYSET, &dimTag, two_val, 1, cgm_surfs );
    if( MB_SUCCESS != result ) return result;

    for( Range::iterator i = cgm_surfs.begin(); i != cgm_surfs.end(); ++i )
    {
        int surf_id;
        result = MBI->tag_get_data( idTag, &( *i ), 1, &surf_id );
        if( MB_SUCCESS != result ) return result;
        if( debug ) std::cout << "surf_id=" << surf_id << std::endl;

        // Find the meshed surface set with the same id
        Range cub_surf;
        const Tag tags[]             = { idTag, dimTag };
        const void* const tag_vals[] = { &surf_id, &two };
        result = MBI->get_entities_by_type_and_tag( cub_file_set, MBENTITYSET, tags, tag_vals, 2, cub_surf );
        if( MB_SUCCESS != result ) return result;
        if( 1 != cub_surf.size() )
        {
            std::cout << "  Surface " << surf_id << ": no meshed representation found, using CAD representation instead"
                      << std::endl;
            continue;
        }

        // Get tris that represent the quads of the cub surf
        Range quads;
        result = MBI->get_entities_by_type( cub_surf.front(), MBQUAD, quads );
        if( MB_SUCCESS != result ) return result;

        Range cub_tris;
        result = make_tris_from_quads( MBI, quads, cub_tris );
        if( MB_SUCCESS != result ) return result;

        // Remove the tris from the cgm surf. Don't forget to remove them from the
        // cgm_file_set because it is not TRACKING.
        Range cgm_tris;
        result = MBI->get_entities_by_type( *i, MBTRI, cgm_tris );
        if( MB_SUCCESS != result ) return result;
        result = MBI->remove_entities( *i, cgm_tris );
        if( MB_SUCCESS != result ) return result;
        result = MBI->remove_entities( cgm_file_set, cgm_tris );
        if( MB_SUCCESS != result ) return result;
        result = MBI->delete_entities( cgm_tris );
        if( MB_SUCCESS != result ) return result;

        // Add the cub_tris to the cgm_surf
        result = MBI->add_entities( *i, cub_tris );
        if( MB_SUCCESS != result ) return result;
    }

    return MB_SUCCESS;
}
ErrorCode summarize_cell_volume_change ( Interface MBI,
const EntityHandle  cgm_file_set,
const Tag  categoryTag,
const Tag  dimTag,
const Tag  sizeTag,
const Tag  nameTag,
const Tag  idTag,
const bool  conserve_mass,
const bool  debug 
)

Definition at line 59 of file cub2h5m.cpp.

References moab::Interface::add_entities(), moab::Range::begin(), CATEGORY_TAG_SIZE, moab::Interface::create_meshset(), moab::Range::end(), ErrorCode, moab::Interface::get_entities_by_type_and_tag(), get_group_names(), MB_SUCCESS, MBENTITYSET, moab::GeomQueryTool::measure_volume(), MESHSET_SET, moab::Interface::remove_entities(), moab::Interface::tag_get_data(), moab::Interface::tag_set_data(), and tokenize().

Referenced by main().

{
    // get groups
    ErrorCode rval;
    const char group_category[CATEGORY_TAG_SIZE] = { "Group\0" };
    const void* const group_val[]                = { &group_category };
    Range groups;
    rval = MBI->get_entities_by_type_and_tag( 0, MBENTITYSET, &categoryTag, group_val, 1, groups );
    if( MB_SUCCESS != rval ) return rval;

    // Get the maximum group id. This is so that new groups do not have
    // duplicate ids.
    int max_grp_id = -1;
    for( Range::const_iterator i = groups.begin(); i != groups.end(); ++i )
    {
        int grp_id;
        rval = MBI->tag_get_data( idTag, &( *i ), 1, &grp_id );
        if( MB_SUCCESS != rval ) return rval;
        if( max_grp_id < grp_id ) max_grp_id = grp_id;
    }
    if( conserve_mass )
    {
        std::cout << "  Altering group densities to conserve mass for each volume." << std::endl;
        std::cout << "    maximum group id=" << max_grp_id << std::endl;
    }

    for( Range::const_iterator i = groups.begin(); i != groups.end(); ++i )
    {
        // get group names
        std::vector< std::string > grp_names;
        rval = get_group_names( MBI, *i, nameTag, grp_names );
        if( MB_SUCCESS != rval ) return MB_FAILURE;

        // determine if it is a material group
        bool material_grp = false;
        int mat_id        = -1;
        double rho        = 0;
        for( std::vector< std::string >::const_iterator j = grp_names.begin(); j != grp_names.end(); ++j )
        {
            if( std::string::npos != ( *j ).find( "mat" ) && std::string::npos != ( *j ).find( "rho" ) )
            {
                material_grp = true;
                std::cout << "  material group: " << *j << std::endl;

                // get the density and material id
                std::vector< std::string > tokens;
                tokenize( *j, tokens, "_" );
                mat_id = atoi( tokens[1].c_str() );
                rho    = atof( tokens[3].c_str() );
            }
        }
        if( !material_grp ) continue;

        // get the volume sets of the material group
        const int three               = 3;
        const void* const three_val[] = { &three };
        Range vols;
        rval = MBI->get_entities_by_type_and_tag( *i, MBENTITYSET, &dimTag, three_val, 1, vols );
        if( MB_SUCCESS != rval ) return rval;

        // for each volume, sum predeformed and deformed volume
        double orig_grp_volume = 0, defo_grp_volume = 0;

        moab::GeomTopoTool gtt  = moab::GeomTopoTool( MBI, false );
        moab::GeomQueryTool gqt = moab::GeomQueryTool( &gtt );
        for( Range::const_iterator j = vols.begin(); j != vols.end(); ++j )
        {
            double defo_size = 0, orig_size = 0;
            rval = gqt.measure_volume( *j, defo_size );
            if( MB_SUCCESS != rval ) return rval;
            defo_grp_volume += defo_size;
            rval = MBI->tag_get_data( sizeTag, &( *j ), 1, &orig_size );
            if( MB_SUCCESS != rval ) return rval;
            orig_grp_volume += orig_size;

            // calculate a new density to conserve mass through the deformation
            if( !conserve_mass ) continue;
            double new_rho = rho * orig_size / defo_size;

            // create a group for the volume with modified density
            EntityHandle new_grp;
            rval = MBI->create_meshset( MESHSET_SET, new_grp );
            if( MB_SUCCESS != rval ) return rval;
            std::stringstream new_name_ss;
            new_name_ss << "mat_" << mat_id << "_rho_" << new_rho << "\0";
            std::string new_name;
            new_name_ss >> new_name;
            rval = MBI->tag_set_data( nameTag, &new_grp, 1, new_name.c_str() );
            if( MB_SUCCESS != rval ) return rval;
            max_grp_id++;
            rval = MBI->tag_set_data( idTag, &new_grp, 1, &max_grp_id );
            if( MB_SUCCESS != rval ) return rval;
            const char group_category2[CATEGORY_TAG_SIZE] = "Group\0";
            rval = MBI->tag_set_data( categoryTag, &new_grp, 1, group_category2 );
            if( MB_SUCCESS != rval ) return rval;

            // add the volume to the new group
            rval = MBI->add_entities( new_grp, &( *j ), 1 );
            if( MB_SUCCESS != rval ) return rval;

            // add the new grp to the cgm_file_set
            rval = MBI->add_entities( cgm_file_set, &new_grp, 1 );
            if( MB_SUCCESS != rval ) return rval;

            // remove the volume from the old group
            rval = MBI->remove_entities( *i, &( *j ), 1 );
            if( MB_SUCCESS != rval ) return rval;
            if( debug ) std::cout << "    new group: " << new_name << " id=" << max_grp_id << std::endl;
        }

        std::cout << "    orig_volume=" << orig_grp_volume << " defo_volume=" << defo_grp_volume
                  << " defo/orig=" << defo_grp_volume / orig_grp_volume << std::endl;
    }

    return MB_SUCCESS;
}
static double tet_volume ( const CartVect v0,
const CartVect v1,
const CartVect v2,
const CartVect v3 
) [inline, static]

Definition at line 585 of file cub2h5m.cpp.

Referenced by measure().

{
    return 1. / 6. * ( ( ( v1 - v0 ) * ( v2 - v0 ) ) % ( v3 - v0 ) );
}
void tokenize ( const std::string &  str,
std::vector< std::string > &  tokens,
const char *  delimiters 
)

Definition at line 25 of file cub2h5m.cpp.

Referenced by summarize_cell_volume_change().

{
    std::string::size_type last = str.find_first_not_of( delimiters, 0 );
    std::string::size_type pos  = str.find_first_of( delimiters, last );
    if( std::string::npos == pos )
        tokens.push_back( str );
    else
        while( std::string::npos != pos && std::string::npos != last )
        {
            tokens.push_back( str.substr( last, pos - last ) );
            last = str.find_first_not_of( delimiters, pos );
            pos  = str.find_first_of( delimiters, last );
            if( std::string::npos == pos ) pos = str.size();
        }
}

Variable Documentation

double DEFAULT_DISTANCE = 0.001

Definition at line 1103 of file cub2h5m.cpp.

double DEFAULT_LEN = 0.0

Definition at line 1104 of file cub2h5m.cpp.

int DEFAULT_NORM = 5

Definition at line 1105 of file cub2h5m.cpp.

Referenced by moab::ReadCGM::set_options().

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines