Merge branch 'master' of git+ssh://scm.gforge.inria.fr//gitroot/scalfmm/scalfmm

CMakeLists.txt

@@ -42,6 +42,7 @@ MESSAGE(STATUS " CXX  ${CMAKE_CXX_COMPILER_ID}" )
 #
 # Options
 OPTION( ScalFMM_USE_BLAS 		"Set to ON to build ScaFMM with BLAS" OFF )
+OPTION( ScalFMM_USE_FFTW 		"Set to ON to build ScaFMM with FFTW" OFF )
 OPTION( ScalFMM_USE_TRACE 		"Set to ON to print trace or use itac trace" OFF )
 OPTION( ScalFMM_BUILD_TESTS 		"Set to ON to build fonctionnalities Tests" OFF )
 OPTION( ScalFMM_BUILD_UTESTS 		"Set to ON to build UTests" OFF )
@@ -138,7 +139,8 @@ endif()
 if( ScalFMM_USE_BLAS )
   OPTION( ScalFMM_USE_MKL_AS_BLAS "Set to ON to use MKL CBLAS" OFF )
   if( ScalFMM_USE_MKL_AS_BLAS )
-     SET(BLAS_LIBRARIES "-L$ENV{MKLROOT}/lib;-lmkl_intel_lp64;-lmkl_sequential;-lmkl_core" CACHE STRING "Set your MKL flags")
+     SET(BLAS_LIBRARIES
+  "-L$ENV{MKLROOT}/lib;-lmkl_intel_lp64;-lmkl_sequential;-lmkl_core" CACHE STRING "Set your MKL flags")
      SET(LAPACK_LIBRARIES "")
   elseif(ScalFMM_USE_EXTERNAL_BLAS) 
          MESSAGE(STATUS "BLAS SET BY EXTERNAL PROGRAM = ${BLAS_LIBRARIES}")
@@ -150,6 +152,18 @@ if( ScalFMM_USE_BLAS )
   MESSAGE(STATUS "SCALFMM_LIBRARIES          = ${SCALFMM_LIBRARIES}")
 endif(ScalFMM_USE_BLAS)
 
+# FFTW option
+if( ScalFMM_USE_FFTW )
+  OPTION( ScalFMM_USE_MKL_AS_FFTW "Set to ON to use MKL FFTW" ON )
+  if( ScalFMM_USE_MKL_AS_FFTW )
+    SET(FFTW_LIBRARIES "-I$ENV{MKLROOT}/include/fftw; -L$ENV{MKLROOT}/lib; -lmkl_intel_lp64; -lmkl_sequential; -lmkl_core; -lpthread; -lm" CACHE STRING "Set your MKL flags")
+  else()
+    SET(FFTW_LIBRARIES "-lfftw3" CACHE STRING "Use LIBFFTW")
+  endif()
+  SET(SCALFMM_LIBRARIES "${SCALFMM_LIBRARIES}; ${FFTW_LIBRARIES}")
+  MESSAGE(STATUS "SCALFMM_LIBRARIES          = ${SCALFMM_LIBRARIES}")
+endif(ScalFMM_USE_FFTW)
+
 # Compile option
 #ADD_DEFINITIONS(-Wall -Wshadow -Wpointer-arith -Wcast-qual -Wconversion -fpic )
 # 

Src/Kernels/Chebyshev/FAbstractChebKernel.hpp

@@ -22,7 +22,7 @@
 
 #include "../../Components/FAbstractKernels.hpp"
 
-#include "FChebP2PKernels.hpp"
+#include "../Interpolation/FInterpP2PKernels.hpp"
 #include "./FChebInterpolator.hpp"
 
 #include "../../Containers/FTreeCoordinate.hpp"

Src/Kernels/Chebyshev/FChebInterpolator.hpp

@@ -17,7 +17,7 @@
 #define FCHEBINTERPOLATOR_HPP
 
 
-#include "./FChebMapping.hpp"
+#include "./../Interpolation/FInterpMapping.hpp"
 #include "./FChebTensor.hpp"
 #include "./FChebRoots.hpp"
 
@@ -296,7 +296,7 @@ class FChebInterpolator : FNoCopyable
 
             // set child info
             FChebTensor<ORDER>::setRelativeChildCenter(child, ChildCenter);
-            FChebTensor<ORDER>::setChebyshevRoots(ChildCenter, ChildWidth, ChildCoords);
+            FChebTensor<ORDER>::setPolynomialsRoots(ChildCenter, ChildWidth, ChildCoords);
 
             // allocate memory
             ChildParentInterpolator[child] = new FReal [3 * ORDER*ORDER];

Src/Kernels/Chebyshev/FChebKernel.hpp

@@ -70,8 +70,8 @@ public:
 			M2LHandler(new M2LHandlerClass(Epsilon))
 	{
 		// read precomputed compressed m2l operators from binary file
-		M2LHandler->ReadFromBinaryFileAndSet();
-		//M2LHandler->ComputeAndCompressAndSet();
+		//M2LHandler->ReadFromBinaryFileAndSet();
+		M2LHandler->ComputeAndCompressAndSet();
 	}
 
 
@@ -209,14 +209,14 @@ public:
                      ContainerClass* const NeighborSourceParticles[27],
                      const int /* size */)
     {
-        DirectInteactionComputer<MatrixKernelClass::Identifier, NVALS>::P2P(TargetParticles,NeighborSourceParticles);
+        DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2P(TargetParticles,NeighborSourceParticles);
     }
 
 
     void P2PRemote(const FTreeCoordinate& /*inPosition*/,
                    ContainerClass* const FRestrict inTargets, const ContainerClass* const FRestrict /*inSources*/,
                    ContainerClass* const inNeighbors[27], const int /*inSize*/){
-        DirectInteactionComputer<MatrixKernelClass::Identifier, NVALS>::P2PRemote(inTargets,inNeighbors,27);
+        DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2PRemote(inTargets,inNeighbors,27);
     }
 
 };

Src/Kernels/Chebyshev/FChebM2LHandler.hpp

@@ -99,8 +99,12 @@ public:
 		// check if aready set
 		if (U||C||B) throw std::runtime_error("Compressed M2L operator already set");
 		rank = ComputeAndCompress(epsilon, U, C, B);
+
+    unsigned long sizeM2L = 343*rank*rank*sizeof(FReal);
+
+
 		// write info
-		std::cout << "Compressed and set M2L operators (" << rank << ") in "
+		std::cout << "Compressed and set M2L operators (" << long(sizeM2L) << " B) in "
 							<< time.tacAndElapsed() << "sec."	<< std::endl;
 	}
 

Src/Kernels/Chebyshev/FChebMatrixKernel.hpp

-// ===================================================================================
-// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Bérenger Bramas, Matthias Messner
-// olivier.coulaud@inria.fr, berenger.bramas@inria.fr
-// This software is a computer program whose purpose is to compute the FMM.
-//
-// This software is governed by the CeCILL-C and LGPL licenses and
-// abiding by the rules of distribution of free software.  
-// 
-// This program is distributed in the hope that it will be useful,
-// but WITHOUT ANY WARRANTY; without even the implied warranty of
-// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
-// GNU General Public and CeCILL-C Licenses for more details.
-// "http://www.cecill.info". 
-// "http://www.gnu.org/licenses".
-// ===================================================================================
-#ifndef FCHEBMATRIXKERNEL_HPP
-#define FCHEBMATRIXKERNEL_HPP
-
-#include "../../Utils/FPoint.hpp"
-#include "../../Utils/FNoCopyable.hpp"
-#include "../../Utils/FMath.hpp"
-#include "../../Utils/FBlas.hpp"
-
-
-// extendable
-enum KERNEL_FUNCCTION_IDENTIFIER {ONE_OVER_R,
-																	ONE_OVER_R_SQUARED,
-																	LEONARD_JONES_POTENTIAL};
-
-// probably not extedable :)
-enum KERNEL_FUNCTION_TYPE {HOMOGENEOUS, NON_HOMOGENEOUS};
-
-
-/**
- * @author Matthias Messner (matthias.messner@inria.fr)
- * @class FChebMatrixKernels
- * Please read the license
- */
-struct FChebAbstractMatrixKernel : FNoCopyable
-{
-    virtual ~FChebAbstractMatrixKernel(){} // to remove warning
-	virtual FReal evaluate(const FPoint&, const FPoint&) const = 0;
-	// I need both functions because required arguments are not always given
-	virtual FReal getScaleFactor(const FReal, const int) const = 0;
-	virtual FReal getScaleFactor(const FReal) const = 0;
-};
-
-
-
-/// One over r
-struct FChebMatrixKernelR : FChebAbstractMatrixKernel
-{
-	static const KERNEL_FUNCTION_TYPE Type = HOMOGENEOUS;
-	static const KERNEL_FUNCCTION_IDENTIFIER Identifier = ONE_OVER_R;
-
-	FChebMatrixKernelR() {}
-
-	FReal evaluate(const FPoint& x, const FPoint& y) const
-	{
-		const FPoint xy(x-y);
-		return FReal(1.) / FMath::Sqrt(xy.getX()*xy.getX() +
-																	 xy.getY()*xy.getY() +
-																	 xy.getZ()*xy.getZ());
-	}
-
-	FReal getScaleFactor(const FReal RootCellWidth, const int TreeLevel) const
-	{
-		const FReal CellWidth(RootCellWidth / FReal(FMath::pow(2, TreeLevel)));
-		return getScaleFactor(CellWidth);
-	}
-
-	FReal getScaleFactor(const FReal CellWidth) const
-	{
-		return FReal(2.) / CellWidth;
-	}
-};
-
-
-
-/// One over r^2
-struct FChebMatrixKernelRR : FChebAbstractMatrixKernel
-{
-	static const KERNEL_FUNCTION_TYPE Type = HOMOGENEOUS;
-	static const KERNEL_FUNCCTION_IDENTIFIER Identifier = ONE_OVER_R_SQUARED;
-
-	FChebMatrixKernelRR() {}
-
-	FReal evaluate(const FPoint& x, const FPoint& y) const
-	{
-		const FPoint xy(x-y);
-		return FReal(1.) / FReal(xy.getX()*xy.getX() +
-														 xy.getY()*xy.getY() +
-														 xy.getZ()*xy.getZ());
-	}
-
-	FReal getScaleFactor(const FReal RootCellWidth, const int TreeLevel) const
-	{
-		const FReal CellWidth(RootCellWidth / FReal(FMath::pow(2, TreeLevel)));
-		return getScaleFactor(CellWidth);
-	}
-
-	FReal getScaleFactor(const FReal CellWidth) const
-	{
-		return FReal(4.) / CellWidth;
-	}
-};
-
-
-
-/// One over r^12 - One over r^6
-struct FChebMatrixKernelLJ : FChebAbstractMatrixKernel
-{
-	static const KERNEL_FUNCTION_TYPE Type = NON_HOMOGENEOUS;
-	static const KERNEL_FUNCCTION_IDENTIFIER Identifier = LEONARD_JONES_POTENTIAL;
-
-	FChebMatrixKernelLJ() {}
-
-	FReal evaluate(const FPoint& x, const FPoint& y) const
-	{
-		const FPoint xy(x-y);
-		const FReal r = xy.norm();
-		const FReal r3 = r*r*r;
-		const FReal one_over_r6 = FReal(1.) / (r3*r3);
-		//return one_over_r6 * one_over_r6;
-		//return one_over_r6;
-		return one_over_r6 * one_over_r6 - one_over_r6;
-	}
-
-	FReal getScaleFactor(const FReal, const int) const
-	{
-		// return 1 because non homogeneous kernel functions cannot be scaled!!!
-		return FReal(1.);
-	}
-
-	FReal getScaleFactor(const FReal) const
-	{
-		// return 1 because non homogeneous kernel functions cannot be scaled!!!
-		return FReal(1.);
-	}
-};
-
-
-
-
-
-/*!  Functor which provides the interface to assemble a matrix based on the
-  number of rows and cols and on the coordinates x and y and the type of the
-  generating matrix-kernel function.
- */
-template <typename MatrixKernel>
-class EntryComputer
-{
-	const MatrixKernel Kernel;
-
-	const unsigned int nx, ny;
-	const FPoint *const px, *const py;
-
-	const FReal *const weights;
-
-public:
-	explicit EntryComputer(const unsigned int _nx, const FPoint *const _px,
-												 const unsigned int _ny, const FPoint *const _py,
-												 const FReal *const _weights = NULL)
-		: Kernel(),	nx(_nx), ny(_ny), px(_px), py(_py), weights(_weights) {}
-
-//	template <typename Point>
-//	void operator()(const unsigned int nx, const Point *const px,
-//									const unsigned int ny, const Point *const py,
-//									FReal *const data) const
-//	{
-//		for (unsigned int j=0; j<ny; ++j)
-//			for (unsigned int i=0; i<nx; ++i)
-//				data[j*nx + i] = Kernel.evaluate(px[i], py[j]);
-//	}
-
-	void operator()(const unsigned int xbeg, const unsigned int xend,
-									const unsigned int ybeg, const unsigned int yend,
-									FReal *const data) const
-	{
-		unsigned int idx = 0;
-		if (weights) {
-			for (unsigned int j=ybeg; j<yend; ++j)
-				for (unsigned int i=xbeg; i<xend; ++i)
-					data[idx++] = weights[i] * weights[j] * Kernel.evaluate(px[i], py[j]);
-		} else {
-			for (unsigned int j=ybeg; j<yend; ++j)
-				for (unsigned int i=xbeg; i<xend; ++i)
-					data[idx++] = Kernel.evaluate(px[i], py[j]);
-		}
-
-		/*
-		// apply weighting matrices
-		if (weights) {
-			if ((xend-xbeg) == (yend-ybeg) && (xend-xbeg) == nx)
-				for (unsigned int n=0; n<nx; ++n) {
-					FBlas::scal(nx, weights[n], data + n,  nx); // scale rows
-					FBlas::scal(nx, weights[n], data + n * nx); // scale cols
-				}
-			else if ((xend-xbeg) == 1 && (yend-ybeg) == ny)
-				for (unsigned int j=0; j<ny; ++j)	data[j] *= weights[j];
-			else if ((yend-ybeg) == 1 && (xend-xbeg) == nx)
-				for (unsigned int i=0; i<nx; ++i)	data[i] *= weights[i];
-		}
-		*/
-
-	}
-};
-
-
-
-
-
-#endif // FCHEBMATRIXKERNEL_HPP
-
-// [--END--]

Src/Kernels/Chebyshev/FChebP2PKernels.hpp

-#ifndef FCHEBP2PKERNELS_HPP
-#define FCHEBP2PKERNELS_HPP
-
-
-#include "../P2P/FP2P.hpp"
-
-template <KERNEL_FUNCCTION_IDENTIFIER Identifier, int NVALS>
-struct DirectInteactionComputer;
-
-///////////////////////////////////////////////////////
-// P2P Wrappers
-///////////////////////////////////////////////////////
-
-/*! Specialization for Laplace potential */
-template <>
-struct DirectInteactionComputer<ONE_OVER_R, 1>
-{
-    template <typename ContainerClass>
-    static void P2P(		 ContainerClass* const FRestrict TargetParticles,
-                     ContainerClass* const NeighborSourceParticles[27]){
-        FP2P::FullMutual(TargetParticles,NeighborSourceParticles,14);
-    }
-
-    template <typename ContainerClass>
-    static void P2PRemote( ContainerClass* const FRestrict inTargets,
-                           ContainerClass* const inNeighbors[27],
-                           const int inSize){
-        FP2P::FullRemote(inTargets,inNeighbors,inSize);
-    }
-};
-
-
-/*! Specialization for Leonard-Jones potential */
-template <>
-struct DirectInteactionComputer<LEONARD_JONES_POTENTIAL, 1>
-{
-    template <typename ContainerClass>
-    static void P2P(		 ContainerClass* const FRestrict TargetParticles,
-                     ContainerClass* const NeighborSourceParticles[27]){
-        FP2P::FullMutualLJ(TargetParticles,NeighborSourceParticles,14);
-    }
-
-    template <typename ContainerClass>
-    static void P2PRemote( ContainerClass* const FRestrict inTargets,
-                           ContainerClass* const inNeighbors[27],
-                           const int inSize){
-        FP2P::FullRemoteLJ(inTargets,inNeighbors,inSize);
-    }
-};
-
-///////////////////////////////////////////////////////
-// In case of multi right hand side
-///////////////////////////////////////////////////////
-
-
-
-template <int NVALS>
-struct DirectInteactionComputer<ONE_OVER_R, NVALS>
-{
-    template <typename ContainerClass>
-    static void P2P(		 ContainerClass* const FRestrict TargetParticles,
-                     ContainerClass* const NeighborSourceParticles[27]){
-        for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
-            FP2P::FullMutual(TargetParticles,NeighborSourceParticles,14);
-        }
-    }
-
-    template <typename ContainerClass>
-    static void P2PRemote( ContainerClass* const FRestrict inTargets,
-                           ContainerClass* const inNeighbors[27],
-                           const int inSize){
-        for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
-            FP2P::FullRemote(inTargets,inNeighbors,inSize);
-        }
-    }
-};
-
-
-/*! Specialization for Leonard-Jones potential */
-template <int NVALS>
-struct DirectInteactionComputer<LEONARD_JONES_POTENTIAL, NVALS>
-{
-    template <typename ContainerClass>
-    static void P2P(		 ContainerClass* const FRestrict TargetParticles,
-                     ContainerClass* const NeighborSourceParticles[27]){
-        for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
-            FP2P::FullMutualLJ(TargetParticles,NeighborSourceParticles,14);
-        }
-    }
-
-    template <typename ContainerClass>
-    static void P2PRemote( ContainerClass* const FRestrict inTargets,
-                           ContainerClass* const inNeighbors[27],
-                           const int inSize){
-        for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
-            FP2P::FullRemoteLJ(inTargets,inNeighbors,inSize);
-        }
-    }
-};
-
-#endif // FCHEBP2PKERNELS_HPP

Src/Kernels/Chebyshev/FChebSymKernel.hpp

@@ -454,14 +454,14 @@ public:
                      ContainerClass* const NeighborSourceParticles[27],
                      const int /* size */)
     {
-        DirectInteactionComputer<MatrixKernelClass::Identifier, NVALS>::P2P(TargetParticles,NeighborSourceParticles);
+        DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2P(TargetParticles,NeighborSourceParticles);
     }
 
 
     void P2PRemote(const FTreeCoordinate& /*inPosition*/,
                    ContainerClass* const FRestrict inTargets, const ContainerClass* const FRestrict /*inSources*/,
                    ContainerClass* const inNeighbors[27], const int /*inSize*/){
-        DirectInteactionComputer<MatrixKernelClass::Identifier, NVALS>::P2PRemote(inTargets,inNeighbors,27);
+        DirectInteractionComputer<MatrixKernelClass::Identifier, NVALS>::P2PRemote(inTargets,inNeighbors,27);
     }
 
 };

Src/Kernels/Chebyshev/FChebSymM2LHandler.hpp

@@ -479,9 +479,12 @@ static void precompute(const MatrixKernelClass *const MatrixKernel, const FReal
 			}
 		}
 	}
-	//	std::cout << "The approximation of the " << counter
-	//					<< " far-field interactions (overall rank " << overall_rank << ") took "
-	//					<< overall_time << "s\n" << std::endl;
+		std::cout << "The approximation of the " << counter
+              << " far-field interactions (overall rank " << overall_rank 
+              << " / " << 16*nnodes 
+              << " , sizeM2L= " << 2*overall_rank*nnodes*sizeof(FReal) << ""
+              << " / " << 16*nnodes*nnodes*sizeof(FReal) << " B"
+              << ") took " << overall_time << "s\n" << std::endl;
 	delete [] U;
 	delete [] WORK;
 	delete [] VT;

Src/Kernels/Chebyshev/FChebTensor.hpp

@@ -4,13 +4,13 @@
 // This software is a computer program whose purpose is to compute the FMM.
 //
 // This software is governed by the CeCILL-C and LGPL licenses and
-// abiding by the rules of distribution of free software.  
-// 
+// abiding by the rules of distribution of free software.
+//
 // This program is distributed in the hope that it will be useful,
 // but WITHOUT ANY WARRANTY; without even the implied warranty of
 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 // GNU General Public and CeCILL-C Licenses for more details.
-// "http://www.cecill.info". 
+// "http://www.cecill.info".
 // "http://www.gnu.org/licenses".
 // ===================================================================================
 #ifndef FCHEBTENSOR_HPP
@@ -19,7 +19,7 @@
 #include "../../Utils/FMath.hpp"
 
 #include "./FChebRoots.hpp"
-#include "./FChebMapping.hpp"
+#include "./../Interpolation/FInterpTensor.hpp"
 
 
 /**
@@ -27,19 +27,6 @@
  * Please read the license
  */
 
-/**
- * @class TensorTraits
- *
- * The class @p TensorTraits gives the number of interpolation nodes per
- * cluster in 3D, depending on the interpolation order.
- *
- * @tparam ORDER interpolation order
- */
-template <int ORDER> struct TensorTraits
-{
-	enum {nnodes = ORDER*ORDER*ORDER};
-};
-
 
 
 /**
@@ -51,123 +38,37 @@ template <int ORDER> struct TensorTraits
  * @tparam ORDER interpolation order \f$\ell\f$
  */
 template <int ORDER>
-class FChebTensor : FNoCopyable
+class FChebTensor : public FInterpTensor<ORDER,FChebRoots<ORDER>>
 {
-	enum {nnodes = TensorTraits<ORDER>::nnodes};
-  typedef FChebRoots<ORDER> BasisType;
-	
-public:
-	/**
-	 * Sets the ids of the coordinates of all \f$\ell^3\f$ interpolation
-	 * nodes
-	 *
-	 * @param[out] NodeIds ids of coordinates of interpolation nodes
-	 */
-	static
-	void setNodeIds(unsigned int NodeIds[nnodes][3])
-	{
-		for (unsigned int n=0; n<nnodes; ++n) {
-			NodeIds[n][0] =  n         % ORDER;
-			NodeIds[n][1] = (n/ ORDER) % ORDER;
-			NodeIds[n][2] =  n/(ORDER  * ORDER);
-		}
-	}
-
-
-	/**
-	 * Sets the roots of the Chebyshev quadrature weights defined as \f$w_i =
-	 * \frac{\pi}{\ell}\sqrt{1-\bar x_i^2}\f$ with the Chebyshev roots \f$\bar
-	 * x\f$.
-	 *
-	 * @param weights[out] the root of the weights \f$\sqrt{w_i}\f$
-	 */
-	static
-	void setRootOfWeights(FReal weights[nnodes])
-	{
-		// weights in 1d
-		FReal weights_1d[ORDER];
-		for (unsigned int o=0; o<ORDER; ++o)
-			weights_1d[o] = FMath::FPi/ORDER * FMath::Sqrt(FReal(1.)-FReal(BasisType::roots[o])*FReal(BasisType::roots[o]));
-		// weights in 3d (tensor structure)
-		unsigned int node_ids[nnodes][3];
-		setNodeIds(node_ids);
-		for (unsigned int n=0; n<nnodes; ++n) {
-			weights[n] = FMath::Sqrt(weights_1d[node_ids[n][0]]*weights_1d[node_ids[n][1]]*weights_1d[node_ids[n][2]]);
-		}
-	}
-
-
-	/**
-	 * Sets the interpolation points in the cluster with @p center and @p width
-	 *
-	 * @param[in] center of cluster
-	 * @param[in] width of cluster
-	 * @param[out] rootPositions coordinates of interpolation points
-	 */
-	static
-	void setRoots(const FPoint& center, const FReal width, FPoint rootPositions[nnodes])
-	{
-		unsigned int node_ids[nnodes][3];
-		setNodeIds(node_ids);
-		const map_loc_glob map(center, width);
-		FPoint localPosition;
+  enum {nnodes = TensorTraits<ORDER>::nnodes};
+  typedef FChebRoots<ORDER> BasisType; 
+  typedef FInterpTensor<ORDER,BasisType> ParentTensor;
+
+ public:
+
+  /**
+   * Sets the roots of the Chebyshev quadrature weights defined as \f$w_i =
+   * \frac{\pi}{\ell}\sqrt{1-\bar x_i^2}\f$ with the Chebyshev roots \f$\bar
+   * x\f$.
+   *
+   * @param weights[out] the root of the weights \f$\sqrt{w_i}\f$
+   */
+  static
+    void setRootOfWeights(FReal weights[nnodes])
+  {
+    // weights in 1d
+    FReal weights_1d[ORDER];
+    for (unsigned int o=0; o<ORDER; ++o)
+      weights_1d[o] = FMath::FPi/ORDER * FMath::Sqrt(FReal(1.)-FReal(BasisType::roots[o])*FReal(BasisType::roots[o]));
+    // weights in 3d (tensor structure)
+    unsigned int node_ids[nnodes][3];
+    ParentTensor::setNodeIds(node_ids);