#ifndef __topology_spa_mesh_h__
#define __topology_spa_mesh_h__

#include "Mesh.hpp"
#include "ReferenceFiniteElement.hpp"
#include "ElementIterator.hpp"

class MeshFEDataAccessor
{
public:
  typedef SpaMesh::Connectivity::TopoIndexVector TopoIndexVector;
  typedef SpaMesh::Connectivity::TopoPointVector TopoPointVector;

  MeshFEDataAccessor(unsigned int Dim) :
    _dim(Dim)
  {
    ;
  }
  
  ~MeshFEDataAccessor() {}
  template <unsigned int Dim>
  void initData(const SpaMesh::Mesh<Dim>& mesh, 
                const int& elem)
  {
    int elemTypeId;
    int spaMeshElemTypeDim = mesh.elemMesh.typeFromLID(elem, elemTypeId);
    mesh.getCoords(mesh.elemMesh.elemNodesConnectivity[spaMeshElemTypeDim], elemTypeId, _elem_coord);
    _spaMeshElemType = static_cast<SpaMesh::ElemTypesMesh>(mesh.elemMesh.elemTypeFromDimType(spaMeshElemTypeDim));
    _num_elem_nodes = _elem_coord.rows();
  }
  template <unsigned int Dim>
  void initDataSide(const SpaMesh::Mesh<Dim>& mesh, 
                    const int& side)
  {
    int sideTypeId;
    int spaMeshElemTypeDim = mesh.faceMesh.typeFromLID(side, sideTypeId);
    mesh.getCoords(mesh.faceMesh.elemNodesConnectivity[spaMeshElemTypeDim], sideTypeId, _elem_coord);
    _spaMeshElemType = static_cast<SpaMesh::ElemTypesMesh>(mesh.faceMesh.elemTypeFromDimType(spaMeshElemTypeDim));
    _num_elem_nodes = _elem_coord.rows();
  }

  template <unsigned int Dim>
  void initData(const SpaMesh::Mesh<Dim>& mesh, 
                const SpaMesh::ElementIterator& elemIt)
  {
    int elemTypeId;
    int spaMeshElemTypeDim = mesh.elemMesh.typeFromLID(elemIt.value(), elemTypeId);
    mesh.getCoords(mesh.elemMesh.elemNodesConnectivity[spaMeshElemTypeDim], elemTypeId, _elem_coord);
    _spaMeshElemType = static_cast<SpaMesh::ElemTypesMesh>(mesh.elemMesh.elemTypeFromDimType(spaMeshElemTypeDim));
    _num_elem_nodes = _elem_coord.rows();
  }
  
  template <unsigned int Dim>
  void initData(const SpaMesh::Mesh<Dim>& mesh, 
                const SpaMesh::NeighborIterator& elemIt)
  {
    mesh.getCoords(mesh.elemMesh.elemNodesConnectivity[elemIt.getType()], elemIt.getTypeID(), _elem_coord);
    _spaMeshElemType = static_cast<SpaMesh::ElemTypesMesh>(mesh.elemMesh.elemTypeFromDimType(elemIt.getType()));
    _num_elem_nodes = _elem_coord.rows();
  }

  template <unsigned int Dim>
  void initData(const SpaMesh::Mesh<Dim>& mesh, 
                const SpaMesh::ElementIterator& elemIt, 
		const SpaMesh::NeighborIterator& neighborIt, 
		const bool referenceFrameSideInfo = true)
  {
    // connectivity for elem and neighbor
    _elem_side_number = neighborIt.index();
    
    int elemTypeId;
    int spaMeshElemTypeDim = mesh.elemMesh.typeFromLID(elemIt.value(), elemTypeId);
    _spaMeshElemType = static_cast<SpaMesh::ElemTypesMesh>(mesh.elemMesh.elemTypeFromDimType(spaMeshElemTypeDim));
    mesh.getCoords(mesh.elemMesh.elemNodesConnectivity[spaMeshElemTypeDim], elemTypeId, _elem_coord);
    _num_elem_nodes = _elem_coord.rows();
    
    int neighborTypeId;
    int spaMeshNeighborTypeDim = mesh.elemMesh.typeFromLID(neighborIt.value(), neighborTypeId);
    _spaMeshNeighborType = static_cast<SpaMesh::ElemTypesMesh>(mesh.elemMesh.elemTypeFromDimType(spaMeshNeighborTypeDim));
    mesh.getCoords(mesh.elemMesh.elemNodesConnectivity[spaMeshNeighborTypeDim], neighborTypeId, _neighbor_coord);
    _num_neighbor_nodes = _neighbor_coord.rows();

    // reference frame node indices and physical frame nodes coordinates for the side
    _elem_side_indices.resize(SpaMesh::MAX_SUB_ENTITY_VERTICES);
    SpaMesh::SpaCN::SubEntityNodeIndices(_spaMeshElemType, _num_elem_nodes, _dim-1, _elem_side_number, _spaMeshSideType, _num_side_nodes, &_elem_side_indices(0));
    _elem_side_indices.conservativeResize(_num_side_nodes);
    int side_conn[_num_side_nodes];
    _side_coord.resize(_num_side_nodes,Eigen::NoChange);
    for (unsigned int i = 0; i < _num_side_nodes; i++)
    {
      side_conn[i] = mesh.elemMesh.elemNodesConnectivity[spaMeshElemTypeDim](_elem_side_indices[i],elemTypeId);
      for (unsigned int j = 0; j < 3; j++)
        _side_coord(i,j) = _elem_coord(_elem_side_indices[i],j);
    }
    
    _neighbor_side_indices.resize(_num_side_nodes);
    if (referenceFrameSideInfo)
    {
      for (unsigned int i = 0; i < _num_side_nodes; i++)
        for(unsigned int j = 0; j < _num_neighbor_nodes; j++)
  	  if (side_conn[i] == mesh.elemMesh.elemNodesConnectivity[spaMeshNeighborTypeDim](j,neighborTypeId))  
            _neighbor_side_indices[i] = j;	
    }
    
     assert (_num_side_verts <= _num_side_nodes < _num_elem_nodes && _num_side_nodes < _num_neighbor_nodes);
  }

  template <unsigned int Dim>
  void initData(const SpaMesh::Mesh<Dim>& mesh, 
                const SpaMesh::ElementIterator& elemIt, 
		const int side_number, 
		const bool referenceFrameSideInfo = true)
  {
    // connectivity for elem and neighbor
    _elem_side_number = side_number;
    int elemTypeId;
    int spaMeshElemTypeDim = mesh.elemMesh.typeFromLID(elemIt.value(), elemTypeId);
    _spaMeshElemType = static_cast<SpaMesh::ElemTypesMesh>(mesh.elemMesh.elemTypeFromDimType(spaMeshElemTypeDim));
    mesh.getCoords(mesh.elemMesh.elemNodesConnectivity[spaMeshElemTypeDim], elemTypeId, _elem_coord);
  
    // reference frame node indices and physical frame nodes coordinates for the side
    _num_elem_nodes = _elem_coord.rows();
    _elem_side_indices.resize(SpaMesh::MAX_SUB_ENTITY_VERTICES);
    SpaMesh::SpaCN::SubEntityNodeIndices(_spaMeshElemType, _num_elem_nodes, _dim-1, _elem_side_number, _spaMeshSideType, _num_side_nodes, &_elem_side_indices(0));
    _elem_side_indices.conservativeResize(_num_side_nodes);
    int side_conn[_num_side_nodes];
    _side_coord.resize(_num_side_nodes,Eigen::NoChange);
    for (unsigned int i = 0; i < _num_side_nodes; i++)
    {
      side_conn[i] = mesh.elemMesh.elemNodesConnectivity[spaMeshElemTypeDim](_elem_side_indices[i],elemTypeId);
      for (unsigned int j = 0; j < 3; j++)
        _side_coord(i,j) = _elem_coord(_elem_side_indices[i],j);
    }
    
    if (referenceFrameSideInfo)
    {
      _num_side_verts = SpaMesh::SpaCN::VerticesPerEntity(_spaMeshSideType);
    }
    assert(_num_side_verts <= _num_side_nodes < _num_elem_nodes);
  }

  template <unsigned int Dim>
  void initData(const SpaMesh::Mesh<Dim>& mesh, 
                const int& elem, 
		const int& neighbor, 
		const SpaMesh::FaceIterator& sideIt, 
		const bool referenceFrameSideInfo = true)
  {
    int elemTypeId;
    int spaMeshElemTypeDim = mesh.elemMesh.typeFromLID(elem, elemTypeId);
    _spaMeshElemType = static_cast<SpaMesh::ElemTypesMesh>(mesh.elemMesh.elemTypeFromDimType(spaMeshElemTypeDim));
    mesh.getCoords(mesh.elemMesh.elemNodesConnectivity[spaMeshElemTypeDim], elemTypeId, _elem_coord);
    _num_elem_nodes = _elem_coord.rows();
    
    int neighborTypeId;
    int spaMeshNeighborTypeDim = mesh.elemMesh.typeFromLID(neighbor, neighborTypeId);
    _spaMeshNeighborType = static_cast<SpaMesh::ElemTypesMesh>(mesh.elemMesh.elemTypeFromDimType(spaMeshNeighborTypeDim));
    mesh.getCoords(mesh.elemMesh.elemNodesConnectivity[spaMeshNeighborTypeDim], neighborTypeId, _neighbor_coord);
    _num_neighbor_nodes = _neighbor_coord.rows();

    int sideTypeId;
    int spaMeshSideTypeDim = mesh.faceMesh.typeFromLID(sideIt.value(), sideTypeId);
    _spaMeshSideType = static_cast<SpaMesh::ElemTypesMesh>(mesh.faceMesh.elemTypeFromDimType(spaMeshSideTypeDim));
    if (_spaMeshSideType == SpaMesh::node)
      _num_side_nodes = 1;
    else 
       _num_side_nodes = mesh.faceMesh.elemNodesConnectivity[spaMeshSideTypeDim].rows();
    _num_side_verts = SpaMesh::SpaCN::VerticesPerEntity(_spaMeshSideType);
    int sense;
    int offset;
    SpaMesh::SpaCN::SideNumber(_spaMeshElemType, 
	                       &mesh.elemMesh.elemNodesConnectivity[spaMeshElemTypeDim](0,elemTypeId), 
			       &mesh.faceMesh.elemNodesConnectivity[spaMeshSideTypeDim](0,sideTypeId),
			       _num_side_verts, _dim-1, _elem_side_number, sense, offset);
    _elem_side_indices.resize(_num_side_nodes);
    SpaMesh::SpaCN::SubEntityNodeIndices(_spaMeshElemType, _num_elem_nodes, _dim-1, _elem_side_number, _spaMeshSideType, _num_side_nodes, &_elem_side_indices(0));
  
    _side_coord.resize(_num_side_nodes,Eigen::NoChange);
    int side_conn[_num_side_nodes];
    for (unsigned int i = 0; i < _num_side_nodes; i++)
    {
      side_conn[i] = mesh.elemMesh.elemNodesConnectivity[spaMeshElemTypeDim](_elem_side_indices[i],elemTypeId);
      for (unsigned int d = 0; d < 3; d++)
        _side_coord(i,d) = _elem_coord(_elem_side_indices[i],d);
    }
  
    _neighbor_side_indices.resize(_num_side_nodes);
    if (referenceFrameSideInfo)
    {
      for (unsigned int i = 0; i < _num_side_nodes; i++)
        for(unsigned int j = 0; j < _num_neighbor_nodes; j++)
  	  if (side_conn[i] == mesh.elemMesh.elemNodesConnectivity[spaMeshNeighborTypeDim](j,neighborTypeId))  
            _neighbor_side_indices[i] = j;	
    }
    assert (_num_side_verts <= _num_side_nodes < _num_elem_nodes && _num_side_nodes < _num_neighbor_nodes);
  }
  
  template <unsigned int Dim>
    void initData(const SpaMesh::Mesh<Dim>& mesh, const int& elem, const int& side, const bool referenceFrameSideInfo = true, const bool faceInfo = false)
  {
    int elemTypeId;
    int spaMeshElemTypeDim = mesh.elemMesh.typeFromLID(elem, elemTypeId);
    _spaMeshElemType = static_cast<SpaMesh::ElemTypesMesh>(mesh.elemMesh.elemTypeFromDimType(spaMeshElemTypeDim));
    mesh.getCoords(mesh.elemMesh.elemNodesConnectivity[spaMeshElemTypeDim], elemTypeId, _elem_coord);
    _num_elem_nodes = _elem_coord.rows();
    
    int sideTypeId;
    int spaMeshSideTypeDim = mesh.faceMesh.typeFromLID(side, sideTypeId);
    _spaMeshSideType = static_cast<SpaMesh::ElemTypesMesh>(mesh.faceMesh.elemTypeFromDimType(spaMeshSideTypeDim));
    if (_spaMeshSideType == SpaMesh::node)
      _num_side_nodes = 1;
    else 
       _num_side_nodes = mesh.faceMesh.elemNodesConnectivity[spaMeshSideTypeDim].rows();
    _num_side_verts = SpaMesh::SpaCN::VerticesPerEntity(_spaMeshSideType);
    int sense;
    int offset;
    SpaMesh::SpaCN::SideNumber(_spaMeshElemType, 
	                       &mesh.elemMesh.elemNodesConnectivity[spaMeshElemTypeDim](0,elemTypeId), 
			       &mesh.faceMesh.elemNodesConnectivity[spaMeshSideTypeDim](0,sideTypeId),
			       _num_side_verts, _dim-1, _elem_side_number, sense, offset);
    _elem_side_indices.resize(_num_side_nodes);
    SpaMesh::SpaCN::SubEntityNodeIndices(_spaMeshElemType, _num_elem_nodes, _dim-1, _elem_side_number, _spaMeshSideType, _num_side_nodes, &_elem_side_indices(0));
  
    _side_coord.resize(_num_side_nodes,Eigen::NoChange);
    for (unsigned int i = 0; i < _num_side_nodes; i++)
      for (unsigned int d = 0; d < 3; d++)
        _side_coord(i,d) = _elem_coord(_elem_side_indices[i],d);
  
    if (faceInfo)
      mesh.getCoords(mesh.faceMesh.elemNodesConnectivity[spaMeshSideTypeDim], sideTypeId, _face_coord);
    
    assert (_num_side_verts <= _num_side_nodes < _num_elem_nodes);
  }

  const TopoPointVector& getElemCoord() const
  {
    return _elem_coord;
  }
  const TopoPointVector& getNeighborCoord() const
  {
    return _neighbor_coord;
  }
  const TopoPointVector& getSideCoord() const
  {
    return _side_coord;
  }
  const TopoPointVector& getFaceCoord() const
  {
    return _face_coord;
  }
  const TopoIndexVector& getElemSideIndices() const
  {
    return _elem_side_indices;
  }
  const TopoIndexVector& getElemSideIndices(int & elem_side_number) const
  {
    elem_side_number = _elem_side_number;
    return _elem_side_indices;
  }
  const TopoIndexVector& getNeighborSideIndices() const
  {
    return _neighbor_side_indices;
  }
  ElemType getElemType() const
  {
    return ReferenceFiniteElement::elemTypeFromEntityType(_spaMeshElemType, _num_elem_nodes);
  }
  ElemType getNeighborType() const
  {
    return ReferenceFiniteElement::elemTypeFromEntityType(_spaMeshNeighborType, _num_neighbor_nodes);
  }
  ElemType getSideType() const
  {
    return ReferenceFiniteElement::elemTypeFromEntityType(_spaMeshSideType, _num_side_nodes);
  }

protected:
  const unsigned int _dim;

  TopoPointVector _elem_coord;
  TopoPointVector _neighbor_coord;
  TopoPointVector _side_coord;
  TopoPointVector _face_coord;

  TopoIndexVector _elem_side_indices;
  TopoIndexVector _neighbor_side_indices;

  int _elem_side_number;
  int _num_elem_nodes;
  int _num_neighbor_nodes;
  int _num_side_nodes;
  int _num_side_verts;

  SpaMesh::ElemTypesMesh _spaMeshElemType; 
  SpaMesh::ElemTypesMesh _spaMeshNeighborType; 
  SpaMesh::ElemTypesMesh _spaMeshSideType;
};
#endif
