diff --git a/Geo/boundaryLayersData.cpp b/Geo/boundaryLayersData.cpp
index af3b5c5f7e5b2c3e8dd6aec1c25b0e60389fd320..f579e417d216043e93c82d969f7c723fdf81552c 100644
--- a/Geo/boundaryLayersData.cpp
+++ b/Geo/boundaryLayersData.cpp
@@ -32,7 +32,7 @@ edgeColumn BoundaryLayerColumns::getColumns(MVertex *v1, MVertex *v2 , int side)
#include "Field.h"
-const int FANSIZE__ = 4;
+const int FANSIZE__ = 5;
/*
@@ -80,26 +80,26 @@ edgeColumn BoundaryLayerColumns::getColumns(MVertex *v1, MVertex *v2 , int side)
if (nbSides == 1){
if (it1 != _fans.end() && it2 == _fans.end() ){
if (aaa(it1->second._e1,e))
- return edgeColumn(getColumn (v1,0),getColumn(v2,0));
+ return edgeColumn(getColumn (v1,0),getColumn(v2,0));
else
- return edgeColumn(getColumn (v1,N1-1),getColumn(v2,0));
+ return edgeColumn(getColumn (v1,N1-1),getColumn(v2,0));
}
if (it2 != _fans.end() && it1 == _fans.end() ){
if (aaa(it2->second._e1,e))
- return edgeColumn(getColumn (v1,0),getColumn(v2,0));
+ return edgeColumn(getColumn (v1,0),getColumn(v2,0));
else
- return edgeColumn(getColumn (v1,0),getColumn(v2,N2-1));
+ return edgeColumn(getColumn (v1,0),getColumn(v2,N2-1));
}
if (it2 != _fans.end() && it1 != _fans.end() ){
int c1, c2;
if (aaa(it1->second._e1,e))
- c1 = 0;
+ c1 = 0;
else
- c1 = N1-1;
+ c1 = N1-1;
if (aaa(it2->second._e1,e))
- c2 = 0;
+ c2 = 0;
else
- c2 = N2-1;
+ c2 = N2-1;
return edgeColumn(getColumn (v1,c1),getColumn(v2,c2));
}
// fan on the right
@@ -169,56 +169,56 @@ static void treat2Connections(GFace *gf, MVertex *_myVert, MEdge &e1, MEdge &e2,
{
std::vector<SVector3> N1,N2;
for (std::multimap<MEdge,SVector3,Less_Edge>::iterator itm =
- _columns->_normals.lower_bound(e1);
+ _columns->_normals.lower_bound(e1);
itm != _columns->_normals.upper_bound(e1); ++itm) N1.push_back(itm->second);
for (std::multimap<MEdge,SVector3,Less_Edge>::iterator itm =
- _columns->_normals.lower_bound(e2);
+ _columns->_normals.lower_bound(e2);
itm != _columns->_normals.upper_bound(e2); ++itm) N2.push_back(itm->second);
if (N1.size() == N2.size()){
for (unsigned int SIDE = 0; SIDE < N1.size() ; SIDE++){
if (!fan){
- SVector3 x = N1[SIDE]*1.01+N2[SIDE];
- x.normalize();
- _dirs.push_back(x);
+ SVector3 x = N1[SIDE]*1.01+N2[SIDE];
+ x.normalize();
+ _dirs.push_back(x);
}
else if (fan){
- // printf("fan \n");
-
- int fanSize = FANSIZE__;
- // if the angle is greater than PI, than reverse the sense
- double alpha1 = atan2(N1[SIDE].y(),N1[SIDE].x());
- double alpha2 = atan2(N2[SIDE].y(),N2[SIDE].x());
- double AMAX = std::max(alpha1,alpha2);
- double AMIN = std::min(alpha1,alpha2);
- MEdge ee[2];
- if (alpha1 > alpha2){
- ee[0] = e2;ee[1] = e1;
- }
- else {
- ee[0] = e1;ee[1] = e2;
- }
- if ( AMAX - AMIN >= M_PI){
- double temp = AMAX;
- AMAX = AMIN + 2*M_PI;
- AMIN = temp;
- MEdge eee0 = ee[0];
- ee[0] = ee[1];ee[1] = eee0;
- }
- _columns->addFan (_myVert,ee[0],ee[1],true);
- for (int i=-1; i<=fanSize; i++){
- double t = (double)(i+1)/ (fanSize+1);
- double alpha = t * AMAX + (1.-t)* AMIN;
- SVector3 x (cos(alpha),sin(alpha),0);
- x.normalize();
- _dirs.push_back(x);
- }
+ // printf("fan \n");
+
+ int fanSize = FANSIZE__;
+ // if the angle is greater than PI, than reverse the sense
+ double alpha1 = atan2(N1[SIDE].y(),N1[SIDE].x());
+ double alpha2 = atan2(N2[SIDE].y(),N2[SIDE].x());
+ double AMAX = std::max(alpha1,alpha2);
+ double AMIN = std::min(alpha1,alpha2);
+ MEdge ee[2];
+ if (alpha1 > alpha2){
+ ee[0] = e2;ee[1] = e1;
+ }
+ else {
+ ee[0] = e1;ee[1] = e2;
+ }
+ if ( AMAX - AMIN >= M_PI){
+ double temp = AMAX;
+ AMAX = AMIN + 2*M_PI;
+ AMIN = temp;
+ MEdge eee0 = ee[0];
+ ee[0] = ee[1];ee[1] = eee0;
+ }
+ _columns->addFan (_myVert,ee[0],ee[1],true);
+ for (int i=-1; i<=fanSize; i++){
+ double t = (double)(i+1)/ (fanSize+1);
+ double alpha = t * AMAX + (1.-t)* AMIN;
+ SVector3 x (cos(alpha),sin(alpha),0);
+ x.normalize();
+ _dirs.push_back(x);
+ }
}
/*
else {
- _dirs.push_back(N1[SIDE]);
- _dirs.push_back(N2[SIDE]);
- }
+ _dirs.push_back(N1[SIDE]);
+ _dirs.push_back(N2[SIDE]);
+ }
*/
}
}
@@ -231,13 +231,13 @@ static void treat3Connections(GFace *gf, MVertex *_myVert, MEdge &e1,
{
std::vector<SVector3> N1,N2,N3;
for (std::multimap<MEdge,SVector3,Less_Edge>::iterator itm =
- _columns->_normals.lower_bound(e1);
+ _columns->_normals.lower_bound(e1);
itm != _columns->_normals.upper_bound(e1); ++itm) N1.push_back(itm->second);
for (std::multimap<MEdge,SVector3,Less_Edge>::iterator itm =
- _columns->_normals.lower_bound(e2);
+ _columns->_normals.lower_bound(e2);
itm != _columns->_normals.upper_bound(e2); ++itm) N2.push_back(itm->second);
for (std::multimap<MEdge,SVector3,Less_Edge>::iterator itm =
- _columns->_normals.lower_bound(e3);
+ _columns->_normals.lower_bound(e3);
itm != _columns->_normals.upper_bound(e3); ++itm) N3.push_back(itm->second);
SVector3 x1,x2;
@@ -319,13 +319,13 @@ static void getEdgesData(GFace *gf,
// check if this edge generates a boundary layer
if (isEdgeOfFaceBL (gf,*ite,blf)){
for(unsigned int i = 0; i< (*ite)->lines.size(); i++){
- MVertex *v1 = (*ite)->lines[i]->getVertex(0);
- MVertex *v2 = (*ite)->lines[i]->getVertex(1);
- allEdges.insert(MEdge(v1,v2));
- _columns->_non_manifold_edges.insert(std::make_pair(v1,v2));
- _columns->_non_manifold_edges.insert(std::make_pair(v2,v1));
- _vertices.insert(v1);
- _vertices.insert(v2);
+ MVertex *v1 = (*ite)->lines[i]->getVertex(0);
+ MVertex *v2 = (*ite)->lines[i]->getVertex(1);
+ allEdges.insert(MEdge(v1,v2));
+ _columns->_non_manifold_edges.insert(std::make_pair(v1,v2));
+ _columns->_non_manifold_edges.insert(std::make_pair(v2,v1));
+ _vertices.insert(v1);
+ _vertices.insert(v2);
}
}
else {
@@ -411,8 +411,15 @@ void getLocalInfoAtNode (MVertex *v, BoundaryLayerField *blf, double &hwall) {
double t;
v->getParameter(0,t);
double hwall_beg = blf->hwall (ge->getBeginVertex()->tag());
- double hwall_end = blf->hwall (ge->getEndVertex()->tag());
- hwall = hwall_beg + (t-t_begin)/(t_end-t_begin) * (hwall_end - hwall_beg);
+ double hwall_end = blf->hwall (ge->getEndVertex()->tag());
+ double x = (t-t_begin)/(t_end-t_begin);
+ double hwallLin = hwall_beg + x * (hwall_end - hwall_beg);
+ double hwall_mid = std::min(hwall_beg, hwall_end);
+ double hwallQuad = hwall_beg * (1-x)*(1-x) +
+ hwall_mid * 2*x*(1-x) +
+ hwall_end * x*x;
+ // we prefer a quadratic growing:
+ hwall = 0*hwallLin + 1*hwallQuad;
}
}
@@ -476,8 +483,8 @@ bool buildAdditionalPoints2D(GFace *gf)
MEdge e1 (*it,_connections[0]);
std::vector<SVector3> N1;
for (std::multimap<MEdge,SVector3,Less_Edge>::iterator itm =
- _columns->_normals.lower_bound(e1);
- itm != _columns->_normals.upper_bound(e1); ++itm) N1.push_back(itm->second);
+ _columns->_normals.lower_bound(e1);
+ itm != _columns->_normals.upper_bound(e1); ++itm) N1.push_back(itm->second);
// one point has only one side and one normal : it has to be at the end of the BL
// then, we have the tangent to the connecting edge
@@ -486,57 +493,58 @@ bool buildAdditionalPoints2D(GFace *gf)
// NO BL BL
if (N1.size() == 1){
- std::vector<MVertex*> Ts;
- for (std::multimap<MVertex*,MVertex*>::iterator itm =
- tangents.lower_bound(*it);
- itm != tangents.upper_bound(*it); ++itm) Ts.push_back(itm->second);
- // end of the BL --> let's add a column that correspond to the
- // model edge that lies after the end of teh BL
- if (Ts.size() == 1){
- // printf("HERE WE ARE IN FACE %d %d\n",gf->tag(),Ts.size());
- // printf("Classif dim %d %d\n",(*it)->onWhat()->dim(),Ts[0]->onWhat()->dim());
- GEdge *ge = dynamic_cast<GEdge*>(Ts[0]->onWhat());
- GVertex *gv = dynamic_cast<GVertex*>((*it)->onWhat());
- if (ge && gv){
- addColumnAtTheEndOfTheBL (ge,gv,_columns,blf);
- }
- }
- else {
- Msg::Error("Impossible BL Configuration -- One Edge -- Tscp.size() = %d",Ts.size());
- }
+ std::vector<MVertex*> Ts;
+ for (std::multimap<MVertex*,MVertex*>::iterator itm =
+ tangents.lower_bound(*it);
+ itm != tangents.upper_bound(*it); ++itm) Ts.push_back(itm->second);
+ // end of the BL --> let's add a column that correspond to the
+ // model edge that lies after the end of teh BL
+ if (Ts.size() == 1){
+ // printf("HERE WE ARE IN FACE %d %d\n",gf->tag(),Ts.size());
+ // printf("Classif dim %d %d\n",(*it)->onWhat()->dim(),Ts[0]->onWhat()->dim());
+ GEdge *ge = dynamic_cast<GEdge*>(Ts[0]->onWhat());
+ GVertex *gv = dynamic_cast<GVertex*>((*it)->onWhat());
+ if (ge && gv){
+ addColumnAtTheEndOfTheBL (ge,gv,_columns,blf);
+ }
+ }
+ else {
+ Msg::Error("Impossible BL Configuration -- One Edge -- Tscp.size() = %d",Ts.size());
+ }
}
else if (N1.size() == 2){
// printf("%g %g --> %g %g \n",e1.getVertex(0)->x(),e1.getVertex(0)->y(),
// e1.getVertex(1)->x(),e1.getVertex(1)->y());
- // printf("N1.size = %d %g %g %g %g\n",N1.size(),N1[0].x(),N1[0].y(),N1[1].x(),N1[1].y());
- SPoint2 p0,p1;
- reparamMeshEdgeOnFace(*it,_connections[0], gf, p0, p1);
-
- int fanSize = FANSIZE__;
- double alpha1 = atan2(N1[0].y(),N1[0].x());
- double alpha2 = atan2(N1[1].y(),N1[1].x());
- double alpha3 = atan2(p1.y()-p0.y(),p1.x()-p0.x());
- double AMAX = std::max(alpha1,alpha2);
- double AMIN = std::min(alpha1,alpha2);
- if (alpha3 > AMAX){
- AMIN += M_PI;
- AMAX += M_PI;
- }
- if ( AMAX - AMIN >= M_PI){
- double temp = AMAX;
- AMAX = AMIN + 2*M_PI;
- AMIN = temp;
- }
- _columns->addFan (*it,e1,e1,true);
- // printf("%g %g --> %g %g\n",N1[0].x(),N1[0].y(),N1[1].x(),N1[1].y());
- for (int i=-1; i<=fanSize; i++){
- double t = (double)(i+1)/ (fanSize+1);
- double alpha = t * AMAX + (1.-t)* AMIN;
- SVector3 x (cos(alpha),sin(alpha),0);
- // printf("%d %g %g %g\n",i,x.x(),x.y(),alpha);
- x.normalize();
- _dirs.push_back(x);
- }
+ // printf("N1.size = %d %g %g %g %g\n",N1.size(),N1[0].x(),N1[0].y(),N1[1].x(),N1[1].y());
+ SPoint2 p0,p1;
+ reparamMeshEdgeOnFace(_connections[0], *it, gf, p0, p1);
+
+ int fanSize = fan ? FANSIZE__ : 0 ;
+ double alpha1 = atan2(N1[0].y(),N1[0].x());
+ double alpha2 = atan2(N1[1].y(),N1[1].x());
+ double alpha3 = atan2(p1.y()-p0.y(),p1.x()-p0.x());
+ double AMAX = std::max(alpha1,alpha2);
+ double AMIN = std::min(alpha1,alpha2);
+ if (alpha3 > AMAX){
+ double temp = AMAX;
+ AMAX = AMIN + 2*M_PI;
+ AMIN = temp;
+ }
+ if (alpha3 < AMIN) {
+ double temp = AMIN;
+ AMIN = AMAX - 2 * M_PI;
+ AMAX = temp;
+ }
+ if (fan) _columns->addFan (*it,e1,e1,true);
+ // printf("%g %g --> %g %g\n",N1[0].x(),N1[0].y(),N1[1].x(),N1[1].y());
+ for (int i=-1; i<=fanSize; i++){
+ double t = (double)(i+1)/ (fanSize+1);
+ double alpha = t * AMAX + (1.-t)* AMIN;
+ SVector3 x (cos(alpha),sin(alpha),0);
+ // printf("%d %g %g %g\n",i,x.x(),x.y(),alpha);
+ x.normalize();
+ _dirs.push_back(x);
+ }
}
}
@@ -552,33 +560,33 @@ bool buildAdditionalPoints2D(GFace *gf)
// < ------------------------------- > //
/* if (endOfTheBL){
- printf("%g %g %d %d %g\n", (*it)->x(), (*it)->y(), DIR, (int)_dirs.size(),
+ printf("%g %g %d %d %g\n", (*it)->x(), (*it)->y(), DIR, (int)_dirs.size(),
dot(n, dirEndOfBL));
}
*/
if (endOfTheBL && dot(n,dirEndOfBL) > .99){
- // printf( "coucou c'est moi\n");
+ // printf( "coucou c'est moi\n");
}
else {
- MVertex *first = *it;
- double hwall;
- getLocalInfoAtNode (first, blf, hwall);
- std::vector<MVertex*> _column;
- SPoint2 par = gf->parFromPoint(SPoint3(first->x(),first->y(),first->z()));
- double L = hwall;
- while(1){
- // printf("L = %g\n",L);
- if (L > blf->thickness) break;
- SPoint2 pnew (par.x() + L * n.x(),
- par.y() + L * n.y());
- GPoint pp = gf->point(pnew);
- MFaceVertex *_current = new MFaceVertex (pp.x(),pp.y(),pp.z(),gf,pnew.x(),pnew.y());
- _current->bl_data = new MVertexBoundaryLayerData;
- _column.push_back(_current);
- int ith = _column.size() ;
- L+= hwall * pow (blf->ratio, ith);
- }
- _columns->addColumn(n,*it, _column /*,_metrics*/);
+ MVertex *first = *it;
+ double hwall;
+ getLocalInfoAtNode (first, blf, hwall);
+ std::vector<MVertex*> _column;
+ SPoint2 par = gf->parFromPoint(SPoint3(first->x(),first->y(),first->z()));
+ double L = hwall;
+ while(1){
+ // printf("L = %g\n",L);
+ if (L > blf->thickness) break;
+ SPoint2 pnew (par.x() + L * n.x(),
+ par.y() + L * n.y());
+ GPoint pp = gf->point(pnew);
+ MFaceVertex *_current = new MFaceVertex (pp.x(),pp.y(),pp.z(),gf,pnew.x(),pnew.y());
+ _current->bl_data = new MVertexBoundaryLayerData;
+ _column.push_back(_current);
+ int ith = _column.size() ;
+ L+= hwall * pow (blf->ratio, ith);
+ }
+ _columns->addColumn(n,*it, _column /*,_metrics*/);
}
}
}
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index 79978c14feda7bd7c9e28bc50bd80f7e1c013bfd..042c91cfc1fa44ad79ea1df9e732c81bf4d9f851 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -451,7 +451,8 @@ static void createPoints(GVertex *gv, GEdge *ge, BoundaryLayerField *blf,
std::vector<MVertex*>& v, const SVector3 &dir)
{
#if defined(HAVE_ANN)
- double L = blf->hwall_n;
+ const double hwall = blf->hwall(gv->tag());
+ double L = hwall;
double LEdge = distance (ge->getBeginVertex()->mesh_vertices[0],
ge->getEndVertex()->mesh_vertices[0]);
while (1){
@@ -459,7 +460,7 @@ static void createPoints(GVertex *gv, GEdge *ge, BoundaryLayerField *blf,
SPoint3 p (gv->x() + dir.x() * L, gv->y() + dir.y() * L, 0.0);
v.push_back(new MEdgeVertex(p.x(), p.y(), p.z(), ge, ge->parFromPoint(p), blf->hfar));
int ith = v.size() ;
- L += blf->hwall_n * pow (blf->ratio, ith);
+ L += hwall * pow (blf->ratio, ith);
}
#endif
}
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index b5675e43a7ae2a4089821b63ca9be2ffa45121f4..8c913f24eb550b59312b2de7ccd309616c326ea9 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -646,7 +646,10 @@ static void addOrRemove(MVertex *v1, MVertex *v2, std::set<MEdge,Less_Edge> & be
}
}
-static void modifyInitialMeshForTakingIntoAccountBoundaryLayers(GFace *gf)
+static void modifyInitialMeshForTakingIntoAccountBoundaryLayers(GFace *gf,
+ std::vector<MQuadrangle*> &blQuads,
+ std::vector<MTriangle*> &blTris,
+ std::set<MVertex*> &verts)
{
if (!buildAdditionalPoints2D (gf))return;
BoundaryLayerColumns* _columns = gf->getColumns();
@@ -654,15 +657,12 @@ static void modifyInitialMeshForTakingIntoAccountBoundaryLayers(GFace *gf)
std::set<MEdge,Less_Edge> bedges;
std::set<MEdge,Less_Edge> removed;
- std::vector<MQuadrangle*> blQuads;
- std::vector<MTriangle*> blTris;
std::list<GEdge*> edges = gf->edges();
std::list<GEdge*> embedded_edges = gf->embeddedEdges();
edges.insert(edges.begin(), embedded_edges.begin(),embedded_edges.end());
std::list<GEdge*>::iterator ite = edges.begin();
FILE *ff2 = Fopen ("tato.pos","w");
if(ff2) fprintf(ff2,"View \" \"{\n");
- std::set<MVertex*> verts;
std::vector<MLine*> _lines;
@@ -702,6 +702,22 @@ static void modifyInitialMeshForTakingIntoAccountBoundaryLayers(GFace *gf)
v22->x(),v22->y(),v22->z(),
v21->x(),v21->y(),v21->z(),l+1,l+1,l+1,l+1);
}
+ if (c1._column.size() != c2._column.size()) {
+ MVertex *v11,*v12,*v;
+ v11 = c1._column[N-1];
+ v12 = c2._column[N-1];
+ v = c1._column.size() > c2._column.size() ?
+ c1._column[N] : c2._column[N];
+ MTriangle *qq = new MTriangle(v11,v12,v);
+ qq->setPartition (N+1);
+ myCol.push_back(qq);
+ blTris.push_back(qq);
+ if(ff2)
+ fprintf(ff2,"ST (%g,%g,%g,%g,%g,%g,%g,%g,%g){%d,%d,%d};\n",
+ v11->x(),v11->y(),v11->z(),
+ v12->x(),v12->y(),v12->z(),
+ v->x(),v->y(),v->z(),N+1,N+1,N+1);
+ }
// int M = std::max(c1._column.size(),c2._column.size());
for (unsigned int l=0;l<myCol.size();l++)_columns->_toFirst[myCol[l]] = myCol[0];
_columns->_elemColumns[myCol[0]] = myCol;
@@ -808,10 +824,6 @@ static void modifyInitialMeshForTakingIntoAccountBoundaryLayers(GFace *gf)
deMeshGFace kil_;
kil_(gf);
meshGenerator(gf, 0, 0, true , false, &hop);
-
- gf->quadrangles = blQuads;
- gf->triangles.insert(gf->triangles.begin(),blTris.begin(),blTris.end());
- gf->mesh_vertices.insert(gf->mesh_vertices.begin(),verts.begin(),verts.end());
}
static bool inside_domain(MElementOctree* octree,double x,double y)
@@ -1532,6 +1544,10 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
// fill the small gmsh structures
BDS2GMSH(m, gf, recoverMap);
+ std::vector<MQuadrangle*> blQuads;
+ std::vector<MTriangle*> blTris;
+ std::set<MVertex*> verts;
+
bool infty = false;
if (gf->getMeshingAlgo() == ALGO_2D_FRONTAL_QUAD ||
gf->getMeshingAlgo() == ALGO_2D_PACK_PRLGRMS ||
@@ -1541,7 +1557,7 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
if (infty)
buildBackGroundMesh (gf);
// BOUNDARY LAYER
- modifyInitialMeshForTakingIntoAccountBoundaryLayers(gf);
+ modifyInitialMeshForTakingIntoAccountBoundaryLayers(gf, blQuads, blTris, verts);
}
// the delaunay algo is based directly on internal gmsh structures
@@ -1582,9 +1598,18 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
delete m;
+
+
+ gf->quadrangles.insert(gf->quadrangles.begin(),blQuads.begin(),blQuads.end());
+ gf->triangles.insert(gf->triangles.begin(),blTris.begin(),blTris.end());
+ gf->mesh_vertices.insert(gf->mesh_vertices.begin(),verts.begin(),verts.end());
+
+
if((CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine) &&
!CTX::instance()->mesh.optimizeLloyd && !onlyInitialMesh && CTX::instance()->mesh.algoRecombine != 2)
recombineIntoQuads(gf);
+
+
computeElementShapes(gf, gf->meshStatistics.worst_element_shape,
gf->meshStatistics.average_element_shape,
@@ -2334,7 +2359,13 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
if (infty)
buildBackGroundMesh (gf, &equivalence, ¶metricCoordinates);
// BOUNDARY LAYER
- modifyInitialMeshForTakingIntoAccountBoundaryLayers(gf);
+ std::vector<MQuadrangle*> blQuads;
+ std::vector<MTriangle*> blTris;
+ std::set<MVertex*> verts;
+ modifyInitialMeshForTakingIntoAccountBoundaryLayers(gf, blQuads, blTris, verts);
+ gf->quadrangles.insert(gf->quadrangles.begin(),blQuads.begin(),blQuads.end());
+ gf->triangles.insert(gf->triangles.begin(),blTris.begin(),blTris.end());
+ gf->mesh_vertices.insert(gf->mesh_vertices.begin(),verts.begin(),verts.end());
if(algoDelaunay2D(gf)){
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index ac8061047be80da751b3e8a601889a273df308d1..785d7f39440205e6b3a3669fcde739ea1d8940a3 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -119,7 +119,7 @@ void buildMeshGenerationDataStructures(GFace *gf,
std::cout << "***************************************** NULL" << std::endl;
throw;
}
-
+
for(unsigned int i = 0;i < gf->triangles.size(); i++)
setLcs(gf->triangles[i], vSizesMap, data);
@@ -478,7 +478,8 @@ static int _removeTwoQuadsNodes(GFace *gf)
std::set<MElement*> touched;
std::set<MVertex*> vtouched;
while (it != adj.end()) {
- if(it->second.size()==2 && it->first->onWhat()->dim() == 2) {
+ MVertex *v = it->first;
+ if(it->second.size()==2 && v->onWhat()->dim() == 2) {
MElement *q1 = it->second[0];
MElement *q2 = it->second[1];
if (q1->getNumVertices() == 4 &&
@@ -486,7 +487,7 @@ static int _removeTwoQuadsNodes(GFace *gf)
touched.find(q1) == touched.end() && touched.find(q2) == touched.end()){
int comm = 0;
for (int i=0;i<4;i++){
- if (q1->getVertex(i) == it->first){
+ if (q1->getVertex(i) == v){
comm = i;
break;
}
@@ -496,8 +497,8 @@ static int _removeTwoQuadsNodes(GFace *gf)
MVertex *v3 = q1->getVertex((comm+3)%4);
MVertex *v4 = 0;
for (int i=0;i<4;i++){
- if (q2->getVertex(i) != v1 && q2->getVertex(i) != v2 &&
- q2->getVertex(i) != v3 && q2->getVertex(i) != it->first){
+ if (q2->getVertex(i) != v1 && q2->getVertex(i) != v3 &&
+ q2->getVertex(i) != v){
v4 = q2->getVertex(i);
break;
}
@@ -518,13 +519,14 @@ static int _removeTwoQuadsNodes(GFace *gf)
touched.insert(q1);
touched.insert(q2);
gf->quadrangles.push_back(q);
- vtouched.insert(it->first);
+ vtouched.insert(v);
}
}
}
it++;
}
std::vector<MQuadrangle*> quadrangles2;
+ quadrangles2.reserve(gf->quadrangles.size() - touched.size());
for(unsigned int i = 0; i < gf->quadrangles.size(); i++){
if(touched.find(gf->quadrangles[i]) == touched.end()){
quadrangles2.push_back(gf->quadrangles[i]);
@@ -536,6 +538,7 @@ static int _removeTwoQuadsNodes(GFace *gf)
gf->quadrangles = quadrangles2;
std::vector<MVertex*> mesh_vertices2;
+ mesh_vertices2.reserve(gf->mesh_vertices.size() - vtouched.size());
for(unsigned int i = 0; i < gf->mesh_vertices.size(); i++){
if(vtouched.find(gf->mesh_vertices[i]) == vtouched.end()){
mesh_vertices2.push_back(gf->mesh_vertices[i]);
@@ -707,7 +710,7 @@ static int _removeDiamonds(GFace *gf)
v2->onWhat()->dim() == 2 &&
v3->onWhat()->dim() == 2 &&
v4->onWhat()->dim() == 2 &&
- it1->second.size() ==3 && it3->second.size() == 3 &&
+ it1->second.size() == 3 && it3->second.size() == 3 &&
_tryToCollapseThatVertex (gf, it1->second, it3->second,
q, v1, v3)){
touched.insert(v1);
@@ -847,14 +850,11 @@ void getAllBoundaryLayerVertices (GFace *gf, std::set<MVertex*> &vs){
vs.clear();
BoundaryLayerColumns* _columns = gf->getColumns();
if (!_columns)return;
- for ( std::map<MElement*,std::vector<MElement*> >::iterator it = _columns->_elemColumns.begin();
- it != _columns->_elemColumns.end();it++){
- std::vector<MElement *> &e = it->second;
- for (unsigned int i=0;i<e.size();i++){
- for (int j=0;j<e[i]->getNumVertices();j++){
- vs.insert(e[i]->getVertex(j));
- }
- }
+ for ( std::multimap<MVertex*, BoundaryLayerData>::iterator it = _columns->_data.begin();
+ it != _columns->_data.end();it++){
+ BoundaryLayerData &data = it->second;
+ for (size_t i = 0;i<data._column.size();i++)
+ vs.insert(data._column[i]);
}
}
@@ -1187,6 +1187,7 @@ static int _recombineIntoQuads(GFace *gf, double minqual, bool cubicGraph = 1)
std::sort(pairs.begin(),pairs.end());
std::set<MElement*> touched;
+
if(CTX::instance()->mesh.algoRecombine != 0){
#if defined(HAVE_BLOSSOM)
int ncount = gf->triangles.size();
@@ -1233,6 +1234,7 @@ static int _recombineIntoQuads(GFace *gf, double minqual, bool cubicGraph = 1)
}
}
+
double matzeit = 0.0;
char MATCHFILE[256];
sprintf(MATCHFILE,".face.match");
@@ -1294,6 +1296,7 @@ static int _recombineIntoQuads(GFace *gf, double minqual, bool cubicGraph = 1)
#endif
}
+
std::vector<RecombineTriangle>::iterator itp = pairs.begin();
while(itp != pairs.end()){
// recombine if difference between max quad angle and right
@@ -1302,6 +1305,7 @@ static int _recombineIntoQuads(GFace *gf, double minqual, bool cubicGraph = 1)
if(itp->angle < gf->meshAttributes.recombineAngle){
MElement *t1 = itp->t1;
MElement *t2 = itp->t2;
+
if(touched.find(t1) == touched.end() &&
touched.find(t2) == touched.end()){
touched.insert(t1);
@@ -1339,6 +1343,7 @@ static int _recombineIntoQuads(GFace *gf, double minqual, bool cubicGraph = 1)
}
gf->triangles = triangles2;
+
if(CTX::instance()->mesh.algoRecombine != 1){
quadsToTriangles(gf, minqual);
}
@@ -1374,7 +1379,7 @@ void recombineIntoQuads(GFace *gf,
double t1 = Cpu();
bool haveParam = true;
- bool saveAll = false; //CTX::instance()->mesh.saveAll;
+ bool saveAll = CTX::instance()->mesh.saveAll;
if(gf->geomType() == GEntity::DiscreteSurface && !gf->getCompound())
haveParam = false;
@@ -1382,10 +1387,10 @@ void recombineIntoQuads(GFace *gf,
int success = _recombineIntoQuads(gf, minqual);
if (saveAll) gf->model()->writeMSH("raw.msh");
+
if(haveParam && nodeRepositioning){
RelocateVertices (gf,CTX::instance()->mesh.nbSmoothing);
}
-
// blossom-quad algo
if(success && CTX::instance()->mesh.algoRecombine != 0){
if(topologicalOpti){