From feae43c6fec35f85a96b76c4538b86e6210397dd Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Tue, 11 Feb 2014 10:41:21 +0000
Subject: [PATCH] 3D Delaunnay-ization of a set of points is now fully gmsh
based
---
Mesh/meshGRegionDelaunayInsertion.cpp | 159 +++++++++++++++-
Mesh/meshGRegionDelaunayInsertion.h | 2 +
Numeric/CMakeLists.txt | 2 +
Numeric/HilbertCurve.cpp | 264 ++++++++++++++++++++++++++
Numeric/HilbertCurve.h | 4 +
5 files changed, 422 insertions(+), 9 deletions(-)
create mode 100644 Numeric/HilbertCurve.cpp
create mode 100644 Numeric/HilbertCurve.h
diff --git a/Mesh/meshGRegionDelaunayInsertion.cpp b/Mesh/meshGRegionDelaunayInsertion.cpp
index c34c4a8e0f..481854e1ba 100644
--- a/Mesh/meshGRegionDelaunayInsertion.cpp
+++ b/Mesh/meshGRegionDelaunayInsertion.cpp
@@ -18,6 +18,7 @@
#include "MTriangle.h"
#include "Numeric.h"
#include "Context.h"
+#include "HilbertCurve.h"
int MTet4::radiusNorm = 2;
static double LIMIT_ = 1;
@@ -330,15 +331,6 @@ void nonrecurFindCavity(std::list<faceXtet> & shell,
MVertex *v ,
MTet4 *t)
{
- // Msg::Info("tet %d %d %d %d",t->tet()->getVertex(0)->getNum(),
- // t->tet()->getVertex(1)->getNum(),
- // t->tet()->getVertex(2)->getNum(),
- // t->tet()->getVertex(3)->getNum());
-
- // invariant : this one has to be inserted in the cavity
- // consider this tet deleted
- // remove its reference to its neighbors
-
std::stack<MTet4*> _stack;
_stack.push(t);
while(!_stack.empty()){
@@ -1612,3 +1604,152 @@ void bowyerWatsonFrontalLayers(GRegion *gr, bool hex)
MTet4::radiusNorm = 2;
LIMIT_ = 1;
}
+
+
+///// do a 3D delaunay mesh assuming a set of vertices
+
+static void initialCube (std::vector<MVertex*> &v,
+ MVertex *box[8],
+ std::vector<MTet4*> &t){
+ SBoundingBox3d bbox ;
+ for (size_t i=0;i<v.size();i++){
+ MVertex *pv = v[i];
+ bbox += SPoint3(pv->x(),pv->y(),pv->z());
+ }
+ bbox *= 1.3;
+ box[0] = new MVertex (bbox.min().x(),bbox.min().y(),bbox.min().z());
+ box[1] = new MVertex (bbox.max().x(),bbox.min().y(),bbox.min().z());
+ box[2] = new MVertex (bbox.max().x(),bbox.max().y(),bbox.min().z());
+ box[3] = new MVertex (bbox.min().x(),bbox.max().y(),bbox.min().z());
+ box[4] = new MVertex (bbox.min().x(),bbox.min().y(),bbox.max().z());
+ box[5] = new MVertex (bbox.max().x(),bbox.min().y(),bbox.max().z());
+ box[6] = new MVertex (bbox.max().x(),bbox.max().y(),bbox.max().z());
+ box[7] = new MVertex (bbox.min().x(),bbox.max().y(),bbox.max().z());
+ std::vector<MTetrahedron*> t_box;
+ MTetrahedron *t0 = new MTetrahedron (box[2],box[7],box[3],box[1]);
+ MTetrahedron *t1 = new MTetrahedron (box[0],box[7],box[1],box[3]);
+ MTetrahedron *t2 = new MTetrahedron (box[6],box[1],box[7],box[2]);
+ MTetrahedron *t3 = new MTetrahedron (box[0],box[1],box[7],box[4]);
+ MTetrahedron *t4 = new MTetrahedron (box[1],box[4],box[5],box[7]);
+ MTetrahedron *t5 = new MTetrahedron (box[1],box[7],box[5],box[6]);
+ t.push_back(new MTet4(t0,0.0));
+ t.push_back(new MTet4(t1,0.0));
+ t.push_back(new MTet4(t2,0.0));
+ t.push_back(new MTet4(t3,0.0));
+ t.push_back(new MTet4(t4,0.0));
+ t.push_back(new MTet4(t5,0.0));
+ connectTets(t);
+}
+
+int intersect_line_triangle(double X[3], double Y[3], double Z[3] ,
+ double P[3], double N[3], double);
+
+static MTet4* search4Tet (MTet4 *t, MVertex *v, int _size) {
+ if (t->inCircumSphere(v)) return t;
+ int ITER = 0;
+ SPoint3 p2 (v->x(),v->y(),v->z());
+ while (1){
+ SPoint3 p1 = t->tet()->barycenter();
+ int found = -1;
+ for (int i = 0; i < 4; i++){
+ MTet4 *neigh = t->getNeigh(i);
+ if (neigh){
+ faceXtet fxt (t, i);
+ double X[3] = {fxt.v[0]->x(),fxt.v[1]->x(),fxt.v[2]->x()};
+ double Y[3] = {fxt.v[0]->y(),fxt.v[1]->y(),fxt.v[2]->y()};
+ double Z[3] = {fxt.v[0]->z(),fxt.v[1]->z(),fxt.v[2]->z()};
+ // printf("%d %d\n",i,intersect_line_triangle(X,Y,Z,p1,p2-p1));
+ if (intersect_line_triangle(X,Y,Z,p1,p1-p2, 1.e-12) > 0) {
+ found = i;
+ break;
+ }
+ }
+ }
+ if (found < 0){
+ return 0;
+ }
+ t = t->getNeigh(found);
+ if (t->inCircumSphere(v)) {
+ // printf("FOUND\n");
+ return t;
+ }
+ if (ITER++ > _size) break;
+ }
+ return 0;
+}
+
+
+MTet4 * getTetToBreak (MVertex *v, std::vector<MTet4*> &t, int &NB_GLOBAL_SEARCH){
+ // last inserted is used as starting point
+ // we know it is not deleted
+ MTet4 *start = t[t.size() - 1];
+ start = search4Tet (start,v,(int)t.size());
+ if (start)return start;
+ NB_GLOBAL_SEARCH++;
+ for (size_t i = 0;i<t.size();i++){
+ if (!t[i]->isDeleted() && t[i]->inCircumSphere(v))return t[i];
+ }
+ return 0;
+}
+
+bool tetOnBox (MTetrahedron *t, MVertex *box[8]){
+ for (size_t i = 0;i<4;i++)
+ for (size_t j = 0;j<8;j++)
+ if (t->getVertex(i) == box[j])return true;
+ return false;
+}
+
+void delaunayMeshIn3D(std::vector<MVertex*> &v, std::vector<MTetrahedron*> &result)
+{
+ std::vector<MTet4*> t;
+ MVertex *box[8];
+ initialCube (v,box,t);
+
+ int NB_GLOBAL_SEARCH = 0;
+ SortHilbert(v);
+ clock_t t1 = clock();
+
+ for (size_t i=0;i<v.size();i++){
+ MVertex *pv = v[i];
+ MTet4 * found = getTetToBreak (pv,t,NB_GLOBAL_SEARCH);
+ std::list<faceXtet> shell;
+ std::list<MTet4*> cavity;
+ recurFindCavity(shell, cavity, pv, found);
+ std::vector<MTet4*> extended_cavity;
+ std::list<faceXtet>::iterator it = shell.begin();
+ while (it != shell.end()){
+ MTetrahedron *tr = new MTetrahedron(it->getVertex(0),
+ it->getVertex(1),
+ it->getVertex(2), pv);
+ MTet4 *t4 = new MTet4(tr, 0.0);
+ t.push_back(t4);
+ extended_cavity.push_back(t4);
+ MTet4 *otherSide = it->t1->getNeigh(it->i1);
+ if (otherSide)
+ extended_cavity.push_back(otherSide);
+ ++it;
+ }
+ // printf("connecting %i tets\n",extended_cavity.size());
+ connectTets_vector(extended_cavity.begin(),extended_cavity.end());
+ }
+
+ clock_t t2 = clock();
+ printf("%d global searches among %d CPU = %g\n",NB_GLOBAL_SEARCH,v.size(),(double)(t2-t1)/CLOCKS_PER_SEC);
+
+
+ // FILE *f = fopen ("tet.pos","w");
+ // fprintf(f,"View \"\"{\n");
+ for (size_t i = 0;i<t.size();i++){
+ if (t[i]->isDeleted() || tetOnBox (t[i]->tet(),box)) delete t[i]->tet();
+ else {
+ result.push_back(t[i]->tet());
+ // t[i]->tet()->writePOS (f, false,false,true,false,false,false);
+ }
+ delete t[i];
+ }
+ // fprintf(f,"};\n");
+ // fclose(f);
+}
+
+
+
diff --git a/Mesh/meshGRegionDelaunayInsertion.h b/Mesh/meshGRegionDelaunayInsertion.h
index 15a467cd8a..9360b51439 100644
--- a/Mesh/meshGRegionDelaunayInsertion.h
+++ b/Mesh/meshGRegionDelaunayInsertion.h
@@ -173,6 +173,8 @@ class MTet4
void connectTets(std::list<MTet4*> &);
void connectTets(std::vector<MTet4*> &);
+// IN --> Vertices ---- OUT --> Tets
+void delaunayMeshIn3D(std::vector<MVertex*> &, std::vector<MTetrahedron*> &);
void insertVerticesInRegion(GRegion *gr, int maxVert = 2000000000, bool _classify = true);
void bowyerWatsonFrontalLayers(GRegion *gr, bool hex);
GRegion *getRegionFromBoundingFaces(GModel *model,
diff --git a/Numeric/CMakeLists.txt b/Numeric/CMakeLists.txt
index a3382e8d92..ae0f52d721 100644
--- a/Numeric/CMakeLists.txt
+++ b/Numeric/CMakeLists.txt
@@ -7,6 +7,7 @@ set(SRC
Numeric.cpp
fullMatrix.cpp
BasisFactory.cpp
+ discreteFrechetDistance.cpp
nodalBasis.cpp
polynomialBasis.cpp
pyramidalBasis.cpp
@@ -27,6 +28,7 @@ set(SRC
GaussQuadraturePyr.cpp
GaussLegendreSimplex.cpp
GaussJacobi1D.cpp
+ HilbertCurve.cpp
robustPredicates.cpp
mathEvaluator.cpp
Iso.cpp
diff --git a/Numeric/HilbertCurve.cpp b/Numeric/HilbertCurve.cpp
new file mode 100644
index 0000000000..a8f11067d6
--- /dev/null
+++ b/Numeric/HilbertCurve.cpp
@@ -0,0 +1,264 @@
+#include "SBoundingBox3d.h"
+#include "MVertex.h"
+
+
+struct HilbertSort
+{
+// The code for generating table transgc
+// from: http://graphics.stanford.edu/~seander/bithacks.html.
+ int transgc[8][3][8];
+ int tsb1mod3[8];
+ int maxDepth;
+ int Limit;
+ SBoundingBox3d bbox ;
+ void ComputeGrayCode(int n);
+ int Split(MVertex** vertices,
+ int arraysize,int GrayCode0,int GrayCode1,
+ double BoundingBoxXmin, double BoundingBoxXmax,
+ double BoundingBoxYmin, double BoundingBoxYmax,
+ double BoundingBoxZmin, double BoundingBoxZmax);
+ void Sort(MVertex** vertices, int arraysize, int e, int d,
+ double BoundingBoxXmin, double BoundingBoxXmax, double BoundingBoxYmin, double BoundingBoxYmax,
+ double BoundingBoxZmin, double BoundingBoxZmax, int depth);
+ HilbertSort (int m = 0, int l=1) : maxDepth(m),Limit(l)
+ {
+ ComputeGrayCode(3);
+ }
+ void MultiscaleSortHilbert(MVertex** vertices, int arraysize,
+ int threshold, double ratio, int *depth)
+ {
+ int middle;
+
+ middle = 0;
+ if (arraysize >= threshold) {
+ (*depth)++;
+ middle = arraysize * ratio;
+ MultiscaleSortHilbert(vertices, middle, threshold, ratio, depth);
+ }
+ Sort (&(vertices[middle]),arraysize - middle,0,0,
+ bbox.min().x(),bbox.max().x(),
+ bbox.min().y(),bbox.max().y(),
+ bbox.min().z(),bbox.max().z(),0);
+ }
+
+ void Apply (std::vector<MVertex*> &v)
+ {
+ for (size_t i=0;i<v.size();i++){
+ MVertex *pv = v[i];
+ bbox += SPoint3(pv->x(),pv->y(),pv->z());
+ }
+ bbox *= 1.01;
+ MVertex**pv = &v[0];
+ int depth;
+ MultiscaleSortHilbert(pv, v.size(), 128, 0.125,&depth);
+ }
+};
+
+
+void HilbertSort::ComputeGrayCode(int n)
+{
+ int gc[8], N, mask, travel_bit;
+ int e, d, f, k, g;
+ int v, c;
+ int i;
+
+ N = (n == 2) ? 4 : 8;
+ mask = (n == 2) ? 3 : 7;
+
+ // Generate the Gray code sequence.
+ for (i = 0; i < N; i++) {
+ gc[i] = i ^ (i >> 1);
+ }
+
+ for (e = 0; e < N; e++) {
+ for (d = 0; d < n; d++) {
+ // Calculate the end point (f).
+ f = e ^ (1 << d); // Toggle the d-th bit of 'e'.
+ // travel_bit = 2**p, the bit we want to travel.
+ travel_bit = e ^ f;
+ for (i = 0; i < N; i++) {
+ // // Rotate gc[i] left by (p + 1) % n bits.
+ k = gc[i] * (travel_bit * 2);
+ g = ((k | (k / N)) & mask);
+ // Calculate the permuted Gray code by xor with the start point (e).
+ transgc[e][d][i] = (g ^ e);
+ }
+ // assert(transgc[e][d][0] == e);
+ // assert(transgc[e][d][N - 1] == f);
+ } // d
+ } // e
+
+ // Count the consecutive '1' bits (trailing) on the right.
+ tsb1mod3[0] = 0;
+ for (i = 1; i < N; i++) {
+ v = ~i; // Count the 0s.
+ v = (v ^ (v - 1)) >> 1; // Set v's trailing 0s to 1s and zero rest
+ for (c = 0; v; c++) {
+ v >>= 1;
+ }
+ tsb1mod3[i] = c % n;
+ }
+}
+
+
+int HilbertSort::Split(MVertex** vertices,
+ int arraysize,int GrayCode0,int GrayCode1,
+ double BoundingBoxXmin, double BoundingBoxXmax, double BoundingBoxYmin, double BoundingBoxYmax,
+ double BoundingBoxZmin, double BoundingBoxZmax)
+{
+ MVertex* swapvert;
+ int axis, d;
+ double split;
+ int i, j;
+
+
+ // Find the current splitting axis. 'axis' is a value 0, or 1, or 2, which
+ // correspoding to x-, or y- or z-axis.
+ axis = (GrayCode0 ^ GrayCode1) >> 1;
+
+ // Calulate the split position along the axis.
+ if (axis == 0) {
+ split = 0.5 * (BoundingBoxXmin + BoundingBoxXmax);
+ } else if (axis == 1) {
+ split = 0.5 * (BoundingBoxYmin + BoundingBoxYmax);
+ } else { // == 2
+ split = 0.5 * (BoundingBoxZmin + BoundingBoxZmax);
+ }
+
+ // Find the direction (+1 or -1) of the axis. If 'd' is +1, the direction
+ // of the axis is to the positive of the axis, otherwise, it is -1.
+ d = ((GrayCode0 & (1<<axis)) == 0) ? 1 : -1;
+
+
+ // Partition the vertices into left- and right-arrays such that left points
+ // have Hilbert indices lower than the right points.
+ i = 0;
+ j = arraysize - 1;
+
+ // Partition the vertices into left- and right-arrays.
+ if (d > 0) {
+ do {
+ for (; i < arraysize; i++) {
+ if (vertices[i]->point()[axis] >= split) break;
+ }
+ for (; j >= 0; j--) {
+ if (vertices[j]->point()[axis] < split) break;
+ }
+ // Is the partition finished?
+ if (i == (j + 1)) break;
+ // Swap i-th and j-th vertices.
+ swapvert = vertices[i];
+ vertices[i] = vertices[j];
+ vertices[j] = swapvert;
+ // Continue patitioning the array;
+ } while (true);
+ } else {
+ do {
+ for (; i < arraysize; i++) {
+ if (vertices[i]->point()[axis] <= split) break;
+ }
+ for (; j >= 0; j--) {
+ if (vertices[j]->point()[axis] > split) break;
+ }
+ // Is the partition finished?
+ if (i == (j + 1)) break;
+ // Swap i-th and j-th vertices.
+ swapvert = vertices[i];
+ vertices[i] = vertices[j];
+ vertices[j] = swapvert;
+ // Continue patitioning the array;
+ } while (true);
+ }
+
+ return i;
+}
+
+// The sorting code is inspired by Tetgen 1.5
+
+void HilbertSort::Sort(MVertex** vertices, int arraysize, int e, int d,
+ double BoundingBoxXmin, double BoundingBoxXmax, double BoundingBoxYmin, double BoundingBoxYmax,
+ double BoundingBoxZmin, double BoundingBoxZmax, int depth)
+{
+ double x1, x2, y1, y2, z1, z2;
+ int p[9], w, e_w, d_w, k, ei, di;
+ int n = 3, mask = 7;
+
+ p[0] = 0;
+ p[8] = arraysize;
+
+ p[4] = Split(vertices, p[8], transgc[e][d][3], transgc[e][d][4],
+ BoundingBoxXmin, BoundingBoxXmax, BoundingBoxYmin, BoundingBoxYmax, BoundingBoxZmin, BoundingBoxZmax);
+ p[2] = Split(vertices, p[4], transgc[e][d][1], transgc[e][d][2],
+ BoundingBoxXmin, BoundingBoxXmax, BoundingBoxYmin, BoundingBoxYmax, BoundingBoxZmin, BoundingBoxZmax);
+ p[1] = Split(vertices, p[2], transgc[e][d][0], transgc[e][d][1],
+ BoundingBoxXmin, BoundingBoxXmax, BoundingBoxYmin, BoundingBoxYmax, BoundingBoxZmin, BoundingBoxZmax);
+ p[3] = Split(&(vertices[p[2]]), p[4] - p[2],
+ transgc[e][d][2], transgc[e][d][3],
+ BoundingBoxXmin, BoundingBoxXmax, BoundingBoxYmin, BoundingBoxYmax, BoundingBoxZmin, BoundingBoxZmax) + p[2];
+ p[6] = Split(&(vertices[p[4]]), p[8] - p[4],
+ transgc[e][d][5], transgc[e][d][6],
+ BoundingBoxXmin, BoundingBoxXmax, BoundingBoxYmin, BoundingBoxYmax, BoundingBoxZmin, BoundingBoxZmax) + p[4];
+ p[5] = Split(&(vertices[p[4]]), p[6] - p[4],
+ transgc[e][d][4], transgc[e][d][5],
+ BoundingBoxXmin, BoundingBoxXmax, BoundingBoxYmin, BoundingBoxYmax, BoundingBoxZmin, BoundingBoxZmax) + p[4];
+ p[7] = Split(&(vertices[p[6]]), p[8] - p[6],
+ transgc[e][d][6], transgc[e][d][7],
+ BoundingBoxXmin, BoundingBoxXmax, BoundingBoxYmin, BoundingBoxYmax, BoundingBoxZmin, BoundingBoxZmax) + p[6];
+
+ if (maxDepth > 0) {
+ if ((depth + 1) == maxDepth) {
+ return;
+ }
+ }
+
+ // Recursively sort the points in sub-boxes.
+ for (w = 0; w < 8; w++) {
+ if ((p[w+1] - p[w]) > Limit) {
+ if (w == 0) {
+ e_w = 0;
+ } else {
+ k = 2 * ((w - 1) / 2);
+ e_w = k ^ (k >> 1);
+ }
+ k = e_w;
+ e_w = ((k << (d+1)) & mask) | ((k >> (n-d-1)) & mask);
+ ei = e ^ e_w;
+ if (w == 0) {
+ d_w = 0;
+ } else {
+ d_w = ((w % 2) == 0) ? tsb1mod3[w - 1] : tsb1mod3[w];
+ }
+ di = (d + d_w + 1) % n;
+ if (transgc[e][d][w] & 1) {
+ x1 = 0.5 * (BoundingBoxXmin + BoundingBoxXmax);
+ x2 = BoundingBoxXmax;
+ } else {
+ x1 = BoundingBoxXmin;
+ x2 = 0.5 * (BoundingBoxXmin + BoundingBoxXmax);
+ }
+ if (transgc[e][d][w] & 2) { // y-axis
+ y1 = 0.5 * (BoundingBoxYmin + BoundingBoxYmax);
+ y2 = BoundingBoxYmax;
+ } else {
+ y1 = BoundingBoxYmin;
+ y2 = 0.5 * (BoundingBoxYmin + BoundingBoxYmax);
+ }
+ if (transgc[e][d][w] & 4) { // z-axis
+ z1 = 0.5 * (BoundingBoxZmin + BoundingBoxZmax);
+ z2 = BoundingBoxZmax;
+ } else {
+ z1 = BoundingBoxZmin;
+ z2 = 0.5 * (BoundingBoxZmin + BoundingBoxZmax);
+ }
+ Sort(&(vertices[p[w]]), p[w+1] - p[w], ei, di,
+ x1, x2, y1, y2, z1, z2, depth+1);
+ }
+ }
+}
+
+
+
+void SortHilbert (std::vector<MVertex*>& v){
+ HilbertSort h;
+ h.Apply(v);
+}
diff --git a/Numeric/HilbertCurve.h b/Numeric/HilbertCurve.h
new file mode 100644
index 0000000000..b413f587ab
--- /dev/null
+++ b/Numeric/HilbertCurve.h
@@ -0,0 +1,4 @@
+#ifndef _HILBERT_CURVE_
+#define _HILBERT_CURVE_
+void SortHilbert (std::vector<MVertex*>&);
+#endif
--
GitLab