Commit fed7d3b4 authored by COULAUD Olivier's avatar COULAUD Olivier

Fix compilation warnings and add update comments

Add
  - CutOffKernel which uses the matrixkernel class in the P2P Operator
  - CutOffAlgorithm to explain how we can use scalFMM with only cutOff method. The Octree dos not contain any local or multipole values.
parent a9aa9aee
// ===================================================================================
// Copyright ScalFmm 2016 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.
// An extension to the license is given to allow static linking of scalfmm
// inside a proprietary application (no matter its license).
// See the main license file for more details.
//
// 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".
// ===================================================================================
// ==== CMAKE =====
// @FUSE_FFT
// @FUSE_BLAS
// ==== Git =====
// @SCALFMM_PRIVATE
// ================
#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <string>
#include "ScalFmmConfig.h"
#include "Utils/FGlobal.hpp"
#include "Utils/FParameters.hpp"
#include "Utils/FParameterNames.hpp"
// Our particles loader
#include "Files/FFmaGenericLoader.hpp"
//
// Cells without local and multipole values
#include "Components/FBasicCell.hpp"
// Leaves
#include "Components/FSimpleLeaf.hpp"
#include "Kernels/P2P/FP2PParticleContainerIndexed.hpp"
#include "Kernels/Interpolation/FInterpMatrixKernel.hpp" // For the predefined Matrix kernels
#include "Kernels/Interpolation/FCutOffKernel.hpp" // Generic CutOff kernel defined with Matrix Kernel
// The classical Octree
#include "Containers/FOctree.hpp"
#ifdef _OPENMP
#include "Core/FFmmAlgorithmThread.hpp"
#else
#include "Core/FFmmAlgorithm.hpp"
#endif
#include <memory>
/**
* This program runs the only the P2P Algorithm with the cutoff kernel (1/r) and compares the results with a direct computation.
*/
/// \file CutOffAlhorithm.cpp
//!
//! \brief This program runs the only the P2P Algorithm with the cutoff kernel (1/r) and compares the results with a direct computation.
//! \authors O. Coulaud
//!
//! This code is a short example to use ScalFMM with only P2P algorithm for the 1/r kernel
// Simply create particles and try the kernels
int main(int argc, char* argv[])
{
FHelpDescribeAndExit(argc, argv,
"Driver for CutOff kernel (1/r kernel).",
FParameterDefinitions::InputFile, FParameterDefinitions::OctreeHeight,
FParameterDefinitions::OctreeSubHeight, FParameterDefinitions::InputFile,FParameterDefinitions::OutputFile,
FParameterDefinitions::NbThreads);
const std::string defaultFile(SCALFMMDataPath+"unitCubeXYZQ100.bfma" );
const std::string filename = FParameters::getStr(argc,argv,FParameterDefinitions::InputFile.options, defaultFile.c_str());
const unsigned int TreeHeight = FParameters::getValue(argc, argv, FParameterDefinitions::OctreeHeight.options, 2);
const unsigned int SubTreeHeight = FParameters::getValue(argc, argv, FParameterDefinitions::OctreeSubHeight.options, 1);
const unsigned int NbThreads = FParameters::getValue(argc, argv, FParameterDefinitions::NbThreads.options, 4);
#ifdef _OPENMP
omp_set_num_threads(NbThreads);
std::cout << "\n>> Using " << omp_get_max_threads() << " threads.\n" << std::endl;
#else
std::cout << "\n>> Sequential version.\n" << std::endl;
#endif
//
std::cout << "Parameters "<< std::endl
<< " Octree Depth "<< TreeHeight <<std::endl
<< " SubOctree depth " << SubTreeHeight <<std::endl
<< " Input file name: " <<filename <<std::endl
<< " Thread number: " << NbThreads <<std::endl
<<std::endl;
//
// init timer
FTic time;
// open particle file
////////////////////////////////////////////////////////////////////
//
typedef double FReal;
FFmaGenericLoader<FReal> loader(filename);
//
////////////////////////////////////////////////////////////////////
// begin CutOff kernel
// accuracy
// typedefs
typedef FP2PParticleContainerIndexed<FReal> ContainerClass;
typedef FSimpleLeaf<FReal, ContainerClass > LeafClass;
typedef FBasicCell CellClass;
typedef FOctree<FReal, CellClass,ContainerClass,LeafClass> OctreeClass;
//
typedef FInterpMatrixKernelR<FReal> MatrixKernelClass;
const MatrixKernelClass MatrixKernel;
typedef FCutOffKernel<FReal,CellClass,ContainerClass,MatrixKernelClass> KernelClass;
//
#ifdef _OPENMP
typedef FFmmAlgorithmThread<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
#else
typedef FFmmAlgorithm<OctreeClass,CellClass,ContainerClass,KernelClass,LeafClass> FmmClass;
#endif
// init oct-tree
OctreeClass tree(TreeHeight, SubTreeHeight, loader.getBoxWidth(), loader.getCenterOfBox());
{ // -----------------------------------------------------
std::cout << "Creating & Inserting " << loader.getNumberOfParticles()
<< " particles ..." << std::endl;
std::cout << "\tHeight : " << TreeHeight << " \t sub-height : " << SubTreeHeight << std::endl;
time.tic();
//
FPoint<FReal> position;
FReal physicalValue = 0.0;
//
for(FSize idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
//
// Read particle per particle from file
loader.fillParticle(&position,&physicalValue);
//
// put particle in octree
tree.insert(position, idxPart, physicalValue);
}
time.tac();
std::cout << "Done " << "(@Creating and Inserting Particles = "
<< time.elapsed() << " s) ." << std::endl;
} // -----------------------------------------------------
{ // -----------------------------------------------------
std::cout << " CutOff Algorithm ... " << std::endl;
time.tic();
//
std::unique_ptr<KernelClass> kernels(new KernelClass(TreeHeight, loader.getBoxWidth(), loader.getCenterOfBox(),&MatrixKernel));
//
FmmClass algo(&tree, kernels.get());
//
algo.execute( FFmmP2P); // Here the call of the FMM algorithm
//
time.tac();
std::cout << "Timers Far Field \n"
<< "P2M " << algo.getTime(FAlgorithmTimers::P2MTimer) << " seconds\n"
<< "M2M " << algo.getTime(FAlgorithmTimers::M2MTimer) << " seconds\n"
<< "M2L " << algo.getTime(FAlgorithmTimers::M2LTimer) << " seconds\n"
<< "L2L " << algo.getTime(FAlgorithmTimers::L2LTimer) << " seconds\n"
<< "P2P and L2P " << algo.getTime(FAlgorithmTimers::NearTimer) << " seconds\n"
<< std::endl;
std::cout << "Done " << "(@Algorithm = " << time.elapsed() << " s) ." << std::endl;
}
// -----------------------------------------------------
//
// Some output
//
//
{ // -----------------------------------------------------
FSize N1=0, N2= loader.getNumberOfParticles()/2, N3= loader.getNumberOfParticles() -1; ;
FReal energy =0.0 ;
//
// Loop over all leaves
//
std::cout <<std::endl<<" &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& "<<std::endl;
std::cout << std::scientific;
std::cout.precision(10) ;
tree.forEachLeaf([&](LeafClass* leaf){
const FReal*const posX = leaf->getTargets()->getPositions()[0];
const FReal*const posY = leaf->getTargets()->getPositions()[1];
const FReal*const posZ = leaf->getTargets()->getPositions()[2];
const FReal*const potentials = leaf->getTargets()->getPotentials();
const FReal*const forcesX = leaf->getTargets()->getForcesX();
const FReal*const forcesY = leaf->getTargets()->getForcesY();
const FReal*const forcesZ = leaf->getTargets()->getForcesZ();
const FSize nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();
const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
const FSize indexPartOrig = indexes[idxPart];
if ((indexPartOrig == N1) || (indexPartOrig == N2) || (indexPartOrig == N3) ) {
std::cout << "Index "<< indexPartOrig <<" potential " << potentials[idxPart]
<< " Pos "<<posX[idxPart]<<" "<<posY[idxPart]<<" "<<posZ[idxPart]
<< " Forces: " << forcesX[idxPart] << " " << forcesY[idxPart] << " "<< forcesZ[idxPart] <<std::endl;
}
energy += potentials[idxPart]*physicalValues[idxPart] ;
}
});
std::cout <<std::endl<<"Energy: "<< energy<<std::endl;
std::cout <<std::endl<<" &&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& "<<std::endl<<std::endl;
}
// -----------------------------------------------------
if(FParameters::existParameter(argc, argv, FParameterDefinitions::OutputFile.options)){
std::string name(FParameters::getStr(argc,argv,FParameterDefinitions::OutputFile.options, "output.fma"));
FFmaGenericWriter<FReal> writer(name) ;
//
FSize NbPoints = loader.getNumberOfParticles();
FReal * particles ;
particles = new FReal[8*NbPoints] ;
memset(particles,0,8*NbPoints*sizeof(FReal));
FSize j = 0 ;
tree.forEachLeaf([&](LeafClass* leaf){
//
// Input
const FReal*const posX = leaf->getTargets()->getPositions()[0];
const FReal*const posY = leaf->getTargets()->getPositions()[1];
const FReal*const posZ = leaf->getTargets()->getPositions()[2];
const FReal*const physicalValues = leaf->getTargets()->getPhysicalValues();
const FVector<FSize>& indexes = leaf->getTargets()->getIndexes();
//
// Computed data
const FReal*const potentials = leaf->getTargets()->getPotentials();
const FReal*const forcesX = leaf->getTargets()->getForcesX();
const FReal*const forcesY = leaf->getTargets()->getForcesY();
const FReal*const forcesZ = leaf->getTargets()->getForcesZ();
//
const FSize nbParticlesInLeaf = leaf->getTargets()->getNbParticles();
for(FSize idxPart = 0 ; idxPart < nbParticlesInLeaf ; ++idxPart){
j = 8*indexes[idxPart];
particles[j] = posX[idxPart] ;
particles[j+1] = posY[idxPart] ;
particles[j+2] = posZ[idxPart] ;
particles[j+3] = physicalValues[idxPart] ;
particles[j+4] = potentials[idxPart] ;
particles[j+5] = forcesX[idxPart] ;
particles[j+6] = forcesY[idxPart] ;
particles[j+7] = forcesZ[idxPart] ;
}
});
writer.writeHeader( loader.getCenterOfBox(), loader.getBoxWidth() , NbPoints, sizeof(FReal), 8) ;
writer.writeArrayOfReal(particles, 8 , NbPoints);
delete[] particles;
//
std::string name1( "output.fma");
//
FFmaGenericWriter<FReal> writer1(name1) ;
writer1.writeDistributionOfParticlesFromOctree(&tree,NbPoints) ;
}
return 0;
}
......@@ -52,118 +52,118 @@
// Simply create particles and try the kernels
int main(int argc, char ** argv){
FHelpDescribeAndExit(argc, argv,
">> This executable has to be used to compute interaction either for periodic or non periodic system.\n"
">> Example -fin filenameIN.{fma or bfma) -fout filenameOUT{fma or bfma) \n"
">> Default input file : Data/unitCubeXYZQ20k.fma\n"
">> Only our FMA (.bma, .bfma) is allowed as input.\n"
">> Output file with extension (default output.bfma).",
FParameterDefinitions::InputFile, FParameterDefinitions::OutputFile,
FParameterDefinitions::EnabledVerbose);
//////////////////////////////////////////////////////////////
typedef double FReal;
const std::string defaultFile(/*SCALFMMDataPath+*/"../Data/unitCubeXYZQ20k.fma");
const std::string filenameIn(FParameters::getStr(argc,argv,FParameterDefinitions::InputFile.options, defaultFile.c_str()));
const std::string filenameOut(FParameters::getStr(argc,argv,FParameterDefinitions::OutputFile.options, "output.bfma"));
//
FTic counter;
// -----------------------------------------------------
// LOADER
// -----------------------------------------------------
// ---------------------------------------------------------------------------------
// Read particles in the Octree
// ---------------------------------------------------------------------------------
std::cout << "Opening : " << filenameIn << "\n";
//
FFmaGenericLoader<FReal> loader(filenameIn);
//
FSize nbParticles = static_cast<int>(loader.getNumberOfParticles());
std::cout << "Read " << nbParticles << " particles ..." << std::endl;
double BoxWith=loader.getBoxWidth();
FPoint<FReal> Centre(loader.getCenterOfBox().getX(), loader.getCenterOfBox().getY() , loader.getCenterOfBox().getZ());
std::cout << "\tWidth : " <<BoxWith << " \t center x : " << loader.getCenterOfBox().getX()
<< " y : " << loader.getCenterOfBox().getY() << " z : " << loader.getCenterOfBox().getZ() << std::endl;
FHelpDescribeAndExit(argc, argv,
">> This executable has to be used to compute interaction either for periodic or non periodic system.\n"
">> Example -fin filenameIN.{fma or bfma) -fout filenameOUT{fma or bfma) \n"
">> Default input file : Data/unitCubeXYZQ20k.fma\n"
">> Only our FMA (.bma, .bfma) is allowed as input.\n"
">> Output file with extension (default output.bfma).",
FParameterDefinitions::InputFile, FParameterDefinitions::OutputFile,
FParameterDefinitions::EnabledVerbose);
//////////////////////////////////////////////////////////////
typedef double FReal;
const std::string defaultFile(/*SCALFMMDataPath+*/"../Data/unitCubeXYZQ20k.fma");
const std::string filenameIn(FParameters::getStr(argc,argv,FParameterDefinitions::InputFile.options, defaultFile.c_str()));
const std::string filenameOut(FParameters::getStr(argc,argv,FParameterDefinitions::OutputFile.options, "output.bfma"));
//
FTic counter;
// -----------------------------------------------------
// LOADER
// -----------------------------------------------------
// ---------------------------------------------------------------------------------
// Read particles in the Octree
// ---------------------------------------------------------------------------------
std::cout << "Opening : " << filenameIn << "\n";
//
FFmaGenericLoader<FReal> loader(filenameIn);
//
FSize nbParticles = static_cast<int>(loader.getNumberOfParticles());
std::cout << "Read " << nbParticles << " particles ..." << std::endl;
double BoxWith=loader.getBoxWidth();
FPoint<FReal> Centre(loader.getCenterOfBox().getX(), loader.getCenterOfBox().getY() , loader.getCenterOfBox().getZ());
std::cout << "\tWidth : " <<BoxWith << " \t center x : " << loader.getCenterOfBox().getX()
<< " y : " << loader.getCenterOfBox().getY() << " z : " << loader.getCenterOfBox().getZ() << std::endl;
counter.tic();
//
FmaRWParticle<FReal, 4,8> * particles = new FmaRWParticle<FReal, 4,8>[nbParticles];
memset(particles, 0, sizeof(FmaRWParticle<FReal, 4,8>) * nbParticles) ;
//
double totalCharge = 0.0;
//
// int nbDataToRead = particles[0].getReadDataNumber();
for(int idx = 0 ; idx<nbParticles ; ++idx){
//
loader.fillParticle(particles[idx].getPtrFirstData(), particles[idx].getReadDataNumber());
// loader.fillParticle(particles[idx].getPtrFirstData(), nbDataToRead); // OK
// loader.fillParticle(particles[idx]); // OK
// std::cout << idx <<" "<< particles[idx].getPosition() << " "<<particles[idx].getPhysicalValue() << " "<<particles[idx].getPotential()
// <<" " << particles[idx].getForces()[0]<<" " <<particles[idx].getForces()[1]<<" " <<particles[idx].getForces()[2]<<" " <<std::endl;
//
totalCharge += particles[idx].getPhysicalValue() ;
}
counter.tac();
std::cout << std::endl;
std::cout << "Total Charge = "<< totalCharge <<std::endl;
std::cout << std::endl;
std::cout << "Done " << "(@ reading Particles " << counter.elapsed() << " s)." << std::endl;
//
// ----------------------------------------------------------------------------------------------------------
// COMPUTATION
// ----------------------------------------------------------------------------------------------------------
// interaction kernel evaluator
typedef FInterpMatrixKernelR<FReal> MatrixKernelClass;
const MatrixKernelClass MatrixKernel;
FReal denergy = 0.0;
//
// computation
//
{
printf("Compute :\n");
counter.tic();
#pragma omp parallel shared(nbParticles, particles,denergy)
{
#pragma omp for
for(int idxTarget = 0 ; idxTarget < nbParticles ; ++idxTarget){
//
// compute with all other except itself
//
// Compute force and potential between particles[idxTarget] and particles inside the box
//
for(int idxOther = 0; idxOther < nbParticles ; ++idxOther){
if( idxOther != idxTarget ){
FP2P::NonMutualParticles(particles[idxTarget].getPosition().getX(), particles[idxTarget].getPosition().getY(),
particles[idxTarget].getPosition().getZ(),particles[idxTarget].getPhysicalValue(),
&particles[idxTarget].setForces()[0],&particles[idxTarget].setForces()[1],
&particles[idxTarget].setForces()[2],particles[idxTarget].setPotential(),
particles[idxOther].getPosition().getX(), particles[idxOther].getPosition().getY(),
particles[idxOther].getPosition().getZ(),particles[idxOther].getPhysicalValue(),
&MatrixKernel);
}
}
} // end for
// Compute the energy
#pragma omp for reduction(+:denergy)
for(int idx = 0 ; idx < nbParticles ; ++idx){
denergy += particles[idx].getPotential()*(particles[idx].getPhysicalValue()) ;
}
} // end pragma parallel
//
FmaRWParticle<FReal, 4,8> * particles = new FmaRWParticle<FReal, 4,8>[nbParticles];
memset(particles, 0, sizeof(FmaRWParticle<FReal, 4,8>) * nbParticles) ;
//
double totalCharge = 0.0;
//
// int nbDataToRead = particles[0].getReadDataNumber();
for(int idx = 0 ; idx<nbParticles ; ++idx){
//
loader.fillParticle(particles[idx].getPtrFirstData(), particles[idx].getReadDataNumber());
// loader.fillParticle(particles[idx].getPtrFirstData(), nbDataToRead); // OK
// loader.fillParticle(particles[idx]); // OK
// std::cout << idx <<" "<< particles[idx].getPosition() << " "<<particles[idx].getPhysicalValue() << " "<<particles[idx].getPotential()
// <<" " << particles[idx].getForces()[0]<<" " <<particles[idx].getForces()[1]<<" " <<particles[idx].getForces()[2]<<" " <<std::endl;
//
totalCharge += particles[idx].getPhysicalValue() ;
}
denergy *= 0.5 ;
counter.tac();
std::cout << std::endl;
std::cout << "Total Charge = "<< totalCharge <<std::endl;
std::cout << std::endl;
std::cout << "Done " << "(@ reading Particles " << counter.elapsed() << " s)." << std::endl;
//
// ----------------------------------------------------------------------------------------------------------
// COMPUTATION
// ----------------------------------------------------------------------------------------------------------
// interaction kernel evaluator
typedef FInterpMatrixKernelR<FReal> MatrixKernelClass;
const MatrixKernelClass MatrixKernel;
FReal denergy = 0.0;
//
// computation
printf("Energy = %.14e\n",denergy);
//
{
printf("Compute :\n");
counter.tic();
#pragma omp parallel shared(nbParticles, particles,denergy)
{
#pragma omp for
for(int idxTarget = 0 ; idxTarget < nbParticles ; ++idxTarget){
//
// compute with all other except itself
//
// Compute force and potential between particles[idxTarget] and particles inside the box
//
for(int idxOther = 0; idxOther < nbParticles ; ++idxOther){
if( idxOther != idxTarget ){
FP2P::NonMutualParticles(particles[idxTarget].getPosition().getX(), particles[idxTarget].getPosition().getY(),
particles[idxTarget].getPosition().getZ(),particles[idxTarget].getPhysicalValue(),
&particles[idxTarget].setForces()[0],&particles[idxTarget].setForces()[1],
&particles[idxTarget].setForces()[2],particles[idxTarget].setPotential(),
particles[idxOther].getPosition().getX(), particles[idxOther].getPosition().getY(),
particles[idxOther].getPosition().getZ(),particles[idxOther].getPhysicalValue(),
&MatrixKernel);
}
}
} // end for
// Compute the energy
#pragma omp for reduction(+:denergy)
for(int idx = 0 ; idx < nbParticles ; ++idx){
denergy += particles[idx].getPotential()*(particles[idx].getPhysicalValue()) ;
}
} // end pragma parallel
//
denergy *= 0.5 ;
counter.tac();
//
printf("Energy = %.14e\n",denergy);
//
std::cout << "Done " << "(@ Direct computation done = " << counter.elapsed() << " s)." << std::endl;
std::cout << "\n"<< "END "
<< "-------------------------------------------------------------------------"
<< std::endl << std::endl ;
} // END
std::cout << "Done " << "(@ Direct computation done = " << counter.elapsed() << " s)." << std::endl;
std::cout << "\n"<< "END "
<< "-------------------------------------------------------------------------"
<< std::endl << std::endl ;
} // END
//
// ----------------------------------------------------------------
......@@ -171,34 +171,34 @@ int main(int argc, char ** argv){
//
//
std::cout << "Generate " << filenameOut <<" for output file" << std::endl;
//
std::cout << " nbParticles: " << nbParticles <<" " << sizeof(nbParticles) <<std::endl;
std::cout << " denergy: " << denergy <<" " << sizeof(denergy) <<std::endl;
std::cout << " Box size: " << loader.getBoxWidth() << " " << sizeof(loader.getBoxWidth())<<std::endl;
//
FFmaGenericWriter<FReal> writer(filenameOut) ;
writer.writeHeader(Centre,BoxWith, nbParticles,*particles) ;
writer.writeArrayOfParticles(particles, nbParticles);
//
// end generate
// -----------------------------------------------------
//
if(FParameters::existParameter(argc, argv, FParameterDefinitions::EnabledVerbose.options)){
denergy = 0 ;
for(int idx = 0 ; idx < nbParticles ; ++idx){
std::cout << ">> index " << idx << std::endl;
std::cout << " x " << particles[idx].getPosition().getX() << " y " << particles[idx].getPosition().getY() << " z " << particles[idx].getPosition().getZ() << std::endl;
std::cout << " Q " << particles[idx].getPhysicalValue() << " V " << particles[idx].getPotential() << std::endl;
std::cout << " fx " << particles[idx].getForces()[0] << " fy " << particles[idx].getForces()[1] << " fz " << particles[idx].getForces()[2] << std::endl;
std::cout << "\n";
denergy += particles[idx].getPotential()*particles[idx].getPhysicalValue() ;
}
std::cout << "Generate " << filenameOut <<" for output file" << std::endl;
//
std::cout << " nbParticles: " << nbParticles <<" " << sizeof(nbParticles) <<std::endl;
std::cout << " denergy: " << denergy <<" " << sizeof(denergy) <<std::endl;
std::cout << " Box size: " << loader.getBoxWidth() << " " << sizeof(loader.getBoxWidth())<<std::endl;
//
FFmaGenericWriter<FReal> writer(filenameOut) ;
writer.writeHeader(Centre,BoxWith, nbParticles,*particles) ;
writer.writeArrayOfParticles(particles, nbParticles);
//
// end generate
// -----------------------------------------------------
//
if(FParameters::existParameter(argc, argv, FParameterDefinitions::EnabledVerbose.options)){
denergy = 0 ;
for(int idx = 0 ; idx < nbParticles ; ++idx){
std::cout << ">> index " << idx << std::endl;
std::cout << " x " << particles[idx].getPosition().getX() << " y " << particles[idx].getPosition().getY() << " z " << particles[idx].getPosition().getZ() << std::endl;
std::cout << " Q " << particles[idx].getPhysicalValue() << " V " << particles[idx].getPotential() << std::endl;
std::cout << " fx " << particles[idx].getForces()[0] << " fy " << particles[idx].getForces()[1] << " fz " << particles[idx].getForces()[2] << std::endl;
std::cout << "\n";
denergy += particles[idx].getPotential()*particles[idx].getPhysicalValue() ;
}
std::cout << " ENERGY " << denergy << std::endl;
//
delete[] particles;
return 0;
}
std::cout << " ENERGY " << denergy << std::endl;
//
delete[] particles;
return 0;
}
......
// ===================================================================================
// Copyright ScalFmm 2016 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.
// An extension to the license is given to allow static linking of scalfmm
// inside a proprietary application (no matter its license).
// See the main license file for more details.
//
// 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".
// ===================================================================================
// Keep in private GIT
#ifndef FCutOffKernel_HPP
#define FCutOffKernel_HPP
#include "Utils/FGlobal.hpp"
#include "Components/FAbstractKernels.hpp"