#include "libmesh.h"
#include "mesh.h"
#include "mesh_generation.h"
#include "mesh_refinement.h"
#include "elem_range.h"
#include "error_vector.h"
#include "boundary_info.h"

void print_info(Mesh & mesh)
{
  // print out
  ConstElemRange range(mesh.active_local_elements_begin(), mesh.active_local_elements_end(), 1);
  for (ConstElemRange::const_iterator el = range.begin() ; el != range.end(); ++el)
  {
    const Elem *e = (*el);
    std::cout << "elem id = " << e->id() << ": ";
    
    for (unsigned int i = 0; i < e->n_nodes(); i++)
    {
      const Point & pt = e->point(i);
      std::cout << "  (" << pt(0) << ", " << pt(1) << ", " << pt(2) << ")";
    }
    std::cout << std::endl;

    for (unsigned int side = 0; side < e->n_sides(); side++)
    {
      std::cout << "  side = " << side << ": ";

      std::cout << "  ";
      short int bnd_id = mesh.boundary_info->boundary_id (e, side);
      std::cout << " " << bnd_id << std::endl;
    }
  }
}

int main(int argc, char * argv[])
{
  LibMeshInit init(argc, argv);

  Mesh mesh;

  // generate mesh
  int nx = 2, ny = 2; 
  MeshTools::Generation::build_square(mesh, nx, ny, 0, 1, 0, 1, QUAD4);

  print_info(mesh);  

  // refine me
  std::cout << "----------------------------------------" << std::endl;
    
  MeshRefinement mr(mesh);
  mr.refine_fraction() = 1.0;
  mr.coarsen_fraction() = 0.0;
  mr.max_h_level() = 10;

  ErrorVector ev(mesh.n_local_elem(), 0.5);
  mr.flag_elements_by_error_fraction (ev);
  mr.refine_and_coarsen_elements();

  mesh.boundary_info->build_node_list_from_side_list();
  
  print_info(mesh);  

  return 0;
}