diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index e5d2d9329cff7b87005f2e5a0aeebd37ff8a94ac..3546ef6515ec281167a1f05a47a284b004e70640 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1423,224 +1423,6 @@ static bool buildConsecutiveListOfVertices(GFace *gf, GEdgeLoop &gel,
return true;
}
-
-static void printParamGrid(GFace *gf, std::vector<MVertex*> vert1, std::vector<MVertex*> vert2,
- std::vector<MVertex*> e01, std::vector<MVertex*> e10,
- std::vector<MVertex*> e23, std::vector<MVertex*> e32,
- std::vector<MVertex*> e02, std::vector<MVertex*> e13,
- std::vector<MQuadrangle*> quads)
-{
-
- std::vector<SPoint2> p1,p2;
- for (unsigned int i = 0; i< vert1.size(); i++){
- SPoint2 pi;
- reparamMeshVertexOnFace(vert1[i], gf, pi);
- p1.push_back(pi);
- }
- for (unsigned int j = 0; j< vert2.size(); j++){
- SPoint2 pj;
- reparamMeshVertexOnFace(vert2[j], gf, pj);
- p2.push_back(pj);
- }
-
-
- char name[234];
- sprintf(name,"paramGrid_%d.pos", gf->tag());
- FILE *f = fopen(name,"w");
- fprintf(f,"View \"%s\" {\n",name);
-
- for (unsigned int i = 0; i < p1.size(); i++)
- fprintf(f,"SP(%g,%g,%g) {%d};\n", p1[i].x(), p1[i].y(), 0.0, i);
-
- for (unsigned int j = 0; j < p2.size(); j++)
- fprintf(f,"SP(%g,%g,%g) {%d};\n", p2[j].x(), p2[j].y(), 0.0, 100+j);
-
- fprintf(f,"};\n");
- fclose(f);
-
- char name2[234];
- sprintf(name2,"paramEdges_%d.pos", gf->tag());
- FILE *f2 = fopen(name2,"w");
- fprintf(f2,"View \"%s\" {\n",name2);
-
- for (unsigned int i = 0; i < e01.size(); i++){
- SPoint2 pi; reparamMeshVertexOnFace(e01[i], gf, pi);
- fprintf(f2,"SP(%g,%g,%g) {%d};\n", pi.x(), pi.y(), 0.0, 1);
- }
-
- for (unsigned int i = 0; i < e10.size(); i++){
- SPoint2 pi; reparamMeshVertexOnFace(e10[i], gf, pi);
- fprintf(f2,"SP(%g,%g,%g) {%d};\n", pi.x(), pi.y(), 0.0, 10);
- }
-
- for (unsigned int i = 0; i < e23.size(); i++){
- SPoint2 pi; reparamMeshVertexOnFace(e23[i], gf, pi);
- fprintf(f2,"SP(%g,%g,%g) {%d};\n", pi.x(), pi.y(), 0.0, 23);
- }
-
- for (unsigned int i = 0; i < e32.size(); i++){
- SPoint2 pi; reparamMeshVertexOnFace(e32[i], gf, pi);
- fprintf(f2,"SP(%g,%g,%g) {%d};\n", pi.x(), pi.y(), 0.0, 32);
- }
-
- for (unsigned int i = 0; i < e02.size(); i++){
- SPoint2 pi; reparamMeshVertexOnFace(e02[i], gf, pi);
- fprintf(f2,"SP(%g,%g,%g) {%d};\n", pi.x(), pi.y(), 0.0, 2);
- }
-
- for (unsigned int i = 0; i < e13.size(); i++){
- SPoint2 pi; reparamMeshVertexOnFace(e13[i], gf, pi);
- fprintf(f2,"SP(%g,%g,%g) {%d};\n", pi.x(), pi.y(), 0.0, 13);
- }
- fprintf(f2,"};\n");
- fclose(f2);
-
-
- char name3[234];
- sprintf(name3,"quad_%d.pos", gf->tag());
- FILE *f3 = fopen(name3,"w");
- fprintf(f3,"View \"%s\" {\n",name3);
- for (unsigned int i = 0; i < quads.size(); i++){
- quads[i]->writePOS(f3,true,false,false,false,false,false);
- }
- fprintf(f3,"};\n");
- fclose(f3);
-
-
- return;
-
-}
-static void createRegularTwoCircleGrid (Centerline *center, GFace *gf)
-{
-#if defined(HAVE_ANN)
- std::list<GEdge*> bedges = gf->edges();
- std::list<GEdge*>::iterator itb = bedges.begin();
- std::list<GEdge*>::iterator ite = bedges.end(); ite--;
- GEdge *ge1 = (*itb)->getCompound();
- GEdge *ge2 = (*ite)->getCompound();
- int N = ge1->mesh_vertices.size() + 1;
- int N2 = ge2->mesh_vertices.size() + 1;
- if (N != N2 || N%2 != 0){
- Msg::Error("You should have a pair number of points specified in centerline field \n");
- }
-
- //vert1 is the outer circle
- //vert2 is the inner circle
- // - - - -
- // - -
- // - - -
- // v0- v2- -v3 -v1
- // - - -
- // - -
- // - - - -
-
- std::vector<MVertex*> vert1, vert2, vert_temp;
- vert1.push_back(ge1->getBeginVertex()->mesh_vertices[0]);
- vert2.push_back(ge2->getBeginVertex()->mesh_vertices[0]);
- for(int i = 0; i < N-1; i++) {
- vert1.push_back(ge1->mesh_vertices[i]);
- vert2.push_back(ge2->mesh_vertices[i]);
- }
- SPoint2 pv1; reparamMeshVertexOnFace(vert1[0], gf, pv1);
- SPoint2 pv2; reparamMeshVertexOnFace(vert2[0], gf, pv2);
- SPoint2 pv1b; reparamMeshVertexOnFace(vert1[1], gf, pv1b);
- SPoint2 pv2b; reparamMeshVertexOnFace(vert2[1], gf, pv2b);
- SPoint2 pv1c; reparamMeshVertexOnFace(vert1[2], gf, pv1c);
- SPoint2 pv2c; reparamMeshVertexOnFace(vert2[2], gf, pv2c);
- SVector3 vec1(pv1b.x()-pv1.x(),pv1b.y()-pv1.y() , 0.0);
- SVector3 vec1b(pv1c.x()-pv1b.x(),pv1c.y()-pv1b.y() ,0.0);
- SVector3 vec2(pv2b.x()-pv2.x(),pv2b.y()-pv2.y() ,0.0);
- SVector3 vec2b(pv2c.x()-pv2b.x(),pv2c.y()-pv2b.y() , 0.0);
- int sign2 = (dot(crossprod(vec1,vec1b),crossprod(vec2,vec2b)) < 0) ? -1 : +1 ;
- double n1 = sqrt(pv1.x()*pv1.x()+pv1.y()*pv1.y());
- double n2 = sqrt(pv2.x()*pv2.x()+pv2.y()*pv2.y());
- if (n2 > n1) {
- vert_temp = vert1;
- vert1.clear();
- vert1 = vert2;
- vert2.clear();
- vert2 = vert_temp;
- }
-
- ANNpointArray nodes = annAllocPts(N, 3);
- ANNidxArray index = new ANNidx[1];
- ANNdistArray dist = new ANNdist[1];
- for (int ind = 0; ind < N; ind++){
- SPoint2 pp2; reparamMeshVertexOnFace(vert2[ind], gf, pp2);
- nodes[ind][0] = pp2.x();
- nodes[ind][1] = pp2.y();
- nodes[ind][2] = 0.0;
- }
- ANNkd_tree *kdtree = new ANNkd_tree(nodes, N, 3);
- SPoint2 pp1; reparamMeshVertexOnFace(vert1[0], gf, pp1);
- double xyz[3] = {pp1.x(), pp1.y(),0.0};
- kdtree->annkSearch(xyz, 1, index, dist);
- int close_ind = index[0];
- double length = sqrt((vert1[0]->x()-vert2[close_ind]->x())*(vert1[0]->x()-vert2[close_ind]->x())+
- (vert1[0]->y()-vert2[close_ind]->y())*(vert1[0]->y()-vert2[close_ind]->y())+
- (vert1[0]->z()-vert2[close_ind]->z())*(vert1[0]->z()-vert2[close_ind]->z()));
-
- double rad = center->operator()(vert1[0]->x(), vert1[0]->y(), vert1[0]->z());
- double lc = 2*M_PI*rad/N;
- int M = length/lc;
- printf("length =%g arc=%g \n", length, 2*M_PI*rad);
-
- delete kdtree;
- delete[]index;
- delete[]dist;
- annDeallocPts(nodes);
-
- MVertex *v0 = vert1[0]; SPoint2 p0; reparamMeshVertexOnFace(v0, gf, p0);
- MVertex *v1 = vert1[N/2]; SPoint2 p1; reparamMeshVertexOnFace(v1, gf, p1);
- MVertex *v2 = vert2[close_ind]; SPoint2 p2; reparamMeshVertexOnFace(v2, gf, p2);
- MVertex *v3 = vert2[(close_ind+N/2)%N]; SPoint2 p3; reparamMeshVertexOnFace(v3, gf, p3);
-
- printf("grid N = %d M = %d\n", N , M);
- printf("v2 =%d v3=%d \n", close_ind, (close_ind+N/2)%N);
- std::vector<MVertex*> e01,e10,e23,e32;//edges without first and last vertex
- for (int i=1;i<N/2-1;i++) e01.push_back(vert1[i]);
- for (int i=N/2+1;i<N-1;i++) e10.push_back(vert1[i]);
- for (int i=1;i<N/2-1;i++) e23.push_back(vert2[(close_ind+i)%N]);
- for (int i=N/2+1;i<N-1;i++) e32.push_back(vert2[(close_ind+i)%N]);
-
- std::vector<MVertex*> e02 = saturateEdge (gf,p0,p2,M);
- std::vector<MVertex*> e13 = saturateEdge (gf,p1,p3,M);
-
- std::vector<MQuadrangle*> q;
- std::vector<MVertex*> e_inner1 = e23;
- std::vector<MVertex*> e_inner2 = e32;
- if (sign2 <0){
- e_inner1 = e32;
- e_inner2 = e23;
- }
-
- createRegularMesh (gf,
- v0, p0,
- e01, +1,
- v1,p1,
- e13, +1,
- v3,p3,
- e_inner1, -sign2,
- v2,p2,
- e02, -1,
- q);
-
- createRegularMesh (gf,
- v0,p0,
- e02, +1,
- v2,p2,
- e_inner2, -sign2,
- v3,p3,
- e13, -1,
- v1,p1,
- e10, +1,
- q);
-
- printParamGrid(gf, vert1, vert2, e01,e10,e23,e32,e02,e13, q);
-#endif;
-}
-
->>>>>>> .r12521
static bool meshGeneratorElliptic(GFace *gf, bool debug = true)
{
#if defined(HAVE_ANN)
diff --git a/Mesh/meshGFaceElliptic.cpp b/Mesh/meshGFaceElliptic.cpp
index 9ac94b28829d8b900cd814dbfcdf315b963feb91..a26a58f797aa1d7e25a3bc12338a937a5716e00d 100644
--- a/Mesh/meshGFaceElliptic.cpp
+++ b/Mesh/meshGFaceElliptic.cpp
@@ -496,8 +496,8 @@ bool createRegularTwoCircleGrid (Centerline *center, GFace *gf)
pe_inner2 = pe23;
}
- std::vector<MQuadrangle*> quads, smoothQuad;
- std::vector<MVertex*> newv, smoothv;
+ std::vector<MQuadrangle*> quads, quadS1, quadS2;
+ std::vector<MVertex*> newv, newv1, newv2;
fullMatrix<SPoint2> uv1, uv2;
createRegularGrid (gf,
v0, p0,
@@ -510,8 +510,9 @@ bool createRegularTwoCircleGrid (Centerline *center, GFace *gf)
e02, pe02, -1,
quads, newv, uv1);
- printQuads(gf, uv1, quads, 0);
- ellipticSmoother(uv1, gf, smoothQuad, smoothv);
+ //updateFaceQuads(gf, quads, newv);
+ //printQuads(gf, uv1, quads, 0);
+ ellipticSmoother(uv1, gf, quadS1, newv1);
createRegularGrid (gf,
v0,p0,
@@ -523,11 +524,14 @@ bool createRegularTwoCircleGrid (Centerline *center, GFace *gf)
v1,p1,
e10, pe10, +1,
quads, newv, uv2);
+
+ ellipticSmoother(uv2, gf, quadS2, newv2);
- //updateFaceQuads(gf, quads, newv);
- //ellipticSmoother
-
+ quads.clear();
+ for (int i= 0; i< quadS1.size(); i++) quads.push_back(quadS1[i]);
+ for (int i= 0; i< quadS2.size(); i++) quads.push_back(quadS2[i]);
printParamGrid(gf, vert1, vert2, e01,e10,e23,e32,e02,e13, quads);
+
return true;
}