[MOAB-dev] r3006 - MOAB/trunk/tools/dagmc
wilsonp at engr.wisc.edu
wilsonp at engr.wisc.edu
Tue Jul 14 11:49:39 CDT 2009
Author: wilsonp
Date: 2009-07-14 11:49:39 -0500 (Tue, 14 Jul 2009)
New Revision: 3006
Modified:
MOAB/trunk/tools/dagmc/DagMC.cpp
MOAB/trunk/tools/dagmc/DagMC.hpp
MOAB/trunk/tools/dagmc/main.cc
Log:
MAJOR CHANGE to DagMC interface and user interaction
* replace cgm2moab function with ReadCGM in both DagMC & cgm2moab app
* allow default values for tags in get_tag
* separate load_file_and_init into two functions in preparation for
parallel version
* load_file (now based on ReadCGM)
* init_OBBTree
* separate tight connection with MCNP in code base
* rename and cleanup DagMC settings
* separate write_log into two functions
* parse_metadata
* write_mcnp (this can be replicated for other codes)
* change parsing of group names
* change some keywords to include '.' instead of '_' so that
'_' can be used as delimiter for tokenizer
Modified: MOAB/trunk/tools/dagmc/DagMC.cpp
===================================================================
--- MOAB/trunk/tools/dagmc/DagMC.cpp 2009-07-14 14:57:15 UTC (rev 3005)
+++ MOAB/trunk/tools/dagmc/DagMC.cpp 2009-07-14 16:49:39 UTC (rev 3006)
@@ -4,6 +4,7 @@
#include "MBRange.hpp"
#include "MBCore.hpp"
#include "MBGeomUtil.hpp"
+#include "FileOptions.hpp"
#ifdef CGM
#include "InitCGMA.hpp"
@@ -31,11 +32,25 @@
#define MB_OBB_TREE_TAG_NAME "OBB_TREE"
#define CATEGORY_TAG_LENGTH 32
+#define MAT_GROUP 0
+#define BC_SPEC 1
+#define BC_WHITE 2
+#define IMP_ZERO 3
+#define TALLY_GROUP 4
+
+#define MATERIAL_TAG_NAME "DAGMC_MATERIAL_ID"
+#define DENSITY_TAG_NAME "DAGMC_MATERIAL_DENSITY"
+#define BC_TAG_NAME "DAGMC_BOUNDARY_CONDITION"
+#define IMP_TAG_NAME "DAGMC_IMPORTANCE"
+#define TALLY_TAG_NAME "DAGMC_TALLY"
+
+
+#define MBI moab_instance()
+
DagMC *DagMC::instance_ = NULL;
const bool debug = false;
-
unsigned int DagMC::interface_revision()
{
unsigned int result = 0;
@@ -44,13 +59,17 @@
}
+void DagMC::create_instance(MBInterface *mb_impl)
+{
+ if (NULL == mb_impl) mb_impl = new MBCore();
+ instance_ = new DagMC(mb_impl);
+}
+
+
DagMC::DagMC(MBInterface *mb_impl)
- : mbImpl(mb_impl), obbTree(mb_impl),
- discardDistTol(1e-8),
- addDistTol(1e-6),
- facetingTolerance(0.001),
- moabMCNPSourceCell(0), moabMCNPUseDistLimit(false)
+ : mbImpl(mb_impl), obbTree(mb_impl), is_geom(false)
{
+ // This is the correct place to uniquely define default values for the dagmc settings
options[0] = Option( "source_cell", "source cell ID, or zero if unknown", "0" );
options[1] = Option( "discard_distance_tolerance", "positive real value", "0.001" );
options[2] = Option( "use_distance_limit", "one to enable distance limit optimization, zero otherwise", "0" );
@@ -58,6 +77,9 @@
options[4] = Option( "faceting_tolerance", "graphics faceting tolerance", "0.001" );
options[5] = Option( "add_distance_tolerance", "positive real value", "0.001" );
+ // call parse settings to initialize default values for settings from options
+ parse_settings();
+
memset( specReflectName, 0, NAME_TAG_SIZE );
strcpy( specReflectName, "spec_reflect" );
memset( whiteReflectName, 0, NAME_TAG_SIZE );
@@ -67,7 +89,7 @@
distanceLimit = std::numeric_limits<double>::max();
- nameTag = get_tag(NAME_TAG_NAME, NAME_TAG_SIZE, MB_TAG_SPARSE, MB_TYPE_OPAQUE, false);
+ nameTag = get_tag(NAME_TAG_NAME, NAME_TAG_SIZE, MB_TAG_SPARSE, MB_TYPE_OPAQUE, NULL, false);
idTag = get_tag( GLOBAL_ID_TAG_NAME, sizeof(int), MB_TAG_DENSE, MB_TYPE_INTEGER );
@@ -78,6 +100,17 @@
// get sense of surfaces wrt volumes
senseTag = get_tag( "GEOM_SENSE_2", 2*sizeof(MBEntityHandle), MB_TAG_DENSE, MB_TYPE_HANDLE );
+ int matid = 0;
+ const void *def_matid = &matid;
+ matTag = get_tag(MATERIAL_TAG_NAME, sizeof(int), MB_TAG_DENSE, MB_TYPE_INTEGER, def_matid );
+ densTag = get_tag(DENSITY_TAG_NAME, sizeof(double), MB_TAG_DENSE, MB_TYPE_DOUBLE );
+ bcTag = get_tag(BC_TAG_NAME, sizeof(int), MB_TAG_SPARSE, MB_TYPE_INTEGER);
+
+ double imp_one = 1;
+ const void *def_imp = &imp_one;
+ impTag = get_tag(IMP_TAG_NAME, sizeof(double), MB_TAG_DENSE, MB_TYPE_DOUBLE, def_imp );
+ tallyTag = get_tag(TALLY_TAG_NAME, sizeof(int), MB_TAG_SPARSE, MB_TYPE_INTEGER);
+
}
MBErrorCode DagMC::ray_fire(const MBEntityHandle vol, const MBEntityHandle last_surf_hit,
@@ -183,24 +216,80 @@
return MB_SUCCESS;
}
-void DagMC::create_instance(MBInterface *mb_impl)
-{
- if (NULL == mb_impl) mb_impl = new MBCore();
- instance_ = new DagMC(mb_impl);
-}
+// **** HANDLING FILES TO INPUT DAGMC SETTINGS ****
+// there is a movement away from using a settings file and, instead, incorporating this
+// information into the native input file of the calling application
-void DagMC::write_settings( FILE* fptr, bool desc ) {
+void DagMC::read_settings( const char* filename )
+{
int num_opt = sizeof(options) / sizeof(options[0]);
- if (desc) for (int i = 0; i < num_opt; ++i)
- fprintf( fptr, "%s = %s # %s\n", options[i].name.c_str(), options[i].value.c_str(), options[i].desc.c_str() );
- else for (int i = 0; i < num_opt; ++i)
- fprintf( fptr, "%s = %s\n", options[i].name.c_str(), options[i].value.c_str() );
+ FILE* file;
+
+ if (filename && (file = fopen( filename, "r" ))) {
+ int line = 0;
+ char buffer[256];
+ bool havenl = true;
+ while ( fgets( buffer, sizeof(buffer), file ) ) {
+ if (havenl) {
+ havenl = false;
+ ++line;
+ }
+ int len = strlen(buffer);
+ if (len && buffer[len-1] == '\n')
+ havenl = true;
+
+ // chop off any thing after the '#' comman indicator
+ char* c = strchr( buffer, '#' );
+ if (c) *c = '\0';
+ // skip leading white space
+ char* p = buffer;
+ while (isspace(*p)) ++p;
+ // if empty line, done
+ if (!*p)
+ continue;
+ // find '='
+ c = strchr( p, '=' );
+ if (!c) {
+ fprintf(stderr, "Invalid option at line %d of config file '%s'\n", line, filename );
+ exit(2);
+ }
+ // get rid of white space around '='
+ char* v = c+1;
+ *c = ' ';
+ while (c > p && isspace(*c)) {
+ *c = '\0';
+ --c;
+ }
+ while (isspace(*v))
+ ++v;
+
+ // search for option
+ bool found = false;
+ for (int i = 0; i < num_opt; ++i) {
+ if (options[i].name == p) {
+ found = true;
+ options[i].value = v;
+ options[i].user_set = true;
+ }
+ }
+ if (!found)
+ fprintf(stderr,"Warning: unknown option at line %d of '%s': %s\n", line, filename, p );
+ } // while( fgets() )
+
+ fclose(file);
+ } // file(file)
+
+ // always do this, even if no file was found.
+ // it will either read the default values, or re-parse
+ // the values previously read from some other file.
+ parse_settings();
}
+ // this method sets the default values from options[] if none were read from file
void DagMC::parse_settings() {
- moabMCNPSourceCell = atoi( options[0].value.c_str() );
- if (moabMCNPSourceCell < 0) {
- std::cerr << "Invalid source_cell = " << moabMCNPSourceCell << std::endl;
+ sourceCell = atoi( options[0].value.c_str() );
+ if (sourceCell < 0) {
+ std::cerr << "Invalid source_cell = " << sourceCell << std::endl;
exit(2);
}
@@ -216,7 +305,7 @@
exit(2);
}
- moabMCNPUseDistLimit = !!atoi( options[2].value.c_str() );
+ useDistLimit = !!atoi( options[2].value.c_str() );
useCAD = !!atoi( options[3].value.c_str() );
#ifndef HAVE_CGM_FIRE_RAY
@@ -235,16 +324,28 @@
}
}
+
+void DagMC::write_settings( FILE* fptr, bool desc ) {
+ int num_opt = sizeof(options) / sizeof(options[0]);
+ if (desc) for (int i = 0; i < num_opt; ++i)
+ fprintf( fptr, "%s = %s # %s\n", options[i].name.c_str(), options[i].value.c_str(), options[i].desc.c_str() );
+ else for (int i = 0; i < num_opt; ++i)
+ fprintf( fptr, "%s = %s\n", options[i].name.c_str(), options[i].value.c_str() );
+}
+
+// **** PASS SETTINGS FROM CALLING APPLICATION ****
+// This is the prefered way to set DAGMC settings
+
void DagMC::set_settings(int source_cell, int use_cad, int use_dist_limit,
double add_distance_tolerance,
double discard_distance_tolerance) {
- moabMCNPSourceCell = source_cell;
- if (moabMCNPSourceCell < 0) {
- std::cerr << "Invalid source_cell = " << moabMCNPSourceCell << std::endl;
+ sourceCell = source_cell;
+ if (sourceCell < 0) {
+ std::cerr << "Invalid source_cell = " << sourceCell << std::endl;
exit(2);
}
- std::cout << "Set Source Cell = " << moabMCNPSourceCell << std::endl;
+ std::cout << "Set Source Cell = " << sourceCell << std::endl;
addDistTol = add_distance_tolerance;
if (addDistTol <= 0 || addDistTol > 1) {
@@ -262,9 +363,9 @@
std::cout << "Set distance tolerance 2 = " << discardDistTol << std::endl;
- moabMCNPUseDistLimit = !!(use_dist_limit);
+ useDistLimit = !!(use_dist_limit);
- std::cout << "Turned " << (moabMCNPUseDistLimit?"ON":"OFF") << " distance limit." << std::endl;
+ std::cout << "Turned " << (useDistLimit?"ON":"OFF") << " distance limit." << std::endl;
useCAD = !!(use_cad);
std::cout << "Turned " << (useCAD?"ON":"OFF") << " ray firing on full CAD model." << std::endl;
@@ -280,71 +381,8 @@
}
-void DagMC::read_settings( const char* filename )
-{
- int num_opt = sizeof(options) / sizeof(options[0]);
- FILE* file;
-
- if (filename && (file = fopen( filename, "r" ))) {
- int line = 0;
- char buffer[256];
- bool havenl = true;
- while ( fgets( buffer, sizeof(buffer), file ) ) {
- if (havenl) {
- havenl = false;
- ++line;
- }
- int len = strlen(buffer);
- if (len && buffer[len-1] == '\n')
- havenl = true;
-
- // chop off any thing after the '#' comman indicator
- char* c = strchr( buffer, '#' );
- if (c) *c = '\0';
- // skip leading white space
- char* p = buffer;
- while (isspace(*p)) ++p;
- // if empty line, done
- if (!*p)
- continue;
- // find '='
- c = strchr( p, '=' );
- if (!c) {
- fprintf(stderr, "Invalid option at line %d of config file '%s'\n", line, filename );
- exit(2);
- }
- // get rid of white space around '='
- char* v = c+1;
- *c = ' ';
- while (c > p && isspace(*c)) {
- *c = '\0';
- --c;
- }
- while (isspace(*v))
- ++v;
-
- // search for option
- bool found = false;
- for (int i = 0; i < num_opt; ++i) {
- if (options[i].name == p) {
- found = true;
- options[i].value = v;
- options[i].user_set = true;
- }
- }
- if (!found)
- fprintf(stderr,"Warning: unknown option at line %d of '%s': %s\n", line, filename, p );
- } // while( fgets() )
-
- fclose(file);
- } // file(file)
-
- // always do this, even if no file was found.
- // it will either read the default values, or re-parse
- // the values previously read from some other file.
- parse_settings();
-}
+
// point_in_volume_slow, including poly_solid_angle helper subroutine
// are adapted from "Point in Polyhedron Testing Using Spherical Polygons", Paulo Cezar
// Pinto Carvalho and Paulo Roma Cavalcanti, _Graphics Gems V_, pg. 42. Original algorithm
@@ -360,7 +398,7 @@
// Get connectivity
const MBEntityHandle* conn;
int len;
- rval = moab_instance()->get_connectivity( face, conn, len, true );
+ rval = MBI->get_connectivity( face, conn, len, true );
if (MB_SUCCESS != rval)
return rval;
@@ -374,7 +412,7 @@
}
// get coordinates
- rval = moab_instance()->get_coords( conn, len, coords->array() );
+ rval = MBI->get_coords( conn, len, coords->array() );
if (MB_SUCCESS != rval)
return rval;
@@ -421,7 +459,7 @@
double sum = 0.0;
const MBCartVect point(x,y,z);
- rval = moab_instance()->get_child_meshsets( volume, surfs );
+ rval = MBI->get_child_meshsets( volume, surfs );
if (MB_SUCCESS != rval)
return rval;
@@ -436,7 +474,7 @@
double surf_area = 0.0, face_area;
faces.clear();
- rval = moab_instance()->get_entities_by_dimension( surfs[i], 2, faces );
+ rval = MBI->get_entities_by_dimension( surfs[i], 2, faces );
if (MB_SUCCESS != rval)
return rval;
@@ -562,10 +600,10 @@
int len, sense;
MBCartVect coords[3], normal;
if (tris.size() == 1) {
- rval = moab_instance()->get_connectivity( tris.front(), conn, len );
+ rval = MBI->get_connectivity( tris.front(), conn, len );
if (MB_SUCCESS != rval || len != 3)
return MB_FAILURE;
- rval = moab_instance()->get_coords( conn, len, coords[0].array() );
+ rval = MBI->get_coords( conn, len, coords[0].array() );
if (MB_SUCCESS != rval) return rval;
rval = surface_sense( volume, surfs.front(), sense );
@@ -593,10 +631,10 @@
for (unsigned i = 0; i < tris.size(); ++i) {
some_above = some_below = false;
- rval = moab_instance()->get_connectivity( tris[i], conn, len );
+ rval = MBI->get_connectivity( tris[i], conn, len );
if (MB_SUCCESS != rval || len != 3)
return MB_FAILURE;
- rval = moab_instance()->get_coords( conn, len, coords[0].array() );
+ rval = MBI->get_coords( conn, len, coords[0].array() );
if (MB_SUCCESS != rval) return rval;
rval = surface_sense( volume, surfs[i], sense );
@@ -613,10 +651,10 @@
if (j == i)
continue;
- rval = moab_instance()->get_connectivity( tris[j], conn, len );
+ rval = MBI->get_connectivity( tris[j], conn, len );
if (MB_SUCCESS != rval || len != 3)
return MB_FAILURE;
- rval = moab_instance()->get_coords( conn, len, coords[0].array() );
+ rval = MBI->get_coords( conn, len, coords[0].array() );
if (MB_SUCCESS != rval)
return rval;
@@ -695,7 +733,7 @@
}
// get surfaces from volume
- rval = moab_instance()->get_child_meshsets( volume, surfaces );
+ rval = MBI->get_child_meshsets( volume, surfaces );
if (MB_SUCCESS != rval) return rval;
// get surface senses
@@ -714,7 +752,7 @@
// get triangles in surface
MBRange triangles;
- rval = moab_instance()->get_entities_by_dimension( surfaces[i], 2, triangles );
+ rval = MBI->get_entities_by_dimension( surfaces[i], 2, triangles );
if (MB_SUCCESS != rval)
return rval;
if (!triangles.all_of_type(MBTRI)) {
@@ -728,10 +766,10 @@
int len;
MBCartVect coords[3];
for (MBRange::iterator j = triangles.begin(); j != triangles.end(); ++j) {
- rval = moab_instance()->get_connectivity( *j, conn, len, true );
+ rval = MBI->get_connectivity( *j, conn, len, true );
if (MB_SUCCESS != rval) return rval;
assert(3 == len);
- rval = moab_instance()->get_coords( conn, 3, coords[0].array() );
+ rval = MBI->get_coords( conn, 3, coords[0].array() );
if (MB_SUCCESS != rval) return rval;
coords[1] -= coords[0];
@@ -750,7 +788,7 @@
{
// get triangles in surface
MBRange triangles;
- MBErrorCode rval = moab_instance()->get_entities_by_dimension( surface, 2, triangles );
+ MBErrorCode rval = MBI->get_entities_by_dimension( surface, 2, triangles );
if (MB_SUCCESS != rval)
return rval;
if (!triangles.all_of_type(MBTRI)) {
@@ -764,10 +802,10 @@
int len;
MBCartVect coords[3];
for (MBRange::iterator j = triangles.begin(); j != triangles.end(); ++j) {
- rval = moab_instance()->get_connectivity( *j, conn, len, true );
+ rval = MBI->get_connectivity( *j, conn, len, true );
if (MB_SUCCESS != rval) return rval;
assert(3 == len);
- rval = moab_instance()->get_coords( conn, 3, coords[0].array() );
+ rval = MBI->get_coords( conn, 3, coords[0].array() );
if (MB_SUCCESS != rval) return rval;
coords[1] -= coords[0];
@@ -793,7 +831,7 @@
volume = (MBEntityHandle) 0;
std::vector<MBEntityHandle> surf_volumes( 2*num_surfaces );
- MBErrorCode rval = moab_instance()->tag_get_data( sense_tag(), surfaces, num_surfaces, &surf_volumes[0] );
+ MBErrorCode rval = MBI->tag_get_data( sense_tag(), surfaces, num_surfaces, &surf_volumes[0] );
if (MB_SUCCESS != rval) return rval;
const MBEntityHandle* end = surfaces + num_surfaces;
@@ -828,7 +866,7 @@
// get sense of surfaces wrt volumes
MBEntityHandle surf_volumes[2];
- MBErrorCode rval = moab_instance()->tag_get_data( sense_tag(), &surface, 1, surf_volumes );
+ MBErrorCode rval = MBI->tag_get_data( sense_tag(), &surface, 1, surf_volumes );
if (MB_SUCCESS != rval) return rval;
if (surf_volumes[0] == volume)
@@ -841,93 +879,56 @@
return MB_SUCCESS;
}
-MBErrorCode DagMC::load_file_and_init(const char* cfile,
- const int clen,
- const char* ffile,
- const int flen,
- const double facet_tolerance)
+MBErrorCode DagMC::load_file(const char* cfile,
+ const double facet_tolerance)
{
- // Always do this first to make sure we have at least the
- // default values if noone ever calls read_settings(). If
- // read_settings() has already been called, it will just
- // re-parse the same values.
- read_settings(0);
-
+ MBErrorCode rval;
+ MBEntityHandle file_set;
+
+ // override default value of facetingTolerance with passed value
if (facet_tolerance > 0 )
facetingTolerance = facet_tolerance;
- MBErrorCode rval;
-
- std::string scfile = cfile;
+ char facetTolStr[16];
- // save whether it's geometry or mesh, building obb later depends
- // on this setting
- bool is_geom = false;
-
- if (scfile.find(".sat",0) > (unsigned int) 0 &&
- scfile.find(".sat",0) <= scfile.length()) {
+ sprintf(facetTolStr,"%g",facetingTolerance);
- is_geom = true;
+ // what if we are using default faceting tolerance???
+ char options[120] = "CGM_ATTRIBS=yes;FACET_DISTANCE_TOLERANCE=";
+ strcat(options,facetTolStr);
-#ifdef CGM
- // initialize cgm
- InitCGMA::initialize_cgma("ACIS");
- CGMApp::instance()->attrib_manager()->set_all_auto_read_flags(true);
- CGMApp::instance()->attrib_manager()->set_all_auto_actuate_flags(true);
- CubitStatus s = GeometryQueryTool::instance()->import_solid_model( scfile.c_str(), "ACIS_SAT");
- if (CUBIT_SUCCESS != s) {
- std::cerr << "Failed to read '" << scfile << "' of type 'ACIS_SAT'"
- << std::endl;
- exit(2);
- }
- if (!cgm2moab(moab_instance(), facetingTolerance)) {
- std::cerr << "Internal error copying geometry" << std::endl;
- exit(5);
- }
-
- if (!useCAD) {
- // ok, we're done with CGM
- CGMApp::instance()->shutdown();
- }
-#else
- std::cerr << "CGM not supported in this version." << std::endl;
-#endif
+ rval = MBI->load_file(cfile, file_set, options, NULL, 0, 0);
+ if (MB_SUCCESS != rval) {
+ std::cerr << "Couldn't read file " << cfile << std::endl;
+ return rval;
}
- else if (scfile.find(".cub",0) > (unsigned int) 0 &&
- scfile.find(".cub",0) <= scfile.length()) {
- std::cerr << "Error: can't read .cub files yet." << std::endl;
- exit(7);
- }
- else {
- rval = moab_instance()->load_mesh(scfile.c_str());
- if (MB_SUCCESS != rval) {
- std::cerr << "Couldn't read file " << scfile << std::endl;
- return rval;
- }
- else if (useCAD) {
- std::cerr << "useCAD parameter turned on but geometry read from "
- << "a mesh file;" << std::endl
- << " turning that parameter off for you." << std::endl;
- useCAD = 0;
- }
- else if (options[4].user_set == true) {
- std::cerr << "Warning: user-set tolerance won't apply to pre-generated mesh file."
- << std::endl;
- }
- }
+
+ // How do we know what type of file we read for downstream issues:
+ // * useCad
+ // * is_geom in build_indices
+
+ return MB_SUCCESS;
+
+}
+
+MBErrorCode DagMC::init_OBBTree()
+{
+
+ MBErrorCode rval;
+
MBRange surfs, vols;
const int three = 3;
const void* const three_val[] = {&three};
- rval =
- moab_instance()->get_entities_by_type_and_tag( 0, MBENTITYSET, &geomTag,
- three_val, 1, vols );
+ rval = MBI->get_entities_by_type_and_tag( 0, MBENTITYSET, &geomTag,
+ three_val, 1, vols );
if (MB_SUCCESS != rval)
return rval;
const int two = 2;
const void* const two_val[] = {&two};
- rval = moab_instance()->get_entities_by_type_and_tag( 0, MBENTITYSET, &geomTag, two_val, 1, surfs );
+ rval = MBI->get_entities_by_type_and_tag( 0, MBENTITYSET, &geomTag,
+ two_val, 1, surfs );
if (MB_SUCCESS != rval)
return rval;
@@ -950,23 +951,32 @@
}
// build the various index vectors used for efficiency
- rval = build_indices(surfs, vols, is_geom);
+ rval = build_indices(surfs, vols);
if (MB_SUCCESS != rval) {
std::cerr << "Failed to build surface/volume indices." << std::endl;
return rval;
}
- // finally, if user entered a facet file name, export facets to
- // that file
+ return MB_SUCCESS;
+}
+
+
+MBErrorCode DagMC::write_mesh(const char* ffile,
+ const int flen)
+{
+ MBErrorCode rval;
+
+ // write out a mesh file if requested
if (ffile && 0 < flen) {
- rval = moab_instance()->write_mesh(ffile);
+ rval = MBI->write_mesh(ffile);
if (MB_SUCCESS != rval) {
std::cerr << "Failed to write mesh to " << ffile << "." << std::endl;
return rval;
}
}
-
+
return MB_SUCCESS;
+
}
bool DagMC::have_obb_tree()
@@ -986,7 +996,7 @@
MBErrorCode rval;
MBRange results;
- rval = moab_instance()->get_entities_by_type_and_tag( 0, MBENTITYSET, tags, vals, 2, results );
+ rval = MBI->get_entities_by_type_and_tag( 0, MBENTITYSET, tags, vals, 2, results );
if (MB_SUCCESS != rval || results.empty())
return 0;
@@ -1000,356 +1010,337 @@
return 0;
int result = 0;
- moab_instance()->tag_get_data( idTag, &h, 1, &result );
+ MBI->tag_get_data( idTag, &h, 1, &result );
return result;
}
-void DagMC::write_log(std::string ifile, const bool overwrite)
+// functionality in write_log
+//
+// 1. for each group
+// a) check for material assignment -> get list of member vols and set mat/dens array
+// b) check for b.c. -> get list of surfaces and set bc array (incl. periodic???)
+// c) check for graveyard -> set importance
+// d) check for tallies
+// * record tally ID
+// 2. write cell cards
+// 3. write surf cards
+// 4. if tallies
+// a) loop through all groups to find tallies
+// * parse tally type, etc
+// * get list of cells or surfaces
+// * write tally cards
+
+
+
+MBErrorCode DagMC::write_mcnp(std::string ifile, const bool overwrite)
{
+ std::map<std::string,int> tallyKeywords;
+
+ tallyKeywords["surf.current"] = 1;
+ tallyKeywords["surf.flux"] = 2;
+ tallyKeywords["cell.flux"] = 4;
+ tallyKeywords["cell.heating"] = 6;
+ tallyKeywords["cell.fission"] = 7;
+ tallyKeywords["pulse.height"] = 8;
+
+ std::vector<MBEntityHandle>::iterator iter;
+
if (!overwrite) {
// if not overwriting, test for file, and if it exists, just return;
// this replaces the old inquir function
std::ifstream testfile;
testfile.open(ifile.c_str());
- if (!testfile.fail()) return;
+ if (!testfile.fail()) return MB_SUCCESS;
}
const char *cifile = ifile.c_str();
std::ofstream cgmfile;
cgmfile.open(cifile);
- int ncells;
- int nsurfs;
- int ngroups;
- int maxital;
-// Get the number of surfaces and cells for dynamic memory
- ncells = vol_handles().size()-1;
- nsurfs = surf_handles().size()-1;
- ngroups = group_handles().size()-1;
- MBRange vols_save, surfs_save;
- std::copy(vol_handles().begin()+1, vol_handles().end(),
- mb_range_inserter(vols_save));
- std::copy(surf_handles().begin()+1, surf_handles().end(),
- mb_range_inserter(surfs_save));
-
- std::cout << "ncells: " << ncells << std::endl;
- std::cout << "nsurfs: " << nsurfs << std::endl;
- std::cout << "ngroups: " << ngroups << std::endl;
- std::vector<int> mats(ncells);
- std::vector<double> dens(ncells);
- std::vector<int> imp(ncells);
- std::vector<int> surflist(nsurfs);
- std::vector<int> surfbcs(nsurfs);
- std::vector<std::string> gruplist;
-// Zero out material and density info
- for (int ii=1; ii <=ncells; ii++) {
- mats[ii-1] = 0;
- dens[ii-1] = 0.0;
- imp[ii-1] = 1;
+ std::cout << " # cells: " << vol_handles().size()-1 << std::endl;
+ std::cout << " # surfs: " << surf_handles().size()-1 << std::endl;
+ std::cout << " # groups: " << group_handles().size()-1 << std::endl;
+ std::cout << "# tallies: " << tallyList.size() << std::endl;
+
+ MBRange grp_sets, grp_ents;
+ MBErrorCode rval;
+
+ MBRange surfs, vols;
+
+
+ // write cell information (skip first entry)
+ for (iter = vol_handles().begin()+1; iter != vol_handles().end(); iter++) {
+
+ int cellid,matid;
+ double imp, density;
+
+ MBI->tag_get_data( idTag , &(*iter), 1, &cellid );
+ MBI->tag_get_data(impTag , &(*iter), 1, &imp );
+ MBI->tag_get_data(matTag , &(*iter), 1, &matid );
+ if (0 == matid) {
+ cgmfile << cellid << " " << matid << " imp:n=" << imp ;
+ } else {
+ MBI->tag_get_data(densTag, &(*iter), 1, &density );
+ cgmfile << cellid << " " << matid << " " << density << " imp:n=" << imp ;
+ }
+ if (*iter == impl_compl_handle)
+ cgmfile << " $ implicit complement ";
+ cgmfile << std::endl;
}
- maxital = -1;
- for (int ii=1; ii <=nsurfs; ii++)
- surfbcs[ii-1] = 0;
-// Read in each group name, check if a material specification
-// or surface boundary condition
- // get all the groups
+
+ // skip a line
+ cgmfile << std::endl;
+
+ // write surface info (skip first entry)
+ for (iter = surf_handles().begin()+1; iter != surf_handles().end(); iter++ ) {
+
+ int surfid, bc_id;
+
+ MBI->tag_get_data( idTag, &(*iter), 1, &surfid );
+ MBI->tag_get_data( bcTag, &(*iter), 1, &bc_id );
+
+ if (BC_SPEC == bc_id)
+ cgmfile << "*";
+ if (BC_WHITE == bc_id)
+ cgmfile << "+";
+ cgmfile << surfid << std::endl;
+
+ }
+
+ // add a final blank line
+ cgmfile << std::endl;
+
std::vector<std::string> grp_names;
- for (int ii=1; ii <=ngroups; ii++) {
- // get group & names
- MBEntityHandle group = group_handles()[ii];
+ std::vector<std::string> tokens;
+
+ // write tally info
+ for (std::vector<int>::iterator grp = tallyList.begin(); grp != tallyList.end(); grp++) {
+
+ std::stringstream tallyCard;
+
+ MBEntityHandle group = group_handles()[*grp];
grp_names.clear();
bool success = get_group_names(group, grp_names);
assert(success);
if (grp_names.empty()) continue;
+
+ // get sets associated with this group
+ rval = MBI->get_entities_by_type(group, MBENTITYSET, grp_sets);
+ if (MB_SUCCESS != rval) continue;
- // get sets in the group
- MBRange grp_sets;
- MBErrorCode result = moab_instance()->get_entities_by_type(group, MBENTITYSET, grp_sets);
- if (MB_SUCCESS != result) continue;
-
- if (grp_names[0].find("mat",0) == 0 && (grp_names[0].find("rho",0) > 0)) {
- // Extract material ID and density
- int cpos = grp_names[0].find("mat",0) + 4;
- std::string matid_str;
- while ((int(grp_names[0][cpos]) >= 48) && (int(grp_names[0][cpos]) <= 57)) {
- matid_str += grp_names[0][cpos];
- cpos++;
+ tokens.clear();
+ tokenize(grp_names[0],tokens,"_");
+
+ // Get user number for tally
+ int talNum = atoi(tokens[1].c_str());
+ char talMod = ' ';
+
+ // Get tally modifier
+ if ( 'e' == tokens[2][0] )
+ talMod = '*';
+ if ( 'q' == tokens[2][0] )
+ talMod = '+';
+ if (talMod != ' ')
+ tokens[2].erase(0,1); // strip special character for modified tallies
+
+ // Get tally type
+ if ( tallyKeywords.count(tokens[2]) == 0) {
+ std::cout << "Invalid tally defined in geometry: " << grp_names[0] << std::endl;
+ continue;
+ }
+
+ // Create MCNP number for tally by combining user number and MCNP type code
+ talNum = talNum*10 + tallyKeywords[tokens[2]];
+
+ // get tally particles
+ std::string partList = "n";
+ if (tokens.size()>3)
+ {
+ partList = tokens[3].substr(0,1);
+ for (std::string::iterator pos = tokens[3].begin()+1; pos != tokens[3].end(); pos++)
+ {
+ partList += ",";
+ partList += *pos;
+ }
}
- cpos = grp_names[0].find("rho",0) + 4;
- std::string density_str;
- while (((int(grp_names[0][cpos]) >= 48) && (int(grp_names[0][cpos]) <= 57)) ||
- (int(grp_names[0][cpos]) == 46) || (int(grp_names[0][cpos]) == 45) ||
- (int(grp_names[0][cpos]) == 69) || (int(grp_names[0][cpos]) == 101) ||
- (int(grp_names[0][cpos]) == 68) || (int(grp_names[0][cpos]) == 100)) {
- density_str += grp_names[0][cpos];
- cpos++;
- }
- const char *cmatid = matid_str.c_str();
- const char *cdensity = density_str.c_str();
- int matid = atoi(cmatid);
- double density = atof(cdensity);
+
+ // write tally type info
+ tallyCard << talMod << "f" << talNum << ":" << partList ;
+
+ int dim;
+ if (tallyKeywords[tokens[2]] < 3)
+ dim = 2;
+ else
+ dim = 3;
- // check for material definition in complement
- if (grp_names[0].find("comp",0) < 20) {
+ // get entities of correct dimension
+ const void* dim_val[] = {&dim};
+ MBRange ents;
+ rval = MBI->get_entities_by_type_and_tag( 0, MBENTITYSET, &geomTag, dim_val, 1, ents );
+ if (MB_SUCCESS != rval)
+ return rval;
+ MBRange grp_ents = intersect(grp_sets,ents);
- std::cout<<"write_log: complement material and density specified" << std::endl;
- mats[ncells-1]=matid;
- dens[ncells-1]=density;
+ for (MBRange::iterator ent = grp_ents.begin(); ent != grp_ents.end(); ent++)
+ tallyCard << " " << get_entity_id(*ent);
- }
- else {
+ tallyCard << " T";
- // Get the volumes in the current group
- MBRange grp_vols = intersect( grp_sets, vols_save);
- for (MBRange::iterator rit = grp_vols.begin();
- rit != grp_vols.end(); rit++) {
- MBEntityHandle vol = *rit;
- int cell_num = index_by_handle(vol);
- mats[cell_num-1]=matid;
- dens[cell_num-1]=density;
- }
+ std::string tallyCardStr = tallyCard.str();
+ tallyCard.str("");
+ // write tallyCardStr with 80 char limit per line
+
+ while ( tallyCardStr.length() > 72 )
+ {
+ std::string::size_type pos = tallyCardStr.rfind(' ',72);
+ cgmfile << tallyCardStr.substr(0,pos) << " &" << std::endl;
+ tallyCardStr.erase(0,pos);
+
+ cgmfile << " ";
}
+ cgmfile << tallyCardStr << std::endl;
+
+ }
+
+ cgmfile.close();
+ return MB_SUCCESS;
+}
+
+
+
+// desired functionality
+//
+// 1. for each group
+// a) check for material assignment -> get list of member vols and tag mat & dens
+// b) check for b.c. -> get list of surfaces and tag b.c.
+// c) check for graveyard -> tag imp
+
+MBErrorCode DagMC::parse_metadata()
+{
+ std::vector<std::string> grp_names;
+ std::vector<std::string> tokens;
+ std::map<std::string,int> keywords;
+
+ MBEntityHandle group;
+ MBRange grp_sets, grp_ents;
+ MBErrorCode rval;
+
+ MBRange vols, surfs;
+
+ std::copy(vol_handles().begin()+1, vol_handles().end(),
+ mb_range_inserter(vols));
+ std::copy(surf_handles().begin()+1, surf_handles().end(),
+ mb_range_inserter(surfs));
+
+ keywords["mat"] = MAT_GROUP;
+ keywords["spec.reflect"] = BC_SPEC;
+ keywords["white.reflect"] = BC_WHITE;
+ keywords["graveyard"] = IMP_ZERO;
+ keywords["outside.world"] = IMP_ZERO;
+ keywords["rest.of.world"] = IMP_ZERO;
+ keywords["tally"] = TALLY_GROUP;
+
+ for (unsigned int grp=0; grp < group_handles().size(); grp++) {
+
+ // get group handle & names
+ group = group_handles()[grp];
+ grp_names.clear();
+ bool success = get_group_names(group, grp_names);
+ assert(success);
+ if (grp_names.empty()) continue;
+
+ // get sets associated with this group
+ grp_sets.clear();
+ rval = MBI->get_entities_by_type(group, MBENTITYSET, grp_sets);
+ if (MB_SUCCESS != rval) continue;
+
+ tokens.clear();
+ tokenize(grp_names[0],tokens,"_");
+
+ // is this a keyword?
+ if (keywords.count(tokens[0]) == 0 ) {
+ std::cout << "Ignoring group with name: " << grp_names[0] << std:: endl;
+ continue;
}
- else if ((grp_names[0].find("spec_reflect",0)==0) ||
- (grp_names[0].find("white_reflect",0)==0)) {
- // Check whether group is reflecing or white
- int bc_id;
- if (grp_names[0].find("spec_reflect",0)==0)
- bc_id = 1;
- else
- bc_id = 2;
- // Get the surfaces in the current group
- MBRange grp_faces = intersect( grp_sets, surfs_save);
- for (MBRange::iterator rit = grp_faces.begin(); rit != grp_faces.end(); rit++) {
- MBEntityHandle surf = *rit;
- int surf_num = index_by_handle(surf);
- surfbcs[surf_num-1]=bc_id;
+
+ int matid, bc_id;
+ std::vector<int> matid_list, bc_id_list;
+ double density, imp_zero = 0;
+ std::vector<double> density_list, imp_zero_list;
+
+ // actions for each keyword
+ switch (keywords[tokens[0]]) {
+ case MAT_GROUP:
+ // missing density
+ if ( tokens[2].compare("rho") != 0 ) {
+ std::cout << "Specified material with no density " << grp_names[0] << std::endl;
+ continue;
}
- }
- else if (grp_names[0].find("periodic",0)==0) {
- // Get the periodic surface id
- int bc_id;
- std::string bc_str;
- int cpos = 9;
- bc_str.clear();
- while ((int(grp_names[0][cpos]) >= 48) && (int(grp_names[0][cpos]) <= 57)) {
- bc_str += grp_names[0][cpos];
- cpos++;
+
+ matid = atoi(tokens[1].c_str());
+ density = atof(tokens[3].c_str());
+
+ // check to see if this is the implicit complement
+ if ( tokens.size() > 4 && tokens[4].compare("comp") == 0 ) {
+ std::cout<<"parse_metadata: complement material and density specified" << std::endl;
+ // TAG IMPL COMPL WITH MAT & DENS
+ MBI->tag_set_data(matTag ,&impl_compl_handle,1,&matid);
+ MBI->tag_set_data(densTag,&impl_compl_handle,1,&density);
+ } else {
+ // apply tag to range of volumes
+ grp_ents.clear();
+ grp_ents = intersect(grp_sets.vols);
+ matid_list.assign(grp_ents.size(),matid);
+ density_list.assign(grp_ents.size(),density);
+ MBI->tag_set_data(matTag ,grp_ents,&(*matid_list.begin()));
+ MBI->tag_set_data(densTag,grp_ents,&(*density_list.begin()));
}
- const char *cbc_id = bc_str.c_str();
- bc_id = -1*atoi(cbc_id);
- // Get the surfaces in the current group
- MBRange grp_faces = intersect( grp_sets, surfs_save);
- for (MBRange::iterator rit = grp_faces.begin(); rit != grp_faces.end(); rit++) {
- MBEntityHandle surf = *rit;
- int surf_num = index_by_handle(surf);
- surfbcs[surf_num-1]=bc_id;
- }
+ break;
+ case BC_SPEC:
+ case BC_WHITE:
+ bc_id = keywords[tokens[0]];
+ grp_ents.clear();
+ grp_ents = intersect(grp_sets.surfs);
+ bc_id_list.assign(grp_ents.size(),bc_id);
+ MBI->tag_set_data(bcTag ,grp_ents, &(*bc_id_list.begin()));
+ break;
+ case IMP_ZERO:
+ grp_ents.clear();
+ grp_ents = intersect(grp_sets,vols);
+ imp_zero_list.assign(grp_ents.size(),imp_zero);
+ MBI->tag_set_data(impTag ,grp_ents, &(*imp_zero_list.begin()));
+ break;
+ case TALLY_GROUP:
+ // make list of groups that are tallies
+ grp_ents.clear();
+ grp_ents = intersect(grp_sets,vols);
+ tallyList.push_back(grp);
+ break;
}
-// Set graveyard to zero importance
- else if ((grp_names[0].find("graveyard",0)==0)||(grp_names[0].find("outside_world",0)==0)
- ||(grp_names[0].find("rest_of_world",0)==0)) {
- // Get the volumes in the current group
- MBRange grp_vols = intersect( grp_sets, vols_save);
- for (MBRange::iterator rit = grp_vols.begin(); rit != grp_vols.end(); rit++) {
- MBEntityHandle vol = *rit;
- int vol_num = index_by_handle(vol);
- imp[vol_num-1]=0;
- }
- }
-// Check for tally groups, assign maximum CUBIT index
- else if (grp_names[0].find("tally_",0)==0) {
- std::string stal;
- int cpos = 6;
- while ((int(grp_names[0][cpos]) >= 48) && (int(grp_names[0][cpos]) <= 57)) {
- stal += grp_names[0][cpos];
- cpos++;
- }
- const char *ctal = stal.c_str();
- if (atoi(ctal) > maxital)
- maxital = atoi(ctal);
- }
+
}
-// Generate list of surfaces
- for (int ii=1; ii <= nsurfs; ii++) {
- MBEntityHandle surf = surf_handles()[ii];
- surflist[ii-1] = get_entity_id(surf);
-// cout << "surf#: " << ii << " = " << surflist[ii-1] << endl;
- }
- // write in cell information
- for (int ii=1; ii <= ncells; ii++) {
- if (mats[ii-1] == 0) {
- cgmfile << id_by_index(3, ii) << " " << mats[ii-1] << " " <<
- "imp:n=" << imp[ii-1] << std::endl;
+ return MB_SUCCESS;
+}
+
+void DagMC::tokenize( const std::string& str,
+ std::vector<std::string>& tokens,
+ const char* delimiters )
+{
+ 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();
}
- else {
- cgmfile << id_by_index(3, ii) << " " << mats[ii-1] << " " <<
- dens[ii-1] << " " << "imp:n=" << imp[ii-1] << std::endl;
- }
- }
- // skip a line
- cgmfile << std::endl;
- // write out surface information (need boundary conditions)
- for (int ii=1; ii <= nsurfs; ii++) {
- if (surfbcs[ii-1] == 1)
- cgmfile << "*" << surflist[ii-1] << std::endl;
- else if (surfbcs[ii-1] == 2)
- cgmfile << "+" << surflist[ii-1] << std::endl;
- else if (surfbcs[ii-1] < 0)
- cgmfile << surflist[ii-1] << " " << surfbcs[ii-1] << std::endl;
- else
- cgmfile << surflist[ii-1] << std::endl;
- }
- // add a final blank line
- cgmfile << " ";
- // if tallies were specified, write those out at end of .cgm file
- if (maxital >= 0) {
- cgmfile << std::endl;
- // loop through each possible CUBIT tally index, write information
- for (int ii=0; ii <= maxital; ii++) {
- // loop through each group name, check to see if the desired tally type
- for (int jj=1; jj <= ngroups; jj++) {
- MBEntityHandle group = group_handles()[jj];
- grp_names.clear();
- bool success = get_group_names(group, grp_names);
- assert(success);
- if (!success || grp_names.empty()) continue;
- std::string talstr;
- talstr = "tally_" + itos(ii) + "_";
- const char *ctalstr = talstr.c_str();
- if (grp_names[0].find(ctalstr,0)==0) {
- unsigned int cpos = talstr.length();
- int etal = 0;
- // check if next character is "e" or "q" denoting energy or charge
- // need to also check if next character conforms to acceptable argument
- if ((int(grp_names[0][cpos]) == 101) &&
- ((int(grp_names[0][cpos+1]) == 115) ||
- (int(grp_names[0][cpos+1]) == 99) || (int(grp_names[0][cpos+1]) == 112))) {
- etal = 1;
- cpos++;
- }
- else if ((int(grp_names[0][cpos]) == 113) && (int(grp_names[0][cpos+1]) == 112)) {
- etal = -1;
- cpos++;
- }
- // check to obtain tally type
- int ital = 0;
- int tidx = 0;
- int tklen = 0;
- if (grp_names[0].find("surf_current",cpos)==cpos) {
- tidx = 1;
- ital = 10*ii + tidx;
- tklen = 13;
- }
- else if (grp_names[0].find("surf_flux",cpos)==cpos) {
- tidx = 2;
- ital = 10*ii + tidx;
- tklen = 10;
- }
- else if (grp_names[0].find("cell_flux",cpos)==cpos) {
- tidx = 4;
- ital = 10*ii + tidx;
- tklen = 10;
- }
- else if (grp_names[0].find("cell_heating",cpos)==cpos) {
- tidx = 6;
- ital = 10*ii + tidx;
- tklen = 13;
- }
- else if (grp_names[0].find("cell_fission",cpos)==cpos) {
- tidx = 7;
- ital = 10*ii + tidx;
- tklen = 13;
- }
- else if (grp_names[0].find("pulse_height",cpos)==cpos) {
- tidx = 8;
- ital = 10*ii + tidx;
- tklen = 13;
- }
- if (tidx > 0) {
- int ipos = 1;
- // place the tally header in the .cgm file
- if (etal == 1) {
- cgmfile << "*";
- ipos++;
- }
- else if ((etal == -1) && (tidx == 8)) {
- cgmfile << "+";
- ipos++;
- }
- cgmfile << "f" << ital << ":";
- std::stringstream spos;
- spos << ital;
- ipos += spos.str().length() + 2;
- // get the particles the tally covers (default is n)
- cpos += tklen;
- std::string talptls;
- while((cpos < grp_names[0].length()) && (int(grp_names[0][cpos]) >= 65) &&
- (int(grp_names[0][cpos]) <= 122) && (int(grp_names[0][cpos]) != 95)) {
- talptls += grp_names[0][cpos];
- cpos++;
- }
- // read the tally particle string and place in .cgm file
- if (talptls.empty()) {
- cgmfile << "n" << " ";
- ipos += 2;
- }
- else {
- for(unsigned int kk=1; kk <= talptls.length(); kk++) {
- cgmfile << talptls[kk-1];
- if (kk < talptls.length())
- cgmfile << ",";
- else
- cgmfile << " ";
- ipos += 2*talptls.length();
- }
- }
- // begin read of surfaces in current group
- MBRange grp_sets;
- MBErrorCode result = moab_instance()->get_entities_by_type(group, MBENTITYSET, grp_sets);
- if (MB_SUCCESS != result) continue;
- if ((tidx == 1) || (tidx == 2)) {
- MBRange grp_faces = intersect( grp_sets, surfs_save);
- for (MBRange::iterator rit = grp_faces.begin(); rit != grp_faces.end(); rit++) {
- int jtal = get_entity_id(*rit);
- std::stringstream spos1;
- spos1 << jtal;
- ipos += spos1.str().length() + 1;
- if (ipos > 78) {
- cgmfile << "&" << std::endl << " ";
- ipos = spos1.str().length() + 7;
- }
- cgmfile << jtal << " ";
- }
- }
- else if ((tidx == 4) || ((tidx >= 6) && (tidx <= 8))) {
- MBRange grp_vols = intersect( grp_sets, vols_save);
- for (MBRange::iterator rit = grp_vols.begin(); rit != grp_vols.end(); rit++) {
- int jtal = get_entity_id(*rit);
- std::stringstream spos1;
- spos1 << jtal;
- ipos += spos1.str().length() + 1;
- if (ipos > 78) {
- cgmfile << "&" << std::endl << " ";
- ipos = spos1.str().length() + 7;
- }
- cgmfile << jtal << " ";
- }
- }
- // set default behavior to provide total over all tally regions
- if (ipos+2 > 78) {
- cgmfile << "&" << std::endl << " ";
- }
- cgmfile << " T" << std::endl;
- }
- }
- }
- }
- }
- cgmfile.close();
}
std::string DagMC::itos(int ival) {
@@ -1364,7 +1355,7 @@
// get names
char name0[NAME_TAG_SIZE];
std::fill(name0, name0+NAME_TAG_SIZE, '\0');
- MBErrorCode result = moab_instance()->tag_get_data(name_tag(), &group_set, 1,
+ MBErrorCode result = MBI->tag_get_data(name_tag(), &group_set, 1,
&name0);
if (MB_SUCCESS != result && MB_TAG_NOT_FOUND != result) return false;
@@ -1378,7 +1369,7 @@
false);
if (0 == this_tag) break;
std::fill(name0, name0+NAME_TAG_SIZE, '\0');
- result = moab_instance()->tag_get_data(this_tag, &group_set, 1, &name0);
+ result = MBI->tag_get_data(this_tag, &group_set, 1, &name0);
if (MB_SUCCESS != result && MB_TAG_NOT_FOUND != result) return false;
if (MB_TAG_NOT_FOUND == result) break;
else grp_names.push_back(std::string(name0));
@@ -1387,12 +1378,13 @@
return true;
}
-MBTag DagMC::get_tag( const char* name, int size, MBTagType store, MBDataType type,
- bool create_if_missing)
+MBTag DagMC::get_tag( const char* name, int size, MBTagType store,
+ MBDataType type, const void* def_value,
+ bool create_if_missing)
{
MBTag retval = 0;
- MBErrorCode result = moab_instance()->tag_create(name, size, store, type,
- retval, NULL, create_if_missing);
+ MBErrorCode result = MBI->tag_create(name, size, store, type,
+ retval, def_value, create_if_missing);
if (create_if_missing && MB_SUCCESS != result)
std::cerr << "Couldn't find nor create tag named " << name << std::endl;
@@ -1402,9 +1394,9 @@
int DagMC::get_entity_id(MBEntityHandle this_ent)
{
int id = 0;
- MBErrorCode result = moab_instance()->tag_get_data(idTag, &this_ent, 1, &id);
+ MBErrorCode result = MBI->tag_get_data(idTag, &this_ent, 1, &id);
if (MB_TAG_NOT_FOUND == result)
- id = moab_instance()->id_from_handle(this_ent);
+ id = MBI->id_from_handle(this_ent);
return id;
}
@@ -1501,8 +1493,7 @@
#endif
}
-MBErrorCode DagMC::build_indices(MBRange &surfs, MBRange &vols,
- bool is_geom)
+MBErrorCode DagMC::build_indices(MBRange &surfs, MBRange &vols)
{
MBErrorCode rval = MB_SUCCESS;
@@ -1538,7 +1529,7 @@
for (MBRange::iterator rit = vols.begin(); rit != vols.end(); rit++) {
entIndices[*rit-setOffset] = idx++;
int result=0;
- moab_instance()->tag_get_data( idTag, &*rit, 1, &result );
+ MBI->tag_get_data( idTag, &*rit, 1, &result );
max_id = (max_id > result ? max_id : result);
}
// add implicit complement to entity index
@@ -1546,7 +1537,7 @@
// assign ID to implicit complement
max_id++;
- moab_instance()->tag_set_data(idTag, &impl_compl_handle, 1, &max_id);
+ MBI->tag_set_data(idTag, &impl_compl_handle, 1, &max_id);
@@ -1556,7 +1547,7 @@
// get geometry entities by id and cache in this vector
std::vector<int> ids;
ids.resize(surfs.size());
- rval = moab_instance()->tag_get_data(id_tag(), surfs, &ids[0]);
+ rval = MBI->tag_get_data(id_tag(), surfs, &ids[0]);
if (MB_SUCCESS != rval) return MB_FAILURE;
int i = 0;
MBRange::iterator rit = surfs.begin();
@@ -1567,7 +1558,7 @@
geomEntities[*rit - setOffset] = this_surf;
}
ids.resize(vols.size());
- rval = moab_instance()->tag_get_data(id_tag(), vols, &ids[0]);
+ rval = MBI->tag_get_data(id_tag(), vols, &ids[0]);
if (MB_SUCCESS != rval) return MB_FAILURE;
i = 0;
rit = vols.begin();
@@ -1588,8 +1579,8 @@
sprintf(group_category, "%s", "Group");
const void* const group_val[] = {&group_category};
MBRange groups;
- rval = moab_instance()->get_entities_by_type_and_tag(0, MBENTITYSET, &category_tag,
- group_val, 1, groups);
+ rval = MBI->get_entities_by_type_and_tag(0, MBENTITYSET, &category_tag,
+ group_val, 1, groups);
if (MB_SUCCESS != rval)
return rval;
group_handles().resize(groups.size()+1);
@@ -1599,7 +1590,7 @@
// populate root sets vector
std::vector<MBEntityHandle> rsets;
rsets.resize(surfs.size());
- rval = moab_instance()->tag_get_data(obb_tag(), surfs, &rsets[0]);
+ rval = MBI->tag_get_data(obb_tag(), surfs, &rsets[0]);
if (MB_SUCCESS != rval) return MB_FAILURE;
MBRange::iterator rit;
int i;
@@ -1607,14 +1598,14 @@
rootSets[*rit-setOffset] = rsets[i];
rsets.resize(vols.size());
- rval = moab_instance()->tag_get_data(obb_tag(), vols, &rsets[0]);
+ rval = MBI->tag_get_data(obb_tag(), vols, &rsets[0]);
if (MB_SUCCESS != rval) return MB_FAILURE;
for (i = 0, rit = vols.begin(); rit != vols.end(); rit++, i++)
rootSets[*rit-setOffset] = rsets[i];
// add implicit complement root set
rsets.resize(1);
- rval = moab_instance()->tag_get_data(obb_tag(), &impl_compl_handle, 1,
+ rval = MBI->tag_get_data(obb_tag(), &impl_compl_handle, 1,
&rsets[0]);
if (MB_SUCCESS != rval) return MB_FAILURE;
rootSets[impl_compl_handle-setOffset] = rsets[0];
@@ -1629,7 +1620,7 @@
for (MBRange::iterator i = surfs.begin(); i != surfs.end(); ++i) {
MBEntityHandle root;
MBRange tris;
- rval = moab_instance()->get_entities_by_dimension( *i, 2, tris );
+ rval = MBI->get_entities_by_dimension( *i, 2, tris );
if (MB_SUCCESS != rval)
return rval;
if (tris.empty())
@@ -1637,10 +1628,10 @@
rval = obbTree.build( tris, root );
if (MB_SUCCESS != rval)
return rval;
- rval = moab_instance()->add_entities( root, &*i, 1 );
+ rval = MBI->add_entities( root, &*i, 1 );
if (MB_SUCCESS != rval)
return rval;
- rval = moab_instance()->tag_set_data( obbTag, &*i, 1, &root );
+ rval = MBI->tag_set_data( obbTag, &*i, 1, &root );
if (MB_SUCCESS != rval)
return rval;
}
@@ -1648,7 +1639,7 @@
for (MBRange::iterator i = vols.begin(); i != vols.end(); ++i) {
// get all surfaces in volume
MBRange tmp_surfs;
- rval = moab_instance()->get_child_meshsets( *i, tmp_surfs );
+ rval = MBI->get_child_meshsets( *i, tmp_surfs );
if (MB_SUCCESS != rval)
return rval;
@@ -1667,7 +1658,7 @@
if (!sense)
continue;
- rval = moab_instance()->tag_get_data( obbTag, &*j, 1, &root );
+ rval = MBI->tag_get_data( obbTag, &*j, 1, &root );
if (MB_SUCCESS != rval || !root)
return MB_FAILURE;
trees.insert( root );
@@ -1678,7 +1669,7 @@
if (MB_SUCCESS != rval)
return rval;
- rval = moab_instance()->tag_set_data( obbTag, &*i, 1, &root );
+ rval = MBI->tag_set_data( obbTag, &*i, 1, &root );
if (MB_SUCCESS != rval)
return rval;
@@ -1705,13 +1696,13 @@
parent_vols.clear();
// get parents of each surface
- rval = moab_instance()->get_parent_meshsets( *surf_i, parent_vols );
+ rval = MBI->get_parent_meshsets( *surf_i, parent_vols );
if (MB_SUCCESS != rval)
return rval;
// if only one parent, get the OBB root for this surface
if (parent_vols.size() == 1 ) {
- rval = moab_instance()->tag_get_data( obbTag, &*surf_i, 1, &surf_obb_root );
+ rval = MBI->tag_get_data( obbTag, &*surf_i, 1, &surf_obb_root );
if (MB_SUCCESS != rval)
return rval;
if (!surf_obb_root)
@@ -1728,7 +1719,7 @@
return rval;
// tag the implicit complement handle with the handle for its own OBB tree
- rval = moab_instance()->tag_set_data( obbTag, &impl_compl_handle, 1, &comp_root );
+ rval = MBI->tag_set_data( obbTag, &impl_compl_handle, 1, &comp_root );
if (MB_SUCCESS != rval)
return rval;
@@ -1774,24 +1765,27 @@
MBErrorCode rval = mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET,
&nameTag, tagdata, 1,
entities );
+ // query error
if (MB_SUCCESS != rval) {
std::cerr << "Unable to query for implicit complement." << std::endl;
return rval;
}
+ // found too many
if (entities.size() > 1) {
std::cerr << "Too many implicit complement sets." << std::endl;
return MB_MULTIPLE_ENTITIES_FOUND;
}
+ // found none
if (entities.empty()) {
- rval= moab_instance()->create_meshset(MESHSET_SET,impl_compl_handle);
+ rval= MBI->create_meshset(MESHSET_SET,impl_compl_handle);
if (MB_SUCCESS != rval) {
std::cerr << "Failed to create mesh set for implicit complement." << std::endl;
return rval;
}
// tag this entity with name for implicit complement
- rval = moab_instance()->tag_set_data(nameTag,&impl_compl_handle,1,&implComplName);
+ rval = MBI->tag_set_data(nameTag,&impl_compl_handle,1,&implComplName);
if (MB_SUCCESS != rval) {
std::cerr << "Failed to tag new entity as implicit complement." << std::endl;
}
@@ -1799,9 +1793,11 @@
return rval;
} else {
+ // found a single implicit complement
impl_compl_handle = entities.front();
return MB_SUCCESS;
}
}
+
Modified: MOAB/trunk/tools/dagmc/DagMC.hpp
===================================================================
--- MOAB/trunk/tools/dagmc/DagMC.hpp 2009-07-14 14:57:15 UTC (rev 3005)
+++ MOAB/trunk/tools/dagmc/DagMC.hpp 2009-07-14 16:49:39 UTC (rev 3006)
@@ -74,13 +74,19 @@
// {-1 -> reversed, 0 -> both, 1 -> forward}
MBErrorCode surface_sense( MBEntityHandle volume, MBEntityHandle surface, int& sense_out );
- // load mesh and initialize
- MBErrorCode load_file_and_init(const char* cfile,
- const int clen,
- const char* ffile,
- const int flen,
- const double facet_tolerance = 0);
+ // load mesh
+ MBErrorCode load_file(const char* cfile,
+ const double facet_tolerance = 0);
+ MBErrorCode init_OBBtree();
+
+ MBErrorCode write_mesh(const char* ffile,
+ const int flen);
+
+
+ // initialize data structures and OBB tree
+ MBErrorCode init_OBBTree();
+
// map between MBEntityHandle, base-1 index, and GLOBAL_ID
MBEntityHandle entity_by_index( int dimension, int index );
MBEntityHandle entity_by_id( int dimension, int id );
@@ -92,10 +98,12 @@
// if dimension == 2).
int num_entities( int dimension );
- // read/write settings
+ // read/write settings from file
void read_settings( const char* filename );
void write_settings( FILE* filp, bool with_description = true );
void parse_settings();
+
+ // pass settings from calling application
void set_settings(int source_cell, int use_cad, int use_dist_limit,
double add_distance_tolerance,
double discard_distance_tolerance);
@@ -125,7 +133,11 @@
std::vector<MBEntityHandle>& group_handles() {return entHandles[4];}
- void write_log(std::string ifile, const bool overwrite = true);
+ MBErrorCode write_mcnp(std::string ifile, const bool overwrite = true);
+ MBErrorCode parse_metadata();
+ void tokenize( const std::string& str,
+ std::vector<std::string>& tokens,
+ const char* delimiters );
MBErrorCode poly_solid_angle( MBEntityHandle face, const MBCartVect& point, double& area );
@@ -154,9 +166,9 @@
double faceting_tolerance() {return facetingTolerance;}
- int source_cell() {return moabMCNPSourceCell;}
+ int source_cell() {return sourceCell;}
- bool use_dist_limit() {return moabMCNPUseDistLimit;}
+ bool use_dist_limit() {return useDistLimit;}
bool use_cad() {return useCAD;}
@@ -181,15 +193,14 @@
MBErrorCode get_impl_compl();
MBTag get_tag( const char* name, int size, MBTagType store, MBDataType type,
- bool create_if_missing = true);
+ const void* def_value = NULL, bool create_if_missing = true);
bool get_group_names(MBEntityHandle group_set, std::vector<std::string> &grp_names);
static void create_instance(MBInterface *mb_impl = NULL);
// build internal index vectors that speed up handle-by-id, etc.
- MBErrorCode build_indices(MBRange &surfs, MBRange &vols,
- bool is_geom);
+ MBErrorCode build_indices(MBRange &surfs, MBRange &vols);
// build obb structure
MBErrorCode build_obbs(MBRange &surfs, MBRange &vols);
@@ -204,12 +215,16 @@
bool user_set;
};
+ std::vector<int> tallyList;
+
std::string itos(int ival);
MBInterface *mbImpl;
MBOrientedBoxTreeTool obbTree;
MBTag obbTag, geomTag, idTag, nameTag, senseTag;
+ // metadata
+ MBTag matTag, densTag, bcTag, impTag, tallyTag;
Option options[6];
@@ -224,9 +239,10 @@
double discardDistTol;
double addDistTol;
double facetingTolerance;
- int moabMCNPSourceCell;
- bool moabMCNPUseDistLimit;
+ int sourceCell;
+ bool useDistLimit;
bool useCAD;
+ bool is_geom;
double distanceLimit;
Modified: MOAB/trunk/tools/dagmc/main.cc
===================================================================
--- MOAB/trunk/tools/dagmc/main.cc 2009-07-14 14:57:15 UTC (rev 3005)
+++ MOAB/trunk/tools/dagmc/main.cc 2009-07-14 16:49:39 UTC (rev 3006)
@@ -42,6 +42,8 @@
const char* input_name = 0;
const char* output_name = 0;
+ char options[256];
+
double dist_tol = 0.001, len_tol = 0.0;
int norm_tol = 5;
bool actuate_attribs = true;
@@ -112,21 +114,7 @@
if (!output_name)
usage(argv[0]);
-
-
- // Initialize CGM
- InitCGMA::initialize_cgma();
- if (actuate_attribs) {
- CGMApp::instance()->attrib_manager()->set_all_auto_read_flags( actuate_attribs );
- CGMApp::instance()->attrib_manager()->set_all_auto_actuate_flags( actuate_attribs );
- }
-
-
- // Intitalize MOAB
- MBCore moab;
- MBInterface* iface = &moab;
-
- // Get CGM file type
+
if (!file_type) {
file_type = get_geom_file_type( input_name );
if (!file_type) {
@@ -135,8 +123,17 @@
}
}
- // If CUB file, extract ACIS geometry
- CubitStatus s;
+
+ sprintf(options,"CGM_ATTRIBS=%s;FACET_DISTANCE_TOLERANCE=%g;FACET_NORMAL_TOLERANCE=%d;MAX_FACET_EDGE_LENGTH=%g;",
+ actuate_attribs?"yes":"no",dist_tol,norm_tol,len_tol);
+
+ // Intitalize MOAB
+ MBCore moab;
+ MBInterface* iface = &moab;
+ MBEntityHandle file_set;
+ MBErrorCode rval;
+
+// If CUB file, extract ACIS geometry
if (!strcmp( file_type, "CUBIT" )) {
char* temp_name = tempnam( 0, "cgm" );
if (!temp_name) {
@@ -158,36 +155,29 @@
exit(2);
}
- int rval = cub_file_type( cub_file, tmp_file, CUB_FILE_ACIS );
+ int crval = cub_file_type( cub_file, tmp_file, CUB_FILE_ACIS );
fclose( cub_file );
fclose( tmp_file );
- if (rval) {
+ if (crval) {
remove( temp_name );
free( temp_name );
exit(2);
}
-
- s = GeometryQueryTool::instance()->import_solid_model( temp_name, "ACIS_SAT" );
+
+ rval = iface->load_file(temp_name,file_set,options,NULL,0,0);
remove( temp_name );
free( temp_name );
}
else {
- s = GeometryQueryTool::instance()->import_solid_model( input_name, file_type );
+ rval = iface->load_file(input_name,file_set,options,NULL,0,0);
}
- if (CUBIT_SUCCESS != s) {
+ if (MB_SUCCESS != rval) {
std::cerr << "Failed to read '" << input_name << "' of type '" << file_type << "'" << std::endl;
exit(2);
}
-
- // copy geometry facets into mesh database
- if (!cgm2moab(iface, dist_tol, norm_tol,
- len_tol, actuate_attribs)) {
- std::cerr << "Internal error copying geometry" << std::endl;
- exit(5);
- }
-
+
// write mesh database
- MBErrorCode rval = iface->write_mesh( output_name );
+ rval = iface->write_mesh( output_name );
if (MB_SUCCESS != rval) {
std::cerr << "Failed to write '" << output_name << "'" << std::endl;
exit(2);
More information about the moab-dev
mailing list