| #ifndef _box_utils_hpp_ |
| #define _box_utils_hpp_ |
| |
| //@HEADER |
| // ************************************************************************ |
| // |
| // MiniFE: Simple Finite Element Assembly and Solve |
| // Copyright (2006-2013) Sandia Corporation |
| // |
| // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive |
| // license for use of this work by or on behalf of the U.S. Government. |
| // |
| // This library is free software; you can redistribute it and/or modify |
| // it under the terms of the GNU Lesser General Public License as |
| // published by the Free Software Foundation; either version 2.1 of the |
| // License, or (at your option) any later version. |
| // |
| // This library is distributed in the hope that it will be useful, but |
| // WITHOUT ANY WARRANTY; without even the implied warranty of |
| // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU |
| // Lesser General Public License for more details. |
| // |
| // You should have received a copy of the GNU Lesser General Public |
| // License along with this library; if not, write to the Free Software |
| // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 |
| // USA |
| // |
| // ************************************************************************ |
| //@HEADER |
| |
| #include <vector> |
| #include <set> |
| #include <map> |
| |
| #ifdef HAVE_MPI |
| #include <mpi.h> |
| #endif |
| |
| #include <TypeTraits.hpp> |
| #include <Box.hpp> |
| |
| namespace miniFE { |
| |
| inline void copy_box(const Box& from_box, Box& to_box) |
| { |
| for(int i=0; i<3; ++i) { |
| to_box[i][0] = from_box[i][0]; |
| to_box[i][1] = from_box[i][1]; |
| } |
| } |
| |
| template<typename GlobalOrdinal> |
| #ifdef __CUDACC__ |
| __host__ __device__ __inline__ |
| #endif |
| void get_int_coords(GlobalOrdinal ID, int nx, int ny, int nz, |
| int& x, int& y, int& z) |
| { |
| z = ID/(nx*ny); |
| y = (ID%(nx*ny))/nx; |
| x = ID%nx; |
| } |
| |
| template<typename GlobalOrdinal,typename Scalar> |
| #ifdef __CUDACC__ |
| __host__ __device__ __inline__ |
| #endif |
| void get_coords(GlobalOrdinal ID, int nx, int ny, int nz, |
| Scalar& x, Scalar& y, Scalar& z) |
| { |
| const int xdiv = nx>1 ? nx-1 : 1; |
| const int ydiv = ny>1 ? ny-1 : 1; |
| const int zdiv = nz>1 ? nz-1 : 1; |
| |
| //This code assumes that ID is 0-based. |
| // |
| //compute coordinates that lie on (or in) the unit cube. |
| //that's why we're dividing by nz,ny,nx: |
| z = (1.0*(ID/(nx*ny)))/zdiv; |
| y = 1.0*((ID%(nx*ny))/nx)/ydiv; |
| x = 1.0*(ID%nx)/xdiv; |
| } |
| |
| template<typename GlobalOrdinal> |
| GlobalOrdinal get_num_ids(const Box& box) |
| { |
| int nx = box[0][1] - box[0][0]; |
| int ny = box[1][1] - box[1][0]; |
| int nz = box[2][1] - box[2][0]; |
| GlobalOrdinal tmp = nx*ny; |
| tmp *= nz; |
| return tmp; |
| } |
| |
| template<typename GlobalOrdinal> |
| #ifdef __CUDACC__ |
| __host__ __device__ __inline__ |
| #endif |
| GlobalOrdinal get_id(int nx, int ny, int nz, |
| int x, int y, int z) |
| { |
| if (x<0 || y<0 || z<0) return -1; |
| if (x>=nx || y>=ny || z>=nz) return -1; |
| |
| //form x + nx*y + nx*ny*z: |
| |
| GlobalOrdinal tmp = nx*ny; |
| tmp *= z; |
| tmp = x + nx * y + tmp; |
| return tmp; |
| } |
| |
| template<typename GlobalOrdinal> |
| void get_ids(int nx, int ny, int nz, |
| const Box& box, |
| std::vector<GlobalOrdinal>& ids, |
| bool include_ghost_layer=false) |
| { |
| ids.clear(); |
| int minz = box[2][0]; |
| int maxz = box[2][1]; |
| int miny = box[1][0]; |
| int maxy = box[1][1]; |
| int minx = box[0][0]; |
| int maxx = box[0][1]; |
| |
| if (include_ghost_layer) { |
| if (minz > 0) minz--; |
| if (miny > 0) miny--; |
| if (minx > 0) minx--; |
| if (maxz < nz) maxz++; |
| if (maxy < ny) maxy++; |
| if (maxx < nx) maxx++; |
| } |
| |
| size_t ids_size = ((maxz - minz) * (maxy - miny)) * (maxx - minx); |
| ids.reserve(ids_size); |
| |
| for(int z=minz; z<maxz; ++z) { |
| for(int y=miny; y<maxy; ++y) { |
| for(int x=minx; x<maxx; ++x) { |
| ids.push_back(get_id<GlobalOrdinal>(nx, ny, nz, x, y, z)); |
| } |
| } |
| } |
| } |
| |
| template<typename GlobalOrdinal> |
| void get_ghost_ids(int nx, int ny, int nz, |
| const Box& box, |
| std::vector<GlobalOrdinal>& ids) |
| { |
| ids.clear(); |
| int minz,maxz,miny,maxy,minx,maxx; |
| int orig_minz = minz = box[2][0]; |
| int orig_maxz = maxz = box[2][1]; |
| int orig_miny = miny = box[1][0]; |
| int orig_maxy = maxy = box[1][1]; |
| int orig_minx = minx = box[0][0]; |
| int orig_maxx = maxx = box[0][1]; |
| |
| if (minz > 0) minz--; |
| if (miny > 0) miny--; |
| if (minx > 0) minx--; |
| if (maxz < nz) maxz++; |
| if (maxy < ny) maxy++; |
| if (maxx < nx) maxx++; |
| |
| for(int z=minz; z<maxz; ++z) { |
| for(int y=miny; y<maxy; ++y) { |
| for(int x=minx; x<maxx; ++x) { |
| bool x_in_ghost_layer = (x < orig_minx) || (x >= orig_maxx); |
| bool y_in_ghost_layer = (y < orig_miny) || (y >= orig_maxy); |
| bool z_in_ghost_layer = (z < orig_minz) || (z >= orig_maxz); |
| //we are in the ghost layer if any one of x,y,z are in the ghost layer |
| if (!x_in_ghost_layer && !y_in_ghost_layer && !z_in_ghost_layer) continue; |
| ids.push_back(get_id<GlobalOrdinal>(nx, ny, nz, x, y, z)); |
| } |
| } |
| } |
| } |
| |
| inline void print_box(int myproc, const char* name, const Box& box, |
| const char* name2, const Box& box2) |
| { |
| std::cout << "proc " << myproc << " "<<name |
| <<" ("<<box[0][0]<<","<<box[0][1]<<") " |
| <<" ("<<box[1][0]<<","<<box[1][1]<<") " |
| <<" ("<<box[2][0]<<","<<box[2][1]<<") " |
| <<name2 |
| <<" ("<<box2[0][0]<<","<<box2[0][1]<<") " |
| <<" ("<<box2[1][0]<<","<<box2[1][1]<<") " |
| <<" ("<<box2[2][0]<<","<<box2[2][1]<<") "<<std::endl; |
| } |
| |
| bool is_neighbor(const Box& box1, const Box& box2) |
| { |
| //neighbors in the x dimension if: |
| bool x_neighbor = (box1[0][1] == box2[0][0]) || (box1[0][0] == box2[0][1]) || // min matches max |
| (box1[0][0] == box2[0][0]) || (box1[0][1] == box2[0][1]) || // mins or maxs match |
| (box1[0][0] > box2[0][0] && box1[0][1] < box2[0][1]) || // range contains other |
| (box2[0][0] > box1[0][0] && box2[0][1] < box1[0][1]) || // range contains other |
| (box1[0][0] > box2[0][0] && box1[0][0] < box2[0][1]) || // min contained in rng |
| (box2[0][0] > box1[0][0] && box2[0][0] < box1[0][1]); // min contained in rng |
| if (!x_neighbor) { |
| x_neighbor = (box1[0][1] == box2[0][0]-1) || (box1[0][0] == box2[0][1]+1); |
| } |
| |
| bool y_neighbor = (box1[1][1] == box2[1][0]) || (box1[1][0] == box2[1][1]) || // min matches max |
| (box1[1][0] == box2[1][0]) || (box1[1][1] == box2[1][1]) || // mins or maxs match |
| (box1[1][0] > box2[1][0] && box1[1][1] < box2[1][1]) || // range contains other |
| (box2[1][0] > box1[1][0] && box2[1][1] < box1[1][1]) || // range contains other |
| (box1[1][0] > box2[1][0] && box1[1][0] < box2[1][1]) || // min contained in rng |
| (box2[1][0] > box1[1][0] && box2[1][0] < box1[1][1]); // min contained in rng |
| if (!y_neighbor) { |
| y_neighbor = (box1[1][1] == box2[1][0]-1) || (box1[1][0] == box2[1][1]+1); |
| } |
| |
| bool z_neighbor = (box1[2][1] == box2[2][0]) || (box1[2][0] == box2[2][1]) || // min matches max |
| (box1[2][0] == box2[2][0]) || (box1[2][1] == box2[2][1]) || // mins or maxs match |
| (box1[2][0] > box2[2][0] && box1[2][1] < box2[2][1]) || // range contains other |
| (box2[2][0] > box1[2][0] && box2[2][1] < box1[2][1]) || // range contains other |
| (box1[2][0] > box2[2][0] && box1[2][0] < box2[2][1]) || // min contained in rng |
| (box2[2][0] > box1[2][0] && box2[2][0] < box1[2][1]); // min contained in rng |
| if (!z_neighbor) { |
| z_neighbor = (box1[2][1] == box2[2][0]-1) || (box1[2][0] == box2[2][1]+1); |
| } |
| |
| return x_neighbor && y_neighbor && z_neighbor; |
| } |
| |
| template<typename GlobalOrdinal> |
| void create_map_id_to_row(int global_nx, int global_ny, int global_nz, |
| const Box& box, |
| std::map<GlobalOrdinal,GlobalOrdinal>& id_to_row) |
| { |
| GlobalOrdinal num_my_ids = get_num_ids<GlobalOrdinal>(box); |
| |
| typename std::vector<GlobalOrdinal> all_ids; |
| bool include_ghost_layer = false; |
| get_ids(global_nx, global_ny, global_nz, box, all_ids, include_ghost_layer); |
| |
| GlobalOrdinal my_first_row = 0; |
| typename std::vector<GlobalOrdinal> global_offsets; |
| std::vector<int> all_boxes; |
| int numprocs = 1, myproc = 0; |
| #ifdef HAVE_MPI |
| MPI_Comm_size(MPI_COMM_WORLD, &numprocs); |
| MPI_Comm_rank(MPI_COMM_WORLD, &myproc); |
| |
| GlobalOrdinal local_num_ids = num_my_ids; |
| global_offsets.resize(numprocs); |
| MPI_Datatype mpi_dtype = TypeTraits<GlobalOrdinal>::mpi_type(); |
| MPI_Allgather(&local_num_ids, 1, mpi_dtype, &global_offsets[0], 1, mpi_dtype, MPI_COMM_WORLD); |
| GlobalOrdinal offset = 0; |
| for(int i=0; i<numprocs; ++i) { |
| GlobalOrdinal tmp = global_offsets[i]; |
| global_offsets[i] = offset; |
| offset += tmp; |
| } |
| |
| my_first_row = global_offsets[myproc]; |
| |
| all_boxes.resize(6*numprocs); |
| int* local_box_ranges = const_cast<int*>(&box.ranges[0]); |
| MPI_Allgather(local_box_ranges, 6, MPI_INT, &all_boxes[0], 6, MPI_INT, MPI_COMM_WORLD); |
| #endif |
| |
| if (all_ids.size() > 0) { |
| id_to_row.insert(std::make_pair(all_ids[0], my_first_row)); |
| } |
| |
| for(size_t i=1; i<all_ids.size(); ++i) { |
| if (all_ids[i] != all_ids[i-1]+1) { |
| id_to_row.insert(std::make_pair(all_ids[i], my_first_row+i)); |
| } |
| } |
| |
| // int num_neighbors = 0; |
| for(int i=0; i<numprocs; ++i) { |
| if (i == myproc) continue; |
| Box box_i; |
| for(int r=0; r<6; ++r) box_i.ranges[r] = all_boxes[i*6 + r]; |
| // bool neighbor= is_neighbor(box, box_i); |
| //if(myproc==2) { |
| // std::cout<<"i: "<<i<<" "<<neighbor<<" "; |
| // print_box(myproc, " ", box, " ", box_i); |
| //} |
| if (!is_neighbor(box, box_i)) { |
| // if (myproc==50) { |
| // std::cout<<"box ("<<box[0][0]<<","<<box[0][1]<<" - "<<box[1][0]<<","<<box[1][1]<<" - "<<box[2][0]<<","<<box[2][1]<<")"<<std::endl<<" and ("<<box_i[0][0]<<","<<box_i[0][1]<<" - "<<box_i[1][0]<<","<<box_i[1][1]<<" - "<<box_i[2][0]<<","<<box_i[2][1]<<") not neighbors."<<std::endl; |
| // } |
| continue; |
| } |
| // ++num_neighbors; |
| |
| get_ids(global_nx, global_ny, global_nz, box_i, all_ids, include_ghost_layer); |
| |
| GlobalOrdinal first_row = global_offsets[i]; |
| if (all_ids.size() > 0) { |
| id_to_row.insert(std::make_pair(all_ids[0], first_row)); |
| } |
| for(size_t j=1; j<all_ids.size(); ++j) { |
| if (all_ids[j] != all_ids[j-1]+1) { |
| id_to_row.insert(std::make_pair(all_ids[j], first_row+j)); |
| } |
| } |
| } |
| |
| //std::cout<<"proc "<<myproc<<": num_neighbors: "<<num_neighbors<<", id_to_row.size(): "<<id_to_row.size()<<std::endl; |
| //typename std::map<GlobalOrdinal,GlobalOrdinal>::iterator iter = id_to_row.begin(), end = id_to_row.end(); |
| //for(; iter!=end; ++iter) { |
| // std::cout<<"proc "<<myproc<<": "<<iter->first<<" :: "<<iter->second<<std::endl; |
| //} |
| } |
| |
| }//namespace miniFE |
| |
| #endif |
| |