#include <getfem/getfem_mesh.h>
#include <getfem/bgeot_mesh_structure.h>
#include <getfem/getfem_export.h>
#include <getfem/getfem_regular_meshes.h>

int main(int argc, char *argv[]) {
  try {
     getfem::mesh mymesh;

     std::vector<getfem::size_type> nsubdiv(2);
     nsubdiv[0] = 4; nsubdiv[1] = 3;
    getfem::regular_unit_mesh(mymesh, nsubdiv,  
                      bgeot::parallelepiped_geotrans(2,1));
    getfem::vtk_export exp("square.vtk", true);
    exp.exporting(mymesh);
    exp.write_mesh_quality(mymesh);

    if (argc < 2 ) {
      std::cerr << "Usage " << argv[0] << " slice_angle\n";
      exit(-1);
    }
    getfem::stored_mesh_slice sl;
    getfem::base_node x0(0.5,0.5);
    double angle = M_PI*atof(argv[1])/180;
    getfem::base_node n(-sin(angle), cos(angle));
    sl.build(mymesh, getfem::slicer_half_space(x0, n, 0), 1);
    sl.merge_nodes();

    std::cout << sl.nb_simplexes(1) << " " << sl.nb_points() << "\n";
    getfem::vtk_export expsl("sqrare_slice.vtk", true);
    expsl.exporting(sl);
    expsl.write_mesh_quality(mymesh);
  } GMM_STANDARD_CATCH_ERROR;
  return 0;
}