[MOAB-dev] r1843 - in MOAB/trunk: . tools/mcnpmit
wilsonp at mcs.anl.gov
wilsonp at mcs.anl.gov
Mon May 26 10:07:56 CDT 2008
Author: wilsonp
Date: 2008-05-26 10:07:55 -0500 (Mon, 26 May 2008)
New Revision: 1843
Modified:
MOAB/trunk/ReadIDEAS.cpp
MOAB/trunk/ReadIDEAS.hpp
MOAB/trunk/tools/mcnpmit/main.cpp
MOAB/trunk/tools/mcnpmit/mcnpmit.cpp
Log:
Bug-fix for ReadIDEAS: remove off-by-one error on test for empty string
Feature extension for ReadIDEAS: read material and physical property table data and tag elements with integer
Bug-fix in mcnpmit: update tags to be typed: MB_TYPE_DOUBLE
Feature extension in mcnpmit:
Instead of following unique behavior only for .unv (IDEAS Unviersal format), follow unique
behavior only for .qnv (centroid data). This allows the use of any file format supported by MOAB
to specifiy CFD mesh.
Bug-fix in mcnpmit: reverse order of loops when reading MCNP data
Modified: MOAB/trunk/ReadIDEAS.cpp
===================================================================
--- MOAB/trunk/ReadIDEAS.cpp 2008-05-23 18:57:33 UTC (rev 1842)
+++ MOAB/trunk/ReadIDEAS.cpp 2008-05-26 15:07:55 UTC (rev 1843)
@@ -79,7 +79,7 @@
il = std::strtol(line, &ctmp, 10);
if (il == -1) {
- s = ctmp+1;
+ s = ctmp;
if (s.empty()) end_of_block++;
}
else end_of_block = 0;
@@ -91,6 +91,7 @@
}
+
void ReadIDEAS::create_vertices() {
// Read two lines: first has some data, second has coordinates
@@ -114,8 +115,8 @@
il1 = std::strtol(line1, &ctmp1, 10);
il2 = std::strtol(line2, &ctmp2, 10);
if ((il1 == -1) && (il2 == -1)) {
- s1 = ctmp1+1;
- s2 = ctmp2+1;
+ s1 = ctmp1;
+ s2 = ctmp2;
if ((s1.empty()) && (s2.empty())) break;
}
num_verts++;
@@ -163,6 +164,10 @@
MBRange verts;
rval = MBI->get_entities_by_type( mesh_handle, MBVERTEX, verts, true);
MBEntityHandle vstart = *(verts.begin());
+
+ MBTag mat_prop_tag, phys_prop_tag;
+ rval = MBI->tag_create( MAT_PROP_TABLE_TAG , sizeof(int), MB_TAG_DENSE, mat_prop_tag, 0);
+ rval = MBI->tag_create( PHYS_PROP_TABLE_TAG , sizeof(int), MB_TAG_DENSE, phys_prop_tag, 0);
while (! file.eof() ) {
@@ -173,11 +178,15 @@
il1 = std::strtol(line1, &ctmp1, 10);
il2 = std::strtol(line2, &ctmp2, 10);
if ((il1 == -1) && (il2 == -1)) {
- s1 = ctmp1+1;
- s2 = ctmp2+1;
+ s1 = ctmp1;
+ s2 = ctmp2;
if ((s1.empty()) && (s2.empty())) break;
}
+ // Get property tables out of 1st line
+ int phys_table = strtol(line1+21, &ctmp1, 10);
+ int mat_table = strtol(line1+31, &ctmp1, 10);
+
// Get the connectivity out of the 2nd line
connect[0] = vstart + strtol( line2, &ctmp2, 10) - 1;
connect[1] = vstart + strtol( ctmp2+1, &ctmp2, 10) - 1;
@@ -187,10 +196,16 @@
// Make the element
rval = MBI->create_element(MBTET, connect, 4, handle);
assert( MB_SUCCESS == rval );
-
+
+ rval = MBI->tag_set_data(phys_prop_tag,&handle,1,&phys_table);
+ assert( MB_SUCCESS == rval);
+
+ rval = MBI->tag_set_data(mat_prop_tag,&handle,1,&mat_table);
+ assert( MB_SUCCESS == rval);
+
rval = MBI->add_entities( mesh_handle, &handle, 1);
assert( MB_SUCCESS == rval );
-
+
}
}
Modified: MOAB/trunk/ReadIDEAS.hpp
===================================================================
--- MOAB/trunk/ReadIDEAS.hpp 2008-05-23 18:57:33 UTC (rev 1842)
+++ MOAB/trunk/ReadIDEAS.hpp 2008-05-26 15:07:55 UTC (rev 1843)
@@ -8,6 +8,9 @@
#define VERTEX_LIST 2411
#define MAKE_TETRAHEDRA 2412
+#define MAT_PROP_TABLE_TAG "mat_prop_table"
+#define PHYS_PROP_TABLE_TAG "phys_prop_table"
+
class MBReadUtilIface;
class ReadIDEAS : public MBReaderIface
@@ -17,12 +20,12 @@
static MBReaderIface* factory( MBInterface* );
-MBErrorCode load_file(const char* fname,
- MBEntityHandle& meshset,
- const FileOptions&,
- const int* material_set_list,
- int num_material_sets );
-
+ MBErrorCode load_file(const char* fname,
+ MBEntityHandle& meshset,
+ const FileOptions&,
+ const int* material_set_list,
+ int num_material_sets );
+
//! Constructor
ReadIDEAS(MBInterface* impl = NULL);
@@ -30,21 +33,21 @@
virtual ~ReadIDEAS() {}
protected:
-
+
void skip_header();
void create_vertices();
void create_tetrahedral_elements();
-
+
private:
-
+
std::ifstream file;
-
+
// Read mesh interface
MBReadUtilIface* readMeshIface;
-
+
// MOAB Interface
MBInterface* MBI;
-
+
// Handle for the mesh
MBEntityHandle mesh_handle;
Modified: MOAB/trunk/tools/mcnpmit/main.cpp
===================================================================
--- MOAB/trunk/tools/mcnpmit/main.cpp 2008-05-23 18:57:33 UTC (rev 1842)
+++ MOAB/trunk/tools/mcnpmit/main.cpp 2008-05-26 15:07:55 UTC (rev 1843)
@@ -29,7 +29,7 @@
std::string output_filename;
bool skip_build = false;
-bool read_unv = false;
+bool read_qnv = false;
clock_t start_time, load_time, build_time, interp_time;
@@ -45,12 +45,14 @@
// 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;
@@ -73,6 +75,10 @@
}
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 (MB_SUCCESS == MBI->get_last_error(message) && !message.empty())
+ std::cerr << "Error message: " << message << std::endl;
return 1;
}
@@ -96,7 +102,12 @@
}
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 (MB_SUCCESS == MBI->get_last_error(message) && !message.empty())
+ std::cerr << "Error message: " << message << std::endl;
return 1;
+ return 1;
}
}
@@ -123,12 +134,12 @@
std::string s;
char line[10000];
- // Used only when reading a .unv file to get vertex info
+ // Used only when reading a mesh file to get vertex info
double *cfd_coords;
MBRange::iterator cfd_iter;
MBEntityHandle meshset;
- if (!read_unv) {
+ if (read_qnv) {
cadfile.open( CAD_filename.c_str() );
cadfile.getline(line, 10000);
cadfile.getline(line, 10000);
@@ -147,7 +158,7 @@
MBresult = MBI->get_coords( cfd_verts , cfd_coords );
cfd_iter = cfd_verts.begin();
- MBresult = MBI->tag_create("heating_tag", sizeof(double), MB_TAG_DENSE, cfd_heating_tag, 0);
+ MBresult = MBI->tag_create("heating_tag", sizeof(double), MB_TAG_DENSE, MB_TYPE_DOUBLE, cfd_heating_tag, 0);
std::cout << std::endl << "Read in mesh with query points." << std::endl << std::endl;
@@ -176,10 +187,15 @@
// double davg = 0.0;
// unsigned int nmax = 0, nmin = 1000000000 ;
+ int status_freq = int(num_pts/100);
+
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_unv) {
+ if (read_qnv) {
cadfile.getline(line, 10000);
nl = std::strtol(line, &ctmp, 10); n = (unsigned int) nl;
@@ -239,7 +255,7 @@
<< testpt[2] << ","
<< taldata << std::endl;
- if (read_unv) {
+ if (!read_qnv) {
MBresult = MBI->tag_set_data(cfd_heating_tag, &(*cfd_iter), 1, &taldata);
cfd_iter++;
}
@@ -270,10 +286,9 @@
interp_time = clock() - build_time;
- if (read_unv) {
- std::string cfd_mesh_fname = CAD_filename;
- cfd_mesh_fname.erase( CAD_filename.find(".unv"), CAD_filename.length());
- MBresult = MBI->write_mesh( (cfd_mesh_fname + ".h5m").c_str(), &meshset, 1);
+ 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);
}
@@ -313,9 +328,9 @@
str = argv[2];
CAD_filename = str;
- itmp = str.find(".unv");
+ itmp = str.find(".qnv");
if ((itmp > 0) && (itmp < str.length()))
- read_unv = true;
+ read_qnv = true;
// Set the output filename
str = argv[3];
Modified: MOAB/trunk/tools/mcnpmit/mcnpmit.cpp
===================================================================
--- MOAB/trunk/tools/mcnpmit/mcnpmit.cpp 2008-05-23 18:57:33 UTC (rev 1842)
+++ MOAB/trunk/tools/mcnpmit/mcnpmit.cpp 2008-05-26 15:07:55 UTC (rev 1843)
@@ -152,9 +152,9 @@
// vstart = MCNP_vertices.front();
vstart = *(vert_handles.begin());
+ for (int i=0; i < nv[0]-1; i++) {
+ for (int j=0; j < nv[1]-1; j++) {
for (int k=0; k < nv[2]-1; k++) {
- for (int j=0; j < nv[1]-1; j++) {
- for (int i=0; i < nv[0]-1; i++) {
vijk = vstart + (i + j*nv[0] + k*nv[0]*nv[1]);
//std::cout << vijk << std::endl;
@@ -312,8 +312,8 @@
MBErrorCode rval;
- rval = MBI->tag_create(TALLY_TAG, sizeof(double), MB_TAG_DENSE, tally_tag, 0);
- rval = MBI->tag_create(ERROR_TAG, sizeof(double), MB_TAG_DENSE, relerr_tag, 0);
+ rval = MBI->tag_create(TALLY_TAG, sizeof(double), MB_TAG_DENSE, MB_TYPE_DOUBLE, tally_tag, 0);
+ rval = MBI->tag_create(ERROR_TAG, sizeof(double), MB_TAG_DENSE, MB_TYPE_DOUBLE, relerr_tag, 0);
return MCNP_SUCCESS;
More information about the moab-dev
mailing list