/**
 * \file
 * \copyright
 * Copyright (c) 2012-2020, OpenGeoSys Community (http://www.opengeosys.org)
 *            Distributed under a Modified BSD License.
 *              See accompanying file LICENSE.txt or
 *              http://www.opengeosys.org/project/license
 *
 */

#include "MeshGenerator.h"

#include <memory>
#include <numeric>

#include "MeshLib/Elements/Hex.h"
#include "MeshLib/Elements/Line.h"
#include "MeshLib/Elements/Pyramid.h"
#include "MeshLib/Elements/Quad.h"
#include "MeshLib/Elements/Tet.h"
#include "MeshLib/Elements/Tri.h"
#include "MeshLib/MeshEditing/AddLayerToMesh.h"
#include "MeshLib/MeshEditing/RemoveMeshComponents.h"
#include "MeshLib/Node.h"

namespace MeshLib
{

std::vector<MeshLib::Node*> MeshGenerator::generateRegularNodes(
    const std::vector<const std::vector<double>*> &vec_xyz_coords,
    const MathLib::Point3d& origin)
{
    auto const shift_coordinates = [](auto const& in, auto& out,
                                      auto const& shift) {
        std::transform(in.begin(), in.end(), std::back_inserter(out),
                       [&shift](auto const& v) { return v + shift; });
    };
    std::array<std::vector<double>, 3> coords;
    for (std::size_t i = 0; i < 3; ++i)
    {
        coords[i].reserve(vec_xyz_coords[i]->size());
        shift_coordinates(*vec_xyz_coords[i], coords[i], origin[i]);
    }

    std::vector<Node*> nodes;
    nodes.reserve(coords[0].size() * coords[1].size() * coords[2].size());

    for (auto const z : coords[2])
    {
        for (auto const y : coords[1])
        {
            std::transform(
                coords[0].begin(), coords[0].end(), std::back_inserter(nodes),
                [&y, &z](double const& x) { return new Node(x, y, z); });
        }
    }
    return nodes;
}

std::vector<MeshLib::Node*> MeshGenerator::generateRegularNodes(
    const std::vector<double> &vec_x_coords,
    const MathLib::Point3d& origin)
{
    std::vector<const std::vector<double>*> vec_xyz_coords;
    vec_xyz_coords.push_back(&vec_x_coords);
    std::vector<double> dummy(1,0.0);
    for (unsigned i = vec_xyz_coords.size() - 1; i < 3u; i++)
    {
        vec_xyz_coords.push_back(&dummy);
    }
    return generateRegularNodes(vec_xyz_coords, origin);
}

std::vector<MeshLib::Node*> MeshGenerator::generateRegularNodes(
    std::vector<double> &vec_x_coords,
    std::vector<double> &vec_y_coords,
    const MathLib::Point3d& origin)
{
    std::vector<const std::vector<double>*> vec_xyz_coords;
    vec_xyz_coords.push_back(&vec_x_coords);
    vec_xyz_coords.push_back(&vec_y_coords);
    std::vector<double> dummy(1,0.0);
    for (unsigned i = vec_xyz_coords.size() - 1; i < 3u; i++)
    {
        vec_xyz_coords.push_back(&dummy);
    }
    return generateRegularNodes(vec_xyz_coords, origin);
}

std::vector<MeshLib::Node*> MeshGenerator::generateRegularNodes(
    std::vector<double> &vec_x_coords,
    std::vector<double> &vec_y_coords,
    std::vector<double> &vec_z_coords,
    const MathLib::Point3d& origin)
{
    std::vector<const std::vector<double>*> vec_xyz_coords;
    vec_xyz_coords.push_back(&vec_x_coords);
    vec_xyz_coords.push_back(&vec_y_coords);
    vec_xyz_coords.push_back(&vec_z_coords);
    return generateRegularNodes(vec_xyz_coords, origin);
}

std::vector<MeshLib::Node*> MeshGenerator::generateRegularNodes(
    const std::array<unsigned,3> &n_cells,
    const std::array<double,3> &cell_size,
    const MathLib::Point3d& origin)
{
    std::vector<Node*> nodes;
    nodes.reserve((n_cells[0]+1)*(n_cells[1]+1)*(n_cells[2]+1));

    for (std::size_t i = 0; i < n_cells[2]+1; i++)
    {
        const double z (origin[2] + cell_size[2] * i);
        for (std::size_t j = 0; j < n_cells[1]+1; j++)
        {
            const double y (origin[1] + cell_size[1] * j);
            for (std::size_t k = 0; k < n_cells[0]+1; k++)
            {
                nodes.push_back (new Node(origin[0] + cell_size[0] * k, y, z));
            }
        }
    }
    return nodes;
}

namespace MeshGenerator
{
std::vector<MeshLib::Node*> generateRegularPyramidTopNodes(
    std::vector<double> const& x_coords,
    std::vector<double> const& y_coords,
    std::vector<double> const& z_coords,
    const MathLib::Point3d& origin)
{
    std::vector<Node*> nodes;
    nodes.reserve((x_coords.size() - 1) * (y_coords.size() - 1) *
                  (z_coords.size() - 1));

    auto const n_x_coords = x_coords.size() - 1;
    auto const n_y_coords = y_coords.size() - 1;
    auto const n_z_coords = z_coords.size() - 1;

    for (std::size_t i = 0; i < n_z_coords; i++)
    {
        const double z((z_coords[i] + z_coords[i + 1]) / 2 + origin[2]);
        for (std::size_t j = 0; j < n_y_coords; j++)
        {
            const double y((y_coords[j] + y_coords[j + 1]) / 2 + origin[1]);
            for (std::size_t k = 0; k < n_x_coords; k++)
            {
                const double x((x_coords[k] + x_coords[k + 1]) / 2 + origin[0]);
                nodes.push_back(new Node(x, y, z));
            }
        }
    }
    return nodes;
}
}  // end namespace MeshGenerator

Mesh* MeshGenerator::generateLineMesh(
    const double length,
    const std::size_t subdivision,
    const MathLib::Point3d& origin,
    std::string   const& mesh_name)
{
    return generateLineMesh(subdivision, length/subdivision, origin, mesh_name);
}

Mesh* MeshGenerator::generateLineMesh(
    const unsigned n_cells,
    const double   cell_size,
    MathLib::Point3d const& origin,
    std::string   const& mesh_name)
{
    return generateLineMesh(BaseLib::UniformSubdivision(n_cells*cell_size, n_cells), origin, mesh_name);
}

Mesh* MeshGenerator::generateLineMesh(
    const BaseLib::ISubdivision &div,
    MathLib::Point3d const& origin,
    std::string const& mesh_name)
{
    const std::vector<double> vec_x(div());
    std::vector<Node*> nodes(generateRegularNodes(vec_x, origin));

    //elements
    const std::size_t n_cells = nodes.size()-1;
    std::vector<Element*> elements;
    elements.reserve(n_cells);

    for (std::size_t i = 0; i < n_cells; i++)
    {
        elements.push_back(new Line({nodes[i], nodes[i + 1]}));
    }

    return new Mesh(mesh_name, nodes, elements);
}

Mesh* MeshGenerator::generateRegularQuadMesh(
    const double length,
    const std::size_t subdivision,
    const MathLib::Point3d& origin,
    std::string const& mesh_name)
{
    return generateRegularQuadMesh(subdivision, subdivision,
        length/subdivision, length/subdivision, origin, mesh_name);
}

Mesh* MeshGenerator::generateRegularQuadMesh(
    const double x_length,
    const double y_length,
    const std::size_t x_subdivision,
    const std::size_t y_subdivision,
    const MathLib::Point3d& origin,
    std::string const& mesh_name)
{
    return generateRegularQuadMesh(x_subdivision, y_subdivision,
        x_length/x_subdivision, y_length/y_subdivision, origin, mesh_name);
}

Mesh* MeshGenerator::generateRegularQuadMesh(
    const unsigned n_x_cells,
    const unsigned n_y_cells,
    const double cell_size,
    MathLib::Point3d const& origin,
    std::string const& mesh_name)
{
    return generateRegularQuadMesh(n_x_cells, n_y_cells, cell_size, cell_size, origin, mesh_name);
}

Mesh* MeshGenerator::generateRegularQuadMesh(
    const unsigned n_x_cells,
    const unsigned n_y_cells,
    const double cell_size_x,
    const double cell_size_y,
    MathLib::Point3d const& origin,
    std::string const& mesh_name)
{
    return generateRegularQuadMesh(BaseLib::UniformSubdivision(n_x_cells*cell_size_x, n_x_cells),
        BaseLib::UniformSubdivision(n_y_cells*cell_size_y, n_y_cells), origin, mesh_name);
}

Mesh* MeshGenerator::generateRegularQuadMesh(
    const BaseLib::ISubdivision &div_x,
    const BaseLib::ISubdivision &div_y,
    MathLib::Point3d const& origin,
    std::string const& mesh_name)
{
    std::vector<double> vec_x(div_x());
    std::vector<double> vec_y(div_y());
    std::vector<Node*> nodes(generateRegularNodes(vec_x, vec_y, origin));
    const unsigned n_x_nodes (vec_x.size());

    //elements
    std::vector<Element*> elements;
    const unsigned n_x_cells (vec_x.size()-1);
    const unsigned n_y_cells (vec_y.size()-1);
    elements.reserve(n_x_cells * n_y_cells);

    for (std::size_t j = 0; j < n_y_cells; j++)
    {
        const std::size_t offset_y1 = j * n_x_nodes;
        const std::size_t offset_y2 = (j + 1) * n_x_nodes;
        for (std::size_t k = 0; k < n_x_cells; k++)
        {
            elements.push_back(
                new Quad({nodes[offset_y1 + k], nodes[offset_y1 + k + 1],
                          nodes[offset_y2 + k + 1], nodes[offset_y2 + k]}));
        }
    }

    return new Mesh(mesh_name, nodes, elements);
}

Mesh* MeshGenerator::generateRegularHexMesh(
    const double length,
    const std::size_t subdivision,
    const MathLib::Point3d& origin,
    std::string const& mesh_name)
{
    return MeshGenerator::generateRegularHexMesh(subdivision, subdivision,
        subdivision, length/subdivision, origin, mesh_name);
}

Mesh* MeshGenerator::generateRegularHexMesh(
    const double x_length,
    const double y_length,
    const double z_length,
    const std::size_t x_subdivision,
    const std::size_t y_subdivision,
    const std::size_t z_subdivision,
    const MathLib::Point3d& origin,
    std::string const& mesh_name)
{
    return MeshGenerator::generateRegularHexMesh(x_subdivision, y_subdivision, z_subdivision,
        x_length/x_subdivision, y_length/y_subdivision, z_length/z_subdivision, origin, mesh_name);
}

Mesh* MeshGenerator::generateRegularHexMesh(
    const unsigned n_x_cells,
    const unsigned n_y_cells,
    const unsigned n_z_cells,
    const double   cell_size,
    MathLib::Point3d const& origin,
    std::string   const& mesh_name)
{
    return MeshGenerator::generateRegularHexMesh(n_x_cells, n_y_cells, n_z_cells,
        cell_size, cell_size, cell_size, origin, mesh_name);
}

Mesh* MeshGenerator::generateRegularHexMesh(
    const unsigned n_x_cells,
    const unsigned n_y_cells,
    const unsigned n_z_cells,
    const double   cell_size_x,
    const double   cell_size_y,
    const double   cell_size_z,
    MathLib::Point3d const& origin,
    std::string   const& mesh_name)
{
    return generateRegularHexMesh(
            BaseLib::UniformSubdivision(n_x_cells*cell_size_x, n_x_cells),
            BaseLib::UniformSubdivision(n_y_cells*cell_size_y, n_y_cells),
            BaseLib::UniformSubdivision(n_z_cells*cell_size_z, n_z_cells),
            origin, mesh_name);
}

Mesh* MeshGenerator::generateRegularHexMesh(
    const BaseLib::ISubdivision &div_x,
    const BaseLib::ISubdivision &div_y,
    const BaseLib::ISubdivision &div_z,
    MathLib::Point3d const& origin,
    std::string const& mesh_name)
{
    std::vector<double> vec_x(div_x());
    std::vector<double> vec_y(div_y());
    std::vector<double> vec_z(div_z());
    std::vector<Node*> nodes(generateRegularNodes(vec_x, vec_y, vec_z, origin));

    const unsigned n_x_nodes (vec_x.size());
    const unsigned n_y_nodes (vec_y.size());
    const unsigned n_x_cells (vec_x.size()-1);
    const unsigned n_y_cells (vec_y.size()-1);
    const unsigned n_z_cells (vec_z.size()-1);

    //elements
    std::vector<Element*> elements;
    elements.reserve(n_x_cells * n_y_cells * n_z_cells);

    for (std::size_t i = 0; i < n_z_cells; i++)
    {
        const std::size_t offset_z1 = i * n_x_nodes * n_y_nodes; // bottom
        const std::size_t offset_z2 = (i + 1) * n_x_nodes * n_y_nodes; // top
        for (std::size_t j = 0; j < n_y_cells; j++)
        {
            const std::size_t offset_y1 = j * n_x_nodes;
            const std::size_t offset_y2 = (j + 1) * n_x_nodes;
            for (std::size_t k = 0; k < n_x_cells; k++)
            {
                elements.push_back(
                    new Hex({// bottom
                             nodes[offset_z1 + offset_y1 + k],
                             nodes[offset_z1 + offset_y1 + k + 1],
                             nodes[offset_z1 + offset_y2 + k + 1],
                             nodes[offset_z1 + offset_y2 + k],
                             // top
                             nodes[offset_z2 + offset_y1 + k],
                             nodes[offset_z2 + offset_y1 + k + 1],
                             nodes[offset_z2 + offset_y2 + k + 1],
                             nodes[offset_z2 + offset_y2 + k]}));
            }
        }
    }

    return new Mesh(mesh_name, nodes, elements);
}

Mesh* MeshGenerator::generateRegularPyramidMesh(
    const BaseLib::ISubdivision& div_x,
    const BaseLib::ISubdivision& div_y,
    const BaseLib::ISubdivision& div_z,
    MathLib::Point3d const& origin,
    std::string const& mesh_name)
{
    std::vector<double> vec_x(div_x());
    std::vector<double> vec_y(div_y());
    std::vector<double> vec_z(div_z());
    std::vector<Node*> nodes(generateRegularNodes(vec_x, vec_y, vec_z, origin));
    std::vector<Node*> const top_nodes(
        generateRegularPyramidTopNodes(vec_x, vec_y, vec_z, origin));

    nodes.insert(nodes.end(), top_nodes.begin(), top_nodes.end());

    const unsigned n_x_nodes(vec_x.size());
    const unsigned n_y_nodes(vec_y.size());
    const unsigned n_z_nodes(vec_z.size());
    const unsigned n_x_cells(vec_x.size() - 1);
    const unsigned n_y_cells(vec_y.size() - 1);
    const unsigned n_z_cells(vec_z.size() - 1);

    // elements
    std::vector<Element*> elements;
    auto const top_node_offset(n_x_nodes * n_y_nodes * n_z_nodes);
    elements.reserve(n_x_cells * n_y_cells * n_z_cells);

    for (std::size_t i = 0; i < n_z_cells; i++)
    {
        const std::size_t offset_z1 = i * n_x_nodes * n_y_nodes;  // bottom
        const std::size_t offset_z2 = (i + 1) * n_x_nodes * n_y_nodes;  // top
        for (std::size_t j = 0; j < n_y_cells; j++)
        {
            const std::size_t offset_y1 = j * n_x_nodes;
            const std::size_t offset_y2 = (j + 1) * n_x_nodes;
            for (std::size_t k = 0; k < n_x_cells; k++)
            {
                // generate 6 pyramids within the virtual hexahedron cell
                int const pyramid_top_index(i * n_x_cells * n_y_cells +
                                            j * n_x_cells + k +
                                            top_node_offset);
                elements.push_back(
                    new Pyramid{{// bottom 'hexahedron' face
                                 nodes[offset_z1 + offset_y1 + k],
                                 nodes[offset_z1 + offset_y1 + k + 1],
                                 nodes[offset_z1 + offset_y2 + k + 1],
                                 nodes[offset_z1 + offset_y2 + k],
                                 // top
                                 nodes[pyramid_top_index]}});
                elements.push_back(
                    new Pyramid{{// top 'hexahedron' face
                                 nodes[offset_z2 + offset_y1 + k + 1],
                                 nodes[offset_z2 + offset_y1 + k],
                                 nodes[offset_z2 + offset_y2 + k],
                                 nodes[offset_z2 + offset_y2 + k + 1],
                                 // top of pyramid directed towards the bottom
                                 nodes[pyramid_top_index]}});
                elements.push_back(
                    new Pyramid{{// right 'hexahedron' face
                                 nodes[offset_z1 + offset_y1 + k + 1],
                                 nodes[offset_z2 + offset_y1 + k + 1],
                                 nodes[offset_z2 + offset_y2 + k + 1],
                                 nodes[offset_z1 + offset_y2 + k + 1],
                                 // top of pyramid directed towards the bottom
                                 nodes[pyramid_top_index]}});
                elements.push_back(
                    new Pyramid{{// left 'hexahedron' face
                                 nodes[offset_z2 + offset_y1 + k],
                                 nodes[offset_z1 + offset_y1 + k],
                                 nodes[offset_z1 + offset_y2 + k],
                                 nodes[offset_z2 + offset_y2 + k],
                                 // top of pyramid directed towards the bottom
                                 nodes[pyramid_top_index]}});
                elements.push_back(
                    new Pyramid{{// front 'hexahedron' face
                                 nodes[offset_z2 + offset_y1 + k],
                                 nodes[offset_z2 + offset_y1 + k + 1],
                                 nodes[offset_z1 + offset_y1 + k + 1],
                                 nodes[offset_z1 + offset_y1 + k],
                                 // top of pyramid directed towards the bottom
                                 nodes[pyramid_top_index]}});
                elements.push_back(
                    new Pyramid{{// back 'hexahedron' face
                                 nodes[offset_z1 + offset_y2 + k],
                                 nodes[offset_z1 + offset_y2 + k + 1],
                                 nodes[offset_z2 + offset_y2 + k + 1],
                                 nodes[offset_z2 + offset_y2 + k],
                                 // top of pyramid directed towards the bottom
                                 nodes[pyramid_top_index]}});
            }
        }
    }
    return new Mesh(mesh_name, nodes, elements);
}

Mesh* MeshGenerator::generateRegularTriMesh(
    const double length,
    const std::size_t subdivision,
    const MathLib::Point3d& origin,
    std::string const& mesh_name)
{
    return generateRegularTriMesh(subdivision, subdivision, length/subdivision, origin, mesh_name);
}

Mesh* MeshGenerator::generateRegularTriMesh(
    const double x_length,
    const double y_length,
    const std::size_t x_subdivision,
    const std::size_t y_subdivision,
    const MathLib::Point3d& origin,
    std::string const& mesh_name)
{
    return generateRegularTriMesh(x_subdivision, y_subdivision, x_length/x_subdivision, y_length/y_subdivision, origin, mesh_name);
}

Mesh* MeshGenerator::generateRegularTriMesh(
    const unsigned n_x_cells,
    const unsigned n_y_cells,
    const double cell_size,
    MathLib::Point3d const& origin,
    std::string const& mesh_name)
{
    return generateRegularTriMesh(n_x_cells, n_y_cells, cell_size, cell_size, origin, mesh_name);
}

Mesh* MeshGenerator::generateRegularTriMesh(
    const unsigned n_x_cells,
    const unsigned n_y_cells,
    const double   cell_size_x,
    const double   cell_size_y,
    MathLib::Point3d const& origin,
    std::string   const& mesh_name)
{
    return generateRegularTriMesh(BaseLib::UniformSubdivision(n_x_cells*cell_size_x, n_x_cells),
        BaseLib::UniformSubdivision(n_y_cells*cell_size_y, n_y_cells), origin, mesh_name);
}

Mesh* MeshGenerator::generateRegularTriMesh(
    const BaseLib::ISubdivision &div_x,
    const BaseLib::ISubdivision &div_y,
    MathLib::Point3d const& origin,
    std::string const& mesh_name)
{
    std::vector<double> vec_x(div_x());
    std::vector<double> vec_y(div_y());
    std::vector<Node*> nodes(generateRegularNodes(vec_x, vec_y, origin));
    const unsigned n_x_nodes (vec_x.size());
    const unsigned n_x_cells (vec_x.size()-1);
    const unsigned n_y_cells (vec_y.size()-1);

    //elements
    std::vector<Element*> elements;
    elements.reserve(n_x_cells * n_y_cells * 2);

    for (std::size_t j = 0; j < n_y_cells; j++)
    {
        const std::size_t offset_y1 = j * n_x_nodes;
        const std::size_t offset_y2 = (j + 1) * n_x_nodes;
        for (std::size_t k = 0; k < n_x_cells; k++)
        {
            elements.push_back(
                new Tri({nodes[offset_y1 + k], nodes[offset_y2 + k + 1],
                         nodes[offset_y2 + k]}));

            elements.push_back(
                new Tri({nodes[offset_y1 + k], nodes[offset_y1 + k + 1],
                         nodes[offset_y2 + k + 1]}));
        }
    }

    return new Mesh(mesh_name, nodes, elements);
}

Mesh* MeshGenerator::generateRegularPrismMesh(
    const double x_length,
    const double y_length,
    const double z_length,
    const std::size_t x_subdivision,
    const std::size_t y_subdivision,
    const std::size_t z_subdivision,
    const MathLib::Point3d& origin,
    std::string const& mesh_name)
{
    return generateRegularPrismMesh(x_subdivision, y_subdivision, z_subdivision,
        x_length/x_subdivision, y_length/y_subdivision, z_length/z_subdivision,
        origin, mesh_name);
}

Mesh* MeshGenerator::generateRegularPrismMesh(
    const unsigned n_x_cells,
    const unsigned n_y_cells,
    const unsigned n_z_cells,
    const double cell_size,
    MathLib::Point3d const& origin,
    std::string const& mesh_name)
{
    return generateRegularPrismMesh(n_x_cells, n_y_cells, n_z_cells,
        cell_size, cell_size, cell_size, origin, mesh_name);
}

Mesh* MeshGenerator::generateRegularPrismMesh(
    const unsigned n_x_cells,
    const unsigned n_y_cells,
    const unsigned n_z_cells,
    const double   cell_size_x,
    const double   cell_size_y,
    const double   cell_size_z,
    MathLib::Point3d const& origin,
    std::string   const& mesh_name)
{
    std::unique_ptr<MeshLib::Mesh> mesh (
        generateRegularTriMesh(n_x_cells, n_y_cells, cell_size_x, cell_size_y, origin, mesh_name));
    std::size_t const n_tris (mesh->getNumberOfElements());
    bool const copy_material_ids = false;
    for (std::size_t i = 0; i < n_z_cells; ++i)
    {
        mesh.reset(MeshLib::addTopLayerToMesh(*mesh, cell_size_z, mesh_name,
                                              copy_material_ids));
    }
    std::vector<std::size_t> elem_ids (n_tris);
    std::iota(elem_ids.begin(), elem_ids.end(), 0);
    return MeshLib::removeElements(*mesh, elem_ids, mesh_name);
}

Mesh* MeshGenerator::generateRegularTetMesh(
    const double x_length,
    const double y_length,
    const double z_length,
    const std::size_t x_subdivision,
    const std::size_t y_subdivision,
    const std::size_t z_subdivision,
    const MathLib::Point3d& origin,
    std::string const& mesh_name)
{
    return generateRegularTetMesh(x_subdivision, y_subdivision, z_subdivision,
        x_length/x_subdivision, y_length/y_subdivision, z_length/z_subdivision,
        origin, mesh_name);
}

Mesh* MeshGenerator::generateRegularTetMesh(
    const unsigned n_x_cells,
    const unsigned n_y_cells,
    const unsigned n_z_cells,
    const double   cell_size_x,
    const double   cell_size_y,
    const double   cell_size_z,
    MathLib::Point3d const& origin,
    std::string   const& mesh_name)
{
    return generateRegularTetMesh(
            BaseLib::UniformSubdivision(n_x_cells*cell_size_x, n_x_cells),
            BaseLib::UniformSubdivision(n_y_cells*cell_size_y, n_y_cells),
            BaseLib::UniformSubdivision(n_z_cells*cell_size_z, n_z_cells),
            origin, mesh_name);
}

Mesh* MeshGenerator::generateRegularTetMesh(
    const BaseLib::ISubdivision &div_x,
    const BaseLib::ISubdivision &div_y,
    const BaseLib::ISubdivision &div_z,
    MathLib::Point3d const& origin,
    std::string const& mesh_name)
{
    std::vector<double> vec_x(div_x());
    std::vector<double> vec_y(div_y());
    std::vector<double> vec_z(div_z());
    std::vector<Node*> nodes(generateRegularNodes(vec_x, vec_y, vec_z, origin));

    const unsigned n_x_nodes (vec_x.size());
    const unsigned n_y_nodes (vec_y.size());
    const unsigned n_x_cells (vec_x.size()-1);
    const unsigned n_y_cells (vec_y.size()-1);
    const unsigned n_z_cells (vec_z.size()-1);

    //elements
    std::vector<Element*> elements;
    elements.reserve(n_x_cells * n_y_cells * n_z_cells * 6);

    for (std::size_t i = 0; i < n_z_cells; i++)
    {
        const std::size_t offset_z1 = i * n_x_nodes * n_y_nodes; // bottom
        const std::size_t offset_z2 = (i + 1) * n_x_nodes * n_y_nodes; // top
        for (std::size_t j = 0; j < n_y_cells; j++)
        {
            const std::size_t offset_y1 = j * n_x_nodes;
            const std::size_t offset_y2 = (j + 1) * n_x_nodes;
            for (std::size_t k = 0; k < n_x_cells; k++)
            {
                // tet 1
                elements.push_back(
                    new Tet({// bottom
                             nodes[offset_z1 + offset_y1 + k],
                             nodes[offset_z1 + offset_y2 + k + 1],
                             nodes[offset_z1 + offset_y2 + k],
                             // top
                             nodes[offset_z2 + offset_y1 + k]}));
                // tet 2
                elements.push_back(
                    new Tet({// bottom
                             nodes[offset_z1 + offset_y2 + k + 1],
                             nodes[offset_z1 + offset_y2 + k],
                             // top
                             nodes[offset_z2 + offset_y1 + k],
                             nodes[offset_z2 + offset_y2 + k + 1]}));
                // tet 3
                elements.push_back(
                    new Tet({// bottom
                             nodes[offset_z1 + offset_y2 + k],
                             // top
                             nodes[offset_z2 + offset_y1 + k],
                             nodes[offset_z2 + offset_y2 + k + 1],
                             nodes[offset_z2 + offset_y2 + k]}));
                // tet 4
                elements.push_back(
                    new Tet({// bottom
                             nodes[offset_z1 + offset_y1 + k],
                             nodes[offset_z1 + offset_y1 + k + 1],
                             nodes[offset_z1 + offset_y2 + k + 1],
                             // top
                             nodes[offset_z2 + offset_y1 + k + 1]}));
                // tet 5
                elements.push_back(
                    new Tet({// bottom
                             nodes[offset_z1 + offset_y1 + k],
                             nodes[offset_z1 + offset_y2 + k + 1],
                             // top
                             nodes[offset_z2 + offset_y1 + k],
                             nodes[offset_z2 + offset_y1 + k + 1]}));
                // tet 6
                elements.push_back(
                    new Tet({// bottom
                             nodes[offset_z1 + offset_y2 + k + 1],
                             // top
                             nodes[offset_z2 + offset_y1 + k],
                             nodes[offset_z2 + offset_y1 + k + 1],
                             nodes[offset_z2 + offset_y2 + k + 1]}));
            }
        }
    }

    return new Mesh(mesh_name, nodes, elements);
}

MeshLib::Mesh*
MeshGenerator::createSurfaceMesh(std::string const& mesh_name,
    MathLib::Point3d const& ll, MathLib::Point3d const& ur,
    std::array<std::size_t, 2> const& n_steps,
    const std::function<double(double,double)>& f)
{
    std::array<double, 2> step_size{{
        (ur[0]-ll[0])/(n_steps[0]-1), (ur[1]-ll[1])/(n_steps[1]-1)}};

    std::vector<MeshLib::Node*> nodes;
    for (std::size_t j(0); j<n_steps[1]; ++j) {
        for (std::size_t i(0); i<n_steps[0]; ++i) {
            std::size_t const id = i + j * n_steps[1];
            double const x = ll[0] + i * step_size[0];
            double const y = ll[1] + j * step_size[1];

            nodes.push_back(new MeshLib::Node({x, y, f(x, y)}, id));
        }
    }

    std::vector<MeshLib::Element*> sfc_eles;
    for (std::size_t j(0); j<n_steps[1]-1; ++j) {
        for (std::size_t i(0); i<n_steps[0]-1; ++i) {
            std::size_t id_ll(i+j*n_steps[0]);
            std::size_t id_lr(i+1+j*n_steps[0]);
            std::size_t id_ul(i+(j+1)*n_steps[0]);
            std::size_t id_ur(i+1+(j+1)*n_steps[0]);
            sfc_eles.push_back(
                new MeshLib::Tri({nodes[id_ll], nodes[id_lr], nodes[id_ur]}));
            sfc_eles.push_back(
                new MeshLib::Tri({nodes[id_ll], nodes[id_ur], nodes[id_ul]}));
        }
    }

    return new MeshLib::Mesh(mesh_name, nodes, sfc_eles);
}

}  // namespace MeshLib