Mesh Oriented datABase  (version 5.4.1)
Array-based unstructured mesh datastructure
main.cpp File Reference
#include <iostream>
#include <fstream>
#include <cmath>
#include <ctime>
#include <cstdlib>
#include <cassert>
#include "mcnpmit.hpp"
#include "moab/CartVect.hpp"
#include "moab/Core.hpp"
#include "MBTagConventions.hpp"
#include "moab/AdaptiveKDTree.hpp"
#include "moab/GeomUtil.hpp"
#include "moab/FileOptions.hpp"
#include "../tools/mbcoupler/ElemUtil.hpp"
+ Include dependency graph for main.cpp:

Go to the source code of this file.

Defines

#define MBI   mb_instance()

Functions

McnpDatamc_instance ()
moab::Interfacemb_instance ()
MCNPError read_files (int, char **)
MCNPError next_double (std::string, double &, int &)
MCNPError next_int (std::string, int &, int &)
int main (int argc, char **argv)

Variables

moab::Tag coord_tag
moab::Tag rotation_tag
moab::Tag cfd_heating_tag
moab::Tag cfd_error_tag
std::string h5m_filename
std::string CAD_filename
std::string output_filename
bool skip_build = false
bool read_qnv = false
clock_t start_time
clock_t load_time
clock_t build_time
clock_t interp_time

Define Documentation

#define MBI   mb_instance()

Definition at line 18 of file main.cpp.

Referenced by main().


Function Documentation

int main ( int  argc,
char **  argv 
)

Definition at line 37 of file main.cpp.

References moab::Range::begin(), box_max(), box_min(), build_time, moab::AdaptiveKDTree::build_tree(), CAD_filename, CARTESIAN, cfd_error_tag, cfd_heating_tag, moab::Range::clear(), coord_tag, CYLINDRICAL, moab::Range::end(), ErrorCode, moab::Tree::find_all_trees(), h5m_filename, interp_time, load_time, MB_CHK_ERR, MB_SUCCESS, MB_TAG_CREAT, MB_TAG_DENSE, MB_TYPE_DOUBLE, MB_TYPE_INTEGER, MBHEX, MBI, MBVERTEX, MCNP, MCNP_FAILURE, MCNP_SUCCESS, next_int(), nl, output_filename, moab::ElemUtil::point_in_trilinear_hex(), moab::AdaptiveKDTree::point_search(), read_files(), read_qnv, rotation_tag, moab::Range::size(), skip_build, and start_time.

{
    MCNPError result;

    start_time = clock();

    // Read in file names from command line
    result = read_files( argc, argv );
    if( result == MCNP_FAILURE ) return 1;

    result = MCNP->initialize_tags();

    // Parse the MCNP input file
    if( !skip_build )
    {

        result = MCNP->read_mcnpfile( skip_build );
        if( result == MCNP_FAILURE )
        {
            std::cout << "Failure reading MCNP file!" << std::endl;
            return 1;
        }
    }

    load_time = clock() - start_time;

    // Make the KD-Tree
    moab::ErrorCode MBresult;
    moab::AdaptiveKDTree kdtree( MBI );
    moab::EntityHandle root;

    MBI->tag_get_handle( "CoordTag", 1, moab::MB_TYPE_INTEGER, coord_tag, moab::MB_TAG_DENSE | moab::MB_TAG_CREAT );
    MBI->tag_get_handle( "RotationTag", 16, moab::MB_TYPE_DOUBLE, rotation_tag,
                         moab::MB_TAG_DENSE | moab::MB_TAG_CREAT );

    if( skip_build )
    {
        MBresult = MBI->load_mesh( h5m_filename.c_str() );

        if( moab::MB_SUCCESS == MBresult )
        {
            std::cout << std::endl << "Read in mesh from h5m file." << std::endl << std::endl;
            std::cout << "Querying mesh file..." << std::endl;
        }
        else
        {
            std::cout << "Failure reading h5m file!" << std::endl;
            std::cerr << "Error code: " << MBI->get_error_string( MBresult ) << " (" << MBresult << ")" << std::endl;
            std::string message;
            if( moab::MB_SUCCESS == MBI->get_last_error( message ) && !message.empty() )
                std::cerr << "Error message: " << message << std::endl;
            return 1;
        }

        moab::Range tmprange;
        kdtree.find_all_trees( tmprange );
        root = tmprange[0];
    }
    else
    {
        std::cout << "Building KD-Tree..." << std::endl;
        moab::FileOptions opts( "CANDIDATE_PLANE_SET=SUBDIVISION" );
        MBresult = kdtree.build_tree( MCNP->elem_handles, &root, &opts );
        if( MBresult == moab::MB_SUCCESS )
        {

            MBI->tag_set_data( coord_tag, &root, 1, &( MCNP->coord_system ) );
            MBI->tag_set_data( rotation_tag, &root, 1, &( MCNP->rotation_matrix ) );

            std::cout << "KD-Tree has been built successfully!" << std::endl << std::endl;
            MBresult = MBI->write_mesh( ( MCNP->MCNP_filename + ".h5m" ).c_str() );

            std::cout << "Querying mesh file..." << std::endl;
        }
        else
        {
            std::cout << "Error building KD-Tree!" << std::endl << std::endl;
            std::cerr << "Error code: " << MBI->get_error_string( MBresult ) << " (" << MBresult << ")" << std::endl;
            std::string message;
            if( moab::MB_SUCCESS == MBI->get_last_error( message ) && !message.empty() )
                std::cerr << "Error message: " << message << std::endl;
            return 1;
        }
    }

    int coord_sys;
    double rmatrix[16];

    MBresult = MBI->tag_get_data( coord_tag, &root, 1, &coord_sys );MB_CHK_ERR( MBresult );
    MBresult = MBI->tag_get_data( rotation_tag, &root, 1, &rmatrix );MB_CHK_ERR( MBresult );

    build_time = clock() - load_time;

    // Read the CAD mesh data and query the tree
    std::ifstream cadfile;
    std::ofstream outfile;

    outfile.open( output_filename.c_str() );

    int num_pts;
    int n;
    long int nl = 0;
    char* ctmp;
    int elems_read = 0;
    int p          = 0;
    char line[10000];

    // Used only when reading a mesh file to get vertex info
    double* cfd_coords = NULL;
    moab::Range::iterator cfd_iter;
    moab::EntityHandle meshset;

    if( read_qnv )
    {
        cadfile.open( CAD_filename.c_str() );
        cadfile.getline( line, 10000 );
        cadfile.getline( line, 10000 );
        result = next_int( line, num_pts, p );
    }
    else
    {

        meshset  = 0;
        MBresult = MBI->load_file( CAD_filename.c_str(), &meshset );MB_CHK_ERR( MBresult );
        assert( 0 != meshset );

        moab::Range cfd_verts;
        MBresult = MBI->get_entities_by_type( meshset, moab::MBVERTEX, cfd_verts, true );MB_CHK_ERR( MBresult );
        num_pts = cfd_verts.size();

        cfd_coords = new double[3 * num_pts];
        MBresult   = MBI->get_coords( cfd_verts, cfd_coords );MB_CHK_ERR( MBresult );

        cfd_iter = cfd_verts.begin();
        MBresult = MBI->tag_get_handle( "heating_tag", 1, moab::MB_TYPE_DOUBLE, cfd_heating_tag,
                                        moab::MB_TAG_DENSE | moab::MB_TAG_CREAT );MB_CHK_ERR( MBresult );
        MBresult = MBI->tag_get_handle( "error_tag", 1, moab::MB_TYPE_DOUBLE, cfd_error_tag,
                                        moab::MB_TAG_DENSE | moab::MB_TAG_CREAT );MB_CHK_ERR( MBresult );

        std::cout << std::endl << "Read in mesh with query points." << std::endl << std::endl;
    }

    double testpt[3];
    double transformed_pt[3];
    double taldata;
    double errdata;

    moab::CartVect testvc;

    bool found = false;

    // MBRange verts;
    std::vector< moab::EntityHandle > verts;
    moab::Range range;
    moab::CartVect box_max, box_min;

    moab::CartVect hexverts[8];
    moab::CartVect tmp_cartvect;
    std::vector< double > coords;

    double tal_sum = 0.0, err_sum = 0.0, tal_sum_sqr = 0.0, err_sum_sqr = 0.0;

    //  double davg = 0.0;
    //  unsigned int    nmax = 0, nmin = 1000000000 ;

    for( unsigned int i = 0; i < (unsigned int)num_pts; i++ )
    {

        // if (i%status_freq == 0)
        //  std::cerr << "Completed " << i/status_freq << "%" << std::endl;

        // Grab the coordinates to query
        if( read_qnv )
        {
            cadfile.getline( line, 10000 );

            nl        = std::strtol( line, &ctmp, 10 );
            n         = (unsigned int)nl;
            testpt[0] = std::strtod( ctmp + 1, &ctmp );
            testpt[1] = std::strtod( ctmp + 1, &ctmp );
            testpt[2] = std::strtod( ctmp + 1, NULL );
        }
        else
        {
            testpt[0] = cfd_coords[3 * i];
            testpt[1] = cfd_coords[3 * i + 1];
            testpt[2] = cfd_coords[3 * i + 2];
            n         = i + 1;
        }

        result = MCNP->transform_point( testpt, transformed_pt, coord_sys, rmatrix );

        testvc[0] = transformed_pt[0];
        testvc[1] = transformed_pt[1];
        testvc[2] = transformed_pt[2];

        // Find the leaf containing the point
        moab::EntityHandle tree_node;
        MBresult = kdtree.point_search( transformed_pt, tree_node );
        if( moab::MB_SUCCESS != MBresult )
        {
            double x = 0.0, y = 0.0, z = 0.0;
            if( CARTESIAN == coord_sys )
            {
                x = testvc[0];
                y = testvc[1];
                z = testvc[2];
            }
            else if( CYLINDRICAL == coord_sys )
            {
                x = testvc[0] * cos( 2 * M_PI * testvc[2] );
                y = testvc[0] * sin( 2 * M_PI * testvc[2] );
                z = testvc[1];
            }
            else
            {
                std::cout << "MOAB WARNING: Unhandled error code during point search in KdTree, "
                             "ErrorCode = "
                          << MBresult << " and Coord xyz=" << x << " " << y << " " << z << std::endl;
            }
            std::cout << "No leaf found, MCNP coord xyz=" << x << " " << y << " " << z << std::endl;
            ++cfd_iter;
            continue;
        }

        range.clear();
        MBresult = MBI->get_entities_by_type( tree_node, moab::MBHEX, range );MB_CHK_ERR( MBresult );

        // davg += (double) range.size();
        // if (range.size() > nmax) nmax = range.size();
        // if (range.size() < nmin) nmin = range.size();

        for( moab::Range::iterator rit = range.begin(); rit != range.end(); ++rit )
        {
            verts.clear();
            const moab::EntityHandle* connect;
            int num_connect;
            MBresult = MBI->get_connectivity( *rit, connect, num_connect, true );MB_CHK_ERR( MBresult );

            coords.resize( 3 * num_connect );
            MBresult = MBI->get_coords( connect, num_connect, &coords[0] );MB_CHK_ERR( MBresult );

            for( unsigned int j = 0; j < (unsigned int)num_connect; j++ )
            {
                hexverts[j][0] = coords[3 * j];
                hexverts[j][1] = coords[3 * j + 1];
                hexverts[j][2] = coords[3 * j + 2];
            }

            if( moab::ElemUtil::point_in_trilinear_hex( hexverts, testvc, 1.e-6 ) )
            {
                MBresult = MBI->tag_get_data( MCNP->tally_tag, &( *rit ), 1, &taldata );MB_CHK_ERR( MBresult );
                MBresult = MBI->tag_get_data( MCNP->relerr_tag, &( *rit ), 1, &errdata );MB_CHK_ERR( MBresult );

                outfile << n << "," << testpt[0] << "," << testpt[1] << "," << testpt[2] << "," << taldata << ","
                        << errdata << std::endl;

                if( !read_qnv )
                {
                    MBresult = MBI->tag_set_data( cfd_heating_tag, &( *cfd_iter ), 1, &taldata );MB_CHK_ERR( MBresult );
                    MBresult = MBI->tag_set_data( cfd_error_tag, &( *cfd_iter ), 1, &errdata );MB_CHK_ERR( MBresult );
                }

                found = true;
                elems_read++;

                tal_sum     = tal_sum + taldata;
                err_sum     = err_sum + errdata;
                tal_sum_sqr = tal_sum_sqr + taldata * taldata;
                err_sum_sqr = err_sum_sqr + errdata * errdata;

                break;
            }
        }

        if( !read_qnv ) ++cfd_iter;

        if( !found )
        {
            std::cout << n << " " << testvc << std::endl;
        }
        found = false;
    }

    cadfile.close();
    outfile.close();

    if( result == MCNP_SUCCESS )
    {
        std::cout << "Success! " << elems_read << " elements interpolated." << std::endl << std::endl;
    }
    else
    {
        std::cout << "Failure during query! " << elems_read << " elements interpolated." << std::endl << std::endl;
    }

    double tal_std_dev = sqrt( ( 1.0 / elems_read ) * ( tal_sum_sqr - ( 1.0 / elems_read ) * tal_sum * tal_sum ) );
    double err_std_dev = sqrt( ( 1.0 / elems_read ) * ( err_sum_sqr - ( 1.0 / elems_read ) * err_sum * err_sum ) );

    std::cout << "Tally Mean:               " << tal_sum / elems_read << std::endl;
    std::cout << "Tally Standard Deviation: " << tal_std_dev << std::endl;
    std::cout << "Error Mean:               " << err_sum / elems_read << std::endl;
    std::cout << "Error Standard Deviation: " << err_std_dev << std::endl;

    interp_time = clock() - build_time;

    if( !read_qnv )
    {
        std::string out_mesh_fname = output_filename;
        MBresult                   = MBI->write_mesh( ( out_mesh_fname + ".h5m" ).c_str(), &meshset, 1 );
        // MBresult = MBI->write_file( (cfd_mesh_fname + ".vtk").c_str(), "vtk", NULL, &meshset, 1,
        // &cfd_heating_tag, 1);
    }

    std::cout << "Time to read in file:     " << (double)load_time / CLOCKS_PER_SEC << std::endl;
    std::cout << "Time to build kd-tree:    " << (double)build_time / CLOCKS_PER_SEC << std::endl;
    std::cout << "Time to interpolate data: " << (double)interp_time / CLOCKS_PER_SEC << std::endl;

    return 0;
}

Definition at line 453 of file main.cpp.

{
    static moab::Core inst;
    return &inst;
}

Definition at line 447 of file main.cpp.

{
    static McnpData inst;
    return &inst;
}
MCNPError next_double ( std::string  s,
double &  d,
int &  p 
)

Definition at line 399 of file main.cpp.

References DONE, and MCNP_SUCCESS.

{

    unsigned int slen = s.length();
    unsigned int j;

    for( unsigned int i = p; i < slen; i++ )
    {
        if( ( ( s[i] >= 48 ) && ( s[i] <= 57 ) ) || ( s[i] == 45 ) )
        {

            j = s.find( ",", i );
            if( j > slen ) j = slen;

            d = std::atof( s.substr( i, j - i ).c_str() );
            p = j + 1;

            return MCNP_SUCCESS;
        }
    }

    return DONE;
}
MCNPError next_int ( std::string  s,
int &  k,
int &  p 
)

Definition at line 423 of file main.cpp.

References DONE, and MCNP_SUCCESS.

Referenced by main().

{

    unsigned int slen = s.length();
    unsigned int j;

    for( unsigned int i = p; i < slen; i++ )
    {
        if( ( ( s[i] >= 48 ) && ( s[i] <= 57 ) ) || ( s[i] == 45 ) )
        {

            j = s.find( ",", i );
            if( j > slen ) j = slen;

            k = std::atoi( s.substr( i, j - i ).c_str() );
            p = j + 1;

            return MCNP_SUCCESS;
        }
    }

    return DONE;
}
MCNPError read_files ( int  argc,
char **  argv 
)

Definition at line 358 of file main.cpp.

References CAD_filename, h5m_filename, MCNP, MCNP_FAILURE, output_filename, read_qnv, and skip_build.

Referenced by main().

{
    MCNPError result = MCNP_FAILURE;

    // Check to see if appropriate command lines specified
    if( argc < 3 )
    {
        std::cout << "Source and Target mesh filenames NOT specified!";
        std::cout << std::endl;
        return MCNP_FAILURE;
    }

    // Set the MCNP or H5M filename
    std::string str;
    str = argv[1];

    unsigned int itmp = str.find( ".h5m" );
    if( ( itmp > 0 ) && ( itmp < str.length() ) )
    {
        skip_build   = true;
        h5m_filename = str;
    }
    else
    {
        result = MCNP->set_filename( str );
    }

    // Set the CAD filename
    str          = argv[2];
    CAD_filename = str;

    itmp = str.find( ".qnv" );
    if( ( itmp > 0 ) && ( itmp < str.length() ) ) read_qnv = true;

    // Set the output filename
    str             = argv[3];
    output_filename = str;

    return result;
}

Variable Documentation

clock_t build_time

Definition at line 35 of file main.cpp.

Referenced by main().

std::string CAD_filename

Definition at line 29 of file main.cpp.

Referenced by main(), and read_files().

Definition at line 26 of file main.cpp.

Referenced by main().

Definition at line 26 of file main.cpp.

Referenced by main().

Definition at line 26 of file main.cpp.

Referenced by main().

std::string h5m_filename

Definition at line 28 of file main.cpp.

Referenced by main(), and read_files().

clock_t interp_time

Definition at line 35 of file main.cpp.

Referenced by main().

clock_t load_time

Definition at line 35 of file main.cpp.

Referenced by main().

std::string output_filename

Definition at line 30 of file main.cpp.

Referenced by main(), read_files(), and TestMeshRefiner().

bool read_qnv = false

Definition at line 33 of file main.cpp.

Referenced by main(), and read_files().

Definition at line 26 of file main.cpp.

Referenced by main().

bool skip_build = false

Definition at line 32 of file main.cpp.

Referenced by main(), and read_files().

clock_t start_time

Definition at line 35 of file main.cpp.

Referenced by main().

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines