Skip to content
Snippets Groups Projects
Commit 590ce8da authored by Christophe Geuzaine's avatar Christophe Geuzaine
Browse files

fix compile

parent 46901577
No related branches found
No related tags found
No related merge requests found
...@@ -11,8 +11,8 @@ ...@@ -11,8 +11,8 @@
#include <vector> #include <vector>
#include <algorithm> #include <algorithm>
#include <math.h> #include <math.h>
#include <sys/time.h>
#include "SBoundingBox3d.h" #include "SBoundingBox3d.h"
#include "OS.h"
#include "delaunay3d_private.h" #include "delaunay3d_private.h"
#include "delaunay3d.h" #include "delaunay3d.h"
#include "MVertex.h" #include "MVertex.h"
...@@ -305,24 +305,7 @@ double walltime( double *t0 ) ...@@ -305,24 +305,7 @@ double walltime( double *t0 )
#ifdef _OPENMP #ifdef _OPENMP
return omp_get_wtime(); return omp_get_wtime();
#else #else
double mic, time; return GetTimeInSeconds();
double mega = 0.000001;
struct timeval tp;
struct timezone tzp;
static long base_sec = 0;
static long base_usec = 0;
(void) gettimeofday(&tp,&tzp);
if (base_sec == 0)
{
base_sec = tp.tv_sec;
base_usec = tp.tv_usec;
}
time = (double) (tp.tv_sec - base_sec);
mic = (double) (tp.tv_usec - base_usec);
time = (time + mic * mega) - *t0;
return(time);
#endif #endif
} }
...@@ -642,7 +625,7 @@ void delaunayTrgl (const unsigned int numThreads, ...@@ -642,7 +625,7 @@ void delaunayTrgl (const unsigned int numThreads,
#pragma omp barrier #pragma omp barrier
cavitySize[myThread] = 0; cavitySize[myThread] = 0;
std::vector<Tet*> t(NPTS_AT_ONCE); std::vector<Tet*> t(NPTS_AT_ONCE);
// clock_t c1 = clock(); // double c1 = Cpu();
for (unsigned int K=0; K< NPTS_AT_ONCE; K++) { for (unsigned int K=0; K< NPTS_AT_ONCE; K++) {
vToAdd[K] = iPGlob < locSizeK[K] ? &allocatedVerts[K][iPGlob] : NULL; vToAdd[K] = iPGlob < locSizeK[K] ? &allocatedVerts[K][iPGlob] : NULL;
if(vToAdd[K]){ if(vToAdd[K]){
...@@ -681,7 +664,7 @@ void delaunayTrgl (const unsigned int numThreads, ...@@ -681,7 +664,7 @@ void delaunayTrgl (const unsigned int numThreads,
} }
} }
} }
// clock_t c1 = clock(); // double c1 = Cpu();
for (unsigned int K=0; K< NPTS_AT_ONCE; K++) { for (unsigned int K=0; K< NPTS_AT_ONCE; K++) {
if(vToAdd[K]){ if(vToAdd[K]){
cavityContainer &cavityK = cavity[K]; cavityContainer &cavityK = cavity[K];
...@@ -719,7 +702,7 @@ void delaunayTrgl (const unsigned int numThreads, ...@@ -719,7 +702,7 @@ void delaunayTrgl (const unsigned int numThreads,
else ok[K] = canWeProcessCavity (cavity[K], myThread, K); else ok[K] = canWeProcessCavity (cavity[K], myThread, K);
} }
// clock_t ck = clock(); // double ck = Cpu();
// std::set<Tet*> touched; // std::set<Tet*> touched;
for (unsigned int K=0; K< NPTS_AT_ONCE; K++) { for (unsigned int K=0; K< NPTS_AT_ONCE; K++) {
if (ok[K]){ if (ok[K]){
......
...@@ -7,6 +7,7 @@ ...@@ -7,6 +7,7 @@
#include <algorithm> #include <algorithm>
#include <math.h> #include <math.h>
#include "GmshMessage.h" #include "GmshMessage.h"
#include "OS.h"
#include "SPoint3.h" #include "SPoint3.h"
#include "delaunay3d_private.h" #include "delaunay3d_private.h"
#include "delaunay3d.h" #include "delaunay3d.h"
...@@ -39,10 +40,10 @@ struct edgeContainer ...@@ -39,10 +40,10 @@ struct edgeContainer
std::vector<Edge> &v = _hash[h %_hash.size()]; std::vector<Edge> &v = _hash[h %_hash.size()];
AVGSEARCH+=v.size(); AVGSEARCH+=v.size();
for (unsigned int i=0; i< v.size();i++)if (e == v[i]) {return false;} for (unsigned int i=0; i< v.size();i++)if (e == v[i]) {return false;}
return true; return true;
} }
bool addNewEdge (const Edge &e) bool addNewEdge (const Edge &e)
{ {
size_t h = (size_t) e.first >> 3; size_t h = (size_t) e.first >> 3;
std::vector<Edge> &v = _hash[h %_hash.size()]; std::vector<Edge> &v = _hash[h %_hash.size()];
...@@ -59,14 +60,14 @@ struct IPT { ...@@ -59,14 +60,14 @@ struct IPT {
_x1(x1),_x2(x2),_x3(x3),_x4(x4){}; _x1(x1),_x2(x2),_x3(x3),_x4(x4){};
}; };
double adaptiveTrapezoidalRule (SPoint3 p1 , SPoint3 p2 , double adaptiveTrapezoidalRule (SPoint3 p1 , SPoint3 p2 ,
double lc1 , double lc2 , double lc1 , double lc2 ,
double (*f)(const SPoint3 &p, void *), double (*f)(const SPoint3 &p, void *),
void *data, std::vector< IPT > & _result, void *data, std::vector< IPT > & _result,
double &dl, double epsilon = 1.e-6) double &dl, double epsilon = 1.e-6)
{ {
std::stack< IPT > _stack; std::stack< IPT > _stack;
_result.clear(); _result.clear();
// local parameters on the edge // local parameters on the edge
double t1 = 0.0; double t1 = 0.0;
double t2 = 1.0; double t2 = 1.0;
...@@ -98,7 +99,7 @@ double adaptiveTrapezoidalRule (SPoint3 p1 , SPoint3 p2 , ...@@ -98,7 +99,7 @@ double adaptiveTrapezoidalRule (SPoint3 p1 , SPoint3 p2 ,
SPoint3 pmid = p1 + dp*t12; SPoint3 pmid = p1 + dp*t12;
double fmid = 0.5* (f1+f2);//f(pmid,data); double fmid = 0.5* (f1+f2);//f(pmid,data);
double dt = t2-t1; double dt = t2-t1;
// average should be compared to mid value // average should be compared to mid value
double f12 = 0.5* (f1+f2); double f12 = 0.5* (f1+f2);
if (fabs (f12 - 0.5*(f1+f2)) > epsilon*dt ) { if (fabs (f12 - 0.5*(f1+f2)) > epsilon*dt ) {
IPT left (t1,t12,f1,f12); IPT left (t1,t12,f1,f12);
...@@ -109,7 +110,7 @@ double adaptiveTrapezoidalRule (SPoint3 p1 , SPoint3 p2 , ...@@ -109,7 +110,7 @@ double adaptiveTrapezoidalRule (SPoint3 p1 , SPoint3 p2 ,
else { else {
_result.push_back (pp); _result.push_back (pp);
// compute the integral using trapezoidal rule on both sides // compute the integral using trapezoidal rule on both sides
totalValue += 1./((0.5*f12+0.25*(f1+f2)))*dt; totalValue += 1./((0.5*f12+0.25*(f1+f2)))*dt;
} }
} }
// take into account the real length of the edge // take into account the real length of the edge
...@@ -119,16 +120,16 @@ double adaptiveTrapezoidalRule (SPoint3 p1 , SPoint3 p2 , ...@@ -119,16 +120,16 @@ double adaptiveTrapezoidalRule (SPoint3 p1 , SPoint3 p2 ,
} }
void saturateEdge (Edge &e, std::vector<Vertex*> &S, double (*f)(const SPoint3 &p, void *), void *data) { void saturateEdge (Edge &e, std::vector<Vertex*> &S, double (*f)(const SPoint3 &p, void *), void *data) {
std::vector< IPT > _result; std::vector< IPT > _result;
double dl; double dl;
SPoint3 p1 = e.first->point(); SPoint3 p1 = e.first->point();
SPoint3 p2 = e.second->point(); SPoint3 p2 = e.second->point();
double dN = adaptiveTrapezoidalRule (p1,p2,e.first->lc(), e.second->lc(), f,data,_result, dl); double dN = adaptiveTrapezoidalRule (p1,p2,e.first->lc(), e.second->lc(), f,data,_result, dl);
const int N = (int) (dN+0.1); const int N = (int) (dN+0.1);
double interval = dN/N; double interval = dN/N;
double L = 0.0; double L = 0.0;
// printf("edge length %g %d intervals of size %g\n",dl,N,interval); // printf("edge length %g %d intervals of size %g\n",dl,N,interval);
for (unsigned int i=0; i< _result.size() ; i++) { for (unsigned int i=0; i< _result.size() ; i++) {
...@@ -145,8 +146,8 @@ void saturateEdge (Edge &e, std::vector<Vertex*> &S, double (*f)(const SPoint3 & ...@@ -145,8 +146,8 @@ void saturateEdge (Edge &e, std::vector<Vertex*> &S, double (*f)(const SPoint3 &
break; break;
} }
else { else {
SPoint3 p = p1 * (1.-t) + p2*t; SPoint3 p = p1 * (1.-t) + p2*t;
double lc = e.first->lc() * (1.-t) + e.second->lc()*t; double lc = e.first->lc() * (1.-t) + e.second->lc()*t;
const double dx = 1.e-16 * (double) rand() / RAND_MAX; const double dx = 1.e-16 * (double) rand() / RAND_MAX;
const double dy = 1.e-16 * (double) rand() / RAND_MAX; const double dy = 1.e-16 * (double) rand() / RAND_MAX;
const double dz = 1.e-16 * (double) rand() / RAND_MAX; const double dz = 1.e-16 * (double) rand() / RAND_MAX;
...@@ -155,16 +156,16 @@ void saturateEdge (Edge &e, std::vector<Vertex*> &S, double (*f)(const SPoint3 & ...@@ -155,16 +156,16 @@ void saturateEdge (Edge &e, std::vector<Vertex*> &S, double (*f)(const SPoint3 &
} }
} }
} }
// printf("%d points added\n",S.size()); // printf("%d points added\n",S.size());
// exit(1); // exit(1);
} }
void saturateEdges ( edgeContainer &ec, void saturateEdges ( edgeContainer &ec,
std::vector<Tet*> &T, std::vector<Tet*> &T,
int nbThreads, int nbThreads,
std::vector<Vertex*> &S, std::vector<Vertex*> &S,
double (*f)(const SPoint3 &p, void *), void *data) { double (*f)(const SPoint3 &p, void *), void *data) {
AVGSEARCH= 0; AVGSEARCH= 0;
int counter = 0; int counter = 0;
...@@ -175,25 +176,25 @@ void saturateEdges ( edgeContainer &ec, ...@@ -175,25 +176,25 @@ void saturateEdges ( edgeContainer &ec,
Edge ed = T[i]->getEdge(iEdge); Edge ed = T[i]->getEdge(iEdge);
bool isNew = ec.addNewEdge(ed); bool isNew = ec.addNewEdge(ed);
counter++; counter++;
if (isNew){ if (isNew){
saturateEdge (ed, S, f, data); saturateEdge (ed, S, f, data);
} }
} }
} }
} }
// printf("a total of %d Tets counter %d AVG %d\n",T.size(),counter,AVGSEARCH/counter); // printf("a total of %d Tets counter %d AVG %d\n",T.size(),counter,AVGSEARCH/counter);
} }
///////////////////////// F I L T E R I N G //////////////////////////////////////////////////// ///////////////////////// F I L T E R I N G ////////////////////////////////////////////////////
class volumePointWithExclusionRegion { class volumePointWithExclusionRegion {
public : public :
Vertex *_v; Vertex *_v;
volumePointWithExclusionRegion (Vertex *v) : _v(v){ volumePointWithExclusionRegion (Vertex *v) : _v(v){
} }
bool inExclusionZone (volumePointWithExclusionRegion *p) bool inExclusionZone (volumePointWithExclusionRegion *p)
{ {
double d = sqrt ((p->_v->x() - _v->x()) * (p->_v->x() - _v->x())+ double d = sqrt ((p->_v->x() - _v->x()) * (p->_v->x() - _v->x())+
(p->_v->y() - _v->y()) * (p->_v->y() - _v->y())+ (p->_v->y() - _v->y()) * (p->_v->y() - _v->y())+
(p->_v->z() - _v->z()) * (p->_v->z() - _v->z())); (p->_v->z() - _v->z()) * (p->_v->z() - _v->z()));
...@@ -213,7 +214,7 @@ public : ...@@ -213,7 +214,7 @@ public :
struct my_wrapper_3D { struct my_wrapper_3D {
bool _tooclose; bool _tooclose;
volumePointWithExclusionRegion *_p; volumePointWithExclusionRegion *_p;
my_wrapper_3D (volumePointWithExclusionRegion *sp) : my_wrapper_3D (volumePointWithExclusionRegion *sp) :
_tooclose (false), _p(sp) {} _tooclose (false), _p(sp) {}
}; };
...@@ -231,7 +232,7 @@ bool rtree_callback(volumePointWithExclusionRegion *neighbour,void* point) ...@@ -231,7 +232,7 @@ bool rtree_callback(volumePointWithExclusionRegion *neighbour,void* point)
class vertexFilter { class vertexFilter {
RTree<volumePointWithExclusionRegion*,double,3,double> _rtree; RTree<volumePointWithExclusionRegion*,double,3,double> _rtree;
public: public:
void insert (Vertex * v) { void insert (Vertex * v) {
volumePointWithExclusionRegion *sp = new volumePointWithExclusionRegion (v); volumePointWithExclusionRegion *sp = new volumePointWithExclusionRegion (v);
double _min[3],_max[3]; double _min[3],_max[3];
...@@ -249,15 +250,15 @@ public: ...@@ -249,15 +250,15 @@ public:
} }
}; };
void filterVertices (const int numThreads, void filterVertices (const int numThreads,
vertexFilter &_filter, vertexFilter &_filter,
std::vector<Vertex*> &S, std::vector<Vertex*> &S,
std::vector<Vertex*> &add, std::vector<Vertex*> &add,
double (*f)(const SPoint3 &p, void *), double (*f)(const SPoint3 &p, void *),
void *data) { void *data) {
// printf("before : %d points to add\n",add.size()); // printf("before : %d points to add\n",add.size());
std::vector<Vertex*> _add = add; std::vector<Vertex*> _add = add;
for (unsigned int i=0;i<S.size();i++){ for (unsigned int i=0;i<S.size();i++){
SPoint3 p (S[i]->x(),S[i]->y(),S[i]->z()); SPoint3 p (S[i]->x(),S[i]->y(),S[i]->z());
double l = f (p, data); double l = f (p, data);
...@@ -273,7 +274,7 @@ void filterVertices (const int numThreads, ...@@ -273,7 +274,7 @@ void filterVertices (const int numThreads,
} }
else else
delete _add[i]; delete _add[i];
} }
// printf("after filter : %d points to add\n",add.size()); // printf("after filter : %d points to add\n",add.size());
} }
...@@ -310,8 +311,8 @@ void computeAdjacencies (Tet *t, int iFace, connSet &faceToTet){ ...@@ -310,8 +311,8 @@ void computeAdjacencies (Tet *t, int iFace, connSet &faceToTet){
} }
void edgeBasedRefinement (const int numThreads, void edgeBasedRefinement (const int numThreads,
const int nptsatonce, const int nptsatonce,
GRegion *gr) { GRegion *gr) {
// fill up old Datastructures // fill up old Datastructures
...@@ -388,40 +389,38 @@ void edgeBasedRefinement (const int numThreads, ...@@ -388,40 +389,38 @@ void edgeBasedRefinement (const int numThreads,
Msg::Info("----------------------------------- SATUR FILTR SORTH DELNY TIME TETS"); Msg::Info("----------------------------------- SATUR FILTR SORTH DELNY TIME TETS");
clock_t __t__ = clock(); double __t__ = Cpu();
while(1){ while(1){
char name[256]; char name[256];
// sprintf(name,"beforeRefinement%d.pos",iter); // sprintf(name,"beforeRefinement%d.pos",iter);
// print (name,_tets); // print (name,_tets);
// printf("ITER %3d\n",iter); // printf("ITER %3d\n",iter);
std::vector<Vertex*> add; std::vector<Vertex*> add;
clock_t t1 = clock(); double t1 = Cpu();
saturateEdges (ec, _tets, numThreads, add, _fx, NULL); saturateEdges (ec, _tets, numThreads, add, _fx, NULL);
// printf("%d points to add\n",add.size()); // printf("%d points to add\n",add.size());
//sprintf(name,"Points%d.pos",iter); //sprintf(name,"Points%d.pos",iter);
// _print (name,add); // _print (name,add);
clock_t t2 = clock(); double t2 = Cpu();
filterVertices (numThreads, _filter, _vertices, add, _fx, NULL); filterVertices (numThreads, _filter, _vertices, add, _fx, NULL);
clock_t t3 = clock(); double t3 = Cpu();
if (add.empty())break; if (add.empty())break;
std::vector<int> indices; std::vector<int> indices;
SortHilbert(add, indices); SortHilbert(add, indices);
clock_t t4 = clock(); double t4 = Cpu();
// sprintf(name,"PointsFiltered%d.pos",iter); // sprintf(name,"PointsFiltered%d.pos",iter);
// _print (name,add); // _print (name,add);
delaunayTrgl (1,1,add.size(), _tets, &add); delaunayTrgl (1,1,add.size(), _tets, &add);
clock_t t5 = clock(); double t5 = Cpu();
Msg::Info("IT %3d %6d points added, timings %5.2f %5.2f %5.2f %5.2f %5.2f %5d",iter,add.size(), Msg::Info("IT %3d %6d points added, timings %5.2f %5.2f %5.2f %5.2f %5.2f %5d",iter,add.size(),
(float)(t2-t1)/CLOCKS_PER_SEC, (t2-t1),
(float)(t3-t2)/CLOCKS_PER_SEC, (t3-t2),
(float)(t4-t3)/CLOCKS_PER_SEC, (t4-t3),
(float)(t5-t4)/CLOCKS_PER_SEC, (t5-t4),
(float)(t5-__t__)/CLOCKS_PER_SEC, (t5-__t__),
_tets.size()); _tets.size());
iter++; iter++;
} }
print ("afterRefinement.pos",_tets); print ("afterRefinement.pos",_tets);
} }
} }
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment