Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
#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 argc, char *argv[])
{
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;
}