From 6be38de8b3a48e451cae3b9c15bdd4747a0cdc48 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Wed, 23 Dec 2015 16:34:07 +0000
Subject: [PATCH] 3d re-writing in progress
---
Mesh/delaunay3d.cpp | 263 +++++++++++++++++++----------
Mesh/delaunay3d_private.h | 9 +-
Mesh/delaunay_refinement.cpp | 92 +++++++---
Mesh/meshGRegion.cpp | 9 +-
Mesh/meshGRegionRelocateVertex.cpp | 23 +--
Mesh/meshGRegionRelocateVertex.h | 2 +-
6 files changed, 268 insertions(+), 130 deletions(-)
diff --git a/Mesh/delaunay3d.cpp b/Mesh/delaunay3d.cpp
index ce3ebac442..9c29cab96d 100644
--- a/Mesh/delaunay3d.cpp
+++ b/Mesh/delaunay3d.cpp
@@ -322,13 +322,12 @@ void computeAdjacencies (Tet *t, int iFace, connContainer &faceToTet){
}
}
-void delaunayCavity (Tet *t,
- Tet *prev,
- Vertex *v,
- cavityContainer &cavity,
- connContainer &bnd,
- int thread, int iPnt){
-
+static void delaunayCavity2 (Tet *t,
+ Tet *prev,
+ Vertex *v,
+ cavityContainer &cavity,
+ connContainer &bnd,
+ int thread, int iPnt, int PP){
t->set(thread, iPnt); // Mark the triangle
cavity.push_back(t);
for (int iNeigh=0; iNeigh<4 ; iNeigh++){
@@ -337,15 +336,40 @@ void delaunayCavity (Tet *t,
bnd.push_back(conn(t->getFace(iNeigh),iNeigh,NULL));
}
else if (neigh == prev){
+ if (PP == 40)printf("coucou\n");
}
else if (!neigh->inSphere(v,thread)){
- // look if v sees face t->getFace(iNeigh)
Face f = t->getFace(iNeigh);
bnd.push_back(conn(f,iNeigh,neigh));
- neigh->set(thread, iPnt);
+ // neigh->set(thread, iPnt);
}
else if (!(neigh->isSet(thread, iPnt))) {
- delaunayCavity (neigh, t, v, cavity,bnd,thread, iPnt);
+ delaunayCavity2 (neigh, t, v, cavity,bnd,thread, iPnt,PP);
+ }
+ }
+}
+
+static void delaunayCavity (Tet *t,
+ Tet *prev,
+ Vertex *v,
+ cavityContainer &cavity,
+ connContainer &bnd,
+ int thread, int iPnt, int PP){
+
+ t->set(thread, iPnt); // Mark the triangle
+ cavity.push_back(t);
+ for (int iNeigh=0; iNeigh<4 ; iNeigh++){
+ Tet *neigh = t->T[iNeigh];
+ if (neigh == NULL){
+ bnd.push_back(conn(t->getFace(iNeigh),iNeigh,NULL));
+ }
+ else if (neigh->inSphere(v,thread)){
+ if (!(neigh->isSet(thread, iPnt))) {
+ delaunayCavity (neigh, t, v, cavity,bnd,thread, iPnt,PP);
+ }
+ }
+ else {
+ bnd.push_back(conn( t->getFace(iNeigh),iNeigh,neigh));
}
}
}
@@ -366,7 +390,7 @@ bool straddling_segment_intersects_triangle(Vertex *p1, Vertex *p2,
}
void print (Vertex *v){
- printf("%3d ",v->getNum());
+ printf("%p ",v);
}
void print (Tet *t, bool recur = true){
@@ -424,64 +448,45 @@ Tet* walk (Tet *t, Vertex *v, int maxx, double &totSearch, int thread){
}
-Tet* walk1 (Tet *t, Vertex *v, int maxx, double &totSearch, int thread){
- int NUMSEARCH = 0;
- int STARTER = 0;
-
- maxx = std::max(100,maxx);
- while (1){
- if (NUMSEARCH++ > 9*maxx) {
- if (t->T[0] && t->T[0]->T[0]) t = t->T[0]->T[0];
- else if (t->T[1] && t->T[1]->T[0]) t = t->T[1]->T[0];
- else if (t->T[2] && t->T[2]->T[0]) t = t->T[2]->T[0];
- if (NUMSEARCH > 10*maxx) {
- printf("trying to find a path from %p to %p: datastructure broken (%d searches, maxx %d)\n",v,t,NUMSEARCH,maxx);
- NUMSEARCH=0;
- return NULL;
- }
- }
- if (t == NULL) {
- Msg::Warning("search went through the boundary of the mesh without finding a tet");
- return NULL; // we should NEVER return here
- }
- if (t->inSphere(v,thread)) {totSearch += NUMSEARCH; return t;}
- Vertex c = t->centroid();
- Tet *oldt=t;
- for (int iNeigh=STARTER; iNeigh<(STARTER+4) ; iNeigh++){
- const int NEIGH = iNeigh %4;
- Face f = t->getFace (NEIGH);
- bool inters = straddling_segment_intersects_triangle(&c, v, f.V[0], f.V[1], f.V[2]);
- if (inters){
- t = t->T[NEIGH];
- break;
- }
- }
- if (t==oldt){
- printf("trying to find a path from %p to %p, ",v,t);
- printf("strange : no intersection\n");
- for (int iNeigh=0; iNeigh<20 ; iNeigh++){
- int rr = rand()%4;
- if (t->T[rr])t =t->T[rr];
- }
- }
- STARTER++;
+void __print (const char *name, connContainer &conn){
+ FILE *f = fopen(name,"w");
+ fprintf(f,"View \"\"{\n");
+
+ 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(),
+ conn[i].f.V[1]->x(),conn[i].f.V[1]->y(),conn[i].f.V[1]->z(),
+ conn[i].f.V[2]->x(),conn[i].f.V[2]->y(),conn[i].f.V[2]->z(),1.,1.,1.);
}
+ fprintf(f,"};\n");
+ fclose(f);
}
-void print (const char *name, std::vector<Tet*> &T){
+void __print (const char *name, std::vector<Tet*> &T, 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<T.size();i++){
if (T[i]->V[0]){
// double val = robustPredicates::orient3d((double*)T[i]->V[0],(double*)T[i]->V[1],(double*)T[i]->V[2],(double*)T[i]->V[3]);
-
- fprintf(f,"SS(%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g,%g};\n",
- T[i]->V[0]->x(),T[i]->V[0]->y(),T[i]->V[0]->z(),
- T[i]->V[1]->x(),T[i]->V[1]->y(),T[i]->V[1]->z(),
- T[i]->V[2]->x(),T[i]->V[2]->y(),T[i]->V[2]->z(),
- T[i]->V[3]->x(),T[i]->V[3]->y(),T[i]->V[3]->z(),
- T[i]->V[0]->lc(),T[i]->V[1]->lc(),T[i]->V[2]->lc(),T[i]->V[3]->lc());
+ if (!v){
+ fprintf(f,"SS(%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g,%g};\n",
+ T[i]->V[0]->x(),T[i]->V[0]->y(),T[i]->V[0]->z(),
+ T[i]->V[1]->x(),T[i]->V[1]->y(),T[i]->V[1]->z(),
+ T[i]->V[2]->x(),T[i]->V[2]->y(),T[i]->V[2]->z(),
+ T[i]->V[3]->x(),T[i]->V[3]->y(),T[i]->V[3]->z(),
+ T[i]->V[0]->lc(),T[i]->V[1]->lc(),T[i]->V[2]->lc(),T[i]->V[3]->lc());
+ }
+ else{
+ fprintf(f,"SS(%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){%d,%d,%d,%d};\n",
+ T[i]->V[0]->x(),T[i]->V[0]->y(),T[i]->V[0]->z(),
+ T[i]->V[1]->x(),T[i]->V[1]->y(),T[i]->V[1]->z(),
+ T[i]->V[2]->x(),T[i]->V[2]->y(),T[i]->V[2]->z(),
+ T[i]->V[3]->x(),T[i]->V[3]->y(),T[i]->V[3]->z(),
+ T[i]->V[0]->getNum(),T[i]->V[1]->getNum(),T[i]->V[2]->getNum(),T[i]->V[3]->getNum());
+ }
}
}
fprintf(f,"};\n");
@@ -550,11 +555,63 @@ static Tet* stoopidSearch (Tet *t, Vertex *v, int thread){
return NULL;
}
+static void stoopidCavity (Tet *t, Vertex *v, int thread, std::vector<Tet*> &cavity){
+ cavity.clear();
+ std::stack<Tet*> _stack;
+ _stack.push(t);
+ std::set<Tet*> all;
+ while(!_stack.empty()){
+ t = _stack.top();
+ all.insert(t);
+ _stack.pop();
+ if (t->inSphere(v,thread)){cavity.push_back(t);}
+ for (int i=0;i<4;i++){
+ if (t->T[i]){
+ if (all.find(t->T[i]) == all.end()){
+ all.insert(t->T[i]);
+ _stack.push(t->T[i]);
+ }
+ }
+ }
+ }
+}
+
+void checkCavity (std::vector<Tet*> &c, Vertex *v){
+ std::set<Face> ff;
+ for (unsigned int i=0; i<c.size(); i++) {
+ for (int j=0;j<4;j++){
+ std::set<Face>::iterator it = ff.find(c[i]->getFace(j));
+ if (it != ff.end())ff.erase(it);
+ else ff.insert(c[i]->getFace(j));
+ }
+ }
+ printf("%d free faces\n",ff.size());
+
+ std::vector<Tet*> bad;
+ for (unsigned int i=0; i<c.size(); i++) {
+ for (int j=0;j<4;j++){
+ std::set<Face>::iterator it = ff.find(c[i]->getFace(j));
+ if (it != ff.end()){
+ Tet *neigh = c[i]->T[j];
+ if (neigh){
+ int val = neigh->inSphere(v,0);
+ if (val){printf("arghhhh %p %p set %d\n",neigh,c[i], neigh->isSet(0,0));
+ bad.push_back(neigh);
+ // bad.push_back(c[i]);
+ }
+ }
+ std::vector<Tet*>::iterator it2= std::find(c.begin(),c.end(),neigh);
+ if (it2 != c.end())printf("datastructure broken\n");
+ }
+ }
+ }
+ __print("bad.pos",bad);
+}
-void delaunayTrgl (const unsigned int numThreads,
- const unsigned int NPTS_AT_ONCE,
- unsigned int Npts,
- std::vector<Tet*> &T,
+void delaunayTrgl (const unsigned int numThreads,
+ const unsigned int NPTS_AT_ONCE,
+ unsigned int Npts,
+ std::vector<Tet*> &T,
std::vector<Vertex*> assignTo[]){
double totSearchGlob=0;
double totCavityGlob=0;
@@ -664,14 +721,35 @@ void delaunayTrgl (const unsigned int numThreads,
}
}
// double c1 = Cpu();
+
+ // FIXME TEST
+ // int NBNULL;
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);
+ delaunayCavity(t[K], NULL, vToAdd[K], cavityK, bndK, myThread, K,iPGlob);
+ 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);
+ //// FIXME TEST
+ // NBNULL = 0;
+ // for (unsigned int i=0; i<cavityK.size(); i++) {
+ // for (int j=0;j<4;j++)if(!cavityK[i]->T[j])NBNULL++;
+ // }
+ /* if (iPGlob==40){
+ char name[245];
+ printf("cavity for point %d size %d shell %d NBNULL %d\n",
+ iPGlob,cavityK.size(),bndK.size(),NBNULL);
+ sprintf(name,"Cavity%d.pos",iPGlob);
+ __print(name,cavityK,vToAdd[K]);
+ sprintf(name,"Bnd%d.pos",iPGlob);
+ __print(name,bndK);
+ checkCavity(cavityK,vToAdd[K]);
+ }
+ */
// verify the cavity
for (unsigned int i=0; i< bndK.size(); i++) {
double val = robustPredicates::orient3d((double*)bndK[i].f.V[0],
@@ -679,10 +757,20 @@ void delaunayTrgl (const unsigned int numThreads,
(double*)bndK[i].f.V[2],
(double*)vToAdd[K]);
if (val <= 0 ) {
- // char name[245];
- // sprintf(name,"invalidCavity%d.pos",invalidCavities [myThread]);
- // print(name,cavityK);
- // printf("val = %22.15E\n",val);
+ // char name[245];
+ // sprintf(name,"invalidCavity%d.pos",invalidCavities [myThread]);
+ // __print(name,cavityK,vToAdd[K]);
+ // printf("val = %22.15E\n",val);
+ // std::vector<Tet*> temp;
+ // stoopidCavity (t[K], vToAdd[K], myThread, temp);
+ // sprintf(name,"validCavity%d.pos",invalidCavities [myThread]);
+ // __print(name,temp,vToAdd[K]);
+ // printf("---> cavity %d vs. %d elements in cavity\n",temp.size(),cavityK.size());
+ // sprintf(name,"Bnd%d.pos",invalidCavities [myThread]);
+ // __print(name,bndK);
+ // for (unsigned int i=0; i<temp.size(); i++) {
+ // print(temp[i]);
+ // }
vToAdd[K] = NULL;
invalidCavities [myThread]++;
break;
@@ -712,6 +800,8 @@ void delaunayTrgl (const unsigned int numThreads,
const unsigned int bSize = bndK.size();
totCavity+= cSize;
Choice[K] = cavityK[0];
+ // FIXME TEST
+ std::vector<Tet*> _news;
for (unsigned int i=0; i<bSize; i++) {
// reuse memory slots of invalid elements
Tet *t;
@@ -723,17 +813,8 @@ void delaunayTrgl (const unsigned int numThreads,
allocatedTable.push_back(t);
newCounter = allocatedTable.size();
}
- // if (newCounter >= allocatedSize){
- // Msg::Fatal("JF : write the reallocation now !");
- // }
if (i < cSize && t->V[0]->_thread != myThread)cashMisses[myThread]++;
- int sign = t->setVertices (bndK[i].f.V[0], bndK[i].f.V[1], bndK[i].f.V[2], vToAdd[K]);
- if (sign <= 0){
- // A coplanar hull face. We need to test if this hull face is
- // Delaunay or not. We test if the adjacent tet (not faked)
- // of this hull face is Delaunay or not.
- Msg::Fatal ("JF: Fix Invalid Cavity %d",sign);
- }
+ t->setVerticesNoTest (bndK[i].f.V[0], bndK[i].f.V[1], bndK[i].f.V[2], vToAdd[K]);
Tet *neigh = bndK[i].t;
t->T[0] = neigh;
t->T[1] = t->T[2] = t->T[3] = NULL;
@@ -747,7 +828,19 @@ void delaunayTrgl (const unsigned int numThreads,
computeAdjacencies (t,1,faceToTet);
computeAdjacencies (t,2,faceToTet);
computeAdjacencies (t,3,faceToTet);
+ _news.push_back(t);
}
+ // FIXME TEST
+ // int NBNULL2 = 0;
+ // for (unsigned int i=0; i<_news.size(); i++) {
+ // for (int j=0;j<4;j++)if(!_news[i]->T[j])NBNULL2++;
+ // }
+ // if (NBNULL != NBNULL2){
+ // char name[245];
+ // sprintf(name,"invalidConnexion%d.pos",iPGlob);
+ // __print(name,_news,vToAdd[K]);
+ // printf("%d %d\n",NBNULL,NBNULL2);
+ // }
for (unsigned int i=bSize; i<cSize; i++) {
counterMiss++;
cavityK[i]->V[0] = cavityK[i]->V[1] = cavityK[i]->V[2] = cavityK[i]->V[3] = NULL;
@@ -774,7 +867,7 @@ void delaunayTrgl (const unsigned int numThreads,
#pragma omp barrier
}
- printf("%d invalid cavities\n",invalidCavities[0]);
+ if (invalidCavities[0])Msg::Warning("%d invalid cavities",invalidCavities[0]);
#ifdef _VERBOSE
double _t2 = walltime(&_t);
@@ -906,7 +999,7 @@ void delaunayTriangulation (const int numThreads,
// _t = 0;
// double _t2 = walltime(&_t);
// printf("WALL CLOCK TIME %12.5E\n",_t2-_t1);
- print ("finalTetrahedrization.pos",T);
+ __print ("finalTetrahedrization.pos",T);
}
@@ -930,9 +1023,9 @@ void delaunayTriangulation (const int numThreads,
for (unsigned int i=0;i<N;i++){
MVertex *mv = S[i];
// FIXME : should be zero !!!!
- double dx = d*1.e-16 * (double)rand() / RAND_MAX;
- double dy = d*1.e-16 * (double)rand() / RAND_MAX;
- double dz = d*1.e-16 * (double)rand() / RAND_MAX;
+ double dx = d*1.e-14 * (double)rand() / RAND_MAX;
+ double dy = d*1.e-14 * (double)rand() / RAND_MAX;
+ double dz = d*1.e-14 * (double)rand() / RAND_MAX;
mv->x() += dx;
mv->y() += dy;
mv->z() += dz;
diff --git a/Mesh/delaunay3d_private.h b/Mesh/delaunay3d_private.h
index 343b701db2..4abd71e548 100644
--- a/Mesh/delaunay3d_private.h
+++ b/Mesh/delaunay3d_private.h
@@ -161,6 +161,11 @@ struct Tet {
T[0] = T[1] = T[2] = T[3] = NULL;
setAllDeleted();
}
+ int setVerticesNoTest (Vertex *v0, Vertex *v1, Vertex *v2, Vertex *v3){
+ _modified=true;
+ V[0] = v0; V[1] = v1; V[2] = v2; V[3] = v3;
+ return 1;
+ }
int setVertices (Vertex *v0, Vertex *v1, Vertex *v2, Vertex *v3){
_modified=true;
double val = robustPredicates::orient3d((double*)v0,(double*)v1,(double*)v2,(double*)v3);
@@ -169,10 +174,12 @@ struct Tet {
return 1;
}
else if (val < 0){
+ // throw;
V[0] = v1; V[1] = v0; V[2] = v2; V[3] = v3;
return -1;
}
else {
+ // throw;
return 0;
}
}
@@ -232,7 +239,7 @@ typedef std::vector<conn> connContainer;
void SortHilbert (std::vector<Vertex*>& v, std::vector<int> &indices);
void computeAdjacencies (Tet *t, int iFace, connContainer &faceToTet);
-void print (const char *name, std::vector<Tet*> &T);
+void __print (const char *name, std::vector<Tet*> &T, Vertex *v = 0);
void delaunayTrgl (const unsigned int numThreads,
const unsigned int NPTS_AT_ONCE,
unsigned int Npts,
diff --git a/Mesh/delaunay_refinement.cpp b/Mesh/delaunay_refinement.cpp
index 857db4ab39..76f0d4de74 100644
--- a/Mesh/delaunay_refinement.cpp
+++ b/Mesh/delaunay_refinement.cpp
@@ -35,14 +35,6 @@ struct edgeContainer
_hash2.insert(e);
return true;
}
- bool isNewEdge (const Edge &e) {
- size_t h = (size_t) e.first >> 3;
- std::vector<Edge> &v = _hash[h %_hash.size()];
- AVGSEARCH+=v.size();
- for (unsigned int i=0; i< v.size();i++)if (e == v[i]) {return false;}
- return true;
- }
-
bool addNewEdge (const Edge &e)
{
size_t h = (size_t) e.first >> 3;
@@ -148,9 +140,9 @@ void saturateEdge (Edge &e, std::vector<Vertex*> &S, double (*f)(const SPoint3 &
else {
SPoint3 p = p1 * (1.-t) + p2*t;
double lc = e.first->lc() * (1.-t) + e.second->lc()*t;
- const double dx = 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 dx = 1.e-10 * (double) rand() / RAND_MAX;
+ const double dy = 1.e-10 * (double) rand() / RAND_MAX;
+ const double dz = 1.e-10 * (double) rand() / RAND_MAX;
S.push_back(new Vertex(p.x()+dx,p.y()+dy,p.z()+dz,lc));
L += interval;
}
@@ -198,16 +190,16 @@ public :
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->z() - _v->z()) * (p->_v->z() - _v->z()));
- return d < .7*p->_v->lc();//_l;
+ return d < .8*p->_v->lc();//_l;
}
void minmax (double _min[3], double _max[3]) const
{
- _min[0] = _v->x() - 1.1*_v->lc();
- _min[1] = _v->y() - 1.1*_v->lc();
- _min[2] = _v->z() - 1.1*_v->lc();
- _max[0] = _v->x() + 1.1*_v->lc();
- _max[1] = _v->y() + 1.1*_v->lc();
- _max[2] = _v->z() + 1.1*_v->lc();
+ _min[0] = _v->x() - 1.01*_v->lc();
+ _min[1] = _v->y() - 1.01*_v->lc();
+ _min[2] = _v->z() - 1.01*_v->lc();
+ _max[0] = _v->x() + 1.01*_v->lc();
+ _max[1] = _v->y() + 1.01*_v->lc();
+ _max[2] = _v->z() + 1.01*_v->lc();
}
};
@@ -317,13 +309,14 @@ void edgeBasedRefinement (const int numThreads,
// fill up old Datastructures
- std::vector<MTetrahedron*> T = gr->tetrahedra;
std::vector<Tet*> _tets;
- _tets.resize(T.size());
std::vector<Vertex *> _vertices;
edgeContainer ec;
+ std::map<Vertex*,MVertex*> _ma;
{
+ std::vector<MTetrahedron*> &T = gr->tetrahedra;
+ _tets.resize(T.size());
std::set<MVertex *> all;
for (int i=0;i<T.size();i++){
for (int j=0;j<4;j++){
@@ -338,6 +331,7 @@ void edgeBasedRefinement (const int numThreads,
mv->setIndex(counter);
Vertex *v = new Vertex (mv->x(),mv->y(),mv->z(),1.e22, counter);
_vertices[counter] = v;
+ _ma[v] = mv;
counter++;
}
connSet faceToTet;
@@ -352,9 +346,11 @@ void edgeBasedRefinement (const int numThreads,
computeAdjacencies (t,2,faceToTet);
computeAdjacencies (t,3,faceToTet);
_tets[i] = t;
+ delete T[i];
}
+ T.clear();
}
-
+
// do not allow to saturate boundary edges
{
for (unsigned int i=0;i<_tets.size();i++) {
@@ -377,12 +373,13 @@ void edgeBasedRefinement (const int numThreads,
for (unsigned int i=0;i<_tets.size();i++) {
for (int j=0;j<6;j++){
Edge e = _tets[i]->getEdge(j);
- if(e.first->lc() == 1.e22){printf("coucou\n");e.first->lc() = e.second->lc();}
- else if(e.second->lc() == 1.e22){printf("coucou\n");e.second->lc() = e.first->lc();}
+ if(e.first->lc() == 1.e22){/*printf("coucou\n");*/e.first->lc() = e.second->lc();}
+ else if(e.second->lc() == 1.e22){/*printf("coucou\n");*/e.second->lc() = e.first->lc();}
}
}
}
+ std::vector<Vertex*> add_all;
{
vertexFilter _filter;
int iter = 1;
@@ -410,8 +407,9 @@ void edgeBasedRefinement (const int numThreads,
double t4 = Cpu();
// sprintf(name,"PointsFiltered%d.pos",iter);
// _print (name,add);
- delaunayTrgl (1,1,add.size(), _tets, &add);
- double t5 = Cpu();
+ delaunayTrgl (1,1,add.size(), _tets, &add);
+ add_all.insert (add_all.end(), add.begin(), add.end());
+ clock_t t5 = clock();
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),
@@ -421,6 +419,48 @@ void edgeBasedRefinement (const int numThreads,
_tets.size());
iter++;
}
- print ("afterRefinement.pos",_tets);
+ // print ("afterRefinement.pos",_tets);
+ }
+
+ std::map<Edge,double> _sizes;
+ for (unsigned int i=0; i<_tets.size();i++){
+ MVertex *mvs[4];
+
+
+ if (_tets[i]->V[0]){
+ for (int j=0;j<6;j++){
+ Edge e = _tets[i]->getEdge(j);
+ std::map<Edge,double>::iterator it = _sizes.find(e);
+ if (it == _sizes.end()){
+ double l = sqrt ((e.first->x() - e.second->x()) * (e.first->x() - e.second->x()) +
+ (e.first->y() - e.second->y()) * (e.first->y() - e.second->y()) +
+ (e.first->z() - e.second->z()) * (e.first->z() - e.second->z()));
+ _sizes[e]= 2* l / (e.first->lc() + e.second->lc());
+ }
+ }
+ for (int j=0;j<4;j++){
+ Vertex *v = _tets[i]->V[j];
+ std::map<Vertex*,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);
+ _ma[v] = mv;
+ mvs[j] = mv;
+ }
+ else mvs[j] = it->second;
+ }
+ gr->tetrahedra.push_back(new MTetrahedron(mvs[0],mvs[1],mvs[2],mvs[3]));
+ }
+ delete _tets[i];
}
+
+ std::map<Edge,double>::iterator it = _sizes.begin();
+ double sum = 0;
+ for (; it !=_sizes.end();++it){
+ double d = it->second;
+ double tau = d < 1 ? d - 1 : 1./d - 1;
+ sum += tau;
+ }
+ Msg::Info("MESH EFFICIENCY : %22.15E",exp (sum / _sizes.size()));
+
}
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 20e13256c1..b826362385 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -628,17 +628,14 @@ static void MeshDelaunayVolumeNewCode(std::vector<GRegion*> ®ions) {
// restore the initial set of faces
gr->set(faces);
-
-
-#if 1
+
void edgeBasedRefinement (const int numThreads,
const int nptsatonce,
GRegion *gr);
+ // just to remove tets that are not to be meshed
insertVerticesInRegion(gr,0);
edgeBasedRefinement (1,1,gr);
-#else
- insertVerticesInRegion(gr,200000000);
-#endif
+ // RelocateVertices (regions,-1);
}
diff --git a/Mesh/meshGRegionRelocateVertex.cpp b/Mesh/meshGRegionRelocateVertex.cpp
index 1c004f84d7..10085b9596 100644
--- a/Mesh/meshGRegionRelocateVertex.cpp
+++ b/Mesh/meshGRegionRelocateVertex.cpp
@@ -34,15 +34,14 @@ static double objective_function (double xi, MVertex *ver,
#define sqrt5 2.236067977499789696
-static int Stopping_Rule(double x0, double x1)
+static int Stopping_Rule(double x0, double x1, double tol)
{
- double tolerance = 1.e-2;
- return ( fabs( x1 - x0 ) < tolerance ) ? 1 : 0;
+ return ( fabs( x1 - x0 ) < tol) ? 1 : 0;
}
double Maximize_Quality_Golden_Section( MVertex *ver,
double xTarget, double yTarget, double zTarget,
- const std::vector<MElement*> < )
+ const std::vector<MElement*> < , double tol)
{
static const double lambda = 0.5 * (sqrt5 - 1.0);
@@ -55,18 +54,20 @@ double Maximize_Quality_Golden_Section( MVertex *ver,
double fx1 = objective_function (x1, ver, xTarget, yTarget, zTarget , lt );
double fx2 = objective_function (x2, ver, xTarget, yTarget, zTarget , lt );
- while ( ! Stopping_Rule( a, b ) ) {
+ if (tol < 0.0)return fx1 > fx2 ? x1 : x2;
+
+ while ( ! Stopping_Rule( a, b , tol) ) {
// printf("GOLDEN : %g %g (%12.5E,%12.5E)\n",a,b,fa,fb);
if (fx1 < fx2) {
a = x1;
- if ( Stopping_Rule( a, b ) ) break;
+ if ( Stopping_Rule( a, b , tol) ) break;
x1 = x2;
fx1 = fx2;
x2 = b - mu * (b - a);
fx2 = objective_function (x2, ver, xTarget, yTarget, zTarget , lt );
} else {
b = x2;
- if ( Stopping_Rule( a, b ) ) break;
+ if ( Stopping_Rule( a, b , tol) ) break;
x2 = x1;
fx2 = fx1;
x1 = a + mu * (b - a);
@@ -78,7 +79,7 @@ double Maximize_Quality_Golden_Section( MVertex *ver,
}
static void _relocateVertex(MVertex *ver,
- const std::vector<MElement*> <)
+ const std::vector<MElement*> <, double tol)
{
if(ver->onWhat()->dim() != 3) return;
double x = 0, y=0, z=0;
@@ -96,13 +97,13 @@ static void _relocateVertex(MVertex *ver,
N += lt[i]->getNumVertices();
}
- double xi = Maximize_Quality_Golden_Section( ver, x/N, y/N, z/N, lt );
+ double xi = Maximize_Quality_Golden_Section( ver, x/N, y/N, z/N, lt , tol);
ver->x() = (1.-xi) * ver->x() + xi * x/N;
ver->y() = (1.-xi) * ver->y() + xi * y/N;
ver->z() = (1.-xi) * ver->z() + xi * z/N;
}
-void RelocateVertices (std::vector<GRegion*> ®ions) {
+void RelocateVertices (std::vector<GRegion*> ®ions, double tol) {
for(unsigned int k = 0; k < regions.size(); k++){
v2t_cont adj;
buildVertexToElement(regions[k]->tetrahedra, adj);
@@ -111,7 +112,7 @@ void RelocateVertices (std::vector<GRegion*> ®ions) {
buildVertexToElement(regions[k]->hexahedra, adj);
v2t_cont::iterator it = adj.begin();
while (it != adj.end()){
- _relocateVertex( it->first, it->second);
+ _relocateVertex( it->first, it->second, tol);
++it;
}
}
diff --git a/Mesh/meshGRegionRelocateVertex.h b/Mesh/meshGRegionRelocateVertex.h
index 678403aade..5a1b49e429 100644
--- a/Mesh/meshGRegionRelocateVertex.h
+++ b/Mesh/meshGRegionRelocateVertex.h
@@ -2,5 +2,5 @@
#define _MESHGREGIONRELOCATEVERTEX_
#include <vector>
class GRegion;
-void RelocateVertices (std::vector<GRegion*> ®ions);
+void RelocateVertices (std::vector<GRegion*> ®ions, double tol = 1.e-2);
#endif
--
GitLab