diff --git a/Mesh/delaunay3d.cpp b/Mesh/delaunay3d.cpp
index c97434a7456f85be18ae4c3a3784b48886d9f250..dea924898fe661474f67b61ae61b5fd7a35f6194 100644
--- a/Mesh/delaunay3d.cpp
+++ b/Mesh/delaunay3d.cpp
@@ -30,6 +30,9 @@ const int MEASURE_BARR = 0;
#include <numa.h>
#endif
+int Tet::in_sphere_counter = 0;
+
+
struct HilbertSortB
{
// The code for generating table transgc
@@ -333,9 +336,8 @@ static void delaunayCavity2 (Tet *t,
else if (neigh == prev){
}
else if (!neigh->inSphere(v,thread)){
- Face f = t->getFace(iNeigh);
- bnd.push_back(conn(f,iNeigh,neigh));
- // neigh->set(thread, iPnt);
+ bnd.push_back(conn(t->getFace(iNeigh),iNeigh,neigh));
+ neigh->set(thread, iPnt);
}
else if (!(neigh->isSet(thread, iPnt))) {
delaunayCavity2 (neigh, t, v, cavity,bnd,thread, iPnt);
@@ -399,7 +401,6 @@ Tet* walk (Tet *t, Vertex *v, int maxx, double &totSearch, int thread){
}
-
void __print (const char *name, connContainer &conn){
FILE *f = fopen(name,"w");
fprintf(f,"View \"\"{\n");
@@ -517,7 +518,7 @@ static Tet* randomTet (int thread, tetContainer &allocator ){
}
-#define _VERBOSE 1
+//#define _VERBOSE 1
void delaunayTrgl (const unsigned int numThreads,
const unsigned int NPTS_AT_ONCE,
unsigned int Npts,
@@ -528,6 +529,8 @@ void delaunayTrgl (const unsigned int numThreads,
double totCavityGlob=0;
#endif
+ double t1,t2=0,t3=0,t4=0;
+
// checkLocalDelaunayness(allocator, 0, "initial");
std::vector<int> invalidCavities(numThreads);
@@ -588,6 +591,7 @@ void delaunayTrgl (const unsigned int numThreads,
std::vector<Tet*> t(NPTS_AT_ONCE);
// double c1 = Cpu();
// FIND SEEDS
+ t1 = Cpu();
for (unsigned int K=0; K< NPTS_AT_ONCE; K++) {
vToAdd[K] = iPGlob < locSizeK[K] ? &allocatedVerts[K][iPGlob] : NULL;
if(vToAdd[K]){
@@ -602,24 +606,24 @@ void delaunayTrgl (const unsigned int numThreads,
}
}
}
+ t2+= Cpu() - t1;
// double c1 = Cpu();
// BUILD CAVITIES
+ t1 = Cpu();
for (unsigned int K=0; K< NPTS_AT_ONCE; K++) {
if(vToAdd[K]){
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);
+ 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);
+ delaunayCavity2(t[K], NULL, vToAdd[K], cavityK, bndK, myThread, K);
// verify the cavity
for (unsigned int i=0; i< bndK.size(); i++) {
- double val = robustPredicates::orient3d((double*)bndK[i].f.V[0],
- (double*)bndK[i].f.V[1],
- (double*)bndK[i].f.V[2],
- (double*)vToAdd[K]);
+ const double val = robustPredicates::orient3d((double*)bndK[i].f.V[0],
+ (double*)bndK[i].f.V[1],
+ (double*)bndK[i].f.V[2],
+ (double*)vToAdd[K]);
if (val <= 0 ) {
vToAdd[K] = NULL;
invalidCavities [myThread]++;
@@ -629,6 +633,7 @@ void delaunayTrgl (const unsigned int numThreads,
}
}
+ t3 += Cpu() - t1;
#pragma omp barrier
@@ -636,7 +641,8 @@ void delaunayTrgl (const unsigned int numThreads,
if (!vToAdd[K])ok[K]=false;
else ok[K] = canWeProcessCavity (cavity[K], myThread, K);
}
-
+ t1 = Cpu();
+
for (unsigned int K=0; K< NPTS_AT_ONCE; K++) {
if (ok[K]){
cavityContainer &cavityK = cavity[K];
@@ -670,6 +676,7 @@ void delaunayTrgl (const unsigned int numThreads,
}
}
}
+ t4 += Cpu() - t1;
}
#ifdef _VERBOSE
#pragma omp critical
@@ -679,6 +686,11 @@ void delaunayTrgl (const unsigned int numThreads,
}
#endif
#pragma omp barrier
+ // clear last cavity
+ for (unsigned int K=0; K< NPTS_AT_ONCE; K++) {
+ for (unsigned int i=0; i<cavity[K].size(); i++)cavity[K][i]->unset(myThread,K);
+ for (unsigned int i=0; i<bnd[K].size(); i++)if(bnd[K][i].t)bnd[K][i].t->unset(myThread,K);
+ }
}
if (invalidCavities[0])Msg::Warning("%d invalid cavities",invalidCavities[0]);
@@ -686,12 +698,14 @@ void delaunayTrgl (const unsigned int numThreads,
//checkLocalDelaunayness(T,"final");
+ printf(" %12.5E %12.5E %12.5E tot %12.5E \n",t2,t3,t4,t2+t3+t4);
+
#ifdef _VERBOSE
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)cacheMisses[i]);
+ printf("cache misses: ");
+ for (unsigned int i=0;i<numThreads;i++){
+ printf("%4ud ",(int)cacheMisses[i]);
}
printf("\n");
@@ -815,6 +829,8 @@ void delaunayTriangulation (const int numThreads,
}
double d = 1*sqrt ( maxx*maxx + maxy*maxy + maxz*maxz );
+ tetContainer allocator (1, S.size() * 10);
+
for (unsigned int i=0;i<N;i++){
MVertex *mv = S[i];
// FIXME : should be zero !!!!
@@ -824,6 +840,7 @@ void delaunayTriangulation (const int numThreads,
mv->x() += dx;
mv->y() += dy;
mv->z() += dz;
+ // Vertex *v = new Vertex (mv->x(),mv->y(),mv->z(),1.e22,i+1);
Vertex *v = new Vertex (mv->x(),mv->y(),mv->z(),1.e22,i+1);
_vertices.push_back(v);
_temp [v->getNum()] = mv;
@@ -832,7 +849,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/delaunay3d_private.h b/Mesh/delaunay3d_private.h
index f1a38f2a54520641ab8651a3b5c1e87d8da23e68..def190ba3e0a872a298663d2f74b86047acb7149 100644
--- a/Mesh/delaunay3d_private.h
+++ b/Mesh/delaunay3d_private.h
@@ -36,7 +36,7 @@ public :
inline double &z() {return _x[2];}
inline double &lc() {return _lc;}
inline operator double *() { return _x; }
- Vertex (double X, double Y, double Z, double lc, int num = 0) : _num(num), _thread(0){
+ Vertex (double X=0, double Y=0, double Z=0, double lc=0, int num = 0) : _num(num), _thread(0){
_x[0]=X; _x[1]=Y;_x[2]=Z; _lc = lc;}
Vertex operator + (const Vertex &other){
return Vertex (x()+other.x(),y()+other.y(),z()+other.z(),other.lc() + _lc);
@@ -68,7 +68,7 @@ inline double orientationTestFast(double *pa, double *pb, double *pc, double *pd
inline bool inSphereTest_s (Vertex *va, Vertex *vb, Vertex *vc, Vertex *vd , Vertex *ve){
double val = robustPredicates::insphere ((double*)va,(double*)vb,(double*)vc,(double*)vd,(double*)ve);
if (val == 0.0){
- Msg::Debug("symbolic perturbation needed vol %22.15E",orientationTestFast((double*)va,(double*)vb,(double*)vc,(double*)vd));
+ Msg::Info("symbolic perturbation needed vol %22.15E",orientationTestFast((double*)va,(double*)vb,(double*)vc,(double*)vd));
int count;
// symbolic perturbation
Vertex *pt[5] = {va,vb,vc,vd,ve};
@@ -160,7 +160,8 @@ struct Tet {
Tet *T[4];
Vertex *V[4];
CHECKTYPE _bitset [MAX_NUM_THREADS_];
- char _modified;
+ bool _modified;
+ static int in_sphere_counter;
Tet () : _modified(true){
V[0] = V[1] = V[2] = V[3] = NULL;
T[0] = T[1] = T[2] = T[3] = NULL;
@@ -169,6 +170,7 @@ struct Tet {
int setVerticesNoTest (Vertex *v0, Vertex *v1, Vertex *v2, Vertex *v3){
_modified=true;
V[0] = v0; V[1] = v1; V[2] = v2; V[3] = v3;
+ // for (int i=0;i<4;i++)_copy[i] = *V[i];
return 1;
}
int setVertices (Vertex *v0, Vertex *v1, Vertex *v2, Vertex *v3){
@@ -176,11 +178,13 @@ struct Tet {
double val = robustPredicates::orient3d((double*)v0,(double*)v1,(double*)v2,(double*)v3);
V[0] = v0; V[1] = v1; V[2] = v2; V[3] = v3;
if (val > 0){
+ // for (int i=0;i<4;i++)_copy[i] = *V[i];
return 1;
}
else if (val < 0){
// throw;
V[0] = v1; V[1] = v0; V[2] = v2; V[3] = v3;
+ // for (int i=0;i<4;i++)_copy[i] = *V[i];
return -1;
}
else {
@@ -220,13 +224,10 @@ struct Tet {
return Edge (std::min(V[edg[k][0]],V[edg[k][1]]),
std::max(V[edg[k][0]],V[edg[k][1]]));
}
- inline bool inSphere (Vertex *vd, int thread) const{
+ inline bool inSphere (Vertex *vd, int thread) {
+ in_sphere_counter++;
return inSphereTest_s (V[0],V[1],V[2],V[3],vd);
}
- Vertex centroid () const {
- return (*V[0]+*V[1]+*V[2]+*V[3])*0.25;
- }
-
};
struct conn {
@@ -246,32 +247,32 @@ struct conn {
// tetrahedra (one per thread)
-struct aBunchOfTets {
- std::vector<Tet*> _all;
+template <class T> class aBunchOfStuff {
+public:
+ std::vector<T*> _all;
unsigned int _current;
unsigned int _nbAlloc;
unsigned int size () {
return _current + (_all.size()-1) * _nbAlloc;
}
- Tet* operator () (unsigned int i){
+ T* operator () (unsigned int i){
unsigned int _array = i / _nbAlloc;
unsigned int _offset = i % _nbAlloc;
// printf("%d %d %d\n",i,_array,_offset);
return _all [_array] + _offset;
}
- aBunchOfTets (unsigned int s) : _current(0), _nbAlloc (s) {
- _all.push_back(new Tet [_nbAlloc]);
+ aBunchOfStuff (unsigned int s) : _current(0), _nbAlloc (s) {
+ _all.push_back(new T [_nbAlloc]);
}
- ~aBunchOfTets(){
+ ~aBunchOfStuff(){
for (unsigned int i=0;i<_all.size();i++){
delete [] _all[i];
}
}
- Tet* newTet() {
-
+ T* newStuff() {
if (_current == _nbAlloc){
- _all.push_back(new Tet [_nbAlloc]);
- printf("REALLOCATION %d\n",_all.size());
+ _all.push_back(new T [_nbAlloc]);
+ // printf("REALLOCATION %d\n",_all.size());
_current = 0;
}
_current++;
@@ -281,21 +282,37 @@ struct aBunchOfTets {
// tetAllocator owns the tets that have been allocated by itself
class tetContainer {
- std::vector<aBunchOfTets*> _perThread;
-public:
+ // std::vector<aBunchOfStuff<Vertex> *> _perThreadV;
+ std::vector<aBunchOfStuff<Tet> *> _perThread;
+ public:
unsigned int size(int thread) const {return _perThread[thread]->size();}
- Tet * operator () (int thread, int j) const {return (*_perThread[thread]) (j);}
+ // unsigned int sizeV(int thread) const {return _perThreadV[thread]->size();}
+ inline Tet * operator () (int thread, int j) const {return (*_perThread[thread]) (j);}
+ // Vertex * getVertex(int thread, int j) const {return (*_perThreadV[thread]) (j);}
tetContainer (int nbThreads, int preallocSizePerThread) {
// FIXME !!!
if (nbThreads != 1) throw;
- _perThread.push_back(new aBunchOfTets (preallocSizePerThread) );
- }
- Tet * newTet(int thread) {
- return _perThread[thread]->newTet();
- }
+ _perThread.push_back (new aBunchOfStuff<Tet> (preallocSizePerThread) );
+ // _perThreadV.push_back(new aBunchOfStuff<Vertex> (preallocSizePerThread) );
+ }
+ inline Tet * newTet(int thread) {
+ return _perThread[thread]->newStuff();
+ }
+ /*
+ Vertex * newVertex(int thread, double x, double y, double z, double lc, int num) {
+ Vertex *v = _perThreadV[thread]->newStuff();
+ v->x() = x;
+ v->y() = y;
+ v->z() = z;
+ v->lc() = lc;
+ v->setNum(num);
+ return v;
+ }
+ */
~tetContainer(){
// FIXME !!!
- delete _perThread[0];
+ delete _perThread [0];
+ // delete _perThreadV[0];
}
};
diff --git a/Mesh/delaunay_refinement.cpp b/Mesh/delaunay_refinement.cpp
index 0628eb7745225cd9d9efc4408b69316e33b43d27..7f7fedf6641ed2052332cbb760b499f9507fcf03 100644
--- a/Mesh/delaunay_refinement.cpp
+++ b/Mesh/delaunay_refinement.cpp
@@ -327,6 +327,16 @@ void edgeBasedRefinement (const int numThreads,
}
}
+
+ FILE *f = fopen ("pts_init.dat","w");
+ fprintf(f,"%d\n",all.size());
+ for (std::set<MVertex*>::iterator it = all.begin();it !=all.end(); ++it){
+ MVertex *mv = *it;
+ fprintf(f,"%12.5E %12.5E %12.5E\n",mv->x(),mv->y(),mv->z());
+ }
+ fclose(f);
+
+
_vertices.resize(all.size());
int counter=0;
for (std::set<MVertex*>::iterator it = all.begin();it !=all.end(); ++it){
@@ -393,6 +403,7 @@ void edgeBasedRefinement (const int numThreads,
Msg::Info("----------------------------------- SATUR FILTR SORTH DELNY TIME TETS");
double __t__ = Cpu();
+ Tet::in_sphere_counter = 0;
while(1){
std::vector<Vertex*> add;
double t1 = Cpu();
@@ -408,8 +419,16 @@ void edgeBasedRefinement (const int numThreads,
SortHilbert(add, indices);
double t4 = Cpu();
delaunayTrgl (1,1,add.size(), &add,allocator);
- add_all.insert (add_all.end(), add.begin(), add.end());
double t5 = Cpu();
+ add_all.insert (add_all.end(), add.begin(), add.end());
+ // if (iter == 1){
+ // FILE *f = fopen ("pts.dat","w");
+ // fprintf(f,"%d\n",add.size());
+ // for (unsigned int i=0;i<add.size();i++){
+ // fprintf(f,"%12.5E %12.5E %12.5E\n",add[i]->x(),add[i]->y(),add[i]->z());
+ // }
+ // fclose(f);
+ // }
Msg::Info("IT %3d %6d points added, timings %5.2f %5.2f %5.2f %5.2f %5.2f %5d",iter,add.size(),
(t2-t1),
(t3-t2),
@@ -417,10 +436,13 @@ void edgeBasedRefinement (const int numThreads,
(t5-t4),
(t5-__t__),
allocator.size(0));
+ printf("%d calls to inSphere\n",Tet::in_sphere_counter);
iter++;
}
}
+
+
for (unsigned int i=0; i< allocator.size(0);i++){
Tet *tt = allocator (0,i);
MVertex *mvs[4];
diff --git a/Numeric/robustPredicates.cpp b/Numeric/robustPredicates.cpp
index e2441a484b033174a313f3f0692aabd61a429d05..846c3a60b1ddb42eb19e66b8e60bd540103d4c5e 100644
--- a/Numeric/robustPredicates.cpp
+++ b/Numeric/robustPredicates.cpp
@@ -2353,7 +2353,7 @@ REAL orient3d(REAL *pa, REAL *pb, REAL *pc, REAL *pd)
if (det < -o3dstaticfilter) return det;
}
-
+
permanent = (Absolute(bdxcdy) + Absolute(cdxbdy)) * Absolute(adz)
+ (Absolute(cdxady) + Absolute(adxcdy)) * Absolute(bdz)
+ (Absolute(adxbdy) + Absolute(bdxady)) * Absolute(cdz);
@@ -2362,7 +2362,7 @@ REAL orient3d(REAL *pa, REAL *pb, REAL *pc, REAL *pd)
return det;
}
- return orient3dadapt(pa, pb, pc, pd, permanent);
+ return orient3dadapt(pa, pb, pc, pd, permanent);
}
#endif // ifdef INEXACT_GEOM_PRED