Commit 0cce46c3 authored by Matthias Messner's avatar Matthias Messner
Browse files

added tensor product interpolation for M2M and L2L for the Chebyshev scheme;

Chebyshev roots are now [-1,1] instead of [1,-1]
parent ff588f62
......@@ -33,6 +33,10 @@ class FChebInterpolator : FNoCopyable
unsigned int node_ids[nnodes][3];
FReal* ChildParentInterpolator[8];
// permutations (only needed in the tensor product interpolation case)
unsigned int perm[3][nnodes];
/**
* Initialize the child - parent - interpolator, it is basically the matrix
* S which is precomputed and reused for all M2M and L2L operations, ie for
......@@ -63,6 +67,51 @@ class FChebInterpolator : FNoCopyable
}
}
/**
* Initialize the child - parent - interpolator, it is basically the matrix
* S which is precomputed and reused for all M2M and L2L operations, ie for
* all non leaf inter/anterpolations.
*/
void initTensorM2MandL2L()
{
FPoint ParentRoots[nnodes];
FReal ChildCoords[3][ORDER];
const FReal ParentWidth(2.);
const FPoint ParentCenter(0., 0., 0.);
FChebTensor<ORDER>::setRoots(ParentCenter, ParentWidth, ParentRoots);
FPoint ChildCenter;
const FReal ChildWidth(1.);
// loop: child cells
for (unsigned int child=0; child<8; ++child) {
// set child info
FChebTensor<ORDER>::setRelativeChildCenter(child, ChildCenter);
FChebTensor<ORDER>::setChebyshevRoots(ChildCenter, ChildWidth, ChildCoords);
// allocate memory
ChildParentInterpolator[child] = new FReal [3 * ORDER*ORDER];
assembleInterpolator(ORDER, ChildCoords[0], ChildParentInterpolator[child]);
assembleInterpolator(ORDER, ChildCoords[1], ChildParentInterpolator[child] + 1 * ORDER*ORDER);
assembleInterpolator(ORDER, ChildCoords[2], ChildParentInterpolator[child] + 2 * ORDER*ORDER);
}
// init permutations
for (unsigned int i=0; i<ORDER; ++i) {
for (unsigned int j=0; j<ORDER; ++j) {
for (unsigned int k=0; k<ORDER; ++k) {
const unsigned int index = k*ORDER*ORDER + j*ORDER + i;
perm[0][index] = k*ORDER*ORDER + j*ORDER + i;
perm[1][index] = i*ORDER*ORDER + k*ORDER + j;
perm[2][index] = j*ORDER*ORDER + i*ORDER + k;
}
}
}
}
public:
......@@ -82,7 +131,10 @@ public:
// initialize interpolation operator for non M2M and L2L (non leaf
// operations)
this -> initM2MandL2L();
//this -> initM2MandL2L();
this -> initTensorM2MandL2L();
}
......@@ -145,6 +197,33 @@ public:
}
void assembleInterpolator(const unsigned int M, const FReal *const x, FReal *const S) const
{
// values of chebyshev polynomials of source particle: T_o(x_i)
FReal T_of_x[ORDER];
// loop: local points (mapped in [-1,1])
for (unsigned int m=0; m<M; ++m) {
// evaluate chebyshev polynomials at local points
for (unsigned int o=1; o<ORDER; ++o)
T_of_x[o] = BasisType::T(o, x[m]);
for (unsigned int n=0; n<ORDER; ++n) {
S[n*M + m] = FReal(1.) / ORDER;
for (unsigned int o=1; o<ORDER; ++o)
S[n*M + m] += FReal(2.) / ORDER * T_of_x[o] * T_of_roots[o][n];
}
}
}
/**
* Particle to moment: application of \f$S_\ell(y,\bar y_n)\f$
......@@ -188,6 +267,7 @@ public:
ContainerClass *const localParticles) const;
/*
void applyM2M(const unsigned int ChildIndex,
const FReal *const ChildExpansion,
FReal *const ParentExpansion) const
......@@ -205,6 +285,55 @@ public:
ChildParentInterpolator[ChildIndex],
const_cast<FReal*>(ParentExpansion), ChildExpansion);
}
*/
void applyM2M(const unsigned int ChildIndex,
const FReal *const ChildExpansion,
FReal *const ParentExpansion) const
{
FReal Exp[nnodes], PermExp[nnodes];
FBlas::gemtm(ORDER, ORDER, ORDER*ORDER, FReal(1.),
ChildParentInterpolator[ChildIndex], ORDER,
const_cast<FReal*>(ChildExpansion), ORDER, PermExp, ORDER);
for (unsigned int n=0; n<nnodes; ++n) Exp[n] = PermExp[perm[1][n]];
FBlas::gemtm(ORDER, ORDER, ORDER*ORDER, FReal(1.),
ChildParentInterpolator[ChildIndex] + 2 * ORDER*ORDER, ORDER,
Exp, ORDER, PermExp, ORDER);
for (unsigned int n=0; n<nnodes; ++n) Exp[perm[1][n]] = PermExp[perm[2][n]];
FBlas::gemtm(ORDER, ORDER, ORDER*ORDER, FReal(1.),
ChildParentInterpolator[ChildIndex] + 1 * ORDER*ORDER, ORDER,
Exp, ORDER, PermExp, ORDER);
for (unsigned int n=0; n<nnodes; ++n) ParentExpansion[perm[2][n]] += PermExp[n];
}
void applyL2L(const unsigned int ChildIndex,
const FReal *const ParentExpansion,
FReal *const ChildExpansion) const
{
FReal Exp[nnodes], PermExp[nnodes];
FBlas::gemm(ORDER, ORDER, ORDER*ORDER, FReal(1.),
ChildParentInterpolator[ChildIndex], ORDER,
const_cast<FReal*>(ParentExpansion), ORDER, PermExp, ORDER);
for (unsigned int n=0; n<nnodes; ++n) Exp[n] = PermExp[perm[1][n]];
FBlas::gemm(ORDER, ORDER, ORDER*ORDER, FReal(1.),
ChildParentInterpolator[ChildIndex] + 2 * ORDER*ORDER, ORDER,
Exp, ORDER, PermExp, ORDER);
for (unsigned int n=0; n<nnodes; ++n) Exp[perm[1][n]] = PermExp[perm[2][n]];
FBlas::gemm(ORDER, ORDER, ORDER*ORDER, FReal(1.),
ChildParentInterpolator[ChildIndex] + 1 * ORDER*ORDER, ORDER,
Exp, ORDER, PermExp, ORDER);
for (unsigned int n=0; n<nnodes; ++n) ChildExpansion[perm[2][n]] += PermExp[n];
}
};
......@@ -232,46 +361,45 @@ inline void FChebInterpolator<ORDER>::applyP2M(const FPoint& center,
const map_glob_loc map(center, width);
FPoint localPosition;
FReal T_of_x[ORDER][3];
FReal S[3],c1;
FReal S[3], c1;
//
FReal xpx,ypy,zpz ;
c1 = FReal(8.) / nnodes ;
c1 = FReal(8.) / nnodes ; // 1 flop
// loop over source particles
typename ContainerClass::ConstBasicIterator iter(*sourceParticles);
while(iter.hasNotFinished()){
// map global position to [-1,1]
map(iter.data().getPosition(), localPosition);
map(iter.data().getPosition(), localPosition); // 15 flops
// evaluate chebyshev polynomials of source particle: T_o(x_i)
T_of_x[0][0] = FReal(1.); T_of_x[1][0] = localPosition.getX();
T_of_x[0][1] = FReal(1.); T_of_x[1][1] = localPosition.getY();
T_of_x[0][2] = FReal(1.); T_of_x[1][2] = localPosition.getZ();
xpx = FReal(2.) * localPosition.getX() ;
ypy = FReal(2.) * localPosition.getY() ;
zpz = FReal(2.) * localPosition.getZ() ;
xpx = FReal(2.) * localPosition.getX() ; // 1 flop
ypy = FReal(2.) * localPosition.getY() ; // 1 flop
zpz = FReal(2.) * localPosition.getZ() ; // 1 flop
for (unsigned int o=2; o<ORDER; ++o) {
T_of_x[o][0] = xpx * T_of_x[o-1][0] - T_of_x[o-2][0];
T_of_x[o][1] = ypy * T_of_x[o-1][1] - T_of_x[o-2][1];
T_of_x[o][2] = zpz * T_of_x[o-1][2] - T_of_x[o-2][2];
T_of_x[o][0] = xpx * T_of_x[o-1][0] - T_of_x[o-2][0]; // 2 flops
T_of_x[o][1] = ypy * T_of_x[o-1][1] - T_of_x[o-2][1]; // 2 flops
T_of_x[o][2] = zpz * T_of_x[o-1][2] - T_of_x[o-2][2]; // 2 flops
}
// anterpolate
const FReal sourceValue = iter.data().getPhysicalValue();
for (unsigned int n=0; n<nnodes; ++n) {
const unsigned int j[3] = {node_ids[n][0], node_ids[n][1], node_ids[n][2]};
S[0] = FReal(0.5) + T_of_x[1][0] * T_of_roots[1][j[0]];
S[1] = FReal(0.5) + T_of_x[1][1] * T_of_roots[1][j[1]];
S[2] = FReal(0.5) + T_of_x[1][2] * T_of_roots[1][j[2]];
S[0] = FReal(0.5) + T_of_x[1][0] * T_of_roots[1][j[0]]; // 2 flops
S[1] = FReal(0.5) + T_of_x[1][1] * T_of_roots[1][j[1]]; // 2 flops
S[2] = FReal(0.5) + T_of_x[1][2] * T_of_roots[1][j[2]]; // 2 flops
for (unsigned int o=2; o<ORDER; ++o) {
S[0] += T_of_x[o][0] * T_of_roots[o][j[0]];
S[1] += T_of_x[o][1] * T_of_roots[o][j[1]];
S[2] += T_of_x[o][2] * T_of_roots[o][j[2]];
S[0] += T_of_x[o][0] * T_of_roots[o][j[0]]; // 2 flops
S[1] += T_of_x[o][1] * T_of_roots[o][j[1]]; // 2 flops
S[2] += T_of_x[o][2] * T_of_roots[o][j[2]]; // 2 flops
}
// gather contributions
//
multipoleExpansion[n] += c1 * S[0] * S[1] * S[2] * sourceValue;
multipoleExpansion[n] += c1 * S[0] * S[1] * S[2] * sourceValue; // 4 flops
}
// increment source iterator
iter.gotoNext();
......@@ -489,7 +617,7 @@ inline void FChebInterpolator<ORDER>::applyL2PTotal(const FPoint& center,
// setup local to global mapping
const map_glob_loc map(center, width);
FPoint Jacobian;
map.computeJacobian(Jacobian);
map.computeJacobian(Jacobian); // 6 flops
const FReal jacobian[3] = {Jacobian.getX(), Jacobian.getY(), Jacobian.getZ()};
FPoint localPosition;
FReal T_of_x[ORDER][3];
......@@ -497,19 +625,19 @@ inline void FChebInterpolator<ORDER>::applyL2PTotal(const FPoint& center,
FReal P[6];
//
FReal xpx,ypy,zpz ;
FReal c1 = FReal(8.0) / nnodes ;
FReal c1 = FReal(8.0) / nnodes; // 1 flop
//
typename ContainerClass::BasicIterator iter(*localParticles);
while(iter.hasNotFinished()){
// map global position to [-1,1]
map(iter.data().getPosition(), localPosition);
map(iter.data().getPosition(), localPosition); // 15 flops
// evaluate chebyshev polynomials of source particle
// T_0(x_i) and T_1(x_i)
xpx = FReal(2.) * localPosition.getX() ;
ypy = FReal(2.) * localPosition.getY() ;
zpz = FReal(2.) * localPosition.getZ() ;
xpx = FReal(2.) * localPosition.getX(); // 1 flop
ypy = FReal(2.) * localPosition.getY(); // 1 flop
zpz = FReal(2.) * localPosition.getZ(); // 1 flop
//
T_of_x[0][0] = FReal(1.); T_of_x[1][0] = localPosition.getX();
T_of_x[0][1] = FReal(1.); T_of_x[1][1] = localPosition.getY();
......@@ -530,20 +658,20 @@ inline void FChebInterpolator<ORDER>::applyL2PTotal(const FPoint& center,
// U_of_x[o][0] = FReal(2.)*localPosition.getX()*U_of_x[o-1][0] - U_of_x[o-2][0];
// U_of_x[o][1] = FReal(2.)*localPosition.getY()*U_of_x[o-1][1] - U_of_x[o-2][1];
// U_of_x[o][2] = FReal(2.)*localPosition.getZ()*U_of_x[o-1][2] - U_of_x[o-2][2];
T_of_x[o][0] = xpx*T_of_x[o-1][0] - T_of_x[o-2][0];
T_of_x[o][1] = ypy*T_of_x[o-1][1] - T_of_x[o-2][1];
T_of_x[o][2] = zpz*T_of_x[o-1][2] - T_of_x[o-2][2];
T_of_x[o][0] = xpx * T_of_x[o-1][0] - T_of_x[o-2][0]; // 2 flops
T_of_x[o][1] = ypy * T_of_x[o-1][1] - T_of_x[o-2][1]; // 2 flops
T_of_x[o][2] = zpz * T_of_x[o-1][2] - T_of_x[o-2][2]; // 2 flops
// U_o(x_i)
U_of_x[o][0] = xpx*U_of_x[o-1][0] - U_of_x[o-2][0];
U_of_x[o][1] = ypy*U_of_x[o-1][1] - U_of_x[o-2][1];
U_of_x[o][2] = zpz*U_of_x[o-1][2] - U_of_x[o-2][2];
U_of_x[o][0] = xpx * U_of_x[o-1][0] - U_of_x[o-2][0]; // 2 flops
U_of_x[o][1] = ypy * U_of_x[o-1][1] - U_of_x[o-2][1]; // 2 flops
U_of_x[o][2] = zpz * U_of_x[o-1][2] - U_of_x[o-2][2]; // 2 flops
}
// scale, because dT_o/dx = oU_{o-1}
for (unsigned int o=2; o<ORDER; ++o) {
U_of_x[o-1][0] *= FReal(o);
U_of_x[o-1][1] *= FReal(o);
U_of_x[o-1][2] *= FReal(o);
U_of_x[o-1][0] *= FReal(o); // 1 flops
U_of_x[o-1][1] *= FReal(o); // 1 flops
U_of_x[o-1][2] *= FReal(o); // 1 flops
}
// apply P and increment forces
......@@ -557,38 +685,38 @@ inline void FChebInterpolator<ORDER>::applyL2PTotal(const FPoint& center,
// tensor indices of chebyshev nodes
const unsigned int j[3] = {node_ids[n][0], node_ids[n][1], node_ids[n][2]};
//
P[0] = FReal(0.5) + T_of_x[1][0] * T_of_roots[1][j[0]];
P[1] = FReal(0.5) + T_of_x[1][1] * T_of_roots[1][j[1]];
P[2] = FReal(0.5) + T_of_x[1][2] * T_of_roots[1][j[2]];
P[3] = U_of_x[0][0] * T_of_roots[1][j[0]];
P[4] = U_of_x[0][1] * T_of_roots[1][j[1]];
P[5] = U_of_x[0][2] * T_of_roots[1][j[2]];
P[0] = FReal(0.5) + T_of_x[1][0] * T_of_roots[1][j[0]]; // 2 flops
P[1] = FReal(0.5) + T_of_x[1][1] * T_of_roots[1][j[1]]; // 2 flops
P[2] = FReal(0.5) + T_of_x[1][2] * T_of_roots[1][j[2]]; // 2 flops
P[3] = U_of_x[0][0] * T_of_roots[1][j[0]]; // 1 flop
P[4] = U_of_x[0][1] * T_of_roots[1][j[1]]; // 1 flop
P[5] = U_of_x[0][2] * T_of_roots[1][j[2]]; // 1 flop
for (unsigned int o=2; o<ORDER; ++o) {
P[0] += T_of_x[o ][0] * T_of_roots[o][j[0]];
P[1] += T_of_x[o ][1] * T_of_roots[o][j[1]];
P[2] += T_of_x[o ][2] * T_of_roots[o][j[2]];
P[3] += U_of_x[o-1][0] * T_of_roots[o][j[0]];
P[4] += U_of_x[o-1][1] * T_of_roots[o][j[1]];
P[5] += U_of_x[o-1][2] * T_of_roots[o][j[2]];
P[0] += T_of_x[o ][0] * T_of_roots[o][j[0]]; // 2 flop
P[1] += T_of_x[o ][1] * T_of_roots[o][j[1]]; // 2 flop
P[2] += T_of_x[o ][2] * T_of_roots[o][j[2]]; // 2 flop
P[3] += U_of_x[o-1][0] * T_of_roots[o][j[0]]; // 2 flop
P[4] += U_of_x[o-1][1] * T_of_roots[o][j[1]]; // 2 flop
P[5] += U_of_x[o-1][2] * T_of_roots[o][j[2]]; // 2 flop
}
//
potential += P[0] * P[1] * P[2] * localExpansion[n];
forces[0] += P[3] * P[1] * P[2] * localExpansion[n];
forces[1] += P[0] * P[4] * P[2] * localExpansion[n];
forces[2] += P[0] * P[1] * P[5] * localExpansion[n];
potential += P[0] * P[1] * P[2] * localExpansion[n]; // 4 flops
forces[0] += P[3] * P[1] * P[2] * localExpansion[n]; // 4 flops
forces[1] += P[0] * P[4] * P[2] * localExpansion[n]; // 4 flops
forces[2] += P[0] * P[1] * P[5] * localExpansion[n]; // 4 flops
}
//
potential *= c1 ;
forces[0] *= jacobian[0] *c1;
forces[1] *= jacobian[1] *c1;
forces[2] *= jacobian[2] *c1;
potential *= c1 ; // 1 flop
forces[0] *= jacobian[0] *c1; // 2 flops
forces[1] *= jacobian[1] *c1; // 2 flops
forces[2] *= jacobian[2] *c1; // 2 flops
// set computed potential
iter.data().incPotential(potential);
iter.data().incPotential(potential); // 1 flop
// set computed forces
iter.data().incForces(forces[0] * iter.data().getPhysicalValue(),
forces[1] * iter.data().getPhysicalValue(),
forces[2] * iter.data().getPhysicalValue());
forces[1] * iter.data().getPhysicalValue(),
forces[2] * iter.data().getPhysicalValue()); // 6 flops
// increment iterator
iter.gotoNext();
......
......@@ -94,20 +94,17 @@ public:
*/
void operator()(const FPoint& globPos, FPoint& loclPos) const
{
loclPos.setX((FReal(2.)*globPos.getX()-b.getX()-a.getX()) /
(b.getX()-a.getX()));
loclPos.setY((FReal(2.)*globPos.getY()-b.getY()-a.getY()) /
(b.getY()-a.getY()));
loclPos.setZ((FReal(2.)*globPos.getZ()-b.getZ()-a.getZ()) /
(b.getZ()-a.getZ()));
loclPos.setX((FReal(2.)*globPos.getX()-b.getX()-a.getX()) / (b.getX()-a.getX())); // 5 flops
loclPos.setY((FReal(2.)*globPos.getY()-b.getY()-a.getY()) / (b.getY()-a.getY())); // 5 flops
loclPos.setZ((FReal(2.)*globPos.getZ()-b.getZ()-a.getZ()) / (b.getZ()-a.getZ())); // 5 flops
}
// jacobian = 2 / (b - a);
void computeJacobian(FPoint& jacobian) const
{
jacobian.setX(FReal(2.) / (b.getX() - a.getX()));
jacobian.setY(FReal(2.) / (b.getY() - a.getY()));
jacobian.setZ(FReal(2.) / (b.getZ() - a.getZ()));
jacobian.setX(FReal(2.) / (b.getX() - a.getX())); // 2 flops
jacobian.setY(FReal(2.) / (b.getY() - a.getY())); // 2 flops
jacobian.setZ(FReal(2.) / (b.getZ() - a.getZ())); // 2 flops
}
};
......
......@@ -81,119 +81,119 @@ struct FChebRoots : FNoCopyable
// order 2
template<> const double FChebRoots<2>::roots[] = {0.707106781186548,
-0.707106781186547};
template<> const double FChebRoots<2>::roots[] = {-0.707106781186548,
0.707106781186547};
// order 3
template<> const double FChebRoots<3>::roots[] = {8.66025403784439e-01,
0.0,
-8.66025403784438e-01};
template<> const double FChebRoots<3>::roots[] = {-8.66025403784439e-01,
0.0,
8.66025403784438e-01};
// order 4
template<> const double FChebRoots<4>::roots[] = {0.923879532511287,
0.382683432365090,
-0.382683432365090,
-0.923879532511287};
template<> const double FChebRoots<4>::roots[] = {-0.923879532511287,
-0.382683432365090,
0.382683432365090,
0.923879532511287};
// order 5
template<> const double FChebRoots<5>::roots[] = {9.51056516295154e-01,
5.87785252292473e-01,
0.0,
-5.87785252292473e-01,
-9.51056516295154e-01};
template<> const double FChebRoots<5>::roots[] = {-9.51056516295154e-01,
-5.87785252292473e-01,
0.0,
5.87785252292473e-01,
9.51056516295154e-01};
// order 6
template<> const double FChebRoots<6>::roots[] = {0.965925826289068,
0.707106781186548,
0.258819045102521,
-0.258819045102521,
-0.707106781186547,
-0.965925826289068};
template<> const double FChebRoots<6>::roots[] = {-0.965925826289068,
-0.707106781186548,
-0.258819045102521,
0.258819045102521,
0.707106781186547,
0.965925826289068};
// order 7
template<> const double FChebRoots<7>::roots[] = {9.74927912181824e-01,
7.81831482468030e-01,
4.33883739117558e-01,
0.0,
-4.33883739117558e-01,
-7.81831482468030e-01,
-9.74927912181824e-01};
template<> const double FChebRoots<7>::roots[] = {-9.74927912181824e-01,
-7.81831482468030e-01,
-4.33883739117558e-01,
0.0,
4.33883739117558e-01,
7.81831482468030e-01,
9.74927912181824e-01};
// order 8
template<> const double FChebRoots<8>::roots[] = {0.980785280403230,
0.831469612302545,
0.555570233019602,
0.195090322016128,
-0.195090322016128,
-0.555570233019602,
-0.831469612302545,
-0.980785280403230};
template<> const double FChebRoots<8>::roots[] = {-0.980785280403230,
-0.831469612302545,
-0.555570233019602,
-0.195090322016128,
0.195090322016128,
0.555570233019602,
0.831469612302545,
0.980785280403230};
// order 9
template<> const double FChebRoots<9>::roots[] = {9.84807753012208e-01,
8.66025403784439e-01,
6.42787609686539e-01,
3.42020143325669e-01,
0.0,
-3.42020143325669e-01,
-6.42787609686539e-01,
-8.66025403784438e-01,
-9.84807753012208e-01,};
template<> const double FChebRoots<9>::roots[] = {-9.84807753012208e-01,
-8.66025403784439e-01,
-6.42787609686539e-01,
-3.42020143325669e-01,
0.0,
3.42020143325669e-01,
6.42787609686539e-01,
8.66025403784438e-01,
9.84807753012208e-01,};
// order 10
template<> const double FChebRoots<10>::roots[] = {0.987688340595138,
0.891006524188368,
0.707106781186548,
0.453990499739547,
0.156434465040231,
-0.156434465040231,
-0.453990499739547,
-0.707106781186547,
-0.891006524188368,
-0.987688340595138};
template<> const double FChebRoots<10>::roots[] = {-0.987688340595138,
-0.891006524188368,
-0.707106781186548,
-0.453990499739547,
-0.156434465040231,
0.156434465040231,
0.453990499739547,
0.707106781186547,
0.891006524188368,
0.987688340595138};
// order 11
template<> const double FChebRoots<11>::roots[] = {9.89821441880933e-01,
9.09631995354518e-01,
7.55749574354258e-01,
5.40640817455598e-01,
2.81732556841430e-01,
0.0,
-2.81732556841430e-01,
-5.40640817455597e-01,
-7.55749574354258e-01,
-9.09631995354518e-01,
-9.89821441880933e-01};
template<> const double FChebRoots<11>::roots[] = {-9.89821441880933e-01,
-9.09631995354518e-01,
-7.55749574354258e-01,
-5.40640817455598e-01,
-2.81732556841430e-01,
0.0,
2.81732556841430e-01,
5.40640817455597e-01,
7.55749574354258e-01,
9.09631995354518e-01,
9.89821441880933e-01};
// order 12
template<> const double FChebRoots<12>::roots[] = {0.991444861373810,
0.923879532511287,
0.793353340291235,
0.608761429008721,
0.382683432365090,
0.130526192220052,
-0.130526192220052,
-0.382683432365090,
-0.608761429008721,
-0.793353340291235,
-0.923879532511287,
-0.991444861373810};
template<> const double FChebRoots<12>::roots[] = {-0.991444861373810,
-0.923879532511287,
-0.793353340291235,
-0.608761429008721,
-0.382683432365090,
-0.130526192220052,
0.130526192220052,
0.382683432365090,
0.608761429008721,
0.793353340291235,
0.923879532511287,
0.991444861373810};
// order 13
template<> const double FChebRoots<13>::roots[] = {9.92708874098054e-01,
9.35016242685415e-01,