// SpaFEDTe, a Template based C++ library for creating 
// Discontinuous Finite Element Spaces,
// Copyright (C) 2012 Lorenzo Alessio Botti

/* 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 3.0 of the License, or (at your option) any later version. */

/* This software 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 software; if not, a copy of the full */
/* GNU Lesser General Public License can be found at */
/* http://www.gnu.org/licenses/ */

#ifndef __dof_map_h__
#define __dof_map_h__

#include "Eigen/Dense"
#include "Mesh.hpp"
#include "ElementIterator.hpp"

class DofMap
{
public:
  typedef Eigen::Matrix<std::string, 10, 1> VarNamesType;
  typedef Eigen::Matrix<double, 10, 1> VarScalingsType;

  DofMap():
      _n_local_elements(0),
      _n_global_elements(0),
      _n_vars(0),
      _p_order(0),
      _element_dofs(0),
      _element_dofs_col(0)
  {
    this->setVarNames("a");
     _var_scalings.setOnes();
  }
  
  DofMap(unsigned int n_vars, unsigned int p_order, unsigned int n_dofs):
      _n_local_elements(0),
      _n_global_elements(0),
      _n_vars(n_vars),
      _p_order(p_order),
      _element_dofs(n_dofs),
      _element_dofs_col(n_dofs)
  {
    this->setVarNames("a");
     _var_scalings.setOnes();
  }
  
  DofMap(unsigned int n_vars, unsigned int p_order, unsigned int n_dofs_row, unsigned int n_dofs_col):
      _n_local_elements(0),
      _n_global_elements(0),
      _n_vars(n_vars),
      _p_order(p_order),
      _element_dofs(n_dofs_row),
      _element_dofs_col(n_dofs_col)
  {
    this->setVarNames("a");
     _var_scalings.setOnes();
  }

  ~DofMap() {}
  
  const Eigen::VectorXi& getLocalNonZeroBlocks() const
  {
    return _local_nonzero_blocks;
  }
  
  const Eigen::VectorXi& getNonLocalNonZeroBlocks() const
  {
    return _nonlocal_nonzero_blocks;
  }

  const int& getLocalElementNumber() const
  {
    return _n_local_elements;
  }
  
  const int& getGhostElementNumber() const
  {
    return _n_ghost_elements;
  }
  
  const int& getGlobalElementNumber() const
  {
    return _n_global_elements;
  }

  Eigen::VectorXi& getLocalNonZeroBlocks()
  {
    return _local_nonzero_blocks;
  }
  
  Eigen::VectorXi& getNonLocalNonZeroBlocks()
  {
    return _nonlocal_nonzero_blocks;
  }
  
  Eigen::VectorXi& getGhostBlockGlobalIndices()
  {
    return _ghost_block_global_indices;
  }
  
  const Eigen::VectorXi& getGhostBlockGlobalIndices() const
  {
    return _ghost_block_global_indices;
  }

  void setLocalElementNumber(const int& local_elements_number)
  {
    _n_local_elements = local_elements_number;
  }
  
  void setGhostElementNumber(const int& ghost_elements_number)
  {
    _n_ghost_elements = ghost_elements_number;
  }
  
  void setGlobalElementNumber(const int& global_elements_number)
  {
    _n_global_elements = global_elements_number;
  }

  const unsigned int& getPOrder() const
  {
     return _p_order;
  }

  const unsigned int& getNVars() const
  {
     return _n_vars;
  }

  const unsigned int& getNDofs() const
  {
     return _element_dofs;
  }
  
  const unsigned int& getNDofsCol() const
  {
     return _element_dofs_col;
  }
  
  void setPOrder(const unsigned int& p_order) 
  {
     _p_order = p_order;
  }

  void setNVars(const unsigned int& n_vars) 
  {
     _n_vars = n_vars;
  }

  void setNDofs(const unsigned int& element_dofs) 
  {
     _element_dofs = element_dofs;
  }

  const VarNamesType& getVarNames() const
  {
     return _var_names;
  }
  
  const VarScalingsType& getVarScalings() const
  {
     return _var_scalings;
  }
  
  void setVarNames(std::string one, 
                   std::string two = "b", 
		   std::string three = "c",
		   std::string four = "d",
		   std::string five = "e",
		   std::string six = "f",
		   std::string seven = "g",
		   std::string eight = "h",
		   std::string nine = "i",
		   std::string ten = "l")
  {
     _var_names[0] = one;
     _var_names[1] = two;
     _var_names[2] = three;
     _var_names[3] = four;
     _var_names[4] = five;
     _var_names[5] = six;
     _var_names[6] = seven;
     _var_names[7] = eight;
     _var_names[8] = nine;
     _var_names[9] = ten;
  }
  
  void setVarScalings(double one, 
                      double two = 1., 
		      double three = 1.,
		      double four = 1.,
		      double five = 1.,
		      double six = 1.,
		      double seven = 1.,
		      double eight = 1.,
		      double nine = 1.,
		      double ten = 1.)
  {
     _var_scalings[0] = one;
     _var_scalings[1] = two;
     _var_scalings[2] = three;
     _var_scalings[3] = four;
     _var_scalings[4] = five;
     _var_scalings[5] = six;
     _var_scalings[6] = seven;
     _var_scalings[7] = eight;
     _var_scalings[8] = nine;
     _var_scalings[9] = ten;
  }

  template<unsigned int Dim>
  void computeMatrixSparsity(const SpaMesh::Mesh<Dim>& mesh)
  {
    _n_local_elements = mesh.elemCard();
    _n_global_elements = mesh.elemCard();
    _n_ghost_elements = 0;
    _local_nonzero_blocks = Eigen::VectorXi::Constant(_n_local_elements,1);
    _nonlocal_nonzero_blocks = Eigen::VectorXi::Zero(_n_local_elements);

    for (int elem = 0; elem < mesh.elemCard(); elem++)
    {
      SpaMesh::Connectivity::ElemConnVector neighborsVec = mesh.elemNeighbors(elem);
      unsigned int neighborsNumber = 0;
      for (unsigned int i = 0; i < neighborsVec.size(); i++)
      {
	if (neighborsVec[i] == elem)
	  continue;
	neighborsNumber++;
      }
      _local_nonzero_blocks[elem] += neighborsNumber;
    }
  }
  template<unsigned int Dim>
  void computeMatrixSparsityHHO(const SpaMesh::Mesh<Dim>& mesh)
  {
    _n_local_elements = mesh.faceCard();
    _n_global_elements = mesh.faceCard();
    _n_ghost_elements = 0;
    _local_nonzero_blocks = Eigen::VectorXi::Constant(_n_local_elements,1);
    _nonlocal_nonzero_blocks = Eigen::VectorXi::Zero(_n_local_elements);

    for (SpaMesh::FaceIterator faceIt(mesh); faceIt; ++faceIt)
    {
      SpaMesh::Connectivity::NeighborsPair neighborsPair = mesh.faceMesh.sideElemsConnectivity[faceIt.getType()].col(faceIt.getTypeID());
      int neighborsPairSize = 2;
      if (neighborsPair[0] == neighborsPair[1])
	neighborsPairSize = 1;
      int faceFound = 0;
      int facesNumber = 0;
      for (unsigned int n = 0; n < neighborsPairSize; n++)
      {
        SpaMesh::Connectivity::ElemConnVector neighborsVec = mesh.elemNeighbors(neighborsPair[n]);
        SpaMesh::Connectivity::ElemConnVector facesVec = mesh.elemFaces(neighborsPair[n]);
        assert(faceVec.size() == neighborsVec.size());
        unsigned int neighborsNumber = 0;
        for (unsigned int i = 0; i < facesVec.size(); i++)
        {
          if (facesVec[i] == faceIt.value())
          {
	    faceFound++;
            continue;
          }
          facesNumber++;
        }
      }
      assert(faceFound == neighborsPairSize);
      _local_nonzero_blocks[faceIt.value()] += facesNumber;
    }
  }

private:

  //matrix blocks count
  Eigen::VectorXi _local_nonzero_blocks;
  Eigen::VectorXi _nonlocal_nonzero_blocks;
  
  //vector ghost blocks indices
  Eigen::VectorXi _ghost_block_global_indices;

  int _n_local_elements;
  int _n_ghost_elements;
  int _n_global_elements;
  unsigned int _n_vars;
  unsigned int _p_order;
  unsigned int _element_dofs;
  unsigned int _element_dofs_col;
  VarNamesType _var_names;
  VarScalingsType _var_scalings;
};
#endif
