/**
 * \copyright
 * Copyright (c) 2012-2021, 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 "NumLib/Fem/CoordinatesMapping/NaturalNodeCoordinates.h"

#include <gtest/gtest.h>

#include <vector>

#include "MeshLib/ElementCoordinatesMappingLocal.h"
#include "NumLib/Fem/CoordinatesMapping/ShapeMatrices.h"
#include "NumLib/Fem/FiniteElement/C0IsoparametricElements.h"
#include "NumLib/Fem/InitShapeMatrices.h"
#include "NumLib/Fem/Integration/GaussLegendreIntegrationPolicy.h"
#include "NumLib/Fem/ShapeMatrixPolicy.h"

using Node = MeshLib::Node;

template <typename ShapeFunction, int GlobalDim>
bool test(MeshLib::Element const& element)
{
    using ShapeMatricesType = ShapeMatrixPolicyType<ShapeFunction, GlobalDim>;

    // Test if the cast is possible.
    bool const dynamic_cast_failed =
        dynamic_cast<typename ShapeFunction::MeshElement const*>(&element) ==
        nullptr;
    EXPECT_FALSE(dynamic_cast_failed);
    if (dynamic_cast_failed)
    {
        return false;
    }

    // For each node evaluate the shape functions at natural coordinates and
    // test if only the corresponding shape function has value 1 and all others
    // must return 0.
    int const number_nodes = element.getNumberOfNodes();
    std::vector<MathLib::Point3d> points;
    points.reserve(number_nodes);
    for (int n = 0; n < number_nodes; ++n)
    {
        // Evaluate shape matrices at natural coordinates.
        // Compute only N, because for pyramid the detJ becomes zero at the
        // apex, and we only use N in the following test anyway.
        points.emplace_back(
            NumLib::NaturalCoordinates<
                typename ShapeFunction::MeshElement>::coordinates[n]);
    }
    auto const shape_matrices =
        NumLib::computeShapeMatrices<ShapeFunction, ShapeMatricesType,
                                     GlobalDim, NumLib::ShapeMatrixType::N>(
            element, false /*is_axially_symmetric*/, points);

    bool result = true;
    for (int n = 0; n < number_nodes; ++n)
    {
        auto const& N = shape_matrices[n].N;
        for (int p = 0; p < static_cast<int>(ShapeFunction::NPOINTS); ++p)
        {
            if (p == n)
            {
                if (N[p] != 1)
                {
                    EXPECT_EQ(1, N[p]) << "for n = " << n << ", p = " << p
                                       << " and dimension " << GlobalDim;
                    result = false;
                }
            }
            else
            {
                if (N[p] != 0)
                {
                    EXPECT_EQ(0, N[p]) << "for n = " << n << ", p = " << p
                                       << " and dimension " << GlobalDim;
                    result = false;
                }
            }
        }
    }
    return result;
}

TEST(NumLibNaturalCoordinates, Line2)
{
    using ShapeFunction = NumLib::ShapeLine2;
    using MeshElement = ShapeFunction::MeshElement;

    auto const& xs = NumLib::NaturalCoordinates<MeshElement>::coordinates;
    std::array<Node, MeshElement::n_all_nodes> ns{{Node{xs[0]}, Node{xs[1]}}};
    std::array<Node*, MeshElement::n_all_nodes> nodes{{&ns[0], &ns[1]}};
    MeshElement element(nodes);

    bool const test_result1 = test<ShapeFunction, 1>(element);
    ASSERT_TRUE(test_result1) << "Test failed for dimension 1.";

    bool const test_result2 = test<ShapeFunction, 2>(element);
    ASSERT_TRUE(test_result2) << "Test failed for dimension 2.";

    bool const test_result3 = test<ShapeFunction, 3>(element);
    ASSERT_TRUE(test_result3) << "Test failed for dimension 3.";
}

TEST(NumLibNaturalCoordinates, Line3)
{
    using ShapeFunction = NumLib::ShapeLine3;
    using MeshElement = ShapeFunction::MeshElement;

    auto const& xs = NumLib::NaturalCoordinates<MeshElement>::coordinates;
    std::array<Node, MeshElement::n_all_nodes> ns{
        {Node{xs[0]}, Node{xs[1]}, Node{xs[2]}}};
    std::array<Node*, MeshElement::n_all_nodes> nodes{{&ns[0], &ns[1], &ns[2]}};
    MeshElement element(nodes);

    bool const test_result1 = test<ShapeFunction, 1>(element);
    ASSERT_TRUE(test_result1) << "Test failed for dimension 1.";

    bool const test_result2 = test<ShapeFunction, 2>(element);
    ASSERT_TRUE(test_result2) << "Test failed for dimension 2.";

    bool const test_result3 = test<ShapeFunction, 3>(element);
    ASSERT_TRUE(test_result3) << "Test failed for dimension 3.";
}

TEST(NumLibNaturalCoordinates, Tri3)
{
    using ShapeFunction = NumLib::ShapeTri3;
    using MeshElement = ShapeFunction::MeshElement;

    auto const& xs = NumLib::NaturalCoordinates<MeshElement>::coordinates;
    std::array<Node, MeshElement::n_all_nodes> ns{
        {Node{xs[0]}, Node{xs[1]}, Node{xs[2]}}};
    std::array<Node*, MeshElement::n_all_nodes> nodes{{&ns[0], &ns[1], &ns[2]}};
    MeshElement element(nodes);

    bool const test_result2 = test<ShapeFunction, 2>(element);
    ASSERT_TRUE(test_result2) << "Test failed for dimension 2.";

    bool const test_result3 = test<ShapeFunction, 3>(element);
    ASSERT_TRUE(test_result3) << "Test failed for dimension 3.";
}

TEST(NumLibNaturalCoordinates, Tri6)
{
    using ShapeFunction = NumLib::ShapeTri6;
    using MeshElement = ShapeFunction::MeshElement;

    auto const& xs = NumLib::NaturalCoordinates<MeshElement>::coordinates;
    std::array<Node, MeshElement::n_all_nodes> ns{{Node{xs[0]}, Node{xs[1]},
                                                   Node{xs[2]}, Node{xs[3]},
                                                   Node{xs[4]}, Node{xs[5]}}};
    std::array<Node*, MeshElement::n_all_nodes> nodes{
        {&ns[0], &ns[1], &ns[2], &ns[3], &ns[4], &ns[5]}};
    MeshElement element(nodes);

    bool const test_result2 = test<ShapeFunction, 2>(element);
    ASSERT_TRUE(test_result2) << "Test failed for dimension 2.";

    bool const test_result3 = test<ShapeFunction, 3>(element);
    ASSERT_TRUE(test_result3) << "Test failed for dimension 3.";
}

TEST(NumLibNaturalCoordinates, Quad4)
{
    using ShapeFunction = NumLib::ShapeQuad4;
    using MeshElement = ShapeFunction::MeshElement;

    auto const& xs = NumLib::NaturalCoordinates<MeshElement>::coordinates;
    std::array<Node, MeshElement::n_all_nodes> ns{
        {Node{xs[0]}, Node{xs[1]}, Node{xs[2]}, Node{xs[3]}}};
    std::array<Node*, MeshElement::n_all_nodes> nodes{
        {&ns[0], &ns[1], &ns[2], &ns[3]}};
    MeshElement element(nodes);

    bool const test_result2 = test<ShapeFunction, 2>(element);
    ASSERT_TRUE(test_result2) << "Test failed for dimension 2.";

    bool const test_result3 = test<ShapeFunction, 3>(element);
    ASSERT_TRUE(test_result3) << "Test failed for dimension 3.";
}

TEST(NumLibNaturalCoordinates, Quad8)
{
    using ShapeFunction = NumLib::ShapeQuad8;
    using MeshElement = ShapeFunction::MeshElement;

    auto const& xs = NumLib::NaturalCoordinates<MeshElement>::coordinates;
    std::array<Node, MeshElement::n_all_nodes> ns{
        {Node{xs[0]}, Node{xs[1]}, Node{xs[2]}, Node{xs[3]}, Node{xs[4]},
         Node{xs[5]}, Node{xs[6]}, Node{xs[7]}}};
    std::array<Node*, MeshElement::n_all_nodes> nodes{
        {&ns[0], &ns[1], &ns[2], &ns[3], &ns[4], &ns[5], &ns[6], &ns[7]}};
    MeshElement element(nodes);

    bool const test_result2 = test<ShapeFunction, 2>(element);
    ASSERT_TRUE(test_result2) << "Test failed for dimension 2.";

    bool const test_result3 = test<ShapeFunction, 3>(element);
    ASSERT_TRUE(test_result3) << "Test failed for dimension 3.";
}

TEST(NumLibNaturalCoordinates, Quad9)
{
    using ShapeFunction = NumLib::ShapeQuad9;
    using MeshElement = ShapeFunction::MeshElement;

    auto const& xs = NumLib::NaturalCoordinates<MeshElement>::coordinates;
    std::array<Node, MeshElement::n_all_nodes> ns{
        {Node{xs[0]}, Node{xs[1]}, Node{xs[2]}, Node{xs[3]}, Node{xs[4]},
         Node{xs[5]}, Node{xs[6]}, Node{xs[7]}, Node{xs[8]}}};
    std::array<Node*, MeshElement::n_all_nodes> nodes{{&ns[0], &ns[1], &ns[2],
                                                       &ns[3], &ns[4], &ns[5],
                                                       &ns[6], &ns[7], &ns[8]}};
    MeshElement element(nodes);

    bool const test_result2 = test<ShapeFunction, 2>(element);
    ASSERT_TRUE(test_result2) << "Test failed for dimension 2.";

    bool const test_result3 = test<ShapeFunction, 3>(element);
    ASSERT_TRUE(test_result3) << "Test failed for dimension 3.";
}

TEST(NumLibNaturalCoordinates, Tet4)
{
    using ShapeFunction = NumLib::ShapeTet4;
    using MeshElement = ShapeFunction::MeshElement;

    auto const& xs = NumLib::NaturalCoordinates<MeshElement>::coordinates;
    std::array<Node, MeshElement::n_all_nodes> ns{
        {Node{xs[0]}, Node{xs[1]}, Node{xs[2]}, Node{xs[3]}}};
    std::array<Node*, MeshElement::n_all_nodes> nodes{
        {&ns[0], &ns[1], &ns[2], &ns[3]}};
    MeshElement element(nodes);

    bool const test_result3 = test<ShapeFunction, 3>(element);
    ASSERT_TRUE(test_result3) << "Test failed for dimension 3.";
}

TEST(NumLibNaturalCoordinates, Tet10)
{
    using ShapeFunction = NumLib::ShapeTet10;
    using MeshElement = ShapeFunction::MeshElement;

    auto const& xs = NumLib::NaturalCoordinates<MeshElement>::coordinates;
    std::array<Node, MeshElement::n_all_nodes> ns{
        {Node{xs[0]}, Node{xs[1]}, Node{xs[2]}, Node{xs[3]}, Node{xs[4]},
         Node{xs[5]}, Node{xs[6]}, Node{xs[7]}, Node{xs[8]}, Node{xs[9]}}};
    std::array<Node*, MeshElement::n_all_nodes> nodes{
        {&ns[0], &ns[1], &ns[2], &ns[3], &ns[4], &ns[5], &ns[6], &ns[7], &ns[8],
         &ns[9]}};
    MeshElement element(nodes);

    bool const test_result3 = test<ShapeFunction, 3>(element);
    ASSERT_TRUE(test_result3) << "Test failed for dimension 3.";
}

TEST(NumLibNaturalCoordinates, Prism6)
{
    using ShapeFunction = NumLib::ShapePrism6;
    using MeshElement = ShapeFunction::MeshElement;

    auto const& xs = NumLib::NaturalCoordinates<MeshElement>::coordinates;
    std::array<Node, MeshElement::n_all_nodes> ns{{Node{xs[0]}, Node{xs[1]},
                                                   Node{xs[2]}, Node{xs[3]},
                                                   Node{xs[4]}, Node{xs[5]}}};
    std::array<Node*, MeshElement::n_all_nodes> nodes{
        {&ns[0], &ns[1], &ns[2], &ns[3], &ns[4], &ns[5]}};
    MeshElement element(nodes);

    bool const test_result3 = test<ShapeFunction, 3>(element);
    ASSERT_TRUE(test_result3) << "Test failed for dimension 3.";
}

TEST(NumLibNaturalCoordinates, Prism15)
{
    using ShapeFunction = NumLib::ShapePrism15;
    using MeshElement = ShapeFunction::MeshElement;

    auto const& xs = NumLib::NaturalCoordinates<MeshElement>::coordinates;
    std::array<Node, MeshElement::n_all_nodes> ns{
        {Node{xs[0]}, Node{xs[1]}, Node{xs[2]}, Node{xs[3]}, Node{xs[4]},
         Node{xs[5]}, Node{xs[6]}, Node{xs[7]}, Node{xs[8]}, Node{xs[9]},
         Node{xs[10]}, Node{xs[11]}, Node{xs[12]}, Node{xs[13]}, Node{xs[14]}}};
    std::array<Node*, MeshElement::n_all_nodes> nodes{
        {&ns[0], &ns[1], &ns[2], &ns[3], &ns[4], &ns[5], &ns[6], &ns[7], &ns[8],
         &ns[9], &ns[10], &ns[11], &ns[12], &ns[13], &ns[14]}};
    MeshElement element(nodes);

    bool const test_result3 = test<ShapeFunction, 3>(element);
    ASSERT_TRUE(test_result3) << "Test failed for dimension 3.";
}

TEST(NumLibNaturalCoordinates, Pyramid5)
{
    using ShapeFunction = NumLib::ShapePyra5;
    using MeshElement = ShapeFunction::MeshElement;

    auto const& xs = NumLib::NaturalCoordinates<MeshElement>::coordinates;
    std::array<Node, MeshElement::n_all_nodes> ns{
        {Node{xs[0]}, Node{xs[1]}, Node{xs[2]}, Node{xs[3]}, Node{xs[4]}}};
    std::array<Node*, MeshElement::n_all_nodes> nodes{
        {&ns[0], &ns[1], &ns[2], &ns[3], &ns[4]}};
    MeshElement element(nodes);

    bool const test_result3 = test<ShapeFunction, 3>(element);
    ASSERT_TRUE(test_result3) << "Test failed for dimension 3.";
}

TEST(NumLibNaturalCoordinates, Pyramid13)
{
    using ShapeFunction = NumLib::ShapePyra13;
    using MeshElement = ShapeFunction::MeshElement;

    auto const& xs = NumLib::NaturalCoordinates<MeshElement>::coordinates;
    std::array<Node, MeshElement::n_all_nodes> ns{
        {Node{xs[0]}, Node{xs[1]}, Node{xs[2]}, Node{xs[3]}, Node{xs[4]},
         Node{xs[5]}, Node{xs[6]}, Node{xs[7]}, Node{xs[8]}, Node{xs[9]},
         Node{xs[10]}, Node{xs[11]}, Node{xs[12]}}};

    std::array<Node*, MeshElement::n_all_nodes> nodes{
        {&ns[0], &ns[1], &ns[2], &ns[3], &ns[4], &ns[5], &ns[6], &ns[7], &ns[8],
         &ns[9], &ns[10], &ns[11], &ns[12]}};
    MeshElement element(nodes);

    bool const test_result3 = test<ShapeFunction, 3>(element);
    ASSERT_TRUE(test_result3) << "Test failed for dimension 3.";
}

TEST(NumLibNaturalCoordinates, Hex8)
{
    using ShapeFunction = NumLib::ShapeHex8;
    using MeshElement = ShapeFunction::MeshElement;

    auto const& xs = NumLib::NaturalCoordinates<MeshElement>::coordinates;
    std::array<Node, MeshElement::n_all_nodes> ns{
        {Node{xs[0]}, Node{xs[1]}, Node{xs[2]}, Node{xs[3]}, Node{xs[4]},
         Node{xs[5]}, Node{xs[6]}, Node{xs[7]}}};
    std::array<Node*, MeshElement::n_all_nodes> nodes{
        {&ns[0], &ns[1], &ns[2], &ns[3], &ns[4], &ns[5], &ns[6], &ns[7]}};
    MeshElement element(nodes);

    bool const test_result3 = test<ShapeFunction, 3>(element);
    ASSERT_TRUE(test_result3) << "Test failed for dimension 3.";
}

TEST(NumLibNaturalCoordinates, Hex20)
{
    using ShapeFunction = NumLib::ShapeHex20;
    using MeshElement = ShapeFunction::MeshElement;

    auto const& xs = NumLib::NaturalCoordinates<MeshElement>::coordinates;
    std::array<Node, MeshElement::n_all_nodes> ns{
        {Node{xs[0]},  Node{xs[1]},  Node{xs[2]},  Node{xs[3]},  Node{xs[4]},
         Node{xs[5]},  Node{xs[6]},  Node{xs[7]},  Node{xs[8]},  Node{xs[9]},
         Node{xs[10]}, Node{xs[11]}, Node{xs[12]}, Node{xs[13]}, Node{xs[14]},
         Node{xs[15]}, Node{xs[16]}, Node{xs[17]}, Node{xs[18]}, Node{xs[19]}}};

    std::array<Node*, MeshElement::n_all_nodes> nodes{
        {&ns[0],  &ns[1],  &ns[2],  &ns[3],  &ns[4],  &ns[5],  &ns[6],
         &ns[7],  &ns[8],  &ns[9],  &ns[10], &ns[11], &ns[12], &ns[13],
         &ns[14], &ns[15], &ns[16], &ns[17], &ns[18], &ns[19]}};

    MeshElement element(nodes);

    bool const test_result3 = test<ShapeFunction, 3>(element);
    ASSERT_TRUE(test_result3) << "Test failed for dimension 3.";
}