StdDifferentiation1D.cpp

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

#if WITHSOLUTION
    for (size_t i = 0; i < nQuadPoints; ++i)
    {
        for (size_t j = 0; j < nQuadPoints; ++j)
        {
            quadDerivs[i] += (*derivMatrix)(i, j) * pow(quadZeros[j], 7);
        }
    }
#endif

    // 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).
}