#include <cstdio>
#include <cstdlib>

#include <LibUtilities/Foundations/ManagerAccess.h>
#include <LibUtilities/Polylib/Polylib.h>
#include <LibUtilities/LinearAlgebra/NekTypeDefs.hpp>

using namespace Nektar;
using namespace std;

#define WITHSOLUTION 1

int main(int, char **)
{
    cout << "=========================================================" << endl;
    cout << "|     DIFFERENTIATION IN 2D ELEMENT in Local Region     |" << endl;
    cout << "=========================================================" << endl;
    cout << endl;
    cout << "Differentiate the function f(x1,x2) = x1^7 * x2^9 " << endl;
    cout << "in a local quadrilateral element:" << endl;

    // Specify the number of quadrature points in both directions
    int nQuadPointsDir1 = 8;
    int nQuadPointsDir2 = 10;

    // Specify the type of quadrature points in both directions
    LibUtilities::PointsType quadPointsTypeDir1 =
        LibUtilities::eGaussLobattoLegendre;
    LibUtilities::PointsType quadPointsTypeDir2 =
        LibUtilities::eGaussLobattoLegendre;

    // Declare variables (of type Array) to hold the quadrature zeros
    // and the values of the derivatives at the quadrature points in both
    // directions
    Array<OneD, NekDouble> quadZerosDir1(nQuadPointsDir1);
    Array<OneD, NekDouble> quadZerosDir2(nQuadPointsDir2);
    Array<TwoD, NekDouble> quadDerivsDir1(nQuadPointsDir1, nQuadPointsDir2);
    Array<TwoD, NekDouble> quadDerivsDir2(nQuadPointsDir1, nQuadPointsDir2);

    // Declare pointers (to type NekMatrix<NekDouble>) to hold the
    // differentiation matrices
    DNekMatSharedPtr derivMatrixDir1;
    DNekMatSharedPtr derivMatrixDir2;

    // Calculate the GLL-quadrature zeros and the differentiation
    // matrices in both directions. This is done in 2 steps.

    // Step 1: Declare the PointsKeys which uniquely defines the
    // quadrature points
    const LibUtilities::PointsKey quadPointsKeyDir1(nQuadPointsDir1,
                                                    quadPointsTypeDir1);
    const LibUtilities::PointsKey quadPointsKeyDir2(nQuadPointsDir2,
                                                    quadPointsTypeDir2);

    // Step 2: Using this key, the quadrature zeros and differentiation
    // matrices can now be retrieved through the PointsManager
    quadZerosDir1   = LibUtilities::PointsManager()[quadPointsKeyDir1]->GetZ();
    derivMatrixDir1 = LibUtilities::PointsManager()[quadPointsKeyDir1]->GetD();

    quadZerosDir2   = LibUtilities::PointsManager()[quadPointsKeyDir2]->GetZ();
    derivMatrixDir2 = LibUtilities::PointsManager()[quadPointsKeyDir2]->GetD();

    // The local (straight-sided) quadrilateral element has the
    // following vertices:
    //
    // - Vertex A: (x1_A,x2_A) = (0,-1)
    // - Vertex B: (x1_A,x2_A) = (1,-1)
    // - Vertex C: (x1_A,x2_A) = (1,1)
    // - Vertex D: (x1_A,x2_A) = (0,0)
    //
    NekDouble x1_A = 0.0;
    NekDouble x2_A = -1.0;

    NekDouble x1_B = 1.0;
    NekDouble x2_B = -1.0;

    NekDouble x1_C = 1.0;
    NekDouble x2_C = 1.0;

    NekDouble x1_D = 0.0;
    NekDouble x2_D = 0.0;

    // Differentiate the function f(x1,x2) = x1^7 * x2^9 in a local
    // quadrilateral element (defined above). Use Gauss-Lobatto-Legendre
    // quadrature in both directions.
    //
    // Your code can be based on the previous exercise. However,
    // as we are calculating the derivatives of a function defined in
    // a local element rather than in a reference element, we have
    // to take into account the geometry of the element.
    //
    // Therefore, the implementation should be altered in two ways:
    //
    // (1) The quadrature zeros should be transformed to local
    //     coordinates to evaluate the function f(x1,x2)
    //
    // (2) Take into account the Jacobian matrix of the transformation
    //     between local and reference coordinates when evaluating
    //     the derivative. (Evaluate the expression for the Jacobian
    //     matrix analytically rather than using numerical
    //     differentiation and invert it)
    //
    // Store the solution in the matrices 'quadDerivsDir1' and 'quadDerivsDir2'

    // Apply the Gaussian quadrature technique to differentiate the
    // function f(x1,x2) = x1^7 * x2^9 in the standard
    // quadrilateral.  To do so, edit the (double) loop which performs
    // the summation.
    NekDouble x1_master, x2_master, x1_slave, x2_slave;
    NekDouble dx1dxi1, dx1dxi2, dx2dxi1, dx2dxi2;
    NekDouble dxi1dx1, dxi1dx2, dxi2dx1, dxi2dx2;
    NekDouble jacobian;
    NekDouble error = 0;

    for (size_t i = 0; i < nQuadPointsDir1; ++i)
    {
        for (size_t j = 0; j < nQuadPointsDir2; ++j)
        {
            // Compute the local coordinates of the quadrature point
            x1_master =
                x1_A * 0.25 * (1 - quadZerosDir1[i]) * (1 - quadZerosDir2[j]) +
                x1_B * 0.25 * (1 + quadZerosDir1[i]) * (1 - quadZerosDir2[j]) +
                x1_D * 0.25 * (1 - quadZerosDir1[i]) * (1 + quadZerosDir2[j]) +
                x1_C * 0.25 * (1 + quadZerosDir1[i]) * (1 + quadZerosDir2[j]);
            x2_master =
                x2_A * 0.25 * (1 - quadZerosDir1[i]) * (1 - quadZerosDir2[j]) +
                x2_B * 0.25 * (1 + quadZerosDir1[i]) * (1 - quadZerosDir2[j]) +
                x2_D * 0.25 * (1 - quadZerosDir1[i]) * (1 + quadZerosDir2[j]) +
                x2_C * 0.25 * (1 + quadZerosDir1[i]) * (1 + quadZerosDir2[j]);

            // Fill the Jacobian matrix analytically with the dx?dxi? terms
#if WITHSOLUTION
            dx1dxi1 = 0.25 * (1 - quadZerosDir2[j]) * (x1_B - x1_A) +
                      0.25 * (1 + quadZerosDir2[j]) * (x1_C - x1_D);
            dx2dxi2 = 0.25 * (1 - quadZerosDir1[i]) * (x2_D - x2_A) +
                      0.25 * (1 + quadZerosDir1[i]) * (x2_C - x2_B);
            dx1dxi2 = 0.25 * (1 - quadZerosDir1[i]) * (x1_D - x1_A) +
                      0.25 * (1 + quadZerosDir1[i]) * (x1_C - x1_B);
            dx2dxi1 = 0.25 * (1 - quadZerosDir2[j]) * (x2_B - x2_A) +
                      0.25 * (1 + quadZerosDir2[j]) * (x2_C - x2_D);
#endif

            // Compute the Jacobian determinant
            jacobian = dx1dxi1 * dx2dxi2 - dx1dxi2 * dx2dxi1;

            // Invert the Jacobian matrix to obtain the dxi?dx? terms
#if WITHSOLUTION
            dxi1dx1 = dx2dxi2 / jacobian;
            dxi2dx2 = dx1dxi1 / jacobian;
            dxi1dx2 = -dx1dxi2 / jacobian;
            dxi2dx1 = -dx2dxi1 / jacobian;
#endif

            for (size_t k = 0; k < nQuadPointsDir1; ++k)
            {
                // Compute the local coordinates of the quadrature point
                x1_slave = x1_A * 0.25 * (1 - quadZerosDir1[k]) *
                               (1 - quadZerosDir2[j]) +
                           x1_B * 0.25 * (1 + quadZerosDir1[k]) *
                               (1 - quadZerosDir2[j]) +
                           x1_D * 0.25 * (1 - quadZerosDir1[k]) *
                               (1 + quadZerosDir2[j]) +
                           x1_C * 0.25 * (1 + quadZerosDir1[k]) *
                               (1 + quadZerosDir2[j]);

                x2_slave = x2_A * 0.25 * (1 - quadZerosDir1[k]) *
                               (1 - quadZerosDir2[j]) +
                           x2_B * 0.25 * (1 + quadZerosDir1[k]) *
                               (1 - quadZerosDir2[j]) +
                           x2_D * 0.25 * (1 - quadZerosDir1[k]) *
                               (1 + quadZerosDir2[j]) +
                           x2_C * 0.25 * (1 + quadZerosDir1[k]) *
                               (1 + quadZerosDir2[j]);

                // Add the contribution of this quadrature point to
                // quadDerivsDir1[i][j] and quadDerivsDir2[i][j]
#if WITHSOLUTION
                quadDerivsDir1[i][j] += (*derivMatrixDir1)(i, k) *
                                        pow(x1_slave, 7) * pow(x2_slave, 9) *
                                        dxi1dx1;
                quadDerivsDir2[i][j] += (*derivMatrixDir1)(i, k) *
                                        pow(x1_slave, 7) * pow(x2_slave, 9) *
                                        dxi1dx2;
#endif
            }
            for (size_t k = 0; k < nQuadPointsDir2; ++k)
            {
                // Compute the local coordinates of the quadrature point
                x1_slave = x1_A * 0.25 * (1 - quadZerosDir1[i]) *
                               (1 - quadZerosDir2[k]) +
                           x1_B * 0.25 * (1 + quadZerosDir1[i]) *
                               (1 - quadZerosDir2[k]) +
                           x1_D * 0.25 * (1 - quadZerosDir1[i]) *
                               (1 + quadZerosDir2[k]) +
                           x1_C * 0.25 * (1 + quadZerosDir1[i]) *
                               (1 + quadZerosDir2[k]);

                x2_slave = x2_A * 0.25 * (1 - quadZerosDir1[i]) *
                               (1 - quadZerosDir2[k]) +
                           x2_B * 0.25 * (1 + quadZerosDir1[i]) *
                               (1 - quadZerosDir2[k]) +
                           x2_D * 0.25 * (1 - quadZerosDir1[i]) *
                               (1 + quadZerosDir2[k]) +
                           x2_C * 0.25 * (1 + quadZerosDir1[i]) *
                               (1 + quadZerosDir2[k]);

                // Add the contribution of this quadrature point to
                // quadDerivsDir1[i][j] and quadDerivsDir2[i][j]
#if WITHSOLUTION
                quadDerivsDir2[i][j] += (*derivMatrixDir2)(j, k) *
                                        pow(x1_slave, 7) * pow(x2_slave, 9) *
                                        dxi2dx2;
                quadDerivsDir1[i][j] += (*derivMatrixDir2)(j, k) *
                                        pow(x1_slave, 7) * pow(x2_slave, 9) *
                                        dxi2dx1;
#endif
            }

            error += fabs(quadDerivsDir1[i][j] -
                          7 * pow(x1_master, 6) * pow(x2_master, 9));
            error += fabs(quadDerivsDir2[i][j] -
                          9 * pow(x1_master, 7) * pow(x2_master, 8));
        }
    }

    // Display the averag error
    cout << "\t q1 = " << nQuadPointsDir1 << ", q2 = " << nQuadPointsDir2;
    cout << ": Average Error = ";
    cout << error / (nQuadPointsDir1 * nQuadPointsDir2) << endl;
    cout << endl;
}