StdDifferentiation1D.cpp

#include <cstdio>
#include <cstdlib>

#include <LibUtilities/Polylib/Polylib.h>
#include <LibUtilities/Foundations/ManagerAccess.h>

using namespace Nektar;
using namespace std;

typedef boost::shared_ptr<NekMatrix<NekDouble> > MatrixSharedPtrType;

#define WITHSOLUTION 1

int main(int argc, char *argv[])
{
    cout << "======================================================" <<endl;
    cout << "|      DIFFERENTIATION ON A 1D STANDARD REGION       |" <<endl;
    cout << "======================================================" <<endl;
    
    cout << "Differentiate the function f(xi) = xi^7 on 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
	MatrixSharedPtrType 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 on 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).
}