diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index 62d046b4590d66a7e0b02cdb22ed3019a892f256..3b659a7ef9a2660769aad6d1157fc22ec8c2d5e2 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -975,7 +975,7 @@ StringXNumber MeshOptions_Number[] = {
"Only display elements whose longest edge is greater than RadiusInf" },
{ F|O, "RadiusSup" , opt_mesh_radius_sup , 0.0 ,
"Only display elements whose longest edge is smaller than RadiusSup" },
- { F|O, "RandomFactor" , opt_mesh_rand_factor , 1.e-4 ,
+ { F|O, "RandomFactor" , opt_mesh_rand_factor , 1.e-10 ,
"Random factor used in 2D and 3D meshing algorithm (test other values when the algorithm fails)" },
{ F|O, "RecombineAlgo" , opt_mesh_recombine_algo , 1 ,
"Recombine algorithm (1=mixed triangles-quadrangles, 2=all quadrangles)" },
diff --git a/Geo/GRegion.cpp b/Geo/GRegion.cpp
index a5cd03c822920e6697401797f83415529705878a..9273ded37ee46d2ce1d2532ffd3203e8f4330905 100644
--- a/Geo/GRegion.cpp
+++ b/Geo/GRegion.cpp
@@ -1,4 +1,4 @@
-// $Id: GRegion.cpp,v 1.10 2006-08-26 13:34:46 geuzaine Exp $
+// $Id: GRegion.cpp,v 1.11 2006-09-11 15:23:54 remacle Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -110,3 +110,22 @@ void GRegion::deleteMeshPartitions()
for(unsigned int i = 0; i < pyramids.size(); i++)
pyramids[i]->setPartition(0);
}
+std::list<GEdge*> GRegion::edges() const
+{
+ std::list<GEdge*> e;
+ std::list<GFace *>::const_iterator it = l_faces.begin();
+ while (it != l_faces.end())
+ {
+ std::list<GEdge*> e2;
+ e2 = (*it)->edges();
+ std::list<GEdge*>::const_iterator it2 = e2.begin();
+ while (it2 != e2.end())
+ {
+ if (std::find(e.begin(),e.end(),*it2) == e.end())
+ e.push_back(*it2);
+ ++it2;
+ }
+ ++it;
+ }
+ return e;
+}
diff --git a/Geo/GRegion.h b/Geo/GRegion.h
index ebc66decf3c11c7c1d036d49b6bc95c2315d6e3c..2e0fe1d484d3db13d7172c7fdc991d835b2a7ef2 100644
--- a/Geo/GRegion.h
+++ b/Geo/GRegion.h
@@ -36,6 +36,9 @@ class GRegion : public GEntity {
virtual int dim() const {return 3;}
virtual void setVisibility(char val, bool recursive=false);
virtual std::list<GFace*> faces() const{return l_faces;}
+ virtual std::list<GEdge*> edges() const;
+ void set(std::list<GFace*> &f) {l_faces= f;}
+
// The bounding box
virtual SBoundingBox3d bounds() const;
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 38d5d31157c30da7f3690d3a7ea99b25f17222b2..4cea2e14deba5fa83882cb8f2e20f8d0a4bda995 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -995,4 +995,5 @@ class MPyramid14 : public MPyramid {
}
};
+
#endif
diff --git a/Mesh/2D_DivAndConq.cpp b/Mesh/2D_DivAndConq.cpp
index bdb0131b99b0b87e1b0f381f178bc34dd708c8be..26b55470a4f4988ecbbdf1257be5e7edc5dc6344 100644
--- a/Mesh/2D_DivAndConq.cpp
+++ b/Mesh/2D_DivAndConq.cpp
@@ -1,4 +1,4 @@
-// $Id: 2D_DivAndConq.cpp,v 1.21 2006-01-06 00:34:25 geuzaine Exp $
+// $Id: 2D_DivAndConq.cpp,v 1.22 2006-09-11 15:23:54 remacle Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -183,7 +183,7 @@ Segment UpperCommonTangent(DT vl, DT vr)
/* return 1 if the point k is NOT in the circumcircle of
triangle hij */
-
+#ifndef HAVE_TETGEN
int Qtest(PointNumero h, PointNumero i, PointNumero j, PointNumero k)
{
double xc, yc, rcarre, distca;
@@ -209,6 +209,28 @@ int Qtest(PointNumero h, PointNumero i, PointNumero j, PointNumero k)
else
return (0); /* point not in circle, because no circle ! */
}
+#else
+#include "tetgen.h"
+int Qtest(PointNumero h, PointNumero i, PointNumero j, PointNumero k)
+{
+ exactinit() ;
+ if((h == i) && (h == j) && (h == k)) {
+ Msg(GERROR, "3 identical points in Qtest");
+ return (0); /* returning 1 will cause looping for ever */
+ }
+
+ REAL pa[2] = {(REAL)pPointArray[h].where.h,(REAL)pPointArray[h].where.v};
+ REAL pb[2] = {(REAL)pPointArray[i].where.h,(REAL)pPointArray[i].where.v};
+ REAL pc[2] = {(REAL)pPointArray[j].where.h,(REAL)pPointArray[j].where.v};
+ REAL pd[2] = {(REAL)pPointArray[k].where.h,(REAL)pPointArray[k].where.v};
+
+ REAL result = incircle ( pa, pb, pc, pd ) * orient2d (pa,pb,pc);
+ // Msg(INFO, "Qtest %12.5E",result);
+
+ return (result < 0)?1:0;
+}
+#endif
+
int merge(DT vl, DT vr)
{
diff --git a/Mesh/Makefile b/Mesh/Makefile
index 78b0288b0e6140d0417641fff53df7bf1d28973e..229ed8216404e6c66f00d52cc8941acbf0e288d1 100644
--- a/Mesh/Makefile
+++ b/Mesh/Makefile
@@ -1,4 +1,4 @@
-# $Id: Makefile,v 1.130 2006-09-07 04:54:55 geuzaine Exp $
+# $Id: Makefile,v 1.131 2006-09-11 15:23:54 remacle Exp $
#
# Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
#
@@ -66,6 +66,7 @@ SRC = 1D_Mesh.cpp \
meshGEdge.cpp \
meshGFace.cpp \
meshGFaceTransfinite.cpp \
+ meshGFaceDelaunayInstertion.cpp \
meshGRegion.cpp \
Nurbs.cpp \
Interpolation.cpp \
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index dcbc2c7585363f423802fe8c045bcda67bc7a414..8c1067414dfc0a7e65685bf877ac5c0c45a882d4 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFace.cpp,v 1.12 2006-09-08 02:39:43 geuzaine Exp $
+// $Id: meshGFace.cpp,v 1.13 2006-09-11 15:23:54 remacle Exp $
//
// Copyright (C) 1997-2006 C. Geuzaine, J.-F. Remacle
//
@@ -443,7 +443,7 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
}
}
IT++;
- // Msg(INFO," %d split %d swap %d collapse %d smooth",nb_split,nb_swap,nb_collaps,nb_smooth);
+ Msg(INFO," iter %d minL %g maxL %g %d split %d swap %d collapse %d smooth",IT,minL,maxL,nb_split,nb_swap,nb_collaps,nb_smooth);
m.cleanup();
}
}
@@ -895,7 +895,7 @@ void gmsh2DMeshGenerator ( GFace *gf )
// goto hhh;
// start mesh generation
- RefineMesh (gf,*m,20);
+ RefineMesh (gf,*m,15);
// OptimizeMesh (gf,*m,2);
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 92d678ffe16a262866f3eb4f1203a3235af5e1b2..72023179e2fb97ff6563dc80cec2d7500ad23ee3 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -1,10 +1,15 @@
#include "meshGRegion.h"
+#include "GModel.h"
#include "GRegion.h"
#include "GFace.h"
+#include "GEdge.h"
#include "BDS.h"
#include "Message.h"
#include <vector>
+
+
+
void getAllBoundingVertices (GRegion *gr, int importVolumeMesh, std::set<MVertex*> &allBoundingVertices )
{
std::list<GFace*> faces = gr->faces();
@@ -19,9 +24,9 @@ void getAllBoundingVertices (GRegion *gr, int importVolumeMesh, std::set<MVertex
for (unsigned int i = 0; i< gf->triangles.size(); i++)
{
MTriangle *t = gf->triangles[i];
- allBoundingVertices.insert (t->getVertex(0));
- allBoundingVertices.insert (t->getVertex(1));
- allBoundingVertices.insert (t->getVertex(2));
+ for (int k=0;k<3;k++)
+ if (allBoundingVertices.find(t->getVertex(k)) == allBoundingVertices.end())
+ allBoundingVertices.insert (t->getVertex(k));
}
++it;
}
@@ -66,8 +71,6 @@ void buildTetgenStructure ( GRegion *gr, tetgenio &in, std::vector<MVertex*> &
++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];
@@ -84,7 +87,8 @@ void buildTetgenStructure ( GRegion *gr, tetgenio &in, std::vector<MVertex*> &
tetgenio::init(f);
f->numberofholes = 0;
f->numberofpolygons = 1;
- tetgenio::polygon *p = f->polygonlist = new tetgenio::polygon[f->numberofpolygons];
+ f->polygonlist = new tetgenio::polygon[f->numberofpolygons];
+ tetgenio::polygon *p = &f->polygonlist[0];
tetgenio::init(p);
p->numberofvertices = 3;
p->vertexlist = new int[p->numberofvertices];
@@ -97,17 +101,75 @@ void buildTetgenStructure ( GRegion *gr, tetgenio &in, std::vector<MVertex*> &
}
++it;
}
+
+
}
void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out, std::vector<MVertex*> & numberedV)
{
- for (int i = 0; i < out.numberofpoints; i++)
+ int I = numberedV.size() + 1;
+ for (int i = numberedV.size(); i < out.numberofpoints; i++)
{
- MVertex *v = new MVertex (out.pointlist[i * 3 + 0],out.pointlist[i * 3 + 1],out.pointlist[i * 3 + 2]);
+ MVertex *v = new MVertex (out.pointlist[i * 3 + 0],
+ out.pointlist[i * 3 + 1],
+ out.pointlist[i * 3 + 2], gr);
gr->mesh_vertices.push_back(v);
numberedV.push_back(v);
}
+ Msg(INFO,"%d points %d edges and %d faces in the final mesh",out.numberofpoints,out.numberofedges,out.numberoftrifaces);
+ // Tetgen modifies both surface & edge mesh. So, we recreate all mes
+
+ {
+ std::list<GFace*> faces = gr->faces();
+ std::list<GFace*>::iterator it = faces.begin();
+ while (it != faces.end())
+ {
+ GFace *gf = (*it);
+ for (int i=0;i<gf->triangles.size();i++)
+ {
+ delete gf->triangles[i];
+ }
+ gf->triangles.clear();
+ ++it;
+ }
+ }
+
+ // re create 1D mesh
+ // 2 be done ...
+ for (int i = 0; i < out.numberofedges; i++)
+ {
+ MVertex *v[2];
+ v[0] = numberedV[out.edgelist[i * 2 + 0] -1];
+ v[1] = numberedV[out.edgelist[i * 2 + 1] -1];
+ }
+
+ // re-create the triangular meshes
+ for (int i = 0; i < out.numberoftrifaces; i++)
+ {
+ MVertex *v[3];
+ v[0] = numberedV[out.trifacelist[i * 3 + 0] -1];
+ v[1] = numberedV[out.trifacelist[i * 3 + 1] -1];
+ v[2] = numberedV[out.trifacelist[i * 3 + 2] -1];
+ GFace *gf = gr->model()->faceByTag ( out.trifacemarkerlist[i] );
+ for (int j=0;j<3;j++)
+ {
+ if ( v[j]->onWhat()->dim() == 3 )
+ {
+ v[j]->onWhat()->mesh_vertices.erase(std::find(v[j]->onWhat()->mesh_vertices.begin(),
+ v[j]->onWhat()->mesh_vertices.end(),
+ v[j]));
+ // SPoint2 pp = gf->parFromPoint(SPoint3 (v[j]->x(),v[j]->y(),v[j]->z()));
+ SPoint2 pp (0,0);
+ MFaceVertex *v1b = new MFaceVertex (v[j]->x(),v[j]->y(),v[j]->z(),gf,pp.x(),pp.y() );
+ numberedV[out.trifacelist[i * 3 + j] -1] = v1b;
+ delete v[j];
+ v[j] = v1b;
+ }
+ }
+ MTriangle *t = new MTriangle(v[0],v[1],v[2]);
+ gf->triangles.push_back(t);
+ }
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];
@@ -116,10 +178,10 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out, std::vector<MV
MTetrahedron *t = new MTetrahedron(v1,v2,v3,v4);
gr->tetrahedra.push_back(t);
}
-
-
-}
+
+}
+
#endif
@@ -200,7 +262,7 @@ void TransferVolumeMesh(GRegion *gr, Ng_Mesh * _ngmesh, std::vector<MVertex*> &n
{
double tmp[3];
Ng_GetPoint(_ngmesh, i+1, tmp);
- MVertex *v = new MVertex (tmp[0],tmp[1],tmp[2]);
+ MVertex *v = new MVertex (tmp[0],tmp[1],tmp[2],gr);
numberedV.push_back(v);
gr->mesh_vertices.push_back(v);
}
@@ -243,7 +305,7 @@ int intersect_line_triangle ( double X[3],
{
double mat[3][3], det;
double b[3], res[3];
- const double eps_prec = 1.e-5;
+ const double eps_prec = 1.e-9;
mat[0][0] = X[1] - X[0];
mat[0][1] = X[2] - X[0];
@@ -362,23 +424,35 @@ void meshGRegion :: operator() (GRegion *gr)
// Send a messsage to the GMSH environment
Msg(STATUS2, "Meshing volume %d", gr->tag());
- // 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
+ // put all the faces in the same model
+ GModel::riter rit = gr->model()->firstRegion() ;
+ if (gr != *rit)return;
+ std::list<GFace*> faces = gr->faces();
+ std::list<GFace*> allFaces;
+ GModel::fiter fit = gr->model()->firstFace() ;
+ while (fit != gr->model()->lastFace())
+ {
+ allFaces.push_back(*fit);
+ ++fit;
+ }
+ gr->set ( allFaces );
+ // mesh with tetgen, possibly changing the mesh on boundaries
tetgenio in, out;
std::vector<MVertex*> numberedV;
char opts[128];
buildTetgenStructure ( gr, in, numberedV);
- sprintf(opts, "pqa%f%c", (float)CTX.mesh.quality,
+ sprintf(opts, "pe%c",
(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);
+ // restore the initial set of faces
+ gr->set ( faces );
+
#endif
}
@@ -387,6 +461,8 @@ void meshGRegion :: operator() (GRegion *gr)
#if !defined(HAVE_NETGEN)
Msg(GERROR, "Netgen is not compiled in this version of Gmsh");
#else
+ // orient the triangles of with respect to this region
+ meshNormalsPointOutOfTheRegion (gr);
std::vector<MVertex*> numberedV;
Ng_Mesh * _ngmesh = buildNetgenStructure (gr, 0, numberedV);
Ng_Meshing_Parameters mp;
diff --git a/contrib/Tetgen/tetgen.h b/contrib/Tetgen/tetgen.h
index e67f87d14a42fc09d7447a7f529cba6e981bd2e9..896c3f2279614d83acf0ac35c4d489734ed599ad 100644
--- a/contrib/Tetgen/tetgen.h
+++ b/contrib/Tetgen/tetgen.h
@@ -520,7 +520,9 @@ class tetgenbehavior {
///////////////////////////////////////////////////////////////////////////////
REAL exactinit();
+REAL incircle(REAL *pa, REAL *pb, REAL *pc, REAL *pd);
REAL orient3d(REAL *pa, REAL *pb, REAL *pc, REAL *pd);
+REAL orient2d(REAL *pa, REAL *pb, REAL *pc);
REAL insphere(REAL *pa, REAL *pb, REAL *pc, REAL *pd, REAL *pe);
///////////////////////////////////////////////////////////////////////////////