Commit 1c2a0c02 authored by Dave Moxey's avatar Dave Moxey
Browse files

Remove unnecessary StdRegions demos

parent db37229b
......@@ -38,47 +38,6 @@ ADD_NEKTAR_EXECUTABLE(StdProject_Diff3D demos StdProject_Diff3DSource)
TARGET_LINK_LIBRARIES(StdProject_Diff3D ${LinkLibraries})
SET_LAPACK_LINK_LIBRARIES(StdProject_Diff3D)
SET(StdTetExpDemoSource StdTetExpDemo.cpp)
ADD_NEKTAR_EXECUTABLE(StdTetExpDemo demos StdTetExpDemoSource)
TARGET_LINK_LIBRARIES(StdTetExpDemo ${LinkLibraries})
SET_LAPACK_LINK_LIBRARIES(StdTetExpDemo)
SET(StdHexExpDemoSource StdHexExpDemo.cpp)
ADD_NEKTAR_EXECUTABLE(StdHexExpDemo demos StdHexExpDemoSource)
TARGET_LINK_LIBRARIES(StdHexExpDemo ${LinkLibraries})
SET_LAPACK_LINK_LIBRARIES(StdHexExpDemo)
SET(StdTetExp_Diff3DSource StdTetExp_Diff3D.cpp)
ADD_NEKTAR_EXECUTABLE(StdTetExp_Diff3D demos StdTetExp_Diff3DSource)
TARGET_LINK_LIBRARIES(StdTetExp_Diff3D ${LinkLibraries})
SET_LAPACK_LINK_LIBRARIES(StdTetExp_Diff3D)
SET(StdHexExp_Diff3DSource StdHexExp_Diff3D.cpp)
ADD_NEKTAR_EXECUTABLE(StdHexExp_Diff3D demos StdHexExp_Diff3DSource)
TARGET_LINK_LIBRARIES(StdHexExp_Diff3D ${LinkLibraries})
SET_LAPACK_LINK_LIBRARIES(StdHexExp_Diff3D)
SET(StdPrismExpDemoSource StdPrismExpDemo.cpp)
ADD_NEKTAR_EXECUTABLE(StdPrismExpDemo demos StdPrismExpDemoSource)
TARGET_LINK_LIBRARIES(StdPrismExpDemo ${LinkLibraries})
SET_LAPACK_LINK_LIBRARIES(StdPrismExpDemo)
SET(StdPrismExp_Diff3DSource StdPrismExp_Diff3D.cpp)
ADD_NEKTAR_EXECUTABLE(StdPrismExp_Diff3D demos StdPrismExp_Diff3DSource)
TARGET_LINK_LIBRARIES(StdPrismExp_Diff3D ${LinkLibraries})
SET_LAPACK_LINK_LIBRARIES(StdPrismExp_Diff3D)
SET(StdPyrExpDemoSource StdPyrExpDemo.cpp)
ADD_NEKTAR_EXECUTABLE(StdPyrExpDemo demos StdPyrExpDemoSource)
TARGET_LINK_LIBRARIES(StdPyrExpDemo ${LinkLibraries})
SET_LAPACK_LINK_LIBRARIES(StdPyrExpDemo)
SET(StdPyrExp_Diff3DSource StdPyrExp_Diff3D.cpp)
ADD_NEKTAR_EXECUTABLE(StdPyrExp_Diff3D demos StdPyrExp_Diff3DSource)
TARGET_LINK_LIBRARIES(StdPyrExp_Diff3D ${LinkLibraries})
SET_LAPACK_LINK_LIBRARIES(StdPyrExp_Diff3D)
# Generate list of available subdirectories
FILE(GLOB dir_list "*")
LIST(REMOVE_ITEM dir_list "${CMAKE_SOURCE_DIR}/.svn")
......@@ -88,7 +47,6 @@ FOREACH(dir ${dir_list})
ENDIF(IS_DIRECTORY ${dir} AND EXISTS ${dir}/CMakeLists.txt)
ENDFOREACH(dir ${dir_list})
ADD_NEKTAR_TEST(StdProject1D_Seg_Orth_P6_Q7)
ADD_NEKTAR_TEST(StdProject1D_Seg_Mod_P6_Q7)
ADD_NEKTAR_TEST(StdProject1D_Single_Mode_Fourier_P2_Q2)
......@@ -126,4 +84,4 @@ ADD_NEKTAR_TEST(StdProject_Diff3D_Hex_Orth_P6_Q7)
ADD_NEKTAR_TEST(StdProject_Diff3D_Hex_Mod_P6_Q7)
ADD_NEKTAR_TEST(StdProject_Diff3D_Hex_Lagrange_P6_Q7)
ADD_NEKTAR_TEST(StdProject_Diff3D_Hex_Legendre_P6_Q7)
ADD_NEKTAR_TEST(StdProject_Diff3D_Hex_Chebyshev_P6_Q7)
\ No newline at end of file
ADD_NEKTAR_TEST(StdProject_Diff3D_Hex_Chebyshev_P6_Q7)
#include <iostream>
#include <cstdlib>
#include <cmath>
using namespace std;
#include <StdRegions/StdHexExp.h>
#include <StdRegions/StdRegions.hpp>
#include <LibUtilities/Foundations/Foundations.hpp>
#include <LibUtilities/Foundations/Basis.h>
#include <LibUtilities/BasicConst/NektarUnivTypeDefs.hpp>
using namespace Nektar;
using namespace Nektar::LibUtilities;
using namespace Nektar::StdRegions;
/// Defines a solution which excites all modes in a StdHex expansion.
NekDouble Hex_sol(NekDouble x, NekDouble y, NekDouble z,
int order1, int order2, int order3,
LibUtilities::BasisType bType_x,
LibUtilities::BasisType bType_y,
LibUtilities::BasisType bType_z);
/// modification to deal with exact solution. Return 1 if integer < 0
static double pow_loc(const double val, const int i)
{
return (i < 0)? 1.0: pow(val,i);
}
/// This routine projects a polynomial or trigonmetric functions which
/// has energy in all mdoes of the expansions and reports and error
int main(int argc, char *argv[]) {
if( argc != 10 ) {
cerr << "Usage: HexDemo Type_x Type_y Type_z numModes_x numModes_y "
<< "numModes_z Qx Qy Qz" << endl;
cerr << "Where type is an interger value which dictates the basis as:"
<< endl;
cerr << "\t Ortho_A = 1\n";
cerr << "\t Ortho_B = 2\n";
cerr << "\t Ortho_C = 3\n";
cerr << "\t Modified_A = 4\n";
cerr << "\t Modified_B = 5\n";
cerr << "\t Modified_C = 6\n";
cerr << "\t Fourier = 7\n";
cerr << "\t Lagrange = 8\n";
cerr << "\t Legendre = 9\n";
cerr << "\t Chebyshev = 10\n";
cerr << "\n\n" << "Example: " << argv[0] << " 1 2 3 3 3 3 5 5 5"
<< endl << endl;
exit(1);
}
LibUtilities::ShapeType regionShape = LibUtilities::eHexahedron;
int bType_x_val = atoi(argv[1]);
int bType_y_val = atoi(argv[2]);
int bType_z_val = atoi(argv[3]);
BasisType bType_x = static_cast<BasisType>( bType_x_val );
BasisType bType_y = static_cast<BasisType>( bType_y_val );
BasisType bType_z = static_cast<BasisType>( bType_z_val );
if( (bType_x == eOrtho_B) || (bType_x == eModified_B) ) {
NEKERROR(ErrorUtil::efatal,
"Basis 1 cannot be of type Ortho_B or Modified_B");
}
if( (bType_x == eOrtho_C) || (bType_x == eModified_C) ) {
NEKERROR(ErrorUtil::efatal,
"Basis 1 cannot be of type Ortho_C or Modified_C");
}
if( (bType_y == eOrtho_C) || (bType_y == eModified_C) ) {
NEKERROR(ErrorUtil::efatal,
"Basis 2 cannot be of type Ortho_C or Modified_C");
}
int xModes = atoi(argv[4]);
int yModes = atoi(argv[5]);
int zModes = atoi(argv[6]);
int Qx = atoi(argv[7]);
int Qy = atoi(argv[8]);
int Qz = atoi(argv[9]);
int P = xModes - 1;
int Q = yModes - 1;
int R = zModes - 1;
LibUtilities::PointsType Qtype_x, Qtype_y, Qtype_z;
Array<OneD, NekDouble> x( Qx * Qy * Qz );
Array<OneD, NekDouble> y( Qx * Qy * Qz );
Array<OneD, NekDouble> z( Qx * Qy * Qz );
Array<OneD, NekDouble> solution( Qx * Qy * Qz );
Array<OneD, NekDouble> diff_solution( Qx * Qy * Qz );
if(bType_x != LibUtilities::eFourier)
{
Qtype_x = LibUtilities::eGaussLobattoLegendre;
}
else
{
Qtype_x = LibUtilities::eFourierEvenlySpaced;
}
if(bType_y != LibUtilities::eFourier)
{
Qtype_y = LibUtilities::eGaussLobattoLegendre;
}
else
{
Qtype_y = LibUtilities::eFourierEvenlySpaced;
}
if(bType_z != LibUtilities::eFourier)
{
Qtype_z = LibUtilities::eGaussLobattoLegendre;
}
else
{
Qtype_z = LibUtilities::eFourierEvenlySpaced;
}
//-----------------------------------------------
// Define a 3D expansion based on basis definition
StdRegions::StdExpansion *she;
const LibUtilities::PointsKey pointsKey_x( Qx, Qtype_x );
const LibUtilities::PointsKey pointsKey_y( Qy, Qtype_y );
const LibUtilities::PointsKey pointsKey_z( Qz, Qtype_z );
const LibUtilities::BasisKey basisKey_x( bType_x, xModes, pointsKey_x );
const LibUtilities::BasisKey basisKey_y( bType_y, yModes, pointsKey_y );
const LibUtilities::BasisKey basisKey_z( bType_z, zModes, pointsKey_z );
she = new StdRegions::StdHexExp( basisKey_x, basisKey_y, basisKey_z );
she->GetCoords(x,y,z);
//----------------------------------------------
// Define solution to be projected
for(int n = 0; n < Qx * Qy * Qz; ++n) {
solution[n] = Hex_sol( x[n], y[n], z[n], P, Q, R,
bType_x, bType_y, bType_z );
}
//----------------------------------------------
Array<OneD, NekDouble> phys (Qx*Qy*Qz);
Array<OneD, NekDouble> coeffs(xModes*yModes*zModes);
//---------------------------------------------
// Project onto Expansion
she->FwdTrans( solution, coeffs );
//---------------------------------------------
//-------------------------------------------
// Backward Transform Solution to get projected values
she->BwdTrans(coeffs, phys);
//-------------------------------------------
//--------------------------------------------
// Calculate L_p error
cout << "L infinity error: " << she->Linf(phys,solution) << endl;
cout << "L 2 error: " << she->L2 (phys,solution) << endl;
//--------------------------------------------
//-------------------------------------------
// Evaulate solution at x = y = z = 0 and print error
Array<OneD, NekDouble> t = Array<OneD, NekDouble>(3);
t[0] = -0.39;
t[1] = -0.25;
t[2] = 0.5;
if( regionShape == LibUtilities::eHexahedron ) {
solution[0] = Hex_sol( t[0], t[1], t[2], P, Q, R,
bType_x, bType_y, bType_z );
}
NekDouble numericSolution = she->PhysEvaluate(t,phys);
cout << "Interpolation difference from actual solution at x = ( "
<< t[0] << ", " << t[1] << ", " << t[2] << " ): "
<< numericSolution - solution[0] << endl;
//-------------------------------------------
cout << "numericSolution = " << numericSolution << endl;
cout << "exact solution = " << solution[0] << endl;
return 0;
}
NekDouble Hex_sol(NekDouble x, NekDouble y, NekDouble z,
int order1, int order2, int order3,
LibUtilities::BasisType btype1,
LibUtilities::BasisType btype2,
LibUtilities::BasisType btype3 ) {
int i,j,k;
NekDouble sol = 0.0;
int Nx = (btype1 == LibUtilities::eFourier ? order1/2 : order1);
int Ny = (btype2 == LibUtilities::eFourier ? order2/2 : order2);
int Nz = (btype3 == LibUtilities::eFourier ? order3/2 : order3);
bool Fx = (btype1 == LibUtilities::eFourier);
bool Fy = (btype2 == LibUtilities::eFourier);
bool Fz = (btype3 == LibUtilities::eFourier);
NekDouble a;
for (i = 0; i < Nx; ++i)
{
for (j = 0; j < Ny; ++j)
{
for (k = 0; k < Nz; ++k)
{
a = (Fx ? sin(M_PI*i*x) + cos(M_PI*i*x) : pow_loc(x,i));
a *= (Fy ? sin(M_PI*j*y) + cos(M_PI*j*y) : pow_loc(y,j));
a *= (Fz ? sin(M_PI*k*z) + cos(M_PI*k*z) : pow_loc(z,k));
sol += a;
}
}
}
return sol;
}
#include <iostream>
#include <cstdlib>
#include <cmath>
using namespace std;
#include <StdRegions/StdHexExp.h>
#include <StdRegions/StdRegions.hpp>
#include <LibUtilities/Foundations/Foundations.hpp>
#include <LibUtilities/Foundations/Basis.h>
#include <LibUtilities/BasicConst/NektarUnivTypeDefs.hpp>
using namespace Nektar;
using namespace Nektar::LibUtilities;
using namespace Nektar::StdRegions;
/// Defines a solution which excites all modes in a Hex expansion.
NekDouble Hex_sol(NekDouble x, NekDouble y, NekDouble z,
int order1, int order2, int order3,
LibUtilities::BasisType bType_x,
LibUtilities::BasisType bType_y,
LibUtilities::BasisType bType_z);
/// Derivative of solution on Hex expansion.
NekDouble Hex_Diff_Sol(NekDouble x, NekDouble y, NekDouble z,
int P, int Q, int R,
LibUtilities::BasisType bType_x,
LibUtilities::BasisType bType_y,
LibUtilities::BasisType bType_z, int dir);
/// modification to deal with exact solution. Return 1 if integer < 0
static double pow_loc(const double val, const int i)
{
return (i < 0)? 1.0: pow(val,i);
}
/// This routine projects a polynomial or trigonmetric functions which
/// has energy in all mdoes of the expansions and reports and error
int main(int argc, char *argv[]) {
if( argc != 10 ) {
cerr << "Usage: HexDemo Type_x Type_y Type_z numModes_x numModes_y "
"numModes_z Qx Qy Qz" << endl;
cerr << "Where type is an interger value which dictates the basis as:"
<< endl;
cerr << "\t Ortho_A = 1\n";
cerr << "\t Ortho_B = 2\n";
cerr << "\t Ortho_C = 3\n";
cerr << "\t Modified_A = 4\n";
cerr << "\t Modified_B = 5\n";
cerr << "\t Modified_C = 6\n";
cerr << "\t Fourier = 7\n";
cerr << "\t Lagrange = 8\n";
cerr << "\t Legendre = 9\n";
cerr << "\t Chebyshev = 10\n";
cerr << "\n\n" << "Example: " << argv[0] << " 1 2 3 3 3 3 5 5 5"
<< endl << endl;
exit(1);
}
int bType_x_val = atoi(argv[1]);
int bType_y_val = atoi(argv[2]);
int bType_z_val = atoi(argv[3]);
BasisType bType_x = static_cast<BasisType>( bType_x_val );
BasisType bType_y = static_cast<BasisType>( bType_y_val );
BasisType bType_z = static_cast<BasisType>( bType_z_val );
if( (bType_x == eOrtho_B) || (bType_x == eModified_B) ) {
NEKERROR(ErrorUtil::efatal,
"Basis 1 cannot be of type Ortho_B or Modified_B");
}
if( (bType_x == eOrtho_C) || (bType_x == eModified_C) ) {
NEKERROR(ErrorUtil::efatal,
"Basis 1 cannot be of type Ortho_C or Modified_C");
}
if( (bType_y == eOrtho_C) || (bType_y == eModified_C) ) {
NEKERROR(ErrorUtil::efatal,
"Basis 2 cannot be of type Ortho_C or Modified_C");
}
int xModes = atoi(argv[4]);
int yModes = atoi(argv[5]);
int zModes = atoi(argv[6]);
int Qx = atoi(argv[7]);
int Qy = atoi(argv[8]);
int Qz = atoi(argv[9]);
int P = xModes - 1;
int Q = yModes - 1;
int R = zModes - 1;
Array<OneD, NekDouble> solution( Qx * Qy * Qz );
Array<OneD, NekDouble> dx(Qx*Qy*Qz,0.0);
Array<OneD, NekDouble> dy(Qx*Qy*Qz,0.0);
Array<OneD, NekDouble> dz(Qx*Qy*Qz,0.0);
Array<OneD, NekDouble> diff_solution_x( Qx * Qy * Qz, 0.0 );
Array<OneD, NekDouble> diff_solution_y( Qx * Qy * Qz, 0.0 );
Array<OneD, NekDouble> diff_solution_z( Qx * Qy * Qz, 0.0 );
Array<OneD, NekDouble> derivatives( Qx * Qy * Qz, 0.0 );
LibUtilities::PointsType Qtype_x, Qtype_y, Qtype_z;
Array<OneD, NekDouble> x( Qx * Qy * Qz );
Array<OneD, NekDouble> y( Qx * Qy * Qz );
Array<OneD, NekDouble> z( Qx * Qy * Qz );
if(bType_x != LibUtilities::eFourier)
{
Qtype_x = LibUtilities::eGaussLobattoLegendre;
}
else
{
Qtype_x = LibUtilities::eFourierEvenlySpaced;
}
if(bType_y != LibUtilities::eFourier)
{
Qtype_y = LibUtilities::eGaussLobattoLegendre;
}
else
{
Qtype_y = LibUtilities::eFourierEvenlySpaced;
}
if(bType_z != LibUtilities::eFourier)
{
Qtype_z = LibUtilities::eGaussLobattoLegendre;
}
else
{
Qtype_z = LibUtilities::eFourierEvenlySpaced;
}
//-----------------------------------------------
// Define a 3D expansion based on basis definition
StdRegions::StdExpansion *she;
const LibUtilities::PointsKey pointsKey_x( Qx, Qtype_x );
const LibUtilities::PointsKey pointsKey_y( Qy, Qtype_y );
const LibUtilities::PointsKey pointsKey_z( Qz, Qtype_z );
const LibUtilities::BasisKey basisKey_x( bType_x, xModes, pointsKey_x );
const LibUtilities::BasisKey basisKey_y( bType_y, yModes, pointsKey_y );
const LibUtilities::BasisKey basisKey_z( bType_z, zModes, pointsKey_z );
she = new StdRegions::StdHexExp( basisKey_x, basisKey_y, basisKey_z );
she->GetCoords(x,y,z);
//----------------------------------------------
// Define solution to be projected
for(int n = 0; n < Qx * Qy * Qz; ++n) {
solution[n] = Hex_sol( x[n], y[n], z[n], P, Q, R,
bType_x, bType_y, bType_z );
}
//----------------------------------------------
//----------------------------------------------
// Define the derivative solution
for(int n = 0; n < Qx*Qy*Qz; ++n){
diff_solution_x[n] = Hex_Diff_Sol( x[n], y[n], z[n], P, Q, R,
bType_x, bType_y, bType_z, 1);
diff_solution_y[n] = Hex_Diff_Sol( x[n], y[n], z[n], P, Q, R,
bType_x, bType_y, bType_z, 2);
diff_solution_z[n] = Hex_Diff_Sol( x[n], y[n], z[n], P, Q, R,
bType_x, bType_y, bType_z, 3);
}
//---------------------------------------------
Array<OneD, NekDouble> phys (Qx*Qy*Qz);
Array<OneD, NekDouble> coeffs(xModes*yModes*zModes);
//---------------------------------------------
// Project onto Expansion
she->FwdTrans( solution, coeffs);
//---------------------------------------------
//-------------------------------------------
// Backward Transform Solution to get projected values
she->BwdTrans( coeffs, phys );
//-------------------------------------------
//--------------------------------------------
// Calculate L_p error
cout << "L infinity error: " << she->Linf(phys,solution) << endl;
cout << "L 2 error: " << she->L2 (phys,solution) << endl;
//--------------------------------------------
//-------------------------------------------
// Evaulate solution at x = y = z = 0 and print error
Array<OneD, NekDouble> t = Array<OneD, NekDouble>(3);
t[0] = -0.39;
t[1] = -0.25;
t[2] = 0.5;
solution[0] = Hex_sol( t[0], t[1], t[2], P, Q, R,
bType_x, bType_y, bType_z );
NekDouble numericSolution = she->PhysEvaluate(t,phys);
cout << "Interpolation difference from actual solution at x = ( "
<< t[0] << ", " << t[1] << ", " << t[2] << " ): "
<< numericSolution - solution[0] << endl;
cout << "numericSolution = " << numericSolution << endl;
cout << "exact solution = " << solution[0] << endl;
//-------------------------------------------
//--------------------------------------------
// Taking the physical derivative and putting them into dx, dy, dz.
she->PhysDeriv( phys, dx, dy, dz );
//--------------------------------------------
double error_x = 0, error_y=0, error_z=0;
for( int n = 0; n < Qx*Qy*Qz; ++n )
{
error_x += fabs(diff_solution_x[n] - dx[n]);
error_y += fabs(diff_solution_y[n] - dy[n]);
error_z += fabs(diff_solution_z[n] - dz[n]);
}
cout << "L 1 error of derivatives X = " << error_x << endl;
cout << "L 1 error of derivatives Y = " << error_y << endl;
cout << "L 1 error of derivatives Z = " << error_z << endl;
return 0;
}
NekDouble Hex_sol(NekDouble x, NekDouble y, NekDouble z,
int order1, int order2, int order3,
LibUtilities::BasisType btype1,
LibUtilities::BasisType btype2,
LibUtilities::BasisType btype3 )
{
int i,j,k;
NekDouble sol = 0.0;
int Nx = (btype1 == LibUtilities::eFourier ? order1/2 : order1);
int Ny = (btype2 == LibUtilities::eFourier ? order2/2 : order2);
int Nz = (btype3 == LibUtilities::eFourier ? order3/2 : order3);
bool Fx = (btype1 == LibUtilities::eFourier);
bool Fy = (btype2 == LibUtilities::eFourier);
bool Fz = (btype3 == LibUtilities::eFourier);
NekDouble a;
for (i = 0; i < Nx; ++i)
{
for (j = 0; j < Ny; ++j)
{
for (k = 0; k < Nz; ++k)
{
a = (Fx ? sin(M_PI*i*x) + cos(M_PI*i*x) : pow_loc(x,i));
a *= (Fy ? sin(M_PI*j*y) + cos(M_PI*j*y) : pow_loc(y,j));
a *= (Fz ? sin(M_PI*k*z) + cos(M_PI*k*z) : pow_loc(z,k));
sol += a;
}
}
}
return sol;
}
NekDouble Hex_Diff_Sol(NekDouble x, NekDouble y, NekDouble z,
int order1, int order2, int order3,
LibUtilities::BasisType btype1,
LibUtilities::BasisType btype2,
LibUtilities::BasisType btype3, int dir)
{
int i,j,k;
NekDouble sol = 0.0;
int Nx = (btype1 == LibUtilities::eFourier ? order1/2 : order1);
int Ny = (btype2 == LibUtilities::eFourier ? order2/2 : order2);
int Nz = (btype3 == LibUtilities::eFourier ? order3/2 : order3);
bool Fx = (btype1 == LibUtilities::eFourier);
bool Fy = (btype2 == LibUtilities::eFourier);
bool Fz = (btype3 == LibUtilities::eFourier);
NekDouble a;
for (i = 0; i < Nx; ++i)
{
for (j = 0; j < Ny; ++j)
{
for (k = 0; k < Nz; ++k)
{
a = i*(Fx ? M_PI*(cos(M_PI*i*x) - sin(M_PI*i*x))
: pow_loc(x,i-1));
a *= (Fy ? (sin(M_PI*j*y) + cos(M_PI*j*y)) : pow_loc(y,j));
a *= (Fz ? (sin(M_PI*k*z) + cos(M_PI*k*z)) : pow_loc(z,k));
sol += a;
a = (Fx ? (sin(M_PI*i*x) + cos(M_PI*i*x)) : pow_loc(x,i));
a *= j*(Fy ? M_PI*(cos(M_PI*j*y) - sin(M_PI*j*y))
: pow_loc(y,j-1));
a *= (Fz ? (sin(M_PI*k*z) + cos(M_PI*k*z)) : pow_loc(z,k));
sol += a;
a = (Fx ? (sin(M_PI*i*x) + cos(M_PI*i*x)) : pow_loc(x,i));
a *= (Fy ? (sin(M_PI*j*y) + cos(M_PI*j*y)) : pow_loc(y,j));
a *= k*(Fz ? M_PI*(cos(M_PI*k*z) - sin(M_PI*k*z))
: pow_loc(z,k-1));
sol += a;