diff --git a/Mesh/delaunay3d.cpp b/Mesh/delaunay3d.cpp
index a89dfef0995a6251ffca179a16bda7e69777c206..ae37b17ba2578cc30cc4d61dfaa50c12a17a5c38 100644
--- a/Mesh/delaunay3d.cpp
+++ b/Mesh/delaunay3d.cpp
@@ -31,7 +31,6 @@
//int Tet::in_sphere_counter = 0;
-
struct HilbertSortB
{
// The code for generating table transgc from:
@@ -42,12 +41,12 @@ struct HilbertSortB
int Limit;
SBoundingBox3d bbox ;
void ComputeGrayCode(int n);
- int Split(Vertex** vertices,
+ int Split(Vert** vertices,
int arraysize,int GrayCode0,int GrayCode1,
double BoundingBoxXmin, double BoundingBoxXmax,
double BoundingBoxYmin, double BoundingBoxYmax,
double BoundingBoxZmin, double BoundingBoxZmax);
- void Sort(Vertex** vertices, int arraysize, int e, int d,
+ void Sort(Vert** vertices, int arraysize, int e, int d,
double BoundingBoxXmin, double BoundingBoxXmax,
double BoundingBoxYmin, double BoundingBoxYmax,
double BoundingBoxZmin, double BoundingBoxZmax, int depth);
@@ -55,7 +54,7 @@ struct HilbertSortB
{
ComputeGrayCode(3);
}
- void MultiscaleSortHilbert(Vertex** vertices, int arraysize,
+ void MultiscaleSortHilbert(Vert** vertices, int arraysize,
int threshold, double ratio, int *depth,
std::vector<int> &indices)
{
@@ -74,14 +73,14 @@ struct HilbertSortB
bbox.min().y(),bbox.max().y(),
bbox.min().z(),bbox.max().z(),0);
}
- void Apply (std::vector<Vertex*> &v, std::vector<int> &indices)
+ void Apply (std::vector<Vert*> &v, std::vector<int> &indices)
{
for (size_t i=0;i<v.size();i++){
- Vertex *pv = v[i];
+ Vert *pv = v[i];
bbox += SPoint3(pv->x(),pv->y(),pv->z());
}
bbox *= 1.01;
- Vertex**pv = &v[0];
+ Vert**pv = &v[0];
int depth;
indices.clear();
MultiscaleSortHilbert(pv, (int)v.size(), 64, .125,&depth,indices);
@@ -136,13 +135,13 @@ void HilbertSortB::ComputeGrayCode(int n)
}
-int HilbertSortB::Split(Vertex** vertices,
+int HilbertSortB::Split(Vert** vertices,
int arraysize,int GrayCode0,int GrayCode1,
double BoundingBoxXmin, double BoundingBoxXmax,
double BoundingBoxYmin, double BoundingBoxYmax,
double BoundingBoxZmin, double BoundingBoxZmax)
{
- Vertex* swapvert;
+ Vert* swapvert;
int axis, d;
double split;
int i, j;
@@ -210,7 +209,7 @@ int HilbertSortB::Split(Vertex** vertices,
// The sorting code is inspired by Tetgen 1.5
-void HilbertSortB::Sort(Vertex** vertices, int arraysize, int e, int d,
+void HilbertSortB::Sort(Vert** vertices, int arraysize, int e, int d,
double BoundingBoxXmin, double BoundingBoxXmax,
double BoundingBoxYmin, double BoundingBoxYmax,
double BoundingBoxZmin, double BoundingBoxZmax, int depth)
@@ -301,7 +300,7 @@ void HilbertSortB::Sort(Vertex** vertices, int arraysize, int e, int d,
}
}
-void SortHilbert (std::vector<Vertex*>& v, std::vector<int> &indices)
+void SortHilbert (std::vector<Vert*>& v, std::vector<int> &indices)
{
HilbertSortB h(1000);
h.Apply(v, indices);
@@ -354,10 +353,10 @@ static void computeNeighboringTetsOfACavity(const std::vector<Tet*> &cavity,
}
static bool buildEdgeCavity(Tet *t, int iLocalEdge,
- Vertex **v1, Vertex **v2,
+ Vert **v1, Vert **v2,
std::vector<Tet*> &cavity,
std::vector<Tet*> &outside,
- std::vector<Vertex*> &ring)
+ std::vector<Vert*> &ring)
{
cavity.clear();
ring.clear();
@@ -366,13 +365,13 @@ static bool buildEdgeCavity(Tet *t, int iLocalEdge,
*v2 = t->V[edges[iLocalEdge][1]];
// the 5 - i th edge contains the other 2 points of the tet
- Vertex *lastinring = t->V[edges[5 - iLocalEdge][0]];
+ Vert *lastinring = t->V[edges[5 - iLocalEdge][0]];
ring.push_back(lastinring);
cavity.push_back(t);
while (1){
- Vertex *ov1 = t->V[edges[5 - iLocalEdge][0]];
- Vertex *ov2 = t->V[edges[5 - iLocalEdge][1]];
+ Vert *ov1 = t->V[edges[5 - iLocalEdge][0]];
+ Vert *ov2 = t->V[edges[5 - iLocalEdge][1]];
int K = ov1 == lastinring ? 1 : 0;
lastinring = ov1 == lastinring ? ov2 : ov1;
// look in the 2 faces sharing this edge the one that has vertex
@@ -395,8 +394,8 @@ static bool buildEdgeCavity(Tet *t, int iLocalEdge,
cavity.push_back(t);
iLocalEdge = -1;
for (int i = 0; i < 6; i++){
- Vertex *a = t->V[edges[i][0]];
- Vertex *b = t->V[edges[i][1]];
+ Vert *a = t->V[edges[i][0]];
+ Vert *b = t->V[edges[i][1]];
if ((a == *v1 && b == *v2) || (a == *v2 && b == *v1)){
iLocalEdge = i;
break;
@@ -417,8 +416,8 @@ bool edgeSwap(Tet *tet, int iLocalEdge, tetContainer &T, int myThread)
{
std::vector<Tet*> cavity;
std::vector<Tet*> outside;
- std::vector<Vertex*> ring;
- Vertex *v1, *v2;
+ std::vector<Vert*> ring;
+ Vert *v1, *v2;
bool closed = buildEdgeCavity(tet, iLocalEdge, &v1, &v2, cavity, outside, ring);
@@ -436,10 +435,10 @@ bool edgeSwap(Tet *tet, int iLocalEdge, tetContainer &T, int myThread)
double tetQualityRef = 1.e22;
for (unsigned int i = 0; i < cavity.size(); i++){
- Vertex *vx0 = cavity[i]->V[0];
- Vertex *vx1 = cavity[i]->V[1];
- Vertex *vx2 = cavity[i]->V[2];
- Vertex *vx3 = cavity[i]->V[3];
+ Vert *vx0 = cavity[i]->V[0];
+ Vert *vx1 = cavity[i]->V[1];
+ Vert *vx2 = cavity[i]->V[2];
+ Vert *vx3 = cavity[i]->V[3];
double volume = robustPredicates::orient3d((double*)vx0, (double*)vx1, (double*)vx2, (double*)vx3);
const double a =
SQR (vx0->x() - vx1->x()) + SQR (vx0->y() - vx1->y()) + SQR (vx0->z() - vx1->z()) +
@@ -504,9 +503,9 @@ bool edgeSwap(Tet *tet, int iLocalEdge, tetContainer &T, int myThread)
int p1 = sp.triangles[iT][0];
int p2 = sp.triangles[iT][1];
int p3 = sp.triangles[iT][2];
- Vertex *pv1 = ring[p1];
- Vertex *pv2 = ring[p2];
- Vertex *pv3 = ring[p3];
+ Vert *pv1 = ring[p1];
+ Vert *pv2 = ring[p2];
+ Vert *pv3 = ring[p3];
Tet *nt1,*nt2;
if (counter < cavity.size()) nt1 = cavity[counter++];
@@ -544,7 +543,7 @@ static void edgeSwapPass (tetContainer &T)
Finite Elements in Analysis and Design 25 (1997) 297-317
*/
-static void starShapeness (Vertex *v, connContainer &bndK,
+static void starShapeness (Vert *v, connContainer &bndK,
std::vector<unsigned int> &_negatives,
double threshold)
{
@@ -561,7 +560,7 @@ static void starShapeness (Vertex *v, connContainer &bndK,
}
}
-static Tet* tetContainsV (Vertex *v, cavityContainer &cavity)
+static Tet* tetContainsV (Vert *v, cavityContainer &cavity)
{
for (unsigned int i=0; i< cavity.size(); i++) {
unsigned int count = 0;
@@ -681,7 +680,7 @@ static Tet *tetInsideCavityWithFAce (Face &f, cavityContainer &cavity)
}
static bool fixDelaunayCavity (double threshold,
- Vertex *v,
+ Vert *v,
cavityContainer &cavity,
connContainer &bndK,
int myThread, int K,
@@ -723,7 +722,7 @@ static bool fixDelaunayCavity (double threshold,
static void delaunayCavity2 (Tet *t,
Tet *prev,
- Vertex *v,
+ Vert *v,
cavityContainer &cavity,
connContainer &bnd,
int thread, int iPnt)
@@ -747,7 +746,7 @@ static void delaunayCavity2 (Tet *t,
}
}
-Tet* walk (Tet *t, Vertex *v, int maxx, double &totSearch, int thread)
+Tet* walk (Tet *t, Vert *v, int maxx, double &totSearch, int thread)
{
while (1){
totSearch++;
@@ -781,7 +780,7 @@ Tet* walk (Tet *t, Vertex *v, int maxx, double &totSearch, int thread)
Msg::Fatal("Jump-and-Walk Failed (No neighbor)");
}
-void __print (const char *name, connContainer &conn, Vertex *v)
+void __print (const char *name, connContainer &conn, Vert *v)
{
FILE *f = fopen(name,"w");
fprintf(f,"View \"\"{\n");
@@ -798,7 +797,7 @@ void __print (const char *name, connContainer &conn, Vertex *v)
fclose(f);
}
-void __print (const char *name, int thread, tetContainer &T, Vertex *v){
+void __print (const char *name, int thread, tetContainer &T, Vert *v){
FILE *f = fopen(name,"w");
fprintf(f,"View \"\"{\n");
if (v)fprintf(f,"SP(%g,%g,%g){%d};\n",v->x(),v->y(),v->z(),v->getNum());
@@ -829,9 +828,9 @@ void __print (const char *name, int thread, tetContainer &T, Vertex *v){
fclose(f);
}
-void print (std::vector<Vertex*> &V, std::vector<Tet*> &T)
+void print (std::vector<Vert*> &V, std::vector<Tet*> &T)
{
- std::map<Vertex*,int> nums;
+ std::map<Vert*,int> nums;
for (unsigned int i=0;i<V.size();i++){
nums[V[i]] = i;
}
@@ -842,7 +841,7 @@ void print (std::vector<Vertex*> &V, std::vector<Tet*> &T)
}
}
-void print (const char *name, std::vector<Vertex*> &T)
+void print (const char *name, std::vector<Vert*> &T)
{
FILE *f = fopen(name,"w");
fprintf(f,"View \"\"{\n");
@@ -910,7 +909,7 @@ static Tet* randomTet (int thread, tetContainer &allocator)
void delaunayTrgl (const unsigned int numThreads,
const unsigned int NPTS_AT_ONCE,
unsigned int Npts,
- std::vector<Vertex*> assignTo[],
+ std::vector<Vert*> assignTo[],
tetContainer &allocator,
double threshold)
{
@@ -957,15 +956,15 @@ void delaunayTrgl (const unsigned int numThreads,
invalidCavities [myThread] = 0;
unsigned int locSize=0;
std::vector<unsigned int> locSizeK(NPTS_AT_ONCE);
- std::vector<Vertex*> allocatedVerts(NPTS_AT_ONCE);
+ std::vector<Vert*> allocatedVerts(NPTS_AT_ONCE);
for (unsigned int K=0;K<NPTS_AT_ONCE;K++){
locSizeK[K] = assignTo[K+myThread*NPTS_AT_ONCE].size();
locSize += locSizeK[K];
#if defined(_HAVE_NUMA)
- allocatedVerts [K] = (Vertex*)numa_alloc_local (locSizeK[K]*sizeof(Vertex));
+ allocatedVerts [K] = (Vert*)numa_alloc_local (locSizeK[K]*sizeof(Vert));
#else
- // allocatedVerts [K] = (Vertex*)calloc (locSizeK[K],sizeof(Vertex));
- allocatedVerts [K] = new Vertex [locSizeK[K]];
+ // allocatedVerts [K] = (Vert*)calloc (locSizeK[K],sizeof(Vert));
+ allocatedVerts [K] = new Vert [locSizeK[K]];
#endif
for (unsigned int iP=0 ; iP < locSizeK[K] ; iP++){
allocatedVerts[K][iP] = *(assignTo[K+myThread*NPTS_AT_ONCE][iP]);
@@ -973,7 +972,7 @@ void delaunayTrgl (const unsigned int numThreads,
}
}
- std::vector<Vertex*> vToAdd(NPTS_AT_ONCE);
+ std::vector<Vert*> vToAdd(NPTS_AT_ONCE);
#if defined(_OPENMP)
#pragma omp barrier
@@ -1110,26 +1109,26 @@ void delaunayTrgl (const unsigned int numThreads,
}
-static void initialCube (std::vector<Vertex*> &v,
- Vertex *box[8],
+static void initialCube (std::vector<Vert*> &v,
+ Vert *box[8],
tetContainer & allocator)
{
SBoundingBox3d bbox ;
// bbox += SPoint3(0,0,0);
// bbox += SPoint3(1,1,1);
for (size_t i=0;i<v.size();i++){
- Vertex *pv = v[i];
+ Vert *pv = v[i];
bbox += SPoint3(pv->x(),pv->y(),pv->z());
}
bbox *= 1.3;
- box[0] = new Vertex (bbox.min().x(),bbox.min().y(),bbox.min().z(),bbox.diag());
- box[1] = new Vertex (bbox.max().x(),bbox.min().y(),bbox.min().z(),bbox.diag());
- box[2] = new Vertex (bbox.max().x(),bbox.max().y(),bbox.min().z(),bbox.diag());
- box[3] = new Vertex (bbox.min().x(),bbox.max().y(),bbox.min().z(),bbox.diag());
- box[4] = new Vertex (bbox.min().x(),bbox.min().y(),bbox.max().z(),bbox.diag());
- box[5] = new Vertex (bbox.max().x(),bbox.min().y(),bbox.max().z(),bbox.diag());
- box[6] = new Vertex (bbox.max().x(),bbox.max().y(),bbox.max().z(),bbox.diag());
- box[7] = new Vertex (bbox.min().x(),bbox.max().y(),bbox.max().z(),bbox.diag());
+ box[0] = new Vert (bbox.min().x(),bbox.min().y(),bbox.min().z(),bbox.diag());
+ box[1] = new Vert (bbox.max().x(),bbox.min().y(),bbox.min().z(),bbox.diag());
+ box[2] = new Vert (bbox.max().x(),bbox.max().y(),bbox.min().z(),bbox.diag());
+ box[3] = new Vert (bbox.min().x(),bbox.max().y(),bbox.min().z(),bbox.diag());
+ box[4] = new Vert (bbox.min().x(),bbox.min().y(),bbox.max().z(),bbox.diag());
+ box[5] = new Vert (bbox.max().x(),bbox.min().y(),bbox.max().z(),bbox.diag());
+ box[6] = new Vert (bbox.max().x(),bbox.max().y(),bbox.max().z(),bbox.diag());
+ box[7] = new Vert (bbox.min().x(),bbox.max().y(),bbox.max().z(),bbox.diag());
/* Tet *t0 = new Tet (box[2],box[7],box[3],box[1]);
Tet *t1 = new Tet (box[0],box[7],box[1],box[3]);
@@ -1159,8 +1158,8 @@ static void initialCube (std::vector<Vertex*> &v,
void delaunayTriangulation (const int numThreads,
const int nptsatonce,
- std::vector<Vertex*> &S,
- Vertex *box[8],
+ std::vector<Vert*> &S,
+ Vert *box[8],
tetContainer & allocator)
{
int N = S.size();
@@ -1175,8 +1174,8 @@ void delaunayTriangulation (const int numThreads,
// int blockSize = (nbBlocks * (N / nbBlocks))/nbBlocks;
- std::vector<Vertex*> assignTo0[1];
- std::vector<std::vector<Vertex*> > assignTo(nbBlocks);
+ std::vector<Vert*> assignTo0[1];
+ std::vector<std::vector<Vert*> > assignTo(nbBlocks);
for (unsigned int i=1;i<indices.size();i++){
int start = indices[i-1];
@@ -1216,7 +1215,7 @@ void delaunayTriangulation (const int numThreads,
std::vector<MTetrahedron*> &T)
{
std::vector<MVertex*> _temp;
- std::vector<Vertex*> _vertices;
+ std::vector<Vert*> _vertices;
unsigned int N = S.size();
_temp.resize(N+1+8);
double maxx=0, maxy=0,maxz=0;
@@ -1239,8 +1238,8 @@ 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);
+ // Vert *v = new Vert (mv->x(),mv->y(),mv->z(),1.e22,i+1);
+ Vert *v = new Vert (mv->x(),mv->y(),mv->z(),1.e22,i+1);
_vertices.push_back(v);
_temp [v->getNum()] = mv;
}
@@ -1249,13 +1248,13 @@ void delaunayTriangulation (const int numThreads,
// FIXME numThreads
- Vertex *box[8];
+ Vert *box[8];
delaunayTriangulation (numThreads, nptsatonce, _vertices, box, allocator);
//__print ("finalTetrahedrization.pos",0, allocator);
MVertex *VV[8];
for (int i=0;i<8;i++){
- Vertex *v = box[i];
+ Vert *v = box[i];
v->setNum(N+i+1);
VV[i]= new MVertex (v->x(),v->y(),v->z(),NULL,N+(i+1));
_temp [v->getNum()] = VV[i];
diff --git a/Mesh/delaunay3d.h b/Mesh/delaunay3d.h
index e20921b739cba0bc79d20d5a88941a38c13f20b1..58c578d4e49996e9e16093b7de7e1bed46529628 100644
--- a/Mesh/delaunay3d.h
+++ b/Mesh/delaunay3d.h
@@ -8,9 +8,10 @@
class MVertex;
class MTetrahedron;
-void delaunayTriangulation (const int numThreads,
- const int nptsatonce,
- std::vector<MVertex*> &S,
- std::vector<MTetrahedron*> &T);
+
+void delaunayTriangulation(const int numThreads,
+ const int nptsatonce,
+ std::vector<MVertex*> &S,
+ std::vector<MTetrahedron*> &T);
#endif
diff --git a/Mesh/delaunay3d_private.h b/Mesh/delaunay3d_private.h
index 077b66be8492f3e9d53fdf901965068939122ce1..bfab084c0743dd343bb735cb1a2ea3d61ae7f906 100644
--- a/Mesh/delaunay3d_private.h
+++ b/Mesh/delaunay3d_private.h
@@ -21,8 +21,7 @@
typedef unsigned char CHECKTYPE;
-
-struct Vertex {
+struct Vert {
private :
double _x[3];
double _lc;
@@ -40,18 +39,22 @@ public :
inline double &z() {return _x[2];}
inline double &lc() {return _lc;}
inline operator double *() { return _x; }
- 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);
+ Vert (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 double &other){
- return Vertex (x()*other,y()*other,z()*other,_lc*other);
+ Vert operator + (const Vert &other)
+ {
+ return Vert (x()+other.x(),y()+other.y(),z()+other.z(),other.lc() + _lc);
+ }
+ Vert operator * (const double &other)
+ {
+ return Vert(x()*other, y()*other, z()*other, _lc*other);
}
inline SPoint3 point() const { return SPoint3(x(), y(), z()); }
};
-
inline double orientationTestFast(double *pa, double *pb, double *pc, double *pd)
{
const double adx = pa[0] - pd[0];
@@ -69,13 +72,16 @@ inline double orientationTestFast(double *pa, double *pb, double *pc, double *pd
+ cdx * (ady * bdz - adz * bdy);
}
-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);
+inline bool inSphereTest_s (Vert *va, Vert *vb, Vert *vc, Vert *vd , Vert *ve)
+{
+ double val = robustPredicates::insphere((double*)va, (double*)vb, (double*)vc,
+ (double*)vd, (double*)ve);
if (val == 0.0){
- Msg::Info("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};
+ Vert *pt[5] = {va,vb,vc,vd,ve};
int swaps = 0;
int n = 5;
do {
@@ -83,20 +89,22 @@ inline bool inSphereTest_s (Vertex *va, Vertex *vb, Vertex *vc, Vertex *vd , Ver
n = n - 1;
for (int i = 0; i < n; i++) {
if (pt[i] > pt[i+1]) {
- Vertex *swappt = pt[i]; pt[i] = pt[i+1]; pt[i+1] = swappt;
+ Vert *swappt = pt[i]; pt[i] = pt[i+1]; pt[i+1] = swappt;
count++;
}
}
swaps += count;
} while (count > 0);
- double oriA = robustPredicates::orient3d ((double*)pt[1], (double*)pt[2], (double*)pt[3], (double*)pt[4]);
+ double oriA = robustPredicates::orient3d((double*)pt[1], (double*)pt[2],
+ (double*)pt[3], (double*)pt[4]);
if (oriA != 0.0) {
// Flip the sign if there are odd number of swaps.
if ((swaps % 2) != 0) oriA = -oriA;
val = oriA;
}
else {
- double oriB = -robustPredicates::orient3d ((double*)pt[0], (double*)pt[2], (double*)pt[3], (double*)pt[4]);
+ double oriB = -robustPredicates::orient3d((double*)pt[0], (double*)pt[2],
+ (double*)pt[3], (double*)pt[4]);
if (oriB == 0.0) {
Msg::Fatal("Symbolic perturbation failed in icCircle Predicate");
}
@@ -108,18 +116,19 @@ inline bool inSphereTest_s (Vertex *va, Vertex *vb, Vertex *vc, Vertex *vd , Ver
return val > 0 ? 1 : 0;
}
-inline double orientationTest (Vertex *va, Vertex *vb, Vertex *vc, Vertex *vd){
+inline double orientationTest (Vert *va, Vert *vb, Vert *vc, Vert *vd)
+{
return robustPredicates::orient3d ((double*)va,(double*)vb,(double*)vc,(double*)vd);
}
-inline double orientationTestFast (Vertex *va, Vertex *vb, Vertex *vc, Vertex *vd){
+inline double orientationTestFast (Vert *va, Vert *vb, Vert *vc, Vert *vd)
+{
return orientationTestFast ((double*)va,(double*)vb,(double*)vc,(double*)vd);
}
-
struct Edge {
- Vertex *first,*second;
- Edge (Vertex *v1, Vertex *v2) : first(v1), second(v2)
+ Vert *first, *second;
+ Edge (Vert *v1, Vert *v2) : first(v1), second(v2)
{
}
bool operator == (const Edge &e) const{
@@ -133,24 +142,25 @@ struct Edge {
}
};
-//typedef std::pair<Vertex*, Vertex*> Edge;
-
struct Face {
- Vertex *v[3];
- Vertex *V[3];
- Face (Vertex *v1, Vertex *v2, Vertex *v3){
+ Vert *v[3];
+ Vert *V[3];
+ Face (Vert *v1, Vert *v2, Vert *v3)
+ {
V[0] = v[0] = v1;
V[1] = v[1] = v2;
V[2] = v[2] = v3;
-#define cswap(a,b) do { if(a > b) { Vertex* tmp = a; a = b; b = tmp; } } while(0)
+#define cswap(a,b) do { if(a > b) { Vert* tmp = a; a = b; b = tmp; } } while(0)
cswap(v[0], v[1]);
cswap(v[1], v[2]);
cswap(v[0], v[1]);
}
- bool operator == (const Face &other) const {
+ bool operator == (const Face &other) const
+ {
return v[0] == other.v[0] && v[1] == other.v[1] && v[2] == other.v[2];
}
- bool operator < (const Face &other) const {
+ bool operator < (const Face &other) const
+ {
if (v[0] < other.v[0])return true;
if (v[0] > other.v[0])return false;
if (v[1] < other.v[1])return true;
@@ -161,9 +171,9 @@ struct Face {
};
struct Tet {
- Tet *T[4];
- Vertex *V[4];
- CHECKTYPE _bitset [MAX_NUM_THREADS_];
+ Tet *T[4];
+ Vert *V[4];
+ CHECKTYPE _bitset[MAX_NUM_THREADS_];
bool _modified;
// static int in_sphere_counter;
Tet () : _modified(true){
@@ -171,66 +181,76 @@ struct Tet {
T[0] = T[1] = T[2] = T[3] = NULL;
setAllDeleted();
}
- int setVerticesNoTest (Vertex *v0, Vertex *v1, Vertex *v2, Vertex *v3){
+ int setVerticesNoTest (Vert *v0, Vert *v1, Vert *v2, Vert *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){
+ int setVertices (Vert *v0, Vert *v1, Vert *v2, Vert *v3)
+ {
_modified=true;
- double val = robustPredicates::orient3d((double*)v0,(double*)v1,(double*)v2,(double*)v3);
+ 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];
+ // for (int i=0;i<4;i++)_copy[i] = *V[i];
return 1;
}
else if (val < 0){
- // throw;
+ // throw;
V[0] = v1; V[1] = v0; V[2] = v2; V[3] = v3;
- // for (int i=0;i<4;i++)_copy[i] = *V[i];
+ // for (int i=0;i<4;i++)_copy[i] = *V[i];
return -1;
}
else {
- // throw;
+ // throw;
return 0;
}
}
- Tet (Vertex *v0, Vertex *v1, Vertex *v2, Vertex *v3)
+ Tet(Vert *v0, Vert *v1, Vert *v2, Vert *v3)
{
setVertices (v0,v1,v2,v3);
T[0] = T[1] = T[2] = T[3] = NULL;
setAllDeleted();
}
- void setAllDeleted (){
- for (int i=0;i<MAX_NUM_THREADS_;i++) _bitset [i] = 0x0;
+ void setAllDeleted()
+ {
+ for (int i = 0; i < MAX_NUM_THREADS_; i++) _bitset [i] = 0x0;
}
- inline void unset ( int thread, int iPnt ){
- _bitset[thread] &= ~(1 << iPnt);
+ inline void unset(int thread, int iPnt)
+ {
+ _bitset[thread] &= ~(1 << iPnt);
}
- inline void set ( int thread, int iPnt ){
- _bitset [thread] |= (1 << iPnt);
+ inline void set(int thread, int iPnt)
+ {
+ _bitset [thread] |= (1 << iPnt);
}
- inline CHECKTYPE isSet ( int thread, int iPnt ) const{
+ inline CHECKTYPE isSet(int thread, int iPnt) const
+ {
return _bitset[thread] & (1 << iPnt);
}
- inline Face getFace (int k) const {
+ inline Face getFace(int k) const
+ {
const int fac[4][3] = {{0,1,2},{1,3,2},{2,3,0},{1,0,3}};
- return Face (V[fac[k][0]],V[fac[k][1]],V[fac[k][2]]);
+ return Face(V[fac[k][0]], V[fac[k][1]], V[fac[k][2]]);
}
- inline Vertex* getOppositeVertex (int k) const {
+ inline Vert* getOppositeVertex(int k) const
+ {
const int o[4] = {3,0,1,2};
return V[o[k]];
}
-
- inline Edge getEdge (int k) const {
+ inline Edge getEdge (int k) const
+ {
const int edg[6][2] = {{0,1},{0,2},{0,3},{1,2},{1,3},{2,3}};
- return Edge (std::min(V[edg[k][0]],V[edg[k][1]]),
- std::max(V[edg[k][0]],V[edg[k][1]]));
+ 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) {
+ inline bool inSphere (Vert *vd, int thread)
+ {
// in_sphere_counter++;
- return inSphereTest_s (V[0],V[1],V[2],V[3],vd);
+ return inSphereTest_s(V[0], V[1], V[2], V[3], vd);
}
};
@@ -240,11 +260,14 @@ struct conn {
Tet *t;
conn () : f(0,0,0), i(0), t(0){}
conn (Face _f, int _i, Tet *_t) : f(_f), i(_i), t(_t)
- { }
- inline bool operator == (const conn & c) const{
+ {
+ }
+ inline bool operator == (const conn & c) const
+ {
return f == c.f;
}
- inline bool operator < (const conn & c) const{
+ inline bool operator < (const conn & c) const
+ {
return f < c.f;
}
};
@@ -256,27 +279,32 @@ public:
std::vector<T*> _all;
unsigned int _current;
unsigned int _nbAlloc;
- unsigned int size () {
+ unsigned int size ()
+ {
return _current + (_all.size()-1) * _nbAlloc;
}
- inline T* operator () (unsigned int i){
+ inline T* operator () (unsigned int i)
+ {
const unsigned int _array = i / _nbAlloc;
const unsigned int _offset = i % _nbAlloc;
- // printf("%d %d %d\n",i,_array,_offset);
+ // printf("%d %d %d\n",i,_array,_offset);
return _all [_array] + _offset;
}
- aBunchOfStuff (unsigned int s) : _current(0), _nbAlloc (s) {
+ aBunchOfStuff (unsigned int s) : _current(0), _nbAlloc (s)
+ {
_all.push_back(new T [_nbAlloc]);
}
- ~aBunchOfStuff(){
+ ~aBunchOfStuff()
+ {
for (unsigned int i=0;i<_all.size();i++){
delete [] _all[i];
}
}
- T* newStuff() {
+ T* newStuff()
+ {
if (_current == _nbAlloc){
_all.push_back(new T [_nbAlloc]);
- // printf("REALLOCATION %d\n",_all.size());
+ // printf("REALLOCATION %d\n",_all.size());
_current = 0;
}
_current++;
@@ -288,9 +316,13 @@ public:
class tetContainer {
std::vector<aBunchOfStuff<Tet> *> _perThread;
public:
- unsigned int size(int thread) const {return _perThread[thread]->size();}
- inline Tet * operator () (int thread, int j) const {return (*_perThread[thread]) (j);}
- tetContainer (int nbThreads, int preallocSizePerThread) {
+ unsigned int size(int thread) const { return _perThread[thread]->size(); }
+ inline Tet * operator () (int thread, int j) const
+ {
+ return (*_perThread[thread])(j);
+ }
+ tetContainer (int nbThreads, int preallocSizePerThread)
+ {
// FIXME !!!
if (nbThreads != 1) throw;
_perThread.resize(nbThreads);
@@ -303,28 +335,30 @@ class tetContainer {
#else
int myThread = 0;
#endif
- _perThread [myThread] = new aBunchOfStuff<Tet> (preallocSizePerThread) ;
+ _perThread [myThread] = new aBunchOfStuff<Tet>(preallocSizePerThread) ;
}
}
- inline Tet * newTet(int thread) {
- return _perThread[thread]->newStuff();
+ inline Tet * newTet(int thread)
+ {
+ return _perThread[thread]->newStuff();
}
- ~tetContainer(){
+ ~tetContainer()
+ {
delete _perThread [0];
}
};
typedef std::vector<Tet*> cavityContainer;
-typedef std::vector<conn> connContainer;
+typedef std::vector<conn> connContainer;
-void SortHilbert (std::vector<Vertex*>& v, std::vector<int> &indices);
+void SortHilbert(std::vector<Vert*>& v, std::vector<int> &indices);
void computeAdjacencies (Tet *t, int iFace, connContainer &faceToTet);
-void __print (const char *name, int thread, tetContainer &T, Vertex *v = 0);
-void delaunayTrgl (const unsigned int numThreads,
- const unsigned int NPTS_AT_ONCE,
- unsigned int Npts,
- std::vector<Vertex*> assignTo[],
- tetContainer &allocator, double threshold = 0.0);
+void __print(const char *name, int thread, tetContainer &T, Vert *v = 0);
+void delaunayTrgl(const unsigned int numThreads,
+ const unsigned int NPTS_AT_ONCE,
+ unsigned int Npts,
+ std::vector<Vert*> assignTo[],
+ tetContainer &allocator, double threshold = 0.0);
bool edgeSwap(Tet *tet, int iLocalEdge, tetContainer &T, int myThread);
#endif
diff --git a/Mesh/delaunay_refinement.cpp b/Mesh/delaunay_refinement.cpp
index 4d53dd74d1eb005fd31f2791ab561c909b983e17..a9c259e0a4223b72c08d37fea53898acdb4b075b 100644
--- a/Mesh/delaunay_refinement.cpp
+++ b/Mesh/delaunay_refinement.cpp
@@ -32,10 +32,12 @@ struct edgeContainer
{
std::set< Edge > _hash2;
std::vector<std::vector<Edge> > _hash;
- edgeContainer (unsigned int N = 1000000) {
+ edgeContainer(unsigned int N = 1000000)
+ {
_hash.resize(N);
}
- bool addNewEdge2 (const Edge &e) {
+ bool addNewEdge2 (const Edge &e)
+ {
std::set< Edge >::iterator it = _hash2.find(e);
if (it != _hash2.end())return false;
_hash2.insert(e);
@@ -54,15 +56,18 @@ struct edgeContainer
struct IPT {
double _x1,_x2,_x3,_x4;
- IPT(double x1, double x2, double x3, double x4) :
- _x1(x1),_x2(x2),_x3(x3),_x4(x4){};
+ IPT(double x1, double x2, double x3, double x4)
+ : _x1(x1), _x2(x2), _x3(x3), _x4(x4)
+ {
+ };
};
-double adaptiveTrapezoidalRule (SPoint3 p1 , SPoint3 p2 ,
- double lc1 , double lc2 ,
- double (*f)(const SPoint3 &p, void *),
- void *data, std::vector< IPT > & _result,
- double &dl, std::stack<IPT> &_stack, double epsilon = 1.e-5)
+double adaptiveTrapezoidalRule(SPoint3 p1 , SPoint3 p2 ,
+ double lc1 , double lc2 ,
+ double (*f)(const SPoint3 &p, void *),
+ void *data, std::vector< IPT > & _result,
+ double &dl, std::stack<IPT> &_stack,
+ double epsilon = 1.e-5)
{
// _stack.clear();
_result.clear();
@@ -117,20 +122,22 @@ double adaptiveTrapezoidalRule (SPoint3 p1 , SPoint3 p2 ,
}
-void saturateEdge (Edge &e, std::vector<Vertex*> &S,
- double (*f)(const SPoint3 &p, void *),
- void *data, std::stack<IPT> &temp)
+void saturateEdge(Edge &e, std::vector<Vert*> &S,
+ double (*f)(const SPoint3 &p, void *),
+ void *data, std::stack<IPT> &temp)
{
std::vector< IPT > _result;
double dl;
SPoint3 p1 = e.first->point();
SPoint3 p2 = e.second->point();
- const double dN = adaptiveTrapezoidalRule (p1,p2,e.first->lc(), e.second->lc(), f,data,_result, dl, temp);
+ const double dN = adaptiveTrapezoidalRule(p1,p2,e.first->lc(), e.second->lc(),
+ f, data, _result, dl, temp);
const int N = (int) (dN+0.1);
const double interval = dN/N;
double L = 0.0;
- // printf("edge length %g %d intervals of size %g (%d results)\n",dl,N,interval,_result.size());
+ // printf("edge length %g %d intervals of size %g (%d results)\n",
+ // dl,N,interval,_result.size());
const unsigned int Nr = _result.size();
for (unsigned int i=0; i< Nr ; i++) {
const IPT & rr = _result[i];
@@ -140,7 +147,7 @@ void saturateEdge (Edge &e, std::vector<Vertex*> &S,
const double f2 = rr._x4;
const double dL = 2.*(t2-t1) * dl / (f1+f2);
- // printf("%g --> %g for %g --> %g\n",L,dL,t1,t2);
+ // printf("%g --> %g for %g --> %g\n",L,dL,t1,t2);
double L0 = L;
while (1) {
const double t = t1 + (L+interval-L0)*(t2-t1) / dL;
@@ -148,13 +155,13 @@ void saturateEdge (Edge &e, std::vector<Vertex*> &S,
break;
}
else {
- // printf("%g ",t);
+ // printf("%g ",t);
SPoint3 p = p1 * (1.-t) + p2*t;
double lc = e.first->lc() * (1.-t) + e.second->lc()*t;
const double dx = 0;//1.e-12 * (double) rand() / RAND_MAX;
const double dy = 0;//1.e-12 * (double) rand() / RAND_MAX;
const double dz = 0;//1.e-12 * (double) rand() / RAND_MAX;
- S.push_back(new Vertex(p.x()+dx,p.y()+dy,p.z()+dz,lc));
+ S.push_back(new Vert(p.x()+dx,p.y()+dy,p.z()+dz,lc));
L += interval;
}
}
@@ -169,7 +176,7 @@ void saturateEdge (Edge &e, std::vector<Vertex*> &S,
void saturateEdges(edgeContainer &ec,
tetContainer &T,
int nbThreads,
- std::vector<Vertex*> &S,
+ std::vector<Vert*> &S,
double (*f)(const SPoint3 &p, void *), void *data)
{
std::stack<IPT> temp;
@@ -191,14 +198,14 @@ void saturateEdges(edgeContainer &ec,
}
}
-///////////////////////// F I L T E R I N G ////////////////////////////////////////////////////
+// F I L T E R I N G
#define SQR(X) (X)*(X)
class volumePointWithExclusionRegion {
public :
- Vertex *_v;
- volumePointWithExclusionRegion (Vertex *v) : _v(v) {}
+ Vert *_v;
+ volumePointWithExclusionRegion (Vert *v) : _v(v) {}
inline bool inExclusionZone (volumePointWithExclusionRegion *p) const
{
@@ -243,7 +250,7 @@ bool rtree_callback(volumePointWithExclusionRegion *neighbour,void* point)
class vertexFilter {
RTree<volumePointWithExclusionRegion*,double,3,double> _rtree;
public:
- void insert (Vertex * v) {
+ void insert (Vert * v) {
volumePointWithExclusionRegion *sp = new volumePointWithExclusionRegion (v);
double _min[3],_max[3];
sp->minmax(_min,_max);
@@ -260,18 +267,18 @@ public:
}
};
-void filterVertices (const int numThreads,
- vertexFilter &_filter,
- std::vector<Vertex*> &add,
- double (*f)(const SPoint3 &p, void *),
- void *data)
+void filterVertices(const int numThreads,
+ vertexFilter &_filter,
+ std::vector<Vert*> &add,
+ double (*f)(const SPoint3 &p, void *),
+ void *data)
{
std::vector<int> indices;
SortHilbert(add, indices);
- std::vector<Vertex*> _add=add;
+ std::vector<Vert*> _add=add;
- // std::vector<Vertex*> _add;
- // Vertex *current = add[0];
+ // std::vector<Vert*> _add;
+ // Vert *current = add[0];
// printf("before %d\n",add.size());
// for (unsigned int i=1;i<add.size();i++){
// const double d = sqrt (SQR(add[i]->x()-current->x()) +
@@ -306,7 +313,7 @@ double _fx (const SPoint3 &p, void *)
}
/*
-static void _print (const char *name, std::vector<Vertex*> &T)
+static void _print (const char *name, std::vector<Vert*> &T)
{
FILE *f = fopen(name,"w");
fprintf(f,"View \"\"{\n");
@@ -341,19 +348,18 @@ bool edgeSwaps(tetContainer &T, int myThread)
return false;
}
-
-void edgeBasedRefinement (const int numThreads,
- const int nptsatonce,
- GRegion *gr)
+void edgeBasedRefinement(const int numThreads,
+ const int nptsatonce,
+ GRegion *gr)
{
// fill up old Datastructures
tetContainer allocator (numThreads,1000000);
SBoundingBox3d bb;
- std::vector<Vertex *> _vertices;
+ std::vector<Vert *> _vertices;
edgeContainer ec;
- std::map<Vertex*,MVertex*> _ma;
+ std::map<Vert*,MVertex*> _ma;
{
std::vector<MTetrahedron*> &T = gr->tetrahedra;
@@ -364,7 +370,6 @@ 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){
@@ -373,13 +378,12 @@ void edgeBasedRefinement (const int numThreads,
// }
// fclose(f);
-
_vertices.resize(all.size());
int counter=0;
for (std::set<MVertex*>::iterator it = all.begin();it !=all.end(); ++it){
MVertex *mv = *it;
mv->setIndex(counter);
- Vertex *v = new Vertex (mv->x(),mv->y(),mv->z(),1.e22, counter);
+ Vert *v = new Vert (mv->x(),mv->y(),mv->z(),1.e22, counter);
_vertices[counter] = v;
bb += SPoint3(v->x(),v->y(),v->z());
_ma[v] = mv;
@@ -395,7 +399,8 @@ void edgeBasedRefinement (const int numThreads,
int i1 = tt->getVertex(1)->getIndex();
int i2 = tt->getVertex(2)->getIndex();
int i3 = tt->getVertex(3)->getIndex();
- Tet *t = allocator.newTet(0) ; t->setVertices (_vertices[i0],_vertices[i1],_vertices[i2],_vertices[i3]);
+ Tet *t = allocator.newTet(0) ;
+ t->setVertices(_vertices[i0],_vertices[i1],_vertices[i2],_vertices[i3]);
computeAdjacencies (t,0,faceToTet);
computeAdjacencies (t,1,faceToTet);
computeAdjacencies (t,2,faceToTet);
@@ -414,8 +419,8 @@ void edgeBasedRefinement (const int numThreads,
if (!tt->T[j]){
Face f = tt->getFace(j);
for (int k=0;k<3;k++){
- Vertex *vi = f.V[k];
- Vertex *vj = f.V[(k+1)%3];
+ Vert *vi = f.V[k];
+ Vert *vj = f.V[(k+1)%3];
double l = sqrt ((vi->x()-vj->x())*(vi->x()-vj->x())+
(vi->y()-vj->y())*(vi->y()-vj->y())+
(vi->z()-vj->z())*(vi->z()-vj->z()));
@@ -436,7 +441,7 @@ void edgeBasedRefinement (const int numThreads,
}
}
- std::vector<Vertex*> add_all;
+ std::vector<Vert*> add_all;
{
// vertexFilter _filter (bb, 20);
vertexFilter _filter;
@@ -451,7 +456,7 @@ void edgeBasedRefinement (const int numThreads,
double __t__ = Cpu();
// Tet::in_sphere_counter = 0;
while(1){
- std::vector<Vertex*> add;
+ std::vector<Vert*> add;
double t1 = Cpu();
saturateEdges (ec, allocator, numThreads, add, _fx, NULL);
double t2 = Cpu();
@@ -467,7 +472,8 @@ void edgeBasedRefinement (const int numThreads,
delaunayTrgl (1,1,add.size(), &add,allocator,1.e-28);
double t5 = Cpu();
add_all.insert (add_all.end(), add.begin(), add.end());
- Msg::Info("IT %3d %8d points added, timings %5.2f %5.2f %5.2f %5.2f %5.2f %5d",iter,add.size(),
+ Msg::Info("IT %3d %8d points added, timings %5.2f %5.2f %5.2f %5.2f %5.2f %5d",
+ iter,add.size(),
(t2-t1),
(t3-t2),
(t4-t3),
@@ -485,8 +491,8 @@ void edgeBasedRefinement (const int numThreads,
MVertex *mvs[4];
if (tt->V[0]){
for (int j=0;j<4;j++){
- Vertex *v = tt->V[j];
- std::map<Vertex*,MVertex*>::iterator it = _ma.find(v);
+ Vert *v = tt->V[j];
+ std::map<Vert*,MVertex*>::iterator it = _ma.find(v);
if (it == _ma.end()){
MVertex *mv = new MVertex (v->x(),v->y(),v->z(),gr);
gr->mesh_vertices.push_back(mv);
@@ -525,6 +531,7 @@ void edgeBasedRefinement (const int numThreads,
if (d > 2.)nbBad++;
sum += tau;
}
- Msg::Info("MESH EFFICIENCY : %22.15E %6d edges among %d are out of range",exp (sum / _sizes.size()),nbBad,_sizes.size());
+ Msg::Info("MESH EFFICIENCY : %22.15E %6d edges among %d are out of range",
+ exp (sum / _sizes.size()),nbBad,_sizes.size());
}
}
diff --git a/Mesh/delaunay_refinement.h b/Mesh/delaunay_refinement.h
index 1911d165aa0a1baed46b2460e256660bf484301e..b9cfbcb5c9364dcea6a72d82141232a6a9c67f8b 100644
--- a/Mesh/delaunay_refinement.h
+++ b/Mesh/delaunay_refinement.h
@@ -9,10 +9,10 @@
#include "SPoint3.h"
#include <vector>
class Tet;
-class Vertex;
+class Vert;
void delaunayRefinement (const int numThreads,
const int nptsatonce,
- std::vector<Vertex*> &S,
+ std::vector<Vert*> &S,
std::vector<Tet*> &T,
double (*f)(const SPoint3 &p, void *),
void *data);