46 #ifndef _ZOLTAN2_COORDCOMMGRAPH_HPP_ 47 #define _ZOLTAN2_COORDCOMMGRAPH_HPP_ 56 #include "Teuchos_CommHelpers.hpp" 57 #include "Teuchos_Comm.hpp" 58 #include "Teuchos_ArrayViewDecl.hpp" 59 #include "Teuchos_RCPDecl.hpp" 64 #define Z2_ABS(x) ((x) >= 0 ? (x) : -(x)) 69 template <
typename scalar_t,
typename part_t>
83 part_t *minHashIndices;
84 part_t *maxHashIndices;
87 std::vector <part_t> *gridIndices;
89 std::set <part_t> neighbors;
97 maxScalar (
std::numeric_limits<scalar_t>::max()),
98 _EPSILON(
std::numeric_limits<scalar_t>::
epsilon()),
101 gridIndices(0), neighbors()
103 lmins =
new scalar_t [dim];
104 lmaxs =
new scalar_t [dim];
106 minHashIndices =
new part_t [dim];
107 maxHashIndices =
new part_t [dim];
108 gridIndices =
new std::vector <part_t> ();
109 for (
int i = 0; i < dim; ++i){
110 lmins[i] = -this->maxScalar;
111 lmaxs[i] = this->maxScalar;
121 maxScalar (
std::numeric_limits<scalar_t>::max()),
122 _EPSILON(
std::numeric_limits<scalar_t>::
epsilon()),
125 gridIndices(0), neighbors()
127 lmins =
new scalar_t [dim];
128 lmaxs =
new scalar_t [dim];
129 minHashIndices =
new part_t [dim];
130 maxHashIndices =
new part_t [dim];
131 gridIndices =
new std::vector <part_t> ();
132 for (
int i = 0; i < dim; ++i){
146 maxScalar (
std::numeric_limits<scalar_t>::max()),
147 _EPSILON(
std::numeric_limits<scalar_t>::
epsilon()),
150 gridIndices(0), neighbors()
153 lmins =
new scalar_t [dim];
154 lmaxs =
new scalar_t [dim];
155 minHashIndices =
new part_t [dim];
156 maxHashIndices =
new part_t [dim];
157 gridIndices =
new std::vector <part_t> ();
158 scalar_t *othermins = other.
getlmins();
159 scalar_t *othermaxs = other.
getlmaxs();
160 for (
int i = 0; i < dim; ++i){
161 lmins[i] = othermins[i];
162 lmaxs[i] = othermaxs[i];
168 delete []this->lmins;
169 delete [] this->lmaxs;
170 delete []this->minHashIndices;
171 delete [] this->maxHashIndices;
205 for (
int i = 0; i < this->dim; i++)
206 centroid[i] = 0.5 * (this->lmaxs[i] + this->lmins[i]);
213 return this->gridIndices;
219 return &(this->neighbors);
225 if (pointdim != this->dim)
226 throw std::logic_error(
"dim of point must match dim of box");
227 for (
int i = 0; i < pointdim; i++) {
228 if (point[i] < this->lmins[i])
return false;
229 if (point[i] > this->lmaxs[i])
return false;
237 if (cdim != this->dim)
238 throw std::logic_error(
"dim of given box must match dim of box");
242 for (
int i = 0; i < cdim; i++) {
243 if (!((lower[i] >= this->lmins[i] && lower[i] <= this->lmaxs[i])
245 || (upper[i] >= this->lmins[i] && upper[i] <= this->lmaxs[i])
247 || (lower[i] < this->lmins[i] && upper[i] > this->lmaxs[i]))) {
266 for (
int i = 0; i < dim; ++i){
268 if (omins[i] - this->lmaxs[i] > _EPSILON ||
269 this->lmins[i] - omaxs[i] > _EPSILON ) {
272 else if (
Z2_ABS(omins[i] - this->lmaxs[i]) < _EPSILON ||
273 Z2_ABS(this->lmins[i] - omaxs[i]) < _EPSILON ){
283 std::cout <<
"something is wrong: equality:" 284 << equality << std::endl;
293 neighbors.insert(nIndex);
299 if (neighbors.end() != neighbors.find(nIndex)){
310 scalar_t *minMaxBoundaries,
311 scalar_t *sliceSizes,
312 part_t numSlicePerDim
314 for (
int j = 0; j < dim; ++j){
315 scalar_t distance = (lmins[j] - minMaxBoundaries[j]);
317 if (distance > _EPSILON && sliceSizes[j] > _EPSILON){
318 minInd =
static_cast<part_t
>(floor((lmins[j] - minMaxBoundaries[j])/ sliceSizes[j]));
321 if(minInd >= numSlicePerDim){
322 minInd = numSlicePerDim - 1;
325 distance = (lmaxs[j] - minMaxBoundaries[j]);
326 if (distance > _EPSILON && sliceSizes[j] > _EPSILON){
327 maxInd =
static_cast<part_t
>(ceil((lmaxs[j] - minMaxBoundaries[j])/ sliceSizes[j]));
329 if(maxInd >= numSlicePerDim){
330 maxInd = numSlicePerDim - 1;
335 minHashIndices[j] = minInd;
336 maxHashIndices[j] = maxInd;
338 std::vector <part_t> *in =
new std::vector <part_t> ();
340 std::vector <part_t> *out =
new std::vector <part_t> ();
342 for (
int j = 0; j < dim; ++j){
344 part_t minInd = minHashIndices[j];
345 part_t maxInd = maxHashIndices[j];
348 part_t pScale = part_t(pow (
float(numSlicePerDim),
int(dim - j -1)));
350 part_t inSize = in->size();
352 for (part_t k = minInd; k <= maxInd; ++k){
353 for (part_t i = 0; i < inSize; ++i){
354 out->push_back((*in)[i] + k * pScale);
358 std::vector <part_t> *tmp = in;
363 std::vector <part_t> *tmp = in;
375 for(
int i = 0; i < this->dim; ++i){
376 std::cout <<
"\tbox:" << this->pID <<
" dim:" << i <<
" min:" << lmins[i] <<
" max:" << lmaxs[i] << std::endl;
384 lmaxs[dimInd] = newBoundary;
387 lmins[dimInd] = newBoundary;
394 int numCorners = (int(1)<<dim);
395 scalar_t *corner1 =
new scalar_t [dim];
396 scalar_t *corner2 =
new scalar_t [dim];
398 for (
int i = 0; i < dim; ++i){
410 mm << lmins[i] <<
" ";
414 for (
int i = 0; i < dim; ++i){
428 mm << lmaxs[i] <<
" ";
433 for (
int j = 0; j < numCorners; ++j){
434 std::vector <int> neighborCorners;
435 for (
int i = 0; i < dim; ++i){
436 if(
int(j & (
int(1)<<i)) == 0){
437 corner1[i] = lmins[i];
440 corner1[i] = lmaxs[i];
442 if (j % (
int(1)<<(i + 1)) >= (int(1)<<i)){
443 int c1 = j - (int(1)<<i);
446 neighborCorners.push_back(c1);
451 int c1 = j + (int(1)<<i);
452 if (c1 < (
int(1) << dim)) {
453 neighborCorners.push_back(c1);
458 for (
int m = 0; m < int (neighborCorners.size()); ++m){
460 int n = neighborCorners[m];
462 for (
int i = 0; i < dim; ++i){
463 if(
int(n & (
int(1)<<i)) == 0){
464 corner2[i] = lmins[i];
467 corner2[i] = lmaxs[i];
471 std::string arrowline =
"set arrow from ";
472 for (
int i = 0; i < dim - 1; ++i){
474 Teuchos::toString<scalar_t>(corner1[i]) +
",";
477 Teuchos::toString<scalar_t>(corner1[dim -1]) +
" to ";
479 for (
int i = 0; i < dim - 1; ++i){
481 Teuchos::toString<scalar_t>(corner2[i]) +
",";
484 Teuchos::toString<scalar_t>(corner2[dim -1]) +
501 template <
typename scalar_t,
typename part_t>
505 const RCP < std::vector <Zoltan2::coordinateModelPartBox <scalar_t, part_t> > > pBoxes;
508 scalar_t *minMaxBoundaries;
510 scalar_t *maxMinBoundaries;
512 scalar_t *sliceSizes;
516 part_t numSlicePerDim;
520 std::vector <std::vector <part_t> > grids;
522 ArrayRCP <part_t> comXAdj;
523 ArrayRCP <part_t> comAdj;
530 part_t ntasks_,
int dim_):
533 maxMinBoundaries(0), sliceSizes(0),
536 numSlicePerDim(part_t(pow(double(ntasks_), 1.0 / dim))),
542 minMaxBoundaries =
new scalar_t[dim];
543 maxMinBoundaries =
new scalar_t[dim];
544 sliceSizes =
new scalar_t[dim];
547 if (numSlicePerDim == 0) numSlicePerDim = 1;
549 numGrids = part_t(pow(
float(numSlicePerDim),
int(dim)));
552 std::vector <std::vector <part_t> > grids_ (numGrids);
553 this->grids = grids_;
562 ArrayRCP <part_t> tmpComXadj(ntasks_+1);
563 ArrayRCP <part_t> tmpComAdj(nCount);
564 comXAdj = tmpComXadj;
575 delete []minMaxBoundaries;
576 delete []maxMinBoundaries;
588 for(part_t i = 0; i < this->nTasks; ++i){
589 std::set<part_t> *neigbors = (*pBoxes)[i].getNeighbors();
591 part_t s = neigbors->size();
593 comXAdj[i+1] = comXAdj[i] + s;
594 typedef typename std::set<part_t> mySet;
595 typedef typename mySet::iterator myIT;
597 for (it=neigbors->begin(); it!=neigbors->end(); ++it)
599 comAdj[adjIndex++] = *it;
611 ArrayRCP <part_t> &comXAdj_,
612 ArrayRCP <part_t> &comAdj_){
613 comXAdj_ = this->comXAdj;
614 comAdj_ = this->comAdj;
622 for(part_t i = 0; i < this->nTasks; ++i){
623 std::vector <part_t> *gridIndices =(*pBoxes)[i].getGridIndices();
624 part_t gridCount = gridIndices->size();
626 for (part_t j = 0; j < gridCount; ++j){
627 part_t grid = (*gridIndices)[j];
628 part_t boxCount = grids[grid].size();
629 for (part_t k = 0; k < boxCount; ++k){
630 part_t boxIndex = grids[grid][k];
632 if((!(*pBoxes)[i].isAlreadyNeighbor(boxIndex))&& (*pBoxes)[i].isNeighborWith((*pBoxes)[boxIndex])){
634 (*pBoxes)[i].addNeighbor(boxIndex);
635 (*pBoxes)[boxIndex].addNeighbor(i);
652 for(part_t i = 0; i < this->nTasks; ++i){
653 (*pBoxes)[i].setMinMaxHashIndices(minMaxBoundaries, sliceSizes, numSlicePerDim);
655 std::vector <part_t> *gridIndices =(*pBoxes)[i].getGridIndices();
657 part_t gridCount = gridIndices->size();
659 for (part_t j = 0; j < gridCount; ++j){
660 part_t grid = (*gridIndices)[j];
663 (grids)[grid].push_back(i);
684 scalar_t *mins = (*pBoxes)[0].getlmins();
685 scalar_t *maxs = (*pBoxes)[0].getlmaxs();
687 for (
int j = 0; j < dim; ++j){
688 minMaxBoundaries[j] = maxs[j];
689 maxMinBoundaries[j] = mins[j];
692 for (part_t i = 1; i < nTasks; ++i){
694 mins = (*pBoxes)[i].getlmins();
695 maxs = (*pBoxes)[i].getlmaxs();
697 for (
int j = 0; j < dim; ++j){
699 if (minMaxBoundaries[j] > maxs[j]){
700 minMaxBoundaries[j] = maxs[j];
702 if (maxMinBoundaries[j] < mins[j]){
703 maxMinBoundaries[j] = mins[j];
709 for (
int j = 0; j < dim; ++j){
710 sliceSizes[j] = (maxMinBoundaries[j] - minMaxBoundaries[j]) / numSlicePerDim;
711 if (sliceSizes[j] < 0) sliceSizes[j] = 0;
int getDim() const
function to set the dimension
GridHash Class, Hashing Class for part boxes.
std::set< part_t > * getNeighbors()
function to get the indices of the neighboring parts.
coordinateModelPartBox(part_t pid, int dim_)
Constructor.
~GridHash()
GridHash Class, Destructor.
void setpId(part_t pid)
function to set the part id
void getAdjArrays(ArrayRCP< part_t > &comXAdj_, ArrayRCP< part_t > &comAdj_)
GridHash Class, returns the adj arrays.
void fillAdjArrays()
GridHash Class, Function to fill adj arrays.
bool boxesOverlap(int cdim, scalar_t *lower, scalar_t *upper) const
function to test whether this box overlaps a given box
std::vector< part_t > * getGridIndices()
function to get the indices of the buckets that the part is inserted to
coordinateModelPartBox Class, represents the boundaries of the box which is a result of a geometric p...
scalar_t * getlmaxs() const
function to get maximum values along all dimensions
bool pointInBox(int pointdim, scalar_t *point) const
function to test whether a point is in the box
void computeCentroid(scalar_t *centroid) const
compute the centroid of the box
bool isNeighborWith(const coordinateModelPartBox< scalar_t, part_t > &other) const
function to check if two boxes are neighbors.
void writeGnuPlot(std::ofstream &file, std::ofstream &mm)
function for visualization.
void print()
function to print the boundaries.
void updateMinMax(scalar_t newBoundary, int isMax, int dimInd)
function to update the boundary of the box.
GridHash(const RCP< std::vector< Zoltan2::coordinateModelPartBox< scalar_t, part_t > > > &pBoxes_, part_t ntasks_, int dim_)
GridHash Class, Constructor.
part_t getpId() const
function to get the part id
part_t calculateNeighbors()
GridHash Class, For each box compares the adjacency against the boxes that are in the same buckets...
void getMinMaxBoundaries()
GridHash Class, calculates the minimum of maximum box boundaries, and maxium of minimum box boundarie...
~coordinateModelPartBox()
Destructor.
void setMinMaxHashIndices(scalar_t *minMaxBoundaries, scalar_t *sliceSizes, part_t numSlicePerDim)
function to obtain the min and max hash values along all dimensions.
bool isAlreadyNeighbor(part_t nIndex)
function to check if a given part is already in the neighbor list.
void insertToHash()
GridHash Class, For each box calculates the buckets which it should be inserted to.
coordinateModelPartBox(const coordinateModelPartBox< scalar_t, part_t > &other)
Copy Constructor deep copy of the maximum and minimum boundaries.
void addNeighbor(part_t nIndex)
function to add a new neighbor to the neighbor list.
scalar_t * getlmins() const
function to get minimum values along all dimensions
coordinateModelPartBox(part_t pid, int dim_, scalar_t *lmi, scalar_t *lma)
Constructor deep copy of the maximum and minimum boundaries.