diff --git a/Geo/GFace.h b/Geo/GFace.h
index 2daf19571c65f1b128305a1c5c8c65fafcc466bf..67447ca4e9bf611d20950bb3a98c515817d8f83c 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -131,8 +131,17 @@ class GFace : public GEntity
std::vector<MQuadrangle*> quadrangles;
struct {
+ // do we recombine the triangles of the mesh ?
int recombine;
+ // what is the treshold angle for recombination
double recombineAngle;
+ // is this surface meshed using a transfinite interpolation
+ int Method;
+ // these are the 3 corners of the interpolation
+ std::vector<GVertex* > corners;
+ // all diagonals of the triangulation are left (1), right (2)
+ // or alternated (3)
+ int transfiniteArrangement;
} meshAttributes ;
};
diff --git a/Geo/gmshFace.cpp b/Geo/gmshFace.cpp
index 332f886bd60228ef1bb70a7a3e0022bc324f5d11..6622a37dfa343a10c39763693d3764f242cd3f1f 100644
--- a/Geo/gmshFace.cpp
+++ b/Geo/gmshFace.cpp
@@ -1,4 +1,4 @@
-// $Id: gmshFace.cpp,v 1.17 2006-08-29 10:39:48 remacle Exp $
+// $Id: gmshFace.cpp,v 1.18 2006-09-05 21:37:59 remacle Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -71,6 +71,20 @@ gmshFace::gmshFace(GModel *m, Surface *face)
meshAttributes.recombine = s->Recombine;
meshAttributes.recombineAngle = s->RecombineAngle;
+ meshAttributes.Method = s->Method;
+ if (s->Method == TRANSFINI)
+ {
+ meshAttributes.Method = s->Method;
+ meshAttributes.transfiniteArrangement = s->Recombine_Dir;
+ for (int i=0;i< List_Nbr(s->TrsfPoints);i++)
+ {
+ Vertex *corn;
+ List_Read(s->TrsfPoints, i, &corn);
+ GVertex *gv = m->vertexByTag(corn->Num);
+ if(!gv) throw;
+ meshAttributes.corners.push_back(gv);
+ }
+ }
}
gmshFace::gmshFace(GModel *m, int num)
diff --git a/Mesh/Makefile b/Mesh/Makefile
index 76fe131e9839459e579602b360317ed1ee2cd2d5..0e79b693dcc9c3053cf7cc75a34fd0a24270e66a 100644
--- a/Mesh/Makefile
+++ b/Mesh/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.128 2006-09-03 07:44:10 geuzaine Exp $
+# $Id: Makefile,v 1.129 2006-09-05 21:37:59 remacle Exp $
#
# Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
#
@@ -65,6 +65,7 @@ SRC = 1D_Mesh.cpp \
Metric.cpp \
meshGEdge.cpp \
meshGFace.cpp \
+ meshGFaceTransfinite.cpp \
meshGRegion.cpp \
Nurbs.cpp \
Interpolation.cpp \
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index d7481e218ce1bb5c1a893d2ef430904a1262875c..4f204b1f336f0bc50a0fe9c5490fb69c460e02e4 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFace.cpp,v 1.10 2006-09-01 10:10:05 remacle Exp $
+// $Id: meshGFace.cpp,v 1.11 2006-09-05 21:37:59 remacle Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -103,21 +103,18 @@ public :
}
};
-class fromParametricToCartesian
+fromParametricToCartesian::fromParametricToCartesian ( GFace *_gf )
+ : gf(_gf)
{
- GFace *gf;
-public :
- fromParametricToCartesian ( GFace *_gf )
- : gf(_gf){}
- void operator () (MVertex * v)
- {
- GPoint coords = gf->point (SPoint2(v->x(),v->y()));
- v->x() = coords.x();
- v->y() = coords.y();
- v->z() = coords.z();
- }
-};
+}
+void fromParametricToCartesian::operator () (MVertex * v)
+{
+ GPoint coords = gf->point (SPoint2(v->x(),v->y()));
+ v->x() = coords.x();
+ v->y() = coords.y();
+ v->z() = coords.z();
+}
void computeEdgeLoops (const GFace *gf,
std::vector<MVertex*> & all_mvertices,
@@ -362,7 +359,7 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
++it;
}
- Msg(STATUS2," %d triangles : conv %g -> %g (%g sec)",m.triangles.size(),maxL,1.5,(double)(clock()-t1)/CLOCKS_PER_SEC);
+ // Msg(STATUS2," %d triangles : conv %g -> %g (%g sec)",m.triangles.size(),maxL,1.5,(double)(clock()-t1)/CLOCKS_PER_SEC);
if ((minL > 0.4 && maxL < 1.5) || IT > NIT)break;
@@ -1010,6 +1007,9 @@ void meshGFace :: operator() (GFace *gf)
// Only apply this technique to unknown surfaces or planar surfaces
// when it is unknown, try your best ...
+ if (MeshTransfiniteSurface(gf))return;
+
+
std::vector<MVertex*> points;
std::vector<int> indices;
// compute loops on the fly
diff --git a/Mesh/meshGFace.h b/Mesh/meshGFace.h
index a3c4191f045c2d2da3e6b1d6b312f287807e709e..cfd1401e6d79b238159092aa0dcc8ab81e9f3f9d 100644
--- a/Mesh/meshGFace.h
+++ b/Mesh/meshGFace.h
@@ -20,7 +20,25 @@
//
// Please report all bugs and problems to <gmsh@geuz.org>.
+#include <vector>
+class MVertex;
class GFace;
+// compute edge loops of the face, all_mvertices
+// are the vertices of the
+void computeEdgeLoops (const GFace *gf,
+ std::vector<MVertex*> & all_mvertices,
+ std::vector<int> &indices);
+
+int MeshTransfiniteSurface( GFace *gf);
+
+class fromParametricToCartesian
+{
+ GFace *gf;
+public :
+ fromParametricToCartesian ( GFace *_gf ) ;
+ void operator () (MVertex * v);
+};
+
// Create the mesh of the face
class meshGFace
{
diff --git a/Mesh/meshGFaceTransfinite.cpp b/Mesh/meshGFaceTransfinite.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..121effb6849610c39b23cf2ae599b54017042ea7
--- /dev/null
+++ b/Mesh/meshGFaceTransfinite.cpp
@@ -0,0 +1,226 @@
+// $Id: meshGFaceTransfinite.cpp,v 1.1 2006-09-05 21:37:59 remacle Exp $
+//
+// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
+//
+// This program is free software; you can redistribute it and/or modify
+// it under the terms of the GNU General Public License as published by
+// the Free Software Foundation; either version 2 of the License, or
+// (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program; if not, write to the Free Software
+// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
+// USA.
+//
+// Please report all bugs and problems to <gmsh@geuz.org>.
+
+#include <map>
+#include "meshGFace.h"
+#include "GVertex.h"
+#include "GEdge.h"
+#include "GFace.h"
+#include "MVertex.h"
+#include "MElement.h"
+#include "Context.h"
+#include "Message.h"
+
+#define TRAN_QUA(c1,c2,c3,c4,s1,s2,s3,s4,u,v) \
+ (1.-u)*c4+u*c2+(1.-v)*c1+v*c3-((1.-u)*(1.-v)*s1+u*(1.-v)*s2+u*v*s3+(1.-u)*v*s4)
+
+int MeshTransfiniteSurface( GFace *gf)
+{
+ if(gf->meshAttributes.Method != TRANSFINI)
+ return 0;
+ std::vector <MVertex *> corners;
+ std::vector <MVertex *> d_vertices;
+ std::vector <int> indices;
+
+ for (int i=0;i<gf->meshAttributes.corners.size();i++)
+ corners.push_back(gf->meshAttributes.corners[i]->mesh_vertices[0]);
+
+ computeEdgeLoops (gf,d_vertices, indices);
+
+ if (corners.size () != 3 && corners.size () != 4)
+ Msg (GERROR,"Surface %d is transfinite but has %d corners",gf->tag(),corners.size());
+ if (indices.size () != 2)
+ Msg (GERROR,"Surface %d is transfinite but has %d holes",gf->tag(),indices.size()-2);
+
+ std::vector <MVertex *> m_vertices;
+ int I;
+ for (I=0;I<d_vertices.size();I++)
+ if(d_vertices[I] == corners[0])break;
+ for (int j=0;j<d_vertices.size();j++)
+ m_vertices.push_back(d_vertices[(I+j)%d_vertices.size()]);
+
+ int iCorner = 0;
+ int N[4];
+ std::vector<double> U;
+ std::vector<double> V;
+
+ std::map<std::pair<int,int> , MVertex*> tab;
+
+ for (int i=0;i<m_vertices.size();i++)
+ {
+ MVertex *v = m_vertices[i];
+ if (v==corners[0]||v==corners[1]||v==corners[2]|| (corners.size()==4 && v==corners[3]))
+ {
+ N[iCorner++] = i;
+ if (iCorner > 4)
+ Msg (GERROR,"Surface %d transfinite parameters are incoherent",gf->tag());
+ }
+ SPoint2 param = gf->parFromPoint (SPoint3(v->x(),v->y(),v->z()));
+ U.push_back(param.x());
+ V.push_back(param.y());
+ }
+
+ int N1 = N[0];
+ int N2 = N[1];
+ int N3 = N[2];
+ int N4 = N[3];
+
+ int L = N2-N1;
+ int H = N3-N2;
+
+ std::vector<double> lengths_i;
+ std::vector<double> lengths_j;
+
+ double L_i=0;
+ double L_j=0;
+
+ lengths_i.push_back( 0. );
+ lengths_j.push_back( 0. );
+ for (int i=0;i<L;i++)
+ {
+ MVertex *v1 = m_vertices[i];
+ MVertex *v2 = m_vertices[i+1];
+ double l = sqrt ((v1->x()-v2->x())*(v1->x()-v2->x())+
+ (v1->y()-v2->y())*(v1->y()-v2->y())+
+ (v1->z()-v2->z())*(v1->z()-v2->z()));
+ L_i+= l;
+ lengths_i.push_back( L_i );
+ }
+ for (int i=L;i<L+H;i++)
+ {
+ MVertex *v1 = m_vertices[i];
+ MVertex *v2 = m_vertices[i+1];
+ double l = sqrt ((v1->x()-v2->x())*(v1->x()-v2->x())+
+ (v1->y()-v2->y())*(v1->y()-v2->y())+
+ (v1->z()-v2->z())*(v1->z()-v2->z()));
+ L_j+= l;
+ lengths_j.push_back( L_j );
+ }
+
+
+ //Msg (INFO,"L %d H %d -- %d -- %d %d %d %d",L,H,m_vertices.size(),N1,N2,N3,N4);
+
+ /*
+ 2L+H
+ +------------+ L+H
+ | |
+ | |
+ | |
+2L+2H+2 | |
+ +------------+
+ 0 L
+
+ */
+
+ tab[std::make_pair(0,0)] = m_vertices[0];
+ tab[std::make_pair(L,0)] = m_vertices[L];
+ tab[std::make_pair(L,H)] = m_vertices[L+H];
+ tab[std::make_pair(0,H)] = m_vertices[2*L+H];
+
+ for (int i=1;i<L;i++)
+ {
+ tab[std::make_pair(i,0)] = m_vertices[i];
+ tab[std::make_pair(i,H)] = m_vertices[2*L+H-i];
+ }
+ for (int i=1;i<H;i++)
+ {
+ tab[std::make_pair(L,i)] = m_vertices[L+i];
+ tab[std::make_pair(0,i)] = m_vertices[2*L+2*H-i];
+ }
+
+ double UC1 = U[N1];
+ double UC2 = U[N2];
+ double UC3 = U[N3];
+ double UC4 = U[N4];
+ double VC1 = V[N1];
+ double VC2 = V[N2];
+ double VC3 = V[N3];
+ double VC4 = V[N4];
+
+ for(int i = 1; i<L; i++)
+ {
+ double u = lengths_i[i]/L_i;
+ for(int j = 1; j < H; j++)
+ {
+ double v = lengths_j[j]/L_j;
+ int iP1 = N1+i;
+ int iP2 = N2+j;
+ int iP3 = N4-i;
+ int iP4 = (N4+(N3-N2)-j)%m_vertices.size();
+
+ double Up = TRAN_QUA ( U[iP1], U[iP2], U[iP3], U[iP4] , UC1, UC2, UC3, UC4, u, v );
+ double Vp = TRAN_QUA ( V[iP1], V[iP2], V[iP3], V[iP4] , VC1, VC2, VC3, VC4, u, v );
+
+ GPoint gp = gf->point (SPoint2(Up,Vp));
+ MFaceVertex *newv = new MFaceVertex ( gp.x(),gp.y(),gp.z(), gf, Up, Vp );
+ gf->mesh_vertices.push_back (newv);
+ tab[std::make_pair(i,j)] = newv;
+ }
+ }
+
+
+ for (int IT = 0;IT< 5;IT++)
+ {
+ for(int i = 1; i<L; i++)
+ {
+ for(int j = 1; j < H; j++)
+ {
+ MVertex *v1 = tab[std::make_pair(i,j)];
+ MVertex *v2 = tab[std::make_pair(i+1,j)];
+ MVertex *v3 = tab[std::make_pair(i-1,j)];
+ MVertex *v4 = tab[std::make_pair(i,j+1)];
+ MVertex *v5 = tab[std::make_pair(i,j-1)];
+
+ v1->x() = 0.25 * (v2->x() +v3->x() +v4->x() +v5->x());
+ v1->y() = 0.25 * (v2->y() +v3->y() +v4->y() +v5->y());
+ v1->z() = 0.25 * (v2->z() +v3->z() +v4->z() +v5->z());
+ }
+ }
+ }
+
+ for(int i = 0; i < L ; i++)
+ {
+ for(int j = 0; j < H; j++)
+ {
+ MVertex *v1 = tab[std::make_pair(i,j)];
+ MVertex *v2 = tab[std::make_pair(i+1,j)];
+ MVertex *v3 = tab[std::make_pair(i+1,j+1)];
+ MVertex *v4 = tab[std::make_pair(i,j+1)];
+ if (gf->meshAttributes.recombine)
+ gf->quadrangles.push_back(new MQuadrangle (v1,v2,v3,v4));
+ else if ( gf->meshAttributes.transfiniteArrangement == 1 ||
+ (gf->meshAttributes.transfiniteArrangement == 0 &&
+ (( i % 2 == 0 && j % 2 == 1 ) ||
+ ( i % 2 == 1 && j % 2 == 0 ) ) ) )
+ {
+ gf->triangles.push_back(new MTriangle (v1,v2,v3));
+ gf->triangles.push_back(new MTriangle (v3,v4,v1));
+ }
+ else
+ {
+ gf->triangles.push_back(new MTriangle (v1,v2,v4));
+ gf->triangles.push_back(new MTriangle (v4,v2,v3));
+ }
+ }
+ }
+ return 1;
+}
+
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 42ea7dc5f595b07a0160215f385c15065261a2f0..1c27b806c9bb443656baf554d6946a0ebc230843 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -204,6 +204,12 @@ void meshNormalsPointOutOfTheRegion (GRegion *gr)
double v2[3] = {X[2]-X[1],Y[2]-Y[1],Z[2]-Z[1]};
double N[3] ;
prodve ( v1, v2 , N );
+ norme (v1);
+ norme (v2);
+ norme(N);
+ N[0]+= 0.05*v1[0]-0.05*v2[0];
+ N[1]+= 0.03*v1[1]-0.15*v2[1];
+ N[2]+= 0.05*v1[2]-0.05*v2[2];
norme(N);
std::list<GFace*>::iterator it_b = faces.begin();
while (it_b != faces.end())
diff --git a/benchmarks/2d/transfinite_b.geo b/benchmarks/2d/transfinite_b.geo
new file mode 100644
index 0000000000000000000000000000000000000000..8efc8c45351a05dde8707afd9f3250c640ebb243
--- /dev/null
+++ b/benchmarks/2d/transfinite_b.geo
@@ -0,0 +1,19 @@
+Point(1) = {0,0,0,0.1};
+Point(2) = {1,0,0,0.1};
+Point(3) = {0,1,0,0.1};
+Point(4) = {3,0,0,0.1};
+Point(5) = {3,3,0,0.1};
+Point(6) = {0,3,0,0.1};
+Line(1) = {3,6};
+Line(2) = {6,5};
+Line(3) = {5,4};
+Line(4) = {4,2};
+Circle(5) = {2,1,3};
+Line Loop(6) = {3,4,5,1,2};
+Plane Surface(7) = {6};
+Transfinite Line {1} = 20 Using Progression 1.2;
+Transfinite Line {4} = 20 Using Progression 1./1.2;
+Transfinite Line {3,2} = 10 Using Progression 1;
+Transfinite Line {5} = 19 Using Progression 1;
+Transfinite Surface {7} = {3,6,4,2} Alternated;
+Recombine Surface {7};