diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 2591aa0a3bd0e5276a850033b0a1080646f8019f..b4f3e016545516f1da76dfe4f0724cee4b301c0a 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -355,6 +355,8 @@ bool GFaceCompound::trivial() const
}
return false;
}
+
+
//For the conformal map the linear system cannot guarantee there is no folding
//of triangles
bool GFaceCompound::checkFolding(std::vector<MVertex*> &ordered) const
@@ -792,7 +794,7 @@ GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
: GFace(m, tag), _compound(compound), _U0(U0), _U1(U1), _V0(V0), _V1(V1), oct(0),
_lsys(lsys),_mapping(mpg), _allowPartition(allowPartition)
{
-
+
if (!_lsys) {
#if defined(HAVE_PETSC) && !defined(HAVE_TAUCS)
_lsys = new linearSystemPETSc<double>;
diff --git a/Geo/GFaceCompound.h b/Geo/GFaceCompound.h
index e1263ac04cc97b202eb6917889253d92d7771ec8..eeab87d002a610459cf2a44f661814b0bbdb2424 100644
--- a/Geo/GFaceCompound.h
+++ b/Geo/GFaceCompound.h
@@ -134,7 +134,7 @@ template<class scalar> class linearSystem;
class GFaceCompound : public GFace {
public:
typedef enum {HARMONIC=1,CONFORMAL=2, CONVEXCOMBINATION=3, MULTISCALE=4} typeOfMapping;
- GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
+ GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
std::list<GEdge*> &U0, std::list<GEdge*> &U1,
std::list<GEdge*> &V0, std::list<GEdge*> &V1,
linearSystem<double>* lsys =0,
@@ -153,6 +153,7 @@ class GFaceCompound : public GFace {
SVector3 *dudu, SVector3 *dvdv, SVector3 *dudv) const{}
virtual SPoint2 getCoordinates(MVertex *v) const { return SPoint2(); }
void parametrize() const {}
+
};
#endif
diff --git a/Geo/GModelIO_Geo.cpp b/Geo/GModelIO_Geo.cpp
index 37873357044aa9641956f81a740afbbaf9337fbf..8c5d0c15a6665ad62f1786363a1eb2cdcee84573 100644
--- a/Geo/GModelIO_Geo.cpp
+++ b/Geo/GModelIO_Geo.cpp
@@ -182,10 +182,15 @@ int GModel::importGEOInternals()
}
int allowPartition = 1;
if (abs(s->TypeOfMapping) != 1) allowPartition = 0;
+
f = new GFaceCompound(this, s->Num, comp,
b[0], b[1], b[2], b[3], 0,
s->TypeOfMapping > 0 ? GFaceCompound::HARMONIC :
GFaceCompound::CONFORMAL, allowPartition);
+ f->meshAttributes.recombine = s->Recombine;
+ f->meshAttributes.recombineAngle = s->RecombineAngle;
+ f->meshAttributes.Method = s->Method;
+ f->meshAttributes.extrude = s->Extrude;
add(f);
}
else if(!f){
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index 854b569d89d1cb21298dbe08fca008298817f685..712018fa3564e056bb204c69febaddc3db33a7cd 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -407,7 +407,7 @@ static void Mesh2D(GModel *m)
// and curve meshes) is global as it depends on a smooth normal
// field generated from the surface mesh of the source surfaces
if(!Mesh2DWithBoundaryLayers(m)){
- //std::for_each(m->firstFace(), m->lastFace(), meshGFace());
+
std::set<GFace*> classFaces;
std::set<GFace*> compFaces;
for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){
@@ -419,6 +419,19 @@ static void Mesh2D(GModel *m)
std::for_each(classFaces.begin(), classFaces.end(), meshGFace());
std::for_each(compFaces.begin(), compFaces.end(), meshGFace());
+ //lloyd optimization
+ if (CTX::instance()->mesh.optimize > 0 ){
+ for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it){
+ GFace *gf = *it;
+ int recombine = gf->meshAttributes.recombine;
+ Msg::Info("LLoyd optimisation for face %d", gf->tag());
+ gf->lloyd(40,recombine);
+
+ if(recombine) recombineIntoQuads(gf);
+
+ }
+ }
+
int nIter = 0;
while(1){
meshGFace mesher;
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 7c2c4959075336a198ee124e9dd2f2df6fe568fe..2bbc4c2c33878fc3718972f8d9cbab2305392053 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -669,18 +669,25 @@ static bool meshGenerator(GFace *gf, int RECUR_ITER,
outputScalarField(m->triangles, name,1);
}
if(CTX::instance()->mesh.remove4triangles)
- removeFourTrianglesNodes(gf,false);
+ removeFourTrianglesNodes(gf,false);
// delete the mesh
delete m;
- if(gf->meshAttributes.recombine)
- recombineIntoQuads(gf);
+ printf("to recombine gf->meshAttributes.recombine=%d \n", gf->meshAttributes.recombine);
+ if(gf->meshAttributes.recombine && CTX::instance()->mesh.optimize == 0){
+ printf("in recombine \n");
+ recombineIntoQuads(gf);
+ }
+
computeElementShapes(gf, gf->meshStatistics.worst_element_shape,
gf->meshStatistics.average_element_shape,
gf->meshStatistics.best_element_shape,
gf->meshStatistics.nbTriangle,
gf->meshStatistics.nbGoodQuality);
+
+
+
return true;
}
@@ -1218,9 +1225,9 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
// delete the mesh
delete m;
- if(gf->meshAttributes.recombine)
- recombineIntoQuads(gf);
-
+ if(gf->meshAttributes.recombine && CTX::instance()->mesh.optimize == 0)
+ recombineIntoQuads(gf);
+
computeElementShapes(gf, gf->meshStatistics.worst_element_shape,
gf->meshStatistics.average_element_shape,
gf->meshStatistics.best_element_shape,
@@ -1405,7 +1412,7 @@ void partitionAndRemesh(GFaceCompound *gf)
int num_gfc = numf + NF + i ;
f_compound.push_back(pf);
Msg::Info("Parametrize Compound Surface (%d) = %d discrete face", num_gfc, pf->tag() );
- GFaceCompound *gfc = new GFaceCompound(gf->model(), num_gfc, f_compound,
+ GFaceCompound *gfc = new GFaceCompound(gf->model(), num_gfc, f_compound,
b[0], b[1], b[2], b[3], 0, gf->getTypeOfMapping() );
gf->model()->add(gfc);
diff --git a/Mesh/meshGFaceLloyd.cpp b/Mesh/meshGFaceLloyd.cpp
index 48fd0c8e1958675fb90379d90825667ab1786b0f..a22b7de4ea40c44ea7ca079a5668f66a08bb00f0 100644
--- a/Mesh/meshGFaceLloyd.cpp
+++ b/Mesh/meshGFaceLloyd.cpp
@@ -58,8 +58,8 @@ void lloydAlgorithm::operator () ( GFace * gf) {
// compute the Voronoi diagram
triangulator.Voronoi();
- // printf("hullSize = %d\n",triangulator.hullSize());
- // triangulator.makePosView("LloydInit.pos");
+ //printf("hullSize = %d\n",triangulator.hullSize());
+ triangulator.makePosView("LloydInit.pos");
// now do the Lloyd iterations
int ITER = 0;
@@ -104,8 +104,9 @@ void lloydAlgorithm::operator () ( GFace * gf) {
}
}
- Msg::Debug("GFace %d Lloyd iteration %d energy %g",gf->tag(),ITER++,ENERGY);
+ Msg::Info("GFace %d Lloyd iteration %d energy %g",gf->tag(),ITER++,ENERGY);
if (ITER > ITER_MAX)break;
+
// compute the Voronoi diagram
triangulator.Voronoi();
}
@@ -133,7 +134,7 @@ void lloydAlgorithm::operator () ( GFace * gf) {
// remesh the face with all those vertices in
meshGFace mesher;
mesher(gf);
- // assign thos vertices to the face internal vertices
+ // assign those vertices to the face internal vertices
gf->mesh_vertices.insert(gf->mesh_vertices.begin(),
gf->_additional_vertices.begin(),
gf->_additional_vertices.end());
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index 319d09c79817060af4d2f34f5e8e7067de604b2f..0d6279373241e387e66c0a162aa0af26494b8042 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -1011,6 +1011,7 @@ struct RecombineTriangle
static void _recombineIntoQuads(GFace *gf)
{
+
std::set<MVertex*> emb_edgeverts;
{
std::list<GEdge*> emb_edges = gf->embeddedEdges();
@@ -1091,6 +1092,7 @@ static void _recombineIntoQuads(GFace *gf)
void recombineIntoQuads(GFace *gf)
{
+
_recombineIntoQuads(gf);
laplaceSmoothing(gf);
_recombineIntoQuads(gf);
diff --git a/Parser/Gmsh.tab.cpp b/Parser/Gmsh.tab.cpp
index 7b12c48ebe88fd092b92a7a443c86cb62f9e7047..602857aa04d5d331ec437d58d2db8188b4414d6b 100644
--- a/Parser/Gmsh.tab.cpp
+++ b/Parser/Gmsh.tab.cpp
@@ -5405,7 +5405,7 @@ yyreduce:
yymsg(0, "Surface %d already exists", num);
}
else{
- Surface *s = Create_Surface(num, MSH_SURF_COMPOUND);
+ Surface *s = Create_Surface(num, MSH_SURF_COMPOUND);
for(int i = 0; i < List_Nbr((yyvsp[(7) - (9)].l)); i++){
s->compound.push_back((int)*(double*)List_Pointer((yyvsp[(7) - (9)].l), i));
s->TypeOfMapping = (yyvsp[(8) - (9)].i);
diff --git a/Parser/Gmsh.y b/Parser/Gmsh.y
index 943e3042cc35c01f13589960668b5fe7304a2323..36a560ddb4a93e603dd872d7bdd20bb3d164bbef 100644
--- a/Parser/Gmsh.y
+++ b/Parser/Gmsh.y
@@ -1462,7 +1462,7 @@ Shape :
yymsg(0, "Surface %d already exists", num);
}
else{
- Surface *s = Create_Surface(num, MSH_SURF_COMPOUND);
+ Surface *s = Create_Surface(num, MSH_SURF_COMPOUND);
for(int i = 0; i < List_Nbr($7); i++){
s->compound.push_back((int)*(double*)List_Pointer($7, i));
s->TypeOfMapping = $8;