Skip to content
Snippets Groups Projects
Commit 22cbc692 authored by Jean-François Remacle's avatar Jean-François Remacle
Browse files

*** empty log message ***

parent 0dc30a0a
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
Geo/MElement.h
+ 18
− 0
View file @ 22cbc692
...@@ -347,6 +347,24 @@ class MTriangle6 : public MTriangle { ...@@ -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 { class MQuadrangle : public MElement {
protected: protected:
MVertex *_v[4]; MVertex *_v[4];
......
Mesh/meshGFace.cpp
+ 16
− 7
View file @ 22cbc692
// $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 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
// //
...@@ -314,7 +314,8 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT) ...@@ -314,7 +314,8 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
double lone = NewGetLc ( *it); double lone = NewGetLc ( *it);
const double coord = 0.5; 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)) if ((*it)->numfaces() == 2 && (lone > 1.3))
{ {
BDS_Point *mid ; BDS_Point *mid ;
...@@ -364,7 +365,7 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT) ...@@ -364,7 +365,7 @@ void RefineMesh ( GFace *gf, BDS_Mesh &m , const int NIT)
{ {
if (NN2++ >= NN1)break; if (NN2++ >= NN1)break;
double lone = NewGetLc ( *it); 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 ) if (!(*it)->deleted && (*it)->numfaces() == 2 && lone < 0.6 )
{ {
...@@ -920,6 +921,9 @@ bool buildConsecutiveListOfVertices ( GFace *gf, ...@@ -920,6 +921,9 @@ bool buildConsecutiveListOfVertices ( GFace *gf,
std::vector<SPoint2> mesh1d_seam; std::vector<SPoint2> mesh1d_seam;
bool seam = ges.ge->isSeam(gf); 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); Range<double> range = ges.ge->parBounds(0);
MVertex *here = ges.ge->getBeginVertex()->mesh_vertices[0]; MVertex *here = ges.ge->getBeginVertex()->mesh_vertices[0];
...@@ -1196,6 +1200,9 @@ bool gmsh2DMeshGeneratorPeriodic ( GFace *gf ) ...@@ -1196,6 +1200,9 @@ bool gmsh2DMeshGeneratorPeriodic ( GFace *gf )
free (doc.points); free (doc.points);
free (doc.delaunay); free (doc.delaunay);
char name[245];
sprintf(name,"param%d.pos",gf->tag());
outputScalarField(m->triangles, name,1);
// Recover the boundary edges // Recover the boundary edges
...@@ -1332,7 +1339,12 @@ bool gmsh2DMeshGeneratorPeriodic ( GFace *gf ) ...@@ -1332,7 +1339,12 @@ bool gmsh2DMeshGeneratorPeriodic ( GFace *gf )
MVertex *v2 = recover_map[n[1]]; MVertex *v2 = recover_map[n[1]];
MVertex *v3 = recover_map[n[2]]; MVertex *v3 = recover_map[n[2]];
if (!n[3]) if (!n[3])
{
// 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) ); gf->triangles.push_back(new MTriangle (v1,v2,v3) );
}
else else
{ {
MVertex *v4 = recover_map[n[3]]; MVertex *v4 = recover_map[n[3]];
...@@ -1345,9 +1357,6 @@ bool gmsh2DMeshGeneratorPeriodic ( GFace *gf ) ...@@ -1345,9 +1357,6 @@ bool gmsh2DMeshGeneratorPeriodic ( GFace *gf )
// delete the mesh // 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()); // sprintf(name,"real%d.pos",gf->tag());
//outputScalarField(m->triangles, name,0); //outputScalarField(m->triangles, name,0);
......
Mesh/meshGRegion.cpp
+ 34
− 11
View file @ 22cbc692
// $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 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
// //
...@@ -111,7 +111,9 @@ void buildTetgenStructure(GRegion *gr, tetgenio &in, std::vector<MVertex*> &numb ...@@ -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) std::vector<MVertex*> &numberedV)
{ {
for(int i = numberedV.size(); i < out.numberofpoints; i++){ for(int i = numberedV.size(); i < out.numberofpoints; i++){
...@@ -439,10 +441,13 @@ void meshGRegion::operator() (GRegion *gr) ...@@ -439,10 +441,13 @@ void meshGRegion::operator() (GRegion *gr)
Msg(STATUS2, "Meshing volume %d", gr->tag()); Msg(STATUS2, "Meshing volume %d", gr->tag());
// destroy the mesh if it exists // destroy the mesh if it exists
if(gr->meshAttributes.Method == TRANSFINI)
{
deMeshGRegion dem; deMeshGRegion dem;
dem(gr); dem(gr);
MeshTransfiniteVolume(gr);
if(MeshTransfiniteVolume(gr)) return; return;
}
std::list<GFace*> faces = gr->faces(); std::list<GFace*> faces = gr->faces();
...@@ -458,9 +463,15 @@ void meshGRegion::operator() (GRegion *gr) ...@@ -458,9 +463,15 @@ void meshGRegion::operator() (GRegion *gr)
#if !defined(HAVE_TETGEN) #if !defined(HAVE_TETGEN)
Msg(GERROR, "Tetgen is not compiled in this version of Gmsh"); Msg(GERROR, "Tetgen is not compiled in this version of Gmsh");
#else #else
// put all the faces in the same model // delete the mesh for all regions
GModel::riter rit = gr->model()->firstRegion() ; GModel::riter rit = gr->model()->firstRegion() ;
if (gr != *rit)return; 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; std::list<GFace*> allFaces;
GModel::fiter fit = gr->model()->firstFace() ; GModel::fiter fit = gr->model()->firstFace() ;
while (fit != gr->model()->lastFace()){ while (fit != gr->model()->lastFace()){
...@@ -476,14 +487,24 @@ void meshGRegion::operator() (GRegion *gr) ...@@ -476,14 +487,24 @@ void meshGRegion::operator() (GRegion *gr)
sprintf(opts, "pe%c", (CTX.verbosity < 3) ? 'Q': (CTX.verbosity > 6)? 'V': '\0'); sprintf(opts, "pe%c", (CTX.verbosity < 3) ? 'Q': (CTX.verbosity > 6)? 'V': '\0');
tetrahedralize(opts, &in, &out); tetrahedralize(opts, &in, &out);
TransferTetgenMesh(gr, in, out, numberedV); 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 // restore the initial set of faces
// region--we need to recategorize them into each separate region gr->set(faces);
// now do insertion of points // now do insertion of points
insertVerticesInRegion(gr); insertVerticesInRegion(gr);
// restore the initial set of faces
gr->set(faces);
// meshNormalsPointOutOfTheRegion(gr); // meshNormalsPointOutOfTheRegion(gr);
#endif #endif
} }
...@@ -492,6 +513,8 @@ void meshGRegion::operator() (GRegion *gr) ...@@ -492,6 +513,8 @@ void meshGRegion::operator() (GRegion *gr)
#if !defined(HAVE_NETGEN) #if !defined(HAVE_NETGEN)
Msg(GERROR, "Netgen is not compiled in this version of Gmsh"); Msg(GERROR, "Netgen is not compiled in this version of Gmsh");
#else #else
deMeshGRegion dem;
dem(gr);
// orient the triangles of with respect to this region // orient the triangles of with respect to this region
meshNormalsPointOutOfTheRegion(gr); meshNormalsPointOutOfTheRegion(gr);
std::vector<MVertex*> numberedV; std::vector<MVertex*> numberedV;
......
Mesh/meshGRegionDelaunayInsertion.cpp
+ 103
− 10
View file @ 22cbc692
// $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 // Copyright (C) 1997-2007 C. Geuzaine, J.-F. Remacle
// //
...@@ -21,6 +21,7 @@ ...@@ -21,6 +21,7 @@
#include "BackgroundMesh.h" #include "BackgroundMesh.h"
#include "meshGRegionDelaunayInsertion.h" #include "meshGRegionDelaunayInsertion.h"
#include "GModel.h"
#include "GRegion.h" #include "GRegion.h"
#include "Numeric.h" #include "Numeric.h"
#include "Message.h" #include "Message.h"
...@@ -38,7 +39,7 @@ int MTet4::inCircumSphere ( const double *p ) const ...@@ -38,7 +39,7 @@ int MTet4::inCircumSphere ( const double *p ) const
return (result > 0) ? 1 : 0; 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 struct faceXtet
{ {
...@@ -119,7 +120,7 @@ void recurFindCavity ( std::list<faceXtet> & shell, ...@@ -119,7 +120,7 @@ void recurFindCavity ( std::list<faceXtet> & shell,
else if (!neigh->isDeleted() ) else if (!neigh->isDeleted() )
{ {
int circ = neigh->inCircumSphere ( v ); int circ = neigh->inCircumSphere ( v );
if (circ) if (circ && (neigh->onWhat() == t->onWhat()))
recurFindCavity ( shell, cavity,v , neigh); recurFindCavity ( shell, cavity,v , neigh);
else else
shell.push_back ( faceXtet ( t, i ) ); shell.push_back ( faceXtet ( t, i ) );
...@@ -203,7 +204,7 @@ bool insertVertex (MVertex *v , ...@@ -203,7 +204,7 @@ bool insertVertex (MVertex *v ,
while (it != shell.end()) while (it != shell.end())
{ {
MTetrahedron *t = new MTetrahedron ( it->v[0], MTetrahedron *tr = new MTetrahedron ( it->v[0],
it->v[1], it->v[1],
it->v[2], it->v[2],
v); v);
...@@ -219,7 +220,8 @@ bool insertVertex (MVertex *v , ...@@ -219,7 +220,8 @@ bool insertVertex (MVertex *v ,
// it->v[2]->y(), // it->v[2]->y(),
// it->v[2]->z()); // it->v[2]->z());
MTet4 *t4 = new MTet4 ( t , vSizes); MTet4 *t4 = new MTet4 ( tr , vSizes);
t4->setOnWhat(t->onWhat());
newTets[k++]=t4; newTets[k++]=t4;
// all new tets are pushed front in order to // all new tets are pushed front in order to
// ba able to destroy them if the cavity is not // ba able to destroy them if the cavity is not
...@@ -297,6 +299,80 @@ static void setLcs ( MTetrahedron *t, std::map<MVertex*,double> &vSizes) ...@@ -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) void insertVerticesInRegion (GRegion *gr)
{ {
...@@ -319,6 +395,23 @@ void insertVerticesInRegion (GRegion *gr) ...@@ -319,6 +395,23 @@ void insertVerticesInRegion (GRegion *gr)
gr->tetrahedra.clear(); gr->tetrahedra.clear();
connectTets ( allTets.begin(), allTets.end() ); 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()); Msg(DEBUG,"All %d tets were connected",allTets.size());
// here the classification should be done // here the classification should be done
...@@ -389,7 +482,7 @@ void insertVerticesInRegion (GRegion *gr) ...@@ -389,7 +482,7 @@ void insertVerticesInRegion (GRegion *gr)
} }
else else
{ {
gr->tetrahedra.push_back(worst->tet()); worst->onWhat()->tetrahedra.push_back(worst->tet());
} }
delete worst; delete worst;
allTets.erase(allTets.begin()); allTets.erase(allTets.begin());
......
Mesh/meshGRegionDelaunayInsertion.h
+ 6
− 1
View file @ 22cbc692
...@@ -32,13 +32,18 @@ class MTet4 ...@@ -32,13 +32,18 @@ class MTet4
double circum_radius; double circum_radius;
MTetrahedron *base; MTetrahedron *base;
MTet4 *neigh[4]; MTet4 *neigh[4];
GRegion *gr;
public : public :
inline GRegion * onWhat () const {return gr;}
inline void setOnWhat (GRegion *g) {gr=g;}
bool isDeleted () const {return deleted;} bool isDeleted () const {return deleted;}
void forceRadius (double r){circum_radius=r;} void forceRadius (double r){circum_radius=r;}
double getRadius ()const {return circum_radius;} 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; neigh[0] = neigh[1] = neigh[2] = neigh[3] = 0;
double center[3]; double center[3];
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment