#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 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 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; }