[MOAB-dev] r1840 - MOAB/trunk
kraftche at mcs.anl.gov
kraftche at mcs.anl.gov
Fri May 23 11:38:10 CDT 2008
Author: kraftche
Date: 2008-05-23 11:38:10 -0500 (Fri, 23 May 2008)
New Revision: 1840
Added:
MOAB/trunk/MBBSPTree.cpp
MOAB/trunk/MBBSPTree.hpp
MOAB/trunk/bsp_tree_test.cpp
Modified:
MOAB/trunk/Makefile.am
MOAB/trunk/configure.in
Log:
BSP Tree
Added: MOAB/trunk/MBBSPTree.cpp
===================================================================
--- MOAB/trunk/MBBSPTree.cpp (rev 0)
+++ MOAB/trunk/MBBSPTree.cpp 2008-05-23 16:38:10 UTC (rev 1840)
@@ -0,0 +1,823 @@
+/*
+ * MOAB, a Mesh-Oriented datABase, is a software component for creating,
+ * storing and accessing finite element mesh data.
+ *
+ * Copyright 2008 Sandia Corporation. Under the terms of Contract
+ * DE-AC04-94AL85000 with Sandia Coroporation, the U.S. Government
+ * retains certain rights in this software.
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ */
+
+/**\file MBBSPTree.cpp
+ *\author Jason Kraftcheck (kraftche at cae.wisc.edu)
+ *\date 2008-05-13
+ */
+
+#include "MBBSPTree.hpp"
+#include "MBGeomUtil.hpp"
+#include "MBRange.hpp"
+#include "MBInternals.hpp"
+
+#include <assert.h>
+#include <algorithm>
+#include <limits>
+
+#ifdef _MSC_VER
+# include <float.h>
+# define finite(A) _finite(A)
+#elif defined(HAVE_IEEEFP_H)
+# include <ieeefp.h>
+#endif
+
+#define MB_BSP_TREE_DEFAULT_TAG_NAME "BSPTree"
+
+MBBSPTree::MBBSPTree( MBInterface* mb,
+ const char* tagname,
+ unsigned set_flags )
+ : mbInstance(mb), meshSetFlags(set_flags)
+{
+ MBErrorCode rval;
+
+ if (!tagname)
+ tagname = MB_BSP_TREE_DEFAULT_TAG_NAME;
+
+ std::string rootname(tagname);
+ rootname += "_box";
+
+ rval = mb->tag_create( tagname, 4*sizeof(double), MB_TAG_DENSE, MB_TYPE_DOUBLE, planeTag, 0, true );
+ if (MB_SUCCESS != rval)
+ planeTag = 0;
+
+ rval = mb->tag_create( rootname.c_str(), 6*sizeof(double), MB_TAG_SPARSE, MB_TYPE_DOUBLE, rootTag, 0, true );
+ if (MB_SUCCESS != rval)
+ rootTag = 0;
+}
+
+MBErrorCode MBBSPTree::set_split_plane( MBEntityHandle node, const Plane& p )
+{
+ // check for unit-length normal
+ const double lensqr = p.norm[0]*p.norm[0]
+ + p.norm[1]*p.norm[1]
+ + p.norm[2]*p.norm[2];
+ if (fabs(lensqr - 1.0) < std::numeric_limits<double>::epsilon())
+ return moab()->tag_set_data( planeTag, &node, 1, &p );
+
+ const double inv_len = 1.0/sqrt(lensqr);
+ Plane p2(p);
+ p2.norm[0] *= inv_len;
+ p2.norm[1] *= inv_len;
+ p2.norm[2] *= inv_len;
+ p2.coeff *= inv_len;
+
+ // check for zero-length normal
+ if (!finite(p2.norm[0]+p2.norm[1]+p2.norm[2]+p2.coeff))
+ return MB_FAILURE;
+
+ // store plane
+ return moab()->tag_set_data( planeTag, &node, 1, &p2 );
+}
+
+MBErrorCode MBBSPTree::set_tree_box( MBEntityHandle root_handle,
+ const double box_min[3],
+ const double box_max[3] )
+{
+ const double box[6] = { box_min[0], box_min[1], box_min[2],
+ box_max[0], box_max[1], box_max[2] };
+ return moab()->tag_set_data( rootTag, &root_handle, 1, box );
+}
+
+MBErrorCode MBBSPTree::get_tree_box( MBEntityHandle root_handle,
+ double box_min_out[3],
+ double box_max_out[3] )
+{
+ double box[6];
+ MBErrorCode rval = moab()->tag_get_data( rootTag, &root_handle, 1, box );
+ box_min_out[0] = box[0]; box_min_out[1] = box[1]; box_min_out[2] = box[2];
+ box_max_out[0] = box[3]; box_max_out[1] = box[4]; box_max_out[2] = box[5];
+ return rval;
+}
+
+MBErrorCode MBBSPTree::create_tree( MBEntityHandle& root_handle )
+{
+ const double min[3] = { -HUGE_VAL, -HUGE_VAL, -HUGE_VAL };
+ const double max[3] = { HUGE_VAL, HUGE_VAL, HUGE_VAL };
+ return create_tree( min, max, root_handle );
+}
+
+MBErrorCode MBBSPTree::create_tree( const double box_min[3],
+ const double box_max[3],
+ MBEntityHandle& root_handle )
+{
+ MBErrorCode rval = moab()->create_meshset( meshSetFlags, root_handle );
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ rval = set_tree_box( root_handle, box_min, box_max );
+ if (MB_SUCCESS != rval) {
+ moab()->delete_entities( &root_handle, 1 );
+ root_handle = 0;
+ return rval;
+ }
+
+ return MB_SUCCESS;
+}
+
+MBErrorCode MBBSPTree::delete_tree( MBEntityHandle root_handle )
+{
+ MBErrorCode rval;
+
+ std::vector<MBEntityHandle> children, dead_sets, current_sets;
+ current_sets.push_back( root_handle );
+ while (!current_sets.empty()) {
+ MBEntityHandle set = current_sets.back();
+ current_sets.pop_back();
+ dead_sets.push_back( set );
+ rval = moab()->get_child_meshsets( set, children );
+ if (MB_SUCCESS != rval)
+ return rval;
+ std::copy( children.begin(), children.end(), std::back_inserter(current_sets) );
+ children.clear();
+ }
+
+ rval = moab()->tag_delete_data( rootTag, &root_handle, 1 );
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ return moab()->delete_entities( &dead_sets[0], dead_sets.size() );
+}
+
+MBErrorCode MBBSPTree::find_all_trees( MBRange& results )
+{
+ return moab()->get_entities_by_type_and_tag( 0, MBENTITYSET,
+ &rootTag, 0, 1,
+ results );
+}
+
+MBErrorCode MBBSPTree::get_tree_iterator( MBEntityHandle root,
+ MBBSPTreeIter& iter )
+{
+ MBErrorCode rval = iter.initialize( this, root );
+ if (MB_SUCCESS != rval)
+ return rval;
+ return iter.step_to_first_leaf( MBBSPTreeIter::LEFT );
+}
+
+MBErrorCode MBBSPTree::get_tree_end_iterator( MBEntityHandle root,
+ MBBSPTreeIter& iter )
+{
+ MBErrorCode rval = iter.initialize( this, root );
+ if (MB_SUCCESS != rval)
+ return rval;
+ return iter.step_to_first_leaf( MBBSPTreeIter::RIGHT );
+}
+
+MBErrorCode MBBSPTree::split_leaf( MBBSPTreeIter& leaf,
+ Plane plane,
+ MBEntityHandle& left,
+ MBEntityHandle& right )
+{
+ MBErrorCode rval;
+
+ rval = moab()->create_meshset( meshSetFlags, left );
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ rval = moab()->create_meshset( meshSetFlags, right );
+ if (MB_SUCCESS != rval) {
+ moab()->delete_entities( &left, 1 );
+ return rval;
+ }
+
+ if (MB_SUCCESS != set_split_plane( leaf.handle(), plane ) ||
+ MB_SUCCESS != moab()->add_child_meshset( leaf.handle(), left ) ||
+ MB_SUCCESS != moab()->add_child_meshset( leaf.handle(), right) ||
+ MB_SUCCESS != leaf.step_to_first_leaf(MBBSPTreeIter::LEFT)) {
+ MBEntityHandle children[] = { left, right };
+ moab()->delete_entities( children, 2 );
+ return MB_FAILURE;
+ }
+
+ return MB_SUCCESS;
+}
+
+MBErrorCode MBBSPTree::split_leaf( MBBSPTreeIter& leaf, Plane plane )
+{
+ MBEntityHandle left, right;
+ return split_leaf( leaf, plane, left, right );
+}
+
+MBErrorCode MBBSPTree::split_leaf( MBBSPTreeIter& leaf,
+ Plane plane,
+ const MBRange& left_entities,
+ const MBRange& right_entities )
+{
+ MBEntityHandle left, right, parent = leaf.handle();
+ MBErrorCode rval = split_leaf( leaf, plane, left, right );
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ if (MB_SUCCESS == moab()->add_entities( left, left_entities ) &&
+ MB_SUCCESS == moab()->add_entities(right,right_entities ) &&
+ MB_SUCCESS == moab()->clear_meshset( &parent, 1 ))
+ return MB_SUCCESS;
+
+ moab()->remove_child_meshset( parent, left );
+ moab()->remove_child_meshset( parent, right );
+ MBEntityHandle children[] = { left, right };
+ moab()->delete_entities( children, 2 );
+ return MB_FAILURE;
+}
+
+MBErrorCode MBBSPTree::split_leaf( MBBSPTreeIter& leaf, Plane plane,
+ const std::vector<MBEntityHandle>& left_entities,
+ const std::vector<MBEntityHandle>& right_entities )
+{
+ MBEntityHandle left, right, parent = leaf.handle();
+ MBErrorCode rval = split_leaf( leaf, plane, left, right );
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ if (MB_SUCCESS == moab()->add_entities( left, &left_entities[0], left_entities.size() ) &&
+ MB_SUCCESS == moab()->add_entities(right,&right_entities[0],right_entities.size() ) &&
+ MB_SUCCESS == moab()->clear_meshset( &parent, 1 ))
+ return MB_SUCCESS;
+
+ moab()->remove_child_meshset( parent, left );
+ moab()->remove_child_meshset( parent, right );
+ MBEntityHandle children[] = { left, right };
+ moab()->delete_entities( children, 2 );
+ return MB_FAILURE;
+}
+
+MBErrorCode MBBSPTree::merge_leaf( MBBSPTreeIter& iter )
+{
+ MBErrorCode rval;
+ if (iter.depth() == 1) // at root
+ return MB_FAILURE;
+
+ // Move iter to parent
+ iter.up();
+
+ // Get all entities from children and put them in parent
+ MBEntityHandle parent = iter.handle();
+ moab()->remove_child_meshset( parent, iter.childVect[0] );
+ moab()->remove_child_meshset( parent, iter.childVect[1] );
+ std::vector<MBEntityHandle> stack( iter.childVect );
+
+ MBRange range;
+ while (!stack.empty()) {
+ MBEntityHandle h = stack.back();
+ stack.pop_back();
+ range.clear();
+ rval = moab()->get_entities_by_handle( h, range );
+ if (MB_SUCCESS != rval)
+ return rval;
+ rval = moab()->add_entities( parent, range );
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ iter.childVect.clear();
+ moab()->get_child_meshsets( h, iter.childVect );
+ if (!iter.childVect.empty()) {
+ moab()->remove_child_meshset( h, iter.childVect[0] );
+ moab()->remove_child_meshset( h, iter.childVect[1] );
+ stack.push_back( iter.childVect[0] );
+ stack.push_back( iter.childVect[1] );
+ }
+
+ rval = moab()->delete_entities( &h, 1 );
+ if (MB_SUCCESS != rval)
+ return rval;
+ }
+
+ return MB_SUCCESS;
+}
+
+
+
+MBErrorCode MBBSPTreeIter::initialize( MBBSPTree* tool,
+ MBEntityHandle root,
+ const double* point )
+{
+ treeTool = tool;
+ mStack.clear();
+ mStack.push_back( root );
+ return MB_SUCCESS;
+}
+
+
+MBErrorCode MBBSPTreeIter::step_to_first_leaf( Direction direction )
+{
+ MBErrorCode rval;
+ for (;;) {
+ childVect.clear();
+ rval = tool()->moab()->get_child_meshsets( mStack.back(), childVect );
+ if (MB_SUCCESS != rval)
+ return rval;
+ if (childVect.empty()) // leaf
+ break;
+
+ mStack.push_back( childVect[direction] );
+ }
+ return MB_SUCCESS;
+}
+
+MBErrorCode MBBSPTreeIter::step( Direction direction )
+{
+ MBEntityHandle node, parent;
+ MBErrorCode rval;
+ const Direction opposite = static_cast<Direction>(1-direction);
+
+ // If stack is empty, then either this iterator is uninitialized
+ // or we reached the end of the iteration (and return
+ // MB_ENTITY_NOT_FOUND) already.
+ if (mStack.empty())
+ return MB_FAILURE;
+
+ // Pop the current node from the stack.
+ // The stack should then contain the parent of the current node.
+ // If the stack is empty after this pop, then we've reached the end.
+ node = mStack.back();
+ mStack.pop_back();
+
+ while(!mStack.empty()) {
+ // Get data for parent entity
+ parent = mStack.back();
+ childVect.clear();
+ rval = tool()->moab()->get_child_meshsets( parent, childVect );
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ // If we're at the left child
+ if (childVect[opposite] == node) {
+ // push right child on stack
+ mStack.push_back( childVect[direction] );
+ // descend to left-most leaf of the right child
+ return step_to_first_leaf(opposite);
+ }
+
+ // The current node is the right child of the parent,
+ // continue up the tree.
+ assert( childVect[direction] == node );
+ node = parent;
+ mStack.pop_back();
+ }
+
+ return MB_ENTITY_NOT_FOUND;
+}
+
+MBErrorCode MBBSPTreeIter::up()
+{
+ if (mStack.size() < 2)
+ return MB_ENTITY_NOT_FOUND;
+ mStack.pop_back();
+ return MB_SUCCESS;
+}
+
+MBErrorCode MBBSPTreeIter::down( const MBBSPTree::Plane& plane, Direction dir )
+{
+ childVect.clear();
+ MBErrorCode rval = tool()->moab()->get_child_meshsets( mStack.back(), childVect );
+ if (MB_SUCCESS != rval)
+ return rval;
+ if (childVect.empty())
+ return MB_ENTITY_NOT_FOUND;
+
+ mStack.push_back( childVect[dir] );
+ return MB_SUCCESS;
+}
+
+MBErrorCode MBBSPTreeIter::get_parent_split_plane( MBBSPTree::Plane& plane ) const
+{
+ if (mStack.size() < 2) // at tree root
+ return MB_ENTITY_NOT_FOUND;
+
+ MBEntityHandle parent = mStack[mStack.size()-2];
+ return tool()->get_split_plane( parent, plane );
+}
+
+MBErrorCode MBBSPTreeBoxIter::initialize( MBBSPTree* tool_ptr,
+ MBEntityHandle root,
+ const double* point )
+{
+ MBErrorCode rval = MBBSPTreeIter::initialize( tool_ptr, root );
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ double box_min[3], box_max[3];
+ tool()->get_tree_box( root, box_min, box_max );
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ if (point) {
+ if (point[0] < box_min[0] || point[0] > box_max[0]
+ || point[1] < box_min[1] || point[1] > box_max[1]
+ || point[2] < box_min[2] || point[2] > box_max[2])
+ return MB_ENTITY_NOT_FOUND;
+ }
+
+ const double* ranges[] = { box_min, box_max };
+ for (int z = 0; z < 2; ++z) {
+ leafCoords[4*z ][0] = box_min[0];
+ leafCoords[4*z ][1] = box_min[1];
+ leafCoords[4*z ][2] = ranges[z][2];
+
+ leafCoords[4*z+1][0] = box_max[0];
+ leafCoords[4*z+1][1] = box_min[1];
+ leafCoords[4*z+1][2] = ranges[z][2];
+
+ leafCoords[4*z+2][0] = box_max[0];
+ leafCoords[4*z+2][1] = box_max[1];
+ leafCoords[4*z+2][2] = ranges[z][2];
+
+ leafCoords[4*z+3][0] = box_min[0];
+ leafCoords[4*z+3][1] = box_max[1];
+ leafCoords[4*z+3][2] = ranges[z][2];
+ }
+ stackData.resize(1);
+ return MB_SUCCESS;
+}
+
+MBBSPTreeBoxIter::SideBits
+MBBSPTreeBoxIter::side_above_plane( const double hex_coords[8][3],
+ const MBBSPTree::Plane& plane )
+{
+ unsigned result = 0;
+ for (unsigned i = 0; i < 8u; ++i)
+ result |= plane.above(hex_coords[i]) << i;
+ return (MBBSPTreeBoxIter::SideBits)result;
+}
+
+MBBSPTreeBoxIter::SideBits
+MBBSPTreeBoxIter::side_on_plane( const double hex_coords[8][3],
+ const MBBSPTree::Plane& plane )
+{
+ unsigned result = 0;
+ for (unsigned i = 0; i < 8u; ++i) {
+ bool on = plane.distance(hex_coords[i]) <= MBBSPTree::epsilon();
+ result |= on << i;
+ }
+ return (MBBSPTreeBoxIter::SideBits)result;
+}
+
+/** \brief Clip an edge using a plane
+ *
+ * Given an edge from keep_end_coords to cut_end_coords,
+ * cut the edge using the passed plane, such that cut_end_coords
+ * is updated with a new location on the plane, and old_coords_out
+ * contains the original value of cut_end_coords.
+ */
+static inline
+void plane_cut_edge( double old_coords_out[3],
+ const double keep_end_coords[3],
+ double cut_end_coords[3],
+ const MBBSPTree::Plane& plane )
+{
+ const MBCartVect start( keep_end_coords ), end( cut_end_coords );
+ const MBCartVect norm( plane.norm );
+ MBCartVect xsect_point;
+
+ const MBCartVect m = end - start;
+ const double t = -(norm % start + plane.coeff) / (norm % m);
+ assert( t > 0.0 && t < 1.0 );
+ xsect_point = start + t * m;
+
+ end.get( old_coords_out );
+ xsect_point.get( cut_end_coords );
+}
+
+/** Given the corners of a hexahedron in corners_input and a
+ * plane, cut the hex with the plane, updating corners_input
+ * and storing the original,cut-off side of the hex in cut_face_out.
+ *
+ * The portion of the hex below the plane is retained. cut_face_out
+ * will contain the side of the hex that is entirely above the plane.
+ *\return MB_FAILURE if plane/hex intersection is not a quadrilateral.
+ */
+static MBErrorCode plane_cut_box( double cut_face_out[4][3],
+ double corners_inout[8][3],
+ const MBBSPTree::Plane& plane )
+{
+ switch (MBBSPTreeBoxIter::side_above_plane( corners_inout, plane )) {
+ case MBBSPTreeBoxIter::B0154:
+ plane_cut_edge( cut_face_out[0], corners_inout[3], corners_inout[0], plane );
+ plane_cut_edge( cut_face_out[1], corners_inout[2], corners_inout[1], plane );
+ plane_cut_edge( cut_face_out[2], corners_inout[6], corners_inout[5], plane );
+ plane_cut_edge( cut_face_out[3], corners_inout[7], corners_inout[4], plane );
+ break;
+ case MBBSPTreeBoxIter::B1265:
+ plane_cut_edge( cut_face_out[0], corners_inout[0], corners_inout[1], plane );
+ plane_cut_edge( cut_face_out[1], corners_inout[3], corners_inout[2], plane );
+ plane_cut_edge( cut_face_out[2], corners_inout[7], corners_inout[6], plane );
+ plane_cut_edge( cut_face_out[3], corners_inout[4], corners_inout[5], plane );
+ break;
+ case MBBSPTreeBoxIter::B2376:
+ plane_cut_edge( cut_face_out[0], corners_inout[1], corners_inout[2], plane );
+ plane_cut_edge( cut_face_out[1], corners_inout[0], corners_inout[3], plane );
+ plane_cut_edge( cut_face_out[2], corners_inout[4], corners_inout[7], plane );
+ plane_cut_edge( cut_face_out[3], corners_inout[5], corners_inout[6], plane );
+ break;
+ case MBBSPTreeBoxIter::B3047:
+ plane_cut_edge( cut_face_out[0], corners_inout[2], corners_inout[3], plane );
+ plane_cut_edge( cut_face_out[1], corners_inout[1], corners_inout[0], plane );
+ plane_cut_edge( cut_face_out[2], corners_inout[5], corners_inout[4], plane );
+ plane_cut_edge( cut_face_out[3], corners_inout[6], corners_inout[7], plane );
+ break;
+ case MBBSPTreeBoxIter::B3210:
+ plane_cut_edge( cut_face_out[0], corners_inout[7], corners_inout[3], plane );
+ plane_cut_edge( cut_face_out[1], corners_inout[6], corners_inout[2], plane );
+ plane_cut_edge( cut_face_out[2], corners_inout[5], corners_inout[1], plane );
+ plane_cut_edge( cut_face_out[3], corners_inout[4], corners_inout[0], plane );
+ break;
+ case MBBSPTreeBoxIter::B4567:
+ plane_cut_edge( cut_face_out[0], corners_inout[0], corners_inout[4], plane );
+ plane_cut_edge( cut_face_out[1], corners_inout[1], corners_inout[5], plane );
+ plane_cut_edge( cut_face_out[2], corners_inout[2], corners_inout[6], plane );
+ plane_cut_edge( cut_face_out[3], corners_inout[3], corners_inout[7], plane );
+ break;
+ default:
+ return MB_FAILURE; // child is not a box
+ }
+
+ return MB_SUCCESS;
+}
+
+static inline
+void copy_coords( double dest[3], const double source[3] )
+{
+ dest[0] = source[0];
+ dest[1] = source[1];
+ dest[2] = source[2];
+}
+
+/** reverse of plane_cut_box */
+static inline
+MBErrorCode plane_uncut_box( const double cut_face_in[4][3],
+ double corners_inout[8][3],
+ const MBBSPTree::Plane& plane )
+{
+ switch (MBBSPTreeBoxIter::side_on_plane( corners_inout, plane )) {
+ case MBBSPTreeBoxIter::B0154:
+ copy_coords( corners_inout[0], cut_face_in[0] );
+ copy_coords( corners_inout[1], cut_face_in[1] );
+ copy_coords( corners_inout[5], cut_face_in[2] );
+ copy_coords( corners_inout[4], cut_face_in[3] );
+ break;
+ case MBBSPTreeBoxIter::B1265:
+ copy_coords( corners_inout[1], cut_face_in[0] );
+ copy_coords( corners_inout[2], cut_face_in[1] );
+ copy_coords( corners_inout[6], cut_face_in[2] );
+ copy_coords( corners_inout[5], cut_face_in[3] );
+ break;
+ case MBBSPTreeBoxIter::B2376:
+ copy_coords( corners_inout[2], cut_face_in[0] );
+ copy_coords( corners_inout[3], cut_face_in[1] );
+ copy_coords( corners_inout[7], cut_face_in[2] );
+ copy_coords( corners_inout[6], cut_face_in[3] );
+ break;
+ case MBBSPTreeBoxIter::B3047:
+ copy_coords( corners_inout[3], cut_face_in[0] );
+ copy_coords( corners_inout[0], cut_face_in[1] );
+ copy_coords( corners_inout[4], cut_face_in[2] );
+ copy_coords( corners_inout[7], cut_face_in[3] );
+ break;
+ case MBBSPTreeBoxIter::B3210:
+ copy_coords( corners_inout[3], cut_face_in[0] );
+ copy_coords( corners_inout[2], cut_face_in[1] );
+ copy_coords( corners_inout[1], cut_face_in[2] );
+ copy_coords( corners_inout[0], cut_face_in[3] );
+ break;
+ case MBBSPTreeBoxIter::B4567:
+ copy_coords( corners_inout[4], cut_face_in[0] );
+ copy_coords( corners_inout[5], cut_face_in[1] );
+ copy_coords( corners_inout[6], cut_face_in[2] );
+ copy_coords( corners_inout[7], cut_face_in[3] );
+ break;
+ default:
+ return MB_FAILURE; // child is not a box
+ }
+
+ return MB_SUCCESS;
+}
+
+MBErrorCode MBBSPTreeBoxIter::step_to_first_leaf( Direction direction )
+{
+ MBErrorCode rval;
+ MBBSPTree::Plane plane;
+ Corners clipped_corners;
+
+ for (;;) {
+ childVect.clear();
+ rval = tool()->moab()->get_child_meshsets( mStack.back(), childVect );
+ if (MB_SUCCESS != rval)
+ return rval;
+ if (childVect.empty()) // leaf
+ break;
+
+ rval = tool()->get_split_plane( mStack.back(), plane );
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ if (direction == RIGHT)
+ plane.flip();
+ rval = plane_cut_box( clipped_corners.coords, leafCoords, plane );
+ if (MB_SUCCESS != rval)
+ return rval;
+ mStack.push_back( childVect[direction] );
+ stackData.push_back( clipped_corners );
+ }
+ return MB_SUCCESS;
+}
+
+MBErrorCode MBBSPTreeBoxIter::up()
+{
+ MBErrorCode rval;
+ if (mStack.size() == 1)
+ return MB_ENTITY_NOT_FOUND;
+
+ MBEntityHandle node = mStack.back();
+ Corners clipped_face = stackData.back();
+ mStack.pop_back();
+ stackData.pop_back();
+
+ MBBSPTree::Plane plane;
+ rval = tool()->get_split_plane( mStack.back(), plane );
+ if (MB_SUCCESS != rval) {
+ mStack.push_back( node );
+ stackData.push_back( clipped_face );
+ return rval;
+ }
+
+ rval = plane_uncut_box( clipped_face.coords, leafCoords, plane );
+ if (MB_SUCCESS != rval) {
+ mStack.push_back( node );
+ stackData.push_back( clipped_face );
+ return rval;
+ }
+
+ return MB_SUCCESS;
+}
+
+MBErrorCode MBBSPTreeBoxIter::down( const MBBSPTree::Plane& plane_ref, Direction direction )
+{
+ childVect.clear();
+ MBErrorCode rval = tool()->moab()->get_child_meshsets( mStack.back(), childVect );
+ if (MB_SUCCESS != rval)
+ return rval;
+ if (childVect.empty())
+ return MB_ENTITY_NOT_FOUND;
+
+ MBBSPTree::Plane plane(plane_ref);
+ if (direction == RIGHT)
+ plane.flip();
+
+ Corners clipped_face;
+ rval = plane_cut_box( clipped_face.coords, leafCoords, plane );
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ mStack.push_back( childVect[direction] );
+ stackData.push_back( clipped_face );
+ return MB_SUCCESS;
+}
+
+MBErrorCode MBBSPTreeBoxIter::step( Direction direction )
+{
+ MBEntityHandle node, parent;
+ Corners clipped_face;
+ MBErrorCode rval;
+ MBBSPTree::Plane plane;
+ const Direction opposite = static_cast<Direction>(1-direction);
+
+ // If stack is empty, then either this iterator is uninitialized
+ // or we reached the end of the iteration (and return
+ // MB_ENTITY_NOT_FOUND) already.
+ if (mStack.empty())
+ return MB_FAILURE;
+
+ // Pop the current node from the stack.
+ // The stack should then contain the parent of the current node.
+ // If the stack is empty after this pop, then we've reached the end.
+ node = mStack.back();
+ mStack.pop_back();
+ clipped_face = stackData.back();
+ stackData.pop_back();
+
+ while(!mStack.empty()) {
+ // Get data for parent entity
+ parent = mStack.back();
+ childVect.clear();
+ rval = tool()->moab()->get_child_meshsets( parent, childVect );
+ if (MB_SUCCESS != rval)
+ return rval;
+ rval = tool()->get_split_plane( parent, plane );
+ if (MB_SUCCESS != rval)
+ return rval;
+ if (direction == LEFT)
+ plane.flip();
+
+ // If we're at the left child
+ if (childVect[opposite] == node) {
+ // change from box of left child to box of parent
+ plane_uncut_box( clipped_face.coords, leafCoords, plane );
+ // change from box of parent to box of right child
+ plane.flip();
+ plane_cut_box( clipped_face.coords, leafCoords, plane );
+ // push right child on stack
+ mStack.push_back( childVect[direction] );
+ stackData.push_back( clipped_face );
+ // descend to left-most leaf of the right child
+ return step_to_first_leaf(opposite);
+ }
+
+ // The current node is the right child of the parent,
+ // continue up the tree.
+ assert( childVect[direction] == node );
+ plane.flip();
+ plane_uncut_box( clipped_face.coords, leafCoords, plane );
+ node = parent;
+ clipped_face = stackData.back();
+ mStack.pop_back();
+ stackData.pop_back();
+ }
+
+ return MB_ENTITY_NOT_FOUND;
+}
+
+MBErrorCode MBBSPTreeBoxIter::get_box_corners( double coords[8][3] ) const
+{
+ memcpy( coords, leafCoords, 24*sizeof(double) );
+ return MB_SUCCESS;
+}
+
+MBErrorCode MBBSPTreeBoxIter::sibling_side( SideBits& side_out ) const
+{
+ if (mStack.size() < 2) // at tree root
+ return MB_ENTITY_NOT_FOUND;
+
+ MBEntityHandle parent = mStack[mStack.size()-2];
+ MBBSPTree::Plane plane;
+ MBErrorCode rval = tool()->get_split_plane( parent, plane );
+ if (MB_SUCCESS != rval)
+ return MB_FAILURE;
+
+ side_out = side_on_plane( leafCoords, plane );
+ return MB_SUCCESS;
+}
+
+MBErrorCode MBBSPTree::leaf_containing_point( MBEntityHandle tree_root,
+ const double point[3],
+ MBEntityHandle& leaf_out )
+{
+ std::vector<MBEntityHandle> children;
+ Plane plane;
+ MBEntityHandle node = tree_root;
+ MBErrorCode rval = moab()->get_child_meshsets( node, children );
+ if (MB_SUCCESS != rval)
+ return rval;
+ while (!children.empty()) {
+ rval = get_split_plane( node, plane );
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ node = children[plane.above(point)];
+ children.clear();
+ rval = moab()->get_child_meshsets( node, children );
+ if (MB_SUCCESS != rval)
+ return rval;
+ }
+ leaf_out = node;
+ return MB_SUCCESS;
+}
+
+MBErrorCode MBBSPTree::leaf_containing_point( MBEntityHandle root,
+ const double point[3],
+ MBBSPTreeIter& iter )
+{
+ MBErrorCode rval;
+
+ rval = iter.initialize( this, root, point );
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ for (;;) {
+ iter.childVect.clear();
+ rval = moab()->get_child_meshsets( iter.handle(), iter.childVect );
+ if (MB_SUCCESS != rval || iter.childVect.empty())
+ return rval;
+
+ Plane plane;
+ rval = get_split_plane( iter.handle(), plane );
+ if (MB_SUCCESS != rval)
+ return rval;
+
+ rval = iter.down( plane, (MBBSPTreeIter::Direction)(plane.above( point )) );
+ if (MB_SUCCESS != rval)
+ return rval;
+ }
+}
Added: MOAB/trunk/MBBSPTree.hpp
===================================================================
--- MOAB/trunk/MBBSPTree.hpp (rev 0)
+++ MOAB/trunk/MBBSPTree.hpp 2008-05-23 16:38:10 UTC (rev 1840)
@@ -0,0 +1,318 @@
+/*
+ * MOAB, a Mesh-Oriented datABase, is a software component for creating,
+ * storing and accessing finite element mesh data.
+ *
+ * Copyright 2004 Sandia Corporation. Under the terms of Contract
+ * DE-AC04-94AL85000 with Sandia Coroporation, the U.S. Government
+ * retains certain rights in this software.
+ *
+ * This library is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public
+ * License as published by the Free Software Foundation; either
+ * version 2.1 of the License, or (at your option) any later version.
+ *
+ */
+
+/**\file MBAdaptiveKDTree.hpp
+ *\author Jason Kraftcheck (kraftche at cae.wisc.edu)
+ *\date 2008-05-12
+ */
+
+#ifndef MB_BSP_TREE_HPP
+#define MB_BSP_TREE_HPP
+
+#include "MBTypes.h"
+#include "MBInterface.hpp"
+
+#include <math.h>
+#include <string>
+#include <vector>
+
+class MBBSPTreeBoxIter;
+class MBBSPTreeIter;
+class MBRange;
+
+class MBBSPTree
+{
+private:
+ MBInterface* mbInstance;
+ MBTag planeTag, rootTag;
+ unsigned meshSetFlags;
+
+public:
+
+ static double epsilon() { return 1e-6; }
+
+ MBBSPTree( MBInterface* iface,
+ const char* tagname = 0,
+ unsigned meshset_creation_flags = MESHSET_SET );
+
+ /**\brief struct to store a plane
+ *
+ * If plane is defined as Ax+By+Cz+D=0, then
+ * norm={A,B,C} and coeff=-D.
+ */
+ struct Plane {
+ Plane() {}
+ Plane( double n[3], double d ) : coeff(d)
+ { norm[0] = n[0]; norm[1] = n[1]; norm[2] = n[2]; }
+ /** a x + b y + c z + d = 0 */
+ Plane( double a, double b, double c, double d ) : coeff(d)
+ { norm[0] = a; norm[1] = b; norm[2] = c; }
+
+ double norm[3]; //!< Unit normal of plane
+ double coeff; //!< norm[0]*x + norm[1]*y + norm[2]*z + coeff = 0
+
+ /** Signed distance from point to plane:
+ * absolute value is distance from point to plane
+ * positive if 'above' plane, negative if 'below'
+ * Note: assumes unit-length normal.
+ */
+ double signed_distance( const double point[3] ) const
+ { return point[0]*norm[0]+point[1]*norm[1]+point[2]*norm[2] + coeff; }
+
+ /** return true if point is below the plane */
+ bool below( const double point[3] ) const
+ { return signed_distance(point) <= 0.0; }
+
+ /** return true if point is above the plane */
+ bool above( const double point[3] ) const
+ { return signed_distance(point) >= 0.0; }
+
+ double distance( const double point[3] ) const
+ { return fabs( signed_distance(point) ); }
+
+ /** reverse plane normal */
+ void flip()
+ { norm[0] = -norm[0]; norm[1] = -norm[1]; norm[2] = -norm[2]; coeff = -coeff; }
+ };
+
+ //! Get split plane for tree node
+ MBErrorCode get_split_plane( MBEntityHandle node, Plane& plane )
+ { return moab()->tag_get_data( planeTag, &node, 1, &plane ); }
+
+ //! Set split plane for tree node
+ MBErrorCode set_split_plane( MBEntityHandle node, const Plane& plane );
+
+ //! Get bounding box for entire tree
+ MBErrorCode get_tree_box( MBEntityHandle root_node,
+ double box_min_out[3],
+ double box_max_out[3] );
+
+ //! Set bounding box for entire tree
+ MBErrorCode set_tree_box( MBEntityHandle root_node,
+ const double box_min[3],
+ const double box_max[3] );
+
+ //! Create tree root node
+ MBErrorCode create_tree( const double box_min[3],
+ const double box_max[3],
+ MBEntityHandle& root_handle );
+
+ //! Create tree root node
+ MBErrorCode create_tree( MBEntityHandle& root_handle );
+
+ //! Find all tree roots
+ MBErrorCode find_all_trees( MBRange& results );
+
+ //! Destroy a tree
+ MBErrorCode delete_tree( MBEntityHandle root_handle );
+
+ MBInterface* moab() { return mbInstance; }
+
+ //! Get iterator for tree
+ MBErrorCode get_tree_iterator( MBEntityHandle tree_root,
+ MBBSPTreeIter& result );
+
+ //! Get iterator at right-most ('last') leaf.
+ MBErrorCode get_tree_end_iterator( MBEntityHandle tree_root,
+ MBBSPTreeIter& result );
+
+ //! Split leaf of tree
+ //! Updates iterator location to point to first new leaf node.
+ MBErrorCode split_leaf( MBBSPTreeIter& leaf, Plane plane );
+
+ //! Split leaf of tree
+ //! Updates iterator location to point to first new leaf node.
+ MBErrorCode split_leaf( MBBSPTreeIter& leaf,
+ Plane plane,
+ MBEntityHandle& left_child,
+ MBEntityHandle& right_child );
+
+ //! Split leaf of tree
+ //! Updates iterator location to point to first new leaf node.
+ MBErrorCode split_leaf( MBBSPTreeIter& leaf,
+ Plane plane,
+ const MBRange& left_entities,
+ const MBRange& right_entities );
+
+ //! Split leaf of tree
+ //! Updates iterator location to point to first new leaf node.
+ MBErrorCode split_leaf( MBBSPTreeIter& leaf,
+ Plane plane,
+ const std::vector<MBEntityHandle>& left_entities,
+ const std::vector<MBEntityHandle>& right_entities );
+
+ //! Merge the leaf pointed to by the current iterator with it's
+ //! sibling. If the sibling is not a leaf, multiple merges may
+ //! be done.
+ MBErrorCode merge_leaf( MBBSPTreeIter& iter );
+
+ //! Get leaf contianing input position.
+ //!
+ //! Does not take into account global bounding box of tree.
+ //! - Therefore there is always one leaf containing the point.
+ //! - If caller wants to account for global bounding box, then
+ //! caller can test against that box and not call this method
+ //! at all if the point is outside the box, as there is no leaf
+ //! containing the point in that case.
+ MBErrorCode leaf_containing_point( MBEntityHandle tree_root,
+ const double point[3],
+ MBEntityHandle& leaf_out );
+
+ //! Get iterator at leaf containing input position.
+ //!
+ //! Returns MB_ENTITY_NOT_FOUND if point is not within
+ //! bounding box of tree.
+ MBErrorCode leaf_containing_point( MBEntityHandle tree_root,
+ const double xyz[3],
+ MBBSPTreeIter& result );
+};
+
+class MBBSPTreeIter
+{
+public:
+
+ enum Direction { LEFT = 0, RIGHT = 1 };
+
+private:
+ friend class MBBSPTree;
+
+ MBBSPTree* treeTool;
+
+protected:
+ std::vector<MBEntityHandle> mStack;
+ mutable std::vector<MBEntityHandle> childVect;
+
+ virtual MBErrorCode step_to_first_leaf( Direction direction );
+
+ virtual MBErrorCode up();
+ virtual MBErrorCode down( const MBBSPTree::Plane& plane, Direction direction );
+
+ virtual MBErrorCode initialize( MBBSPTree* tool,
+ MBEntityHandle root,
+ const double* point = 0 );
+
+public:
+
+ MBBSPTreeIter() : treeTool(0), childVect(2) {}
+ virtual ~MBBSPTreeIter() {}
+
+
+ MBBSPTree* tool() const
+ { return treeTool; }
+
+ //! Get handle for current leaf
+ MBEntityHandle handle() const
+ { return mStack.back(); }
+
+ //! Get depth in tree. root is at depth of 1.
+ unsigned depth() const
+ { return mStack.size(); }
+
+ //! Advance the iterator either left or right in the tree
+ //! Note: stepping past the end of the tree will invalidate
+ //! the iterator. It will *not* be work step the
+ //! other direction.
+ virtual MBErrorCode step( Direction direction );
+
+ //! Advance to next leaf
+ //! Returns MB_ENTITY_NOT_FOUND if at end.
+ //! Note: steping past the end of the tree will invalidate
+ //! the iterator. Calling back() will not work.
+ MBErrorCode step() { return step(RIGHT); }
+
+ //! Move back to previous leaf
+ //! Returns MB_ENTITY_NOT_FOUND if at beginning.
+ //! Note: steping past the start of the tree will invalidate
+ //! the iterator. Calling step() will not work.
+ MBErrorCode back() { return step(LEFT); }
+
+
+ //! Get split plane that separates this node from its immediate sibling.
+ MBErrorCode get_parent_split_plane( MBBSPTree::Plane& plane ) const;
+};
+
+class MBBSPTreeBoxIter : public MBBSPTreeIter
+{
+ private:
+
+ double leafCoords[8][3];
+ struct Corners { double coords[4][3]; };
+ std::vector<Corners> stackData;
+
+ protected:
+
+ virtual MBErrorCode step_to_first_leaf( Direction direction );
+
+ virtual MBErrorCode up();
+ virtual MBErrorCode down( const MBBSPTree::Plane& plane, Direction direction );
+
+ virtual MBErrorCode initialize( MBBSPTree* tool,
+ MBEntityHandle root,
+ const double* point = 0 );
+ public:
+
+
+ //! Faces of a hex : corner bitmap
+ enum SideBits { B0154 = 0x33, //!< Face defined by corners {0,1,5,4}: 1st Exodus side
+ B1265 = 0x66, //!< Face defined by corners {1,2,6,5}: 2nd Exodus side
+ B2376 = 0xCC, //!< Face defined by corners {2,3,7,6}: 3rd Exodus side
+ B3047 = 0x99, //!< Face defined by corners {3,0,4,7}: 4th Exodus side
+ B3210 = 0x0F, //!< Face defined by corners {3,2,1,0}: 5th Exodus side
+ B4567 = 0xF0 //!< Face defined by corners {4,5,6,7}: 6th Exodus side
+ };
+
+ static SideBits side_above_plane( const double hex_coords[8][3],
+ const MBBSPTree::Plane& plane );
+
+ static SideBits side_on_plane( const double hex_coords[8][3],
+ const MBBSPTree::Plane& plane );
+
+ static SideBits opposite_face( const SideBits& bits )
+ { return (SideBits)(~bits); }
+
+ //! Advance the iterator either left or right in the tree
+ //! Note: stepping past the end of the tree will invalidate
+ //! the iterator. It will *not* be work step the
+ //! other direction.
+ virtual MBErrorCode step( Direction direction );
+
+ //! Advance to next leaf
+ //! Returns MB_ENTITY_NOT_FOUND if at end.
+ //! Note: steping past the end of the tree will invalidate
+ //! the iterator. Calling back() will not work.
+ MBErrorCode step() { return MBBSPTreeIter::step(); }
+
+ //! Move back to previous leaf
+ //! Returns MB_ENTITY_NOT_FOUND if at beginning.
+ //! Note: steping past the start of the tree will invalidate
+ //! the iterator. Calling step() will not work.
+ MBErrorCode back() { return MBBSPTreeIter::back(); }
+
+ //! Get coordinates of box corners, in Exodus II hexahedral ordering.
+ MBErrorCode get_box_corners( double coords[8][3] ) const;
+
+ //! Return the side of the box bounding this tree node
+ //! that is shared with the immediately adjacent sibling
+ //! (the tree node that shares a common parent node with
+ //! this node in the binary tree.)
+ //!
+ //!\return MB_ENTITY_NOT FOUND if root node.
+ //! MB_FAILURE if internal error.
+ //! MB_SUCCESS otherwise.
+ MBErrorCode sibling_side( SideBits& side_out ) const;
+
+};
+
+#endif
Modified: MOAB/trunk/Makefile.am
===================================================================
--- MOAB/trunk/Makefile.am 2008-05-22 21:52:58 UTC (rev 1839)
+++ MOAB/trunk/Makefile.am 2008-05-23 16:38:10 UTC (rev 1840)
@@ -28,6 +28,7 @@
obb_test \
vtk_test \
adaptive_kd_tree_tests \
+ bsp_tree_test \
file_options_test \
kd_tree_test \
var_len_test var_len_test_no_template \
@@ -119,6 +120,8 @@
MBAxisBox.hpp \
MBBits.cpp \
MBBits.hpp \
+ MBBSPTree.cpp \
+ MBBSPTree.hpp \
MBCN.cpp \
MBCNArrays.hpp \
MBCartVect.cpp \
@@ -333,6 +336,9 @@
kd_tree_test_SOURCES = kd_tree_test.cpp
kd_tree_test_LDADD = $(top_builddir)/libMOAB.la
+bsp_tree_test_SOURCES = bsp_tree_test.cpp
+bsp_tree_test_LDADD = $(top_builddir)/libMOAB.la
+
mbparallelcomm_test_SOURCES = mbparallelcomm_test.cpp
mbparallelcomm_test_LDADD = $(top_builddir)/libMOAB.la
mbparallelcomm_test_DEPENDENCIES = $(mbparallelcomm_test_LDADD)
Added: MOAB/trunk/bsp_tree_test.cpp
===================================================================
--- MOAB/trunk/bsp_tree_test.cpp (rev 0)
+++ MOAB/trunk/bsp_tree_test.cpp 2008-05-23 16:38:10 UTC (rev 1840)
@@ -0,0 +1,1103 @@
+#include "MBCore.hpp"
+#include "TestUtil.hpp"
+#include "MBBSPTree.hpp"
+
+void test_iterator();
+void test_box_iterator();
+void test_tree_create();
+void test_box_tree_create();
+void test_leaf_containing_point_bounded_tree();
+void test_leaf_containing_point_unbounded_tree();
+void test_merge_leaf();
+
+int main()
+{
+ int failures = 0;
+
+ failures += RUN_TEST( test_iterator );
+ failures += RUN_TEST( test_box_iterator );
+ failures += RUN_TEST( test_tree_create );
+ failures += RUN_TEST( test_box_tree_create );
+ failures += RUN_TEST( test_leaf_containing_point_bounded_tree );
+ failures += RUN_TEST( test_leaf_containing_point_unbounded_tree );
+ failures += RUN_TEST( test_merge_leaf );
+
+ return failures;
+}
+
+
+// Make CHECK_EQUAL macro work for planes
+void check_equal( const MBBSPTree::Plane& p1,
+ const MBBSPTree::Plane& p2,
+ const char* exp1,
+ const char* exp2,
+ int line,
+ const char* file )
+{
+ if (fabs(p1.norm[0]-p2.norm[0]) > 1e-6 ||
+ fabs(p1.norm[1]-p2.norm[1]) > 1e-6 ||
+ fabs(p1.norm[2]-p2.norm[2]) > 1e-6 ||
+ fabs(p1.coeff -p2.coeff ) > 1e-6) {
+ printf( "Equality Test Failed: %s == %s\n", exp1, exp2 );
+ printf( " at line %d of '%s'\n", line, file );
+ printf( " Expected: %f x + %f y + %f z + %f = 0\n",
+ p1.norm[0], p1.norm[1], p1.norm[2], p1.coeff );
+ printf( " Actual : %f x + %f y + %f z + %f = 0\n",
+ p2.norm[0], p2.norm[1], p2.norm[2], p2.coeff );
+ printf( "\n" );
+ FLAG_ERROR;
+ }
+}
+
+void test_iterator()
+{
+ MBCore moab;
+ MBBSPTree tool( &moab );
+ MBErrorCode rval;
+ MBEntityHandle root;
+ MBBSPTreeIter iter;
+
+ // create a depth-1 tree
+ rval = tool.create_tree( root );
+ CHECK_ERR(rval);
+ rval = tool.get_tree_iterator( root, iter );
+ CHECK_ERR(rval);
+
+ // check initial state of iterator
+ CHECK_EQUAL( &tool, iter.tool() );
+ CHECK_EQUAL( root, iter.handle() );
+ CHECK_EQUAL( 1u, iter.depth() );
+
+ // should fail if at root
+ MBBSPTree::Plane p;
+ rval = iter.get_parent_split_plane( p );
+ CHECK_EQUAL( MB_ENTITY_NOT_FOUND, rval );
+
+ // check that step past end returns correct value
+ rval = iter.step();
+ CHECK_EQUAL( MB_ENTITY_NOT_FOUND, rval );
+
+ // reset iterator
+ rval = tool.get_tree_iterator( root, iter );
+ CHECK_ERR(rval);
+
+ // check that step past start returns correct value
+ rval = iter.back();
+ CHECK_EQUAL( MB_ENTITY_NOT_FOUND, rval );
+
+ // reset iterator
+ rval = tool.get_tree_iterator( root, iter );
+ CHECK_ERR(rval);
+
+ // insert a single split plane
+ rval = tool.split_leaf( iter, MBBSPTree::Plane(2,0,0,0) );
+ CHECK_ERR(rval);
+
+ // check initial location
+ CHECK_EQUAL( 2u, iter.depth() );
+ CHECK( iter.handle() != root );
+
+ // create new iterators at left and right ends
+ MBBSPTreeIter left, right;
+ rval = tool.get_tree_iterator( root, left );
+ CHECK_ERR(rval);
+ rval = tool.get_tree_end_iterator( root, right );
+ CHECK_ERR(rval);
+
+ // compare post-split iterator to left
+ CHECK_EQUAL( left.depth(), iter.depth() );
+ CHECK_EQUAL( left.handle(), iter.handle() );
+
+ // step to other leaf
+ rval = iter.step();
+ CHECK_ERR(rval);
+
+ // check location
+ CHECK_EQUAL( 2u, iter.depth() );
+ CHECK( iter.handle() != root );
+
+ // compare stepped iterator to right
+ CHECK_EQUAL( right.depth(), iter.depth() );
+ CHECK_EQUAL( right.handle(), iter.handle() );
+
+ // step to back to first leaf
+ rval = iter.back();
+ CHECK_ERR(rval);
+
+ // compare to left
+ CHECK_EQUAL( left.depth(), iter.depth() );
+ CHECK_EQUAL( left.handle(), iter.handle() );
+
+ // check plane
+ // should have unit normal
+ left.get_parent_split_plane( p );
+ CHECK_EQUAL( MBBSPTree::Plane(1,0,0,0), p );
+ p.norm[0] = 11;
+ right.get_parent_split_plane( p );
+ CHECK_EQUAL( MBBSPTree::Plane(1,0,0,0), p );
+
+ // check step past ends
+ rval = left.back();
+ CHECK_EQUAL( MB_ENTITY_NOT_FOUND, rval );
+ rval = right.step();
+ CHECK_EQUAL( MB_ENTITY_NOT_FOUND, rval );
+}
+
+bool compare_hexes( const double expected[8][3],
+ const double actual[8][3],
+ double epsilon )
+{
+ // for each of three possible rotations
+ const int rotation_maps[3][8] = { { 0, 1, 2, 3, 4, 5, 6, 7 },
+ { 3, 2, 6, 7, 0, 1, 5, 4 },
+ { 7, 3, 2, 6, 4, 0, 1, 5 } };
+ for (int i = 0; i < 3; ++i) {
+ // compare for rotationts about axis from face 4 to 5
+ for (int j = 0; j < 4; ++j) {
+ bool match = true;
+ // for each face vertex
+ for (int k = 0; k < 4 && match; ++k) {
+ // for each coordinate
+ for (int d = 0; d < 3; ++d) {
+ if (fabs(expected[(j+k)%4 ][d] - actual[rotation_maps[i][k ]][d]) > epsilon
+ || fabs(expected[(j+k)%4+4][d] - actual[rotation_maps[i][k+4]][d]) > epsilon) {
+ match = false;
+ break;
+ }
+ }
+ }
+
+ if (match)
+ return true;
+ }
+ }
+
+ printf("Hex vertex copmarison failed.\n");
+ printf(" Exected Actual\n");
+ for (int v = 0; v < 8; ++v) {
+ printf("% 9f % 9f % 9f % 9f % 9f % 9f\n",
+ expected[v][0], expected[v][1], expected[v][2],
+ actual[v][0], actual[v][1], actual[v][2]);
+ }
+
+
+ return false;
+}
+
+static void aabox_corners( const double min[3],
+ const double max[3],
+ double corners[8][3] )
+{
+ const double expt[8][3] = { { min[0], min[1], min[2] },
+ { max[0], min[1], min[2] },
+ { max[0], max[1], min[2] },
+ { min[0], max[1], min[2] },
+ { min[0], min[1], max[2] },
+ { max[0], min[1], max[2] },
+ { max[0], max[1], max[2] },
+ { min[0], max[1], max[2] } };
+ memcpy( corners, expt, sizeof(expt) );
+}
+
+static void aabox_corners( double min_x, double min_y, double min_z,
+ double max_x, double max_y, double max_z,
+ double corners[8][3] )
+{
+ double min[] = { min_x, min_y, min_z };
+ double max[] = { max_x, max_y, max_z };
+ aabox_corners( min, max, corners );
+}
+
+
+void test_box_iterator()
+{
+ MBCore moab;
+ MBBSPTree tool( &moab );
+ MBErrorCode rval;
+ MBEntityHandle root;
+ MBBSPTreeBoxIter iter;
+ const double min[3] = { -1, -2, -3 };
+ const double max[3] = { 3, 2, 1 };
+
+ // create a depth-1 tree
+ rval = tool.create_tree( min, max, root );
+ CHECK_ERR(rval);
+ rval = tool.get_tree_iterator( root, iter );
+ CHECK_ERR(rval);
+
+ // check initial state of iterator
+ CHECK_EQUAL( &tool, iter.tool() );
+ CHECK_EQUAL( root, iter.handle() );
+ CHECK_EQUAL( 1u, iter.depth() );
+
+ // check initial corner values
+ double corners[8][3], expt[8][3];
+ aabox_corners( min, max, expt );
+ rval = iter.get_box_corners( corners );
+ CHECK_ERR( rval );
+ CHECK( compare_hexes( expt, corners, 1e-6 ) );
+
+ // should fail if at root
+ MBBSPTree::Plane p;
+ rval = iter.get_parent_split_plane( p );
+ CHECK_EQUAL( MB_ENTITY_NOT_FOUND, rval );
+
+ // check that step past end returns correct value
+ rval = iter.step();
+ CHECK_EQUAL( MB_ENTITY_NOT_FOUND, rval );
+
+ // reset iterator
+ rval = tool.get_tree_iterator( root, iter );
+ CHECK_ERR(rval);
+
+ // check that step past start returns correct value
+ rval = iter.back();
+ CHECK_EQUAL( MB_ENTITY_NOT_FOUND, rval );
+
+ // reset iterator
+ rval = tool.get_tree_iterator( root, iter );
+ CHECK_ERR(rval);
+
+ // insert a single split plane
+ rval = tool.split_leaf( iter, MBBSPTree::Plane(2,0,0,0) );
+ CHECK_ERR(rval);
+
+ // check initial location
+ CHECK_EQUAL( 2u, iter.depth() );
+ CHECK( iter.handle() != root );
+
+ // create new iterators at left and right ends
+ MBBSPTreeIter left, right;
+ rval = tool.get_tree_iterator( root, left );
+ CHECK_ERR(rval);
+ rval = tool.get_tree_end_iterator( root, right );
+ CHECK_ERR(rval);
+
+ // compare post-split iterator to left
+ CHECK_EQUAL( left.depth(), iter.depth() );
+ CHECK_EQUAL( left.handle(), iter.handle() );
+
+ // check box
+ aabox_corners( min, max, expt );
+ for (int i = 0; i < 8; ++i)
+ if (expt[i][0] > 0)
+ expt[i][0] = 0;
+ rval = iter.get_box_corners( corners );
+ CHECK_ERR( rval );
+ CHECK( compare_hexes( expt, corners, 1e-6 ) );
+
+ // step to other leaf
+ rval = iter.step();
+ CHECK_ERR(rval);
+
+ // check location
+ CHECK_EQUAL( 2u, iter.depth() );
+ CHECK( iter.handle() != root );
+
+ // compare stepped iterator to right
+ CHECK_EQUAL( right.depth(), iter.depth() );
+ CHECK_EQUAL( right.handle(), iter.handle() );
+
+ // check box
+ aabox_corners( min, max, expt );
+ for (int i = 0; i < 8; ++i)
+ if (expt[i][0] < 0)
+ expt[i][0] = 0;
+ rval = iter.get_box_corners( corners );
+ CHECK_ERR( rval );
+ CHECK( compare_hexes( expt, corners, 1e-6 ) );
+
+ // step to back to first leaf
+ rval = iter.back();
+ CHECK_ERR(rval);
+
+ // compare to left
+ CHECK_EQUAL( left.depth(), iter.depth() );
+ CHECK_EQUAL( left.handle(), iter.handle() );
+
+ // check box
+ aabox_corners( min, max, expt );
+ for (int i = 0; i < 8; ++i)
+ if (expt[i][0] > 0)
+ expt[i][0] = 0;
+ rval = iter.get_box_corners( corners );
+ CHECK_ERR( rval );
+ CHECK( compare_hexes( expt, corners, 1e-6 ) );
+
+ // check plane
+ // should have unit normal
+ left.get_parent_split_plane( p );
+ CHECK_EQUAL( MBBSPTree::Plane(1,0,0,0), p );
+ p.norm[0] = 11;
+ right.get_parent_split_plane( p );
+ CHECK_EQUAL( MBBSPTree::Plane(1,0,0,0), p );
+
+ // check step past ends
+ rval = left.back();
+ CHECK_EQUAL( MB_ENTITY_NOT_FOUND, rval );
+ rval = right.step();
+ CHECK_EQUAL( MB_ENTITY_NOT_FOUND, rval );
+}
+
+void test_tree_create()
+{
+ MBCore moab;
+ MBBSPTree tool( &moab );
+ MBErrorCode rval;
+ MBEntityHandle root;
+ MBBSPTreeIter iter;
+ MBBSPTree::Plane p;
+
+ // create a depth-1 tree
+ rval = tool.create_tree( root );
+ CHECK_ERR(rval);
+ rval = tool.get_tree_iterator( root, iter );
+ CHECK_ERR(rval);
+
+ // check initial state of iterator
+ CHECK_EQUAL( &tool, iter.tool() );
+ CHECK_EQUAL( root, iter.handle() );
+ CHECK_EQUAL( 1u, iter.depth() );
+
+ // split with z=0
+ rval = tool.split_leaf( iter, MBBSPTree::Plane(0,0,0.5,0) );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( 2u, iter.depth() );
+
+ // check that parent split plane is correct
+ rval = iter.get_parent_split_plane( p );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( MBBSPTree::Plane(0,0,1,0), p );
+
+ // split lower leaf with diagonal plane
+ rval = tool.split_leaf( iter, MBBSPTree::Plane(1,1,0,0) );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( 3u, iter.depth() );
+
+ const double r = 0.5*sqrt(2.0);
+
+ // check that parent split plane is correct
+ rval = iter.get_parent_split_plane( p );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( MBBSPTree::Plane(r,r,0,0), p );
+
+ // step to upper leaf
+ rval = iter.step();
+ CHECK_ERR(rval);
+
+ // split upper leaf with diagonal plane
+ rval = tool.split_leaf( iter, MBBSPTree::Plane(-1,1,0,0) );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( 4u, iter.depth() );
+
+ // check that parent split plane is correct
+ rval = iter.get_parent_split_plane( p );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( MBBSPTree::Plane(-r,r,0,0), p );
+
+ // iterate over four leaves
+
+ // first leaf
+ rval = tool.get_tree_iterator( root, iter );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( 3u, iter.depth() );
+ rval = iter.get_parent_split_plane( p );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( MBBSPTree::Plane(r,r,0,0), p );
+
+ // second leaf
+ rval = iter.step();
+ CHECK_ERR(rval);
+ CHECK_EQUAL( 4u, iter.depth() );
+ rval = iter.get_parent_split_plane( p );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( MBBSPTree::Plane(-r,r,0,0), p );
+
+ // third leaf
+ rval = iter.step();
+ CHECK_ERR(rval);
+ CHECK_EQUAL( 4u, iter.depth() );
+ rval = iter.get_parent_split_plane( p );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( MBBSPTree::Plane(-r,r,0,0), p );
+
+ // fourth leaf
+ rval = iter.step();
+ CHECK_ERR(rval);
+ CHECK_EQUAL( 2u, iter.depth() );
+ rval = iter.get_parent_split_plane( p );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( MBBSPTree::Plane(0,0,1,0), p );
+
+ // no more leaves
+ rval = iter.step();
+ CHECK_EQUAL( MB_ENTITY_NOT_FOUND, rval );
+}
+
+void test_box_tree_create()
+{
+ MBCore moab;
+ MBBSPTree tool( &moab );
+ MBErrorCode rval;
+ MBEntityHandle root;
+ MBBSPTreeBoxIter iter;
+ MBBSPTree::Plane p;
+ const double min[3] = { -2, -2, -2 };
+ const double max[3] = { 2, 2, 2 };
+
+ // create a depth-1 tree
+ rval = tool.create_tree( min, max, root );
+ CHECK_ERR(rval);
+ rval = tool.get_tree_iterator( root, iter );
+ CHECK_ERR(rval);
+
+ // check initial state of iterator
+ CHECK_EQUAL( &tool, iter.tool() );
+ CHECK_EQUAL( root, iter.handle() );
+ CHECK_EQUAL( 1u, iter.depth() );
+
+ // check initial corner values
+ double corners[8][3], expt[8][3];
+ aabox_corners( min, max, expt );
+ rval = iter.get_box_corners( corners );
+ CHECK_ERR( rval );
+ CHECK( compare_hexes( expt, corners, 1e-6 ) );
+
+
+ // Try splits that should fail because they
+ // do not intersect the leaf at all
+ rval = tool.split_leaf( iter, MBBSPTree::Plane(0,0,1, 4) );
+ CHECK_EQUAL( MB_FAILURE, rval );
+ rval = tool.split_leaf( iter, MBBSPTree::Plane(0,0,1,-4) );
+ CHECK_EQUAL( MB_FAILURE, rval );
+ rval = tool.split_leaf( iter, MBBSPTree::Plane(0,1,0, 4) );
+ CHECK_EQUAL( MB_FAILURE, rval );
+ rval = tool.split_leaf( iter, MBBSPTree::Plane(0,1,0,-4) );
+ CHECK_EQUAL( MB_FAILURE, rval );
+ rval = tool.split_leaf( iter, MBBSPTree::Plane(1,0,0, 4) );
+ CHECK_EQUAL( MB_FAILURE, rval );
+ rval = tool.split_leaf( iter, MBBSPTree::Plane(1,0,0,-4) );
+ CHECK_EQUAL( MB_FAILURE, rval );
+ rval = tool.split_leaf( iter, MBBSPTree::Plane(-1,-1,-1, 7) );
+ CHECK_EQUAL( MB_FAILURE, rval );
+ rval = tool.split_leaf( iter, MBBSPTree::Plane(-1,-1,-1,-7) );
+ CHECK_EQUAL( MB_FAILURE, rval );
+ rval = tool.split_leaf( iter, MBBSPTree::Plane( 1,-1,-1, 7) );
+ CHECK_EQUAL( MB_FAILURE, rval );
+ rval = tool.split_leaf( iter, MBBSPTree::Plane( 1,-1,-1,-7) );
+ CHECK_EQUAL( MB_FAILURE, rval );
+ rval = tool.split_leaf( iter, MBBSPTree::Plane( 1, 1,-1, 7) );
+ CHECK_EQUAL( MB_FAILURE, rval );
+ rval = tool.split_leaf( iter, MBBSPTree::Plane( 1, 1,-1,-7) );
+ CHECK_EQUAL( MB_FAILURE, rval );
+ rval = tool.split_leaf( iter, MBBSPTree::Plane(-1, 1,-1, 7) );
+ CHECK_EQUAL( MB_FAILURE, rval );
+ rval = tool.split_leaf( iter, MBBSPTree::Plane(-1, 1,-1,-7) );
+ CHECK_EQUAL( MB_FAILURE, rval );
+
+
+ // Try a split that should fail because the
+ // resulting leaf would not be a hexahedron.
+ // Clip each corner twice using planes with opposing normals
+ rval = tool.split_leaf( iter, MBBSPTree::Plane(-1,-1,-1,-4) );
+ CHECK_EQUAL( MB_FAILURE, rval );
+ rval = tool.split_leaf( iter, MBBSPTree::Plane( 1, 1, 1, 4) );
+ CHECK_EQUAL( MB_FAILURE, rval );
+ rval = tool.split_leaf( iter, MBBSPTree::Plane( 1,-1,-1,-4) );
+ CHECK_EQUAL( MB_FAILURE, rval );
+ rval = tool.split_leaf( iter, MBBSPTree::Plane(-1, 1, 1, 4) );
+ CHECK_EQUAL( MB_FAILURE, rval );
+ rval = tool.split_leaf( iter, MBBSPTree::Plane( 1, 1,-1,-4) );
+ CHECK_EQUAL( MB_FAILURE, rval );
+ rval = tool.split_leaf( iter, MBBSPTree::Plane(-1,-1, 1, 4) );
+ CHECK_EQUAL( MB_FAILURE, rval );
+ rval = tool.split_leaf( iter, MBBSPTree::Plane(-1, 1,-1,-4) );
+ CHECK_EQUAL( MB_FAILURE, rval );
+ rval = tool.split_leaf( iter, MBBSPTree::Plane( 1,-1, 1, 4) );
+ CHECK_EQUAL( MB_FAILURE, rval );
+ rval = tool.split_leaf( iter, MBBSPTree::Plane(-1,-1, 1,-4) );
+ CHECK_EQUAL( MB_FAILURE, rval );
+ rval = tool.split_leaf( iter, MBBSPTree::Plane( 1, 1,-1, 4) );
+ CHECK_EQUAL( MB_FAILURE, rval );
+ rval = tool.split_leaf( iter, MBBSPTree::Plane( 1,-1, 1,-4) );
+ CHECK_EQUAL( MB_FAILURE, rval );
+ rval = tool.split_leaf( iter, MBBSPTree::Plane(-1, 1,-1, 4) );
+ CHECK_EQUAL( MB_FAILURE, rval );
+ rval = tool.split_leaf( iter, MBBSPTree::Plane( 1, 1, 1,-4) );
+ CHECK_EQUAL( MB_FAILURE, rval );
+ rval = tool.split_leaf( iter, MBBSPTree::Plane(-1,-1,-1, 4) );
+ CHECK_EQUAL( MB_FAILURE, rval );
+ rval = tool.split_leaf( iter, MBBSPTree::Plane(-1, 1, 1,-4) );
+ CHECK_EQUAL( MB_FAILURE, rval );
+ rval = tool.split_leaf( iter, MBBSPTree::Plane( 1,-1,-1, 4) );
+ CHECK_EQUAL( MB_FAILURE, rval );
+ // Clip each edge
+ rval = tool.split_leaf( iter, MBBSPTree::Plane(-1,-1, 0,-2) );
+ CHECK_EQUAL( MB_FAILURE, rval );
+ rval = tool.split_leaf( iter, MBBSPTree::Plane( 1,-1, 0,-2) );
+ CHECK_EQUAL( MB_FAILURE, rval );
+ rval = tool.split_leaf( iter, MBBSPTree::Plane( 1, 1, 0,-2) );
+ CHECK_EQUAL( MB_FAILURE, rval );
+ rval = tool.split_leaf( iter, MBBSPTree::Plane(-1, 1, 0,-2) );
+ CHECK_EQUAL( MB_FAILURE, rval );
+ rval = tool.split_leaf( iter, MBBSPTree::Plane( 0,-1,-1,-2) );
+ CHECK_EQUAL( MB_FAILURE, rval );
+ rval = tool.split_leaf( iter, MBBSPTree::Plane( 0, 1,-1,-2) );
+ CHECK_EQUAL( MB_FAILURE, rval );
+ rval = tool.split_leaf( iter, MBBSPTree::Plane( 0, 1, 1,-2) );
+ CHECK_EQUAL( MB_FAILURE, rval );
+ rval = tool.split_leaf( iter, MBBSPTree::Plane( 0,-1, 1,-2) );
+ CHECK_EQUAL( MB_FAILURE, rval );
+ rval = tool.split_leaf( iter, MBBSPTree::Plane(-1, 0,-1,-2) );
+ CHECK_EQUAL( MB_FAILURE, rval );
+ rval = tool.split_leaf( iter, MBBSPTree::Plane( 1, 0,-1,-2) );
+ CHECK_EQUAL( MB_FAILURE, rval );
+ rval = tool.split_leaf( iter, MBBSPTree::Plane( 1, 0, 1,-2) );
+ CHECK_EQUAL( MB_FAILURE, rval );
+ rval = tool.split_leaf( iter, MBBSPTree::Plane(-1, 0, 1,-2) );
+ CHECK_EQUAL( MB_FAILURE, rval );
+
+ // verify that iterator is still valid after many failed splits
+ CHECK_EQUAL( &tool, iter.tool() );
+ CHECK_EQUAL( root, iter.handle() );
+ CHECK_EQUAL( 1u, iter.depth() );
+ rval = iter.get_box_corners( corners );
+ CHECK_ERR( rval );
+ CHECK( compare_hexes( expt, corners, 1e-6 ) );
+
+
+ // split with z=0
+ rval = tool.split_leaf( iter, MBBSPTree::Plane(0,0,0.5,0) );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( 2u, iter.depth() );
+
+ // check that parent split plane is correct
+ rval = iter.get_parent_split_plane( p );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( MBBSPTree::Plane(0,0,1,0), p );
+
+ // check that box corners are correct
+ aabox_corners( min, max, expt );
+ for (unsigned i = 0; i < 8; ++i)
+ if (expt[i][2] > 0.0)
+ expt[i][2] = 0.0;
+ rval = iter.get_box_corners( corners );
+ CHECK_ERR( rval );
+ CHECK( compare_hexes( expt, corners, 1e-6 ) );
+
+
+ // split with z=-1 and normal in opposite direction
+ rval = tool.split_leaf( iter, MBBSPTree::Plane(0,0,-2,-2) );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( 3u, iter.depth() );
+ for (unsigned i = 0; i < 8; ++i)
+ if (expt[i][2] < -1.0)
+ expt[i][2] = -1.0;
+ rval = iter.get_box_corners( corners );
+ CHECK_ERR( rval );
+ CHECK( compare_hexes( expt, corners, 1e-6 ) );
+
+ // step to next leaf (z < -1)
+ rval = iter.step();
+ CHECK_ERR(rval);
+ CHECK_EQUAL( 3u, iter.depth() );
+ aabox_corners( min, max, expt );
+ for (unsigned i = 0; i < 8; ++i)
+ if (expt[i][2] > -1.0)
+ expt[i][2] = -1.0;
+ rval = iter.get_box_corners( corners );
+ CHECK_ERR( rval );
+ CHECK( compare_hexes( expt, corners, 1e-6 ) );
+
+
+ // split at x = -1
+ rval = tool.split_leaf( iter, MBBSPTree::Plane(-0.1,0,0,-0.1) );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( 4u, iter.depth() );
+
+ // check that parent split plane is correct
+ rval = iter.get_parent_split_plane( p );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( MBBSPTree::Plane(-1,0,0,-1), p );
+
+ // check that leaf box is correct
+ aabox_corners( -1, -2, -2, 2, 2, -1, expt );
+ rval = iter.get_box_corners( corners );
+ CHECK_ERR( rval );
+ CHECK( compare_hexes( expt, corners, 1e-6 ) );
+
+
+ // split at x = 1
+ rval = tool.split_leaf( iter, MBBSPTree::Plane(5,0,0,-5) );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( 5u, iter.depth() );
+
+ // check that parent split plane is correct
+ rval = iter.get_parent_split_plane( p );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( MBBSPTree::Plane(1,0,0,-1), p );
+
+ // check that leaf box is correct
+ aabox_corners( -1, -2, -2, 1, 2, -1, expt );
+ rval = iter.get_box_corners( corners );
+ CHECK_ERR( rval );
+ CHECK( compare_hexes( expt, corners, 1e-6 ) );
+
+
+ // split at y = -1
+ rval = tool.split_leaf( iter, MBBSPTree::Plane(0,-1,0,-1) );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( 6u, iter.depth() );
+
+ // check that parent split plane is correct
+ rval = iter.get_parent_split_plane( p );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( MBBSPTree::Plane(0,-1,0,-1), p );
+
+ // check that leaf box is correct
+ aabox_corners( -1, -1, -2, 1, 2, -1, expt );
+ rval = iter.get_box_corners( corners );
+ CHECK_ERR( rval );
+ CHECK( compare_hexes( expt, corners, 1e-6 ) );
+
+
+ // split at y = 1
+ rval = tool.split_leaf( iter, MBBSPTree::Plane(0,1,0,-1) );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( 7u, iter.depth() );
+
+ // check that parent split plane is correct
+ rval = iter.get_parent_split_plane( p );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( MBBSPTree::Plane(0,1,0,-1), p );
+
+ // check that leaf box is correct
+ aabox_corners( -1, -1, -2, 1, 1, -1, expt );
+ rval = iter.get_box_corners( corners );
+ CHECK_ERR( rval );
+ CHECK( compare_hexes( expt, corners, 1e-6 ) );
+
+
+
+ // iterate over tree, verifying
+ rval = tool.get_tree_iterator( root, iter );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( 3u, iter.depth() );
+ rval = iter.get_parent_split_plane( p );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( MBBSPTree::Plane(0,0,-1,-1), p );
+ aabox_corners( -2, -2, -1, 2, 2, 0, expt );
+ rval = iter.get_box_corners( corners );
+ CHECK_ERR( rval );
+ CHECK( compare_hexes( expt, corners, 1e-6 ) );
+
+ rval = iter.step();
+ CHECK_ERR(rval);
+ CHECK_EQUAL( 7u, iter.depth() );
+ rval = iter.get_parent_split_plane( p );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( MBBSPTree::Plane(0,1,0,-1), p );
+ aabox_corners( -1, -1, -2, 1, 1, -1, expt );
+ rval = iter.get_box_corners( corners );
+ CHECK_ERR( rval );
+ CHECK( compare_hexes( expt, corners, 1e-6 ) );
+
+ rval = iter.step();
+ CHECK_ERR(rval);
+ CHECK_EQUAL( 7u, iter.depth() );
+ rval = iter.get_parent_split_plane( p );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( MBBSPTree::Plane(0,1,0,-1), p );
+ aabox_corners( -1, 1, -2, 1, 2, -1, expt );
+ rval = iter.get_box_corners( corners );
+ CHECK_ERR( rval );
+ CHECK( compare_hexes( expt, corners, 1e-6 ) );
+
+ rval = iter.step();
+ CHECK_ERR(rval);
+ CHECK_EQUAL( 6u, iter.depth() );
+ rval = iter.get_parent_split_plane( p );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( MBBSPTree::Plane(0,-1,0,-1), p );
+ aabox_corners( -1, -2, -2, 1, -1, -1, expt );
+ rval = iter.get_box_corners( corners );
+ CHECK_ERR( rval );
+ CHECK( compare_hexes( expt, corners, 1e-6 ) );
+
+ rval = iter.step();
+ CHECK_ERR(rval);
+ CHECK_EQUAL( 5u, iter.depth() );
+ rval = iter.get_parent_split_plane( p );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( MBBSPTree::Plane(1,0,0,-1), p );
+ aabox_corners( 1, -2, -2, 2, 2, -1, expt );
+ rval = iter.get_box_corners( corners );
+ CHECK_ERR( rval );
+ CHECK( compare_hexes( expt, corners, 1e-6 ) );
+
+ rval = iter.step();
+ CHECK_ERR(rval);
+ CHECK_EQUAL( 4u, iter.depth() );
+ rval = iter.get_parent_split_plane( p );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( MBBSPTree::Plane(-1,0,0,-1), p );
+ aabox_corners( -2, -2, -2, -1, 2, -1, expt );
+ rval = iter.get_box_corners( corners );
+ CHECK_ERR( rval );
+ CHECK( compare_hexes( expt, corners, 1e-6 ) );
+
+ rval = iter.step();
+ CHECK_ERR(rval);
+ CHECK_EQUAL( 2u, iter.depth() );
+ rval = iter.get_parent_split_plane( p );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( MBBSPTree::Plane(0,0,1,0), p );
+ aabox_corners( -2, -2, 0, 2, 2, 2, expt );
+ rval = iter.get_box_corners( corners );
+ CHECK_ERR( rval );
+ CHECK( compare_hexes( expt, corners, 1e-6 ) );
+
+ // no more leaves
+ rval = iter.step();
+ CHECK_EQUAL( MB_ENTITY_NOT_FOUND, rval );
+}
+
+void test_leaf_containing_point_bounded_tree()
+{
+ MBCore moab;
+ MBBSPTree tool( &moab );
+ MBErrorCode rval;
+ MBEntityHandle root;
+ MBBSPTreeIter iter;
+ MBBSPTreeBoxIter b_iter;
+ MBBSPTree::Plane p;
+ MBEntityHandle h;
+ double expected[8][3], corners[8][3];
+
+
+/* Build Tree
+
+ 1.0 +---------+--------------+
+ | A | |
+ | | |
+ 0.7 +---------+ C |
+ | | |
+ | | |
+ | B | |
+ | +--------------+ 0.3
+ | | D |
+ | | |
+ 0.0 +---------+--------------+
+ 0.0 0.4 1.0 */
+
+ const MBBSPTree::Plane X_split(1.0, 0.0, 0.0,-0.4);
+ const MBBSPTree::Plane AB_split(0.0,-1.0, 0.0, 0.7);
+ const MBBSPTree::Plane CD_split(0.0,-1.0, 0.0, 0.3);
+
+
+ const double min[3] = { 0, 0, 0 };
+ const double max[3] = { 1, 1, 1 };
+ rval = tool.create_tree( min, max, root );
+ CHECK_ERR(rval);
+ rval = tool.get_tree_iterator( root, iter );
+ CHECK_ERR(rval);
+
+ rval = tool.split_leaf( iter, X_split );
+ CHECK_ERR(rval);
+
+ rval = tool.split_leaf( iter, AB_split );
+ CHECK_ERR(rval);
+ const MBEntityHandle A = iter.handle();
+
+ rval = iter.step();
+ CHECK_ERR(rval);
+ const MBEntityHandle B = iter.handle();
+
+ rval = iter.step();
+ CHECK_ERR(rval);
+ rval = tool.split_leaf( iter, CD_split );
+ CHECK_ERR(rval);
+ const MBEntityHandle C = iter.handle();
+
+ rval = iter.step();
+ CHECK_ERR(rval);
+ const MBEntityHandle D = iter.handle();
+
+
+ // Test queries inside tree
+
+ const double A_point[] = { 0.1, 0.8, 0.5 };
+ rval = tool.leaf_containing_point( root, A_point, h );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( A, h );
+ rval = tool.leaf_containing_point( root, A_point, iter );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( A, iter.handle() );
+ rval = iter.get_parent_split_plane( p );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( AB_split, p );
+ rval = tool.leaf_containing_point( root, A_point, b_iter );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( A, b_iter.handle() );
+ aabox_corners( 0.0, 0.7, 0.0, 0.4, 1.0, 1.0, expected );
+ rval = b_iter.get_box_corners( corners );
+ CHECK_ERR( rval );
+ CHECK( compare_hexes( expected, corners, 1e-6 ) );
+
+ const double B_point[] = { 0.3, 0.1, 0.6 };
+ rval = tool.leaf_containing_point( root, B_point, h );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( B, h );
+ rval = tool.leaf_containing_point( root, B_point, iter );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( B, iter.handle() );
+ rval = iter.get_parent_split_plane( p );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( AB_split, p );
+ rval = tool.leaf_containing_point( root, B_point, b_iter );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( B, b_iter.handle() );
+ aabox_corners( 0.0, 0.0, 0.0, 0.4, 0.7, 1.0, expected );
+ rval = b_iter.get_box_corners( corners );
+ CHECK_ERR( rval );
+ CHECK( compare_hexes( expected, corners, 1e-6 ) );
+
+ const double C_point[] = { 0.9, 0.9, 0.1 };
+ rval = tool.leaf_containing_point( root, C_point, h );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( C, h );
+ rval = tool.leaf_containing_point( root, C_point, iter );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( C, iter.handle() );
+ rval = iter.get_parent_split_plane( p );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( CD_split, p );
+ rval = tool.leaf_containing_point( root, C_point, b_iter );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( C, b_iter.handle() );
+ aabox_corners( 0.4, 0.3, 0.0, 1.0, 1.0, 1.0, expected );
+ rval = b_iter.get_box_corners( corners );
+ CHECK_ERR( rval );
+ CHECK( compare_hexes( expected, corners, 1e-6 ) );
+
+ const double D_point[] = { 0.5, 0.2, 0.9 };
+ rval = tool.leaf_containing_point( root, D_point, h );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( D, h );
+ rval = tool.leaf_containing_point( root, D_point, iter );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( D, iter.handle() );
+ rval = iter.get_parent_split_plane( p );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( CD_split, p );
+ rval = tool.leaf_containing_point( root, D_point, b_iter );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( D, b_iter.handle() );
+ aabox_corners( 0.4, 0.0, 0.0, 1.0, 0.3, 1.0, expected );
+ rval = b_iter.get_box_corners( corners );
+ CHECK_ERR( rval );
+ CHECK( compare_hexes( expected, corners, 1e-6 ) );
+
+
+ // Try a couple points outside of the tree
+
+ const double above_pt[] = { 0.5, 0.5, 2.0 };
+ rval = tool.leaf_containing_point( root, above_pt, b_iter );
+ CHECK_EQUAL( MB_ENTITY_NOT_FOUND, rval );
+
+ const double left_pt[] = { -1.0, 0.5, 0.5 };
+ rval = tool.leaf_containing_point( root, left_pt, b_iter );
+ CHECK_EQUAL( MB_ENTITY_NOT_FOUND, rval );
+}
+
+void test_leaf_containing_point_unbounded_tree()
+{
+ MBCore moab;
+ MBBSPTree tool( &moab );
+ MBErrorCode rval;
+ MBEntityHandle root;
+ MBBSPTreeIter iter;
+ MBBSPTree::Plane p;
+ MBEntityHandle h;
+
+
+/* Build Tree
+
+ \ |
+ \ C o (0,4,0)
+ \ |
+ \ |
+ \ |
+ B \ | D
+ \|
+ ________o (0,0,0)
+ \
+ \
+ A \
+ o (1,-2,0)
+ \
+ */
+ const MBBSPTree::Plane X_split( 2.0, 1.0, 0.0, 0.0);
+ const MBBSPTree::Plane AB_split( 0.0, 1.0, 0.0, 0.0);
+ const MBBSPTree::Plane CD_split( 1.0, 0.0, 0.0, 0.0);
+
+
+ rval = tool.create_tree( root );
+ CHECK_ERR(rval);
+ rval = tool.get_tree_iterator( root, iter );
+ CHECK_ERR(rval);
+
+ rval = tool.split_leaf( iter, X_split );
+ CHECK_ERR(rval);
+
+ rval = tool.split_leaf( iter, AB_split );
+ CHECK_ERR(rval);
+ const MBEntityHandle A = iter.handle();
+
+ rval = iter.step();
+ CHECK_ERR(rval);
+ const MBEntityHandle B = iter.handle();
+
+ rval = iter.step();
+ CHECK_ERR(rval);
+ rval = tool.split_leaf( iter, CD_split );
+ CHECK_ERR(rval);
+ const MBEntityHandle C = iter.handle();
+
+ rval = iter.step();
+ CHECK_ERR(rval);
+ const MBEntityHandle D = iter.handle();
+
+
+ // Test queries inside tree
+
+ const double A_point[] = { -1000, -1000, -1000 };
+ rval = tool.leaf_containing_point( root, A_point, h );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( A, h );
+ rval = tool.leaf_containing_point( root, A_point, iter );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( A, iter.handle() );
+ rval = iter.get_parent_split_plane( p );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( AB_split, p );
+
+ const double B_point[] = { -3, 1, 100 };
+ rval = tool.leaf_containing_point( root, B_point, h );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( B, h );
+ rval = tool.leaf_containing_point( root, B_point, iter );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( B, iter.handle() );
+ rval = iter.get_parent_split_plane( p );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( AB_split, p );
+
+ const double C_point[] = { -10, 500, 0 };
+ rval = tool.leaf_containing_point( root, C_point, h );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( C, h );
+ rval = tool.leaf_containing_point( root, C_point, iter );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( C, iter.handle() );
+ rval = iter.get_parent_split_plane( p );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( CD_split, p );
+
+ const double D_point[] = { 10, 500, 0 };
+ rval = tool.leaf_containing_point( root, D_point, h );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( D, h );
+ rval = tool.leaf_containing_point( root, D_point, iter );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( D, iter.handle() );
+ rval = iter.get_parent_split_plane( p );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( CD_split, p );
+}
+
+void test_merge_leaf()
+{
+ MBCore moab;
+ MBBSPTree tool( &moab );
+ MBErrorCode rval;
+ MBEntityHandle root;
+ MBBSPTreeBoxIter iter;
+ MBBSPTree::Plane p;
+ double expected[8][3], corners[8][3];
+
+
+/* Build Tree
+
+ 1.0 +---------+--------------+
+ | A | |
+ | | |
+ 0.7 +---------+ C |
+ | | |
+ | | |
+ | B | |
+ | +--------------+ 0.3
+ | | D |
+ | | |
+ 0.0 +---------+--------------+
+ 0.0 0.4 1.0 */
+
+ const MBBSPTree::Plane X_split(1.0, 0.0, 0.0,-0.4);
+ const MBBSPTree::Plane AB_split(0.0,-1.0, 0.0, 0.7);
+ const MBBSPTree::Plane CD_split(0.0,-1.0, 0.0, 0.3);
+
+ const double min[3] = { 0, 0, 0 };
+ const double max[3] = { 1, 1, 1 };
+ rval = tool.create_tree( min, max, root );
+ CHECK_ERR(rval);
+ rval = tool.get_tree_iterator( root, iter );
+ CHECK_ERR(rval);
+ rval = tool.split_leaf( iter, X_split );
+ CHECK_ERR(rval);
+ const MBEntityHandle AB = iter.handle();
+ rval = tool.split_leaf( iter, AB_split );
+ CHECK_ERR(rval);
+ rval = iter.step();
+ CHECK_ERR(rval);
+ rval = iter.step();
+ CHECK_ERR(rval);
+ const MBEntityHandle CD = iter.handle();
+ rval = tool.split_leaf( iter, CD_split );
+ CHECK_ERR(rval);
+ rval = iter.step();
+ CHECK_ERR(rval);
+
+ rval = tool.get_tree_iterator( root, iter );
+ CHECK_ERR(rval);
+ rval = tool.merge_leaf( iter );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( AB, iter.handle() );
+ CHECK_EQUAL( 2u, iter.depth() );
+ rval = iter.get_parent_split_plane( p );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( X_split, p );
+ aabox_corners( 0.0, 0.0, 0.0, 0.4, 1.0, 1.0, expected );
+ rval = iter.get_box_corners( corners );
+ CHECK_ERR( rval );
+ CHECK( compare_hexes( expected, corners, 1e-6 ) );
+
+ rval = iter.step();
+ CHECK_ERR(rval);
+ rval = iter.step();
+ CHECK_ERR(rval);
+ rval = tool.merge_leaf( iter );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( CD, iter.handle() );
+ CHECK_EQUAL( 2u, iter.depth() );
+ rval = iter.get_parent_split_plane( p );
+ CHECK_ERR(rval);
+ CHECK_EQUAL( X_split, p );
+ aabox_corners( 0.4, 0.0, 0.0, 1.0, 1.0, 1.0, expected );
+ rval = iter.get_box_corners( corners );
+ CHECK_ERR( rval );
+ CHECK( compare_hexes( expected, corners, 1e-6 ) );
+
+ rval = iter.step();
+ CHECK_EQUAL( MB_ENTITY_NOT_FOUND, rval );
+}
+
Modified: MOAB/trunk/configure.in
===================================================================
--- MOAB/trunk/configure.in 2008-05-22 21:52:58 UTC (rev 1839)
+++ MOAB/trunk/configure.in 2008-05-23 16:38:10 UTC (rev 1840)
@@ -96,6 +96,12 @@
fi
################################################################################
+# System Headers
+################################################################################
+AC_CHECK_HEADER([ieeefp.h],[DEFINES="$DEFINES -DHAVE_IEEEFP_H"] )
+
+
+################################################################################
# BOOST OPTIONS
################################################################################
More information about the moab-dev
mailing list