[MOAB-dev] commit/MOAB: 2 new changesets
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Mon Jan 13 19:09:13 CST 2014
2 new commits in MOAB:
https://bitbucket.org/fathomteam/moab/commits/293272f5da39/
Changeset: 293272f5da39
Branch: None
User: tautges
Date: 2014-01-13 21:10:10
Summary: Changing is_inside in ElemEvaluation stuff to be int instead of bool (should have remembered that from Effective STL)
Adding BoundBox::get().
Initial addition of code for parallel spatial location using intermediate point location. Compiles but not tested yet.
Eliminating many compiler warnings from test/oldinc having to do with poorly-chosen variable names.
Passes make check in parallel.
Affected #: 29 files
diff --git a/.gitignore b/.gitignore
index d842c2b..482fcd4 100644
--- a/.gitignore
+++ b/.gitignore
@@ -155,6 +155,7 @@ test/io/gmsh_test
test/io/ideas_test
test/io/nastran_test
test/io/read_cgm_test
+test/io/read_mpas_nc
test/io/read_nc
test/io/read_ucd_nc
test/io/readutil_test
@@ -169,6 +170,7 @@ test/lloyd_smoother_test
test/mbcn_test
test/mbfacet_test
test/mbground_test
+test/mergemesh_test
test/mesh_set_test
test/moab_test
test/obb/obb_test
@@ -188,6 +190,7 @@ test/parallel/parallel_unit_tests
test/parallel/parallel_write_test
test/parallel/par_coupler_test
test/parallel/par_intx_sph
+test/parallel/par_spatial_locator_test
test/parallel/parmerge
test/parallel/partcheck
test/parallel/pcomm_serial
@@ -208,6 +211,8 @@ test/perf/seqperf
test/perf/tstt_perf_binding
test/perf/point_location/elem_eval_time
test/perf/point_location/point_location
+test/perf/point_location/sploc_searching_perf
+test/perf/point_location/tree_searching_perf
test/range_test
test/reorder_test
test/scdseq_test
@@ -243,6 +248,11 @@ tools/mbcslam/lagr.h5m
tools/mbcslam/spectral.vtk
tools/mbcslam/spec_visu_test
tools/mbcslam/spherical_area_test
+tools/mbcslam/cslam_par_test
+tools/mbcslam/create_dp
+tools/mbcslam/diffusion
+tools/mbcslam/intx_imesh
+tools/mbcslam/proj1
tools/mbdepth
tools/mbgsets
tools/mbmem
diff --git a/src/BVHTree.cpp b/src/BVHTree.cpp
index b65eae4..d54e24d 100644
--- a/src/BVHTree.cpp
+++ b/src/BVHTree.cpp
@@ -443,7 +443,7 @@ namespace moab
const TreeNode &node = myTree[index];
treeStats.nodesVisited++;
CartVect params;
- bool is_inside;
+ int is_inside;
ErrorCode rval = MB_SUCCESS;
if(node.dim == 3){
treeStats.leavesVisited++;
diff --git a/src/LocalDiscretization/ElemEvaluator.cpp b/src/LocalDiscretization/ElemEvaluator.cpp
index 4ca3943..209cd11 100644
--- a/src/LocalDiscretization/ElemEvaluator.cpp
+++ b/src/LocalDiscretization/ElemEvaluator.cpp
@@ -18,7 +18,7 @@ namespace moab {
ErrorCode EvalSet::evaluate_reverse(EvalFcn eval, JacobianFcn jacob, InsideFcn inside_f,
const double *posn, const double *verts, const int nverts,
const int ndim, const double iter_tol, const double inside_tol,
- double *work, double *params, bool *inside) {
+ double *work, double *params, int *inside) {
// TODO: should differentiate between epsilons used for
// Newton Raphson iteration, and epsilons used for curved boundary geometry errors
// right now, fix the tolerance used for NR
@@ -39,7 +39,7 @@ namespace moab {
// residual is diff between old and new pos; need to minimize that
CartVect res = new_pos - *cvposn;
Matrix3 J;
- bool dum, *tmp_inside = (inside ? inside : &dum);
+ int dum, *tmp_inside = (inside ? inside : &dum);
int iters=0;
// while |res| larger than tol
@@ -77,7 +77,7 @@ namespace moab {
return MB_SUCCESS;
}// Map::evaluate_reverse()
- bool EvalSet::inside_function(const double *params, const int ndims, const double tol)
+ int EvalSet::inside_function(const double *params, const int ndims, const double tol)
{
if (params[0] >= -1-tol && params[0] <= 1+tol &&
(ndims < 2 || (params[1] >= -1-tol && params[1] <= 1+tol)) &&
@@ -118,7 +118,7 @@ namespace moab {
const double inside_tol, EntityHandle &containing_ent,
double *params, unsigned int *num_evals)
{
- bool is_inside;
+ int is_inside;
ErrorCode rval = MB_SUCCESS;
unsigned int nevals = 0;
Range::iterator i;
diff --git a/src/LocalDiscretization/LinearHex.cpp b/src/LocalDiscretization/LinearHex.cpp
index 5e730af..e9044a8 100644
--- a/src/LocalDiscretization/LinearHex.cpp
+++ b/src/LocalDiscretization/LinearHex.cpp
@@ -97,14 +97,14 @@ namespace moab
ErrorCode LinearHex::reverseEvalFcn(EvalFcn eval, JacobianFcn jacob, InsideFcn ins,
const double *posn, const double *verts, const int nverts, const int ndim,
const double iter_tol, const double inside_tol, double *work,
- double *params, bool *is_inside)
+ double *params, int *is_inside)
{
assert(posn && verts);
return EvalSet::evaluate_reverse(eval, jacob, ins, posn, verts, nverts, ndim, iter_tol, inside_tol, work,
params, is_inside);
}
- bool LinearHex::insideFcn(const double *params, const int ndim, const double tol)
+ int LinearHex::insideFcn(const double *params, const int ndim, const double tol)
{
return EvalSet::inside_function(params, ndim, tol);
}
diff --git a/src/LocalDiscretization/LinearQuad.cpp b/src/LocalDiscretization/LinearQuad.cpp
index 63c9959..667bf0b 100644
--- a/src/LocalDiscretization/LinearQuad.cpp
+++ b/src/LocalDiscretization/LinearQuad.cpp
@@ -76,13 +76,13 @@ namespace moab
ErrorCode LinearQuad::reverseEvalFcn(EvalFcn eval, JacobianFcn jacob, InsideFcn ins,
const double *posn, const double *verts, const int nverts, const int ndim,
const double iter_tol, const double inside_tol, double *work,
- double *params, bool *is_inside)
+ double *params, int *is_inside)
{
return EvalSet::evaluate_reverse(eval, jacob, ins, posn, verts, nverts, ndim, iter_tol, inside_tol, work,
params, is_inside);
}
- bool LinearQuad::insideFcn(const double *params, const int ndim, const double tol)
+ int LinearQuad::insideFcn(const double *params, const int ndim, const double tol)
{
return EvalSet::inside_function(params, ndim, tol);
}
diff --git a/src/LocalDiscretization/LinearTet.cpp b/src/LocalDiscretization/LinearTet.cpp
index 4965cfc..6d590e7 100644
--- a/src/LocalDiscretization/LinearTet.cpp
+++ b/src/LocalDiscretization/LinearTet.cpp
@@ -75,14 +75,14 @@ namespace moab
ErrorCode LinearTet::reverseEvalFcn(EvalFcn eval, JacobianFcn jacob, InsideFcn ins,
const double *posn, const double *verts, const int nverts, const int ndim,
const double iter_tol, const double inside_tol, double *work,
- double *params, bool *is_inside)
+ double *params, int *is_inside)
{
assert(posn && verts);
return evaluate_reverse(eval, jacob, ins, posn, verts, nverts, ndim, iter_tol, inside_tol,
work, params, is_inside);
}
- bool LinearTet::insideFcn(const double *params, const int , const double tol)
+ int LinearTet::insideFcn(const double *params, const int , const double tol)
{
return (params[0] >= -1.0-tol && params[1] >= -1.0-tol && params[2] >= -1.0-tol &&
params[0] + params[1] + params[2] <= 1.0+tol);
@@ -92,7 +92,7 @@ namespace moab
ErrorCode LinearTet::evaluate_reverse(EvalFcn eval, JacobianFcn jacob, InsideFcn inside_f,
const double *posn, const double *verts, const int nverts,
const int ndim, const double iter_tol, const double inside_tol,
- double *work, double *params, bool *inside) {
+ double *work, double *params, int *inside) {
// TODO: should differentiate between epsilons used for
// Newton Raphson iteration, and epsilons used for curved boundary geometry errors
// right now, fix the tolerance used for NR
diff --git a/src/LocalDiscretization/LinearTri.cpp b/src/LocalDiscretization/LinearTri.cpp
index 8fe9d23..e6f7f1c 100644
--- a/src/LocalDiscretization/LinearTri.cpp
+++ b/src/LocalDiscretization/LinearTri.cpp
@@ -72,14 +72,14 @@ namespace moab
ErrorCode LinearTri::reverseEvalFcn(EvalFcn eval, JacobianFcn jacob, InsideFcn ins,
const double *posn, const double *verts, const int nverts, const int ndim,
const double iter_tol, const double inside_tol, double *work,
- double *params, bool *is_inside)
+ double *params, int *is_inside)
{
assert(posn && verts);
return evaluate_reverse(eval, jacob, ins, posn, verts, nverts, ndim, iter_tol, inside_tol, work,
params, is_inside);
}
- bool LinearTri::insideFcn(const double *params, const int , const double tol)
+ int LinearTri::insideFcn(const double *params, const int , const double tol)
{
return (params[0] >= -1.0-tol && params[1] >= -1.0-tol &&
params[0] + params[1] <= 1.0+tol);
@@ -89,7 +89,7 @@ namespace moab
ErrorCode LinearTri::evaluate_reverse(EvalFcn eval, JacobianFcn jacob, InsideFcn inside_f,
const double *posn, const double *verts, const int nverts,
const int ndim, const double iter_tol, const double inside_tol,
- double *work, double *params, bool *inside) {
+ double *work, double *params, int *inside) {
// TODO: should differentiate between epsilons used for
// Newton Raphson iteration, and epsilons used for curved boundary geometry errors
// right now, fix the tolerance used for NR
diff --git a/src/LocalDiscretization/QuadraticHex.cpp b/src/LocalDiscretization/QuadraticHex.cpp
index 5c621f6..e8d2ad2 100644
--- a/src/LocalDiscretization/QuadraticHex.cpp
+++ b/src/LocalDiscretization/QuadraticHex.cpp
@@ -108,14 +108,14 @@ namespace moab
ErrorCode QuadraticHex::reverseEvalFcn(EvalFcn eval, JacobianFcn jacob, InsideFcn ins,
const double *posn, const double *verts, const int nverts, const int ndim,
const double iter_tol, const double inside_tol, double *work,
- double *params, bool *is_inside)
+ double *params, int *is_inside)
{
assert(posn && verts);
return EvalSet::evaluate_reverse(eval, jacob, ins, posn, verts, nverts, ndim, iter_tol, inside_tol,
work, params, is_inside);
}
- bool QuadraticHex::insideFcn(const double *params, const int ndim, const double tol)
+ int QuadraticHex::insideFcn(const double *params, const int ndim, const double tol)
{
return EvalSet::inside_function(params, ndim, tol);
}
diff --git a/src/LocalDiscretization/SpectralHex.cpp b/src/LocalDiscretization/SpectralHex.cpp
index 31f5391..be383fa 100644
--- a/src/LocalDiscretization/SpectralHex.cpp
+++ b/src/LocalDiscretization/SpectralHex.cpp
@@ -180,7 +180,7 @@ void SpectralHex::integrate_vector(const double *field_values, int num_tuples, d
//std::cout << "volume: " << volume << "\n";
}
- bool SpectralHex::insideFcn(const double *params, const int ndim, const double tol)
+ int SpectralHex::insideFcn(const double *params, const int ndim, const double tol)
{
return EvalSet::inside(params, ndim, tol);
}
diff --git a/src/LocalDiscretization/SpectralQuad.cpp b/src/LocalDiscretization/SpectralQuad.cpp
index 5c11d13..d7355e6 100644
--- a/src/LocalDiscretization/SpectralQuad.cpp
+++ b/src/LocalDiscretization/SpectralQuad.cpp
@@ -96,7 +96,7 @@ CartVect SpectralQuad::evalFcn(const double *params, const double *field, const
// replicate the functionality of hex_findpt
bool SpectralQuad::reverseEvalFcn(const double *posn, const double *verts, const int nverts, const int ndim,
const double iter_tol, const double inside_tol, double *work,
- double *params, bool *is_inside)
+ double *params, int *is_inside)
{
params = init;
@@ -147,10 +147,10 @@ void SpectralQuad:: integrate_vector(const double *field, const double *verts, c
// not implemented
}
- bool SpectralQuad::insideFcn(const double *params, const int ndim, const double tol)
- {
- return EvalSet::inside(params, ndim, tol);
- }
+int SpectralQuad::insideFcn(const double *params, const int ndim, const double tol)
+{
+ return EvalSet::inside(params, ndim, tol);
+}
// something we don't do for spectral hex; we do it here because
// we do not store the position of gl points in a tag yet
diff --git a/src/LocalDiscretization/moab/ElemEvaluator.hpp b/src/LocalDiscretization/moab/ElemEvaluator.hpp
index 29d1051..866388c 100644
--- a/src/LocalDiscretization/moab/ElemEvaluator.hpp
+++ b/src/LocalDiscretization/moab/ElemEvaluator.hpp
@@ -21,12 +21,12 @@ namespace moab {
typedef ErrorCode (*InitFcn)(const double *verts, const int nverts, double *&work);
- typedef bool (*InsideFcn)(const double *verts, const int ndims, const double tol);
+ typedef int (*InsideFcn)(const double *verts, const int ndims, const double tol);
typedef ErrorCode (*ReverseEvalFcn)(EvalFcn eval, JacobianFcn jacob, InsideFcn ins,
const double *posn, const double *verts, const int nverts, const int ndim,
const double iter_tol, const double inside_tol,
- double *work, double *params, bool *is_inside);
+ double *work, double *params, int *is_inside);
class EvalSet
{
@@ -78,9 +78,9 @@ namespace moab {
static ErrorCode evaluate_reverse(EvalFcn eval, JacobianFcn jacob, InsideFcn inside_f,
const double *posn, const double *verts, const int nverts,
const int ndim, const double iter_tol, const double inside_tol,
- double *work, double *params, bool *inside);
+ double *work, double *params, int *inside);
/** \brief Common function that returns true if params is in [-1,1]^ndims */
- static bool inside_function(const double *params, const int ndims, const double tol);
+ static int inside_function(const double *params, const int ndims, const double tol);
};
/** \brief Given an entity handle, get an appropriate eval set, based on type & #vertices */
@@ -141,7 +141,7 @@ namespace moab {
* (in most cases, within [-1]*(dim)
*/
ErrorCode reverse_eval(const double *posn, double iter_tol, double inside_tol, double *params,
- bool *is_inside = NULL) const;
+ int *is_inside = NULL) const;
/** \brief Evaluate the jacobian of the cached entity at a given parametric location
* \param params Parameters at which to evaluate jacobian
@@ -158,7 +158,7 @@ namespace moab {
* \param params Parameters at which to query the element
* \param tol Tolerance, usually 10^-6 or so
*/
- bool inside(const double *params, const double tol) const;
+ int inside(const double *params, const double tol) const;
/** \brief Given a list of entities, return the entity the point is in, or none
* This function reverse-evaluates the entities, returning the first entity containing the point.
@@ -451,7 +451,7 @@ namespace moab {
}
inline ErrorCode ElemEvaluator::reverse_eval(const double *posn, const double iter_tol, const double inside_tol,
- double *params, bool *ins) const
+ double *params, int *ins) const
{
assert(entHandle && MBMAXTYPE != entType);
return (*evalSets[entType].reverseEvalFcn)(evalSets[entType].evalFcn, evalSets[entType].jacobianFcn, evalSets[entType].insideFcn,
@@ -493,7 +493,7 @@ namespace moab {
else return find_containing_entity(entities, point, iter_tol, inside_tol, containing_ent, params, num_evals);
}
- inline bool ElemEvaluator::inside(const double *params, const double tol) const
+ inline int ElemEvaluator::inside(const double *params, const double tol) const
{
return (*evalSets[entType].insideFcn)(params, entDim, tol);
}
diff --git a/src/LocalDiscretization/moab/LinearHex.hpp b/src/LocalDiscretization/moab/LinearHex.hpp
index 03ed33f..32d62f8 100644
--- a/src/LocalDiscretization/moab/LinearHex.hpp
+++ b/src/LocalDiscretization/moab/LinearHex.hpp
@@ -18,7 +18,7 @@ public:
static ErrorCode reverseEvalFcn(EvalFcn eval, JacobianFcn jacob, InsideFcn ins,
const double *posn, const double *verts, const int nverts, const int ndim,
const double iter_tol, const double inside_tol, double *work,
- double *params, bool *is_inside);
+ double *params, int *is_inside);
/** \brief Evaluate the jacobian at a specified parametric position */
static ErrorCode jacobianFcn(const double *params, const double *verts, const int nverts, const int ndim,
@@ -29,7 +29,7 @@ public:
double *work, double *result);
/** \brief Function that returns whether or not the parameters are inside the natural space of the element */
- static bool insideFcn(const double *params, const int ndim, const double tol);
+ static int insideFcn(const double *params, const int ndim, const double tol);
static EvalSet eval_set()
{
diff --git a/src/LocalDiscretization/moab/LinearQuad.hpp b/src/LocalDiscretization/moab/LinearQuad.hpp
index 2eb3027..897e3b4 100644
--- a/src/LocalDiscretization/moab/LinearQuad.hpp
+++ b/src/LocalDiscretization/moab/LinearQuad.hpp
@@ -18,7 +18,7 @@ public:
static ErrorCode reverseEvalFcn(EvalFcn eval, JacobianFcn jacob, InsideFcn ins,
const double *posn, const double *verts, const int nverts, const int ndim,
const double iter_tol, const double inside_tol, double *work,
- double *params, bool *is_inside);
+ double *params, int *is_inside);
/** \brief Evaluate the jacobian at a specified parametric position */
static ErrorCode jacobianFcn(const double *params, const double *verts, const int nverts, const int ndim,
@@ -29,7 +29,7 @@ public:
double *work, double *result);
/** \brief Function that returns whether or not the parameters are inside the natural space of the element */
- static bool insideFcn(const double *params, const int ndim, const double tol);
+ static int insideFcn(const double *params, const int ndim, const double tol);
static EvalSet eval_set()
{
diff --git a/src/LocalDiscretization/moab/LinearTet.hpp b/src/LocalDiscretization/moab/LinearTet.hpp
index f621b2f..713192a 100644
--- a/src/LocalDiscretization/moab/LinearTet.hpp
+++ b/src/LocalDiscretization/moab/LinearTet.hpp
@@ -18,7 +18,7 @@ public:
static ErrorCode reverseEvalFcn(EvalFcn eval, JacobianFcn jacob, InsideFcn ins,
const double *posn, const double *verts, const int nverts, const int ndim,
const double iter_tol, const double inside_tol, double *work,
- double *params, bool *is_inside);
+ double *params, int *is_inside);
/** \brief Evaluate the jacobian at a specified parametric position */
static ErrorCode jacobianFcn(const double *params, const double *verts, const int nverts, const int ndim,
@@ -32,12 +32,12 @@ public:
static ErrorCode initFcn(const double *verts, const int nverts, double *&work);
/** \brief Function that returns whether or not the parameters are inside the natural space of the element */
- static bool insideFcn(const double *params, const int ndim, const double tol);
+ static int insideFcn(const double *params, const int ndim, const double tol);
static ErrorCode evaluate_reverse(EvalFcn eval, JacobianFcn jacob, InsideFcn inside_f,
const double *posn, const double *verts, const int nverts,
const int ndim, const double iter_tol, const double inside_tol, double *work,
- double *params, bool *inside);
+ double *params, int *inside);
static EvalSet eval_set()
{
diff --git a/src/LocalDiscretization/moab/LinearTri.hpp b/src/LocalDiscretization/moab/LinearTri.hpp
index 93820a7..3ab4eca 100644
--- a/src/LocalDiscretization/moab/LinearTri.hpp
+++ b/src/LocalDiscretization/moab/LinearTri.hpp
@@ -18,7 +18,7 @@ public:
static ErrorCode reverseEvalFcn(EvalFcn eval, JacobianFcn jacob, InsideFcn ins,
const double *posn, const double *verts, const int nverts, const int ndim,
const double iter_tol, const double inside_tol, double *work,
- double *params, bool *is_inside);
+ double *params, int *is_inside);
/** \brief Evaluate the jacobian at a specified parametric position */
static ErrorCode jacobianFcn(const double *params, const double *verts, const int nverts, const int ndim,
@@ -32,12 +32,12 @@ public:
static ErrorCode initFcn(const double *verts, const int nverts, double *&work);
/** \brief Function that returns whether or not the parameters are inside the natural space of the element */
- static bool insideFcn(const double *params, const int ndim, const double tol);
+ static int insideFcn(const double *params, const int ndim, const double tol);
static ErrorCode evaluate_reverse(EvalFcn eval, JacobianFcn jacob, InsideFcn inside_f,
const double *posn, const double *verts, const int nverts,
const int ndim, const double iter_tol, const double inside_tol, double *work,
- double *params, bool *inside);
+ double *params, int *inside);
static EvalSet eval_set()
{
diff --git a/src/LocalDiscretization/moab/QuadraticHex.hpp b/src/LocalDiscretization/moab/QuadraticHex.hpp
index 0e5cfb8..ee0dcc8 100644
--- a/src/LocalDiscretization/moab/QuadraticHex.hpp
+++ b/src/LocalDiscretization/moab/QuadraticHex.hpp
@@ -18,7 +18,7 @@ public:
static ErrorCode reverseEvalFcn(EvalFcn eval, JacobianFcn jacob, InsideFcn ins,
const double *posn, const double *verts, const int nverts, const int ndim,
const double iter_tol, const double inside_tol, double *work,
- double *params, bool *is_inside);
+ double *params, int *is_inside);
/** \brief Evaluate the jacobian at a specified parametric position */
static ErrorCode jacobianFcn(const double *params, const double *verts, const int nverts, const int ndim,
@@ -29,7 +29,7 @@ public:
double *work, double *result);
/** \brief Function that returns whether or not the parameters are inside the natural space of the element */
- static bool insideFcn(const double *params, const int ndim, const double tol);
+ static int insideFcn(const double *params, const int ndim, const double tol);
static EvalSet eval_set()
{
diff --git a/src/LocalDiscretization/moab/SpectralHex.hpp b/src/LocalDiscretization/moab/SpectralHex.hpp
index 648cc70..9e60c63 100644
--- a/src/LocalDiscretization/moab/SpectralHex.hpp
+++ b/src/LocalDiscretization/moab/SpectralHex.hpp
@@ -20,7 +20,7 @@ public:
static ErrorCode reverseEvalFcn(EvalFcn eval, JacobianFcn jacob, InsideFcn ins,
const double *posn, const double *verts, const int nverts, const int ndim,
const double iter_tol, const double inside_tol, double *work,
- double *params, bool *is_inside);
+ double *params, int *is_inside);
/** \brief Evaluate the jacobian at a specified parametric position */
static ErrorCode jacobianFcn(const double *params, const double *verts, const int nverts, const int ndim,
@@ -34,7 +34,7 @@ public:
static ErrorCode initFcn(const double *verts, const int nverts, double *&work);
/** \brief Function that returns whether or not the parameters are inside the natural space of the element */
- static bool insideFcn(const double *params, const int ndim, const double tol);
+ static int insideFcn(const double *params, const int ndim, const double tol);
static EvalSet eval_set()
{
diff --git a/src/LocalDiscretization/moab/SpectralQuad.hpp b/src/LocalDiscretization/moab/SpectralQuad.hpp
index a72f9b4..4390d23 100644
--- a/src/LocalDiscretization/moab/SpectralQuad.hpp
+++ b/src/LocalDiscretization/moab/SpectralQuad.hpp
@@ -20,7 +20,7 @@ public:
static ErrorCode reverseEvalFcn(EvalFcn eval, JacobianFcn jacob, InsideFcn ins,
const double *posn, const double *verts, const int nverts, const int ndim,
const double iter_tol, const double inside_tol, double *work,
- double *params, bool *is_inside);
+ double *params, int *is_inside);
/** \brief Evaluate the jacobian at a specified parametric position */
static ErrorCode jacobianFcn(const double *params, const double *verts, const int nverts, const int ndim,
@@ -33,7 +33,7 @@ public:
static ErrorCode initFcn(const double *verts, const int nverts, double *&work);
/** \brief Function that returns whether or not the parameters are inside the natural space of the element */
- static bool insideFcn(const double *params, const int ndim, const double tol);
+ static int insideFcn(const double *params, const int ndim, const double tol);
static EvalSet eval_set()
{
diff --git a/src/SpatialLocator.cpp b/src/SpatialLocator.cpp
index f765ab9..593c5e0 100644
--- a/src/SpatialLocator.cpp
+++ b/src/SpatialLocator.cpp
@@ -3,6 +3,13 @@
#include "moab/ElemEvaluator.hpp"
#include "moab/AdaptiveKDTree.hpp"
+// include ScdInterface for box partitioning
+#include "moab/ScdInterface.hpp"
+
+#ifdef USE_MPI
+#include "moab/ParallelComm.hpp"
+#endif
+
bool debug = true;
namespace moab
@@ -34,21 +41,256 @@ namespace moab
}
#ifdef USE_MPI
- ErrorCode SpatialLocator::par_locate_points(Range &/*vertices*/,
- const double /*rel_iter_tol*/, const double /*abs_iter_tol*/,
- const double /*inside_tol*/)
+ ErrorCode SpatialLocator::initialize_global_box(ParallelComm *pc)
{
- return MB_UNSUPPORTED_OPERATION;
+ if (!pc) return MB_FAILURE;
+
+ BoundBox box, gbox;
+ ErrorCode rval = myTree->get_bounding_box(box);
+
+ //step 2
+ // get the global bounding box
+ double sendbuffer[6];
+ double rcvbuffer[6];
+
+ box.get(sendbuffer); //fill sendbuffer with local box, max values in [0:2] min values in [3:5]
+ sendbuffer[3] *= -1;
+ sendbuffer[4] *= -1; //negate Xmin,Ymin,Zmin to get their minimum using MPI_MAX
+ sendbuffer[5] *= -1; //to avoid calling MPI_Allreduce again with MPI_MIN
+
+ int mpi_err = MPI_Allreduce(sendbuffer, rcvbuffer, 6, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
+ if (MPI_SUCCESS != mpi_err) return MB_FAILURE;
+
+ rcvbuffer[3] *= -1;
+ rcvbuffer[4] *= -1; //negate Xmin,Ymin,Zmin again to get original values
+ rcvbuffer[5] *= -1;
+
+ globalBox.update_max(&rcvbuffer[0]); //saving values in globalBox
+ globalBox.update_min(&rcvbuffer[3]);
+
+ // compute the alternate decomposition; use ScdInterface::compute_partition_sqijk for this
+ ScdParData spd;
+ spd.partMethod = ScdParData::SQIJK;
+ spd.gPeriodic[0] = spd.gPeriodic[1] = spd.gPeriodic[2] = 0;
+ double lg = log10((box.bMax - box.bMin).length());
+ double mfactor = pow(10.0, 6 - lg);
+ int ldims[3], lper[3];
+ double dgijk[6];
+ box.get(dgijk);
+ for (int i = 0; i < 6; i++) spd.gDims[i] = dgijk[i] * mfactor;
+ rval = ScdInterface::compute_partition(pc->size(), pc->rank(), spd,
+ ldims, lper, regNums);
+ if (MB_SUCCESS != rval) return rval;
+ // all we're really interested in is regNums[i], #procs in each direction
+
+ for (int i = 0; i < 3; i++)
+ regDeltaXYZ[i] = (globalBox.bMax[i] - globalBox.bMin[i])/double(regNums[i]); //size of each region
+
+ return MB_SUCCESS;
+ }
+
+//this function sets up the TupleList TLreg_o containing the registration messages
+// and sends it
+ ErrorCode SpatialLocator::register_with_forwarders(ParallelComm *pc, TupleList &TLreg_o)
+ {
+
+ //step 4
+ //set up TLreg_o
+ TLreg_o.initialize(1,0,0,6,0);
+ // TLreg_o (int destProc, real Xmin, Ymin, Zmin, Xmax, Ymax, Zmax)
+
+ int dest;
+ double boxtosend[6];
+
+ BoundBox box;
+ ErrorCode result = myTree->get_bounding_box(box);
+ if (result != MB_SUCCESS)
+ return result;
+
+ box.get(boxtosend);
+
+ //iterate over all regions overlapping with my bounding box using the computerd corner IDs
+ for (int k = cornerIJK[2]; k <= cornerIJK[5]; k++) {
+ for (int j = cornerIJK[1]; j <= cornerIJK[4]; j++) {
+ for (int i = cornerIJK[0]; i <= cornerIJK[3]; i++) {
+ dest = k * regNums[0]*regNums[1] + j * regNums[0] + i;
+ TLreg_o.push_back(&dest, NULL, NULL, boxtosend);
+ }
+ }
+ }
+
+ //step 5
+ //send TLreg_o, receive TLrequests_i
+ if (pc) pc->proc_config().crystal_router()->gs_transfer(1, TLreg_o, 0);
+
+ //step 6
+ //Read registration requests from TLreg_o and add to list of procs to forward to
+ //get number of tuples sent to me
+
+ //read tuples and fill processor list;
+ int NN = TLreg_o.get_n();
+ for (int i=0; i < NN; i++)
+ //TLreg_o is now TLrequests_i
+ srcProcBoxes[TLreg_o.vi_rd[i]] = BoundBox(TLreg_o.vr_rd+6*i);
+
+ return MB_SUCCESS;
}
- ErrorCode SpatialLocator::par_locate_points(const double */*pos*/, int /*num_points*/,
+ ErrorCode SpatialLocator::par_locate_points(ParallelComm */*pc*/,
+ Range &/*vertices*/,
const double /*rel_iter_tol*/, const double /*abs_iter_tol*/,
const double /*inside_tol*/)
{
return MB_UNSUPPORTED_OPERATION;
}
-#endif
+
+ bool is_initialized(int i) {return (i == -1);}
+
+ ErrorCode SpatialLocator::par_locate_points(ParallelComm *pc,
+ const double *pos, int num_points,
+ const double rel_iter_tol, const double abs_iter_tol,
+ const double inside_tol)
+ {
+ ErrorCode rval;
+ //TUpleList used for communication
+ TupleList TLreg_o; //TLregister_outbound
+ TupleList TLquery_o; //TLquery_outbound
+ TupleList TLforward_o; //TLforward_outbound
+ TupleList TLsearch_results_o; //TLsearch_results_outbound
+
+
+ // steps 1-2 - initialize the alternative decomposition box from global box
+ rval = initialize_global_box(pc);
+
+ //step 3 - compute the IDs of the regions which contain each corner of local box
+ rval = compute_corner_ijks();
+ if (rval != MB_SUCCESS) return rval;
+
+ //steps 4-6 - set up TLreg_o, gs_transfer, gather registrations
+ rval = register_with_forwarders(pc, TLreg_o);
+ if (rval != MB_SUCCESS) return rval;
+
+ // actual parallel point location using intermediate partition
+
+ // target_pts: TL(to_proc, tgt_index, x, y, z): tuples sent to source mesh procs representing pts to be located
+ // source_pts: TL(from_proc, tgt_index, src_index): results of source mesh proc point location, ready to send
+ // back to tgt procs; src_index of -1 indicates point not located (arguably not useful...)
+
+ unsigned int my_rank = (pc? pc->proc_config().proc_rank() : 0);
+
+ //TLquery_o: Tuples sent to forwarder proc
+ //TL (toProc, OriginalSourceProc, targetIndex, X,Y,Z)
+
+ //TLforw_req_i: Tuples to forward to corresponding procs (forwarding requests)
+ //TL (sourceProc, OriginalSourceProc, targetIndex, X,Y,Z)
+
+ TLquery_o.initialize(3,0,0,3,0);
+
+ int iargs[3];
+
+ for (int pnt=0; pnt < 3*num_points; pnt+=3)
+ {
+ int forw_id = proc_from_point(pos+pnt); //get ID of proc resonsible of the region the proc is in
+
+ iargs[0] = forw_id; //toProc
+ iargs[1] = my_rank; //originalSourceProc
+ iargs[2] = pnt/3; //targetIndex
+
+ TLquery_o.push_back(iargs, NULL, NULL, const_cast<double*>(pos+pnt));
+ }
+
+ //send point search queries to forwarders
+ if (pc)
+ pc->proc_config().crystal_router()->gs_transfer(1, TLquery_o, 0);
+
+ //now read forwarding requests and forward to corresponding procs
+ //TLquery_o is now TLforw_req_i
+
+ //TLforward_o: query messages forwarded to corresponding procs
+ //TL (toProc, OriginalSourceProc, targetIndex, X,Y,Z)
+
+ TLforward_o.initialize(3,0,0,3,0);
+
+ int NN = TLquery_o.get_n();
+
+ for (int i=0; i < NN; i++) {
+ iargs[1] = TLquery_o.vi_rd[3*i+1]; //get OriginalSourceProc
+ iargs[2] = TLquery_o.vi_rd[3*i+2]; //targetIndex
+ CartVect tmp_pnt(TLquery_o.vr_rd+3*i);
+
+ //compare coordinates to list of bounding boxes
+ for (std::map<int, BoundBox>::iterator mit = srcProcBoxes.begin(); mit != srcProcBoxes.end(); mit++) {
+ if ((*mit).second.contains_point(tmp_pnt.array(), abs_iter_tol)) {
+ iargs[0] = (*mit).first;
+ TLforward_o.push_back(iargs, NULL, NULL, tmp_pnt.array());
+ }
+ }
+
+ }
+
+ if (pc)
+ pc->proc_config().crystal_router()->gs_transfer(1, TLforward_o, 0);
+
+ //step 12
+ //now read Point Search requests
+ //TLforward_o is now TLsearch_req_i
+ //TLsearch_req_i: (sourceProc, OriginalSourceProc, targetIndex, X,Y,Z)
+
+ NN = TLforward_o.get_n();
+
+ //TLsearch_results_o
+ //TL: (OriginalSourceProc, targetIndex, sourceIndex, U,V,W);
+ TLsearch_results_o.initialize(3,0,0,3,0);
+
+ //step 13 is done in test_local_box
+
+ std::vector<double> params(3*NN);
+ std::vector<int> is_inside(NN, 0);
+ std::vector<EntityHandle> ents(NN, 0);
+
+ rval = locate_points(TLforward_o.vr_rd, TLforward_o.get_n(),
+ &ents[0], ¶ms[0], &is_inside[0],
+ rel_iter_tol, abs_iter_tol, inside_tol);
+ if (MB_SUCCESS != rval)
+ return rval;
+ for (int i = 0; i < NN; i++) {
+ if (is_inside[i]) {
+ iargs[0] = TLforward_o.vi_rd[3*i+1];
+ iargs[1] = TLforward_o.vi_rd[3*i+2];
+ iargs[2] = locTable.get_n();
+ TLsearch_results_o.push_back(iargs, NULL, NULL, NULL);
+ locTable.push_back(const_cast<int*>(&TLforward_o.vi_rd[3*i+1]), NULL, &ents[i], ¶ms[3*i]);
+ }
+ }
+
+ //step 14: send TLsearch_results_o and receive TLloc_i
+ if (pc)
+ pc->proc_config().crystal_router()->gs_transfer(1, TLsearch_results_o, 0);
+
+
+ // store proc/index tuples in parLocTable
+ parLocTable.initialize(2, 0, 0, 0, num_points);
+ parLocTable.enableWriteAccess();
+ std::fill(parLocTable.vi_wr, parLocTable.vi_wr + 2*num_points, -1);
+
+ for (unsigned int i = 0; i < TLsearch_results_o.get_n(); i++) {
+ int idx = TLsearch_results_o.vi_rd[3*i+1];
+ parLocTable.vi_wr[2*idx] = TLsearch_results_o.vi_rd[3*i];
+ parLocTable.vi_wr[2*idx+1] = TLsearch_results_o.vi_rd[3*i+2];
+ }
+
+ if (debug) {
+ int num_found = num_points - 0.5 *
+ std::count_if(parLocTable.vi_wr, parLocTable.vi_wr + 2*num_points, is_initialized);
+ std::cout << "Points found = " << num_found << " out of " << num_points << "." << std::endl;
+ }
+
+ return MB_SUCCESS;
+ }
+
+#endif
+
ErrorCode SpatialLocator::locate_points(Range &verts,
const double rel_iter_tol, const double abs_iter_tol,
const double inside_tol)
@@ -82,7 +324,7 @@ namespace moab
}
ErrorCode SpatialLocator::locate_points(Range &verts,
- EntityHandle *ents, double *params, bool *is_inside,
+ EntityHandle *ents, double *params, int *is_inside,
const double rel_iter_tol, const double abs_iter_tol,
const double inside_tol)
{
@@ -94,7 +336,7 @@ namespace moab
}
ErrorCode SpatialLocator::locate_points(const double *pos, int num_points,
- EntityHandle *ents, double *params, bool *is_inside,
+ EntityHandle *ents, double *params, int *is_inside,
const double rel_iter_tol, const double abs_iter_tol,
const double inside_tol)
{
@@ -154,8 +396,8 @@ namespace moab
if(rval != MB_SUCCESS) return rval;
// loop over the range_leaf
- bool tmp_inside;
- bool *is_ptr = (is_inside ? is_inside+i : &tmp_inside);
+ int tmp_inside;
+ int *is_ptr = (is_inside ? is_inside+i : &tmp_inside);
*is_ptr = false;
EntityHandle ent = 0;
for(Range::iterator rit = range_leaf.begin(); rit != range_leaf.end(); rit++)
diff --git a/src/moab/BoundBox.hpp b/src/moab/BoundBox.hpp
index 5c51341..3e525ac 100644
--- a/src/moab/BoundBox.hpp
+++ b/src/moab/BoundBox.hpp
@@ -27,6 +27,7 @@ namespace moab {
void update_min(const double *coords);
void update_max(const BoundBox &other_box);
void update_max(const double *coords);
+ ErrorCode get(double *coords);
/** \brief Return the diagonal length of this box
*/
@@ -123,6 +124,13 @@ namespace moab {
bMax[2] = std::max(bMax[2], coords[2]);
}
+ inline ErrorCode BoundBox::get(double *coords)
+ {
+ bMin.get(coords);
+ bMax.get(coords+3);
+ return MB_SUCCESS;
+ }
+
inline void BoundBox::compute_center(CartVect ¢er){
center = 0.5 * (bMin + bMax);
}
diff --git a/src/moab/SpatialLocator.hpp b/src/moab/SpatialLocator.hpp
index bb639d0..00a7857 100644
--- a/src/moab/SpatialLocator.hpp
+++ b/src/moab/SpatialLocator.hpp
@@ -34,15 +34,18 @@
#include "moab/Tree.hpp"
#include "moab/Range.hpp"
#include "moab/TupleList.hpp"
+#include "moab/BoundBox.hpp"
#include <string>
#include <vector>
+#include <map>
#include <math.h>
namespace moab {
class Interface;
class ElemEvaluator;
+ class ParallelComm;
class SpatialLocator
{
@@ -61,13 +64,13 @@ namespace moab {
/* locate a set of vertices, Range variant */
ErrorCode locate_points(Range &vertices,
- EntityHandle *ents, double *params, bool *is_inside = NULL,
+ EntityHandle *ents, double *params, int *is_inside = NULL,
const double rel_iter_tol = 1.0e-10, const double abs_iter_tol = 1.0e-10,
const double inside_tol = 1.0e-6);
/* locate a set of points */
ErrorCode locate_points(const double *pos, int num_points,
- EntityHandle *ents, double *params, bool *is_inside = NULL,
+ EntityHandle *ents, double *params, int *is_inside = NULL,
const double rel_iter_tol = 1.0e-10, const double abs_iter_tol = 1.0e-10,
const double inside_tol = 1.0e-6);
@@ -107,20 +110,26 @@ namespace moab {
* variables (TupleList) locTable and parLocTable (see comments on locTable and parLocTable for
* structure of those tuples).
*/
- ErrorCode par_locate_points(Range &vertices,
- const double rel_iter_tol = 1.0e-10, const double abs_iter_tol = 1.0e-10,
- const double inside_tol = 1.0e-6);
+ ErrorCode par_locate_points(ParallelComm *pc,
+ Range &vertices,
+ const double rel_iter_tol = 1.0e-10, const double abs_iter_tol = 1.0e-10,
+ const double inside_tol = 1.0e-6);
/* locate a set of points, storing results on TupleList in this class
* Locate a set of points, storing the detailed results in member
* variables (TupleList) locTable and parLocTable (see comments on locTable and parLocTable for
* structure of those tuples).
*/
- ErrorCode par_locate_points(const double *pos, int num_points,
- const double rel_iter_tol = 1.0e-10, const double abs_iter_tol = 1.0e-10,
- const double inside_tol = 1.0e-6);
+ ErrorCode par_locate_points(ParallelComm *pc,
+ const double *pos, int num_points,
+ const double rel_iter_tol = 1.0e-10, const double abs_iter_tol = 1.0e-10,
+ const double inside_tol = 1.0e-6);
#endif
+ /** \brief Return the MOAB interface associated with this locator
+ */
+ Interface* moab() { return mbImpl; }
+
/* return the tree */
Tree *get_tree() {return myTree;}
@@ -153,10 +162,38 @@ namespace moab {
/* locate a point */
ErrorCode locate_point(const double *pos,
- EntityHandle &ent, double *params, bool *is_inside = NULL,
+ EntityHandle &ent, double *params, int *is_inside = NULL,
const double rel_iter_tol = 1.0e-10, const double abs_iter_tol = 1.0e-10,
const double inside_tol = 1.0e-6);
+#ifdef USE_MPI
+ /* MPI_ReduceAll source mesh bounding boxes to get global source mesh bounding box
+ */
+ ErrorCode initialize_global_box(ParallelComm *pc);
+
+ /* compute the ijk location of this proc's target mesh bounding box based on intermediate decomp
+ * extents/divisions
+ */
+ inline ErrorCode compute_corner_ijks();
+
+ /* for a given point in space, compute its ijk location in the intermediate decomposition
+ */
+ inline ErrorCode get_point_ijk(const CartVect &point, int *ijk) const;
+
+ /* given an ijk location in the intermediate partition, return the proc rank for that location
+ */
+ inline int proc_from_ijk(const int *ijk) const;
+
+ /* given a point in space, return the proc responsible for that point from the intermediate decomp
+ */
+ inline int proc_from_point(const double *pos) const;
+
+ /* register my source mesh with intermediate decomposition procs
+ */
+ ErrorCode register_with_forwarders(ParallelComm *pc, TupleList &TLreg_o);
+
+#endif
+
/* MOAB instance */
Interface* mbImpl;
@@ -197,6 +234,28 @@ namespace moab {
* The indexing of parLocTable corresponds to that of the points/entities passed in.
*/
TupleList parLocTable;
+
+ /* \brief Global bounding box, used in parallel spatial location
+ */
+ BoundBox globalBox;
+
+ /* \brief Regional delta xyz, used in parallel spatial location
+ */
+ CartVect regDeltaXYZ;
+
+ /* \brief Number of regions in each of 3 directions
+ */
+ int regNums[3];
+
+ /* \brief Region numbers for my box's min/max corners, used in parallel spatial location
+ */
+ int cornerIJK[6];
+
+ /* \brief Map from source processor to bounding box of that proc's source mesh
+ *
+ */
+ std::map<int, BoundBox> srcProcBoxes;
+
};
inline SpatialLocator::~SpatialLocator()
@@ -205,7 +264,7 @@ namespace moab {
}
inline ErrorCode SpatialLocator::locate_point(const double *pos,
- EntityHandle &ent, double *params, bool *is_inside,
+ EntityHandle &ent, double *params, int *is_inside,
const double rel_iter_tol, const double abs_iter_tol,
const double inside_tol)
{
@@ -217,6 +276,40 @@ namespace moab {
return myTree->get_bounding_box(box);
}
+ inline ErrorCode SpatialLocator::compute_corner_ijks() {
+
+ BoundBox box;
+ ErrorCode rval = myTree->get_bounding_box(box);
+ if (rval != MB_SUCCESS) return rval;
+
+ rval = get_point_ijk(box.bMin, cornerIJK);
+ if (MB_SUCCESS != rval) return rval;
+ rval = get_point_ijk(box.bMax, cornerIJK+3);
+ return rval;
+ }
+
+ inline ErrorCode SpatialLocator::get_point_ijk(const CartVect &point, int *ijk) const
+ {
+ for (int i = 0; i < 3; i++)
+ ijk[i] = (regNums[i] > 1 ? (point[i] - globalBox.bMin[i]) / regDeltaXYZ[i] : 0);
+
+ return MB_SUCCESS;
+ }
+
+ inline int SpatialLocator::proc_from_ijk(const int *ijk) const
+ {
+ return ijk[2] * regNums[0]*regNums[1] + ijk[1] * regNums[0] + ijk[0];
+ }
+
+ inline int SpatialLocator::proc_from_point(const double *pos) const
+ {
+ int ijk[3];
+ ErrorCode rval = get_point_ijk(CartVect(pos), ijk);
+ if (MB_SUCCESS != rval) return -1;
+
+ return ijk[2] * regNums[0]*regNums[1] + ijk[1] * regNums[0] + ijk[0];
+ }
+
} // namespace moab
#endif
diff --git a/test/elem_eval_test.cpp b/test/elem_eval_test.cpp
index beb8704..212931f 100644
--- a/test/elem_eval_test.cpp
+++ b/test/elem_eval_test.cpp
@@ -50,7 +50,7 @@ void test_eval(ElemEvaluator &ee, bool test_integrate)
{
CartVect params, posn, params2;
- bool is_inside;
+ int is_inside;
Matrix3 jacob;
ErrorCode rval;
int ent_dim = ee.get_moab()->dimension_from_handle(ee.get_ent_handle());
diff --git a/test/oldinc/test_oldinc.cpp b/test/oldinc/test_oldinc.cpp
index 59fa57b..fc4b026 100644
--- a/test/oldinc/test_oldinc.cpp
+++ b/test/oldinc/test_oldinc.cpp
@@ -8,12 +8,12 @@ typedef MBProcConfig Foo3;
#include "MBEntityHandle.h"
MBEntityHandle handle = 0;
#include "MBTypes.h"
-MBEntityType type = MBVERTEX;
-MBTag tag = 0;
+MBEntityType tp = MBVERTEX;
+MBTag tg = 0;
MBErrorCode code = MB_SUCCESS;
-int i = MB_VARIABLE_LENGTH;
+int it = MB_VARIABLE_LENGTH;
MBTagType tagtype = MB_TAG_DENSE;
-MBDataType datatype = MB_TYPE_HANDLE;
+MBDataType dttype = MB_TYPE_HANDLE;
MBEntitySetProperty prop = MESHSET_SET;
// test version info
diff --git a/test/parallel/Makefile.am b/test/parallel/Makefile.am
index 6dd3207..ca35879 100644
--- a/test/parallel/Makefile.am
+++ b/test/parallel/Makefile.am
@@ -35,6 +35,7 @@ TESTS = pcomm_unit \
uber_parallel_test \
scdtest \
pcomm_serial \
+ par_spatial_locator_test \
$(NETCDF_TESTS) \
$(HDF5_TESTS) \
$(COUPLER_TESTS) \
@@ -89,6 +90,7 @@ scdpart_SOURCES = scdpart.cpp
read_nc_par_SOURCES = ../io/read_nc.cpp
ucdtrvpart_SOURCES = ucdtrvpart.cpp
mpastrvpart_SOURCES = mpastrvpart.cpp
+par_spatial_locator_test_SOURCES = par_spatial_locator_test.cpp
if ENABLE_mbcoupler
par_coupler_test_SOURCES = par_coupler_test.cpp
diff --git a/test/parallel/par_spatial_locator_test.cpp b/test/parallel/par_spatial_locator_test.cpp
new file mode 100644
index 0000000..40cb98d
--- /dev/null
+++ b/test/parallel/par_spatial_locator_test.cpp
@@ -0,0 +1,187 @@
+#include "moab/Core.hpp"
+#include "moab/SpatialLocator.hpp"
+#include "moab/Tree.hpp"
+#include "moab/HomXform.hpp"
+#include "moab/ScdInterface.hpp"
+#include "moab/CartVect.hpp"
+#include "moab/BVHTree.hpp"
+#include "moab/ProgOptions.hpp"
+#include "moab/CpuTimer.hpp"
+#include "moab/ParallelComm.hpp"
+#include "moab_mpi.h"
+
+#include "TestUtil.hpp"
+
+#include <cstdlib>
+#include <sstream>
+#include <string>
+
+using namespace moab;
+
+void test_kd_tree();
+void test_bvh_tree();
+void test_locator(SpatialLocator *sl);
+
+ErrorCode create_hex_mesh(Interface &mb, Range &elems, int n, int dim);
+ErrorCode load_file(Interface &mb, std::string &fn, Range &elems);
+
+int max_depth = 30;
+int npoints = 1000;
+int leaf = 6;
+bool print_tree = false;
+int ints = 10;
+std::string fname;
+
+int main(int argc, char **argv)
+{
+ int fail = MPI_Init(&argc, &argv);
+ if (fail) return fail;
+
+ ProgOptions po("spatial_locator_test options" );
+ po.addOpt<std::string>( "file,f", "Input filename", &fname);
+ po.addOpt<int>( "ints,i", "Number of intervals on each side of scd mesh", &ints);
+ po.addOpt<int>( "leaf,l", "Maximum number of elements per leaf", &leaf);
+ po.addOpt<int>( "max_depth,m", "Maximum depth of tree", &max_depth);
+ po.addOpt<int>( "npoints,n", "Number of query points", &npoints);
+ po.addOpt<void>( "print,p", "Print tree details", &print_tree);
+ po.parseCommandLine(argc, argv);
+
+ RUN_TEST(test_kd_tree);
+ RUN_TEST(test_bvh_tree);
+
+ fail = MPI_Finalize();
+ if (fail) return fail;
+
+ return 0;
+}
+
+void test_kd_tree()
+{
+ ErrorCode rval;
+ Core mb;
+
+ // create a simple mesh to test
+ Range elems;
+ if (fname.empty()) {
+ rval = create_hex_mesh(mb, elems, ints, 3); CHECK_ERR(rval);
+ }
+ else {
+ rval = load_file(mb, fname, elems); CHECK_ERR(rval);
+ }
+
+ // initialize spatial locator with the elements and the default tree type
+ SpatialLocator *sl = new SpatialLocator(&mb, elems);
+
+ test_locator(sl);
+
+ // destroy spatial locator, and tree along with it
+ delete sl;
+}
+
+void test_bvh_tree()
+{
+ ErrorCode rval;
+ Core mb;
+
+ // create a simple mesh to test
+ Range elems;
+ if (fname.empty()) {
+ rval = create_hex_mesh(mb, elems, ints, 3); CHECK_ERR(rval);
+ }
+ else {
+ rval = load_file(mb, fname, elems); CHECK_ERR(rval);
+ }
+
+ // initialize spatial locator with the elements and a BVH tree
+ BVHTree bvh(&mb);
+ std::ostringstream opts;
+ opts << "MAX_DEPTH=" << max_depth << ";MAX_PER_LEAF=" << leaf;
+ FileOptions fo(opts.str().c_str());
+ rval = bvh.parse_options(fo);
+ SpatialLocator *sl = new SpatialLocator(&mb, elems, &bvh);
+ test_locator(sl);
+
+ // destroy spatial locator, and tree along with it
+ delete sl;
+}
+
+bool is_neg(int is_neg)
+{
+ return (is_neg < 0);
+}
+
+void test_locator(SpatialLocator *sl)
+{
+ CartVect box_del;
+ BoundBox box;
+ ErrorCode rval = sl->get_bounding_box(box); CHECK_ERR(rval);
+ box_del = box.bMax - box.bMin;
+
+ double denom = 1.0 / (double)RAND_MAX;
+ std::vector<CartVect> test_pts(npoints);
+
+ for (int i = 0; i < npoints; i++) {
+ // generate a small number of random point to test
+ double rx = (double)rand() * denom, ry = (double)rand() * denom, rz = (double)rand() * denom;
+ test_pts[i] = box.bMin + CartVect(rx*box_del[0], ry*box_del[1], rz*box_del[2]);
+ }
+
+ // call spatial locator to locate points
+ ParallelComm *pc = ParallelComm::get_pcomm(sl->moab(), 0);
+ CHECK(pc != NULL);
+
+ rval = sl->par_locate_points(pc, test_pts[0].array(), npoints); CHECK_ERR(rval);
+ if (pc->rank() == 0) {
+ int num_out = std::count_if(sl->par_loc_table().vi_rd, sl->par_loc_table().vi_rd+2*npoints, is_neg);
+ num_out /= 2;
+
+ std::cout << "Number of points inside an element = " << npoints-num_out << " out of " << npoints << std::endl;
+ std::cout << "Traversal stats:" << std::endl;
+ sl->get_tree()->tree_stats().output();
+
+ if (print_tree) {
+ std::cout << "Tree information: " << std::endl;
+ rval = sl->get_tree()->print();
+ CHECK_ERR(rval);
+ }
+ }
+}
+
+ErrorCode create_hex_mesh(Interface &mb, Range &elems, int n, int dim)
+{
+ ScdInterface *scdi;
+ ErrorCode rval = mb.query_interface(scdi); CHECK_ERR(rval);
+ ScdParData spd;
+ spd.gDims[0] = n-1;
+ spd.partMethod = ScdParData::SQIJK;
+ if (dim > 1) spd.gDims[1] = n-1;
+ if (dim > 2) spd.gDims[2] = n-1;
+ ScdBox *new_box;
+ rval = scdi->construct_box(HomCoord(0, 0, 0), HomCoord(0, 0, 0),
+ NULL, 0, new_box, NULL, &spd, false, 0); CHECK_ERR(rval);
+
+ rval = mb.get_entities_by_dimension(0, dim, elems); CHECK_ERR(rval);
+
+ return rval;
+}
+
+ErrorCode load_file(Interface &mb, std::string &fn, Range &elems)
+{
+ std::string options;
+ ErrorCode rval;
+
+ options = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;";
+ rval = mb.load_file(fn.c_str(), NULL, options.c_str());
+ if (MB_SUCCESS != rval) return rval;
+
+ rval = mb.get_entities_by_dimension(0, 3, elems);
+ if (MB_SUCCESS != rval) return rval;
+
+ if (elems.empty()) {
+ rval = mb.get_entities_by_dimension(0, 2, elems);
+ if (MB_SUCCESS != rval) return rval;
+ }
+
+ return MB_SUCCESS;
+}
+
diff --git a/test/perf/point_location/elem_eval_time.cpp b/test/perf/point_location/elem_eval_time.cpp
index ef86762..2f573b2 100644
--- a/test/perf/point_location/elem_eval_time.cpp
+++ b/test/perf/point_location/elem_eval_time.cpp
@@ -155,7 +155,7 @@ ErrorCode time_reverse_eval(Interface *mbi, int method, Range &elems,
ElemEvaluator eeval(mbi);
eeval.set_eval_set(*elems.begin());
eeval.set_tag_handle(0, 0); // indicates coordinates as the field to evaluate
- bool ins;
+ int ins;
for (rit = elems.begin(), i = 0; rit != elems.end(); rit++, i++) {
eeval.set_ent_handle(*rit);
rval = eeval.reverse_eval(coords[i].array(), 1.0e-10, 1.0e-6, params[i].array(), &ins);
diff --git a/test/perf/point_location/sploc_searching_perf.cpp b/test/perf/point_location/sploc_searching_perf.cpp
index 158504d..ad17241 100644
--- a/test/perf/point_location/sploc_searching_perf.cpp
+++ b/test/perf/point_location/sploc_searching_perf.cpp
@@ -148,7 +148,7 @@ ErrorCode test_locator(SpatialLocator &sl, int npoints, double rtol, double &cpu
std::vector<CartVect> test_pts(npoints), test_res(npoints);
std::vector<EntityHandle> ents(npoints);
- bool *is_in = new bool[npoints];
+ int *is_in = new int[npoints];
double denom = 1.0 / (double)RAND_MAX;
for (int i = 0; i < npoints; i++) {
diff --git a/test/perf/point_location/tree_searching_perf.cpp b/test/perf/point_location/tree_searching_perf.cpp
index 9e35d9b..230e911 100644
--- a/test/perf/point_location/tree_searching_perf.cpp
+++ b/test/perf/point_location/tree_searching_perf.cpp
@@ -139,7 +139,7 @@ ErrorCode test_locator(SpatialLocator &sl, int npoints, double &cpu_time, double
std::vector<CartVect> test_pts(npoints), test_res(npoints);
std::vector<EntityHandle> ents(npoints);
- bool *is_in = new bool[npoints];
+ int *is_in = new int[npoints];
double denom = 1.0 / (double)RAND_MAX;
for (int i = 0; i < npoints; i++) {
diff --git a/test/spatial_locator_test.cpp b/test/spatial_locator_test.cpp
index 6ffa109..2ede04f 100644
--- a/test/spatial_locator_test.cpp
+++ b/test/spatial_locator_test.cpp
@@ -108,7 +108,7 @@ void test_locator(SpatialLocator *sl)
box_del = box.bMax - box.bMin;
double denom = 1.0 / (double)RAND_MAX;
- bool is_in;
+ int is_in;
EntityHandle ent;
for (int i = 0; i < npoints; i++) {
// generate a small number of random point to test
https://bitbucket.org/fathomteam/moab/commits/7ad97741f8d3/
Changeset: 7ad97741f8d3
Branch: tautges/par_spatial_locator
User: tautges
Date: 2014-01-14 02:09:00
Summary: Some debugging for parallel spatial location based on intermediate decomposition. Works in parallel now.
Affected #: 3 files
diff --git a/src/SpatialLocator.cpp b/src/SpatialLocator.cpp
index 593c5e0..ae124dc 100644
--- a/src/SpatialLocator.cpp
+++ b/src/SpatialLocator.cpp
@@ -26,6 +26,9 @@ namespace moab
myDim = mbImpl->dimension_from_handle(*elems.rbegin());
ErrorCode rval = myTree->build_tree(myElems);
if (MB_SUCCESS != rval) throw rval;
+
+ rval = myTree->get_bounding_box(localBox);
+ if (MB_SUCCESS != rval) throw rval;
}
}
@@ -41,45 +44,44 @@ namespace moab
}
#ifdef USE_MPI
- ErrorCode SpatialLocator::initialize_global_box(ParallelComm *pc)
+ ErrorCode SpatialLocator::initialize_intermediate_partition(ParallelComm *pc)
{
if (!pc) return MB_FAILURE;
- BoundBox box, gbox;
- ErrorCode rval = myTree->get_bounding_box(box);
+ BoundBox gbox;
//step 2
// get the global bounding box
double sendbuffer[6];
double rcvbuffer[6];
- box.get(sendbuffer); //fill sendbuffer with local box, max values in [0:2] min values in [3:5]
- sendbuffer[3] *= -1;
- sendbuffer[4] *= -1; //negate Xmin,Ymin,Zmin to get their minimum using MPI_MAX
- sendbuffer[5] *= -1; //to avoid calling MPI_Allreduce again with MPI_MIN
+ localBox.get(sendbuffer); //fill sendbuffer with local box, max values in [0:2] min values in [3:5]
+ sendbuffer[0] *= -1;
+ sendbuffer[1] *= -1; //negate Xmin,Ymin,Zmin to get their minimum using MPI_MAX
+ sendbuffer[2] *= -1; //to avoid calling MPI_Allreduce again with MPI_MIN
int mpi_err = MPI_Allreduce(sendbuffer, rcvbuffer, 6, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
if (MPI_SUCCESS != mpi_err) return MB_FAILURE;
- rcvbuffer[3] *= -1;
- rcvbuffer[4] *= -1; //negate Xmin,Ymin,Zmin again to get original values
- rcvbuffer[5] *= -1;
+ rcvbuffer[0] *= -1;
+ rcvbuffer[1] *= -1; //negate Xmin,Ymin,Zmin again to get original values
+ rcvbuffer[2] *= -1;
- globalBox.update_max(&rcvbuffer[0]); //saving values in globalBox
- globalBox.update_min(&rcvbuffer[3]);
+ globalBox.update_max(&rcvbuffer[3]); //saving values in globalBox
+ globalBox.update_min(&rcvbuffer[0]);
// compute the alternate decomposition; use ScdInterface::compute_partition_sqijk for this
ScdParData spd;
spd.partMethod = ScdParData::SQIJK;
spd.gPeriodic[0] = spd.gPeriodic[1] = spd.gPeriodic[2] = 0;
- double lg = log10((box.bMax - box.bMin).length());
+ double lg = log10((localBox.bMax - localBox.bMin).length());
double mfactor = pow(10.0, 6 - lg);
int ldims[3], lper[3];
double dgijk[6];
- box.get(dgijk);
+ localBox.get(dgijk);
for (int i = 0; i < 6; i++) spd.gDims[i] = dgijk[i] * mfactor;
- rval = ScdInterface::compute_partition(pc->size(), pc->rank(), spd,
- ldims, lper, regNums);
+ ErrorCode rval = ScdInterface::compute_partition(pc->size(), pc->rank(), spd,
+ ldims, lper, regNums);
if (MB_SUCCESS != rval) return rval;
// all we're really interested in is regNums[i], #procs in each direction
@@ -91,8 +93,16 @@ namespace moab
//this function sets up the TupleList TLreg_o containing the registration messages
// and sends it
- ErrorCode SpatialLocator::register_with_forwarders(ParallelComm *pc, TupleList &TLreg_o)
+ ErrorCode SpatialLocator::register_src_with_intermediate_procs(ParallelComm *pc, double abs_iter_tol, TupleList &TLreg_o)
{
+ int corner_ijk[6];
+
+ // step 3: compute ijks of local box corners in intermediate partition
+ // get corner ijk values for my box
+ ErrorCode rval = get_point_ijk(localBox.bMin-CartVect(abs_iter_tol), abs_iter_tol, corner_ijk);
+ if (MB_SUCCESS != rval) return rval;
+ rval = get_point_ijk(localBox.bMax+CartVect(abs_iter_tol), abs_iter_tol, corner_ijk+3);
+ if (MB_SUCCESS != rval) return rval;
//step 4
//set up TLreg_o
@@ -102,17 +112,12 @@ namespace moab
int dest;
double boxtosend[6];
- BoundBox box;
- ErrorCode result = myTree->get_bounding_box(box);
- if (result != MB_SUCCESS)
- return result;
-
- box.get(boxtosend);
+ localBox.get(boxtosend);
//iterate over all regions overlapping with my bounding box using the computerd corner IDs
- for (int k = cornerIJK[2]; k <= cornerIJK[5]; k++) {
- for (int j = cornerIJK[1]; j <= cornerIJK[4]; j++) {
- for (int i = cornerIJK[0]; i <= cornerIJK[3]; i++) {
+ for (int k = corner_ijk[2]; k <= corner_ijk[5]; k++) {
+ for (int j = corner_ijk[1]; j <= corner_ijk[4]; j++) {
+ for (int i = corner_ijk[0]; i <= corner_ijk[3]; i++) {
dest = k * regNums[0]*regNums[1] + j * regNums[0] + i;
TLreg_o.push_back(&dest, NULL, NULL, boxtosend);
}
@@ -144,7 +149,7 @@ namespace moab
return MB_UNSUPPORTED_OPERATION;
}
- bool is_initialized(int i) {return (i == -1);}
+ bool is_neg(int i) {return (i == -1);}
ErrorCode SpatialLocator::par_locate_points(ParallelComm *pc,
const double *pos, int num_points,
@@ -160,14 +165,10 @@ namespace moab
// steps 1-2 - initialize the alternative decomposition box from global box
- rval = initialize_global_box(pc);
+ rval = initialize_intermediate_partition(pc);
- //step 3 - compute the IDs of the regions which contain each corner of local box
- rval = compute_corner_ijks();
- if (rval != MB_SUCCESS) return rval;
-
- //steps 4-6 - set up TLreg_o, gs_transfer, gather registrations
- rval = register_with_forwarders(pc, TLreg_o);
+ //steps 3-6 - set up TLreg_o, gs_transfer, gather registrations
+ rval = register_src_with_intermediate_procs(pc, abs_iter_tol, TLreg_o);
if (rval != MB_SUCCESS) return rval;
// actual parallel point location using intermediate partition
@@ -190,7 +191,7 @@ namespace moab
for (int pnt=0; pnt < 3*num_points; pnt+=3)
{
- int forw_id = proc_from_point(pos+pnt); //get ID of proc resonsible of the region the proc is in
+ int forw_id = proc_from_point(pos+pnt, abs_iter_tol); //get ID of proc resonsible of the region the proc is in
iargs[0] = forw_id; //toProc
iargs[1] = my_rank; //originalSourceProc
@@ -240,7 +241,7 @@ namespace moab
//TLsearch_results_o
//TL: (OriginalSourceProc, targetIndex, sourceIndex, U,V,W);
- TLsearch_results_o.initialize(3,0,0,3,0);
+ TLsearch_results_o.initialize(3,0,0,0,0);
//step 13 is done in test_local_box
@@ -254,6 +255,8 @@ namespace moab
if (MB_SUCCESS != rval)
return rval;
+ locTable.initialize(1, 0, 1, 3, 0);
+ locTable.enableWriteAccess();
for (int i = 0; i < NN; i++) {
if (is_inside[i]) {
iargs[0] = TLforward_o.vi_rd[3*i+1];
@@ -263,6 +266,7 @@ namespace moab
locTable.push_back(const_cast<int*>(&TLforward_o.vi_rd[3*i+1]), NULL, &ents[i], ¶ms[3*i]);
}
}
+ locTable.disableWriteAccess();
//step 14: send TLsearch_results_o and receive TLloc_i
if (pc)
@@ -282,8 +286,9 @@ namespace moab
if (debug) {
int num_found = num_points - 0.5 *
- std::count_if(parLocTable.vi_wr, parLocTable.vi_wr + 2*num_points, is_initialized);
- std::cout << "Points found = " << num_found << " out of " << num_points << "." << std::endl;
+ std::count_if(parLocTable.vi_wr, parLocTable.vi_wr + 2*num_points, is_neg);
+ std::cout << "Points found = " << num_found << "/" << num_points
+ << " (" << 100.0*((double)num_found/num_points) << "%)" << std::endl;
}
return MB_SUCCESS;
@@ -343,9 +348,7 @@ namespace moab
double tmp_abs_iter_tol = abs_iter_tol;
if (rel_iter_tol && !tmp_abs_iter_tol) {
// relative epsilon given, translate to absolute epsilon using box dimensions
- BoundBox box;
- myTree->get_bounding_box(box);
- tmp_abs_iter_tol = rel_iter_tol * box.diagonal_length();
+ tmp_abs_iter_tol = rel_iter_tol * localBox.diagonal_length();
}
EntityHandle closest_leaf;
diff --git a/src/moab/SpatialLocator.hpp b/src/moab/SpatialLocator.hpp
index 00a7857..0307b1c 100644
--- a/src/moab/SpatialLocator.hpp
+++ b/src/moab/SpatialLocator.hpp
@@ -60,7 +60,10 @@ namespace moab {
ErrorCode add_elems(Range &elems);
/* get bounding box of this locator */
- ErrorCode get_bounding_box(BoundBox &box);
+ BoundBox &local_box() {return localBox;}
+
+ /* get bounding box of this locator */
+ const BoundBox &local_box() const {return localBox;}
/* locate a set of vertices, Range variant */
ErrorCode locate_points(Range &vertices,
@@ -169,28 +172,25 @@ namespace moab {
#ifdef USE_MPI
/* MPI_ReduceAll source mesh bounding boxes to get global source mesh bounding box
*/
- ErrorCode initialize_global_box(ParallelComm *pc);
-
- /* compute the ijk location of this proc's target mesh bounding box based on intermediate decomp
- * extents/divisions
- */
- inline ErrorCode compute_corner_ijks();
+ ErrorCode initialize_intermediate_partition(ParallelComm *pc);
- /* for a given point in space, compute its ijk location in the intermediate decomposition
+ /* for a given point in space, compute its ijk location in the intermediate decomposition; tolerance is
+ * used only to test proximity to global box extent, not for local box test
*/
- inline ErrorCode get_point_ijk(const CartVect &point, int *ijk) const;
+ inline ErrorCode get_point_ijk(const CartVect &point, const double abs_iter_tol, int *ijk) const;
/* given an ijk location in the intermediate partition, return the proc rank for that location
*/
inline int proc_from_ijk(const int *ijk) const;
- /* given a point in space, return the proc responsible for that point from the intermediate decomp
+ /* given a point in space, return the proc responsible for that point from the intermediate decomp; no tolerances
+ * applied here, so first proc in lexicographic ijk ordering is returned
*/
- inline int proc_from_point(const double *pos) const;
+ inline int proc_from_point(const double *pos, const double abs_iter_tol) const;
/* register my source mesh with intermediate decomposition procs
*/
- ErrorCode register_with_forwarders(ParallelComm *pc, TupleList &TLreg_o);
+ ErrorCode register_src_with_intermediate_procs(ParallelComm *pc, double abs_iter_tol, TupleList &TLreg_o);
#endif
@@ -235,6 +235,10 @@ namespace moab {
*/
TupleList parLocTable;
+ /* \brief Local bounding box, duplicated from tree
+ */
+ BoundBox localBox;
+
/* \brief Global bounding box, used in parallel spatial location
*/
BoundBox globalBox;
@@ -247,10 +251,6 @@ namespace moab {
*/
int regNums[3];
- /* \brief Region numbers for my box's min/max corners, used in parallel spatial location
- */
- int cornerIJK[6];
-
/* \brief Map from source processor to bounding box of that proc's source mesh
*
*/
@@ -271,29 +271,19 @@ namespace moab {
return locate_points(pos, 1, &ent, params, is_inside, rel_iter_tol, abs_iter_tol, inside_tol);
}
- inline ErrorCode SpatialLocator::get_bounding_box(BoundBox &box)
- {
- return myTree->get_bounding_box(box);
- }
-
- inline ErrorCode SpatialLocator::compute_corner_ijks() {
-
- BoundBox box;
- ErrorCode rval = myTree->get_bounding_box(box);
- if (rval != MB_SUCCESS) return rval;
-
- rval = get_point_ijk(box.bMin, cornerIJK);
- if (MB_SUCCESS != rval) return rval;
- rval = get_point_ijk(box.bMax, cornerIJK+3);
- return rval;
- }
-
- inline ErrorCode SpatialLocator::get_point_ijk(const CartVect &point, int *ijk) const
+ inline ErrorCode SpatialLocator::get_point_ijk(const CartVect &point, const double abs_iter_tol, int *ijk) const
{
- for (int i = 0; i < 3; i++)
- ijk[i] = (regNums[i] > 1 ? (point[i] - globalBox.bMin[i]) / regDeltaXYZ[i] : 0);
+ for (int i = 0; i < 3; i++) {
+ if (point[i] < globalBox.bMin[i]-abs_iter_tol || point[i] > globalBox.bMax[i]+abs_iter_tol)
+ ijk[i] = -1;
+ else {
+ ijk[i] = point[i] - globalBox.bMin[i] / regDeltaXYZ[i];
+ if (ijk[i] >= regNums[i] && point[i] <= globalBox.bMax[i]+abs_iter_tol)
+ ijk[i] = regNums[i]-1;
+ }
+ }
- return MB_SUCCESS;
+ return (ijk[0] >= 0 && ijk[1] >= 0 && ijk[2] >= 0 ? MB_SUCCESS : MB_FAILURE);;
}
inline int SpatialLocator::proc_from_ijk(const int *ijk) const
@@ -301,10 +291,10 @@ namespace moab {
return ijk[2] * regNums[0]*regNums[1] + ijk[1] * regNums[0] + ijk[0];
}
- inline int SpatialLocator::proc_from_point(const double *pos) const
+ inline int SpatialLocator::proc_from_point(const double *pos, const double abs_iter_tol) const
{
int ijk[3];
- ErrorCode rval = get_point_ijk(CartVect(pos), ijk);
+ ErrorCode rval = get_point_ijk(CartVect(pos), abs_iter_tol, ijk);
if (MB_SUCCESS != rval) return -1;
return ijk[2] * regNums[0]*regNums[1] + ijk[1] * regNums[0] + ijk[0];
diff --git a/test/parallel/par_spatial_locator_test.cpp b/test/parallel/par_spatial_locator_test.cpp
index 40cb98d..334f958 100644
--- a/test/parallel/par_spatial_locator_test.cpp
+++ b/test/parallel/par_spatial_locator_test.cpp
@@ -28,6 +28,7 @@ ErrorCode load_file(Interface &mb, std::string &fn, Range &elems);
int max_depth = 30;
int npoints = 1000;
int leaf = 6;
+int tree = -1;
bool print_tree = false;
int ints = 10;
std::string fname;
@@ -44,9 +45,12 @@ int main(int argc, char **argv)
po.addOpt<int>( "max_depth,m", "Maximum depth of tree", &max_depth);
po.addOpt<int>( "npoints,n", "Number of query points", &npoints);
po.addOpt<void>( "print,p", "Print tree details", &print_tree);
+ po.addOpt<int>( "tree,t", "Tree type (-1=all (default), 0=AdaptiveKD, 1=BVH", &tree);
po.parseCommandLine(argc, argv);
- RUN_TEST(test_kd_tree);
+ if (-1 == tree || 0 == tree)
+ RUN_TEST(test_kd_tree);
+ if (-1 == tree || 1 == tree)
RUN_TEST(test_bvh_tree);
fail = MPI_Finalize();
@@ -113,8 +117,7 @@ bool is_neg(int is_neg)
void test_locator(SpatialLocator *sl)
{
CartVect box_del;
- BoundBox box;
- ErrorCode rval = sl->get_bounding_box(box); CHECK_ERR(rval);
+ BoundBox box = sl->local_box();
box_del = box.bMax - box.bMin;
double denom = 1.0 / (double)RAND_MAX;
@@ -130,12 +133,13 @@ void test_locator(SpatialLocator *sl)
ParallelComm *pc = ParallelComm::get_pcomm(sl->moab(), 0);
CHECK(pc != NULL);
- rval = sl->par_locate_points(pc, test_pts[0].array(), npoints); CHECK_ERR(rval);
+ ErrorCode rval = sl->par_locate_points(pc, test_pts[0].array(), npoints); CHECK_ERR(rval);
if (pc->rank() == 0) {
int num_out = std::count_if(sl->par_loc_table().vi_rd, sl->par_loc_table().vi_rd+2*npoints, is_neg);
num_out /= 2;
- std::cout << "Number of points inside an element = " << npoints-num_out << " out of " << npoints << std::endl;
+ std::cout << "Number of points inside an element = " << npoints-num_out << "/" << npoints
+ << " (" << 100.0*((double)npoints-num_out)/npoints << "%)" << std::endl;
std::cout << "Traversal stats:" << std::endl;
sl->get_tree()->tree_stats().output();
@@ -152,10 +156,11 @@ ErrorCode create_hex_mesh(Interface &mb, Range &elems, int n, int dim)
ScdInterface *scdi;
ErrorCode rval = mb.query_interface(scdi); CHECK_ERR(rval);
ScdParData spd;
- spd.gDims[0] = n-1;
+ spd.gDims[0] = spd.gDims[1] = spd.gDims[2] = 0;
+ spd.gDims[3] = n;
spd.partMethod = ScdParData::SQIJK;
- if (dim > 1) spd.gDims[1] = n-1;
- if (dim > 2) spd.gDims[2] = n-1;
+ if (dim > 1) spd.gDims[4] = n;
+ if (dim > 2) spd.gDims[5] = n;
ScdBox *new_box;
rval = scdi->construct_box(HomCoord(0, 0, 0), HomCoord(0, 0, 0),
NULL, 0, new_box, NULL, &spd, false, 0); CHECK_ERR(rval);
Repository URL: https://bitbucket.org/fathomteam/moab/
--
This is a commit notification from bitbucket.org. You are receiving
this because you have the service enabled, addressing the recipient of
this email.
More information about the moab-dev
mailing list