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];
_______________________________________________ Mailing list: https://launchpad.net/~yade-dev Post to : yade-dev@lists.launchpad.net Unsubscribe : https://launchpad.net/~yade-dev More help : https://help.launchpad.net/ListHelp