Newer
Older
#include <cstdio>
#include <cstdlib>
#include <LibUtilities/Foundations/ManagerAccess.h>
Jeremy Cohen
committed
#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 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
Jeremy Cohen
committed
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();
// 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 matrices 'quadDerivsDir1' and 'quadDerivsDir2'
#if WITHSOLUTION
for (size_t i = 0; i < nQuadPointsDir1; ++i)
{
for (size_t j = 0; j < nQuadPointsDir2; ++j)
{
for (size_t k = 0; k < nQuadPointsDir1; ++k)
{
quadDerivsDir1[i][j] += (*derivMatrixDir1)(i, k) *
pow(quadZerosDir1[k], 7) *
pow(quadZerosDir2[j], 9);
}
for (size_t k = 0; k < nQuadPointsDir2; ++k)
{
quadDerivsDir2[i][j] += (*derivMatrixDir2)(j, k) *
pow(quadZerosDir1[i], 7) *
pow(quadZerosDir2[k], 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));
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;
}