*** empty log message ***

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 ;
 
 };

Geo/gmshFace.cpp

-// $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)

Mesh/Makefile

-# $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 \

Mesh/meshGFace.cpp

-// $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

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 
 {

Mesh/meshGFaceTransfinite.cpp

+// $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;
+}
+  

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())

benchmarks/2d/transfinite_b.geo

+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};