diff --git a/Geo/MElement.h b/Geo/MElement.h
index 0d8bdc168eb54fa49eb29d6fa72e21fad80f1756..248a3d3929235c2d4abe78d9df15dc5eefff70d1 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -347,6 +347,24 @@ class MTriangle6 : public MTriangle {
}
};
+struct compareMTriangleLexicographic
+{
+ bool operator () ( MTriangle *t1 , MTriangle *t2 ) const
+ {
+ MVertex *_v1[3] = {t1->getVertex(0),t1->getVertex(1),t1->getVertex(2)};
+ MVertex *_v2[3] = {t2->getVertex(0),t2->getVertex(1),t2->getVertex(2)};
+ std::sort(_v1,_v1+3);
+ std::sort(_v2,_v2+3);
+ if (_v1[0] < _v2[0]) return true;
+ if (_v1[0] > _v2[0]) return false;
+ if (_v1[1] < _v2[1]) return true;
+ if (_v1[1] > _v2[1]) return false;
+ if (_v1[2] < _v2[2]) return true;
+ return false;
+ }
+};
+
+
class MQuadrangle : public MElement {
protected:
MVertex *_v[4];
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 2ae5a72370c33828d1092a584c9565d5f4fbeefa..6c729257e3297befe70af72db97abb375b88eac5 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGFace.cpp,v 1.48 2007-01-16 11:31:41 geuzaine Exp $
+// $Id: meshGFace.cpp,v 1.49 2007-01-16 14:19:31 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -314,7 +314,8 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
double lone = NewGetLc ( *it);
const double coord = 0.5;
- if (lone < 1.e-12 && computeParametricEdgeLength((*it)->p1,(*it)->p2) > 1.e-5) lone = 2;
+ // take care with seams :
+ if (lone < 1.e-10 && computeParametricEdgeLength((*it)->p1,(*it)->p2) > 1.e-5) lone = 2;
if ((*it)->numfaces() == 2 && (lone > 1.3))
{
BDS_Point *mid ;
@@ -364,7 +365,7 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
{
if (NN2++ >= NN1)break;
double lone = NewGetLc ( *it);
- if (lone < 1.e-12 && computeParametricEdgeLength((*it)->p1,(*it)->p2) > 1.e-5) lone = 2;
+ if (lone < 1.e-10 && computeParametricEdgeLength((*it)->p1,(*it)->p2) > 1.e-5) lone = 2;
if (!(*it)->deleted && (*it)->numfaces() == 2 && lone < 0.6 )
{
@@ -920,6 +921,9 @@ bool buildConsecutiveListOfVertices ( GFace *gf,
std::vector<SPoint2> mesh1d_seam;
bool seam = ges.ge->isSeam(gf);
+
+ printf("face %d edge %d seam %d (%d %d)\n",gf->tag(),ges.ge->tag(),seam,ges.ge->getBeginVertex()->tag(),ges.ge->getEndVertex()->tag());
+
Range<double> range = ges.ge->parBounds(0);
MVertex *here = ges.ge->getBeginVertex()->mesh_vertices[0];
@@ -1196,6 +1200,9 @@ bool gmsh2DMeshGeneratorPeriodic ( GFace *gf )
free (doc.points);
free (doc.delaunay);
+ char name[245];
+ sprintf(name,"param%d.pos",gf->tag());
+ outputScalarField(m->triangles, name,1);
// Recover the boundary edges
@@ -1332,7 +1339,12 @@ bool gmsh2DMeshGeneratorPeriodic ( GFace *gf )
MVertex *v2 = recover_map[n[1]];
MVertex *v3 = recover_map[n[2]];
if (!n[3])
- gf->triangles.push_back(new MTriangle (v1,v2,v3) );
+ {
+ // when a singular point is present, degenerated triangles may be created,
+ // for example on a sphere that contains one pole
+ if (v1 != v2 && v1 != v3 && v2 != v3)
+ gf->triangles.push_back(new MTriangle (v1,v2,v3) );
+ }
else
{
MVertex *v4 = recover_map[n[3]];
@@ -1345,9 +1357,6 @@ bool gmsh2DMeshGeneratorPeriodic ( GFace *gf )
// delete the mesh
- // char name[245];
- // sprintf(name,"param%d.pos",gf->tag());
- //outputScalarField(m->triangles, name,1);
// sprintf(name,"real%d.pos",gf->tag());
//outputScalarField(m->triangles, name,0);
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 2beeccc3e4c80645c33138b44de8fb6f52f73967..0b28eed08c8e51b1ef60b167abb1cb040300b1f1 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGRegion.cpp,v 1.23 2007-01-16 11:31:42 geuzaine Exp $
+// $Id: meshGRegion.cpp,v 1.24 2007-01-16 14:19:31 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -111,7 +111,9 @@ void buildTetgenStructure(GRegion *gr, tetgenio &in, std::vector<MVertex*> &numb
}
}
-void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out,
+void TransferTetgenMesh(GRegion *gr,
+ tetgenio &in,
+ tetgenio &out,
std::vector<MVertex*> &numberedV)
{
for(int i = numberedV.size(); i < out.numberofpoints; i++){
@@ -439,10 +441,13 @@ void meshGRegion::operator() (GRegion *gr)
Msg(STATUS2, "Meshing volume %d", gr->tag());
// destroy the mesh if it exists
- deMeshGRegion dem;
- dem(gr);
-
- if(MeshTransfiniteVolume(gr)) return;
+ if(gr->meshAttributes.Method == TRANSFINI)
+ {
+ deMeshGRegion dem;
+ dem(gr);
+ MeshTransfiniteVolume(gr);
+ return;
+ }
std::list<GFace*> faces = gr->faces();
@@ -458,9 +463,15 @@ void meshGRegion::operator() (GRegion *gr)
#if !defined(HAVE_TETGEN)
Msg(GERROR, "Tetgen is not compiled in this version of Gmsh");
#else
- // put all the faces in the same model
+ // delete the mesh for all regions
GModel::riter rit = gr->model()->firstRegion() ;
if (gr != *rit)return;
+ for (; rit != gr->model()->lastRegion();++rit)
+ {
+ deMeshGRegion dem;
+ dem(*rit);
+ }
+ // put all the faces in the same model
std::list<GFace*> allFaces;
GModel::fiter fit = gr->model()->firstFace() ;
while (fit != gr->model()->lastFace()){
@@ -476,14 +487,24 @@ void meshGRegion::operator() (GRegion *gr)
sprintf(opts, "pe%c", (CTX.verbosity < 3) ? 'Q': (CTX.verbosity > 6)? 'V': '\0');
tetrahedralize(opts, &in, &out);
TransferTetgenMesh(gr, in, out, numberedV);
+ // sort triangles in all model faces in order to be able to search in vectors
+ {
+ std::list<GFace*>::iterator itf = allFaces.begin();
+ while(itf!=allFaces.end())
+ {
+ compareMTriangleLexicographic cmp;
+ std::sort((*itf)->triangles.begin(),
+ (*itf)->triangles.end(),
+ cmp);
+ ++itf;
+ }
+ }
- // FIXME: all the volume nodes+tets now belong to the first
- // region--we need to recategorize them into each separate region
+ // restore the initial set of faces
+ gr->set(faces);
// now do insertion of points
insertVerticesInRegion(gr);
- // restore the initial set of faces
- gr->set(faces);
// meshNormalsPointOutOfTheRegion(gr);
#endif
}
@@ -492,6 +513,8 @@ void meshGRegion::operator() (GRegion *gr)
#if !defined(HAVE_NETGEN)
Msg(GERROR, "Netgen is not compiled in this version of Gmsh");
#else
+ deMeshGRegion dem;
+ dem(gr);
// orient the triangles of with respect to this region
meshNormalsPointOutOfTheRegion(gr);
std::vector<MVertex*> numberedV;
diff --git a/Mesh/meshGRegionDelaunayInsertion.cpp b/Mesh/meshGRegionDelaunayInsertion.cpp
index aec76e24f8adb9bbce9d317346cf2697ec842a09..6684fedca997df6f767a074023b857631f86f755 100644
--- a/Mesh/meshGRegionDelaunayInsertion.cpp
+++ b/Mesh/meshGRegionDelaunayInsertion.cpp
@@ -1,4 +1,4 @@
-// $Id: meshGRegionDelaunayInsertion.cpp,v 1.10 2007-01-12 13:16:59 remacle Exp $
+// $Id: meshGRegionDelaunayInsertion.cpp,v 1.11 2007-01-16 14:19:31 remacle Exp $
//
// Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
//
@@ -21,6 +21,7 @@
#include "BackgroundMesh.h"
#include "meshGRegionDelaunayInsertion.h"
+#include "GModel.h"
#include "GRegion.h"
#include "Numeric.h"
#include "Message.h"
@@ -38,7 +39,7 @@ int MTet4::inCircumSphere ( const double *p ) const
return (result > 0) ? 1 : 0;
}
-int faces[4][3] = {{0,1,2},{0,2,3},{0,1,3},{1,2,3}};
+static int faces[4][3] = {{0,1,2},{0,2,3},{0,1,3},{1,2,3}};
struct faceXtet
{
@@ -116,10 +117,10 @@ void recurFindCavity ( std::list<faceXtet> & shell,
MTet4 *neigh = t->getNeigh(i) ;
if (!neigh)
shell.push_back ( faceXtet ( t, i ) );
- else if (!neigh->isDeleted())
+ else if (!neigh->isDeleted() )
{
int circ = neigh->inCircumSphere ( v );
- if (circ)
+ if (circ && (neigh->onWhat() == t->onWhat()))
recurFindCavity ( shell, cavity,v , neigh);
else
shell.push_back ( faceXtet ( t, i ) );
@@ -203,10 +204,10 @@ bool insertVertex (MVertex *v ,
while (it != shell.end())
{
- MTetrahedron *t = new MTetrahedron ( it->v[0],
- it->v[1],
- it->v[2],
- v);
+ MTetrahedron *tr = new MTetrahedron ( it->v[0],
+ it->v[1],
+ it->v[2],
+ v);
// Msg(INFO,"shell %d %d %d",it->v[0]->getNum(),it->v[1]->getNum(),it->v[2]->getNum());
// fprintf(ff,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g) {0,0,0};\n",
// it->v[0]->x(),
@@ -219,7 +220,8 @@ bool insertVertex (MVertex *v ,
// it->v[2]->y(),
// it->v[2]->z());
- MTet4 *t4 = new MTet4 ( t , vSizes);
+ MTet4 *t4 = new MTet4 ( tr , vSizes);
+ t4->setOnWhat(t->onWhat());
newTets[k++]=t4;
// all new tets are pushed front in order to
// ba able to destroy them if the cavity is not
@@ -297,6 +299,80 @@ static void setLcs ( MTetrahedron *t, std::map<MVertex*,double> &vSizes)
// }
+bool find_triangle_in_model ( GModel *model , MTriangle *tri, GFace **gfound)
+{
+ GModel::fiter fit = model->firstFace() ;
+ compareMTriangleLexicographic cmp;
+ while (fit != model->lastFace()){
+
+ bool found = binary_search((*fit)->triangles.begin(),
+ (*fit)->triangles.end(),
+ tri,cmp);
+ if (found )
+ {
+ *gfound = *fit;
+ return true;
+ }
+ ++fit;
+ }
+ return false;
+
+}
+
+GRegion *getRegionFromBoundingFaces (GModel *model,
+ std::set<GFace *> &faces_bound)
+{
+ GModel::riter git = model->firstRegion() ;
+ while (git != model->lastRegion()){
+
+ std::list <GFace *> _faces = (*git)->faces();
+ if (_faces.size() == faces_bound.size())
+ {
+ bool ok = true;
+ for (std::list <GFace *>::iterator it = _faces.begin();it != _faces.end();++it)
+ {
+ if (faces_bound.find (*it) == faces_bound.end())ok = false;
+ }
+ if (ok)return *git;
+ }
+ ++git;
+ }
+ return 0;
+}
+
+
+void recur_classify ( MTet4 *t ,
+ std::set<MTet4*,compareTet4Ptr> &theRegion,
+ std::set<GFace *> &faces_bound,
+ GRegion *bidon ,
+ GModel *model)
+{
+ if (!t) Msg (GERROR,"a tet is not connected by a boundary face");
+ if (t->onWhat())return;
+ theRegion.insert(t);
+ t->setOnWhat(bidon);
+
+ for (int i=0;i<4;i++)
+ {
+ MTriangle tri ( t->tet()->getVertex ( faces[i][0] ),
+ t->tet()->getVertex ( faces[i][1] ),
+ t->tet()->getVertex ( faces[i][2] ) );
+ GFace *gfound;
+ bool found = find_triangle_in_model ( model , &tri, &gfound);
+
+ // Msg(INFO,"found ? %d",found);
+ if (found)
+ {
+ faces_bound.insert(gfound);
+ }
+ else
+ {
+ recur_classify ( t->getNeigh(i) , theRegion, faces_bound, bidon, model );
+ }
+ }
+
+}
+
void insertVerticesInRegion (GRegion *gr)
{
@@ -319,6 +395,23 @@ void insertVerticesInRegion (GRegion *gr)
gr->tetrahedra.clear();
connectTets ( allTets.begin(), allTets.end() );
+ // classify the tets on the right region
+ Msg (INFO,"reclassifying %d tets",allTets.size());
+ for (std::set<MTet4*,compareTet4Ptr>::iterator it = allTets.begin();it!=allTets.end();++it)
+ {
+ if (!(*it)->onWhat())
+ {
+ std::set<MTet4*,compareTet4Ptr> theRegion;
+ std::set<GFace *> faces_bound;
+ GRegion *bidon;
+ recur_classify ( *it , theRegion, faces_bound, bidon , gr->model());
+ GRegion *myGRegion = getRegionFromBoundingFaces (gr->model() , faces_bound );
+ Msg (INFO,"found a region with %d tets %p",theRegion.size(),myGRegion);
+ for (std::set<MTet4*,compareTet4Ptr>::iterator it2 = theRegion.begin();it2!=theRegion.end();++it2)(*it2)->setOnWhat(myGRegion);
+ }
+ }
+
+
Msg(DEBUG,"All %d tets were connected",allTets.size());
// here the classification should be done
@@ -389,7 +482,7 @@ void insertVerticesInRegion (GRegion *gr)
}
else
{
- gr->tetrahedra.push_back(worst->tet());
+ worst->onWhat()->tetrahedra.push_back(worst->tet());
}
delete worst;
allTets.erase(allTets.begin());
diff --git a/Mesh/meshGRegionDelaunayInsertion.h b/Mesh/meshGRegionDelaunayInsertion.h
index 24a1e821c6a724468e8967ffd33d5947e1299723..313f2609f618314f9bb416ab51234c2b441b10ff 100644
--- a/Mesh/meshGRegionDelaunayInsertion.h
+++ b/Mesh/meshGRegionDelaunayInsertion.h
@@ -32,13 +32,18 @@ class MTet4
double circum_radius;
MTetrahedron *base;
MTet4 *neigh[4];
+ GRegion *gr;
public :
+
+ inline GRegion * onWhat () const {return gr;}
+ inline void setOnWhat (GRegion *g) {gr=g;}
bool isDeleted () const {return deleted;}
void forceRadius (double r){circum_radius=r;}
double getRadius ()const {return circum_radius;}
- MTet4 ( MTetrahedron * t, std::vector<double> & sizes) : deleted(false), base (t)
+
+ MTet4 ( MTetrahedron * t, std::vector<double> & sizes) : deleted(false), base (t), gr(0)
{
neigh[0] = neigh[1] = neigh[2] = neigh[3] = 0;
double center[3];