diff --git a/Geo/GModelIO_OCC.cpp b/Geo/GModelIO_OCC.cpp
index b2b73bd7afc14e5a0cfa24862d44f2c2d462cec9..08d0deaa980e3536fe16482c4e31f7f1d56c8eee 100644
--- a/Geo/GModelIO_OCC.cpp
+++ b/Geo/GModelIO_OCC.cpp
@@ -678,7 +678,6 @@ void OCC_Internals::applyBooleanOperator(TopoDS_Shape tool, const BooleanOperato
break;
case OCC_Internals::Fuse :
{
- printf("coucou\n");
BRepAlgoAPI_Fuse BO (tool, shape);
if (!BO.IsDone()) {
Msg::Error("Fuse operation can not be performed on the given shapes");
diff --git a/Mesh/Generator.cpp b/Mesh/Generator.cpp
index da4b10461f375df7335a9825abbc5c0061ab6bdd..bea223c2ff82004cd7aa4caecaa9004d0beb2c7f 100644
--- a/Mesh/Generator.cpp
+++ b/Mesh/Generator.cpp
@@ -509,7 +509,7 @@ static void Mesh3D(GModel *m)
// global operation, which can require changing the surface mesh!)
SubdivideExtrudedMesh(m);
- // then mesh all the non-delaunay regions
+ // then mesh all the non-delaunay regions (front3D with netgen)
std::vector<GRegion*> delaunay;
std::for_each(m->firstRegion(), m->lastRegion(), meshGRegion(delaunay));
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index 6501751ca361c1c62fbe359e41620d051696e0ec..12f4ea05c9f6b4e4237dd7af113ff44cc650e9b1 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -5,6 +5,7 @@
#include <sstream>
#include <stdlib.h>
+#include <map>
#include "meshGFace.h"
#include "meshGFaceBDS.h"
#include "meshGFaceDelaunayInsertion.h"
@@ -389,6 +390,7 @@ static bool meshGenerator(GFace *gf, int RECUR_ITER,
std::map<BDS_Point*, MVertex*> recoverMap;
std::map<MVertex*, BDS_Point*> recoverMapInv;
std::list<GEdge*> edges = gf->edges();
+ std::list<int> dir = gf->edgeOrientations();
// replace edges by their compounds
// if necessary split compound and remesh parts
@@ -397,36 +399,99 @@ static bool meshGenerator(GFace *gf, int RECUR_ITER,
isMeshed = checkMeshCompound((GFaceCompound*) gf, edges);
if (isMeshed) return true;
}
-
+
+ // build face normal
+ std::list<GEdge*>::iterator ite = edges.begin();
+ std::list<int>::iterator itd = dir.begin();
+ // SVector3 tan1, tan2;
+ // for (int i=0; i<2; i++){
+ // MVertex *v1 = (*ite)->lines[0]->getVertex(0);
+ // MVertex *v2 = (*ite)->lines[0]->getVertex(1);
+ // SPoint2 p1, p2;
+ // if (*itd > 0.0){
+ // //printf("v1=%g %g v2=%g %g \n", v1->x(), v1->y(), v2->x(), v2->y());
+ // reparamMeshVertexOnFace(v1, gf, p1);
+ // reparamMeshVertexOnFace(v2, gf, p2);
+ // //printf("orient OK p1=%g %g p2=%g %g \n", p1.x(), p1.y(), p2.x(), p2.y());
+ // }
+ // else{
+ // //printf("v1=%g %g v2=%g %g \n", v1->x(), v1->y(), v2->x(), v2->y());
+ // reparamMeshVertexOnFace(v2, gf, p1);
+ // reparamMeshVertexOnFace(v1, gf, p2);
+ // //printf("orient -1 p1=%g %g p2=%g %g \n", p1.x(), p1.y(), p2.x(), p2.y());
+ // }
+ // //printf("edge=%d tan=%g %g \n", (*ite)->tag(), p2.x()-p1.x(),p2.y()-p1.y());
+ // if (i==0) tan1 = SVector3(p2.x()-p1.x(),p2.y()-p1.y(),0.0);
+ // else tan2 = SVector3(p2.x()-p1.x(),p2.y()-p1.y(),0.0);
+ // ite++; itd++;
+ // }
+ // tan1.normalize(); tan2.normalize();
+ // SVector3 next = crossprod(tan1,tan2);
+ // next.normalize();
+ // printf("next =%g %g %g \n", next.x(), next.y(), next.z());
+ SVector3 next(0., 0., 1.0);
+
// build a set with all points of the boundaries
std::set<MVertex*> all_vertices;
+ std::map<SPoint2, SVector3> pt2Normal;
+ ite = edges.begin();
+ itd = dir.begin();
+ while(ite != edges.end()){
+ if((*ite)->isSeam(gf)) return false;
+ if(!(*ite)->isMeshDegenerated()){
+ for(unsigned int i = 0; i< (*ite)->lines.size(); i++){
+ MVertex *v1 = (*ite)->lines[i]->getVertex(0);
+ MVertex *v2 = (*ite)->lines[i]->getVertex(1);
+ all_vertices.insert(v1);
+ all_vertices.insert(v2);
+
+ SPoint2 p1, p2;
+ if (*itd > 0.0){
+ reparamMeshVertexOnFace(v1, gf, p1);
+ reparamMeshVertexOnFace(v2, gf, p2);
+ }
+ else{
+ reparamMeshVertexOnFace(v2, gf, p1);
+ reparamMeshVertexOnFace(v1, gf, p2);
+ }
+ SVector3 tan(p2.x()-p1.x(),p2.y()-p1.y(),0.0);
+ tan.normalize();
+ SVector3 ne = crossprod(tan, next);
+ ne.normalize();
+ std::map<SPoint2, SVector3>::iterator it = pt2Normal.find(p1);
+ if (it == pt2Normal.end())
+ pt2Normal.insert(std::make_pair(p1,ne));
+ else{
+ SVector3 n1 = it->second;
+ it->second = n1+ne;
+ }
+ std::map<SPoint2, SVector3>::iterator it2 = pt2Normal.find(p2);
+ if (it2 == pt2Normal.end())
+ pt2Normal.insert(std::make_pair(p2,ne));
+ else{
+ SVector3 n1 = it2->second;
+ it2->second = n1+ne;
+ }
- std::list<GEdge*>::iterator it = edges.begin();
- while(it != edges.end()){
- if((*it)->isSeam(gf)) return false;
- if(!(*it)->isMeshDegenerated()){
- for(unsigned int i = 0; i< (*it)->lines.size(); i++){
- all_vertices.insert((*it)->lines[i]->getVertex(0));
- all_vertices.insert((*it)->lines[i]->getVertex(1));
}
}
else
- printf("edge %d degenerated mesh \n", (*it)->tag());
- ++it;
+ printf("edge %d degenerated mesh \n", (*ite)->tag());
+ ++ite; ++itd;
}
std::list<GEdge*> emb_edges = gf->embeddedEdges();
- it = emb_edges.begin();
- while(it != emb_edges.end()){
- if(!(*it)->isMeshDegenerated()){
- all_vertices.insert((*it)->mesh_vertices.begin(),
- (*it)->mesh_vertices.end() );
- all_vertices.insert((*it)->getBeginVertex()->mesh_vertices.begin(),
- (*it)->getBeginVertex()->mesh_vertices.end());
- all_vertices.insert((*it)->getEndVertex()->mesh_vertices.begin(),
- (*it)->getEndVertex()->mesh_vertices.end());
+ ite = emb_edges.begin();
+ while(ite != emb_edges.end()){
+ if(!(*ite)->isMeshDegenerated()){
+ all_vertices.insert((*ite)->mesh_vertices.begin(),
+ (*ite)->mesh_vertices.end() );
+ all_vertices.insert((*ite)->getBeginVertex()->mesh_vertices.begin(),
+ (*ite)->getBeginVertex()->mesh_vertices.end());
+ all_vertices.insert((*ite)->getEndVertex()->mesh_vertices.begin(),
+ (*ite)->getEndVertex()->mesh_vertices.end());
}
- ++it;
+ ++ite;
}
// add embedded vertices
@@ -485,6 +550,7 @@ static bool meshGenerator(GFace *gf, int RECUR_ITER,
// use a divide & conquer type algorithm to create a triangulation.
// We add to the triangulation a box with 4 points that encloses the
// domain.
+ std::map<SPoint2, SVector3> pt2NormalNEW;
{
DocRecord doc(points.size() + 4);
for(unsigned int i = 0; i < points.size(); i++){
@@ -496,6 +562,15 @@ static bool meshGenerator(GFace *gf, int RECUR_ITER,
doc.points[i].where.v = points[i]->v + YY;
doc.points[i].data = points[i];
doc.points[i].adjacent = NULL;
+
+ SPoint2 p(points[i]->u, points[i]->v);
+ SPoint2 pNEW(points[i]->u+XX, points[i]->v+YY);
+ std::map<SPoint2, SVector3>::iterator it = pt2Normal.find(p);
+ if (it != pt2Normal.end() ){
+ pt2NormalNEW.insert(std::make_pair(pNEW,it->second));
+ }
+ else
+ Msg::Error("no point found \n");
}
// increase the size of the bounding box
@@ -548,31 +623,37 @@ static bool meshGenerator(GFace *gf, int RECUR_ITER,
Msg::Debug("Recovering %d model Edges", edges.size());
std::set<EdgeToRecover> edgesToRecover;
std::set<EdgeToRecover> edgesNotRecovered;
- it = edges.begin();
- while(it != edges.end()){
- if(!(*it)->isMeshDegenerated())
+ ite = edges.begin();
+ while(ite != edges.end()){
+ if(!(*ite)->isMeshDegenerated())
recoverEdge
- (m, *it, recoverMapInv, &edgesToRecover, &edgesNotRecovered, 1);
- ++it;
+ (m, *ite, recoverMapInv, &edgesToRecover, &edgesNotRecovered, 1);
+ ++ite;
}
- it = emb_edges.begin();
- while(it != emb_edges.end()){
- if(!(*it)->isMeshDegenerated())
+ ite = emb_edges.begin();
+ while(ite != emb_edges.end()){
+ if(!(*ite)->isMeshDegenerated())
recoverEdge
- (m, *it, recoverMapInv, &edgesToRecover, &edgesNotRecovered, 1);
- ++it;
+ (m, *ite, recoverMapInv, &edgesToRecover, &edgesNotRecovered, 1);
+ ++ite;
}
Msg::Debug("Recovering %d mesh Edges", edgesToRecover.size());
+
+ if (Msg::GetVerbosity() == 10){
+ doc.Voronoi();
+ doc.makePosView("voronoi.pos", gf);
+ doc.printMedialAxis(pt2NormalNEW, "skeleton.pos", gf);
+ }
// effectively recover the medge
- it = edges.begin();
- while(it != edges.end()){
- if(!(*it)->isMeshDegenerated()){
+ ite = edges.begin();
+ while(ite != edges.end()){
+ if(!(*ite)->isMeshDegenerated()){
recoverEdge
- (m, *it, recoverMapInv, &edgesToRecover, &edgesNotRecovered, 2);
+ (m, *ite, recoverMapInv, &edgesToRecover, &edgesNotRecovered, 2);
}
- ++it;
+ ++ite;
}
if(edgesNotRecovered.size()){
@@ -649,12 +730,12 @@ static bool meshGenerator(GFace *gf, int RECUR_ITER,
}
}
- it = emb_edges.begin();
- while(it != emb_edges.end()){
- if(!(*it)->isMeshDegenerated())
+ ite = emb_edges.begin();
+ while(ite != emb_edges.end()){
+ if(!(*ite)->isMeshDegenerated())
recoverEdge
- (m, *it, recoverMapInv, &edgesToRecover, &edgesNotRecovered, 2);
- ++it;
+ (m, *ite, recoverMapInv, &edgesToRecover, &edgesNotRecovered, 2);
+ ++ite;
}
// compute characteristic lengths at vertices
diff --git a/Mesh/meshGFaceLloyd.cpp b/Mesh/meshGFaceLloyd.cpp
index 2015481ca59bd12a5b3de3710fb7966eaf47a06d..8fc7b948729a4a282da1b85773efb6d67ddceb74 100644
--- a/Mesh/meshGFaceLloyd.cpp
+++ b/Mesh/meshGFaceLloyd.cpp
@@ -14,11 +14,13 @@ void lloydAlgorithm::operator () (GFace *gf)
std::set<MVertex*> all;
// get all the points of the face ...
-
for (unsigned int i = 0; i < gf->getNumMeshElements(); i++){
MElement *e = gf->getMeshElement(i);
for (int j = 0;j<e->getNumVertices(); j++){
- all.insert(e->getVertex(j));
+ MVertex *v = e->getVertex(j);
+ //if (v->onWhat()->dim() < 2){
+ all.insert(v);
+ //}
}
}
@@ -38,7 +40,6 @@ void lloydAlgorithm::operator () (GFace *gf)
// gf->getNumMeshElements(), (int)all.size(), LC2D);
int i = 0;
-
for (std::set<MVertex*>::iterator it = all.begin(); it != all.end(); ++it){
SPoint2 p;
bool success = reparamMeshVertexOnFace(*it, gf, p);
@@ -55,11 +56,12 @@ void lloydAlgorithm::operator () (GFace *gf)
triangulator.y(i) = p.y() + YY;
triangulator.data(i++) = (*it);
}
-
+
// compute the Voronoi diagram
triangulator.Voronoi();
//printf("hullSize = %d\n",triangulator.hullSize());
triangulator.makePosView("LloydInit.pos");
+ //triangulator.printMedialAxis("medialAxis.pos");
// now do the Lloyd iterations
int ITER = 0;
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 2bc09056ceaba5a38e8fdd836b8510d8c33517c2..7cfb9069bb12c61682a9f45741b3aa43b6618e2c 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -272,8 +272,10 @@ void MeshDelaunayVolume(std::vector<GRegion*> ®ions)
std::vector<MVertex*> numberedV;
char opts[128];
buildTetgenStructure(gr, in, numberedV);
- sprintf(opts, "pe%c",
- (Msg::GetVerbosity() < 3) ? 'Q': (Msg::GetVerbosity() > 6)? 'V': '\0');
+ if (Msg::GetVerbosity() == 10)
+ sprintf(opts, "peVv");
+ else
+ sprintf(opts, "pe%c", (Msg::GetVerbosity() < 3) ? 'Q': (Msg::GetVerbosity() > 6)? 'V': '\0');
try{
tetrahedralize(opts, &in, &out);
}
@@ -325,6 +327,7 @@ void MeshDelaunayVolume(std::vector<GRegion*> ®ions)
// now do insertion of points
insertVerticesInRegion(gr);
+
#endif
}
diff --git a/Numeric/DivideAndConquer.cpp b/Numeric/DivideAndConquer.cpp
index 3375e2102735c9fd26a40eab7729497d5fb5507c..b313a24d5afe742ca249aa9f6a4fe9a11977eeba 100644
--- a/Numeric/DivideAndConquer.cpp
+++ b/Numeric/DivideAndConquer.cpp
@@ -19,6 +19,8 @@
#include "DivideAndConquer.h"
#include "Numeric.h"
#include "robustPredicates.h"
+#include "GPoint.h"
+#include "GFace.h"
#define Pred(x) ((x)->prev)
#define Succ(x) ((x)->next)
@@ -564,7 +566,7 @@ void DocRecord::voronoiCell(PointNumero pt, std::vector<SPoint2> &pts) const
*/
-void DocRecord::makePosView(std::string fileName)
+void DocRecord::makePosView(std::string fileName, GFace *gf)
{
FILE *f = fopen(fileName.c_str(),"w");
if (_adjacencies){
@@ -573,15 +575,24 @@ void DocRecord::makePosView(std::string fileName)
std::vector<SPoint2> pts;
double pc[2] = {(double)points[i].where.h, (double)points[i].where.v};
if (!onHull(i)){
- fprintf(f,"SP(%g,%g,%g) {%g};\n",pc[0],pc[1],0.0,(double)i);
+ GPoint p0(pc[0], pc[1], 0.0);
+ if (gf) p0 = gf->point(pc[0], pc[1]);
+ fprintf(f,"SP(%g,%g,%g) {%g};\n",p0.x(),p0.y(),p0.z(),(double)i);
voronoiCell (i,pts);
for (unsigned int j = 0; j < pts.size(); j++){
+ SPoint2 pp1 = pts[j];
+ SPoint2 pp2 = pts[(j+1)%pts.size()];
+ GPoint p1(pp1.x(), pp1.y(), 0.0);
+ GPoint p2(pp2.x(), pp2.y(), 0.0);
+ if (gf) {
+ p1 = gf->point(p1.x(), p1.y());
+ p2 = gf->point(p2.x(), p2.y());
+ }
fprintf(f,"SL(%g,%g,%g,%g,%g,%g) {%g,%g};\n",
- pts[j].x(),pts[j].y(),0.0,
- pts[(j+1)%pts.size()].x(),pts[(j+1)%pts.size()].y(),0.0,
+ p1.x(),p1.y(),p1.z(),p2.x(),p2.y(),p2.z(),
(double)i,(double)i);
}
- }
+ }
else {
fprintf(f,"SP(%g,%g,%g) {%g};\n",pc[0],pc[1],0.0,-(double)i);
}
@@ -591,6 +602,51 @@ void DocRecord::makePosView(std::string fileName)
fclose(f);
}
+void DocRecord::printMedialAxis(std::map<SPoint2, SVector3> &pt2Normal,
+ std::string fileName, GFace *gf)
+{
+
+ FILE *f = fopen(fileName.c_str(),"w");
+ if (_adjacencies){
+ fprintf(f,"View \"medial axis\" {\n");
+ for(PointNumero i = 0; i < numPoints; i++) {
+ std::vector<SPoint2> pts;
+ SPoint2 pc((double)points[i].where.h, (double)points[i].where.v);
+ if (!onHull(i)){
+ std::map<SPoint2, SVector3>::const_iterator it = pt2Normal.find(pc);
+ if (it == pt2Normal.end()) printf("pt not found \n");
+ SVector3 n = it->second;
+ //fprintf(f,"VP(%g,%g,%g) {%g,%g,%g};\n",pc.x(),pc.y(), 0.0, n.x(), n.y(), n.z());
+ GPoint p0(pc[0], pc[1], 0.0);
+ if (gf) p0 = gf->point(pc[0], pc[1]);
+ fprintf(f,"SP(%g,%g,%g) {%g};\n",p0.x(),p0.y(),p0.z(),(double)i);
+ voronoiCell (i,pts);
+ for (unsigned int j = 0; j < pts.size(); j++){
+ SVector3 pp1(pts[j].x(), pts[j].y(), 0.0);
+ SVector3 pp2(pts[(j+1)%pts.size()].x(),pts[(j+1)%pts.size()].y(), 0.0);
+ SVector3 v1(pp1.x()-pc.x(),pp1.y()-pc.y(),0.0);
+ SVector3 v2(pp2.x()-pc.x(),pp2.y()-pc.y(),0.0);
+ GPoint p1(pp1.x(), pp1.y(), 0.0);
+ GPoint p2(pp2.x(), pp2.y(), 0.0);
+ if (gf) {
+ p1 = gf->point(p1.x(), p1.y());
+ p2 = gf->point(p2.x(), p2.y());
+ }
+ if (dot(v1,n) > 0.0 && dot(v2,n) > 0.0){
+ fprintf(f,"SL(%g,%g,%g,%g,%g,%g) {%g,%g};\n",
+ p1.x(),p1.y(),p1.z(),
+ p2.x(),p2.y(),p2.z(),
+ (double)i,(double)i);
+ }
+ }
+ }
+ }
+ fprintf(f,"};\n");
+ }
+ fclose(f);
+
+}
+
void centroidOfOrientedBox(std::vector<SPoint2> &pts, const double &angle,
double &xc, double &yc, double &inertia, double &area)
{
@@ -821,8 +877,8 @@ void DocRecord::registerBindings(binding *b)
cm = cb->addMethod("hullSize", &DocRecord::hullSize);
cm->setDescription("returns the size of the hull");
cm = cb->addMethod("makePosView", &DocRecord::makePosView);
- cm->setDescription("save a .pos file with the voronoi");
- cm->setArgNames("FileName",NULL);
+ cm->setDescription("save a .pos file with the voronoi in GFace");
+ cm->setArgNames("FileName", "GFace", NULL);
cm = cb->addMethod("Lloyd", &DocRecord::Lloyd);
cm->setDescription("do one iteration of Lloyd's algorithm");
cm->setArgNames("Type",NULL);
diff --git a/Numeric/DivideAndConquer.h b/Numeric/DivideAndConquer.h
index aa75d2b79e7818500b89e99e6e758176122ba00a..87eeee10fab48a1d5ae1295a3641d6f259fd6dd0 100644
--- a/Numeric/DivideAndConquer.h
+++ b/Numeric/DivideAndConquer.h
@@ -94,7 +94,8 @@ class DocRecord{
void Voronoi ();
int hullSize() {return _hullSize;}
int onHull(PointNumero i) { return std::binary_search(_hull, _hull+_hullSize, i); }
- void makePosView(std::string);
+ void makePosView(std::string, GFace *gf=NULL);
+ void printMedialAxis(std::map<SPoint2, SVector3> &pt2Normal, std::string, GFace *gf=NULL);
double Lloyd (int);
void voronoiCell (PointNumero pt, std::vector<SPoint2> &pts) const;
static void registerBindings(binding *b);
diff --git a/Plugin/Distance.cpp b/Plugin/Distance.cpp
index fb7199ac2f39a14ebe69374f7ae03fbd1615ba26..32c14b2147d79abe86bfe2cd27ecdf42a5dc84df 100644
--- a/Plugin/Distance.cpp
+++ b/Plugin/Distance.cpp
@@ -359,9 +359,12 @@ PView *GMSH_DistancePlugin::execute(PView *v)
//compute also orthogonal vector to distance field
// A Uortho = -C DIST
//------------------------------------------------
- if (ortho > 0){
+ if (ortho > 0 && _maxDim != 2){
+ Msg::Error("The orthogonal field is only implemented for a 2D distance field !");
+ }
+ else if (ortho > 0 && _maxDim == 2){
#if defined(HAVE_SOLVER)
-
+
#ifdef HAVE_TAUCS
linearSystemCSRTaucs<double> *lsys2 = new linearSystemCSRTaucs<double>;
#else
diff --git a/benchmarks/3d/Torus2.geo b/benchmarks/3d/Torus2.geo
index 245f082097f7e64738fb512962545044168a83d3..227ea82289155fc8d5e5914ce14299966c0d8900 100644
--- a/benchmarks/3d/Torus2.geo
+++ b/benchmarks/3d/Torus2.geo
@@ -1,3 +1,6 @@
+
+Mesh.CharacteristicLengthFactor = 0.1;
+
lc = 0.7;
Point(1) = {0,0,0,lc};
Point(2) = {1,0,0,lc};
@@ -14,7 +17,7 @@ Extrude Surface {6, {0,1,0}, {-2,0,0}, 2*Pi/3};
Extrude Surface {28, {0,1,0}, {-2,0,0}, 2*Pi/3};
Extrude Surface {50, {0,1,0}, {-2,0,0}, 2*Pi/3};
Surface Loop(72) = {45,23,67,71,49,27,15,59,37,41,19,63};
+//Volume(73) = {72};
//Compound Surface(100)={45,23,67,71,49,27,15,59,37,41,19,63};
-//Volume(73) = {72};
diff --git a/benchmarks/3d/Torus_GEO.geo b/benchmarks/3d/Torus_GEO.geo
deleted file mode 100644
index bdda121a002a85fac8b4d8aec41bdd98b5971184..0000000000000000000000000000000000000000
--- a/benchmarks/3d/Torus_GEO.geo
+++ /dev/null
@@ -1,11 +0,0 @@
-Mesh.CharacteristicLengthFactor=0.4;
-
-Merge "Torus2_CLASS.msh";
-CreateTopology;
-
-Compound Line(100)={1};
-Compound Line(200)={2};
-
-Compound Surface(1000)={150};
-Compound Surface(2000)={250};
-
diff --git a/contrib/Tetgen/tetgen.cxx b/contrib/Tetgen/tetgen.cxx
index 23a5449f67ee9dcbe76e2643c2b4e6ec8d4820dc..76fbc8f53f0a5e224aff8d671a6fd781634c7be7 100644
--- a/contrib/Tetgen/tetgen.cxx
+++ b/contrib/Tetgen/tetgen.cxx
@@ -32386,7 +32386,8 @@ void tetgenmesh::outvoronoi(tetgenio* out)
terminatetetgen(1);
}
// Number of voronoi points, 3 dim, no attributes, no marker.
- fprintf(outfile, "%ld 3 0 0\n", tetrahedrons->items);
+ //fprintf(outfile, "%ld 3 0 0\n", tetrahedrons->items);
+ fprintf(outfile, "View \"voronoi nodes\" {\n");//EMI
} else {
// Allocate space for 'vpointlist'.
out->numberofvpoints = (int) tetrahedrons->items;
@@ -32404,6 +32405,8 @@ void tetgenmesh::outvoronoi(tetgenio* out)
tetloop.tet = tetrahedrontraverse();
vpointcount = 0;
index = 0;
+ REAL *vlist;
+ vlist = new REAL[(int)tetrahedrons->items * 3];
while (tetloop.tet != (tetrahedron *) NULL) {
// Calculate the circumcenter.
for (i = 0; i < 4; i++) {
@@ -32412,8 +32415,12 @@ void tetgenmesh::outvoronoi(tetgenio* out)
}
circumsphere(pt[0], pt[1], pt[2], pt[3], ccent, NULL);
if (out == (tetgenio *) NULL) {
- fprintf(outfile, "%4d %16.8e %16.8e %16.8e\n", vpointcount + shift,
- ccent[0], ccent[1], ccent[2]);
+ //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
+ vlist[(vpointcount+shift)*3+0] = ccent[0];//EMI
+ vlist[(vpointcount+shift)*3+1] = ccent[1];//EMI
+ vlist[(vpointcount+shift)*3+2] = ccent[2];//EMI
} else {
out->vpointlist[index++] = ccent[0];
out->vpointlist[index++] = ccent[1];
@@ -32428,7 +32435,8 @@ void tetgenmesh::outvoronoi(tetgenio* out)
* (int *) (dummytet + elemmarkerindex) = -1;
if (out == (tetgenio *) NULL) {
- fprintf(outfile, "# Generated by %s\n", b->commandline);
+ //fprintf(outfile, "# Generated by %s\n", b->commandline);
+ fprintf(outfile,"};\n"); //EMI
fclose(outfile);
}
@@ -32453,7 +32461,8 @@ void tetgenmesh::outvoronoi(tetgenio* out)
terminatetetgen(1);
}
// Number of Voronoi edges, no marker.
- fprintf(outfile, "%ld 0\n", faces);
+ //sfprintf(outfile, "%ld 0\n", faces);
+ fprintf(outfile, "View \"voronoi edges\" {\n");//EMI
} else {
// Allocate space for 'vpointlist'.
out->numberofedges = (int) faces;
@@ -32479,7 +32488,7 @@ void tetgenmesh::outvoronoi(tetgenio* out)
decode(tetloop.tet[i], worktet);
if ((worktet.tet == dummytet) || (tetloop.tet < worktet.tet)) {
if (out == (tetgenio *) NULL) {
- fprintf(outfile, "%4d %4d", vedgecount + shift, end1 + shift);
+ //fprintf(outfile, "%4d %4d", vedgecount + shift, end1 + shift); //EMI
} else {
vedge = &(out->vedgelist[index++]);
vedge->v1 = end1 + shift;
@@ -32502,8 +32511,8 @@ void tetgenmesh::outvoronoi(tetgenio* out)
+ infvec[2] * infvec[2]);
if (L > 0) for (j = 0; j < 3; j++) infvec[j] /= L;
if (out == (tetgenio *) NULL) {
- fprintf(outfile, " -1");
- fprintf(outfile, " %g %g %g\n", infvec[0], infvec[1], infvec[2]);
+ //fprintf(outfile, " -1");//EMI
+ //fprintf(outfile, " %g %g %g\n", infvec[0], infvec[1], infvec[2]);//EMI
} else {
vedge->v2 = -1;
vedge->vnormal[0] = infvec[0];
@@ -32512,7 +32521,11 @@ void tetgenmesh::outvoronoi(tetgenio* out)
}
} else {
if (out == (tetgenio *) NULL) {
- fprintf(outfile, " %4d\n", end2 + shift);
+ //fprintf(outfile, " %4d\n", end2 + shift);//EMI
+ 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);
} else {
vedge->v2 = end2 + shift;
vedge->vnormal[0] = 0.0;
@@ -32532,7 +32545,8 @@ void tetgenmesh::outvoronoi(tetgenio* out)
}
if (out == (tetgenio *) NULL) {
- fprintf(outfile, "# Generated by %s\n", b->commandline);
+ //fprintf(outfile, "# Generated by %s\n", b->commandline);
+ fprintf(outfile,"};\n"); //EMI
fclose(outfile);
}
@@ -34856,7 +34870,8 @@ void tetrahedralize(tetgenbehavior *b, tetgenio *in, tetgenio *out,
}
if (b->voroout) {
- m.outvoronoi(out);
+ m.outvoronoi(NULL);
+ //m.outvoronoi(out);
}
tv[13] = clock();