#include <iostream>
#include <vector>
#include <fstream>
#include <map>
#include <stdio.h>
#include <deal.II/grid/grid_generator.h>
#include <deal.II/grid/grid_out.h>
#include <deal.II/grid/manifold_lib.h>

using namespace dealii;
using namespace std;

//auxiliary function for checking whether vertex/line/center locations of active cells of boundary mesh coincide with vertex/line/center locations of corresponding faces of domain cells (all quantities are computed with respect_manifold=true)
template<unsigned int spaceDim> void checkVertexCoodinates(const typename Triangulation<spaceDim-1,spaceDim>::cell_iterator& cellBoundary, const typename Triangulation<spaceDim,spaceDim>::face_iterator& faceDomain){
	Assert( (cellBoundary->has_children()==faceDomain->has_children()), ExcMessage("Boundary cells and corresponding domain faces not refined to the same level!"));
	if(cellBoundary->has_children())
		for(unsigned int i1=0; i1<cellBoundary->n_children(); ++i1)
			checkVertexCoodinates<spaceDim>(cellBoundary->child(i1), faceDomain->child(i1));
	else{
		printf("Boundary cell %4i, Domain face %4i:\n", cellBoundary->index(), faceDomain->index());

		for(unsigned int vertex=0; vertex<GeometryInfo<spaceDim>::vertices_per_face; ++vertex)
			if(spaceDim==2)
				printf("  Vertex %1i: ( %- 1.0e, %- 1.0e )\n", vertex ,faceDomain->vertex(vertex)[0]-cellBoundary->vertex(vertex)[0], faceDomain->vertex(vertex)[1]-cellBoundary->vertex(vertex)[1]);
			else if(spaceDim==3)
				printf("  Vertex %1i: ( %- 1.0e, %- 1.0e, %- 1.0e )\n", vertex ,faceDomain->vertex(vertex)[0]-cellBoundary->vertex(vertex)[0], faceDomain->vertex(vertex)[1]-cellBoundary->vertex(vertex)[1], faceDomain->vertex(vertex)[2]-cellBoundary->vertex(vertex)[2]);
		for(unsigned int line=0; line<GeometryInfo<spaceDim>::lines_per_face; ++line)
			if(spaceDim==2)
				printf("    Line %1i: ( %- 1.0e, %- 1.0e )\n", line ,faceDomain->line(line)->center(true)[0]-cellBoundary->line(line)->center(true)[0], faceDomain->line(line)->center(true)[1]-cellBoundary->line(line)->center(true)[1]);
			else if(spaceDim==3)
				printf("    Line %1i: ( %- 1.0e, %- 1.0e, %- 1.0e )\n", line ,faceDomain->line(line)->center(true)[0]-cellBoundary->line(line)->center(true)[0], faceDomain->line(line)->center(true)[1]-cellBoundary->line(line)->center(true)[1], faceDomain->line(line)->center(true)[2]-cellBoundary->line(line)->center(true)[2]);
		if(spaceDim==2)
				printf("    Center: ( %- 1.0e, %- 1.0e )\n", faceDomain->center(true)[0]-cellBoundary->center(true)[0], faceDomain->center(true)[1]-cellBoundary->center(true)[1]);
		else if(spaceDim==3)
				printf("    Center: ( %- 1.0e, %- 1.0e, %- 1.0e )\n", faceDomain->center(true)[0]-cellBoundary->center(true)[0], faceDomain->center(true)[1]-cellBoundary->center(true)[1], faceDomain->center(true)[2]-cellBoundary->center(true)[2]);

	}
	return;
}

int main(){

	const unsigned int spaceDim=3;			//spatial dimension
	const unsigned int n_refinements=2;		//number of global refinements (the program behaves as expected for spaceDim=2 and arbitrary n_refinement as well as for spaceDim=3 and n_refinements<=1)

	//triangulate ball
	Triangulation<spaceDim, spaceDim> triaDomain;
	GridGenerator::hyper_ball(triaDomain);

	//extract boundary of ball (essentially copied from source code of hyper_sphere(); hyper_sphere() is not used because it doesn't provide with the mapping between boundary cells and corresponding faces of the domain cells)
	//fortunately, the exterior faces of the ball are of type 0, 2, 4 such that the vertex ordering of domain faces coincides with the vertex ordering of the corresponding boundary cells generated by extract_boundary_mesh())
	Triangulation<spaceDim-1, spaceDim> triaBoundary;
	set<types::boundary_id> boundary_ids;
	boundary_ids.insert(0);
	map<Triangulation<spaceDim-1,spaceDim>::cell_iterator, Triangulation<spaceDim,spaceDim>::face_iterator> boundary_to_domain=GridGenerator::extract_boundary_mesh(triaDomain, triaBoundary, boundary_ids);
	triaBoundary.set_all_manifold_ids(0);
	triaBoundary.set_manifold(0, SphericalManifold<spaceDim-1, spaceDim>());

	//refine
	triaDomain.refine_global(n_refinements);
	triaBoundary.refine_global(n_refinements);

	//check whether vertex/line/center locations of active cells of boundary mesh coincide with vertex/line/center locations of corresponding faces of domain cells (all quantities are computed with respect_manifold=true)
	cout << "Differences of coordinates of corresponding vertices/line centers/face centers of domain mesh and boundary mesh (I'd expect all values to be zero within numerical accuracy)" << endl;
	for(typename map<Triangulation<spaceDim-1,spaceDim>::cell_iterator, Triangulation<spaceDim,spaceDim>::face_iterator>::iterator it=boundary_to_domain.begin(); it!=boundary_to_domain.end(); ++it){
		Triangulation<spaceDim-1, spaceDim>::cell_iterator cellBoundary=it->first;
		Triangulation<spaceDim, spaceDim>::face_iterator faceDomain=it->second;
		checkVertexCoodinates<spaceDim>(cellBoundary, faceDomain);
	}

	//output of triangulations (discrepancy can also be observed graphically - just zoom in)
	ofstream outBoundary("triaBoundary.vtk");
	GridOut gridOutBoundary;
	gridOutBoundary.write_vtk(triaBoundary, outBoundary);
	ofstream outDomain("triaDomain.vtk");
	GridOut gridOutDomain;
	gridOutDomain.write_vtk(triaDomain, outDomain);
}