diff --git a/Mesh/delaunay3d.cpp b/Mesh/delaunay3d.cpp
index 99c7813a3a7114a7489abf33bdc4f27731cfa808..e6d4b668a420f6ce5dab751d1c9d639d68285c97 100644
--- a/Mesh/delaunay3d.cpp
+++ b/Mesh/delaunay3d.cpp
@@ -23,6 +23,7 @@
#include "MVertex.h"
#include "MEdge.h"
#include "MTetrahedron.h"
+#include "qualityMeasures.h"
#include "meshGRegionLocalMeshMod.h"
#if defined(_HAVE_NUMA)
@@ -31,6 +32,22 @@
//int Tet::in_sphere_counter = 0;
+#define SQR(X) (X)*(X)
+
+inline double tetQuality (Vert *vx0,Vert *vx1,Vert *vx2,Vert *vx3, double *volume)
+{
+ return qmTetrahedron::gamma (vx0->x(),vx0->y(),vx0->z(),
+ vx1->x(),vx1->y(),vx1->z(),
+ vx2->x(),vx2->y(),vx2->z(),
+ vx3->x(),vx3->y(),vx3->z(),volume);
+}
+
+double tetQuality (Tet* t, double *volume)
+{
+ return tetQuality(t->V[0],t->V[1],t->V[2],t->V[3], volume);
+}
+
+
struct HilbertSortB
{
// The code for generating table transgc from:
@@ -332,9 +349,12 @@ void computeAdjacencies (Tet *t, int iFace, connContainer &faceToTet)
************************************************************************/
static int edges[6][2] = {{0,1},{0,2},{0,3},{1,2},{1,3},{2,3}};
-static int efaces[6][2] = {{0,2},{0,1},{1,2},{0,3},{2,3},{1,3}};
-static int faces[4][3] = {{0,1,2},{0,2,3},{0,1,3},{1,2,3}};
-//static int vnofaces[4] = {3,1,2,0};
+//static int efaces[6][2] = {{0,2},{0,1},{1,2},{0,3},{2,3},{1,3}};
+//static int faces[4][3] = {{0,1,2},{0,2,3},{0,1,3},{1,2,3}};
+static int efaces[6][2] = {{0,3},{0,2},{2,3},{0,1},{1,3},{1,2}};
+static int faces[4][3] = {{0,1,2},{1,3,2},{2,3,0},{1,0,3}};
+static int fnovert[4] = {1,2,3,0};
+
//static int vFac[4][3] = {{0,1,2},{0,2,3},{0,1,3},{1,2,3}};
static void computeNeighboringTetsOfACavity(const std::vector<Tet*> &cavity,
@@ -352,6 +372,66 @@ static void computeNeighboringTetsOfACavity(const std::vector<Tet*> &cavity,
}
}
+void printtet (const char *c, Tet *t){
+ printf("%s ",c);
+
+ if (t->V[0])
+ printf("%7d %7d %7d %7d \n",
+ t->V[0]->getNum(),
+ t->V[1]->getNum(),
+ t->V[2]->getNum(),
+ t->V[3]->getNum());
+
+}
+
+static bool buildVertexCavity(Vert *v, std::vector<Tet*> &cavity, std::stack<Tet*> &_stack)
+{
+ // _stack.clear();
+
+ Tet *t = v->getT();
+ if (!t){
+ // Msg::Error("Vertex with no tet connected");
+ return false;
+ }
+ // std::stack<Tet*> _stack;
+ _stack.push (t);
+ while (!_stack.empty()){
+ t = _stack.top();
+ _stack.pop();
+ if (std::find(cavity.begin(), cavity.end(),t) == cavity.end()){
+ cavity.push_back(t);
+ int i;
+ for (i=0;i<4;i++){
+ if (t->V[i] == v)break;
+ }
+ if (i == 4){
+ Msg::Error("Error in Vertex Cavity");
+ while (!_stack.empty())_stack.pop();
+ return false;
+ }
+
+ if (cavity.size() > 1000){
+ Msg::Error("Error in Vertex Cavity (too big)");
+ while (!_stack.empty())_stack.pop();
+ return false;
+ }
+
+
+ for (int j=0;j<4;j++){
+ if (j!=fnovert[i]){
+ if (t->T[j])_stack.push(t->T[j]);
+ else {
+ while (!_stack.empty())_stack.pop();
+ return false;
+ }
+ }
+ }
+ }
+ }
+ return true;
+}
+
+
static bool buildEdgeCavity(Tet *t, int iLocalEdge,
Vert **v1, Vert **v2,
std::vector<Tet*> &cavity,
@@ -363,15 +443,20 @@ static bool buildEdgeCavity(Tet *t, int iLocalEdge,
*v1 = t->V[edges[iLocalEdge][0]];
*v2 = t->V[edges[iLocalEdge][1]];
-
+
// the 5 - i th edge contains the other 2 points of the tet
Vert *lastinring = t->V[edges[5 - iLocalEdge][0]];
+
ring.push_back(lastinring);
cavity.push_back(t);
+ // printf("edge %d %d \n",(*v1)->getNum(),(*v2)->getNum());
+ /// printtet ("first", t);
+
while (1){
Vert *ov1 = t->V[edges[5 - iLocalEdge][0]];
Vert *ov2 = t->V[edges[5 - iLocalEdge][1]];
+ // printf("%p %p %p %p\n",t->V[0],t->V[1],t->V[2],t->V[3]);
int K = ov1 == lastinring ? 1 : 0;
lastinring = ov1 == lastinring ? ov2 : ov1;
// look in the 2 faces sharing this edge the one that has vertex
@@ -388,10 +473,14 @@ static bool buildEdgeCavity(Tet *t, int iLocalEdge,
else { Msg::Error("Error of connexion"); return false; }
t=t->T[iFace];
if (!t) return false;
+ // printf("xh\n");
+ // printtet ("next", t);
+ // printf("ams\n");
if (!t->V[0]){ Msg::Error("Weird!!"); return false; }
if (t == cavity[0]) break;
ring.push_back(lastinring);
cavity.push_back(t);
+ if (cavity.size() > 8)return false;
iLocalEdge = -1;
for (int i = 0; i < 6; i++){
Vert *a = t->V[edges[i][0]];
@@ -402,17 +491,25 @@ static bool buildEdgeCavity(Tet *t, int iLocalEdge,
}
}
if (iLocalEdge == -1){
- Msg::Error("loc = %d", iLocalEdge);
+ for (int i = 0; i < 6; i++){
+ Vert *a = t->V[edges[i][0]];
+ Vert *b = t->V[edges[i][1]];
+ printf ("%d %d vs. %d %d\n",a->getNum(),b->getNum(),(*v1)->getNum(),(*v2)->getNum());
+ }
+ // getchar();
+ Msg::Error("loc = %d size %lu", iLocalEdge, cavity.size());
return false;
}
}
+
+ // printf("%d elements of the cavity\n",cavity.size());
computeNeighboringTetsOfACavity (cavity,outside);
+ // for (int i=0;i<outside.size();++i)printtet("outside",outside[i]);
+ // getchar();
return true;
}
-#define SQR(X) (X)*(X)
-
-bool edgeSwap(Tet *tet, int iLocalEdge, tetContainer &T, int myThread)
+static bool edgeSwap(Tet *tet, int iLocalEdge, tetContainer &T, int myThread)
{
std::vector<Tet*> cavity;
std::vector<Tet*> outside;
@@ -422,6 +519,7 @@ bool edgeSwap(Tet *tet, int iLocalEdge, tetContainer &T, int myThread)
bool closed = buildEdgeCavity(tet, iLocalEdge, &v1, &v2, cavity, outside, ring);
if (!closed) return false;
+ // return false;
SwapPattern sp;
switch (ring.size()){
@@ -434,55 +532,38 @@ bool edgeSwap(Tet *tet, int iLocalEdge, tetContainer &T, int myThread)
}
double tetQualityRef = 1.e22;
+ double volumeOld = 0, volume;
for (unsigned int i = 0; i < cavity.size(); i++){
- 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()) +
- SQR (vx0->x() - vx2->x()) + SQR (vx0->y() - vx2->y()) + SQR (vx0->z() - vx2->z()) +
- SQR (vx0->x() - vx3->x()) + SQR (vx0->y() - vx3->y()) + SQR (vx0->z() - vx3->z()) +
- SQR (vx1->x() - vx2->x()) + SQR (vx1->y() - vx2->y()) + SQR (vx1->z() - vx2->z()) +
- SQR (vx1->x() - vx3->x()) + SQR (vx1->y() - vx3->y()) + SQR (vx1->z() - vx3->z()) +
- SQR (vx2->x() - vx3->x()) + SQR (vx2->y() - vx3->y()) + SQR (vx2->z() - vx3->z()) ;
- tetQualityRef = std::min(tetQualityRef, volume / a);
+ tetQualityRef = std::min(tetQualityRef, tetQuality (cavity[i], &volume));
+ volumeOld += volume;
}
// compute qualities of all tets that appear in the patterns
double tetQuality1[100], tetQuality2[100];
+ double volume1[100], volume2[100];
for (int i = 0; i < sp.nbr_triangles; i++){
// FIXME VERIFY ORIENTATION OF TRIANGULAR PATTERNS
int p1 = sp.triangles[i][0];
int p2 = sp.triangles[i][1];
int p3 = sp.triangles[i][2];
- const double volume1 = robustPredicates::orient3d((double*)ring[p1], (double*)ring[p2], (double*)ring[p3], (double*) v1);
- const double volume2 = robustPredicates::orient3d((double*)ring[p1], (double*)ring[p2], (double*)ring[p3], (double*) v2);
- const double a12 = SQR (ring[p1]->x() - ring[p2]->x())+SQR (ring[p1]->y() - ring[p2]->y())+SQR (ring[p1]->z() - ring[p2]->z());
- const double a13 = SQR (ring[p1]->x() - ring[p3]->x())+SQR (ring[p1]->y() - ring[p3]->y())+SQR (ring[p1]->z() - ring[p3]->z());
- const double a23 = SQR (ring[p2]->x() - ring[p3]->x())+SQR (ring[p2]->y() - ring[p3]->y())+SQR (ring[p2]->z() - ring[p3]->z());
- const double a123 = a12 + a13 + a23;
- const double v11 = SQR (ring[p1]->x() - v1->x())+SQR (ring[p1]->y() - v1->y())+SQR (ring[p1]->z() - v1->z());
- const double v12 = SQR (ring[p2]->x() - v1->x())+SQR (ring[p2]->y() - v1->y())+SQR (ring[p2]->z() - v1->z());
- const double v13 = SQR (ring[p3]->x() - v1->x())+SQR (ring[p3]->y() - v1->y())+SQR (ring[p3]->z() - v1->z());
- const double v21 = SQR (ring[p1]->x() - v2->x())+SQR (ring[p1]->y() - v2->y())+SQR (ring[p1]->z() - v2->z());
- const double v22 = SQR (ring[p2]->x() - v2->x())+SQR (ring[p2]->y() - v2->y())+SQR (ring[p2]->z() - v2->z());
- const double v23 = SQR (ring[p3]->x() - v2->x())+SQR (ring[p3]->y() - v2->y())+SQR (ring[p3]->z() - v2->z());
- tetQuality1[i] = volume1 / (a123 + v11 + v12 + v13);
- tetQuality2[i] = volume2 / (a123 + v21 + v22 + v23);
+ tetQuality1[i] = tetQuality (ring[p1],ring[p2],ring[p3],v1,&volume1[i]);
+ tetQuality2[i] = tetQuality (ring[p3],ring[p2],ring[p1],v2,&volume2[i]);
}
-
+
// look for the best triangulation, i.e. the one that maximize the
// minimum element quality
double minQuality[100];
+ double volumes[100];
// for all triangulations
for (int i = 0; i < sp.nbr_trianguls; i++){
+ volumes[i] = 0;
// for all triangles in a triangulation
minQuality[i] = 1.e22;
for (int j = 0; j < sp.nbr_triangles_2; j++){
int iT = sp.trianguls[i][j];
+ // if (ring.size() == 5)printf("triangulation %d triangle %d qualities %g %g\n",i,j,tetQuality1[iT],tetQuality2[iT]);
minQuality[i] = std::min(minQuality[i], std::min (tetQuality1[iT],tetQuality2[iT]));
+ volumes[i] += volume1[iT] + volume2[iT];
}
}
@@ -495,7 +576,12 @@ bool edgeSwap(Tet *tet, int iLocalEdge, tetContainer &T, int myThread)
}
}
+ if (fabs(volumes[iBest] - volumeOld) > 1.e-14*volumeOld){
+ // printf("NOPE --> RING %d BEST %12.5E INITIAL %12.5E VOLS %12.5E %12.5E\n",ring.size(),best, tetQualityRef,volumes[iBest],volumeOld);
+ return false;
+ }
if (best <= tetQualityRef) return false;
+ // printf("RING %d BEST %12.5E INITIAL %12.5E VOLS %12.5E %12.5E\n",ring.size(),best, tetQualityRef,volumes[iBest],volumeOld);
unsigned int counter = 0;
for (int j = 0; j < sp.nbr_triangles_2; j++){
@@ -513,8 +599,10 @@ bool edgeSwap(Tet *tet, int iLocalEdge, tetContainer &T, int myThread)
if (counter < cavity.size()) nt2 = cavity[counter++];
else nt2 = T.newTet(myThread);
- nt1->setVertices(pv1, pv2, pv3, v1);
- nt2->setVertices(pv3, pv2, pv3, v1);
+ nt1->setVertices(pv1, pv2, pv3, v2);
+ nt2->setVertices(pv1, pv3, pv2, v1);
+ nt1->T[0] = nt1->T[1] = nt1->T[2] = nt1->T[3] = NULL;
+ nt2->T[0] = nt2->T[1] = nt2->T[2] = nt2->T[3] = NULL;
outside.push_back(nt1);
outside.push_back(nt2);
@@ -523,18 +611,116 @@ bool edgeSwap(Tet *tet, int iLocalEdge, tetContainer &T, int myThread)
for(unsigned int i = counter; i < cavity.size(); i++) cavity[i]->V[0] = NULL;
connContainer ctnr;
- for(unsigned int i = 0; i < outside.size(); i++)
- for (int j=0;j<4;j++)
+ for(unsigned int i = 0; i < outside.size(); i++){
+ for (int j=0;j<4;j++){
computeAdjacencies (outside[i],j,ctnr);
-
+ }
+ }
+
return true;
}
-/*
-static void edgeSwapPass (tetContainer &T)
+static bool relocateVertex (Vert* v, std::vector<Tet*> &cavity){
+
+ double x=0, y=0, z=0;
+ double minQual = 1.0;
+ double volumeOld = 0, volume;
+ for (unsigned int i=0;i<cavity.size();i++){
+ minQual = std::min(minQual, tetQuality(cavity[i], &volume));
+ volumeOld += volume;
+ for (unsigned int j=0;j<4;j++){
+ x += cavity[i]->V[j]->x();
+ y += cavity[i]->V[j]->y();
+ z += cavity[i]->V[j]->z();
+ }
+ }
+ if (minQual > 0.3)return false;
+ x /= 4.*cavity.size();
+ y /= 4.*cavity.size();
+ z /= 4.*cavity.size();
+ double oldX=v->x(), oldY=v->y(), oldZ=v->z();
+ v->x()=x, v->y()=y, v->z()=z;
+ // double minQualNew = 1.0;
+ double volumeNew = 0;
+ for (unsigned int i=0;i<cavity.size();i++){
+ if ( tetQuality(cavity[i], &volume) < minQual){
+ v->x()=oldX, v->y()=oldY, v->z()=oldZ;
+ return false;
+ }
+ volumeNew += volume;
+ }
+
+ if (fabs(volumeNew-volumeOld) < 1.e-14*volumeOld) {
+ return true;
+ }
+ v->x()=oldX, v->y()=oldY, v->z()=oldZ;
+ return false;
+
+}
+
+static bool relocateVertex (Vert* v, std::stack<Tet*> &_work, std::vector<Tet*> &cavity){
+ cavity.clear();
+ if (buildVertexCavity(v, cavity, _work)){
+ return relocateVertex (v, cavity);
+ }
+ return false;
+}
+
+void vertexRelocationPass (int numThreads, std::vector<Vert*> &v){
+
+ int N = 0;
+ int iter = 0;
+ std::stack<Tet*> _work;
+ std::vector<Tet*> _work2;
+ while (1){
+ for (unsigned int i=0;i<v.size();i++){
+ if (relocateVertex (v[i],_work, _work2))N++;
+ }
+ // printf("done\n");
+ if (!N)break;
+ if (iter++ >= 0)break;
+ }
+
+}
+
+
+void edgeSwapPass (int numThreads, tetContainer &allocator, edgeContainer &embeddedEdges)
{
+ edgeContainer ec;
+ while(1){
+ int N = 0;
+ for (int myThread=0; myThread < numThreads; myThread++) {
+ int NNN = allocator.size(myThread);
+ for (unsigned int i=0;i<NNN;i++){
+ // printf("%d %d\n",NNN,i);
+ Tet *t = allocator(myThread,i);
+ if (t->V[0]){
+ double vol;
+ if (tetQuality (t,&vol) < 0.3){
+ for (int j=0;j<6;j++){
+ int iLocalEdge = j;
+ Edge e (t->V[edges[iLocalEdge][0]], t->V[edges[iLocalEdge][1]]);
+ if (! embeddedEdges.find(e) && !ec.find(e)){
+ ec.addNewEdge(e);
+ bool result = edgeSwap(t, iLocalEdge, allocator, myThread);
+ if (result){
+ N++;
+ j=6;
+ break;
+ }
+ }
+ else{
+ // printf("cannot swop\n");
+ }
+ }
+ }
+ }
+ }
+ }
+ if (!N) break;
+ }
}
-*/
+
/*
Fixing a non star shaped cavity (non delaunay triangulations)
@@ -908,6 +1094,7 @@ static Tet* randomTet (int thread, tetContainer &allocator)
int isCavityCompatibleWithEmbeddedEdges(cavityContainer &cavity,
connContainer &bndK,
edgeContainer &allEmbeddedEdges){
+ // return true;
const unsigned int bSize = bndK.size();
std::vector<Edge> ed;
diff --git a/Mesh/delaunay3d_private.h b/Mesh/delaunay3d_private.h
index 45db65b4c2e7a42add9ef41f1a28f8f4e501fa42..f1b82ffa4e313c961ba99f947861738add12ab70 100644
--- a/Mesh/delaunay3d_private.h
+++ b/Mesh/delaunay3d_private.h
@@ -20,13 +20,17 @@
#endif
typedef unsigned char CHECKTYPE;
+struct Tet;
struct Vert {
private :
double _x[3];
double _lc;
unsigned int _num;
+ Tet *_t;
public :
+ inline void setT(Tet*t) { _t = t;}
+ inline Tet* getT() const { return _t;}
inline unsigned int getNum () const {return _num;}
inline void setNum (unsigned int n) {_num=n;}
unsigned char _thread;
@@ -40,7 +44,7 @@ public :
inline double &lc() {return _lc;}
inline operator double *() { return _x; }
Vert (double X=0, double Y=0, double Z=0, double lc=0, int num = 0)
- : _num(num), _thread(0)
+ : _num(num),_t(NULL), _thread(0)
{
_x[0] = X; _x[1] = Y; _x[2] = Z; _lc = lc;
}
@@ -223,6 +227,7 @@ struct Tet {
{
_modified=true;
V[0] = v0; V[1] = v1; V[2] = v2; V[3] = v3;
+ for (int i=0;i<4;i++)if (V[i])V[i]->setT(this);
// for (int i=0;i<4;i++)_copy[i] = *V[i];
return 1;
}
@@ -232,6 +237,7 @@ struct Tet {
double val = robustPredicates::orient3d((double*)v0, (double*)v1,
(double*)v2, (double*)v3);
V[0] = v0; V[1] = v1; V[2] = v2; V[3] = v3;
+ for (int i=0;i<4;i++)if (V[i])V[i]->setT(this);
if (val > 0){
// for (int i=0;i<4;i++)_copy[i] = *V[i];
return 1;
@@ -397,6 +403,7 @@ void delaunayTrgl(const unsigned int numThreads,
unsigned int Npts,
std::vector<Vert*> assignTo[],
tetContainer &allocator, edgeContainer *embedded = 0);
-bool edgeSwap(Tet *tet, int iLocalEdge, tetContainer &T, int myThread);
+void edgeSwapPass (int numThreads, tetContainer &allocator, edgeContainer &embeddedEdges);
+void vertexRelocationPass (int numThreads, std::vector<Vert*> &v);
#endif
diff --git a/Mesh/delaunay_refinement.cpp b/Mesh/delaunay_refinement.cpp
index 93ed5524402eecacdc2e7cc3216be9b4f202dfed..1b0bd9b5a7999d3ede129e923c3244bb52c5f186 100644
--- a/Mesh/delaunay_refinement.cpp
+++ b/Mesh/delaunay_refinement.cpp
@@ -142,20 +142,10 @@ void saturateEdge(Edge &e, std::vector<Vert*> &S, std::stack<IPT> &temp)
bgm = fields->get(fields->getBackgroundField());
}
- if (0){
- const double length = p1.distance(p2);
- const double lc_avg = (e.first->lc() + e.second->lc())*.5;
- const double l_adim_approx = length / lc_avg;
- if (l_adim_approx > 10.){
- SPoint3 pmid = (p1+p2)*0.5;
- const double lc = GMSHSIZE(pmid, bgm, lc_avg);
- S.push_back(new Vert(pmid.x(),pmid.y(),pmid.z(),lc));
- return;
- }
- }
-
// double _a = Cpu();
+ // printf("%g %g \n",e.first->lc(), e.second->lc());
+
const double dN = adaptiveTrapezoidalRule(p1,p2,e.first->lc(), e.second->lc(),
_result, dl, temp, bgm);
// _C1 += Cpu() - _a;
@@ -449,15 +439,25 @@ void computeMeshSizes (GRegion *gr, std::map<MVertex*, double> &s){
updateSize ((*it)->triangles[i]->getVertex((j+1)%3), s, d);
}
}
- }
+ }
+ // there may be steiner points
+ for (unsigned int i=0;i<gr->tetrahedra.size();i++){
+ for (unsigned int j = 0; j < 4; j++){
+ MVertex *v = gr->tetrahedra[i]->getVertex(j);
+ if (s.find(v) == s.end()){
+ s[v] = 1.0;
+ }
+ }
+ }
}
+extern double tetQuality (Tet* t, double *volume);
+
void edgeBasedRefinement(const int numThreads,
const int nptsatonce,
GRegion *gr)
{
// fill up old Datastructures
-
edgeContainer embeddedEdges (10000);
std::map<MVertex*, double> sizes;
@@ -525,7 +525,6 @@ void edgeBasedRefinement(const int numThreads,
}
}
-
// do not allow to saturate boundary edges
{
for (unsigned int i=0;i< allocator.size(0);i++) {
@@ -550,7 +549,6 @@ void edgeBasedRefinement(const int numThreads,
for (unsigned int i=0;i<_vertices.size();i++){
_filter.insert( _vertices[i] );
}
-
int iter = 1;
Msg::Info("------------------------------------- SATUR FILTR SORTH DELNY TIME TETS");
@@ -561,10 +559,10 @@ void edgeBasedRefinement(const int numThreads,
double t1 = Cpu();
// C_COUNT = 0;
saturateEdges (ec, allocator, numThreads, add);
- // printf("%d calls %d points added",C_COUNT,add.size());
+ // printf("%d points added",add.size());
double t2 = Cpu();
filterVertices (numThreads, _filter, add);
- // printf(" %d points remain\n",add.size());
+ // printf(" %d points remain\n",add.size());
double t3 = Cpu();
if (add.empty())break;
// randomize vertices (EXTREMELY IMPORTANT FOR NOT DETERIORATING PERFORMANCE)
@@ -572,6 +570,7 @@ void edgeBasedRefinement(const int numThreads,
// sort them using BRIO
std::vector<int> indices;
SortHilbert(add, indices);
+ // printf("%d indices\n",indices.size());
double t4 = Cpu();
delaunayTrgl (1,1,add.size(), &add,allocator,embeddedEdges.empty() ? NULL : &embeddedEdges);
double t5 = Cpu();
@@ -588,12 +587,45 @@ void edgeBasedRefinement(const int numThreads,
}
}
-
-
+ std::vector<Vert*> vv;
+ for (int myThread=0; myThread < numThreads; myThread++)
+ for (unsigned int i=0;i<allocator.size(myThread);i++)
+ for (unsigned int j=0;j<4;j++)
+ if (allocator(myThread,i)->V[j])
+ allocator(myThread,i)->V[j]->setNum(0);
+
+ for (int myThread=0; myThread < numThreads; myThread++)
+ for (unsigned int i=0;i<allocator.size(myThread);i++)
+ for (unsigned int j=0;j<4;j++)
+ if (allocator(myThread,i)->V[j] && allocator(myThread,i)->V[j]->getNum() == 0){
+ allocator(myThread,i)->V[j]->setNum(1);
+ std::map<Vert*,MVertex*>::iterator it = _ma.find(allocator(myThread,i)->V[j]);
+ if (it == _ma.end())vv.push_back(allocator(myThread,i)->V[j]);
+ }
+
+ double t6 = Cpu();
+ Msg::Info("Optimizing");
+ edgeSwapPass (numThreads, allocator, embeddedEdges);
+ vertexRelocationPass (numThreads,vv);
+ edgeSwapPass (numThreads, allocator, embeddedEdges);
+ vertexRelocationPass (numThreads,vv);
+ edgeSwapPass (numThreads, allocator, embeddedEdges);
+ vertexRelocationPass (numThreads,vv);
+ double t7 = Cpu();
+
+ Msg::Info("Optimization done (%g seconds)",t7-t6);
+
+ int cat [10] = {0,0,0,0,0,0,0,0,0,0};
+ double MIN = 1.0;
+
for (unsigned int i=0; i< allocator.size(0);i++){
Tet *tt = allocator (0,i);
MVertex *mvs[4];
if (tt->V[0]){
+ double vol;
+ double q = tetQuality (tt,&vol);
+ cat [(int)(q*10)] ++;
+ MIN = std::min(MIN,q);
for (int j=0;j<4;j++){
Vert *v = tt->V[j];
std::map<Vert*,MVertex*>::iterator it = _ma.find(v);
@@ -609,6 +641,11 @@ void edgeBasedRefinement(const int numThreads,
}
}
+ Msg::Info ("Min Tet Quality %22.15E",MIN);
+ for (int i=0;i<10;i++){
+ Msg::Info ("Tet Quality [%4.3f,%4.3f] %8d",0.1*i,0.1*(i+1),cat[i]);
+ }
+
if (Msg::GetVerbosity() == 99) {
std::map<Edge,double> _sizes;
for (unsigned int i=0; i< allocator.size(0);i++){
@@ -638,4 +675,6 @@ void edgeBasedRefinement(const int numThreads,
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/meshGRegionBoundaryRecovery.cpp b/Mesh/meshGRegionBoundaryRecovery.cpp
index 0a5ec9e1b400ac045dbd8f68cce52a8d8ed30710..469b6e638878c1db3a26e3a107e5cb14a829dc2d 100644
--- a/Mesh/meshGRegionBoundaryRecovery.cpp
+++ b/Mesh/meshGRegionBoundaryRecovery.cpp
@@ -107,7 +107,7 @@ bool tetgenmesh::reconstructmesh(void *p)
double t_start = Cpu();
std::vector<MVertex*> _vertices;
-
+ std::map<int,MVertex*> _extras;
// Get the set of vertices from GRegion.
{
std::set<MVertex*, MVertexLessThanNum> all;
@@ -526,6 +526,7 @@ bool tetgenmesh::reconstructmesh(void *p)
clock_t t;
// printf("coucou2\n");
+ Msg::Info("Boundary Recovery...");
recoverboundary(t);
// printf("coucou3\n");
@@ -651,7 +652,7 @@ bool tetgenmesh::reconstructmesh(void *p)
}
v->setIndex(pointmark(pointloop));
_gr->mesh_vertices.push_back(v);
- _vertices.push_back(v);
+ _extras[pointmark(pointloop)] = v;
}
spivot(parentseg, parentsh);
if (parentsh.sh != NULL) {
@@ -675,7 +676,7 @@ bool tetgenmesh::reconstructmesh(void *p)
}
v->setIndex(pointmark(pointloop));
_gr->mesh_vertices.push_back(v);
- _vertices.push_back(v);
+ _extras[pointmark(pointloop)] = v;
}
}
// Record all the GFaces' tag at this segment.
@@ -690,9 +691,8 @@ bool tetgenmesh::reconstructmesh(void *p)
// Create an interior mesh vertex.
MVertex *v = new MVertex(pointloop[0], pointloop[1], pointloop[2], _gr);
v->setIndex(pointmark(pointloop));
- _gr->mesh_vertices.push_back(v);
- _vertices.push_back(v);
- }
+ _extras[pointmark(pointloop)] = v;
+ _gr->mesh_vertices.push_back(v); }
}
else if (pointtype(pointloop) == FREEFACETVERTEX) {
sdecode(point2sh(pointloop), parentsh);
@@ -718,29 +718,34 @@ bool tetgenmesh::reconstructmesh(void *p)
}
v->setIndex(pointmark(pointloop));
_gr->mesh_vertices.push_back(v);
- _vertices.push_back(v);
+ _extras[pointmark(pointloop)] = v;
}
else {
// Create a mesh vertex.
MVertex *v = new MVertex(pointloop[0], pointloop[1], pointloop[2], _gr);
v->setIndex(pointmark(pointloop));
_gr->mesh_vertices.push_back(v);
- _vertices.push_back(v);
+ _extras[pointmark(pointloop)] = v;
}
}
else {
MVertex *v = new MVertex(pointloop[0], pointloop[1], pointloop[2], _gr);
v->setIndex(pointmark(pointloop));
_gr->mesh_vertices.push_back(v);
- _vertices.push_back(v);
+ _extras[pointmark(pointloop)] = v;
}
}
pointloop = pointtraverse();
}
- assert((int)_vertices.size() == points->items);
+ // assert((int)_vertices.size() == points->items);
}
+ if (!_extras.empty())Msg::Info("We add %d steiner points...",_extras.size());
+
+
if (l_edges.size() > 0) {
+
+
// There are Steiner points on segments!
face segloop;
// Re-create the segment mesh in the corresponding GEdges.
@@ -756,7 +761,7 @@ bool tetgenmesh::reconstructmesh(void *p)
}
}
assert(ge != NULL);
- Msg::Info("Steiner points exist on GEdge %d", ge->tag());
+ // Msg::Info("Steiner points exist on GEdge %d", ge->tag());
// Delete the old triangles.
for(unsigned int i = 0; i < ge->lines.size(); i++)
delete ge->lines[i];
@@ -770,8 +775,11 @@ bool tetgenmesh::reconstructmesh(void *p)
if (shellmark(segloop) == etag) {
p[0] = sorg(segloop);
p[1] = sdest(segloop);
- MVertex *v1 = _vertices[pointmark(p[0])];
- MVertex *v2 = _vertices[pointmark(p[1])];
+ MVertex *v1 = pointmark(p[0]) >= _vertices.size() ? _extras[pointmark(p[0])] : _vertices[pointmark(p[0])];
+ MVertex *v2 = pointmark(p[1]) >= _vertices.size() ? _extras[pointmark(p[1])] : _vertices[pointmark(p[1])];
+ // MVertex *v2 = _vertices[pointmark(p[1])];
+ // printf("%d %d %d\n",pointmark(p[0]),pointmark(p[1]),_vertices.size());
+ // printf("%g %g %g %g\n",v1->x(),v1->y(),v2->x(),v2->y());
MLine *t = new MLine(v1, v2);
ge->lines.push_back(t);
}
@@ -814,9 +822,15 @@ bool tetgenmesh::reconstructmesh(void *p)
p[0] = sorg(subloop);
p[1] = sdest(subloop);
p[2] = sapex(subloop);
- MVertex *v1 = _vertices[pointmark(p[0])];
- MVertex *v2 = _vertices[pointmark(p[1])];
- MVertex *v3 = _vertices[pointmark(p[2])];
+ MVertex *v1 = pointmark(p[0]) >= _vertices.size() ? _extras[pointmark(p[0])] : _vertices[pointmark(p[0])];
+ MVertex *v2 = pointmark(p[1]) >= _vertices.size() ? _extras[pointmark(p[1])] : _vertices[pointmark(p[1])];
+ MVertex *v3 = pointmark(p[2]) >= _vertices.size() ? _extras[pointmark(p[2])] : _vertices[pointmark(p[2])];
+ if (pointmark(p[0]) >= _vertices.size())printf("F %d %d\n",v1->getIndex(),pointmark(p[0]));
+ if (pointmark(p[1]) >= _vertices.size())printf("F %d %d\n",v2->getIndex(),pointmark(p[1]));
+ if (pointmark(p[2]) >= _vertices.size())printf("F %d %d\n",v3->getIndex(),pointmark(p[2]));
+ // MVertex *v1 = _vertices[pointmark(p[0])];
+ // MVertex *v2 = _vertices[pointmark(p[1])];
+ // MVertex *v3 = _vertices[pointmark(p[2])];
MTriangle *t = new MTriangle(v1, v2, v3);
gf->triangles.push_back(t);
}
@@ -837,10 +851,20 @@ bool tetgenmesh::reconstructmesh(void *p)
p[2] = apex(tetloop);
p[3] = oppo(tetloop);
- MVertex *v1 = _vertices[pointmark(p[0])];
- MVertex *v2 = _vertices[pointmark(p[1])];
- MVertex *v3 = _vertices[pointmark(p[2])];
- MVertex *v4 = _vertices[pointmark(p[3])];
+ MVertex *v1 = pointmark(p[0]) >= _vertices.size() ? _extras[pointmark(p[0])] : _vertices[pointmark(p[0])];
+ MVertex *v2 = pointmark(p[1]) >= _vertices.size() ? _extras[pointmark(p[1])] : _vertices[pointmark(p[1])];
+ MVertex *v3 = pointmark(p[2]) >= _vertices.size() ? _extras[pointmark(p[2])] : _vertices[pointmark(p[2])];
+ MVertex *v4 = pointmark(p[3]) >= _vertices.size() ? _extras[pointmark(p[3])] : _vertices[pointmark(p[3])];
+
+ // if (pointmark(p[0]) >= _vertices.size())printf("%d %d\n",v1->getIndex(),pointmark(p[0]));
+ // if (pointmark(p[1]) >= _vertices.size())printf("%d %d\n",v2->getIndex(),pointmark(p[1]));
+ // if (pointmark(p[2]) >= _vertices.size())printf("%d %d\n",v3->getIndex(),pointmark(p[2]));
+ // if (pointmark(p[3]) >= _vertices.size())printf("%d %d\n",v4->getIndex(),pointmark(p[3]));
+
+ // MVertex *v1 = _vertices[pointmark(p[0])];
+ // MVertex *v2 = _vertices[pointmark(p[1])];
+ // MVertex *v3 = _vertices[pointmark(p[2])];
+ // MVertex *v4 = _vertices[pointmark(p[3])];
MTetrahedron *t = new MTetrahedron(v1, v2, v3, v4);
_gr->tetrahedra.push_back(t);
tetloop.tet = tetrahedrontraverse();
diff --git a/Mesh/meshGRegionDelaunayInsertion.cpp b/Mesh/meshGRegionDelaunayInsertion.cpp
index 5e153903773efd1ea2f504cf333d77f115ba4079..84a0a1b7dfd88bdc80bf3c82ca8904009d9cac54 100644
--- a/Mesh/meshGRegionDelaunayInsertion.cpp
+++ b/Mesh/meshGRegionDelaunayInsertion.cpp
@@ -193,16 +193,18 @@ struct faceXtet{
}
};
-void connectTets_vector2(std::vector<MTet4*> &t, std::vector<faceXtet> &conn)
+template <class ITER>
+void connectTets_vector2_templ(size_t _size, ITER beg, ITER end, std::vector<faceXtet> &conn)
{
conn.clear();
- conn.reserve(4*t.size());
+ conn.reserve(4*_size);
// unsigned int k = 0;
// printf("COUCOU1 %d tets\n",t.size());
- for (unsigned int i=0;i<t.size();i++){
- if (!t[i]->isDeleted()){
+ for (ITER IT = beg; IT != end; ++IT){
+ MTet4* t = *IT;
+ if (!t->isDeleted()){
for (int j = 0; j < 4; j++){
- conn.push_back(faceXtet(t[i], j));
+ conn.push_back(faceXtet(t, j));
}
}
}
@@ -253,6 +255,8 @@ void connectTets(ITER beg, ITER end, std::set<MFace, Less_Face> *allEmbeddedFace
void connectTets(std::list<MTet4*> &l) { connectTets(l.begin(), l.end()); }
void connectTets(std::vector<MTet4*> &l) { connectTets(l.begin(), l.end()); }
+void connectTets_vector2(std::list<MTet4*> &l, std::vector<faceXtet> &conn) { connectTets_vector2_templ(l.size(), l.begin(), l.end(), conn); }
+void connectTets_vector2(std::vector<MTet4*> &l, std::vector<faceXtet> &conn) { connectTets_vector2_templ(l.size(), l.begin(), l.end(), conn); }
// Ensure the star-shapeness of the delaunay cavity
// We use the visibility criterion : the vertex should be visible
@@ -865,7 +869,7 @@ void optimizeMesh(GRegion *gr, const qmTetrahedron::Measures &qm)
while (1){
// printf("coucou\n");
std::vector<MTet4*> newTets;
- for (CONTAINER::iterator it = allTets.begin(); it != allTets.end(); ++it){
+ /* for (CONTAINER::iterator it = allTets.begin(); it != allTets.end(); ++it){
if (!(*it)->isDeleted()){
double qq = (*it)->getQuality();
if (qq < qMin){
@@ -878,6 +882,7 @@ void optimizeMesh(GRegion *gr, const qmTetrahedron::Measures &qm)
}
}
}
+ */
// printf("coucou\n");
illegals.clear();
@@ -1188,6 +1193,11 @@ void insertVerticesInRegion (GRegion *gr, int maxVert, bool _classify)
}
gr->tetrahedra.clear();
+ {
+ std::vector<faceXtet> conn;
+ // connectTets_vector2_templ(allTets.size(), allTets.begin(), allTets.end() , conn);
+ }
+ // SLOW
connectTets(allTets.begin(), allTets.end());
// classify the tets on the right region
diff --git a/Mesh/meshGRegionLocalMeshMod.cpp b/Mesh/meshGRegionLocalMeshMod.cpp
index 7e0073e4cafbdf566d61db2963614e1855f4bd08..dff9b341c5791e916446ad0107be2a53a546357e 100644
--- a/Mesh/meshGRegionLocalMeshMod.cpp
+++ b/Mesh/meshGRegionLocalMeshMod.cpp
@@ -128,10 +128,22 @@ void BuildSwapPattern3(SwapPattern *sc)
sc->trianguls = trgul ;
}
+/*
+ 0 1
+ +------------+
+ | |
+ | |
+ | |
+ +------------+
+ 3 2
+
+*/
+
+
void BuildSwapPattern4(SwapPattern *sc)
{
static int trgl[][3] =
- { {0,1,2}, {0,2,3}, {0,1,3}, {1,2,3} };
+ { {0,1,2}, {2,3,0}, {1,2,3}, {3,0,1} };
static int trgul[][5] =
{ {0,1,-1,-1,-1}, {2,3,-1,-1,-1} };
@@ -282,6 +294,7 @@ bool edgeSwap(std::vector<MTet4 *> &newTets,
// if there exist no swap that enhance the quality
if (best <= tetQualityRef) return false;
+ // does random swaps if (best < .01) return false;
// we have the best configuration, so we swap
// printf("outside size = %d\n",outside.size());
diff --git a/Mesh/meshGRegionRelocateVertex.cpp b/Mesh/meshGRegionRelocateVertex.cpp
index b7b7185eed13d11a8c6df6b89bd27d79ba4ea119..b4711b3f5a8bf0e36f241a9ea9b6daad17bc891d 100644
--- a/Mesh/meshGRegionRelocateVertex.cpp
+++ b/Mesh/meshGRegionRelocateVertex.cpp
@@ -60,6 +60,34 @@ static double objective_function (double xi, MVertex *ver, GFace *gf,
}
+static double objective_function (double xi, MVertex *ver, GFace *gf,
+ SPoint3 &p1, SPoint3 &p2,
+ const std::vector<MElement*> <){
+ double x = ver->x();
+ double y = ver->y();
+ double z = ver->z();
+ SPoint3 p = p1 * (1.-xi) + p2 * xi;
+
+ double initialGuess[2]={0,0};
+
+ GPoint pp = gf->closestPoint(p,initialGuess);
+ ver->x() = pp.x();
+ ver->y() = pp.y();
+ ver->z() = pp.z();
+ double minQual = 1.0;
+ for(unsigned int i = 0; i < lt.size(); i++){
+ if (lt[i]->getNumVertices () == 4)
+ minQual = std::min((lt[i]->etaShapeMeasure()), minQual);
+ else
+ minQual = std::min(fabs(lt[i]->gammaShapeMeasure()), minQual);
+ }
+ ver->x() = x;
+ ver->y() = y;
+ ver->z() = z;
+ return minQual;
+}
+
+
#define sqrt5 2.236067977499789696
@@ -158,6 +186,55 @@ double Maximize_Quality_Golden_Section( MVertex *ver, GFace *gf,
return a;
}
+
+double Maximize_Quality_Golden_Section( MVertex *ver, GFace *gf,
+ SPoint3 &p1, SPoint3 &p2,
+ const std::vector<MElement*> < ,
+ double tol, double &worst)
+{
+
+ static const double lambda = 0.5 * (sqrt5 - 1.0);
+ static const double mu = 0.5 * (3.0 - sqrt5); // = 1 - lambda
+ double a = 0.0;
+ double b = 1.0;
+
+ worst = objective_function (0.0, ver, gf, p1, p2, lt );
+
+ if (worst > 0.5) return 0.0;
+
+ double x1 = b - lambda * (b - a);
+ double x2 = a + lambda * (b - a);
+ double fx1 = objective_function (x1, ver, gf, p1, p2, lt );
+ double fx2 = objective_function (x2, ver, gf, p1, p2, lt );
+
+ 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 , tol) ) break;
+ x1 = x2;
+ fx1 = fx2;
+ x2 = b - mu * (b - a);
+ fx2 = objective_function (x2, ver, gf, p1, p2, lt );
+ } else {
+ b = x2;
+ if ( Stopping_Rule( a, b , tol) ) break;
+ x2 = x1;
+ fx2 = fx1;
+ x1 = a + mu * (b - a);
+ fx1 = objective_function (x1, ver, gf, p1, p2, lt );
+ }
+ }
+ double final = objective_function (a, ver, gf, p1, p2, lt );
+ if (final < worst) return 0.0;
+ worst = final;
+ // printf("finally : %g %g (%12.5E,%12.5E)\n",a,b,fa,fb);
+ return a;
+}
+
+
void _relocateVertexGolden(MVertex *ver,
const std::vector<MElement*> <, double relax , double tol)
{
@@ -184,6 +261,37 @@ void _relocateVertexGolden(MVertex *ver,
ver->z() = (1.-xi) * ver->z() + xi * z/N;
}
+
+// use real space + projection at the end
+static double _relocateVertex2(GFace* gf, MVertex *ver,
+ const std::vector<MElement*> <,
+ double tol) {
+
+ SPoint3 p1(0,0,0);
+ int counter = 0;
+ for(unsigned int i = 0; i < lt.size(); i++){
+ for (int j=0;j<lt[i]->getNumVertices();j++){
+ MVertex* v = lt[i]->getVertex(j);
+ p1 += SPoint3(v->x(),v->y(),v->z());
+ counter++;
+ }
+ }
+ p1 *= 1./(double)counter;
+ SPoint3 p2(ver->x(),ver->y(),ver->z());
+ double worst;
+ double xi = Maximize_Quality_Golden_Section( ver, gf, p1, p2, lt , tol, worst);
+
+ SPoint3 p = p1*(1-xi) + p2*xi;
+ double initialGuess[2]={0,0};
+ GPoint pp = gf->closestPoint(p,initialGuess);
+ if (!pp.succeeded())return 2.0;
+ ver->x() = pp.x();
+ ver->y() = pp.y();
+ ver->z() = pp.z();
+ return worst;
+
+}
+
static double _relocateVertex(GFace* gf, MVertex *ver,
const std::vector<MElement*> <,
double tol) {
@@ -191,15 +299,19 @@ static double _relocateVertex(GFace* gf, MVertex *ver,
SPoint2 p1(0,0);
SPoint2 p2;
- ver->getParameter(0,p2[0]);
- ver->getParameter(1,p2[1]);
+ if (ver->getParameter(0,p2[0])){
+ ver->getParameter(1,p2[1]);
+ }
+ else {
+ return _relocateVertex2(gf,ver,lt,tol);
+ }
int counter=0;
for(unsigned int i = 0; i < lt.size(); i++){
for (int j=0;j<lt[i]->getNumVertices();j++){
MVertex* v = lt[i]->getVertex(j);
SPoint2 pp;
- reparamMeshVertexOnFace(v, gf, pp);
+ reparamMeshVertexOnFace(v, gf, pp);
counter++;
if (v->onWhat()->dim() == 1) {
GEdge *ge = dynamic_cast<GEdge*> (v->onWhat());
@@ -230,7 +342,7 @@ void getAllBoundaryLayerVertices (GFace *gf, std::set<MVertex*> &vs);
void RelocateVertices (GFace* gf, int niter, double tol) {
std::set<MVertex*> vs;
getAllBoundaryLayerVertices (gf, vs);
-
+
v2t_cont adj;
buildVertexToElement(gf->triangles, adj);
buildVertexToElement(gf->quadrangles, adj);
diff --git a/Mesh/surfaceFiller.cpp b/Mesh/surfaceFiller.cpp
index 65bc9c890752094d2480312ed83bde18809734d2..1697072a807c9f61529100ca3d10e97ef4df8647 100644
--- a/Mesh/surfaceFiller.cpp
+++ b/Mesh/surfaceFiller.cpp
@@ -691,11 +691,11 @@ void packingOfParallelogramsSmoothness(GFace* gf, std::vector<MVertex*> &packed
// add the vertices as additional vertices in the
// surface mesh
- char ccc[256]; sprintf(ccc,"points%d.pos",gf->tag());
- FILE *f = Fopen(ccc,"w");
- if(f) fprintf(f, "View \"\"{\n");
+ // char ccc[256]; sprintf(ccc,"points%d.pos",gf->tag());
+ // FILE *f = Fopen(ccc,"w");
+ // if(f) fprintf(f, "View \"\"{\n");
for (unsigned int i=0;i<vertices.size();i++){
- if(f) vertices[i]->print(f,i);
+ // if(f) vertices[i]->print(f,i);
if(vertices[i]->_v->onWhat() == gf) {
packed.push_back(vertices[i]->_v);
metrics.push_back(vertices[i]->_meshMetric);
@@ -704,10 +704,10 @@ void packingOfParallelogramsSmoothness(GFace* gf, std::vector<MVertex*> &packed
}
delete vertices[i];
}
- if(f){
- fprintf(f,"};");
- fclose(f);
- }
+ // if(f){
+ // fprintf(f,"};");
+ // fclose(f);
+ // }
}