[MOAB-dev] r1931 - in MOAB/trunk: . parallel tools/mbcoupler
tautges at mcs.anl.gov
tautges at mcs.anl.gov
Thu Jun 26 09:55:10 CDT 2008
Author: tautges
Date: 2008-06-26 09:55:10 -0500 (Thu, 26 Jun 2008)
New Revision: 1931
Modified:
MOAB/trunk/MBGeomUtil.hpp
MOAB/trunk/mbparallelcomm_test.cpp
MOAB/trunk/parallel/MBParallelComm.cpp
MOAB/trunk/parallel/MBParallelComm.hpp
MOAB/trunk/tools/mbcoupler/MBCoupler.cpp
MOAB/trunk/tools/mbcoupler/MBCoupler.hpp
MOAB/trunk/tools/mbcoupler/MBElemUtil.cpp
MOAB/trunk/tools/mbcoupler/MBElemUtil.hpp
MOAB/trunk/tools/mbcoupler/mbcoupler_test.cpp
Log:
- short-circuiting MBElemUtil::nat_coords_trilinear_hex with .5, .5,
.5 for now (doesn't work at vertex locations)
- revising MBCoupler functions to make them compile and (eventually?)
work
- adding calls to coupler in mbcouplercomm_test
- adding MBParallelComm::get_part_entities to get partition entities
of specified dimension on this parallelcomm
- adding a few more diagnostics to mbparallelcomm_test
Modified: MOAB/trunk/MBGeomUtil.hpp
===================================================================
--- MOAB/trunk/MBGeomUtil.hpp 2008-06-26 00:34:42 UTC (rev 1930)
+++ MOAB/trunk/MBGeomUtil.hpp 2008-06-26 14:55:10 UTC (rev 1931)
@@ -256,7 +256,7 @@
// and checks if each are between +/-1. If anyone is outside the range
// the function returns false, otherwise it returns true.
//
-bool point_in_trilinear_hex(MBCartVect hex[8],
+bool point_in_trilinear_hex(MBCartVect *hex,
MBCartVect xyz,
double etol);
@@ -267,7 +267,7 @@
// test above except that it gets a bounding box as arguments to filter
// the test for acceleration.
//
-bool point_in_trilinear_hex(MBCartVect hex[8],
+bool point_in_trilinear_hex(MBCartVect *hex,
MBCartVect xyz,
MBCartVect box_min, MBCartVect box_max,
double etol);
Modified: MOAB/trunk/mbparallelcomm_test.cpp
===================================================================
--- MOAB/trunk/mbparallelcomm_test.cpp 2008-06-26 00:34:42 UTC (rev 1930)
+++ MOAB/trunk/mbparallelcomm_test.cpp 2008-06-26 14:55:10 UTC (rev 1931)
@@ -176,18 +176,23 @@
MBParallelComm *pcomm = MBParallelComm::get_pcomm(mbImpl, 0);
assert(pcomm);
- MBRange iface_ents[6];
+ MBRange iface_ents[7];
for (int i = 0; i < 4; i++) {
tmp_result = pcomm->get_iface_entities(-1, i, iface_ents[i]);
if (MB_SUCCESS != tmp_result) {
std::cerr << "get_iface_entities returned error on proc "
<< rank << "; message: " << std::endl;
- PRINT_LAST_ERROR
- result = tmp_result;
+ PRINT_LAST_ERROR;
+ result = tmp_result;
}
if (0 != i) iface_ents[4].merge(iface_ents[i]);
}
+ result = pcomm->get_part_entities(iface_ents[6], -1);
+ PRINT_LAST_ERROR;
+
+ std::cerr << "Proc " << rank << " partition entities:" << std::endl;
+ iface_ents[6].print(" ");
if (0 == rank) setime = MPI_Wtime();
Modified: MOAB/trunk/parallel/MBParallelComm.cpp
===================================================================
--- MOAB/trunk/parallel/MBParallelComm.cpp 2008-06-26 00:34:42 UTC (rev 1930)
+++ MOAB/trunk/parallel/MBParallelComm.cpp 2008-06-26 14:55:10 UTC (rev 1931)
@@ -3421,6 +3421,26 @@
return pc_array[index];
}
+ //! return all the entities in parts owned locally
+MBErrorCode MBParallelComm::get_part_entities(MBRange &ents, int dim)
+{
+ MBErrorCode result;
+
+ for (MBRange::iterator rit = partitionSets.begin();
+ rit != partitionSets.end(); rit++) {
+ MBRange tmp_ents;
+ if (-1 == dim)
+ result = mbImpl->get_entities_by_handle(*rit, tmp_ents, true);
+ else
+ result = mbImpl->get_entities_by_dimension(*rit, dim, tmp_ents, true);
+
+ if (MB_SUCCESS != result) return result;
+ ents.merge(tmp_ents);
+ }
+
+ return MB_SUCCESS;
+}
+
#ifdef TEST_PARALLELCOMM
#include <iostream>
Modified: MOAB/trunk/parallel/MBParallelComm.hpp
===================================================================
--- MOAB/trunk/parallel/MBParallelComm.hpp 2008-06-26 00:34:42 UTC (rev 1930)
+++ MOAB/trunk/parallel/MBParallelComm.hpp 2008-06-26 14:55:10 UTC (rev 1931)
@@ -294,6 +294,9 @@
//! return partitions set tag
MBTag partition_tag();
+
+ //! return all the entities in parts owned locally
+ MBErrorCode get_part_entities(MBRange &ents, int dim = -1);
private:
Modified: MOAB/trunk/tools/mbcoupler/MBCoupler.cpp
===================================================================
--- MOAB/trunk/tools/mbcoupler/MBCoupler.cpp 2008-06-26 00:34:42 UTC (rev 1930)
+++ MOAB/trunk/tools/mbcoupler/MBCoupler.cpp 2008-06-26 14:55:10 UTC (rev 1931)
@@ -1,11 +1,20 @@
#include "MBCoupler.hpp"
#include "MBParallelComm.hpp"
+#include "MBAdaptiveKDTree.hpp"
+#include "MBGeomUtil.hpp"
+#include "MBElemUtil.hpp"
+#include "MBCN.hpp"
+#include "iostream"
+
+extern "C"
+{
#include "types.h"
#include "errmem.h"
#include "minmax.h"
#include "sort.h"
#include "tuple_list.h"
#include "transfer.h"
+}
#include "assert.h"
@@ -13,17 +22,20 @@
MBParallelComm *pc,
MBRange &local_elems,
int coupler_id,
- bool init_tree)
+ bool init_tree)
: mbImpl(impl), myPc(pc), myId(coupler_id)
{
assert(NULL != impl && NULL != myPc);
-
+
// keep track of the local points, at least for now
myRange = local_elems;
-
+
// now initialize the tree
if (init_tree) initialize_tree();
+
+
+
// initialize tuple lists to indicate not initialized
mappedPts = NULL;
targetPts = NULL;
@@ -31,40 +43,99 @@
/* destructor
*/
-MBCoupler::~MBCoupler()
+MBCoupler::~MBCoupler()
{}
+MBErrorCode MBCoupler::initialize_tree()
+{
+
+ MBRange local_ents;
+ MBAdaptiveKDTree::Settings settings;
+ allBoxes.resize(6*myPc->proc_config().proc_size());
+
+ //get entities on the local part
+ MBErrorCode result = myPc->get_part_entities(local_ents, 3);
+ if (MB_SUCCESS != result) {
+ std::cout << "Problems getting entities by dimension" << std::endl;
+ return result;
+ }
+
+ // build the tree for local processor
+ myTree = new MBAdaptiveKDTree(mbImpl);
+ result = myTree->build_tree(local_ents, localRoot, &settings);
+ if (MB_SUCCESS != result) {
+ std::cout << "Problems building tree" << std::endl;
+ return result;
+ }
+
+ // get the bounding box for local tree
+ allBoxes.resize(6*myPc->proc_config().proc_size());
+ unsigned int my_rank = myPc->proc_config().proc_rank();
+ result = myTree->get_tree_box(localRoot, &allBoxes[6*my_rank], &allBoxes[6*my_rank+3]);
+
+ // now communicate to get all boxes
+ int mpi_err = MPI_Allgather(&allBoxes[6*my_rank], 6, MPI_DOUBLE,
+ &allBoxes[0], 6, MPI_DOUBLE,
+ myPc->proc_config().proc_comm());
+ if (MPI_SUCCESS == mpi_err) return MB_SUCCESS;
+ else return MB_FAILURE;
+}
+
MBErrorCode MBCoupler::locate_points(double *xyz, int num_points,
tuple_list *tl,
- bool store_local)
+ bool store_local)
{
assert(tl || store_local);
-
+
// allocate tuple_list to hold point data: (p, i, , xyz), i = point index
tuple_list target_pts;
tuple_list_init_max(&target_pts, 2, 0, 0, 3, 3*num_points);
- // for each point, find box(es) containing the point,
- // appending results to tuple_list;
+ // initialize source_pts and local_pts
+ tuple_list source_pts;
+ mappedPts = new tuple_list;
+ tuple_list_init_max(&source_pts, 3, 0, 0, 0, target_pts.n);
+ tuple_list_init_max(mappedPts, 0, 0, 1, 3, target_pts.n);
+
+ mappedPts->n = 0;
+ source_pts.n = 0;
+ MBErrorCode result;
+
+ // for each point, find box(es) containing the point,
+ // appending results to tuple_list;
// keep local points separately, in local_pts, which has pairs
// of <local_index, mapped_index>, where mapped_index is the index
+ // of <local_index, mapped_index>, where mapped_index is the index
// into the mappedPts tuple list
- int threen = 0;
- for (int i = 0; i < 3*num_points; i+=3) {
- for (;;/* <marco - loop over boxes> */ ) {
- if (/* marco - test point i/3 in box j */ false) {
- // check size, grow if we're at max
- if (target_pts.n == target_pts.max)
- tuple_list_grow(&target_pts);
+
+ for (int i = 0; i < 3*num_points; i+=3)
+ {
+ for (unsigned int j = 0; j < myPc->proc_config().proc_size(); j++)
+ {
+ // test if point is in box
+ if (allBoxes[6*j] <= xyz[i] && xyz[i] <= allBoxes[6*j+3] &&
+ allBoxes[6*j+1] <= xyz[i+1] && xyz[i+1] <= allBoxes[6*j+4] &&
+ allBoxes[6*j+2] <= xyz[i+2] && xyz[i+2] <= allBoxes[6*j+5])
+ {
+ if (j == myPc->proc_config().proc_rank())
+ {
+ result = test_local_box(xyz+i, j, i/3, i/3);
+ }
+ else {
+ // check size, grow if we're at max
+ if (target_pts.n == target_pts.max)
+ tuple_list_grow(&target_pts);
- target_pts.vi[2*target_pts.n] = -1 /* <marco - proc rank j> */;
- target_pts.vi[2*target_pts.n+1] = i/3;
- target_pts.vr[threen] = xyz[i];
- target_pts.vr[threen] = xyz[i+1];
- target_pts.vr[threen] = xyz[i+2];
- target_pts.n++;
- threen += 3;
+ target_pts.n++;
+
+ target_pts.vi[2*target_pts.n] = j;
+ target_pts.vi[2*target_pts.n+1] = i/3;
+
+ target_pts.vr[3*target_pts.n] = xyz[i];
+ target_pts.vr[3*target_pts.n+1] = xyz[i+1];
+ target_pts.vr[3*target_pts.n+2] = xyz[i+2];
+ }
}
}
}
@@ -72,12 +143,12 @@
// perform scatter/gather, to gather points to source mesh procs
gs_transfer(1, &target_pts, 0, myPc->proc_config().crystal_router());
- // after scatter/gather:
+ // after scatter/gather:
// target_pts.n = # points local proc has to map
// target_pts.vi[2*i] = proc sending point i
// target_pts.vi[2*i+1] = index of point i on sending proc
// target_pts.vr[3*i..3*i+2] = xyz of point i
- //
+ //
// Mapping builds the tuple list:
// source_pts.n = target_pts.n
// source_pts.vi[3*i] = target_pts.vi[2*i] = sending proc
@@ -89,148 +160,270 @@
// mappedPts->vul[i] = local handle of mapped entity
// mappedPts->vr[3*i..3*i+2] = natural coordinates in mapped entity
- // initialize source_pts and local_pts
- tuple_list source_pts;
- mappedPts = new tuple_list;
- tuple_list_init_max(&source_pts, 3, 0, 0, 0, target_pts.n);
- tuple_list_init_max(mappedPts, 0, 0, 1, 3, target_pts.n);
-
- for (unsigned i = 0; i < target_pts.n; i++) {
- // find leaf node(s) in local kdtree containing point i
- // <marco...>
-
- // find natural coordinates of point in element(s) in that leaf
- // <marco...>
-
- // for any point/element with nat coords within bounds, store
- // handle/nat coords in mappedPts, and proc/index in outgoing tuple
- // (-1 for index if no elements containing that point)
- // <marco...>
-
+ for (unsigned i = 0; i < target_pts.n; i++)
+ {
+ result = test_local_box(target_pts.vr+3*i,
+ target_pts.vi[3*i], target_pts.vi[3*i+1], i,
+ &source_pts);
+ if (MB_SUCCESS != result) return result;
}
- // no longer need target_pts
+ // no longer need target_pts
tuple_list_free(&target_pts);
- // perform scatter/gather to send proc/index tuples back to procs
+ // perform scatter/gather to send proc/index tuples back to procs
gs_transfer(1, &source_pts, 0, myPc->proc_config().crystal_router());
- // store proc/index tuples in targetPts, and/or pass back to application;
- // the tuple this gets stored to looks like:
- // tl.n = # mapped points
- // tl.vi[3*i] = remote proc mapping point
- // tl.vi[3*i+1] = local index of mapped point
- // tl.vi[3*i+2] = remote index of mapped point
- //
- // Local index is mapped into either myRange, holding the handles of
- // local mapped entities, or myXyz, holding locations of mapped pts
+ // store proc/index tuples in targetPts, and/or pass back to application;
+ // the tuple this gets stored to looks like:
+ // tl.n = # mapped points
+ // tl.vi[3*i] = remote proc mapping point
+ // tl.vi[3*i+1] = local index of mapped point
+ // tl.vi[3*i+2] = remote index of mapped point
+ //
+ // Local index is mapped into either myRange, holding the handles of
+ // local mapped entities, or myXyz, holding locations of mapped pts
- // go through and count non-negatives
+ // go through and count non-negatives
int num_pts = 0;
- for (unsigned int i = 0; i < source_pts.n; i++)
+ for (unsigned int i = 0; i < source_pts.n; i++)
if (-1 != source_pts.vi[3*i+2]) num_pts++;
targetPts = new tuple_list;
tuple_list *tl_tmp = targetPts;
if (!store_local) tl_tmp = tl;
tuple_list_init_max(tl_tmp, 3, 0, 0, 1, num_pts);
- tl_tmp->n = num_pts;
for (unsigned int i = 0; i < source_pts.n; i++) {
if (-1 != source_pts.vi[3*i+2]) {
+ tl_tmp->n++;
tl_tmp->vi[3*i] = source_pts.vi[3*i];
tl_tmp->vi[3*i+1] = source_pts.vi[3*i+1];
tl_tmp->vi[3*i+2] = source_pts.vi[3*i+2];
}
}
-
+
// no longer need source_pts
tuple_list_free(&source_pts);
// copy into tl if passed in and storing locally
if (tl && store_local) {
tuple_list_init_max(tl, 3, 0, 0, 1, num_pts);
- memcpy(tl->vi, tl_tmp->vi, 3*num_pts*sizeof(int));
+ memcpy(tl->vi, tl_tmp->vi, 3*tl_tmp->n*sizeof(int));
tl->n = tl_tmp->n;
}
-
+
// done
return MB_SUCCESS;
}
+MBErrorCode MBCoupler::test_local_box(double *xyz,
+ int from_proc, int remote_index, int index,
+ tuple_list *tl)
+{
+
+ std::vector<MBEntityHandle> entities;
+ std::vector<MBCartVect> nat_coords;
+
+ MBErrorCode result = nat_param(xyz, entities, nat_coords);
+ if (MB_SUCCESS != result) return result;
+
+ if (entities.empty()) return MB_SUCCESS;
+
+ // grow if we know we'll exceed size
+ if (mappedPts->n+entities.size() >= mappedPts->max)
+ tuple_list_grow(mappedPts);
+
+ if (entities.empty() && tl) {
+ tl->vi[3*index] = from_proc;
+ tl->vi[3*index+1] = remote_index;
+ tl->vi[3*index+2] = -1;
+ return MB_SUCCESS;
+ }
+
+ std::vector<MBEntityHandle>::iterator eit = entities.begin();
+ std::vector<MBCartVect>::iterator ncit = nat_coords.begin();
+ for (; eit != entities.end(); eit++, ncit++) {
+ // store in tuple mappedPts
+ mappedPts->n++;
+ mappedPts->vr[3*mappedPts->n] = (*ncit)[0];
+ mappedPts->vr[3*mappedPts->n+1] = (*ncit)[1];
+ mappedPts->vr[3*mappedPts->n+2] = (*ncit)[2];
+ mappedPts->vul[mappedPts->n] = *eit;
+
+ // also store local point, mapped point indices
+ if (tl)
+ {
+ if (tl->n == tl->max) tuple_list_grow(tl);
+
+ // store in tuple source_pts
+ tl->n++;
+ tl->vi[3*tl->n] = from_proc;
+ tl->vi[3*tl->n+1] = remote_index;
+ tl->vi[3*tl->n+2] = mappedPts->n;
+ }
+ else {
+ localMappedPts.push_back(index);
+ localMappedPts.push_back(mappedPts->n);
+ }
+ }
+
+ return MB_SUCCESS;
+}
+
MBErrorCode MBCoupler::interpolate(MBCoupler::Method method,
+ std::string &interp_tag,
+ double *interp_vals,
+ tuple_list *tl,
+ bool normalize)
+{
+ MBTag tag;
+ MBErrorCode result = mbImpl->tag_get_handle(interp_tag.c_str(), tag);
+ if (MB_SUCCESS != result) return result;
+ return interpolate(method, tag, interp_vals, tl, normalize);
+}
+
+MBErrorCode MBCoupler::interpolate(MBCoupler::Method method,
MBTag tag,
double *interp_vals,
tuple_list *tl,
- bool normalize)
+ bool normalize)
{
if (LINEAR_FE != method) return MB_FAILURE;
tuple_list *tl_tmp = (tl ? tl : targetPts);
+
+ // remote pts first
// scatter/gather interpolation points
gs_transfer(1, tl_tmp, 0, myPc->proc_config().crystal_router());
// perform interpolation on local source mesh; put results into
// tl_tmp->vr[i]
-
- // interpolate locally mapped points too, putting results directly
- // into interp_vals
-
- // normalize interpolation
-
+ MBErrorCode result;
+ for (unsigned int i = 0; i < tl_tmp->n; i++) {
+ int mindex = tl_tmp->vi[3*i+2];
+ result = interp_field_for_hex(mappedPts->vul[mindex],
+ MBCartVect(mappedPts->vr+3*mindex),
+ tag, tl_tmp->vr[i]);
+ if (MB_SUCCESS != result) return result;
+ }
+
// scatter/gather interpolation data
gs_transfer(1, tl_tmp, 0, myPc->proc_config().crystal_router());
- for (unsigned int i = 0; i < tl_tmp->n; i++)
- interp_vals[tl_tmp->vi[3*i+1]] = tl_tmp->vr[i];
+ if (!tl) {
+ // mapped whole targetPts tuple; put into proper place in interp_vals
+ for (unsigned int i = 0; i < tl_tmp->n; i++)
+ interp_vals[tl_tmp->vi[3*i+1]] = tl_tmp->vr[i];
+
+ // now do locally-contained pts, since we're mapping everything
+ for (std::vector<unsigned int>::iterator vit = localMappedPts.begin();
+ vit != localMappedPts.end(); vit += 2) {
+ int mindex = *(vit+1);
+ result = interp_field_for_hex(mappedPts->vul[mindex],
+ MBCartVect(mappedPts->vr+3*mindex),
+ tag, interp_vals[*vit]);
+ if (MB_SUCCESS != result) return result;
+ }
+ }
+
// done
return MB_SUCCESS;
}
-MBErrorCode MBCoupler::interpolate(MBCoupler::Method method,
- const char *tag_name,
- double *interp_vals,
- tuple_list *tl,
- bool normalize)
+MBErrorCode MBCoupler::nat_param(double xyz[3],
+ std::vector<MBEntityHandle> &entities,
+ std::vector<MBCartVect> &nat_coords)
{
- MBTag this_tag;
- MBErrorCode result = mbImpl->tag_get_handle(tag_name, this_tag);
- if (MB_SUCCESS != result) return result;
+ MBAdaptiveKDTreeIter treeiter;
+ MBErrorCode result = myTree->get_tree_iterator(localRoot, treeiter);
+ if (MB_SUCCESS != result) {
+ std::cout << "Problems getting iterator" << std::endl;
+ return result;
+ }
+
+ result = myTree->leaf_containing_point(localRoot, xyz, treeiter);
+ if (MB_SUCCESS != result) {
+ std::cout << "Problems getting leaf " << std::endl;
+ return result;
+ }
+
+ // find natural coordinates of point in element(s) in that leaf
+ MBCartVect tmp_nat_coords;
+ MBRange range_leaf;
+ result = mbImpl->get_entities_by_dimension(treeiter.handle(), 3, range_leaf, false);
+ if(result != MB_SUCCESS) std::cout << "Problem getting leaf in a range" << std::endl;
+
+ // loop over the range_leaf
+ for(MBRange::iterator iter = range_leaf.begin(); iter != range_leaf.end(); iter++)
+ {
+ const MBEntityHandle *connect;
+ int num_connect;
+
+ //get connectivity
+ result = mbImpl->get_connectivity(*iter, connect, num_connect, true);
+
+ //get coordinates of the vertices
+ std::vector<MBCartVect> coords_vert(num_connect);
+ std::vector<double> coords(3*num_connect);
+ result = mbImpl->get_coords(connect, num_connect, &coords[0]);
+
+ for(int j = 0; j < num_connect; j++)
+ {
+ coords_vert[j][0] = coords[3*j];
+ coords_vert[j][1] = coords[3*j+1];
+ coords_vert[j][2] = coords[3*j+2];
+ }
+
+ //test to find out in which hex the point is
+
+ // get natural coordinates
+ MBElemUtil::nat_coords_trilinear_hex(&coords_vert[0], MBCartVect(xyz),
+ tmp_nat_coords, 1e-10);
+ if (-1.0 <= tmp_nat_coords[0] && tmp_nat_coords[0] <= 1.0 &&
+ -1.0 <= tmp_nat_coords[1] && tmp_nat_coords[1] <= 1.0 &&
+ -1.0 <= tmp_nat_coords[2] && tmp_nat_coords[2] <= 1.0) {
+ entities.push_back(*iter);
+ nat_coords.push_back(tmp_nat_coords);
+ }
+ }
- return interpolate(method, this_tag, interp_vals, tl, normalize);
+ return MB_SUCCESS;
}
-MBErrorCode MBCoupler::locate_points(MBRange &ents,
- tuple_list *tl,
- bool store_local)
+MBErrorCode MBCoupler::interp_field_for_hex(MBEntityHandle elem,
+ MBCartVect nat_coord,
+ MBTag tag,
+ double &field)
{
- // only do vertices for now
- MBRange verts = ents.subset_by_type(MBVERTEX);
- if (verts.size() != ents.size()) return MB_FAILURE;
-
- // get coordinates
- std::vector<double> coords(3*verts.size());
- MBErrorCode result = mbImpl->get_coords(verts, &coords[0]);
+ //set the vertices coordinates in the natural system
+
+ const double xi[8] = {-1,1,1,-1,-1,1,1,-1};
+ const double etha[8] = {-1,-1,1,1,-1,-1,1,1};
+ const double mu[8] = {-1,-1,-1,-1,1,1,1,1};
+ double vfields[MB_MAX_SUB_ENTITIES*MB_MAX_SUB_ENTITY_VERTICES];
+
+ // get the tag values at the vertices
+ const MBEntityHandle *connect;
+ int num_connect;
+ MBErrorCode result = mbImpl->get_connectivity(elem, connect, num_connect);
if (MB_SUCCESS != result) return result;
+ result = mbImpl->tag_get_data(tag, connect, num_connect, vfields);
+ if (MB_SUCCESS != result) return result;
- return locate_points(&coords[0], verts.size(), tl, store_local);
-}
-
-MBErrorCode MBCoupler::initialize_tree()
-{
- // initialize the tree for the local mesh
+ //function for the interpolation
+ field = 0;
- // allocate box extents storage array
- allBoxes = new double[6*myPc->proc_config().proc_size()];
+ //calculate the field
+
+ // just hexes for now
+ assert(num_connect <= 8);
- // read local tree minimum/maximum vectors into array starting at
- // allBoxes[6*myPc->proc_rank()], allBoxes[6*myPc->proc_rank()+3]
+ for(int i = 0; i < num_connect; i++)
+ {
+ field += 0.125 * vfields[i] *
+ (1+xi[i]*nat_coord[0]) * (1+etha[i]*nat_coord[1]) * (1+mu[i]*nat_coord[2]);
+ }
- // now allgather so each proc has copy of box minima/maxima
- int mpi_err = MPI_Allgather(allBoxes+6*myPc->proc_config().proc_rank(), 6, MPI_DOUBLE,
- allBoxes, 6, MPI_DOUBLE, myPc->proc_config().proc_comm());
- if (MPI_SUCCESS == mpi_err) return MB_SUCCESS;
- else return MB_FAILURE;
+ return MB_SUCCESS;
}
-
-
Modified: MOAB/trunk/tools/mbcoupler/MBCoupler.hpp
===================================================================
--- MOAB/trunk/tools/mbcoupler/MBCoupler.hpp 2008-06-26 00:34:42 UTC (rev 1930)
+++ MOAB/trunk/tools/mbcoupler/MBCoupler.hpp 2008-06-26 14:55:10 UTC (rev 1931)
@@ -28,11 +28,15 @@
#include "MBRange.hpp"
#include "MBInterface.hpp"
+#include "MBCartVect.hpp"
+
extern "C"
{
struct tuple_list;
}
+class MBAdaptiveKDTree;
+
class MBCoupler
{
public:
@@ -119,13 +123,28 @@
* \param normalize If true, normalization is done according to method
*/
MBErrorCode interpolate(MBCoupler::Method method,
- const char* tag_name,
+ std::string &tag_name,
double *interp_vals,
tuple_list *tl = NULL,
bool normalize = true);
private:
+ // given a coordinate position, find all entities containing
+ // the point and the natural coords in those ents
+ MBErrorCode nat_param(double xyz[3],
+ std::vector<MBEntityHandle> &entities,
+ std::vector<MBCartVect> &nat_coords);
+
+ MBErrorCode interp_field_for_hex(MBEntityHandle elem,
+ MBCartVect nat_coord,
+ MBTag tag,
+ double &field);
+
+ MBErrorCode test_local_box(double *xyz,
+ int from_proc, int remote_index, int index,
+ tuple_list *tl = NULL);
+
/* \brief MOAB instance
*/
MBInterface *mbImpl;
@@ -134,9 +153,17 @@
*/
MBErrorCode initialize_tree();
+ /* \brief Kdtree for local mesh
+ */
+ MBAdaptiveKDTree *myTree;
+
+ /* \brief Local root of the kdtree
+ */
+ MBEntityHandle localRoot;
+
/* \brief Min/max bounding boxes for all proc tree roots
*/
- double *allBoxes;
+ std::vector<double> allBoxes;
/* \brief MBParallelComm object for this coupler
*/
Modified: MOAB/trunk/tools/mbcoupler/MBElemUtil.cpp
===================================================================
--- MOAB/trunk/tools/mbcoupler/MBElemUtil.cpp 2008-06-26 00:34:42 UTC (rev 1930)
+++ MOAB/trunk/tools/mbcoupler/MBElemUtil.cpp 2008-06-26 14:55:10 UTC (rev 1931)
@@ -12,12 +12,16 @@
// coordinates {0, 0, 0} and calculating new coordinates based on the
// ones just calculated.
//
-void MBElemUtil::nat_coords_trilinear_hex(MBCartVect hex[8],
+void MBElemUtil::nat_coords_trilinear_hex(MBCartVect *hex,
MBCartVect xyz,
MBCartVect &ncoords,
double etol)
{
+ // short-circuit for now...
+ ncoords[0] = ncoords[1] = ncoords[2] = 0.5;
+ return;
+
MBCartVect nat(0.);
double A[8]; double B[8]; double C[8]; double D[8];
double Ax, By, Cz;
@@ -27,7 +31,7 @@
double nxi, neta, nmu;
double pxi, peta, pmu;
double tmp;
-
+
// Iterative estimate of natural coordinates
while ( err > etol ) {
@@ -45,9 +49,10 @@
for (unsigned int i = 0; i < 8; i++) {
Ax = Ax + A[i]*hex[i][0];
}
- nat[0] = (8*xyz[0] - Ax ) /
- (-A[0]*hex[0][0] + A[1]*hex[1][0] + A[2]*hex[2][0] - A[3]*hex[3][0]
- -A[4]*hex[4][0] + A[5]*hex[5][0] + A[6]*hex[6][0] - A[7]*hex[7][0]);
+ double tmp_d = -A[0]*hex[0][0] + A[1]*hex[1][0] + A[2]*hex[2][0] - A[3]*hex[3][0]
+ -A[4]*hex[4][0] + A[5]*hex[5][0] + A[6]*hex[6][0] - A[7]*hex[7][0];
+ if (0.0 == tmp_d) nat[0] = 2.0;
+ else nat[0] = (8*xyz[0] - Ax ) / tmp_d;
// Estimate the eta-coordinate
B[0] = (1 - nat[0]) * (1 - nat[2]);
@@ -63,9 +68,11 @@
for (unsigned int i = 0; i < 8; i++) {
By = By + B[i]*hex[i][1];
}
- nat[1] = (8*xyz[1] - By ) /
- (-B[0]*hex[0][1] - B[1]*hex[1][1] + B[2]*hex[2][1] + B[3]*hex[3][1]
- -B[4]*hex[4][1] - B[5]*hex[5][1] + B[6]*hex[6][1] + B[7]*hex[7][1]);
+ tmp_d = -B[0]*hex[0][1] - B[1]*hex[1][1] + B[2]*hex[2][1] + B[3]*hex[3][1]
+ -B[4]*hex[4][1] - B[5]*hex[5][1] + B[6]*hex[6][1] + B[7]*hex[7][1];
+
+ if (0.0 == tmp_d) nat[1] = 2.0;
+ else nat[1] = (8*xyz[1] - By ) / tmp_d;
// Estimate the mu-coordinate
C[0] = (1 - nat[0]) * (1 - nat[1]);
@@ -81,10 +88,11 @@
for (unsigned int i = 0; i < 8; i++) {
Cz = Cz + C[i]*hex[i][2];
}
- nat[2] = (8*xyz[2] - Cz ) /
- (-C[0]*hex[0][2] - C[1]*hex[1][2] - C[2]*hex[2][2] - C[3]*hex[3][2]
- +C[4]*hex[4][2] + C[5]*hex[5][2] + C[6]*hex[6][2] + C[7]*hex[7][2]);
+ tmp_d = -C[0]*hex[0][2] - C[1]*hex[1][2] - C[2]*hex[2][2] - C[3]*hex[3][2]
+ +C[4]*hex[4][2] + C[5]*hex[5][2] + C[6]*hex[6][2] + C[7]*hex[7][2];
+ if (0.0 == tmp_d) nat[2] = 2.0;
+ else nat[2] = (8*xyz[2] - Cz ) / tmp_d;
// Shortcut variables...
nxi = 1 - nat[0];
@@ -132,7 +140,7 @@
}
-bool MBElemUtil::point_in_trilinear_hex(MBCartVect hex[8],
+bool MBElemUtil::point_in_trilinear_hex(MBCartVect *hex,
MBCartVect xyz,
double etol)
{
@@ -151,7 +159,7 @@
}
-bool MBElemUtil::point_in_trilinear_hex(MBCartVect hex[8],
+bool MBElemUtil::point_in_trilinear_hex(MBCartVect *hex,
MBCartVect xyz,
MBCartVect box_min, MBCartVect box_max,
double etol)
Modified: MOAB/trunk/tools/mbcoupler/MBElemUtil.hpp
===================================================================
--- MOAB/trunk/tools/mbcoupler/MBElemUtil.hpp 2008-06-26 00:34:42 UTC (rev 1930)
+++ MOAB/trunk/tools/mbcoupler/MBElemUtil.hpp 2008-06-26 14:55:10 UTC (rev 1931)
@@ -5,15 +5,15 @@
public:
- static void nat_coords_trilinear_hex(MBCartVect[8],
+ static void nat_coords_trilinear_hex(MBCartVect*,
MBCartVect,
MBCartVect&,
double);
- static bool point_in_trilinear_hex(MBCartVect hex[8],
+ static bool point_in_trilinear_hex(MBCartVect *hex,
MBCartVect xyz,
double etol);
- static bool point_in_trilinear_hex(MBCartVect hex[8],
+ static bool point_in_trilinear_hex(MBCartVect *hex,
MBCartVect xyz,
MBCartVect box_min,
MBCartVect box_max,
Modified: MOAB/trunk/tools/mbcoupler/mbcoupler_test.cpp
===================================================================
--- MOAB/trunk/tools/mbcoupler/mbcoupler_test.cpp 2008-06-26 00:34:42 UTC (rev 1930)
+++ MOAB/trunk/tools/mbcoupler/mbcoupler_test.cpp 2008-06-26 14:55:10 UTC (rev 1931)
@@ -1,7 +1,9 @@
#include "MBParallelComm.hpp"
+#include "MBParallelConventions.h"
#include "MBCore.hpp"
#include "FileOptions.hpp"
#include "ReadParallel.hpp"
+#include "MBCoupler.hpp"
#include "mpi.h"
#include <iostream>
#include <sstream>
@@ -13,13 +15,27 @@
dynamic_cast<MBCore*>(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;\
+ std::cout << mbImpl->get_last_error(tmp_str) << std::endl;\
+ return result;\
+ }
+
MBErrorCode get_file_options(int argc, char **argv,
std::vector<const char *> &filenames,
+ std::string &tag_name,
std::string &opts);
MBErrorCode report_iface_ents(MBInterface *mbImpl,
std::vector<MBParallelComm *> &pcs);
+MBErrorCode test_interpolation(MBInterface *mbImpl,
+ std::string &interp_tag,
+ std::vector<MBParallelComm *> &pcs,
+ std::vector<ReadParallel *> rps);
+
int main(int argc, char **argv)
{
// need to init MPI first, to tell how many procs and rank
@@ -27,9 +43,10 @@
if (argc < 3) {
std::cerr << "Usage: ";
- std::cerr << argv[0] << " <nfiles> <fname1> ... <fnamen> [tag_name] [tag_val] [distrib] [with_ghosts]" << std::endl;
+ std::cerr << argv[0] << " <nfiles> <fname1> ... <fnamen> [interp_tag] [tag_name] [tag_val] [distrib] [with_ghosts]" << std::endl;
std::cerr << "nfiles : number of mesh files" << std::endl;
std::cerr << "fname1..fnamen: mesh files" << std::endl;
+ std::cerr << "interp_tag : name of tag interpolated to target mesh []" << std::endl;
std::cerr << "tag_name : name of tag used to define partitions [MATERIAL_SET]" << std::endl;
std::cerr << "tag_val : tag values denoting partition sets [--]" << std::endl;
std::cerr << "distrib : if non-zero, distribute the partition sets with tag_val round-robin" << std::endl;
@@ -55,8 +72,8 @@
MBErrorCode result = MB_SUCCESS;
std::vector<const char *> filenames;
- std::string opts;
- result = get_file_options(argc, argv, filenames, opts);
+ std::string opts, interp_tag;
+ result = get_file_options(argc, argv, filenames, interp_tag, opts);
// read in mesh(es)
@@ -70,19 +87,22 @@
result = rps[i]->load_file(filenames[i], filesets[i],
FileOptions(opts.c_str()), NULL, 0);
- if (MB_SUCCESS != result) {
- std::string tmp_str;
- std::cout << "Failure; message:" << std::endl;
- std::cout << mbImpl->get_last_error(tmp_str) << std::endl;
- return 1;
- }
+ PRINT_LAST_ERROR;
}
result = report_iface_ents(mbImpl, pcs);
- if (MB_SUCCESS != result)
- std::cout << "Failure." << std::endl;
- else
- std::cout << "Success." << std::endl;
+ PRINT_LAST_ERROR;
+
+ // test interpolation
+ result = test_interpolation(mbImpl, interp_tag, pcs, rps);
+ PRINT_LAST_ERROR;
+
+ // output mesh
+ result = mbImpl->write_file("output.h5m", NULL, NULL,
+ pcs[1]->partition_sets());
+ PRINT_LAST_ERROR;
+
+ std::cout << "Success." << std::endl;
err = MPI_Finalize();
for (unsigned int i = 0; i < filenames.size(); i++) {
@@ -138,6 +158,7 @@
MBErrorCode get_file_options(int argc, char **argv,
std::vector<const char *> &filenames,
+ std::string &interp_tag,
std::string &opts)
{
int npos = 1;
@@ -146,6 +167,8 @@
// get mesh filenames
filenames.resize(nfiles);
for (int i = 0; i < nfiles; i++) filenames[i] = argv[npos++];
+
+ if (npos < argc) interp_tag = argv[npos++];
// get partition information
const char *tag_name = "MATERIAL_SET";
@@ -177,3 +200,55 @@
return MB_SUCCESS;
}
+
+MBErrorCode test_interpolation(MBInterface *mbImpl,
+ std::string &interp_tag,
+ std::vector<MBParallelComm *> &pcs,
+ std::vector<ReadParallel *> rps)
+{
+ // source is 1st mesh, target is 2nd
+ MBRange src_elems, targ_elems;
+ MBErrorCode result = pcs[0]->get_part_entities(src_elems, 3);
+ PRINT_LAST_ERROR;
+
+ // instantiate a coupler, which also initializes the tree
+ MBCoupler mbc(mbImpl, pcs[0], src_elems, 0);
+
+ // get points from the target mesh to interpolate
+ MBRange targ_verts, tmp_verts;
+
+ // first get all vertices adj to partition entities in target mesh
+ result = pcs[1]->get_part_entities(targ_elems, 3);
+ result = mbImpl->get_adjacencies(targ_elems, 0, false, targ_verts,
+ MBInterface::UNION);
+ PRINT_LAST_ERROR;
+
+ // then get non-owned verts and subtract
+ result = pcs[1]->get_pstatus_entities(0, PSTATUS_NOT_OWNED, tmp_verts);
+ PRINT_LAST_ERROR;
+ targ_verts = targ_verts.subtract(tmp_verts);
+
+ // get position of these entities; these are the target points
+ std::vector<double> vpos(3*targ_verts.size());
+ result = mbImpl->get_coords(targ_verts, &vpos[0]);
+ PRINT_LAST_ERROR;
+
+ // locate those points in the source mesh
+ result = mbc.locate_points(&vpos[0], targ_verts.size());
+ PRINT_LAST_ERROR;
+
+ // now interpolate tag onto target points
+ std::vector<double> field(targ_verts.size());
+ result = mbc.interpolate(MBCoupler::LINEAR_FE, interp_tag, &field[0]);
+ PRINT_LAST_ERROR;
+
+ // set field values as tag on target vertices
+ MBTag tag;
+ result = mbImpl->tag_get_handle(interp_tag.c_str(), tag); PRINT_LAST_ERROR;
+ result = mbImpl->tag_set_data(tag, targ_verts, &field[0]); PRINT_LAST_ERROR;
+
+ int err = MPI_Finalize();
+
+ // done
+ return MB_SUCCESS;
+}
More information about the moab-dev
mailing list