[MOAB-dev] commit/MOAB: 2 new changesets
commits-noreply at bitbucket.org
commits-noreply at bitbucket.org
Tue Mar 11 15:34:35 CDT 2014
2 new commits in MOAB:
https://bitbucket.org/fathomteam/moab/commits/fc3351604b22/
Changeset: fc3351604b22
Branch: None
User: tautges
Date: 2014-02-20 19:44:45
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 a8a8f07..44f8494 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/41dc120f7e2c/
Changeset: 41dc120f7e2c
Branch: deformed-mesh-remap
User: tautges
Date: 2014-03-11 21:31:21
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 44f8494..6f3cf22 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;
}
@@ -365,90 +365,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;
+}
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