[MOAB-dev] commit/MOAB: 8 new changesets
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Thu Feb 13 11:09:52 CST 2014
8 new commits in MOAB:
https://bitbucket.org/fathomteam/moab/commits/9a7a7209f500/
Changeset: 9a7a7209f500
Branch: None
User: tautges
Date: 2014-02-13 18:06:29
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/05548459620a/
Changeset: 05548459620a
Branch: None
User: tautges
Date: 2014-02-13 18:06:29
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 #: 2 files
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;
diff --git a/examples/makefile b/examples/makefile
index 6c1ed9a..26e3592 100644
--- a/examples/makefile
+++ b/examples/makefile
@@ -7,7 +7,7 @@ include ${MOAB_DIR}/lib/iMesh-Defs.inc
# MESH_DIR is the directory containing mesh files that come with MOAB source
MESH_DIR="../MeshFiles/unittest"
-EXAMPLES = HelloMOAB GetEntities SetsNTags StructuredMeshSimple DirectAccessWithHoles DirectAccessNoHoles
+EXAMPLES = HelloMOAB GetEntities SetsNTags StructuredMeshSimple DirectAccessWithHoles DirectAccessNoHoles DeformMeshRemap
PAREXAMPLES = HelloParMOAB ReduceExchangeTags LloydRelaxation
F90EXAMPLES = DirectAccessNoHolesF90 PushParMeshIntoMoabF90
EXOIIEXAMPLES = TestExodusII
https://bitbucket.org/fathomteam/moab/commits/d716d5ca9799/
Changeset: d716d5ca9799
Branch: None
User: tautges
Date: 2014-02-13 18:06:29
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/076554e403f6/
Changeset: 076554e403f6
Branch: None
User: tautges
Date: 2014-02-13 18:06:29
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/c953d6fe171a/
Changeset: c953d6fe171a
Branch: None
User: tautges
Date: 2014-02-13 18:06:29
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/dad32f253cfc/
Changeset: dad32f253cfc
Branch: None
User: tautges
Date: 2014-02-13 18:06:29
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 36016c1..c56c9ef 100644
--- a/src/SpatialLocator.cpp
+++ b/src/SpatialLocator.cpp
@@ -30,7 +30,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/b0b14cc655e6/
Changeset: b0b14cc655e6
Branch: None
User: tautges
Date: 2014-02-13 18:06:29
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/edc9cad082de/
Changeset: edc9cad082de
Branch: deformed-mesh-remap
User: tautges
Date: 2014-02-13 18:08:00
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 #: 9 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 b65eae4..3ad8d72 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 c56c9ef..f487e39 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"
bool debug = true;
@@ -11,10 +12,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);
@@ -22,6 +21,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 5c51341..303039e 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 ¢er);
void update(const BoundBox &other_box);
void update(const double *coords);
@@ -74,10 +74,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/SpatialLocator.hpp b/src/moab/SpatialLocator.hpp
index e8ab5cf..3834588 100644
--- a/src/moab/SpatialLocator.hpp
+++ b/src/moab/SpatialLocator.hpp
@@ -157,6 +157,18 @@ namespace moab {
private:
+ /* locate a point */
+ ErrorCode locate_point(const double *pos,
+ EntityHandle &ent, double *params, bool *is_inside = NULL,
+ const double rel_iter_tol = 1.0e-10, const double abs_iter_tol = 1.0e-10,
+ const double inside_tol = 1.0e-6);
+
+ /** 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;
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
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