[MOAB-dev] r3418 - in MOAB/trunk: . tools/dagmc
bmsmith6 at wisc.edu
bmsmith6 at wisc.edu
Fri Dec 18 14:20:07 CST 2009
Author: bmsmith
Date: 2009-12-18 14:20:06 -0600 (Fri, 18 Dec 2009)
New Revision: 3418
Modified:
MOAB/trunk/MBOrientedBoxTreeTool.cpp
MOAB/trunk/MBOrientedBoxTreeTool.hpp
MOAB/trunk/tools/dagmc/DagMC.cpp
MOAB/trunk/tools/dagmc/DagMC.hpp
MOAB/trunk/tools/dagmc/ray_fire_test.cc
Log:
Pass facets back from ray casting code.
-This is a short term fix in support of mesh-based particle transport
-This might be part of the long term robust transport algorithm
Add a simple test to the ray_fire testing tool.
-Use this to fire one ray at a single volume
Modified: MOAB/trunk/MBOrientedBoxTreeTool.cpp
===================================================================
--- MOAB/trunk/MBOrientedBoxTreeTool.cpp 2009-12-16 19:50:46 UTC (rev 3417)
+++ MOAB/trunk/MBOrientedBoxTreeTool.cpp 2009-12-18 20:20:06 UTC (rev 3418)
@@ -795,11 +795,12 @@
const double tol;
std::vector<double>& intersections;
std::vector<MBEntityHandle>& sets;
-
- MBEntityHandle lastSet;
+ std::vector<MBEntityHandle>& facets;
+
+ MBEntityHandle lastSet;
int lastSetDepth;
- void add_intersection( double t );
+ void add_intersection( double t, MBEntityHandle facet );
public:
RayIntersectSets( MBOrientedBoxTreeTool* tool_ptr,
@@ -809,12 +810,13 @@
double tolerance,
unsigned min_tol_intersections,
std::vector<double>& intersections,
- std::vector<MBEntityHandle>& surfaces )
+ std::vector<MBEntityHandle>& surfaces,
+ std::vector<MBEntityHandle>& facets )
: tool(tool_ptr), minTolInt(min_tol_intersections),
b(ray_point), m(unit_ray_dir),
len(ray_length), tol(tolerance),
intersections(intersections),
- sets(surfaces), lastSet(0)
+ sets(surfaces), facets(facets), lastSet(0)
{ }
virtual MBErrorCode visit( MBEntityHandle node,
@@ -900,12 +902,12 @@
// NOTE: add_intersection may modify the 'len' member, which
// will affect subsequent calls to ray_tri_intersect in
// this loop.
- add_intersection( td );
+ add_intersection( td, *t );
}
return MB_SUCCESS;
}
-void RayIntersectSets::add_intersection( double t )
+void RayIntersectSets::add_intersection( double t, MBEntityHandle facet )
{
// Check if the 'len' pointer is pointing into the intersection
// list. If this is the case, then the list contains, at that
@@ -925,6 +927,7 @@
if (intersections.size() >= minTolInt) {
intersections[len_idx] = t;
sets[len_idx] = lastSet;
+ facets[len_idx] = facet;
// From now on, we want only intersections within the tolerance,
// so update length accordingly
len = &tol;
@@ -933,6 +936,7 @@
else {
intersections.push_back(t);
sets.push_back(lastSet);
+ facets.push_back(facet);
len = &intersections[len_idx];
}
}
@@ -940,6 +944,7 @@
else {
intersections.push_back(t);
sets.push_back(lastSet);
+ facets.push_back(facet);
// If we have all the intersections we want, set
// length such that we will only find further intersections
// within the tolerance
@@ -954,6 +959,7 @@
if (t <= *len) {
intersections[len_idx] = t;
sets[len_idx] = lastSet;
+ facets[len_idx] = facet;
}
}
// Otherwise if we want an intersection outside the tolerance
@@ -961,6 +967,7 @@
else if (intersections.size() < minTolInt) {
intersections.push_back( t );
sets.push_back( lastSet );
+ facets.push_back(facet);
// udpate length. this is currently the closest intersection, so
// only want further intersections that are closer than this one.
len = &intersections.back();
@@ -971,6 +978,7 @@
MBErrorCode MBOrientedBoxTreeTool::ray_intersect_sets(
std::vector<double>& distances_out,
std::vector<MBEntityHandle>& sets_out,
+ std::vector<MBEntityHandle>& facets_out,
MBEntityHandle root_set,
double tolerance,
unsigned min_tolerace_intersections,
@@ -979,7 +987,7 @@
const double* ray_length )
{
RayIntersectSets op( this, ray_point, unit_ray_dir, ray_length, tolerance,
- min_tolerace_intersections, distances_out, sets_out );
+ min_tolerace_intersections, distances_out, sets_out, facets_out );
return preorder_traverse( root_set, op );
}
Modified: MOAB/trunk/MBOrientedBoxTreeTool.hpp
===================================================================
--- MOAB/trunk/MBOrientedBoxTreeTool.hpp 2009-12-16 19:50:46 UTC (rev 3417)
+++ MOAB/trunk/MBOrientedBoxTreeTool.hpp 2009-12-18 20:20:06 UTC (rev 3418)
@@ -188,6 +188,7 @@
*/
MBErrorCode ray_intersect_sets( std::vector<double>& distances_out,
std::vector<MBEntityHandle>& sets_out,
+ std::vector<MBEntityHandle>& facets_out,
MBEntityHandle root_set,
double tolerance,
unsigned min_tolerace_intersections,
Modified: MOAB/trunk/tools/dagmc/DagMC.cpp
===================================================================
--- MOAB/trunk/tools/dagmc/DagMC.cpp 2009-12-16 19:50:46 UTC (rev 3417)
+++ MOAB/trunk/tools/dagmc/DagMC.cpp 2009-12-18 20:20:06 UTC (rev 3418)
@@ -157,6 +157,7 @@
// it for every call
std::vector<double> &distances = distList;
std::vector<MBEntityHandle> &surfaces = surfList;
+ std::vector<MBEntityHandle> &facets = facetList;
// do ray fire
@@ -164,11 +165,13 @@
const double dir[] = {uuu, vvv, www};
distances.clear();
surfaces.clear();
+ facets.clear();
double len = use_dist_limit() ? distance_limit() : huge_val;
unsigned min_tolerance_intersections = 1000;
rval = obbTree.ray_intersect_sets( distances,
- surfaces,
+ surfaces,
+ facets,
root,
add_dist_tol(),
min_tolerance_intersections,
Modified: MOAB/trunk/tools/dagmc/DagMC.hpp
===================================================================
--- MOAB/trunk/tools/dagmc/DagMC.hpp 2009-12-16 19:50:46 UTC (rev 3417)
+++ MOAB/trunk/tools/dagmc/DagMC.hpp 2009-12-18 20:20:06 UTC (rev 3418)
@@ -272,7 +272,7 @@
static DagMC *instance_;
// temporary storage so functions don't have to reallocate vectors
- std::vector<MBEntityHandle> triList, surfList;
+ std::vector<MBEntityHandle> triList, surfList, facetList;
std::vector<double> distList;
};
Modified: MOAB/trunk/tools/dagmc/ray_fire_test.cc
===================================================================
--- MOAB/trunk/tools/dagmc/ray_fire_test.cc 2009-12-16 19:50:46 UTC (rev 3417)
+++ MOAB/trunk/tools/dagmc/ray_fire_test.cc 2009-12-18 20:20:06 UTC (rev 3418)
@@ -52,6 +52,7 @@
std::cerr << "Usage: " << argv[0] << " [-f] <filename> "
<< " <facet_tol> <source_rad> <vol_index> <#calls> [<location_az>] [<direction_az>]" << std::endl;
std::cerr << "-f: report full tree statistics" << std::endl;
+ std::cerr << "-s: perform a simple ray_fire with <filename> <facet_tol> <vol_id> <x y z> <u v w>" << std::endl;
std::cerr << "filename: mesh or geometry filename" << std::endl;
std::cerr << "facet_tol: facet tolerance" << std::endl;
std::cerr << "source_rad: if < 0, ray at source_rad and random angle pointed at origin;" << std::endl;
@@ -63,10 +64,69 @@
return 1;
}
+ // declare variables used for all ray_fires
const double PI = acos(-1.0);
-
+ char *filename;
double facet_tol;
+ double x, y, z, u, v, w, dist;
int ncalls;
+ MBEntityHandle surf = 0, vol = 0;
+
+ // Simply test 1 ray against 1 volume with no timing
+ if(strcmp(argv[2], "-s")) {
+ if(11 > argc) {
+ std::cerr << "Need more input arguments for simple test" << std::endl;
+ return 1;
+ }
+
+ int i=2;
+ filename = argv[i++];
+ facet_tol = atof(argv[i++]);
+ int vol_id = atoi(argv[i++]);
+ x = atof(argv[i++]);
+ y = atof(argv[i++]);
+ z = atof(argv[i++]);
+ u = atof(argv[i++]);
+ v = atof(argv[i++]);
+ w = atof(argv[i++]);
+
+ DagMC& dagmc = *DagMC::instance();
+ rval = dagmc.load_file( filename, facet_tol );
+ if(MB_SUCCESS != rval) {
+ std::cerr << "Failed to load file." << std::endl;
+ return 2;
+ }
+
+ rval = dagmc.init_OBBTree( );
+ if(MB_SUCCESS != rval) {
+ std::cerr << "Failed to initialize DagMC." << std::endl;
+ return 2;
+ }
+
+ vol = dagmc.entity_by_id(3, vol_id);
+ if(0 == vol) {
+ std::cerr << "Problem getting volume." << std::endl;
+ return 2;
+ }
+
+ rval = dagmc.ray_fire( vol, 0, 1, u, v, w, x, y, z, DBL_MAX, dist, surf );
+ if(MB_SUCCESS != rval) {
+ std::cerr << "Problem with ray_fire" << std::endl;
+ return 2;
+ }
+
+ if(0 == surf) {
+ std::cerr << "Particle is lost" << std::endl;
+ return 2;
+ }
+
+ int surf_id = dagmc.id_by_index( 2, dagmc.index_by_handle( surf ) );
+ std::cout << "dist=" << dist << " surf_id=" << surf_id << " new_xyz="
+ << x+dist*u << " " << y+dist*v << " " << z+dist*w << std::endl;
+ return 0;
+
+ } else {
+
bool full = false;
int i = 1;
if (!strcmp(argv[i], "-f")) {
@@ -74,7 +134,7 @@
i++;
}
- char* filename = argv[i++];
+ filename = argv[i++];
facet_tol = atof(argv[i++]);
double rad = atof(argv[i++]);
int vol_idx = atoi(argv[i++]);
@@ -98,7 +158,7 @@
}
- MBEntityHandle vol = dagmc.entity_by_index(3, vol_idx);
+ vol = dagmc.entity_by_index(3, vol_idx);
if (0 == vol) {
std::cerr << "Problem getting first volume." << std::endl;
return 2;
@@ -109,8 +169,6 @@
// initialize random number generator using ttime1
srand((unsigned int) ttime1);
- double x, y, z, u, v, w, dist;
- MBEntityHandle nsurf;
#ifdef DEBUG
double uavg = 0.0, vavg = 0.0, wavg = 0.0;
@@ -134,7 +192,7 @@
#endif
dagmc.ray_fire(vol, 0, 1, u, v, w, x, y, z, DBL_MAX,
- dist, nsurf);
+ dist, surf);
}
get_time_mem(ttime2, utime2, stime2, tmem1);
@@ -184,6 +242,7 @@
#endif
return 0;
+ }
}
void get_time_mem(double &tot_time, double &user_time,
More information about the moab-dev
mailing list