
#ifndef _simple_mesh_description_hpp_
#define _simple_mesh_description_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 <utils.hpp>
#include <set>
#include <map>

namespace miniFE {

template<typename GlobalOrdinal>
class simple_mesh_description {
public:
  simple_mesh_description(const Box& global_box_in, const Box& local_box_in)
  {
   Box local_node_box;
    for(int i=0; i<3; ++i) {
      global_box[i][0] = global_box_in[i][0];
      global_box[i][1] = global_box_in[i][1];
      local_box[i][0] = local_box_in[i][0];
      local_box[i][1] = local_box_in[i][1];
      local_node_box[i][0] = local_box_in[i][0];
      local_node_box[i][1] = local_box_in[i][1];
      //num-owned-nodes == num-elems+1 in this dimension if the elem box is not empty
      //and we are at the high end of the global range in that dimension:
      if (local_box_in[i][1] > local_box_in[i][0] && local_box_in[i][1] == global_box[i][1]) local_node_box[i][1] += 1;
    }

    int max_node_x = global_box[0][1]+1;
    int max_node_y = global_box[1][1]+1;
    int max_node_z = global_box[2][1]+1;
    create_map_id_to_row(max_node_x, max_node_y, max_node_z, local_node_box,
                         map_ids_to_rows);

    //As described in analytic_soln.hpp,
    //we will impose a 0 boundary-condition on faces x=0, y=0, z=0, y=1, z=1
    //we will impose a 1 boundary-condition on face x=1

#ifdef MINIFE_DEBUG
std::cout<<std::endl;
#endif
    const int X=0;
    const int Y=1;
    const int Z=2;

    const int x1 = max_node_x - 1;
    const int y1 = max_node_y - 1;
    const int z1 = max_node_z - 1;

    //if we're on the x=0 face:
    if (global_box[X][0] == local_box[X][0]) {
      int miny = local_node_box[Y][0];
      int minz = local_node_box[Z][0];
      int maxy = local_node_box[Y][1];
      int maxz = local_node_box[Z][1];
      //expand y and z dimensions to include ghost layer
      if (local_node_box[Y][0] > 0) --miny;
      if (local_node_box[Z][0] > 0) --minz;
      if (local_node_box[Y][1] < max_node_y) ++maxy;
      if (local_node_box[Z][1] < max_node_z) ++maxz;

      for(int iz=minz; iz<maxz; ++iz) {
        for(int iy=miny; iy<maxy; ++iy) {
          GlobalOrdinal nodeID = get_id<GlobalOrdinal>(max_node_x, max_node_y, max_node_z,
             0, iy, iz);
#ifdef MINIFE_DEBUG
std::cout<<"x=0 BC, node "<<nodeID<<", (0,"<<iy<<","<<iz<<")"<<std::endl;
#endif
          bc_rows_0.insert(map_id_to_row(nodeID));
        }
      }
    }

    //if we're on the y=0 face:
    if (global_box[Y][0] == local_box[Y][0]) {
      int minx = local_node_box[X][0];
      int minz = local_node_box[Z][0];
      int maxx = local_node_box[X][1];
      int maxz = local_node_box[Z][1];
      //expand x and z dimensions to include ghost layer
      if (local_node_box[X][0] > 0) --minx;
      if (local_node_box[Z][0] > 0) --minz;
      if (local_node_box[X][1] < max_node_x) ++maxx;
      if (local_node_box[Z][1] < max_node_z) ++maxz;

      for(int iz=minz; iz<maxz; ++iz) {
        for(int ix=minx; ix<maxx; ++ix) {
          GlobalOrdinal nodeID = get_id<GlobalOrdinal>(max_node_x, max_node_y, max_node_z,
             ix, 0, iz);
#ifdef MINIFE_DEBUG
std::cout<<"y=0 BC, node "<<nodeID<<", ("<<ix<<",0,"<<iz<<")"<<std::endl;
#endif
          GlobalOrdinal row = map_id_to_row(nodeID);
          if (row < 0) {
            std::cout<<"on the y==0 face (ix="<<ix<<", iz="<<iz<<"), ERROR: found negative row ("<<row<<") for nodeID="<<nodeID<<std::endl;
          }
          bc_rows_0.insert(row);
        }
      }
    }

    //if we're on the z=0 face:
    if (global_box[Z][0] == local_box[Z][0]) {
      int minx = local_node_box[X][0];
      int miny = local_node_box[Y][0];
      int maxx = local_node_box[X][1];
      int maxy = local_node_box[Y][1];
      //expand x and y dimensions to include ghost layer
      if (local_node_box[X][0] > 0) --minx;
      if (local_node_box[Y][0] > 0) --miny;
      if (local_node_box[X][1] < max_node_x) ++maxx;
      if (local_node_box[Y][1] < max_node_y) ++maxy;

      for(int iy=miny; iy<maxy; ++iy) {
        for(int ix=minx; ix<maxx; ++ix) {
          GlobalOrdinal nodeID = get_id<GlobalOrdinal>(max_node_x, max_node_y, max_node_z,
             ix, iy, 0);
#ifdef MINIFE_DEBUG
std::cout<<"z=0 BC, node "<<nodeID<<", ("<<ix<<","<<iy<<",0)"<<std::endl;
#endif
          bc_rows_0.insert(map_id_to_row(nodeID));
        }
      }
    }

    //if we're on the x=1 face:
    if (global_box[X][1] == local_box[X][1]) {
      int minz = local_node_box[Z][0];
      int miny = local_node_box[Y][0];
      int maxz = local_node_box[Z][1];
      int maxy = local_node_box[Y][1];
      //expand z and y dimensions to include ghost layer
      if (local_node_box[Z][0] > 0) --minz;
      if (local_node_box[Y][0] > 0) --miny;
      if (local_node_box[Z][1] < max_node_z) ++maxz;
      if (local_node_box[Y][1] < max_node_y) ++maxy;

      for(int iy=miny; iy<maxy; ++iy) {
        for(int iz=minz; iz<maxz; ++iz) {
          GlobalOrdinal nodeID = get_id<GlobalOrdinal>(max_node_x, max_node_y, max_node_z,
             x1, iy, iz);
          GlobalOrdinal row = map_id_to_row(nodeID);
#ifdef MINIFE_DEBUG
std::cout<<"x=1 BC, node "<<nodeID<<", row "<<row<<", ("<<x1<<","<<iy<<","<<iz<<")"<<std::endl;
#endif
          bc_rows_1.insert(row);
        }
      }
    }

    //if we're on the y=1 face:
    if (global_box[Y][1] == local_box[Y][1]) {
      int minz = local_node_box[Z][0];
      int minx = local_node_box[X][0];
      int maxz = local_node_box[Z][1];
      int maxx = local_node_box[X][1];
      //expand z and x dimensions to include ghost layer
      if (local_node_box[Z][0] > 0) --minz;
      if (local_node_box[X][0] > 0) --minx;
      if (local_node_box[Z][1] < max_node_z) ++maxz;
      if (local_node_box[X][1] < max_node_x) ++maxx;

      for(int ix=minx; ix<maxx; ++ix) {
        for(int iz=minz; iz<maxz; ++iz) {
          GlobalOrdinal nodeID = get_id<GlobalOrdinal>(max_node_x, max_node_y, max_node_z,
             ix, y1, iz);
#ifdef MINIFE_DEBUG
std::cout<<"y=1 BC, node "<<nodeID<<", ("<<ix<<","<<y1<<","<<iz<<")"<<std::endl;
#endif
          bc_rows_0.insert(map_id_to_row(nodeID));
        }
      }
    }

    //if we're on the z=1 face:
    if (global_box[Z][1] == local_box[Z][1]) {
      int miny = local_node_box[Y][0];
      int minx = local_node_box[X][0];
      int maxy = local_node_box[Y][1];
      int maxx = local_node_box[X][1];
      //expand x and y dimensions to include ghost layer
      if (local_node_box[Y][0] > 0) --miny;
      if (local_node_box[X][0] > 0) --minx;
      if (local_node_box[Y][1] < max_node_y) ++maxy;
      if (local_node_box[X][1] < max_node_x) ++maxx;

      for(int ix=minx; ix<maxx; ++ix) {
        for(int iy=miny; iy<maxy; ++iy) {
          GlobalOrdinal nodeID = get_id<GlobalOrdinal>(max_node_x, max_node_y, max_node_z,
             ix, iy, z1);
#ifdef MINIFE_DEBUG
std::cout<<"z=1 BC, node "<<nodeID<<", ("<<ix<<","<<iy<<","<<z1<<")"<<std::endl;
#endif
          bc_rows_0.insert(map_id_to_row(nodeID));
        }
      }
    }

  }

  GlobalOrdinal map_id_to_row(const GlobalOrdinal& id) const
  {
    return find_row_for_id(id, map_ids_to_rows);
  }

  GlobalOrdinal max_row_in_map() const {
    if (map_ids_to_rows.empty()) return 0;
    typename std::map<GlobalOrdinal,GlobalOrdinal>::const_iterator mend = map_ids_to_rows.end();
    --mend;
    return mend->second;
  }
  std::set<GlobalOrdinal> bc_rows_0;
  std::set<GlobalOrdinal> bc_rows_1;
  std::map<GlobalOrdinal,GlobalOrdinal> map_ids_to_rows;
  Box global_box;
  Box local_box;
};//class simple_mesh_description

}//namespace miniFE

#endif