StdDifferentiation2D.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 IN 2D ELEMENT in Standard Region     |"<< endl;
    cout << "==========================================================="<< endl;
    cout << endl;
    cout << "Differentiate the function f(x1,x2) = (x1)^7*(x2)^9" << endl;
    cout << "in the standard quadrilateral element:" << endl;
    
    // Specify the number of quadrature points in both directions
    int nQuadPointsDir1 = 7;
    int nQuadPointsDir2 = 9;
    
    // 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 derivative 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
	MatrixSharedPtrType derivMatrixDir1;
	MatrixSharedPtrType 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();
    
    // Now you have the quadrature zeros and the differentiation matrix,
	// apply the Gaussian quadrature technique to differentiate the function
    // f(x_1,i,x_2,j) = x_1,i^7 * x2,j^9 on the standard
    // quadrilateral.  To do so, write a (double) loop which performs
    // the summation.
    //
    // Store the solution in the variable 'quadDerivsDir1' and 'quadDerivsDir2'
    
#if WITHSOLUTION
	for (size_t i = 0; i < nQuadPointsDir1; ++i)
		for (size_t j = 0; j < nQuadPointsDir1; ++j)
			for (size_t k = 0; k < nQuadPointsDir2; ++k)
				quadDerivsDir1[i][k] += (*derivMatrixDir1)(i, j) * pow(quadZerosDir1[j], 7) * pow(quadZerosDir2[k], 9);
	for (size_t i = 0; i < nQuadPointsDir2; ++i)
		for (size_t j = 0; j < nQuadPointsDir2; ++j)
			for (size_t k = 0; k < nQuadPointsDir1; ++k)
				quadDerivsDir2[k][i] += (*derivMatrixDir2)(i, j) * pow(quadZerosDir1[k], 7) * pow(quadZerosDir2[j], 9);
#endif
	
	// Compute the total error
	NekDouble error = 0.0;
	for (size_t i = 0; i < nQuadPointsDir1; ++i)
		for (size_t j = 0; j < nQuadPointsDir2; ++j)
			error += fabs(quadDerivsDir1[i][j] - 7 * pow(quadZerosDir1[i], 6) * pow(quadZerosDir2[j], 9));
	for (size_t i = 0; i < nQuadPointsDir1; ++i)
		for (size_t j = 0; j < nQuadPointsDir2; ++j)
			error += fabs(quadDerivsDir2[i][j] - 9 * pow(quadZerosDir1[i], 7) * pow(quadZerosDir2[j], 8));

    // Display the output
    cout << "\t q1 = " << nQuadPointsDir1 << ", q2 = " << nQuadPointsDir2 ;
    cout << ": Error = "<< error << endl;
    cout << endl;
}