diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 640344e6a8b3b8d56f5f6b8971210a8117af8ccb..b18f2fceed8a56a324289fc60a87a299db18caa6 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -1331,10 +1331,12 @@ static void _associateEntityWithElementVertices(GEntity *ge, std::vector<T*> &el
void GModel::_createGeometryOfDiscreteEntities(bool force)
{
+ printf("%d vertices\n",getNumVertices());
if (CTX::instance()->meshDiscrete){
createTopologyFromMeshNew ();
exportDiscreteGEOInternals();
}
+ printf("%d vertices\n",getNumVertices());
if (force || CTX::instance()->meshDiscrete){
Msg::Info("Creating the geometry of discrete curves");
@@ -1345,6 +1347,7 @@ void GModel::_createGeometryOfDiscreteEntities(bool force)
}
}
}
+ printf("%d vertices\n",getNumVertices());
if (CTX::instance()->meshDiscrete){
Msg::Info("Creating the geometry of discrete surfaces");
@@ -1355,6 +1358,7 @@ void GModel::_createGeometryOfDiscreteEntities(bool force)
}
}
}
+ printf("%d vertices\n",getNumVertices());
}
diff --git a/Geo/GModelCreateTopologyFromMesh.cpp b/Geo/GModelCreateTopologyFromMesh.cpp
index c22f363c801ebf01f153333be40bf74b1dd7c214..55ded1a21b3c2cc44baabbb08d500351690fb691 100644
--- a/Geo/GModelCreateTopologyFromMesh.cpp
+++ b/Geo/GModelCreateTopologyFromMesh.cpp
@@ -103,7 +103,7 @@ void createTopologyFromMesh1D ( GModel *gm , int &num) {
if (it->second.stuff.size() > 1) {
// it->second.print();
num++;
- discreteVertex *dv = new discreteVertex ( gm , GModel::current()->getGEOInternals()->MaxPointNum++);
+ discreteVertex *dv = new discreteVertex ( gm , ++GModel::current()->getGEOInternals()->MaxPointNum);
gm->add(dv);
MVertex *v = it->first;
MPoint *mp = new MPoint(v);
@@ -135,7 +135,23 @@ void createTopologyFromMesh1D ( GModel *gm , int &num) {
}
}
+ for(GModel::eiter it = gm->firstEdge(); it != gm->lastEdge(); it++) {
+ if (!(*it)->getBeginVertex() && !(*it)->getEndVertex()){
+ num++;
+ discreteVertex *dv = new discreteVertex ( gm , ++GModel::current()->getGEOInternals()->MaxPointNum);
+ gm->add(dv);
+ printf("creating vertex %d for periodic edge %d\n",dv->tag(),(*it)->tag());
+ MVertex *v = (*it)->lines[0]->getVertex(0);
+ MPoint *mp = new MPoint(v);
+ dv->points.push_back(mp);
+ dv->addEdge (*it);
+ (*it)->setBeginVertex(dv);
+ (*it)->setEndVertex(dv);
+ _topology[*it].insert(dv);
+ }
+ }
+
// NICE :-)
}
@@ -231,7 +247,7 @@ void ensureManifoldFaces ( GModel *gm ) {
for(unsigned int i=0;i<f.size(); i++)ensureManifoldFace (f[i]);
}
-
+/// FIXME : 2 TIMES THE SAME MLINE IN EACH CONNECTED PART IF PERIODIC
std::vector<GEdge*> ensureSimplyConnectedEdge ( GEdge *ge ) {
std::vector<GEdge*> _all;
std::set<MLine*> _lines;
@@ -381,6 +397,9 @@ void createTopologyFromMesh2D ( GModel *gm , int & num) {
std::map<GEdge*, std::vector<GEdge*> > _parts;
for(GModel::eiter it = gm->firstEdge(); it != gm->lastEdge(); it++) {
_parts[*it] = ensureSimplyConnectedEdge (*it);
+ // for (int i=0;i<(*it)->lines.size();i++)
+ // printf("%d %d\n",(*it)->lines[i]->getVertex(0)->getNum(),(*it)->lines[i]->getVertex(1)->getNum());
+ // getchar();
}
// create Face 2 Edge topology
@@ -648,13 +667,14 @@ void GModel::createTopologyFromMeshNew ( ) {
getEntities(entities);
std::set<MVertex*> vs;
for(unsigned int i = 0; i < entities.size(); i++){
+ // printf("entity %d of dim %d\n",entities[i]->tag(),entities[i]->dim());
vs.insert (entities[i]->mesh_vertices.begin(), entities[i]->mesh_vertices.end());
entities[i]->mesh_vertices.clear();
}
std::vector<MVertex*> cc;
cc.insert(cc.begin(), vs.begin(), vs.end());
_storeVerticesInEntities(cc);
-
+
// printf("%d vertices\n", getNumVertices());
double t2 = Cpu();
diff --git a/Geo/discreteDiskFace.cpp b/Geo/discreteDiskFace.cpp
index a4dfab09db998b2d609d818c99d753ced3bc0410..a1433ccc66017057254c03198d331463933982de 100644
--- a/Geo/discreteDiskFace.cpp
+++ b/Geo/discreteDiskFace.cpp
@@ -983,7 +983,7 @@ GPoint discreteDiskFace::intersectionWithCircle(const SVector3 &n1, const SVecto
m[0][1] = n.z();
m[1][0] = t.y();
m[1][1] = t.z();
- if (fabs(det2x2(m)) > 1.e-8){
+ if (fabs(det2x2(m)) > 1.e-12){
sys2x2(m,rhs,r);
x0 = SVector3(0,r[0],r[1]);
}
@@ -992,7 +992,7 @@ GPoint discreteDiskFace::intersectionWithCircle(const SVector3 &n1, const SVecto
m[0][1] = n.z();
m[1][0] = t.x();
m[1][1] = t.z();
- if (fabs(det2x2(m)) > 1.e-8){
+ if (fabs(det2x2(m)) > 1.e-12){
sys2x2(m,rhs,r);
x0 = SVector3(r[0],0,r[1]);
}
@@ -1005,7 +1005,7 @@ GPoint discreteDiskFace::intersectionWithCircle(const SVector3 &n1, const SVecto
x0 = SVector3(r[0],r[1],0);
}
else{
- printf("mauvaise pioche\n");
+ // printf("mauvaise pioche\n");
continue;
}
}
diff --git a/Geo/discreteEdge.cpp b/Geo/discreteEdge.cpp
index 40f0272def4ab9943a86403b3929611e82188c29..8af4bdf1e973bf84843effc239b16666d8682a74 100644
--- a/Geo/discreteEdge.cpp
+++ b/Geo/discreteEdge.cpp
@@ -70,6 +70,7 @@ void discreteEdge::orderMLines()
it != segments.end(); ++it){
MVertex *vL = (*it)->getVertex(0);
MVertex *vR = (*it)->getVertex(1);
+
std::map<MVertex*,MLine*>::iterator it1 = boundv.find(vL);
if (it1 == boundv.end()) boundv.insert(std::make_pair(vL,*it));
else boundv.erase(it1);
@@ -78,6 +79,7 @@ void discreteEdge::orderMLines()
else boundv.erase(it2);
}
// find the first MLine and erase it from the list segments
+
MLine *firstLine;
if (boundv.size() == 2){ // non periodic
firstLine = (boundv.begin())->second;
@@ -214,7 +216,7 @@ void discreteEdge::setBoundVertices()
v1 = bound_vertices[1];
}
else if (boundv.size() == 0){
- //printf("closed loop for discrete Edge =%d \n", this->tag());
+ printf("closed loop for discrete Edge =%d \n", this->tag());
GVertex* bound_vertex;
std::vector<MLine*>::const_iterator it = lines.begin();
MVertex* vE = (*it)->getVertex(0);
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index 53f5cb42190b45edaf6960ce7d9fe66ddb1c75ed..907f8abacd58339db56ac6cac1f014999b15eab3 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -704,6 +704,7 @@ bool insertVertexB (std::list<edgeXface> &shell,
double *metric,
MTri3 **oneNewTriangle)
{
+ if (cavity.size() == 1) return false;
if (shell.size() != cavity.size() + 2) return false;
std::list<MTri3*> new_cavity;
@@ -747,7 +748,7 @@ bool insertVertexB (std::list<edgeXface> &shell,
// avoid angles that are too obtuse
double cosv = ((d1*d1+d2*d2-d3*d3)/(2.*d1*d2));
- if ((d1 < LL * .35 || d2 < LL * .35 || cosv < -.9999) && !force) {
+ if ((d1 < LL * .55 || d2 < LL * .55 || cosv < -.9999) && !force) {
onePointIsTooClose = true;
// printf("%12.5E %12.5E %12.5E %12.5E \n",d1,d2,LL,cosv);
}
@@ -827,7 +828,7 @@ bool insertVertex(bool force, GFace *gf, MVertex *v, double *param , MTri3 *t,
else{
recurFindCavityAniso(gf, shell, cavity, metric, param, t, data);
}
-
+
return insertVertexB(shell, cavity, force, gf, v, param , t,
allTets,
activeTets,
@@ -1058,7 +1059,7 @@ void bowyerWatson(GFace *gf, int MAXPNT,
// if (equivalence)_printTris ("before.pos", AllTris.begin(), AllTris.end(), DATA);
int nbSwaps = edgeSwapPass(gf, AllTris, SWCR_DEL, DATA);
// _printTris ("after2.pos", AllTris, Us,Vs);
- Msg::Debug("Delaunization of the initial mesh done (%d swaps)", nbSwaps);
+ // Msg::Debug("Delaunization of the initial mesh done (%d swaps)", nbSwaps);
if(AllTris.empty()){
Msg::Error("No triangles in initial mesh");
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index a6b1963883f8859b595cf8d2b04e84f21e417010..c1e34275a448c496a8e6ec3ce2bf6fdf28e491fe 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -818,7 +818,7 @@ void _relocateVertex(GFace *gf, MVertex *ver,
}
void getAllBoundaryLayerVertices (GFace *gf, std::set<MVertex*> &vs){
- //return;
+ // return;
vs.clear();
BoundaryLayerColumns* _columns = gf->getColumns();
if (!_columns)return;
@@ -836,11 +836,6 @@ void getAllBoundaryLayerVertices (GFace *gf, std::set<MVertex*> &vs){
void laplaceSmoothing(GFace *gf, int niter, bool infinity_norm)
{
if (!niter)return;
- RelocateVertices (gf, niter);
- return;
-
-
- printf("argh\n");
std::set<MVertex*> vs;
getAllBoundaryLayerVertices (gf, vs);
v2t_cont adj;
@@ -877,12 +872,6 @@ bool edgeSwap(std::set<swapquad> &configs, MTri3 *t1, GFace *gf, int iLocalEdge,
MTri3 *t2 = t1->getNeigh(iLocalEdge);
if(!t2) return false;
- double triQualityRef = 0;
- if (cr == SWCR_QUAL){
- triQualityRef=std::min(qmTriangle::gamma(t1->tri()), qmTriangle::gamma(t2->tri()));
- if (triQualityRef > 0.5)return false;
- }
-
MVertex *v1 = t1->tri()->getVertex(iLocalEdge == 0 ? 2 : iLocalEdge - 1);
MVertex *v2 = t1->tri()->getVertex((iLocalEdge) % 3);
MVertex *v3 = t1->tri()->getVertex((iLocalEdge + 1) % 3);
@@ -893,14 +882,12 @@ bool edgeSwap(std::set<swapquad> &configs, MTri3 *t1, GFace *gf, int iLocalEdge,
for(int i = 0; i < 3; i++)
- if(t2->tri()->getVertex(i) != v1 && t2->tri()->getVertex(i) != v2){
+ if(t2->tri()->getVertex(i) != v1 && t2->tri()->getVertex(i) != v2)
v4 = t2->tri()->getVertex(i);
- break;
- }
- // swapquad sq (v1, v2, v3, v4);
- // if(configs.find(sq) != configs.end()) return false;
- // configs.insert(sq);
+ swapquad sq (v1, v2, v3, v4);
+ if(configs.find(sq) != configs.end()) return false;
+ configs.insert(sq);
if (edgeSwapDelProj(v3,v4,v2,v1))return false;
@@ -908,15 +895,15 @@ bool edgeSwap(std::set<swapquad> &configs, MTri3 *t1, GFace *gf, int iLocalEdge,
MTriangle *t1b = new MTriangle(v2, v3, v4);
MTriangle *t2b = new MTriangle(v4, v3, v1);
-
switch(cr){
case SWCR_QUAL:
{
+ const double triQualityRef = std::min(qmTriangle::gamma(t1->tri()),
+ qmTriangle::gamma(t2->tri()));
const double triQuality = std::min(qmTriangle::gamma(t1b),
qmTriangle::gamma(t2b));
-
if (!edgeSwapDelProj(v1,v2,v3,v4)){
- if(triQuality < triQualityRef/* && triQuality < .001*/){
+ if(triQuality < triQualityRef){
delete t1b;
delete t2b;
return false;
@@ -1006,13 +993,47 @@ bool edgeSwap(std::set<swapquad> &configs, MTri3 *t1, GFace *gf, int iLocalEdge,
}
int edgeSwapPass(GFace *gf, std::set<MTri3*, compareTri3Ptr> &allTris,
+ const swapCriterion &cr,bidimMeshData & data)
+{
+ typedef std::set<MTri3*, compareTri3Ptr> CONTAINER;
+
+ int nbSwapTot = 0;
+ std::set<swapquad> configs;
+ for(int iter = 0; iter < 1200; iter++){
+ int nbSwap = 0;
+ std::vector<MTri3*> newTris;
+ for(CONTAINER::iterator it = allTris.begin(); it != allTris.end(); ++it){
+ if(!(*it)->isDeleted()){
+ for(int i = 0; i < 3; i++){
+ if(edgeSwap(configs, *it, gf, i, newTris, cr, data)){
+ nbSwap++;
+ break;
+ }
+ }
+ }
+ else{
+ delete *it;
+ CONTAINER::iterator itb = it;
+ ++it;
+ allTris.erase(itb);
+ if(it == allTris.end()) break;
+ }
+ }
+ allTris.insert(newTris.begin(), newTris.end());
+ nbSwapTot += nbSwap;
+ if(nbSwap == 0) break;
+ }
+ return nbSwapTot;
+}
+
+int edgeSwapPass2(GFace *gf, std::set<MTri3*, compareTri3Ptr> &allTris,
const swapCriterion &cr,bidimMeshData & data)
{
typedef std::set<MTri3*, compareTri3Ptr> CONTAINER;
int nbSwapTot = 0;
std::set<swapquad> configs;
- for(int iter = 0; iter < 3; iter++){
+ for(int iter = 0; iter < 1200; iter++){
int nbSwap = 0;
std::vector<MTri3*> newTris;
for(CONTAINER::iterator it = allTris.begin(); it != allTris.end(); ++it){
@@ -1024,12 +1045,17 @@ int edgeSwapPass(GFace *gf, std::set<MTri3*, compareTri3Ptr> &allTris,
}
}
}
+ else{
+ delete *it;
+ CONTAINER::iterator itb = it;
+ ++it;
+ allTris.erase(itb);
+ if(it == allTris.end()) break;
+ }
}
-
allTris.insert(newTris.begin(), newTris.end());
nbSwapTot += nbSwap;
if(nbSwap == 0) break;
- // printf("%d swaps %d new tris\n",nbSwap, newTris.size());
}
return nbSwapTot;
}