From 6a4cb61a6a55566f8276b3f3a64d54239d9575ca Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Fri, 29 Jan 2016 09:33:19 +0000
Subject: [PATCH] new 3D mesher online
---
Mesh/delaunay3d.cpp | 251 +++++++++++++++++++++++++--
Mesh/delaunay3d_private.h | 2 +-
Mesh/delaunay_refinement.cpp | 25 ++-
Mesh/meshGRegionBoundaryRecovery.cpp | 28 +--
Mesh/meshGRegionBoundaryRecovery.h | 4 +-
benchmarks/3d/Cube-01.geo | 2 +-
6 files changed, 280 insertions(+), 32 deletions(-)
diff --git a/Mesh/delaunay3d.cpp b/Mesh/delaunay3d.cpp
index 7ba803bdd4..f075922b24 100644
--- a/Mesh/delaunay3d.cpp
+++ b/Mesh/delaunay3d.cpp
@@ -319,6 +319,231 @@ void computeAdjacencies (Tet *t, int iFace, connContainer &faceToTet){
faceToTet.erase(it);
}
}
+/*
+ Fixing a non star shaped cavity (non delaunay triangulations)
+ See P.L. George's paper
+ Improvements on Delaunay-based three-dimensional automatic mesh generator
+ Finite Elements in Analysis and Design 25 (1997) 297-317
+*/
+
+static void starShapeness (Vertex *v, connContainer &bndK,
+ std::vector<unsigned int> &_negatives,
+ std::vector<unsigned int> &_nulls,
+ double threshold){
+ _negatives.clear();
+ _nulls.clear();
+ for (unsigned int i=0; i< bndK.size(); i++) {
+ // no symbolic perturbation
+ const double val = robustPredicates::orient3d((double*)bndK[i].f.V[0],
+ (double*)bndK[i].f.V[1],
+ (double*)bndK[i].f.V[2],
+ (double*)v);
+ if (val <= threshold ) {
+ if(val >= 0)_nulls.push_back(i);
+ else _negatives.push_back(i);
+ }
+ }
+}
+
+static Tet* tetContainsV (Vertex *v, cavityContainer &cavity){
+
+ for (unsigned int i=0; i< cavity.size(); i++) {
+ unsigned int count = 0;
+ for (unsigned int j=0;j<4;j++){
+ Face f = cavity[i]->getFace(j);
+ const double val = robustPredicates::orient3d((double*)f.V[0],
+ (double*)f.V[1],
+ (double*)f.V[2],
+ (double*)v);
+ if (val >= 0 ) {
+ count++;
+ }
+ }
+ if (count == 4)return cavity[i];
+ }
+ // printf("AIE\n");
+ return NULL;
+}
+
+static void buildDelaunayBall (cavityContainer &cavity, connContainer &faceToTet){
+ faceToTet.clear();
+ for (unsigned int i=0; i< cavity.size(); i++) {
+ Tet *t = cavity[i];
+ for (unsigned int iFace=0; iFace< 4; iFace++) {
+ Tet *neigh = t->T[iFace];
+ conn c (t->getFace(iFace), iFace, neigh);
+ connContainer::iterator it = std::find(faceToTet.begin(),faceToTet.end(),c);
+ if (it == faceToTet.end()){
+ faceToTet.push_back(c);
+ }
+ else {
+ faceToTet.erase(it);
+ }
+ }
+ }
+}
+
+static void findPossibleTets (cavityContainer &cavity,
+ connContainer &bndK,
+ std::vector<Tet*> &_possible){
+ _possible.clear();
+ std::vector<Vertex*> _all;
+ for (unsigned int i=0; i< bndK.size(); i++) {
+ for (unsigned int j=0; j< 3; j++) {
+ _all.push_back(bndK[i].f.V[j]);
+ }
+ }
+ // std::vector<Vertex*>::iterator last = std::unique(_all.begin(), _all.end());
+ for (unsigned int i=0; i< cavity.size(); i++) {
+ Tet *t = cavity[i];
+ for (unsigned int j=0; j< 4; j++) {
+ Vertex *v = t->V[j];
+ if (std::find (_all.begin(),_all.end(), v) == _all.end()){
+ _possible.push_back(t);
+ break;
+ }
+ }
+ }
+}
+
+static bool removeIsolatedTets(Tet *containsV,
+ cavityContainer &cavity,
+ connContainer &bndK,
+ int myThread, int K){
+ cavityContainer cc;
+ cc.push_back(containsV);
+ std::stack<Tet*> _stack;
+ _stack.push(containsV);
+
+ while (!_stack.empty()){
+ Tet *t = _stack.top();
+ _stack.pop();
+ for (int i=0;i<4;i++){
+ Tet *neigh = t->T[i];
+ if (neigh &&
+ (std::find(cc.begin(),cc.end(),neigh) == cc.end()) &&
+ (std::find(cavity.begin(),cavity.end(),neigh) != cavity.end())){
+ cc.push_back(neigh);
+ _stack.push(neigh);
+ }
+ }
+ }
+ if (cc.size() == cavity.size())return false;
+ // // Msg::Info(" cavity updated (%3ld elements) %3ld isolated tet removed",cavity.size(),cavity.size()-cc.size());
+ cavity = cc;
+ return true;
+}
+
+static Tet *tetInsideCavityWithFAce (Face &f, cavityContainer &cavity){
+ // printf("size of cavity %ld\n",cavity.size());
+ for (unsigned int i=0; i< cavity.size(); i++) {
+ Tet *t = cavity[i];
+ for (unsigned int iFace=0; iFace< 4; iFace++) {
+ if (t->getFace(iFace) == f) return t;
+ }
+ }
+ return NULL;
+}
+
+
+void __print (const char *name, connContainer &conn, Vertex *v);
+static bool fixDelaunayCavity (double threshold,
+ Vertex *v,
+ cavityContainer &cavity,
+ connContainer &bndK,
+ int myThread, int K,
+ std::vector<unsigned int> & _negatives,
+ std::vector<unsigned int> & _nulls){
+
+ static int XXX = 1;
+ starShapeness (v, bndK, _negatives, _nulls, threshold);
+
+#if 0
+ // test coherence ...
+ connContainer bndK2;
+ buildDelaunayBall (cavity,bndK2);
+ if (bndK.size() != bndK2.size())printf("BUG1\n");
+ for (unsigned int i=0;i<bndK.size();i++){
+ connContainer::iterator it = std::find(bndK.begin(),bndK.end(),bndK2[i]);
+ if( it == bndK.end())printf("%d BUG2\n",XXX);
+ if(it->t != bndK2[i].t)printf("%d BUG3\n",XXX);
+ if(it->i != bndK2[i].i)printf("%d BUG4\n",XXX);
+ }
+#endif
+
+
+ if (_negatives.empty() && _nulls.empty())return false;
+
+ // unset all tets of the cavity
+ for (unsigned int i=0; i< cavity.size(); i++)cavity[i]->unset(myThread,K);
+ for (unsigned int i=0; i<bndK.size(); i++)if(bndK[i].t)bndK[i].t->unset(myThread,K);
+
+ Msg::Debug("Fixing cavity %3d (%3ld,%3ld) : %ld negatives and %ld Nulls",XXX,
+ cavity.size(),bndK.size(), _negatives.size(),_nulls.size());
+
+ Tet *containsV = tetContainsV (v,cavity);
+
+ if (! containsV) return true;
+
+ while (1) {
+ if (removeIsolatedTets(containsV, cavity,bndK,myThread,K)){
+ buildDelaunayBall (cavity,bndK);
+ starShapeness (v, bndK, _negatives, _nulls, threshold);
+ }
+
+ for (unsigned int i=0;i<_negatives.size();i++){
+ // FIXME
+ int index = _negatives[i];
+ conn &c = bndK[index];
+ Tet *toRemove = tetInsideCavityWithFAce (c.f, cavity);
+ if (toRemove){
+ std::vector<Tet*>::iterator it = std::find(cavity.begin(), cavity.end(), toRemove);
+ if (it != cavity.end())
+ cavity.erase(it);
+ else
+ Msg::Fatal("Datastructure Broken in %s line %5d",__FILE__,__LINE__);
+ }
+ }
+ for (unsigned int i=0;i<_nulls.size();i++){
+ return true;
+ int index = _nulls[i];
+ conn &c = bndK[index];
+ Tet *toAdd = c.t;
+ if (toAdd && std::find(cavity.begin(), cavity.end(), toAdd) == cavity.end())cavity.push_back(toAdd);
+ else return true;
+ }
+ // rebuild the boundary of the cavity
+ // this step could be done faster by directly
+ // updating bndK in the steps above
+ buildDelaunayBall (cavity,bndK);
+ // if isolated points are present in the cavity
+ // find all tets that contain thos vertices
+ std::vector<Tet*> _possible;
+ if (_nulls.size())findPossibleTets (cavity,bndK,_possible);
+ if (!_possible.empty()){
+ printf("cavity is dead\n");
+ throw;
+ // FIXME
+ return true;
+ std::random_shuffle(_possible.begin(), _possible.end());
+ Tet *toRemove = _possible[0];
+ std::remove(cavity.begin(), cavity.end(), toRemove);
+ buildDelaunayBall (cavity,bndK);
+ starShapeness (v, bndK, _negatives, _nulls, threshold);
+ }
+ else {
+ starShapeness (v, bndK, _negatives, _nulls, threshold);
+ // Msg::Info(" cavity updated (%3ld elements) %3ld negatives",cavity.size(),_negatives.size());
+ if (_negatives.empty() && _nulls.empty()){
+ // sprintf(name,"final%d.pos",XXX++);
+ // __print (name, bndK, v);
+ //Msg::Info(" cavity is fixed");
+ // getchar();
+ return false;
+ }
+ }
+ }
+}
static void delaunayCavity2 (Tet *t,
Tet *prev,
@@ -375,20 +600,23 @@ Tet* walk (Tet *t, Vertex *v, int maxx, double &totSearch, int thread){
if (t == NULL) {
return NULL; // we should NEVER return here
}
- if (t->inSphere(v,thread)) {return t;}
+ // if (t->inSphere(v,thread)) {return t;}
double _min = 0.0;
int NEIGH = -1;
+ int count = 0;
for (int iNeigh=0; iNeigh<4; iNeigh++){
Face f = t->getFace (iNeigh);
double val = robustPredicates::orient3d((double*)f.V[0],
(double*)f.V[1],
(double*)f.V[2],
(double*)v);
+ if (val >=-1.e-12) count++;
if (val < _min){
NEIGH = iNeigh;
_min = val;
}
}
+ if (count == 4 && t->inSphere(v,thread))return t;
if (NEIGH >= 0){
t = t->T[NEIGH];
}
@@ -400,10 +628,12 @@ Tet* walk (Tet *t, Vertex *v, int maxx, double &totSearch, int thread){
}
-void __print (const char *name, connContainer &conn){
+void __print (const char *name, connContainer &conn, Vertex *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());
+
for (unsigned int i=0;i<conn.size();i++){
fprintf(f,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
conn[i].f.V[0]->x(),conn[i].f.V[0]->y(),conn[i].f.V[0]->z(),
@@ -522,7 +752,8 @@ void delaunayTrgl (const unsigned int numThreads,
const unsigned int NPTS_AT_ONCE,
unsigned int Npts,
std::vector<Vertex*> assignTo[],
- tetContainer &allocator){
+ tetContainer &allocator,
+ double threshold){
#ifdef _VERBOSE
double totSearchGlob=0;
double totCavityGlob=0;
@@ -552,6 +783,7 @@ void delaunayTrgl (const unsigned int numThreads,
double totSearch=0;
double totCavity=0;
+ std::vector<unsigned int> _negatives, _nulls;
std::vector<cavityContainer> cavity(NPTS_AT_ONCE);
std::vector<connContainer> bnd(NPTS_AT_ONCE);
std::vector<bool> ok(NPTS_AT_ONCE);
@@ -620,16 +852,9 @@ void delaunayTrgl (const unsigned int numThreads,
delaunayCavity2(t[K], NULL, vToAdd[K], cavityK, bndK, myThread, K);
//delaunayCavity(t[K], vToAdd[K], cavityK, bndK, myThread, K);
// verify the cavity
- for (unsigned int i=0; i< bndK.size(); i++) {
- 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]++;
- break;
- }
+ if (fixDelaunayCavity (threshold, vToAdd[K], cavityK, bndK, myThread, K, _negatives, _nulls)){
+ vToAdd[K] = NULL;
+ invalidCavities [myThread]++;
}
}
}
diff --git a/Mesh/delaunay3d_private.h b/Mesh/delaunay3d_private.h
index 36af19a9bf..a71990cb00 100644
--- a/Mesh/delaunay3d_private.h
+++ b/Mesh/delaunay3d_private.h
@@ -338,6 +338,6 @@ void delaunayTrgl (const unsigned int numThreads,
const unsigned int NPTS_AT_ONCE,
unsigned int Npts,
std::vector<Vertex*> assignTo[],
- tetContainer &allocator);
+ tetContainer &allocator, double threshold = 0.0);
#endif
diff --git a/Mesh/delaunay_refinement.cpp b/Mesh/delaunay_refinement.cpp
index 5bd391a679..803ddc8b39 100644
--- a/Mesh/delaunay_refinement.cpp
+++ b/Mesh/delaunay_refinement.cpp
@@ -263,8 +263,25 @@ void filterVertices (const int numThreads,
std::vector<int> indices;
SortHilbert(add, indices);
+ std::vector<Vertex*> _add=add;
+
+ // std::vector<Vertex*> _add;
+ // Vertex *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()) +
+ // SQR(add[i]->y()-current->y()) +
+ // SQR(add[i]->z()-current->z()) );
+ // if (0.8*current->lc() > d){
+ // delete add[i];
+ // }
+ // else {
+ // current = add[i];
+ // _add.push_back(add[i]);
+ // }
+ // }
+ // printf("after %d\n",_add.size());
- std::vector<Vertex*> _add = add;
add.clear();
for (unsigned int i=0;i<_add.size();i++){
SPoint3 p (_add[i]->x(),_add[i]->y(),_add[i]->z());
@@ -415,7 +432,7 @@ void edgeBasedRefinement (const int numThreads,
int iter = 1;
- Msg::Info("----------------------------------- SATUR FILTR SORTH DELNY TIME TETS");
+ Msg::Info("------------------------------------- SATUR FILTR SORTH DELNY TIME TETS");
double __t__ = Cpu();
// Tet::in_sphere_counter = 0;
@@ -433,10 +450,10 @@ void edgeBasedRefinement (const int numThreads,
std::vector<int> indices;
SortHilbert(add, indices);
double t4 = Cpu();
- delaunayTrgl (1,1,add.size(), &add,allocator);
+ delaunayTrgl (1,1,add.size(), &add,allocator,1.e-24);
double t5 = Cpu();
add_all.insert (add_all.end(), add.begin(), add.end());
- 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 %8d points added, timings %5.2f %5.2f %5.2f %5.2f %5.2f %5d",iter,add.size(),
(t2-t1),
(t3-t2),
(t4-t3),
diff --git a/Mesh/meshGRegionBoundaryRecovery.cpp b/Mesh/meshGRegionBoundaryRecovery.cpp
index 58db969066..91f82a9f23 100644
--- a/Mesh/meshGRegionBoundaryRecovery.cpp
+++ b/Mesh/meshGRegionBoundaryRecovery.cpp
@@ -12894,7 +12894,7 @@ long meshGRegionBoundaryRecovery::lawsonflip3d(flipconstraints *fc)
return totalcount + sliver_peels;
}
-void meshGRegionBoundaryRecovery::recoverdelaunay()
+bool meshGRegionBoundaryRecovery::recoverdelaunay()
{
arraypool *flipqueue, *nextflipqueue, *swapqueue;
triface tetloop, neightet, *parytet;
@@ -12941,13 +12941,13 @@ void meshGRegionBoundaryRecovery::recoverdelaunay()
fc.enqflag = 2;
lawsonflip3d(&fc);
-
+
if (b->verbose > 1) {
printf(" obj (after Lawson) = %.17g\n", tetprism_vol_sum);
}
-
+
if (unflipqueue->objects == 0l) {
- return; // The mesh is Delaunay.
+ return true; // The mesh is Delaunay.
}
fc.unflip = 1; // Unflip if the edge is not flipped.
@@ -13042,7 +13042,9 @@ void meshGRegionBoundaryRecovery::recoverdelaunay()
}
} // while (flipqueue->objects > 0l)
+ bool returnValue = true;
if (flipqueue->objects > 0l) {
+ returnValue = false;
if (b->verbose > 1) {
printf(" %ld non-Delaunay edges remained.\n", flipqueue->objects);
}
@@ -13055,6 +13057,7 @@ void meshGRegionBoundaryRecovery::recoverdelaunay()
b->flipstarsize = bakmaxflipstarsize;
delete flipqueue;
delete nextflipqueue;
+ return returnValue;
}
int meshGRegionBoundaryRecovery::gettetrahedron(point pa, point pb, point pc,
@@ -14450,11 +14453,11 @@ void meshGRegionBoundaryRecovery::jettisonnodes()
///////////////////////////////////////////////////////////////////////////////
-void meshGRegionBoundaryRecovery::reconstructmesh(GRegion *_gr)
+bool meshGRegionBoundaryRecovery::reconstructmesh(GRegion *_gr)
{
+ bool returnValue (true);
double t_start = Cpu();
-
std::vector<MVertex*> _vertices;
// Get the set of vertices from GRegion.
@@ -14512,7 +14515,7 @@ void meshGRegionBoundaryRecovery::reconstructmesh(GRegion *_gr)
longest = sqrt(x * x + y * y + z * z);
if (longest == 0.0) {
Msg::Warning("Error: The point set is trivial.\n");
- return;
+ return true;
}
// Two identical points are distinguished by 'lengthlimit'.
@@ -14854,9 +14857,10 @@ void meshGRegionBoundaryRecovery::reconstructmesh(GRegion *_gr)
suppresssteinerpoints();
}
- recoverdelaunay();
+ returnValue = recoverdelaunay();
- // optimizemesh();
+ // let's trry
+ optimizemesh();
if ((dupverts > 0l) || (unuverts > 0l)) {
// Remove hanging nodes.
@@ -15159,8 +15163,10 @@ void meshGRegionBoundaryRecovery::reconstructmesh(GRegion *_gr)
_gr->tetrahedra.push_back(t);
tetloop.tet = tetrahedrontraverse();
}
-} // mesh output
- Msg::Info("Reconstruct time : %g sec",Cpu()-t_start);
+ } // mesh output
+ if (returnValue)Msg::Info("Reconstruct time : %g sec (mesh is Delaunay)",Cpu()-t_start);
+ else Msg::Info("Reconstruct time : %g sec (mesh is not Delaunay)",Cpu()-t_start);
+ return returnValue;
}
void terminateBoundaryRecovery(void *, int exitcode)
diff --git a/Mesh/meshGRegionBoundaryRecovery.h b/Mesh/meshGRegionBoundaryRecovery.h
index 0023617e71..0bfc189d40 100644
--- a/Mesh/meshGRegionBoundaryRecovery.h
+++ b/Mesh/meshGRegionBoundaryRecovery.h
@@ -743,7 +743,7 @@ class meshGRegionBoundaryRecovery {
// Mesh optimize
long lawsonflip3d(flipconstraints *fc);
- void recoverdelaunay();
+ bool recoverdelaunay();
int gettetrahedron(point, point, point, point, triface *);
long improvequalitybyflips();
int smoothpoint(point smtpt, arraypool*, int ccw, optparameters *opm);
@@ -874,7 +874,7 @@ class meshGRegionBoundaryRecovery {
void unifysubfaces(face *f1, face *f2);
void unifysegments();
void jettisonnodes();
- void reconstructmesh(GRegion *_gr);
+ bool reconstructmesh(GRegion *_gr);
};
void terminateBoundaryRecovery(void *, int exitcode);
diff --git a/benchmarks/3d/Cube-01.geo b/benchmarks/3d/Cube-01.geo
index 239b5af5e4..da9b61afa8 100644
--- a/benchmarks/3d/Cube-01.geo
+++ b/benchmarks/3d/Cube-01.geo
@@ -18,4 +18,4 @@ Line Loop(5) = {2,3,4,1};
Plane Surface(6) = {5};
Extrude Surface { 6, {0,0.0,1} };
//Characteristic Length {10} = lc/100;
-Physical Point(1)={1 };
+//Physical Point(1)={1 };
--
GitLab