[MOAB-dev] commit/MOAB: 2 new changesets
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Fri Oct 11 17:49:57 CDT 2013
2 new commits in MOAB:
https://bitbucket.org/fathomteam/moab/commits/50b351b17830/
Changeset: 50b351b17830
Branch: None
User: tautges
Date: 2013-10-07 23:31:36
Summary: DeformMeshRemap: add some of the data interpolation calls, basically filling
out the algorithm outlined on confluence.
makefile: adding MOAB lib as a dependency, and mbcoupler lib
to libs for DeformMeshRemap
ElemEvaluator: hardwiring the number of tuples in a few cases where we know
coords is the field we're evaluating; also changing the name of a few params from
tag_dim to tagged_ent_dim to be more descriptive
LinearQuad: fixed jacobian scaling so J(2,2) is 1 instead of .25 (otherwise Jacobian for
ideal cube is wrong)
SpatialLocator: added par_locate_points function, with no definition yet; also commented more fully,
to describe use of SpatialLocator and assumptions made, and distinguishing serial vs. parallel use
test/elem_eval_test: added test for linear quad, and corrected test for avg position to divide by elem volume
test/spatial_locator_test: corrected arg order for SpatialLocator::locate_points
tools/mbcoupler/DataCoupler: new class analagous to Coupler that takes advantage of new SpatialLocator functionality
Affected #: 12 files
diff --git a/examples/DeformMeshRemap.cpp b/examples/DeformMeshRemap.cpp
index fac83ae..500d969 100644
--- a/examples/DeformMeshRemap.cpp
+++ b/examples/DeformMeshRemap.cpp
@@ -15,6 +15,7 @@
#include "moab/LloydSmoother.hpp"
#include "moab/ProgOptions.hpp"
#include "MBTagConventions.hpp"
+#include "DataCoupler.hpp"
#include <iostream>
#include <set>
@@ -84,7 +85,7 @@ public:
private:
//! apply a known deformation to the solid elements, putting the results in the xNew tag; also
//! write current coordinates to the xNew tag for fluid elements
- ErrorCode deform_master(Range &fluid_elems, Range &solid_elems);
+ ErrorCode deform_master(Range &fluid_elems, Range &solid_elems, const char *tag_name = NULL);
//! read a file and establish proper ranges
ErrorCode read_file(int m_or_s, string &fname, EntityHandle &seth);
@@ -193,20 +194,46 @@ ErrorCode DeformMeshRemap::execute()
if (MB_SUCCESS != rval) return rval;
// deform the master's solid mesh, put results in a new tag
- rval = deform_master(fluidElems[MASTER], solidElems[MASTER]); RR("");
+ rval = deform_master(fluidElems[MASTER], solidElems[MASTER], "xnew"); RR("");
if (debug) write_and_save(solidElems[MASTER], masterSet, xNew, "deformed.vtk");
-
- // smooth the master mesh
- LloydSmoother *ll = new LloydSmoother(mbImpl, NULL, fluidElems[MASTER], xNew);
- rval = ll->perform_smooth();
- RR("Failed in lloyd smoothing.");
- cout << "Lloyd smoothing required " << ll->num_its() << " iterations." << endl;
- if (debug) write_and_save(fluidElems[MASTER], masterSet, xNew, "smoothed.vtk");
+ { // to isolate the lloyd smoother & delete when done
+ // smooth the master mesh
+ LloydSmoother ll(mbImpl, NULL, fluidElems[MASTER], xNew);
+ rval = ll.perform_smooth();
+ RR("Failed in lloyd smoothing.");
+ cout << "Lloyd smoothing required " << ll.num_its() << " iterations." << endl;
+ }
+
// map new locations to slave
+ // locate slave vertices in master, orig coords; do this with a data coupler, so you can
+ // later interpolate
+ Range src_elems = solidElems[MASTER];
+ src_elems.merge(fluidElems[MASTER]);
+
+ // initialize data coupler on source elements
+ DataCoupler dc_master(mbImpl, NULL, src_elems, 0);
- delete ll;
+ Range tgt_verts, tmp_range = solidElems[SLAVE];
+ tmp_range.merge(fluidElems[SLAVE]);
+ rval = mbImpl->get_adjacencies(tmp_range, 0, false, tgt_verts, Interface::UNION);
+ RR("Failed to get target verts.");
+
+ // locate slave vertices, caching results in dc
+ rval = dc_master.locate_points(tgt_verts, 0.0, 1e-10); RR("Point location of tgt verts failed.");
+
+ // interpolate xNew to slave points
+ rval = dc_master.interpolate((int)DataCoupler::VOLUME, "xnew"); RR("Failed to interpolate target solution.");
+
+ // transfer xNew to coords, for master and slave
+ rval = write_to_coords(fluidElems[MASTER], xNew); RR("Failed writing tag to master fluid verts.");
+ rval = write_to_coords(tgt_verts, xNew); RR("Failed writing tag to slave verts.");
+
+ if (debug) {
+ rval = mbImpl->write_file("smoothed_master.vtk", NULL, NULL, &masterSet, 1);
+ rval = mbImpl->write_file("slave_interp.vtk", NULL, NULL, &slaveSet, 1);
+ }
return MB_SUCCESS;
}
@@ -302,12 +329,12 @@ void deform_func(double *xold, double *xnew)
xnew[1] = xold[1] + yfrac * (1.0 + xfrac) * 0.05;
}
-ErrorCode DeformMeshRemap::deform_master(Range &fluid_elems, Range &solid_elems)
+ErrorCode DeformMeshRemap::deform_master(Range &fluid_elems, Range &solid_elems, const char *tag_name)
{
// deform elements with an analytic function
// create the tag
- ErrorCode rval = mbImpl->tag_get_handle("", 3, MB_TYPE_DOUBLE, xNew, MB_TAG_CREAT|MB_TAG_DENSE);
+ ErrorCode rval = mbImpl->tag_get_handle((tag_name ? tag_name : ""), 3, MB_TYPE_DOUBLE, xNew, MB_TAG_CREAT|MB_TAG_DENSE);
RR("Failed to create xnew tag.");
// get all the vertices and coords in the fluid, set xnew to them
diff --git a/examples/makefile b/examples/makefile
index 87c0132..773c443 100644
--- a/examples/makefile
+++ b/examples/makefile
@@ -14,43 +14,43 @@ EXOIIEXAMPLES = TestExodusII
default: ${EXAMPLES}
-HelloMOAB : HelloMOAB.o
+HelloMOAB : HelloMOAB.o ${MOAB_LIBDIR}/libMOAB.la
${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK}
-GetEntities: GetEntities.o
+GetEntities: GetEntities.o ${MOAB_LIBDIR}/libMOAB.la
${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK}
-SetsNTags: SetsNTags.o
+SetsNTags: SetsNTags.o ${MOAB_LIBDIR}/libMOAB.la
${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK}
-LloydRelaxation: LloydRelaxation.o
+LloydRelaxation: LloydRelaxation.o ${MOAB_LIBDIR}/libMOAB.la
${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK}
-StructuredMeshSimple : StructuredMeshSimple.o
+StructuredMeshSimple : StructuredMeshSimple.o ${MOAB_LIBDIR}/libMOAB.la
${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK}
-DirectAccessWithHoles: DirectAccessWithHoles.o
+DirectAccessWithHoles: DirectAccessWithHoles.o ${MOAB_LIBDIR}/libMOAB.la
${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK}
-DirectAccessNoHoles: DirectAccessNoHoles.o
+DirectAccessNoHoles: DirectAccessNoHoles.o ${MOAB_LIBDIR}/libMOAB.la
${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK}
-DirectAccessNoHolesF90: DirectAccessNoHolesF90.o
+DirectAccessNoHolesF90: DirectAccessNoHolesF90.o ${MOAB_LIBDIR}/libMOAB.la
${MOAB_CXX} -o $@ $< ${IMESH_LIBS}
-ReduceExchangeTags : ReduceExchangeTags.o
+ReduceExchangeTags : ReduceExchangeTags.o ${MOAB_LIBDIR}/libMOAB.la
${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK}
-HelloParMOAB: HelloParMOAB.o
+HelloParMOAB: HelloParMOAB.o ${MOAB_LIBDIR}/libMOAB.la
${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK}
-TestExodusII: TestExodusII.o
+TestExodusII: TestExodusII.o ${MOAB_LIBDIR}/libMOAB.la
${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK}
-point_in_elem_search: point_in_elem_search.o
+point_in_elem_search: point_in_elem_search.o ${MOAB_LIBDIR}/libMOAB.la
-DeformMeshRemap: DeformMeshRemap.o
- ${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK}
+DeformMeshRemap: DeformMeshRemap.o ${MOAB_LIBDIR}/libMOAB.la
+ ${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK} -lmbcoupler ${MOAB_LIBS_LINK}
clean:
rm -rf *.o ${EXAMPLES} ${PAREXAMPLES} ${EXOIIEXAMPLES}
diff --git a/src/LocalDiscretization/ElemEvaluator.cpp b/src/LocalDiscretization/ElemEvaluator.cpp
index b132da9..4446281 100644
--- a/src/LocalDiscretization/ElemEvaluator.cpp
+++ b/src/LocalDiscretization/ElemEvaluator.cpp
@@ -30,7 +30,9 @@ namespace moab {
CartVect new_pos;
// evaluate that first guess to get a new position
- ErrorCode rval = (*eval)(cvparams->array(), verts, ndim, ndim, work, new_pos.array());
+ ErrorCode rval = (*eval)(cvparams->array(), verts, ndim,
+ 3, // hardwire to num_tuples to 3 since the field is coords
+ work, new_pos.array());
if (MB_SUCCESS != rval) return rval;
// residual is diff between old and new pos; need to minimize that
@@ -58,7 +60,9 @@ namespace moab {
*cvparams -= J.inverse(1.0/det) * res;
// get the new forward-evaluated position, and its difference from the target pt
- rval = (*eval)(params, verts, ndim, ndim, work, new_pos.array());
+ rval = (*eval)(params, verts, ndim,
+ 3, // hardwire to num_tuples to 3 since the field is coords
+ work, new_pos.array());
if (MB_SUCCESS != rval) return rval;
res = new_pos - *cvposn;
}
diff --git a/src/LocalDiscretization/LinearQuad.cpp b/src/LocalDiscretization/LinearQuad.cpp
index ddd4918..93237f9 100644
--- a/src/LocalDiscretization/LinearQuad.cpp
+++ b/src/LocalDiscretization/LinearQuad.cpp
@@ -16,23 +16,23 @@ namespace moab
*/
const double LinearQuad::gauss[1][2] = { { 2.0, 0.0 } };
- ErrorCode LinearQuad::jacobianFcn(const double *params, const double *verts, const int /*nverts*/, const int ndim,
+ ErrorCode LinearQuad::jacobianFcn(const double *params, const double *verts, const int /*nverts*/, const int /*ndim*/,
double *, double *result)
{
Matrix3 *J = reinterpret_cast<Matrix3*>(result);
*J = Matrix3(0.0);
for (unsigned i = 0; i < 4; ++i) {
- const double params_p = 1 + params[0]*corner[i][0];
+ const double xi_p = 1 + params[0]*corner[i][0];
const double eta_p = 1 + params[1]*corner[i][1];
- const double dNi_dparams = corner[i][0] * eta_p;
- const double dNi_deta = corner[i][1] * params_p;
- (*J)(0,0) += dNi_dparams * verts[i*ndim+0];
- (*J)(1,0) += dNi_dparams * verts[i*ndim+1];
- (*J)(0,1) += dNi_deta * verts[i*ndim+0];
- (*J)(1,1) += dNi_deta * verts[i*ndim+1];
+ const double dNi_dxi = corner[i][0] * eta_p;
+ const double dNi_deta = corner[i][1] * xi_p;
+ (*J)(0,0) += dNi_dxi * verts[i*3+0];
+ (*J)(1,0) += dNi_dxi * verts[i*3+1];
+ (*J)(0,1) += dNi_deta * verts[i*3+0];
+ (*J)(1,1) += dNi_deta * verts[i*3+1];
}
- (*J)(2,2) = 1.0; /* to make sure the Jacobian determinant is non-zero */
(*J) *= 0.25;
+ (*J)(2,2) = 1.0; /* to make sure the Jacobian determinant is non-zero */
return MB_SUCCESS;
}// LinearQuad::jacobian()
@@ -62,7 +62,7 @@ namespace moab
for(unsigned int j2 = 0; j2 < LinearQuad::gauss_count; ++j2) {
x[1] = LinearQuad::gauss[j2][1];
double w2 = LinearQuad::gauss[j2][0];
- rval = evalFcn(x.array(),field, ndim, num_tuples, NULL, tmp_result);
+ rval = evalFcn(x.array(), field, ndim, num_tuples, NULL, tmp_result);
if (MB_SUCCESS != rval) return rval;
rval = jacobianFcn(x.array(), verts, nverts, ndim, work, J[0]);
if (MB_SUCCESS != rval) return rval;
diff --git a/src/LocalDiscretization/moab/ElemEvaluator.hpp b/src/LocalDiscretization/moab/ElemEvaluator.hpp
index 48d34ac..eac008d 100644
--- a/src/LocalDiscretization/moab/ElemEvaluator.hpp
+++ b/src/LocalDiscretization/moab/ElemEvaluator.hpp
@@ -119,9 +119,9 @@ namespace moab {
* \param impl MOAB instance
* \param ent Entity handle to cache on the evaluator
* \param tag Tag to cache on the evaluator
- * \param tag_dim Tag dimension to cache on the evaluator
+ * \param tagged_ent_dim Dimension of entities to be tagged to cache on the evaluator
*/
- ElemEvaluator(Interface *impl, EntityHandle ent = 0, Tag tag = 0, int tag_dim = -1);
+ ElemEvaluator(Interface *impl, EntityHandle ent = 0, Tag tag = 0, int tagged_ent_dim = -1);
/** \brief Evaluate cached tag at a given parametric location within the cached entity
* If evaluating coordinates, call set_tag(0, 0), which indicates coords instead of a tag.
@@ -233,23 +233,23 @@ namespace moab {
* To designate that vertex coordinates are the desired tag, pass in a tag handle of 0
* and a tag dimension of 0.
* \param tag Tag handle to cache, or 0 to cache vertex positions
- * \param tag_dim Dimension of entities tagged with this tag
+ * \param tagged_ent_dim Dimension of entities tagged with this tag
*/
- inline ErrorCode set_tag_handle(Tag tag, int tag_dim = -1);
+ inline ErrorCode set_tag_handle(Tag tag, int tagged_ent_dim = -1);
/* \brief Set the name of the tag to cache on this evaluator
* To designate that vertex coordinates are the desired tag, pass in "COORDS" as the name
* and a tag dimension of 0.
- * \param tag Tag handle to cache, or 0 to cache vertex positions
- * \param tag_dim Dimension of entities tagged with this tag
+ * \param tag_name Tag handle to cache, or 0 to cache vertex positions
+ * \param tagged_ent_dim Dimension of entities tagged with this tag
*/
- inline ErrorCode set_tag(const char *tag_name, int ent_dim = -1);
+ inline ErrorCode set_tag(const char *tag_name, int tagged_ent_dim = -1);
/* \brief Get the dimension of the entities on which tag is set */
- inline int get_tag_dim() const {return tagDim;};
+ inline int get_tagged_ent_dim() const {return taggedEntDim;};
/* \brief Set the dimension of entities having the tag */
- inline ErrorCode set_tag_dim(int dim);
+ inline ErrorCode set_tagged_ent_dim(int dim);
/* \brief MOAB interface cached on this evaluator */
inline Interface *get_moab() {return mbImpl;}
@@ -287,7 +287,7 @@ namespace moab {
int numTuples;
/** \brief Dimension of entities from which to grab tag */
- int tagDim;
+ int taggedEntDim;
/** \brief Tag space */
std::vector<unsigned char> tagSpace;
@@ -300,13 +300,13 @@ namespace moab {
}; // class ElemEvaluator
- inline ElemEvaluator::ElemEvaluator(Interface *impl, EntityHandle ent, Tag tag, int tag_dim)
+ inline ElemEvaluator::ElemEvaluator(Interface *impl, EntityHandle ent, Tag tag, int tagged_ent_dim)
: mbImpl(impl), entHandle(0), entType(MBMAXTYPE), entDim(-1), numVerts(0),
vertHandles(NULL), tagHandle(0), tagCoords(false), numTuples(0),
- tagDim(0), workSpace(NULL)
+ taggedEntDim(0), workSpace(NULL)
{
if (ent) set_ent_handle(ent);
- if (tag) set_tag_handle(tag, tag_dim);
+ if (tag) set_tag_handle(tag, tagged_ent_dim);
}
inline ErrorCode ElemEvaluator::set_ent_handle(EntityHandle ent)
@@ -329,18 +329,17 @@ namespace moab {
rval = set_tag_handle(tagHandle);
if (MB_SUCCESS != rval) return rval;
}
-
if (evalSets[entType].initFcn) return (*evalSets[entType].initFcn)(vertPos[0].array(), numVerts, workSpace);
return MB_SUCCESS;
}
- inline ErrorCode ElemEvaluator::set_tag_handle(Tag tag, int tag_dim)
+ inline ErrorCode ElemEvaluator::set_tag_handle(Tag tag, int tagged_ent_dim)
{
ErrorCode rval = MB_SUCCESS;
- if (!tag && !tag_dim) {
+ if (!tag && !tagged_ent_dim) {
tagCoords = true;
numTuples = 3;
- tagDim = 0;
+ taggedEntDim = 0;
tagHandle = 0;
return rval;
}
@@ -355,14 +354,14 @@ namespace moab {
tagCoords = false;
}
- tagDim = (-1 == tag_dim ? 0 : tag_dim);
+ taggedEntDim = (-1 == tagged_ent_dim ? 0 : tagged_ent_dim);
if (entHandle) {
- if (0 == tagDim) {
+ if (0 == taggedEntDim) {
rval = mbImpl->tag_get_data(tagHandle, vertHandles, numVerts, &tagSpace[0]);
if (MB_SUCCESS != rval) return rval;
}
- else if (tagDim == entDim) {
+ else if (taggedEntDim == entDim) {
rval = mbImpl->tag_get_data(tagHandle, &entHandle, 1, &tagSpace[0]);
if (MB_SUCCESS != rval) return rval;
}
@@ -371,14 +370,14 @@ namespace moab {
return rval;
}
- inline ErrorCode ElemEvaluator::set_tag(const char *tag_name, int tag_dim)
+ inline ErrorCode ElemEvaluator::set_tag(const char *tag_name, int tagged_ent_dim)
{
ErrorCode rval = MB_SUCCESS;
if (!tag_name || strlen(tag_name) == 0) return MB_FAILURE;
Tag tag;
if (!strcmp(tag_name, "COORDS")) {
tagCoords = true;
- tagDim = 0;
+ taggedEntDim = 0;
numTuples = 3;
tagHandle = 0;
// can return here, because vertex coords already cached when entity handle set
@@ -399,15 +398,15 @@ namespace moab {
tagCoords = false;
}
- tagDim = (-1 == tag_dim ? entDim : tag_dim);
+ taggedEntDim = (-1 == tagged_ent_dim ? entDim : tagged_ent_dim);
}
if (entHandle) {
- if (0 == tagDim) {
+ if (0 == taggedEntDim) {
rval = mbImpl->tag_get_data(tagHandle, vertHandles, numVerts, &tagSpace[0]);
if (MB_SUCCESS != rval) return rval;
}
- else if (tagDim == entDim) {
+ else if (taggedEntDim == entDim) {
rval = mbImpl->tag_get_data(tagHandle, &entHandle, 1, &tagSpace[0]);
if (MB_SUCCESS != rval) return rval;
}
@@ -431,16 +430,16 @@ namespace moab {
assert(entHandle && MBMAXTYPE != entType);
return (*evalSets[entType].evalFcn)(params,
(tagCoords ? (const double*) vertPos[0].array(): (const double*)&tagSpace[0]),
- entDim, (-1 == num_tuples ? numTuples : num_tuples),
- workSpace, result);
+ entDim,
+ (-1 == num_tuples ? numTuples : num_tuples), workSpace, result);
}
inline ErrorCode ElemEvaluator::reverse_eval(const double *posn, const double tol, double *params, bool *ins) const
{
assert(entHandle && MBMAXTYPE != entType);
return (*evalSets[entType].reverseEvalFcn)(evalSets[entType].evalFcn, evalSets[entType].jacobianFcn, evalSets[entType].insideFcn,
- posn, vertPos[0].array(), numVerts, entDim, tol, workSpace,
- params, ins);
+ posn, vertPos[0].array(), numVerts,
+ entDim, tol, workSpace, params, ins);
}
/** \brief Evaluate the jacobian of the cached entity at a given parametric location */
@@ -456,7 +455,7 @@ namespace moab {
assert(entHandle && MBMAXTYPE != entType && (tagCoords || tagHandle));
ErrorCode rval = MB_SUCCESS;
if (!tagCoords) {
- if (0 == tagDim) rval = mbImpl->tag_get_data(tagHandle, vertHandles, numVerts, (void*)&tagSpace[0]);
+ if (0 == taggedEntDim) rval = mbImpl->tag_get_data(tagHandle, vertHandles, numVerts, (void*)&tagSpace[0]);
else rval = mbImpl->tag_get_data(tagHandle, &entHandle, 1, (void*)&tagSpace[0]);
if (MB_SUCCESS != rval) return rval;
}
diff --git a/src/SpatialLocator.cpp b/src/SpatialLocator.cpp
index a543dc9..d3b521e 100644
--- a/src/SpatialLocator.cpp
+++ b/src/SpatialLocator.cpp
@@ -31,10 +31,62 @@ namespace moab
return MB_SUCCESS;
}
+#ifdef USE_MPI
+ ErrorCode SpatialLocator::par_locate_points(Range &/*vertices*/,
+ double /*rel_tol*/, double /*abs_tol*/)
+ {
+ return MB_UNSUPPORTED_OPERATION;
+ }
+
+ ErrorCode SpatialLocator::par_locate_points(const double */*pos*/, int /*num_points*/,
+ double /*rel_tol*/, double /*abs_tol*/)
+ {
+ return MB_UNSUPPORTED_OPERATION;
+ }
+#endif
+
+ ErrorCode SpatialLocator::locate_points(Range &verts,
+ double rel_eps, double abs_eps)
+ {
+ assert(!verts.empty() && mbImpl->type_from_handle(*verts.rbegin()) == MBVERTEX);
+ std::vector<double> pos(3*verts.size());
+ ErrorCode rval = mbImpl->get_coords(verts, &pos[0]);
+ if (MB_SUCCESS != rval) return rval;
+ rval = locate_points(&pos[0], verts.size(), rel_eps, abs_eps);
+ if (MB_SUCCESS != rval) return rval;
+
+ return MB_SUCCESS;
+ }
+
+ ErrorCode SpatialLocator::locate_points(const double *pos, int num_points,
+ double rel_eps, double abs_eps)
+ {
+ // initialize to tuple structure (p_ui, hs_ul, r[3]_d) (see header comments for locTable)
+ locTable.initialize(1, 0, 1, 3, num_points);
+ locTable.enableWriteAccess();
+
+ // pass storage directly into locate_points, since we know those arrays are contiguous
+ ErrorCode rval = locate_points(pos, num_points, locTable.vul_wr, locTable.vr_wr, NULL, rel_eps, abs_eps);
+ std::fill(locTable.vi_wr, locTable.vi_wr+num_points, 0);
+ if (MB_SUCCESS != rval) return rval;
+
+ return MB_SUCCESS;
+ }
+
+ ErrorCode SpatialLocator::locate_points(Range &verts,
+ EntityHandle *ents, double *params, bool *is_inside,
+ double rel_eps, double abs_eps)
+ {
+ assert(!verts.empty() && mbImpl->type_from_handle(*verts.rbegin()) == MBVERTEX);
+ std::vector<double> pos(3*verts.size());
+ ErrorCode rval = mbImpl->get_coords(verts, &pos[0]);
+ if (MB_SUCCESS != rval) return rval;
+ return locate_points(&pos[0], verts.size(), ents, params, is_inside, rel_eps, abs_eps);
+ }
+
ErrorCode SpatialLocator::locate_points(const double *pos, int num_points,
- EntityHandle *ents, double *params,
- double rel_eps, double abs_eps,
- bool *is_inside)
+ EntityHandle *ents, double *params, bool *is_inside,
+ double rel_eps, double abs_eps)
{
if (rel_eps && !abs_eps) {
diff --git a/src/moab/SpatialLocator.hpp b/src/moab/SpatialLocator.hpp
index a8f694f..780c086 100644
--- a/src/moab/SpatialLocator.hpp
+++ b/src/moab/SpatialLocator.hpp
@@ -1,6 +1,30 @@
/**\file SpatialLocator.hpp
* \class moab::SpatialLocator
* \brief Tool to facilitate spatial location of a point in a mesh
+ *
+ * SpatialLocator facilitates searching for points in or performing ray traces on collections of mesh entities
+ * in 2D or 3D. This searching is facilitated by a tree-based decomposition of the mesh. Various types
+ * of trees are implemented in MOAB and can be used by this tool, but by default it uses AdaptiveKDTree
+ * (see child classes of Tree for which others are available). Parallel and serial searching are both
+ * supported.
+ *
+ * SpatialLocator can either cache the search results for points or pass back this information in arguments.
+ * Cached information is kept in locTable, indexed in the same order as search points passed in. This information
+ * consists of the entity handle containing the point and the parametric coordinates inside that element.
+ * Information about the points searched, e.g. the entities from which those points are derived, can be stored
+ * in the calling application if desired.
+ *
+ * In parallel, there is a separation between the proc deciding which points to search for (the "target" proc),
+ * and the proc locating the point in its local mesh (the "source" proc). On the source proc, location
+ * information is cached in locTable, as in the serial case. By default, this location information (handle and
+ * parametric coords) is not returned to the target proc, since it would be of no use there. Instead, the rank
+ * of the source proc locating the point, and the index of that location info in the source proc's locTable, is
+ * returned; this information is stored on the target proc in this class's parLocTable variable. Again,
+ * information about the points searched should be stored in the calling application, if desired.
+ *
+ * This class uses the ElemEvaluator class for specification and evaluation of basis functions (used for computing
+ * parametric coords within an entity). See documentation and examples for that class for usage information.
+ *
*/
#ifndef MOAB_SPATIALLOCATOR_HPP
@@ -9,6 +33,7 @@
#include "moab/Types.hpp"
#include "moab/Tree.hpp"
#include "moab/Range.hpp"
+#include "moab/TupleList.hpp"
#include <string>
#include <vector>
@@ -34,23 +59,83 @@ namespace moab {
/* get bounding box of this locator */
ErrorCode get_bounding_box(BoundBox &box);
+ /* locate a set of vertices, Range variant */
+ ErrorCode locate_points(Range &vertices,
+ EntityHandle *ents, double *params, bool *is_inside = NULL,
+ double rel_tol = 0.0, double abs_tol = 0.0);
+
/* locate a set of points */
ErrorCode locate_points(const double *pos, int num_points,
- EntityHandle *ents, double *params,
- double rel_tol = 0.0, double abs_tol = 0.0,
- bool *is_inside = NULL);
+ EntityHandle *ents, double *params, bool *is_inside = NULL,
+ double rel_tol = 0.0, double abs_tol = 0.0);
- /* locate a point */
- ErrorCode locate_point(const double *pos,
- EntityHandle &ent, double *params,
- double rel_tol = 0.0, double abs_tol = 0.0,
- bool *is_inside = NULL);
+ /* locate a set of vertices or entity centroids, storing results on TupleList in this class
+ * Locate a set of vertices or entity centroids, storing the detailed results in member
+ * variable (TupleList) locTable (see comments on locTable for structure of that tuple).
+ */
+ ErrorCode locate_points(Range &ents,
+ double rel_tol = 0.0, double abs_tol = 0.0);
+
+ /* locate a set of points, storing results on TupleList in this class
+ * Locate a set of points, storing the detailed results in member variable (TupleList) locTable
+ * (see comments on locTable for structure of that tuple).
+ */
+ ErrorCode locate_points(const double *pos, int num_points,
+ double rel_tol = 0.0, double abs_tol = 0.0);
+
+#ifdef USE_MPI
+ /* locate a set of vertices or entity centroids, storing results on TupleList in this class
+ * Locate a set of vertices or entity centroids, 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(Range &vertices,
+ double rel_tol = 0.0, double abs_tol = 0.0);
+
+ /* 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,
+ double rel_tol = 0.0, double abs_tol = 0.0);
+#endif
/* return the tree */
Tree *get_tree() {return myTree;}
+
+ /* get the locTable
+ */
+ TupleList &loc_table() {return locTable;}
+
+ /* get the locTable
+ */
+ const TupleList &loc_table() const {return locTable;}
+
+ /* get the parLocTable
+ */
+ TupleList &par_loc_table() {return parLocTable;}
+
+ /* get the parLocTable
+ */
+ const TupleList &par_loc_table() const {return parLocTable;}
+
+ /* get elemEval */
+ ElemEvaluator *elem_eval() {return elemEval;}
+
+ /* get elemEval */
+ const ElemEvaluator *elem_eval() const {return elemEval;}
+
+ /* set elemEval */
+ void elem_eval(ElemEvaluator *eval) {elemEval = eval;}
private:
+ /* locate a point */
+ ErrorCode locate_point(const double *pos,
+ EntityHandle &ent, double *params, bool *is_inside = NULL,
+ double rel_tol = 0.0, double abs_tol = 0.0);
+
/* MOAB instance */
Interface* mbImpl;
@@ -68,6 +153,29 @@ namespace moab {
/* whether I created the tree or not (determines whether to delete it or not on destruction) */
bool iCreatedTree;
+
+ /* \brief local locations table
+ * This table stores detailed local location data results from locate_points, that is, location data
+ * for points located on the local mesh. Data is stored
+ * in a TupleList, where each tuple consists of (p_i, hs_ul, r[3]_d), where
+ * p_i = (int) proc from which request for this point was made (0 if serial)
+ * hs_ul = (unsigned long) source entity containing the point
+ * r[3]_d = (double) parametric coordinates of the point in hs
+ */
+ TupleList locTable;
+
+ /* \brief parallel locations table
+ * This table stores information about points located on a local or remote processor. For
+ * points located on this processor's local mesh, detailed location data is stored in locTable.
+ * For points located on remote processors, more communication is required to retrieve specific
+ * location data (usually that information isn't useful on this processor).
+ *
+ * The tuple structure of this TupleList is (p_i, ri_i), where:
+ * p_i = (int) processor rank containing this point
+ * ri_i = (int) index into locTable on remote proc containing this point's location information
+ * The indexing of parLocTable corresponds to that of the points/entities passed in.
+ */
+ TupleList parLocTable;
};
inline SpatialLocator::~SpatialLocator()
@@ -76,11 +184,10 @@ namespace moab {
}
inline ErrorCode SpatialLocator::locate_point(const double *pos,
- EntityHandle &ent, double *params,
- double rel_tol, double abs_tol,
- bool *is_inside)
+ EntityHandle &ent, double *params, bool *is_inside,
+ double rel_tol, double abs_tol)
{
- return locate_points(pos, 1, &ent, params, rel_tol, abs_tol, is_inside);
+ return locate_points(pos, 1, &ent, params, is_inside, rel_tol, abs_tol);
}
inline ErrorCode SpatialLocator::get_bounding_box(BoundBox &box)
diff --git a/test/elem_eval_test.cpp b/test/elem_eval_test.cpp
index a5907fd..6f07d7a 100644
--- a/test/elem_eval_test.cpp
+++ b/test/elem_eval_test.cpp
@@ -7,6 +7,7 @@
#include "moab/Core.hpp"
#include "moab/ReadUtilIface.hpp"
#include "moab/ElemEvaluator.hpp"
+#include "moab/LinearQuad.hpp"
#include "moab/LinearHex.hpp"
#include "moab/LinearTet.hpp"
#include "moab/QuadraticHex.hpp"
@@ -21,6 +22,7 @@ std::string TestDir(".");
using namespace moab;
+void test_linear_quad();
void test_linear_hex();
void test_linear_tet();
void test_quadratic_hex();
@@ -49,6 +51,7 @@ void test_eval(ElemEvaluator &ee, bool test_integrate)
bool is_inside;
Matrix3 jacob;
ErrorCode rval;
+ int ent_dim = ee.get_moab()->dimension_from_handle(ee.get_ent_handle());
for (params[0] = -1; params[0] <= 1; params[0] += 0.2) {
for (params[1] = -1; params[1] <= 1; params[1] += 0.2) {
@@ -59,9 +62,11 @@ void test_eval(ElemEvaluator &ee, bool test_integrate)
rval = ee.eval(params.array(), posn.array()); CHECK_ERR(rval);
rval = ee.reverse_eval(posn.array(), EPS1, params2.array(), &is_inside); CHECK_ERR(rval);
- if ((params - params2).length() > 3*EPS1)
- std::cerr << params << std::endl;
- CHECK_REAL_EQUAL(0.0, (params - params2).length(), 3*EPS1);
+ CHECK_REAL_EQUAL(0.0, params[0] - params2[0], 3*EPS1);
+ if (ent_dim > 1)
+ CHECK_REAL_EQUAL(0.0, params[1] - params2[1], 3*EPS1);
+ if (ent_dim > 2)
+ CHECK_REAL_EQUAL(0.0, params[2] - params2[2], 3*EPS1);
// jacobian should be >= 0
rval = ee.jacobian(params.array(), jacob.array()); CHECK_ERR(rval);
@@ -81,8 +86,21 @@ void test_eval(ElemEvaluator &ee, bool test_integrate)
// set that temporary tag on the evaluator so that's what gets integrated
rval = ee.set_tag_handle(tag, 0); CHECK_ERR(rval);
- CartVect integral, avg(0.0);
+ CartVect integral, avg;
rval = ee.integrate(integral.array()); CHECK_ERR(rval);
+
+ // now integrate a const 1-valued function, using direct call to the integrate function
+ std::vector<double> one(ee.get_num_verts(), 1.0);
+ double measure;
+ EvalSet es;
+ EntityHandle eh = ee.get_ent_handle();
+ rval = EvalSet::get_eval_set(ee.get_moab(), eh, es); CHECK_ERR(rval);
+ rval = (*es.integrateFcn)(&one[0], hex_verts[0].array(), ee.get_num_verts(), ee.get_moab()->dimension_from_handle(eh),
+ 1, NULL, &measure); CHECK_ERR(rval);
+ if (measure) integral /= measure;
+
+ // check against avg of entity's vertices' positions
+ rval = ee.get_moab()->get_coords(&eh, 1, avg.array()); CHECK_ERR(rval);
CHECK_REAL_EQUAL(0.0, (avg - integral).length(), EPS1);
// delete the temporary tag
@@ -118,6 +136,7 @@ int main()
{
int failures = 0;
+ failures += RUN_TEST(test_linear_quad);
failures += RUN_TEST(test_linear_hex);
failures += RUN_TEST(test_quadratic_hex);
failures += RUN_TEST(test_linear_tet);
@@ -125,6 +144,23 @@ int main()
return failures;
}
+void test_linear_quad()
+{
+ Core mb;
+ Range verts;
+ ErrorCode rval = mb.create_vertices((double*)hex_verts, 4, verts); CHECK_ERR(rval);
+ EntityHandle quad;
+ std::vector<EntityHandle> connect;
+ std::copy(verts.begin(), verts.end(), std::back_inserter(connect));
+ rval = mb.create_element(MBQUAD, &connect[0], 4, quad); CHECK_ERR(rval);
+
+ ElemEvaluator ee(&mb, quad, 0);
+ ee.set_tag_handle(0, 0);
+ ee.set_eval_set(MBQUAD, LinearQuad::eval_set());
+
+ test_eval(ee, true);
+}
+
void test_linear_hex()
{
Core mb;
diff --git a/test/spatial_locator_test.cpp b/test/spatial_locator_test.cpp
index 31bf08a..6ffa109 100644
--- a/test/spatial_locator_test.cpp
+++ b/test/spatial_locator_test.cpp
@@ -116,7 +116,7 @@ void test_locator(SpatialLocator *sl)
test_pt = box.bMin + CartVect(rx*box_del[0], ry*box_del[1], rz*box_del[2]);
// call spatial locator to locate points
- rval = sl->locate_points(test_pt.array(), 1, &ent, test_res.array(), 0.0, 0.0, &is_in); CHECK_ERR(rval);
+ rval = sl->locate_points(test_pt.array(), 1, &ent, test_res.array(), &is_in); CHECK_ERR(rval);
// verify that the point was found
CHECK_EQUAL(is_in, true);
diff --git a/tools/mbcoupler/DataCoupler.cpp b/tools/mbcoupler/DataCoupler.cpp
new file mode 100644
index 0000000..71f3497
--- /dev/null
+++ b/tools/mbcoupler/DataCoupler.cpp
@@ -0,0 +1,225 @@
+#include "DataCoupler.hpp"
+#include "moab/ParallelComm.hpp"
+#include "moab/Tree.hpp"
+#include "moab/TupleList.hpp"
+#include "moab/SpatialLocator.hpp"
+#include "moab/ElemEvaluator.hpp"
+#include "moab/Error.hpp"
+
+#include "iostream"
+#include <stdio.h>
+#include <algorithm>
+#include <sstream>
+
+#include "assert.h"
+
+namespace moab {
+
+bool debug = false;
+
+DataCoupler::DataCoupler(Interface *impl,
+ ParallelComm *pc,
+ Range &source_ents,
+ int coupler_id,
+ bool init_locator,
+ int dim)
+ : mbImpl(impl), myPcomm(pc), myId(coupler_id), myDim(dim)
+{
+ assert(NULL != mbImpl && (myPcomm || !source_ents.empty()));
+
+ // now initialize the tree
+ if (init_locator) {
+ myLocator = new SpatialLocator(mbImpl, source_ents);
+ myLocator->elem_eval(new ElemEvaluator(mbImpl));
+
+ // initialize element evaluator with the default for the entity types in source_ents;
+ // can be replaced later by application if desired
+ if (!source_ents.empty()) {
+ Range::pair_iterator pit = source_ents.pair_begin();
+ EntityType last_type = MBMAXTYPE;
+ for (; pit != source_ents.pair_end(); pit++) {
+ EntityType this_type = mbImpl->type_from_handle(pit->first);
+ if (last_type == this_type) continue;
+ myLocator->elem_eval()->set_eval_set(pit->first);
+ last_type = this_type;
+ }
+ }
+ }
+
+ if (-1 == dim && !source_ents.empty())
+ dim = mbImpl->dimension_from_handle(*source_ents.rbegin());
+
+ ErrorCode rval = impl->query_interface(mError);
+ if (MB_SUCCESS != rval) throw(rval);
+}
+
+ /* destructor
+ */
+DataCoupler::~DataCoupler()
+{
+ delete myLocator;
+}
+
+ErrorCode DataCoupler::locate_points(Range &targ_ents,
+ double rel_eps,
+ double abs_eps)
+{
+ targetEnts = targ_ents;
+
+ return myLocator->locate_points(targ_ents, rel_eps, abs_eps);
+}
+
+ErrorCode DataCoupler::locate_points(double *xyz, int num_points,
+ double rel_eps,
+ double abs_eps)
+{
+ return myLocator->locate_points(xyz, num_points, rel_eps, abs_eps);
+}
+
+ErrorCode DataCoupler::interpolate(/*DataCoupler::Method*/ int method,
+ const std::string &interp_tag,
+ double *interp_vals,
+ std::vector<int> *point_indices,
+ bool normalize)
+{
+ // tag name input, translate to tag handle and pass down the chain
+
+ // not inlined because of call to set_last_error, class Error isn't in public interface
+ Tag tag;
+ ErrorCode result = mbImpl->tag_get_handle(interp_tag.c_str(), tag);
+ if (MB_SUCCESS != result) {
+ std::ostringstream str;
+ str << "Failed to get handle for interpolation tag \"" << interp_tag << "\"";
+ mError->set_last_error(str.str());
+ return result;
+ }
+ return interpolate(method, tag, interp_vals, point_indices, normalize);
+}
+
+ErrorCode DataCoupler::interpolate(/*DataCoupler::Method*/ int *methods,
+ Tag *tags,
+ int *points_per_method,
+ int num_methods,
+ double *interp_vals,
+ std::vector<int> *point_indices,
+ bool /*normalize*/)
+{
+ // lowest-level interpolate function, does actual interpolation using calls to ElemEvaluator
+
+ ErrorCode result = MB_SUCCESS;
+
+ unsigned int pts_total = 0;
+ for (int i = 0; i < num_methods; i++) pts_total += (points_per_method ? points_per_method[i] : targetEnts.size());
+
+ unsigned int num_indices = (point_indices ? point_indices->size() : targetEnts.size());
+
+ int max_tsize = -1;
+ for (int i = 0; i < num_methods; i++) {
+ int tmp_tsize;
+ result = mbImpl->tag_get_length(tags[i], tmp_tsize);
+ if (MB_SUCCESS != result) return MB_FAILURE;
+ max_tsize = std::max(max_tsize, tmp_tsize);
+ }
+
+ // if tl was passed in non-NULL, just have those points, otherwise have targetPts plus
+ // locally mapped pts
+ if (pts_total != num_indices)
+ return MB_FAILURE;
+
+ if (myPcomm) {
+ // TL to send interpolation indices to target procs
+ // Tuple structure: (pto_i, ridx_i, lidx_i, meth_i, tagidx_i, interp_val[max_tsize]_d)
+ TupleList tinterp;
+ tinterp.initialize(5, 0, 0, max_tsize, num_indices);
+ int t = 0;
+ tinterp.enableWriteAccess();
+ for (int i = 0; i < num_methods; i++) {
+ int num_points = (points_per_method ? points_per_method[i] : targetEnts.size());
+ for (int j = 0; j < num_points; j++) {
+ int idx = (point_indices ? (*point_indices)[j] : j);
+
+ // remote proc/idx from myLocator->parLocTable
+ tinterp.vi_wr[5*t] = myLocator->par_loc_table().vi_rd[2*idx]; // proc
+ tinterp.vi_wr[5*t+1] = myLocator->par_loc_table().vi_rd[2*idx+1]; // remote idx
+
+ // local entity index, tag/method index from my data
+ tinterp.vi_wr[5*t+2] = idx;
+ tinterp.vi_wr[5*t+3] = methods[i];
+ tinterp.vi_wr[5*t+4] = i;
+ tinterp.inc_n();
+ t++;
+ }
+ }
+
+ // scatter/gather interpolation points
+ myPcomm->proc_config().crystal_router()->gs_transfer(1, tinterp, 0);
+
+ // perform interpolation on local source mesh; put results into
+ // tinterp.vr_wr
+
+ for (unsigned int i = 0; i < tinterp.get_n(); i++) {
+ int lidx = tinterp.vi_rd[5*i+1];
+// /*Method*/ int method = (/*Method*/ int)tinterp.vi_rd[5*i+3];
+ Tag tag = tags[tinterp.vi_rd[5*i+4]];
+
+ myLocator->elem_eval()->set_tag_handle(tag);
+ myLocator->elem_eval()->set_ent_handle(myLocator->loc_table().vul_rd[lidx]);
+ result = myLocator->elem_eval()->eval(myLocator->loc_table().vr_rd+3*lidx, tinterp.vr_rd+i*max_tsize);
+ if (MB_SUCCESS != result) return result;
+ }
+
+ // scatter/gather interpolation data
+ myPcomm->proc_config().crystal_router()->gs_transfer(1, tinterp, 0);
+
+ // copy the interpolated field as a unit
+ std::copy(tinterp.vr_rd, tinterp.vr_rd+tinterp.get_n()*max_tsize, interp_vals);
+ }
+ else {
+ std::vector<double> tmp_vals;
+ std::vector<EntityHandle> tmp_ents;
+ double *tmp_dbl = interp_vals;
+ for (int i = 0; i < num_methods; i++) {
+ int num_points = (points_per_method ? points_per_method[i] : targetEnts.size());
+
+ // interpolated data is tsize long, which is either max size (if data passed back to caller in tinterp)
+ // or tag size (if data will be set on entities, in which case it shouldn't have padding)
+ int tsize = max_tsize, tsize_bytes = 0;
+ if (!interp_vals) {
+ tmp_vals.resize(num_points*max_tsize);
+ tmp_dbl = &tmp_vals[0];
+ tmp_ents.resize(num_points);
+ result = mbImpl->tag_get_length(tags[i], tsize);
+ result = mbImpl->tag_get_bytes(tags[i], tsize_bytes);
+ }
+
+ for (int j = 0; j < num_points; j++) {
+ int lidx;
+ if (point_indices) {
+ lidx = (*point_indices)[j];
+ }
+ else {
+ lidx = j;
+ }
+
+ myLocator->elem_eval()->set_tag_handle(tags[i]);
+ myLocator->elem_eval()->set_ent_handle(myLocator->loc_table().vul_rd[lidx]);
+ if (!interp_vals) tmp_ents[j] = targetEnts[lidx]; // could be performance-sensitive, thus the if test
+ result = myLocator->elem_eval()->eval(myLocator->loc_table().vr_rd+3*lidx, tmp_dbl);
+ tmp_dbl += tsize;
+ if (MB_SUCCESS != result) return result;
+ } // for j
+
+ if (!interp_vals) {
+ // set tags on tmp_ents; data is already w/o padding, due to tsize setting above
+ result = mbImpl->tag_set_data(tags[i], &tmp_ents[0], tmp_ents.size(), &tmp_vals[0]);
+ if (MB_SUCCESS != result) return result;
+ }
+
+ } // for i
+ } // if myPcomm
+
+ // done
+ return MB_SUCCESS;
+}
+
+} // namespace_moab
diff --git a/tools/mbcoupler/DataCoupler.hpp b/tools/mbcoupler/DataCoupler.hpp
new file mode 100644
index 0000000..3ad1d74
--- /dev/null
+++ b/tools/mbcoupler/DataCoupler.hpp
@@ -0,0 +1,255 @@
+/**
+ * \class moab::DataCoupler
+ *
+ * \brief This class couples data between meshes.
+ *
+ * The coupler interpolates solution data at a set of points. Data being interpolated resides on a "source"
+ * mesh, in a tag or in vertex coords. Applications calling this coupler send in entities, and receive back
+ * data interpolated at those points. Entities in the source mesh containing those points do not have to reside
+ * on the same processor.
+ *
+ * To use, an application should:
+ * - instantiate this DataCoupler by calling the constructor collectively on all processors in the communicator
+ * - call locate_points, which locates the points to be interpolated and (optionally) caches the results in
+ * this class and SpatialLocator
+ * - call interpolate, which does the interpolation
+ *
+ * Multiple interpolations (of multiple tags, or element-average vs. true interpolation) can be done after
+ * locating the points.
+ *
+ * SpatialLocator is used for the spatial location portion of this work.
+ *
+ * This class is a next-generation implementation of Coupler.
+ */
+#ifndef DATACOUPLER_HPP
+#define DATACOUPLER_HPP
+
+#include "moab/Range.hpp"
+#include "moab/Interface.hpp"
+
+#include <sstream>
+
+namespace moab {
+
+class ParallelComm;
+class SpatialLocator;
+class Error;
+
+class DataCoupler
+{
+public:
+
+ enum IntegType {VOLUME};
+
+ /* constructor
+ * Constructor, which also optionally initializes the coupler
+ * \param pc ParallelComm object to be used with this coupler, representing the union
+ * of processors containing source and target meshes
+ * \param source_elems Elements in the source mesh
+ * \param coupler_id Id of this coupler, should be the same over all procs
+ * \param init_locator If true, initializes a spatial locator inside the constructor
+ * \param dim Dimension of entities to be coupled; if -1, get from source_elems
+ */
+ DataCoupler(Interface *impl,
+ ParallelComm *pc,
+ Range &source_ents,
+ int coupler_id,
+ bool init_locator = true,
+ int dim = -1);
+
+ /* destructor
+ */
+ virtual ~DataCoupler();
+
+ /* \brief Locate points on the source mesh
+ * This is a pass-through function to SpatialLocator::locate_points
+ * \param xyz Point locations (interleaved) being located
+ * \param num_points Number of points in xyz
+ * \param rel_eps Relative tolerance for the non-linear iteration inside a given element
+ * \param abs_eps Absolute tolerance for the non-linear iteration inside a given element
+ * \param loc_results Tuple list containing the results; two types of results are possible,
+ * controlled by value of store_local parameter:
+ * store_local = true: each tuple T[j] consists of (p, i), p = proc with src element containing
+ * point j (or 0 if serial), i = index in stored list (in SpatialLocator) on src proc
+ * store_local = false: each tuple T[j] consists of (p, ht, hs, pr[3]), where ht = target mesh entity
+ * handle, hs = source mesh entity containing point (0 if not found), pr = parameters in hs
+ * \param store_local If true, stores the located points on SpatialLocator
+ */
+ ErrorCode locate_points(double *xyz, int num_points,
+ double rel_eps = 0.0, double abs_eps = 0.0);
+
+ /* \brief Locate points on the source mesh
+ * This is a pass-through function to SpatialLocator::locate_points
+ * \param ents Target entities being located
+ * \param rel_eps Relative tolerance for the non-linear iteration inside a given element
+ * \param abs_eps Absolute tolerance for the non-linear iteration inside a given element
+ * \param loc_results Tuple list containing the results; two types of results are possible,
+ * controlled by value of store_local parameter:
+ * store_local = true: each tuple T[j] consists of (p, i), p = proc with src element containing
+ * point j (or 0 if serial), i = index in stored list (in SpatialLocator) on src proc
+ * store_local = false: each tuple T[j] consists of (p, ht, hs, pr[3]), where ht = target mesh entity
+ * handle, hs = source mesh entity containing point (0 if not found), pr = parameters in hs
+ * \param store_local If true, stores the located points on SpatialLocator
+ */
+ ErrorCode locate_points(Range &ents,
+ double rel_eps = 0.0,
+ double abs_eps = 0.0);
+
+ /* \brief Interpolate data from the source mesh onto points
+ * All entities/points or, if tuple_list is input, only those points
+ * are interpolated from the source mesh. Application should
+ * allocate enough memory in interp_vals to hold interpolation results.
+ *
+ * If normalization is requested, technique used depends on the coupling
+ * method.
+ *
+ * \param method Interpolation/normalization method
+ * \param tag Tag on source mesh holding data to be interpolated
+ * \param interp_vals Memory holding interpolated data; if NULL, data is written to same tag on target ents
+ * \param point_indices If non-NULL, a set of indices of points input to
+ * locate_points at which to interpolate; if NULL, interpolates at all points
+ * input to locate_points
+ * \param normalize If true, normalization is done according to method
+ */
+ ErrorCode interpolate(/*DataCoupler::Method*/ int method,
+ Tag tag,
+ double *interp_vals = NULL,
+ std::vector<int> *point_indices = NULL,
+ bool normalize = true);
+
+ /* \brief Interpolate data from the source mesh onto points
+ * All entities/points or, if tuple_list is input, only those points
+ * are interpolated from the source mesh. Application should
+ * allocate enough memory in interp_vals to hold interpolation results.
+ *
+ * If normalization is requested, technique used depends on the coupling
+ * method.
+ *
+ * \param method Interpolation/normalization method
+ * \param tag_name Tag name on source mesh holding data to be interpolated
+ * \param interp_vals Memory holding interpolated data; if NULL, data is written to same tag on target ents
+ * \param point_indices If non-NULL, a set of indices of points input to
+ * locate_points at which to interpolate; if NULL, interpolates at all points
+ * input to locate_points
+ * \param normalize If true, normalization is done according to method
+ */
+ ErrorCode interpolate(/*DataCoupler::Method*/ int method,
+ const std::string &tag_name,
+ double *interp_vals = NULL,
+ std::vector<int> *point_indices = NULL,
+ bool normalize = true);
+
+ /* \brief Interpolate data from multiple tags
+ * All entities/points or, if tuple_list is input, only those points
+ * are interpolated from the source mesh. Application should
+ * allocate enough memory in interp_vals to hold interpolation results.
+ *
+ * In this variant, multiple tags, possibly with multiple interpolation
+ * methods, are specified. Sum of values in points_per_method should be
+ * the number of points in tl or, if NULL, targetPts.
+ *
+ * If normalization is requested, technique used depends on the coupling
+ * method.
+ *
+ * \param methods Vector of Interpolation/normalization methods
+ * \param tag_names Names of tag being interpolated for each method
+ * \param points_per_method Number of points for each method
+ * \param num_methods Length of vectors in previous 3 arguments
+ * \param interp_vals Memory holding interpolated data; if NULL, data is written to same tag on target ents
+ * \param point_indices If non-NULL, a set of indices of points input to
+ * locate_points at which to interpolate; if NULL, interpolates at all points
+ * input to locate_points
+ * \param normalize If true, normalization is done according to method
+ */
+ ErrorCode interpolate(/*DataCoupler::Method*/ int *methods,
+ const std::string *tag_names,
+ int *points_per_method,
+ int num_methods,
+ double *interp_vals = NULL,
+ std::vector<int> *point_indices = NULL,
+ bool normalize = true);
+
+
+ /* \brief Interpolate data from multiple tags
+ * All entities/points or, if tuple_list is input, only those points
+ * are interpolated from the source mesh. Application should
+ * allocate enough memory in interp_vals to hold interpolation results.
+ *
+ * In this variant, multiple tags, possibly with multiple interpolation
+ * methods, are specified. Sum of values in points_per_method should be
+ * the number of points in tl or, if NULL, targetPts.
+ *
+ * If normalization is requested, technique used depends on the coupling
+ * method.
+ *
+ * \param methods Vector of Interpolation/normalization methods
+ * \param tag_names Names of tag being interpolated for each method
+ * \param points_per_method Number of points for each method
+ * \param num_methods Length of vectors in previous 3 arguments
+ * \param interp_vals Memory holding interpolated data; if NULL, data is written to same tag on target ents
+ * \param point_indices If non-NULL, a set of indices of points input to
+ * locate_points at which to interpolate; if NULL, interpolates at all points
+ * input to locate_points
+ * \param normalize If true, normalization is done according to method
+ */
+ ErrorCode interpolate(/*DataCoupler::Method*/ int *methods,
+ Tag *tag_names,
+ int *points_per_method,
+ int num_methods,
+ double *interp_vals = NULL,
+ std::vector<int> *point_indices = NULL,
+ bool normalize = true);
+
+ /* Get functions */
+ inline SpatialLocator *spatial_locator() {return myLocator;}
+ inline int my_id() const {return myId;}
+ inline const Range &target_ents() const {return targetEnts;}
+ inline Range &target_ents() {return targetEnts;}
+ inline int get_dim() const {return myDim;}
+
+private:
+
+ /* \brief MOAB instance
+ */
+ Interface *mbImpl;
+
+ /* \brief ParallelComm object for this coupler
+ */
+ ParallelComm *myPcomm;
+
+ /* \brief SpatialLocator for local mesh
+ */
+ SpatialLocator *myLocator;
+
+ /* \brief error object used to set last error on interface
+ */
+ Error *mError;
+
+ /* \brief Id of this coupler
+ */
+ int myId;
+
+ /* \brief Range of target entities
+ */
+ Range targetEnts;
+
+ // entity dimension
+ int myDim;
+
+};
+
+ inline ErrorCode DataCoupler::interpolate(/*DataCoupler::Method*/ int method,
+ Tag tag,
+ double *interp_vals,
+ std::vector<int> *point_indices,
+ bool normalize)
+{
+ // no point indices input,
+ int num_pts = (point_indices ? point_indices->size() : targetEnts.size());
+ return interpolate(&method, &tag, &num_pts, 1,
+ interp_vals, point_indices, normalize);
+}
+
+} // namespace moab
+
+#endif
diff --git a/tools/mbcoupler/Makefile.am b/tools/mbcoupler/Makefile.am
index 25a58c2..ea32ff9 100644
--- a/tools/mbcoupler/Makefile.am
+++ b/tools/mbcoupler/Makefile.am
@@ -9,6 +9,7 @@ AM_CPPFLAGS += -DSRCDIR=$(srcdir) \
-I$(top_srcdir)/src \
-I$(top_builddir)/src \
-I$(top_srcdir)/src/parallel \
+ -I$(top_srcdir)/src/LocalDiscretization \
-I$(top_srcdir)/src/moab/point_locater/lotte \
-I$(top_srcdir)/test \
-I$(top_srcdir)/itaps \
@@ -23,10 +24,12 @@ LDADD = libmbcoupler.la $(top_builddir)/src/libMOAB.la $(top_builddir)/itaps/ime
libmbcoupler_la_SOURCES = \
Coupler.cpp \
+ DataCoupler.cpp \
ElemUtil.cpp
libmbcoupler_la_include_HEADERS = \
Coupler.hpp \
+ DataCoupler.hpp \
ElemUtil.hpp
libmbcoupler_la_includedir = $(includedir)
https://bitbucket.org/fathomteam/moab/commits/97620d2ef9c1/
Changeset: 97620d2ef9c1
Branch: deformed-mesh-remap
User: tautges
Date: 2013-10-12 00:49:43
Summary: General debugging and small feature additions.
DeformMeshRemap: moving location in source mesh of slave mesh to before lloyd smoothing.
Adding bounding box calculation for use in relative displacement calculation, using displacements
of about 1% based on communication with LLNL.
SpatialLocator: adding local_num_located and remote_num_located functions, which evaluate contents of locTable and parLocTable for this info.
Also, when eval_set is set, pass it down to the tree too.
elem_eval_test.cpp: adding (commented out for now) test for LinearTri.
Tree, AdaptiveKDTree, BVHTree, various element types: changing tolerance model to using both iteration tolerance, for non-linear
point in element iteration, and inside tolerance, for determining whether parameters are inside.
ElemEvaluator: instead of just asserting on negative jacobian, only do so if we're inside the element.
Also using new tolerance model from above, and incorporating linear tri.
Adding a LinearTri class.
Affected #: 31 files
diff --git a/examples/DeformMeshRemap.cpp b/examples/DeformMeshRemap.cpp
index 500d969..7614934 100644
--- a/examples/DeformMeshRemap.cpp
+++ b/examples/DeformMeshRemap.cpp
@@ -14,6 +14,8 @@
#include "moab/Range.hpp"
#include "moab/LloydSmoother.hpp"
#include "moab/ProgOptions.hpp"
+#include "moab/BoundBox.hpp"
+#include "moab/SpatialLocator.hpp"
#include "MBTagConventions.hpp"
#include "DataCoupler.hpp"
@@ -30,7 +32,7 @@ using namespace std;
ErrorCode read_file(string &fname, EntityHandle &seth,
Range &solids, Range &solid_elems, Range &fluids, Range &fluid_elems);
-void deform_func(double *xold, double *xnew);
+void deform_func(const BoundBox &bbox, double *xold, double *xnew);
ErrorCode deform_master(Range &fluid_elems, Range &solid_elems, Tag &xnew);
ErrorCode smooth_master(int dim, Tag xnew, EntityHandle &master, Range &fluids);
ErrorCode write_to_coords(Range &elems, Tag tagh);
@@ -193,9 +195,30 @@ ErrorCode DeformMeshRemap::execute()
rval = read_file(SLAVE, slaveFileName, slaveSet);
if (MB_SUCCESS != rval) return rval;
+ Range src_elems = solidElems[MASTER];
+ src_elems.merge(fluidElems[MASTER]);
+ // locate slave vertices in master, orig coords; do this with a data coupler, so you can
+ // later interpolate
+ Range tgt_verts, tmp_range = solidElems[SLAVE];
+ tmp_range.merge(fluidElems[SLAVE]);
+ rval = mbImpl->get_adjacencies(tmp_range, 0, false, tgt_verts, Interface::UNION);
+ RR("Failed to get target verts.");
+
+
+ // initialize data coupler on source elements
+ DataCoupler dc_master(mbImpl, NULL, src_elems, 0);
+
+ // locate slave vertices, caching results in dc
+ rval = dc_master.locate_points(tgt_verts); RR("Point location of tgt verts failed.");
+ int num_located = dc_master.spatial_locator()->local_num_located();
+ if (num_located != (int)tgt_verts.size()) {
+ rval = MB_FAILURE;
+ std::cout << "Only " << num_located << " out of " << tgt_verts.size() << " target points successfully located." << std::endl;
+ return rval;
+ }
+
// deform the master's solid mesh, put results in a new tag
rval = deform_master(fluidElems[MASTER], solidElems[MASTER], "xnew"); RR("");
- if (debug) write_and_save(solidElems[MASTER], masterSet, xNew, "deformed.vtk");
{ // to isolate the lloyd smoother & delete when done
@@ -207,22 +230,6 @@ ErrorCode DeformMeshRemap::execute()
}
// map new locations to slave
- // locate slave vertices in master, orig coords; do this with a data coupler, so you can
- // later interpolate
- Range src_elems = solidElems[MASTER];
- src_elems.merge(fluidElems[MASTER]);
-
- // initialize data coupler on source elements
- DataCoupler dc_master(mbImpl, NULL, src_elems, 0);
-
- Range tgt_verts, tmp_range = solidElems[SLAVE];
- tmp_range.merge(fluidElems[SLAVE]);
- rval = mbImpl->get_adjacencies(tmp_range, 0, false, tgt_verts, Interface::UNION);
- RR("Failed to get target verts.");
-
- // locate slave vertices, caching results in dc
- rval = dc_master.locate_points(tgt_verts, 0.0, 1e-10); RR("Point location of tgt verts failed.");
-
// interpolate xNew to slave points
rval = dc_master.interpolate((int)DataCoupler::VOLUME, "xnew"); RR("Failed to interpolate target solution.");
@@ -317,9 +324,10 @@ ErrorCode DeformMeshRemap::write_to_coords(Range &elems, Tag tagh)
return MB_SUCCESS;
}
-void deform_func(double *xold, double *xnew)
+void deform_func(const BoundBox &bbox, double *xold, double *xnew)
{
- const double RODWIDTH = 0.2, RODHEIGHT = 0.5;
+/* Deformation function based on max delx and dely at top of rod
+ const double RODWIDTH = 0.2, RODHEIGHT = 0.5;
// function: origin is at middle base of rod, and is .5 high
// top of rod is (0,.55) on left and (.2,.6) on right
double delx = 0.5*RODWIDTH;
@@ -327,6 +335,13 @@ void deform_func(double *xold, double *xnew)
double xfrac = (xold[0] + .5*RODWIDTH)/RODWIDTH, yfrac = xold[1]/RODHEIGHT;
xnew[0] = xold[0] + yfrac * delx;
xnew[1] = xold[1] + yfrac * (1.0 + xfrac) * 0.05;
+*/
+
+/* Deformation function based on fraction of bounding box dimension in each direction */
+ double frac = 0.01; // taken from approximate relative deformation from LLNL Diablo of XX09 assys
+ CartVect *xo = reinterpret_cast<CartVect*>(xold), *xn = reinterpret_cast<CartVect*>(xnew);
+ CartVect disp = frac * (*xo - bbox.bMin);
+ *xn = *xo + disp;
}
ErrorCode DeformMeshRemap::deform_master(Range &fluid_elems, Range &solid_elems, const char *tag_name)
@@ -346,6 +361,10 @@ ErrorCode DeformMeshRemap::deform_master(Range &fluid_elems, Range &solid_elems,
RR("Failed to get vertex coords.");
rval = mbImpl->tag_set_data(xNew, verts, &coords[0]);
RR("Failed to set xnew tag on fluid verts.");
+
+ // get the bounding box of the solid mesh
+ BoundBox bbox;
+ bbox.update(*mbImpl, solid_elems);
// get all the vertices and coords in the solid
verts.clear();
@@ -356,7 +375,7 @@ ErrorCode DeformMeshRemap::deform_master(Range &fluid_elems, Range &solid_elems,
RR("Failed to get vertex coords.");
unsigned int num_verts = verts.size();
for (unsigned int i = 0; i < num_verts; i++)
- deform_func(&coords[3*i], &coords[3*i]);
+ deform_func(bbox, &coords[3*i], &coords[3*i]);
// set the new tag to those coords
rval = mbImpl->tag_set_data(xNew, verts, &coords[0]);
diff --git a/src/AdaptiveKDTree.cpp b/src/AdaptiveKDTree.cpp
index f66acbc..1324c30 100644
--- a/src/AdaptiveKDTree.cpp
+++ b/src/AdaptiveKDTree.cpp
@@ -1303,7 +1303,8 @@ namespace moab {
ErrorCode AdaptiveKDTree::point_search(const double *point,
EntityHandle& leaf_out,
- double tol,
+ const double iter_tol,
+ const double inside_tol,
bool *multiple_leaves,
EntityHandle *start_node,
CartVect *params)
@@ -1322,7 +1323,7 @@ namespace moab {
treeStats.nodesVisited++;
ErrorCode rval = get_bounding_box(box, &node);
if (MB_SUCCESS != rval) return rval;
- if (!box.contains_point(point, tol)) return MB_SUCCESS;
+ if (!box.contains_point(point, iter_tol)) return MB_SUCCESS;
rval = moab()->get_child_meshsets( node, children );
if (MB_SUCCESS != rval)
@@ -1346,7 +1347,7 @@ namespace moab {
treeStats.leavesVisited++;
if (myEval && params) {
- rval = myEval->find_containing_entity(node, point, tol,
+ rval = myEval->find_containing_entity(node, point, iter_tol, inside_tol,
leaf_out, params->array(), &treeStats.leafObjectTests);
if (MB_SUCCESS != rval) return rval;
}
@@ -1358,7 +1359,8 @@ namespace moab {
ErrorCode AdaptiveKDTree::point_search(const double *point,
AdaptiveKDTreeIter& leaf_it,
- double tol,
+ const double iter_tol,
+ const double /*inside_tol*/,
bool *multiple_leaves,
EntityHandle *start_node)
{
@@ -1372,7 +1374,7 @@ namespace moab {
leaf_it.mBox[1] = boundBox.bMax;
// test that point is inside tree
- if (!boundBox.contains_point(point, tol)) {
+ if (!boundBox.contains_point(point, iter_tol)) {
treeStats.nodesVisited++;
return MB_ENTITY_NOT_FOUND;
}
@@ -1423,7 +1425,8 @@ namespace moab {
ErrorCode AdaptiveKDTree::distance_search(const double from_point[3],
const double distance,
std::vector<EntityHandle>& result_list,
- double tol,
+ const double iter_tol,
+ const double inside_tol,
std::vector<double> *result_dists,
std::vector<CartVect> *result_params,
EntityHandle *tree_root)
@@ -1445,7 +1448,7 @@ namespace moab {
// (zero if inside box)
BoundBox box;
rval = get_bounding_box(box);
- if (MB_SUCCESS == rval && !box.contains_point(from_point, tol)) {
+ if (MB_SUCCESS == rval && !box.contains_point(from_point, iter_tol)) {
treeStats.nodesVisited++;
return MB_SUCCESS;
}
@@ -1484,7 +1487,7 @@ namespace moab {
if (myEval && result_params) {
EntityHandle ent;
CartVect params;
- rval = myEval->find_containing_entity(node.handle, from_point, tol,
+ rval = myEval->find_containing_entity(node.handle, from_point, iter_tol, inside_tol,
ent, params.array(), &treeStats.leafObjectTests);
if (MB_SUCCESS != rval) return rval;
else if (ent) {
@@ -1802,7 +1805,7 @@ namespace moab {
// the same distance from the input point as the current closest
// point is.
CartVect diff = closest_pt - from;
- rval = distance_search(from_coords, sqrt(diff%diff), leaves, 0.0, NULL, NULL, &tree_root);
+ rval = distance_search(from_coords, sqrt(diff%diff), leaves, 1.0e-10, 1.0e-6, NULL, NULL, &tree_root);
if (MB_SUCCESS != rval) return rval;
// Check any close leaves to see if they contain triangles that
@@ -1834,7 +1837,7 @@ namespace moab {
// get leaves of tree that intersect sphere
assert(tree_root);
- rval = distance_search(center, radius, leaves, 0.0, NULL, NULL, &tree_root);
+ rval = distance_search(center, radius, leaves, 1.0e-10, 1.0e-6, NULL, NULL, &tree_root);
if (MB_SUCCESS != rval) return rval;
// search each leaf for triangles intersecting sphere
diff --git a/src/BVHTree.cpp b/src/BVHTree.cpp
index 8903070..b65eae4 100644
--- a/src/BVHTree.cpp
+++ b/src/BVHTree.cpp
@@ -435,7 +435,8 @@ namespace moab
ErrorCode BVHTree::find_point(const std::vector<double> &point,
const unsigned int &index,
- const double tol,
+ const double iter_tol,
+ const double inside_tol,
std::pair<EntityHandle, CartVect> &result)
{
if (index == 0) treeStats.numTraversals++;
@@ -453,7 +454,7 @@ namespace moab
for(Range::iterator i = entities.begin(); i != entities.end(); i++) {
treeStats.leafObjectTests++;
myEval->set_ent_handle(*i);
- myEval->reverse_eval(&point[0], tol, params.array(), &is_inside);
+ myEval->reverse_eval(&point[0], iter_tol, inside_tol, params.array(), &is_inside);
if (is_inside) {
result.first = *i;
result.second = params;
@@ -469,11 +470,11 @@ namespace moab
rval = mbImpl->get_child_meshsets(startSetHandle+index, children);
if (MB_SUCCESS != rval || children.size() != 2) return rval;
- if((node.Lmax+tol) < (node.Rmin-tol)){
- if(point[node.dim] <= (node.Lmax + tol))
- return find_point(point, children[0]-startSetHandle, tol, result);
- else if(point[ node.dim] >= (node.Rmin - tol))
- return find_point(point, children[1]-startSetHandle, tol, result);
+ if((node.Lmax+iter_tol) < (node.Rmin-iter_tol)){
+ if(point[node.dim] <= (node.Lmax + iter_tol))
+ return find_point(point, children[0]-startSetHandle, iter_tol, inside_tol, result);
+ else if(point[ node.dim] >= (node.Rmin - iter_tol))
+ return find_point(point, children[1]-startSetHandle, iter_tol, inside_tol, result);
result = std::make_pair(0, CartVect(&point[0])); //point lies in empty space.
return MB_SUCCESS;
}
@@ -482,11 +483,11 @@ namespace moab
//left of Rmin, you must be on the left
//we can't be sure about the boundaries since the boxes overlap
//this was a typo in the paper which caused pain.
- if(point[node.dim] < (node.Rmin - tol))
- return find_point(point, children[0]-startSetHandle, tol, result);
+ if(point[node.dim] < (node.Rmin - iter_tol))
+ return find_point(point, children[0]-startSetHandle, iter_tol, inside_tol, result);
//if you are on the right Lmax, you must be on the right
- else if(point[ node.dim] > (node.Lmax+tol))
- return find_point(point, children[1]-startSetHandle, tol, result);
+ else if(point[ node.dim] > (node.Lmax+iter_tol))
+ return find_point(point, children[1]-startSetHandle, iter_tol, inside_tol, result);
/* pg5 of paper
* However, instead of always traversing either subtree
@@ -501,22 +502,22 @@ namespace moab
//branch predicted..
//bool dir = (point[ node.dim] - node.Rmin) <=
// (node.Lmax - point[ node.dim]);
- find_point(point, children[0]-startSetHandle, tol, result);
+ find_point(point, children[0]-startSetHandle, iter_tol, inside_tol, result);
if(result.first == 0){
- return find_point(point, children[1]-startSetHandle, tol, result);
+ return find_point(point, children[1]-startSetHandle, iter_tol, inside_tol, result);
}
return MB_SUCCESS;
}
- EntityHandle BVHTree::bruteforce_find(const double *point, const double tol) {
+ EntityHandle BVHTree::bruteforce_find(const double *point, const double iter_tol, const double inside_tol) {
treeStats.numTraversals++;
CartVect params;
for(unsigned int i = 0; i < myTree.size(); i++) {
- if(myTree[i].dim != 3 || !myTree[i].box.contains_point(point, tol)) continue;
+ if(myTree[i].dim != 3 || !myTree[i].box.contains_point(point, iter_tol)) continue;
if (myEval) {
EntityHandle entity = 0;
treeStats.leavesVisited++;
- ErrorCode rval = myEval->find_containing_entity(startSetHandle+i, point, tol,
+ ErrorCode rval = myEval->find_containing_entity(startSetHandle+i, point, iter_tol, inside_tol,
entity, params.array(), &treeStats.leafObjectTests);
if (entity) return entity;
else if (MB_SUCCESS != rval) return 0;
@@ -542,7 +543,8 @@ namespace moab
ErrorCode BVHTree::point_search(const double *point,
EntityHandle& leaf_out,
- double tol,
+ const double iter_tol,
+ const double inside_tol,
bool *multiple_leaves,
EntityHandle *start_node,
CartVect *params)
@@ -572,7 +574,7 @@ namespace moab
// test box of this node
ErrorCode rval = get_bounding_box(box, &this_set);
if (MB_SUCCESS != rval) return rval;
- if (!box.contains_point(point, tol)) continue;
+ if (!box.contains_point(point, iter_tol)) continue;
// else if not a leaf, test children & put on list
else if (myTree[ind].dim != 3) {
@@ -581,7 +583,7 @@ namespace moab
continue;
}
else if (myTree[ind].dim == 3 && myEval && params) {
- rval = myEval->find_containing_entity(startSetHandle+ind, point, tol,
+ rval = myEval->find_containing_entity(startSetHandle+ind, point, iter_tol, inside_tol,
leaf_out, params->array(), &treeStats.leafObjectTests);
if (leaf_out || MB_SUCCESS != rval) return rval;
}
@@ -599,7 +601,8 @@ namespace moab
ErrorCode BVHTree::distance_search(const double from_point[3],
const double distance,
std::vector<EntityHandle>& result_list,
- double params_tol,
+ const double iter_tol,
+ const double inside_tol,
std::vector<double> *result_dists,
std::vector<CartVect> *result_params,
EntityHandle *tree_root)
@@ -653,7 +656,7 @@ namespace moab
if (myEval && result_params) {
EntityHandle ent;
CartVect params;
- rval = myEval->find_containing_entity(startSetHandle+ind, from_point, params_tol,
+ rval = myEval->find_containing_entity(startSetHandle+ind, from_point, iter_tol, inside_tol,
ent, params.array(), &treeStats.leafObjectTests);
if (MB_SUCCESS != rval) return rval;
else if (ent) {
diff --git a/src/LocalDiscretization/ElemEvaluator.cpp b/src/LocalDiscretization/ElemEvaluator.cpp
index 4446281..4ca3943 100644
--- a/src/LocalDiscretization/ElemEvaluator.cpp
+++ b/src/LocalDiscretization/ElemEvaluator.cpp
@@ -6,6 +6,7 @@
// need to include eval set types here to support get_eval_set; alternative would be to have some
// type of registration, but we'd still need static registration for the built-in types
+#include "moab/LinearTri.hpp"
#include "moab/LinearQuad.hpp"
#include "moab/LinearTet.hpp"
#include "moab/LinearHex.hpp"
@@ -16,12 +17,12 @@
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 tol, double *work,
- double *params, bool *inside) {
+ const int ndim, const double iter_tol, const double inside_tol,
+ double *work, double *params, bool *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
- const double error_tol_sqr = tol*tol;
+ const double error_tol_sqr = iter_tol*iter_tol;
CartVect *cvparams = reinterpret_cast<CartVect*>(params);
const CartVect *cvposn = reinterpret_cast<const CartVect*>(posn);
@@ -38,23 +39,26 @@ 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 iters=0;
// while |res| larger than tol
while (res % res > error_tol_sqr) {
if(++iters>10) {
- if (inside) {
- // if we haven't converged but we're outside, that's defined as success
- *inside = (*inside_f)(params, ndim, tol);
- if (!(*inside)) return MB_SUCCESS;
- }
- return MB_FAILURE;
+ // if we haven't converged but we're outside, that's defined as success
+ *tmp_inside = (*inside_f)(params, ndim, inside_tol);
+ if (!(*tmp_inside)) return MB_SUCCESS;
+ else return MB_FAILURE;
}
// get jacobian at current params
rval = (*jacob)(cvparams->array(), verts, nverts, ndim, work, J[0]);
double det = J.determinant();
- assert(det > std::numeric_limits<double>::epsilon());
+ if (det < std::numeric_limits<double>::epsilon()) {
+ *tmp_inside = (*inside_f)(params, ndim, inside_tol);
+ if (!(*tmp_inside)) return MB_SUCCESS;
+ else return MB_FAILURE;
+ }
// new params tries to eliminate residual
*cvparams -= J.inverse(1.0/det) * res;
@@ -68,7 +72,7 @@ namespace moab {
}
if (inside)
- *inside = (*inside_f)(params, ndim, tol);
+ *inside = (*inside_f)(params, ndim, inside_tol);
return MB_SUCCESS;
}// Map::evaluate_reverse()
@@ -89,6 +93,7 @@ namespace moab {
case MBEDGE:
break;
case MBTRI:
+ if (LinearTri::compatible(tp, num_vertices, eval_set)) return MB_SUCCESS;
break;
case MBQUAD:
if (LinearQuad::compatible(tp, num_vertices, eval_set)) return MB_SUCCESS;
@@ -109,9 +114,9 @@ namespace moab {
return MB_NOT_IMPLEMENTED;
}
- ErrorCode ElemEvaluator::find_containing_entity(Range &entities, const double *point, double tol,
- EntityHandle &containing_ent, double *params,
- unsigned int *num_evals)
+ ErrorCode ElemEvaluator::find_containing_entity(Range &entities, const double *point, const double iter_tol,
+ const double inside_tol, EntityHandle &containing_ent,
+ double *params, unsigned int *num_evals)
{
bool is_inside;
ErrorCode rval = MB_SUCCESS;
@@ -120,7 +125,7 @@ namespace moab {
for(i = entities.begin(); i != entities.end(); i++) {
nevals++;
set_ent_handle(*i);
- rval = reverse_eval(point, tol, params, &is_inside);
+ rval = reverse_eval(point, iter_tol, inside_tol, params, &is_inside);
if (MB_SUCCESS != rval) return rval;
if (is_inside) break;
}
diff --git a/src/LocalDiscretization/LinearHex.cpp b/src/LocalDiscretization/LinearHex.cpp
index 14de2d9..5e730af 100644
--- a/src/LocalDiscretization/LinearHex.cpp
+++ b/src/LocalDiscretization/LinearHex.cpp
@@ -96,10 +96,12 @@ 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 tol, double *work, double *params, bool *is_inside)
+ const double iter_tol, const double inside_tol, double *work,
+ double *params, bool *is_inside)
{
assert(posn && verts);
- return EvalSet::evaluate_reverse(eval, jacob, ins, posn, verts, nverts, ndim, tol, work, params, is_inside);
+ 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)
diff --git a/src/LocalDiscretization/LinearQuad.cpp b/src/LocalDiscretization/LinearQuad.cpp
index 93237f9..63c9959 100644
--- a/src/LocalDiscretization/LinearQuad.cpp
+++ b/src/LocalDiscretization/LinearQuad.cpp
@@ -75,9 +75,11 @@ 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 tol, double *work, double *params, bool *is_inside)
+ const double iter_tol, const double inside_tol, double *work,
+ double *params, bool *is_inside)
{
- return EvalSet::evaluate_reverse(eval, jacob, ins, posn, verts, nverts, ndim, tol, work, params, 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)
diff --git a/src/LocalDiscretization/LinearTet.cpp b/src/LocalDiscretization/LinearTet.cpp
index 705dbee..4965cfc 100644
--- a/src/LocalDiscretization/LinearTet.cpp
+++ b/src/LocalDiscretization/LinearTet.cpp
@@ -74,27 +74,29 @@ 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 tol, double *work, double *params, bool *is_inside)
+ const double iter_tol, const double inside_tol, double *work,
+ double *params, bool *is_inside)
{
assert(posn && verts);
- return evaluate_reverse(eval, jacob, ins, posn, verts, nverts, ndim, tol, work, params, is_inside);
+ 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)
{
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);
+ params[0] + params[1] + params[2] <= 1.0+tol);
}
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 tol, double *work,
- double *params, bool *inside) {
+ const int ndim, const double iter_tol, const double inside_tol,
+ double *work, double *params, bool *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
- const double error_tol_sqr = tol*tol;
+ const double error_tol_sqr = iter_tol*iter_tol;
CartVect *cvparams = reinterpret_cast<CartVect*>(params);
const CartVect *cvposn = reinterpret_cast<const CartVect*>(posn);
@@ -138,7 +140,7 @@ namespace moab
}
if (inside)
- *inside = (*inside_f)(params, ndim, tol);
+ *inside = (*inside_f)(params, ndim, inside_tol);
return MB_SUCCESS;
}// Map::evaluate_reverse()
diff --git a/src/LocalDiscretization/LinearTri.cpp b/src/LocalDiscretization/LinearTri.cpp
new file mode 100644
index 0000000..8fe9d23
--- /dev/null
+++ b/src/LocalDiscretization/LinearTri.cpp
@@ -0,0 +1,145 @@
+#include "moab/LinearTri.hpp"
+#include "moab/Forward.hpp"
+#include <algorithm>
+
+namespace moab
+{
+
+ const double LinearTri::corner[3][2] = { {0,0},
+ {1,0},
+ {0,1}};
+
+ ErrorCode LinearTri::initFcn(const double *verts, const int /*nverts*/, double *&work) {
+ // allocate work array as:
+ // work[0..8] = T
+ // work[9..17] = Tinv
+ // work[18] = detT
+ // work[19] = detTinv
+ assert(!work && verts);
+ work = new double[20];
+ Matrix3 *T = reinterpret_cast<Matrix3*>(work),
+ *Tinv = reinterpret_cast<Matrix3*>(work+9);
+ double *detT = work+18, *detTinv = work+19;
+
+ *T = Matrix3(verts[1*3+0]-verts[0*3+0],verts[2*3+0]-verts[0*3+0],0.0,
+ verts[1*3+1]-verts[0*3+1],verts[2*3+1]-verts[0*3+1],0.0,
+ verts[1*3+2]-verts[0*3+2],verts[2*3+2]-verts[0*3+2],1.0);
+ *T *= 0.5;
+ (*T)(2,2) = 1.0;
+
+ *Tinv = T->inverse();
+ *detT = T->determinant();
+ *detTinv = (0.0 == *detT ? HUGE : 1.0 / *detT);
+
+ return MB_SUCCESS;
+ }
+
+ ErrorCode LinearTri::evalFcn(const double *params, const double *field, const int /*ndim*/, const int num_tuples,
+ double */*work*/, double *result) {
+ assert(params && field && num_tuples > 0);
+ // convert to [0,1]
+ double p1 = 0.5 * (1.0 + params[0]),
+ p2 = 0.5 * (1.0 + params[1]),
+ p0 = 1.0 - p1 - p2;
+
+ for (int j = 0; j < num_tuples; j++)
+ result[j] = p0 * field[0*num_tuples+j] + p1 * field[1*num_tuples+j] + p2 * field[2*num_tuples+j];
+
+ return MB_SUCCESS;
+ }
+
+ ErrorCode LinearTri::integrateFcn(const double *field, const double */*verts*/, const int /*nverts*/, const int /*ndim*/, const int num_tuples,
+ double *work, double *result)
+ {
+ assert(field && num_tuples > 0);
+ double tmp = work[18];
+
+ for (int i = 0; i < num_tuples; i++)
+ result[i] = tmp * (field[num_tuples+i] + field[2*num_tuples+i]);
+
+ return MB_SUCCESS;
+ }
+
+ ErrorCode LinearTri::jacobianFcn(const double *, const double *, const int, const int ,
+ double *work, double *result)
+ {
+ // jacobian is cached in work array
+ assert(work);
+ std::copy(work, work+9, result);
+ return MB_SUCCESS;
+ }
+
+ 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)
+ {
+ 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)
+ {
+ return (params[0] >= -1.0-tol && params[1] >= -1.0-tol &&
+ params[0] + params[1] <= 1.0+tol);
+
+ }
+
+ 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) {
+ // 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
+ const double error_tol_sqr = iter_tol*iter_tol;
+ CartVect *cvparams = reinterpret_cast<CartVect*>(params);
+ const CartVect *cvposn = reinterpret_cast<const CartVect*>(posn);
+
+ // find best initial guess to improve convergence
+ CartVect tmp_params[] = {CartVect(-1,-1,-1), CartVect(1,-1,-1), CartVect(-1,1,-1)};
+ double resl = HUGE;
+ CartVect new_pos, tmp_pos;
+ ErrorCode rval;
+ for (unsigned int i = 0; i < 3; i++) {
+ rval = (*eval)(tmp_params[i].array(), verts, ndim, 3, work, tmp_pos.array());
+ if (MB_SUCCESS != rval) return rval;
+ double tmp_resl = (tmp_pos-*cvposn).length_squared();
+ if (tmp_resl < resl) {
+ *cvparams = tmp_params[i];
+ new_pos = tmp_pos;
+ resl = tmp_resl;
+ }
+ }
+
+ // residual is diff between old and new pos; need to minimize that
+ CartVect res = new_pos - *cvposn;
+ Matrix3 J;
+ rval = (*jacob)(cvparams->array(), verts, nverts, ndim, work, J[0]);
+ double det = J.determinant();
+ assert(det > std::numeric_limits<double>::epsilon());
+ Matrix3 Ji = J.inverse(1.0/det);
+
+ int iters=0;
+ // while |res| larger than tol
+ while (res % res > error_tol_sqr) {
+ if(++iters>25)
+ return MB_FAILURE;
+
+ // new params tries to eliminate residual
+ *cvparams -= Ji * res;
+
+ // get the new forward-evaluated position, and its difference from the target pt
+ rval = (*eval)(params, verts, ndim, 3, work, new_pos.array());
+ if (MB_SUCCESS != rval) return rval;
+ res = new_pos - *cvposn;
+ }
+
+ if (inside)
+ *inside = (*inside_f)(params, ndim, inside_tol);
+
+ return MB_SUCCESS;
+ }// Map::evaluate_reverse()
+
+} // namespace moab
diff --git a/src/LocalDiscretization/Makefile.am b/src/LocalDiscretization/Makefile.am
index 7309aef..cb34071 100644
--- a/src/LocalDiscretization/Makefile.am
+++ b/src/LocalDiscretization/Makefile.am
@@ -17,6 +17,7 @@ MOAB_LOCALDISCR_SRCS = \
LinearHex.cpp \
LinearQuad.cpp \
LinearTet.cpp \
+ LinearTri.cpp \
QuadraticHex.cpp
MOAB_LOCALDISCR_HDRS = \
@@ -24,6 +25,7 @@ MOAB_LOCALDISCR_HDRS = \
moab/LinearHex.hpp \
moab/LinearQuad.hpp \
moab/LinearTet.hpp \
+ moab/LinearTri.hpp \
moab/QuadraticHex.hpp
# The list of source files, and any header files that do not need to be installed
diff --git a/src/LocalDiscretization/QuadraticHex.cpp b/src/LocalDiscretization/QuadraticHex.cpp
index ca4758b..5c621f6 100644
--- a/src/LocalDiscretization/QuadraticHex.cpp
+++ b/src/LocalDiscretization/QuadraticHex.cpp
@@ -107,10 +107,12 @@ 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 tol, double *work, double *params, bool *is_inside)
+ const double iter_tol, const double inside_tol, double *work,
+ double *params, bool *is_inside)
{
assert(posn && verts);
- return EvalSet::evaluate_reverse(eval, jacob, ins, posn, verts, nverts, ndim, tol, work, params, is_inside);
+ 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)
diff --git a/src/LocalDiscretization/SpectralHex.cpp b/src/LocalDiscretization/SpectralHex.cpp
index 0e40eba..31f5391 100644
--- a/src/LocalDiscretization/SpectralHex.cpp
+++ b/src/LocalDiscretization/SpectralHex.cpp
@@ -73,7 +73,8 @@ CartVect SpectralHex::evaluate( const CartVect& params ) const
return result;
}
// replicate the functionality of hex_findpt
-bool SpectralHex::evaluate_reverse(CartVect const & xyz, CartVect ¶ms, double tol, const CartVect &init) const
+ bool SpectralHex::evaluate_reverse(CartVect const & xyz, CartVect ¶ms, double iter_tol, const double inside_tol,
+ const CartVect &init) const
{
params = init;
@@ -92,7 +93,7 @@ bool SpectralHex::evaluate_reverse(CartVect const & xyz, CartVect ¶ms, doubl
//copy parametric coords back
params = r;
- return is_inside(params, tol);
+ return is_inside(params, inside_tol);
}
Matrix3 SpectralHex::jacobian(const CartVect& params) const
{
diff --git a/src/LocalDiscretization/SpectralQuad.cpp b/src/LocalDiscretization/SpectralQuad.cpp
index cbd1251..5c11d13 100644
--- a/src/LocalDiscretization/SpectralQuad.cpp
+++ b/src/LocalDiscretization/SpectralQuad.cpp
@@ -95,7 +95,8 @@ 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 tol, double *work, double *params, bool *is_inside)
+ const double iter_tol, const double inside_tol, double *work,
+ double *params, bool *is_inside)
{
params = init;
@@ -114,7 +115,7 @@ bool SpectralQuad::reverseEvalFcn(const double *posn, const double *verts, const
//copy parametric coords back
params = r;
- return insideFcn(params, 2, tol);
+ return insideFcn(params, 2, inside_tol);
}
diff --git a/src/LocalDiscretization/moab/ElemEvaluator.hpp b/src/LocalDiscretization/moab/ElemEvaluator.hpp
index eac008d..29d1051 100644
--- a/src/LocalDiscretization/moab/ElemEvaluator.hpp
+++ b/src/LocalDiscretization/moab/ElemEvaluator.hpp
@@ -25,7 +25,8 @@ namespace moab {
typedef ErrorCode (*ReverseEvalFcn)(EvalFcn eval, JacobianFcn jacob, InsideFcn ins,
const double *posn, const double *verts, const int nverts, const int ndim,
- const double tol, double *work, double *params, bool *is_inside);
+ const double iter_tol, const double inside_tol,
+ double *work, double *params, bool *is_inside);
class EvalSet
{
@@ -76,8 +77,8 @@ namespace moab {
/** \brief Common function to do reverse evaluation based on evaluation and jacobian functions */
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 tol, double *work, double *params,
- bool *inside);
+ const int ndim, const double iter_tol, const double inside_tol,
+ double *work, double *params, bool *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);
};
@@ -133,12 +134,14 @@ namespace moab {
/** \brief Reverse-evaluate the cached entity at a given physical position
* \param posn Position at which to evaluate parameters
- * \param tol Tolerance of reverse evaluation, usually 10^-6 or so
+ * \param iter_tol Tolerance of reverse evaluation non-linear iteration, usually 10^-10 or so
+ * \param inside_tol Tolerance of is_inside evaluation, usually 10^-6 or so
* \param params Result of evaluation
* \param is_inside If non-NULL, returns true of resulting parameters place the point inside the element
* (in most cases, within [-1]*(dim)
*/
- ErrorCode reverse_eval(const double *posn, double tol, double *params, bool *is_inside = NULL) const;
+ ErrorCode reverse_eval(const double *posn, double iter_tol, double inside_tol, double *params,
+ bool *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
@@ -166,14 +169,16 @@ namespace moab {
* object is changed.
* \param entities Entities tested
* \param point Point tested, must have 3 dimensions, even for edge and face entities
- * \param tol Tolerance for is_inside test
+ * \param iter_tol Tolerance for non-linear reverse evaluation
+ * \param inside_tol Tolerance for is_inside test
* \param containing_ent Entity containing the point, returned 0 if no entity
* \param params Parameters of point in containing entity, unchanged if no containing entity
* \param num_evals If non-NULL, incremented each time reverse_eval is called
* \return Returns non-success only if evaulation failed for some reason (point not in element is NOT a
* reason for failure)
*/
- ErrorCode find_containing_entity(Range &entities, const double *point, double tol,
+ ErrorCode find_containing_entity(Range &entities, const double *point,
+ const double iter_tol, const double inside_tol,
EntityHandle &containing_ent, double *params,
unsigned int *num_evals = NULL);
@@ -186,14 +191,16 @@ namespace moab {
* object is changed.
* \param ent_set Entity set containing the entities to be tested
* \param point Point tested, must have 3 dimensions, even for edge and face entities
- * \param tol Tolerance for is_inside test
+ * \param iter_tol Tolerance for non-linear reverse evaluation
+ * \param inside_tol Tolerance for is_inside test
* \param containing_ent Entity containing the point, returned 0 if no entity
* \param params Parameters of point in containing entity, unchanged if no containing entity
* \param num_evals If non-NULL, incremented each time reverse_eval is called
* \return Returns non-success only if evaulation failed for some reason (point not in element is NOT a
* reason for failure)
*/
- ErrorCode find_containing_entity(EntityHandle ent_set, const double *point, double tol,
+ ErrorCode find_containing_entity(EntityHandle ent_set, const double *point,
+ const double iter_tol, const double inside_tol,
EntityHandle &containing_ent, double *params,
unsigned int *num_evals = NULL);
@@ -218,8 +225,12 @@ namespace moab {
inline ErrorCode set_ent_handle(EntityHandle ent);
/** \brief Get entity handle for this ElemEval */
- inline EntityHandle get_ent_handle() const {return entHandle;};
+ inline EntityHandle get_ent_handle() const {return entHandle;}
+ /* \brief Get vertex positions cached on this EE
+ */
+ inline double *get_vert_pos() {return vertPos[0].array();}
+
/* \brief Get the vertex handles cached here */
inline const EntityHandle *get_vert_handles() const {return vertHandles;}
@@ -251,6 +262,11 @@ namespace moab {
/* \brief Set the dimension of entities having the tag */
inline ErrorCode set_tagged_ent_dim(int dim);
+ /* \brief Get work space, sometimes this is useful for evaluating data you don't want to set as tags on entities
+ * Can't be const because most of the functions (evaluate, integrate, etc.) take non-const work space *'s
+ */
+ inline double *get_work_space() {return workSpace;}
+
/* \brief MOAB interface cached on this evaluator */
inline Interface *get_moab() {return mbImpl;}
@@ -434,12 +450,13 @@ namespace moab {
(-1 == num_tuples ? numTuples : num_tuples), workSpace, result);
}
- inline ErrorCode ElemEvaluator::reverse_eval(const double *posn, const double tol, double *params, bool *ins) const
+ inline ErrorCode ElemEvaluator::reverse_eval(const double *posn, const double iter_tol, const double inside_tol,
+ double *params, bool *ins) const
{
assert(entHandle && MBMAXTYPE != entType);
return (*evalSets[entType].reverseEvalFcn)(evalSets[entType].evalFcn, evalSets[entType].jacobianFcn, evalSets[entType].insideFcn,
posn, vertPos[0].array(), numVerts,
- entDim, tol, workSpace, params, ins);
+ entDim, iter_tol, inside_tol, workSpace, params, ins);
}
/** \brief Evaluate the jacobian of the cached entity at a given parametric location */
@@ -464,7 +481,8 @@ namespace moab {
workSpace, result);
}
- inline ErrorCode ElemEvaluator::find_containing_entity(EntityHandle ent_set, const double *point, double tol,
+ inline ErrorCode ElemEvaluator::find_containing_entity(EntityHandle ent_set, const double *point,
+ const double iter_tol, const double inside_tol,
EntityHandle &containing_ent, double *params,
unsigned int *num_evals)
{
@@ -472,7 +490,7 @@ namespace moab {
Range entities;
ErrorCode rval = mbImpl->get_entities_by_handle(ent_set, entities);
if (MB_SUCCESS != rval) return rval;
- else return find_containing_entity(entities, point, tol, containing_ent, params, num_evals);
+ 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
diff --git a/src/LocalDiscretization/moab/LinearHex.hpp b/src/LocalDiscretization/moab/LinearHex.hpp
index e7ca3f4..03ed33f 100644
--- a/src/LocalDiscretization/moab/LinearHex.hpp
+++ b/src/LocalDiscretization/moab/LinearHex.hpp
@@ -17,7 +17,8 @@ public:
/** \brief Reverse-evaluation of parametric coordinates at physical space position */
static ErrorCode reverseEvalFcn(EvalFcn eval, JacobianFcn jacob, InsideFcn ins,
const double *posn, const double *verts, const int nverts, const int ndim,
- const double tol, double *work, double *params, bool *is_inside);
+ const double iter_tol, const double inside_tol, double *work,
+ double *params, bool *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,
diff --git a/src/LocalDiscretization/moab/LinearQuad.hpp b/src/LocalDiscretization/moab/LinearQuad.hpp
index 26882a1..2eb3027 100644
--- a/src/LocalDiscretization/moab/LinearQuad.hpp
+++ b/src/LocalDiscretization/moab/LinearQuad.hpp
@@ -17,7 +17,8 @@ public:
/** \brief Reverse-evaluation of parametric coordinates at physical space position */
static ErrorCode reverseEvalFcn(EvalFcn eval, JacobianFcn jacob, InsideFcn ins,
const double *posn, const double *verts, const int nverts, const int ndim,
- const double tol, double *work, double *params, bool *is_inside);
+ const double iter_tol, const double inside_tol, double *work,
+ double *params, bool *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,
diff --git a/src/LocalDiscretization/moab/LinearTet.hpp b/src/LocalDiscretization/moab/LinearTet.hpp
index 464b1f6..f621b2f 100644
--- a/src/LocalDiscretization/moab/LinearTet.hpp
+++ b/src/LocalDiscretization/moab/LinearTet.hpp
@@ -17,7 +17,8 @@ public:
/** \brief Reverse-evaluation of parametric coordinates at physical space position */
static ErrorCode reverseEvalFcn(EvalFcn eval, JacobianFcn jacob, InsideFcn ins,
const double *posn, const double *verts, const int nverts, const int ndim,
- const double tol, double *work, double *params, bool *is_inside);
+ const double iter_tol, const double inside_tol, double *work,
+ double *params, bool *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,
@@ -35,7 +36,7 @@ public:
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 tol, double *work,
+ const int ndim, const double iter_tol, const double inside_tol, double *work,
double *params, bool *inside);
static EvalSet eval_set()
@@ -45,7 +46,7 @@ public:
static bool compatible(EntityType tp, int numv, EvalSet &eset)
{
- if (tp == MBTET && numv == 4) {
+ if (tp == MBTET && numv >= 4) {
eset = eval_set();
return true;
}
diff --git a/src/LocalDiscretization/moab/LinearTri.hpp b/src/LocalDiscretization/moab/LinearTri.hpp
new file mode 100644
index 0000000..93820a7
--- /dev/null
+++ b/src/LocalDiscretization/moab/LinearTri.hpp
@@ -0,0 +1,63 @@
+#ifndef LINEAR_TRI_HPP
+#define LINEAR_TRI_HPP
+ /**\brief Shape function space for trilinear tetrahedron, obtained by a pushforward of the canonical linear (affine) functions. */
+
+#include "moab/ElemEvaluator.hpp"
+
+namespace moab
+{
+
+class LinearTri
+{
+public:
+ /** \brief Forward-evaluation of field at parametric coordinates */
+ static ErrorCode evalFcn(const double *params, const double *field, const int ndim, const int num_tuples,
+ double *work, double *result);
+
+ /** \brief Reverse-evaluation of parametric coordinates at physical space position */
+ 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);
+
+ /** \brief Evaluate the jacobian at a specified parametric position */
+ static ErrorCode jacobianFcn(const double *params, const double *verts, const int nverts, const int ndim,
+ double *work, double *result);
+
+ /** \brief Forward-evaluation of field at parametric coordinates */
+ static ErrorCode integrateFcn(const double *field, const double *verts, const int nverts, const int ndim, const int num_tuples,
+ double *work, double *result);
+
+ /** \brief Initialize this EvalSet */
+ 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 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);
+
+ static EvalSet eval_set()
+ {
+ return EvalSet(evalFcn, reverseEvalFcn, jacobianFcn, integrateFcn, initFcn, insideFcn);
+ }
+
+ static bool compatible(EntityType tp, int numv, EvalSet &eset)
+ {
+ if (tp == MBTRI && numv >= 3) {
+ eset = eval_set();
+ return true;
+ }
+ else return false;
+ }
+
+protected:
+
+ static const double corner[3][2];
+};// class LinearTri
+
+} // namespace moab
+
+#endif
diff --git a/src/LocalDiscretization/moab/QuadraticHex.hpp b/src/LocalDiscretization/moab/QuadraticHex.hpp
index 0801ad7..0e5cfb8 100644
--- a/src/LocalDiscretization/moab/QuadraticHex.hpp
+++ b/src/LocalDiscretization/moab/QuadraticHex.hpp
@@ -17,7 +17,8 @@ public:
/** \brief Reverse-evaluation of parametric coordinates at physical space position */
static ErrorCode reverseEvalFcn(EvalFcn eval, JacobianFcn jacob, InsideFcn ins,
const double *posn, const double *verts, const int nverts, const int ndim,
- const double tol, double *work, double *params, bool *is_inside);
+ const double iter_tol, const double inside_tol, double *work,
+ double *params, bool *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,
diff --git a/src/LocalDiscretization/moab/SpectralHex.hpp b/src/LocalDiscretization/moab/SpectralHex.hpp
index e6464c3..648cc70 100644
--- a/src/LocalDiscretization/moab/SpectralHex.hpp
+++ b/src/LocalDiscretization/moab/SpectralHex.hpp
@@ -19,7 +19,8 @@ public:
/** \brief Reverse-evaluation of parametric coordinates at physical space position */
static ErrorCode reverseEvalFcn(EvalFcn eval, JacobianFcn jacob, InsideFcn ins,
const double *posn, const double *verts, const int nverts, const int ndim,
- const double tol, double *work, double *params, bool *is_inside);
+ const double iter_tol, const double inside_tol, double *work,
+ double *params, bool *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,
diff --git a/src/LocalDiscretization/moab/SpectralQuad.hpp b/src/LocalDiscretization/moab/SpectralQuad.hpp
index 8716a43..a72f9b4 100644
--- a/src/LocalDiscretization/moab/SpectralQuad.hpp
+++ b/src/LocalDiscretization/moab/SpectralQuad.hpp
@@ -19,7 +19,8 @@ public:
/** \brief Reverse-evaluation of parametric coordinates at physical space position */
static ErrorCode reverseEvalFcn(EvalFcn eval, JacobianFcn jacob, InsideFcn ins,
const double *posn, const double *verts, const int nverts, const int ndim,
- const double tol, double *work, double *params, bool *is_inside);
+ const double iter_tol, const double inside_tol, double *work,
+ double *params, bool *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,
diff --git a/src/MergeMesh.cpp b/src/MergeMesh.cpp
index 4cc8158..72c4f66 100644
--- a/src/MergeMesh.cpp
+++ b/src/MergeMesh.cpp
@@ -171,7 +171,7 @@ moab::ErrorCode MergeMesh::find_merged_to(moab::EntityHandle &tree_root,
// check close-by leaves too
leaves_out.clear();
result = tree.distance_search(from.array(), mergeTol,
- leaves_out, 0.0, NULL, NULL, &tree_root);
+ leaves_out, mergeTol, 1.0e-6, NULL, NULL, &tree_root);
leaf_range2.clear();
for (std::vector<moab::EntityHandle>::iterator vit = leaves_out.begin();
vit != leaves_out.end(); vit++) {
diff --git a/src/SpatialLocator.cpp b/src/SpatialLocator.cpp
index d3b521e..f765ab9 100644
--- a/src/SpatialLocator.cpp
+++ b/src/SpatialLocator.cpp
@@ -3,6 +3,8 @@
#include "moab/ElemEvaluator.hpp"
#include "moab/AdaptiveKDTree.hpp"
+bool debug = true;
+
namespace moab
{
@@ -33,41 +35,47 @@ namespace moab
#ifdef USE_MPI
ErrorCode SpatialLocator::par_locate_points(Range &/*vertices*/,
- double /*rel_tol*/, double /*abs_tol*/)
+ const double /*rel_iter_tol*/, const double /*abs_iter_tol*/,
+ const double /*inside_tol*/)
{
return MB_UNSUPPORTED_OPERATION;
}
ErrorCode SpatialLocator::par_locate_points(const double */*pos*/, int /*num_points*/,
- double /*rel_tol*/, double /*abs_tol*/)
+ const double /*rel_iter_tol*/, const double /*abs_iter_tol*/,
+ const double /*inside_tol*/)
{
return MB_UNSUPPORTED_OPERATION;
}
#endif
ErrorCode SpatialLocator::locate_points(Range &verts,
- double rel_eps, double abs_eps)
+ const double rel_iter_tol, const double abs_iter_tol,
+ const double inside_tol)
{
assert(!verts.empty() && mbImpl->type_from_handle(*verts.rbegin()) == MBVERTEX);
std::vector<double> pos(3*verts.size());
ErrorCode rval = mbImpl->get_coords(verts, &pos[0]);
if (MB_SUCCESS != rval) return rval;
- rval = locate_points(&pos[0], verts.size(), rel_eps, abs_eps);
+ rval = locate_points(&pos[0], verts.size(), rel_iter_tol, abs_iter_tol, inside_tol);
if (MB_SUCCESS != rval) return rval;
return MB_SUCCESS;
}
ErrorCode SpatialLocator::locate_points(const double *pos, int num_points,
- double rel_eps, double abs_eps)
+ const double rel_iter_tol, const double abs_iter_tol,
+ const double inside_tol)
{
// initialize to tuple structure (p_ui, hs_ul, r[3]_d) (see header comments for locTable)
locTable.initialize(1, 0, 1, 3, num_points);
locTable.enableWriteAccess();
// pass storage directly into locate_points, since we know those arrays are contiguous
- ErrorCode rval = locate_points(pos, num_points, locTable.vul_wr, locTable.vr_wr, NULL, rel_eps, abs_eps);
+ ErrorCode rval = locate_points(pos, num_points, locTable.vul_wr, locTable.vr_wr, NULL, rel_iter_tol, abs_iter_tol,
+ inside_tol);
std::fill(locTable.vi_wr, locTable.vi_wr+num_points, 0);
+ locTable.set_n(num_points);
if (MB_SUCCESS != rval) return rval;
return MB_SUCCESS;
@@ -75,25 +83,27 @@ namespace moab
ErrorCode SpatialLocator::locate_points(Range &verts,
EntityHandle *ents, double *params, bool *is_inside,
- double rel_eps, double abs_eps)
+ const double rel_iter_tol, const double abs_iter_tol,
+ const double inside_tol)
{
assert(!verts.empty() && mbImpl->type_from_handle(*verts.rbegin()) == MBVERTEX);
std::vector<double> pos(3*verts.size());
ErrorCode rval = mbImpl->get_coords(verts, &pos[0]);
if (MB_SUCCESS != rval) return rval;
- return locate_points(&pos[0], verts.size(), ents, params, is_inside, rel_eps, abs_eps);
+ return locate_points(&pos[0], verts.size(), ents, params, is_inside, rel_iter_tol, abs_iter_tol, inside_tol);
}
ErrorCode SpatialLocator::locate_points(const double *pos, int num_points,
EntityHandle *ents, double *params, bool *is_inside,
- double rel_eps, double abs_eps)
+ const double rel_iter_tol, const double abs_iter_tol,
+ const double inside_tol)
{
-
- if (rel_eps && !abs_eps) {
+ 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);
- abs_eps = rel_eps * box.diagonal_length();
+ tmp_abs_iter_tol = rel_iter_tol * box.diagonal_length();
}
EntityHandle closest_leaf;
@@ -104,8 +114,8 @@ namespace moab
for (int i = 0; i < num_points; i++) {
int i3 = 3*i;
ents[i] = 0;
- if (abs_eps) {
- rval = myTree->distance_search(pos+i3, abs_eps, leaves, abs_eps, &dists);
+ if (tmp_abs_iter_tol) {
+ rval = myTree->distance_search(pos+i3, tmp_abs_iter_tol, leaves, tmp_abs_iter_tol, inside_tol, &dists);
if (MB_SUCCESS != rval) return rval;
if (!leaves.empty()) {
// get closest leaf
@@ -145,22 +155,71 @@ namespace moab
// loop over the range_leaf
bool tmp_inside;
+ bool *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++)
{
- bool *is_ptr = (is_inside ? is_inside+i : &tmp_inside);
rval = elemEval->set_ent_handle(*rit);
if (MB_SUCCESS != rval) return rval;
- rval = elemEval->reverse_eval(pos+i3, abs_eps, params+i3, is_ptr);
+ rval = elemEval->reverse_eval(pos+i3, tmp_abs_iter_tol, inside_tol, params+i3, is_ptr);
if (MB_SUCCESS != rval) return rval;
if (*is_ptr) {
- ents[i] = *rit;
+ ent = *rit;
break;
}
}
+ if (debug && !ent) {
+ std::cout << "Point " << i << " not found; point: ("
+ << pos[i3] << "," << pos[i3+1] << "," << pos[i3+2] << ")" << std::endl;
+ std::cout << "Source element candidates: " << std::endl;
+ range_leaf.print(" ");
+ for(Range::iterator rit = range_leaf.begin(); rit != range_leaf.end(); rit++)
+ {
+ std::cout << "Candidate " << CN::EntityTypeName(mbImpl->type_from_handle(*rit)) << " " << mbImpl->id_from_handle(*rit) << ": ";
+ rval = elemEval->set_ent_handle(*rit);
+ if (MB_SUCCESS != rval) return rval;
+ rval = elemEval->reverse_eval(pos+i3, tmp_abs_iter_tol, inside_tol, params+i3, is_ptr);
+ if (MB_SUCCESS != rval) return rval;
+ std::cout << "Parameters: (" << params[i3] << "," << params[i3+1] << "," << params[i3+2] << ")"
+ << " inside = " << *is_ptr << std::endl;
+ }
+ }
+ ents[i] = ent;
}
return MB_SUCCESS;
}
+ /* Count the number of located points in locTable
+ * Return the number of entries in locTable that have non-zero entity handles, which
+ * represents the number of points in targetEnts that were inside one element in sourceEnts
+ *
+ */
+ int SpatialLocator::local_num_located()
+ {
+ int num_located = locTable.get_n() - std::count(locTable.vul_rd, locTable.vul_rd+locTable.get_n(), 0);
+ if (num_located != (int)locTable.get_n()) {
+ unsigned long *nl = std::find(locTable.vul_rd, locTable.vul_rd+locTable.get_n(), 0);
+ if (nl) {
+ int idx = nl - locTable.vul_rd;
+ if (idx) {}
+ }
+ }
+ return num_located;
+ }
+
+ /* Count the number of located points in parLocTable
+ * Return the number of entries in parLocTable that have a non-negative index in on a remote
+ * proc in parLocTable, which gives the number of points located in at least one element in a
+ * remote proc's sourceEnts.
+ */
+ int SpatialLocator::remote_num_located()
+ {
+ int located = 0;
+ for (unsigned int i = 0; i < parLocTable.get_n(); i++)
+ if (parLocTable.vi_rd[2*i] != -1) located++;
+ return located;
+ }
} // namespace moab
diff --git a/src/moab/AdaptiveKDTree.hpp b/src/moab/AdaptiveKDTree.hpp
index 4a5ce1b..54d48e5 100644
--- a/src/moab/AdaptiveKDTree.hpp
+++ b/src/moab/AdaptiveKDTree.hpp
@@ -68,7 +68,8 @@ namespace moab {
* containing the point in that case.
* \param point Point to be located in tree
* \param leaf_out Leaf containing point
- * \param tol Tolerance below which a point is "in"
+ * \param iter_tol Tolerance for convergence of point search
+ * \param inside_tol Tolerance for inside element calculation
* \param multiple_leaves Some tree types can have multiple leaves containing a point;
* if non-NULL, this parameter is returned true if multiple leaves contain
* the input point
@@ -77,7 +78,8 @@ namespace moab {
*/
virtual ErrorCode point_search(const double *point,
EntityHandle& leaf_out,
- double tol = 0.0,
+ const double iter_tol = 1.0e-10,
+ const double inside_tol = 1.0e-6,
bool *multiple_leaves = NULL,
EntityHandle *start_node = NULL,
CartVect *params = NULL);
@@ -92,7 +94,8 @@ namespace moab {
* containing the point in that case.
* \param point Point to be located in tree
* \param leaf_it Iterator to leaf containing point
- * \param tol Tolerance below which a point is "in"
+ * \param iter_tol Tolerance for convergence of point search
+ * \param inside_tol Tolerance for inside element calculation
* \param multiple_leaves Some tree types can have multiple leaves containing a point;
* if non-NULL, this parameter is returned true if multiple leaves contain
* the input point
@@ -101,7 +104,8 @@ namespace moab {
*/
ErrorCode point_search(const double *point,
AdaptiveKDTreeIter& leaf_it,
- double tol = 0.0,
+ const double iter_tol = 1.0e-10,
+ const double inside_tol = 1.0e-6,
bool *multiple_leaves = NULL,
EntityHandle *start_node = NULL);
@@ -114,7 +118,8 @@ namespace moab {
* \param point Point to be located in tree
* \param distance Distance within which to query
* \param leaves_out Leaves within distance or containing point
- * \param tol Tolerance below which a point is "in"
+ * \param iter_tol Tolerance for convergence of point search
+ * \param inside_tol Tolerance for inside element calculation
* \param dists_out If non-NULL, will contain distsances to leaves
* \param params_out If non-NULL, will contain parameters of the point in the ents in leaves_out
* \param start_node Start from this tree node (non-NULL) instead of tree root (NULL)
@@ -122,7 +127,8 @@ namespace moab {
virtual ErrorCode distance_search(const double *point,
const double distance,
std::vector<EntityHandle>& leaves_out,
- double tol = 0.0,
+ const double iter_tol = 1.0e-10,
+ const double inside_tol = 1.0e-6,
std::vector<double> *dists_out = NULL,
std::vector<CartVect> *params_out = NULL,
EntityHandle *start_node = NULL);
diff --git a/src/moab/BVHTree.hpp b/src/moab/BVHTree.hpp
index 7f0c7ec..d435ebf 100644
--- a/src/moab/BVHTree.hpp
+++ b/src/moab/BVHTree.hpp
@@ -74,7 +74,8 @@ namespace moab {
* containing the point in that case.
* \param point Point to be located in tree
* \param leaf_out Leaf containing point
- * \param tol Tolerance below which a point is "in"
+ * \param iter_tol Tolerance for convergence of point search
+ * \param inside_tol Tolerance for inside element calculation
* \param multiple_leaves Some tree types can have multiple leaves containing a point;
* if non-NULL, this parameter is returned true if multiple leaves contain
* the input point
@@ -83,7 +84,8 @@ namespace moab {
*/
virtual ErrorCode point_search(const double *point,
EntityHandle& leaf_out,
- double tol = 0.0,
+ const double iter_tol = 1.0e-10,
+ const double inside_tol = 1.0e-6,
bool *multiple_leaves = NULL,
EntityHandle *start_node = NULL,
CartVect *params = NULL);
@@ -97,7 +99,8 @@ namespace moab {
* \param from_point Point to be located in tree
* \param distance Distance within which to query
* \param result_list Leaves within distance or containing point
- * \param tol Tolerance below which a point is "in"
+ * \param iter_tol Tolerance for convergence of point search
+ * \param inside_tol Tolerance for inside element calculation
* \param result_dists If non-NULL, will contain distsances to leaves
* \param result_params If non-NULL, will contain parameters of the point in the ents in leaves_out
* \param tree_root Start from this tree node (non-NULL) instead of tree root (NULL)
@@ -105,7 +108,8 @@ namespace moab {
virtual ErrorCode distance_search(const double from_point[3],
const double distance,
std::vector<EntityHandle>& result_list,
- double tol = 0.0,
+ const double iter_tol = 1.0e-10,
+ const double inside_tol = 1.0e-6,
std::vector<double> *result_dists = NULL,
std::vector<CartVect> *result_params = NULL,
EntityHandle *tree_root = NULL);
@@ -242,10 +246,13 @@ namespace moab {
ErrorCode find_point(const std::vector<double> &point,
const unsigned int &index,
- const double tol,
+ const double iter_tol,
+ const double inside_tol,
std::pair<EntityHandle, CartVect> &result);
- EntityHandle bruteforce_find(const double *point, const double tol);
+ EntityHandle bruteforce_find(const double *point,
+ const double iter_tol = 1.0e-10,
+ const double inside_tol = 1.0e-6);
int local_build_tree(std::vector<Node> &tree_nodes,
HandleDataVec::iterator begin,
diff --git a/src/moab/SpatialLocator.hpp b/src/moab/SpatialLocator.hpp
index 780c086..bb639d0 100644
--- a/src/moab/SpatialLocator.hpp
+++ b/src/moab/SpatialLocator.hpp
@@ -62,26 +62,44 @@ namespace moab {
/* locate a set of vertices, Range variant */
ErrorCode locate_points(Range &vertices,
EntityHandle *ents, double *params, bool *is_inside = NULL,
- double rel_tol = 0.0, double abs_tol = 0.0);
+ 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,
- double rel_tol = 0.0, double abs_tol = 0.0);
+ 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 vertices or entity centroids, storing results on TupleList in this class
* Locate a set of vertices or entity centroids, storing the detailed results in member
* variable (TupleList) locTable (see comments on locTable for structure of that tuple).
*/
ErrorCode locate_points(Range &ents,
- double rel_tol = 0.0, double abs_tol = 0.0);
+ 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 variable (TupleList) locTable
* (see comments on locTable for structure of that tuple).
*/
ErrorCode locate_points(const double *pos, int num_points,
- double rel_tol = 0.0, double abs_tol = 0.0);
+ const double rel_iter_tol = 1.0e-10, const double abs_iter_tol = 1.0e-10,
+ const double inside_tol = 1.0e-6);
+
+ /* Count the number of located points in locTable
+ * Return the number of entries in locTable that have non-zero entity handles, which
+ * represents the number of points in targetEnts that were inside one element in sourceEnts
+ *
+ */
+ int local_num_located();
+
+ /* Count the number of located points in parLocTable
+ * Return the number of entries in parLocTable that have a non-negative index in on a remote
+ * proc in parLocTable, which gives the number of points located in at least one element in a
+ * remote proc's sourceEnts.
+ */
+ int remote_num_located();
#ifdef USE_MPI
/* locate a set of vertices or entity centroids, storing results on TupleList in this class
@@ -90,7 +108,8 @@ namespace moab {
* structure of those tuples).
*/
ErrorCode par_locate_points(Range &vertices,
- double rel_tol = 0.0, double abs_tol = 0.0);
+ 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
@@ -98,7 +117,8 @@ namespace moab {
* structure of those tuples).
*/
ErrorCode par_locate_points(const double *pos, int num_points,
- double rel_tol = 0.0, double abs_tol = 0.0);
+ const double rel_iter_tol = 1.0e-10, const double abs_iter_tol = 1.0e-10,
+ const double inside_tol = 1.0e-6);
#endif
/* return the tree */
@@ -127,14 +147,15 @@ namespace moab {
const ElemEvaluator *elem_eval() const {return elemEval;}
/* set elemEval */
- void elem_eval(ElemEvaluator *eval) {elemEval = eval;}
+ void elem_eval(ElemEvaluator *eval) {elemEval = eval; if (myTree) myTree->set_eval(eval);}
private:
/* locate a point */
ErrorCode locate_point(const double *pos,
EntityHandle &ent, double *params, bool *is_inside = NULL,
- double rel_tol = 0.0, double abs_tol = 0.0);
+ const double rel_iter_tol = 1.0e-10, const double abs_iter_tol = 1.0e-10,
+ const double inside_tol = 1.0e-6);
/* MOAB instance */
Interface* mbImpl;
@@ -185,9 +206,10 @@ namespace moab {
inline ErrorCode SpatialLocator::locate_point(const double *pos,
EntityHandle &ent, double *params, bool *is_inside,
- double rel_tol, double abs_tol)
+ const double rel_iter_tol, const double abs_iter_tol,
+ const double inside_tol)
{
- return locate_points(pos, 1, &ent, params, is_inside, rel_tol, abs_tol);
+ 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)
diff --git a/src/moab/Tree.hpp b/src/moab/Tree.hpp
index c2a117c..4d2694a 100644
--- a/src/moab/Tree.hpp
+++ b/src/moab/Tree.hpp
@@ -100,7 +100,8 @@ namespace moab {
* containing the point in that case.
* \param point Point to be located in tree
* \param leaf_out Leaf containing point
- * \param tol Tolerance below which a point is "in"
+ * \param iter_tol Tolerance for convergence of point search
+ * \param inside_tol Tolerance for inside element calculation
* \param multiple_leaves Some tree types can have multiple leaves containing a point;
* if non-NULL, this parameter is returned true if multiple leaves contain
* the input point
@@ -109,7 +110,8 @@ namespace moab {
*/
virtual ErrorCode point_search(const double *point,
EntityHandle& leaf_out,
- double tol = 0.0,
+ const double iter_tol = 1.0e-10,
+ const double inside_tol = 1.0e-6,
bool *multiple_leaves = NULL,
EntityHandle *start_node = NULL,
CartVect *params = NULL) = 0;
@@ -123,7 +125,8 @@ namespace moab {
* \param point Point to be located in tree
* \param distance Distance within which to query
* \param leaves_out Leaves within distance or containing point
- * \param tol Tolerance below which a point is "in"
+ * \param iter_tol Tolerance for convergence of point search
+ * \param inside_tol Tolerance for inside element calculation
* \param dists_out If non-NULL, will contain distsances to leaves
* \param params_out If non-NULL, will contain parameters of the point in the ents in leaves_out
* \param start_node Start from this tree node (non-NULL) instead of tree root (NULL)
@@ -131,7 +134,8 @@ namespace moab {
virtual ErrorCode distance_search(const double *point,
const double distance,
std::vector<EntityHandle>& leaves_out,
- double params_tol = 0.0,
+ const double iter_tol = 1.0e-10,
+ const double inside_tol = 1.0e-6,
std::vector<double> *dists_out = NULL,
std::vector<CartVect> *params_out = NULL,
EntityHandle *start_node = NULL) = 0;
diff --git a/test/elem_eval_test.cpp b/test/elem_eval_test.cpp
index 6f07d7a..1e669c3 100644
--- a/test/elem_eval_test.cpp
+++ b/test/elem_eval_test.cpp
@@ -7,6 +7,7 @@
#include "moab/Core.hpp"
#include "moab/ReadUtilIface.hpp"
#include "moab/ElemEvaluator.hpp"
+#include "moab/LinearTri.hpp"
#include "moab/LinearQuad.hpp"
#include "moab/LinearHex.hpp"
#include "moab/LinearTet.hpp"
@@ -22,6 +23,7 @@ std::string TestDir(".");
using namespace moab;
+void test_linear_tri();
void test_linear_quad();
void test_linear_hex();
void test_linear_tet();
@@ -78,11 +80,11 @@ void test_eval(ElemEvaluator &ee, bool test_integrate)
if (!test_integrate) return;
- // tag equal to coordinates should integrate to avg position, since volume is 1
+ // tag equal to coordinates should integrate to avg position, test that
Tag tag;
// make a temporary tag and set it on vertices to the vertex positions
rval = ee.get_moab()->tag_get_handle(NULL, 3, MB_TYPE_DOUBLE, tag, MB_TAG_DENSE | MB_TAG_CREAT); CHECK_ERR(rval);
- rval = ee.get_moab()->tag_set_data(tag, ee.get_vert_handles(), ee.get_num_verts(), hex_verts[0].array()); CHECK_ERR(rval);
+ rval = ee.get_moab()->tag_set_data(tag, ee.get_vert_handles(), ee.get_num_verts(), ee.get_vert_pos()); CHECK_ERR(rval);
// set that temporary tag on the evaluator so that's what gets integrated
rval = ee.set_tag_handle(tag, 0); CHECK_ERR(rval);
@@ -95,8 +97,8 @@ void test_eval(ElemEvaluator &ee, bool test_integrate)
EvalSet es;
EntityHandle eh = ee.get_ent_handle();
rval = EvalSet::get_eval_set(ee.get_moab(), eh, es); CHECK_ERR(rval);
- rval = (*es.integrateFcn)(&one[0], hex_verts[0].array(), ee.get_num_verts(), ee.get_moab()->dimension_from_handle(eh),
- 1, NULL, &measure); CHECK_ERR(rval);
+ rval = (*es.integrateFcn)(&one[0], ee.get_vert_pos(), ee.get_num_verts(), ee.get_moab()->dimension_from_handle(eh),
+ 1, ee.get_work_space(), &measure); CHECK_ERR(rval);
if (measure) integral /= measure;
// check against avg of entity's vertices' positions
@@ -136,6 +138,7 @@ int main()
{
int failures = 0;
+// failures += RUN_TEST(test_linear_tri); currently failing linear tri, bad formulation, working on it...
failures += RUN_TEST(test_linear_quad);
failures += RUN_TEST(test_linear_hex);
failures += RUN_TEST(test_quadratic_hex);
@@ -144,6 +147,28 @@ int main()
return failures;
}
+void test_linear_tri()
+{
+ Core mb;
+ Range verts;
+ double tri_verts[] = {-1.0, -1.0, -1.0,
+ 1.0, -1.0, -1.0,
+ -1.0, 1.0, -1.0};
+
+
+ ErrorCode rval = mb.create_vertices(tri_verts, 3, verts); CHECK_ERR(rval);
+ EntityHandle tri;
+ std::vector<EntityHandle> connect;
+ std::copy(verts.begin(), verts.end(), std::back_inserter(connect));
+ rval = mb.create_element(MBTRI, &connect[0], 3, tri); CHECK_ERR(rval);
+
+ ElemEvaluator ee(&mb, tri, 0);
+ ee.set_tag_handle(0, 0);
+ ee.set_eval_set(MBTRI, LinearTri::eval_set());
+
+ test_eval(ee, true);
+}
+
void test_linear_quad()
{
Core mb;
diff --git a/tools/mbcoupler/Coupler.cpp b/tools/mbcoupler/Coupler.cpp
index 71231b8..0a10d14 100644
--- a/tools/mbcoupler/Coupler.cpp
+++ b/tools/mbcoupler/Coupler.cpp
@@ -614,7 +614,7 @@ ErrorCode Coupler::nat_param(double xyz[3],
if (epsilon) {
std::vector<double> dists;
std::vector<EntityHandle> leaves;
- result = myTree->distance_search(xyz, epsilon, leaves, 0.0, &dists, NULL, &localRoot);
+ result = myTree->distance_search(xyz, epsilon, leaves, 1.0e-10, 1.0e-6, &dists, NULL, &localRoot);
if (leaves.empty())
// not found returns success here, with empty list, just like case with no epsilon
return MB_SUCCESS;
@@ -631,7 +631,7 @@ ErrorCode Coupler::nat_param(double xyz[3],
}
}
else {
- result = myTree->point_search(xyz, treeiter, 0.0, NULL, &localRoot);
+ result = myTree->point_search(xyz, treeiter, 1.0e-10, 1.0e-6, NULL, &localRoot);
if(MB_ENTITY_NOT_FOUND==result) //point is outside of myTree's bounding box
return MB_SUCCESS;
else if (MB_SUCCESS != result) {
diff --git a/tools/mbcoupler/DataCoupler.cpp b/tools/mbcoupler/DataCoupler.cpp
index 71f3497..421f45c 100644
--- a/tools/mbcoupler/DataCoupler.cpp
+++ b/tools/mbcoupler/DataCoupler.cpp
@@ -40,7 +40,8 @@ DataCoupler::DataCoupler(Interface *impl,
for (; pit != source_ents.pair_end(); pit++) {
EntityType this_type = mbImpl->type_from_handle(pit->first);
if (last_type == this_type) continue;
- myLocator->elem_eval()->set_eval_set(pit->first);
+ ErrorCode rval = myLocator->elem_eval()->set_eval_set(pit->first);
+ if (MB_SUCCESS != rval) throw(rval);
last_type = this_type;
}
}
@@ -61,19 +62,19 @@ DataCoupler::~DataCoupler()
}
ErrorCode DataCoupler::locate_points(Range &targ_ents,
- double rel_eps,
- double abs_eps)
+ const double rel_iter_tol, const double abs_iter_tol,
+ const double inside_tol)
{
targetEnts = targ_ents;
- return myLocator->locate_points(targ_ents, rel_eps, abs_eps);
+ return myLocator->locate_points(targ_ents, rel_iter_tol, abs_iter_tol, inside_tol);
}
ErrorCode DataCoupler::locate_points(double *xyz, int num_points,
- double rel_eps,
- double abs_eps)
+ const double rel_iter_tol, const double abs_iter_tol,
+ const double inside_tol)
{
- return myLocator->locate_points(xyz, num_points, rel_eps, abs_eps);
+ return myLocator->locate_points(xyz, num_points, rel_iter_tol, abs_iter_tol, inside_tol);
}
ErrorCode DataCoupler::interpolate(/*DataCoupler::Method*/ int method,
diff --git a/tools/mbcoupler/DataCoupler.hpp b/tools/mbcoupler/DataCoupler.hpp
index 3ad1d74..bfe2ba5 100644
--- a/tools/mbcoupler/DataCoupler.hpp
+++ b/tools/mbcoupler/DataCoupler.hpp
@@ -65,8 +65,9 @@ public:
* This is a pass-through function to SpatialLocator::locate_points
* \param xyz Point locations (interleaved) being located
* \param num_points Number of points in xyz
- * \param rel_eps Relative tolerance for the non-linear iteration inside a given element
- * \param abs_eps Absolute tolerance for the non-linear iteration inside a given element
+ * \param rel_iter_tol Relative tolerance for non-linear iteration
+ * \param abs_iter_tol Relative tolerance for non-linear iteration, usually 10^-10 or so
+ * \param inside_tol Tolerance of is_inside evaluation, usually 10^-6 or so
* \param loc_results Tuple list containing the results; two types of results are possible,
* controlled by value of store_local parameter:
* store_local = true: each tuple T[j] consists of (p, i), p = proc with src element containing
@@ -76,13 +77,15 @@ public:
* \param store_local If true, stores the located points on SpatialLocator
*/
ErrorCode locate_points(double *xyz, int num_points,
- double rel_eps = 0.0, double abs_eps = 0.0);
+ const double rel_iter_tol = 1.0e-10, const double abs_iter_tol = 1.0e-10,
+ const double inside_tol = 1.0e-6);
/* \brief Locate points on the source mesh
* This is a pass-through function to SpatialLocator::locate_points
* \param ents Target entities being located
- * \param rel_eps Relative tolerance for the non-linear iteration inside a given element
- * \param abs_eps Absolute tolerance for the non-linear iteration inside a given element
+ * \param rel_iter_tol Relative tolerance for non-linear iteration
+ * \param abs_iter_tol Relative tolerance for non-linear iteration, usually 10^-10 or so
+ * \param inside_tol Tolerance of is_inside evaluation, usually 10^-6 or so
* \param loc_results Tuple list containing the results; two types of results are possible,
* controlled by value of store_local parameter:
* store_local = true: each tuple T[j] consists of (p, i), p = proc with src element containing
@@ -92,8 +95,8 @@ public:
* \param store_local If true, stores the located points on SpatialLocator
*/
ErrorCode locate_points(Range &ents,
- double rel_eps = 0.0,
- double abs_eps = 0.0);
+ const double rel_iter_tol = 1.0e-10, const double abs_iter_tol = 1.0e-10,
+ const double inside_tol = 1.0e-6);
/* \brief Interpolate data from the source mesh onto points
* All entities/points or, if tuple_list is input, only those points
diff --git a/tools/skin.cpp b/tools/skin.cpp
index 761f7d2..25caa43 100644
--- a/tools/skin.cpp
+++ b/tools/skin.cpp
@@ -525,7 +525,7 @@ ErrorCode merge_duplicate_vertices( Interface& moab, const double epsilon )
if (MB_SUCCESS != rval) return rval;
leaves.clear();;
- rval = tree.distance_search(coords, epsilon, leaves, epsilon);
+ rval = tree.distance_search(coords, epsilon, leaves, epsilon, epsilon);
if (MB_SUCCESS != rval) return rval;
Range near;
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