From 97565c65c436821da304f75fd9cd5c58f25441ab Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Thu, 7 Jan 2016 14:37:08 +0000
Subject: [PATCH] random shuffle added in delaunay refinement
---
Mesh/delaunay3d.cpp | 18 +++++++++---------
Mesh/delaunay_refinement.cpp | 6 ++++--
2 files changed, 13 insertions(+), 11 deletions(-)
diff --git a/Mesh/delaunay3d.cpp b/Mesh/delaunay3d.cpp
index 8eaa669da5..d4f341eecc 100644
--- a/Mesh/delaunay3d.cpp
+++ b/Mesh/delaunay3d.cpp
@@ -80,7 +80,6 @@ struct HilbertSortB
bbox *= 1.01;
Vertex**pv = &v[0];
int depth;
- // MultiscaleSortHilbert(pv, (int)v.size(), 64, 0.125,&depth);
indices.clear();
MultiscaleSortHilbert(pv, (int)v.size(), 64, .125,&depth,indices);
indices.push_back(v.size());
@@ -323,7 +322,7 @@ static void delaunayCavity2 (Tet *t,
Vertex *v,
cavityContainer &cavity,
connContainer &bnd,
- int thread, int iPnt, int PP){
+ int thread, int iPnt){
t->set(thread, iPnt); // Mark the triangle
cavity.push_back(t);
for (int iNeigh=0; iNeigh<4 ; iNeigh++){
@@ -332,7 +331,6 @@ static void delaunayCavity2 (Tet *t,
bnd.push_back(conn(t->getFace(iNeigh),iNeigh,NULL));
}
else if (neigh == prev){
- if (PP == 40)printf("coucou\n");
}
else if (!neigh->inSphere(v,thread)){
Face f = t->getFace(iNeigh);
@@ -340,7 +338,7 @@ static void delaunayCavity2 (Tet *t,
// neigh->set(thread, iPnt);
}
else if (!(neigh->isSet(thread, iPnt))) {
- delaunayCavity2 (neigh, t, v, cavity,bnd,thread, iPnt,PP);
+ delaunayCavity2 (neigh, t, v, cavity,bnd,thread, iPnt);
}
}
}
@@ -519,6 +517,7 @@ static Tet* randomTet (int thread, tetContainer &allocator ){
}
+//#define _VERBOSE 1
void delaunayTrgl (const unsigned int numThreads,
const unsigned int NPTS_AT_ONCE,
unsigned int Npts,
@@ -532,7 +531,7 @@ void delaunayTrgl (const unsigned int numThreads,
// checkLocalDelaunayness(T, "initial");
std::vector<int> invalidCavities(numThreads);
- std::vector<int> cashMisses(numThreads, 0);
+ std::vector<int> cacheMisses(numThreads, 0);
unsigned int maxLocSizeK = 0;
for (unsigned int i = 0; i < numThreads * NPTS_AT_ONCE; i++){
@@ -609,9 +608,11 @@ void delaunayTrgl (const unsigned int numThreads,
cavityContainer &cavityK = cavity[K];
connContainer &bndK = bnd[K];
for (unsigned int i=0; i<cavityK.size(); i++)cavityK[i]->unset(myThread,K);
+ // for (unsigned int i=0; i< bndK.size(); i++)if(bndK[i].t)bndK[i].t->unset(myThread,K);
cavityK.clear(); bndK.clear();
delaunayCavity(t[K], NULL, vToAdd[K], cavityK, bndK, myThread, K);
for (unsigned int i=0; i<cavityK.size(); i++)cavityK[i]->unset(myThread,K);
+ // for (unsigned int i=0; i< bndK.size(); i++)if(bndK[i].t)bndK[i].t->unset(myThread,K);
// verify the cavity
for (unsigned int i=0; i< bndK.size(); i++) {
double val = robustPredicates::orient3d((double*)bndK[i].f.V[0],
@@ -654,7 +655,7 @@ void delaunayTrgl (const unsigned int numThreads,
else {
t = allocator.newTet(myThread);
}
- if (i < cSize && t->V[0]->_thread != myThread)cashMisses[myThread]++;
+ if (i < cSize && t->V[0]->_thread != myThread)cacheMisses[myThread]++;
t->setVerticesNoTest (bndK[i].f.V[0], bndK[i].f.V[1], bndK[i].f.V[2], vToAdd[K]);
Tet *neigh = bndK[i].t;
t->T[0] = neigh;
@@ -692,13 +693,11 @@ void delaunayTrgl (const unsigned int numThreads,
#ifdef _VERBOSE
- printf("total time for %d points %12.5E seconds\n",Npts,(double)(_t2-_t1));
- printf("%8d tets generated (%6f/sec)\n",T.size(),(double) T.size()/ ((double)(_t2-_t1)));
printf("average searches per point %12.5E\n",totSearchGlob/Npts);
printf("average size for del cavity %12.5E\n",totCavityGlob/Npts);
printf("cash misses: ");
for (int i=0;i<numThreads;i++){
- printf("%4d ",(int)cashMisses[i]);
+ printf("%4d ",(int)cacheMisses[i]);
}
printf("\n");
@@ -838,6 +837,7 @@ void delaunayTriangulation (const int numThreads,
robustPredicates::exactinit(1,maxx,maxy,maxz);
+ // FIXME numThreads
tetContainer allocator (1, S.size() * 10);
Vertex *box[8];
delaunayTriangulation (numThreads, nptsatonce, _vertices, box, allocator);
diff --git a/Mesh/delaunay_refinement.cpp b/Mesh/delaunay_refinement.cpp
index d5b84701e5..b7598f82a3 100644
--- a/Mesh/delaunay_refinement.cpp
+++ b/Mesh/delaunay_refinement.cpp
@@ -311,7 +311,7 @@ void edgeBasedRefinement (const int numThreads,
// fill up old Datastructures
- tetContainer allocator (numThreads,1000000);
+ tetContainer allocator (numThreads,3000000);
std::vector<Vertex *> _vertices;
@@ -401,6 +401,9 @@ void edgeBasedRefinement (const int numThreads,
filterVertices (numThreads, _filter, _vertices, add, _fx, NULL);
double t3 = Cpu();
if (add.empty())break;
+ // randomize vertices (EXTREMELY IMPORTANT FOR NOT DETERIORATING PERFORMANCE)
+ std::random_shuffle(add.begin(), add.end());
+ // sort them using BRIO
std::vector<int> indices;
SortHilbert(add, indices);
double t4 = Cpu();
@@ -441,7 +444,6 @@ void edgeBasedRefinement (const int numThreads,
std::map<Edge,double> _sizes;
for (unsigned int i=0; i< allocator.size(0);i++){
Tet *tt = allocator (0,i);
- MVertex *mvs[4];
if (tt->V[0]){
for (int j=0;j<6;j++){
Edge e = tt->getEdge(j);
--
GitLab