[MOAB-dev] r1304 - MOAB/trunk/tools/dagmc
tautges at mcs.anl.gov
tautges at mcs.anl.gov
Wed Oct 10 09:34:46 CDT 2007
Author: tautges
Date: 2007-10-10 09:34:46 -0500 (Wed, 10 Oct 2007)
New Revision: 1304
Added:
MOAB/trunk/tools/dagmc/pt_vol_test.cc
Modified:
MOAB/trunk/tools/dagmc/DagMC.cpp
MOAB/trunk/tools/dagmc/DagMC.hpp
MOAB/trunk/tools/dagmc/Makefile.am
MOAB/trunk/tools/dagmc/cgm2moab.cc
MOAB/trunk/tools/dagmc/cgm2moab.hpp
MOAB/trunk/tools/dagmc/test_geom.cc
Log:
1. Adding parsing of useCAD option
2. Adding a facet tolerance option, and passing that to function which
gets faceting from geometry
3. Adding simple test program for pt_in_volume from Paul. To make,
either do a 'make check' in the root directory of MOAB, or in
tools/dagmc, do a 'make pt_vol_test'. Run the code w/o arguments to
find usage.
4. Adding Brandon's implicit complement code.
This has not been tested in MCNPX, so that should be done quickly.
Modified: MOAB/trunk/tools/dagmc/DagMC.cpp
===================================================================
--- MOAB/trunk/tools/dagmc/DagMC.cpp 2007-10-10 14:30:46 UTC (rev 1303)
+++ MOAB/trunk/tools/dagmc/DagMC.cpp 2007-10-10 14:34:46 UTC (rev 1304)
@@ -34,17 +34,21 @@
DagMC::DagMC(MBInterface *mb_impl)
: mbImpl(mb_impl), obbTree(mb_impl),
- moabMCNPTolerance(1e-6), moabMCNPSourceCell(0), moabMCNPUseDistLimit(false)
+ distanceTolerance(1e-6), facetingTolerance(0.001),
+ moabMCNPSourceCell(0), moabMCNPUseDistLimit(false)
{
options[0] = Option( "source_cell", "source cell ID, or zero if unknown", "0" );
options[1] = Option( "distance_tolerance", "positive real value", "0.001" );
options[2] = Option( "use_distance_limit", "one to enable distance limit optimization, zero otherwise", "0" );
options[3] = Option( "use_cad", "one to ray-trace to cad, zero for just facets", "0" );
+ options[4] = Option( "faceting_tolerance", "graphics faceting tolerance", "0.001" );
memset( specReflectName, 0, NAME_TAG_SIZE );
strcpy( specReflectName, "spec_reflect" );
memset( whiteReflectName, 0, NAME_TAG_SIZE );
strcpy( whiteReflectName, "white_reflect" );
+ memset( implComplName, 0, NAME_TAG_SIZE );
+ strcpy( implComplName , "impl_complement" );
distanceLimit = std::numeric_limits<double>::max();
@@ -95,7 +99,7 @@
rval = obbTree.ray_intersect_sets( distances,
surfaces,
root,
- tolerance(),
+ distance_tolerance(),
2,
point, dir,
&len);
@@ -104,7 +108,7 @@
#ifdef CGM
if (useCAD) {
#ifdef CUBIT_CGM
- std::cout << "USE_CAD = 1 not supported with this build of CGM/DagMC;"
+ std::cout << "use_cad = 1 not supported with this build of CGM/DagMC;"
<< std:: endl
<< "need an ACIS-based install of CGM." << std::endl;
return MB_NOT_IMPLEMENTED;
@@ -128,7 +132,8 @@
int smallest = std::min_element( distances.begin(), distances.end() ) - distances.begin();
// If intersected previous surface near start of ray, reject it
- if (surfaces[smallest] == last_surf_hit && distances[smallest] < tolerance()) {
+ if (surfaces[smallest] == last_surf_hit &&
+ distances[smallest] < distance_tolerance()) {
distances.erase( distances.begin() + smallest );
surfaces.erase( surfaces.begin() + smallest );
@@ -176,12 +181,23 @@
std::cerr << "Invalid source_cell = " << moabMCNPSourceCell << std::endl;
exit(2);
}
- moabMCNPTolerance = strtod( options[1].value.c_str(), 0 );
- if (moabMCNPTolerance <= 0 || moabMCNPTolerance > 1) {
- std::cerr << "Invalid distance_tolerance = " << moabMCNPTolerance << std::endl;
+ distanceTolerance = strtod( options[1].value.c_str(), 0 );
+ if (distanceTolerance <= 0 || distanceTolerance > 1) {
+ std::cerr << "Invalid distance_tolerance = " << distanceTolerance << std::endl;
exit(2);
}
+
moabMCNPUseDistLimit = !!atoi( options[2].value.c_str() );
+
+ useCAD = !!atoi( options[3].value.c_str() );
+#ifdef CUBIT_CGM
+ if (useCAD) {
+ std::cout << "Warning: use_cad = 1 not supported with this build of "
+ << "CGM/DagMC;" << std:: endl
+ << "need an ACIS-based install of CGM." << std::endl;
+ }
+#endif
+
}
void DagMC::read_settings( const char* filename )
@@ -376,7 +392,7 @@
int& result )
{
MBErrorCode rval;
- const double epsilon = moabMCNPTolerance;
+ const double epsilon = distanceTolerance;
// Get OBB Tree for volume
assert(volume - setOffset < rootSets.size());
@@ -534,6 +550,12 @@
std::vector<MBEntityHandle> surfaces, surf_volumes;
result = 0.0;
+ // don't try to calculate volume of implicit complement
+ if (volume == impl_compl_handle) {
+ result = 1.0;
+ return MB_SUCCESS;
+ }
+
// get surfaces from volume
rval = moab_instance()->get_child_meshsets( volume, surfaces );
if (MB_SUCCESS != rval) return rval;
@@ -704,7 +726,7 @@
<< std::endl;
exit(2);
}
- if (!cgm2moab(moab_instance())) {
+ if (!cgm2moab(moab_instance(), facetingTolerance)) {
std::cerr << "Internal error copying geometry" << std::endl;
exit(5);
}
@@ -751,6 +773,14 @@
if (MB_SUCCESS != rval)
return rval;
+ // If it doesn't already exist, create implicit complement
+ // Create data structures for implicit complement
+ rval = get_impl_compl();
+ if (MB_SUCCESS != rval) {
+ std::cerr << "Failed to find or create implicit complement handle." << std::endl;
+ return rval;
+ }
+
// Build OBB trees for everything, but only if we only read geometry
// Changed to build obb tree if tree does not already exist. -- JK
if (!have_obb_tree()) {
@@ -901,14 +931,28 @@
const char *cdensity = density_str.c_str();
int matid = atoi(cmatid);
double density = atof(cdensity);
- // Get the volumes in the current group
- MBRange grp_vols = grp_sets.intersect(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;
+
+ // check for material definition in complement
+ if (grp_names[0].find("comp",0) < 20) {
+
+ std::cout<<"write_log: complement material and density specified" << std::endl;
+ mats[ncells-1]=matid;
+ dens[ncells-1]=density;
+
}
+ else {
+
+ // Get the volumes in the current group
+ MBRange grp_vols = grp_sets.intersect(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;
+ }
+
+ }
}
else if ((grp_names[0].find("spec_reflect",0)==0) ||
(grp_names[0].find("white_reflect",0)==0)) {
@@ -1215,7 +1259,7 @@
const double in_pt[] = { xxx, yyy, zzz };
std::vector<MBEntityHandle> &facets = triList;
facets.clear();
- MBErrorCode rval = obbTree.closest_to_location( in_pt, root, tolerance(), facets );
+ MBErrorCode rval = obbTree.closest_to_location( in_pt, root, distance_tolerance(), facets );
assert(MB_SUCCESS == rval);
MBCartVect coords[3], normal(0.0);
@@ -1306,9 +1350,14 @@
// surf/vol offsets are just first handles
setOffset = (*surfs.begin() < *vols.begin() ? *surfs.begin() : *vols.begin());
- unsigned int tmp_offset = (surfs.back() > vols.back() ? surfs.back() : vols.back())
- - setOffset + 1;
- rootSets.resize(tmp_offset);
+ setOffset = (impl_compl_handle < setOffset ? impl_compl_handle : setOffset);
+ // max
+ unsigned int tmp_offset = (surfs.back() > vols.back() ?
+ surfs.back() : vols.back());
+ tmp_offset = (impl_compl_handle > tmp_offset ?
+ impl_compl_handle : tmp_offset);
+ // set size
+ rootSets.resize(tmp_offset - setOffset + 1);
entIndices.resize(rootSets.size());
// store surf/vol handles lists (surf/vol by index) and
@@ -1317,6 +1366,7 @@
std::vector<MBEntityHandle>::iterator iter = entHandles[2].begin();
*(iter++) = 0;
std::copy( surfs.begin(), surfs.end(), iter );
+ entHandles[3].push_back(impl_compl_handle);
int idx = 1;
for (MBRange::iterator rit = surfs.begin(); rit != surfs.end(); rit++)
entIndices[*rit-setOffset] = idx++;
@@ -1326,9 +1376,22 @@
*(iter++) = 0;
std::copy( vols.begin(), vols.end(), iter );
idx = 1;
- for (MBRange::iterator rit = vols.begin(); rit != vols.end(); rit++)
+ int max_id = -1;
+ 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 );
+ max_id = (max_id > result ? max_id : result);
+ }
+ // add implicit complement to entity index
+ entIndices[impl_compl_handle-setOffset] = idx++ ;
+ // assign ID to implicit complement
+ max_id++;
+ moab_instance()->tag_set_data(idTag, &impl_compl_handle, 1, &max_id);
+
+
+
#ifdef CGM
if (is_geom) {
geomEntities.resize(rootSets.size());
@@ -1391,6 +1454,13 @@
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,
+ &rsets[0]);
+ if (MB_SUCCESS != rval) return MB_FAILURE;
+ rootSets[impl_compl_handle-setOffset] = rsets[0];
+
return MB_SUCCESS;
}
@@ -1456,5 +1526,126 @@
}
+ rval = build_obb_impl_compl(surfs);
+ if (MB_SUCCESS != rval) {
+ std::cerr << "Unable to build OBB tree for implicit complement." << std::endl;
+ return rval;
+ }
+
return MB_SUCCESS;
}
+
+MBErrorCode DagMC::build_obb_impl_compl(MBRange &surfs)
+{
+ MBEntityHandle comp_root, surf_obb_root;
+ MBRange comp_tree;
+ MBErrorCode rval;
+ std::vector<MBEntityHandle> parent_vols;
+
+ // search through all surfaces
+ for (MBRange::iterator surf_i = surfs.begin(); surf_i != surfs.end(); ++surf_i) {
+
+ parent_vols.clear();
+ // get parents of each surface
+ rval = moab_instance()->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 );
+ if (MB_SUCCESS != rval)
+ return rval;
+ if (!surf_obb_root)
+ return MB_FAILURE;
+
+ // add obb root to list of obb roots
+ comp_tree.insert( surf_obb_root );
+ }
+ }
+
+ // join surface trees to make OBB tree for implicit complement
+ rval = obbTree.join_trees( comp_tree, comp_root );
+ if (MB_SUCCESS != rval)
+ 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 );
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ return MB_SUCCESS;
+
+}
+
+MBErrorCode DagMC::getobb(MBEntityHandle volume, double minPt[3],
+ double maxPt[3])
+{
+ double center[3], axis1[3], axis2[3], axis3[3];
+
+ // get center point and vectors to OBB faces
+ MBErrorCode rval = getobb(volume, center, axis1, axis2, axis3);
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ // compute min and max verticies
+ for (int i=0; i<3; i++)
+ {
+ minPt[i] = center[i] - (axis1[i] + axis2[i] + axis3[i]);
+ maxPt[i] = center[i] + (axis1[i] + axis2[i] + axis3[i]);
+ }
+ return MB_SUCCESS;
+}
+
+MBErrorCode DagMC::getobb(MBEntityHandle volume, double center[3], double axis1[3],
+ double axis2[3], double axis3[3])
+{
+ //find MBEntityHandle node_set for use in box
+ MBEntityHandle root = rootSets[volume - setOffset];
+
+ // call box to get center and vectors to faces
+ return obbTree.box(root, center, axis1, axis2, axis3);
+
+}
+
+MBErrorCode DagMC::get_impl_compl()
+{
+ MBRange entities;
+ const void* const tagdata[] = {implComplName};
+ MBErrorCode rval = mbImpl->get_entities_by_type_and_tag( 0, MBENTITYSET,
+ &nameTag, tagdata, 1,
+ entities );
+ if (MB_SUCCESS != rval) {
+ std::cerr << "Unable to query for implicit complement." << std::endl;
+ return rval;
+ }
+
+ if (entities.size() > 1) {
+ std::cerr << "Too many implicit complement sets." << std::endl;
+ return MB_MULTIPLE_ENTITIES_FOUND;
+ }
+
+ if (entities.empty()) {
+ rval= moab_instance()->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);
+ if (MB_SUCCESS != rval) {
+ std::cerr << "Failed to tag new entity as implicit complement." << std::endl;
+ }
+
+ return rval;
+
+ } else {
+ impl_compl_handle = entities.front();
+ return MB_SUCCESS;
+ }
+
+
+}
+
+
+
Modified: MOAB/trunk/tools/dagmc/DagMC.hpp
===================================================================
--- MOAB/trunk/tools/dagmc/DagMC.hpp 2007-10-10 14:30:46 UTC (rev 1303)
+++ MOAB/trunk/tools/dagmc/DagMC.hpp 2007-10-10 14:34:46 UTC (rev 1304)
@@ -120,13 +120,23 @@
MBErrorCode get_angle(MBEntityHandle surf,
double xxx, double yyy, double zzz, double *ang);
+
+ // get the corners of the OBB for a given volume
+ MBErrorCode getobb(MBEntityHandle volume, double minPt[3], double maxPt[3]);
+
+ // get the center point and three vectors for the OBB of a given volume
+ MBErrorCode getobb(MBEntityHandle volume, double center[3],
+ double axis1[3], double axis2[3], double axis3[3]);
+
int get_entity_id(MBEntityHandle this_ent);
// Get the instance of MOAB used by functions in this file.
MBInterface* moab_instance() {return mbImpl;}
- double tolerance() {return moabMCNPTolerance;}
+ double distance_tolerance() {return distanceTolerance;}
+ double faceting_tolerance() {return facetingTolerance;}
+
int source_cell() {return moabMCNPSourceCell;}
bool use_dist_limit() {return moabMCNPUseDistLimit;}
@@ -143,6 +153,8 @@
private:
bool have_obb_tree();
+ MBErrorCode get_impl_compl();
+
MBTag get_tag( const char* name, int size, MBTagType store, MBDataType type,
bool create_if_missing = true);
@@ -156,6 +168,7 @@
// build obb structure
MBErrorCode build_obbs(MBRange &surfs, MBRange &vols);
+ MBErrorCode build_obb_impl_compl(MBRange &surfs);
class Option {
public:
@@ -176,10 +189,14 @@
char specReflectName[NAME_TAG_SIZE];
char whiteReflectName[NAME_TAG_SIZE];
+ char implComplName[NAME_TAG_SIZE];
std::vector<MBEntityHandle> entHandles[5];
- double moabMCNPTolerance;
+ MBEntityHandle impl_compl_handle;
+
+ double distanceTolerance;
+ double facetingTolerance;
int moabMCNPSourceCell;
bool moabMCNPUseDistLimit;
bool useCAD;
Modified: MOAB/trunk/tools/dagmc/Makefile.am
===================================================================
--- MOAB/trunk/tools/dagmc/Makefile.am 2007-10-10 14:30:46 UTC (rev 1303)
+++ MOAB/trunk/tools/dagmc/Makefile.am 2007-10-10 14:34:46 UTC (rev 1304)
@@ -34,8 +34,10 @@
cfg_DATA = dagmc.make
TESTS = test_geom
-check_PROGRAMS = $(TESTS)
+check_PROGRAMS = $(TESTS) pt_vol_test
test_geom_SOURCES = test_geom.cc
test_geom_LDADD = libdagmc.la $(top_builddir)/libMOAB.la $(CGM_LIBS_LTFLAGS) $(CGM_LIBS_LINK)
+pt_vol_test_SOURCES = pt_vol_test.cc
+pt_vol_test_LDADD = libdagmc.la $(top_builddir)/libMOAB.la $(CGM_LIBS_LTFLAGS) $(CGM_LIBS_LINK)
Modified: MOAB/trunk/tools/dagmc/cgm2moab.cc
===================================================================
--- MOAB/trunk/tools/dagmc/cgm2moab.cc 2007-10-10 14:30:46 UTC (rev 1303)
+++ MOAB/trunk/tools/dagmc/cgm2moab.cc 2007-10-10 14:34:46 UTC (rev 1304)
@@ -41,7 +41,7 @@
// copy geometry into mesh database
bool cgm2moab(MBInterface* iface,
- double dist_tol,
+ double faceting_tol,
int norm_tol,
double len_tol,
int actuate_attribs)
@@ -285,7 +285,7 @@
GeometryQueryEngine* gqe = curve->get_geometry_query_engine();
data.clean_out();
int count;
- CubitStatus s = gqe->get_graphics( curve, count, &data, dist_tol);
+ CubitStatus s = gqe->get_graphics( curve, count, &data, faceting_tol);
if (CUBIT_SUCCESS != s)
return false;
@@ -371,7 +371,7 @@
data.clean_out();
int nt, np, nf;
CubitStatus s = gqe->get_graphics( surf, nt, np, nf, &data,
- norm_tol, dist_tol, len_tol);
+ norm_tol, faceting_tol, len_tol);
if (CUBIT_SUCCESS != s)
return false;
Modified: MOAB/trunk/tools/dagmc/cgm2moab.hpp
===================================================================
--- MOAB/trunk/tools/dagmc/cgm2moab.hpp 2007-10-10 14:30:46 UTC (rev 1303)
+++ MOAB/trunk/tools/dagmc/cgm2moab.hpp 2007-10-10 14:34:46 UTC (rev 1304)
@@ -14,7 +14,7 @@
#define CGM2MOAB_HPP
bool cgm2moab(MBInterface* iface,
- double dist_tol = 0.001,
+ double faceting_tol = 0.001,
int norm_tol = 5,
double len_tol = 0.0,
int actuate_attribs = true);
Added: MOAB/trunk/tools/dagmc/pt_vol_test.cc
===================================================================
--- MOAB/trunk/tools/dagmc/pt_vol_test.cc (rev 0)
+++ MOAB/trunk/tools/dagmc/pt_vol_test.cc 2007-10-10 14:34:46 UTC (rev 1304)
@@ -0,0 +1,92 @@
+#include "MBInterface.hpp"
+#include "MBCore.hpp"
+#include "DagMC.hpp"
+#include "MBTagConventions.hpp"
+
+#include <vector>
+#include <iostream>
+#include <math.h>
+#include <limits>
+
+#define CHKERR if (MB_SUCCESS != rval) return rval
+
+
+
+MBErrorCode test_pt_volume(DagMC &dagmc, int volID, double xxx, double yyy, double zzz, int &inside)
+{
+ MBErrorCode rval;
+
+ MBEntityHandle vol = dagmc.entity_by_id(3,volID);
+
+ rval = dagmc.point_in_volume( vol, xxx, yyy, zzz, inside);
+ CHKERR;
+
+ return MB_SUCCESS;
+
+}
+
+MBErrorCode test_pt_volume_slow(DagMC &dagmc, int volID, double xxx, double yyy, double zzz, int &inside)
+{
+ MBErrorCode rval;
+
+ MBEntityHandle vol = dagmc.entity_by_id(3,volID);
+
+ rval = dagmc.point_in_volume_slow( vol, xxx, yyy, zzz, inside);
+ CHKERR;
+
+ return MB_SUCCESS;
+
+}
+
+int main( int argc, char* argv[] )
+{
+ MBErrorCode rval;
+
+ if (argc < 6) {
+ std::cerr << "Usage: " << argv[0] << " <mesh_filename> "
+ << " <vol_id> <xxx> <yyy> <zzz> " << std::endl;
+ return 1;
+ }
+
+ char* filename = argv[1];
+ int volID = atoi(argv[2]);
+ double xxx = atof(argv[3]);
+ double yyy = atof(argv[4]);
+ double zzz = atof(argv[5]);
+ int inside;
+
+ std::cout << "Checking pt_in_volume for:" << std::endl
+ << "(x,y,z) = (" << xxx << "," << yyy << "," << zzz << ")"
+ << std::endl << "in volume " << volID
+ << " of geometry " << filename << std::endl;
+
+
+ DagMC& dagmc = *DagMC::instance();
+ rval = dagmc.load_file_and_init( filename, strlen(filename), 0, 0 );
+ if (MB_SUCCESS != rval) {
+ std::cerr << "Failed to initialize DagMC." << std::endl;
+ return 2;
+ }
+
+ int errors = 0;
+
+ rval = test_pt_volume(dagmc,volID,xxx,yyy,zzz,inside);
+ if (MB_SUCCESS != rval) {
+ std::cerr << "Failed to test point in volume [fast]." << std::endl;
+ return 3;
+ }
+
+ std::cout << "[fast] Point is " << (inside?"inside":"outside") << " volume "
+ << volID << std::endl;
+
+ rval = test_pt_volume_slow(dagmc,volID,xxx,yyy,zzz,inside);
+ if (MB_SUCCESS != rval) {
+ std::cerr << "Failed to test point in volume [slow]." << std::endl;
+ return 3;
+ }
+
+ std::cout << "[slow] Point is " << (inside?"inside":"outside") << " volume "
+ << volID << std::endl;
+
+ return errors;
+}
Modified: MOAB/trunk/tools/dagmc/test_geom.cc
===================================================================
--- MOAB/trunk/tools/dagmc/test_geom.cc 2007-10-10 14:30:46 UTC (rev 1303)
+++ MOAB/trunk/tools/dagmc/test_geom.cc 2007-10-10 14:34:46 UTC (rev 1304)
@@ -154,7 +154,6 @@
std::cout << #A << "... " << std::endl; \
if (MB_SUCCESS != A ( dagmc ) ) { \
++errors; \
- std::cout << "<<FAILED>>" << std::endl; \
} \
} \
} while(false)
@@ -180,7 +179,7 @@
}
int errors = 0;
- //RUN_TEST( test_ray_fire );
+ RUN_TEST( test_ray_fire );
RUN_TEST( test_point_in_volume );
RUN_TEST( test_measure_volume );
RUN_TEST( test_measure_area );
More information about the moab-dev
mailing list