diff --git a/Examples/generateDistributions.cpp b/Examples/generateDistributions.cpp
index 57affde649902c03c567e0895d0481bcdd6f376c..cffab37cdd8f1b7344783991e3005d63680a1fec 100644
--- a/Examples/generateDistributions.cpp
+++ b/Examples/generateDistributions.cpp
@@ -102,7 +102,8 @@ int main(int argc, char ** argv){
" -filename name: generic name for files (without extension) and save data\n"
" with following format in name.fma or name.bfma in -bin is set\n"
" -visufmt vtk, vtp, cosmo or cvs format.",
- FParameterDefinitions::InputFile, FParameterDefinitions::NbParticles);
+ FParameterDefinitions::InputFile, FParameterDefinitions::OutputBinFormat, FParameterDefinitions::NbParticles);
+
@@ -164,6 +165,7 @@ int main(int argc, char ** argv){
BoxWith = FMath::Max(A,FMath::Max(B,C) );
FReal halfBW = BoxWith*0.5;
Centre.setPosition(halfBW,halfBW,halfBW);
+ std::cout << "Cuboid "<< A << ":"<< B<<":"<<C<<std::endl;
}
else if(FParameters::existParameter(argc, argv, "-unitSphere")){
unifRandonPointsOnUnitSphere(NbPoints, particles) ;
diff --git a/Src/Containers/FParForEachOctree.hpp b/Src/Containers/FParForEachOctree.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..21b3b1523cdb8651009c14e87d3c278a5b8da44c
--- /dev/null
+++ b/Src/Containers/FParForEachOctree.hpp
@@ -0,0 +1,141 @@
+#ifndef FPARFOREACHOCTREE_HPP
+#define FPARFOREACHOCTREE_HPP
+
+#include "FOctree.hpp"
+
+#include <functional>
+
+namespace FParForEachOctree {
+
+/////////////////////////////////////////////////////////
+// Lambda function to apply to all member
+/////////////////////////////////////////////////////////
+
+/**
+ * @brief forEachLeaf iterate on the leaf and apply the function
+ * @param function
+ */
+template< template< class CellClass, class ContainerClass, class LeafClass, class CellAllocatorClass > class FOctree,
+ class CellClass, class ContainerClass, class LeafClass, class CellAllocatorClass,
+ class FunctionTemplate>
+void forEachLeaf(FOctree<CellClass,ContainerClass,LeafClass,CellAllocatorClass>* tree, FunctionTemplate function){
+ #pragma omp parallel
+ {
+ #pragma omp single
+ {
+ typename FOctree<CellClass,ContainerClass,LeafClass,CellAllocatorClass>::Iterator octreeIterator(tree);
+ octreeIterator.gotoBottomLeft();
+
+ do{
+ LeafClass* lf = octreeIterator.getCurrentLeaf();
+ #pragma omp task firstprivate(lf)
+ {
+ function(lf);
+ }
+ } while(octreeIterator.moveRight());
+
+ #pragma omp taskwait
+ }
+ }
+}
+
+/**
+ * @brief forEachLeaf iterate on the cell and apply the function
+ * @param function
+ */
+template< template< class CellClass, class ContainerClass, class LeafClass, class CellAllocatorClass > class FOctree,
+ class CellClass, class ContainerClass, class LeafClass, class CellAllocatorClass,
+ class FunctionTemplate>
+void forEachCell(FOctree<CellClass,ContainerClass,LeafClass,CellAllocatorClass>* tree, FunctionTemplate function){
+ #pragma omp parallel
+ {
+ #pragma omp single
+ {
+ typename FOctree<CellClass,ContainerClass,LeafClass,CellAllocatorClass>::Iterator octreeIterator(tree);
+ octreeIterator.gotoBottomLeft();
+
+ typename FOctree<CellClass,ContainerClass,LeafClass,CellAllocatorClass>::Iterator avoidGoLeft(octreeIterator);
+
+ for(int idx = tree->getHeight()-1 ; idx >= 1 ; --idx ){
+ do{
+ CellClass* cell = octreeIterator.getCurrentCell();
+ #pragma omp task firstprivate(cell)
+ {
+ function(cell);
+ }
+ } while(octreeIterator.moveRight());
+ avoidGoLeft.moveUp();
+ octreeIterator = avoidGoLeft;
+ }
+
+ #pragma omp taskwait
+ }
+ }
+}
+
+/**
+ * @brief forEachLeaf iterate on the cell and apply the function
+ * @param function
+ */
+template< template< class CellClass, class ContainerClass, class LeafClass, class CellAllocatorClass > class FOctree,
+ class CellClass, class ContainerClass, class LeafClass, class CellAllocatorClass,
+ class FunctionTemplate>
+void forEachCellWithLevel(FOctree<CellClass,ContainerClass,LeafClass,CellAllocatorClass>* tree, FunctionTemplate function){
+ #pragma omp parallel
+ {
+ #pragma omp single
+ {
+ typename FOctree<CellClass,ContainerClass,LeafClass,CellAllocatorClass>::Iterator octreeIterator(tree);
+ octreeIterator.gotoBottomLeft();
+
+ typename FOctree<CellClass,ContainerClass,LeafClass,CellAllocatorClass>::Iterator avoidGoLeft(octreeIterator);
+
+ for(int idx = tree->getHeight()-1 ; idx >= 1 ; --idx ){
+ do{
+ CellClass* cell = octreeIterator.getCurrentCell();
+ #pragma omp task firstprivate(cell, idx)
+ {
+ function(cell,idx);
+ }
+ } while(octreeIterator.moveRight());
+ avoidGoLeft.moveUp();
+ octreeIterator = avoidGoLeft;
+ }
+
+ #pragma omp taskwait
+ }
+ }
+}
+
+/**
+ * @brief forEachLeaf iterate on the cell and apply the function
+ * @param function
+ */
+template< template< class CellClass, class ContainerClass, class LeafClass, class CellAllocatorClass > class FOctree,
+ class CellClass, class ContainerClass, class LeafClass, class CellAllocatorClass,
+ class FunctionTemplate>
+void forEachCellLeaf(FOctree<CellClass,ContainerClass,LeafClass,CellAllocatorClass>* tree, FunctionTemplate function){
+ #pragma omp parallel
+ {
+ #pragma omp single
+ {
+ typename FOctree<CellClass,ContainerClass,LeafClass,CellAllocatorClass>::Iterator octreeIterator(tree);
+ octreeIterator.gotoBottomLeft();
+
+ do{
+ CellClass* cell = octreeIterator.getCurrentCell();
+ LeafClass* lf = octreeIterator.getCurrentLeaf();
+ #pragma omp task firstprivate(cell, lf)
+ {
+ function(cell,lf);
+ }
+ } while(octreeIterator.moveRight());
+
+ #pragma omp taskwait
+ }
+ }
+}
+
+}
+
+#endif // FPARFOREACHOCTREE_HPP
diff --git a/Src/Utils/FParObject.hpp b/Src/Utils/FParObject.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..be47ce30b878e2b2d116af72dcad2bb590147d7d
--- /dev/null
+++ b/Src/Utils/FParObject.hpp
@@ -0,0 +1,197 @@
+#ifndef FPAROBJECT_HPP
+#define FPAROBJECT_HPP
+
+#include <omp.h>
+#include <functional>
+
+/**
+ * This class is a way to hide concurent access to the user.
+ * Especially with the FParForEachOctree functions.
+ * Please refer to the testOctreeParallelFuncteur to see an example.
+ */
+template <class SharedClass>
+class FParObject{
+ const int nbThreads;
+ SharedClass* objects;
+
+ FParObject(const FParObject&) = delete;
+ FParObject& operator=(const FParObject&) = delete;
+
+public:
+
+ template <class... ArgsClass>
+ explicit FParObject(ArgsClass... args) : nbThreads(omp_get_max_threads()), objects(nullptr) {
+ objects = reinterpret_cast<SharedClass*>(new unsigned char[sizeof(SharedClass)*nbThreads]);
+ for(int idxThread = 0 ; idxThread < nbThreads ; ++idxThread){
+ new (&objects[idxThread]) SharedClass(args...);
+ }
+ }
+
+ ~FParObject(){
+ for(int idxThread = 0 ; idxThread < nbThreads ; ++idxThread){
+ objects[idxThread].~SharedClass();
+ }
+ delete[] reinterpret_cast<unsigned char*>(objects);
+ }
+
+ SharedClass& getMine(){
+ return objects[omp_get_thread_num()];
+ }
+
+ const SharedClass& getMine() const {
+ return objects[omp_get_thread_num()];
+ }
+
+ void applyToAll(std::function<void(SharedClass&)> func){
+ for(int idxThread = 0 ; idxThread < nbThreads ; ++idxThread){
+ func(objects[idxThread]);
+ }
+ }
+
+ void applyToAllPar(std::function<void(SharedClass&)> func){
+ #pragma omp parallel for num_threads(nbThreads) schedule(static)
+ for(int idxThread = 0 ; idxThread < nbThreads ; ++idxThread){
+ func(objects[idxThread]);
+ }
+ }
+
+ SharedClass reduce(std::function<SharedClass(const SharedClass&, const SharedClass&)> func){
+ if(nbThreads == 0){
+ return SharedClass();
+ }
+ else if(nbThreads == 1){
+ return objects[0];
+ }
+ else{
+ SharedClass res = func(objects[0], objects[1]);
+ for(int idxThread = 2 ; idxThread < nbThreads ; ++idxThread){
+ res = func(res, objects[idxThread]);
+ }
+ return res;
+ }
+ }
+
+ void set(const SharedClass& val){
+ #pragma omp parallel for num_threads(nbThreads) schedule(static)
+ for(int idxThread = 0 ; idxThread < nbThreads ; ++idxThread){
+ objects[idxThread] = val;
+ }
+ }
+
+ FParObject& operator=(const SharedClass& val){
+ set(val);
+ return *this;
+ }
+};
+
+
+template <class SharedClass>
+class FParArray{
+ const int nbThreads;
+ SharedClass** arrays;
+ size_t arraySize;
+
+ FParArray(const FParArray&) = delete;
+ FParArray& operator=(const FParArray&) = delete;
+
+public:
+
+ explicit FParArray(const size_t inArraySize = 0)
+ : nbThreads(omp_get_max_threads()), arrays(nullptr), arraySize(inArraySize) {
+ arrays = new SharedClass*[nbThreads];
+ #pragma omp parallel for schedule(static)
+ for(int idxThread = 0 ; idxThread < nbThreads ; ++idxThread){
+ arrays[idxThread] = new SharedClass[inArraySize];
+ if(std::is_pod<SharedClass>::value){
+ memset(arrays[idxThread], 0, sizeof(SharedClass)*inArraySize);
+ }
+ }
+ }
+
+ ~FParArray(){
+ for(int idxThread = 0 ; idxThread < nbThreads ; ++idxThread){
+ delete[] arrays[idxThread];
+ }
+ delete[] arrays;
+ }
+
+ void resize(const size_t inArraySize){
+ if(inArraySize != arraySize){
+ for(int idxThread = 0 ; idxThread < nbThreads ; ++idxThread){
+ delete[] arrays[idxThread];
+ }
+
+ arraySize = inArraySize;
+
+ #pragma omp parallel for schedule(static)
+ for(int idxThread = 0 ; idxThread < nbThreads ; ++idxThread){
+ arrays[idxThread] = new SharedClass[inArraySize];
+ if(std::is_pod<SharedClass>::value){
+ memset(arrays[idxThread], 0, sizeof(SharedClass)*inArraySize);
+ }
+ }
+ }
+ }
+
+ SharedClass* getMine(){
+ return arrays[omp_get_thread_num()];
+ }
+
+ const SharedClass* getMine() const {
+ return arrays[omp_get_thread_num()];
+ }
+
+ void applyToAll(std::function<void(SharedClass*)> func){
+ for(int idxThread = 0 ; idxThread < nbThreads ; ++idxThread){
+ func(arrays[idxThread]);
+ }
+ }
+
+ void applyToAllPar(std::function<void(SharedClass*)> func){
+ #pragma omp parallel for num_threads(nbThreads) schedule(static)
+ for(int idxThread = 0 ; idxThread < nbThreads ; ++idxThread){
+ func(arrays[idxThread]);
+ }
+ }
+
+ void reduceArray(SharedClass* resArray, std::function<SharedClass(const SharedClass&, const SharedClass&)> func){
+ if(nbThreads == 0){
+ }
+ else if(nbThreads == 1){
+ if(std::is_pod<SharedClass>::value){
+ memcpy(resArray, arrays[0], sizeof(SharedClass)*arraySize);
+ }
+ else{
+ for(int idxVal = 0 ; idxVal < arraySize ; ++idxVal){
+ resArray = arrays[0][idxVal];
+ }
+ }
+ }
+ else{
+ for(int idxVal = 0 ; idxVal < arraySize ; ++idxVal){
+ resArray[idxVal] = func(arrays[0][idxVal], arrays[1][idxVal]);
+ }
+ for(int idxThread = 2 ; idxThread < nbThreads ; ++idxThread){
+ for(int idxVal = 0 ; idxVal < arraySize ; ++idxVal){
+ resArray[idxVal] = func(resArray[idxVal], arrays[idxThread][idxVal]);
+ }
+ }
+ }
+ }
+
+ void copyArray(const SharedClass* array){
+ #pragma omp parallel for num_threads(nbThreads) schedule(static)
+ for(int idxThread = 0 ; idxThread < nbThreads ; ++idxThread){
+ if(std::is_pod<SharedClass>::value){
+ memcpy(arrays[idxThread], array, sizeof(SharedClass)*arraySize);
+ }
+ else{
+ for(int idxVal = 0 ; idxVal < arraySize ; ++idxVal){
+ arrays[idxThread][idxVal] = array[idxVal];
+ }
+ }
+ }
+ }
+};
+
+#endif // FPAROBJECT_HPP
diff --git a/Tests/Utils/testOctreeParallelFuncteur.cpp b/Tests/Utils/testOctreeParallelFuncteur.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..37e3d3e9c7ce4ba76df7d65ddfb50eb83bec6e69
--- /dev/null
+++ b/Tests/Utils/testOctreeParallelFuncteur.cpp
@@ -0,0 +1,116 @@
+// ===================================================================================
+// 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".
+// ===================================================================================
+
+#include <iostream>
+
+#include <cstdio>
+#include <cstdlib>
+#include <time.h>
+
+#include "../../Src/Utils/FTic.hpp"
+#include "../../Src/Utils/FParameters.hpp"
+
+#include "../../Src/Containers/FOctree.hpp"
+#include "../../Src/Containers/FParForEachOctree.hpp"
+#include "../../Src/Containers/FVector.hpp"
+
+#include "../../Src/Utils/FAssert.hpp"
+#include "../../Src/Utils/FPoint.hpp"
+#include "../../Src/Utils/FParObject.hpp"
+
+#include "../../Src/Components/FBasicParticleContainer.hpp"
+#include "../../Src/Components/FBasicCell.hpp"
+#include "../../Src/Components/FSimpleLeaf.hpp"
+
+#include "../../Src/Files/FRandomLoader.hpp"
+
+#include "../../Src/Utils/FParameterNames.hpp"
+
+/**
+* In this file we show how to use octree's functeur
+*/
+
+int main(int argc, char ** argv){
+ FHelpDescribeAndExit(argc, argv,
+ "Show how to use an octree functeur parallelized (only the code is interesting)",
+ FParameterDefinitions::NbParticles);
+
+ typedef FBasicCell CellClass;
+ typedef FBasicParticleContainer<0> ContainerClass;
+ typedef FSimpleLeaf< ContainerClass > LeafClass;
+ typedef FOctree<CellClass, ContainerClass , LeafClass > OctreeClass;
+ ///////////////////////What we do/////////////////////////////
+ std::cout << ">> This executable is useless to execute.\n";
+ std::cout << ">> It is only interesting to wath the code to understand\n";
+ std::cout << ">> how to use the Octree\n";
+ //////////////////////////////////////////////////////////////
+ const int NbPart = FParameters::getValue(argc,argv,FParameterDefinitions::NbParticles.options, 2000);
+ FTic counter;
+
+ FRandomLoader loader(NbPart, 1, FPoint(0.5,0.5,0.5), 1);
+ OctreeClass tree(10, 3, loader.getBoxWidth(), loader.getCenterOfBox());
+
+ // -----------------------------------------------------
+ std::cout << "Creating and inserting " << NbPart << " particles ..." << std::endl;
+ counter.tic();
+
+ {
+ FPoint particlePosition;
+ for(int idxPart = 0 ; idxPart < loader.getNumberOfParticles() ; ++idxPart){
+ loader.fillParticle(&particlePosition);
+ tree.insert(particlePosition);
+ }
+ }
+
+ counter.tac();
+ std::cout << "Done " << "(" << counter.elapsed() << ")." << std::endl;
+ // -----------------------------------------------------
+
+ // Call a function on each leaf
+ FParObject<long> nbParticles(0);
+ FParForEachOctree::forEachLeaf(&tree, [&](LeafClass* leaf){
+ nbParticles.getMine() += leaf->getSrc()->getNbParticles();
+ });
+ const long totalNbParticles = nbParticles.reduce([](long v1, long v2) -> long {return v1 + v2;});
+ std::cout << "There are " << totalNbParticles << " particles " << std::endl;
+
+ // Call a function on each cell
+ FParObject<long> nbCells;
+ nbCells = 0;
+ FParForEachOctree::forEachCell(&tree, [&nbCells](CellClass* /*cell*/){
+ nbCells.getMine() += 1;
+ });
+ std::cout << "There are " << nbCells.reduce([](long v1, long v2) -> long {return v1 + v2;}) << " cells " << std::endl;
+
+ // To get cell and particles at leaf level
+ FParForEachOctree::forEachCellLeaf(&tree, [&](CellClass* cell, LeafClass* /*leaf*/){
+ cell->resetToInitialState();
+ });
+
+ // If we need an array of long
+ FParArray<long> arrayExample(10);
+ arrayExample.resize(50);
+
+ FParForEachOctree::forEachCellLeaf(&tree, [&](CellClass* /*cell*/, LeafClass* /*leaf*/){
+ arrayExample.getMine()[0] += 10;
+ });
+
+ return 0;
+}
+
+
+
+