#include <cstdio> #include <cstdlib> #include <LibUtilities/Polylib/Polylib.h> #include <LibUtilities/Foundations/ManagerAccess.h> using namespace Nektar; using namespace std; #define WITHSOLUTION 1 int main(int, char **) { cout << "==========================================================="<< endl; cout << "| INTEGRATION ON 2D ELEMENT in Standard Region |"<< endl; cout << "==========================================================="<< endl; cout << endl; cout << "Integrate the function f(x1,x2) = (x1)^12*(x2)^14" << endl; cout << "on the standard quadrilateral element:" << endl; // Specify the number of quadrature points in both directions int nQuadPointsDir1 = 6; int nQuadPointsDir2 = 7; // 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 weights Array<OneD, NekDouble> quadZerosDir1 (nQuadPointsDir1); Array<OneD, NekDouble> quadWeightsDir1(nQuadPointsDir1); Array<OneD, NekDouble> quadZerosDir2 (nQuadPointsDir2); Array<OneD, NekDouble> quadWeightsDir2(nQuadPointsDir2); // Calculate the GLL-quadrature zeros and weights 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 weights // can now be retrieved through the PointsManager quadZerosDir1 = LibUtilities::PointsManager()[quadPointsKeyDir1]->GetZ(); quadWeightsDir1 = LibUtilities::PointsManager()[quadPointsKeyDir1]->GetW(); quadZerosDir2 = LibUtilities::PointsManager()[quadPointsKeyDir2]->GetZ(); quadWeightsDir2 = LibUtilities::PointsManager()[quadPointsKeyDir2]->GetW(); // Now you have the quadrature zeros and weight, apply the // Gaussian quadrature technique to integrate the function // f(x_1,i,x_2,j) = x_1,i^12 * x2,j^14 on the standard // quadrilateral. To do so, write a (double) loop which performs // the summation. // // Store the solution in the variable 'result' NekDouble result = 0.0; #if WITHSOLUTION for(int i = 0; i < nQuadPointsDir1; i++) { for(int j = 0; j < nQuadPointsDir2; j++) { result += pow(quadZerosDir1[i],12)*pow(quadZerosDir2[j],14)* quadWeightsDir1[i]*quadWeightsDir2[j]; } } #endif // Display the output NekDouble exactResult = 4.0/195.0; cout << "\t q1 = " << nQuadPointsDir1 << ", q2 = " << nQuadPointsDir2 ; cout << ": Error = "<< fabs(result-exactResult) << endl; cout << endl; }