diff --git a/Common/CommandLine.cpp b/Common/CommandLine.cpp
index 737de44f56e0a6b0092328778f81c70f86bfdc09..f77e3a63154706da32886212555932f45d092e05 100644
--- a/Common/CommandLine.cpp
+++ b/Common/CommandLine.cpp
@@ -1,4 +1,4 @@
-// $Id: CommandLine.cpp,v 1.79 2006-09-03 07:44:09 geuzaine Exp $
+// $Id: CommandLine.cpp,v 1.80 2006-09-07 16:03:32 remacle Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -428,7 +428,9 @@ void Get_Options(int argc, char *argv[])
else if(!strncmp(argv[i], "netgen", 6))
CTX.mesh.algo3d = FRONTAL_NETGEN;
else if(!strncmp(argv[i], "tetgen", 6))
- CTX.mesh.algo3d = DELAUNAY_TETGEN;
+ {
+ CTX.mesh.algo3d = DELAUNAY_TETGEN;
+ }
else {
fprintf(stderr, ERROR_STR "Unknown mesh algorithm\n");
exit(1);
diff --git a/Mesh/meshGFaceTransfinite.cpp b/Mesh/meshGFaceTransfinite.cpp
index 6a4c36bc392d6d69b75c6cd965bb5d6e75c24bca..1a191019ebd5f8d733b29ed36c590611babd526d 100644
--- a/Mesh/meshGFaceTransfinite.cpp
+++ b/Mesh/meshGFaceTransfinite.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFaceTransfinite.cpp,v 1.2 2006-09-06 10:25:24 remacle Exp $
+// $Id: meshGFaceTransfinite.cpp,v 1.3 2006-09-07 16:03:32 remacle Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -86,6 +86,15 @@ int MeshTransfiniteSurface( GFace *gf)
int L = N2-N1;
int H = N3-N2;
+ int Lb = N4-N3;
+ int Hb = m_vertices.size()-N4;
+
+ if (Lb != L || Hb != H)
+ {
+ Msg(GERROR,"Surface %d cannot be meshed using the transfinite algo",gf->tag());
+ return 0;
+ }
+
std::vector<double> lengths_i;
std::vector<double> lengths_j;
@@ -155,6 +164,8 @@ int MeshTransfiniteSurface( GFace *gf)
double VC3 = V[N3];
double VC4 = V[N4];
+
+ //create points using transfinite interpolation
for(int i = 1; i<L; i++)
{
double u = lengths_i[i]/L_i;
@@ -168,7 +179,7 @@ int MeshTransfiniteSurface( GFace *gf)
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);
@@ -176,51 +187,62 @@ int MeshTransfiniteSurface( GFace *gf)
}
}
- // elliptic smoother
- for (int IT = 0;IT< 15;IT++)
+ if (gf->geomType() == GEntity::Plane) // some work 2 be done to smooth the stuff in parametric coords
+ // parametric coords are available (MFaceVertex) !
{
- for(int i = 1; i<L; i++)
+ // elliptic smoother
+ for (int IT = 0;IT< 10;IT++)
{
- for(int j = 1; j < H; j++)
+ for(int i = 1; i<L; i++)
{
- MVertex *v11 = tab[std::make_pair(i-1,j-1)];
- MVertex *v12 = tab[std::make_pair(i-1,j)];
- MVertex *v13 = tab[std::make_pair(i-1,j+1)];
- MVertex *v21 = tab[std::make_pair(i,j-1)];
- MVertex *v22 = tab[std::make_pair(i,j)];
- MVertex *v23 = tab[std::make_pair(i,j+1)];
- MVertex *v31 = tab[std::make_pair(i+1,j-1)];
- MVertex *v32 = tab[std::make_pair(i+1,j)];
- MVertex *v33 = tab[std::make_pair(i+1,j+1)];
-
- double alpha = 0.25 * (DSQR(v23->x() - v21->x()) +
- DSQR(v23->y() - v21->y()));
- double gamma = 0.25 * (DSQR(v32->x() - v12->x()) +
- DSQR(v32->y() - v12->y()));
- double beta = 0.0625 * ((v32->x() - v12->x()) *
- (v23->x() - v21->x()) +
- (v32->y() - v12->y()) *
- (v23->y() - v21->y()));
-
- v22->x() = 0.5 * (alpha * (v32->x() + v12->x()) +
- gamma * (v23->x() + v21->x()) -
- 2. * beta * (v33->x() - v13->x() -
- v31->x() + v11->x()))
- / (alpha + gamma);
- v22->y() = 0.5 * (alpha * (v32->y() + v12->y()) +
- gamma * (v23->y() + v21->y()) -
- 2. * beta * (v33->y() - v13->y() -
- v31->y() + v11->y()))
- / (alpha + gamma);
- v22->z() = 0.5 * (alpha * (v32->z() + v12->z()) +
- gamma * (v23->z() + v21->z()) -
- 2. * beta * (v33->z() - v13->z() -
- v31->z() + v11->z()))
- / (alpha + gamma);
+ for(int j = 1; j < H; j++)
+ {
+ MVertex *v11 = tab[std::make_pair(i-1,j-1)];
+ MVertex *v12 = tab[std::make_pair(i-1,j)];
+ MVertex *v13 = tab[std::make_pair(i-1,j+1)];
+ MVertex *v21 = tab[std::make_pair(i,j-1)];
+ MVertex *v22 = tab[std::make_pair(i,j)];
+ MVertex *v23 = tab[std::make_pair(i,j+1)];
+ MVertex *v31 = tab[std::make_pair(i+1,j-1)];
+ MVertex *v32 = tab[std::make_pair(i+1,j)];
+ MVertex *v33 = tab[std::make_pair(i+1,j+1)];
+
+ double alpha = 0.25 * (DSQR(v23->x() - v21->x()) +
+ DSQR(v23->y() - v21->y()) +
+ DSQR(v23->z() - v21->z())
+ );
+ double gamma = 0.25 * (DSQR(v32->x() - v12->x()) +
+ DSQR(v32->y() - v12->y()) +
+ DSQR(v32->z() - v12->z())
+ );
+ double beta = 0.0625 * ((v32->x() - v12->x()) *
+ (v23->x() - v21->x()) +
+ (v32->y() - v12->y()) *
+ (v23->y() - v21->y()) +
+ (v32->z() - v12->z()) *
+ (v23->z() - v21->z())
+ );
+
+ v22->x() = 0.5 * (alpha * (v32->x() + v12->x()) +
+ gamma * (v23->x() + v21->x()) -
+ 2. * beta * (v33->x() - v13->x() -
+ v31->x() + v11->x()))
+ / (alpha + gamma);
+ v22->y() = 0.5 * (alpha * (v32->y() + v12->y()) +
+ gamma * (v23->y() + v21->y()) -
+ 2. * beta * (v33->y() - v13->y() -
+ v31->y() + v11->y()))
+ / (alpha + gamma);
+ v22->z() = 0.5 * (alpha * (v32->z() + v12->z()) +
+ gamma * (v23->z() + v21->z()) -
+ 2. * beta * (v33->z() - v13->z() -
+ v31->z() + v11->z()))
+ / (alpha + gamma);
+ }
}
}
}
-
+ // create elements
for(int i = 0; i < L ; i++)
{
for(int j = 0; j < H; j++)
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 1c27b806c9bb443656baf554d6946a0ebc230843..b27532bc9ccaf876bef00dfa4c79cec21fc0eef2 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -5,25 +5,14 @@
#include "Message.h"
#include <vector>
-#if defined(HAVE_NETGEN)
-namespace nglib {
-#include "nglib.h"
-#include "nglib_addon.h"
-}
-using namespace nglib;
-
-Ng_Mesh * buildNetgenStructure (GRegion *gr, int importVolumeMesh, std::vector<MVertex*> & numberedV)
+void getAllBoundingVertices (GRegion *gr, int importVolumeMesh, std::set<MVertex*> &allBoundingVertices )
{
- NgAddOn_Init();
- Ng_Mesh *_ngmesh = Ng_NewMesh();
-
std::list<GFace*> faces = gr->faces();
std::list<GFace*>::iterator it = faces.begin();
- std::set<MVertex*> allBoundingVertices;
-
+
if (importVolumeMesh)
allBoundingVertices.insert ( gr->mesh_vertices.begin() , gr->mesh_vertices.end() );
-
+
while (it != faces.end())
{
GFace *gf = (*it);
@@ -36,6 +25,121 @@ Ng_Mesh * buildNetgenStructure (GRegion *gr, int importVolumeMesh, std::vector<M
}
++it;
}
+}
+
+
+#if defined(HAVE_TETGEN)
+#include "tetgen.h"
+void buildTetgenStructure ( GRegion *gr, tetgenio &in, std::vector<MVertex*> & numberedV)
+{
+
+ std::set<MVertex*> allBoundingVertices;
+ getAllBoundingVertices (gr, 0, allBoundingVertices );
+
+ in.mesh_dim = 3;
+ in.firstnumber = 1;
+
+ in.numberofpoints = allBoundingVertices.size();
+ in.pointlist = new REAL[in.numberofpoints * 3];
+ in.pointmarkerlist = NULL;
+
+ std::set<MVertex*>::iterator itv = allBoundingVertices.begin();
+ int I=1;
+ while (itv != allBoundingVertices.end())
+ {
+ in.pointlist[(I-1)*3 + 0] = (*itv)->x();
+ in.pointlist[(I-1)*3 + 1] = (*itv)->y();
+ in.pointlist[(I-1)*3 + 2] = (*itv)->z();
+ (*itv)->setNum(I++);
+ numberedV.push_back(*itv);
+ ++itv;
+ }
+
+ int nbFace = 0;
+
+ std::list<GFace*> faces = gr->faces();
+ std::list<GFace*>::iterator it = faces.begin();
+ while (it != faces.end())
+ {
+ GFace *gf = (*it);
+ nbFace += gf->triangles.size();
+ ++it;
+ }
+
+ Msg(INFO,"%d model faces -- NumFaces = %d",nbFace,faces.size());
+
+ in.numberoffacets = nbFace;
+ in.facetlist = new tetgenio::facet[in.numberoffacets];
+ in.facetmarkerlist = new int[in.numberoffacets];
+
+ it = faces.begin();
+ I= 0;
+ while (it != faces.end())
+ {
+ GFace *gf = (*it);
+ for (int i = 0; i< gf->triangles.size(); i++)
+ {
+ MTriangle *t = gf->triangles[i];
+ tetgenio::facet *f = &in.facetlist[I];
+ tetgenio::init(f);
+ f->numberofholes = 0;
+ f->numberofpolygons = 1;
+ tetgenio::polygon *p = f->polygonlist = new tetgenio::polygon[f->numberofpolygons];
+ tetgenio::init(p);
+ p->numberofvertices = 3;
+ p->vertexlist = new int[p->numberofvertices];
+ p->vertexlist[0] = t->getVertex(0)->getNum();
+ p->vertexlist[1] = t->getVertex(1)->getNum();
+ p->vertexlist[2] = t->getVertex(2)->getNum();
+ // Msg(INFO,"Faces %d = %d %d %d",I,p->vertexlist[0] ,p->vertexlist[1] ,p->vertexlist[2] );
+ in.facetmarkerlist[I] = gf->tag();
+ ++I;
+ }
+ ++it;
+ }
+}
+
+void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out, std::vector<MVertex*> & numberedV)
+{
+ int I = numberedV.size() + 1;
+ for (int i = 0; i < out.numberofpoints; i++)
+ {
+ MVertex *v = new MVertex (out.pointlist[i * 3 + 0],out.pointlist[i * 3 + 1],out.pointlist[i * 3 + 2]);
+ gr->mesh_vertices.push_back(v);
+ numberedV.push_back(v);
+ }
+
+ for (int i = 0; i < out.numberoftetrahedra; i++) {
+ MVertex *v1 = numberedV[out.tetrahedronlist[i * 4 + 0] -1];
+ MVertex *v2 = numberedV[out.tetrahedronlist[i * 4 + 1] -1];
+ MVertex *v3 = numberedV[out.tetrahedronlist[i * 4 + 2] -1];
+ MVertex *v4 = numberedV[out.tetrahedronlist[i * 4 + 3] -1];
+ MTetrahedron *t = new MTetrahedron(v1,v2,v3,v4);
+ gr->tetrahedra.push_back(t);
+ }
+
+
+}
+
+
+#endif
+
+
+#if defined(HAVE_NETGEN)
+namespace nglib {
+#include "nglib.h"
+#include "nglib_addon.h"
+}
+using namespace nglib;
+
+Ng_Mesh * buildNetgenStructure (GRegion *gr, int importVolumeMesh, std::vector<MVertex*> & numberedV)
+{
+ NgAddOn_Init();
+ Ng_Mesh *_ngmesh = Ng_NewMesh();
+
+ std::set<MVertex*> allBoundingVertices;
+ getAllBoundingVertices (gr, importVolumeMesh, allBoundingVertices );
+
std::set<MVertex*>::iterator itv = allBoundingVertices.begin();
int I=1;
while (itv != allBoundingVertices.end())
@@ -50,7 +154,8 @@ Ng_Mesh * buildNetgenStructure (GRegion *gr, int importVolumeMesh, std::vector<M
++itv;
}
- it = faces.begin();
+ std::list<GFace*> faces = gr->faces();
+ std::list<GFace*>::iterator it = faces.begin();
while (it != faces.end())
{
GFace *gf = (*it);
@@ -261,6 +366,23 @@ void meshGRegion :: operator() (GRegion *gr)
// orient the triangles of with respect to this region
meshNormalsPointOutOfTheRegion (gr);
+ if(CTX.mesh.algo3d == DELAUNAY_TETGEN)
+ {
+#if !defined(HAVE_TETGEN)
+ Msg(GERROR, "Tetgen is not compiled in this version of Gmsh");
+#else
+ tetgenio in, out;
+ std::vector<MVertex*> numberedV;
+ char opts[128];
+ buildTetgenStructure ( gr, in, numberedV);
+ sprintf(opts, "pq1.4Ya%f%c", (float)CTX.mesh.quality,
+ (CTX.verbosity < 3)? 'Q': (CTX.verbosity > 6)? 'V': '\0');
+ Msg(STATUS2, "Meshing with volume constraint %f", (float)CTX.mesh.quality);
+ tetrahedralize(opts, &in, &out);
+ TransferTetgenMesh(gr, in, out, numberedV);
+#endif
+ }
+
if(CTX.mesh.algo3d == FRONTAL_NETGEN)
{
#if !defined(HAVE_NETGEN)
diff --git a/benchmarks/3d/Cube-01.geo b/benchmarks/3d/Cube-01.geo
index 2935afb6f6c8686fac4a34170d5daa033b5d59ea..fdd3497130e5473eae6e4823aa00bd572a20b066 100644
--- a/benchmarks/3d/Cube-01.geo
+++ b/benchmarks/3d/Cube-01.geo
@@ -10,4 +10,4 @@ Line(4) = {1,4};
Line Loop(5) = {2,3,4,1};
Plane Surface(6) = {5};
Extrude Surface { 6, {0,0.0,1} };
-Mesh.Algorithm3D = 4;
+