From 1e7f858527e61c1666ea04010c200e0e548d3627 Mon Sep 17 00:00:00 2001
From: Jean-Francois Remacle <jean-francois.remacle@uclouvain.be>
Date: Fri, 12 Apr 2013 11:12:11 +0000
Subject: [PATCH] Big modif in the 3D mesher (2 be TESTED) Removed the -q
option in initial mesh generation (generates meshes that do not respect size
field)
---
Common/CommandLine.cpp | 2 +
Common/GmshDefines.h | 1 +
Mesh/meshGFace.cpp | 8 +-
Mesh/meshGFaceBoundaryLayers.cpp | 4 +-
Mesh/meshGFaceDelaunayInsertion.cpp | 57 ++++++-
Mesh/meshGFaceDelaunayInsertion.h | 5 +-
Mesh/meshGRegion.cpp | 5 +-
Mesh/meshGRegionDelaunayInsertion.cpp | 237 ++++++++++++++++++++++----
benchmarks/2d/Square-01.geo | 2 +-
benchmarks/2d/conge.geo | 2 +-
10 files changed, 281 insertions(+), 42 deletions(-)
diff --git a/Common/CommandLine.cpp b/Common/CommandLine.cpp
index 2decad5527..56ad301f4e 100644
--- a/Common/CommandLine.cpp
+++ b/Common/CommandLine.cpp
@@ -793,6 +793,8 @@ void GetOptions(int argc, char *argv[])
CTX::instance()->mesh.algo2d = ALGO_2D_FRONTAL_QUAD;
else if(!strncmp(argv[i], "pack", 4))
CTX::instance()->mesh.algo2d = ALGO_2D_PACK_PRLGRMS;
+ else if(!strncmp(argv[i], "ruppert", 7))
+ CTX::instance()->mesh.algo2d = ALGO_2D_RUPPERT;
else if(!strncmp(argv[i], "front2d", 7) || !strncmp(argv[i], "frontal", 7))
CTX::instance()->mesh.algo2d = ALGO_2D_FRONTAL;
else if(!strncmp(argv[i], "bamg",4))
diff --git a/Common/GmshDefines.h b/Common/GmshDefines.h
index ced9182544..66f24c25f5 100644
--- a/Common/GmshDefines.h
+++ b/Common/GmshDefines.h
@@ -221,6 +221,7 @@
#define ALGO_2D_BAMG 7
#define ALGO_2D_FRONTAL_QUAD 8
#define ALGO_2D_PACK_PRLGRMS 9
+#define ALGO_2D_RUPPERT 10
// 3D meshing algorithms (numbers should not be changed)
#define ALGO_3D_DELAUNAY 1
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 548b445be1..2f6210170b 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -455,6 +455,7 @@ static bool algoDelaunay2D(GFace *gf)
if(gf->getMeshingAlgo() == ALGO_2D_DELAUNAY ||
gf->getMeshingAlgo() == ALGO_2D_BAMG ||
gf->getMeshingAlgo() == ALGO_2D_FRONTAL ||
+ gf->getMeshingAlgo() == ALGO_2D_RUPPERT ||
gf->getMeshingAlgo() == ALGO_2D_FRONTAL_QUAD ||
gf->getMeshingAlgo() == ALGO_2D_PACK_PRLGRMS ||
gf->getMeshingAlgo() == ALGO_2D_BAMG)
@@ -674,7 +675,7 @@ void modifyInitialMeshForTakingIntoAccountBoundaryLayers(GFace *gf)
MEdge dv2 (v21,v22);
//avoid convergent errors
- if (dv2.length() < 0.3 * dv.length())break;
+ if (dv2.length() < 0.03 * dv.length())break;
blQuads.push_back(new MQuadrangle(v11,v21,v22,v12));
fprintf(ff2,"SQ (%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){1,1,1,1};\n",
v11->x(),v11->y(),v11->z(),
@@ -1297,6 +1298,8 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
else if(gf->getMeshingAlgo() == ALGO_2D_DELAUNAY ||
gf->getMeshingAlgo() == ALGO_2D_AUTO)
bowyerWatson(gf);
+ else if(gf->getMeshingAlgo() == ALGO_2D_RUPPERT)
+ gmshRuppert(gf,.5);
else {
bowyerWatson(gf,15000);
meshGFaceBamg(gf);
@@ -2024,6 +2027,8 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
else if(gf->getMeshingAlgo() == ALGO_2D_DELAUNAY ||
gf->getMeshingAlgo() == ALGO_2D_AUTO)
bowyerWatson(gf,1000000000, &equivalence, ¶metricCoordinates);
+ else if(gf->getMeshingAlgo() == ALGO_2D_RUPPERT)
+ gmshRuppert(gf,.4,10000, &equivalence, ¶metricCoordinates);
else
meshGFaceBamg(gf);
if (!infty || !(CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine))
@@ -2105,6 +2110,7 @@ void meshGFace::operator() (GFace *gf, bool print)
case ALGO_2D_DELAUNAY : algo = "Delaunay"; break;
case ALGO_2D_MESHADAPT_OLD : algo = "MeshAdapt (old)"; break;
case ALGO_2D_BAMG : algo = "Bamg"; break;
+ case ALGO_2D_RUPPERT : algo = "Ruppert"; break;
case ALGO_2D_PACK_PRLGRMS : algo = "Square Packing"; break;
case ALGO_2D_AUTO :
algo = (gf->geomType() == GEntity::Plane) ? "Delaunay" : "MeshAdapt";
diff --git a/Mesh/meshGFaceBoundaryLayers.cpp b/Mesh/meshGFaceBoundaryLayers.cpp
index 8800821807..fc8850ef30 100644
--- a/Mesh/meshGFaceBoundaryLayers.cpp
+++ b/Mesh/meshGFaceBoundaryLayers.cpp
@@ -332,10 +332,10 @@ BoundaryLayerColumns* buildAdditionalPoints2D (GFace *gf)
break;
}
// printf("%g %g %g \n",current->x(),current->y(),blf->current_distance);
- if (blf->current_closest != catt || blf -> current_distance < _current_distance){
+ if (0 && blf->current_closest != catt || blf -> current_distance < _current_distance){
SVector3 aaa (_close- blf->_closest_point);
if (aaa.norm() > 8*blf->hwall_n || blf -> current_distance < _current_distance){
- // printf("reaching the skelton %d\n", (int) _column.size());
+ printf("reaching the skelton %d %g %g\n", (int) _column.size(), aaa.norm(),blf->hwall_n);
delete _column[_column.size()-1];
_column.erase(--_column.end());
_metrics.erase(--_metrics.end());
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index b1d210bf92..f37ea93d76 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -316,7 +316,10 @@ MTri3::MTri3(MTriangle *t, double lc, SMetric3 *metric, bidimMeshData * data, GF
{base->getVertex(2)->x(), base->getVertex(2)->y(), base->getVertex(2)->z()};
if (!metric){
- if (radiusNorm == 2){
+ if (radiusNorm == 3){
+ circum_radius = 1./base->gammaShapeMeasure();
+ }
+ else if (radiusNorm == 2){
circumCenterXYZ(pa, pb, pc, center);
const double dx = base->getVertex(0)->x() - center[0];
const double dy = base->getVertex(0)->y() - center[1];
@@ -867,6 +870,56 @@ static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it
}
}
+void gmshRuppert(GFace *gf, double minqual, int MAXPNT,
+ std::map<MVertex* , MVertex*>* equivalence,
+ std::map<MVertex*, SPoint2> * parametricCoordinates){
+ MTri3::radiusNorm =3;
+
+ std::set<MTri3*,compareTri3Ptr> AllTris;
+ bidimMeshData DATA(equivalence,parametricCoordinates);
+
+ buildMeshGenerationDataStructures(gf, AllTris, DATA);
+
+ int nbSwaps = edgeSwapPass(gf, AllTris, SWCR_DEL, DATA);
+
+ Msg::Debug("Delaunization of the initial mesh done (%d swaps)", nbSwaps);
+
+ if(AllTris.empty()){
+ Msg::Error("No triangles in initial mesh");
+ return;
+ }
+
+ int ITER = 0;
+ int NBDELETED = 0;
+ while (1){
+ MTri3 *worst = *AllTris.begin();
+ if (worst->isDeleted()){
+ delete worst->tri();
+ delete worst;
+ AllTris.erase(AllTris.begin());
+ NBDELETED ++;
+ }
+ else{
+ double center[2],metric[3],r2;
+ // printf("%12.5E\n",worst->getRadius() );
+ if (1./worst->getRadius() > minqual || (int) DATA.vSizes.size() > MAXPNT) break;
+ circUV(worst->tri(), DATA, center, gf);
+ MTriangle *base = worst->tri();
+ int index0 = DATA.getIndex( base->getVertex(0) );
+ int index1 = DATA.getIndex( base->getVertex(1) );
+ int index2 = DATA.getIndex( base->getVertex(2) );
+ double pa[2] = {(DATA.Us[index0] + DATA.Us[index1] + DATA.Us[index2])/ 3.,
+ (DATA.Vs[index0] + DATA.Vs[index1] + DATA.Vs[index2])/ 3.};
+ buildMetric(gf, pa, metric);
+ circumCenterMetric(worst->tri(), metric, DATA, center, r2);
+ insertAPoint(gf, AllTris.begin(), center, metric, DATA, AllTris);
+ }
+ }
+ MTri3::radiusNorm =2;
+ transferDataStructure(gf, AllTris, DATA);
+}
+
+
void bowyerWatson(GFace *gf, int MAXPNT,
std::map<MVertex* , MVertex*>* equivalence,
std::map<MVertex*, SPoint2> * parametricCoordinates)
@@ -1357,7 +1410,7 @@ void buildBackGroundMesh (GFace *gf,
int CurvControl = CTX::instance()->mesh.lcFromCurvature;
CTX::instance()->mesh.lcFromCurvature = 0;
// Do a background mesh
- bowyerWatson(gf,4000, equivalence,parametricCoordinates);
+ gmshRuppert(gf,0.3,4000, equivalence,parametricCoordinates);
// Re-enable curv control if asked
CTX::instance()->mesh.lcFromCurvature = CurvControl;
// apply this to the BGM
diff --git a/Mesh/meshGFaceDelaunayInsertion.h b/Mesh/meshGFaceDelaunayInsertion.h
index 96590d8c08..a07c34a0e2 100644
--- a/Mesh/meshGFaceDelaunayInsertion.h
+++ b/Mesh/meshGFaceDelaunayInsertion.h
@@ -84,7 +84,7 @@ class MTri3
MTri3 *neigh[3];
public :
- static int radiusNorm; // 2 is euclidian norm, -1 is infinite norm
+ static int radiusNorm; // 2 is euclidian norm, -1 is infinite norm , 3 quality
bool isDeleted() const { return deleted; }
void forceRadius(double r) { circum_radius = r; }
inline double getRadius() const { return circum_radius; }
@@ -133,6 +133,9 @@ class compareTri3Ptr
void connectTriangles(std::list<MTri3*> &);
void connectTriangles(std::vector<MTri3*> &);
void connectTriangles(std::set<MTri3*,compareTri3Ptr> &AllTris);
+void gmshRuppert(GFace *gf, double minqual = 0.2, int MAXPNT= 1000000000,
+ std::map<MVertex* , MVertex*>* equivalence = 0,
+ std::map<MVertex*, SPoint2> * parametricCoordinates = 0);
void bowyerWatson(GFace *gf, int MAXPNT= 1000000000,
std::map<MVertex* , MVertex*>* equivalence= 0,
std::map<MVertex*, SPoint2> * parametricCoordinates= 0);
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 97eadb0633..cbe2ba8c2f 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -564,9 +564,12 @@ void MeshDelaunayVolume(std::vector<GRegion*> ®ions)
(Msg::GetVerbosity() > 6) ? 'V': '\0');
}
else {
- sprintf(opts, "-q%gYpe%c", CTX::instance()->mesh.delaunayQ,
+ sprintf(opts, "-Ype%c",
(Msg::GetVerbosity() < 3) ? 'Q':
(Msg::GetVerbosity() > 6) ? 'V': '\0');
+ /* sprintf(opts, "-q%gYpe%c", CTX::instance()->mesh.delaunayQ,
+ (Msg::GetVerbosity() < 3) ? 'Q':
+ (Msg::GetVerbosity() > 6) ? 'V': '\0');*/
}
try{
tetrahedralize(opts, &in, &out);
diff --git a/Mesh/meshGRegionDelaunayInsertion.cpp b/Mesh/meshGRegionDelaunayInsertion.cpp
index b1298f6367..df4d1d3520 100644
--- a/Mesh/meshGRegionDelaunayInsertion.cpp
+++ b/Mesh/meshGRegionDelaunayInsertion.cpp
@@ -53,7 +53,7 @@ int MTet4::inCircumSphere(const double *p) const
return (result > 0) ? 1 : 0;
}
-static int faces[4][3] = {{0,1,2}, {0,2,3}, {0,1,3}, {1,2,3}};
+static int faces[4][3] = {{0,1,2}, {0,2,3}, {0,3,1}, {1,3,2}};
struct faceXtet{
MVertex *v[3];
@@ -66,7 +66,10 @@ struct faceXtet{
v[2] = t1->tet()->getVertex(faces[iFac][2]);
std::sort(v, v + 3);
}
- inline bool operator < (const faceXtet & other) const
+
+ inline MVertex * getVertex (int i) const { return t1->tet()->getVertex(faces[i1][i]);}
+
+ inline bool operator < (const faceXtet & other) const
{
if (v[0] < other.v[0]) return true;
if (v[0] > other.v[0]) return false;
@@ -81,6 +84,17 @@ struct faceXtet{
v[1] == other.v[1] &&
v[2] == other.v[2] );
}
+ bool visible (MVertex *v){
+ MVertex* v0 = t1->tet()->getVertex(faces[i1][0]);
+ MVertex* v1 = t1->tet()->getVertex(faces[i1][1]);
+ MVertex* v2 = t1->tet()->getVertex(faces[i1][2]);
+ double a[3] = {v0->x(),v0->y(),v0->z()};
+ double b[3] = {v1->x(),v1->y(),v1->z()};
+ double c[3] = {v2->x(),v2->y(),v2->z()};
+ double d[3] = {v->x(),v->y(),v->z()};
+ double o = robustPredicates :: orient3d(a,b,c,d);
+ return o < 0;
+ }
};
template <class ITER>
@@ -161,6 +175,106 @@ static bool isActive(MTet4 *t, double limit_, int &i, std::set<MFace,Less_Face>
void connectTets(std::list<MTet4*> &l) { connectTets(l.begin(), l.end()); }
void connectTets(std::vector<MTet4*> &l) { connectTets(l.begin(), l.end()); }
+// Ensure the star-shapeness of the delaunay cavity
+// We use the visibility criterion : the vertex should be visible
+// by all the facets of the cavity
+
+static void removeFromCavity (std::list<faceXtet> & shell,
+ std::list<MTet4*> & cavity,
+ faceXtet &toRemove)
+{
+ toRemove.t1->setDeleted(false);
+ cavity.erase(std::remove_if(cavity.begin(),cavity.end(),
+ std::bind2nd(std::equal_to<MTet4*>(), toRemove.t1)));
+ for (int i=0;i<4;i++){
+ faceXtet fxt2(toRemove.t1,i);
+ std::list<faceXtet>::iterator it = std::find(shell.begin(),shell.end(),fxt2);
+ if (it == shell.end()){
+ MTet4 *opposite = toRemove.t1->getNeigh(toRemove.i1);
+ for (int j=0;j<4;j++){
+ faceXtet fxt3(opposite,j);
+ if (fxt3 == fxt2){
+ shell.push_back(fxt3);
+ }
+ }
+ }
+ else shell.erase(it);
+ }
+}
+
+static void extendCavity (std::list<faceXtet> & shell,
+ std::list<MTet4*> & cavity,
+ faceXtet &toExtend)
+{
+ MTet4 *t = toExtend.t1;
+ MTet4 *opposite = t->getNeigh(toExtend.i1);
+ for (int i=0;i<4;i++){
+ faceXtet fxt(opposite,i);
+ std::list<faceXtet>::iterator it = std::find(shell.begin(),shell.end(),fxt);
+ if (it == shell.end()) shell.push_back(fxt);
+ else shell.erase(it);
+ }
+ cavity.push_back(opposite);
+ opposite->setDeleted(true);
+}
+
+// if all faces of the tet that are not in the shell see v, then it is ok
+// either to add or to remove t from the shell
+static bool verifyShell (MVertex *v, MTet4*t, std::list<faceXtet> & shell){
+ if (!t)return false;
+ return 1;
+ int NBAD_BEFORE=0,NBAD_AFTER=0;
+ for (int i=0;i<4 ; i++){
+ faceXtet fxt(t,i);
+ bool starShaped = fxt.visible(v);
+ if (!starShaped){
+ std::list<faceXtet>::iterator its = std::find(shell.begin(),shell.end(),fxt);
+ if (its == shell.end())NBAD_AFTER ++;
+ else NBAD_BEFORE++;
+ }
+ }
+ return 1;
+ return (NBAD_AFTER < NBAD_BEFORE);
+}
+
+int makeCavityStarShaped (std::list<faceXtet> & shell,
+ std::list<MTet4*> & cavity,
+ MVertex *v ){
+ std::list<faceXtet> wrong;
+ for (std::list<faceXtet>::iterator it = shell.begin(); it != shell.end() ;++it) {
+ faceXtet &fxt = *it;
+ bool starShaped = fxt.visible(v);
+ if (!starShaped){
+ wrong.push_back(fxt);
+ }
+ }
+ if (wrong.empty()) return 0;
+ // printf ("cavity %p (shell size %d cavity size %d)is not star shaped (%d faces not visible), correcting it\n",
+ // v,shell.size(),cavity.size(),wrong.size());
+
+ // bool doneNothing = true;
+ while (!wrong.empty()){
+ faceXtet &fxt = *(wrong.begin());
+ std::list<faceXtet>::iterator its = std::find(shell.begin(),shell.end(),fxt);
+ if (its != shell.end()){
+ if (fxt.t1->getNeigh(fxt.i1) && fxt.t1->getNeigh(fxt.i1)->onWhat() == fxt.t1->onWhat() && verifyShell(v,fxt.t1->getNeigh(fxt.i1),shell)){
+ extendCavity (shell,cavity,fxt);
+ // doneNothing = false;
+ }
+ else if (verifyShell(v,fxt.t1,shell)){
+ removeFromCavity (shell,cavity,fxt);
+ // doneNothing = false;
+ }
+ else {
+ return -1;
+ }
+ }
+ wrong.erase(wrong.begin());
+ }
+ // printf("after : shell size %d cavity size %d\n",shell.size(),cavity.size());
+ return 1;
+}
+
void recurFindCavity(std::list<faceXtet> & shell,
std::list<MTet4*> & cavity,
MVertex *v ,
@@ -184,11 +298,13 @@ void recurFindCavity(std::list<faceXtet> & shell,
if (!neigh)
shell.push_back(faceXtet(t, i));
else if (!neigh->isDeleted()){
+ faceXtet fxt (t, i);
int circ = neigh->inCircumSphere(v);
if (circ && (neigh->onWhat() == t->onWhat()))
recurFindCavity(shell, cavity, v, neigh);
- else
- shell.push_back(faceXtet(t, i));
+ else{
+ shell.push_back(fxt);
+ }
}
}
// printf("cavity size %d\n",cavity.size());
@@ -234,6 +350,24 @@ void nonrecurFindCavity(std::list<faceXtet> & shell,
// printf("cavity size %d\n",cavity.size());
}
+void printTets (const char *fn, std::list<MTet4*> &cavity, bool force = false )
+{
+ FILE *f = fopen (fn,"w");
+ fprintf(f,"View \"\"{\n");
+ std::list<MTet4*>::iterator ittet = cavity.begin();
+ std::list<MTet4*>::iterator ittete = cavity.end();
+ while (ittet != ittete){
+ MTet4 *tet = *ittet;
+ if (force || !tet->isDeleted()){
+ MTetrahedron *t = tet->tet();
+ t->writePOS (f, false,false,true,false,false,false);
+ }
+ ittet++;
+ }
+ fprintf(f,"};\n");
+ fclose(f);
+}
+
bool insertVertexB(std::list<faceXtet> &shell,
std::list<MTet4*> &cavity,
@@ -247,13 +381,13 @@ bool insertVertexB(std::list<faceXtet> &shell,
{
std::list<MTet4*> new_cavity;
// check that volume is conserved
- double newVolume = 0;
- double oldVolume = 0;
+ double newVolume = 0;
+ double oldVolume = 0;
std::list<MTet4*>::iterator ittet = cavity.begin();
std::list<MTet4*>::iterator ittete = cavity.end();
- while (ittet != ittete){
- oldVolume += fabs((*ittet)->getVolume());
+ while (ittet != ittete){
+ oldVolume += fabs((*ittet)->getVolume());
++ittet;
}
@@ -264,7 +398,7 @@ bool insertVertexB(std::list<faceXtet> &shell,
bool onePointIsTooClose = false;
while (it != shell.end()){
- MTetrahedron *tr = new MTetrahedron(it->v[0], it->v[1], it->v[2], v);
+ MTetrahedron *tr = new MTetrahedron(it->getVertex(0), it->getVertex(1), it->getVertex(2), v);
double lc = .25 * (vSizes[tr->getVertex(0)->getIndex()] +
vSizes[tr->getVertex(1)->getIndex()] +
@@ -301,14 +435,15 @@ bool insertVertexB(std::list<faceXtet> &shell,
if (otherSide)
new_cavity.push_back(otherSide);
// if (!it->t1->isDeleted())throw;
- newVolume += fabs(t4->getVolume());
+ newVolume += fabs(t4->getVolume());
++it;
}
// OK, the cavity is star shaped
// if (onePointIsTooClose)printf("One point is too close\n");
- // if (fabs(oldVolume - newVolume) > 1.e-10 * oldVolume)printf("Volume do not match %22.15E %22.15E %22.15E\n",oldVolume,newVolume,fabs(oldVolume-newVolume)/newVolume);
- if (fabs(oldVolume - newVolume) < 1.e-10 * oldVolume &&
- !onePointIsTooClose){
+ //if (fabs(oldVolume - newVolume) > 1.e-10 * oldVolume)printf("Volume do not match %22.15E %22.15E %22.15E\n",oldVolume,newVolume,fabs(oldVolume-newVolume)/newVolume);
+ // if (!onePointIsTooClose){
+ if (fabs(oldVolume - newVolume) < 1.e-10 * oldVolume &&
+ !onePointIsTooClose){
connectTets_vector(new_cavity.begin(), new_cavity.end());
allTets.insert(newTets, newTets + shell.size());
@@ -328,6 +463,11 @@ bool insertVertexB(std::list<faceXtet> &shell,
}
else { // The cavity is NOT star shaped
// printf("hola %d %g %g %22.15E\n",onePointIsTooClose, oldVolume, newVolume, 100.*fabs(oldVolume - newVolume) / oldVolume);
+ // new_cavity.clear();
+ // for (unsigned int i = 0; i <shell.size(); i++) new_cavity.push_back(newTets[i]);
+ // printTets ("oldCavity.pos",cavity,true);
+ // printTets ("newCavity.pos",new_cavity);
+ // Msg::Fatal("");
for (unsigned int i = 0; i <shell.size(); i++) myFactory.Free(newTets[i]);
delete [] newTets;
ittet = cavity.begin();
@@ -941,6 +1081,19 @@ double tetcircumcenter(double a[3], double b[3], double c[3], double d[3],
return xxx;
}
+static void memoryCleanup(MTet4Factory &myFactory, std::set<MTet4*, compareTet4Ptr> &allTets){
+ int n1 = allTets.size();
+ std::set<MTet4*,compareTet4Ptr>::iterator itd = allTets.begin();
+ while(itd != allTets.end()){
+ if((*itd)->isDeleted()){
+ myFactory.Free((*itd));
+ allTets.erase(itd++);
+ }
+ else
+ itd++;
+ }
+ // Msg::Info("cleaning up the memory %d -> %d", n1, allTets.size());
+}
void insertVerticesInRegion (GRegion *gr)
{
@@ -963,11 +1116,16 @@ void insertVerticesInRegion (GRegion *gr)
it->first->setIndex(NUM++);
vSizes.push_back(it->second);
vSizesBGM.push_back(it->second);
+ // it->first->x() += (double)rand()/RAND_MAX * 1.e-8;
+ // it->first->y() += (double)rand()/RAND_MAX * 1.e-8;
+ // it->first->z() += (double)rand()/RAND_MAX * 1.e-8;
}
}
- for(unsigned int i = 0; i < gr->tetrahedra.size(); i++)
+ for(unsigned int i = 0; i < gr->tetrahedra.size(); i++){
+ gr->tetrahedra[i]->setVolumePositive();
allTets.insert(myFactory.Create(gr->tetrahedra[i], vSizes,vSizesBGM));
+ }
gr->tetrahedra.clear();
connectTets(allTets.begin(), allTets.end());
@@ -1013,9 +1171,11 @@ void insertVerticesInRegion (GRegion *gr)
// here the classification should be done
int ITER = 0;
-
+ int NB_CORRECTION_OF_CAVITY = 0;
int COUNT_MISS_1 = 0;
int COUNT_MISS_2 = 0;
+
+ clock_t t1 = clock();
while(1){
if(allTets.empty()){
Msg::Error("No tetrahedra in region %d %d", gr->tag(), allTets.size());
@@ -1029,8 +1189,8 @@ void insertVerticesInRegion (GRegion *gr)
allTets.erase(allTets.begin());
}
else{
- if(ITER++ %5000 == 0)
- Msg::Info("%d points created -- Worst tet radius is %g (MISSES %d %d)",
+ if(ITER++ % 5000 == 0)
+ Msg::Info("%9d points created -- Worst tet radius is %8.3f (PTS removed %4d %4d)",
vSizes.size(), worst->getRadius(), COUNT_MISS_1,COUNT_MISS_2);
if(worst->getRadius() < 1) break;
double center[3];
@@ -1052,12 +1212,12 @@ void insertVerticesInRegion (GRegion *gr)
tetcircumcenter(pa,pb,pc,pd, center,&uvw[0],&uvw[1],&uvw[2] );
-
//// A TEST !!!
std::list<faceXtet> shell;
std::list<MTet4*> cavity;
MVertex vv (center[0], center[1], center[2], worst->onWhat());
recurFindCavity(shell, cavity, &vv, worst);
+
bool FOUND = false;
for (std::list<MTet4*>::iterator itc = cavity.begin(); itc != cavity.end(); ++itc){
MTetrahedron *toto = (*itc)->tet();
@@ -1071,11 +1231,23 @@ void insertVerticesInRegion (GRegion *gr)
}
/// END TETS
- // worst->tet()->xyz2uvw(center,uvw);
-
if(FOUND){
MVertex *v = new MVertex(center[0], center[1], center[2], worst->onWhat());
v->setIndex(NUM++);
+
+ // printTets ("before.pos", cavity, true);
+ bool starShaped = true;
+ bool correctCavity = false;
+ while (1){
+ int k = makeCavityStarShaped (shell, cavity, v);
+ if (k == -1){starShaped = false ; break;}
+ else if (k == 0) break;
+ else if (k == 1) correctCavity = true;
+ }
+ if (correctCavity && starShaped) NB_CORRECTION_OF_CAVITY ++;
+ // printTets ("after.pos", cavity, true);
+ //}
+
double lc1 =
(1 - uvw[0] - uvw[1] - uvw[2]) * vSizes[worst->tet()->getVertex(0)->getIndex()] +
uvw[0] * vSizes[worst->tet()->getVertex(1)->getIndex()] +
@@ -1086,8 +1258,9 @@ void insertVerticesInRegion (GRegion *gr)
vSizes.push_back(lc1);
vSizesBGM.push_back(lc);
// compute mesh spacing there
- if(!insertVertexB(shell,cavity,v, worst, myFactory, allTets, vSizes,vSizesBGM)){
+ if(!starShaped || !insertVertexB(shell,cavity,v, worst, myFactory, allTets, vSizes,vSizesBGM)){
COUNT_MISS_1++;
+ // printf("coucou 1 %d\n",ITER);
myFactory.changeTetRadius(allTets.begin(), 0.);
delete v;
}
@@ -1095,6 +1268,7 @@ void insertVerticesInRegion (GRegion *gr)
v->onWhat()->mesh_vertices.push_back(v);
}
else{
+ // printf("coucou 2 %d\n",ITER);
// printf("point outside %12.5E %12.5E %12.5E %12.5E %12.5E\n",VV,uvw[0], uvw[1], uvw[2],1-uvw[0]-uvw[1]-uvw[2]);
myFactory.changeTetRadius(allTets.begin(), 0.0);
COUNT_MISS_2++;
@@ -1105,21 +1279,18 @@ void insertVerticesInRegion (GRegion *gr)
// Normally, a tet mesh contains about 6 times more tets than
// vertices. This allows to clean up the set of tets when lots of
// deleted ones are present in the mesh
- if(allTets.size() > 7 * vSizes.size()){
- int n1 = allTets.size();
- std::set<MTet4*,compareTet4Ptr>::iterator itd = allTets.begin();
- while(itd != allTets.end()){
- if((*itd)->isDeleted()){
- myFactory.Free((*itd));
- allTets.erase(itd++);
- }
- else
- itd++;
- }
- Msg::Info("cleaning up the memory %d -> %d", n1, allTets.size());
+ if(allTets.size() > 7 * vSizes.size() && ITER > 20000){
+ memoryCleanup(myFactory, allTets);
}
}
+ memoryCleanup(myFactory, allTets);
+ clock_t t2 = clock();
+ double dt = (double)(t2-t1)/CLOCKS_PER_SEC;
+ int COUNT_MISS = COUNT_MISS_1+COUNT_MISS_2;
+ Msg::Info("3D Point Insertion Terminated : %9d Delaunay cavities modified for star shapeness",NB_CORRECTION_OF_CAVITY);
+ Msg::Info(" : %9d points could not be inserted among %d",COUNT_MISS,vSizes.size() - COUNT_MISS);
+ Msg::Info(" : %9d tetrahedra created in %8.2f sec. (%d tets/sec.)",allTets.size(),dt,(int)(allTets.size()/dt));
// relocate vertices
int nbReloc = 0;
diff --git a/benchmarks/2d/Square-01.geo b/benchmarks/2d/Square-01.geo
index 6053b54017..5078de8447 100644
--- a/benchmarks/2d/Square-01.geo
+++ b/benchmarks/2d/Square-01.geo
@@ -10,4 +10,4 @@ Line(3) = {1,4};
Line(4) = {4,3};
Line Loop(5) = {1,2,3,4};
Plane Surface(6) = {5};
-//Recombine Surface {6};
+Recombine Surface {6};
diff --git a/benchmarks/2d/conge.geo b/benchmarks/2d/conge.geo
index 39841f15cb..bf9543f499 100644
--- a/benchmarks/2d/conge.geo
+++ b/benchmarks/2d/conge.geo
@@ -79,4 +79,4 @@ Physical Surface(26) = {24,22};
Physical Line(27) = {10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 14, 13, 12, 11};
Physical Line(28) = {17, 16, 20, 19, 18, 15};
Recombine Surface {24, 22};
-Mesh.RecombinationAlgorithm=1;
\ No newline at end of file
+Mesh.Algorithm=8;
--
GitLab