Commit 8f9b8cd8 authored by Spencer Sherwin's avatar Spencer Sherwin
Browse files

Introduced a regression check to see if routine to map equispaced points to...

Introduced a regression check to see if routine to map equispaced points to coefficients works correctly.
parent fc436cad
......@@ -36,7 +36,9 @@ v4.3.0
- Add module to enable mean mode of 3DH1D to be extracted (!530)
- Fix bug in C^0 projection (!541))
- Add command line option to set number of homogeneous planes (!540)
- Add module to project set of points to a fld file
- Add support for interpolating to a box of points and fix ability to run interppointstofld module in parallel when using a plane or box option
**MeshConvert:**
- Enable face curvature inside core MeshConvert objects (!511)
- Add linearise processing module to remove all curvature from high order
......
......@@ -20,6 +20,10 @@ SET(StdProject_Diff2DSource StdProject_Diff2D.cpp)
ADD_NEKTAR_EXECUTABLE(StdProject_Diff2D demos StdProject_Diff2DSource)
TARGET_LINK_LIBRARIES(StdProject_Diff2D ${LinkLibraries})
SET(StdEquiToCoeff2DSource StdEquiToCoeff2D.cpp)
ADD_NEKTAR_EXECUTABLE(StdEquiToCoeff2D demos StdEquiToCoeff2DSource)
TARGET_LINK_LIBRARIES(StdEquiToCoeff2D ${LinkLibraries})
SET(StdProject3DSource StdProject3D.cpp)
ADD_NEKTAR_EXECUTABLE(StdProject3D demos StdProject3DSource)
TARGET_LINK_LIBRARIES(StdProject3D ${LinkLibraries})
......
#include <cstdio>
#include <cstdlib>
#include <StdRegions/StdExpansion.h>
#include <StdRegions/StdQuadExp.h>
#include <StdRegions/StdTriExp.h>
#include <StdRegions/StdRegions.hpp>
using namespace Nektar;
using namespace StdRegions;
using namespace std;
NekDouble Tri_sol(NekDouble x, NekDouble y, int order1, int order2);
int main(int argc, char *argv[])
{
int i,j;
int order, nq1,nq2;
LibUtilities::PointsType Qtype1,Qtype2;
LibUtilities::BasisType btype1,btype2;
StdExpansion *E;
if(argc != 4)
{
fprintf(stderr,"Usage: EquiToCoeff2D order nq1 nq2 \n");
exit(1);
}
btype1 = (LibUtilities::BasisType) LibUtilities::eModified_A;
btype2 = (LibUtilities::BasisType) LibUtilities::eModified_B;
order = atoi(argv[1]);
nq1 = atoi(argv[2]);
nq2 = atoi(argv[3]);
Qtype1 = LibUtilities::eGaussLobattoLegendre;
Qtype2 = LibUtilities::eGaussRadauMAlpha1Beta0;
//-----------------------------------------------
// Define a segment expansion based on basis definition
const LibUtilities::PointsKey Pkey1(nq1,Qtype1);
const LibUtilities::PointsKey Pkey2(nq2,Qtype2);
const LibUtilities::BasisKey b1(btype1,order,Pkey1);
const LibUtilities::BasisKey b2(btype2,order,Pkey2);
E = new StdTriExp (b1,b2);
Array<OneD,NekDouble> x = Array<OneD,NekDouble>(nq1*nq2);
Array<OneD,NekDouble> y = Array<OneD,NekDouble>(nq1*nq2);
E->GetCoords(x,y);
Array<OneD,NekDouble> sol = Array<OneD,NekDouble>(nq1*nq2);
//----------------------------------------------
// Define solution to be projected
for(i = 0; i < nq1*nq2; ++i)
{
sol[i] = Tri_sol(x[i],y[i],order,order);
}
//----------------------------------------------
Array<OneD, NekDouble> esol(order*(order+1)/2);
Array<OneD, NekDouble> coeffs(order*(order+1)/2);
//---------------------------------------------
// Interpolate to equispaced points
E->PhysInterpToSimplexEquiSpaced(sol,esol,order);
//---------------------------------------------
//--------------------------------------------
// Project from equispaced points back to coeffs
E->EquiSpacedToCoeffs(esol,coeffs);
//--------------------------------------------
Array<OneD,NekDouble> phys = Array<OneD,NekDouble>(nq1*nq2);
//--------------------------------------------
// Project from equispaced points back to coeffs
E->BwdTrans(coeffs,phys);
//--------------------------------------------
//--------------------------------------------
// Calculate L_inf error
cout << "L infinity error: " << E->Linf(phys,sol) << endl;
cout << "L 2 error: " << E->L2 (phys,sol) << endl;
//--------------------------------------------
}
NekDouble Tri_sol(NekDouble x, NekDouble y, int order1, int order2){
int l,k;
NekDouble sol = 0;
for(k = 0; k < order1; ++k)
{
for(l = 0; l < order2-k; ++l)
{
sol += pow(x,k)*pow(y,l);
}
}
return sol;
}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment