yade-dev team mailing list archive

[Janek Kozicki <janek_listy@xxxxx>]

Bruno Chareyre said:     (by the date of Wed, 29 Nov 2017 16:24:01 +0100)

> 
> 
> On 11/29/2017 03:25 PM, Janek Kozicki wrote:
> > I hope that you have, Bruno, among those tests some related to PFV
> > and CGAL :)
> Hi Janek,
> Yes we have, thanks. :)
> DEM-PFV-check.py depends on pfv+cgal+TesselationWrapper+openblas+cholmod 
> (at least...). You can use it as a test if it helps.
> 
> I finished the core part of the cgal/alpha thing, I'm now a bit more 
> free to look at your patch. Did you try to compile it with older 
> ("current" in my situation) cgal versions yet?
> Is there a residual amount of version-specific code or did you see a way 
> to have one unique code compatible with both cgal versions.

Yes, a single code is compatible with both CGAL versions. I removed
almost all #ifdefs. It has compiled successfully both with CGAL 4.9 and 4.11.

However there is some problem, which I didn't foresee. First I wanted
to wait until yade --check with CGAL 4.11 would finish successfully.
And then run the CGAL 4.9 compilation check. I was doing other stuff,
and I did let yade --check run for 4 hours or maybe more. It is stuck
on checkPolyhedraCrush.py

Initially I wanted to continue cleaning this diff after it would
finish this yade --check test.

But then you have wrote this email, so I have ("quickly" ;) compiled
yade with the same patch and CGAL 4.9 successfully, and to my
surprise yade --check did finish very quickly. While the 4.11 still
didn't complete.

So there is a problem here.

I wanted to clean this patch a bit more - because I didn't check yet
if those #ifdefs are necessary in RegularTriangulation.h, probably
not, I just don't know yet. But I have removed it everywhere else
though.

But since we have this other problem I see that removing #ifdefs from
RegularTriangulation.h can wait a bit. And I am sending you my
today's version of this patch.

If you prefer to move this work to github, just tell me in which
branch we are going to do this work, and I will join you there.

best regards
Janek

diff --git a/lib/triangulation/FlowBoundingSphere.ipp b/lib/triangulation/FlowBoundingSphere.ipp
index 86d21f3..56c7585 100644
--- a/lib/triangulation/FlowBoundingSphere.ipp
+++ b/lib/triangulation/FlowBoundingSphere.ipp
@@ -110,7 +110,7 @@ void FlowBoundingSphere<Tesselation>::averageRelativeCellVelocity()
 				CVector Surfk = cell->info()-cell->neighbor(i)->info();
 				Real area = sqrt ( Surfk.squared_length() );
 				Surfk = Surfk/area;
-                        	CVector branch = cell->vertex ( facetVertices[i][0] )->point() - cell->info();
+                        	CVector branch = cell->vertex ( facetVertices[i][0] )->point().point() - cell->info();
                         	posAvFacet = (Point) cell->info() + ( branch*Surfk ) *Surfk;
 				facetFlowRate = (cell->info().kNorm())[i] * (cell->info().shiftedP() - cell->neighbor (i)->info().shiftedP());
 				totFlowRate += facetFlowRate;
@@ -215,7 +215,7 @@ double FlowBoundingSphere<Tesselation>::getPorePressure (double X, double Y, dou
 {
 	if (noCache && T[!currentTes].Max_id()<=0) return 0;//the engine never solved anything
 	RTriangulation& Tri = T[noCache?(!currentTes):currentTes].Triangulation();
-	CellHandle cell = Tri.locate(Point(X,Y,Z));
+	CellHandle cell = Tri.locate(CGT::Sphere(X,Y,Z));
 	return cell->info().p();
 }
 
@@ -224,7 +224,7 @@ int FlowBoundingSphere<Tesselation>::getCell (double X, double Y, double Z)
 {
 	if (noCache && T[!currentTes].Max_id()<=0) {cout<<"Triangulation does not exist. Sorry."<<endl; return -1;}
 	RTriangulation& Tri = T[noCache?(!currentTes):currentTes].Triangulation();
-	CellHandle cell = Tri.locate(Point(X,Y,Z));
+	CellHandle cell = Tri.locate(CGT::Sphere(X,Y,Z));
 	return cell->info().id;
 }
 
@@ -245,7 +245,7 @@ void FlowBoundingSphere<Tesselation>::measurePressureProfile(double WallUpy, dou
 	int cell=0;
 	for (int i=0; i<captures; i++){
         for (double Z=min(zMin,zMax); Z<=max(zMin,zMax); Z+=std::abs(Rz)) {
-		permeameter = Tri.locate(Point(X, Y, Z));
+		permeameter = Tri.locate(CGT::Sphere(X, Y, Z));
 		pressure+=permeameter->info().p();
 		cell++;
         }
@@ -264,7 +264,7 @@ double FlowBoundingSphere<Tesselation>::averageSlicePressure(double Y)
   double Rz = (zMax-zMin)/30;
   for (double X=xMin; X<=xMax+Ry/10; X=X+Rx) {
 	for (double Z=zMin; Z<=zMax+Ry/10; Z=Z+Rz) {
-	  P_ave+=Tri.locate(Point(X, Y, Z))->info().p();
+	  P_ave+=Tri.locate(CGT::Sphere(X, Y, Z))->info().p();
 	  n++;
 	}
   }
@@ -433,7 +433,7 @@ template <class Tesselation>
 CVector FlowBoundingSphere<Tesselation>::cellBarycenter(CellHandle& cell)
 {
 	CVector center ( 0,0,0 );
-	for ( int k=0;k<4;k++ ) center= center + 0.25* (cell->vertex(k)->point()-CGAL::ORIGIN);
+	for ( int k=0;k<4;k++ ) center= center + 0.25* (cell->vertex(k)->point().point()-CGAL::ORIGIN);
 	return center;
 }
 
@@ -446,17 +446,17 @@ void FlowBoundingSphere<Tesselation>::interpolate(Tesselation& Tes, Tesselation&
 		CellHandle& newCell = *cellIt;
 		if (newCell->info().Pcondition || newCell->info().isGhost) continue;
 		CVector center ( 0,0,0 );
-		if (newCell->info().fictious()==0) for ( int k=0;k<4;k++ ) center= center + 0.25* (Tes.vertex(newCell->vertex(k)->info().id())->point()-CGAL::ORIGIN);
+		if (newCell->info().fictious()==0) for ( int k=0;k<4;k++ ) center= center + 0.25* (Tes.vertex(newCell->vertex(k)->info().id())->point().point()-CGAL::ORIGIN);
 		else {
 			Real boundPos=0; int coord=0;
-			for ( int k=0;k<4;k++ ) if (!newCell->vertex (k)->info().isFictious) center= center+(1./(4.-newCell->info().fictious()))*(Tes.vertex(newCell->vertex(k)->info().id())->point()-CGAL::ORIGIN);
+			for ( int k=0;k<4;k++ ) if (!newCell->vertex (k)->info().isFictious) center= center+(1./(4.-newCell->info().fictious()))*(Tes.vertex(newCell->vertex(k)->info().id())->point().point()-CGAL::ORIGIN);
 			for ( int k=0;k<4;k++ ) if (newCell->vertex (k)->info().isFictious) {
 					coord=boundary (newCell->vertex(k)->info().id()).coordinate;
 					boundPos=boundary (newCell->vertex(k)->info().id()).p[coord];
 					center=CVector(coord==0?boundPos:center[0],coord==1?boundPos:center[1],coord==2?boundPos:center[2]);
 				}
 		}
-                oldCell = Tri.locate(Point(center[0],center[1],center[2]));
+                oldCell = Tri.locate(CGT::Sphere(center[0],center[1],center[2]));
 		newCell->info().getInfo(oldCell->info());
 //                 newCell->info().p() = oldCell->info().shiftedP();
         }
@@ -467,12 +467,12 @@ template <class Tesselation>
 Real FlowBoundingSphere<Tesselation>::checkSphereFacetOverlap(const Sphere& v0, const Sphere& v1, const Sphere& v2)
 {
 	//First, check that v0 projection fall between v1 and v2...
-	Real dist=(v0-v1)*(v2-v1);
+	Real dist=(v0.point()-v1.point())*(v2.point()-v1.point());
 	if (dist<0) return 0;
-	Real v1v2=(v2-v1).squared_length();
+	Real v1v2=(v2.point()-v1.point()).squared_length();
 	if (dist>v1v2) return 0;
 	//... then, check distance
-	Real m=(cross_product(v0-v1,v2-v1)).squared_length()/v1v2;
+	Real m=(cross_product(v0.point()-v1.point(),v2.point()-v1.point())).squared_length()/v1v2;
 	if (m<v0.weight()) {
 		Real d=2*sqrt((v0.weight()-m));
 		Real teta=2*acos(sqrt(m/v0.weight()));
@@ -530,9 +530,9 @@ void FlowBoundingSphere<Tesselation>::computePermeability()
 				Sphere& v1 = W[1]->point();
 				Sphere& v2 = W[2]->point();
 				cell->info().facetSphereCrossSections[j]=CVector(
-				   W[0]->info().isFictious ? 0 : 0.5*v0.weight()*acos((v1-v0)*(v2-v0)/sqrt((v1-v0).squared_length()*(v2-v0).squared_length())),
-				   W[1]->info().isFictious ? 0 : 0.5*v1.weight()*acos((v0-v1)*(v2-v1)/sqrt((v1-v0).squared_length()*(v2-v1).squared_length())),
-				   W[2]->info().isFictious ? 0 : 0.5*v2.weight()*acos((v0-v2)*(v1-v2)/sqrt((v1-v2).squared_length()*(v2-v0).squared_length())));
+				   W[0]->info().isFictious ? 0 : 0.5*v0.weight()*acos((v1.point()-v0.point())*(v2.point()-v0.point())/sqrt((v1.point()-v0.point()).squared_length()*(v2.point()-v0.point()).squared_length())),
+				   W[1]->info().isFictious ? 0 : 0.5*v1.weight()*acos((v0.point()-v1.point())*(v2.point()-v1.point())/sqrt((v1.point()-v0.point()).squared_length()*(v2.point()-v1.point()).squared_length())),
+				   W[2]->info().isFictious ? 0 : 0.5*v2.weight()*acos((v0.point()-v2.point())*(v1.point()-v2.point())/sqrt((v1.point()-v2.point()).squared_length()*(v2.point()-v0.point()).squared_length())));
 				//FIXME: it should be possible to skip completely blocked cells, currently the problem is it segfault for undefined areas
 // 				if (cell->info().blocked) continue;//We don't need permeability for blocked cells, it will be set to zero anyway
 				pass+=1;
@@ -785,7 +785,7 @@ void FlowBoundingSphere<Tesselation>::initializePressure( double pZero )
         
         IPCells.clear();
         for (unsigned int n=0; n<imposedP.size();n++) {
-		CellHandle cell=Tri.locate(imposedP[n].first);
+		CellHandle cell=Tri.locate(CGT::Sphere(imposedP[n].first,0));
 		//check redundancy
 		for (unsigned int kk=0;kk<IPCells.size();kk++){
 			if (cell==IPCells[kk]) cerr<<"Two imposed pressures fall in the same cell."<<endl;
@@ -797,7 +797,7 @@ void FlowBoundingSphere<Tesselation>::initializePressure( double pZero )
 
 	IFCells.clear();
 	for (unsigned int n=0; n<imposedF.size();n++) {
-		CellHandle cell=Tri.locate(imposedF[n].first);
+		CellHandle cell=Tri.locate(CGT::Sphere(imposedF[n].first,0));
 		//check redundancy
 		for (unsigned int kk=0;kk<IPCells.size();kk++){
 			if (cell==IPCells[kk]) cerr<<"Both flux and pressure are imposed in the same cell."<<endl;
diff --git a/lib/triangulation/FlowBoundingSphereLinSolv.ipp b/lib/triangulation/FlowBoundingSphereLinSolv.ipp
index ed88603..6eb78e8 100644
--- a/lib/triangulation/FlowBoundingSphereLinSolv.ipp
+++ b/lib/triangulation/FlowBoundingSphereLinSolv.ipp
@@ -9,7 +9,9 @@
 
 // #define XVIEW
 #include "FlowBoundingSphereLinSolv.hpp"//include after #define XVIEW
-#include "CGAL/constructions/constructions_on_weighted_points_cartesian_3.h"
+#if CGAL_VERSION_NR < CGAL_VERSION_NUMBER(4,11,0)
+	#include "CGAL/constructions/constructions_on_weighted_points_cartesian_3.h"
+#endif
 #include <CGAL/Width_3.h>
 #include <iostream>
 #include <fstream>
diff --git a/lib/triangulation/KinematicLocalisationAnalyser.cpp b/lib/triangulation/KinematicLocalisationAnalyser.cpp
index f4e2189..bc74c54 100644
--- a/lib/triangulation/KinematicLocalisationAnalyser.cpp
+++ b/lib/triangulation/KinematicLocalisationAnalyser.cpp
@@ -447,7 +447,7 @@ NormalDisplacementDistribution(vector<Edge_iterator>& edges, vector<pair<Real,Re
 			ed_it!=ed_end; ++ed_it) {
 		Vh1= (*ed_it)->first->vertex((*ed_it)->second);
 		Vh2= (*ed_it)->first->vertex((*ed_it)->third);
-		branch = Vh1->point()- Vh2->point();
+		branch = Vh1->point().point()- Vh2->point().point();
 		NORMALIZE(branch);
 		if (consecutive)
 			U = TS1->grain(Vh1->info().id()).translation -
@@ -741,10 +741,11 @@ CVector KinematicLocalisationAnalyser::Deplacement(Finite_cells_iterator cell, i
 
 void KinematicLocalisationAnalyser::Grad_u(Finite_cells_iterator cell, int facet, CVector &V, Tenseur3& T)
 {
-	CVector S = cross_product((cell->vertex(l_vertices[facet][1])->point())
-							  - (cell->vertex(l_vertices[facet][0])->point()),
-							  (cell->vertex(l_vertices[facet][2])->point()) -
-							  (cell->vertex(l_vertices[facet][1])->point())) /2.f;
+	CVector S = cross_product((cell->vertex(l_vertices[facet][1])->point().point())
+							  - (cell->vertex(l_vertices[facet][0])->point().point()),
+							  (cell->vertex(l_vertices[facet][2])->point().point()) -
+							  (cell->vertex(l_vertices[facet][1])->point().point())) /2.f;
+
 	Somme(T, V, S);
 }
 
diff --git a/lib/triangulation/Network.ipp b/lib/triangulation/Network.ipp
index fb26b3a..a1e67a0 100644
--- a/lib/triangulation/Network.ipp
+++ b/lib/triangulation/Network.ipp
@@ -1,7 +1,9 @@
 
 #ifdef FLOW_ENGINE
 
-#include "CGAL/constructions/constructions_on_weighted_points_cartesian_3.h"
+#if CGAL_VERSION_NR < CGAL_VERSION_NUMBER(4,11,0)
+	#include "CGAL/constructions/constructions_on_weighted_points_cartesian_3.h"
+#endif
 #include <iostream>
 #include <fstream>
 #include <new>
@@ -66,7 +68,7 @@ double Network<Tesselation>::volumePoreVoronoiFraction (CellHandle& cell, int& j
                 VertexHandle& SV2 = W[1];
                 VertexHandle& SV3 = W[2];
 
-                cell->info().facetSurfaces[j]=0.5*CGAL::cross_product(SV1->point()-SV3->point(),SV2->point()-SV3->point());
+                cell->info().facetSurfaces[j]=0.5*CGAL::cross_product(SV1->point().point()-SV3->point().point(),SV2->point().point()-SV3->point().point());
 		if (cell->info().facetSurfaces[j][0]==0 && cell->info().facetSurfaces[j][1]==0 && cell->info().facetSurfaces[j][2]==0) cerr<<"NULL FACET SURF"<<endl;
                 if (cell->info().facetSurfaces[j]*(p2-p1) > 0) cell->info().facetSurfaces[j] = -1.0*cell->info().facetSurfaces[j];
                 Real Vtot = abs(ONE_THIRD*cell->info().facetSurfaces[j]*(p1-p2));
@@ -74,8 +76,8 @@ double Network<Tesselation>::volumePoreVoronoiFraction (CellHandle& cell, int& j
 		
                 double Vsolid1=0, Vsolid2=0;
                 for (int i=0;i<3;i++) {
-                Vsolid1 += sphericalTriangleVolume(v[permut3[i][0]],v[permut3[i][1]],p1,p2);
-                Vsolid2 += sphericalTriangleVolume(v[permut3[i][0]],v[permut3[i][2]],p1,p2);}
+                Vsolid1 += sphericalTriangleVolume(v[permut3[i][0]],v[permut3[i][1]].point(),p1,p2);
+                Vsolid2 += sphericalTriangleVolume(v[permut3[i][0]],v[permut3[i][2]].point(),p1,p2);}
 
 		VSolidTot += Vsolid1 + Vsolid2;
 		vPoral += Vtot - (Vsolid1 + Vsolid2);
@@ -98,7 +100,7 @@ double Network<Tesselation>::volumeSolidPore (const CellHandle& cell)
 {
   double Vsolid=0;
   for (int i=0;i<4;i++) {
-	if ( !cell->vertex(permut4[i][0])->info().isFictious ) Vsolid += sphericalTriangleVolume( cell->vertex(permut4[i][0])->point(), cell->vertex(permut4[i][1])->point(), cell->vertex(permut4[i][2])-> point(), cell->vertex(permut4[i][3])-> point());
+	if ( !cell->vertex(permut4[i][0])->info().isFictious ) Vsolid += sphericalTriangleVolume( cell->vertex(permut4[i][0])->point(), cell->vertex(permut4[i][1])->point().point(), cell->vertex(permut4[i][2])-> point().point(), cell->vertex(permut4[i][3])-> point().point());
   }
   return Vsolid;
 }
@@ -128,8 +130,8 @@ double Network<Tesselation>::volumeSingleFictiousPore(const VertexHandle& SV1, c
         Sphere& SW2 = SV2->point();
         Sphere& SW3 = SV3->point();
 
-        Real Vsolid1 = sphericalTriangleVolume(SW2, AA, PV1, PV2)+sphericalTriangleVolume(SW2, SW3, PV1, PV2);
-        Real Vsolid2 = sphericalTriangleVolume(SW3, BB, PV1, PV2)+sphericalTriangleVolume(SW3, SW2, PV1, PV2);
+        Real Vsolid1 = sphericalTriangleVolume(SW2, AA, PV1, PV2)+sphericalTriangleVolume(SW2, SW3.point(), PV1, PV2);
+        Real Vsolid2 = sphericalTriangleVolume(SW3, BB, PV1, PV2)+sphericalTriangleVolume(SW3, SW2.point(), PV1, PV2);
 	
 	VSolidTot += Vsolid1 + Vsolid2;
 	vPoral += Vtot - (Vsolid1 + Vsolid2);
@@ -151,7 +153,7 @@ double Network<Tesselation>::volumeDoubleFictiousPore(const VertexHandle& SV1, c
         Point AA(A[0],A[1],A[2]);
         Point BB(B[0],B[1],B[2]);
 
-        facetSurface = CGAL::cross_product(SV3->point()-AA,SV3->point()-BB);
+        facetSurface = CGAL::cross_product(SV3->point().point()-AA,SV3->point().point()-BB);
         if (facetSurface*(PV2-PV1) > 0) facetSurface = -1.0*facetSurface;
         Real Vtot = abs(facetSurface*(PV1-PV2))*ONE_THIRD;
 	Vtotalissimo += Vtot;
@@ -179,7 +181,7 @@ double Network<Tesselation>::fastSphericalTriangleArea(const Sphere& STA1, const
         using namespace CGAL;
         double rayon2 = STA1.weight();
         if (rayon2 == 0.0) return 0.0;
-        return rayon2 * fastSolidAngle(STA1,STA2,STA3,PTA1);
+        return rayon2 * fastSolidAngle(STA1.point(),STA2,STA3,PTA1);
 }
 
 template<class Tesselation>
@@ -273,14 +275,14 @@ double Network<Tesselation>::surfaceSolidThroat(CellHandle cell, int j, bool sli
                 VertexHandle& SV2 = W[1];
                 VertexHandle& SV3 = W[2];
 
-		Ssolid1 = fastSphericalTriangleArea(SV1->point(), SV2->point(), p1, p2);
-                Ssolid1n = fastSphericalTriangleArea(SV1->point(), SV3->point(), p1, p2);
+		Ssolid1 = fastSphericalTriangleArea(SV1->point(), SV2->point().point(), p1, p2);
+                Ssolid1n = fastSphericalTriangleArea(SV1->point(), SV3->point().point(), p1, p2);
                 cell->info().solidSurfaces[j][0]=Ssolid1+Ssolid1n;
-                Ssolid2 = fastSphericalTriangleArea(SV2->point(),SV1->point(),p1, p2);
-                Ssolid2n = fastSphericalTriangleArea(SV2->point(),SV3->point(),p1, p2);
+                Ssolid2 = fastSphericalTriangleArea(SV2->point(),SV1->point().point(),p1, p2);
+                Ssolid2n = fastSphericalTriangleArea(SV2->point(),SV3->point().point(),p1, p2);
                 cell->info().solidSurfaces[j][1]=Ssolid2+Ssolid2n;
-                Ssolid3 = fastSphericalTriangleArea(SV3->point(),SV2->point(),p1, p2);
-                Ssolid3n = fastSphericalTriangleArea(SV3->point(),SV1->point(),p1, p2);
+                Ssolid3 = fastSphericalTriangleArea(SV3->point(),SV2->point().point(),p1, p2);
+                Ssolid3n = fastSphericalTriangleArea(SV3->point(),SV1->point().point(),p1, p2);
                 cell->info().solidSurfaces[j][2]=Ssolid3+Ssolid3n;
 
     }; break;
@@ -292,14 +294,14 @@ double Network<Tesselation>::surfaceSolidThroat(CellHandle cell, int j, bool sli
 		Boundary &bi1 =  boundary(SV1->info().id());
                 Ssolid1 = 0;
 		if (bi1.flowCondition && ! slipBoundary) {
-                        Ssolid1 = abs(0.5*CGAL::cross_product(p1-p2, SV2->point()-SV3->point())[bi1.coordinate]);
+                        Ssolid1 = abs(0.5*CGAL::cross_product(p1-p2, SV2->point().point()-SV3->point().point())[bi1.coordinate]);
                         cell->info().solidSurfaces[j][facetF1]=Ssolid1;
                 }
-                Ssolid2 = fastSphericalTriangleArea(SV2->point(),SV1->point(),p1, p2);
-                Ssolid2n = fastSphericalTriangleArea(SV2->point(),SV3->point(),p1, p2);
+                Ssolid2 = fastSphericalTriangleArea(SV2->point(),SV1->point().point(),p1, p2);
+                Ssolid2n = fastSphericalTriangleArea(SV2->point(),SV3->point().point(),p1, p2);
                 cell->info().solidSurfaces[j][facetRe1]=Ssolid2+Ssolid2n;
-                Ssolid3 = fastSphericalTriangleArea(SV3->point(),SV2->point(),p1, p2);
-                Ssolid3n = fastSphericalTriangleArea(SV3->point(),SV1->point(),p1, p2);
+                Ssolid3 = fastSphericalTriangleArea(SV3->point(),SV2->point().point(),p1, p2);
+                Ssolid3n = fastSphericalTriangleArea(SV3->point(),SV1->point().point(),p1, p2);
                 cell->info().solidSurfaces[j][facetRe2]=Ssolid3+Ssolid3n;
     }; break;
     case (2) : {
@@ -331,7 +333,7 @@ double Network<Tesselation>::surfaceSolidThroat(CellHandle cell, int j, bool sli
                 Ssolid1n = fastSphericalTriangleArea(SV3->point(), BB, p1, p2);
                 cell->info().solidSurfaces[j][facetRe1]=Ssolid1+Ssolid1n;
                 //area vector of triangle (p1,sphere,p2)
-                CVector p1p2v1Surface = 0.5*CGAL::cross_product(p1-p2,SV3->point()-p2);
+                CVector p1p2v1Surface = 0.5*CGAL::cross_product(p1-p2,SV3->point().point()-p2);
                 if (bi1.flowCondition && ! slipBoundary) {
                         //projection on boundary 1
                         Ssolid2 = abs(p1p2v1Surface[bi1.coordinate]);
@@ -377,9 +379,9 @@ double Network<Tesselation>::surfaceSolidThroatInPore(CellHandle cell, int j, bo
                 VertexHandle& SV2 = W[1];
                 VertexHandle& SV3 = W[2];
 
-		Ssolid1 = fastSphericalTriangleArea(SV1->point(), SV2->point(), p1, SV3->point());
-                Ssolid2 = fastSphericalTriangleArea(SV2->point(),SV1->point(),p1, SV3->point());
-                Ssolid3 = fastSphericalTriangleArea(SV3->point(),SV2->point(),p1, SV1->point());
+		Ssolid1 = fastSphericalTriangleArea(SV1->point(),SV2->point().point(),p1, SV3->point().point());
+                Ssolid2 = fastSphericalTriangleArea(SV2->point(),SV1->point().point(),p1, SV3->point().point());
+                Ssolid3 = fastSphericalTriangleArea(SV3->point(),SV2->point().point(),p1, SV1->point().point());
     }; break;
     case (1) : {
 		VertexHandle SV1 = cell->vertex(facetVertices[j][facetF1]);
@@ -388,9 +390,9 @@ double Network<Tesselation>::surfaceSolidThroatInPore(CellHandle cell, int j, bo
 
 		Boundary &bi1 =  boundary(SV1->info().id());
                 Ssolid1 = 0;
-		if (bi1.flowCondition && ! slipBoundary) Ssolid1 = abs(0.5*CGAL::cross_product(p1-SV2->point(), SV2->point()-SV3->point())[bi1.coordinate]);
-                Ssolid2 = fastSphericalTriangleArea(SV2->point(),SV3->point(),p1, SV2->point()+bi1.normal);
-                Ssolid3 = fastSphericalTriangleArea(SV3->point(),SV2->point(),p1, SV3->point()+bi1.normal);
+		if (bi1.flowCondition && ! slipBoundary) Ssolid1 = abs(0.5*CGAL::cross_product(p1-SV2->point().point(), SV2->point().point()-SV3->point().point())[bi1.coordinate]);
+                Ssolid2 = fastSphericalTriangleArea(SV2->point(),SV3->point().point(),p1, SV2->point().point()+bi1.normal);
+                Ssolid3 = fastSphericalTriangleArea(SV3->point(),SV2->point().point(),p1, SV3->point().point()+bi1.normal);
     }; break;
     case (2) : {
 		double A [3], B[3], C[3];
@@ -413,7 +415,7 @@ double Network<Tesselation>::surfaceSolidThroatInPore(CellHandle cell, int j, bo
                 Sphere C1(CC, 0);
                 //FIXME : we are computing triangle area twice here, because its computed in volume_double_fictious already -> optimize
                 Ssolid1 = 0.5*(fastSphericalTriangleArea(SV3->point(), AA, p1, p2)+ fastSphericalTriangleArea(SV3->point(), BB, p1, p2));
-                CVector p1p2v1Surface = 0.5*CGAL::cross_product(p1-p2,SV3->point()-p2);
+                CVector p1p2v1Surface = 0.5*CGAL::cross_product(p1-p2,SV3->point().point()-p2);
                 if (bi1.flowCondition && ! slipBoundary) Ssolid2 = 0.5*abs(p1p2v1Surface[bi1.coordinate]);
                 if (bi2.flowCondition && ! slipBoundary) Ssolid3 = 0.5*abs(p1p2v1Surface[bi2.coordinate]); 
     }; break;
@@ -442,7 +444,7 @@ CVector Network<Tesselation>::surfaceSingleFictiousFacet(VertexHandle fSV1, Vert
 //  const Boundary &bi2 = boundary ( fSV2->info().id() );
         CVector mean_height = (bi1.p[bi1.coordinate]-0.5*(SV3->point()[bi1.coordinate]+SV2->point()[bi1.coordinate]))*bi1.normal;
 
-        return CGAL::cross_product(mean_height,SV3->point()-SV2->point());
+        return CGAL::cross_product(mean_height,SV3->point().point()-SV2->point().point());
 }
 
 template<class Tesselation>
diff --git a/lib/triangulation/PeriodicFlow.hpp b/lib/triangulation/PeriodicFlow.hpp
index 83ce6d6..fb17a0c 100644
--- a/lib/triangulation/PeriodicFlow.hpp
+++ b/lib/triangulation/PeriodicFlow.hpp
@@ -61,11 +61,11 @@ void PeriodicFlow<_Tesselation>::interpolate(Tesselation& Tes, Tesselation& NewT
 		CellHandle& newCell = *cellIt;
 		if (newCell->info().Pcondition || newCell->info().isGhost) continue;
 		CVector center ( 0,0,0 );
-		if (newCell->info().fictious()==0) for ( int k=0;k<4;k++ ) center= center + 0.25* (Tes.vertex(newCell->vertex(k)->info().id())->point()-CGAL::ORIGIN);
+		if (newCell->info().fictious()==0) for ( int k=0;k<4;k++ ) center= center + 0.25* (Tes.vertex(newCell->vertex(k)->info().id())->point().point()-CGAL::ORIGIN);
 		else {
 			Real boundPos=0; int coord=0;
 			for ( int k=0;k<4;k++ ) {
-				if (!newCell->vertex (k)->info().isFictious) center= center+0.3333333333*(Tes.vertex(newCell->vertex(k)->info().id())->point()-CGAL::ORIGIN);
+				if (!newCell->vertex (k)->info().isFictious) center= center+0.3333333333*(Tes.vertex(newCell->vertex(k)->info().id())->point().point()-CGAL::ORIGIN);
 				else {
 					coord=boundary (newCell->vertex(k)->info().id()).coordinate;
 					boundPos=boundary (newCell->vertex(k)->info().id()).p[coord];
@@ -73,7 +73,7 @@ void PeriodicFlow<_Tesselation>::interpolate(Tesselation& Tes, Tesselation& NewT
 			}
 			center=CVector(coord==0?boundPos:center[0],coord==1?boundPos:center[1],coord==2?boundPos:center[2]);
 		}
-                oldCell = Tri.locate(Point(center[0],center[1],center[2]));
+                oldCell = Tri.locate(CGT::Sphere(center[0],center[1],center[2]));
 		//FIXME: should use getInfo
                 newCell->info().p() = oldCell->info().shiftedP();
         }
@@ -206,9 +206,9 @@ void PeriodicFlow<_Tesselation>::computePermeability()
 						cell->info().facetSphereCrossSections[j][jj]=0.5*W[jj]->point().weight()*Wm3::FastInvCos1((W[permut3[jj][1]]->point()-W[permut3[jj][0]]->point())*(W[permut3[jj][2]]->point()-W[permut3[jj][0]]->point()));
 #else
 					cell->info().facetSphereCrossSections[j]=CVector(
-					W[0]->info().isFictious ? 0 : 0.5*v0.weight()*acos((v1-v0)*(v2-v0)/sqrt((v1-v0).squared_length()*(v2-v0).squared_length())),
-					W[1]->info().isFictious ? 0 : 0.5*v1.weight()*acos((v0-v1)*(v2-v1)/sqrt((v1-v0).squared_length()*(v2-v1).squared_length())),
-					W[2]->info().isFictious ? 0 : 0.5*v2.weight()*acos((v0-v2)*(v1-v2)/sqrt((v1-v2).squared_length()*(v2-v0).squared_length())));
+					W[0]->info().isFictious ? 0 : 0.5*v0.weight()*acos((v1.point()-v0.point())*(v2.point()-v0.point())/sqrt((v1.point()-v0.point()).squared_length()*(v2.point()-v0.point()).squared_length())),
+					W[1]->info().isFictious ? 0 : 0.5*v1.weight()*acos((v0.point()-v1.point())*(v2.point()-v1.point())/sqrt((v1.point()-v0.point()).squared_length()*(v2.point()-v1.point()).squared_length())),
+					W[2]->info().isFictious ? 0 : 0.5*v2.weight()*acos((v0.point()-v2.point())*(v1.point()-v2.point())/sqrt((v1.point()-v2.point()).squared_length()*(v2.point()-v0.point()).squared_length())));
 #endif
 					//FIXME: it should be possible to skip completely blocked cells, currently the problem is it segfault for undefined areas
 					//if (cell->info().blocked) continue;//We don't need permeability for blocked cells, it will be set to zero anyway
diff --git a/lib/triangulation/PeriodicFlowLinSolv.ipp b/lib/triangulation/PeriodicFlowLinSolv.ipp
index 46d00b6..fa7ad73 100644
--- a/lib/triangulation/PeriodicFlowLinSolv.ipp
+++ b/lib/triangulation/PeriodicFlowLinSolv.ipp
@@ -8,7 +8,9 @@
 #ifdef FLOW_ENGINE
 
 
-#include "CGAL/constructions/constructions_on_weighted_points_cartesian_3.h"
+#if CGAL_VERSION_NR < CGAL_VERSION_NUMBER(4,11,0)
+	#include "CGAL/constructions/constructions_on_weighted_points_cartesian_3.h"
+#endif
 #include <CGAL/Width_3.h>
 #include <iostream>
 #include <fstream>
diff --git a/lib/triangulation/RegularTriangulation.h b/lib/triangulation/RegularTriangulation.h
index 4f895dc..c655e67 100644
--- a/lib/triangulation/RegularTriangulation.h
+++ b/lib/triangulation/RegularTriangulation.h
@@ -11,7 +11,9 @@
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 #include <CGAL/Cartesian.h>
 #include <CGAL/Regular_triangulation_3.h>
-#include <CGAL/Regular_triangulation_euclidean_traits_3.h>
+#if CGAL_VERSION_NR < CGAL_VERSION_NUMBER(4,11,0)
+	#include <CGAL/Regular_triangulation_euclidean_traits_3.h>
+#endif
 #define ALPHASHAPES
 #ifdef ALPHASHAPES
 	#include <CGAL/Alpha_shape_vertex_base_3.h>
@@ -36,7 +38,12 @@ typedef CGAL::Exact_predicates_inexact_constructions_kernel K;
 //A bit faster, but gives crash eventualy
 // typedef CGAL::Cartesian<double> K;
 
-typedef CGAL::Regular_triangulation_euclidean_traits_3<K>   				Traits;
+#if CGAL_VERSION_NR < CGAL_VERSION_NUMBER(4,11,0)
+typedef CGAL::Regular_triangulation_euclidean_traits_3<K>				Traits;
+#else
+typedef K										Traits;
+#endif
+
 typedef K::Point_3									Point;
 typedef Traits::Vector_3 								CVector;
 typedef Traits::Segment_3								Segment;
@@ -44,8 +51,13 @@ typedef Traits::Segment_3								Segment;
 /** compilation inside yade: check that Real in yade is the same as Real we will define; otherwise it might make things go wrong badly (perhaps) **/
 BOOST_STATIC_ASSERT(sizeof(Traits::RT)==sizeof(Real));
 #endif
+#if CGAL_VERSION_NR < CGAL_VERSION_NUMBER(4,11,0)
 typedef Traits::RT									Real; //Dans cartesian, RT = FT
 typedef Traits::Weighted_point								Sphere;
+#else
+typedef Traits::FT									Real; //Dans cartesian, RT = FT
+typedef Traits::Weighted_point_3							Sphere;
+#endif
 typedef Traits::Plane_3									Plane;
 typedef Traits::Triangle_3								Triangle;
 typedef Traits::Tetrahedron_3								Tetrahedron;
@@ -96,8 +108,15 @@ class TriangulationTypes {
 public:
 typedef vertex_info								Vertex_Info;
 typedef cell_info								Cell_Info;
+#if CGAL_VERSION_NR < CGAL_VERSION_NUMBER(4,11,0)
 typedef CGAL::Triangulation_vertex_base_with_info_3<Vertex_Info, Traits>	Vb_info;
 typedef CGAL::Triangulation_cell_base_with_info_3<Cell_Info, Traits>		Cb_info;
+#else
+typedef CGAL::Regular_triangulation_vertex_base_3<K>				Vb0;
+typedef CGAL::Regular_triangulation_cell_base_3<K>				Rcb;
+typedef CGAL::Triangulation_vertex_base_with_info_3<Vertex_Info, Traits, Vb0>	Vb_info;
+typedef CGAL::Triangulation_cell_base_with_info_3<Cell_Info, Traits, Rcb>	Cb_info;
+#endif
 #ifdef ALPHASHAPES
 typedef CGAL::Alpha_shape_vertex_base_3<Traits,Vb_info> Vb;
 typedef CGAL::Alpha_shape_cell_base_3<Traits,Cb_info>   Fb;
diff --git a/lib/triangulation/Tesselation.h b/lib/triangulation/Tesselation.h
index 722c055..8ebf9bb 100644
--- a/lib/triangulation/Tesselation.h
+++ b/lib/triangulation/Tesselation.h
@@ -11,6 +11,8 @@
 namespace CGT {
 	
 //Since template inheritance does not automatically give access to the members of the base class, this macro can be used to declare all members at once. 
+#ifdef ALPHASHAPES
+
 #define DECLARE_TESSELATION_TYPES(baseType)\
 		typedef typename baseType::RTriangulation		 	RTriangulation;\
 		typedef typename baseType::AlphaShape		 		AlphaShape;\
@@ -35,6 +37,33 @@ namespace CGT {
 		typedef typename baseType::ListPoint				ListPoint;\
 		typedef typename baseType::VCellIterator			VCellIterator;
 
+#else
+
+#define DECLARE_TESSELATION_TYPES(baseType)\
+		typedef typename baseType::RTriangulation		 	RTriangulation;\
+		typedef typename baseType::VertexInfo				VertexInfo;\
+		typedef typename baseType::CellInfo				CellInfo;\
+		typedef typename baseType::VertexIterator			VertexIterator;\
+		typedef typename baseType::VertexHandle				VertexHandle;\
+		typedef typename baseType::FiniteVerticesIterator		FiniteVerticesIterator;\
+		typedef typename baseType::CellIterator				CellIterator;\
+		typedef typename baseType::FiniteCellsIterator			FiniteCellsIterator;\
+		typedef typename baseType::CellCirculator			CellCirculator;\
+		typedef typename baseType::CellHandle				CellHandle;\
+		typedef typename baseType::Facet				Facet;\
+		typedef typename baseType::FacetIterator			FacetIterator;\
+		typedef typename baseType::FacetCirculator			FacetCirculator;\
+		typedef typename baseType::FiniteFacetsIterator			FiniteFacetsIterator;\
+		typedef typename baseType::LocateType				LocateType;\
+		typedef typename baseType::EdgeIterator				EdgeIterator;\
+		typedef typename baseType::FiniteEdgesIterator			FiniteEdgesIterator;\
+		typedef typename baseType::VectorVertex				VectorVertex;\
+		typedef typename baseType::VectorCell				VectorCell;\
+		typedef typename baseType::ListPoint				ListPoint;\
+		typedef typename baseType::VCellIterator			VCellIterator;
+
+#endif
+
 // Classe Tesselation, contient les fonctions permettant de calculer la Tessalisation
 // d'une RTriangulation et de stocker les centres dans chacune de ses cellules
 
@@ -44,7 +73,9 @@ class _Tesselation
 {
 public:
 	typedef typename TT::RTriangulation							RTriangulation;
+#ifdef ALPHASHAPES
 	typedef typename TT::AlphaShape						 		AlphaShape;
+#endif
 	typedef typename TT::Vertex_Info							VertexInfo;
 	typedef typename TT::Cell_Info								CellInfo;
 	typedef typename RTriangulation::Vertex_iterator		 			VertexIterator;
diff --git a/lib/triangulation/Tesselation.ipp b/lib/triangulation/Tesselation.ipp
index a73a45b..32b983e 100644
--- a/lib/triangulation/Tesselation.ipp
+++ b/lib/triangulation/Tesselation.ipp
@@ -126,8 +126,8 @@ typename _Tesselation<TT>::RTriangulation & _Tesselation<TT>::Triangulation ( vo
 template<class TT>
 Real _Tesselation<TT>::Volume ( FiniteCellsIterator cell )
 {
-	return ( Tetrahedron ( cell->vertex ( 0 )->point(), cell->vertex ( 1 )->point(),
-						 cell->vertex ( 2 )->point(), cell->vertex ( 3 )->point() ) ).volume();
+	return ( Tetrahedron ( cell->vertex ( 0 )->point().point(), cell->vertex ( 1 )->point().point(),
+						 cell->vertex ( 2 )->point().point(), cell->vertex ( 3 )->point().point() ) ).volume();
 }
 
 template<class TT>
@@ -230,13 +230,13 @@ void _Tesselation<TT>::AssignPartialVolume ( FiniteEdgesIterator& ed_it )
 		{
 			if ( !isFictious1 )
 			{
-				r = std::abs ( ( Tetrahedron ( ed_it->first->vertex ( ed_it->second )->point(), cell0->info(), cell1->info(), cell2->info() ) ).volume() );
+				r = std::abs ( ( Tetrahedron ( ed_it->first->vertex ( ed_it->second )->point().point(), cell0->info(), cell1->info(), cell2->info() ) ).volume() );
 				( ed_it->first )->vertex ( ed_it->second )->info().v() += r;
 				TotalFiniteVoronoiVolume+=r;
 			}
 			if ( !isFictious2 )
 			{
-				r = std::abs ( ( Tetrahedron ( ed_it->first->vertex ( ed_it->third )->point(), cell0->info(),  cell1->info(), cell2->info() ) ).volume() );
+				r = std::abs ( ( Tetrahedron ( ed_it->first->vertex ( ed_it->third )->point().point(), cell0->info(),  cell1->info(), cell2->info() ) ).volume() );
 				ed_it->first->vertex ( ed_it->third )->info().v() +=r;
 				TotalFiniteVoronoiVolume+=r;
 			}
diff --git a/pkg/dem/TesselationWrapper.cpp b/pkg/dem/TesselationWrapper.cpp
index 0b3700c..440a279 100644
--- a/pkg/dem/TesselationWrapper.cpp
+++ b/pkg/dem/TesselationWrapper.cpp
@@ -29,17 +29,17 @@ struct RTraits_for_spatial_sort : public CGT::SimpleTriangulationTypes::RTriangu
 
 	struct Less_x_3 {
 		bool operator()(const Point_3& p,const Point_3& q) const {
-			return Gt::Less_x_3()(* (p.first),* (q.first));
+			return Gt::Less_x_3()( p.first->point() , q.first->point() );
 		}
 	};
 	struct Less_y_3 {
 		bool operator()(const Point_3& p,const Point_3& q) const {
-			return Gt::Less_y_3()(* (p.first),* (q.first));
+			return Gt::Less_y_3()( p.first->point(), q.first->point());
 		}
 	};
 	struct Less_z_3 {
 		bool operator()(const Point_3& p,const Point_3& q) const {
-			return Gt::Less_z_3()(* (p.first),* (q.first));
+			return Gt::Less_z_3()( p.first->point(), q.first->point());
 		}
 	};
 	Less_x_3  less_x_3_object() const {return Less_x_3();}
diff --git a/pkg/pfv/FlowEngine.ipp.in b/pkg/pfv/FlowEngine.ipp.in
index 40b1e15..eaaa61d 100644
--- a/pkg/pfv/FlowEngine.ipp.in
+++ b/pkg/pfv/FlowEngine.ipp.in
@@ -179,7 +179,7 @@ void TemplateFlowEngine_@TEMPLATE_FLOW_NAME@<_CellInfo,_VertexInfo,_Tesselation,
 template< class _CellInfo, class _VertexInfo, class _Tesselation, class solverT >
 void TemplateFlowEngine_@TEMPLATE_FLOW_NAME@<_CellInfo,_VertexInfo,_Tesselation,solverT>::imposeFlux ( Vector3r pos, Real flux){
         solver->imposedF.push_back ( pair<CGT::Point,Real> ( CGT::Point ( pos[0],pos[1],pos[2] ), flux ) );
-	CellHandle cell=solver->T[solver->currentTes].Triangulation().locate(CGT::Point(pos[0],pos[1],pos[2]));
+	CellHandle cell=solver->T[solver->currentTes].Triangulation().locate(CGT::Sphere(pos[0],pos[1],pos[2]));
 	if (cell->info().isGhost) cerr<<"Imposing pressure in a ghost cell."<<endl;
 	for (unsigned int kk=0;kk<solver->IPCells.size();kk++){
 		if (cell==solver->IPCells[kk]) cerr<<"Both flux and pressure are imposed in the same cell."<<endl;
diff --git a/pkg/pfv/PeriodicFlowEngine.cpp b/pkg/pfv/PeriodicFlowEngine.cpp
index d61c587..86e1382 100644
--- a/pkg/pfv/PeriodicFlowEngine.cpp
+++ b/pkg/pfv/PeriodicFlowEngine.cpp
@@ -335,12 +335,12 @@ void PeriodicFlowEngine::locateCell ( CellHandle baseCell, unsigned int& index,
 	Vector3i period;
 
 	if (baseCell->info().fictious()==0)
-		for ( int k=0;k<4;k++ ) center+= 0.25*makeVector3r (baseCell->vertex(k)->point());
+		for ( int k=0;k<4;k++ ) center+= 0.25*makeVector3r (baseCell->vertex(k)->point().point());
 	else {
 		
 		Real boundPos=0; int coord=0;
 		for ( int k=0;k<4;k++ ) {
-			if ( !baseCell->vertex ( k )->info().isFictious ) center+= 0.3333333333*makeVector3r ( baseCell->vertex ( k )->point() );
+			if ( !baseCell->vertex ( k )->info().isFictious ) center+= 0.3333333333*makeVector3r ( baseCell->vertex ( k )->point().point() );
 			else {
 				coord=flow.boundary ( baseCell->vertex ( k )->info().id() ).coordinate;
 				boundPos=flow.boundary ( baseCell->vertex ( k )->info().id() ).p[coord];}
@@ -354,9 +354,7 @@ void PeriodicFlowEngine::locateCell ( CellHandle baseCell, unsigned int& index,
 			baseInfo.isGhost=false;
 			return;
 		}
-		CellHandle ch= Tri.locate ( CGT::Point ( wdCenter[0],wdCenter[1],wdCenter[2] )
-// 					     ,/*hint*/ v0
-					     );
+		CellHandle ch= Tri.locate ( CGT::Sphere ( wdCenter[0],wdCenter[1],wdCenter[2] ) );
 		baseInfo.period[0]=period[0];
 		baseInfo.period[1]=period[1];
 		baseInfo.period[2]=period[2];