[MOAB-dev] commit/MOAB: 16 new changesets

commits-noreply at bitbucket.org commits-noreply at bitbucket.org
Wed Nov 27 19:40:54 CST 2013


16 new commits in MOAB:

https://bitbucket.org/fathomteam/moab/commits/d43e32145950/
Changeset:   d43e32145950
Branch:      None
User:        tautges
Date:        2013-11-28 01:33:08
Summary:     Adding a Lloyd smoother tool, and a test of it.

Affected #:  6 files

diff --git a/.gitignore b/.gitignore
index e824a6b..d842c2b 100644
--- a/.gitignore
+++ b/.gitignore
@@ -165,6 +165,7 @@ test/io/vtk_test
 test/kd_tree_test
 test/kd_tree_time
 test/kd_tree_tool
+test/lloyd_smoother_test
 test/mbcn_test
 test/mbfacet_test
 test/mbground_test

diff --git a/src/LloydSmoother.cpp b/src/LloydSmoother.cpp
new file mode 100644
index 0000000..49177d0
--- /dev/null
+++ b/src/LloydSmoother.cpp
@@ -0,0 +1,220 @@
+#include "moab/LloydSmoother.hpp"
+#include "moab/Error.hpp"
+#include "moab/Skinner.hpp"
+#include "moab/CN.hpp"
+#include "moab/CartVect.hpp"
+#include "moab/BoundBox.hpp"
+
+#ifdef USE_MPI
+#include "moab/ParallelComm.hpp"
+#endif
+
+#include <iostream>
+
+#define RR(a) if (MB_SUCCESS != rval) {        \
+    errorHandler->set_last_error(a);                \
+    return rval;}
+
+namespace moab {
+  
+LloydSmoother::LloydSmoother(Interface *impl, ParallelComm *pc, Range &elms, Tag ftag,
+                             double at, double rt) 
+        : mbImpl(impl), myPcomm(pc), myElems(elms), fixedTag(ftag), absTol(at), relTol(rt),
+          reportIts(0), numIts(0), iCreatedTag(false)
+{
+  ErrorCode rval = mbImpl->query_interface(errorHandler);
+  if (rval) throw rval;
+}
+
+LloydSmoother::~LloydSmoother() 
+{
+  if (iCreatedTag && fixedTag) {
+    ErrorCode rval = mbImpl->tag_delete(fixedTag);
+    if (rval) throw rval;
+  }
+}
+    
+ErrorCode LloydSmoother::perform_smooth() 
+{
+  ErrorCode rval;
+
+    // first figure out tolerance to use
+  if (0 > absTol) {
+      // no tolerance set - get one relative to bounding box around elements
+    BoundBox bb;
+    rval = bb.update(*mbImpl, myElems);
+    RR("Failed to compute bounding box around elements.");
+    absTol = relTol * bb.diagonal_length();
+  }
+
+    // initialize if we need to
+  rval = initialize();
+  RR("Failed to initialize.");
+
+    // get all vertices
+  Range verts;
+  rval = mbImpl->get_adjacencies(myElems, 0, false, verts, Interface::UNION);
+  RR("Failed to get all vertices.");
+  
+    // perform Lloyd relaxation:
+    // 1. setup: set vertex centroids from vertex coords; filter to owned verts; get fixed tags
+
+    // get all verts coords into tag; don't need to worry about filtering out fixed verts, 
+    // we'll just be setting to their fixed coords
+  std::vector<double> vcentroids(3*verts.size());
+  rval = mbImpl->get_coords(verts, &vcentroids[0]); RR("Failed to get vert coords.");
+
+  Tag centroid;
+  rval = mbImpl->tag_get_handle("", 3, MB_TYPE_DOUBLE, centroid, MB_TAG_CREAT | MB_TAG_DENSE); 
+  RR("Couldn't get tag handle.");
+  rval = mbImpl->tag_set_data(centroid, verts, &vcentroids[0]); RR("Failed setting centroid tag.");
+
+  Range owned_verts, shared_owned_verts;
+#ifdef USE_MPI
+    // filter verts down to owned ones and get fixed tag for them
+  if (myPcomm->size() > 1) {
+    rval = myPcomm->filter_pstatus(verts, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, &owned_verts);
+    RR("Failed to filter on pstatus.");
+      // get shared owned verts, for exchanging tags
+    rval = myPcomm->filter_pstatus(owned_verts, PSTATUS_SHARED, PSTATUS_AND, -1, &shared_owned_verts);
+    RR("Failed to filter for shared owned.");
+      // workaround: if no shared owned verts, put a non-shared one in the list, to prevent exchanging tags
+      // for all shared entities
+    if (shared_owned_verts.empty()) shared_owned_verts.insert(*verts.begin());
+  }
+  else
+    owned_verts = verts;
+#else
+  owned_verts = verts;
+#endif
+  
+  std::vector<unsigned char> fix_tag(owned_verts.size());
+  rval = mbImpl->tag_get_data(fixedTag, owned_verts, &fix_tag[0]); RR("Failed to get fixed tag.");
+
+    // now fill vcentroids array with positions of just owned vertices, since those are the ones
+    // we're actually computing
+  vcentroids.resize(3*owned_verts.size());
+  rval = mbImpl->tag_get_data(centroid, owned_verts, &vcentroids[0]); RR("Failed to get centroid tag.");
+
+    // some declarations for later iterations
+  std::vector<double> fcentroids(3*myElems.size()); // fcentroids for element centroids
+  std::vector<double> ctag(3*CN::MAX_NODES_PER_ELEMENT);  // temporary coordinate storage for verts bounding an element
+  const EntityHandle *conn;  // const ptr & size to elem connectivity
+  int nconn;
+  Range::iterator eit, vit;  // for iterating over elems, verts
+  int e, v;  // for indexing into centroid vectors
+  std::vector<EntityHandle> adj_elems;  // used in vertex iteration
+  
+    // 2. while !converged
+  double resid = DBL_MAX;
+  numIts = 0;
+  while (resid > absTol) {
+    numIts++;
+    resid = 0.0;
+    
+    // 2a. foreach elem: centroid = sum(vertex centroids)/num_verts_in_cell
+    for (eit = myElems.begin(), e = 0; eit != myElems.end(); eit++, e++) {
+        // get verts for this elem
+      rval = mbImpl->get_connectivity(*eit, conn, nconn); RR("Failed to get connectivity.");
+        // get centroid tags for those verts
+      rval = mbImpl->tag_get_data(centroid, conn, nconn, &ctag[0]); RR("Failed to get centroid.");
+      fcentroids[3*e+0] = fcentroids[3*e+1] = fcentroids[3*e+2] = 0.0;
+      for (v = 0; v < nconn; v++) {
+        fcentroids[3*e+0] += ctag[3*v+0];
+        fcentroids[3*e+1] += ctag[3*v+1];
+        fcentroids[3*e+2] += ctag[3*v+2];
+      }
+      for (v = 0; v < 3; v++) fcentroids[3*e+v] /= nconn;
+    }
+    rval = mbImpl->tag_set_data(centroid, myElems, &fcentroids[0]); RR("Failed to set elem centroid.");
+
+      // 2b. foreach owned vertex: 
+    for (vit = owned_verts.begin(), v = 0; vit != owned_verts.end(); vit++, v++) {
+        // if !fixed
+      if (fix_tag[v]) continue;
+        // vertex centroid = sum(cell centroids)/ncells
+      adj_elems.clear();
+      rval = mbImpl->get_adjacencies(&(*vit), 1, 2, false, adj_elems); RR("Failed getting adjs.");
+      rval = mbImpl->tag_get_data(centroid, &adj_elems[0], adj_elems.size(), &fcentroids[0]); 
+      RR("Failed to get elem centroid.");
+      double vnew[] = {0.0, 0.0, 0.0};
+      for (e = 0; e < (int)adj_elems.size(); e++) {
+        vnew[0] += fcentroids[3*e+0];
+        vnew[1] += fcentroids[3*e+1];
+        vnew[2] += fcentroids[3*e+2];
+      }
+      for (e = 0; e < 3; e++) vnew[e] /= adj_elems.size();
+      double delta = (CartVect(vnew)-CartVect(&vcentroids[3*v])).length();
+      resid = (v ? std::max(resid, delta) : delta);
+      for (e = 0; e < 3; e++) vcentroids[3*v+e] = vnew[e];
+    }
+
+      // set the centroid tag; having them only in vcentroids array isn't enough, as vertex centroids are
+      // accessed randomly in loop over faces
+    rval = mbImpl->tag_set_data(centroid, owned_verts, &vcentroids[0]); RR("Failed to set vertex centroid.");
+
+#ifdef USE_MPI
+    // 2c. exchange tags on owned verts
+    if (myPcomm->size() > 1) {
+      rval = pcomm->exchange_tags(centroid, shared_owned_verts); RR("Failed to exchange tags.");
+    }
+#endif
+
+    if (reportIts && !(numIts%reportIts)) {
+      double global_max = resid;
+#ifdef USE_MPI
+        // global reduce for maximum delta, then report it
+      if (myPcomm->size() > 1)
+        MPI_Reduce(&resid, &global_max, 1, MPI_DOUBLE, MPI_MAX, 0, myPcomm->comm());
+      if (!pcomm->rank()) 
+#endif
+        std::cout << "Max residual = " << global_max << std::endl;
+    }
+
+  } // end while
+  
+    // write the tag back onto vertex coordinates
+  rval = mbImpl->set_coords(owned_verts, &vcentroids[0]); RR("Failed to set vertex coords.");
+
+  return MB_SUCCESS;
+}
+
+ErrorCode LloydSmoother::initialize()
+{
+  ErrorCode rval = MB_SUCCESS;
+  
+    // first, check for tag; if we don't have it, make one and mark skin as fixed
+  if (!fixedTag) {
+    unsigned char fixed = 0x0;
+    rval = mbImpl->tag_get_handle("", 1, MB_TYPE_OPAQUE, fixedTag, MB_TAG_DENSE|MB_TAG_CREAT, &fixed);
+    RR("Trouble making fixed tag.");
+    iCreatedTag = true;
+  
+      // get the skin; get facets, because we might need to filter on shared entities
+    Skinner skinner(mbImpl);
+    Range skin, skin_verts;
+    rval = skinner.find_skin(0, myElems, false, skin);
+    RR("Unable to find skin.");
+
+#ifdef USE_MPI
+      // need to do a little extra if we're working in parallel
+    if (myPcomm) {
+        // filter out ghost and interface facets
+      rval = myPcomm->filter_pstatus(skin, PSTATUS_GHOST|PSTATUS_INTERFACE, PSTATUS_NOT);
+      RR("Failed to filter on shared status.");
+    }
+#endif
+      // get the vertices from those entities
+    rval = mbImpl->get_adjacencies(skin, 0, false, skin_verts, Interface::UNION);
+    RR("Trouble getting vertices.");
+    
+      // mark them fixed
+    std::vector<unsigned char> marks(skin.size(), 0x1);
+    rval = mbImpl->tag_set_data(fixedTag, skin_verts, &marks[0]);
+    RR("Unable to set tag on skin.");
+  }
+  
+  return MB_SUCCESS;
+}
+
+}

diff --git a/src/Makefile.am b/src/Makefile.am
index 0eb07c5..3a5e142 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -64,6 +64,7 @@ libMOAB_la_SOURCES = \
   HigherOrderFactory.cpp \
   HomXform.cpp \
   Internals.hpp \
+  LloydSmoother.cpp \
   MBCNArrays.hpp \
   MergeMesh.cpp \
   MeshSet.cpp \

diff --git a/src/moab/LloydSmoother.hpp b/src/moab/LloydSmoother.hpp
new file mode 100644
index 0000000..5ab9771
--- /dev/null
+++ b/src/moab/LloydSmoother.hpp
@@ -0,0 +1,145 @@
+/** @class LloydSmoother.cpp \n
+ * \brief Perform Lloyd relaxation on a mesh and its dual \n
+ *
+ * Briefly, Lloyd relaxation is a technique to smooth out a mesh.  The centroid of each cell is computed from its
+ * vertex positions, then vertices are placed at the average of their connected cells' centroids, and the process
+ * iterates until convergence.
+ *
+ * In the parallel algorithm, an extra ghost layer of cells is exchanged.  This allows us to compute the centroids
+ * for boundary cells on each processor where they appear; this eliminates the need for one round of data exchange
+ * (for those centroids) between processors.  New vertex positions must be sent from owning processors to processors
+ * sharing those vertices.  Convergence is measured as the maximum distance moved by any vertex.  
+ * 
+ */
+
+#ifndef LLOYDSMOOTHER_HPP
+#define LLOYDSMOOTHER_HPP
+
+#include "moab/Interface.hpp"
+
+namespace moab {
+
+class ParallelComm;
+class Error;
+
+class LloydSmoother
+{
+public:
+
+    /* \brief Constructor
+     * Bare constructor, data input to this class through methods.
+     * \param impl The MOAB instance for this smoother
+     */
+  LloydSmoother(Interface *impl);
+  
+    /* \brief Constructor
+     * Convenience constructor, data input directly
+     * \param impl The MOAB instance for this smoother
+     * \param pcomm The ParallelComm instance by which this mesh is parallel
+     * \param elems The mesh to be smoothed
+     * \param fixed_tag The tag marking which vertices are fixed
+     * \param abs_tol Absolute tolerance measuring convergence
+     * \param rel_tol Relative tolerance measuring convergence
+     */
+  LloydSmoother(Interface *impl, ParallelComm *pc, Range &elems, Tag fixed_tag = 0,
+                double abs_tol = -1.0, double rel_tol = 1.0e-6);
+
+    /* \brief Destructor
+     */
+  ~LloydSmoother();
+
+    /* \brief perform smoothing operation
+     */
+  ErrorCode perform_smooth();
+  
+    /* \brief get instance
+     */
+  Interface *mb_impl() {return mbImpl;}
+  
+    /* \brief get/set ParallelComm
+     */
+  ParallelComm *pcomm() {return myPcomm;}
+  
+    /* \brief get/set ParallelComm
+     */
+  void pcomm(ParallelComm *pc) {myPcomm = pc;}
+
+    /* \brief get/set elements
+     */
+  Range &elems() {return myElems;}
+  
+    /* \brief get/set elements
+     */
+  const Range &elems() const {return myElems;}
+  
+    /* \brief get/set fixed tag
+     */
+  Tag fixed_tag() {return fixedTag;}
+
+    /* \brief get/set fixed tag
+     */
+  void fixed_tag(Tag fixed) {fixedTag = fixed;}
+
+    /* \brief get/set tolerance
+     */
+  double abs_tol() {return absTol;}
+  
+    /* \brief get/set tolerance
+     */
+  void abs_tol(double tol) {absTol = tol;}
+  
+    /* \brief get/set tolerance
+     */
+  double rel_tol() {return relTol;}
+  
+    /* \brief get/set tolerance
+     */
+  void rel_tol(double tol) {relTol = tol;}
+
+    /* \brief get/set numIts
+     */
+  int num_its() {return numIts;}
+  void num_its(int num) {numIts = num;}
+        
+    /* \brief get/set reportIts
+     */
+  int report_its() {return reportIts;}
+  void report_its(int num) {reportIts = num;}
+        
+  
+private:
+
+    //- initialize some things in certain cases
+  ErrorCode initialize();
+  
+    //- MOAB instance
+  Interface *mbImpl;
+
+    //- ParallelComm
+  ParallelComm *myPcomm;
+  
+    //- elements to smooth
+  Range myElems;
+  
+    //- tag marking which vertices are fixed, 0 = not fixed, otherwise fixed
+  Tag fixedTag;
+  
+    //- tolerances
+  double absTol, relTol;
+
+    //- error handler for this class
+  Error *errorHandler;
+
+    //- number of iterations between reporting
+  int reportIts;
+
+    //- number of iterations taken during smooth
+  int numIts;
+
+    //- keep track of whether I created the fixed tag
+  bool iCreatedTag;
+};
+    
+} // namespace moab
+
+#endif

diff --git a/test/Makefile.am b/test/Makefile.am
index 99c5211..0aabcce 100644
--- a/test/Makefile.am
+++ b/test/Makefile.am
@@ -48,7 +48,8 @@ TESTS = range_test \
 	coords_connect_iterate \
 	elem_eval_test \
 	spatial_locator_test \
-	test_boundbox
+	test_boundbox \
+	lloyd_smoother_test
 
 if HDF5_FILE
   TESTS += mbfacet_test \
@@ -147,6 +148,7 @@ coords_connect_iterate_SOURCES = coords_connect_iterate.cpp
 coords_connect_iterate_CPPFLAGS = $(AM_CPPFLAGS) $(CPPFLAGS)
 
 test_boundbox_SOURCES = test_boundbox.cpp
+lloyd_smoother_test_SOURCES = lloyd_smoother_test.cpp
 
 if PARALLEL
 moab_test_CPPFLAGS += -I$(top_srcdir)/src/parallel

diff --git a/test/lloyd_smoother_test.cpp b/test/lloyd_smoother_test.cpp
new file mode 100644
index 0000000..8d88112
--- /dev/null
+++ b/test/lloyd_smoother_test.cpp
@@ -0,0 +1,41 @@
+#include "moab/Core.hpp"
+#include "moab/LloydSmoother.hpp"
+#include "TestUtil.hpp"
+
+#ifdef MESHDIR
+std::string TestDir( STRINGIFY(MESHDIR) );
+#else
+std::string TestDir(".");
+#endif
+
+std::string filename = TestDir + "/surfrandomtris-4part.h5m";
+
+using namespace moab;
+
+int main(int argc, char**argv) 
+{
+  if (argc > 1) filename = std::string(argv[1]);
+  Core mb;
+  ErrorCode rval = mb.load_file(filename.c_str());
+  CHECK_ERR(rval);
+  
+  Range elems;
+  rval = mb.get_entities_by_dimension(0, 3, elems);
+  CHECK_ERR(rval);
+  if (elems.empty()) {
+    rval = mb.get_entities_by_dimension(0, 2, elems);
+    CHECK_ERR(rval);
+  }
+  if (elems.empty()) {
+    std::cout << "Mesh must have faces or regions for this test." << std::endl;
+    CHECK_ERR(MB_FAILURE);
+  }
+
+  LloydSmoother ll(&mb, NULL, elems);
+  ll.report_its(10);
+  rval = ll.perform_smooth();
+  CHECK_ERR(rval);
+  
+  std::cout << "Mesh smoothed in " << ll.num_its() << " iterations." << std::endl;
+}
+


https://bitbucket.org/fathomteam/moab/commits/562629a3a6d4/
Changeset:   562629a3a6d4
Branch:      None
User:        tautges
Date:        2013-11-28 01:33:08
Summary:     Adding ability to smooth based on a tag instead of coords.

Affected #:  3 files

diff --git a/src/LloydSmoother.cpp b/src/LloydSmoother.cpp
index 49177d0..09546b5 100644
--- a/src/LloydSmoother.cpp
+++ b/src/LloydSmoother.cpp
@@ -17,9 +17,9 @@
 
 namespace moab {
   
-LloydSmoother::LloydSmoother(Interface *impl, ParallelComm *pc, Range &elms, Tag ftag,
-                             double at, double rt) 
-        : mbImpl(impl), myPcomm(pc), myElems(elms), fixedTag(ftag), absTol(at), relTol(rt),
+LloydSmoother::LloydSmoother(Interface *impl, ParallelComm *pc, Range &elms, Tag ctag,
+                             Tag ftag, double at, double rt) 
+        : mbImpl(impl), myPcomm(pc), myElems(elms), coordsTag(ctag), fixedTag(ftag), absTol(at), relTol(rt),
           reportIts(0), numIts(0), iCreatedTag(false)
 {
   ErrorCode rval = mbImpl->query_interface(errorHandler);
@@ -62,7 +62,12 @@ ErrorCode LloydSmoother::perform_smooth()
     // get all verts coords into tag; don't need to worry about filtering out fixed verts, 
     // we'll just be setting to their fixed coords
   std::vector<double> vcentroids(3*verts.size());
-  rval = mbImpl->get_coords(verts, &vcentroids[0]); RR("Failed to get vert coords.");
+  if (!coordsTag) {
+    rval = mbImpl->get_coords(verts, &vcentroids[0]); RR("Failed to get vert coords.");
+  }
+  else {
+    rval = mbImpl->tag_get_data(coordsTag, verts, &vcentroids[0]); RR("Failed to get vert coords tag values.");
+  }
 
   Tag centroid;
   rval = mbImpl->tag_get_handle("", 3, MB_TYPE_DOUBLE, centroid, MB_TAG_CREAT | MB_TAG_DENSE); 
@@ -174,7 +179,13 @@ ErrorCode LloydSmoother::perform_smooth()
   } // end while
   
     // write the tag back onto vertex coordinates
-  rval = mbImpl->set_coords(owned_verts, &vcentroids[0]); RR("Failed to set vertex coords.");
+  if (!coordsTag) {
+    rval = mbImpl->set_coords(owned_verts, &vcentroids[0]); RR("Failed to set vertex coords.");
+  }
+  else {
+    rval = mbImpl->tag_set_data(coordsTag, owned_verts, &vcentroids[0]); 
+    RR("Failed to set vert coords tag values.");
+  }
 
   return MB_SUCCESS;
 }

diff --git a/src/moab/LloydSmoother.hpp b/src/moab/LloydSmoother.hpp
index 5ab9771..4f8570a 100644
--- a/src/moab/LloydSmoother.hpp
+++ b/src/moab/LloydSmoother.hpp
@@ -35,13 +35,15 @@ public:
     /* \brief Constructor
      * Convenience constructor, data input directly
      * \param impl The MOAB instance for this smoother
-     * \param pcomm The ParallelComm instance by which this mesh is parallel
+     * \param pc The ParallelComm instance by which this mesh is parallel
      * \param elems The mesh to be smoothed
+     * \param cds_tag If specified, this tag is used to get/set coordinates, rather than 
+     *     true vertex coordinates
      * \param fixed_tag The tag marking which vertices are fixed
      * \param abs_tol Absolute tolerance measuring convergence
      * \param rel_tol Relative tolerance measuring convergence
      */
-  LloydSmoother(Interface *impl, ParallelComm *pc, Range &elems, Tag fixed_tag = 0,
+  LloydSmoother(Interface *impl, ParallelComm *pc, Range &elems, Tag cds_tag = 0, Tag fixed_tag = 0,
                 double abs_tol = -1.0, double rel_tol = 1.0e-6);
 
     /* \brief Destructor
@@ -80,6 +82,14 @@ public:
      */
   void fixed_tag(Tag fixed) {fixedTag = fixed;}
 
+    /* \brief get/set coords tag
+     */
+  Tag coords_tag() {return coordsTag;}
+
+    /* \brief get/set coords tag
+     */
+  void coords_tag(Tag coords) {coordsTag = coords;}
+
     /* \brief get/set tolerance
      */
   double abs_tol() {return absTol;}
@@ -121,6 +131,9 @@ private:
     //- elements to smooth
   Range myElems;
   
+    //- tag for coordinates; if zero, true vertex coords are used
+  Tag coordsTag;
+  
     //- tag marking which vertices are fixed, 0 = not fixed, otherwise fixed
   Tag fixedTag;
   

diff --git a/test/lloyd_smoother_test.cpp b/test/lloyd_smoother_test.cpp
index 8d88112..662732b 100644
--- a/test/lloyd_smoother_test.cpp
+++ b/test/lloyd_smoother_test.cpp
@@ -1,5 +1,6 @@
 #include "moab/Core.hpp"
 #include "moab/LloydSmoother.hpp"
+#include "moab/CartVect.hpp"
 #include "TestUtil.hpp"
 
 #ifdef MESHDIR
@@ -31,11 +32,47 @@ int main(int argc, char**argv)
     CHECK_ERR(MB_FAILURE);
   }
 
-  LloydSmoother ll(&mb, NULL, elems);
+    // get the vertex positions and set on an intermediate tag, to test smoothing for
+    // a tag instead of for coords
+  Tag ctag;
+  rval = mb.tag_get_handle("vcentroid", 3, MB_TYPE_DOUBLE, ctag, MB_TAG_CREAT|MB_TAG_DENSE);
+  CHECK_ERR(rval);
+  Range verts;
+  rval = mb.get_entities_by_dimension(0, 0, verts);
+  CHECK_ERR(rval);
+  std::vector<double> coords(3*verts.size());
+  rval = mb.get_coords(verts, &coords[0]);
+  CHECK_ERR(rval);
+  rval = mb.tag_set_data(ctag, verts, &coords[0]);
+  CHECK_ERR(rval);
+
+  LloydSmoother ll(&mb, NULL, elems, ctag);
   ll.report_its(10);
   rval = ll.perform_smooth();
   CHECK_ERR(rval);
-  
   std::cout << "Mesh smoothed in " << ll.num_its() << " iterations." << std::endl;
+
+    // now, set vertex coords to almost their converged positions, then re-smooth; should take fewer
+    // iterations
+  std::vector<double> new_coords(3*verts.size());
+  rval = mb.tag_get_data(ctag, verts, &new_coords[0]);
+  CHECK_ERR(rval);
+  unsigned int i;
+  Range::iterator vit;
+  for (vit = verts.begin(), i = 0; vit != verts.end(); vit++, i+=3) {
+    CartVect old_pos(&coords[i]), new_pos(&new_coords[i]);
+    CartVect almost_pos = old_pos + .99 * (new_pos - old_pos);
+    almost_pos.get(&new_coords[i]);
+  }
+  rval = mb.set_coords(verts, &new_coords[0]);
+  CHECK_ERR(rval);
+
+  LloydSmoother ll2(&mb, NULL, elems);
+  ll2.report_its(10);
+  rval = ll2.perform_smooth();
+  CHECK_ERR(rval);
+  std::cout << "Mesh smoothed in " << ll2.num_its() << " iterations." << std::endl;
+  
+  
 }
 


https://bitbucket.org/fathomteam/moab/commits/0dd47d3c7a84/
Changeset:   0dd47d3c7a84
Branch:      None
User:        tautges
Date:        2013-11-28 01:33:44
Summary:     Adding deformed mesh example.

Affected #:  2 files

diff --git a/examples/DeformMeshRemap.cpp b/examples/DeformMeshRemap.cpp
new file mode 100644
index 0000000..ddb3f3a
--- /dev/null
+++ b/examples/DeformMeshRemap.cpp
@@ -0,0 +1,180 @@
+/** @example DeformMeshRemap.cpp
+ * Description: Account for mesh deformation of a solid due to structural mechanics\n
+ * In this example there are two meshes, a "master" and "slave" mesh.  In the master mesh,
+ * the solid material is deformed, to mimic what happens when a solid heats up and deforms.
+ * The fluid mesh is smoothed to account for those deformations, and tags on the fluid are
+ * remapped to those new positions.  Then mesh positions and state variables are transferred
+ * to the slave mesh, mimicing another mesh used by some other physics.
+ *
+ * To run: ./DeformMeshRemap [<master_meshfile><slave_meshfile>]\n
+ * (default values can run if users don't specify the mesh files)
+ */
+
+#include "moab/Core.hpp"
+#include "moab/Range.hpp"
+#include "moab/LloydSmoother.hpp"
+#include "MBTagConventions.hpp"
+
+#include <iostream>
+#include <assert.h>
+
+using namespace moab;
+using namespace std;
+
+#ifndef MESH_DIR
+#define MESH_DIR "."
+#endif
+
+ErrorCode read_file(string &fname, EntityHandle &seth, 
+                    Range &solids, Range &solid_elems, Range &fluids, Range &fluid_elems);
+void deform_func(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);
+
+string trimeshf = string(MESH_DIR) + string("/rodtri.g");
+string quadmeshf = string(MESH_DIR) + string("/rodquad.g");
+const int SOLID_SETNO = 100, FLUID_SETNO = 200;
+
+Interface *mb;
+#define RR(a) if (MB_SUCCESS != rval) {cout << a << endl; return MB_FAILURE;}
+    
+
+int main(int /*argc*/, char **/*argv*/) {
+
+  EntityHandle master, slave;
+  ErrorCode rval;
+
+  mb = new Core();
+  
+    // read master/slave files and get fluid/solid material sets
+  Range fluids[2], solids[2], solid_elems[2], fluid_elems[2];
+  rval = read_file(quadmeshf, master, solids[0], solid_elems[0], fluids[0], fluid_elems[0]); RR("");
+  rval = read_file(trimeshf, slave, solids[1], solid_elems[1], fluids[1], fluid_elems[1]); RR("");
+
+    // deform the master's solid mesh, put results in a new tag
+  Tag xnew;
+  rval = deform_master(fluid_elems[0], solid_elems[0], xnew); RR("");
+  rval = write_to_coords(solid_elems[0], xnew); RR("");
+  rval = mb->write_file("deformed.vtk", NULL, NULL, &master, 1); RR("");
+  
+    // smooth the master mesh
+  LloydSmoother *ll = new LloydSmoother(mb, NULL, fluid_elems[0], xnew);
+  rval = ll->perform_smooth();
+  RR("Failed in lloyd smoothing.");
+  cout << "Lloyd smoothing required " << ll->num_its() << " iterations." << endl;
+
+  rval = mb->write_file("smoothed.vtk", NULL, NULL, &master, 1); RR("");
+  
+  delete ll;
+  delete mb;
+  
+  return MB_SUCCESS;
+}
+
+ErrorCode write_to_coords(Range &elems, Tag tagh) 
+{
+    // write the tag to coordinates
+  Range verts;
+  ErrorCode rval = mb->get_adjacencies(elems, 0, false, verts, Interface::UNION);
+  RR("Failed to get adj vertices.");
+  std::vector<double> coords(3*verts.size());
+  rval = mb->tag_get_data(tagh, verts, &coords[0]);
+  RR("Failed to get tag data.");
+  rval = mb->set_coords(verts, &coords[0]);
+  RR("Failed to set coordinates.");
+  return MB_SUCCESS;
+}
+
+void deform_func(double *xold, double *xnew) 
+{
+  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;
+  
+  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;
+}
+  
+ErrorCode deform_master(Range &fluid_elems, Range &solid_elems, Tag &xnew) 
+{
+    // deform elements with an analytic function
+
+    // create the tag
+  ErrorCode rval = mb->tag_get_handle("", 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
+  Range verts;
+  rval = mb->get_adjacencies(fluid_elems, 0, false, verts, Interface::UNION);
+  RR("Failed to get vertices.");
+  std::vector<double> coords(3*verts.size(), 0.0);
+  rval = mb->get_coords(verts, &coords[0]);
+  RR("Failed to get vertex coords.");
+  rval = mb->tag_set_data(xnew, verts, &coords[0]);
+  RR("Failed to set xnew tag on fluid verts.");
+  
+    // get all the vertices and coords in the solid
+  verts.clear();
+  rval = mb->get_adjacencies(solid_elems, 0, false, verts, Interface::UNION);
+  RR("Failed to get vertices.");
+  coords.resize(3*verts.size(), 0.0);
+  rval = mb->get_coords(verts, &coords[0]);
+  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]);
+    
+    // set the new tag to those coords
+  rval = mb->tag_set_data(xnew, verts, &coords[0]);
+  RR("Failed to set tag data.");
+  
+  return MB_SUCCESS;
+}
+
+ErrorCode read_file(string &fname, EntityHandle &seth, 
+                    Range &solids, Range &solid_elems, Range &fluids, Range &fluid_elems)
+{
+    // create meshset
+  ErrorCode rval = mb->create_meshset(0, seth);
+  RR("Couldn't create master/slave set.");
+  rval = mb->load_file(fname.c_str(), &seth);
+  RR("Couldn't load master/slave mesh.");
+
+    // get material sets for solid/fluid
+  Tag tagh;
+  rval = mb->tag_get_handle(MATERIAL_SET_TAG_NAME, tagh); RR("Couldn't get material set tag name.");
+  const void *setno_ptr = &SOLID_SETNO;
+  rval = mb->get_entities_by_type_and_tag(seth, MBENTITYSET, &tagh, &setno_ptr, 1, solids);
+  if (solids.empty()) rval = MB_FAILURE;
+  RR("Couldn't get any solid sets.");
+
+    // get solid entities, and dimension
+  Range tmp_range;
+  for (Range::iterator rit = solids.begin(); rit != solids.end(); rit++) {
+    rval = mb->get_entities_by_handle(*rit, tmp_range, true);
+    RR("Failed to get entities in solid.");
+  }
+  int dim = mb->dimension_from_handle(*tmp_range.rbegin());
+  assert(dim > 0 && dim < 4);
+  
+  solid_elems = tmp_range.subset_by_dimension(dim);
+
+  setno_ptr = &FLUID_SETNO;
+  rval = mb->get_entities_by_type_and_tag(seth, MBENTITYSET, &tagh, &setno_ptr, 1, fluids);
+  if (fluids.empty()) rval = MB_FAILURE;
+  RR("Couldn't get any fluid sets.");
+  
+  for (Range::iterator rit = fluids.begin(); rit != fluids.end(); rit++) {
+    rval = mb->get_entities_by_dimension(*rit, dim, fluid_elems, true);
+    RR("Failed to get entities in fluid.");
+  }
+  if (mb->dimension_from_handle(*fluid_elems.begin()) != dim) {
+    rval = MB_FAILURE;
+    RR("Fluid and solid elements must be same dimension.");
+  }
+  
+  return MB_SUCCESS;
+}

diff --git a/examples/makefile b/examples/makefile
index 5ae5af0..50643cb 100644
--- a/examples/makefile
+++ b/examples/makefile
@@ -4,7 +4,7 @@ include ${MOAB_DIR}/lib/iMesh-Defs.inc
 
 .SUFFIXES: .o .cpp .F90
 
-# MESH_DIR is the top-level MOAB source directory, used to locate mesh files that come with MOAB source
+# MESH_DIR is the directory containing mesh files that come with MOAB source
 MESH_DIR="../MeshFiles/unittest"
 
 EXAMPLES = HelloMOAB GetEntities SetsNTags StructuredMeshSimple DirectAccessWithHoles DirectAccessNoHoles 
@@ -53,6 +53,9 @@ point_in_elem_search: point_in_elem_search.o
 PushParMeshIntoMoabF90: PushParMeshIntoMoabF90.o
 	${MOAB_CXX} -o $@ $< ${IMESH_LIBS} -lgfortran -L/usr/lib/openmpi/lib -lmpi_f90 -lmpi_f77 -lmpi -lopen-rte -lopen-pal -ldl 
 
+DeformMeshRemap: DeformMeshRemap.o
+	${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK}
+
 clean:
 	rm -rf *.o ${EXAMPLES} ${PAREXAMPLES} ${EXOIIEXAMPLES}
 


https://bitbucket.org/fathomteam/moab/commits/0e1d2f088c4f/
Changeset:   0e1d2f088c4f
Branch:      None
User:        tautges
Date:        2013-11-28 01:33:44
Summary:     Adding option processing to DeformMeshRemap, and making
parallel optional in LloydRelaxation example.

Affected #:  2 files

diff --git a/examples/DeformMeshRemap.cpp b/examples/DeformMeshRemap.cpp
index ddb3f3a..55ac3cc 100644
--- a/examples/DeformMeshRemap.cpp
+++ b/examples/DeformMeshRemap.cpp
@@ -13,6 +13,7 @@
 #include "moab/Core.hpp"
 #include "moab/Range.hpp"
 #include "moab/LloydSmoother.hpp"
+#include "moab/ProgOptions.hpp"
 #include "MBTagConventions.hpp"
 
 #include <iostream>
@@ -32,39 +33,48 @@ 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);
 
-string trimeshf = string(MESH_DIR) + string("/rodtri.g");
-string quadmeshf = string(MESH_DIR) + string("/rodquad.g");
 const int SOLID_SETNO = 100, FLUID_SETNO = 200;
 
 Interface *mb;
 #define RR(a) if (MB_SUCCESS != rval) {cout << a << endl; return MB_FAILURE;}
     
 
-int main(int /*argc*/, char **/*argv*/) {
+int main(int argc, char **argv) {
 
   EntityHandle master, slave;
   ErrorCode rval;
 
+  ProgOptions po("Deformed mesh options");
+  po.addOpt<std::string> ("master,m", "Specify the master meshfile name" );
+  po.addOpt<std::string> ("slave,s", "Specify the slave meshfile name" );
+  po.parseCommandLine(argc, argv);
+  std::string foo;
+  string masterf, slavef;
+  if(!po.getOpt("master", &masterf))
+    masterf = string(MESH_DIR) + string("/rodquad.g");
+  if(!po.getOpt("slave", &slavef))
+    slavef = string(MESH_DIR) + string("/rodtri.g");
+
   mb = new Core();
   
     // read master/slave files and get fluid/solid material sets
   Range fluids[2], solids[2], solid_elems[2], fluid_elems[2];
-  rval = read_file(quadmeshf, master, solids[0], solid_elems[0], fluids[0], fluid_elems[0]); RR("");
-  rval = read_file(trimeshf, slave, solids[1], solid_elems[1], fluids[1], fluid_elems[1]); RR("");
+  rval = read_file(masterf, master, solids[0], solid_elems[0], fluids[0], fluid_elems[0]); RR("");
+  rval = read_file(slavef, slave, solids[1], solid_elems[1], fluids[1], fluid_elems[1]); RR("");
 
     // deform the master's solid mesh, put results in a new tag
   Tag xnew;
   rval = deform_master(fluid_elems[0], solid_elems[0], xnew); RR("");
-  rval = write_to_coords(solid_elems[0], xnew); RR("");
-  rval = mb->write_file("deformed.vtk", NULL, NULL, &master, 1); RR("");
+  if (debug) write_and_save(solid_elems[0], master, xnew, "deformed.vtk");
   
     // smooth the master mesh
   LloydSmoother *ll = new LloydSmoother(mb, NULL, fluid_elems[0], xnew);
   rval = ll->perform_smooth();
   RR("Failed in lloyd smoothing.");
   cout << "Lloyd smoothing required " << ll->num_its() << " iterations." << endl;
+  if (debug) write_and_save(fluid_elems[0], master, xnew, "smoothed.vtk");
 
-  rval = mb->write_file("smoothed.vtk", NULL, NULL, &master, 1); RR("");
+    // map new locations to slave
   
   delete ll;
   delete mb;
@@ -72,6 +82,13 @@ int main(int /*argc*/, char **/*argv*/) {
   return MB_SUCCESS;
 }
 
+ErrorCode write_and_save(Range &ents, EntithHanlde seth, Tag tagh, const char *filename) 
+{
+  rval = write_to_coords(ents, tagh); RR("");
+  rval = mb->write_file("deformed.vtk", NULL, NULL, &seth, 1); RR("");
+  return rval;
+}
+  
 ErrorCode write_to_coords(Range &elems, Tag tagh) 
 {
     // write the tag to coordinates

diff --git a/examples/LloydRelaxation.cpp b/examples/LloydRelaxation.cpp
index 0cf9d1e..34002cd 100644
--- a/examples/LloydRelaxation.cpp
+++ b/examples/LloydRelaxation.cpp
@@ -14,8 +14,10 @@
  * in the current directory (H5M format must be used since the file is written in parallel).
  */
 
-#include "moab/ParallelComm.hpp"
-#include "MBParallelConventions.h"
+#ifdef USE_MPI
+#  include "moab/ParallelComm.hpp"
+#  include "MBParallelConventions.h"
+#endif
 #include "moab/Core.hpp"
 #include "moab/Skinner.hpp"
 #include "moab/CN.hpp"
@@ -30,7 +32,7 @@ string test_file_name = string(MESH_DIR) + string("/surfrandomtris-4part.h5m");
 
 #define RC if (MB_SUCCESS != rval) return rval
 
-ErrorCode perform_lloyd_relaxation(ParallelComm *pc, Range &verts, Range &cells, Tag fixed, 
+ErrorCode perform_lloyd_relaxation(Interface *mb, Range &verts, Range &cells, Tag fixed, 
                                    int num_its, int report_its);
 
 int main(int argc, char **argv)
@@ -38,7 +40,9 @@ int main(int argc, char **argv)
   int num_its = 10;
   int report_its = 1;
 
+#ifdef USE_MPI
   MPI_Init(&argc, &argv);
+#endif
 
     // need option handling here for input filename
   if (argc > 1){
@@ -49,10 +53,13 @@ int main(int argc, char **argv)
   // get MOAB and ParallelComm instances
   Interface *mb = new Core;
   if (NULL == mb) return 1;
+  int nprocs = 1;
+  
+#ifdef USE_MPI
   // get the ParallelComm instance
-  ParallelComm* pcomm = new ParallelComm(mb, MPI_COMM_WORLD);
-  int nprocs = pcomm->size();
-
+  ParallelComm *pcomm = new ParallelComm(mb, MPI_COMM_WORLD);
+  nprocs = pcomm->size();
+#endif
   string options;
   if (nprocs > 1) // if reading in parallel, need to tell it how
     options = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_SHARED_ENTS;PARALLEL_GHOSTS=2.0.1;DEBUG_IO=0;DEBUG_PIO=0";
@@ -76,34 +83,45 @@ int main(int argc, char **argv)
     // ok to mark non-owned skin vertices too, I won't move those anyway
     // use MOAB's skinner class to find the skin
   Skinner skinner(mb);
-  rval = skinner.find_skin(faces, true, skin_verts); RC; // 'true' param indicates we want vertices back, not faces
+  rval = skinner.find_skin(0, faces, true, skin_verts); RC; // 'true' param indicates we want vertices back, not faces
 
   std::vector<int> fix_tag(skin_verts.size(), 1); // initialized to 1 to indicate fixed
   rval = mb->tag_set_data(fixed, skin_verts, &fix_tag[0]); RC;
 
     // now perform the Lloyd relaxation
-  rval = perform_lloyd_relaxation(pcomm, verts, faces, fixed, num_its, report_its); RC;
+  rval = perform_lloyd_relaxation(mb, verts, faces, fixed, num_its, report_its); RC;
 
     // delete fixed tag, since we created it here
   rval = mb->tag_delete(fixed); RC;
   
     // output file, using parallel write
-  rval = mb->write_file("lloydfinal.h5m", NULL, "PARALLEL=WRITE_PART"); RC;
+
+#ifdef USE_MPI
+  options = "PARALLEL=WRITE_PART";
+#endif
+
+  rval = mb->write_file("lloydfinal.h5m", NULL, options.c_str()); RC;
 
     // delete MOAB instance
   delete mb;
   
+#ifdef USE_MPI
   MPI_Finalize();
+#endif
 
   return 0;
 }
 
-ErrorCode perform_lloyd_relaxation(ParallelComm *pcomm, Range &verts, Range &faces, Tag fixed, 
+ErrorCode perform_lloyd_relaxation(Interface *mb, Range &verts, Range &faces, Tag fixed, 
                                    int num_its, int report_its) 
 {
   ErrorCode rval;
-  Interface *mb = pcomm->get_moab();
-  int nprocs = pcomm->size();
+  int nprocs = 1;
+
+#ifdef USE_MPI
+  ParallelComm *pcomm = ParallelComm::get_pcomm(mb, 0);
+  nprocs = pcomm->size();
+#endif
   
     // perform Lloyd relaxation:
     // 1. setup: set vertex centroids from vertex coords; filter to owned verts; get fixed tags
@@ -120,8 +138,10 @@ ErrorCode perform_lloyd_relaxation(ParallelComm *pcomm, Range &verts, Range &fac
     // filter verts down to owned ones and get fixed tag for them
   Range owned_verts, shared_owned_verts;
   if (nprocs > 1) {
+#ifdef USE_MPI
     rval = pcomm->filter_pstatus(verts, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, &owned_verts);
     if (rval != MB_SUCCESS) return rval;
+#endif
   }
   else
     owned_verts = verts;
@@ -133,11 +153,13 @@ ErrorCode perform_lloyd_relaxation(ParallelComm *pcomm, Range &verts, Range &fac
   vcentroids.resize(3*owned_verts.size());
   rval = mb->tag_get_data(centroid, owned_verts, &vcentroids[0]); RC;
 
+#ifdef USE_MPI
     // get shared owned verts, for exchanging tags
   rval = pcomm->get_shared_entities(-1, shared_owned_verts, 0, false, true); RC;
     // workaround: if no shared owned verts, put a non-shared one in the list, to prevent exchanging tags
     // for all shared entities
   if (shared_owned_verts.empty()) shared_owned_verts.insert(*verts.begin());
+#endif
   
     // some declarations for later iterations
   std::vector<double> fcentroids(3*faces.size()); // fcentroids for face centroids
@@ -195,16 +217,22 @@ ErrorCode perform_lloyd_relaxation(ParallelComm *pcomm, Range &verts, Range &fac
 
     // 2c. exchange tags on owned verts
     if (nprocs > 1) {
+#ifdef USE_MPI
       rval = pcomm->exchange_tags(centroid, shared_owned_verts); RC;
+#endif
     }
 
 
     if (!(nit%report_its)) {
         // global reduce for maximum delta, then report it
       double global_max = mxdelta;
+      int myrank = 0;
+#ifdef USE_MPI
       if (nprocs > 1)
         MPI_Reduce(&mxdelta, &global_max, 1, MPI_DOUBLE, MPI_MAX, 0, pcomm->comm());
-      if (1 == nprocs || !pcomm->rank()) 
+      myrank = pcomm->rank();
+#endif
+      if (1 == nprocs || !myrank) 
         cout << "Max delta = " << global_max << endl;
     }
   }


https://bitbucket.org/fathomteam/moab/commits/91e2c8e9c401/
Changeset:   91e2c8e9c401
Branch:      None
User:        tautges
Date:        2013-11-28 01:33:44
Summary:     Make a DeformMeshRemap class and implement through that.

Affected #:  1 file

diff --git a/examples/DeformMeshRemap.cpp b/examples/DeformMeshRemap.cpp
index 55ac3cc..fac83ae 100644
--- a/examples/DeformMeshRemap.cpp
+++ b/examples/DeformMeshRemap.cpp
@@ -17,6 +17,7 @@
 #include "MBTagConventions.hpp"
 
 #include <iostream>
+#include <set>
 #include <assert.h>
 
 using namespace moab;
@@ -37,11 +38,212 @@ const int SOLID_SETNO = 100, FLUID_SETNO = 200;
 
 Interface *mb;
 #define RR(a) if (MB_SUCCESS != rval) {cout << a << endl; return MB_FAILURE;}
-    
+
+const bool debug = true;
+
+class DeformMeshRemap 
+{
+public:
+
+    //! enumerator for solid/fluid, master/slave
+  enum {MASTER=0, SLAVE, SOLID, FLUID};
+  
+    //! constructor
+  DeformMeshRemap(Interface *impl) : mbImpl(impl), masterSet(0), slaveSet(0), xNew(0), xNewName("xnew") {}
+  
+    //! destructor
+  ~DeformMeshRemap();
+
+    //! execute the deformed mesh process
+  ErrorCode execute();
+  
+    //! add a set number
+  ErrorCode add_set_no(int fluid_or_solid, int set_no);
+  
+    //! remove a set number
+  ErrorCode remove_set_no(int fluid_or_solid, int set_no);
+  
+    //! get the set numbers
+  ErrorCode get_set_nos(int fluid_or_solid, std::set<int> &set_nos) const;
+
+    //! get the xNew tag handle
+  inline Tag x_new() const {return xNew;}
+
+    //! get the tag name
+  std::string x_new_name() const {return xNewName;}
+  
+    //! set the tag name
+  void x_new_name(const std::string &name) {xNewName = name;}
+
+    //! get/set the file name
+  std::string get_file_name(int m_or_s) const;
+  
+    //! get/set the file name
+  void set_file_name(int m_or_s, const std::string &name);
+  
+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);
+
+    //! read a file and establish proper ranges
+  ErrorCode read_file(int m_or_s, string &fname, EntityHandle &seth);
+
+    //! write the input tag to the coordinates for the vertices in the input elems
+  ErrorCode write_to_coords(Range &elems, Tag tagh);
+
+    //! write the tag to the vertices, then save to the specified file
+  ErrorCode write_and_save(Range &ents, EntityHandle seth, Tag tagh, const char *filename);
+
+    //! moab interface
+  Interface *mbImpl;
+  
+    //! material set numbers for fluid materials
+  std::set<int> fluidSetNos;
+
+    //! material set numbers for solid materials
+  std::set<int> solidSetNos;
+
+    //! sets defining master/slave meshes
+  EntityHandle masterSet, slaveSet;
+
+    //! sets in master/slave meshes
+  Range fluidSets[2], solidSets[2];
+  
+    //! elements in master/slave meshes
+  Range fluidElems[2], solidElems[2];
+  
+    //! filenames for master/slave meshes
+  std::string masterFileName, slaveFileName;
+
+    //! tag used for new positions
+  Tag xNew;
+  
+    //! tag name used for new positions
+  std::string xNewName;
+};
+
+  //! add a set number
+inline ErrorCode DeformMeshRemap::add_set_no(int f_or_s, int set_no) 
+{
+  std::set<int> *this_set;
+  switch (f_or_s) {
+    case FLUID:
+        this_set = &fluidSetNos; break;
+    case SOLID:
+        this_set = &solidSetNos; break;
+    default:
+        assert(false && "f_or_s should be FLUID or SOLID.");
+        return MB_FAILURE;
+  }
+
+  this_set->insert(set_no);
+  
+  return MB_SUCCESS;
+}
+  
+  //! remove a set number
+inline ErrorCode DeformMeshRemap::remove_set_no(int f_or_s, int set_no) 
+{
+  std::set<int> *this_set;
+  switch (f_or_s) {
+    case FLUID:
+        this_set = &fluidSetNos; break;
+    case SOLID:
+        this_set = &solidSetNos; break;
+    default:
+        assert(false && "f_or_s should be FLUID or SOLID.");
+        return MB_FAILURE;
+  }
+  std::set<int>::iterator sit = this_set->find(set_no);
+  if (sit != this_set->end()) {
+    this_set->erase(*sit);
+    return MB_SUCCESS;
+  }
+
+  return MB_FAILURE;
+}
+  
+  //! get the set numbers
+inline ErrorCode DeformMeshRemap::get_set_nos(int f_or_s, std::set<int> &set_nos) const
+{
+  const std::set<int> *this_set;
+  switch (f_or_s) {
+    case FLUID:
+        this_set = &fluidSetNos; break;
+    case SOLID:
+        this_set = &solidSetNos; break;
+    default:
+        assert(false && "f_or_s should be FLUID or SOLID.");
+        return MB_FAILURE;
+  }
+
+  set_nos = *this_set;
+  
+  return MB_SUCCESS;
+}
+
+ErrorCode DeformMeshRemap::execute() 
+{
+    // read master/slave files and get fluid/solid material sets
+  ErrorCode rval = read_file(MASTER, masterFileName, masterSet);
+  if (MB_SUCCESS != rval) return rval;
+  
+  rval = read_file(SLAVE, slaveFileName, slaveSet);
+  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("");
+  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");
+
+    // map new locations to slave
+  
+  delete ll;
+
+  return MB_SUCCESS;
+}
+
+std::string DeformMeshRemap::get_file_name(int m_or_s) const
+{
+  switch (m_or_s) {
+    case MASTER:
+        return masterFileName;
+    case SLAVE:
+        return slaveFileName;
+    default:
+        assert(false && "m_or_s should be MASTER or SLAVE.");
+        return std::string();
+  }
+}
+  
+void DeformMeshRemap::set_file_name(int m_or_s, const std::string &name) 
+{
+  switch (m_or_s) {
+    case MASTER:
+        masterFileName = name; break;
+    case SLAVE:
+        slaveFileName = name; break;
+    default:
+        assert(false && "m_or_s should be MASTER or SLAVE.");
+  }
+}
+
+DeformMeshRemap::~DeformMeshRemap() 
+{
+    // delete the tag
+  mbImpl->tag_delete(xNew);
+}
 
 int main(int argc, char **argv) {
 
-  EntityHandle master, slave;
   ErrorCode rval;
 
   ProgOptions po("Deformed mesh options");
@@ -56,49 +258,34 @@ int main(int argc, char **argv) {
     slavef = string(MESH_DIR) + string("/rodtri.g");
 
   mb = new Core();
-  
-    // read master/slave files and get fluid/solid material sets
-  Range fluids[2], solids[2], solid_elems[2], fluid_elems[2];
-  rval = read_file(masterf, master, solids[0], solid_elems[0], fluids[0], fluid_elems[0]); RR("");
-  rval = read_file(slavef, slave, solids[1], solid_elems[1], fluids[1], fluid_elems[1]); RR("");
 
-    // deform the master's solid mesh, put results in a new tag
-  Tag xnew;
-  rval = deform_master(fluid_elems[0], solid_elems[0], xnew); RR("");
-  if (debug) write_and_save(solid_elems[0], master, xnew, "deformed.vtk");
+  DeformMeshRemap dfr(mb);
+  dfr.set_file_name(DeformMeshRemap::MASTER, masterf);
+  dfr.set_file_name(DeformMeshRemap::SLAVE, slavef);
+  rval = dfr.add_set_no(DeformMeshRemap::SOLID, SOLID_SETNO); RR("Failed to add solid set no.");
+  rval = dfr.add_set_no(DeformMeshRemap::FLUID, FLUID_SETNO); RR("Failed to add fluid set no.");
   
-    // smooth the master mesh
-  LloydSmoother *ll = new LloydSmoother(mb, NULL, fluid_elems[0], xnew);
-  rval = ll->perform_smooth();
-  RR("Failed in lloyd smoothing.");
-  cout << "Lloyd smoothing required " << ll->num_its() << " iterations." << endl;
-  if (debug) write_and_save(fluid_elems[0], master, xnew, "smoothed.vtk");
-
-    // map new locations to slave
-  
-  delete ll;
-  delete mb;
-  
-  return MB_SUCCESS;
+  rval = dfr.execute();
+  return rval;
 }
 
-ErrorCode write_and_save(Range &ents, EntithHanlde seth, Tag tagh, const char *filename) 
+ErrorCode DeformMeshRemap::write_and_save(Range &ents, EntityHandle seth, Tag tagh, const char *filename) 
 {
-  rval = write_to_coords(ents, tagh); RR("");
-  rval = mb->write_file("deformed.vtk", NULL, NULL, &seth, 1); RR("");
+  ErrorCode rval = write_to_coords(ents, tagh); RR("");
+  rval = mbImpl->write_file(filename, NULL, NULL, &seth, 1); RR("");
   return rval;
 }
   
-ErrorCode write_to_coords(Range &elems, Tag tagh) 
+ErrorCode DeformMeshRemap::write_to_coords(Range &elems, Tag tagh) 
 {
     // write the tag to coordinates
   Range verts;
-  ErrorCode rval = mb->get_adjacencies(elems, 0, false, verts, Interface::UNION);
+  ErrorCode rval = mbImpl->get_adjacencies(elems, 0, false, verts, Interface::UNION);
   RR("Failed to get adj vertices.");
   std::vector<double> coords(3*verts.size());
-  rval = mb->tag_get_data(tagh, verts, &coords[0]);
+  rval = mbImpl->tag_get_data(tagh, verts, &coords[0]);
   RR("Failed to get tag data.");
-  rval = mb->set_coords(verts, &coords[0]);
+  rval = mbImpl->set_coords(verts, &coords[0]);
   RR("Failed to set coordinates.");
   return MB_SUCCESS;
 }
@@ -115,83 +302,92 @@ void deform_func(double *xold, double *xnew)
   xnew[1] = xold[1] + yfrac * (1.0 + xfrac) * 0.05;
 }
   
-ErrorCode deform_master(Range &fluid_elems, Range &solid_elems, Tag &xnew) 
+ErrorCode DeformMeshRemap::deform_master(Range &fluid_elems, Range &solid_elems) 
 {
     // deform elements with an analytic function
 
     // create the tag
-  ErrorCode rval = mb->tag_get_handle("", 3, MB_TYPE_DOUBLE, xnew, MB_TAG_CREAT|MB_TAG_DENSE);
+  ErrorCode rval = mbImpl->tag_get_handle("", 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
   Range verts;
-  rval = mb->get_adjacencies(fluid_elems, 0, false, verts, Interface::UNION);
+  rval = mbImpl->get_adjacencies(fluid_elems, 0, false, verts, Interface::UNION);
   RR("Failed to get vertices.");
   std::vector<double> coords(3*verts.size(), 0.0);
-  rval = mb->get_coords(verts, &coords[0]);
+  rval = mbImpl->get_coords(verts, &coords[0]);
   RR("Failed to get vertex coords.");
-  rval = mb->tag_set_data(xnew, verts, &coords[0]);
+  rval = mbImpl->tag_set_data(xNew, verts, &coords[0]);
   RR("Failed to set xnew tag on fluid verts.");
   
     // get all the vertices and coords in the solid
   verts.clear();
-  rval = mb->get_adjacencies(solid_elems, 0, false, verts, Interface::UNION);
+  rval = mbImpl->get_adjacencies(solid_elems, 0, false, verts, Interface::UNION);
   RR("Failed to get vertices.");
   coords.resize(3*verts.size(), 0.0);
-  rval = mb->get_coords(verts, &coords[0]);
+  rval = mbImpl->get_coords(verts, &coords[0]);
   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]);
     
     // set the new tag to those coords
-  rval = mb->tag_set_data(xnew, verts, &coords[0]);
+  rval = mbImpl->tag_set_data(xNew, verts, &coords[0]);
   RR("Failed to set tag data.");
   
   return MB_SUCCESS;
 }
 
-ErrorCode read_file(string &fname, EntityHandle &seth, 
-                    Range &solids, Range &solid_elems, Range &fluids, Range &fluid_elems)
+ErrorCode DeformMeshRemap::read_file(int m_or_s, string &fname, EntityHandle &seth)
 {
     // create meshset
-  ErrorCode rval = mb->create_meshset(0, seth);
+  ErrorCode rval = mbImpl->create_meshset(0, seth);
   RR("Couldn't create master/slave set.");
-  rval = mb->load_file(fname.c_str(), &seth);
+  rval = mbImpl->load_file(fname.c_str(), &seth);
   RR("Couldn't load master/slave mesh.");
 
     // get material sets for solid/fluid
   Tag tagh;
-  rval = mb->tag_get_handle(MATERIAL_SET_TAG_NAME, tagh); RR("Couldn't get material set tag name.");
-  const void *setno_ptr = &SOLID_SETNO;
-  rval = mb->get_entities_by_type_and_tag(seth, MBENTITYSET, &tagh, &setno_ptr, 1, solids);
-  if (solids.empty()) rval = MB_FAILURE;
-  RR("Couldn't get any solid sets.");
+  rval = mbImpl->tag_get_handle(MATERIAL_SET_TAG_NAME, tagh); RR("Couldn't get material set tag name.");
+  for (std::set<int>::iterator sit = solidSetNos.begin(); sit != solidSetNos.end(); sit++) {
+    Range sets;
+    int set_no = *sit;
+    const void *setno_ptr = &set_no;
+    rval = mbImpl->get_entities_by_type_and_tag(seth, MBENTITYSET, &tagh, &setno_ptr, 1, sets);
+    if (sets.empty()) rval = MB_FAILURE;
+    RR("Couldn't get any solid sets.");
+    solidSets[m_or_s].merge(sets);
+  }
 
     // get solid entities, and dimension
   Range tmp_range;
-  for (Range::iterator rit = solids.begin(); rit != solids.end(); rit++) {
-    rval = mb->get_entities_by_handle(*rit, tmp_range, true);
+  for (Range::iterator rit = solidSets[m_or_s].begin(); rit != solidSets[m_or_s].end(); rit++) {
+    rval = mbImpl->get_entities_by_handle(*rit, tmp_range, true);
     RR("Failed to get entities in solid.");
   }
-  int dim = mb->dimension_from_handle(*tmp_range.rbegin());
+  int dim = mbImpl->dimension_from_handle(*tmp_range.rbegin());
   assert(dim > 0 && dim < 4);
   
-  solid_elems = tmp_range.subset_by_dimension(dim);
+  solidElems[m_or_s] = tmp_range.subset_by_dimension(dim);
 
-  setno_ptr = &FLUID_SETNO;
-  rval = mb->get_entities_by_type_and_tag(seth, MBENTITYSET, &tagh, &setno_ptr, 1, fluids);
-  if (fluids.empty()) rval = MB_FAILURE;
-  RR("Couldn't get any fluid sets.");
-  
-  for (Range::iterator rit = fluids.begin(); rit != fluids.end(); rit++) {
-    rval = mb->get_entities_by_dimension(*rit, dim, fluid_elems, true);
-    RR("Failed to get entities in fluid.");
+  for (std::set<int>::iterator sit = fluidSetNos.begin(); sit != fluidSetNos.end(); sit++) {
+    Range sets;
+    int set_no = *sit;
+    const void *setno_ptr = &set_no;
+    rval = mbImpl->get_entities_by_type_and_tag(seth, MBENTITYSET, &tagh, &setno_ptr, 1, sets);
+    if (sets.empty()) rval = MB_FAILURE;
+    RR("Couldn't get any fluid sets.");
+    fluidSets[m_or_s].merge(sets);
   }
-  if (mb->dimension_from_handle(*fluid_elems.begin()) != dim) {
-    rval = MB_FAILURE;
-    RR("Fluid and solid elements must be same dimension.");
+
+    // get fluid entities, and dimension
+  tmp_range.clear();
+  for (Range::iterator rit = fluidSets[m_or_s].begin(); rit != fluidSets[m_or_s].end(); rit++) {
+    rval = mbImpl->get_entities_by_handle(*rit, tmp_range, true);
+    RR("Failed to get entities in fluid.");
   }
   
+  fluidElems[m_or_s] = tmp_range.subset_by_dimension(dim);
+  
   return MB_SUCCESS;
 }


https://bitbucket.org/fathomteam/moab/commits/4d437a52f00a/
Changeset:   4d437a52f00a
Branch:      None
User:        tautges
Date:        2013-11-28 01:35:14
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 50643cb..a6e3d1f 100644
--- a/examples/makefile
+++ b/examples/makefile
@@ -14,47 +14,47 @@ 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
 	${MOAB_CXX} -o $@ $< ${MOAB_LIBS_LINK}
 
 PushParMeshIntoMoabF90: PushParMeshIntoMoabF90.o
 	${MOAB_CXX} -o $@ $< ${IMESH_LIBS} -lgfortran -L/usr/lib/openmpi/lib -lmpi_f90 -lmpi_f77 -lmpi -lopen-rte -lopen-pal -ldl 
 
-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/12c6587bc6e3/
Changeset:   12c6587bc6e3
Branch:      None
User:        tautges
Date:        2013-11-28 01:35:14
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 &params, double tol, const CartVect &init) const
+    bool SpectralHex::evaluate_reverse(CartVect const & xyz, CartVect &params, 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 &params, 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 7293e17..4978f79 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 e8e4c2b..c77bdbc 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 8d968c5..29dbb3d 100644
--- a/tools/mbcoupler/Coupler.cpp
+++ b/tools/mbcoupler/Coupler.cpp
@@ -654,7 +654,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;
@@ -671,7 +671,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;


https://bitbucket.org/fathomteam/moab/commits/9ec798a8f09c/
Changeset:   9ec798a8f09c
Branch:      None
User:        tautges
Date:        2013-11-28 01:35:15
Summary:     LloydSmoother: proper handling of communicator.
test/perf/point_location files: adding tolerance args for two-tolerance changes.

Affected #:  5 files

diff --git a/src/LloydSmoother.cpp b/src/LloydSmoother.cpp
index 09546b5..d60b615 100644
--- a/src/LloydSmoother.cpp
+++ b/src/LloydSmoother.cpp
@@ -7,6 +7,7 @@
 
 #ifdef USE_MPI
 #include "moab/ParallelComm.hpp"
+#include "MBParallelConventions.h"
 #endif
 
 #include <iostream>
@@ -161,7 +162,7 @@ ErrorCode LloydSmoother::perform_smooth()
 #ifdef USE_MPI
     // 2c. exchange tags on owned verts
     if (myPcomm->size() > 1) {
-      rval = pcomm->exchange_tags(centroid, shared_owned_verts); RR("Failed to exchange tags.");
+      rval = myPcomm->exchange_tags(centroid, shared_owned_verts); RR("Failed to exchange tags.");
     }
 #endif
 
@@ -171,7 +172,7 @@ ErrorCode LloydSmoother::perform_smooth()
         // global reduce for maximum delta, then report it
       if (myPcomm->size() > 1)
         MPI_Reduce(&resid, &global_max, 1, MPI_DOUBLE, MPI_MAX, 0, myPcomm->comm());
-      if (!pcomm->rank()) 
+      if (!myPcomm->rank()) 
 #endif
         std::cout << "Max residual = " << global_max << std::endl;
     }

diff --git a/test/perf/point_location/elem_eval_time.cpp b/test/perf/point_location/elem_eval_time.cpp
index 186bd81..ef86762 100644
--- a/test/perf/point_location/elem_eval_time.cpp
+++ b/test/perf/point_location/elem_eval_time.cpp
@@ -158,7 +158,7 @@ ErrorCode time_reverse_eval(Interface *mbi, int method, Range &elems,
     bool ins;
     for (rit = elems.begin(), i = 0; rit != elems.end(); rit++, i++) {
       eeval.set_ent_handle(*rit);
-      rval = eeval.reverse_eval(coords[i].array(), 1.0e-6, params[i].array(), &ins);
+      rval = eeval.reverse_eval(coords[i].array(), 1.0e-10, 1.0e-6, params[i].array(), &ins);
       assert(ins);
 #ifndef NDEBUG
       if (MB_SUCCESS != rval) return rval;

diff --git a/test/perf/point_location/point_location.cpp b/test/perf/point_location/point_location.cpp
index df4132e..4f5ddcc 100644
--- a/test/perf/point_location/point_location.cpp
+++ b/test/perf/point_location/point_location.cpp
@@ -290,7 +290,7 @@ void do_kdtree_test( Interface& mb, int tree_depth, int elem_per_leaf,
   int len;
   for (long i = 0; i < num_test; ++i) {
     const size_t idx = (size_t)i % points.size();
-    rval = tool.point_search( points[idx].array(), leaf, 0.0, NULL, &root ); CHK(rval);
+    rval = tool.point_search( points[idx].array(), leaf, 1.0e-10, 1.0e-6, NULL, &root ); CHK(rval);
     hexes.clear();
     rval = mb.get_entities_by_handle( leaf, hexes ); CHK(rval);
     for (j = hexes.begin(); j != hexes.end(); ++j) {

diff --git a/test/perf/point_location/sploc_searching_perf.cpp b/test/perf/point_location/sploc_searching_perf.cpp
index bf09a8e..158504d 100644
--- a/test/perf/point_location/sploc_searching_perf.cpp
+++ b/test/perf/point_location/sploc_searching_perf.cpp
@@ -160,7 +160,7 @@ ErrorCode test_locator(SpatialLocator &sl, int npoints, double rtol, double &cpu
   CpuTimer ct;
   
     // call spatial locator to locate points
-  rval = sl.locate_points(test_pts[0].array(), npoints, &ents[0], test_res[0].array(), rtol, 0.0, &is_in[0]);
+  rval = sl.locate_points(test_pts[0].array(), npoints, &ents[0], test_res[0].array(), &is_in[0], rtol, 0.0);
   if (MB_SUCCESS != rval) return rval;
 
   cpu_time = ct.time_elapsed();

diff --git a/test/perf/point_location/tree_searching_perf.cpp b/test/perf/point_location/tree_searching_perf.cpp
index 6d1badb..9e35d9b 100644
--- a/test/perf/point_location/tree_searching_perf.cpp
+++ b/test/perf/point_location/tree_searching_perf.cpp
@@ -32,9 +32,10 @@ int main(int argc, char **argv)
 #endif
 
   int npoints = 100, dim = 3;
-  int dints = 1, dleafs = 1, ddeps = 1;
+  int dints = 1, dleafs = 1, ddeps = 1, csints = 0;
   
-  ProgOptions po("tree_searching_perf options" );
+  ProgOptions po;
+  po.addOpt<int>( "candidateplaneset,c", "Candidate plane set (0=SUBDIVISION,1=SUBDIV_SNAP,2=VERTEX_MEDIAN,3=VERTEX_SAMPLE", &csints);
   po.addOpt<int>( "ints,i", "Number of doublings of intervals on each side of scd mesh", &dints);
   po.addOpt<int>( "leaf,l", "Number of doublings of maximum number of elements per leaf", &dleafs);
   po.addOpt<int>( "max_depth,m", "Number of 5-intervals on maximum depth of tree", &ddeps);
@@ -78,7 +79,7 @@ int main(int argc, char **argv)
       for (std::vector<int>::iterator leafs_it = leafs.begin(); leafs_it != leafs.end(); leafs_it++) {
   
           // iteration: tree type
-        for (int tree_tp = 0; tree_tp < 2; tree_tp++) {
+        for (int tree_tp = 1; tree_tp < 2; tree_tp++) {
             // create tree
           Tree *tree;
           if (0 == tree_tp)
@@ -88,6 +89,11 @@ int main(int argc, char **argv)
 
           std::ostringstream opts;
           opts << "MAX_DEPTH=" << *dep_it << ";MAX_PER_LEAF=" << *leafs_it;
+          if (csints) {
+            if (opts.str().length() > 0) 
+              opts << ";";
+            opts << "PLANE_SET=" << csints;
+          }
           FileOptions fo(opts.str().c_str());
           rval = tree->parse_options(fo);
           SpatialLocator sl(&mb, elems, tree);
@@ -145,7 +151,7 @@ ErrorCode test_locator(SpatialLocator &sl, int npoints, double &cpu_time, double
   CpuTimer ct;
   
     // call spatial locator to locate points
-  rval = sl.locate_points(test_pts[0].array(), npoints, &ents[0], test_res[0].array(), 0.0, 0.0, &is_in[0]);
+  rval = sl.locate_points(test_pts[0].array(), npoints, &ents[0], test_res[0].array(), &is_in[0]);
   if (MB_SUCCESS != rval) return rval;
 
   cpu_time = ct.time_elapsed();


https://bitbucket.org/fathomteam/moab/commits/aadd1febe26b/
Changeset:   aadd1febe26b
Branch:      None
User:        tautges
Date:        2013-11-28 01:45:59
Summary:     Adding a Lloyd smoother tool, and a test of it.

Affected #:  3 files

diff --git a/src/LloydSmoother.cpp b/src/LloydSmoother.cpp
index d60b615..f4892e9 100644
--- a/src/LloydSmoother.cpp
+++ b/src/LloydSmoother.cpp
@@ -172,7 +172,7 @@ ErrorCode LloydSmoother::perform_smooth()
         // global reduce for maximum delta, then report it
       if (myPcomm->size() > 1)
         MPI_Reduce(&resid, &global_max, 1, MPI_DOUBLE, MPI_MAX, 0, myPcomm->comm());
-      if (!myPcomm->rank()) 
+      if (!myPcomm || !myPcomm->rank()) 
 #endif
         std::cout << "Max residual = " << global_max << std::endl;
     }

diff --git a/src/Makefile.am b/src/Makefile.am
index 3a5e142..0d474b0 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -166,6 +166,7 @@ nobase_libMOAB_la_include_HEADERS = \
   moab/Forward.hpp \
   moab/GeomUtil.hpp \
   moab/Interface.hpp \
+  moab/LloydSmoother.hpp \
   moab/point_locater/tree/common_tree.hpp \
   moab/point_locater/tree/element_tree.hpp \
   moab/point_locater/tree/bvh_tree.hpp \

diff --git a/test/lloyd_smoother_test.cpp b/test/lloyd_smoother_test.cpp
index 662732b..9f4fd5a 100644
--- a/test/lloyd_smoother_test.cpp
+++ b/test/lloyd_smoother_test.cpp
@@ -1,6 +1,5 @@
 #include "moab/Core.hpp"
 #include "moab/LloydSmoother.hpp"
-#include "moab/CartVect.hpp"
 #include "TestUtil.hpp"
 
 #ifdef MESHDIR
@@ -72,7 +71,5 @@ int main(int argc, char**argv)
   rval = ll2.perform_smooth();
   CHECK_ERR(rval);
   std::cout << "Mesh smoothed in " << ll2.num_its() << " iterations." << std::endl;
-  
-  
 }
 


https://bitbucket.org/fathomteam/moab/commits/05ec91e72d11/
Changeset:   05ec91e72d11
Branch:      None
User:        tautges
Date:        2013-11-28 01:46:28
Summary:     Adding ability to smooth based on a tag instead of coords.

Affected #:  1 file

diff --git a/test/lloyd_smoother_test.cpp b/test/lloyd_smoother_test.cpp
index 9f4fd5a..8804d48 100644
--- a/test/lloyd_smoother_test.cpp
+++ b/test/lloyd_smoother_test.cpp
@@ -1,5 +1,6 @@
 #include "moab/Core.hpp"
 #include "moab/LloydSmoother.hpp"
+#include "moab/CartVect.hpp"
 #include "TestUtil.hpp"
 
 #ifdef MESHDIR


https://bitbucket.org/fathomteam/moab/commits/b65d9a64b639/
Changeset:   b65d9a64b639
Branch:      None
User:        tautges
Date:        2013-11-28 01:54:46
Summary:     Adding a Lloyd smoother tool, and a test of it.

Affected #:  2 files

diff --git a/src/Makefile.am b/src/Makefile.am
index 0d474b0..c7f48a6 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -154,6 +154,7 @@ nobase_libMOAB_la_include_HEADERS = \
   moab/Core.hpp \
   moab/CpuTimer.hpp \
   moab/DualTool.hpp \
+  moab/EigenDecomp.hpp \
   moab/Error.hpp \
   moab/GeomTopoTool.hpp \
   moab/HigherOrderFactory.hpp \

diff --git a/test/lloyd_smoother_test.cpp b/test/lloyd_smoother_test.cpp
index 8804d48..df1f2b3 100644
--- a/test/lloyd_smoother_test.cpp
+++ b/test/lloyd_smoother_test.cpp
@@ -1,6 +1,5 @@
 #include "moab/Core.hpp"
 #include "moab/LloydSmoother.hpp"
-#include "moab/CartVect.hpp"
 #include "TestUtil.hpp"
 
 #ifdef MESHDIR
@@ -32,21 +31,7 @@ int main(int argc, char**argv)
     CHECK_ERR(MB_FAILURE);
   }
 
-    // get the vertex positions and set on an intermediate tag, to test smoothing for
-    // a tag instead of for coords
-  Tag ctag;
-  rval = mb.tag_get_handle("vcentroid", 3, MB_TYPE_DOUBLE, ctag, MB_TAG_CREAT|MB_TAG_DENSE);
-  CHECK_ERR(rval);
-  Range verts;
-  rval = mb.get_entities_by_dimension(0, 0, verts);
-  CHECK_ERR(rval);
-  std::vector<double> coords(3*verts.size());
-  rval = mb.get_coords(verts, &coords[0]);
-  CHECK_ERR(rval);
-  rval = mb.tag_set_data(ctag, verts, &coords[0]);
-  CHECK_ERR(rval);
-
-  LloydSmoother ll(&mb, NULL, elems, ctag);
+  LloydSmoother ll(&mb, NULL, elems);
   ll.report_its(10);
   rval = ll.perform_smooth();
   CHECK_ERR(rval);


https://bitbucket.org/fathomteam/moab/commits/a44e3309825b/
Changeset:   a44e3309825b
Branch:      None
User:        tautges
Date:        2013-11-28 01:55:05
Summary:     Adding ability to smooth based on a tag instead of coords.

Affected #:  1 file

diff --git a/test/lloyd_smoother_test.cpp b/test/lloyd_smoother_test.cpp
index df1f2b3..8804d48 100644
--- a/test/lloyd_smoother_test.cpp
+++ b/test/lloyd_smoother_test.cpp
@@ -1,5 +1,6 @@
 #include "moab/Core.hpp"
 #include "moab/LloydSmoother.hpp"
+#include "moab/CartVect.hpp"
 #include "TestUtil.hpp"
 
 #ifdef MESHDIR
@@ -31,7 +32,21 @@ int main(int argc, char**argv)
     CHECK_ERR(MB_FAILURE);
   }
 
-  LloydSmoother ll(&mb, NULL, elems);
+    // get the vertex positions and set on an intermediate tag, to test smoothing for
+    // a tag instead of for coords
+  Tag ctag;
+  rval = mb.tag_get_handle("vcentroid", 3, MB_TYPE_DOUBLE, ctag, MB_TAG_CREAT|MB_TAG_DENSE);
+  CHECK_ERR(rval);
+  Range verts;
+  rval = mb.get_entities_by_dimension(0, 0, verts);
+  CHECK_ERR(rval);
+  std::vector<double> coords(3*verts.size());
+  rval = mb.get_coords(verts, &coords[0]);
+  CHECK_ERR(rval);
+  rval = mb.tag_set_data(ctag, verts, &coords[0]);
+  CHECK_ERR(rval);
+
+  LloydSmoother ll(&mb, NULL, elems, ctag);
   ll.report_its(10);
   rval = ll.perform_smooth();
   CHECK_ERR(rval);


https://bitbucket.org/fathomteam/moab/commits/c6066c16c556/
Changeset:   c6066c16c556
Branch:      None
User:        tautges
Date:        2013-11-28 01:55:06
Summary:     Handle case of parallel compile but 1 proc.

Affected #:  1 file

diff --git a/src/LloydSmoother.cpp b/src/LloydSmoother.cpp
index f4892e9..b1614d2 100644
--- a/src/LloydSmoother.cpp
+++ b/src/LloydSmoother.cpp
@@ -78,7 +78,7 @@ ErrorCode LloydSmoother::perform_smooth()
   Range owned_verts, shared_owned_verts;
 #ifdef USE_MPI
     // filter verts down to owned ones and get fixed tag for them
-  if (myPcomm->size() > 1) {
+  if (myPcomm && myPcomm->size() > 1) {
     rval = myPcomm->filter_pstatus(verts, PSTATUS_NOT_OWNED, PSTATUS_NOT, -1, &owned_verts);
     RR("Failed to filter on pstatus.");
       // get shared owned verts, for exchanging tags
@@ -161,7 +161,7 @@ ErrorCode LloydSmoother::perform_smooth()
 
 #ifdef USE_MPI
     // 2c. exchange tags on owned verts
-    if (myPcomm->size() > 1) {
+    if (myPcomm && myPcomm->size() > 1) {
       rval = myPcomm->exchange_tags(centroid, shared_owned_verts); RR("Failed to exchange tags.");
     }
 #endif
@@ -170,7 +170,7 @@ ErrorCode LloydSmoother::perform_smooth()
       double global_max = resid;
 #ifdef USE_MPI
         // global reduce for maximum delta, then report it
-      if (myPcomm->size() > 1)
+      if (myPcomm && myPcomm->size() > 1)
         MPI_Reduce(&resid, &global_max, 1, MPI_DOUBLE, MPI_MAX, 0, myPcomm->comm());
       if (!myPcomm || !myPcomm->rank()) 
 #endif


https://bitbucket.org/fathomteam/moab/commits/72f58109faeb/
Changeset:   72f58109faeb
Branch:      None
User:        tautges
Date:        2013-11-28 01:55:06
Summary:     Add support for parallel, and a few better comments for DeformMeshRemap class.

Affected #:  1 file

diff --git a/examples/DeformMeshRemap.cpp b/examples/DeformMeshRemap.cpp
index 7614934..5442d1d 100644
--- a/examples/DeformMeshRemap.cpp
+++ b/examples/DeformMeshRemap.cpp
@@ -19,8 +19,13 @@
 #include "MBTagConventions.hpp"
 #include "DataCoupler.hpp"
 
+#ifdef USE_MPI
+#  include "moab/ParallelComm.hpp"
+#endif
+
 #include <iostream>
 #include <set>
+#include <sstream>
 #include <assert.h>
 
 using namespace moab;
@@ -52,7 +57,10 @@ public:
   enum {MASTER=0, SLAVE, SOLID, FLUID};
   
     //! constructor
-  DeformMeshRemap(Interface *impl) : mbImpl(impl), masterSet(0), slaveSet(0), xNew(0), xNewName("xnew") {}
+    //! if master is NULL, the MOAB part is run in serial; 
+    //! if slave is NULL but the master isn't, the slave is copied from the master
+    //! Create communicators using moab::ParallelComm::get_pcomm
+  DeformMeshRemap(Interface *impl, ParallelComm *master = NULL, ParallelComm *slave = NULL);
   
     //! destructor
   ~DeformMeshRemap();
@@ -100,6 +108,11 @@ private:
 
     //! moab interface
   Interface *mbImpl;
+
+#ifdef USE_MPI
+    //! ParallelComm for master, slave meshes
+  ParallelComm *pcMaster, *pcSlave;
+#endif
   
     //! material set numbers for fluid materials
   std::set<int> fluidSetNos;
@@ -238,8 +251,18 @@ ErrorCode DeformMeshRemap::execute()
   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);
+    std::string str;
+#ifdef USE_MPI
+    if (pcMaster && pcMaster->size() > 1) 
+      str = "PARALLEL=WRITE_PART";
+#endif
+    rval = mbImpl->write_file("smoothed_master.vtk", NULL, str.c_str(), &masterSet, 1);
+#ifdef USE_MPI
+    str.clear();
+    if (pcSlave && pcSlave->size() > 1) 
+      str = "PARALLEL=WRITE_PART";
+#endif
+    rval = mbImpl->write_file("slave_interp.vtk", NULL, str.c_str(), &slaveSet, 1);
   }
 
   return MB_SUCCESS;
@@ -270,6 +293,13 @@ void DeformMeshRemap::set_file_name(int m_or_s, const std::string &name)
   }
 }
 
+DeformMeshRemap::DeformMeshRemap(Interface *impl, ParallelComm *master, ParallelComm *slave)  
+        : mbImpl(impl), pcMaster(master), pcSlave(slave), masterSet(0), slaveSet(0), xNew(0), xNewName("xnew") 
+{
+  if (!pcSlave && pcMaster)
+    pcSlave = pcMaster;
+}
+  
 DeformMeshRemap::~DeformMeshRemap() 
 {
     // delete the tag
@@ -293,13 +323,23 @@ int main(int argc, char **argv) {
 
   mb = new Core();
 
-  DeformMeshRemap dfr(mb);
-  dfr.set_file_name(DeformMeshRemap::MASTER, masterf);
-  dfr.set_file_name(DeformMeshRemap::SLAVE, slavef);
-  rval = dfr.add_set_no(DeformMeshRemap::SOLID, SOLID_SETNO); RR("Failed to add solid set no.");
-  rval = dfr.add_set_no(DeformMeshRemap::FLUID, FLUID_SETNO); RR("Failed to add fluid set no.");
+  DeformMeshRemap *dfr;
+#ifdef USE_MPI
+  ParallelComm *pc = new ParallelComm(mb, MPI_COMM_WORLD);
+  dfr = new DeformMeshRemap(mb, pc);
+#else  
+  dfr = new DeformMeshRemap(mb);
+#endif
+  dfr->set_file_name(DeformMeshRemap::MASTER, masterf);
+  dfr->set_file_name(DeformMeshRemap::SLAVE, slavef);
+  rval = dfr->add_set_no(DeformMeshRemap::SOLID, SOLID_SETNO); RR("Failed to add solid set no.");
+  rval = dfr->add_set_no(DeformMeshRemap::FLUID, FLUID_SETNO); RR("Failed to add fluid set no.");
+  
+  rval = dfr->execute();
+  
+  delete dfr;
+  delete mb;
   
-  rval = dfr.execute();
   return rval;
 }
 
@@ -389,7 +429,16 @@ ErrorCode DeformMeshRemap::read_file(int m_or_s, string &fname, EntityHandle &se
     // create meshset
   ErrorCode rval = mbImpl->create_meshset(0, seth);
   RR("Couldn't create master/slave set.");
-  rval = mbImpl->load_file(fname.c_str(), &seth);
+  ostringstream ostr;
+#ifdef USE_MPI
+  ParallelComm *pc = (m_or_s == MASTER ? pcMaster : pcSlave);
+  if (pc && pc->size() > 1) {
+    if (debug) ostr << "DEBUG_IO=1;CPUTIME;";
+    ostr << "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARALLEL_RESOLVE_SHARED_ENTS;"
+         << "PARALLEL_GHOSTS=2.0.1;PARALLEL_COMM=" << pc->get_id();
+  }
+#endif  
+  rval = mbImpl->load_file(fname.c_str(), &seth, ostr.str().c_str());
   RR("Couldn't load master/slave mesh.");
 
     // get material sets for solid/fluid


https://bitbucket.org/fathomteam/moab/commits/cf13c67fd862/
Changeset:   cf13c67fd862
Branch:      None
User:        tautges
Date:        2013-11-28 01:55:06
Summary:     Small fixes revealed by make check.

LloydSmoother.cpp: test for NULL myPcomm before referencing it.
Makefile.am: remove EigenDecomp.hpp from Makefile.am.
others: add extra tolerance value to arg lists.

Affected #:  4 files

diff --git a/src/Makefile.am b/src/Makefile.am
index c7f48a6..0d474b0 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -154,7 +154,6 @@ nobase_libMOAB_la_include_HEADERS = \
   moab/Core.hpp \
   moab/CpuTimer.hpp \
   moab/DualTool.hpp \
-  moab/EigenDecomp.hpp \
   moab/Error.hpp \
   moab/GeomTopoTool.hpp \
   moab/HigherOrderFactory.hpp \

diff --git a/test/elem_eval_test.cpp b/test/elem_eval_test.cpp
index 1e669c3..beb8704 100644
--- a/test/elem_eval_test.cpp
+++ b/test/elem_eval_test.cpp
@@ -63,7 +63,7 @@ void test_eval(ElemEvaluator &ee, bool test_integrate)
         if (!ee.inside(params.array(), EPS1)) continue;
         
         rval = ee.eval(params.array(), posn.array()); CHECK_ERR(rval);
-        rval = ee.reverse_eval(posn.array(), EPS1, params2.array(), &is_inside); CHECK_ERR(rval);
+        rval = ee.reverse_eval(posn.array(), EPS1, EPS1, params2.array(), &is_inside); CHECK_ERR(rval);
         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);

diff --git a/test/kd_tree_test.cpp b/test/kd_tree_test.cpp
index 8f479d0..fd5389f 100644
--- a/test/kd_tree_test.cpp
+++ b/test/kd_tree_test.cpp
@@ -377,9 +377,9 @@ void test_point_search()
   }
   
     // compare leaf search to iterator search
-  rval = tool.point_search(right.array(), leaf, 0.0, NULL, const_cast<EntityHandle*>(&root));
+  rval = tool.point_search(right.array(), leaf, 0.0, 0.0, NULL, const_cast<EntityHandle*>(&root));
   CHECK_ERR(rval);
-  rval = tool.point_search(right.array(), iter, 0.0, NULL, const_cast<EntityHandle*>(&root));
+  rval = tool.point_search(right.array(), iter, 0.0, 0.0, NULL, const_cast<EntityHandle*>(&root));
   CHECK_ERR(rval);
   assert( iter.handle() == leaf );
   

diff --git a/test/kd_tree_time.cpp b/test/kd_tree_time.cpp
index af0e33d..9c34430 100644
--- a/test/kd_tree_time.cpp
+++ b/test/kd_tree_time.cpp
@@ -135,7 +135,7 @@ int main(int argc, char* argv[])
     if (query_triangles)
       rval = tool.closest_triangle(root, coords, pt, leaf);
     else
-      rval = tool.point_search(coords, leaf, 0.0, NULL, &root);
+      rval = tool.point_search(coords, leaf, 0.0, 0.0, NULL, &root);
     if (MB_SUCCESS != rval) {
       fprintf(stderr, "Failure (ErrorCode == %d) for point %d (%f, %f, %f)\n",
         (int)rval, (int)(count % length), coords[0], coords[1], coords[2]);


https://bitbucket.org/fathomteam/moab/commits/0eb03d434a42/
Changeset:   0eb03d434a42
Branch:      deformed-mesh-remap
User:        tautges
Date:        2013-11-28 02:40:23
Summary:     Merge branch 'deformed-mesh-remap' of bitbucket.org:fathomteam/moab into deformed-mesh-remap

Conflicts:
	examples/DeformMeshRemap.cpp
	examples/makefile
	src/LloydSmoother.cpp
	test/Makefile.am
	test/lloyd_smoother_test.cpp

Affected #:  0 files

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