#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 2D ELEMENT in Local Region |" << endl;
cout << "=========================================================" << endl;
cout << endl;
cout << "Differentiate the function f(x1,x2) = x1^7 * x2^9 " << endl;
cout << "in a local quadrilateral element:" << endl;
// Specify the number of quadrature points in both directions
int nQuadPointsDir1 = 8;
int nQuadPointsDir2 = 10;
// 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 derivatives 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
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();
// The local (straight-sided) quadrilateral element has the
// following vertices:
//
// - Vertex A: (x1_A,x2_A) = (0,-1)
// - Vertex B: (x1_A,x2_A) = (1,-1)
// - Vertex C: (x1_A,x2_A) = (1,1)
// - Vertex D: (x1_A,x2_A) = (0,0)
//
NekDouble x1_A = 0.0;
NekDouble x2_A = -1.0;
NekDouble x1_B = 1.0;
NekDouble x2_B = -1.0;
NekDouble x1_C = 1.0;
NekDouble x2_C = 1.0;
NekDouble x1_D = 0.0;
NekDouble x2_D = 0.0;
// Differentiate the function f(x1,x2) = x1^7 * x2^9 in a local
// quadrilateral element (defined above). Use Gauss-Lobatto-Legendre
// quadrature in both directions.
//
// Your code can be based on the previous exercise. However,
// as we are calculating the derivatives of a function defined in
// a local element rather than in a reference element, we have
// to take into account the geometry of the element.
//
// Therefore, the implementation should be altered in two ways:
//
// (1) The quadrature zeros should be transformed to local
// coordinates to evaluate the function f(x1,x2)
//
// (2) Take into account the Jacobian matrix of the transformation
// between local and reference coordinates when evaluating
// the derivative. (Evaluate the expression for the Jacobian
// matrix analytically rather than using numerical
// differentiation and invert it)
//
// Store the solution in the matrices 'quadDerivsDir1' and 'quadDerivsDir2'
// Apply the Gaussian quadrature technique to differentiate the
// function f(x1,x2) = x1^7 * x2^9 in the standard
// quadrilateral. To do so, edit the (double) loop which performs
// the summation.
NekDouble x1_master, x2_master, x1_slave, x2_slave;
NekDouble dx1dxi1, dx1dxi2, dx2dxi1, dx2dxi2;
NekDouble dxi1dx1, dxi1dx2, dxi2dx1, dxi2dx2;
NekDouble jacobian;
NekDouble error = 0;
for (size_t i = 0; i < nQuadPointsDir1; ++i)
{
for (size_t j = 0; j < nQuadPointsDir2; ++j)
{
// Compute the local coordinates of the quadrature point
x1_master =
x1_A * 0.25 * (1 - quadZerosDir1[i]) * (1 - quadZerosDir2[j]) +
x1_B * 0.25 * (1 + quadZerosDir1[i]) * (1 - quadZerosDir2[j]) +
x1_D * 0.25 * (1 - quadZerosDir1[i]) * (1 + quadZerosDir2[j]) +
x1_C * 0.25 * (1 + quadZerosDir1[i]) * (1 + quadZerosDir2[j]);
x2_master =
x2_A * 0.25 * (1 - quadZerosDir1[i]) * (1 - quadZerosDir2[j]) +
x2_B * 0.25 * (1 + quadZerosDir1[i]) * (1 - quadZerosDir2[j]) +
x2_D * 0.25 * (1 - quadZerosDir1[i]) * (1 + quadZerosDir2[j]) +
x2_C * 0.25 * (1 + quadZerosDir1[i]) * (1 + quadZerosDir2[j]);
// Fill the Jacobian matrix analytically with the dx?dxi? terms
#if WITHSOLUTION
dx1dxi1 = 0.25 * (1 - quadZerosDir2[j]) * (x1_B - x1_A) +
0.25 * (1 + quadZerosDir2[j]) * (x1_C - x1_D);
dx2dxi2 = 0.25 * (1 - quadZerosDir1[i]) * (x2_D - x2_A) +
0.25 * (1 + quadZerosDir1[i]) * (x2_C - x2_B);
dx1dxi2 = 0.25 * (1 - quadZerosDir1[i]) * (x1_D - x1_A) +
0.25 * (1 + quadZerosDir1[i]) * (x1_C - x1_B);
dx2dxi1 = 0.25 * (1 - quadZerosDir2[j]) * (x2_B - x2_A) +
0.25 * (1 + quadZerosDir2[j]) * (x2_C - x2_D);
#endif
// Compute the Jacobian determinant
jacobian = dx1dxi1 * dx2dxi2 - dx1dxi2 * dx2dxi1;
// Invert the Jacobian matrix to obtain the dxi?dx? terms
#if WITHSOLUTION
dxi1dx1 = dx2dxi2 / jacobian;
dxi2dx2 = dx1dxi1 / jacobian;
dxi1dx2 = -dx1dxi2 / jacobian;
dxi2dx1 = -dx2dxi1 / jacobian;
#endif
for (size_t k = 0; k < nQuadPointsDir1; ++k)
{
// Compute the local coordinates of the quadrature point
x1_slave = x1_A * 0.25 * (1 - quadZerosDir1[k]) *
(1 - quadZerosDir2[j]) +
x1_B * 0.25 * (1 + quadZerosDir1[k]) *
(1 - quadZerosDir2[j]) +
x1_D * 0.25 * (1 - quadZerosDir1[k]) *
(1 + quadZerosDir2[j]) +
x1_C * 0.25 * (1 + quadZerosDir1[k]) *
(1 + quadZerosDir2[j]);
x2_slave = x2_A * 0.25 * (1 - quadZerosDir1[k]) *
(1 - quadZerosDir2[j]) +
x2_B * 0.25 * (1 + quadZerosDir1[k]) *
(1 - quadZerosDir2[j]) +
x2_D * 0.25 * (1 - quadZerosDir1[k]) *
(1 + quadZerosDir2[j]) +
x2_C * 0.25 * (1 + quadZerosDir1[k]) *
(1 + quadZerosDir2[j]);
// Add the contribution of this quadrature point to
// quadDerivsDir1[i][j] and quadDerivsDir2[i][j]
#if WITHSOLUTION
quadDerivsDir1[i][j] += (*derivMatrixDir1)(i, k) *
pow(x1_slave, 7) * pow(x2_slave, 9) *
dxi1dx1;
quadDerivsDir2[i][j] += (*derivMatrixDir1)(i, k) *
pow(x1_slave, 7) * pow(x2_slave, 9) *
dxi1dx2;
#endif
}
for (size_t k = 0; k < nQuadPointsDir2; ++k)
{
// Compute the local coordinates of the quadrature point
x1_slave = x1_A * 0.25 * (1 - quadZerosDir1[i]) *
(1 - quadZerosDir2[k]) +
x1_B * 0.25 * (1 + quadZerosDir1[i]) *
(1 - quadZerosDir2[k]) +
x1_D * 0.25 * (1 - quadZerosDir1[i]) *
(1 + quadZerosDir2[k]) +
x1_C * 0.25 * (1 + quadZerosDir1[i]) *
(1 + quadZerosDir2[k]);
x2_slave = x2_A * 0.25 * (1 - quadZerosDir1[i]) *
(1 - quadZerosDir2[k]) +
x2_B * 0.25 * (1 + quadZerosDir1[i]) *
(1 - quadZerosDir2[k]) +
x2_D * 0.25 * (1 - quadZerosDir1[i]) *
(1 + quadZerosDir2[k]) +
x2_C * 0.25 * (1 + quadZerosDir1[i]) *
(1 + quadZerosDir2[k]);
// Add the contribution of this quadrature point to
// quadDerivsDir1[i][j] and quadDerivsDir2[i][j]
#if WITHSOLUTION
quadDerivsDir2[i][j] += (*derivMatrixDir2)(j, k) *
pow(x1_slave, 7) * pow(x2_slave, 9) *
dxi2dx2;
quadDerivsDir1[i][j] += (*derivMatrixDir2)(j, k) *
pow(x1_slave, 7) * pow(x2_slave, 9) *
dxi2dx1;
#endif
}
error += fabs(quadDerivsDir1[i][j] -
7 * pow(x1_master, 6) * pow(x2_master, 9));
error += fabs(quadDerivsDir2[i][j] -
9 * pow(x1_master, 7) * pow(x2_master, 8));
}
}
// Display the averag error
cout << "\t q1 = " << nQuadPointsDir1 << ", q2 = " << nQuadPointsDir2;
cout << ": Average Error = ";
cout << error / (nQuadPointsDir1 * nQuadPointsDir2) << endl;
cout << endl;
}