diff --git a/Addons/HMat/Src/Clustering/FConnexClustering.hpp b/Addons/HMat/Src/Clustering/FConnexClustering.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..0810206c2a51d2d76c4b87b90984e9df52d1362c
--- /dev/null
+++ b/Addons/HMat/Src/Clustering/FConnexClustering.hpp
@@ -0,0 +1,112 @@
+#ifndef FCONNEXCLUSTERING_HPP
+#define FCONNEXCLUSTERING_HPP
+
+// @SCALFMM_PRIVATE
+
+#include "./Utils/FGlobal.hpp"
+#include "./Utils/FAssert.hpp"
+#include "./Utils/FMath.hpp"
+#include "./Containers/FBoolArray.hpp"
+
+#include "FClusterTree.hpp"
+
+#include <stack>
+#include <vector>
+#include <functional>
+#include <queue>
+#include <limits>
+#include <memory>
+
+
+
+template <class FReal>
+class FConnexClustering {
+protected:
+ const int dim;
+
+ int* permsNewToOrig;
+ int* permsOrigToNew;
+ int* partitions;
+ int* partitionsOffset;
+ int nbPartitions;
+
+public:
+ FConnexClustering(const int inDim, const FReal inDistMat[], const FReal thresh)
+ : dim(inDim),
+ permsNewToOrig(new int[dim]),
+ permsOrigToNew(new int[dim]),
+ partitions(new int[dim]),
+ partitionsOffset(new int[dim+1]),
+ nbPartitions (0){
+
+ std::unique_ptr<int[]> partitionsMappin(new int[dim]);
+
+ for(int idx = 0 ; idx < dim ; ++idx){
+ partitionsMappin[idx] = -1;
+ }
+
+ partitionsOffset[0] = 0;
+ partitions[0] = 0;
+
+ for(int idx = 0 ; idx < dim ; ++idx){
+ if(partitionsMappin[idx] == -1){
+ FAssertLF(nbPartitions < dim);
+
+ partitionsOffset[nbPartitions+1] = partitionsOffset[nbPartitions]+1;
+ partitionsMappin[idx] = nbPartitions;
+
+ int idxPartitionElement = partitionsOffset[nbPartitions];
+ permsNewToOrig[idxPartitionElement] = idx;
+
+ while(idxPartitionElement < partitionsOffset[nbPartitions+1]){
+ FAssertLF(idxPartitionElement < dim);
+
+ for(int idxOther = 0 ; idxOther < dim ; ++idxOther){
+ if(partitionsMappin[idxOther] == -1
+ && inDistMat[permsNewToOrig[idxPartitionElement]*dim + idxOther] < thresh){
+ partitionsMappin[idxOther] = nbPartitions;
+ permsNewToOrig[partitionsOffset[nbPartitions+1]] = idxOther;
+ permsOrigToNew[idxOther] = partitionsOffset[nbPartitions+1];
+ partitionsOffset[nbPartitions+1] += 1;
+ FAssertLF(partitionsOffset[nbPartitions+1] <= dim);
+ }
+ }
+
+ idxPartitionElement += 1;
+ }
+
+ partitions[nbPartitions] = partitionsOffset[nbPartitions+1]-partitionsOffset[nbPartitions];
+ nbPartitions += 1;
+ }
+ }
+
+ FAssertLF(partitionsOffset[nbPartitions] == dim);
+ }
+
+ ~FConnexClustering(){
+ delete[] permsNewToOrig;
+ delete[] permsOrigToNew;
+ delete[] partitionsOffset;
+ delete[] partitions;
+ }
+
+
+ void fillPermutations(int* inPermuts, int* invPermuts = nullptr) const {
+ memcpy(inPermuts, permsOrigToNew, sizeof(int)*dim);
+ if(invPermuts){
+ memcpy(invPermuts, permsNewToOrig, sizeof(int)*dim);
+ }
+ }
+
+ int getNbPartitions() const{
+ return nbPartitions;
+ }
+
+ void getPartitions(const int inNbPartitions, int inNbIdxInPartitions[]) const{
+ FAssertLF(nbPartitions == inNbPartitions);
+ memcpy(inNbIdxInPartitions, partitions, sizeof(int)*nbPartitions);
+ }
+};
+
+#endif // FCONNEXCLUSTERING_HPP
+
diff --git a/Addons/HMat/Tests/testConnectClustering.cpp b/Addons/HMat/Tests/testConnectClustering.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..c888cf3a882534fdaf4022bf6da3ec4fd022d050
--- /dev/null
+++ b/Addons/HMat/Tests/testConnectClustering.cpp
@@ -0,0 +1,158 @@
+
+// @SCALFMM_PRIVATE
+
+#include "../Src/Clustering/FConnexClustering.hpp"
+#include "../Src/Utils/FMatrixIO.hpp"
+
+#include "../Src/Containers/FPartitionsMapping.hpp"
+#include "../Src/Utils/FSvgRect.hpp"
+#include "../Src/Viewers/FDenseBlockWrapper.hpp"
+#include "../Src/Blocks/FDenseBlock.hpp"
+
+#include "Utils/FParameters.hpp"
+#include "Utils/FParameterNames.hpp"
+
+#include "Files/FGenerateDistribution.hpp"
+
+#include <memory>
+
+// FUSE_SCOTCH
+
+template <class FReal, class PartitionerClass>
+void TestPartitions(const PartitionerClass& partitioner,
+ const char outputdir[], const int idFile,
+ const int dim,
+ const FReal coords[], const char config[]){
+
+ std::unique_ptr<int[]> permutations(new int[dim]);
+ std::unique_ptr<int[]> invpermutations(new int[dim]);
+ partitioner.fillPermutations(permutations.get(), invpermutations.get());
+
+ const int nbPartitions = partitioner.getNbPartitions();
+ std::unique_ptr<int[]> partitions(new int[nbPartitions]);
+ partitioner.getPartitions(nbPartitions, partitions.get());
+
+ {
+ char coordfilename[1024];
+ sprintf(coordfilename, "%s/%s-coord-%d.csv", outputdir, config, idFile);
+ std::cout << "\t save coord to " << coordfilename << "\n";
+ FILE* fcoord = fopen(coordfilename, "w");
+
+ int offsetParticles = 0;
+ for(int idxPartition = 0 ; idxPartition < nbPartitions ; ++idxPartition){
+ for(int idxPart = 0 ; idxPart < partitions[idxPartition] ; ++idxPart){
+ const int idxUnk = invpermutations[idxPart+offsetParticles];
+ fprintf(fcoord, "%e,%e,%e,%d\n",
+ coords[idxUnk*4 + 0],
+ coords[idxUnk*4 + 1],
+ coords[idxUnk*4 + 2],
+ idxPartition);
+ }
+ offsetParticles += partitions[idxPartition];
+ }
+ FAssertLF(offsetParticles == dim);
+
+ fclose(fcoord);
+ }
+ {
+ typedef FDenseBlock<FReal> CellClass;
+ typedef FPartitionsMapping<FReal, CellClass> GridClass;
+
+ GridClass bissection(dim, partitions.get(), nbPartitions);
+
+ char svgfilename[1024];
+ sprintf(svgfilename, "%s/%s-%d.svg", outputdir, config, idFile);
+ std::cout << "\t save svg to " << svgfilename << "\n";
+ FSvgRect output(svgfilename, dim);
+
+ bissection.forAllBlocksDescriptor([&](const FBlockDescriptor& info){
+ output.addRectWithLegend(info.col, info.row, info.nbCols, info.nbRows, info.level);
+ });
+ }
+}
+
+int main(int argc, char** argv){
+ static const FParameterNames SvgOutParam = {
+ {"-fout", "--out", "-out"} ,
+ "Svg output directory."
+ };
+ static const FParameterNames DimParam = {
+ {"-N", "-nb", "-dim"} ,
+ "Dim of the matrix."
+ };
+
+ FHelpDescribeAndExit(argc, argv,"Test the bisection.",SvgOutParam,DimParam,FParameterDefinitions::OctreeHeight,
+ FParameterDefinitions::NbParticles);
+
+ const int nbParticles = (FParameters::getValue(argc, argv, FParameterDefinitions::NbParticles.options, 3000) & ~1);
+ const char* outputdir = FParameters::getStr(argc, argv, SvgOutParam.options, "/tmp/");
+
+ typedef double FReal;
+
+ /////////////////////////////////////////////////////////////////
+
+ std::cout << "Create spheres\n";
+
+ const FReal radius = 0.5;
+ const FReal distanceBetweenSpheres = 0.4;
+ const int nbPointsSphere = nbParticles/2;
+
+ std::cout << "\tradius = " << radius << "\n";
+ std::cout << "\tdistanceBetweenSpheres = " << distanceBetweenSpheres << "\n";
+ std::cout << "\tnbPointsSphere = " << nbPointsSphere << "\n";
+
+ std::unique_ptr<FReal[]> spherePoints(new FReal[nbParticles*4]);
+ unifRandonPointsOnSphere(nbPointsSphere, radius, spherePoints.get()) ;
+
+ for(int idxPoint = 0 ; idxPoint < nbPointsSphere ; ++idxPoint){
+ spherePoints[(idxPoint+nbPointsSphere)*4 + 0] = spherePoints[idxPoint*4 + 0] + radius*2 + distanceBetweenSpheres;
+ spherePoints[(idxPoint+nbPointsSphere)*4 + 1] = spherePoints[idxPoint*4 + 1] + radius*2 + distanceBetweenSpheres;
+ spherePoints[(idxPoint+nbPointsSphere)*4 + 2] = spherePoints[idxPoint*4 + 2] + radius*2 + distanceBetweenSpheres;
+ }
+
+ /////////////////////////////////////////////////////////////////
+
+ std::cout << "Compute distance\n";
+
+ const int dim = nbParticles;
+ std::unique_ptr<FReal[]> distances(new FReal[nbParticles*nbParticles]);
+ for(int idxCol = 0 ; idxCol < nbParticles ; ++idxCol){
+ for(int idxRow = 0 ; idxRow < nbParticles ; ++idxRow){
+ const FReal diffx = spherePoints[idxCol*4 + 0] - spherePoints[idxRow*4 + 0];
+ const FReal diffy = spherePoints[idxCol*4 + 1] - spherePoints[idxRow*4 + 1];
+ const FReal diffz = spherePoints[idxCol*4 + 2] - spherePoints[idxRow*4 + 2];
+ distances[idxCol*nbParticles+idxRow] = FMath::Sqrt((diffx*diffx) + (diffy*diffy) + (diffz*diffz));
+ }
+ }
+
+ /////////////////////////////////////////////////////////////////
+
+ FReal distMin = std::numeric_limits<FReal>::max();
+ FReal distMax = std::numeric_limits<FReal>::min();
+
+ for(int idxRow = 0 ; idxRow < dim ; ++idxRow){
+ for(int idxCol = idxRow ; idxCol < dim ; ++idxCol){
+ distMin = FMath::Min(distMin, distances[idxCol*nbParticles+idxRow]);
+ distMax = FMath::Max(distMax, distances[idxCol*nbParticles+idxRow]);
+ }
+ }
+ std::cout << "Dist min " << distMin << "\n";
+ std::cout << "Dist max " << distMax << "\n";
+
+ /////////////////////////////////////////////////////////////////
+ const FReal stepThresh = (distMax-distMin)/10;
+ for(FReal thresh = distMin ; thresh <= distMax ; thresh += stepThresh){
+ std::cout << "Test FConnexClustering for thresh " << thresh << "\n";
+ FConnexClustering<FReal> partitioner(dim, distances.get(), thresh);
+ std::cout << "\tGot " << partitioner.getNbPartitions() << " partitions\n";
+ TestPartitions<FReal, FConnexClustering<FReal>>(partitioner,outputdir,
+ int((thresh-distMin)/stepThresh),
+ dim, spherePoints.get(), "FConnexClustering");
+ }
+
+ return 0;
+}
+
+
+
+