diff --git a/Geo/GFace.h b/Geo/GFace.h
index 00aee5a100a99577dff8b6b3ceae1d0dc568a4bf..f7ae16463a37c39dadabf5217e4021d978037f18 100644
--- a/Geo/GFace.h
+++ b/Geo/GFace.h
@@ -50,7 +50,7 @@ class GFace : public GEntity
mean_plane meanPlane;
std::list<GEdge *> embedded_edges;
std::list<GVertex *> embedded_vertices;
- GFaceCompound *compound; // this model edge belongs to a compound
+ GFaceCompound *compound; // this model ede belongs to a compound
// replace edges (gor gluing) for specific modelers, we have to
// re-create internal data
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index 04465bd228b91b4fe6f6fa0300ba878eb82b7971..e5e3e5ccb54bb45794b594448e410bef1b0365bc 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -1116,6 +1116,8 @@ int GModel::removeDuplicateMeshVertices(double tolerance)
void GModel::createTopologyFromMesh()
{
+ removeDuplicateMeshVertices(CTX::instance()->geom.tolerance);
+
Msg::Info("Creating topology from mesh");
double t1 = Cpu();
@@ -1265,7 +1267,7 @@ void GModel::createTopologyFromFaces(std::vector<discreteFace*> &discFaces)
it != discFaces.end(); it++){
(*it)->findEdges(map_edges);
}
-
+
//return if no boundary edges (torus, sphere, ...)
if (map_edges.empty()) return;
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index fac8ccfa4aacc75e3be6b1e796adf51707148abb..9a81119a761cda740f80909bb5d70f29b041baf2 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -788,11 +788,16 @@ static bool meshGenerator(GFace *gf, int RECUR_ITER,
}
}
+
if (Msg::GetVerbosity() == 10){
+ GEdge *ge = new discreteEdge(gf->model(), 1000,0,0);
MElementOctree octree(gf->model());
doc.Voronoi();
doc.makePosView("voronoi.pos", gf);
- doc.printMedialAxis(octree.getInternalOctree(), "skeleton.pos", gf);
+ doc.printMedialAxis(octree.getInternalOctree(), "skeleton.pos", gf, ge);
+ //todo add corners with lines with closest point
+ ge->addPhysicalEntity(1000);
+ gf->model()->add(ge);
}
gf->triangles.clear();
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index ab61ddceb25172830b27a01cca2f306711d5052a..434ab77d296478ba9812d57c9335601aba756b24 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -261,8 +261,9 @@ void MeshDelaunayVolume(std::vector<GRegion*> ®ions)
allFacesSet.insert(f.begin(), f.end());
}
std::list<GFace*> allFaces;
- for(std::set<GFace*>::iterator it = allFacesSet.begin(); it != allFacesSet.end(); it++)
+ for(std::set<GFace*>::iterator it = allFacesSet.begin(); it != allFacesSet.end(); it++){
allFaces.push_back(*it);
+ }
gr->set(allFaces);
// mesh with tetgen, possibly changing the mesh on boundaries (leave
@@ -314,7 +315,8 @@ void MeshDelaunayVolume(std::vector<GRegion*> ®ions)
TransferTetgenMesh(gr, in, out, numberedV);
}
- // sort triangles in all model faces in order to be able to search in vectors
+
+ // 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;
diff --git a/Numeric/DivideAndConquer.cpp b/Numeric/DivideAndConquer.cpp
index 6e159f8977ab80f39daf64e9b7bec0721b551276..76dc26fd55b75bee2ea9ef41727709a2fd7aecb7 100644
--- a/Numeric/DivideAndConquer.cpp
+++ b/Numeric/DivideAndConquer.cpp
@@ -21,6 +21,7 @@
#include "robustPredicates.h"
#include "GPoint.h"
#include "GFace.h"
+#include "MLine.h"
#define Pred(x) ((x)->prev)
#define Succ(x) ((x)->next)
@@ -601,7 +602,7 @@ void DocRecord::makePosView(std::string fileName, GFace *gf)
fclose(f);
}
-void DocRecord::printMedialAxis(Octree *_octree, std::string fileName, GFace *gf)
+void DocRecord::printMedialAxis(Octree *_octree, std::string fileName, GFace *gf, GEdge *ge)
{
FILE *f = fopen(fileName.c_str(),"w");
@@ -631,6 +632,9 @@ void DocRecord::printMedialAxis(Octree *_octree, std::string fileName, GFace *gf
MElement *m1 = (MElement*)Octree_Search(P1, _octree);
MElement *m2 = (MElement*)Octree_Search(P2, _octree);
if (m1 && m2){
+ MVertex *v0 = new MVertex(p1.x(), p1.y(), p1.z());
+ MVertex *v1 = new MVertex(p2.x(), p2.y(), p2.z());
+ ge->lines.push_back(new MLine(v0, v1));
fprintf(f,"SL(%g,%g,%g,%g,%g,%g) {%g,%g};\n",
p1.x(),p1.y(),p1.z(),
p2.x(),p2.y(),p2.z(),
diff --git a/Numeric/DivideAndConquer.h b/Numeric/DivideAndConquer.h
index e59e05e5e64192f4cdf849cc39753f8898468bdd..223426b60dc0ad044ab221e65077db4e1a56fe67 100644
--- a/Numeric/DivideAndConquer.h
+++ b/Numeric/DivideAndConquer.h
@@ -97,7 +97,7 @@ class DocRecord{
int hullSize() {return _hullSize;}
int onHull(PointNumero i) { return std::binary_search(_hull, _hull+_hullSize, i); }
void makePosView(std::string, GFace *gf=NULL);
- void printMedialAxis(Octree *_octree, std::string, GFace *gf=NULL);
+ void printMedialAxis(Octree *_octree, std::string, GFace *gf=NULL, GEdge *ge=NULL);
double Lloyd (int);
void voronoiCell (PointNumero pt, std::vector<SPoint2> &pts) const;
static void registerBindings(binding *b);
diff --git a/contrib/Tetgen/tetgen.cxx b/contrib/Tetgen/tetgen.cxx
index 76fbc8f53f0a5e224aff8d671a6fd781634c7be7..32bf180d1ef92e58d10f7419244e92b32b1a0139 100644
--- a/contrib/Tetgen/tetgen.cxx
+++ b/contrib/Tetgen/tetgen.cxx
@@ -32406,21 +32406,25 @@ void tetgenmesh::outvoronoi(tetgenio* out)
vpointcount = 0;
index = 0;
REAL *vlist;
+ REAL *vlist2;
vlist = new REAL[(int)tetrahedrons->items * 3];
+ vlist2 = new REAL[(int)tetrahedrons->items ];
while (tetloop.tet != (tetrahedron *) NULL) {
// Calculate the circumcenter.
for (i = 0; i < 4; i++) {
pt[i] = (point) tetloop.tet[4 + i];
setpoint2tet(pt[i], encode(tetloop));
}
- circumsphere(pt[0], pt[1], pt[2], pt[3], ccent, NULL);
+ REAL radius;
+ circumsphere(pt[0], pt[1], pt[2], pt[3], ccent, &radius);
if (out == (tetgenio *) NULL) {
//fprintf(outfile, "%4d %16.8e %16.8e %16.8e\n", vpointcount + shift,
// ccent[0], ccent[1], ccent[2]);
- fprintf(outfile,"SP(%g,%g,%g) {%g};\n", ccent[0], ccent[1], ccent[2], (double)vpointcount + shift);//EMI
+ fprintf(outfile,"SP(%g,%g,%g) {%g};\n", ccent[0], ccent[1], ccent[2], (double)radius);//EMI
vlist[(vpointcount+shift)*3+0] = ccent[0];//EMI
vlist[(vpointcount+shift)*3+1] = ccent[1];//EMI
vlist[(vpointcount+shift)*3+2] = ccent[2];//EMI
+ vlist2[(vpointcount+shift)] = radius;//EMI
} else {
out->vpointlist[index++] = ccent[0];
out->vpointlist[index++] = ccent[1];
@@ -32525,7 +32529,7 @@ void tetgenmesh::outvoronoi(tetgenio* out)
fprintf(outfile,"SL(%g,%g,%g,%g,%g,%g) {%g,%g};\n",
vlist[(end1+shift)*3+0], vlist[(end1+shift)*3+1], vlist[(end1+shift)*3+2],
vlist[(end2+shift)*3+0], vlist[(end2+shift)*3+1], vlist[(end2+shift)*3+2],
- 1.0,1.0);
+ vlist2[(end1+shift)], vlist2[(end2+shift)]);
} else {
vedge->v2 = end2 + shift;
vedge->vnormal[0] = 0.0;
@@ -34567,7 +34571,7 @@ void tetrahedralize(tetgenbehavior *b, tetgenio *in, tetgenio *out,
clock_t tv[14];
tv[0] = clock();
-
+
m.b = b;
m.in = in;
m.macheps = exactinit();