Commit b683706a authored by BRAMAS Berenger's avatar BRAMAS Berenger
Browse files

Remove the toRemove directory

parent 814cf454
// ===================================================================================
// Copyright ScalFmm 2011 INRIA, Olivier Coulaud, Berenger Bramas
// 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 FADAPTCELL_HPP
#define FADAPTCELL_HPP
#include <cstddef>
#include <iostream>
#include <vector>
//
#include "Components/FBasicCell.hpp"
#include "Containers/FVector.hpp"
/**
* @author Olivier Coulaud (Olivier.Coulaud@inria.fr)
* @class FAdaptCell*
* @brief This class defines adaptive cell.
*
* ToDO
*
*/
//class FAdaptCell;
template <class CellClass, class LeafClass>
class FAdaptCell : public FBasicCell {
public:
class FExtACell {
public:
FAdaptCell<CellClass, LeafClass> * cell ;
int level ; //< Level in the octree of cell
FExtACell() { cell=nullptr ; level =-1 ; } ;
int getLevel() { return level;}
int getLevel() const { return level;}
};
protected:
// Global Index of the cell in the octree (This id is unique)
long int gID;
//! Number of particles inside the cell
int nbP ;
//! iAdaptive cell True, false
bool adaptive ;
bool sminMCriteria;
//
// Lists
//
FExtACell adaptiveParent ; //< adaptive Parent of the cell
FVector<FExtACell> adaptiveChild ; //< list of adaptive child of the cell
FVector<LeafClass*> leaves ; //< list of leaf child of the cell
//
//
CellClass * trueFMMCell ; //<a pointer on the cell that contains Multipole and local values
public:
FAdaptCell(): gID(-1), nbP(0), adaptive(false),sminMCriteria(false),trueFMMCell(new CellClass) {
}
/** Default destructor */
virtual ~FAdaptCell(){
}
//! add nbPart in cell
void addPart(const int n){
this->nbP += n ;
}
//! Return the number of particles inside the cell
const int getnbPart() const{
return this->nbP ;
}
//!return true if the sminM criteria is satisfied
bool isSminMCriteria (){
return this->sminMCriteria;
}
//!return true if the sminM criteria is satisfied
bool isSminMCriteria () const{
return this->sminMCriteria;
}
//! Set if the cell is adaptive or not
void setSminMCriteria(const bool bb=true){
this->sminMCriteria=bb ;
}
//! Return the global Id of the cell in the octree
const long int getGlobalId(){
return this->gID ;
}
const long int getGlobalId( ) const{
return this->gID ;
}
//! Set he global Id of the cell in the octree to id
void setGlobalId(const long int & id){
this->gID = id; ;
}
//! Set if the cell is adaptive or not
void setCelladaptive(const bool bb=true){
this->adaptive = bb;
}
//!return true if the cell is adaptive
const bool isadaptive() const{
return this->adaptive;
}
//
//! Add the adaptive parent of the cell
void addadaptiveParent( const FExtACell &ac) {
this->adaptiveParent.cell = ac.cell ;
this->adaptiveParent.level = ac.level ;
// std::cout << " addadaptiveParent " << *(this) <<std::endl;
}
//
//! Add the adaptive child of the cell
void addadaptiveParent( FAdaptCell<CellClass, LeafClass> * cell ,const int idxLevel) {
this->adaptiveParent.cell = cell ; this->adaptiveParent.level = idxLevel ;
}
FExtACell getadaptiveFather() {
return this->adaptiveParent ;
}
FExtACell getadaptiveFather() const {
return this->adaptiveParent ;
}
//
//! Add the adaptive child of the cell
void addadaptiveChild( FAdaptCell<CellClass, LeafClass> * cell ,const int idxLevel) {
FExtACell AC ;
AC.cell = cell ; AC.level = idxLevel ;
this->adaptiveChild.push(AC);
}
//
//! Add the adaptive child of the cell
void addadaptiveChild( const FExtACell &AC){
this->adaptiveChild.push(AC) ;
}
int sizeofadaptiveChild() const{
return this->adaptiveChild.getSize();
}
FVector<FExtACell>* getadaptiveChild(){
return &this->adaptiveChild;
}
const FVector<FExtACell>* getadaptiveChild() const{
return &this->adaptiveChild;
}
FExtACell* getAdaptiveChild(const int i) {
return &this->adaptiveChild[0];
}
const FVector<FExtACell> getAdaptiveChild() const{
return this->adaptiveChild;
}
//
//
CellClass* getKernelCell(){
return trueFMMCell ;
}
//
//! Add the adaptive child of the cell
void addLeafptr( LeafClass * leaf) {
this->leaves.push(leaf) ;
}
//! Return the number of leaves
int getLeavesSize() const{
return this->leaves.getSize();
}
//! Return the number of leaves
LeafClass* getLeaf(const int i ) {
return this->leaves[i];
}
//! Return the number of leaves
const LeafClass* getLeaf(const int i ) const {
return this->leaves[i];
}
/** Make it like the beginning */
void resetToInitialState(){
//this->dataDown = 0;
//this->dataUp = 0;
this->nbP = 0 ;
}
/////////////////////////////////////////////////
/** Save the current cell in a buffer */
template <class BufferWriterClass>
void save(BufferWriterClass& buffer) const{
FBasicCell::save(buffer);
//buffer << dataDown << dataUp << nbP;
}
/** Restore the current cell from a buffer */
template <class BufferReaderClass>
void restore(BufferReaderClass& buffer){
FBasicCell::restore(buffer);
// buffer >> dataDown >> dataUp >> nbP;
}
/** Get Multipole */
FReal* getMultipole(const int inRhs){
return this->trueFMMCell->getMultipole(inRhs);
}
/** Get Local */
FReal* getLocal(const int inRhs){
return this->trueFMMCell->getLocal(inRhs);
}
/** Get Multipole */
const FReal* getMultipole(const int inRhs) const
{ return this->trueFMMCell->getMultipole(inRhs);
}
/** Get Local */
const FReal* getLocal(const int inRhs) const{
return this->trueFMMCell->getLocal(inRhs);
}
// void resetToInitialState(){
// return this->trueFMMCell->resetToInitialState();
// }
/////////////////////////////////////////////////
/** Serialize only up data in a buffer */
template <class BufferWriterClass>
void serializeUp(BufferWriterClass& buffer) const {
//buffer << this->dataUp;
}
/** Deserialize only up data in a buffer */
template <class BufferReaderClass>
void deserializeUp(BufferReaderClass& buffer){
//buffer >> this->dataUp;
}
/** Serialize only down data in a buffer */
template <class BufferWriterClass>
void serializeDown(BufferWriterClass& buffer) const {
//buffer << this->dataDown;
}
/** Deserialize only up data in a buffer */
template <class BufferReaderClass>
void deserializeDown(BufferReaderClass& buffer){
//buffer >> this->dataDown;
}
/**
* Operator stream FAdaptCell to std::ostream
*
* @param[in,out] output where to write the adaptive cell
* @param[in] inPosition the cell to write out
* @return the output for multiple << operators
*/
template <class StreamClass>
friend StreamClass& operator<<(StreamClass& output, const FAdaptCell<CellClass,LeafClass>& cell){
output << "( Cell Id " << cell.getGlobalId() << " Adaptive " << std::boolalpha << cell.isadaptive()
<< " sminM " << cell.isSminMCriteria()<< " "<< cell.getnbPart() ;
if(cell.getLeavesSize() >0){
output << " LF={" ;
for (int i=0; i <cell.getLeavesSize() ; ++i){
output << cell.getLeaf(i) << " ";
}
output << "}" ;
}
output <<" CA={ ";
const FVector<FExtACell> * v =cell.getadaptiveChild() ;
if (cell.sizeofadaptiveChild()> 0 ){
for (int i=0; i < v->getSize() ; ++i){
output << v->operator [](i).cell->getGlobalId() << " ";
}
}
output << "} " ;
if(cell.getadaptiveFather().cell){
output << " FA={" << (cell.getadaptiveFather()).cell->getGlobalId() << "} " ;
}
else
{
output << " FA={} " ;
}
output << " )" <<std::endl;
return output; // for multiple << operators.
}
};
#endif //FADAPTCELL_HPP
#ifndef FADAPTCHEBSYMKERNEL_HPP
#define FADAPTCHEBSYMKERNEL_HPP
// ===================================================================================
// Copyright ScalFmm 2011 INRIA,
// 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".
// ===================================================================================
#include "Utils/FGlobal.hpp"
#include "Utils/FPoint.hpp"
#include "Adaptive/FAdaptiveCell.hpp"
#include "Adaptive/FAdaptiveKernelWrapper.hpp"
#include "Adaptive/FAbstractAdaptiveKernel.hpp"
#include "Kernels/Chebyshev/FChebSymKernel.hpp"
class FTreeCoordinate;
// ==== CMAKE =====
// @FUSE_BLAS
// ================
// for verbosity only!!!
//#define COUNT_BLOCKED_INTERACTIONS
// if timings should be logged
//#define LOG_TIMINGS
/**
* @author O. Coulaud
* @class FAdaptChebSymKernel
* @brief
* Please read the license
*
* This kernels implement the Chebyshev interpolation based FMM operators
* exploiting the symmetries in the far-field. It implements all interfaces
* (P2P, P2M, M2M, M2L, L2L, L2P) which are required by the FFmmAlgorithm and
* FFmmAlgorithmThread.
*
* @tparam CellClass Type of cell
* @tparam ContainerClass Type of container to store particles
* @tparam MatrixKernelClass Type of matrix kernel function
* @tparam ORDER Chebyshev interpolation order
*/
template< class CellClass, class ContainerClass, class MatrixKernelClass, int ORDER, int NVALS = 1>
class FAdaptiveChebSymKernel : FChebSymKernel<FReal,CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS>
, public FAbstractAdaptiveKernel<CellClass, ContainerClass> {
//
typedef FChebSymKernel<FReal,CellClass, ContainerClass, MatrixKernelClass, ORDER, NVALS> KernelBaseClass;
enum {order = ORDER,
nnodes = TensorTraits<ORDER>::nnodes};
const MatrixKernelClass *const MatrixKernel;
public:
using KernelBaseClass::P2M;
using KernelBaseClass::M2M;
using KernelBaseClass::M2L;
using KernelBaseClass::finishedLevelM2L;
using KernelBaseClass::L2L;
using KernelBaseClass::L2P;
using KernelBaseClass::P2P;
using KernelBaseClass::P2PRemote;
// /**
// * The constructor initializes all constant attributes and it reads the
// * precomputed and compressed M2L operators from a binary file (an
// * runtime_error is thrown if the required file is not valid).
// */
FAdaptiveChebSymKernel(const int inTreeHeight, const FReal inBoxWidth,
const FPoint<FReal>& inBoxCenter, const MatrixKernelClass *const inMatrixKernel) : KernelBaseClass(inTreeHeight, inBoxWidth, inBoxCenter, inMatrixKernel), MatrixKernel(inMatrixKernel)
{}
// /** Copy constructor */
FAdaptiveChebSymKernel(const FAdaptiveChebSymKernel& other)
: KernelBaseClass(other), MatrixKernel(other.MatrixKernel)
{ }
//
// /** Destructor */
~FAdaptiveChebSymKernel()
{
//this->~KernelBaseClass() ;
}
void P2M(CellClass* const pole, const int cellLevel, const ContainerClass* const particles) override {
const FPoint<FReal> CellCenter(KernelBaseClass::getCellCenter(pole->getCoordinate(),cellLevel));
const FReal BoxWidth = KernelBaseClass::BoxWidth / FMath::pow(2.0,cellLevel);
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// 1) apply Sy
KernelBaseClass::Interpolator->applyP2M(CellCenter, BoxWidth,
pole->getMultipole(idxRhs), particles);
}
}
void M2M(CellClass* const pole, const int poleLevel, const CellClass* const subCell, const int subCellLevel) override {
const FPoint<FReal> subCellCenter(KernelBaseClass::getCellCenter(subCell->getCoordinate(),subCellLevel));
const FReal subCellWidth(KernelBaseClass::BoxWidth / FReal(FMath::pow(2.0,subCellLevel)));
const FPoint<FReal> poleCellCenter(KernelBaseClass::getCellCenter(pole->getCoordinate(),poleLevel));
const FReal poleCellWidth(KernelBaseClass::BoxWidth / FReal(FMath::pow(2.0, poleLevel)));
////////////////////////////////////////////////////////////////////////////
/// p^6 version
// allocate memory
FReal* subChildParentInterpolator = new FReal [nnodes * nnodes];
// set child info
FPoint<FReal> ChildRoots[nnodes], localChildRoots[nnodes];
FChebTensor<FReal,ORDER>::setRoots(subCellCenter, subCellWidth, ChildRoots);
// map global position of roots to local position in parent cell
const map_glob_loc map(poleCellCenter, poleCellWidth);
for (unsigned int n=0; n<nnodes; ++n)
map(ChildRoots[n], localChildRoots[n]);
// assemble child - parent - interpolator
KernelBaseClass::Interpolator->assembleInterpolator(nnodes, localChildRoots, subChildParentInterpolator);
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// 1) apply Sy (using tensor product M2M with an interpolator computed on the fly)
// Do NOT reset multipole expansion !
//FBlas::scal(nnodes, FReal(0.), pole->getMultipole(idxRhs));
/// p^6 version
FBlas::gemtva(nnodes, nnodes, FReal(1.),
subChildParentInterpolator,
const_cast<FReal*>(subCell->getMultipole(idxRhs)), pole->getMultipole(idxRhs));
}// NVALS
}
void P2L(CellClass* const local, const int localLevel, const ContainerClass* const particles) override {
// Target cell: local
const FReal localCellWidth(KernelBaseClass::BoxWidth / FReal(FMath::pow(2.0, localLevel)));
const FPoint<FReal> localCellCenter(KernelBaseClass::getCellCenter(local->getCoordinate(),localLevel));
// interpolation points of target (X) cell
FPoint<FReal> X[nnodes];
FChebTensor<FReal,order>::setRoots(localCellCenter, localCellWidth, X);
// read positions
const FReal*const positionsX = particles->getPositions()[0];
const FReal*const positionsY = particles->getPositions()[1];
const FReal*const positionsZ = particles->getPositions()[2];
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// read physicalValue
const FReal*const physicalValues = particles->getPhysicalValues();
// apply P2L
for (unsigned int idxPart=0; idxPart<particles->getNbParticles(); ++idxPart){
const FPoint<FReal> y = FPoint<FReal>(positionsX[idxPart],
positionsY[idxPart],
positionsZ[idxPart]);
for (unsigned int m=0; m<nnodes; ++m)
local->getLocal(idxRhs)[m]+=MatrixKernel->evaluate(X[m], y) * physicalValues[idxPart];
}
}// NVALS
}
void M2L(CellClass* const local, const int localLevel, const CellClass* const pole, const int poleLevel) override {
// Source cell: pole
const FReal poleCellWidth(KernelBaseClass::BoxWidth / FReal(FMath::pow(2.0, poleLevel)));
const FPoint<FReal> poleCellCenter(KernelBaseClass::getCellCenter(pole->getCoordinate(),poleLevel));
// Target cell: local
const FReal localCellWidth(KernelBaseClass::BoxWidth / FReal(FMath::pow(2.0, localLevel)));
const FPoint<FReal> localCellCenter(KernelBaseClass::getCellCenter(local->getCoordinate(),localLevel));
// interpolation points of source (Y) and target (X) cell
FPoint<FReal> X[nnodes], Y[nnodes];
FChebTensor<FReal,order>::setRoots(poleCellCenter, poleCellWidth, Y);
FChebTensor<FReal,order>::setRoots(localCellCenter, localCellWidth, X);
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// Dense M2L
const FReal *const MultipoleExpansion = pole->getMultipole(idxRhs);
for (unsigned int m=0; m<nnodes; ++m)
for (unsigned int n=0; n<nnodes; ++n){
local->getLocal(idxRhs)[m]+=MatrixKernel->evaluate(X[m], Y[n]) * MultipoleExpansion[n];
}
}
}
void M2P(const CellClass* const pole, const int poleLevel, ContainerClass* const particles) override {
// Source cell: pole
const FReal poleCellWidth(KernelBaseClass::BoxWidth / FReal(FMath::pow(2.0, poleLevel)));
const FPoint<FReal> poleCellCenter(KernelBaseClass::getCellCenter(pole->getCoordinate(),poleLevel));
// interpolation points of source (Y) cell
FPoint<FReal> Y[nnodes];
FChebTensor<FReal,order>::setRoots(poleCellCenter, poleCellWidth, Y);
// read positions
const FReal*const positionsX = particles->getPositions()[0];
const FReal*const positionsY = particles->getPositions()[1];
const FReal*const positionsZ = particles->getPositions()[2];
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
const FReal *const MultipoleExpansion = pole->getMultipole(idxRhs);
// apply M2P
for (unsigned int idxPart=0; idxPart<particles->getNbParticles(); ++idxPart){
const FPoint<FReal> x = FPoint<FReal>(positionsX[idxPart],positionsY[idxPart],positionsZ[idxPart]);
FReal targetValue=0.;
{
for (unsigned int n=0; n<nnodes; ++n)
targetValue +=
MatrixKernel->evaluate(x, Y[n]) * MultipoleExpansion[/*idxLhs*nnodes+*/n];
}
// get potential
FReal*const potentials = particles->getPotentials(/*idxPot*/);
// add contribution to potential
potentials[idxPart] += (targetValue);
}// Particles
}// NVALS
}
void L2L(const CellClass* const local, const int localLevel, CellClass* const subCell, const int subCellLevel) override {
const FPoint<FReal> subCellCenter(KernelBaseClass::getCellCenter(subCell->getCoordinate(),subCellLevel));
const FReal subCellWidth(KernelBaseClass::BoxWidth / FReal(FMath::pow(2.0,subCellLevel)));
const FPoint<FReal> localCenter(KernelBaseClass::getCellCenter(local->getCoordinate(),localLevel));
const FReal localWidth(KernelBaseClass::BoxWidth / FReal(FMath::pow(2.0,localLevel)));
////////////////////////////////////////////////////////////////////////////
/// p^6 version
// allocate memory
FReal* subChildParentInterpolator = new FReal [nnodes * nnodes];
// set child info
FPoint<FReal> ChildRoots[nnodes], localChildRoots[nnodes];
FChebTensor<FReal,ORDER>::setRoots(subCellCenter, subCellWidth, ChildRoots);
// map global position of roots to local position in parent cell
const map_glob_loc map(localCenter, localWidth);
for (unsigned int n=0; n<nnodes; ++n)
map(ChildRoots[n], localChildRoots[n]);
// assemble child - parent - interpolator
KernelBaseClass::Interpolator->assembleInterpolator(nnodes, localChildRoots, subChildParentInterpolator);
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){
// 2) apply Sx
/// p^6 version
FBlas::gemva(nnodes, nnodes, FReal(1.),
subChildParentInterpolator,
const_cast<FReal*>(local->getLocal(idxRhs)), subCell->getLocal(idxRhs));
}// NVALS
}
void L2P(const CellClass* const local, const int cellLevel, ContainerClass* const particles) override {
const FPoint<FReal> CellCenter(KernelBaseClass::getCellCenter(local->getCoordinate(),cellLevel));
const FReal BoxWidth = KernelBaseClass::BoxWidth / FMath::pow(2.0,cellLevel);
for(int idxRhs = 0 ; idxRhs < NVALS ; ++idxRhs){