diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index abf33469adee1d778df61a6ebbc4b2d6351bd2c5..2988a385269503c94e1c8515cd50a78fc21bbeca 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -1362,6 +1362,7 @@ void GModel::createTopologyFromFaces(std::vector<discreteFace*> &discFaces)
int numE = getMaxElementaryNumber(1) + 1;
discreteEdge *e = new discreteEdge(this, numE, 0, 0);
//printf("*** Created discreteEdge %d \n", numE);
+
add(e);
discEdges.push_back(e);
diff --git a/Geo/GeoInterpolation.cpp b/Geo/GeoInterpolation.cpp
index 47ca545b7a9633e227d8d955d455ed9e7c843552..c600e06584216f6342d584709b4f49bb90e68c9e 100644
--- a/Geo/GeoInterpolation.cpp
+++ b/Geo/GeoInterpolation.cpp
@@ -355,11 +355,11 @@ Vertex InterpolateCurve(Curve *c, double u, int derivee)
break;
case MSH_SEGM_BND_LAYER:
- Msg::Error("Cannot interpolate boundary layer curve");
+ Msg::Debug("Cannot interpolate boundary layer curve");
break;
case MSH_SEGM_DISCRETE:
- Msg::Error("Cannot interpolate discrete curve");
+ Msg::Debug("Cannot interpolate discrete curve");
break;
default:
diff --git a/Geo/MTetrahedron.cpp b/Geo/MTetrahedron.cpp
index 2612ac62ed253b2dd57f5299024b4107a9ff4a3c..145fbd6fa25de081852069422f92c5acea542025 100644
--- a/Geo/MTetrahedron.cpp
+++ b/Geo/MTetrahedron.cpp
@@ -28,6 +28,20 @@ SPoint3 MTetrahedron::circumcenter()
#endif
}
+double MTetrahedron::getCircumRadius()
+{
+#if defined(HAVE_MESH)
+ SPoint3 center = circumcenter();
+ const double dx = getVertex(0)->x() - center.x();
+ const double dy = getVertex(0)->y() - center.y();
+ const double dz = getVertex(0)->z() - center.z();
+ double circum_radius = sqrt(dx * dx + dy * dy + dz * dz);
+ return circum_radius;
+#else
+ return 0.0;
+#endif
+}
+
double MTetrahedron::distoShapeMeasure()
{
#if defined(HAVE_MESH)
diff --git a/Geo/MTetrahedron.h b/Geo/MTetrahedron.h
index e0877afd6f0655d2cb67f62cb307a5556135f617..dde4c80ff144da3d9416a37c328697afc392e741 100644
--- a/Geo/MTetrahedron.h
+++ b/Geo/MTetrahedron.h
@@ -130,6 +130,7 @@ class MTetrahedron : public MElement {
virtual double gammaShapeMeasure();
virtual double distoShapeMeasure();
virtual double getInnerRadius();
+ virtual double getCircumRadius();
virtual double etaShapeMeasure();
void xyz2uvw(double xyz[3], double uvw[3]);
virtual const polynomialBasis* getFunctionSpace(int o=-1) const;
diff --git a/Graphics/drawMesh.cpp b/Graphics/drawMesh.cpp
index 2855e4dd8da28150d2c31184d48b77798dbdce7d..108836ac4d5889b0481073e2b022428eedb9ce61 100644
--- a/Graphics/drawMesh.cpp
+++ b/Graphics/drawMesh.cpp
@@ -945,6 +945,10 @@ class drawMeshGRegion {
drawBarycentricDual(r->polyhedra);
}
+ if(CTX::instance()->mesh.voronoi) {
+ if(CTX::instance()->mesh.tetrahedra) drawVoronoiDual(r->tetrahedra);
+ }
+
if(select) {
glPopName();
glPopName();
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 434ab77d296478ba9812d57c9335601aba756b24..4224df3fd80512491d0fae64e173278642a1b6fd 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -22,6 +22,206 @@
#include "Context.h"
#include "GFaceCompound.h"
+#if defined(HAVE_ANN)
+#include "ANN/ANN.h"
+#endif
+
+void voronoiDual(GRegion *gr)
+{
+
+ std::vector<MTetrahedron*> elements = gr->tetrahedra;
+ std::list<GFace*> allFaces = gr->faces();
+
+ std::set<SPoint3> candidates;
+ std::set<SPoint3> seeds;
+
+ printf("computing box and barycenters\n");
+ double Dmax = -1.0;
+ std::list<GFace*>::iterator itf = allFaces.begin();
+ for(; itf != allFaces.end(); itf ++){
+ SBoundingBox3d bbox = (*itf)->bounds();
+ SVector3 dd(bbox.max(), bbox.min());
+ //if( (*itf)->tag()==5904 || (*itf)->tag()==5902 || (*itf)->tag()==5906) {
+ //if( (*itf)->tag()==6 || (*itf)->tag()==28){
+ if( (*itf)->tag() !=1 ) {
+ Dmax = std::max(Dmax,norm(dd));
+ seeds.insert(bbox.center());
+ }
+ }
+ printf("Dmax =%g \n", Dmax);
+
+ printf("printing nodes \n");
+ FILE *outfile;
+ outfile = fopen("nodes.pos", "w");
+ fprintf(outfile, "View \"voronoi nodes\" {\n");
+ for(unsigned int i = 0; i < elements.size(); i++){
+ MTetrahedron *ele = elements[i];
+ SPoint3 pc = ele->circumcenter();
+ double x[3] = {pc.x(), pc.y(), pc.z()};
+ double uvw[3];
+ ele->xyz2uvw(x, uvw);
+ if(ele->isInside(uvw[0], uvw[1], uvw[2])){
+ double radius = ele->getCircumRadius();
+ if(radius > Dmax/5.) {
+ candidates.insert(pc);
+ fprintf(outfile,"SP(%g,%g,%g) {%g};\n",
+ pc.x(), pc.y(), pc.z(), radius);
+ }
+ }
+ }
+ fprintf(outfile,"};\n");
+ fclose(outfile);
+
+ // printf("building maps \n");
+ // std::map<MVertex*, std::set<MTetrahedron*> > node2Tet;
+ // std::map<MFace, std::vector<MTetrahedron*> , Less_Face> face2Tet;
+ // for(unsigned int i = 0; i < elements.size(); i++){
+ // MTetrahedron *ele = elements[i];
+ // for (int j=0; j<4; j++){
+ // MVertex *v = ele->getVertex(j);
+ // std::map<MVertex*, std::set<MTetrahedron*> >::iterator itmap = node2Tet.find(v);
+ // if (itmap == node2Tet.end()){
+ // std::set<MTetrahedron*> oneTet;
+ // oneTet.insert(ele);
+ // node2Tet.insert(std::make_pair(v, oneTet));
+ // }
+ // else{
+ // std::set<MTetrahedron*> allTets = itmap->second;
+ // allTets.insert(allTets.begin(), ele);
+ // itmap->second = allTets;
+ // }
+ // }
+ // for (int j=0; j<4; j++){
+ // MFace f = ele->getFace(j);
+ // std::map<MFace, std::vector<MTetrahedron*>, Less_Face >::iterator itmap = face2Tet.find(f);
+ // if (itmap == face2Tet.end()){
+ // std::vector<MTetrahedron*> oneTet;
+ // oneTet.push_back(ele);
+ // face2Tet.insert(std::make_pair(f, oneTet));
+ // }
+ // else{
+ // std::vector<MTetrahedron*> allTets = itmap->second;
+ // allTets.insert(allTets.begin(), ele);
+ // itmap->second = allTets;
+ // }
+ // }
+ // }
+
+ //print only poles
+ // FILE *outfile;
+ // outfile = fopen("nodes.pos", "w");
+ // fprintf(outfile, "View \"voronoi nodes\" {\n");
+ // std::map<MVertex*, std::set<MTetrahedron*> >::iterator itmap = node2Tet.begin();
+ // for(; itmap != node2Tet.end(); itmap++){
+ // MVertex *v = itmap->first;
+ // std::set<MTetrahedron*> allTets = itmap->second;
+ // std::set<MTetrahedron*>::const_iterator it = allTets.begin();
+ // MTetrahedron *poleTet = *it;
+ // double maxRadius = poleTet->getCircumRadius();
+ // for(; it != allTets.end(); it++){
+ // double radius = (*it)->getCircumRadius();
+ // if (radius > maxRadius){
+ // maxRadius = radius;
+ // poleTet = *it;
+ // }
+ // }
+ // SPoint3 pc = poleTet->circumcenter();
+ // fprintf(outfile,"SP(%g,%g,%g) {%g};\n",
+ // pc.x(), pc.y(), pc.z(), maxRadius);
+ // candidates.insert(pc);
+ // }
+ // fprintf(outfile,"};\n");
+ // fclose(outfile);
+
+ printf("Ann computation of neighbours and writing edges\n");
+ #if defined(HAVE_ANN)
+// FILE *outfile2;
+// outfile2 = fopen("edges.pos", "w");
+// fprintf(outfile2, "View \"voronoi edges\" {\n");
+
+// ANNkd_tree *_kdtree;
+// ANNpointArray _zeronodes;
+// ANNidxArray _index;
+// ANNdistArray _dist;
+
+// std::set<SPoint3>::iterator itseed = seeds.begin();
+// SPoint3 beginPt=*itseed;
+// seeds.erase(itseed);
+// itseed = seeds.begin();
+// for(; itseed != seeds.end(); itseed++){
+// printf("seed =%g %g %g \n", (*itseed).x(), (*itseed).y(), (*itseed).z());
+// candidates.insert(*itseed);
+// }
+// printf("begin seed =%g %g %g \n", beginPt.x(), beginPt.y(), beginPt.z());
+
+// double color = 1.;
+// while(candidates.size()>0){
+
+// _zeronodes = annAllocPts(candidates.size(), 3);
+// std::set<SPoint3>::iterator itset = candidates.begin();
+// int i=0;
+// for(; itset != candidates.end(); itset++){
+// _zeronodes[i][0] = (*itset).x();
+// _zeronodes[i][1] = (*itset).y();
+// _zeronodes[i][2] = (*itset).z();
+// i++;
+// }
+// _kdtree = new ANNkd_tree(_zeronodes, candidates.size(), 3);
+// _index = new ANNidx[1];
+// _dist = new ANNdist[1];
+
+// double xyz[3] = {beginPt.x(), beginPt.y(), beginPt.z()};
+// _kdtree->annkSearch(xyz, 1, _index, _dist);
+// SPoint3 endPt( _zeronodes[_index[0]][0], _zeronodes[_index[0]][1], _zeronodes[_index[0]][2]);
+// fprintf(outfile2,"SL(%g,%g,%g,%g,%g,%g) {%g,%g};\n",
+// beginPt.x(), beginPt.y(), beginPt.z(),
+// endPt.x(), endPt.y(), endPt.z(),
+// color, color);
+
+// std::set<SPoint3>::iterator itse=seeds.find(endPt);
+// std::set<SPoint3>::iterator its=candidates.find (endPt);
+// if(itse != seeds.end()){
+// seeds.erase(itse);
+// beginPt = *(seeds.begin());
+// std::set<SPoint3>::iterator itsee=candidates.find(beginPt);
+// candidates.erase(itsee);
+// color=color*2.;
+// }
+// else beginPt=endPt;
+
+// if(its != candidates.end()) {
+// candidates.erase(its);
+// }
+
+// delete _kdtree;
+// annDeallocPts(_zeronodes);
+// delete [] _index;
+// delete [] _dist;
+// }
+
+// fprintf(outfile2,"};\n");
+// fclose(outfile2);
+#endif
+
+ //print scalar lines
+ // std::map<MFace, std::vector<MTetrahedron*> , Less_Face >::iterator itmap2 = face2Tet.begin();
+ // for(; itmap2 != face2Tet.end(); itmap2++){
+ // std::vector<MTetrahedron*> allTets = itmap2->second;
+ // if (allTets.size() !=2 ) continue;
+ // SPoint3 pc1 = allTets[0]->circumcenter();
+ // SPoint3 pc2 = allTets[1]->circumcenter();
+ // std::set<SPoint3>::const_iterator it1 = candidates.find(pc1);
+ // std::set<SPoint3>::const_iterator it2 = candidates.find(pc2);
+ // if( it1 != candidates.end() || it2 != candidates.end())
+ // fprintf(outfile2,"SL(%g,%g,%g,%g,%g,%g) {%g,%g};\n",
+ // pc1.x(), pc1.y(), pc1.z(), pc2.x(), pc2.y(), pc2.z(),
+ // allTets[0]->getCircumRadius(),allTets[1]->getCircumRadius());
+ // }
+
+
+}
+
+
void getAllBoundingVertices(GRegion *gr, std::set<MVertex*> &allBoundingVertices)
{
std::list<GFace*> faces = gr->faces();
@@ -135,7 +335,7 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out,
gf->deleteVertexArrays();
++it;
}
-
+
// TODO: re-create 1D mesh
for(int i = 0; i < out.numberofedges; i++){
MVertex *v[2];
@@ -145,6 +345,7 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out,
//implement here the 1D mesh ...
}
+ bool needParam = (CTX::instance()->mesh.order > 1 && CTX::instance()->mesh.secondOrderExperimental);
// re-create the triangular meshes FIXME: this can lead to hanging
// nodes for non manifold geometries (single surface connected to
// volume)
@@ -156,43 +357,45 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out,
GFace *gf = gr->model()->getFaceByTag(out.trifacemarkerlist[i]);
double guess[2] = {0, 0};
- int Count = 0;
- for(int j = 0; j < 3; j++){
- if(!v[j]->onWhat()){
- Msg::Error("Uncategorized vertex %d", v[j]->getNum());
- }
- else if(v[j]->onWhat()->dim() == 2){
- double uu,vv;
- v[j]->getParameter(0, uu);
- v[j]->getParameter(1, vv);
- guess[0] += uu;
- guess[1] += vv;
- Count++;
+ if (needParam) {
+ int Count = 0;
+ for(int j = 0; j < 3; j++){
+ if(!v[j]->onWhat()){
+ Msg::Error("Uncategorized vertex %d", v[j]->getNum());
+ }
+ else if(v[j]->onWhat()->dim() == 2){
+ double uu,vv;
+ v[j]->getParameter(0, uu);
+ v[j]->getParameter(1, vv);
+ guess[0] += uu;
+ guess[1] += vv;
+ Count++;
+ }
+ else if(v[j]->onWhat()->dim() == 1){
+ GEdge *ge = (GEdge*)v[j]->onWhat();
+ double UU;
+ v[j]->getParameter(0, UU);
+ SPoint2 param;
+ param = ge->reparamOnFace(gf, UU, 1);
+ guess[0] += param.x();
+ guess[1] += param.y();
+ Count++;
+ }
+ else if(v[j]->onWhat()->dim() == 0){
+ SPoint2 param;
+ GVertex *gv = (GVertex*)v[j]->onWhat();
+ param = gv->reparamOnFace(gf,1);
+ guess[0] += param.x();
+ guess[1] += param.y();
+ Count++;
+ }
}
- else if(v[j]->onWhat()->dim() == 1){
- GEdge *ge = (GEdge*)v[j]->onWhat();
- double UU;
- v[j]->getParameter(0, UU);
- SPoint2 param;
- param = ge->reparamOnFace(gf, UU, 1);
- guess[0] += param.x();
- guess[1] += param.y();
- Count++;
+ if(Count != 0){
+ guess[0] /= Count;
+ guess[1] /= Count;
}
- else if(v[j]->onWhat()->dim() == 0){
- SPoint2 param;
- GVertex *gv = (GVertex*)v[j]->onWhat();
- param = gv->reparamOnFace(gf,1);
- guess[0] += param.x();
- guess[1] += param.y();
- Count++;
- }
- }
- if(Count != 0){
- guess[0] /= Count;
- guess[1] /= Count;
}
-
+
for(int j = 0; j < 3; j++){
if(v[j]->onWhat()->dim() == 3){
v[j]->onWhat()->mesh_vertices.erase
@@ -200,12 +403,11 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out,
v[j]->onWhat()->mesh_vertices.end(),
v[j]));
MVertex *v1b;
- if(CTX::instance()->mesh.order > 1 &&
- CTX::instance()->mesh.secondOrderExperimental){
+ if(needParam){
// parametric coordinates should be set for the vertex !
// (this is not 100 % safe yet, so we reserve that operation
// for high order meshes)
- GPoint gp = gf->closestPoint(SPoint3(v[j]->x(), v[j]->y(), v[j]->z()),guess);
+ GPoint gp = gf->closestPoint(SPoint3(v[j]->x(), v[j]->y(), v[j]->z()), guess);
if(gp.g()){
v1b = new MFaceVertex(gp.x(), gp.y(), gp.z(), gf, gp.u(), gp.v());
}
@@ -228,6 +430,7 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out,
MTriangle *t = new MTriangle(v[0], v[1], v[2]);
gf->triangles.push_back(t);
}
+
for(int i = 0; i < out.numberoftetrahedra; i++){
MVertex *v1 = numberedV[out.tetrahedronlist[i * 4 + 0] - 1];
MVertex *v2 = numberedV[out.tetrahedronlist[i * 4 + 1] - 1];
@@ -236,6 +439,8 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out,
MTetrahedron *t = new MTetrahedron(v1, v2, v3, v4);
gr->tetrahedra.push_back(t);
}
+
+
}
#endif
@@ -273,10 +478,8 @@ void MeshDelaunayVolume(std::vector<GRegion*> ®ions)
std::vector<MVertex*> numberedV;
char opts[128];
buildTetgenStructure(gr, in, numberedV);
- if (Msg::GetVerbosity() == 10)
- sprintf(opts, "peVvS0s");
- else
- sprintf(opts, "pe%c", (Msg::GetVerbosity() < 3) ? 'Q': (Msg::GetVerbosity() > 6)? 'V': '\0');
+ //if (Msg::GetVerbosity() == 20) sprintf(opts, "peVvS0");
+ sprintf(opts, "pe%c", (Msg::GetVerbosity() < 3) ? 'Q': (Msg::GetVerbosity() > 6)? 'V': '\0');
try{
tetrahedralize(opts, &in, &out);
}
@@ -312,10 +515,10 @@ void MeshDelaunayVolume(std::vector<GRegion*> ®ions)
gr->set(faces);
return;
}
+ printf("transfer tetgen \n");
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()){
@@ -327,6 +530,9 @@ void MeshDelaunayVolume(std::vector<GRegion*> ®ions)
// restore the initial set of faces
gr->set(faces);
+ //EMI VORONOI FOR CENRTERLINE
+ if (Msg::GetVerbosity() == 20) voronoiDual(gr);
+
// now do insertion of points
insertVerticesInRegion(gr);
@@ -374,7 +580,6 @@ Ng_Mesh *buildNetgenStructure(GRegion *gr, bool importVolumeMesh,
Ng_AddPoint(ngmesh, tmp);
}
}
-
std::list<GFace*> faces = gr->faces();
std::list<GFace*>::iterator it = faces.begin();
while(it != faces.end()){
diff --git a/benchmarks/3d/Cube-01.geo b/benchmarks/3d/Cube-01.geo
index be4d5d1231858a77221a4726130f83b83ff941b1..db9342fbcbfd857412d16cfa5d1209618977baad 100644
--- a/benchmarks/3d/Cube-01.geo
+++ b/benchmarks/3d/Cube-01.geo
@@ -1,4 +1,7 @@
-lc = 0.1;
+//Mesh.Dual = 1;
+Mesh.Voronoi=1;
+
+lc = 0.5;
Point(1) = {0.0,0.0,0.0,lc};
Point(2) = {1,0.0,0.0,lc};
Point(3) = {1,1,0.0,lc};
diff --git a/benchmarks/3d/Torus.geo b/benchmarks/3d/Torus.geo
index 56489b2a6ae6d4ab27912aae7f6a5f6fdcee510d..e8ef42bd7f589b16d2968ae165ddf90e4c6c6dda 100644
--- a/benchmarks/3d/Torus.geo
+++ b/benchmarks/3d/Torus.geo
@@ -1,4 +1,4 @@
-lc = 15.5;
+lc = 0.1;
Point(1) = {2.0,0.0,0.0,lc};
Point(2) = {2.0,1,0.0,lc};
Point(3) = {1,0,0.0,lc};
diff --git a/contrib/Tetgen/tetgen.cxx b/contrib/Tetgen/tetgen.cxx
index 32bf180d1ef92e58d10f7419244e92b32b1a0139..cd74509640b9c660d4692e5e2ee91432bede203a 100644
--- a/contrib/Tetgen/tetgen.cxx
+++ b/contrib/Tetgen/tetgen.cxx
@@ -31864,7 +31864,7 @@ void tetgenmesh::outsubfaces(tetgenio* out)
if (out == (tetgenio *) NULL) {
printf("Writing %s.\n", facefilename);
} else {
- printf("Writing faces.\n");
+ printf("Writing subfaces.\n");
}
}
@@ -32357,6 +32357,11 @@ void tetgenmesh::outvoronoi(tetgenio* out)
int end1, end2;
int hitbdry, i, j, k;
+
+ //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ //NODES
+ //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
// Output Voronoi vertices to .v.node file.
if (out == (tetgenio *) NULL) {
strcpy(outfilename, b->outfilename);
@@ -32406,9 +32411,10 @@ 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 ];
+ REAL *radlist;
+ vlist = new REAL[(int)tetrahedrons->items * 3]; //for coords of circumcenters
+ radlist = new REAL[(int)tetrahedrons->items ]; //for radii of circumcenters
+
while (tetloop.tet != (tetrahedron *) NULL) {
// Calculate the circumcenter.
for (i = 0; i < 4; i++) {
@@ -32421,10 +32427,12 @@ void tetgenmesh::outvoronoi(tetgenio* out)
//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)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
+ int indexCenter = (vpointcount+shift);//EMI
+ vlist[indexCenter*3+0] = ccent[0];//EMI
+ vlist[indexCenter*3+1] = ccent[1];//EMI
+ vlist[indexCenter*3+2] = ccent[2];//EMI
+ radlist[indexCenter] = radius;//EMI
+
} else {
out->vpointlist[index++] = ccent[0];
out->vpointlist[index++] = ccent[1];
@@ -32435,6 +32443,7 @@ void tetgenmesh::outvoronoi(tetgenio* out)
vpointcount++;
tetloop.tet = tetrahedrontraverse();
}
+
// Set the outside element marker.
* (int *) (dummytet + elemmarkerindex) = -1;
@@ -32444,6 +32453,10 @@ void tetgenmesh::outvoronoi(tetgenio* out)
fclose(outfile);
}
+ //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ //EDGES
+ //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
// Output Voronoi edges to .v.edge file.
if (out == (tetgenio *) NULL) {
strcpy(outfilename, b->outfilename);
@@ -32526,10 +32539,12 @@ void tetgenmesh::outvoronoi(tetgenio* out)
} else {
if (out == (tetgenio *) NULL) {
//fprintf(outfile, " %4d\n", end2 + shift);//EMI
+ int v1 = (end1+shift);
+ int v2 = (end2+shift);
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],
- vlist2[(end1+shift)], vlist2[(end2+shift)]);
+ vlist[v1*3+0], vlist[v1*3+1], vlist[v1*3+2],
+ vlist[v2*3+0], vlist[v2*3+1], vlist[v2*3+2],
+ radlist[v1], radlist[v2]);
} else {
vedge->v2 = end2 + shift;
vedge->vnormal[0] = 0.0;
@@ -32554,6 +32569,10 @@ void tetgenmesh::outvoronoi(tetgenio* out)
fclose(outfile);
}
+ //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ // FACES
+ //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
// Output Voronoi faces to .v.face file.
if (out == (tetgenio *) NULL) {
strcpy(outfilename, b->outfilename);
@@ -32682,6 +32701,10 @@ void tetgenmesh::outvoronoi(tetgenio* out)
fclose(outfile);
}
+ //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+ // CELLS
+ //%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
// Output Voronoi cells to .v.cell file.
if (out == (tetgenio *) NULL) {
strcpy(outfilename, b->outfilename);