#ifndef OERSTED_FINITEELEMENTMESH_H
#define OERSTED_FINITEELEMENTMESH_H

#include <memory>
#include <vector>

#include "Mesh.hpp"

#include "Boundary.h"
#include "MatrixGroup.h"
#include "Node.h"
#include "Triangle.h"
#include "Region.h"

template<size_t Dimension, size_t Order>
class FiniteElementMesh {};

template<size_t Order>
class FiniteElementMesh<2,Order> {
public:
    FiniteElementMesh() {};

    FiniteElementMesh(std::vector<XY> nodes, std::vector<Triangle<Order>> tris, std::vector<std::shared_ptr<Region<2>>> r, std::vector<std::shared_ptr<Boundary<2>>> b) : Nodes(nodes), Triangles(tris), Regions(r), Boundaries(b) {};

    FiniteElementMesh(Mesh &m);

    size_t size_nodes() const { return Nodes.size(); };

    size_t size_elements() const { return Triangles.size(); };

    std::vector<Triangle<Order>> const &triangles() const { return Triangles; };

    Triangle<Order> const &triangle(size_t i) const { return Triangles[i]; };

    auto &boundaries() { return Boundaries; };

    auto &boundary(size_t i) { return Boundaries[i]; };

    auto const &boundaries() const { return Boundaries; };

    auto const &nodes() const { return Nodes; };

    auto const &node(size_t i) const { return Nodes[i]; };

    auto &regions() { return Regions; };

    auto &region(size_t i) { return Regions[i]; };

    auto const &regions() const { return Regions; };

    auto const &region(size_t i) const { return Regions[i]; };

    template<size_t Q>
    DiagonalMatrixGroup<TriangleQuadrature<Q>::size> determinant() const {
        DiagonalMatrixGroup<TriangleQuadrature<Q>::size> mat(Triangles.size());

        for (size_t i = 0; i != Triangles.size(); ++i) {
            Triangles[i].determinant<Q>(mat, Nodes);
        }

        return mat;
    }

    template<size_t Q>
    SparseMatrixGroup<TriangleQuadrature<Q>::size> basis() const {
        SparseMatrixGroup<TriangleQuadrature<Q>::size> mat(Nodes.size(), Triangles.size(), Triangle<Order>::NumNodes);

        for (size_t i = 0; i != Triangles.size(); ++i) {
            Triangles[i].basis<Q>(mat, Nodes);
        }

        return mat;
    }

    template<size_t Q>
    DerivativeMatrixGroup<TriangleQuadrature<Q>::size> derivative() const {
        DerivativeMatrixGroup<TriangleQuadrature<Q>::size> df(Nodes.size(), Triangles.size(), Triangle<Order>::NumNodes);

        for (size_t i = 0; i!= Triangles.size();++i) {
            Triangles[i].derivative<Q>(df, Nodes);
        }

        return df;
    }

    template<size_t Q>
    std::vector<std::vector<XY>> quadrature_points() const {
        std::vector<std::vector<XY>> qp;
        qp.reserve(Triangles.size());

        for(size_t i = 0; i!= Triangles.size();++i){
            qp.push_back(Triangles[i].quadrature_points<Q>(Nodes));
        }

        return qp;
    }

protected:
    std::vector<XY> Nodes;
    std::vector<Triangle<Order>> Triangles;

    std::vector<std::shared_ptr<Region<2>>> Regions; // Contains vector of size_t referencing Triangles (and later Quadrilaterals)
    std::vector<std::shared_ptr<Boundary<2>>> Boundaries;
};

#endif //OERSTED_FINITEELEMENTMESH_H