testing.h

#ifndef TETHEX_TESTING_H
#define TETHEX_TESTING_H

#include "config.h"

#if defined(TESTING)
#include "gtest/gtest.h"
#include "tethex.h"
#include <fstream>

using namespace tethex;


//---------------------------------------------------------
//
// Common auxiliary functions
//
//---------------------------------------------------------

void check_vertices(const Mesh &mesh,
                    const std::vector<Point> &vertices)
{
  EXPECT_EQ(mesh.get_n_vertices(), vertices.size());
  for (unsigned int i = 0; i < vertices.size(); ++i)
    for (unsigned int j = 0; j < Point::n_coord; ++j)
      EXPECT_DOUBLE_EQ(mesh.get_vertex(i).get_coord(j), vertices[i].get_coord(j));
}



void check_points(const Mesh &mesh,
                  const std::vector<PhysPoint> &points)
{
  EXPECT_EQ(mesh.get_n_points(), points.size());
  for (unsigned int i = 0; i < points.size(); ++i)
    for (unsigned int j = 0; j < PhysPoint::n_vertices; ++j)
      EXPECT_EQ(mesh.get_point(i).get_vertex(j), points[i].get_vertex(j));
}



void check_edges(const Mesh &mesh,
                 const std::vector<Line> &edges)
{
  EXPECT_EQ(mesh.get_n_edges(), edges.size());
  for (unsigned int i = 0; i < edges.size(); ++i)
    for (unsigned int j = 0; j < Line::n_vertices; ++j)
      EXPECT_EQ(mesh.get_edge(i).get_vertex(j), edges[i].get_vertex(j));
}



void check_lines(const Mesh &mesh,
                 const std::vector<Line> &lines)
{
  EXPECT_EQ(mesh.get_n_lines(), lines.size());
  for (unsigned int i = 0; i < lines.size(); ++i)
    for (unsigned int j = 0; j < Line::n_vertices; ++j)
      EXPECT_EQ(mesh.get_line(i).get_vertex(j), lines[i].get_vertex(j));
}



void check_faces(const Mesh &mesh,
                 const std::vector<Triangle> &faces)
{
  EXPECT_EQ(mesh.get_n_faces(), faces.size());
  for (unsigned int i = 0; i < faces.size(); ++i)
    for (unsigned int j = 0; j < Triangle::n_vertices; ++j)
    EXPECT_EQ(mesh.get_face(i).get_vertex(j), faces[i].get_vertex(j));
}



void check_triangle_vertices(const Mesh &mesh,
                             const std::vector<Triangle> &triangles)
{
  EXPECT_EQ(mesh.get_n_triangles(), triangles.size());
  for (unsigned int i = 0; i < triangles.size(); ++i)
    for (unsigned int j = 0; j < Triangle::n_vertices; ++j)
    EXPECT_EQ(mesh.get_triangle(i).get_vertex(j), triangles[i].get_vertex(j));
}



void check_triangle_edges(const Mesh &mesh,
                          const std::vector<Triangle> &triangles)
{
  EXPECT_EQ(mesh.get_n_triangles(), triangles.size());
  for (unsigned int i = 0; i < triangles.size(); ++i)
    for (unsigned int j = 0; j < Triangle::n_edges; ++j)
    EXPECT_EQ(mesh.get_triangle(i).get_edge(j), triangles[i].get_edge(j));
}



void check_quadrangle_vertices(const Mesh &mesh,
                               const std::vector<Quadrangle> &quadrangles)
{
  EXPECT_EQ(mesh.get_n_quadrangles(), quadrangles.size());
  for (unsigned int i = 0; i < quadrangles.size(); ++i)
    for (unsigned int j = 0; j < Quadrangle::n_vertices; ++j)
    EXPECT_EQ(mesh.get_quadrangle(i).get_vertex(j), quadrangles[i].get_vertex(j));
}



void check_tetrahedron_vertices(const Mesh &mesh,
                               const std::vector<Tetrahedron> &tetrahedra)
{
  EXPECT_EQ(mesh.get_n_tetrahedra(), tetrahedra.size());
  for (unsigned int i = 0; i < tetrahedra.size(); ++i)
    for (unsigned int j = 0; j < Tetrahedron::n_vertices; ++j)
    EXPECT_EQ(mesh.get_tetrahedron(i).get_vertex(j), tetrahedra[i].get_vertex(j));
}



void check_tetrahedron_edges(const Mesh &mesh,
                             const std::vector<Tetrahedron> &tetrahedra)
{
  EXPECT_EQ(mesh.get_n_tetrahedra(), tetrahedra.size());
  for (unsigned int i = 0; i < tetrahedra.size(); ++i)
    for (unsigned int j = 0; j < Tetrahedron::n_edges; ++j)
    EXPECT_EQ(mesh.get_tetrahedron(i).get_edge(j), tetrahedra[i].get_edge(j));
}



void check_tetrahedron_faces(const Mesh &mesh,
                             const std::vector<Tetrahedron> &tetrahedra)
{
  EXPECT_EQ(mesh.get_n_tetrahedra(), tetrahedra.size());
  for (unsigned int i = 0; i < tetrahedra.size(); ++i)
    for (unsigned int j = 0; j < Tetrahedron::n_faces; ++j)
    EXPECT_EQ(mesh.get_tetrahedron(i).get_face(j), tetrahedra[i].get_face(j));
}



void check_hexahedron_vertices(const Mesh &mesh,
                               const std::vector<Hexahedron> &hexahedra)
{
  EXPECT_EQ(mesh.get_n_hexahedra(), hexahedra.size());
  for (unsigned int i = 0; i < hexahedra.size(); ++i)
    for (unsigned int j = 0; j < Hexahedron::n_vertices; ++j)
    EXPECT_EQ(mesh.get_hexahedron(i).get_vertex(j), hexahedra[i].get_vertex(j));
}





//---------------------------------------------------------
//
// Common procedure 2D
//
//---------------------------------------------------------

void test_2D(const std::string &input_file,
             const std::string &output_file,
             std::vector<Point> &vertices,
             const std::vector<Point> &new_vertices,
             const std::vector<PhysPoint> points,
             const std::vector<Line> &edges,
             const std::vector<Line> &lines,
             const std::vector<Line> &new_lines,
             const std::vector<Triangle> &triangles,
             const std::vector<Quadrangle> &quadrangles)
{
  // read the mesh
  Mesh mesh;
  try
  {
    mesh.read(input_file);
  }
  catch (std::runtime_error &e)
  {
    std::cout << "\n=== EXCEPTION === \n";
    std::cout << e.what() << std::endl << std::endl;;
  }

  // check that mesh is right
  if (!vertices.empty())
    check_vertices(mesh, vertices);

  if (!lines.empty())
    check_lines(mesh, lines);

  if (!triangles.empty())
    check_triangle_vertices(mesh, triangles);

  // convert triangles into quadrangles
  try
  {
    mesh.convert();
  }
  catch (std::runtime_error &e)
  {
    std::cout << "\n=== EXCEPTION === \n";
    std::cout << e.what() << std::endl << std::endl;;
  }

  // check new vertices
  if (!new_vertices.empty())
  {
    vertices.insert(vertices.end(), new_vertices.begin(), new_vertices.end());
    check_vertices(mesh, vertices);
  }

  // during conversion we generate edges
  if (!edges.empty())
    check_edges(mesh, edges);

  // if we don't delete simplices,
  // we can compare edges of triangles
#if !defined(DELETE_SIMPLICES)
  if (!triangles.empty())
    check_triangle_edges(mesh, triangles);
#endif

  // now we compare quadrangles...
  if (!quadrangles.empty())
    check_quadrangle_vertices(mesh, quadrangles);

  // ... boundary lines ...
  if (!new_lines.empty())
    check_lines(mesh, new_lines);

  // ... and points
  if (!points.empty())
    check_points(mesh, points);

  // write the results for double checking
  if (!output_file.empty())
  {
    try
    {
      mesh.write(output_file);
    }
    catch (std::runtime_error &e)
    {
      std::cout << "\n=== EXCEPTION === \n";
      std::cout << e.what() << std::endl << std::endl;;
    }
  }
}



//---------------------------------------------------------
//
// Common procedure 3D
//
//---------------------------------------------------------

void test_3D(const std::string &input_file,
             const std::string &output_file,
             std::vector<Point> &vertices,
             const std::vector<Point> &new_vertices,
             const std::vector<PhysPoint> points,
             const std::vector<Line> &edges,
             const std::vector<Line> &lines,
             const std::vector<Line> &new_lines,
             const std::vector<Triangle> &faces,
             const std::vector<Triangle> &triangles,
             const std::vector<Quadrangle> &quadrangles,
             const std::vector<Tetrahedron> &tetrahedra,
             const std::vector<Hexahedron> &hexahedra)
{
  // read the mesh
  Mesh mesh;
  try
  {
    mesh.read(input_file);
  }
  catch (std::runtime_error &e)
  {
    std::cout << "\n=== EXCEPTION === \n";
    std::cout << e.what() << std::endl << std::endl;;
  }

  // check that mesh is right
  if (!vertices.empty())
    check_vertices(mesh, vertices);

  if (!lines.empty())
    check_lines(mesh, lines);

  if (!triangles.empty())
    check_triangle_vertices(mesh, triangles);

  if (!tetrahedra.empty())
    check_tetrahedron_vertices(mesh, tetrahedra);

  // convert tetrahedra into hexahedra
  try
  {
    mesh.convert();
  }
  catch (std::runtime_error &e)
  {
    std::cout << "\n=== EXCEPTION === \n";
    std::cout << e.what() << std::endl << std::endl;;
  }

  // check new vertices
  if (!new_vertices.empty())
  {
    vertices.insert(vertices.end(), new_vertices.begin(), new_vertices.end());
    check_vertices(mesh, vertices);
  }

  // during conversion we generate edges
  if (!edges.empty())
    check_edges(mesh, edges);

  // and faces
  if (!faces.empty())
    check_faces(mesh, faces);

  // if we don't delete simplices,
  // we can compare edges and faces of tetrahedra and
  // edges of boundary triangles
#if !defined(DELETE_SIMPLICES)
  if (!tetrahedra.empty())
  {
    check_tetrahedron_edges(mesh, tetrahedra);
    check_tetrahedron_faces(mesh, tetrahedra);
  }
#endif

  // now we compare hexahedra...
  if (!hexahedra.empty())
    check_hexahedron_vertices(mesh, hexahedra);

  // ... boundary quadrangles...
  if (!quadrangles.empty())
    check_quadrangle_vertices(mesh, quadrangles);

  // ... boundary lines ...
  if (!new_lines.empty())
    check_lines(mesh, new_lines);

  // ... and points
  if (!points.empty())
    check_points(mesh, points);

  // write the results for double checking
  if (!output_file.empty())
  {
    try
    {
      mesh.write(output_file);
    }
    catch (std::runtime_error &e)
    {
      std::cout << "\n=== EXCEPTION === \n";
      std::cout << e.what() << std::endl << std::endl;;
    }
  }
}





//---------------------------------------------------------
//
// Testing meshes declarations
//
//---------------------------------------------------------

void mesh_2d_1(std::string &input_file, std::string &output_file);
void mesh_2d_2(std::string &input_file, std::string &output_file);
void mesh_3d_1(std::string &input_file, std::string &output_file);
void mesh_3d_2(std::string &input_file, std::string &output_file);




//---------------------------------------------------------
//
// Tests 2D
//
//---------------------------------------------------------

TEST(Mesh_2D, OneElement)
{
  std::string input_file, output_file;
  mesh_2d_1(input_file, output_file);

  // vertices and mesh elements from the file
  std::vector<Point> vertices;
  vertices.push_back(Point(0, 0, 0));
  vertices.push_back(Point(1, 0, 0));
  vertices.push_back(Point(0, 1, 0));

  std::vector<Line> edges;
  edges.push_back(Line(0, 1, 11));
  edges.push_back(Line(0, 2, 11));
  edges.push_back(Line(1, 2, 11));

  std::vector<Line> lines;
  std::vector<Line> new_lines;

  std::vector<Triangle> triangles;
  triangles.push_back(Triangle(0, 1, 2, 11));
  triangles[0].set_edge(0, 0);
  triangles[0].set_edge(1, 1);
  triangles[0].set_edge(2, 2);

  std::vector<Quadrangle> quadrangles;
  quadrangles.push_back(Quadrangle(0, 3, 6, 4, 11));
  quadrangles.push_back(Quadrangle(1, 5, 6, 3, 11));
  quadrangles.push_back(Quadrangle(2, 4, 6, 5, 11));

  // new vertices created after conversion from triangles to quadrangles
  std::vector<Point> new_vertices;
  new_vertices.push_back(Point(1./2., 0, 0));
  new_vertices.push_back(Point(0, 1./2., 0));
  new_vertices.push_back(Point(1./2., 1./2., 0));
  new_vertices.push_back(Point(1./3., 1./3., 0));

  // physical points
  std::vector<PhysPoint> points;
  points.push_back(PhysPoint(0, 41));
  points.push_back(PhysPoint(2, 42));

  // testing
  test_2D(input_file,
          output_file,
          vertices,
          new_vertices,
          points,
          edges,
          lines,
          new_lines,
          triangles,
          quadrangles);

  // once again for converted mesh
  new_vertices.clear(); // vertices already include new_vertices
  edges.clear(); // we don't generate edges for quadrangles
  lines.clear(); // there are no original lines
#if defined(DELETE_SIMPLICES)
  triangles.clear(); // there are no triangles anymore
#endif
  test_2D(output_file,
          "",
          vertices,
          new_vertices,
          points,
          edges,
          lines,
          new_lines,
          triangles,
          quadrangles);
}




TEST(Mesh_2D, TwoElements)
{
  std::string input_file, output_file;
  mesh_2d_2(input_file, output_file);

  // vertices and mesh elements from the file
  std::vector<Point> vertices;
  vertices.push_back(Point(0, 0, 0));
  vertices.push_back(Point(1, 0, 0));
  vertices.push_back(Point(0, 1, 0));
  vertices.push_back(Point(1, 1, 0));

  std::vector<Line> edges;
  edges.push_back(Line(0, 1, 11));
  edges.push_back(Line(0, 2, 11));
  edges.push_back(Line(1, 2, 11));
  edges.push_back(Line(1, 3, 11));
  edges.push_back(Line(2, 3, 11));

  std::vector<Line> lines;
  lines.push_back(Line(0, 1, 21));
  lines.push_back(Line(0, 2, 22));
  lines.push_back(Line(1, 3, 23));
  lines.push_back(Line(2, 3, 24));

  std::vector<Line> new_lines;
  new_lines.push_back(Line(0, 4, 21));
  new_lines.push_back(Line(0, 5, 22));
  new_lines.push_back(Line(1, 7, 23));
  new_lines.push_back(Line(2, 8, 24));
  new_lines.push_back(Line(4, 1, 21));
  new_lines.push_back(Line(5, 2, 22));
  new_lines.push_back(Line(7, 3, 23));
  new_lines.push_back(Line(8, 3, 24));

  std::vector<Triangle> triangles;
  triangles.push_back(Triangle(0, 1, 2, 11));
  triangles[0].set_edge(0, 0);
  triangles[0].set_edge(1, 1);
  triangles[0].set_edge(2, 2);

  triangles.push_back(Triangle(1, 3, 2, 12));
  triangles[1].set_edge(0, 3);
  triangles[1].set_edge(1, 4);
  triangles[1].set_edge(2, 2);

  std::vector<Quadrangle> quadrangles;
  quadrangles.push_back(Quadrangle(0, 4, 9, 5, 11));
  quadrangles.push_back(Quadrangle(1, 6, 9, 4, 11));
  quadrangles.push_back(Quadrangle(2, 5, 9, 6, 11));
  quadrangles.push_back(Quadrangle(1, 7, 10, 6, 12));
  quadrangles.push_back(Quadrangle(3, 8, 10, 7, 12));
  quadrangles.push_back(Quadrangle(2, 6, 10, 8, 12));

  // new vertices created after conversion from triangles to quadrangles
  std::vector<Point> new_vertices;
  new_vertices.push_back(Point(1./2., 0, 0));
  new_vertices.push_back(Point(0, 1./2., 0));
  new_vertices.push_back(Point(1./2., 1./2., 0));
  new_vertices.push_back(Point(1, 1./2., 0));
  new_vertices.push_back(Point(1./2., 1, 0));
  new_vertices.push_back(Point(1./3., 1./3., 0));
  new_vertices.push_back(Point(2./3., 2./3., 0));

  // physical points
  std::vector<PhysPoint> points;

  // testing
  test_2D(input_file,
          output_file,
          vertices,
          new_vertices,
          points,
          edges,
          lines,
          new_lines,
          triangles,
          quadrangles);

  // once again for converted mesh
  new_vertices.clear(); // vertices already include new_vertices
  edges.clear(); // we don't generate edges for quadrangles
  lines.clear(); // there are no original lines
#if defined(DELETE_SIMPLICES)
  triangles.clear(); // there are no triangles anymore
#endif
  test_2D(output_file,
          "",
          vertices,
          new_vertices,
          points,
          edges,
          lines,
          new_lines,
          triangles,
          quadrangles);
}





//---------------------------------------------------------
//
// Tests 3D
//
//---------------------------------------------------------


TEST(Mesh_3D, OneElement)
{
  std::string input_file, output_file;
  mesh_3d_1(input_file, output_file);

  // vertices and mesh elements from the file
  std::vector<Point> vertices;
  vertices.push_back(Point(0, 0, 0));
  vertices.push_back(Point(1, 0, 0));
  vertices.push_back(Point(0, 1, 0));
  vertices.push_back(Point(0, 0, 1));

  // boundary elements are empty
  std::vector<Line> lines;
  std::vector<Line> new_lines;
  std::vector<Triangle> triangles;
  std::vector<Quadrangle> quadrangles;

  // mesh elements
  std::vector<Line> edges;
  edges.push_back(Line(0, 1, 11));
  edges.push_back(Line(0, 2, 11));
  edges.push_back(Line(1, 2, 11));
  edges.push_back(Line(0, 3, 11));
  edges.push_back(Line(1, 3, 11));
  edges.push_back(Line(2, 3, 11));

  std::vector<Triangle> faces;
  faces.push_back(Triangle(0, 1, 2, 11));
  faces.push_back(Triangle(0, 1, 3, 11));
  faces.push_back(Triangle(0, 2, 3, 11));
  faces.push_back(Triangle(1, 2, 3, 11));

  std::vector<Tetrahedron> tetrahedra;
  tetrahedra.push_back(Tetrahedron(0, 1, 2, 3, 11));
  tetrahedra[0].set_edge(0, 0);
  tetrahedra[0].set_edge(1, 1);
  tetrahedra[0].set_edge(2, 2);
  tetrahedra[0].set_edge(3, 3);
  tetrahedra[0].set_edge(4, 4);
  tetrahedra[0].set_edge(5, 5);
  tetrahedra[0].set_face(0, 0);
  tetrahedra[0].set_face(1, 1);
  tetrahedra[0].set_face(2, 2);
  tetrahedra[0].set_face(3, 3);

  // new vertices created after conversion from triangles to quadrangles
  std::vector<Point> new_vertices;
  new_vertices.push_back(Point(1./2., 0, 0));
  new_vertices.push_back(Point(0, 1./2., 0));
  new_vertices.push_back(Point(1./2., 1./2., 0));
  new_vertices.push_back(Point(0, 0, 1./2.));
  new_vertices.push_back(Point(1./2., 0, 1./2.));
  new_vertices.push_back(Point(0, 1./2., 1./2.));
  new_vertices.push_back(Point(1./3., 1./3., 0));
  new_vertices.push_back(Point(1./3., 0, 1./3.));
  new_vertices.push_back(Point(0, 1./3., 1./3.));
  new_vertices.push_back(Point(1./3., 1./3., 1./3.));
  new_vertices.push_back(Point(1./4., 1./4., 1./4.));

  std::vector<Hexahedron> hexahedra;
  hexahedra.push_back(Hexahedron(7, 11, 14, 12, 0, 4, 10, 5, 11));
  hexahedra.push_back(Hexahedron(1, 4, 10, 6, 8, 11, 14, 13, 11));
  hexahedra.push_back(Hexahedron(9, 12, 14, 13, 2, 5, 10, 6, 11));
  hexahedra.push_back(Hexahedron(3, 7, 11, 8, 9, 12, 14, 13, 11));

  // physical points
  std::vector<PhysPoint> points;
  points.push_back(PhysPoint(0, 41));
  points.push_back(PhysPoint(2, 42));
  points.push_back(PhysPoint(3, 43));

  // testing
  test_3D(input_file,
          output_file,
          vertices,
          new_vertices,
          points,
          edges,
          lines,
          new_lines,
          faces,
          triangles,
          quadrangles,
          tetrahedra,
          hexahedra);

  // once again for converted mesh
  new_vertices.clear(); // vertices already include new_vertices
  edges.clear(); // we don't generate edges for tetrahedra
  lines.clear(); // there are no original lines
  faces.clear(); // we don't generate faces for tetrahedra
#if defined(DELETE_SIMPLICES)
  triangles.clear(); // there are no triangles anymore
  tetrahedra.clear(); // and tetrahedra as well
#endif
  test_3D(output_file,
          "",
          vertices,
          new_vertices,
          points,
          edges,
          lines,
          new_lines,
          faces,
          triangles,
          quadrangles,
          tetrahedra,
          hexahedra);
}



TEST(Mesh_3D, TwoElements)
{
  std::string input_file, output_file;
  mesh_3d_2(input_file, output_file);

  // vertices and mesh elements from the file
  std::vector<Point> vertices;
  vertices.push_back(Point(0, 0, 0));
  vertices.push_back(Point(1, 0, 0));
  vertices.push_back(Point(0, 1, 0));
  vertices.push_back(Point(0, 0, 1));
  vertices.push_back(Point(0, 0, -1));

  std::vector<Line> edges;
  edges.push_back(Line(0, 1, 11));
  edges.push_back(Line(0, 2, 11));
  edges.push_back(Line(1, 2, 11));
  edges.push_back(Line(0, 3, 11));
  edges.push_back(Line(1, 3, 11));
  edges.push_back(Line(2, 3, 11));
  edges.push_back(Line(0, 4, 12));
  edges.push_back(Line(1, 4, 12));
  edges.push_back(Line(2, 4, 12));

  std::vector<Line> lines;
  lines.push_back(Line(0, 3, 31));
  lines.push_back(Line(4, 0, 32));

  std::vector<Line> new_lines;
  new_lines.push_back(Line(0, 8, 31));
  new_lines.push_back(Line(4, 11, 32));
  new_lines.push_back(Line(8, 3, 31));
  new_lines.push_back(Line(11, 0, 32));

  std::vector<Triangle> triangles;
  triangles.push_back(Triangle(0, 1, 3, 21));
  triangles.push_back(Triangle(4, 2, 0, 22));

  std::vector<Quadrangle> quadrangles;
  quadrangles.push_back(Quadrangle(0, 5, 15, 8, 21));
  quadrangles.push_back(Quadrangle(1, 5, 15, 9, 21));
  quadrangles.push_back(Quadrangle(3, 8, 15, 9, 21));
  quadrangles.push_back(Quadrangle(4, 13, 19, 11, 22));
  quadrangles.push_back(Quadrangle(2, 13, 19, 6, 22));
  quadrangles.push_back(Quadrangle(0, 11, 19, 6, 22));

  // new vertices created after conversion from triangles to quadrangles
  std::vector<Point> new_vertices;
  new_vertices.push_back(Point(1./2., 0, 0));
  new_vertices.push_back(Point(0, 1./2., 0));
  new_vertices.push_back(Point(1./2., 1./2., 0));
  new_vertices.push_back(Point(0, 0, 1./2.));
  new_vertices.push_back(Point(1./2., 0, 1./2.));
  new_vertices.push_back(Point(0, 1./2., 1./2.));
  new_vertices.push_back(Point(0, 0, -1./2.));
  new_vertices.push_back(Point(1./2., 0, -1./2.));
  new_vertices.push_back(Point(0, 1./2., -1./2.));
  new_vertices.push_back(Point(1./3., 1./3., 0));
  new_vertices.push_back(Point(1./3., 0, 1./3.));
  new_vertices.push_back(Point(0, 1./3., 1./3.));
  new_vertices.push_back(Point(1./3., 1./3., 1./3.));
  new_vertices.push_back(Point(1./3., 0, -1./3.));
  new_vertices.push_back(Point(0, 1./3., -1./3.));
  new_vertices.push_back(Point(1./3., 1./3., -1./3.));
  new_vertices.push_back(Point(1./4., 1./4., 1./4.));
  new_vertices.push_back(Point(1./4., 1./4., -1./4.));

  std::vector<Triangle> faces;
  faces.push_back(Triangle(0, 1, 2, 11));
  faces.push_back(Triangle(0, 1, 3, 11));
  faces.push_back(Triangle(0, 2, 3, 11));
  faces.push_back(Triangle(1, 2, 3, 11));
  faces.push_back(Triangle(0, 1, 4, 12));
  faces.push_back(Triangle(0, 2, 4, 12));
  faces.push_back(Triangle(1, 2, 4, 12));

  std::vector<Tetrahedron> tetrahedra;
  tetrahedra.push_back(Tetrahedron(0, 1, 2, 3, 11));
  tetrahedra[0].set_edge(0, 0);
  tetrahedra[0].set_edge(1, 1);
  tetrahedra[0].set_edge(2, 2);
  tetrahedra[0].set_edge(3, 3);
  tetrahedra[0].set_edge(4, 4);
  tetrahedra[0].set_edge(5, 5);
  tetrahedra[0].set_face(0, 0);
  tetrahedra[0].set_face(1, 1);
  tetrahedra[0].set_face(2, 2);
  tetrahedra[0].set_face(3, 3);

  tetrahedra.push_back(Tetrahedron(0, 1, 2, 4, 12));
  tetrahedra[1].set_edge(0, 0);
  tetrahedra[1].set_edge(1, 1);
  tetrahedra[1].set_edge(2, 2);
  tetrahedra[1].set_edge(3, 6);
  tetrahedra[1].set_edge(4, 7);
  tetrahedra[1].set_edge(5, 8);
  tetrahedra[1].set_face(0, 0);
  tetrahedra[1].set_face(1, 4);
  tetrahedra[1].set_face(2, 5);
  tetrahedra[1].set_face(3, 6);

  std::vector<Hexahedron> hexahedra;

  // physical points
  std::vector<PhysPoint> points;

  // testing
  test_3D(input_file,
          output_file,
          vertices,
          new_vertices,
          points,
          edges,
          lines,
          new_lines,
          faces,
          triangles,
          quadrangles,
          tetrahedra,
          hexahedra);

  // once again for converted mesh
  new_vertices.clear(); // vertices already include new_vertices
  edges.clear(); // we don't generate edges for tetrahedra
  lines.clear(); // there are no original lines
  faces.clear(); // we don't generate faces for tetrahedra
#if defined(DELETE_SIMPLICES)
  triangles.clear(); // there are no triangles anymore
  tetrahedra.clear(); // and tetrahedra as well
#endif
  test_3D(output_file,
          "",
          vertices,
          new_vertices,
          points,
          edges,
          lines,
          new_lines,
          faces,
          triangles,
          quadrangles,
          tetrahedra,
          hexahedra);
}



//---------------------------------------------------------
//
// Testing meshes definitions
//
//---------------------------------------------------------

void mesh_2d_1(std::string &input_file, std::string &output_file)
{
  input_file = "mesh_2d_1.msh";
  output_file = "mesh_2d_1_hex.msh";
  std::ofstream out(input_file.c_str());
  require(out, "File " + input_file + " cannot be opened!");
  out << "$MeshFormat\n2.2 0 8\n$EndMeshFormat\n$Nodes\n";
  out << "3\n";
  out << "1 0 0 0\n";
  out << "2 1 0 0\n";
  out << "3 0 1 0\n";
  out << "$EndNodes\n$Elements\n";
  out << "3\n";
  out << "1 2 2 11 0 1 2 3\n";
  out << "2 15 2 41 0 1\n";
  out << "3 15 2 42 0 3\n";
  out << "$EndElements\n";
  out.close();
}




void mesh_2d_2(std::string &input_file, std::string &output_file)
{
  input_file = "mesh_2d_2.msh";
  output_file = "mesh_2d_2_hex.msh";
  std::ofstream out(input_file.c_str());
  require(out, "File " + input_file + " cannot be opened!");
  out << "$MeshFormat\n2.2 0 8\n$EndMeshFormat\n$Nodes\n";
  out << "4\n";
  out << "1 0 0 0\n";
  out << "2 1 0 0\n";
  out << "3 0 1 0\n";
  out << "4 1 1 0\n";
  out << "$EndNodes\n$Elements\n";
  out << "6\n";
  out << "1 2 2 11 0 1 2 3\n";
  out << "2 2 2 12 0 2 4 3\n";
  out << "3 1 2 21 0 1 2\n";
  out << "4 1 2 22 0 1 3\n";
  out << "5 1 2 23 0 2 4\n";
  out << "6 1 2 24 0 3 4\n";
  //out << "3 15 2 21 0 1\n"; // physical point
  out << "$EndElements\n";
  out.close();
}




void mesh_3d_1(std::string &input_file, std::string &output_file)
{
  input_file = "mesh_3d_1.msh";
  output_file = "mesh_3d_1_hex.msh";
  std::ofstream out(input_file.c_str());
  require(out, "File " + input_file + " cannot be opened!");
  out << "$MeshFormat\n2.2 0 8\n$EndMeshFormat\n$Nodes\n";
  out << "4\n";
  out << "1 0 0 0\n";
  out << "2 1 0 0\n";
  out << "3 0 1 0\n";
  out << "4 0 0 1\n";
  out << "$EndNodes\n$Elements\n";
  out << "4\n";
  out << "1 4 2 11 0 1 2 3 4\n";
  out << "2 15 2 41 0 1\n";
  out << "3 15 2 42 0 3\n";
  out << "4 15 2 43 0 4\n";
  out << "$EndElements\n";
  out.close();
}




void mesh_3d_2(std::string &input_file, std::string &output_file)
{
  input_file = "mesh_3d_2.msh";
  output_file = "mesh_3d_2_hex.msh";
  std::ofstream out(input_file.c_str());
  require(out, "File " + input_file + " cannot be opened!");
  out << "$MeshFormat\n2.2 0 8\n$EndMeshFormat\n$Nodes\n";
  out << "5\n";
  out << "1 0 0 0\n";
  out << "2 1 0 0\n";
  out << "3 0 1 0\n";
  out << "4 0 0 1\n";
  out << "5 0 0 -1\n";
  out << "$EndNodes\n$Elements\n";
  out << "6\n";
  out << "1 4 2 11 0 1 2 3 4\n";
  out << "2 4 2 12 0 1 2 3 5\n";
  out << "3 2 2 21 0 1 2 4\n";
  out << "4 2 2 22 0 5 3 1\n";
  out << "5 1 2 31 0 1 4\n";
  out << "6 1 2 32 0 5 1\n";
  out << "$EndElements\n";
  out.close();
}

#endif // TESTING

#endif // TETHEX_TESTING_H