#include <deal.II/grid/tria.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_refinement.h>
#include <deal.II/grid/tria_boundary_lib.h>
#include <deal.II/grid/grid_out.h>


// This again is C++:
#include <fstream>
#include <iostream>
using namespace dealii;
using namespace std;


int main(){
    const int dim=3;

    Triangulation<dim>   triangulation;
    Triangulation<dim-1,dim> triangulation_surface;
    //GridGenerator::hyper_cube(triangulation,-100,100,true);

    GridGenerator::cylinder(triangulation,100,200);

    static const CylinderBoundary<dim> outer_cylinder (100,0);

    triangulation.set_boundary(0,outer_cylinder);
    triangulation.refine_global (2);

// just output:
    {
        GridOut grid_out;
        ofstream file("gnu_grid.gpl");
        grid_out.write_gnuplot(triangulation,file);
        file.close();
    }

    static const CylinderBoundary<dim-1,dim> surface_cyl(100,0);
    triangulation_surface.set_boundary(0,surface_cyl);

    GridGenerator::extract_boundary_mesh(triangulation,triangulation_surface);

    {
        GridOut grid_out;
        ofstream file("gnu_surface.gpl");
        grid_out.write_gnuplot(triangulation_surface,file);
        file.close();
    }



    std::cout<<triangulation_surface.n_used_vertices()<<std::endl;
    std::cout<<triangulation_surface.n_active_cells(2)<<std::endl;

}