diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index ebdaa804e5ee63dde13ebc672dd91f8d61556482..68bcd9234143605001da311ff6249b2cbe396e0a 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -2556,8 +2556,13 @@ int GModel::readVTK(const std::string &name, bool bigEndian)
if(!strcmp(buffer, "BINARY")) binary = true;
if(fscanf(fp, "%s %s", buffer, buffer2) != 2) return 0;
- if(strcmp(buffer, "DATASET") || strcmp(buffer2, "UNSTRUCTURED_GRID")){
- Msg::Error("VTK reader can only read unstructured datasets");
+
+ bool unstructured = false;
+ if( !strcmp(buffer, "DATASET") && !strcmp(buffer2, "UNSTRUCTURED_GRID") ) unstructured = true;
+
+ if( (strcmp(buffer, "DATASET") && strcmp(buffer2, "UNSTRUCTURED_GRID")) ||
+ (strcmp(buffer, "DATASET") && strcmp(buffer2, "POLYDATA")) ){
+ Msg::Error("VTK reader can only read unstructured or polydata datasets");
return 0;
}
@@ -2603,11 +2608,14 @@ int GModel::readVTK(const std::string &name, bool bigEndian)
// read mesh elements
int numElements, totalNumInt;
if(fscanf(fp, "%s %d %d\n", buffer, &numElements, &totalNumInt) != 3) return 0;
- if(strcmp(buffer, "CELLS") || !numElements){
- Msg::Warning("No cells in dataset");
+
+ if( !strcmp(buffer, "CELLS") && numElements>0 ) Msg::Info("Reading %d cells", numElements);
+ else if (!strcmp(buffer, "POLYGONS") && numElements>0 ) Msg::Info("Reading %d polygons", numElements);
+ else{
+ Msg::Warning("No cells or polygons in dataset");
return 0;
}
- Msg::Info("Reading %d cells", numElements);
+
std::vector<std::vector<MVertex*> > cells(numElements);
for(unsigned int i = 0; i < cells.size(); i++){
int numVerts, n[100];
@@ -2630,33 +2638,50 @@ int GModel::readVTK(const std::string &name, bool bigEndian)
Msg::Error("Bad vertex index");
}
}
- if(fscanf(fp, "%s %d\n", buffer, &numElements) != 2) return 0;
- if(strcmp(buffer, "CELL_TYPES") || numElements != (int)cells.size()){
- Msg::Error("No or invalid number of cells types");
- return 0;
- }
+
std::map<int, std::vector<MElement*> > elements[8];
- for(unsigned int i = 0; i < cells.size(); i++){
- int type;
- if(binary){
- if(fread(&type, sizeof(int), 1, fp) != 1) return 0;
- if(!bigEndian) SwapBytes((char*)&type, sizeof(int), 1);
+ if (unstructured){
+ if(fscanf(fp, "%s %d\n", buffer, &numElements) != 2 ) return 0;
+ if(strcmp(buffer, "CELL_TYPES") || numElements != (int)cells.size()){
+ Msg::Error("No or invalid number of cells types");
+ return 0;
}
- else{
- if(fscanf(fp, "%d", &type) != 1) return 0;
- }
- switch(type){
- case 1: elements[0][1].push_back(new MPoint(cells[i])); break;
- case 3: elements[1][1].push_back(new MLine(cells[i])); break;
- case 5: elements[2][1].push_back(new MTriangle(cells[i])); break;
- case 9: elements[3][1].push_back(new MQuadrangle(cells[i])); break;
- case 10: elements[4][1].push_back(new MTetrahedron(cells[i])); break;
- case 12: elements[5][1].push_back(new MHexahedron(cells[i])); break;
- case 13: elements[6][1].push_back(new MPrism(cells[i])); break;
- case 14: elements[7][1].push_back(new MPyramid(cells[i])); break;
- default:
- Msg::Error("Unknown type of cell %d", type);
- break;
+ for(unsigned int i = 0; i < cells.size(); i++){
+ int type;
+ if(binary){
+ if(fread(&type, sizeof(int), 1, fp) != 1) return 0;
+ if(!bigEndian) SwapBytes((char*)&type, sizeof(int), 1);
+ }
+ else{
+ if(fscanf(fp, "%d", &type) != 1) return 0;
+ }
+ switch(type){
+ case 1: elements[0][1].push_back(new MPoint(cells[i])); break;
+ case 3: elements[1][1].push_back(new MLine(cells[i])); break;
+ case 5: elements[2][1].push_back(new MTriangle(cells[i])); break;
+ case 9: elements[3][1].push_back(new MQuadrangle(cells[i])); break;
+ case 10: elements[4][1].push_back(new MTetrahedron(cells[i])); break;
+ case 12: elements[5][1].push_back(new MHexahedron(cells[i])); break;
+ case 13: elements[6][1].push_back(new MPrism(cells[i])); break;
+ case 14: elements[7][1].push_back(new MPyramid(cells[i])); break;
+ default:
+ Msg::Error("Unknown type of cell %d", type);
+ break;
+ }
+ }
+ }
+ else{
+ for(unsigned int i = 0; i < cells.size(); i++){
+ int nbNodes = (int)cells[i].size();
+ switch(nbNodes){
+ case 1: elements[0][1].push_back(new MPoint(cells[i])); break;
+ case 2: elements[1][1].push_back(new MLine(cells[i])); break;
+ case 3: elements[2][1].push_back(new MTriangle(cells[i])); break;
+ case 4: elements[3][1].push_back(new MQuadrangle(cells[i])); break;
+ default:
+ Msg::Error("Unknown type of mesh element with %d nodes", nbNodes);
+ break;
+ }
}
}
diff --git a/Geo/MEdge.h b/Geo/MEdge.h
index cc7eb24b6eca274cb6c0fee0b491b10b20d9b1c4..ad517128cbe0dcda0c63cfc0f6d460d83355a907 100644
--- a/Geo/MEdge.h
+++ b/Geo/MEdge.h
@@ -53,6 +53,13 @@ class MEdge {
t.normalize();
return t;
}
+ double length() const
+ {
+ SVector3 t(_v[1]->x() - _v[0]->x(),
+ _v[1]->y() - _v[0]->y(),
+ _v[1]->z() - _v[0]->z());
+ return t.norm();
+ }
SVector3 normal() const
{
// this computes one of the normals to the edge
diff --git a/Geo/MTetrahedron.h b/Geo/MTetrahedron.h
index dde4c80ff144da3d9416a37c328697afc392e741..cdbde40e8b96e6500133810721fa193f41be13df 100644
--- a/Geo/MTetrahedron.h
+++ b/Geo/MTetrahedron.h
@@ -70,6 +70,14 @@ class MTetrahedron : public MElement {
{
return MEdge(_v[edges_tetra(num, 0)], _v[edges_tetra(num, 1)]);
}
+ virtual double getMaxEdgeLength(){
+ double maxe = 0.0;
+ for(int i = 0; i < 6; i++){
+ double loce = getEdge(i).length();
+ maxe = std::max(maxe,loce);
+ }
+ return maxe;
+ }
virtual int getNumEdgesRep(){ return 6; }
virtual void getEdgeRep(int num, double *x, double *y, double *z, SVector3 *n)
{
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 06005de6fad87f41dc1828c1c255b136181a36d4..ddc2525f09f1c8b1e6924661afdc59275988108f 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -21,17 +21,17 @@
#include "BDS.h"
#include "Context.h"
#include "GFaceCompound.h"
+#include "SmoothData.h"
#if defined(HAVE_ANN)
#include "ANN/ANN.h"
#endif
-void printVoronoi(GRegion *gr)
+void printVoronoi(GRegion *gr, std::set<SPoint3> &candidates)
{
std::vector<MTetrahedron*> elements = gr->tetrahedra;
std::list<GFace*> allFaces = gr->faces();
- std::set<SPoint3> candidates;
//building maps
std::map<MVertex*, std::set<MTetrahedron*> > node2Tet;
@@ -68,10 +68,11 @@ void printVoronoi(GRegion *gr)
}
}
- //print voronoi nodes
+ //print voronoi poles
FILE *outfile;
+ smooth_normals *snorm = gr->model()->normals;
outfile = fopen("nodes.pos", "w");
- fprintf(outfile, "View \"voronoi nodes\" {\n");
+ fprintf(outfile, "View \"Voronoi poles\" {\n");
std::map<MVertex*, std::set<MTetrahedron*> >::iterator itmap = node2Tet.begin();
for(; itmap != node2Tet.end(); itmap++){
MVertex *v = itmap->first;
@@ -79,6 +80,7 @@ void printVoronoi(GRegion *gr)
std::set<MTetrahedron*>::const_iterator it = allTets.begin();
MTetrahedron *poleTet = *it;
double maxRadius = poleTet->getCircumRadius();
+ double maxEdgeLength = poleTet->getMaxEdgeLength();
for(; it != allTets.end(); it++){
double radius = (*it)->getCircumRadius();
if (radius > maxRadius){
@@ -86,12 +88,23 @@ void printVoronoi(GRegion *gr)
poleTet = *it;
}
}
- SPoint3 pc = poleTet->circumcenter();
- fprintf(outfile,"SP(%g,%g,%g) {%g};\n",
- pc.x(), pc.y(), pc.z(), maxRadius);
- candidates.insert(pc);
- //}//uncomment this
-
+ if (v->onWhat()->dim() == 2 && maxRadius < maxEdgeLength){
+ SPoint3 pc = poleTet->circumcenter();
+ // double nx,ny,nz;
+ // SVector3 vN = snorm->get(v->x(), v->y(), v->z(), nx,ny,nz);
+ // SVector3 pcv(pc.x()-nx, pc.y()-ny,pc.z()-nz);
+ // printf("nx=%g ny=%g nz=%g dot=%g \n", nx,ny,nz, dot(vN, pcv));
+ // if ( dot(vN, pcv) > 0.0 )
+ double x[3] = {pc.x(), pc.y(), pc.z()};
+ double uvw[3];
+ poleTet->xyz2uvw(x, uvw);
+ bool inside = poleTet->isInside(uvw[0], uvw[1], uvw[2]);
+ if (inside){
+ fprintf(outfile,"SP(%g,%g,%g) {%g};\n",
+ pc.x(), pc.y(), pc.z(), maxRadius);
+ candidates.insert(pc);
+ }
+ }
}
fprintf(outfile,"};\n");
fclose(outfile);
@@ -99,26 +112,26 @@ void printVoronoi(GRegion *gr)
//print scalar lines
FILE *outfile2;
outfile2 = fopen("edges.pos", "w");
- fprintf(outfile2, "View \"voronoi edges\" {\n");
+ fprintf(outfile2, "View \"Voronoi edges\" {\n");
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;
+ 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());
- }
+ //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());
+ }
fprintf(outfile2,"};\n");
fclose(outfile2);
}
-
-void skeletonFromVoronoi(GRegion *gr)
+//------------------------------------------------------------------------
+void skeletonFromVoronoi(GRegion *gr, std::set<SPoint3> &voronoiPoles)
{
std::vector<MTetrahedron*> elements = gr->tetrahedra;
@@ -152,7 +165,11 @@ void skeletonFromVoronoi(GRegion *gr)
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])){
+
+ std::set<SPoint3>::const_iterator it2 = voronoiPoles.find(pc);
+
+ if(ele->isInside(uvw[0], uvw[1], uvw[2]) &&
+ it2 != voronoiPoles.end()){
double radius = ele->getCircumRadius();
if(radius > Dmax/10.) {
candidates.insert(pc);
@@ -233,11 +250,8 @@ void skeletonFromVoronoi(GRegion *gr)
fclose(outfile2);
#endif
-
-
-
}
-
+//------------------------------------------------------------------------
void getAllBoundingVertices(GRegion *gr, std::set<MVertex*> &allBoundingVertices)
{
@@ -255,11 +269,9 @@ void getAllBoundingVertices(GRegion *gr, std::set<MVertex*> &allBoundingVertices
++it;
}
}
-
+//------------------------------------------------------------------------
#if defined(HAVE_TETGEN)
-
#include "tetgen.h"
-
void buildTetgenStructure(GRegion *gr, tetgenio &in, std::vector<MVertex*> &numberedV)
{
std::set<MVertex*> allBoundingVertices;
@@ -455,9 +467,9 @@ void TransferTetgenMesh(GRegion *gr, tetgenio &in, tetgenio &out,
gr->tetrahedra.push_back(t);
}
}
-
#endif
+//------------------------------------------------------------------------
void MeshDelaunayVolume(std::vector<GRegion*> ®ions)
{
if(regions.empty()) return;
@@ -491,7 +503,7 @@ void MeshDelaunayVolume(std::vector<GRegion*> ®ions)
std::vector<MVertex*> numberedV;
char opts[128];
buildTetgenStructure(gr, in, numberedV);
- //if (Msg::GetVerbosity() == 20) sprintf(opts, "peVvS0");
+ //if (Msg::GetVerbosity() == 20) sprintf(opts, "peVvS0");
sprintf(opts, "pe%c", (Msg::GetVerbosity() < 3) ? 'Q':
(Msg::GetVerbosity() > 6) ? 'V': '\0');
try{
@@ -545,8 +557,9 @@ void MeshDelaunayVolume(std::vector<GRegion*> ®ions)
//EMI VORONOI FOR CENRTERLINE
if (Msg::GetVerbosity() == 20) {
- printVoronoi(gr);
- skeletonFromVoronoi(gr);
+ std::set<SPoint3> candidates;
+ printVoronoi(gr, candidates);
+ skeletonFromVoronoi(gr, candidates);
}
// now do insertion of points