[MOAB-dev] r1222 - in MOAB/trunk: . test/obb tools/dagmc
kraftche at mcs.anl.gov
kraftche at mcs.anl.gov
Mon Jul 23 14:20:35 CDT 2007
Author: kraftche
Date: 2007-07-23 14:20:35 -0500 (Mon, 23 Jul 2007)
New Revision: 1222
Modified:
MOAB/trunk/MBOrientedBoxTreeTool.cpp
MOAB/trunk/MBOrientedBoxTreeTool.hpp
MOAB/trunk/test/obb/obb_test.cpp
MOAB/trunk/tools/dagmc/DagMC.cpp
Log:
Change algorithm for closest-point and point-in-volume.
o Old closest point function did:
- find closest point on boundary
- calculate distance d between input point and closest point
- return all points within d+epsilon of input point
o Old closest point function very susceptable to rounding error
if input point is far away (d+epsilon == d for large values of d).
o New closest point function just returns the closest single point
o New addtional function to intersect sphere with triangles
o New point-in-volume test does:
- get closest point
- get any triangles within epsilon of closest point by doing
sphere-intersect-trinagles where center == closest point
and radius == epsilon.
This method is better because
o eliminates potential source of rounding error
o greatly simplifies point-in-volume because the result set
of triangles are all topologically adjacent - no need to
search larger set of tris from old closest-point method
for a subset that are adjacent
This method might be slower because of the second tree traversal.
However, it may not be because the closest-point query is a little
faster and the sphere-intersect-triangles query is much quicker
than either closest-point query.
Fixes bug from Paul (tube_1e-6.h5m).
Modified: MOAB/trunk/MBOrientedBoxTreeTool.cpp
===================================================================
--- MOAB/trunk/MBOrientedBoxTreeTool.cpp 2007-07-20 18:24:11 UTC (rev 1221)
+++ MOAB/trunk/MBOrientedBoxTreeTool.cpp 2007-07-23 19:20:35 UTC (rev 1222)
@@ -480,7 +480,107 @@
return MB_SUCCESS;
}
+/********************** General Sphere/Triangel Intersection ***************/
+MBErrorCode MBOrientedBoxTreeTool::sphere_intersect_triangles(
+ const double* center_v,
+ double radius,
+ MBEntityHandle tree_root,
+ std::vector<MBEntityHandle>& facets_out,
+ std::vector<MBEntityHandle>* sets_out )
+{
+ const double radsqr = radius * radius;
+ MBOrientedBox b;
+ MBErrorCode rval;
+ MBRange sets;
+ const MBCartVect center(center_v);
+ MBCartVect closest, coords[3];
+ const MBEntityHandle* conn;
+ int num_conn;
+#ifndef MB_OBB_USE_VECTOR_QUERIES
+ MBRange tris;
+ MBRange::const_iterator t;
+#else
+ std::vector<MBEntityHandle> tris;
+ std::vector<MBEntityHandle>::const_iterator t;
+#endif
+
+ std::vector<MBEntityHandle> stack, children;
+ stack.reserve(60);
+ stack.push_back(tree_root);
+ stack.push_back(0);
+ while (!stack.empty()) {
+ MBEntityHandle surf = stack.back(); stack.pop_back();
+ MBEntityHandle node = stack.back(); stack.pop_back();
+ if (!surf && sets_out) {
+ rval = get_moab_instance()->get_entities_by_type( node, MBENTITYSET, sets );
+ if (!sets.empty())
+ surf = sets.front();
+ }
+
+ // check if sphere intersects box
+ rval = box( node, b );
+ if (MB_SUCCESS != rval)
+ return rval;
+ b.closest_location_in_box( center, closest );
+ closest -= center;
+ if ((closest % closest) > radsqr)
+ continue;
+
+ // push child boxes on stack
+ children.clear();
+ rval = instance->get_child_meshsets( node, children );
+ if (MB_SUCCESS != rval)
+ return rval;
+ if (!children.empty()) {
+ assert(children.size() == 2);
+ stack.push_back( children[0] );
+ stack.push_back( surf );
+ stack.push_back( children[1] );
+ stack.push_back( surf );
+ continue;
+ }
+
+ // if leaf, intersect sphere with triangles
+#ifndef MB_OBB_USE_VECTOR_QUERIES
+# ifdef MB_OBB_USE_TYPE_QUERIES
+ rval = get_moab_instance()->get_entities_by_type( node, MBTRI, tris );
+# else
+ rval = get_moab_instance()->get_entities_by_handle( node, tris );
+# endif
+#else
+ rval = get_moab_instance()->get_entities_by_handle( node, tris );
+#endif
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ for (t = tris.begin(); t != tris.end(); ++t) {
+ rval = get_moab_instance()->get_connectivity( *t, conn, num_conn, true );
+ if (MB_SUCCESS != rval)
+ return rval;
+ if (num_conn != 3)
+ continue;
+
+ rval = get_moab_instance()->get_coords( conn, num_conn, coords[0].array() );
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ MBGeomUtil::closest_location_on_tri( center, coords, closest );
+ closest -= center;
+ if ((closest % closest) <= radsqr &&
+ std::find(facets_out.begin(),facets_out.end(),*t) == facets_out.end()) {
+ facets_out.push_back( *t );
+ if (sets_out)
+ sets_out->push_back( surf );
+ }
+ }
+ }
+
+ return MB_SUCCESS;
+}
+
+
+
/********************** General Ray/Tree and Ray/Triangel Intersection ***************/
@@ -843,53 +943,134 @@
/********************** Closest Point code ***************/
+struct MBOBBTreeCPFrame {
+ MBOBBTreeCPFrame( double d, MBEntityHandle n, MBEntityHandle s )
+ : dist_sqr(d), node(n), mset(s) {}
+ double dist_sqr;
+ MBEntityHandle node;
+ MBEntityHandle mset;
+};
+
MBErrorCode MBOrientedBoxTreeTool::closest_to_location(
const double* point,
MBEntityHandle root,
- double tolerance,
double* point_out,
MBEntityHandle& facet_out,
MBEntityHandle* set_out )
{
- std::vector<MBEntityHandle> facets(1), sets;
- MBErrorCode rval = closest_to_location( point, root, tolerance, facets, set_out ? &sets : 0 );
- if (MB_SUCCESS != rval)
- return rval;
- if (facets.empty())
- return MB_ENTITY_NOT_FOUND;
- facet_out = facets.front();
- if (set_out)
- *set_out = sets.front();
+ MBErrorCode rval;
+ const MBCartVect loc( point );
+ double smallest_dist_sqr = std::numeric_limits<double>::max();
- const MBEntityHandle* conn;
- int len;
- rval = instance->get_connectivity( facet_out, conn, len, true );
- if (MB_SUCCESS != rval)
- return rval;
+ MBEntityHandle current_set = 0;
+ MBRange sets;
+ std::vector<MBEntityHandle> children(2);
+ std::vector<double> coords;
+ std::vector<MBOBBTreeCPFrame> stack;
+
+ stack.push_back( MBOBBTreeCPFrame(0.0, root, current_set) );
- std::vector<MBCartVect> coords( len );
- rval = instance->get_coords( conn, len, coords[0].array() );
- if (MB_SUCCESS != rval)
- return rval;
+ while( !stack.empty() ) {
+
+ // pop from top of stack
+ MBEntityHandle node = stack.back().node;
+ double dist_sqr = stack.back().dist_sqr;
+ current_set = stack.back().mset;
+ stack.pop_back();
+
+ // If current best result is closer than the box, skip this tree node.
+ if (dist_sqr > smallest_dist_sqr)
+ continue;
+
+ // Check if this node has a set associated with it
+ if (set_out && !current_set) {
+ sets.clear();
+ rval = instance->get_entities_by_type( node, MBENTITYSET, sets );
+ if (MB_SUCCESS != rval)
+ return rval;
+ if (!sets.empty()) {
+ if (sets.size() != 1)
+ return MB_MULTIPLE_ENTITIES_FOUND;
+ current_set = sets.front();
+ }
+ }
+
+ // Get child boxes
+ children.clear();
+ rval = instance->get_child_meshsets( node, children );
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ // if not a leaf node
+ if (!children.empty()) {
+ if (children.size() != 2)
+ return MB_MULTIPLE_ENTITIES_FOUND;
+
+ // get boxes
+ MBOrientedBox box1, box2;
+ rval = box( children[0], box1 );
+ if (MB_SUCCESS != rval) return rval;
+ rval = box( children[1], box2 );
+ if (MB_SUCCESS != rval) return rval;
+
+ // get distance from each box
+ MBCartVect pt1, pt2;
+ box1.closest_location_in_box( loc, pt1 );
+ box2.closest_location_in_box( loc, pt2 );
+ pt1 -= loc;
+ pt2 -= loc;
+ const double dsqr1 = pt1 % pt1;
+ const double dsqr2 = pt2 % pt2;
+
+ // push children on tree such that closer one is on top
+ if (dsqr1 < dsqr2) {
+ stack.push_back( MBOBBTreeCPFrame(dsqr2, children[1], current_set ) );
+ stack.push_back( MBOBBTreeCPFrame(dsqr1, children[0], current_set ) );
+ }
+ else {
+ stack.push_back( MBOBBTreeCPFrame(dsqr1, children[0], current_set ) );
+ stack.push_back( MBOBBTreeCPFrame(dsqr2, children[1], current_set ) );
+ }
+ }
+ else { // LEAF NODE
+ MBRange facets;
+ rval = instance->get_entities_by_dimension( node, 2, facets );
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ const MBEntityHandle* conn;
+ int len;
+ MBCartVect tmp, diff;
+ for (MBRange::iterator i = facets.begin(); i != facets.end(); ++i) {
+ rval = instance->get_connectivity( *i, conn, len, true );
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ coords.resize(3*len);
+ rval = instance->get_coords( conn, len, &coords[0] );
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ if (len == 3)
+ MBGeomUtil::closest_location_on_tri( loc, (MBCartVect*)(&coords[0]), tmp );
+ else
+ MBGeomUtil::closest_location_on_polygon( loc, (MBCartVect*)(&coords[0]), len, tmp );
+
+ diff = tmp - loc;
+ dist_sqr = diff % diff;
+ if (dist_sqr < smallest_dist_sqr) {
+ smallest_dist_sqr = dist_sqr;
+ facet_out = *i;
+ tmp.get( point_out );
+ if (set_out)
+ *set_out = current_set;
+ }
+ }
+ } // LEAF NODE
+ }
- if (len == 3)
- MBGeomUtil::closest_location_on_tri( *reinterpret_cast<const MBCartVect*>(point),
- &coords[0],
- *reinterpret_cast<MBCartVect*>(point_out) );
- else
- MBGeomUtil::closest_location_on_polygon( *reinterpret_cast<const MBCartVect*>(point),
- &coords[0], len,
- *reinterpret_cast<MBCartVect*>(point_out) );
return MB_SUCCESS;
}
-
-struct MBOBBTreeCPFrame {
- MBOBBTreeCPFrame( double d, MBEntityHandle n, MBEntityHandle s )
- : dist_sqr(d), node(n), mset(s) {}
- double dist_sqr;
- MBEntityHandle node;
- MBEntityHandle mset;
-};
MBErrorCode MBOrientedBoxTreeTool::closest_to_location( const double* point,
MBEntityHandle root,
Modified: MOAB/trunk/MBOrientedBoxTreeTool.hpp
===================================================================
--- MOAB/trunk/MBOrientedBoxTreeTool.hpp 2007-07-20 18:24:11 UTC (rev 1221)
+++ MOAB/trunk/MBOrientedBoxTreeTool.hpp 2007-07-23 19:20:35 UTC (rev 1222)
@@ -203,7 +203,6 @@
*/
MBErrorCode closest_to_location( const double* point,
MBEntityHandle tree_root,
- double tolerance,
double* point_out,
MBEntityHandle& facet_out,
MBEntityHandle* set_out = 0);
@@ -221,6 +220,23 @@
std::vector<MBEntityHandle>& facets_out,
std::vector<MBEntityHandle>* sets_out = 0 );
+ /**\brief Find facets intersected by a sphere
+ *
+ * Find facets intersected by a spherical volume.
+ *\param center Sphere center
+ *\param radius Sphere radius
+ *\param tree_root Root of OBB tree
+ *\param facets_out List of triangles intersecting sphere
+ *\param sets_out If not null, sets owning facets are appended to this
+ * list in an order corresponding to the entries in
+ * facets_out.
+ */
+ MBErrorCode sphere_intersect_triangles( const double* center,
+ double radius,
+ MBEntityHandle tree_root,
+ std::vector<MBEntityHandle>& facets_out,
+ std::vector<MBEntityHandle>* sets_out = 0 );
+
/**\brief Get oriented box at node in tree
*
* Get the oriented box for a node in an oriented bounding box tree.
Modified: MOAB/trunk/test/obb/obb_test.cpp
===================================================================
--- MOAB/trunk/test/obb/obb_test.cpp 2007-07-20 18:24:11 UTC (rev 1221)
+++ MOAB/trunk/test/obb/obb_test.cpp 2007-07-23 19:20:35 UTC (rev 1222)
@@ -1295,7 +1295,6 @@
// find closest point usnig tree
rval = tool.closest_to_location( points[i].array(),
root_set,
- tolerance,
t_result.array(),
t_tri,
set_ptr );
Modified: MOAB/trunk/tools/dagmc/DagMC.cpp
===================================================================
--- MOAB/trunk/tools/dagmc/DagMC.cpp 2007-07-20 18:24:11 UTC (rev 1221)
+++ MOAB/trunk/tools/dagmc/DagMC.cpp 2007-07-23 19:20:35 UTC (rev 1222)
@@ -371,164 +371,60 @@
double x, double y, double z,
int& result )
{
+ MBErrorCode rval;
const double epsilon = moabMCNPTolerance;
// Get OBB Tree for volume
assert(volume - setOffset < rootSets.size());
MBEntityHandle root = rootSets[volume - setOffset];
- // Get closest triangles in volume
+ // Get closest point in triangulation
+ const MBCartVect point(x,y,z);
+ MBEntityHandle facet;
+ MBCartVect closest, diff;
+ rval = obbTree.closest_to_location( point.array(), root, closest.array(), facet );
+ if (MB_SUCCESS != rval) return rval;
+
+ // Check for on-boundary case
+ diff = closest - point;
+ if (diff%diff <= epsilon*epsilon) {
+ result = -1;
+ return MB_SUCCESS;
+ }
+
+ // Get triangles at closest point
static std::vector<MBEntityHandle> tris, surfs;
tris.clear();
surfs.clear();
- const MBCartVect point(x,y,z);
- MBErrorCode rval = obbTree.closest_to_location( point.array(), root, epsilon, tris, &surfs );
+ rval = obbTree.sphere_intersect_triangles( closest.array(), epsilon, root, tris, &surfs );
if (MB_SUCCESS != rval) return rval;
- // Get sense of each surface wrt volume
- std::vector<int> senses( surfs.size() );
- rval = surface_sense( volume, surfs.size(), &surfs[0], &senses[0] );
- if (MB_SUCCESS != rval) return rval;
-
- // Get closest point on each triangle, and group triangles
- // by intersection type
- std::vector<int> edge_tris, vertex_tris;
-
- // Get corrected normal, closest point and topolgocal closest for each triangle
- std::vector<MBCartVect> normals( tris.size() ), closest( tris.size() );
- std::vector<int> topo( tris.size() );
+ // One-triangle case : determine if point is above or below triangle
const MBEntityHandle* conn;
- int len;
- MBCartVect coords[3];
- for (unsigned i = 0; i < tris.size(); ++i) {
- // volume tree shouldn't contain non-manifold surfaces
- assert(senses[i]);
-
- rval = moab_instance()->get_connectivity( tris[i], conn, len );
+ int len, sense;
+ MBCartVect coords[3], normal;
+ if (tris.size() == 1) {
+ rval = moab_instance()->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() );
- if (MB_SUCCESS != rval)
- return rval;
+ if (MB_SUCCESS != rval) return rval;
- MBGeomUtil::closest_location_on_tri( point, coords, epsilon, closest[i], topo[i] );
+ rval = surface_sense( volume, 1, &surfs.front(), &sense );
+ if (MB_SUCCESS != rval) return rval;
+
+ // get triangle normal
coords[1] -= coords[0];
coords[2] -= coords[0];
- normals[i] = senses[i] * coords[1] * coords[2];
- }
-
- // If within tolerance of triangle, then return that point is on boundary
- if ((point - closest[0]).length() < epsilon) {
- result = -1;
+ normal = sense * coords[1] * coords[2];
+ // compare relative sense
+ result = (normal % (point - closest) < 0.0);
return MB_SUCCESS;
}
- // if only one, then return
- if (tris.size() == 1) {
- // Ignore the possibility of being in the plane of the
- // triangle. If parallel to one triangle, should have
- // at least one other triangle.
- result = (normals[0] % (closest[0] - point) > 0.0);
- return MB_SUCCESS;
- }
+ // many triangle case : determine if closest point is convexity or concavity
- // If we found any triangle for which the closest point was in
- // the interior, just use that one
- std::vector<int>::iterator ii = std::find( topo.begin(), topo.end(), 6 );
- if (ii != topo.end()) {
- int idx = ii - topo.begin();
- result = (normals[idx] % (closest[idx] - point) > 0.0);
- return MB_SUCCESS;
- }
- // Otherwise if we found any triangle for which the
- // closest point was at an edge, find it's neighbor and use the pair
- for (ii = topo.begin(); ii != topo.end() && *ii < 3; ++ii);
- if (ii != topo.end()) // found one
- {
- unsigned idx1 = ii - topo.end();
- int idx2 = -1;
- MBEntityHandle edge[2];
- moab_instance()->get_connectivity( tris[idx1], conn, len, true );
- edge[0] = conn[*ii - 3];
- edge[1] = conn[(*ii - 2) % 3];
- for (unsigned i = 0; i < tris.size(); ++i) {
- if (i == idx1)
- continue;
- moab_instance()->get_connectivity( tris[i], conn, len, true );
- int j;
- for (j = 0; j < 3; ++j)
- if (conn[j] == edge[0])
- if (conn[(j+1)%3] == edge[1] || conn[(j+2)%3] == edge[1])
- break;
- if (j < 3) {
- idx2 = i;
- break;
- }
- }
-
- if (idx2 < 0) { // if no adjacent triangle
- assert(false);
- result = (normals[idx1] % (closest[idx1] - point) > 0.0);
- return MB_SUCCESS;
- }
-
- // Check if inside or outside of both facets first.
- // This avoids potential rounding/accuracy issues
- // when testing for convexity/concavity if the facets
- // are near co-planer, and handles the special case
- // where both facets are exactly coplanar.
- const double dot1 = normals[idx1] % (closest[idx1] - point);
- const double dot2 = normals[idx2] % (closest[idx2] - point);
- const bool inside1 = dot1 > -std::numeric_limits<double>::epsilon();
- const bool inside2 = dot2 > -std::numeric_limits<double>::epsilon();
- const bool outside1 = dot1 < std::numeric_limits<double>::epsilon();
- const bool outside2 = dot2 < std::numeric_limits<double>::epsilon();
- if (inside1 && inside2)
- result = 1;
- else if (outside1 && outside2)
- result = 0;
- else {
- // Edge is at a concavity if the projection of the normal of
- // triangle A into the plane of triangle B points into triangle
- // B. It is a convexity if it points out of triangle B (from the
- // shared edge). This test can be reduced to checking the sign of
- // dot product of the normal of triangle A with any vector pointing
- // into triangle B from the shared edge.
-
- // Get triangle coordinates
- moab_instance()->get_connectivity( tris[idx2], conn, len, true );
- moab_instance()->get_coords( conn, len, coords[0].array() );
- // Get a vector from edge end point to opposite vertex
- const MBCartVect vect2 = coords[(topo[idx2] - 1) % 3] - coords[topo[idx2]];
- result = (vect2 % normals[idx1] >= 0.0); // true if a concavity
- result = inside1 || inside2;
- }
- return MB_SUCCESS;
- }
-
- // Otherwise pick a vertex that the point was closest to remove
- // any triangles from lists that are not adjacent to that vertex.
- assert( topo[0] < 3 ); // if here, must be closest to vertex
- rval = moab_instance()->get_connectivity( tris[0], conn, len );
- if (MB_SUCCESS != rval || len != 3)
- return MB_FAILURE;
- const MBEntityHandle closest_vert = conn[topo[0]];
- unsigned w = 1;
- for (unsigned r = 1; r < tris.size(); ++r) {
- assert( topo[r] < 3 ); // if here, must be closest to vertex
- moab_instance()->get_connectivity( tris[r], conn, len );
- if (conn[topo[r]] == closest_vert) {
- tris[w] = tris[r];
- surfs[w] = surfs[r];
- normals[w] = normals[r];
- ++w;
- }
- }
- tris.resize(w);
- surfs.resize(w);
- normals.resize(w);
-
// use algorithm from:
// (referance)
// to determine inside vs. outside.
@@ -536,20 +432,40 @@
// first find single triangle for which all other triangles are to
// one side.
bool some_above, some_below;
- MBCartVect vertcoords, closestcoords;
- moab_instance()->get_coords( &closest_vert, 1, closestcoords.array() );
for (unsigned i = 0; i < tris.size(); ++i) {
some_above = some_below = false;
+
+ rval = moab_instance()->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() );
+ if (MB_SUCCESS != rval) return rval;
+
+ rval = surface_sense( volume, 1, &surfs.front(), &sense );
+ if (MB_SUCCESS != rval) return rval;
+
+ // get triangle normal
+ coords[1] -= coords[0];
+ coords[2] -= coords[0];
+ normal = sense * coords[1] * coords[2];
+ closest = coords[0];
+ const double norm_len_sqr = normal % normal;
+
for (unsigned j = 0; j < tris.size(); ++j) {
if (j == i)
continue;
- moab_instance()->get_connectivity( tris[j], conn, len );
+ rval = moab_instance()->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() );
+ if (MB_SUCCESS != rval)
+ return rval;
+
for (int k = 0; k < len; ++k) {
- moab_instance()->get_coords( conn+k, 1, vertcoords.array() );
- const double dot = normals[i] % (closestcoords - vertcoords);
- if (dot * dot / (normals[i] % normals[i]) > epsilon * epsilon) {
- if (dot < 0.0)
+ const double dot = normal % (coords[k] - closest);
+ if (dot * dot / norm_len_sqr > epsilon * epsilon) {
+ if (dot > 0.0)
some_above = true;
else
some_below = true;
@@ -560,7 +476,7 @@
// All triangles are roughly co-planar: if we're inside of
// one then we're inside of all.
if (!some_above && !some_below) {
- result = (normals[i] % (closestcoords - point)) > 0.0;
+ result = (normal % (closest - point)) > 0.0;
return MB_SUCCESS;
}
// All other triangles to one side of this triangle:
@@ -586,9 +502,6 @@
// detemine distance to nearest surface
MBErrorCode DagMC::closest_to_location( MBEntityHandle volume, double* coords, double& result)
{
-
- const double epsilon = moabMCNPTolerance;
-
// Get OBB Tree for volume
assert(volume - setOffset < rootSets.size());
MBEntityHandle root = rootSets[volume - setOffset];
@@ -597,7 +510,7 @@
const MBCartVect point(coords);
MBCartVect nearest;
MBEntityHandle facet_out;
- MBErrorCode rval = obbTree.closest_to_location( point.array(), root, epsilon, nearest.array(), facet_out );
+ MBErrorCode rval = obbTree.closest_to_location( point.array(), root, nearest.array(), facet_out );
if (MB_SUCCESS != rval) return rval;
// calculate distance between point and nearest facet
More information about the moab-dev
mailing list