MOAB: Mesh Oriented datABase
(version 5.4.1)
|
#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"
Go to the source code of this file.
Defines | |
#define | MBI mb_instance() |
Functions | |
McnpData * | mc_instance () |
moab::Interface * | mb_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 MBI mb_instance() |
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, outfile, 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; }
moab::Interface * mb_instance | ( | ) |
Definition at line 453 of file main.cpp.
{ static moab::Core inst; return &inst; }
McnpData * mc_instance | ( | ) |
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; }
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; }
clock_t build_time |
Definition at line 35 of file main.cpp.
Referenced by do_kdtree_test(), do_linear_test(), and 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(), and test_global_value_common().
std::string h5m_filename |
Definition at line 28 of file main.cpp.
Referenced by main(), and read_files().
clock_t interp_time |
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().
bool skip_build = false |
Definition at line 32 of file main.cpp.
Referenced by main(), and read_files().
clock_t start_time |