//Include all deal.II header file
#include <deal.II/base/quadrature.h>
#include <deal.II/base/function.h>
#include <deal.II/base/logstream.h>
#include <deal.II/base/timer.h>
#include <deal.II/base/table.h>
#include <deal.II/lac/vector.h>
#include <deal.II/lac/full_matrix.h>
#include <deal.II/lac/constraint_matrix.h>
#include <deal.II/fe/fe_system.h>
#include <deal.II/fe/fe_q.h>
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/mapping_q1.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/dofs/dof_tools.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/numerics/data_out.h>
#include <deal.II/numerics/vector_tools.h>
#include <deal.II/base/tensor_function.h>
#include <deal.II/base/point.h>
#include <deal.II/base/conditional_ostream.h>
#include <deal.II/base/utilities.h>
#include <deal.II/lac/lapack_full_matrix.h>
#include <deal.II/numerics/matrix_tools.h>
#include <deal.II/numerics/error_estimator.h>
#include <deal.II/base/parameter_handler.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/grid/tria_boundary_lib.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/grid_in.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/dofs/dof_renumbering.h>
#include <deal.II/dofs/dof_accessor.h>
#include <deal.II/lac/solver_cg.h>
#include <deal.II/lac/solver_gmres.h>
#include <deal.II/lac/solver_bicgstab.h>
#include <deal.II/lac/precondition.h>
#include <deal.II/distributed/tria.h>
#include <deal.II/lac/parallel_vector.h>
#include <deal.II/matrix_free/matrix_free.h>
#include <deal.II/matrix_free/fe_evaluation.h>
#include <deal.II/lac/slepc_solver.h>
#include <deal.II/base/config.h>
#include <deal.II/base/smartpointer.h>
//Include generic C++ headers
#include <fstream>
#include <iostream>


using namespace dealii;
int main (int argc, char *argv[])
{
  Utilities::MPI::MPI_InitFinalize mpi_initialization(argc, argv);

  const double L=20;
  Triangulation<3,3> triangulation;  
  GridGenerator::hyper_cube (triangulation, -L, L);

  const double Ltemp=L;//FIXME: Size of constraint matrix changes between Ltemp=L, and Ltemp=-L;
  //mark faces
  typename Triangulation<3,3>::active_cell_iterator cell = triangulation.begin_active(), endc = triangulation.end();
  for(; cell!=endc; ++cell) 
  {
      for(unsigned int f=0; f < GeometryInfo<3>::faces_per_cell; ++f)
	{
	  const Point<3> face_center = cell->face(f)->center();
	  if(cell->face(f)->at_boundary())
	    {
	      if (std::abs(face_center[0]-Ltemp)<1.0e-5)
		cell->face(f)->set_boundary_id(1);
	      else if (std::abs(face_center[0]+Ltemp)<1.0e-5)
		cell->face(f)->set_boundary_id(2);
	      else if (std::abs(face_center[1]-Ltemp)<1.0e-5)
		cell->face(f)->set_boundary_id(3);
	      else if (std::abs(face_center[1]+Ltemp)<1.0e-5)
		cell->face(f)->set_boundary_id(4);
	      else if (std::abs(face_center[2]-Ltemp)<1.0e-5)
		cell->face(f)->set_boundary_id(5);
	      else if (std::abs(face_center[2]+Ltemp)<1.0e-5)
		cell->face(f)->set_boundary_id(6);
	    }
	}
  }
 
  std::vector<GridTools::PeriodicFacePair<typename Triangulation<3,3>::cell_iterator> > periodicity_vector;
  for (int i = 0; i < 3; ++i)
  {
      GridTools::collect_periodic_faces(triangulation, /*b_id1*/ 2*i+1, /*b_id2*/ 2*i+2,/*direction*/ i, periodicity_vector);
  }
  triangulation.add_periodicity(periodicity_vector);

  //two level mesh refinement 
  
  triangulation.refine_global(1);

  typename Triangulation<3,3>::active_cell_iterator cellBegin = triangulation.begin_active();
  cellBegin->set_refine_flag();
  triangulation.execute_coarsening_and_refinement(); 
  
  std::cout << "number of elements: "
	<< triangulation.n_global_active_cells()
	<< std::endl;   
  /////
  
  FESystem<3> FE(FE_Q<3>(QGaussLobatto<1>(2)), 1); //linear shape function
  DoFHandler<3> dofHandler (triangulation);
  dofHandler.distribute_dofs(FE);

  ///creating constraint matrix
  ConstraintMatrix constraints;
  constraints.clear();
  std::cout<< "Adding hanging node constraints... "<< std::endl;
  DoFTools::make_hanging_node_constraints(dofHandler, constraints);
  std::cout<< "Adding periodicity constraints... "<< std::endl;
  std::vector<GridTools::PeriodicFacePair<typename DoFHandler<3>::cell_iterator> > periodicity_vectorDof;
  for (int i = 0; i < 3; ++i)
    {
      GridTools::collect_periodic_faces(dofHandler, /*b_id1*/ 2*i+1, /*b_id2*/ 2*i+2,/*direction*/ i, periodicity_vectorDof);
    }
  DoFTools::make_periodicity_constraints<DoFHandler<3> >(periodicity_vectorDof, constraints);
  constraints.close();
  std::cout<< "size of constraint matrix: "<< constraints.n_constraints()<<std::endl;
  std::cout<< "printing constraint matrix ...."<<std::endl;
  constraints.print(std::cout);

  /////write mesh
  std::cout<<"writing mesh .." <<std::endl;
  DataOut<3> data_out;
  data_out.attach_dof_handler(dofHandler);
  data_out.build_patches ();
  std::ofstream output("mesh.vtu");
  data_out.write_vtu (output);  
}