#ifndef _Hex8_box_utils_hpp_
#define _Hex8_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 <stdexcept>

#include <box_utils.hpp>
#include <ElemData.hpp>
#include <simple_mesh_description.hpp>
#include <Hex8.hpp>

namespace miniFE {


template<typename GlobalOrdinal>
void get_hex8_node_ids(int nx, int ny,
                       GlobalOrdinal node0,
                       GlobalOrdinal* elem_node_ids)
{
//Given box dimensions nx and ny, and a starting node
//(local-node-0 for a hex8), compute the other nodes
//of the hex8 using the exodus ordering convention.
  elem_node_ids[0] = node0;
  elem_node_ids[1] = node0 + 1;
  elem_node_ids[2] = node0 + nx + 1;
  elem_node_ids[3] = node0 + nx;
  elem_node_ids[4] = node0 +     nx*ny;
  elem_node_ids[5] = node0 + 1 + nx*ny;
  elem_node_ids[6] = node0 + nx + nx*ny + 1;
  elem_node_ids[7] = node0 + nx + nx*ny;
}

template<typename Scalar>
void get_hex8_node_coords_3d(Scalar x, Scalar y, Scalar z,
                             Scalar hx, Scalar hy, Scalar hz,
                             Scalar* elem_node_coords)
{
  //Input: x,y,z are the coordinates of local-node 0 for a Hex8.
  //'hx', 'hy', 'hz' are the lengths of the sides of the element
  //in each direction.

  elem_node_coords[0] = x;
  elem_node_coords[1] = y;
  elem_node_coords[2] = z;

  elem_node_coords[3] = x + hx;
  elem_node_coords[4] = y;
  elem_node_coords[5] = z;

  elem_node_coords[6] = x + hx;
  elem_node_coords[7] = y + hy;
  elem_node_coords[8] = z;

  elem_node_coords[9]  = x;
  elem_node_coords[10] = y + hy;
  elem_node_coords[11] = z;

  elem_node_coords[12] = x;
  elem_node_coords[13] = y;
  elem_node_coords[14] = z + hz;

  elem_node_coords[15] = x + hx;
  elem_node_coords[16] = y;
  elem_node_coords[17] = z + hz;

  elem_node_coords[18] = x + hx;
  elem_node_coords[19] = y + hy;
  elem_node_coords[20] = z + hz;

  elem_node_coords[21] = x;
  elem_node_coords[22] = y + hy;
  elem_node_coords[23] = z + hz;
}

template<typename GlobalOrdinal, typename Scalar>
void
get_elem_nodes_and_coords(const simple_mesh_description<GlobalOrdinal>& mesh,
                          GlobalOrdinal elemID,
                          GlobalOrdinal* node_ords, Scalar* node_coords)
{
  int global_nodes_x = mesh.global_box[0][1]+1;
  int global_nodes_y = mesh.global_box[1][1]+1;
  int global_nodes_z = mesh.global_box[2][1]+1;
 
  if (elemID < 0) {
    //I don't think this can happen, but check for the sake of paranoia...
    throw std::runtime_error("get_elem_nodes_and_coords ERROR, negative elemID");
  }

  int elem_int_x, elem_int_y, elem_int_z;
  get_int_coords(elemID, global_nodes_x-1, global_nodes_y-1, global_nodes_z-1,
             elem_int_x, elem_int_y, elem_int_z);
  GlobalOrdinal nodeID = get_id<GlobalOrdinal>(global_nodes_x, global_nodes_y, global_nodes_z, elem_int_x, elem_int_y, elem_int_z);

#ifdef MINIFE_DEBUG_VERBOSE
  std::cout<<"\nelemID: "<<elemID<<", nodeID: "<<nodeID<<std::endl;
#endif
  get_hex8_node_ids(global_nodes_x, global_nodes_y, nodeID, node_ords);

  //Map node-IDs to rows because each processor may have a non-contiguous block of
  //node-ids, but needs a contiguous block of row-numbers:
#ifdef MINIFE_DEBUG_VERBOSE
  std::cout<<"elem "<<elemID<<" nodes: ";
#endif
  for(int i=0; i<Hex8::numNodesPerElem; ++i) {
#ifdef MINIFE_DEBUG_VERBOSE
    std::cout<<node_ords[i]<<" ";
#endif
    node_ords[i] = mesh.map_id_to_row(node_ords[i]);
  }
#ifdef MINIFE_DEBUG_VERBOSE
  std::cout << std::endl;
#endif

  int global_elems_x = mesh.global_box[0][1];
  int global_elems_y = mesh.global_box[1][1];
  int global_elems_z = mesh.global_box[2][1];
 
  Scalar ix,iy,iz;
  get_coords<GlobalOrdinal,Scalar>(nodeID, global_nodes_x,global_nodes_y,global_nodes_z,
                            ix,iy,iz);
  Scalar hx = 1.0/global_elems_x;
  Scalar hy = 1.0/global_elems_y;
  Scalar hz = 1.0/global_elems_z;
  get_hex8_node_coords_3d(ix, iy, iz, hx, hy, hz, node_coords);
#ifdef MINIFE_DEBUG_VERBOSE
  int offset = 0;
  for(int i=0; i<Hex8::numNodesPerElem; ++i) {
    std::cout << "("<<node_coords[offset++]<<","<<node_coords[offset++]<<","<<node_coords[offset++]<<")";
  }
  std::cout << std::endl;
#endif
}

template<typename GlobalOrdinal, typename Scalar>
void
get_elem_nodes_and_coords(const simple_mesh_description<GlobalOrdinal>& mesh,
                          GlobalOrdinal elemID,
                          ElemData<GlobalOrdinal,Scalar>& elem_data)
{
  get_elem_nodes_and_coords(mesh, elemID, elem_data.elem_node_ids, elem_data.elem_node_coords);
}

}//namespace miniFE

#endif