#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;

int main(int, char **)
{
    cout << "======================================================" << endl;
    cout << "|      DIFFERENTIATION IN A 1D STANDARD REGION       |" << endl;
    cout << "======================================================" << endl;

    cout << "Differentiate the function f(xi) = xi^7 in the " << endl;
    cout << "standard segment xi=[-1,1] using quadrature points" << endl;

    // Specify the number of quadrature points
    int nQuadPoints = 7;

    // Specify the type of quadrature points. This is done using the proper
    // Nektar++ syntax.
    LibUtilities::PointsType quadPointsType = LibUtilities::eGaussGaussLegendre;

    // Declare variables (of type Array) to hold the quadrature zeros
    // and the values of the derivative at the quadrature points
    Array<OneD, NekDouble> quadZeros(nQuadPoints);
    Array<OneD, NekDouble> quadDerivs(nQuadPoints);

    // Declare a pointer (to type NekMatrix<NekDouble>) to hold the
    // differentiation matrix
    DNekMatSharedPtr derivMatrix;

    // Calculate the quadrature zeros and the differentiation matrix.
    // This is done in 2 steps.

    // Step 1: Declare a PointsKey which uniquely defines the
    // quadrature points
    const LibUtilities::PointsKey quadPointsKey(nQuadPoints, quadPointsType);

    // Step 2: Using this key, the quadrature zeros and the differentiation
    // matrix can now be retrieved through the PointsManager in namespace
    // LibUtilities
    quadZeros   = (LibUtilities::PointsManager()[quadPointsKey])->GetZ();
    derivMatrix = (LibUtilities::PointsManager()[quadPointsKey])->GetD();

    // Now you have the quadrature zeros and the differentiation matrix,
    // apply the Gaussian quadrature technique to differentiate the function
    // f(xi) = xi^7 in the standard segment xi=[-1,1].  To do so, write a
    // loop which performs the summation for each quadrature point.
    //
    // In C++, a loop is implemented as:
    //     for(i = min; i < max; i++)
    //     {
    //         // do something
    //     }
    //
    // The function f(xi) can be evaluated using the command
    // 'pow(xi,7)' of the cmath standard library
    //
    // Store the solution in the array 'quadDerivs'

    //                 ==> Write your code here <==

    // Compute the total error
    NekDouble error = 0.0;
	for (size_t i = 0; i < nQuadPoints; ++i)
	{
		error += fabs(quadDerivs[i] - 7 * pow(quadZeros[i], 6));
	}

    // Display the output
    cout << "\t Q = " << nQuadPoints << ": Error = " << error << endl;

    // Now evaluate the derivatives for a quadrature order of Q = Q_max
    // where Q_max is the number of quadrature points required for an
    // exact evaluation of the derivative (calculate this value
    // analytically).  Check that the error should then be zero (up to
    // numerical precision).
}