#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