Newer
Older
#include <cstdio>
#include <cstdlib>
#include <LibUtilities/Polylib/Polylib.h>
#include <LibUtilities/Foundations/ManagerAccess.h>
using namespace Nektar;
using namespace std;
typedef boost::shared_ptr<NekMatrix<NekDouble> > MatrixSharedPtrType;
#define WITHSOLUTION 1
int main(int argc, char *argv[])
{
cout << "==========================================================="<< endl;
cout << "| DIFFERENTIATION ON 2D ELEMENT in Standard Region |"<< endl;
cout << "==========================================================="<< endl;
cout << endl;
cout << "Differentiate the function f(x1,x2) = (x1)^7*(x2)^9" << endl;
cout << "on the standard 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 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
MatrixSharedPtrType derivMatrixDir1;
MatrixSharedPtrType 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 variable 'quadDerivsDir1' and 'quadDerivsDir2'
#if WITHSOLUTION
for (size_t i = 0; i < nQuadPointsDir1; ++i)
for (size_t j = 0; j < nQuadPointsDir1; ++j)
for (size_t k = 0; k < nQuadPointsDir2; ++k)
quadDerivsDir1[i][k] += (*derivMatrixDir1)(i, j) * pow(quadZerosDir1[j], 7) * pow(quadZerosDir2[k], 9);
for (size_t i = 0; i < nQuadPointsDir2; ++i)
for (size_t j = 0; j < nQuadPointsDir2; ++j)
for (size_t k = 0; k < nQuadPointsDir1; ++k)
quadDerivsDir2[k][i] += (*derivMatrixDir2)(i, j) * pow(quadZerosDir1[k], 7) * pow(quadZerosDir2[j], 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));
for (size_t i = 0; i < nQuadPointsDir1; ++i)
for (size_t j = 0; j < nQuadPointsDir2; ++j)
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;
}