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

commits-noreply at bitbucket.org commits-noreply at bitbucket.org
Wed Mar 12 11:06:16 CDT 2014


24 new commits in MOAB:

https://bitbucket.org/fathomteam/moab/commits/e0c608b42e90/
Changeset:   e0c608b42e90
Branch:      None
User:        tautges
Date:        2014-03-12 16:51:41
Summary:     Changes to allow specification of existing displacement tag(s), which can be either single name corresponding to xyz[3]
or 3 names x, y, z.  Handle either case in deform_master.

Affected #:  1 file

diff --git a/examples/DeformMeshRemap.cpp b/examples/DeformMeshRemap.cpp
index 5442d1d..23ad554 100644
--- a/examples/DeformMeshRemap.cpp
+++ b/examples/DeformMeshRemap.cpp
@@ -91,6 +91,10 @@ public:
   
     //! get/set the file name
   void set_file_name(int m_or_s, const std::string &name);
+
+    //! get/set the x displacement tag names
+  std::string xdisp_name(int idx = 0);
+  void xdisp_name(const std::string &nm, int idx = 0);
   
 private:
     //! apply a known deformation to the solid elements, putting the results in the xNew tag; also
@@ -132,9 +136,15 @@ private:
     //! filenames for master/slave meshes
   std::string masterFileName, slaveFileName;
 
+    //! tag from file, might be 3
+  Tag xDisp[3];
+
     //! tag used for new positions
   Tag xNew;
   
+    //! tag name used to read disps from file
+  std::string xDispNames[3];
+  
     //! tag name used for new positions
   std::string xNewName;
 };
@@ -199,6 +209,16 @@ inline ErrorCode DeformMeshRemap::get_set_nos(int f_or_s, std::set<int> &set_nos
   return MB_SUCCESS;
 }
 
+inline std::string DeformMeshRemap::xdisp_name(int idx) 
+{
+  return xDispNames[idx];
+}
+
+void DeformMeshRemap::xdisp_name(const std::string &nm, int idx) 
+{
+  xDispNames[idx] = nm;
+}
+
 ErrorCode DeformMeshRemap::execute() 
 {
     // read master/slave files and get fluid/solid material sets
@@ -296,14 +316,14 @@ 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") 
 {
+  xDisp[0] = xDisp[1] = xDisp[2] = 0;
+
   if (!pcSlave && pcMaster)
     pcSlave = pcMaster;
 }
   
 DeformMeshRemap::~DeformMeshRemap() 
 {
-    // delete the tag
-  mbImpl->tag_delete(xNew);
 }
 
 int main(int argc, char **argv) {
@@ -312,14 +332,14 @@ int main(int argc, char **argv) {
 
   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.addOpt<std::string> ("worker,w", "Specify the slave/worker meshfile name" );
+  po.addOpt<std::string> ("d1,", "Tag name for displacement x or xyz" );
+  po.addOpt<std::string> ("d2,", "Tag name for displacement y" );
+  po.addOpt<std::string> ("d3,", "Tag name for displacement z" );
+  po.addOpt<int> ("fluid,f", "Specify fluid material set number(s).");
+  po.addOpt<int> ("solid,s", "Specify fluid material set number(s).");
+  
   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();
 
@@ -330,10 +350,43 @@ int main(int argc, char **argv) {
 #else  
   dfr = new DeformMeshRemap(mb);
 #endif
+
+
+  std::string masterf, slavef;
+  if(!po.getOpt("master", &masterf))
+    masterf = string(MESH_DIR) + string("/rodquad.g");
   dfr->set_file_name(DeformMeshRemap::MASTER, masterf);
+
+  if(!po.getOpt("worker", &slavef))
+    slavef = string(MESH_DIR) + string("/rodtri.g");
   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.");
+
+  std::vector<int> set_nos;
+  po.getOptAllArgs("fluid", set_nos);
+  if (set_nos.empty()) {
+    rval = dfr->add_set_no(DeformMeshRemap::FLUID, FLUID_SETNO); RR("Failed to add solid set no.");
+  }
+  else {
+    for (std::vector<int>::iterator vit = set_nos.begin(); vit != set_nos.end(); vit++) 
+      dfr->add_set_no(DeformMeshRemap::FLUID, *vit);
+  }
+  set_nos.clear();
+  
+  po.getOptAllArgs("solid", set_nos);
+  if (set_nos.empty()) {
+    rval = dfr->add_set_no(DeformMeshRemap::SOLID, SOLID_SETNO); RR("Failed to add solid set no.");
+  }
+  else {
+    for (std::vector<int>::iterator vit = set_nos.begin(); vit != set_nos.end(); vit++) 
+      dfr->add_set_no(DeformMeshRemap::SOLID, *vit);
+  }
+
+  std::string tnames[3];
+  po.getOpt("d1", &tnames[0]);
+  po.getOpt("d2", &tnames[1]);
+  po.getOpt("d3", &tnames[2]);
+  for (int i = 0; i < 3; i++) 
+    if (!tnames[i].empty()) dfr->xdisp_name(tnames[i], i);
   
   rval = dfr->execute();
   
@@ -387,40 +440,71 @@ void deform_func(const BoundBox &bbox, double *xold, double *xnew)
 ErrorCode DeformMeshRemap::deform_master(Range &fluid_elems, Range &solid_elems, const char *tag_name) 
 {
     // deform elements with an analytic function
+  ErrorCode rval;
 
-    // create the tag
-  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
-  Range verts;
-  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 = mbImpl->get_coords(verts, &coords[0]);
-  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();
+  Range verts;
   rval = mbImpl->get_adjacencies(solid_elems, 0, false, verts, Interface::UNION);
   RR("Failed to get vertices.");
-  coords.resize(3*verts.size(), 0.0);
+  std::vector<double> coords(3*verts.size());
   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(bbox, &coords[3*i], &coords[3*i]);
+
+    // get or create the tag
+
+  if (!xDispNames[0].empty() && !xDispNames[1].empty() && !xDispNames[2].empty()) {
+      // 3 tags, specifying xyz individual data, integrate into one tag
+    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.");
+    std::vector<double> disps(num_verts);
+    for (int i = 0; i < 3; i++) {
+      rval = mbImpl->tag_get_handle(xDispNames[0].c_str(), 1, MB_TYPE_DOUBLE, xDisp[i]);
+      RR("Failed to get xDisp tag.");
+      rval = mbImpl->tag_get_data(xDisp[i], verts, &disps[0]);
+      RR("Failed to get xDisp tag values.");
+      for (unsigned int j = 0; j < num_verts; j++)
+        coords[3*j+i] = disps[j];
+    }
+  }
+  else if (!xDispNames[0].empty()) {
+    rval = mbImpl->tag_get_handle(xDispNames[0].c_str(), 3, MB_TYPE_DOUBLE, xDisp[0]);
+    RR("Failed to get first xDisp tag.");
+    xNew = xDisp[0];
+    std::vector<double> disps(3*num_verts);
+    rval = mbImpl->tag_get_data(xDisp[0], verts, &disps[0]);
+    for (unsigned int j = 0; j < 3*num_verts; j++)
+      coords[j] += disps[j];
+  }
+  else {
+      // get the bounding box of the solid mesh
+    BoundBox bbox;
+    bbox.update(*mbImpl, solid_elems);
+  
+    for (unsigned int j = 0; j < num_verts; j++)
+      deform_func(bbox, &coords[3*j], &coords[3*j]);
+  }
+
+  if (!xNew) {
+    rval = mbImpl->tag_get_handle((tag_name ? tag_name : ""), 3, MB_TYPE_DOUBLE, 
+                                  xDisp[0], MB_TAG_CREAT|MB_TAG_DENSE);
+    RR("Failed to get xNew tag.");
+  }
     
     // set the new tag to those coords
   rval = mbImpl->tag_set_data(xNew, verts, &coords[0]);
   RR("Failed to set tag data.");
   
+    // get all the vertices and coords in the fluid, set xnew to them
+  verts.clear();
+  rval = mbImpl->get_adjacencies(fluid_elems, 0, false, verts, Interface::UNION);
+  RR("Failed to get vertices.");
+  coords.resize(3*verts.size());
+  rval = mbImpl->get_coords(verts, &coords[0]);
+  RR("Failed to get vertex coords.");
+  rval = mbImpl->tag_set_data(xNew, verts, &coords[0]);
+  RR("Failed to set xnew tag on fluid verts.");
+    
   return MB_SUCCESS;
 }
 


https://bitbucket.org/fathomteam/moab/commits/37f2b44c6262/
Changeset:   37f2b44c6262
Branch:      None
User:        tautges
Date:        2014-03-12 16:53:02
Summary:     - add ability to specify master and slave mesh solid/fluid set numbers separately, so set numbers don't need to agree
- add ability to specify only fluid or solid sets, with other sets derived from that and all material sets in a file
- add DeformMeshRemap to list of examples in makefile

Affected #:  1 file

diff --git a/examples/DeformMeshRemap.cpp b/examples/DeformMeshRemap.cpp
index 23ad554..b92d994 100644
--- a/examples/DeformMeshRemap.cpp
+++ b/examples/DeformMeshRemap.cpp
@@ -69,13 +69,13 @@ public:
   ErrorCode execute();
   
     //! add a set number
-  ErrorCode add_set_no(int fluid_or_solid, int set_no);
+  ErrorCode add_set_no(int m_or_s, int fluid_or_solid, int set_no);
   
     //! remove a set number
-  ErrorCode remove_set_no(int fluid_or_solid, int set_no);
+  ErrorCode remove_set_no(int m_or_s, 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;
+  ErrorCode get_set_nos(int m_or_s, int fluid_or_solid, std::set<int> &set_nos) const;
 
     //! get the xNew tag handle
   inline Tag x_new() const {return xNew;}
@@ -110,6 +110,9 @@ private:
     //! 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);
 
+    //! find fluid/solid sets from complement of solid/fluid sets 
+  ErrorCode find_other_sets(int m_or_s, EntityHandle file_set);
+  
     //! moab interface
   Interface *mbImpl;
 
@@ -118,11 +121,11 @@ private:
   ParallelComm *pcMaster, *pcSlave;
 #endif
   
-    //! material set numbers for fluid materials
-  std::set<int> fluidSetNos;
+    //! material set numbers for fluid materials, for master/slave
+  std::set<int> fluidSetNos[2];
 
-    //! material set numbers for solid materials
-  std::set<int> solidSetNos;
+    //! material set numbers for solid materials, for master/slave
+  std::set<int> solidSetNos[2];
 
     //! sets defining master/slave meshes
   EntityHandle masterSet, slaveSet;
@@ -150,14 +153,17 @@ private:
 };
 
   //! add a set number
-inline ErrorCode DeformMeshRemap::add_set_no(int f_or_s, int set_no) 
+inline ErrorCode DeformMeshRemap::add_set_no(int m_or_s, int f_or_s, int set_no) 
 {
   std::set<int> *this_set;
+  assert ((m_or_s == MASTER || m_or_s == SLAVE) && "m_or_s should be MASTER or SLAVE.");
+  if (m_or_s != MASTER && m_or_s != SLAVE) return MB_INDEX_OUT_OF_RANGE;
+  
   switch (f_or_s) {
     case FLUID:
-        this_set = &fluidSetNos; break;
+        this_set = &fluidSetNos[m_or_s]; break;
     case SOLID:
-        this_set = &solidSetNos; break;
+        this_set = &solidSetNos[m_or_s]; break;
     default:
         assert(false && "f_or_s should be FLUID or SOLID.");
         return MB_FAILURE;
@@ -169,14 +175,16 @@ inline ErrorCode DeformMeshRemap::add_set_no(int f_or_s, int set_no)
 }
   
   //! remove a set number
-inline ErrorCode DeformMeshRemap::remove_set_no(int f_or_s, int set_no) 
+inline ErrorCode DeformMeshRemap::remove_set_no(int m_or_s, int f_or_s, int set_no) 
 {
   std::set<int> *this_set;
+  assert ((m_or_s == MASTER || m_or_s == SLAVE) && "m_or_s should be MASTER or SLAVE.");
+  if (m_or_s != MASTER && m_or_s != SLAVE) return MB_INDEX_OUT_OF_RANGE;
   switch (f_or_s) {
     case FLUID:
-        this_set = &fluidSetNos; break;
+        this_set = &fluidSetNos[m_or_s]; break;
     case SOLID:
-        this_set = &solidSetNos; break;
+        this_set = &solidSetNos[m_or_s]; break;
     default:
         assert(false && "f_or_s should be FLUID or SOLID.");
         return MB_FAILURE;
@@ -191,14 +199,16 @@ inline ErrorCode DeformMeshRemap::remove_set_no(int f_or_s, int set_no)
 }
   
   //! get the set numbers
-inline ErrorCode DeformMeshRemap::get_set_nos(int f_or_s, std::set<int> &set_nos) const
+inline ErrorCode DeformMeshRemap::get_set_nos(int m_or_s, int f_or_s, std::set<int> &set_nos) const
 {
   const std::set<int> *this_set;
+  assert ((m_or_s == MASTER || m_or_s == SLAVE) && "m_or_s should be MASTER or SLAVE.");
+  if (m_or_s != MASTER && m_or_s != SLAVE) return MB_INDEX_OUT_OF_RANGE;
   switch (f_or_s) {
     case FLUID:
-        this_set = &fluidSetNos; break;
+        this_set = &fluidSetNos[m_or_s]; break;
     case SOLID:
-        this_set = &solidSetNos; break;
+        this_set = &solidSetNos[m_or_s]; break;
     default:
         assert(false && "f_or_s should be FLUID or SOLID.");
         return MB_FAILURE;
@@ -224,10 +234,18 @@ 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;
+
+  if (*solidSetNos[MASTER].begin() == -1 || *fluidSetNos[MASTER].begin() == -1) {
+    rval = find_other_sets(MASTER, masterSet); RR("Failed to find other sets in master mesh.");
+  }
   
   rval = read_file(SLAVE, slaveFileName, slaveSet);
   if (MB_SUCCESS != rval) return rval;
 
+  if (*solidSetNos[SLAVE].begin() == -1 || *fluidSetNos[SLAVE].begin() == -1) {
+    rval = find_other_sets(SLAVE, slaveSet); RR("Failed to find other sets in slave mesh.");
+  }
+  
   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
@@ -288,6 +306,50 @@ ErrorCode DeformMeshRemap::execute()
   return MB_SUCCESS;
 }
 
+ErrorCode DeformMeshRemap::find_other_sets(int m_or_s, EntityHandle file_set) 
+{
+    // solid or fluid sets are missing; find the other
+  Range *filled_sets = NULL, *unfilled_sets = NULL, *unfilled_elems = NULL;
+  
+  if (fluidSets[m_or_s].empty() && !solidSets[m_or_s].empty()) {
+    unfilled_sets = &fluidSets[m_or_s];
+    filled_sets = &solidSets[m_or_s];
+    unfilled_elems = &fluidElems[m_or_s];
+  }
+  else if (!fluidSets[m_or_s].empty() && solidSets[m_or_s].empty()) {
+    filled_sets = &fluidSets[m_or_s];
+    unfilled_sets = &solidSets[m_or_s];
+    unfilled_elems = &solidElems[m_or_s];
+  }
+  
+    // ok, we know the filled sets, now fill the unfilled sets, and the elems from those
+  Tag tagh;
+  ErrorCode rval = mbImpl->tag_get_handle(MATERIAL_SET_TAG_NAME, tagh); RR("Couldn't get material set tag name.");
+  Range matsets;
+  rval = mbImpl->get_entities_by_type_and_tag(file_set, MBENTITYSET, &tagh, NULL, 1, matsets);
+  if (matsets.empty()) rval = MB_FAILURE;
+  RR("Couldn't get any material sets.");
+  *unfilled_sets = subtract(matsets, *filled_sets);
+  if (unfilled_sets->empty()) {
+    rval = MB_FAILURE;
+    RR("Failed to find any unfilled material sets.");
+  }
+  Range tmp_range;
+  for (Range::iterator rit = unfilled_sets->begin(); rit != unfilled_sets->end(); rit++) {
+    rval = mbImpl->get_entities_by_handle(*rit, tmp_range, true);
+    RR("Failed to get entities in unfilled set.");
+  }
+  int dim = mbImpl->dimension_from_handle(*tmp_range.rbegin());
+  assert(dim > 0 && dim < 4);  
+  *unfilled_elems = tmp_range.subset_by_dimension(dim);
+  if (unfilled_elems->empty()) {
+    rval = MB_FAILURE;
+    RR("Failed to find any unfilled set entities.");
+  }
+  
+  return MB_SUCCESS;
+}
+
 std::string DeformMeshRemap::get_file_name(int m_or_s) const
 {
   switch (m_or_s) {
@@ -336,8 +398,10 @@ int main(int argc, char **argv) {
   po.addOpt<std::string> ("d1,", "Tag name for displacement x or xyz" );
   po.addOpt<std::string> ("d2,", "Tag name for displacement y" );
   po.addOpt<std::string> ("d3,", "Tag name for displacement z" );
-  po.addOpt<int> ("fluid,f", "Specify fluid material set number(s).");
-  po.addOpt<int> ("solid,s", "Specify fluid material set number(s).");
+  po.addOpt<int> ("fm,", "Specify master fluid material set number(s). If -1 specified, fluid sets derived from complement of solid sets.");
+  po.addOpt<int> ("fs,", "Specify master solid material set number(s). If -1 specified, fluid sets derived from complement of solid sets.");
+  po.addOpt<int> ("sm,", "Specify slave fluid material set number(s). If -1 specified, fluid sets derived from complement of solid sets.");
+  po.addOpt<int> ("ss,", "Specify slave solid material set number(s). If -1 specified, fluid sets derived from complement of solid sets.");
   
   po.parseCommandLine(argc, argv);
 
@@ -362,23 +426,44 @@ int main(int argc, char **argv) {
   dfr->set_file_name(DeformMeshRemap::SLAVE, slavef);
 
   std::vector<int> set_nos;
-  po.getOptAllArgs("fluid", set_nos);
+  po.getOptAllArgs("fm", set_nos);
   if (set_nos.empty()) {
-    rval = dfr->add_set_no(DeformMeshRemap::FLUID, FLUID_SETNO); RR("Failed to add solid set no.");
+    rval = dfr->add_set_no(DeformMeshRemap::MASTER, DeformMeshRemap::FLUID, FLUID_SETNO); 
+    RR("Failed to add solid set no.");
   }
   else {
     for (std::vector<int>::iterator vit = set_nos.begin(); vit != set_nos.end(); vit++) 
-      dfr->add_set_no(DeformMeshRemap::FLUID, *vit);
+      dfr->add_set_no(DeformMeshRemap::MASTER, DeformMeshRemap::FLUID, *vit);
   }
   set_nos.clear();
   
-  po.getOptAllArgs("solid", set_nos);
+  po.getOptAllArgs("fs", set_nos);
   if (set_nos.empty()) {
-    rval = dfr->add_set_no(DeformMeshRemap::SOLID, SOLID_SETNO); RR("Failed to add solid set no.");
+    rval = dfr->add_set_no(DeformMeshRemap::SLAVE, DeformMeshRemap::FLUID, FLUID_SETNO); 
+    RR("Failed to add solid set no.");
   }
   else {
     for (std::vector<int>::iterator vit = set_nos.begin(); vit != set_nos.end(); vit++) 
-      dfr->add_set_no(DeformMeshRemap::SOLID, *vit);
+      dfr->add_set_no(DeformMeshRemap::SLAVE, DeformMeshRemap::FLUID, *vit);
+  }
+  set_nos.clear();
+  
+  po.getOptAllArgs("sm", set_nos);
+  if (set_nos.empty()) {
+    rval = dfr->add_set_no(DeformMeshRemap::MASTER, DeformMeshRemap::SOLID, SOLID_SETNO); RR("Failed to add solid set no.");
+  }
+  else {
+    for (std::vector<int>::iterator vit = set_nos.begin(); vit != set_nos.end(); vit++) 
+      dfr->add_set_no(DeformMeshRemap::MASTER, DeformMeshRemap::SOLID, *vit);
+  }
+
+  po.getOptAllArgs("ss", set_nos);
+  if (set_nos.empty()) {
+    rval = dfr->add_set_no(DeformMeshRemap::SLAVE, DeformMeshRemap::SOLID, SOLID_SETNO); RR("Failed to add solid set no.");
+  }
+  else {
+    for (std::vector<int>::iterator vit = set_nos.begin(); vit != set_nos.end(); vit++) 
+      dfr->add_set_no(DeformMeshRemap::SLAVE, DeformMeshRemap::SOLID, *vit);
   }
 
   std::string tnames[3];
@@ -489,6 +574,7 @@ ErrorCode DeformMeshRemap::deform_master(Range &fluid_elems, Range &solid_elems,
     rval = mbImpl->tag_get_handle((tag_name ? tag_name : ""), 3, MB_TYPE_DOUBLE, 
                                   xDisp[0], MB_TAG_CREAT|MB_TAG_DENSE);
     RR("Failed to get xNew tag.");
+    xNew = xDisp[0];
   }
     
     // set the new tag to those coords
@@ -525,10 +611,12 @@ ErrorCode DeformMeshRemap::read_file(int m_or_s, string &fname, EntityHandle &se
   rval = mbImpl->load_file(fname.c_str(), &seth, ostr.str().c_str());
   RR("Couldn't load master/slave mesh.");
 
+  if (*solidSetNos[m_or_s].begin() == -1 || *fluidSetNos[m_or_s].begin() == -1) return MB_SUCCESS;
+  
     // get material sets for solid/fluid
   Tag tagh;
   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++) {
+  for (std::set<int>::iterator sit = solidSetNos[m_or_s].begin(); sit != solidSetNos[m_or_s].end(); sit++) {
     Range sets;
     int set_no = *sit;
     const void *setno_ptr = &set_no;
@@ -549,7 +637,7 @@ ErrorCode DeformMeshRemap::read_file(int m_or_s, string &fname, EntityHandle &se
   
   solidElems[m_or_s] = tmp_range.subset_by_dimension(dim);
 
-  for (std::set<int>::iterator sit = fluidSetNos.begin(); sit != fluidSetNos.end(); sit++) {
+  for (std::set<int>::iterator sit = fluidSetNos[m_or_s].begin(); sit != fluidSetNos[m_or_s].end(); sit++) {
     Range sets;
     int set_no = *sit;
     const void *setno_ptr = &set_no;


https://bitbucket.org/fathomteam/moab/commits/a93cd6825b14/
Changeset:   a93cd6825b14
Branch:      None
User:        tautges
Date:        2014-03-12 16:53:02
Summary:     Adding the ability to compute on master only.

Affected #:  1 file

diff --git a/examples/DeformMeshRemap.cpp b/examples/DeformMeshRemap.cpp
index b92d994..ac3ca39 100644
--- a/examples/DeformMeshRemap.cpp
+++ b/examples/DeformMeshRemap.cpp
@@ -238,36 +238,42 @@ ErrorCode DeformMeshRemap::execute()
   if (*solidSetNos[MASTER].begin() == -1 || *fluidSetNos[MASTER].begin() == -1) {
     rval = find_other_sets(MASTER, masterSet); RR("Failed to find other sets in master mesh.");
   }
-  
-  rval = read_file(SLAVE, slaveFileName, slaveSet);
-  if (MB_SUCCESS != rval) return rval;
 
-  if (*solidSetNos[SLAVE].begin() == -1 || *fluidSetNos[SLAVE].begin() == -1) {
-    rval = find_other_sets(SLAVE, slaveSet); RR("Failed to find other sets in slave mesh.");
+  bool have_slave = !(slaveFileName == "none");
+  if (have_slave) {
+    rval = read_file(SLAVE, slaveFileName, slaveSet);
+    if (MB_SUCCESS != rval) return rval;
+
+    if (*solidSetNos[SLAVE].begin() == -1 || *fluidSetNos[SLAVE].begin() == -1) {
+      rval = find_other_sets(SLAVE, slaveSet); RR("Failed to find other sets in slave mesh.");
+    }
   }
   
   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;
-  }
 
+  Range tgt_verts;
+  if (have_slave) {
+      // locate slave vertices in master, orig coords; do this with a data coupler, so you can
+      // later interpolate
+    Range 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); 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("");
 
@@ -279,14 +285,17 @@ ErrorCode DeformMeshRemap::execute()
     RR("Failed in lloyd smoothing.");
     cout << "Lloyd smoothing required " << ll.num_its() << " iterations." << endl;
   }
-  
-    // map new locations to slave
-    // 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
+    // transfer xNew to coords, for master
   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 (have_slave) {
+      // map new locations to slave
+      // interpolate xNew to slave points
+    rval = dc_master.interpolate((int)DataCoupler::VOLUME, "xnew"); RR("Failed to interpolate target solution.");
+      // transfer xNew to coords, for slave
+    rval = write_to_coords(tgt_verts, xNew); RR("Failed writing tag to slave verts.");
+  }
 
   if (debug) {
     std::string str;
@@ -295,13 +304,16 @@ ErrorCode DeformMeshRemap::execute()
       str = "PARALLEL=WRITE_PART";
 #endif
     rval = mbImpl->write_file("smoothed_master.vtk", NULL, str.c_str(), &masterSet, 1);
+
+    if (have_slave) {
 #ifdef USE_MPI
-    str.clear();
-    if (pcSlave && pcSlave->size() > 1) 
-      str = "PARALLEL=WRITE_PART";
+      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);
-  }
+      rval = mbImpl->write_file("slave_interp.vtk", NULL, str.c_str(), &slaveSet, 1);
+    } // if have_slave
+  } // if debug
 
   return MB_SUCCESS;
 }
@@ -394,7 +406,7 @@ int main(int argc, char **argv) {
 
   ProgOptions po("Deformed mesh options");
   po.addOpt<std::string> ("master,m", "Specify the master meshfile name" );
-  po.addOpt<std::string> ("worker,w", "Specify the slave/worker meshfile name" );
+  po.addOpt<std::string> ("worker,w", "Specify the slave/worker meshfile name, or 'none' (no quotes) if master only" );
   po.addOpt<std::string> ("d1,", "Tag name for displacement x or xyz" );
   po.addOpt<std::string> ("d2,", "Tag name for displacement y" );
   po.addOpt<std::string> ("d3,", "Tag name for displacement z" );


https://bitbucket.org/fathomteam/moab/commits/8d9ccfdea33f/
Changeset:   8d9ccfdea33f
Branch:      None
User:        tautges
Date:        2014-03-12 16:53:02
Summary:     Better error handling & debug printing.
In Tree, handle legacy tag on tree root by deleting and re-creating (legacy files have old-sized tag).

Affected #:  2 files

diff --git a/examples/DeformMeshRemap.cpp b/examples/DeformMeshRemap.cpp
index ac3ca39..afb4d9d 100644
--- a/examples/DeformMeshRemap.cpp
+++ b/examples/DeformMeshRemap.cpp
@@ -19,6 +19,10 @@
 #include "MBTagConventions.hpp"
 #include "DataCoupler.hpp"
 
+#define IS_BUILDING_MB
+#include "moab/Error.hpp"
+#undef IS_BUILDING_MB
+
 #ifdef USE_MPI
 #  include "moab/ParallelComm.hpp"
 #endif
@@ -116,6 +120,8 @@ private:
     //! moab interface
   Interface *mbImpl;
 
+  Error *mError;
+  
 #ifdef USE_MPI
     //! ParallelComm for master, slave meshes
   ParallelComm *pcMaster, *pcSlave;
@@ -235,7 +241,7 @@ ErrorCode DeformMeshRemap::execute()
   ErrorCode rval = read_file(MASTER, masterFileName, masterSet);
   if (MB_SUCCESS != rval) return rval;
 
-  if (*solidSetNos[MASTER].begin() == -1 || *fluidSetNos[MASTER].begin() == -1) {
+  if (solidSetNos[MASTER].empty() || fluidSetNos[MASTER].empty()) {
     rval = find_other_sets(MASTER, masterSet); RR("Failed to find other sets in master mesh.");
   }
 
@@ -244,7 +250,7 @@ ErrorCode DeformMeshRemap::execute()
     rval = read_file(SLAVE, slaveFileName, slaveSet);
     if (MB_SUCCESS != rval) return rval;
 
-    if (*solidSetNos[SLAVE].begin() == -1 || *fluidSetNos[SLAVE].begin() == -1) {
+    if (solidSetNos[SLAVE].empty() || fluidSetNos[SLAVE].empty()) {
       rval = find_other_sets(SLAVE, slaveSet); RR("Failed to find other sets in slave mesh.");
     }
   }
@@ -318,50 +324,6 @@ ErrorCode DeformMeshRemap::execute()
   return MB_SUCCESS;
 }
 
-ErrorCode DeformMeshRemap::find_other_sets(int m_or_s, EntityHandle file_set) 
-{
-    // solid or fluid sets are missing; find the other
-  Range *filled_sets = NULL, *unfilled_sets = NULL, *unfilled_elems = NULL;
-  
-  if (fluidSets[m_or_s].empty() && !solidSets[m_or_s].empty()) {
-    unfilled_sets = &fluidSets[m_or_s];
-    filled_sets = &solidSets[m_or_s];
-    unfilled_elems = &fluidElems[m_or_s];
-  }
-  else if (!fluidSets[m_or_s].empty() && solidSets[m_or_s].empty()) {
-    filled_sets = &fluidSets[m_or_s];
-    unfilled_sets = &solidSets[m_or_s];
-    unfilled_elems = &solidElems[m_or_s];
-  }
-  
-    // ok, we know the filled sets, now fill the unfilled sets, and the elems from those
-  Tag tagh;
-  ErrorCode rval = mbImpl->tag_get_handle(MATERIAL_SET_TAG_NAME, tagh); RR("Couldn't get material set tag name.");
-  Range matsets;
-  rval = mbImpl->get_entities_by_type_and_tag(file_set, MBENTITYSET, &tagh, NULL, 1, matsets);
-  if (matsets.empty()) rval = MB_FAILURE;
-  RR("Couldn't get any material sets.");
-  *unfilled_sets = subtract(matsets, *filled_sets);
-  if (unfilled_sets->empty()) {
-    rval = MB_FAILURE;
-    RR("Failed to find any unfilled material sets.");
-  }
-  Range tmp_range;
-  for (Range::iterator rit = unfilled_sets->begin(); rit != unfilled_sets->end(); rit++) {
-    rval = mbImpl->get_entities_by_handle(*rit, tmp_range, true);
-    RR("Failed to get entities in unfilled set.");
-  }
-  int dim = mbImpl->dimension_from_handle(*tmp_range.rbegin());
-  assert(dim > 0 && dim < 4);  
-  *unfilled_elems = tmp_range.subset_by_dimension(dim);
-  if (unfilled_elems->empty()) {
-    rval = MB_FAILURE;
-    RR("Failed to find any unfilled set entities.");
-  }
-  
-  return MB_SUCCESS;
-}
-
 std::string DeformMeshRemap::get_file_name(int m_or_s) const
 {
   switch (m_or_s) {
@@ -390,6 +352,8 @@ 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") 
 {
+  mbImpl->query_interface(mError);
+  
   xDisp[0] = xDisp[1] = xDisp[2] = 0;
 
   if (!pcSlave && pcMaster)
@@ -398,6 +362,7 @@ DeformMeshRemap::DeformMeshRemap(Interface *impl, ParallelComm *master, Parallel
   
 DeformMeshRemap::~DeformMeshRemap() 
 {
+  mbImpl->release_interface(mError);
 }
 
 int main(int argc, char **argv) {
@@ -436,47 +401,29 @@ int main(int argc, char **argv) {
   if(!po.getOpt("worker", &slavef))
     slavef = string(MESH_DIR) + string("/rodtri.g");
   dfr->set_file_name(DeformMeshRemap::SLAVE, slavef);
+  if (slavef.empty()) {
+    std::cerr << "Empty slave file name; if no slave, use filename 'none' (no quotes)." << std::endl;
+    return 1;
+  }
 
   std::vector<int> set_nos;
   po.getOptAllArgs("fm", set_nos);
-  if (set_nos.empty()) {
-    rval = dfr->add_set_no(DeformMeshRemap::MASTER, DeformMeshRemap::FLUID, FLUID_SETNO); 
-    RR("Failed to add solid set no.");
-  }
-  else {
-    for (std::vector<int>::iterator vit = set_nos.begin(); vit != set_nos.end(); vit++) 
-      dfr->add_set_no(DeformMeshRemap::MASTER, DeformMeshRemap::FLUID, *vit);
-  }
+  for (std::vector<int>::iterator vit = set_nos.begin(); vit != set_nos.end(); vit++) 
+    dfr->add_set_no(DeformMeshRemap::MASTER, DeformMeshRemap::FLUID, *vit);
   set_nos.clear();
   
   po.getOptAllArgs("fs", set_nos);
-  if (set_nos.empty()) {
-    rval = dfr->add_set_no(DeformMeshRemap::SLAVE, DeformMeshRemap::FLUID, FLUID_SETNO); 
-    RR("Failed to add solid set no.");
-  }
-  else {
-    for (std::vector<int>::iterator vit = set_nos.begin(); vit != set_nos.end(); vit++) 
-      dfr->add_set_no(DeformMeshRemap::SLAVE, DeformMeshRemap::FLUID, *vit);
-  }
+  for (std::vector<int>::iterator vit = set_nos.begin(); vit != set_nos.end(); vit++) 
+    dfr->add_set_no(DeformMeshRemap::SLAVE, DeformMeshRemap::FLUID, *vit);
   set_nos.clear();
   
   po.getOptAllArgs("sm", set_nos);
-  if (set_nos.empty()) {
-    rval = dfr->add_set_no(DeformMeshRemap::MASTER, DeformMeshRemap::SOLID, SOLID_SETNO); RR("Failed to add solid set no.");
-  }
-  else {
-    for (std::vector<int>::iterator vit = set_nos.begin(); vit != set_nos.end(); vit++) 
-      dfr->add_set_no(DeformMeshRemap::MASTER, DeformMeshRemap::SOLID, *vit);
-  }
+  for (std::vector<int>::iterator vit = set_nos.begin(); vit != set_nos.end(); vit++) 
+    dfr->add_set_no(DeformMeshRemap::MASTER, DeformMeshRemap::SOLID, *vit);
 
   po.getOptAllArgs("ss", set_nos);
-  if (set_nos.empty()) {
-    rval = dfr->add_set_no(DeformMeshRemap::SLAVE, DeformMeshRemap::SOLID, SOLID_SETNO); RR("Failed to add solid set no.");
-  }
-  else {
-    for (std::vector<int>::iterator vit = set_nos.begin(); vit != set_nos.end(); vit++) 
-      dfr->add_set_no(DeformMeshRemap::SLAVE, DeformMeshRemap::SOLID, *vit);
-  }
+  for (std::vector<int>::iterator vit = set_nos.begin(); vit != set_nos.end(); vit++) 
+    dfr->add_set_no(DeformMeshRemap::SLAVE, DeformMeshRemap::SOLID, *vit);
 
   std::string tnames[3];
   po.getOpt("d1", &tnames[0]);
@@ -632,10 +579,13 @@ ErrorCode DeformMeshRemap::read_file(int m_or_s, string &fname, EntityHandle &se
     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);
+    ErrorCode tmp_rval = mbImpl->get_entities_by_type_and_tag(seth, MBENTITYSET, &tagh, &setno_ptr, 1, sets);
+    if (sets.empty() || MB_SUCCESS != tmp_rval) {
+      rval = MB_FAILURE;
+      mError->set_last_error("Couldn't find solid set #%d.\n", *sit);
+    }
+    else
+      solidSets[m_or_s].merge(sets);
   }
 
     // get solid entities, and dimension
@@ -644,19 +594,26 @@ ErrorCode DeformMeshRemap::read_file(int m_or_s, string &fname, EntityHandle &se
     rval = mbImpl->get_entities_by_handle(*rit, tmp_range, true);
     RR("Failed to get entities in solid.");
   }
-  int dim = mbImpl->dimension_from_handle(*tmp_range.rbegin());
-  assert(dim > 0 && dim < 4);
-  
-  solidElems[m_or_s] = tmp_range.subset_by_dimension(dim);
+  if (!tmp_range.empty()) {
+    int dim = mbImpl->dimension_from_handle(*tmp_range.rbegin());
+    assert(dim > 0 && dim < 4);
+    solidElems[m_or_s] = tmp_range.subset_by_dimension(dim);
+  }
+
+  if (debug)
+    std::cout << "Read " << solidElems[m_or_s].size() << " solid elements from " << solidSets[m_or_s].size() << " sets in " << (m_or_s == MASTER ? "master" : "slave") << " mesh." << std::endl;    
 
   for (std::set<int>::iterator sit = fluidSetNos[m_or_s].begin(); sit != fluidSetNos[m_or_s].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);
+    ErrorCode tmp_rval = mbImpl->get_entities_by_type_and_tag(seth, MBENTITYSET, &tagh, &setno_ptr, 1, sets);
+    if (sets.empty() || MB_SUCCESS != tmp_rval) {
+      rval = MB_FAILURE;
+      mError->set_last_error("Couldn't find fluid set #%d.\n", *sit);
+    }
+    else
+      fluidSets[m_or_s].merge(sets);
   }
 
     // get fluid entities, and dimension
@@ -665,8 +622,66 @@ ErrorCode DeformMeshRemap::read_file(int m_or_s, string &fname, EntityHandle &se
     rval = mbImpl->get_entities_by_handle(*rit, tmp_range, true);
     RR("Failed to get entities in fluid.");
   }
+  if (!tmp_range.empty()) {
+    int dim = mbImpl->dimension_from_handle(*tmp_range.rbegin());
+    assert(dim > 0 && dim < 4);
+    fluidElems[m_or_s] = tmp_range.subset_by_dimension(dim);
+  }
   
-  fluidElems[m_or_s] = tmp_range.subset_by_dimension(dim);
+  if (debug)
+    std::cout << "Read " << fluidElems[m_or_s].size() << " fluid elements from " << fluidSets[m_or_s].size() << " sets in " << (m_or_s == MASTER ? "master" : "slave") << " mesh." << std::endl;
+
+  return rval;
+}
+
+ErrorCode DeformMeshRemap::find_other_sets(int m_or_s, EntityHandle file_set) 
+{
+    // solid or fluid sets are missing; find the other
+  Range *filled_sets = NULL, *unfilled_sets = NULL, *unfilled_elems = NULL;
+  
+  if (fluidSets[m_or_s].empty() && !solidSets[m_or_s].empty()) {
+    unfilled_sets = &fluidSets[m_or_s];
+    filled_sets = &solidSets[m_or_s];
+    unfilled_elems = &fluidElems[m_or_s];
+    if (debug)
+      std::cout << "Finding unspecified fluid elements in " << (m_or_s == MASTER ? "master" : "slave") << " mesh...";
+  }
+  else if (!fluidSets[m_or_s].empty() && solidSets[m_or_s].empty()) {
+    filled_sets = &fluidSets[m_or_s];
+    unfilled_sets = &solidSets[m_or_s];
+    unfilled_elems = &solidElems[m_or_s];
+    if (debug)
+      std::cout << "Finding unspecified solid elements in " << (m_or_s == MASTER ? "master" : "slave") << " mesh...";
+  }
+  
+    // ok, we know the filled sets, now fill the unfilled sets, and the elems from those
+  Tag tagh;
+  ErrorCode rval = mbImpl->tag_get_handle(MATERIAL_SET_TAG_NAME, tagh); RR("Couldn't get material set tag name.");
+  Range matsets;
+  rval = mbImpl->get_entities_by_type_and_tag(file_set, MBENTITYSET, &tagh, NULL, 1, matsets);
+  if (matsets.empty()) rval = MB_FAILURE;
+  RR("Couldn't get any material sets.");
+  *unfilled_sets = subtract(matsets, *filled_sets);
+  if (unfilled_sets->empty()) {
+    rval = MB_FAILURE;
+    RR("Failed to find any unfilled material sets.");
+  }
+  Range tmp_range;
+  for (Range::iterator rit = unfilled_sets->begin(); rit != unfilled_sets->end(); rit++) {
+    rval = mbImpl->get_entities_by_handle(*rit, tmp_range, true);
+    RR("Failed to get entities in unfilled set.");
+  }
+  int dim = mbImpl->dimension_from_handle(*tmp_range.rbegin());
+  assert(dim > 0 && dim < 4);  
+  *unfilled_elems = tmp_range.subset_by_dimension(dim);
+  if (unfilled_elems->empty()) {
+    rval = MB_FAILURE;
+    RR("Failed to find any unfilled set entities.");
+  }
+
+  if (debug) 
+    std::cout << "found " << unfilled_sets->size() << " sets and " << unfilled_elems->size() << " elements." << std::endl;
   
   return MB_SUCCESS;
 }
+

diff --git a/src/moab/Tree.hpp b/src/moab/Tree.hpp
index 4d2694a..67b67f9 100644
--- a/src/moab/Tree.hpp
+++ b/src/moab/Tree.hpp
@@ -264,6 +264,14 @@ namespace moab {
       if (!boxTag && create_if_missing) {
         assert(boxTagName.length() > 0);
         ErrorCode rval = moab()->tag_get_handle(boxTagName.c_str(), 6, MB_TYPE_DOUBLE, boxTag, MB_TAG_CREAT | MB_TAG_SPARSE);
+        if (MB_INVALID_SIZE == rval) {
+            // delete the tag and get it again, legacy file...
+          rval = moab()->tag_delete(boxTag);
+          if (MB_SUCCESS != rval) return 0;
+          boxTag = 0;
+          return get_box_tag(true);
+        }
+        
         if (MB_SUCCESS != rval) return 0;
       }
       


https://bitbucket.org/fathomteam/moab/commits/863998bd5433/
Changeset:   863998bd5433
Branch:      None
User:        tautges
Date:        2014-03-12 16:53:02
Summary:     Adding progress-type debug messages.
Was assuming faces in Lloyd smoother, now getting dimension from elements.

Affected #:  2 files

diff --git a/examples/DeformMeshRemap.cpp b/examples/DeformMeshRemap.cpp
index afb4d9d..6357df9 100644
--- a/examples/DeformMeshRemap.cpp
+++ b/examples/DeformMeshRemap.cpp
@@ -254,6 +254,8 @@ ErrorCode DeformMeshRemap::execute()
       rval = find_other_sets(SLAVE, slaveSet); RR("Failed to find other sets in slave mesh.");
     }
   }
+ 
+  if (debug) std::cout << "Constructing data coupler/search tree on master mesh..." << std::endl;
   
   Range src_elems = solidElems[MASTER];
   src_elems.merge(fluidElems[MASTER]);
@@ -271,6 +273,7 @@ ErrorCode DeformMeshRemap::execute()
     RR("Failed to get target verts.");
 
       // locate slave vertices, caching results in dc
+    if (debug) std::cout << "Locating slave vertices in master mesh..." << std::endl;
     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()) {
@@ -281,11 +284,13 @@ ErrorCode DeformMeshRemap::execute()
   }
   
     // deform the master's solid mesh, put results in a new tag
+  if (debug) std::cout << "Deforming fluid elements in master mesh..." << std::endl;
   rval = deform_master(fluidElems[MASTER], solidElems[MASTER], "xnew"); RR("");
 
   { // to isolate the lloyd smoother & delete when done
 
       // smooth the master mesh
+    if (debug) std::cout << "Smoothing fluid elements in master mesh..." << std::endl;
     LloydSmoother ll(mbImpl, NULL, fluidElems[MASTER], xNew);
     rval = ll.perform_smooth();
     RR("Failed in lloyd smoothing.");
@@ -293,13 +298,16 @@ ErrorCode DeformMeshRemap::execute()
   }
 
     // transfer xNew to coords, for master
+  if (debug) std::cout << "Transferring coords tag to vertex coordinates in master mesh..." << std::endl;
   rval = write_to_coords(fluidElems[MASTER], xNew); RR("Failed writing tag to master fluid verts.");
 
   if (have_slave) {
       // map new locations to slave
       // interpolate xNew to slave points
+    if (debug) std::cout << "Interpolating new coordinates to slave vertices..." << std::endl;
     rval = dc_master.interpolate((int)DataCoupler::VOLUME, "xnew"); RR("Failed to interpolate target solution.");
       // transfer xNew to coords, for slave
+    if (debug) std::cout << "Transferring coords tag to vertex coordinates in slave mesh..." << std::endl;
     rval = write_to_coords(tgt_verts, xNew); RR("Failed writing tag to slave verts.");
   }
 
@@ -309,6 +317,7 @@ ErrorCode DeformMeshRemap::execute()
     if (pcMaster && pcMaster->size() > 1) 
       str = "PARALLEL=WRITE_PART";
 #endif
+    if (debug) std::cout << "Writing smoothed_master.vtk..." << std::endl;
     rval = mbImpl->write_file("smoothed_master.vtk", NULL, str.c_str(), &masterSet, 1);
 
     if (have_slave) {
@@ -317,6 +326,7 @@ ErrorCode DeformMeshRemap::execute()
       if (pcSlave && pcSlave->size() > 1) 
         str = "PARALLEL=WRITE_PART";
 #endif
+      if (debug) std::cout << "Writing slave_interp.vtk..." << std::endl;
       rval = mbImpl->write_file("slave_interp.vtk", NULL, str.c_str(), &slaveSet, 1);
     } // if have_slave
   } // if debug

diff --git a/src/LloydSmoother.cpp b/src/LloydSmoother.cpp
index b1614d2..e9cce82 100644
--- a/src/LloydSmoother.cpp
+++ b/src/LloydSmoother.cpp
@@ -39,6 +39,17 @@ ErrorCode LloydSmoother::perform_smooth()
 {
   ErrorCode rval;
 
+  if (myElems.empty()) {
+    rval = MB_FAILURE;
+    RR("No elements specified to Lloyd smoother.");
+  }
+  else if (mbImpl->dimension_from_handle(*myElems.begin()) != mbImpl->dimension_from_handle(*myElems.rbegin())) {
+    rval = MB_FAILURE;
+    RR("Elements of unequal dimension specified to Lloyd smoother.");
+  }    
+
+  int dim = mbImpl->dimension_from_handle(*myElems.begin());
+  
     // first figure out tolerance to use
   if (0 > absTol) {
       // no tolerance set - get one relative to bounding box around elements
@@ -140,7 +151,7 @@ ErrorCode LloydSmoother::perform_smooth()
       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->get_adjacencies(&(*vit), 1, dim, 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};


https://bitbucket.org/fathomteam/moab/commits/bd2eb5c6475b/
Changeset:   bd2eb5c6475b
Branch:      None
User:        tautges
Date:        2014-03-12 16:53:02
Summary:     DeformMeshRemap: more debugging output, and handle coord displacements properly (add to prev coords)
LloydSmoother: fix a bug with fixed tag entity handles size
SpatialLocator: if an element range is specified, build the tree too

Affected #:  3 files

diff --git a/examples/DeformMeshRemap.cpp b/examples/DeformMeshRemap.cpp
index 6357df9..f0d55fa 100644
--- a/examples/DeformMeshRemap.cpp
+++ b/examples/DeformMeshRemap.cpp
@@ -12,6 +12,7 @@
 
 #include "moab/Core.hpp"
 #include "moab/Range.hpp"
+#include "moab/Skinner.hpp"
 #include "moab/LloydSmoother.hpp"
 #include "moab/ProgOptions.hpp"
 #include "moab/BoundBox.hpp"
@@ -288,6 +289,22 @@ ErrorCode DeformMeshRemap::execute()
   rval = deform_master(fluidElems[MASTER], solidElems[MASTER], "xnew"); RR("");
 
   { // to isolate the lloyd smoother & delete when done
+    if (debug) {
+        // output the skin of smoothed elems, as a check
+        // get the skin; get facets, because we might need to filter on shared entities
+      Skinner skinner(mbImpl);
+      Range skin;
+      rval = skinner.find_skin(0, fluidElems[MASTER], false, skin); RR("Unable to find skin.");
+      EntityHandle skin_set;
+      std::cout << "Writing skin_mesh.g and fluid_mesh.g." << std::endl;
+      rval = mbImpl->create_meshset(MESHSET_SET, skin_set); RR("Failed to create skin set.");
+      rval = mbImpl->add_entities(skin_set, skin); RR("Failed to add skin entities to set.");
+      rval = mbImpl->write_file("skin_mesh.vtk", NULL, NULL, &skin_set, 1); RR("Failure to write skin set.");
+      rval = mbImpl->remove_entities(skin_set, skin); RR("Failed to remove skin entities from set.");
+      rval = mbImpl->add_entities(skin_set, fluidElems[MASTER]); RR("Failed to add fluid entities to set.");
+      rval = mbImpl->write_file("fluid_mesh.vtk", NULL, NULL, &skin_set, 1); RR("Failure to write fluid set.");
+      rval = mbImpl->delete_entities(&skin_set, 1); RR("Failed to delete skin set.");
+    }
 
       // smooth the master mesh
     if (debug) std::cout << "Smoothing fluid elements in master mesh..." << std::endl;
@@ -385,10 +402,10 @@ int main(int argc, char **argv) {
   po.addOpt<std::string> ("d1,", "Tag name for displacement x or xyz" );
   po.addOpt<std::string> ("d2,", "Tag name for displacement y" );
   po.addOpt<std::string> ("d3,", "Tag name for displacement z" );
-  po.addOpt<int> ("fm,", "Specify master fluid material set number(s). If -1 specified, fluid sets derived from complement of solid sets.");
-  po.addOpt<int> ("fs,", "Specify master solid material set number(s). If -1 specified, fluid sets derived from complement of solid sets.");
-  po.addOpt<int> ("sm,", "Specify slave fluid material set number(s). If -1 specified, fluid sets derived from complement of solid sets.");
-  po.addOpt<int> ("ss,", "Specify slave solid material set number(s). If -1 specified, fluid sets derived from complement of solid sets.");
+  po.addOpt<int> ("fm,", "Specify master fluid material set number(s). If none specified, fluid sets derived from complement of solid sets.");
+  po.addOpt<int> ("fs,", "Specify master solid material set number(s). If none specified, solid sets derived from complement of fluid sets.");
+  po.addOpt<int> ("sm,", "Specify slave fluid material set number(s). If none specified, fluid sets derived from complement of solid sets.");
+  po.addOpt<int> ("ss,", "Specify slave solid material set number(s). If none specified, solid sets derived from complement of fluid sets.");
   
   po.parseCommandLine(argc, argv);
 
@@ -500,7 +517,7 @@ ErrorCode DeformMeshRemap::deform_master(Range &fluid_elems, Range &solid_elems,
   Range verts;
   rval = mbImpl->get_adjacencies(solid_elems, 0, false, verts, Interface::UNION);
   RR("Failed to get vertices.");
-  std::vector<double> coords(3*verts.size());
+  std::vector<double> coords(3*verts.size()), new_coords(3*verts.size());
   rval = mbImpl->get_coords(verts, &coords[0]);
   RR("Failed to get vertex coords.");
   unsigned int num_verts = verts.size();
@@ -512,13 +529,14 @@ ErrorCode DeformMeshRemap::deform_master(Range &fluid_elems, Range &solid_elems,
     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.");
     std::vector<double> disps(num_verts);
+    
     for (int i = 0; i < 3; i++) {
       rval = mbImpl->tag_get_handle(xDispNames[0].c_str(), 1, MB_TYPE_DOUBLE, xDisp[i]);
       RR("Failed to get xDisp tag.");
       rval = mbImpl->tag_get_data(xDisp[i], verts, &disps[0]);
       RR("Failed to get xDisp tag values.");
       for (unsigned int j = 0; j < num_verts; j++)
-        coords[3*j+i] = disps[j];
+        new_coords[3*j+i] = coords[3*j+i] + disps[j];
     }
   }
   else if (!xDispNames[0].empty()) {
@@ -528,7 +546,7 @@ ErrorCode DeformMeshRemap::deform_master(Range &fluid_elems, Range &solid_elems,
     std::vector<double> disps(3*num_verts);
     rval = mbImpl->tag_get_data(xDisp[0], verts, &disps[0]);
     for (unsigned int j = 0; j < 3*num_verts; j++)
-      coords[j] += disps[j];
+      new_coords[j] = coords[j] + disps[j];
   }
   else {
       // get the bounding box of the solid mesh
@@ -536,9 +554,19 @@ ErrorCode DeformMeshRemap::deform_master(Range &fluid_elems, Range &solid_elems,
     bbox.update(*mbImpl, solid_elems);
   
     for (unsigned int j = 0; j < num_verts; j++)
-      deform_func(bbox, &coords[3*j], &coords[3*j]);
+      deform_func(bbox, &coords[3*j], &new_coords[3*j]);
   }
 
+  if (debug) {
+    double len = 0.0;
+    for (unsigned int i = 0; i < num_verts; i++) {
+      CartVect dx = CartVect(&new_coords[3*i]) - CartVect(&coords[3*i]);
+      double tmp_len = dx.length_squared();
+      if (tmp_len > len) len = tmp_len;
+    }
+    std::cout << "Max displacement = " << len << std::endl;
+  }
+  
   if (!xNew) {
     rval = mbImpl->tag_get_handle((tag_name ? tag_name : ""), 3, MB_TYPE_DOUBLE, 
                                   xDisp[0], MB_TAG_CREAT|MB_TAG_DENSE);

diff --git a/src/LloydSmoother.cpp b/src/LloydSmoother.cpp
index e9cce82..ba8169a 100644
--- a/src/LloydSmoother.cpp
+++ b/src/LloydSmoother.cpp
@@ -232,7 +232,7 @@ ErrorCode LloydSmoother::initialize()
     RR("Trouble getting vertices.");
     
       // mark them fixed
-    std::vector<unsigned char> marks(skin.size(), 0x1);
+    std::vector<unsigned char> marks(skin_verts.size(), 0x1);
     rval = mbImpl->tag_set_data(fixedTag, skin_verts, &marks[0]);
     RR("Unable to set tag on skin.");
   }

diff --git a/src/SpatialLocator.cpp b/src/SpatialLocator.cpp
index f976ae7..b2dc8e3 100644
--- a/src/SpatialLocator.cpp
+++ b/src/SpatialLocator.cpp
@@ -40,7 +40,9 @@ namespace moab
   
       myDim = mbImpl->dimension_from_handle(*elems.begin());
       myElems = elems;
-      return MB_SUCCESS;
+
+      ErrorCode rval = myTree->build_tree(myElems);
+      return rval;
     }
     
 #ifdef USE_MPI


https://bitbucket.org/fathomteam/moab/commits/7d59aad7b01d/
Changeset:   7d59aad7b01d
Branch:      None
User:        tautges
Date:        2014-03-12 16:53:02
Summary:     Adding % change to max displacement debug output.

Affected #:  1 file

diff --git a/examples/DeformMeshRemap.cpp b/examples/DeformMeshRemap.cpp
index f0d55fa..667904b 100644
--- a/examples/DeformMeshRemap.cpp
+++ b/examples/DeformMeshRemap.cpp
@@ -564,7 +564,13 @@ ErrorCode DeformMeshRemap::deform_master(Range &fluid_elems, Range &solid_elems,
       double tmp_len = dx.length_squared();
       if (tmp_len > len) len = tmp_len;
     }
-    std::cout << "Max displacement = " << len << std::endl;
+    Range tmp_elems(fluid_elems);
+    tmp_elems.merge(solid_elems);
+    BoundBox box;
+    box.update(*mbImpl, tmp_elems);
+    double max_len = std::max(box.bMax[2]-box.bMin[2], std::max(box.bMax[1]-box.bMin[1], box.bMax[0]-box.bMin[0]));
+    
+    std::cout << "Max displacement = " << len << " (" << 100.0 * len / max_len << "% of max box length)" << std::endl;
   }
   
   if (!xNew) {


https://bitbucket.org/fathomteam/moab/commits/cd923766bcc7/
Changeset:   cd923766bcc7
Branch:      None
User:        tautges
Date:        2014-03-12 16:53:41
Summary:     DeformMeshRemap: add some tree statistics output
AdaptiveKDTree: make several functions non-static so they show up in profiling; also use BoundBox a bit more
BVHTree: account for BoundBox change
SpatialLocator: choose tree type based on contents of myElems, and move tree construction to a separate function
BoundBox: change name of contains_box to intersects_box, and fix logic that was broken there
TreeStats.hpp: improve statistics output
test_boundbox: fix tests for changes to BoundBox

Affected #:  8 files

diff --git a/examples/DeformMeshRemap.cpp b/examples/DeformMeshRemap.cpp
index 667904b..f752b4e 100644
--- a/examples/DeformMeshRemap.cpp
+++ b/examples/DeformMeshRemap.cpp
@@ -348,6 +348,9 @@ ErrorCode DeformMeshRemap::execute()
     } // if have_slave
   } // if debug
 
+  if (debug) 
+    dc_master.spatial_locator()->get_tree()->tree_stats().print();
+
   return MB_SUCCESS;
 }
 

diff --git a/src/AdaptiveKDTree.cpp b/src/AdaptiveKDTree.cpp
index 1324c30..8aafdd1 100644
--- a/src/AdaptiveKDTree.cpp
+++ b/src/AdaptiveKDTree.cpp
@@ -176,7 +176,6 @@ namespace moab {
         else {
           rval = iter.step();
           if (MB_ENTITY_NOT_FOUND == rval) {
-            treeStats.reset();
             rval = treeStats.compute_stats(mbImpl, myRoot);
             treeStats.initTime = cp.time_elapsed();
             return rval;  // at end
@@ -807,8 +806,7 @@ namespace moab {
                                           t_enter, t_exit );
     }
 
-    static ErrorCode intersect_children_with_elems(
-        AdaptiveKDTree* tool,
+    ErrorCode AdaptiveKDTree::intersect_children_with_elems(
         const Range& elems,
         AdaptiveKDTree::Plane plane,
         double eps,
@@ -823,15 +821,16 @@ namespace moab {
       right_tris.clear();
       both_tris.clear();
       CartVect coords[16];
-      Interface *const moab = tool->moab();
   
         // get extents of boxes for left and right sides
-      CartVect right_min( box_min ), left_max( box_max );
-      right_min[plane.norm] = left_max[plane.norm] = plane.coord;
-      const CartVect left_cen = 0.5*(left_max + box_min);
-      const CartVect left_dim = 0.5*(left_max - box_min);
-      const CartVect right_cen = 0.5*(box_max + right_min);
-      const CartVect right_dim = 0.5*(box_max - right_min);
+      BoundBox left_box(box_min, box_max), right_box(box_min, box_max);
+      right_box.bMin = box_min;
+      left_box.bMax = box_max;
+      right_box.bMin[plane.norm] = left_box.bMax[plane.norm] = plane.coord;
+      const CartVect left_cen = 0.5*(left_box.bMax + box_min);
+      const CartVect left_dim = 0.5*(left_box.bMax - box_min);
+      const CartVect right_cen = 0.5*(box_max + right_box.bMin);
+      const CartVect right_dim = 0.5*(box_max - right_box.bMin);
       const CartVect dim = box_max - box_min;
       const double max_tol = std::max(dim[0], std::max(dim[1], dim[2]))/10;
   
@@ -851,8 +850,8 @@ namespace moab {
   
         // vertices
       for (i = elems.begin(); i != elem_begin; ++i) {
-        tool->tree_stats().leafObjectTests++;
-        rval = moab->get_coords( &*i, 1, coords[0].array() );
+        tree_stats().leafObjectTests++;
+        rval = moab()->get_coords( &*i, 1, coords[0].array() );
         if (MB_SUCCESS != rval)
           return rval;
     
@@ -873,13 +872,13 @@ namespace moab {
         // non-polyhedron elements
       std::vector<EntityHandle> dum_vector;
       for (i = elem_begin; i != poly_begin; ++i) {
-        tool->tree_stats().leafObjectTests++;
-        rval = moab->get_connectivity( *i, conn, count, true, &dum_vector);
+        tree_stats().leafObjectTests++;
+        rval = moab()->get_connectivity( *i, conn, count, true, &dum_vector);
         if (MB_SUCCESS != rval) 
           return rval;
         if (count > (int)(sizeof(coords)/sizeof(coords[0])))
           return MB_FAILURE;
-        rval = moab->get_coords( &conn[0], count, coords[0].array() );
+        rval = moab()->get_coords( &conn[0], count, coords[0].array() );
         if (MB_SUCCESS != rval) return rval;
     
         bool lo = false, ro = false;
@@ -894,12 +893,16 @@ namespace moab {
           // identified that leaf, then we're done.  If triangle is on both
           // sides of plane, do more precise test to ensure that it is really
           // in both.
+//        BoundBox box;
+//        box.update(*moab(), *i);
         if (lo && ro) {
           double tol = eps;
           lo = ro = false;
           while (!lo && !ro && tol <= max_tol) {
-            lo = GeomUtil::box_elem_overlap( coords, TYPE_FROM_HANDLE(*i), left_cen, left_dim+CartVect(tol) );
-            ro = GeomUtil::box_elem_overlap( coords, TYPE_FROM_HANDLE(*i),right_cen,right_dim+CartVect(tol) );
+            tree_stats().boxElemTests+= 2;
+            lo = GeomUtil::box_elem_overlap( coords, TYPE_FROM_HANDLE(*i), left_cen, left_dim+CartVect(tol));
+            ro = GeomUtil::box_elem_overlap( coords, TYPE_FROM_HANDLE(*i), right_cen, right_dim+CartVect(tol));
+            
             tol *= 10.0;
           }
         }
@@ -913,20 +916,20 @@ namespace moab {
   
         // polyhedra
       for (i = poly_begin; i != set_begin; ++i) {
-        tool->tree_stats().leafObjectTests++;
-        rval = moab->get_connectivity( *i, conn, count, true );
+        tree_stats().leafObjectTests++;
+        rval = moab()->get_connectivity( *i, conn, count, true );
         if (MB_SUCCESS != rval) 
           return rval;
       
           // just check the bounding box of the polyhedron
         bool lo = false, ro = false;
         for (int j = 0; j < count; ++j) {
-          rval = moab->get_connectivity( conn[j], conn2, count2, true );
+          rval = moab()->get_connectivity( conn[j], conn2, count2, true );
           if (MB_SUCCESS != rval)
             return rval;
       
           for (int k = 0; k < count2; ++k) {
-            rval = moab->get_coords( conn2 + k, 1, coords[0].array() );
+            rval = moab()->get_coords( conn2 + k, 1, coords[0].array() );
             if (MB_SUCCESS != rval)
               return rval;
             if (coords[0][plane.norm] <= plane.coord)
@@ -947,8 +950,8 @@ namespace moab {
         // sets
       BoundBox tbox;
       for (i = set_begin; i != elems.end(); ++i) {
-        tool->tree_stats().leafObjectTests++;
-        rval = tbox.update(*tool->moab(), *i);
+        tree_stats().leafObjectTests++;
+        rval = tbox.update(*moab(), *i);
         if (MB_SUCCESS != rval)
           return rval;
     
@@ -1005,8 +1008,7 @@ namespace moab {
           AdaptiveKDTree::Plane plane = { box_min[axis] + (p/(1.0+plane_count)) * diff[axis], axis };
           Range left, right, both;
           double val;
-          r = intersect_children_with_elems( iter.tool(),
-                                             entities, plane, eps,
+          r = intersect_children_with_elems( entities, plane, eps,
                                              box_min, box_max,
                                              left, right, both, 
                                              val );
@@ -1057,32 +1059,38 @@ namespace moab {
       if (MB_SUCCESS != r)
         return r;
 
-      tmp_data.resize( vertices.size() );
+      unsigned int nverts = vertices.size();
+      tmp_data.resize( 3*nverts);
+      r = iter.tool()->moab()->get_coords( vertices, &tmp_data[0], &tmp_data[nverts], &tmp_data[2*nverts] );
+      if (MB_SUCCESS != r)
+        return r;
+  
       for (int axis = 0; axis < 3; ++axis) {
         int plane_count = num_planes;
+
+          // if num_planes results in width < eps, reset the plane count
         if ((num_planes+1)*eps >= diff[axis])
           plane_count = (int)(diff[axis] / eps) - 1;
 
-        double *ptrs[] = { 0, 0, 0 };
-        ptrs[axis] = &tmp_data[0];
-        r = iter.tool()->moab()->get_coords( vertices, ptrs[0], ptrs[1], ptrs[2] );
-        if (MB_SUCCESS != r)
-          return r;
-  
         for (int p = 1; p <= plane_count; ++p) {
+
+            // coord of this plane on axis
           double coord = box_min[axis] + (p/(1.0+plane_count)) * diff[axis];
-          double closest_coord = tmp_data[0];
-          for (unsigned i = 1; i < tmp_data.size(); ++i) 
-            if (fabs(coord-tmp_data[i]) < fabs(coord-closest_coord))
-              closest_coord = tmp_data[i];
+
+            // find closest vertex coordinate to this plane position
+          unsigned int istrt = axis*nverts;
+          double closest_coord = tmp_data[istrt];
+          for (unsigned i = 1; i < nverts; ++i) 
+            if (fabs(coord-tmp_data[istrt+i]) < fabs(coord-closest_coord))
+              closest_coord = tmp_data[istrt+i];
           if (closest_coord - box_min[axis] <= eps || box_max[axis] - closest_coord <= eps)
             continue;
           
+            // seprate elems into left/right/both, and compute separating metric
           AdaptiveKDTree::Plane plane = { closest_coord, axis };
           Range left, right, both;
           double val;
-          r = intersect_children_with_elems( iter.tool(),
-                                             entities, plane, eps,
+          r = intersect_children_with_elems( entities, plane, eps,
                                              box_min, box_max,
                                              left, right, both, 
                                              val );
@@ -1160,8 +1168,7 @@ namespace moab {
           AdaptiveKDTree::Plane plane = { *citer, axis };
           Range left, right, both;
           double val;
-          r = intersect_children_with_elems( iter.tool(),
-                                             entities, plane, eps,
+          r = intersect_children_with_elems( entities, plane, eps,
                                              box_min, box_max,
                                              left, right, both, 
                                              val );
@@ -1276,8 +1283,7 @@ namespace moab {
           AdaptiveKDTree::Plane plane = { coords[indices[p]], axis };
           Range left, right, both;
           double val;
-          r = intersect_children_with_elems( iter.tool(),
-                                             entities, plane, eps,
+          r = intersect_children_with_elems( entities, plane, eps,
                                              box_min, box_max,
                                              left, right, both, 
                                              val );

diff --git a/src/BVHTree.cpp b/src/BVHTree.cpp
index d54e24d..e2592fb 100644
--- a/src/BVHTree.cpp
+++ b/src/BVHTree.cpp
@@ -48,7 +48,7 @@ namespace moab
 
 #ifndef NDEBUG
       for(std::vector<HandleData>::const_iterator i = handle_data_vec.begin(); i != handle_data_vec.end(); ++i) {
-        if(!boundBox.contains_box(i->myBox, 0)){
+        if(!boundBox.intersects_box(i->myBox, 0)){
           std::cerr << "BB:" << boundBox << "EB:" << i->myBox << std::endl;
           return MB_FAILURE;
         }
@@ -168,16 +168,16 @@ namespace moab
         for (unsigned int dim = 0; dim < 3; ++dim){
           const unsigned int index = Bucket::bucket_index(splitsPerDir, box, interval, dim);
           Bucket &bucket = buckets[dim][index];
-          if(!bucket.boundingBox.contains_box(box))
+          if(!bucket.boundingBox.intersects_box(box))
             std::cerr << "Buckets not covering elements!" << std::endl;
         }
       }
-      if(!elt_union.contains_box(interval) ){
+      if(!elt_union.intersects_box(interval) ){
         std::cout << "element union: " << std::endl << elt_union; 
         std::cout << "intervals: " << std::endl << interval;
         std::cout << "union of elts does not contain original box!" << std::endl;
       }
-      if (!interval.contains_box(elt_union) ){
+      if (!interval.intersects_box(elt_union) ){
         std::cout << "original box does not contain union of elts" << std::endl;
         std::cout << interval << std::endl << elt_union << std::endl;
       }
@@ -193,9 +193,9 @@ namespace moab
         for(unsigned int i = 0; i < nonempty.size(); ++i)
           test_box.update(buckets_[nonempty[i]].boundingBox);
 
-        if(!test_box.contains_box(interval) )
+        if(!test_box.intersects_box(interval) )
           std::cout << "union of buckets in dimension: " << d << "does not contain original box!" << std::endl;
-        if (!interval.contains_box(test_box) ) {
+        if (!interval.intersects_box(test_box) ) {
           std::cout << "original box does " << "not contain union of buckets" 
                     << "in dimension: " << d << std::endl;
           std::cout << interval << std::endl << test_box << std::endl;
@@ -234,7 +234,7 @@ namespace moab
             s != splits_end; ++s, ++left_end) {
           BoundBox test_box = s->leftBox;
           test_box.update(s->rightBox);
-          if(!data.boundingBox.contains_box(test_box)) {
+          if(!data.boundingBox.intersects_box(test_box)) {
             std::cout << "nr: " << s->nr << std::endl;
             std::cout << "Test box: " << std::endl << 
                 test_box;
@@ -247,7 +247,7 @@ namespace moab
             std::cout << "Split boxes larger than bb" 
                       << std::endl;
           }
-          if(!test_box.contains_box(data.boundingBox)) {
+          if(!test_box.intersects_box(data.boundingBox)) {
             std::cout << "bb larger than union of split boxes" << std::endl;
           }         	
         }
@@ -288,7 +288,7 @@ namespace moab
       data.Lmax = data.leftBox.bMax[data.dim];
 #ifndef NDEBUG
       BoundBox test_box(data.rightBox);
-      if(!data.boundingBox.contains_box(test_box)) {
+      if(!data.boundingBox.intersects_box(test_box)) {
         std::cerr << "MEDIAN: BB Does not contain splits" << std::endl;
         std::cerr << "test_box:         " << test_box << std::endl;
         std::cerr << "data.boundingBox: " << data.boundingBox << std::endl;
@@ -335,7 +335,7 @@ namespace moab
           else {
             right_box.update(i->myBox);
           }
-          if(!right_box.contains_box(i->myBox)) {
+          if(!right_box.intersects_box(i->myBox)) {
             if(!issue) {
               std::cerr << "Bounding right box issue!" 
                         << std::endl;
@@ -352,7 +352,7 @@ namespace moab
           else {
             left_box.update(i->myBox);
           }
-          if(!data.leftBox.contains_box(i->myBox)) {
+          if(!data.leftBox.intersects_box(i->myBox)) {
             if(!issue) {
               std::cerr << "Bounding left box issue!" 
                         << std::endl;
@@ -367,13 +367,13 @@ namespace moab
           }
         }
       }
-      if(!left_box.contains_box(data.leftBox)) 
+      if(!left_box.intersects_box(data.leftBox)) 
         std::cout << "left elts do not contain left box" << std::endl;
-      if(!data.leftBox.contains_box(left_box)) 
+      if(!data.leftBox.intersects_box(left_box)) 
         std::cout << "left box does not contain left elts" << std::endl;
-      if(!right_box.contains_box(data.rightBox))
+      if(!right_box.intersects_box(data.rightBox))
         std::cout << "right elts do not contain right box" << std::endl;
-      if(!data.rightBox.contains_box(right_box))
+      if(!data.rightBox.intersects_box(right_box))
         std::cout << "right box do not contain right elts" << std::endl;
 
       if(count_left != data.nl || count_right != data.nr) {
@@ -400,7 +400,7 @@ namespace moab
     {
 #ifndef NDEBUG
       for(HandleDataVec::const_iterator i = begin; i != end; ++i) {
-        if(!box.contains_box(i->myBox, 0)) {
+        if(!box.intersects_box(i->myBox, 0)) {
           std::cerr << "depth: " << depth << std::endl;
           std::cerr << "BB:" << box << "EB:" << i->myBox << std::endl;
           std::exit(-1);

diff --git a/src/SpatialLocator.cpp b/src/SpatialLocator.cpp
index b2dc8e3..2903155 100644
--- a/src/SpatialLocator.cpp
+++ b/src/SpatialLocator.cpp
@@ -2,6 +2,7 @@
 #include "moab/Interface.hpp"
 #include "moab/ElemEvaluator.hpp"
 #include "moab/AdaptiveKDTree.hpp"
+#include "moab/BVHTree.hpp"
 
 // include ScdInterface for box partitioning
 #include "moab/ScdInterface.hpp"
@@ -18,10 +19,8 @@ namespace moab
     SpatialLocator::SpatialLocator(Interface *impl, Range &elems, Tree *tree, ElemEvaluator *eval) 
             : mbImpl(impl), myElems(elems), myDim(-1), myTree(tree), elemEval(eval), iCreatedTree(false)
     {
-      if (!myTree) {
-        myTree = new AdaptiveKDTree(impl);
-        iCreatedTree = true;
-      }
+      create_tree();
+      
       if (!elems.empty()) {
         myDim = mbImpl->dimension_from_handle(*elems.rbegin());
         ErrorCode rval = myTree->build_tree(myElems);
@@ -32,6 +31,20 @@ namespace moab
       }
     }
 
+    void SpatialLocator::create_tree() 
+    {
+      if (myTree) return;
+      
+      if (myElems.empty() || mbImpl->type_from_handle(*myElems.rbegin()) == MBVERTEX) 
+          // create a kdtree if only vertices
+        myTree = new AdaptiveKDTree(mbImpl);
+      else
+          // otherwise a BVHtree, since it performs better for elements
+        myTree = new BVHTree(mbImpl);
+
+      iCreatedTree = true;
+    }
+
     ErrorCode SpatialLocator::add_elems(Range &elems) 
     {
       if (elems.empty() ||

diff --git a/src/moab/AdaptiveKDTree.hpp b/src/moab/AdaptiveKDTree.hpp
index c77bdbc..450a418 100644
--- a/src/moab/AdaptiveKDTree.hpp
+++ b/src/moab/AdaptiveKDTree.hpp
@@ -251,7 +251,6 @@ namespace moab {
       virtual ErrorCode print();
       
   private:
-
       friend class AdaptiveKDTreeIter;
 
         /** \brief Parse options for tree creation
@@ -273,7 +272,18 @@ namespace moab {
                           DataType type, int count, void* default_val, Tag& tag_handle,
                           std::vector<Tag>& created_tags );
   
-      static ErrorCode best_subdivision_snap_plane( int num_planes,
+      ErrorCode intersect_children_with_elems(
+          const Range& elems,
+          AdaptiveKDTree::Plane plane,
+          double eps,
+          CartVect box_min,
+          CartVect box_max,
+          Range& left_tris,
+          Range& right_tris,
+          Range& both_tris,
+          double& metric_value );
+
+      ErrorCode best_subdivision_snap_plane( int num_planes,
                                                     const AdaptiveKDTreeIter& iter,
                                                     Range& best_left,
                                                     Range& best_right,
@@ -282,7 +292,7 @@ namespace moab {
                                                     std::vector<double>& tmp_data,
                                                     double eps );
   
-      static ErrorCode best_subdivision_plane( int num_planes,
+      ErrorCode best_subdivision_plane( int num_planes,
                                                const AdaptiveKDTreeIter& iter,
                                                Range& best_left,
                                                Range& best_right,
@@ -290,7 +300,7 @@ namespace moab {
                                                AdaptiveKDTree::Plane& best_plane,
                                                double eps );
   
-      static ErrorCode best_vertex_median_plane( int num_planes,
+      ErrorCode best_vertex_median_plane( int num_planes,
                                                  const AdaptiveKDTreeIter& iter,
                                                  Range& best_left,
                                                  Range& best_right,
@@ -299,7 +309,7 @@ namespace moab {
                                                  std::vector<double>& coords,
                                                  double eps);
   
-      static ErrorCode best_vertex_sample_plane( int num_planes,
+      ErrorCode best_vertex_sample_plane( int num_planes,
                                                  const AdaptiveKDTreeIter& iter,
                                                  Range& best_left,
                                                  Range& best_right,

diff --git a/src/moab/BoundBox.hpp b/src/moab/BoundBox.hpp
index 3e525ac..c74707a 100644
--- a/src/moab/BoundBox.hpp
+++ b/src/moab/BoundBox.hpp
@@ -17,7 +17,7 @@ namespace moab {
       ~BoundBox() {}
 
       bool contains_point(const double *point, const double tol = 0.0) const;
-      bool contains_box(const BoundBox &b, const double tol = 0.0) const;
+      bool intersects_box(const BoundBox &b, const double tol = 0.0) const;
       void compute_center(CartVect &center);
       void update(const BoundBox &other_box);
       void update(const double *coords);
@@ -75,10 +75,10 @@ namespace moab {
       else return true;
     }
 
-    inline bool BoundBox::contains_box(const BoundBox &b, const double tol) const {
-      if (b.bMin[0] < bMin[0]-tol || b.bMax[0] > bMax[0]+tol ||
-          b.bMin[1] < bMin[1]-tol || b.bMax[1] > bMax[1]+tol ||
-          b.bMin[2] < bMin[2]-tol || b.bMax[2] > bMax[2]+tol) 
+    inline bool BoundBox::intersects_box(const BoundBox &b, const double tol) const {
+      if (b.bMax[0] < bMin[0]-tol || b.bMin[0] > bMax[0]+tol ||
+          b.bMax[1] < bMin[1]-tol || b.bMin[1] > bMax[1]+tol ||
+          b.bMax[2] < bMin[2]-tol || b.bMin[2] > bMax[2]+tol) 
         return false;
 
       else return true;

diff --git a/src/moab/TreeStats.hpp b/src/moab/TreeStats.hpp
index be619a5..4ae685b 100644
--- a/src/moab/TreeStats.hpp
+++ b/src/moab/TreeStats.hpp
@@ -59,6 +59,7 @@ namespace moab
       unsigned int leavesVisited;
       unsigned int numTraversals;
       unsigned int leafObjectTests;
+      unsigned int boxElemTests;
 
   private:
       ErrorCode traverse(Interface *impl, EntityHandle node, unsigned int &depth);
@@ -68,6 +69,8 @@ namespace moab
     inline ErrorCode TreeStats::compute_stats(Interface *impl, EntityHandle root_node) 
     {
       maxDepth = 0;
+      numNodes = 0;
+      numLeaves = 0;
       return traverse(impl, root_node, maxDepth);
     }
       
@@ -111,21 +114,27 @@ namespace moab
       leavesVisited = 0;
       numTraversals = 0;
       leafObjectTests = 0;
+      boxElemTests = 0;
     }
     
     inline void TreeStats::print() const {
       std::cout << "Tree initialization time = " << initTime << " seconds" << std::endl;
       
+      std::cout << "Num nodes         = " << numNodes << std::endl;
+      std::cout << "Num leaves        = " << numLeaves << std::endl;
+      std::cout << "Max depth         = " << maxDepth << std::endl << std::endl;
+
       std::cout << "NodesVisited      = " << nodesVisited << std::endl;
       std::cout << "LeavesVisited     = " << leavesVisited << std::endl;
       std::cout << "Num Traversals    = " << numTraversals << std::endl;
       std::cout << "Leaf Object Tests = " << leafObjectTests << std::endl;
+      std::cout << "Box-Element Tests = " << boxElemTests << std::endl;
     }
 
     inline void TreeStats::output() const 
     {
-      std::cout << initTime << " " << nodesVisited << " " << leavesVisited << " " << numTraversals << " " << leafObjectTests
-                << " # initTime, nodesVisited, leavesVisited, numTraversals, leafObjectTests" << std::endl;
+      std::cout << initTime << " " << nodesVisited << " " << leavesVisited << " " << numTraversals << " " << leafObjectTests << " " << boxElemTests
+                << " # initTime, nodesVisited, leavesVisited, numTraversals, leafObjectTests, boxElemTests" << std::endl;
     }
 }
 

diff --git a/test/test_boundbox.cpp b/test/test_boundbox.cpp
index 9fe0357..b089972 100644
--- a/test/test_boundbox.cpp
+++ b/test/test_boundbox.cpp
@@ -35,7 +35,7 @@ void test_bound_box()
     // test for contains point failure
   bool result = box.contains_point(vec.array(), tol);
   CHECK(!result);
-  result = box.contains_box(other_box, tol);
+  result = box.intersects_box(other_box, tol);
   CHECK(!result);
   
   box = other_box;
@@ -44,7 +44,7 @@ void test_bound_box()
     // test for success
   result = box.contains_point(&vals[0], tol);
   CHECK(result);
-  result = box.contains_box(other_box, tol);
+  result = box.intersects_box(other_box, tol);
   CHECK(result);
   
     // check update functions
@@ -54,9 +54,9 @@ void test_bound_box()
   CHECK(result);
   result = box.contains_point((three*1.1).array(), tol);
   CHECK(!result);
-  result = box.contains_box(BoundBox(three, three), tol);
+  result = box.intersects_box(BoundBox(three, three), tol);
   CHECK(result);
-  result = box.contains_box(BoundBox(three, 1.1*three), tol);
+  result = box.intersects_box(BoundBox(1.1*three, 3.0*three), tol);
   CHECK(!result);
   
   CartVect negthree(-3.0);
@@ -65,9 +65,9 @@ void test_bound_box()
   CHECK(result);
   result = box.contains_point((negthree*1.1).array(), tol);
   CHECK(!result);
-  result = box.contains_box(BoundBox(negthree, negthree), tol);
+  result = box.intersects_box(BoundBox(negthree, negthree), tol);
   CHECK(result);
-  result = box.contains_box(BoundBox(1.1*negthree, negthree), tol);
+  result = box.intersects_box(BoundBox(3.0*negthree, 1.1*negthree), tol);
   CHECK(!result);
   
   for (int i = 0; i < 3; i++) {
@@ -78,9 +78,9 @@ void test_bound_box()
   CHECK(result);
   result = box.contains_point((CartVect(&vals[0])*1.1).array(), tol);
   CHECK(!result);
-  result = box.contains_box(BoundBox(&vals[0]), tol);
+  result = box.intersects_box(BoundBox(&vals[0]), tol);
   CHECK(result);
-  result = box.contains_box(BoundBox(1.1*CartVect(&vals[0]), CartVect(&vals[3])), tol);
+  result = box.intersects_box(BoundBox(1.2*CartVect(&vals[0]), 1.1*CartVect(&vals[0])), tol);
   CHECK(!result);
   
     // check length functions


https://bitbucket.org/fathomteam/moab/commits/da1daa3012a9/
Changeset:   da1daa3012a9
Branch:      None
User:        tautges
Date:        2014-03-12 16:54:58
Summary:     Changes to allow specification of existing displacement tag(s), which can be either single name corresponding to xyz[3]
or 3 names x, y, z.  Handle either case in deform_master.

Affected #:  1 file

diff --git a/examples/DeformMeshRemap.cpp b/examples/DeformMeshRemap.cpp
index f752b4e..109fd3b 100644
--- a/examples/DeformMeshRemap.cpp
+++ b/examples/DeformMeshRemap.cpp
@@ -532,7 +532,6 @@ ErrorCode DeformMeshRemap::deform_master(Range &fluid_elems, Range &solid_elems,
     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.");
     std::vector<double> disps(num_verts);
-    
     for (int i = 0; i < 3; i++) {
       rval = mbImpl->tag_get_handle(xDispNames[0].c_str(), 1, MB_TYPE_DOUBLE, xDisp[i]);
       RR("Failed to get xDisp tag.");


https://bitbucket.org/fathomteam/moab/commits/173e849acfb2/
Changeset:   173e849acfb2
Branch:      None
User:        tautges
Date:        2014-03-12 16:56:11
Summary:     - add ability to specify master and slave mesh solid/fluid set numbers separately, so set numbers don't need to agree
- add ability to specify only fluid or solid sets, with other sets derived from that and all material sets in a file
- add DeformMeshRemap to list of examples in makefile

Affected #:  1 file

diff --git a/examples/DeformMeshRemap.cpp b/examples/DeformMeshRemap.cpp
index 109fd3b..0292cf7 100644
--- a/examples/DeformMeshRemap.cpp
+++ b/examples/DeformMeshRemap.cpp
@@ -354,6 +354,50 @@ ErrorCode DeformMeshRemap::execute()
   return MB_SUCCESS;
 }
 
+ErrorCode DeformMeshRemap::find_other_sets(int m_or_s, EntityHandle file_set) 
+{
+    // solid or fluid sets are missing; find the other
+  Range *filled_sets = NULL, *unfilled_sets = NULL, *unfilled_elems = NULL;
+  
+  if (fluidSets[m_or_s].empty() && !solidSets[m_or_s].empty()) {
+    unfilled_sets = &fluidSets[m_or_s];
+    filled_sets = &solidSets[m_or_s];
+    unfilled_elems = &fluidElems[m_or_s];
+  }
+  else if (!fluidSets[m_or_s].empty() && solidSets[m_or_s].empty()) {
+    filled_sets = &fluidSets[m_or_s];
+    unfilled_sets = &solidSets[m_or_s];
+    unfilled_elems = &solidElems[m_or_s];
+  }
+  
+    // ok, we know the filled sets, now fill the unfilled sets, and the elems from those
+  Tag tagh;
+  ErrorCode rval = mbImpl->tag_get_handle(MATERIAL_SET_TAG_NAME, tagh); RR("Couldn't get material set tag name.");
+  Range matsets;
+  rval = mbImpl->get_entities_by_type_and_tag(file_set, MBENTITYSET, &tagh, NULL, 1, matsets);
+  if (matsets.empty()) rval = MB_FAILURE;
+  RR("Couldn't get any material sets.");
+  *unfilled_sets = subtract(matsets, *filled_sets);
+  if (unfilled_sets->empty()) {
+    rval = MB_FAILURE;
+    RR("Failed to find any unfilled material sets.");
+  }
+  Range tmp_range;
+  for (Range::iterator rit = unfilled_sets->begin(); rit != unfilled_sets->end(); rit++) {
+    rval = mbImpl->get_entities_by_handle(*rit, tmp_range, true);
+    RR("Failed to get entities in unfilled set.");
+  }
+  int dim = mbImpl->dimension_from_handle(*tmp_range.rbegin());
+  assert(dim > 0 && dim < 4);  
+  *unfilled_elems = tmp_range.subset_by_dimension(dim);
+  if (unfilled_elems->empty()) {
+    rval = MB_FAILURE;
+    RR("Failed to find any unfilled set entities.");
+  }
+  
+  return MB_SUCCESS;
+}
+
 std::string DeformMeshRemap::get_file_name(int m_or_s) const
 {
   switch (m_or_s) {


https://bitbucket.org/fathomteam/moab/commits/cf1c8106890e/
Changeset:   cf1c8106890e
Branch:      None
User:        tautges
Date:        2014-03-12 16:57:08
Summary:     Better error handling & debug printing.
In Tree, handle legacy tag on tree root by deleting and re-creating (legacy files have old-sized tag).

Affected #:  1 file

diff --git a/examples/DeformMeshRemap.cpp b/examples/DeformMeshRemap.cpp
index 0292cf7..109fd3b 100644
--- a/examples/DeformMeshRemap.cpp
+++ b/examples/DeformMeshRemap.cpp
@@ -354,50 +354,6 @@ ErrorCode DeformMeshRemap::execute()
   return MB_SUCCESS;
 }
 
-ErrorCode DeformMeshRemap::find_other_sets(int m_or_s, EntityHandle file_set) 
-{
-    // solid or fluid sets are missing; find the other
-  Range *filled_sets = NULL, *unfilled_sets = NULL, *unfilled_elems = NULL;
-  
-  if (fluidSets[m_or_s].empty() && !solidSets[m_or_s].empty()) {
-    unfilled_sets = &fluidSets[m_or_s];
-    filled_sets = &solidSets[m_or_s];
-    unfilled_elems = &fluidElems[m_or_s];
-  }
-  else if (!fluidSets[m_or_s].empty() && solidSets[m_or_s].empty()) {
-    filled_sets = &fluidSets[m_or_s];
-    unfilled_sets = &solidSets[m_or_s];
-    unfilled_elems = &solidElems[m_or_s];
-  }
-  
-    // ok, we know the filled sets, now fill the unfilled sets, and the elems from those
-  Tag tagh;
-  ErrorCode rval = mbImpl->tag_get_handle(MATERIAL_SET_TAG_NAME, tagh); RR("Couldn't get material set tag name.");
-  Range matsets;
-  rval = mbImpl->get_entities_by_type_and_tag(file_set, MBENTITYSET, &tagh, NULL, 1, matsets);
-  if (matsets.empty()) rval = MB_FAILURE;
-  RR("Couldn't get any material sets.");
-  *unfilled_sets = subtract(matsets, *filled_sets);
-  if (unfilled_sets->empty()) {
-    rval = MB_FAILURE;
-    RR("Failed to find any unfilled material sets.");
-  }
-  Range tmp_range;
-  for (Range::iterator rit = unfilled_sets->begin(); rit != unfilled_sets->end(); rit++) {
-    rval = mbImpl->get_entities_by_handle(*rit, tmp_range, true);
-    RR("Failed to get entities in unfilled set.");
-  }
-  int dim = mbImpl->dimension_from_handle(*tmp_range.rbegin());
-  assert(dim > 0 && dim < 4);  
-  *unfilled_elems = tmp_range.subset_by_dimension(dim);
-  if (unfilled_elems->empty()) {
-    rval = MB_FAILURE;
-    RR("Failed to find any unfilled set entities.");
-  }
-  
-  return MB_SUCCESS;
-}
-
 std::string DeformMeshRemap::get_file_name(int m_or_s) const
 {
   switch (m_or_s) {


https://bitbucket.org/fathomteam/moab/commits/96b0903bde31/
Changeset:   96b0903bde31
Branch:      None
User:        tautges
Date:        2014-03-12 16:57:32
Summary:     DeformMeshRemap: more debugging output, and handle coord displacements properly (add to prev coords)
LloydSmoother: fix a bug with fixed tag entity handles size
SpatialLocator: if an element range is specified, build the tree too

Affected #:  1 file

diff --git a/examples/DeformMeshRemap.cpp b/examples/DeformMeshRemap.cpp
index 109fd3b..f752b4e 100644
--- a/examples/DeformMeshRemap.cpp
+++ b/examples/DeformMeshRemap.cpp
@@ -532,6 +532,7 @@ ErrorCode DeformMeshRemap::deform_master(Range &fluid_elems, Range &solid_elems,
     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.");
     std::vector<double> disps(num_verts);
+    
     for (int i = 0; i < 3; i++) {
       rval = mbImpl->tag_get_handle(xDispNames[0].c_str(), 1, MB_TYPE_DOUBLE, xDisp[i]);
       RR("Failed to get xDisp tag.");


https://bitbucket.org/fathomteam/moab/commits/661a70a4dda5/
Changeset:   661a70a4dda5
Branch:      None
User:        tautges
Date:        2014-03-12 16:57:32
Summary:     DeformMeshRemap: add some tree statistics output
AdaptiveKDTree: make several functions non-static so they show up in profiling; also use BoundBox a bit more
BVHTree: account for BoundBox change
SpatialLocator: choose tree type based on contents of myElems, and move tree construction to a separate function
BoundBox: change name of contains_box to intersects_box, and fix logic that was broken there
TreeStats.hpp: improve statistics output
test_boundbox: fix tests for changes to BoundBox

Affected #:  1 file

diff --git a/src/moab/SpatialLocator.hpp b/src/moab/SpatialLocator.hpp
index 3d29c55..2690058 100644
--- a/src/moab/SpatialLocator.hpp
+++ b/src/moab/SpatialLocator.hpp
@@ -193,6 +193,12 @@ namespace moab {
       ErrorCode register_src_with_intermediate_procs(ParallelComm *pc, double abs_iter_tol, TupleList &TLreg_o);
       
 #endif
+
+        /** Create a tree
+         * Tree type depends on what's in myElems: if empty or all vertices, creates a kdtree,
+         * otherwise creates a BVHTree.
+         */
+      void create_tree();
       
         /* MOAB instance */
       Interface* mbImpl;


https://bitbucket.org/fathomteam/moab/commits/e0f9f1e46ca6/
Changeset:   e0f9f1e46ca6
Branch:      None
User:        tautges
Date:        2014-03-12 16:57:32
Summary:     DeformMeshRemap.cpp: Write more intermediate files if compiled with local debug flag true; change some of the functions for moving coords to tags and vice versa to facilitate this.
SpatialLocator.cpp: Back out change that automatically chose BVH tree for non-vertex entities, there's a bug in point location there currently...

Affected #:  2 files

diff --git a/examples/DeformMeshRemap.cpp b/examples/DeformMeshRemap.cpp
index f752b4e..a8a0e10 100644
--- a/examples/DeformMeshRemap.cpp
+++ b/examples/DeformMeshRemap.cpp
@@ -110,10 +110,13 @@ private:
   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);
+    //! if non-zero tmp_tag is input, save coords to tmp_tag before over-writing with tag value
+  ErrorCode write_to_coords(Range &elems, Tag tagh, Tag tmp_tag = 0);
 
     //! 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);
+    //! if restore_coords is true, coords are restored to their initial state after file is written
+  ErrorCode write_and_save(Range &ents, EntityHandle seth, Tag tagh, const char *filename,
+                           bool restore_coords = false);
 
     //! find fluid/solid sets from complement of solid/fluid sets 
   ErrorCode find_other_sets(int m_or_s, EntityHandle file_set);
@@ -262,7 +265,7 @@ ErrorCode DeformMeshRemap::execute()
   src_elems.merge(fluidElems[MASTER]);
 
     // initialize data coupler on source elements
-  DataCoupler dc_master(mbImpl, NULL, src_elems, 0);
+  DataCoupler dc_master(mbImpl, src_elems, 0, NULL);
 
   Range tgt_verts;
   if (have_slave) {
@@ -316,6 +319,7 @@ ErrorCode DeformMeshRemap::execute()
 
     // transfer xNew to coords, for master
   if (debug) std::cout << "Transferring coords tag to vertex coordinates in master mesh..." << std::endl;
+  rval = write_to_coords(solidElems[MASTER], xNew); RR("Failed writing tag to master fluid verts.");
   rval = write_to_coords(fluidElems[MASTER], xNew); RR("Failed writing tag to master fluid verts.");
 
   if (have_slave) {
@@ -334,8 +338,8 @@ ErrorCode DeformMeshRemap::execute()
     if (pcMaster && pcMaster->size() > 1) 
       str = "PARALLEL=WRITE_PART";
 #endif
-    if (debug) std::cout << "Writing smoothed_master.vtk..." << std::endl;
-    rval = mbImpl->write_file("smoothed_master.vtk", NULL, str.c_str(), &masterSet, 1);
+    if (debug) std::cout << "Writing smoothed_master.h5m..." << std::endl;
+    rval = mbImpl->write_file("smoothed_master.h5m", NULL, str.c_str(), &masterSet, 1);
 
     if (have_slave) {
 #ifdef USE_MPI
@@ -343,8 +347,8 @@ ErrorCode DeformMeshRemap::execute()
       if (pcSlave && pcSlave->size() > 1) 
         str = "PARALLEL=WRITE_PART";
 #endif
-      if (debug) std::cout << "Writing slave_interp.vtk..." << std::endl;
-      rval = mbImpl->write_file("slave_interp.vtk", NULL, str.c_str(), &slaveSet, 1);
+      if (debug) std::cout << "Writing slave_interp.h5m..." << std::endl;
+      rval = mbImpl->write_file("slave_interp.h5m", NULL, str.c_str(), &slaveSet, 1);
     } // if have_slave
   } // if debug
 
@@ -470,20 +474,38 @@ int main(int argc, char **argv) {
   return rval;
 }
 
-ErrorCode DeformMeshRemap::write_and_save(Range &ents, EntityHandle seth, Tag tagh, const char *filename) 
+ErrorCode DeformMeshRemap::write_and_save(Range &ents, EntityHandle seth, Tag tagh, const char *filename,
+                                          bool restore_coords) 
 {
-  ErrorCode rval = write_to_coords(ents, tagh); RR("");
+  Tag tmp_tag = 0;
+  ErrorCode rval;
+  if (restore_coords)
+    rval = mbImpl->tag_get_handle("", 3, MB_TYPE_DOUBLE, tmp_tag, MB_TAG_CREAT|MB_TAG_DENSE);
+
+  rval = write_to_coords(ents, tagh, tmp_tag); RR("");
   rval = mbImpl->write_file(filename, NULL, NULL, &seth, 1); RR("");
+  if (restore_coords) {
+    rval = write_to_coords(ents, tmp_tag); RR("");
+    rval = mbImpl->tag_delete(tmp_tag);
+  }
+  
   return rval;
 }
   
-ErrorCode DeformMeshRemap::write_to_coords(Range &elems, Tag tagh) 
+ErrorCode DeformMeshRemap::write_to_coords(Range &elems, Tag tagh, Tag tmp_tag) 
 {
     // write the tag to coordinates
   Range verts;
   ErrorCode rval = mbImpl->get_adjacencies(elems, 0, false, verts, Interface::UNION);
   RR("Failed to get adj vertices.");
   std::vector<double> coords(3*verts.size());
+
+  if (tmp_tag) {
+      // save the coords to tmp_tag first
+    rval = mbImpl->get_coords(verts, &coords[0]); RR("Failed to get tmp copy of coords.");
+    rval = mbImpl->tag_set_data(tmp_tag, verts, &coords[0]); RR("Failed to save tmp copy of coords.");
+  }
+  
   rval = mbImpl->tag_get_data(tagh, verts, &coords[0]);
   RR("Failed to get tag data.");
   rval = mbImpl->set_coords(verts, &coords[0]);
@@ -503,7 +525,6 @@ void deform_func(const BoundBox &bbox, double *xold, double *xnew)
   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);
@@ -517,13 +538,13 @@ ErrorCode DeformMeshRemap::deform_master(Range &fluid_elems, Range &solid_elems,
   ErrorCode rval;
 
     // get all the vertices and coords in the solid
-  Range verts;
-  rval = mbImpl->get_adjacencies(solid_elems, 0, false, verts, Interface::UNION);
+  Range solid_verts, fluid_verts;
+  rval = mbImpl->get_adjacencies(solid_elems, 0, false, solid_verts, Interface::UNION);
   RR("Failed to get vertices.");
-  std::vector<double> coords(3*verts.size()), new_coords(3*verts.size());
-  rval = mbImpl->get_coords(verts, &coords[0]);
+  std::vector<double> coords(3*solid_verts.size()), new_coords(3*solid_verts.size());
+  rval = mbImpl->get_coords(solid_verts, &coords[0]);
   RR("Failed to get vertex coords.");
-  unsigned int num_verts = verts.size();
+  unsigned int num_verts = solid_verts.size();
 
     // get or create the tag
 
@@ -536,7 +557,7 @@ ErrorCode DeformMeshRemap::deform_master(Range &fluid_elems, Range &solid_elems,
     for (int i = 0; i < 3; i++) {
       rval = mbImpl->tag_get_handle(xDispNames[0].c_str(), 1, MB_TYPE_DOUBLE, xDisp[i]);
       RR("Failed to get xDisp tag.");
-      rval = mbImpl->tag_get_data(xDisp[i], verts, &disps[0]);
+      rval = mbImpl->tag_get_data(xDisp[i], solid_verts, &disps[0]);
       RR("Failed to get xDisp tag values.");
       for (unsigned int j = 0; j < num_verts; j++)
         new_coords[3*j+i] = coords[3*j+i] + disps[j];
@@ -547,7 +568,7 @@ ErrorCode DeformMeshRemap::deform_master(Range &fluid_elems, Range &solid_elems,
     RR("Failed to get first xDisp tag.");
     xNew = xDisp[0];
     std::vector<double> disps(3*num_verts);
-    rval = mbImpl->tag_get_data(xDisp[0], verts, &disps[0]);
+    rval = mbImpl->tag_get_data(xDisp[0], solid_verts, &disps[0]);
     for (unsigned int j = 0; j < 3*num_verts; j++)
       new_coords[j] = coords[j] + disps[j];
   }
@@ -584,18 +605,26 @@ ErrorCode DeformMeshRemap::deform_master(Range &fluid_elems, Range &solid_elems,
   }
     
     // set the new tag to those coords
-  rval = mbImpl->tag_set_data(xNew, verts, &coords[0]);
+  rval = mbImpl->tag_set_data(xNew, solid_verts, &new_coords[0]);
   RR("Failed to set tag data.");
   
     // get all the vertices and coords in the fluid, set xnew to them
-  verts.clear();
-  rval = mbImpl->get_adjacencies(fluid_elems, 0, false, verts, Interface::UNION);
+  rval = mbImpl->get_adjacencies(fluid_elems, 0, false, fluid_verts, Interface::UNION);
   RR("Failed to get vertices.");
-  coords.resize(3*verts.size());
-  rval = mbImpl->get_coords(verts, &coords[0]);
+  fluid_verts = subtract(fluid_verts, solid_verts);
+  
+  if (coords.size() < 3*fluid_verts.size()) coords.resize(3*fluid_verts.size());
+  rval = mbImpl->get_coords(fluid_verts, &coords[0]);
   RR("Failed to get vertex coords.");
-  rval = mbImpl->tag_set_data(xNew, verts, &coords[0]);
+  rval = mbImpl->tag_set_data(xNew, fluid_verts, &coords[0]);
   RR("Failed to set xnew tag on fluid verts.");
+ 
+  if (debug) {
+      // save deformed mesh coords to new file for visualizing
+    Range tmp_range(fluidElems[MASTER]);
+    tmp_range.merge(solidElems[MASTER]);
+    rval = write_and_save(tmp_range, masterSet, xNew, "deformed_master.h5m", true); RR("");
+  }
     
   return MB_SUCCESS;
 }

diff --git a/src/SpatialLocator.cpp b/src/SpatialLocator.cpp
index 2903155..63c2fdc 100644
--- a/src/SpatialLocator.cpp
+++ b/src/SpatialLocator.cpp
@@ -35,12 +35,12 @@ namespace moab
     {
       if (myTree) return;
       
-      if (myElems.empty() || mbImpl->type_from_handle(*myElems.rbegin()) == MBVERTEX) 
+//      if (myElems.empty() || mbImpl->type_from_handle(*myElems.rbegin()) == MBVERTEX) 
           // create a kdtree if only vertices
         myTree = new AdaptiveKDTree(mbImpl);
-      else
+//      else
           // otherwise a BVHtree, since it performs better for elements
-        myTree = new BVHTree(mbImpl);
+//        myTree = new BVHTree(mbImpl);
 
       iCreatedTree = true;
     }


https://bitbucket.org/fathomteam/moab/commits/b7b7922b5372/
Changeset:   b7b7922b5372
Branch:      None
User:        tautges
Date:        2014-03-12 16:57:32
Summary:     Various smallish bug fixes and minor changes, improving the use of SpatialLocator.

ElemEvaluator: in set_ent_handle, if there's no evalSet for that entity type, call EvalSet::get_eval_set
SpatialLocator: set the evaluator on the tree if it's not the same as what's on SpatialLocator; also, replace find_point logic
  (that used distance_search for some reason) with call to Tree::point_search.
AdaptiveKDTree: move parse_options to public function, as it is in Tree
BVHTree: fixed bug, where BVHTree had its own myEval shadowing the one declared in Tree
Tree: fixed bug, initialize myEval to NULL in constructor
spatial_locator_test: test trees with and without evaluators
DataCoupler: improving comments a bit, and move ParallelComm down in arg list for constructor; also, some
  improved comments
datacoupler_test: adding test


Passes make check in parallel.

Affected #:  10 files

diff --git a/src/LocalDiscretization/moab/ElemEvaluator.hpp b/src/LocalDiscretization/moab/ElemEvaluator.hpp
index 866388c..6a1efd0 100644
--- a/src/LocalDiscretization/moab/ElemEvaluator.hpp
+++ b/src/LocalDiscretization/moab/ElemEvaluator.hpp
@@ -339,8 +339,13 @@ namespace moab {
       std::vector<EntityHandle> dum_vec;
       ErrorCode rval = mbImpl->get_connectivity(ent, vertHandles, numVerts, false, &dum_vec);
       if (MB_SUCCESS != rval) return rval;
+
+      if (!evalSets[entType].evalFcn)
+        EvalSet::get_eval_set(entType, numVerts, evalSets[entType]);
+
       rval = mbImpl->get_coords(vertHandles, numVerts, vertPos[0].array());
       if (MB_SUCCESS != rval) return rval;
+
       if (tagHandle) {
         rval = set_tag_handle(tagHandle);
         if (MB_SUCCESS != rval) return rval;

diff --git a/src/SpatialLocator.cpp b/src/SpatialLocator.cpp
index 63c2fdc..702399e 100644
--- a/src/SpatialLocator.cpp
+++ b/src/SpatialLocator.cpp
@@ -35,12 +35,12 @@ namespace moab
     {
       if (myTree) return;
       
-//      if (myElems.empty() || mbImpl->type_from_handle(*myElems.rbegin()) == MBVERTEX) 
+      if (myElems.empty() || mbImpl->type_from_handle(*myElems.rbegin()) == MBVERTEX) 
           // create a kdtree if only vertices
         myTree = new AdaptiveKDTree(mbImpl);
-//      else
+      else
           // otherwise a BVHtree, since it performs better for elements
-//        myTree = new BVHTree(mbImpl);
+        myTree = new BVHTree(mbImpl);
 
       iCreatedTree = true;
     }
@@ -367,90 +367,29 @@ namespace moab
           // relative epsilon given, translate to absolute epsilon using box dimensions
         tmp_abs_iter_tol = rel_iter_tol * localBox.diagonal_length();
       }
-  
-      EntityHandle closest_leaf;
-      std::vector<double> dists;
-      std::vector<EntityHandle> leaves;
-      ErrorCode rval = MB_SUCCESS;
 
+      if (elemEval && myTree->get_eval() != elemEval)
+        myTree->set_eval(elemEval);
+      
+      ErrorCode rval = MB_SUCCESS;
       for (int i = 0; i < num_points; i++) {
         int i3 = 3*i;
-        ents[i] = 0;
-        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
-            double min_dist = *dists.begin();
-            closest_leaf = *leaves.begin();
-            std::vector<EntityHandle>::iterator vit = leaves.begin()+1;
-            std::vector<double>::iterator dit = dists.begin()+1;
-            for (; vit != leaves.end() && min_dist; vit++, dit++) {
-              if (*dit < min_dist) {
-                min_dist = *dit;
-                closest_leaf = *vit;
-              }
-            }
-            dists.clear();
-            leaves.clear();
-          }
-        }
-        else {
-          rval = myTree->point_search(pos+i3, closest_leaf);
-          if (MB_ENTITY_NOT_FOUND == rval) closest_leaf = 0;
-          else if (MB_SUCCESS != rval) return rval;
-        }
-
-          // if no ElemEvaluator, just return the box
-        if (!elemEval) {
-          ents[i] = closest_leaf;
-          params[i3] = params[i3+1] = params[i3+2] = -2;
-          if (is_inside && closest_leaf) is_inside[i] = true;
+        ErrorCode tmp_rval = myTree->point_search(pos+i3, ents[i], abs_iter_tol, inside_tol, NULL, NULL, 
+                                                  (CartVect*)(params+i3));
+        if (MB_SUCCESS != tmp_rval) {
+          rval = tmp_rval;
           continue;
         }
-    
-          // find natural coordinates of point in element(s) in that leaf
-        CartVect tmp_nat_coords; 
-        Range range_leaf;
-        rval = mbImpl->get_entities_by_dimension(closest_leaf, myDim, range_leaf, false);
-        if(rval != MB_SUCCESS) return rval;
-
-          // loop over the range_leaf
-        int tmp_inside;
-        int *is_ptr = (is_inside ? is_inside+i : &tmp_inside);      
-        *is_ptr = false;
-        EntityHandle ent = 0;
-        for(Range::iterator rit = range_leaf.begin(); rit != range_leaf.end(); rit++)
-        {
-          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;
-          if (*is_ptr) {
-            ent = *rit;
-            break;
-          }
-        }
-        if (debug && !ent) {
+
+        if (debug && !ents[i]) {
           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;
+        if (is_inside) is_inside[i] = (ents[i] ? true : false);
+      }
+      
+      return rval;
     }
     
         /* Count the number of located points in locTable

diff --git a/src/moab/AdaptiveKDTree.hpp b/src/moab/AdaptiveKDTree.hpp
index 450a418..e4c5078 100644
--- a/src/moab/AdaptiveKDTree.hpp
+++ b/src/moab/AdaptiveKDTree.hpp
@@ -38,6 +38,13 @@ namespace moab {
 
       ~AdaptiveKDTree();
 
+        /** \brief Parse options for tree creation
+         * \param options Options passed in by application
+         * \return Failure is returned if any options were passed in and not interpreted; could mean
+         * inappropriate options for a particular tree type
+         */
+      ErrorCode parse_options(FileOptions &options);
+
         /** Build the tree
          * Build a tree with the entities input.  If a non-NULL tree_root_set pointer is input, 
          * use the pointed-to set as the root of this tree (*tree_root_set!=0) otherwise construct 
@@ -253,13 +260,6 @@ namespace moab {
   private:
       friend class AdaptiveKDTreeIter;
 
-        /** \brief Parse options for tree creation
-         * \param options Options passed in by application
-         * \return Failure is returned if any options were passed in and not interpreted; could mean
-         * inappropriate options for a particular tree type
-         */
-      ErrorCode parse_options(FileOptions &options);
-
       ErrorCode init();
   
         /**\brief find a triangle near the input point */

diff --git a/src/moab/BVHTree.hpp b/src/moab/BVHTree.hpp
index d435ebf..354d3a1 100644
--- a/src/moab/BVHTree.hpp
+++ b/src/moab/BVHTree.hpp
@@ -272,7 +272,6 @@ namespace moab {
       
       Range entityHandles;
       std::vector<TreeNode> myTree;
-      ElemEvaluator *myEval;
       int splitsPerDir;
       EntityHandle startSetHandle;
       static const char *treeName;
@@ -305,7 +304,7 @@ namespace moab {
     }
 
     inline BVHTree::BVHTree(Interface *impl) : 
-            Tree(impl), myEval(NULL), splitsPerDir(3), startSetHandle(0) {boxTagName = treeName;}
+            Tree(impl), splitsPerDir(3), startSetHandle(0) {boxTagName = treeName;}
 
     inline unsigned int BVHTree::set_interval(BoundBox &interval, 
                                               std::vector<Bucket>::const_iterator begin, 

diff --git a/src/moab/Tree.hpp b/src/moab/Tree.hpp
index 67b67f9..0a995b4 100644
--- a/src/moab/Tree.hpp
+++ b/src/moab/Tree.hpp
@@ -236,7 +236,7 @@ namespace moab {
 
     inline Tree::Tree(Interface* iface) 
             : mbImpl(iface), maxPerLeaf(6), maxDepth(30), treeDepth(-1), minWidth(1.0e-10),
-              meshsetFlags(0), cleanUp(true), myRoot(0), boxTag(0)
+              meshsetFlags(0), cleanUp(true), myRoot(0), boxTag(0), myEval(0)
     {}
 
     inline Tree::~Tree() 

diff --git a/test/spatial_locator_test.cpp b/test/spatial_locator_test.cpp
index 0171bf9..404e668 100644
--- a/test/spatial_locator_test.cpp
+++ b/test/spatial_locator_test.cpp
@@ -4,7 +4,9 @@
 #include "moab/HomXform.hpp"
 #include "moab/ScdInterface.hpp"
 #include "moab/CartVect.hpp"
+#include "moab/AdaptiveKDTree.hpp"
 #include "moab/BVHTree.hpp"
+#include "moab/ElemEvaluator.hpp"
 #include "moab/ProgOptions.hpp"
 #include "moab/CpuTimer.hpp"
 
@@ -69,11 +71,21 @@ void test_kd_tree()
   Range elems;
   rval = create_hex_mesh(mb, elems, ints, 3); CHECK_ERR(rval);
 
-    // initialize spatial locator with the elements and the default tree type
-  SpatialLocator *sl = new SpatialLocator(&mb, elems);
+    // initialize spatial locator with the elements and KDtree
+  AdaptiveKDTree kd(&mb);
+  std::ostringstream opts;
+  opts << "MAX_DEPTH=" << max_depth << ";MAX_PER_LEAF=" << leaf;
+  FileOptions fo(opts.str().c_str());
+  rval = kd.parse_options(fo);
+  SpatialLocator *sl = new SpatialLocator(&mb, elems, &kd);
 
   test_locator(sl);
 
+    // test with an evaluator
+  ElemEvaluator eval(&mb);
+  kd.set_eval(&eval);
+  test_locator(sl);
+
     // destroy spatial locator, and tree along with it
   delete sl;
 }
@@ -96,6 +108,11 @@ void test_bvh_tree()
   SpatialLocator *sl = new SpatialLocator(&mb, elems, &bvh);
   test_locator(sl);
   
+    // test with an evaluator
+  ElemEvaluator eval(&mb);
+  bvh.set_eval(&eval);
+  test_locator(sl);
+
     // destroy spatial locator, and tree along with it
   delete sl;
 }

diff --git a/tools/mbcoupler/DataCoupler.cpp b/tools/mbcoupler/DataCoupler.cpp
index a5a3ec1..06a8564 100644
--- a/tools/mbcoupler/DataCoupler.cpp
+++ b/tools/mbcoupler/DataCoupler.cpp
@@ -1,24 +1,27 @@
 #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"
 
+#ifdef USE_MPI
+#include "moab/ParallelComm.hpp"
+#endif
+
 #include "iostream"
-#include <stdio.h>
 #include <algorithm>
 #include <sstream>
 
+#include <stdio.h>
 #include "assert.h"
 
 namespace moab {
 
 DataCoupler::DataCoupler(Interface *impl,
-                         ParallelComm *pc,
                          Range &source_ents,
                          int coupler_id,
+                         ParallelComm *pc,
                          bool init_locator,
                          int dim)
         : mbImpl(impl), myPcomm(pc), myId(coupler_id), myDim(dim)
@@ -64,7 +67,12 @@ ErrorCode DataCoupler::locate_points(Range &targ_ents,
                                      const double inside_tol)
 {
   targetEnts = targ_ents;
-  
+
+#ifdef USE_MPI
+  if (myPcomm && myPcomm->size() > 1) 
+    return myLocator->par_locate_points(myPcomm, targ_ents, rel_iter_tol, abs_iter_tol, inside_tol);
+#endif
+
   return myLocator->locate_points(targ_ents, rel_iter_tol, abs_iter_tol, inside_tol);
 }
 
@@ -72,6 +80,11 @@ ErrorCode DataCoupler::locate_points(double *xyz, int num_points,
                                      const double rel_iter_tol, const double abs_iter_tol,
                                      const double inside_tol)
 {
+#ifdef USE_MPI
+  if (myPcomm && myPcomm->size() > 1) 
+    return myLocator->par_locate_points(myPcomm, xyz, num_points, rel_iter_tol, abs_iter_tol, inside_tol);
+#endif
+
   return myLocator->locate_points(xyz, num_points, rel_iter_tol, abs_iter_tol, inside_tol);
 }
 
@@ -95,7 +108,7 @@ ErrorCode DataCoupler::interpolate(/*DataCoupler::Method*/ int method,
   return interpolate(method, tag, interp_vals, point_indices, normalize);
 }
   
-ErrorCode DataCoupler::interpolate(/*DataCoupler::Method*/ int *methods,
+ErrorCode DataCoupler::interpolate(/*DataCoupler::Method*/ int */* methods */,
                                    Tag *tags,
                                    int *points_per_method,
                                    int num_methods,
@@ -111,7 +124,12 @@ ErrorCode DataCoupler::interpolate(/*DataCoupler::Method*/ int *methods,
   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());
+    // # points and indices should be identical
+  if (pts_total != num_indices)
+    return MB_FAILURE;
 
+    // since each tuple contains one interpolated tag, if we're interpolating multiple tags, every tuple
+    // needs to be able to store up to max tag size
   int max_tsize = -1;
   for (int i = 0; i < num_methods; i++) {
     int tmp_tsize;
@@ -120,75 +138,82 @@ ErrorCode DataCoupler::interpolate(/*DataCoupler::Method*/ int *methods,
     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());
+  if (myPcomm && myPcomm->size() > 1) {
+      // build up the tuple list to distribute from my target points; assume that
+      // all procs use same method/tag input
+    TupleList TLob; // TLob structure: (pto_i, ridx_i, lidx_i, meth_tagidx_i)
+    TLob.initialize(4, 0, 0, 0, num_indices);
+    int tn = 0; // the tuple number we're currently on
+    TLob.enableWriteAccess();
+    for (int m = 0; m < num_methods; m++) {
+      int num_points = (points_per_method ? points_per_method[m] : targetEnts.size());
       for (int j = 0; j < num_points; j++) {
-        int idx = (point_indices ? (*point_indices)[j] : j);
+        int idx = (point_indices ? (*point_indices)[j] : j); // the index in my targetEnts for this interpolation point
 
           // 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
+        TLob.vi_wr[4*tn]   = myLocator->par_loc_table().vi_rd[2*idx]; // proc
+        TLob.vi_wr[4*tn+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++;
+        TLob.vi_wr[4*tn+2] = idx;
+        TLob.vi_wr[4*tn+3] = m;
+        TLob.inc_n();
+        tn++;
       }
     }
 
       // 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]];
-
+    myPcomm->proc_config().crystal_router()->gs_transfer(1, TLob, 0);
+
+      // perform interpolation on local source mesh; put results into TLinterp
+    TupleList TLinterp; // TLinterp structure: (pto_i, ridx_i, vals[max_tsize]_d)
+    TLinterp.initialize(2, 0, 0, max_tsize, TLob.get_n());
+    TLinterp.set_n(TLob.get_n()); // set the size right away
+    TLinterp.enableWriteAccess();
+
+    for (unsigned int i = 0; i < TLob.get_n(); i++) {
+      int lidx = TLob.vi_rd[4*i+1]; // index into myLocator's local table
+        // method and tag indexed with same index
+//      /*Method*/ int method = methods[TLob.vi_rd[4*i+3]];
+      Tag tag = tags[TLob.vi_rd[4*i+3]];
+        // copy proc/remote index from TLob
+      TLinterp.vi_wr[2*i] = TLob.vi_rd[4*i];
+      TLinterp.vi_wr[2*i+1] = TLob.vi_rd[4*i+2]; 
+
+        // set up the evaluator for the tag and entity, then interpolate, putting result in TLinterp
       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);
+      result = myLocator->elem_eval()->eval(myLocator->loc_table().vr_rd+3*lidx, TLinterp.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);
+    myPcomm->proc_config().crystal_router()->gs_transfer(1, TLinterp, 0);
 
       // copy the interpolated field as a unit
-    std::copy(tinterp.vr_rd, tinterp.vr_rd+tinterp.get_n()*max_tsize, interp_vals);
+    std::copy(TLinterp.vr_rd, TLinterp.vr_rd+TLinterp.get_n()*max_tsize, interp_vals);
   }
   else {
+      // if called serially, interpolate directly from local locations/entities,
+      // into either interp_vals or by setting data directly on tags on those entities
     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)
+        // interpolated data is tsize long, which is either max size (if data passed back to caller in interp_vals)
         // or tag size (if data will be set on entities, in which case it shouldn't have padding)
+        // set sizes here, inside loop over methods, to reduce memory usage
       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);
+        if (MB_SUCCESS != result) return result;
         result = mbImpl->tag_get_bytes(tags[i], tsize_bytes);
+        if (MB_SUCCESS != result) return result;
       }
       
       for (int j = 0; j < num_points; j++) {
@@ -206,7 +231,7 @@ ErrorCode DataCoupler::interpolate(/*DataCoupler::Method*/ int *methods,
         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
+      } // for j over points
 
       if (!interp_vals) {
           // set tags on tmp_ents; data is already w/o padding, due to tsize setting above
@@ -214,7 +239,7 @@ ErrorCode DataCoupler::interpolate(/*DataCoupler::Method*/ int *methods,
         if (MB_SUCCESS != result) return result;
       }
 
-    } // for i
+    } // for i over methods
   } // if myPcomm
   
       // done

diff --git a/tools/mbcoupler/DataCoupler.hpp b/tools/mbcoupler/DataCoupler.hpp
index bfe2ba5..a9d6f5e 100644
--- a/tools/mbcoupler/DataCoupler.hpp
+++ b/tools/mbcoupler/DataCoupler.hpp
@@ -39,21 +39,23 @@ class DataCoupler
 {
 public:
 
+  enum Method {CONSTANT, LINEAR_FE, QUADRATIC_FE, SPECTRAL} ;
+
   enum IntegType {VOLUME};
 
     /* constructor
      * Constructor, which also optionally initializes the coupler
+     * \param source_ents Elements in the source mesh
+     * \param coupler_id Id of this coupler, should be the same over all procs
      * \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,
+              ParallelComm *pc = NULL,
               bool init_locator = true,
               int dim = -1);
 

diff --git a/tools/mbcoupler/Makefile.am b/tools/mbcoupler/Makefile.am
index ea32ff9..8114a52 100644
--- a/tools/mbcoupler/Makefile.am
+++ b/tools/mbcoupler/Makefile.am
@@ -40,7 +40,8 @@ libmbcoupler_la_includedir = $(includedir)
 # check that only libraries are going in $(libdir)
 cfgdir = $(libdir)
 
-TESTS = elem_util_test element_test
+TESTS = elem_util_test element_test datacoupler_test
+datacoupler_test_SOURCES = datacoupler_test.cpp
 check_PROGRAMS = $(TESTS) 
 elem_util_test_SOURCES = ElemUtilTest.cpp
 element_test_SOURCES = ElementTest.cpp

diff --git a/tools/mbcoupler/datacoupler_test.cpp b/tools/mbcoupler/datacoupler_test.cpp
new file mode 100644
index 0000000..b39ca0e
--- /dev/null
+++ b/tools/mbcoupler/datacoupler_test.cpp
@@ -0,0 +1,682 @@
+#include "moab/Core.hpp"
+#include "moab/CpuTimer.hpp"
+#include "DataCoupler.hpp"
+#include "ElemUtil.hpp"
+#include "iMesh.h"
+#include "MBiMesh.hpp"
+
+#include <iostream>
+#include <iomanip>
+#include <sstream>
+#include <assert.h>
+
+#ifdef USE_MPI
+#include "moab/ParallelComm.hpp"
+#include "MBParallelConventions.h"
+#endif
+
+using namespace moab;
+
+#define STRINGIFY_(A) #A
+#define STRINGIFY(A) STRINGIFY_(A)
+#ifdef MESHDIR
+std::string TestDir( STRINGIFY(MESHDIR) );
+#else
+std::string TestDir(".");
+#endif
+
+#define RRA(a) if (MB_SUCCESS != result) {\
+      std::string tmp_str; mbImpl->get_last_error(tmp_str);\
+      tmp_str.append("\n"); tmp_str.append(a);\
+      dynamic_cast<Core*>(mbImpl)->get_error_handler()->set_last_error(tmp_str.c_str()); \
+      return result;}
+
+#define PRINT_LAST_ERROR \
+    if (MB_SUCCESS != result) {\
+      std::string tmp_str;\
+      std::cout << "Failure; message:" << std::endl;\
+      mbImpl->get_last_error(tmp_str);\
+      std::cout << tmp_str << std::endl;\
+      MPI_Abort(MPI_COMM_WORLD, result);        \
+      return result;\
+    }
+
+#define PRINT_LAST_ERR \
+    if (iBase_SUCCESS != err) {\
+      std::string tmp_str;\
+      std::cout << "Failure; message:" << std::endl;\
+      mbImpl->get_last_error(tmp_str);\
+      std::cout << tmp_str << std::endl;\
+      MPI_Abort(MPI_COMM_WORLD, result);        \
+      return result;\
+    }
+
+void print_usage(char **argv);
+
+ErrorCode get_file_options(int argc, char **argv, 
+                           std::vector<std::string> &meshFiles,
+                           DataCoupler::Method &method,
+                           std::string &interpTag,
+                           std::string &gNormTag,
+                           std::string &ssNormTag,
+                           std::vector<const char *> &ssTagNames,
+                           std::vector<const char *> &ssTagValues,
+                           std::string &readOpts,
+                           std::string &outFile,
+                           std::string &writeOpts,
+                           std::string &dbgFile,
+                           bool &help,
+                           double & epsilon);
+
+// ErrorCode get_file_options(int argc, char **argv, 
+//                              std::vector<const char *> &filenames,
+//                              std::string &tag_name,
+//                              std::string &out_fname,
+//                              std::string &opts);
+
+#ifdef USE_MPI
+ErrorCode report_iface_ents(Interface *mbImpl,
+                              std::vector<ParallelComm *> &pcs,
+                              bool print_results);
+#endif
+
+ErrorCode test_interpolation(Interface *mbImpl, 
+                             DataCoupler::Method method,
+                             std::string &interpTag,
+                             std::string &gNormTag,
+                             std::string &ssNormTag,
+                             std::vector<const char *> &ssTagNames,
+                             std::vector<const char *> &ssTagValues,
+                             iBase_EntitySetHandle *roots,
+                             std::vector<ParallelComm *> &pcs,
+                             double &instant_time,
+                             double &pointloc_time,
+                             double &interp_time,
+                             double &gnorm_time,
+                             double &ssnorm_time,
+                             double & toler);
+
+void reduceMax(double &v){
+  double buf;
+
+  MPI_Barrier(MPI_COMM_WORLD);
+  MPI_Allreduce( &v, &buf, 1, MPI_DOUBLE, MPI_MAX, MPI_COMM_WORLD);
+
+  v=buf;
+}
+
+
+int main(int argc, char **argv) 
+{
+
+  int err;
+  
+#ifdef USE_MPI
+    // need to init MPI first, to tell how many procs and rank
+  err = MPI_Init(&argc, &argv);
+#endif
+
+  std::vector<const char *> ssTagNames, ssTagValues;
+  std::vector<std::string> meshFiles;
+  std::string interpTag, gNormTag, ssNormTag, readOpts, outFile, writeOpts, dbgFile;
+  DataCoupler::Method method = DataCoupler::CONSTANT;
+
+  ErrorCode result = MB_SUCCESS;
+  bool help = false;
+  double toler = 5.e-10;
+  result = get_file_options(argc, argv, meshFiles, method, interpTag,
+                            gNormTag, ssNormTag, ssTagNames, ssTagValues,
+                            readOpts, outFile, writeOpts, dbgFile, help, toler);
+
+  if (result != MB_SUCCESS || help) {
+    print_usage(argv);
+#ifdef USE_MPI
+    err = MPI_Finalize();
+#endif
+    return 1;
+  }
+
+  int nprocs = 1, rank = 0;
+
+#ifdef USE_MPI  
+  err = MPI_Comm_size(MPI_COMM_WORLD, &nprocs);
+  err = MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+  std::vector<ParallelComm *> pcs(meshFiles.size()); 
+#endif
+
+    // redirect stdout and stderr if dbgFile is not null
+  if (!dbgFile.empty()) {
+    std::stringstream dfname;
+    dfname << dbgFile << rank << ".txt";
+    std::freopen(dfname.str().c_str(), "a", stdout);
+    std::freopen(dfname.str().c_str(), "a", stderr);
+  }
+
+  // create MOAB instance based on that
+  Interface *mbImpl = new Core();
+  if (NULL == mbImpl) return 1;
+
+  iMesh_Instance iMeshInst = reinterpret_cast<iMesh_Instance>( new MBiMesh(mbImpl) );
+  
+    // read in mesh(es)
+
+    // Create root sets for each mesh using the iMesh API.  Then pass these
+    // to the load_file functions to be populated.
+  iBase_EntitySetHandle *roots = (iBase_EntitySetHandle *) malloc(sizeof(iBase_EntitySetHandle) * meshFiles.size());
+
+  for (unsigned int i = 0; i < meshFiles.size(); i++) {
+    std::string newReadopts;
+    std::ostringstream extraOpt;
+#ifdef USE_MPI
+    pcs[i] = new ParallelComm(mbImpl, MPI_COMM_WORLD);
+    int index = pcs[i]->get_id();
+    extraOpt  << ";PARALLEL_COMM=" << index;
+    newReadopts = readOpts+extraOpt.str();
+#endif
+    
+    iMesh_createEntSet(iMeshInst, 0, &(roots[i]), &err);
+    result = mbImpl->load_file( meshFiles[i].c_str(), (EntityHandle *)&roots[i], newReadopts.c_str() );
+    PRINT_LAST_ERROR;
+  }
+
+#ifdef USE_MPI
+  result = report_iface_ents(mbImpl, pcs, true);
+  PRINT_LAST_ERROR;
+#endif
+
+  double instant_time=0.0, pointloc_time=0.0, interp_time=0.0, gnorm_time=0.0, ssnorm_time=0.0;
+    // test interpolation and global normalization and subset normalization
+
+  result = test_interpolation(mbImpl, method, interpTag, gNormTag, ssNormTag, 
+                              ssTagNames, ssTagValues, roots, pcs, 
+                              instant_time, pointloc_time, interp_time, 
+                              gnorm_time, ssnorm_time, toler);
+  PRINT_LAST_ERROR;
+
+  
+  reduceMax(instant_time);  
+  reduceMax(pointloc_time);
+  reduceMax(interp_time);
+
+  if(0==rank)
+    printf("\nMax time : %g %g %g (inst loc interp -- %d procs )\n", instant_time, 
+	   pointloc_time, interp_time, nprocs);
+
+    // output mesh
+  if (!outFile.empty()) {
+    Range partSets;
+      // only save the target mesh
+    partSets.insert((EntityHandle)roots[1]);
+    std::string newwriteOpts = writeOpts;
+    std::ostringstream extraOpt;
+#ifdef USE_MPI
+    extraOpt  << ";PARALLEL_COMM=" << 1;
+    newwriteOpts += extraOpt.str();
+#endif
+    result = mbImpl->write_file(outFile.c_str(), NULL, newwriteOpts.c_str(), partSets);
+    PRINT_LAST_ERROR;
+    std::cout << "Wrote " << outFile << std::endl;
+    std::cout << "mbcoupler_test complete." << std::endl;
+  }
+
+#ifdef USE_MPI
+  for (unsigned int i = 0; i < meshFiles.size(); i++) {
+    delete pcs[i];
+  }
+#endif
+
+  delete mbImpl;
+  //may be leaking iMeshInst, don't care since it's end of program. Remove above deletes?
+
+#ifdef USE_MPI  
+  err = MPI_Finalize();
+#endif
+
+  return 0;
+}
+
+#ifdef USE_MPI
+ErrorCode report_iface_ents(Interface *mbImpl,
+                              std::vector<ParallelComm *> &pcs,
+                              const bool print_results) 
+{
+  Range iface_ents[6];
+  ErrorCode result = MB_SUCCESS, tmp_result;
+  
+    // now figure out which vertices are shared
+  for (unsigned int p = 0; p < pcs.size(); p++) {
+    for (int i = 0; i < 4; i++) {
+      tmp_result = pcs[p]->get_iface_entities(-1, i, iface_ents[i]);
+      
+      if (MB_SUCCESS != tmp_result) {
+        std::cerr << "get_iface_entities returned error on proc " 
+                  << pcs[p]->proc_config().proc_rank() << "; message: " << std::endl;
+        std::string last_error;
+        result = mbImpl->get_last_error(last_error);
+        if (last_error.empty()) std::cerr << "(none)" << std::endl;
+        else std::cerr << last_error << std::endl;
+        result = tmp_result;
+      }
+      if (0 != i) iface_ents[4].merge(iface_ents[i]);
+    }
+  }
+
+    // report # iface entities
+  result = mbImpl->get_adjacencies(iface_ents[4], 0, false, iface_ents[5], 
+                                   Interface::UNION);
+
+  int rank;
+  MPI_Comm_rank(MPI_COMM_WORLD, &rank);
+
+  if (print_results || iface_ents[0].size() != iface_ents[5].size()) {
+    std::cerr << "Proc " << rank << " iface entities: " << std::endl;
+    for (int i = 0; i < 4; i++)
+      std::cerr << "    " << iface_ents[i].size() << " "
+                << i << "d iface entities." << std::endl;
+    std::cerr << "    (" << iface_ents[5].size() 
+              << " verts adj to other iface ents)" << std::endl;
+  }
+  
+  return result;
+}
+#endif
+
+// Print usage
+void print_usage(char **argv) {
+  std::cerr << "Usage: ";
+  std::cerr << argv[0] << " -meshes <source_mesh><target_mesh> -itag <interp_tag> [-gnorm <gnorm_tag>] [-ssnorm <ssnorm_tag><ssnorm_selection>] [-ropts <roptions>] [-outfile <out_file> [-wopts <woptions>]] [-dbgout [<dbg_file>]]" << std::endl;
+  std::cerr << "    -meshes" << std::endl;
+  std::cerr << "        Read in mesh files <source_mesh> and <target_mesh>." << std::endl;
+  std::cerr << "    -itag" << std::endl;
+  std::cerr << "        Interpolate tag <interp_tag> from source mesh to target mesh." << std::endl;
+  std::cerr << "    -gnorm" << std::endl;
+  std::cerr << "        Normalize the value of tag <gnorm_tag> over then entire mesh and save to" << std::endl;
+  std::cerr << "        tag \"<gnorm_tag>_normf\" on the mesh set.  Do this for all meshes." << std::endl;
+  std::cerr << "    -ssnorm" << std::endl;
+  std::cerr << "        Normalize the value of tag <ssnorm_tag> over subsets of a mesh and save to" << std::endl;
+  std::cerr << "        tag \"<ssnorm_tag>_normf\" on the Entity Set for each subset.  Subsets are selected" << std::endl;
+  std::cerr << "        using criteria in <ssnorm_selection>.  Do this for all meshes." << std::endl;
+  std::cerr << "    -ropts" << std::endl;
+  std::cerr << "        Read in the mesh files using options in <roptions>." << std::endl;
+  std::cerr << "    -outfile" << std::endl;
+  std::cerr << "        Write out target mesh to <out_file>." << std::endl;
+  std::cerr << "    -wopts" << std::endl;
+  std::cerr << "        Write out mesh files using options in <woptions>." << std::endl;
+  std::cerr << "    -dbgout" << std::endl;
+  std::cerr << "        Write stdout and stderr streams to the file \'<dbg_file>.txt\'." << std::endl;
+  std::cerr << "    -eps" << std::endl;
+  std::cerr << "        epsilon" << std::endl;
+  std::cerr << "    -meth <method> (0=CONSTANT, 1=LINEAR_FE, 2=QUADRATIC_FE, 3=SPECTRAL)" << std::endl;
+}
+
+// Check first character for a '-'.
+// Return true if one is found.  False otherwise.
+bool check_for_flag(const char *str) {
+  if (str[0] == '-')
+    return true;
+  else
+    return false;
+}
+
+// New get_file_options() function with added possibilities for mbcoupler_test.
+ErrorCode get_file_options(int argc, char **argv, 
+                           std::vector<std::string> &meshFiles,
+                           DataCoupler::Method &method,
+                           std::string &interpTag,
+                           std::string &gNormTag,
+                           std::string &ssNormTag,
+                           std::vector<const char *> &ssTagNames,
+                           std::vector<const char *> &ssTagValues,
+                           std::string &readOpts,
+                           std::string &outFile,
+                           std::string &writeOpts,
+                           std::string &dbgFile,
+                           bool &help,
+                           double & epsilon)
+{
+  // Initialize some of the outputs to null values indicating not present
+  // in the argument list.
+  gNormTag = "";
+  ssNormTag = "";
+  readOpts = "PARALLEL=READ_PART;PARTITION=PARALLEL_PARTITION;PARTITION_DISTRIBUTE;PARALLEL_RESOLVE_SHARED_ENTS;PARALLEL_GHOSTS=3.0.1;CPUTIME";
+  outFile = "";
+  writeOpts = "PARALLEL=WRITE_PART;CPUTIME";
+  dbgFile = "";
+  std::string defaultDbgFile = argv[0];  // The executable name will be the default debug output file.
+
+  // These will indicate if we've gotten our required parameters at the end of parsing.
+  bool haveMeshes = false;
+  bool haveInterpTag = false;
+
+  // Loop over the values in argv pulling out an parsing each one
+  int npos = 1;
+
+  if (argc > 1 && argv[1] == std::string("-h")) {
+    help = true;
+    return MB_SUCCESS;
+  }
+
+  while (npos < argc) {
+    if (argv[npos] == std::string("-meshes")) {
+      // Parse out the mesh filenames
+      npos++;
+      int numFiles = 2;
+      meshFiles.resize(numFiles);
+      for (int i = 0; i < numFiles; i++) {
+        if ((npos < argc) && (!check_for_flag(argv[npos])))
+          meshFiles[i] = argv[npos++];
+        else {
+          std::cerr << "    ERROR - missing correct number of mesh filenames" << std::endl;
+          return MB_FAILURE;
+        }
+      }
+
+      haveMeshes = true;
+    }
+    else if (argv[npos] == std::string("-itag")) {
+      // Parse out the interpolation tag
+      npos++;
+      if ((npos < argc) && (!check_for_flag(argv[npos])))
+        interpTag = argv[npos++];
+      else {
+        std::cerr << "    ERROR - missing <interp_tag>" << std::endl;
+        return MB_FAILURE;
+      }
+
+      haveInterpTag = true;
+    }
+    else if (argv[npos] == std::string("-meth")) {
+      // Parse out the interpolation tag
+      npos++;
+      if (argv[npos][0] == '0') method = DataCoupler::CONSTANT;
+      else if (argv[npos][0] == '1') method = DataCoupler::LINEAR_FE;
+      else if (argv[npos][0] == '2') method = DataCoupler::QUADRATIC_FE;
+      else if (argv[npos][0] == '3') method = DataCoupler::SPECTRAL;
+      else {
+        std::cerr << "    ERROR - unrecognized method number " << method << std::endl;
+        return MB_FAILURE;
+      }
+      npos++;
+    }
+    else if (argv[npos] == std::string("-eps")) {
+      // Parse out the tolerance
+      npos++;
+      if ((npos < argc) && (!check_for_flag(argv[npos])))
+        epsilon = atof(argv[npos++]);
+      else {
+        std::cerr << "    ERROR - missing <epsilon>" << std::endl;
+        return MB_FAILURE;
+      }
+    }
+    else if (argv[npos] == std::string("-gnorm")) {
+      // Parse out the global normalization tag
+      npos++;
+      if ((npos < argc) && (!check_for_flag(argv[npos])))
+        gNormTag = argv[npos++];
+      else {
+        std::cerr << "    ERROR - missing <gnorm_tag>" << std::endl;
+        return MB_FAILURE;
+      }
+    }
+    else if (argv[npos] == std::string("-ssnorm")) {
+      // Parse out the subset normalization tag and selection criteria
+      npos++;
+      if ((npos < argc) && (!check_for_flag(argv[npos])))
+        ssNormTag = argv[npos++];
+      else {
+        std::cerr << "    ERROR - missing <ssnorm_tag>" << std::endl;
+        return MB_FAILURE;
+      }
+
+      if ((npos < argc) && (!check_for_flag(argv[npos]))) {
+        char* opts = argv[npos++];
+        char sep1[1] = {';'};
+        char sep2[1] = {'='};
+        bool end_vals_seen = false;
+        std::vector<char *> tmpTagOpts;
+
+        // first get the options
+        for (char* i = strtok(opts, sep1); i; i = strtok(0, sep1)) {
+          tmpTagOpts.push_back(i);
+        }
+
+        // parse out the name and val or just name.
+        for (unsigned int j = 0; j < tmpTagOpts.size(); j++) {
+          char* e = strtok(tmpTagOpts[j], sep2);
+          ssTagNames.push_back(e);
+          e = strtok(0, sep2);
+          if (e != NULL) {
+            // We have a value
+            if (end_vals_seen) {
+              // ERROR we should not have a value after none are seen
+              std::cerr << "    ERROR - new value seen after end of values in <ssnorm_selection>" << std::endl;
+              return MB_FAILURE;
+            }
+            // Otherwise get the value string from e and convert it to an int
+            int *valp = new int;
+            *valp = atoi(e);
+            ssTagValues.push_back((const char *) valp);
+          }
+          else {
+            // Otherwise there is no '=' so push a null on the list
+            end_vals_seen = true;
+            ssTagValues.push_back((const char *) 0);
+          }
+        }
+      }
+      else {
+        std::cerr << "    ERROR - missing <ssnorm_selection>" << std::endl;
+        return MB_FAILURE;
+      }
+    }
+    else if (argv[npos] == std::string("-ropts")) {
+      // Parse out the mesh file read options
+      npos++;
+      if ((npos < argc) && (!check_for_flag(argv[npos])))
+        readOpts = argv[npos++];
+      else {
+        std::cerr << "    ERROR - missing <roptions>" << std::endl;
+        return MB_FAILURE;
+      }
+    }
+    else if (argv[npos] == std::string("-outfile")) {
+      // Parse out the output file name
+      npos++;
+      if ((npos < argc) && (!check_for_flag(argv[npos])))
+        outFile = argv[npos++];
+      else {
+        std::cerr << "    ERROR - missing <out_file>" << std::endl;
+        return MB_FAILURE;
+      }
+    }
+    else if (argv[npos] == std::string("-wopts")) {
+      // Parse out the output file write options
+      npos++;
+      if ((npos < argc) && (!check_for_flag(argv[npos])))
+        writeOpts = argv[npos++];
+      else {
+        std::cerr << "    ERROR - missing <woptions>" << std::endl;
+        return MB_FAILURE;
+      }
+    }
+    else if (argv[npos] == std::string("-dbgout")) {
+      // Parse out the debug output file name.
+      // If no name then use the default.
+      npos++;
+      if ((npos < argc) && (!check_for_flag(argv[npos])))
+        dbgFile = argv[npos++];
+      else
+        dbgFile = defaultDbgFile;
+    }
+    else {
+      // Unrecognized parameter.  Skip it and move along.
+      std::cerr << "    ERROR - Unrecognized parameter:" << argv[npos] << std::endl;
+      std::cerr << "            Skipping..." << std::endl;
+      npos++;
+    }
+  }
+
+  if (!haveMeshes) {
+    meshFiles.resize(2);
+    meshFiles[0] = std::string(TestDir + "/64bricks_1khex.h5m");
+    meshFiles[1] = std::string(TestDir + "/64bricks_12ktet.h5m");
+    std::cout << "Mesh files not entered; using default files " 
+              << meshFiles[0] << " and " << meshFiles[1] << std::endl;
+  }
+  
+  if (!haveInterpTag) {
+    interpTag = "vertex_field";
+    std::cout << "Interpolation field name not given, using default of " << interpTag << std::endl;
+  }
+
+#ifdef HDF5_FILE
+  if (1 == argc) {
+    std::cout << "No arguments given; using output file dum.h5m." << std::endl;
+    outFile = "dum.h5m";
+  }
+#endif
+    
+  return MB_SUCCESS;
+}
+
+
+// End new get_file_options()
+
+
+ErrorCode test_interpolation(Interface *mbImpl, 
+                             DataCoupler::Method method,
+                             std::string &interpTag,
+                             std::string &gNormTag,
+                             std::string &ssNormTag,
+                             std::vector<const char *> &ssTagNames,
+                             std::vector<const char *> &ssTagValues,
+                             iBase_EntitySetHandle *roots,
+                             std::vector<ParallelComm *> &pcs,
+                             double &instant_time,
+                             double &pointloc_time,
+                             double &interp_time,
+                             double &gnorm_time,
+                             double &ssnorm_time,
+                             double & toler)
+{
+  assert(method >= DataCoupler::CONSTANT && method <= DataCoupler::SPECTRAL);
+
+    // source is 1st mesh, target is 2nd
+  Range src_elems, targ_elems, targ_verts;
+  ErrorCode result = pcs[0]->get_part_entities(src_elems, 3);
+  PRINT_LAST_ERROR;
+
+  CpuTimer timer;
+
+    // instantiate a coupler, which also initializes the tree
+  DataCoupler dc(mbImpl, src_elems, 0, pcs[0]);
+
+  // initialize spectral elements, if they exist
+  bool specSou=false, specTar = false;
+//  result =  mbc.initialize_spectral_elements((EntityHandle)roots[0], (EntityHandle)roots[1], specSou, specTar);
+
+  instant_time = timer.time_since_birth();
+
+    // get points from the target mesh to interpolate
+  // we have to treat differently the case when the target is a spectral mesh
+  // in that case, the points of interest are the GL points, not the vertex nodes
+  std::vector<double> vpos; // this will have the positions we are interested in
+  int numPointsOfInterest = 0;
+#ifdef USE_MPI    
+  result = pcs[1]->get_part_entities(targ_elems, 3);
+#else
+  result = mbImpl->get_entities_by_dimension((EntityHandle)roots[1], 3);
+#endif
+  PRINT_LAST_ERROR;
+
+    // first get all vertices adj to partition entities in target mesh
+  if (DataCoupler::CONSTANT == method) 
+    targ_verts = targ_elems;
+  else
+    result = mbImpl->get_adjacencies(targ_elems, 0, false, targ_verts,
+                                     Interface::UNION);
+  PRINT_LAST_ERROR;
+
+#ifdef USE_MPI  
+    // then get non-owned verts and subtract
+  Range tmp_verts;
+  result = pcs[1]->get_pstatus_entities(0, PSTATUS_NOT_OWNED, tmp_verts);
+  PRINT_LAST_ERROR;
+  targ_verts = subtract( targ_verts, tmp_verts);
+#endif
+    // get position of these entities; these are the target points
+  numPointsOfInterest = (int)targ_verts.size();
+  vpos.resize(3*targ_verts.size());
+  result = mbImpl->get_coords(targ_verts, &vpos[0]);
+  PRINT_LAST_ERROR;
+
+    // locate those points in the source mesh
+#ifdef USE_MPI
+  std::cout << "rank "<< pcs[0]->proc_config().proc_rank();
+#endif
+  std::cout << " points of interest: " << numPointsOfInterest << "\n";
+  result = dc.locate_points(&vpos[0], numPointsOfInterest, toler);
+  PRINT_LAST_ERROR;
+
+  pointloc_time = timer.time_elapsed();
+
+    // now interpolate tag onto target points
+  std::vector<double> field(numPointsOfInterest);
+
+  result = dc.interpolate(method, interpTag, &field[0]);
+  PRINT_LAST_ERROR;
+  
+  interp_time = timer.time_elapsed();
+
+/*
+    // do global normalization if specified
+  if (!gNormTag.empty()) {
+    int err;
+
+    // Normalize the source mesh
+    err = dc.normalize_mesh(roots[0], gNormTag.c_str(), DataCoupler::VOLUME, 4);
+    PRINT_LAST_ERR;
+
+    // Normalize the target mesh
+    err = dc.normalize_mesh(roots[1], gNormTag.c_str(), DataCoupler::VOLUME, 4);
+    PRINT_LAST_ERR;
+  }
+
+  gnorm_time = timer.time_elapsed();
+
+    // do subset normalization if specified
+
+  if (!ssNormTag.empty()) {
+    int err;
+
+    err = dc.normalize_subset(roots[0], 
+                               ssNormTag.c_str(), 
+                               &ssTagNames[0], 
+                               ssTagNames.size(), 
+                               &ssTagValues[0], 
+                              DataCoupler::VOLUME, 
+                               4);
+    PRINT_LAST_ERR;
+
+    err = dc.normalize_subset(roots[1], 
+                               ssNormTag.c_str(), 
+                               &ssTagNames[0], 
+                               ssTagNames.size(), 
+                               &ssTagValues[0], 
+                               DataCoupler::VOLUME, 
+                               4);
+    PRINT_LAST_ERR;
+  }
+
+  ssnorm_time = timer.time_elapsed();
+*/
+    // set field values as tag on target vertices
+    // use original tag
+  Tag tag;
+  result = mbImpl->tag_get_handle(interpTag.c_str(), 1, MB_TYPE_DOUBLE, tag); PRINT_LAST_ERROR;
+  result = mbImpl->tag_set_data(tag, targ_verts, &field[0]);
+  PRINT_LAST_ERROR;
+
+    // done
+  return MB_SUCCESS;
+}


https://bitbucket.org/fathomteam/moab/commits/cfbad175c3f5/
Changeset:   cfbad175c3f5
Branch:      None
User:        tautges
Date:        2014-03-12 16:59:30
Summary:     - add ability to specify master and slave mesh solid/fluid set numbers separately, so set numbers don't need to agree
- add ability to specify only fluid or solid sets, with other sets derived from that and all material sets in a file
- add DeformMeshRemap to list of examples in makefile

Affected #:  1 file

diff --git a/examples/DeformMeshRemap.cpp b/examples/DeformMeshRemap.cpp
index a8a0e10..c5e2452 100644
--- a/examples/DeformMeshRemap.cpp
+++ b/examples/DeformMeshRemap.cpp
@@ -358,6 +358,50 @@ ErrorCode DeformMeshRemap::execute()
   return MB_SUCCESS;
 }
 
+ErrorCode DeformMeshRemap::find_other_sets(int m_or_s, EntityHandle file_set) 
+{
+    // solid or fluid sets are missing; find the other
+  Range *filled_sets = NULL, *unfilled_sets = NULL, *unfilled_elems = NULL;
+  
+  if (fluidSets[m_or_s].empty() && !solidSets[m_or_s].empty()) {
+    unfilled_sets = &fluidSets[m_or_s];
+    filled_sets = &solidSets[m_or_s];
+    unfilled_elems = &fluidElems[m_or_s];
+  }
+  else if (!fluidSets[m_or_s].empty() && solidSets[m_or_s].empty()) {
+    filled_sets = &fluidSets[m_or_s];
+    unfilled_sets = &solidSets[m_or_s];
+    unfilled_elems = &solidElems[m_or_s];
+  }
+  
+    // ok, we know the filled sets, now fill the unfilled sets, and the elems from those
+  Tag tagh;
+  ErrorCode rval = mbImpl->tag_get_handle(MATERIAL_SET_TAG_NAME, tagh); RR("Couldn't get material set tag name.");
+  Range matsets;
+  rval = mbImpl->get_entities_by_type_and_tag(file_set, MBENTITYSET, &tagh, NULL, 1, matsets);
+  if (matsets.empty()) rval = MB_FAILURE;
+  RR("Couldn't get any material sets.");
+  *unfilled_sets = subtract(matsets, *filled_sets);
+  if (unfilled_sets->empty()) {
+    rval = MB_FAILURE;
+    RR("Failed to find any unfilled material sets.");
+  }
+  Range tmp_range;
+  for (Range::iterator rit = unfilled_sets->begin(); rit != unfilled_sets->end(); rit++) {
+    rval = mbImpl->get_entities_by_handle(*rit, tmp_range, true);
+    RR("Failed to get entities in unfilled set.");
+  }
+  int dim = mbImpl->dimension_from_handle(*tmp_range.rbegin());
+  assert(dim > 0 && dim < 4);  
+  *unfilled_elems = tmp_range.subset_by_dimension(dim);
+  if (unfilled_elems->empty()) {
+    rval = MB_FAILURE;
+    RR("Failed to find any unfilled set entities.");
+  }
+  
+  return MB_SUCCESS;
+}
+
 std::string DeformMeshRemap::get_file_name(int m_or_s) const
 {
   switch (m_or_s) {


https://bitbucket.org/fathomteam/moab/commits/321ee0341127/
Changeset:   321ee0341127
Branch:      None
User:        tautges
Date:        2014-03-12 17:00:26
Summary:     Better error handling & debug printing.
In Tree, handle legacy tag on tree root by deleting and re-creating (legacy files have old-sized tag).

Affected #:  1 file

diff --git a/examples/DeformMeshRemap.cpp b/examples/DeformMeshRemap.cpp
index c5e2452..a8a0e10 100644
--- a/examples/DeformMeshRemap.cpp
+++ b/examples/DeformMeshRemap.cpp
@@ -358,50 +358,6 @@ ErrorCode DeformMeshRemap::execute()
   return MB_SUCCESS;
 }
 
-ErrorCode DeformMeshRemap::find_other_sets(int m_or_s, EntityHandle file_set) 
-{
-    // solid or fluid sets are missing; find the other
-  Range *filled_sets = NULL, *unfilled_sets = NULL, *unfilled_elems = NULL;
-  
-  if (fluidSets[m_or_s].empty() && !solidSets[m_or_s].empty()) {
-    unfilled_sets = &fluidSets[m_or_s];
-    filled_sets = &solidSets[m_or_s];
-    unfilled_elems = &fluidElems[m_or_s];
-  }
-  else if (!fluidSets[m_or_s].empty() && solidSets[m_or_s].empty()) {
-    filled_sets = &fluidSets[m_or_s];
-    unfilled_sets = &solidSets[m_or_s];
-    unfilled_elems = &solidElems[m_or_s];
-  }
-  
-    // ok, we know the filled sets, now fill the unfilled sets, and the elems from those
-  Tag tagh;
-  ErrorCode rval = mbImpl->tag_get_handle(MATERIAL_SET_TAG_NAME, tagh); RR("Couldn't get material set tag name.");
-  Range matsets;
-  rval = mbImpl->get_entities_by_type_and_tag(file_set, MBENTITYSET, &tagh, NULL, 1, matsets);
-  if (matsets.empty()) rval = MB_FAILURE;
-  RR("Couldn't get any material sets.");
-  *unfilled_sets = subtract(matsets, *filled_sets);
-  if (unfilled_sets->empty()) {
-    rval = MB_FAILURE;
-    RR("Failed to find any unfilled material sets.");
-  }
-  Range tmp_range;
-  for (Range::iterator rit = unfilled_sets->begin(); rit != unfilled_sets->end(); rit++) {
-    rval = mbImpl->get_entities_by_handle(*rit, tmp_range, true);
-    RR("Failed to get entities in unfilled set.");
-  }
-  int dim = mbImpl->dimension_from_handle(*tmp_range.rbegin());
-  assert(dim > 0 && dim < 4);  
-  *unfilled_elems = tmp_range.subset_by_dimension(dim);
-  if (unfilled_elems->empty()) {
-    rval = MB_FAILURE;
-    RR("Failed to find any unfilled set entities.");
-  }
-  
-  return MB_SUCCESS;
-}
-
 std::string DeformMeshRemap::get_file_name(int m_or_s) const
 {
   switch (m_or_s) {


https://bitbucket.org/fathomteam/moab/commits/7c68d61d2b8f/
Changeset:   7c68d61d2b8f
Branch:      None
User:        tautges
Date:        2014-03-12 17:01:17
Summary:     Changes to allow specification of existing displacement tag(s), which can be either single name corresponding to xyz[3]
or 3 names x, y, z.  Handle either case in deform_master.

Affected #:  1 file

diff --git a/examples/DeformMeshRemap.cpp b/examples/DeformMeshRemap.cpp
index a8a0e10..344565a 100644
--- a/examples/DeformMeshRemap.cpp
+++ b/examples/DeformMeshRemap.cpp
@@ -387,7 +387,6 @@ DeformMeshRemap::DeformMeshRemap(Interface *impl, ParallelComm *master, Parallel
         : mbImpl(impl), pcMaster(master), pcSlave(slave), masterSet(0), slaveSet(0), xNew(0), xNewName("xnew") 
 {
   mbImpl->query_interface(mError);
-  
   xDisp[0] = xDisp[1] = xDisp[2] = 0;
 
   if (!pcSlave && pcMaster)
@@ -553,7 +552,7 @@ ErrorCode DeformMeshRemap::deform_master(Range &fluid_elems, Range &solid_elems,
     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.");
     std::vector<double> disps(num_verts);
-    
+
     for (int i = 0; i < 3; i++) {
       rval = mbImpl->tag_get_handle(xDispNames[0].c_str(), 1, MB_TYPE_DOUBLE, xDisp[i]);
       RR("Failed to get xDisp tag.");


https://bitbucket.org/fathomteam/moab/commits/fbad7648f20f/
Changeset:   fbad7648f20f
Branch:      None
User:        tautges
Date:        2014-03-12 17:01:17
Summary:     - add ability to specify master and slave mesh solid/fluid set numbers separately, so set numbers don't need to agree
- add ability to specify only fluid or solid sets, with other sets derived from that and all material sets in a file
- add DeformMeshRemap to list of examples in makefile

Affected #:  1 file

diff --git a/examples/DeformMeshRemap.cpp b/examples/DeformMeshRemap.cpp
index 344565a..59e4469 100644
--- a/examples/DeformMeshRemap.cpp
+++ b/examples/DeformMeshRemap.cpp
@@ -358,6 +358,50 @@ ErrorCode DeformMeshRemap::execute()
   return MB_SUCCESS;
 }
 
+ErrorCode DeformMeshRemap::find_other_sets(int m_or_s, EntityHandle file_set) 
+{
+    // solid or fluid sets are missing; find the other
+  Range *filled_sets = NULL, *unfilled_sets = NULL, *unfilled_elems = NULL;
+  
+  if (fluidSets[m_or_s].empty() && !solidSets[m_or_s].empty()) {
+    unfilled_sets = &fluidSets[m_or_s];
+    filled_sets = &solidSets[m_or_s];
+    unfilled_elems = &fluidElems[m_or_s];
+  }
+  else if (!fluidSets[m_or_s].empty() && solidSets[m_or_s].empty()) {
+    filled_sets = &fluidSets[m_or_s];
+    unfilled_sets = &solidSets[m_or_s];
+    unfilled_elems = &solidElems[m_or_s];
+  }
+  
+    // ok, we know the filled sets, now fill the unfilled sets, and the elems from those
+  Tag tagh;
+  ErrorCode rval = mbImpl->tag_get_handle(MATERIAL_SET_TAG_NAME, tagh); RR("Couldn't get material set tag name.");
+  Range matsets;
+  rval = mbImpl->get_entities_by_type_and_tag(file_set, MBENTITYSET, &tagh, NULL, 1, matsets);
+  if (matsets.empty()) rval = MB_FAILURE;
+  RR("Couldn't get any material sets.");
+  *unfilled_sets = subtract(matsets, *filled_sets);
+  if (unfilled_sets->empty()) {
+    rval = MB_FAILURE;
+    RR("Failed to find any unfilled material sets.");
+  }
+  Range tmp_range;
+  for (Range::iterator rit = unfilled_sets->begin(); rit != unfilled_sets->end(); rit++) {
+    rval = mbImpl->get_entities_by_handle(*rit, tmp_range, true);
+    RR("Failed to get entities in unfilled set.");
+  }
+  int dim = mbImpl->dimension_from_handle(*tmp_range.rbegin());
+  assert(dim > 0 && dim < 4);  
+  *unfilled_elems = tmp_range.subset_by_dimension(dim);
+  if (unfilled_elems->empty()) {
+    rval = MB_FAILURE;
+    RR("Failed to find any unfilled set entities.");
+  }
+  
+  return MB_SUCCESS;
+}
+
 std::string DeformMeshRemap::get_file_name(int m_or_s) const
 {
   switch (m_or_s) {


https://bitbucket.org/fathomteam/moab/commits/8c76f8fc316e/
Changeset:   8c76f8fc316e
Branch:      None
User:        tautges
Date:        2014-03-12 17:01:17
Summary:     Better error handling & debug printing.
In Tree, handle legacy tag on tree root by deleting and re-creating (legacy files have old-sized tag).

Affected #:  1 file

diff --git a/examples/DeformMeshRemap.cpp b/examples/DeformMeshRemap.cpp
index 59e4469..344565a 100644
--- a/examples/DeformMeshRemap.cpp
+++ b/examples/DeformMeshRemap.cpp
@@ -358,50 +358,6 @@ ErrorCode DeformMeshRemap::execute()
   return MB_SUCCESS;
 }
 
-ErrorCode DeformMeshRemap::find_other_sets(int m_or_s, EntityHandle file_set) 
-{
-    // solid or fluid sets are missing; find the other
-  Range *filled_sets = NULL, *unfilled_sets = NULL, *unfilled_elems = NULL;
-  
-  if (fluidSets[m_or_s].empty() && !solidSets[m_or_s].empty()) {
-    unfilled_sets = &fluidSets[m_or_s];
-    filled_sets = &solidSets[m_or_s];
-    unfilled_elems = &fluidElems[m_or_s];
-  }
-  else if (!fluidSets[m_or_s].empty() && solidSets[m_or_s].empty()) {
-    filled_sets = &fluidSets[m_or_s];
-    unfilled_sets = &solidSets[m_or_s];
-    unfilled_elems = &solidElems[m_or_s];
-  }
-  
-    // ok, we know the filled sets, now fill the unfilled sets, and the elems from those
-  Tag tagh;
-  ErrorCode rval = mbImpl->tag_get_handle(MATERIAL_SET_TAG_NAME, tagh); RR("Couldn't get material set tag name.");
-  Range matsets;
-  rval = mbImpl->get_entities_by_type_and_tag(file_set, MBENTITYSET, &tagh, NULL, 1, matsets);
-  if (matsets.empty()) rval = MB_FAILURE;
-  RR("Couldn't get any material sets.");
-  *unfilled_sets = subtract(matsets, *filled_sets);
-  if (unfilled_sets->empty()) {
-    rval = MB_FAILURE;
-    RR("Failed to find any unfilled material sets.");
-  }
-  Range tmp_range;
-  for (Range::iterator rit = unfilled_sets->begin(); rit != unfilled_sets->end(); rit++) {
-    rval = mbImpl->get_entities_by_handle(*rit, tmp_range, true);
-    RR("Failed to get entities in unfilled set.");
-  }
-  int dim = mbImpl->dimension_from_handle(*tmp_range.rbegin());
-  assert(dim > 0 && dim < 4);  
-  *unfilled_elems = tmp_range.subset_by_dimension(dim);
-  if (unfilled_elems->empty()) {
-    rval = MB_FAILURE;
-    RR("Failed to find any unfilled set entities.");
-  }
-  
-  return MB_SUCCESS;
-}
-
 std::string DeformMeshRemap::get_file_name(int m_or_s) const
 {
   switch (m_or_s) {


https://bitbucket.org/fathomteam/moab/commits/12bd6aeace15/
Changeset:   12bd6aeace15
Branch:      None
User:        tautges
Date:        2014-03-12 17:01:17
Summary:     DeformMeshRemap: more debugging output, and handle coord displacements properly (add to prev coords)
LloydSmoother: fix a bug with fixed tag entity handles size
SpatialLocator: if an element range is specified, build the tree too

Affected #:  1 file

diff --git a/examples/DeformMeshRemap.cpp b/examples/DeformMeshRemap.cpp
index 344565a..fee03b8 100644
--- a/examples/DeformMeshRemap.cpp
+++ b/examples/DeformMeshRemap.cpp
@@ -552,7 +552,6 @@ ErrorCode DeformMeshRemap::deform_master(Range &fluid_elems, Range &solid_elems,
     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.");
     std::vector<double> disps(num_verts);
-
     for (int i = 0; i < 3; i++) {
       rval = mbImpl->tag_get_handle(xDispNames[0].c_str(), 1, MB_TYPE_DOUBLE, xDisp[i]);
       RR("Failed to get xDisp tag.");


https://bitbucket.org/fathomteam/moab/commits/854a59211ad5/
Changeset:   854a59211ad5
Branch:      None
User:        tautges
Date:        2014-03-12 17:01:17
Summary:     DeformMeshRemap.cpp: Write more intermediate files if compiled with local debug flag true; change some of the functions for moving coords to tags and vice versa to facilitate this.
SpatialLocator.cpp: Back out change that automatically chose BVH tree for non-vertex entities, there's a bug in point location there currently...

Affected #:  1 file

diff --git a/src/SpatialLocator.cpp b/src/SpatialLocator.cpp
index 702399e..6046fb2 100644
--- a/src/SpatialLocator.cpp
+++ b/src/SpatialLocator.cpp
@@ -35,12 +35,12 @@ namespace moab
     {
       if (myTree) return;
       
-      if (myElems.empty() || mbImpl->type_from_handle(*myElems.rbegin()) == MBVERTEX) 
+//      if (myElems.empty() || mbImpl->type_from_handle(*myElems.rbegin()) == MBVERTEX) 
           // create a kdtree if only vertices
         myTree = new AdaptiveKDTree(mbImpl);
-      else
+//      else
           // otherwise a BVHtree, since it performs better for elements
-        myTree = new BVHTree(mbImpl);
+//        myTree = new BVHTree(mbImpl);
 
       iCreatedTree = true;
     }


https://bitbucket.org/fathomteam/moab/commits/1fa59bf61d07/
Changeset:   1fa59bf61d07
Branch:      None
User:        tautges
Date:        2014-03-12 17:01:17
Summary:     Various smallish bug fixes and minor changes, improving the use of SpatialLocator.

ElemEvaluator: in set_ent_handle, if there's no evalSet for that entity type, call EvalSet::get_eval_set
SpatialLocator: set the evaluator on the tree if it's not the same as what's on SpatialLocator; also, replace find_point logic
  (that used distance_search for some reason) with call to Tree::point_search.
AdaptiveKDTree: move parse_options to public function, as it is in Tree
BVHTree: fixed bug, where BVHTree had its own myEval shadowing the one declared in Tree
Tree: fixed bug, initialize myEval to NULL in constructor
spatial_locator_test: test trees with and without evaluators
DataCoupler: improving comments a bit, and move ParallelComm down in arg list for constructor; also, some
  improved comments
datacoupler_test: adding test


Passes make check in parallel.

Affected #:  1 file
Diff not available.

https://bitbucket.org/fathomteam/moab/commits/9c9c716a8e6b/
Changeset:   9c9c716a8e6b
Branch:      deformed-mesh-remap
User:        tautges
Date:        2014-03-12 17:05:39
Summary:     Merge branch 'deformed-mesh-remap' of bitbucket.org:fathomteam/moab into deformed-mesh-remap

Conflicts:
	examples/DeformMeshRemap.cpp

Affected #:  1 file
Diff not available.

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