diff --git a/Common/CommandLine.cpp b/Common/CommandLine.cpp
index c99ac805b0a3db56e8bcff9667d75445059a73e7..7fb7b2c816f579b4d623b1108014e7d39bedcbb2 100644
--- a/Common/CommandLine.cpp
+++ b/Common/CommandLine.cpp
@@ -518,6 +518,8 @@ void GetOptions(int argc, char *argv[])
CTX::instance()->mesh.fileFormat = FORMAT_UNV;
else if(!strcmp(argv[i], "vrml"))
CTX::instance()->mesh.fileFormat = FORMAT_VRML;
+ else if(!strcmp(argv[i], "ply2"))
+ CTX::instance()->mesh.fileFormat = FORMAT_PLY2;
else if(!strcmp(argv[i], "stl"))
CTX::instance()->mesh.fileFormat = FORMAT_STL;
else if(!strcmp(argv[i], "mesh"))
diff --git a/Common/CreateFile.cpp b/Common/CreateFile.cpp
index c359a1d0a030c72a027da7b920ecac3f5cfef1f4..29bb6d9a0013e23beace59ad6b84dd3d111a69ee 100644
--- a/Common/CreateFile.cpp
+++ b/Common/CreateFile.cpp
@@ -48,6 +48,7 @@ int GuessFileFormatFromFileName(std::string fileName)
else if(ext == ".p3d") return FORMAT_P3D;
else if(ext == ".wrl") return FORMAT_VRML;
else if(ext == ".vrml") return FORMAT_VRML;
+ else if(ext == ".ply2") return FORMAT_PLY2;
else if(ext == ".gif") return FORMAT_GIF;
else if(ext == ".jpg") return FORMAT_JPEG;
else if(ext == ".jpeg") return FORMAT_JPEG;
@@ -90,6 +91,7 @@ std::string GetDefaultFileName(int format)
case FORMAT_INP: name += ".inp"; break;
case FORMAT_P3D: name += ".p3d"; break;
case FORMAT_VRML: name += ".wrl"; break;
+ case FORMAT_PLY2: name += ".ply2"; break;
case FORMAT_GIF: name += ".gif"; break;
case FORMAT_JPEG: name += ".jpg"; break;
case FORMAT_MPEG: name += ".mpg"; break;
@@ -202,6 +204,10 @@ void CreateOutputFile(std::string fileName, int format)
(fileName, CTX::instance()->mesh.saveAll, CTX::instance()->mesh.scalingFactor);
break;
+ case FORMAT_PLY2:
+ GModel::current()->writePLY2(fileName);
+ break;
+
case FORMAT_UNV:
GModel::current()->writeUNV
(fileName, CTX::instance()->mesh.saveAll, CTX::instance()->mesh.saveGroupsOfNodes,
diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h
index ceda6e2fbc63ba030c8559f408f3a36b18f4fbca..f188b5deb490f8b791b028568cb0d6b78d6b5b01 100644
--- a/Common/DefaultOptions.h
+++ b/Common/DefaultOptions.h
@@ -1095,7 +1095,7 @@ StringXNumber MeshOptions_Number[] = {
{ F|O, "Format" , opt_mesh_file_format , FORMAT_AUTO ,
"Mesh output format (1=msh, 2=unv, 10=automatic, 19=vrml, 27=stl, 30=mesh, 31=bdf, "
- "32=cgns, 33=med)" },
+ "32=cgns, 33=med, 40=ply2)" },
{ F|O, "Hexahedra" , opt_mesh_hexahedra , 1. ,
"Display mesh hexahedra?" },
diff --git a/Common/GmshDefines.h b/Common/GmshDefines.h
index ff39bdc25d41ae795950d6022be4eb0f5c148e44..6baf3c1854d4bd34e878c1c68a11b68236a9e46b 100644
--- a/Common/GmshDefines.h
+++ b/Common/GmshDefines.h
@@ -43,6 +43,7 @@
#define FORMAT_IGES 37
#define FORMAT_IR3 38
#define FORMAT_INP 39
+#define FORMAT_PLY2 40
// Element types
#define TYPE_PNT 1
diff --git a/Common/OpenFile.cpp b/Common/OpenFile.cpp
index 00218141cfbb4efd2e77342cc9cbecccff398c33..d0098daf3f31f08ac23fea9d8c6322331a3bfda6 100644
--- a/Common/OpenFile.cpp
+++ b/Common/OpenFile.cpp
@@ -337,6 +337,9 @@ int MergeFile(std::string fileName, bool warnIfMissing)
status = binding::instance()->readFile(fileName.c_str());
}
#endif
+ else if(ext == ".ply2"){
+ status = GModel::current()->readPLY2(fileName);
+ }
else {
CTX::instance()->geom.draw = 1;
if(!strncmp(header, "$PTS", 4) || !strncmp(header, "$NO", 3) ||
diff --git a/Fltk/menuWindow.cpp b/Fltk/menuWindow.cpp
index c4c7db394f457527cb1262ed4baa568ca902dd41..306071679230c9e9bbb781b4848641b66bb46982 100644
--- a/Fltk/menuWindow.cpp
+++ b/Fltk/menuWindow.cpp
@@ -127,6 +127,7 @@ static const char *input_formats =
"STL Surface Mesh" TT "*.stl" NN
"VTK Mesh" TT "*.vtk" NN
"VRML Surface Mesh" TT "*.{wrl,vrml}" NN
+ "PLY2 Surface Mesh" TT "*.{ply2}" NN
SEPARATOR_IN
"BMP" TT "*.bmp" NN
#if defined(HAVE_LIBJPEG)
@@ -271,6 +272,8 @@ static int _save_stl(const char *name){ return genericMeshFileDialog
(name, "STL Options", FORMAT_STL, true, false); }
static int _save_vrml(const char *name){ return genericMeshFileDialog
(name, "VRML Options", FORMAT_VRML, false, false); }
+static int _save_ply2(const char *name){ return genericMeshFileDialog
+ (name, "PLY2 Options", FORMAT_PLY2, false, false); }
static int _save_eps(const char *name){ return gl2psFileDialog
(name, "EPS Options", FORMAT_EPS); }
static int _save_gif(const char *name){ return gifFileDialog(name); }
@@ -309,6 +312,7 @@ static int _save_auto(const char *name)
case FORMAT_IR3 : return _save_ir3(name);
case FORMAT_STL : return _save_stl(name);
case FORMAT_VRML : return _save_vrml(name);
+ case FORMAT_PLY2 : return _save_ply2(name);
case FORMAT_EPS : return _save_eps(name);
case FORMAT_GIF : return _save_gif(name);
case FORMAT_JPEG : return _save_jpeg(name);
@@ -357,6 +361,7 @@ static void file_save_as_cb(Fl_Widget *w, void *data)
{"STL Surface Mesh" TT "*.stl", _save_stl},
{"VRML Surface Mesh" TT "*.wrl", _save_vrml},
{"VTK Mesh" TT "*.vtk", _save_vtk},
+ {"PLY2 Mesh" TT "*.ply2", _save_ply2},
SEPARATOR_OUT
{"Encapsulated PostScript" TT "*.eps", _save_eps},
{"GIF" TT "*.gif", _save_gif},
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 292deca07ec9518b91fea889f70450f80f4c2b9a..76648d43e0d4cf488042531eb539c6c881d1ff11 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -81,6 +81,70 @@ static int intersection_segments (SPoint3 &p1, SPoint3 &p2,
}
}
+
+//--------------------------------------------------------------
+static bool orderVertices(const std::list<GEdge*> &e, std::vector<MVertex*> &l,
+ std::vector<double> &coord)
+{
+ l.clear();
+ coord.clear();
+
+ std::list<GEdge*>::const_iterator it = e.begin();
+ std::list<MLine*> temp;
+ double tot_length = 0;
+ for( ; it != e.end(); ++it ){
+ for(unsigned int i = 0; i < (*it)->lines.size(); i++ ){
+ temp.push_back((*it)->lines[i]);
+ MVertex *v0 = (*it)->lines[i]->getVertex(0);
+ MVertex *v1 = (*it)->lines[i]->getVertex(1);
+ const double length = sqrt((v0->x() - v1->x()) * (v0->x() - v1->x()) +
+ (v0->y() - v1->y()) * (v0->y() - v1->y()) +
+ (v0->z() - v1->z()) * (v0->z() - v1->z()));
+ tot_length += length;
+ }
+ }
+
+ MVertex *first_v = (*temp.begin())->getVertex(0);
+ MVertex *current_v = (*temp.begin())->getVertex(1);
+
+ l.push_back(first_v);
+ coord.push_back(0.0);
+ temp.erase(temp.begin());
+
+ while(temp.size()){
+ bool found = false;
+ for(std::list<MLine*>::iterator itl = temp.begin(); itl != temp.end(); ++itl){
+ MLine *ll = *itl;
+ MVertex *v0 = ll->getVertex(0);
+ MVertex *v1 = ll->getVertex(1);
+ if(v0 == current_v){
+ found = true;
+ l.push_back(current_v);
+ current_v = v1;
+ temp.erase(itl);
+ const double length = sqrt((v0->x() - v1->x()) * (v0->x() - v1->x()) +
+ (v0->y() - v1->y()) * (v0->y() - v1->y()) +
+ (v0->z() - v1->z()) * (v0->z() - v1->z()));
+ coord.push_back(coord[coord.size()-1] + length / tot_length);
+ break;
+ }
+ else if(v1 == current_v){
+ found = true;
+ l.push_back(current_v);
+ current_v = v0;
+ temp.erase(itl);
+ const double length = sqrt((v0->x() - v1->x()) * (v0->x() - v1->x()) +
+ (v0->y() - v1->y()) * (v0->y() - v1->y()) +
+ (v0->z() - v1->z()) * (v0->z() - v1->z()));
+ coord.push_back(coord[coord.size()-1] + length / tot_length);
+ break;
+ }
+ }
+ if(!found) return false;
+ }
+ return true;
+}
+
//--------------------------------------------------------------
static void computeCGKernelPolygon(std::map<MVertex*,SPoint3> &coordinates,
std::vector<MVertex*> &cavV, double &ucg, double &vcg)
@@ -269,12 +333,18 @@ void GFaceCompound::fillNeumannBCS() const
{
fillTris.clear();
+ fillNodes.clear();
- //close neuman bcs
+ //closed interior loops
for(std::list<std::list<GEdge*> >::const_iterator iloop = _interior_loops.begin();
iloop != _interior_loops.end(); iloop++){
+ std::list<MTriangle*> loopfillTris;
std::list<GEdge*> loop = *iloop;
if (loop != _U0 ){
+ std::vector<MVertex*> ordered;
+ std::vector<double> coords;
+ bool success = orderVertices(*iloop, ordered, coords);
+
//--- center of Neumann interior loop
int nb = 0;
double x=0.;
@@ -282,65 +352,116 @@ void GFaceCompound::fillNeumannBCS() const
double z=0.;
//EMI- TODO FIND KERNEL OF POLYGON AND PLACE AT CG KERNEL !
//IF NO KERNEL -> DO NOT FILL TRIS
- for (std::list<GEdge*>::iterator ite = loop.begin(); ite != loop.end(); ite++){
- for (unsigned int k= 0; k< (*ite)->getNumMeshElements(); k++){
- MVertex *v0 = (*ite)->getMeshElement(k)->getVertex(0);
- MVertex *v1 = (*ite)->getMeshElement(k)->getVertex(1);
- x += .5*(v0->x() + v1->x());
- y += .5*(v0->y() + v1->y());
- z += .5*(v0->z() + v1->z());
- nb++;
- }
+ for(unsigned int i = 0; i < ordered.size(); ++i){
+ MVertex *v0 = ordered[i];
+ MVertex *v1 = (i==ordered.size()-1) ? ordered[0]: ordered[i+1];
+ x += .5*(v0->x() + v1->x());
+ y += .5*(v0->y() + v1->y());
+ z += .5*(v0->z() + v1->z());
+ nb++;
}
x/=nb; y/=nb; z/=nb;
MVertex *c = new MVertex(x, y, z);
//--- create new triangles
- for (std::list<GEdge*>::iterator ite = loop.begin(); ite != loop.end(); ite++){
- for (unsigned int i= 0; i< (*ite)->getNumMeshElements(); i++){
- MVertex *v0 = (*ite)->getMeshElement(i)->getVertex(0);
- MVertex *v1 = (*ite)->getMeshElement(i)->getVertex(1);
+ for(unsigned int i = 0; i < ordered.size(); ++i){
+ MVertex *v0 = ordered[i];
+ MVertex *v1 = (i==ordered.size()-1) ? ordered[0]: ordered[i+1];
-// fillTris.push_back(new MTriangle(v0,v1, c));
-
- MVertex *v2 = new MVertex(.5*(v0->x()+c->x()), .5*(v0->y()+c->y()), .5*(v0->z()+c->z()));
- MVertex *v3 = new MVertex(.5*(v1->x()+c->x()), .5*(v1->y()+c->y()), .5*(v1->z()+c->z()));
- fillTris.push_back(new MTriangle(v0,v2,v3));
- fillTris.push_back(new MTriangle(v2,c, v3));
- fillTris.push_back(new MTriangle(v0,v3, v1)) ;
-
-// MVertex *v2 = new MVertex(.66*v0->x()+.33*c->x(), .66*v0->y()+.33*c->y(), .66*v0->z()+.33*c->z());
-// MVertex *v3 = new MVertex(.66*v1->x()+.33*c->x(), .66*v1->y()+.33*c->y(), .66*v1->z()+.33*c->z());
-// MVertex *v4 = new MVertex(.33*v0->x()+.66*c->x(), .33*v0->y()+.66*c->y(), .33*v0->z()+.66*c->z());
-// MVertex *v5 = new MVertex(.33*v1->x()+.66*c->x(), .33*v1->y()+.66*c->y(), .33*v1->z()+.66*c->z());
-// fillTris.push_back(new MTriangle(v0,v2,v3));
-// fillTris.push_back(new MTriangle(v2,v5,v3));
-// fillTris.push_back(new MTriangle(v2,v4,v5));
-// fillTris.push_back(new MTriangle(v4,c,v5));
-// fillTris.push_back(new MTriangle(v0,v3,v1));
-
- }
+ //loopfillTris.push_back(new MTriangle(v0,v1, c));
+
+ // MVertex *v2 = new MVertex(.5*(v0->x()+c->x()), .5*(v0->y()+c->y()), .5*(v0->z()+c->z()));
+ // MVertex *v3 = new MVertex(.5*(v1->x()+c->x()), .5*(v1->y()+c->y()), .5*(v1->z()+c->z()));
+ // fillNodes.insert(c); fillNodes.insert(v2); fillNodes.insert(v3);
+ // loopfillTris.push_back(new MTriangle(v0,v2,v3));
+ // loopfillTris.push_back(new MTriangle(v2,c, v3));
+ // loopfillTris.push_back(new MTriangle(v0,v3, v1));
+
+ double k = 1/3.; double kk = 2/3.;
+ MVertex *v2 = new MVertex(kk*v0->x()+k*c->x(), kk*v0->y()+k*c->y(),kk*v0->z()+k*c->z());
+ MVertex *v3 = new MVertex(kk*v1->x()+k*c->x(), kk*v1->y()+k*c->y(),kk*v1->z()+k*c->z());
+ MVertex *v4 = new MVertex(k*v0->x()+kk*c->x(), k*v0->y()+kk*c->y(),k*v0->z()+kk*c->z());
+ MVertex *v5 = new MVertex(k*v1->x()+kk*c->x(), k*v1->y()+kk*c->y(),k*v1->z()+kk*c->z());
+ fillNodes.insert(c); fillNodes.insert(v2); fillNodes.insert(v3);
+ fillNodes.insert(v4); fillNodes.insert(v5);
+ loopfillTris.push_back(new MTriangle(v0,v2,v3));
+ loopfillTris.push_back(new MTriangle(v2,v5,v3));
+ loopfillTris.push_back(new MTriangle(v2,v4,v5));
+ loopfillTris.push_back(new MTriangle(v4,c,v5));
+ loopfillTris.push_back(new MTriangle(v0,v3,v1));
}
+
}
- }
-
- if (fillTris.size() > 0){
- FILE * ftri = fopen("fillTris.pos","a");
- fprintf(ftri,"View \"\"{\n");
- for (std::list<MTriangle*>::iterator it2 = fillTris.begin(); it2 !=fillTris.end(); it2++ ){
- MTriangle *t = (*it2);
- fprintf(ftri,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
- t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z(),
- t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z(),
- t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->z(),
- 1., 1., 1.);
+ //check normal orientations of loopfillTris
+ bool invertTris = false;
+ std::map<MEdge, std::set<MTriangle*>, Less_Edge > edge2tris;
+ for(std::list<MTriangle*>::iterator t= loopfillTris.begin(); t!=loopfillTris.end(); t++){
+ for (int j = 0; j < 3; j++){
+ MEdge me = (*t)->getEdge(j);
+ std::map<MEdge, std::set<MTriangle*, std::less<MTriangle*> >, Less_Edge >::iterator it = edge2tris.find(me);
+ if (it == edge2tris.end()) {
+ std::set<MTriangle*, std::less<MTriangle*> > mySet;
+ mySet.insert(*t);
+ edge2tris.insert(std::make_pair(me, mySet));
+ }
+ else{
+ std::set<MTriangle*, std::less<MTriangle*> > mySet = it->second;
+ mySet.insert(*t);
+ it->second = mySet;
+ }
+ }
+ }
+ int iE, si, iE2, si2;
+ std::list<GFace*>::const_iterator itf = _compound.begin();
+ for( ; itf != _compound.end(); ++itf){
+ for(unsigned int i = 0; i < (*itf)->triangles.size(); ++i){
+ MTriangle *t = (*itf)->triangles[i];
+ for (int j = 0; j < 3; j++){
+ MEdge me = t->getEdge(j);
+ std::map<MEdge, std::set<MTriangle*>, Less_Edge >::iterator it = edge2tris.find(me);
+ if(it != edge2tris.end()){
+ t->getEdgeInfo(me, iE,si);
+ MTriangle* t2 = *((it->second).begin());
+ t2->getEdgeInfo(me,iE2,si2);
+ if(si == si2) {
+ invertTris = true;
+ break;
+ }
+ }
+ }
+ }
+ }
+ if (invertTris){
+ for (std::list<MTriangle*>::iterator it = loopfillTris.begin(); it !=loopfillTris.end(); it++ )
+ (*it)->revert();
+ }
+
+ fillTris.insert(fillTris.begin(),loopfillTris.begin(),loopfillTris.end());
+
+ }
+
+ //printing
+ if( !CTX::instance()->mesh.saveAll){
+ if (fillTris.size() > 0){
+ char name[256];
+ std::list<GFace*>::const_iterator itf = _compound.begin();
+ sprintf(name, "fillTris-%d.pos", (*itf)->tag());
+ FILE * ftri = fopen(name,"w");
+ fprintf(ftri,"View \"\"{\n");
+ for (std::list<MTriangle*>::iterator it2 = fillTris.begin(); it2 !=fillTris.end(); it2++ ){
+ MTriangle *t = (*it2);
+ fprintf(ftri,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
+ t->getVertex(0)->x(), t->getVertex(0)->y(), t->getVertex(0)->z(),
+ t->getVertex(1)->x(), t->getVertex(1)->y(), t->getVertex(1)->z(),
+ t->getVertex(2)->x(), t->getVertex(2)->y(), t->getVertex(2)->z(),
+ 1., 1., 1.);
+ }
+ fprintf(ftri,"};\n");
+ fclose(ftri);
}
- fprintf(ftri,"};\n");
- fclose(ftri);
}
-
}
@@ -358,31 +479,42 @@ bool GFaceCompound::trivial() const
}
-//For the conformal map the linear system cannot guarantee there is no folding
+//For the conformal map the linear system cannot guarantee there is no overlapping
//of triangles
-bool GFaceCompound::checkFolding(std::vector<MVertex*> &ordered) const
+bool GFaceCompound::checkOverlap(std::vector<MVertex*> &ordered) const
{
- bool has_no_folding = true;
-
- for(unsigned int i = 0; i < ordered.size()-1; ++i){
- SPoint3 p1 = coordinates[ordered[i]];
- SPoint3 p2 = coordinates[ordered[i+1]];
- int maxSize = (i==0) ? ordered.size()-2: ordered.size()-1;
- for(int k = i+2; k < maxSize; ++k){
- SPoint3 q1 = coordinates[ordered[k]];
- SPoint3 q2 = coordinates[ordered[k+1]];
- double x[2];
- int inters = intersection_segments (p1,p2,q1,q2,x);
- if (inters > 0) has_no_folding = false;
+ bool has_no_overlap = true;
+
+ for(std::list<std::list<GEdge*> >::const_iterator iloop = _interior_loops.begin();
+ iloop != _interior_loops.end(); iloop++){
+ std::vector<MVertex*> ordered;
+ std::vector<double> coords;
+ bool success = orderVertices(*iloop, ordered, coords);
+
+ for(unsigned int i = 0; i < ordered.size()-1; ++i){
+ SPoint3 p1 = coordinates[ordered[i]];
+ SPoint3 p2 = coordinates[ordered[i+1]];
+ int maxSize = (i==0) ? ordered.size()-2: ordered.size()-1;
+ for(int k = i+2; k < maxSize; ++k){
+ SPoint3 q1 = coordinates[ordered[k]];
+ SPoint3 q2 = coordinates[ordered[k+1]];
+ double x[2];
+ int inters = intersection_segments (p1,p2,q1,q2,x);
+ if (inters > 0){
+ has_no_overlap = false;
+ break;
+ }
+ }
}
+
}
- if ( !has_no_folding ) {
- Msg::Warning("$$$ Folding for compound face %d", this->tag());
+ if ( !has_no_overlap ) {
+ Msg::Warning("$$$ Overlap for compound face %d", this->tag());
}
- return has_no_folding;
+ return has_no_overlap;
}
@@ -424,15 +556,15 @@ bool GFaceCompound::checkOrientation(int iter) const
}
}
- int iterMax = 10;
+ int iterMax = 15;
if(!oriented && iter < iterMax){
- if (iter == 0) Msg::Warning("--- Parametrization is flipping : applying cavity checks.");
+ if (iter == 0) Msg::Info("--- Parametrization is flipping : applying cavity checks.");
Msg::Debug("--- Cavity Check - iter %d -",iter);
one2OneMap();
return checkOrientation(iter+1);
}
else if (oriented && iter < iterMax){
- Msg::Info("Parametrization has no flips :-)");
+ //Msg::Info("Parametrization is bijective (no flips)");
//printStuff();
}
@@ -518,11 +650,11 @@ bool GFaceCompound::parametrize() const
else if (_mapping == CONFORMAL){
Msg::Debug("Parametrizing surface %d with 'conformal map'", tag());
fillNeumannBCS();
- bool withoutFolding = parametrize_conformal_spectral() ;
- //bool withoutFolding = parametrize_conformal();
- if ( withoutFolding == false ){
+ bool withoutOverlap = parametrize_conformal_spectral() ;
+ //bool withoutOverlap = parametrize_conformal();
+ if ( withoutOverlap == false || !checkOrientation(0) ){
Msg::Warning("$$$ Parametrization switched to harmonic map");
- parametrize(ITERU,HARMONIC);
+ parametrize(ITERU,HARMONIC);
parametrize(ITERV,HARMONIC);
//buildOct(); exit(1);
}
@@ -531,7 +663,7 @@ bool GFaceCompound::parametrize() const
buildOct();
if (!checkOrientation(0)){
- Msg::Info("--- Parametrization switched to convex combination map");
+ Msg::Info("### Parametrization switched to convex combination map");
coordinates.clear();
Octree_Delete(oct);
fillNeumannBCS();
@@ -544,8 +676,7 @@ bool GFaceCompound::parametrize() const
double AR = checkAspectRatio();
if (floor(AR) > AR_MAX){
- Msg::Warning("Geometrical aspect ratio too high AR=%d ", (int)AR);
- //exit(1);
+ Msg::Warning("Geometrical aspect ratio is high AR=%d ", (int)AR);
paramOK = true; //false;
}
@@ -834,68 +965,6 @@ GFaceCompound::~GFaceCompound()
if (_lsys)delete _lsys;
}
-// order vertices of a closed loop
-static bool orderVertices(const std::list<GEdge*> &e, std::vector<MVertex*> &l,
- std::vector<double> &coord)
-{
- l.clear();
- coord.clear();
-
- std::list<GEdge*>::const_iterator it = e.begin();
- std::list<MLine*> temp;
- double tot_length = 0;
- for( ; it != e.end(); ++it ){
- for(unsigned int i = 0; i < (*it)->lines.size(); i++ ){
- temp.push_back((*it)->lines[i]);
- MVertex *v0 = (*it)->lines[i]->getVertex(0);
- MVertex *v1 = (*it)->lines[i]->getVertex(1);
- const double length = sqrt((v0->x() - v1->x()) * (v0->x() - v1->x()) +
- (v0->y() - v1->y()) * (v0->y() - v1->y()) +
- (v0->z() - v1->z()) * (v0->z() - v1->z()));
- tot_length += length;
- }
- }
-
- MVertex *first_v = (*temp.begin())->getVertex(0);
- MVertex *current_v = (*temp.begin())->getVertex(1);
-
- l.push_back(first_v);
- coord.push_back(0.0);
- temp.erase(temp.begin());
-
- while(temp.size()){
- bool found = false;
- for(std::list<MLine*>::iterator itl = temp.begin(); itl != temp.end(); ++itl){
- MLine *ll = *itl;
- MVertex *v0 = ll->getVertex(0);
- MVertex *v1 = ll->getVertex(1);
- if(v0 == current_v){
- found = true;
- l.push_back(current_v);
- current_v = v1;
- temp.erase(itl);
- const double length = sqrt((v0->x() - v1->x()) * (v0->x() - v1->x()) +
- (v0->y() - v1->y()) * (v0->y() - v1->y()) +
- (v0->z() - v1->z()) * (v0->z() - v1->z()));
- coord.push_back(coord[coord.size()-1] + length / tot_length);
- break;
- }
- else if(v1 == current_v){
- found = true;
- l.push_back(current_v);
- current_v = v0;
- temp.erase(itl);
- const double length = sqrt((v0->x() - v1->x()) * (v0->x() - v1->x()) +
- (v0->y() - v1->y()) * (v0->y() - v1->y()) +
- (v0->z() - v1->z()) * (v0->z() - v1->z()));
- coord.push_back(coord[coord.size()-1] + length / tot_length);
- break;
- }
- }
- if(!found) return false;
- }
- return true;
-}
SPoint2 GFaceCompound::getCoordinates(MVertex *v) const
{
@@ -1058,7 +1127,7 @@ void GFaceCompound::parametrize(iterationStep step, typeOfMapping tom, double al
}
}
- for (std::list<MTriangle*>::iterator it2 = fillTris.begin(); it2 !=fillTris.end(); it2++ ){
+ for (std::list<MTriangle*>::iterator it2 = fillTris.begin(); it2 !=fillTris.end(); it2++ ){
MTriangle *t = (*it2);
myAssembler.numberVertex(t->getVertex(0), 0, 1);
myAssembler.numberVertex(t->getVertex(1), 0, 1);
@@ -1160,6 +1229,12 @@ bool GFaceCompound::parametrize_conformal_spectral() const
myAssembler.numberVertex(v, 0, 2);
}
+ for(std::set<MVertex *>::iterator itv = fillNodes.begin(); itv !=fillNodes.end() ; ++itv){
+ MVertex *v = *itv;
+ myAssembler.numberVertex(v, 0, 1);
+ myAssembler.numberVertex(v, 0, 2);
+ }
+
simpleFunction<double> ONE(1.0);
simpleFunction<double> MONE(-1.0 );
laplaceTerm laplace1(model(), 1, &ONE);
@@ -1176,6 +1251,15 @@ bool GFaceCompound::parametrize_conformal_spectral() const
cross21.addToMatrix(myAssembler, &se);
}
}
+
+ for (std::list<MTriangle*>::iterator it2 = fillTris.begin(); it2 !=fillTris.end(); it2++ ){
+ SElement se((*it2));
+ laplace1.addToMatrix(myAssembler, &se);
+ laplace2.addToMatrix(myAssembler, &se);
+ cross12.addToMatrix(myAssembler, &se);
+ cross21.addToMatrix(myAssembler, &se);
+ }
+
double epsilon = 1.e-7;
for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
MVertex *v = *itv;
@@ -1184,6 +1268,13 @@ bool GFaceCompound::parametrize_conformal_spectral() const
myAssembler.assemble(v, 0, 2, v, 0, 2, epsilon);
}
}
+ for(std::set<MVertex *>::iterator itv = fillNodes.begin(); itv !=fillNodes.end() ; ++itv){
+ MVertex *v = *itv;
+ if (std::find(ordered.begin(), ordered.end(), v) == ordered.end() ){
+ myAssembler.assemble(v, 0, 1, v, 0, 1, epsilon);
+ myAssembler.assemble(v, 0, 2, v, 0, 2, epsilon);
+ }
+ }
//-------------------------------
myAssembler.setCurrentMatrix("B");
@@ -1200,6 +1291,11 @@ bool GFaceCompound::parametrize_conformal_spectral() const
myAssembler.assemble(v, 0, 2, v, 0, 2, 1.0);
}
}
+ for(std::set<MVertex *>::iterator itv = fillNodes.begin(); itv !=fillNodes.end() ; ++itv){
+ MVertex *v = *itv;
+ myAssembler.assemble(v, 0, 1, v, 0, 1, small);
+ myAssembler.assemble(v, 0, 2, v, 0, 2, small);
+ }
//mettre max NC contraintes par bord
int NB = ordered.size();
@@ -1243,7 +1339,7 @@ bool GFaceCompound::parametrize_conformal_spectral() const
lsysA->clear();
lsysB->clear();
- return checkFolding(ordered);
+ return checkOverlap(ordered);
}
else return false;
@@ -1344,8 +1440,8 @@ bool GFaceCompound::parametrize_conformal() const
_lsys->clear();
- //check for folding
- return checkFolding(ordered);
+ //check for overlapping triangles
+ return checkOverlap(ordered);
}
@@ -1840,15 +1936,20 @@ bool GFaceCompound::checkTopology() const
Msg::Info("--- Split surface %d in %d parts with Multilevel Mesh partitioner", tag(), nbSplit);
}
}
- else if (G == 0 && AR > AR_MAX){
- correctTopo = false;
- nbSplit = -2;
- Msg::Warning("Wrong topology: Genus=%d, Nb boundaries=%d, AR=%d", G, Nb, AR);
- if (_allowPartition){
- Msg::Info("-----------------------------------------------------------");
- Msg::Info("--- Split surface %d in 2 parts with Laplacian Mesh partitioner", tag());
- }
- }
+ else if (G == 0 && AR > AR_MAX){
+ correctTopo = false;
+ Msg::Warning("Wrong topology: Genus=%d, Nb boundaries=%d, AR=%d", G, Nb, AR);
+ if (_allowPartition == 1){
+ nbSplit = -2;
+ Msg::Info("-----------------------------------------------------------");
+ Msg::Info("--- Split surface %d in 2 parts with Laplacian Mesh partitioner", tag());
+ }
+ else if (_allowPartition == 2){
+ nbSplit = 2;
+ Msg::Info("-----------------------------------------------------------");
+ Msg::Info("--- Split surface %d in %d parts with Multilevel Mesh partitioner", tag(), nbSplit);
+ }
+ }
else{
Msg::Debug("Correct topology: Genus=%d and Nb boundaries=%d, AR=%g", G, Nb, H/D);
}
@@ -1865,7 +1966,7 @@ double GFaceCompound::checkAspectRatio() const
if(allNodes.empty()) buildAllNodes();
- double limit = 1.e-17;
+ double limit = 1.e-20;
double areaMin = 1.e20;
double area3D = 0.0;
int nb = 0;
@@ -2029,9 +2130,7 @@ void GFaceCompound::printStuff() const
sprintf(name5, "XYZV-%d.pos", (*it)->tag());
sprintf(name6, "XYZC-%d.pos", (*it)->tag());
- sprintf(name7, "UVME-%d.pos", (*it)->tag());
- sprintf(name8, "UVMF-%d.pos", (*it)->tag());
- sprintf(name9, "UVMG-%d.pos", (*it)->tag());
+ sprintf(name7, "UVM-%d.pos", (*it)->tag());
FILE * uva = fopen(name0,"w");
FILE * uvx = fopen(name1,"w");
@@ -2040,9 +2139,7 @@ void GFaceCompound::printStuff() const
FILE * xyzu = fopen(name4,"w");
FILE * xyzv = fopen(name5,"w");
FILE * xyzc = fopen(name6,"w");
- FILE * uvme = fopen(name7,"w");
- FILE * uvmf = fopen(name8,"w");
- FILE * uvmg = fopen(name9,"w");
+ FILE * uvm = fopen(name7,"w");
fprintf(uva,"View \"\"{\n");
@@ -2052,10 +2149,7 @@ void GFaceCompound::printStuff() const
fprintf(xyzu,"View \"\"{\n");
fprintf(xyzv,"View \"\"{\n");
fprintf(xyzc,"View \"\"{\n");
- fprintf(uvme,"View \"\"{\n");
- fprintf(uvmf,"View \"\"{\n");
- fprintf(uvmg,"View \"\"{\n");
-
+ fprintf(uvm,"View \"\"{\n");
for( ; it != _compound.end() ; ++it){
for(unsigned int i = 0; i < (*it)->triangles.size(); ++i){
@@ -2108,27 +2202,35 @@ void GFaceCompound::printStuff() const
double metric2f = dot(der2.second()*(1./norm(der2.second())), der2.first()*(1./norm(der2.first())));
double metric2g = dot(der2.second(), der2.second());
+ double mat0[2][2], eig0[2];
+ double mat1[2][2], eig1[2];
+ double mat2[2][2], eig2[2];
+ mat0[0][0] = metric0e; mat0[0][1] = metric0f;
+ mat0[1][0] = metric0f; mat0[1][1] = metric0g;
+ eigenvalue2x2(mat0, eig0);
+ mat1[0][0] = metric1e; mat1[0][1] = metric1f;
+ mat1[1][0] = metric1f; mat1[1][1] = metric1g;
+ eigenvalue2x2(mat1, eig1);
+ mat2[0][0] = metric2e; mat2[0][1] = metric2f;
+ mat2[1][0] = metric2f; mat2[1][1] = metric2g;
+ eigenvalue2x2(mat2, eig2);
+
+ double disp0 = sqrt(.5*(eig0[0]*eig0[0]+ (eig0[1]*eig0[1])));
+ double disp1 = sqrt(.5*(eig1[0]*eig1[0]+ (eig1[1]*eig1[1])));
+ double disp2 = sqrt(.5*(eig2[0]*eig2[0]+ (eig2[1]*eig2[1])));
+ double mdisp = .333*(disp0+disp1+disp2);
+
fprintf(uva,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
it0->second.x(), it0->second.y(), 0.0,
it1->second.x(), it1->second.y(), 0.0,
it2->second.x(), it2->second.y(), 0.0,
area, area, area);
- fprintf(uvme,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
+ fprintf(uvm,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
it0->second.x(), it0->second.y(), 0.0,
it1->second.x(), it1->second.y(), 0.0,
it2->second.x(), it2->second.y(), 0.0,
- 0.3*(metric0e+metric1e+metric2e), 0.3*(metric0e+metric1e+metric2e), 0.3*(metric0e+metric1e+metric2e) );
- fprintf(uvmf,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
- it0->second.x(), it0->second.y(), 0.0,
- it1->second.x(), it1->second.y(), 0.0,
- it2->second.x(), it2->second.y(), 0.0,
- 0.3*(metric0f+metric1f+metric2f), 0.3*(metric0f+metric1f+metric2f), 0.3*(metric0f+metric1f+metric2f) );
- fprintf(uvmg,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g};\n",
- it0->second.x(), it0->second.y(), 0.0,
- it1->second.x(), it1->second.y(), 0.0,
- it2->second.x(), it2->second.y(), 0.0,
- 0.3*(metric0g+metric1g+metric2g), 0.3*(metric0g+metric1g+metric2g), 0.3*(metric0g+metric1g+metric2g) );
+ mdisp, mdisp, mdisp);
fprintf(uvx,"ST(%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E){%22.15E,%22.15E,%22.15E};\n",
it0->second.x(), it0->second.y(), 0.0,
@@ -2161,48 +2263,8 @@ void GFaceCompound::printStuff() const
fclose(xyzv);
fprintf(xyzc,"};\n");
fclose(xyzc);
- fprintf(uvme,"};\n");
- fclose(uvme);
- fprintf(uvmf,"};\n");
- fclose(uvmf);
- fprintf(uvmg,"};\n");
- fclose(uvmg);
-
-
- //debug cecile rbf
- // it = _compound.begin();
- // char nameM[256], nameF[256];
- // sprintf(nameM, "mappedMesh-%d.msh", (*it)->tag());
- // sprintf(nameF, "XYZfunction-%d.txt", (*it)->tag());
- // FILE * myF = fopen(nameM,"w");
- // FILE * myF2 = fopen(nameF,"w");
- // fprintf(myF,"$MeshFormat\n");
- // fprintf(myF,"2.2 0 8\n");
- // fprintf(myF,"$EndMeshFormat\n");
- // fprintf(myF,"$Nodes\n");
- // fprintf(myF,"%d\n", (int)allNodes.size());
- // for(std::set<MVertex *>::iterator itv = allNodes.begin(); itv !=allNodes.end() ; ++itv){
- // std::map<MVertex*,SPoint3>::const_iterator it0 = coordinates.find(*itv);
- // fprintf(myF,"%d %g %g %g \n", (*itv)->getNum(), it0->second.x(), it0->second.y(), 0.0);
- // fprintf(myF2,"%d %g %g %g \n", (*itv)->getNum(), (*itv)->x(), (*itv)->y(), (*itv)->z());
- // }
- // fprintf(myF,"$EndNodes\n");
- // fprintf(myF,"$Elements\n");
- // int nbTris = 0;
- // for( ; it != _compound.end() ; ++it) nbTris += (*it)->triangles.size();
- // fprintf(myF, "%d \n", nbTris);
- // int k = 1;
- // for(it = _compound.begin(); it != _compound.end() ; ++it){
- // for(unsigned int i = 0; i < (*it)->triangles.size(); ++i){
- // MTriangle *t = (*it)->triangles[i];
- // fprintf(myF,"%d 2 2 0 1 %d %d %d \n", k, t->getVertex(0)->getNum(), t->getVertex(1)->getNum(), t->getVertex(2)->getNum());
- // k++;
- // }
- // }
- // fprintf(myF,"$EndElements\n");
- // fclose(myF);
- // fclose(myF2);
-
+ fprintf(uvm,"};\n");
+ fclose(uvm);;
}
diff --git a/Geo/GFaceCompound.h b/Geo/GFaceCompound.h
index 1e43c0fd7cbbb6aebd2b076161dcc484cacfb1f9..7f4c347f01dc240628901d19332109e7a30cce4f 100644
--- a/Geo/GFaceCompound.h
+++ b/Geo/GFaceCompound.h
@@ -70,6 +70,7 @@ class GFaceCompound : public GFace {
mutable std::map<SPoint3,SPoint3 > _coordPoints;
mutable std::map<MVertex*, SVector3> _normals;
mutable std::list<MTriangle*> fillTris;
+ mutable std::set<MVertex*> fillNodes;
void buildOct() const ;
void buildAllNodes() const;
void parametrize(iterationStep, typeOfMapping, double alpha=0.) const;
@@ -77,7 +78,7 @@ class GFaceCompound : public GFace {
bool parametrize_conformal_spectral() const;
void compute_distance() const;
bool checkOrientation(int iter) const;
- bool checkFolding(std::vector<MVertex*> &ordered) const;
+ bool checkOverlap(std::vector<MVertex*> &ordered) const;
void one2OneMap() const;
double checkAspectRatio() const;
void computeNormals () const;
diff --git a/Geo/GModel.h b/Geo/GModel.h
index bfc416c9fd01fd06c2488f0b276e8e38714cffb0..35b9a7131df97c77c069dd2247c795664a99dff7 100644
--- a/Geo/GModel.h
+++ b/Geo/GModel.h
@@ -488,6 +488,10 @@ class GModel
int writeSTL(const std::string &name, bool binary=false,
bool saveAll=false, double scalingFactor=1.0);
+ // PLY2 format (asciii text format)
+ int readPLY2(const std::string &name);
+ int writePLY2(const std::string &name);
+
// Inventor/VRML format
int readVRML(const std::string &name);
int writeVRML(const std::string &name,
diff --git a/Geo/GModelIO_Geo.cpp b/Geo/GModelIO_Geo.cpp
index 8c5d0c15a6665ad62f1786363a1eb2cdcee84573..f9c2c42ab997d7ad3f10306cb659e5920e348dc5 100644
--- a/Geo/GModelIO_Geo.cpp
+++ b/Geo/GModelIO_Geo.cpp
@@ -180,8 +180,9 @@ int GModel::importGEOInternals()
if(ge) b[j].push_back(ge);
}
}
- int allowPartition = 1;
- if (abs(s->TypeOfMapping) != 1) allowPartition = 0;
+ int allowPartition = 1; //allow
+ if (abs(s->TypeOfMapping) == 2) allowPartition = 0;//not allow
+ if (abs(s->TypeOfMapping) == 3) allowPartition = 2;//only metis
f = new GFaceCompound(this, s->Num, comp,
b[0], b[1], b[2], b[3], 0,
diff --git a/Geo/GModelIO_Mesh.cpp b/Geo/GModelIO_Mesh.cpp
index 138dd33a0200da3978c192792b5b4d85a5339e45..ebdaa804e5ee63dde13ebc672dd91f8d61556482 100644
--- a/Geo/GModelIO_Mesh.cpp
+++ b/Geo/GModelIO_Mesh.cpp
@@ -1208,6 +1208,69 @@ static int readElementsVRML(FILE *fp, std::vector<MVertex*> &vertexVector, int r
return 1;
}
+int GModel::readPLY2(const std::string &name)
+{
+ FILE *fp = fopen(name.c_str(), "r");
+ if(!fp){
+ Msg::Error("Unable to open file '%s'", name.c_str());
+ return 0;
+ }
+
+ std::vector<MVertex*> vertexVector;
+ std::map<int, std::vector<MElement*> > elements[5];
+
+ char buffer[256], str[256];
+ int elementary = getMaxElementaryNumber(-1) + 1;
+ int nbv, nbf;
+ while(!feof(fp)) {
+ if(!fgets(buffer, sizeof(buffer), fp)) break;
+ if(buffer[0] != '#'){ // skip comments
+ sscanf(buffer, "%d", &nbv);
+ if(!fgets(buffer, sizeof(buffer), fp)) break;
+ sscanf(buffer, "%d", &nbf);
+ Msg::Info("%d vertices", nbv);
+ Msg::Info("%d triangles", nbf);
+ vertexVector.resize(nbv);
+ for(int i = 0; i < nbv; i++) {
+ if(!fgets(buffer, sizeof(buffer), fp)) break;
+ double x, y, z;
+ int nb = sscanf(buffer, "%lf %lf %lf", &x, &y, &z);
+ if (nb !=3){
+ if(!fgets(buffer, sizeof(buffer), fp)) break;
+ sscanf(buffer, "%lf", &y);
+ if(!fgets(buffer, sizeof(buffer), fp)) break;
+ sscanf(buffer, "%lf", &z);
+ }
+ vertexVector[i] = new MVertex(x, y, z);
+ }
+ for(int i = 0; i < nbf; i++) {
+ if(!fgets(buffer, sizeof(buffer), fp)) break;
+ int n[3], nbe;
+ int nb = sscanf(buffer, "%d %d %d %d", &nbe, &n[0], &n[1], &n[2]);
+ if (nb !=4){
+ if(!fgets(buffer, sizeof(buffer), fp)) break;
+ sscanf(buffer, "%d", &n[0]);
+ if(!fgets(buffer, sizeof(buffer), fp)) break;
+ sscanf(buffer, "%d", &n[1]);
+ if(!fgets(buffer, sizeof(buffer), fp)) break;
+ sscanf(buffer, "%d", &n[2]);
+ }
+ std::vector<MVertex*> vertices;
+ if(!getVertices(3, n, vertexVector, vertices)) return 0;
+ elements[0][elementary].push_back(new MTriangle(vertices));
+ }
+ }
+ }
+
+ for(int i = 0; i < (int)(sizeof(elements) / sizeof(elements[0])); i++)
+ _storeElementsInEntities(elements[i]);
+ _associateEntityWithMeshVertices();
+ _storeVerticesInEntities(vertexVector);
+
+ fclose(fp);
+ return 1;
+}
+
int GModel::readVRML(const std::string &name)
{
FILE *fp = fopen(name.c_str(), "r");
@@ -1279,6 +1342,38 @@ int GModel::readVRML(const std::string &name)
return 1;
}
+int GModel::writePLY2(const std::string &name)
+{
+ FILE *fp = fopen(name.c_str(), "w");
+ if(!fp){
+ Msg::Error("Unable to open file '%s'", name.c_str());
+ return 0;
+ }
+
+ int numVertices = indexMeshVertices(true);
+ int numTriangles = 0.0;
+ for(fiter it = firstFace(); it != lastFace(); ++it){
+ numTriangles += (*it)->triangles.size();
+ }
+
+ fprintf(fp, "%d\n", numVertices);
+ fprintf(fp, "%d\n", numTriangles);
+
+ std::vector<GEntity*> entities;
+ getEntities(entities);
+ for(unsigned int i = 0; i < entities.size(); i++)
+ for(unsigned int j = 0; j < entities[i]->mesh_vertices.size(); j++)
+ entities[i]->mesh_vertices[j]->writePLY2(fp);
+
+ for(fiter it = firstFace(); it != lastFace(); ++it){
+ for(unsigned int i = 0; i < (*it)->triangles.size(); i++)
+ (*it)->triangles[i]->writePLY2(fp);
+ }
+
+ fclose(fp);
+ return 1;
+}
+
int GModel::writeVRML(const std::string &name, bool saveAll, double scalingFactor)
{
FILE *fp = fopen(name.c_str(), "w");
diff --git a/Geo/Geo.h b/Geo/Geo.h
index fe5d0bc2ac391a93f88b8692c07adbc8056489d3..c40272fd46a483fb9ffa55165460a2a2a2727762 100644
--- a/Geo/Geo.h
+++ b/Geo/Geo.h
@@ -148,8 +148,9 @@ class Surface{
int Method;
int Recombine;
int Recombine_Dir; // -1 is left, +1 is right, 0 is alternated
- int TypeOfMapping; // +1 is harmonic, -1 is conformal, +2 harmonic_NoSplit,
- // -2 conformal_NoSplit
+ int TypeOfMapping; // +1 is harmonic, -1 is conformal,
+ //+2 harmonic_NoSplit, -2 conformal_NoSplit,
+ //+3 harmonic_SplitMetis, -3 conformal_SplitMetis,
double RecombineAngle;
int TransfiniteSmoothing;
List_T *Generatrices;
diff --git a/Geo/MElement.cpp b/Geo/MElement.cpp
index 6bea2b0079bfc37d6b8be3893a9905b58ec79ce1..48ee94d70d2205d924087dc043f7f98433cb584b 100644
--- a/Geo/MElement.cpp
+++ b/Geo/MElement.cpp
@@ -646,6 +646,14 @@ void MElement::writeSTL(FILE *fp, bool binary, double scalingFactor)
}
}
}
+void MElement::writePLY2(FILE *fp)
+{
+ setVolumePositive();
+ fprintf(fp, "3 ");
+ for(int i = 0; i < getNumVertices(); i++)
+ fprintf(fp, " %d", getVertex(i)->getIndex() - 1);
+ fprintf(fp, "\n");
+}
void MElement::writeVRML(FILE *fp)
{
diff --git a/Geo/MElement.h b/Geo/MElement.h
index 7f0992c0af1877fe9c8420f19cc4d8e92b6e3360..4b0bcb558599130b7cb7037c559309a3483e1373 100644
--- a/Geo/MElement.h
+++ b/Geo/MElement.h
@@ -274,6 +274,7 @@ class MElement
bool printDisto,double scalingFactor=1.0, int elementary=1);
virtual void writeSTL(FILE *fp, bool binary=false, double scalingFactor=1.0);
virtual void writeVRML(FILE *fp);
+ virtual void writePLY2(FILE *fp);
virtual void writeUNV(FILE *fp, int num=0, int elementary=1, int physical=1);
virtual void writeVTK(FILE *fp, bool binary=false, bool bigEndian=false);
virtual void writeMESH(FILE *fp, int elementTagType=1, int elementary=1,
diff --git a/Geo/MVertex.cpp b/Geo/MVertex.cpp
index ac1a6a256bd377dbe70c4c1140a8619e18306485..bfdedec319d9ee6bfbea6a8ddd4e9a1e8cc8ef93 100644
--- a/Geo/MVertex.cpp
+++ b/Geo/MVertex.cpp
@@ -108,6 +108,13 @@ void MVertex::writeMSH(FILE *fp, bool binary, bool saveParametric, double scalin
fprintf(fp, "\n");
}
}
+void MVertex::writePLY2(FILE *fp)
+{
+ if(_index < 0) return; // negative index vertices are never saved
+
+ fprintf(fp, "%.16g %.16g %.16g\n",
+ x(), y(), z());
+}
void MVertex::writeVRML(FILE *fp, double scalingFactor)
{
diff --git a/Geo/MVertex.h b/Geo/MVertex.h
index 40c8d4fd934d8e875bcf8475e066226f8e340d53..a568d7aa7eafe3a7e278386904abc240a1c54cb7 100644
--- a/Geo/MVertex.h
+++ b/Geo/MVertex.h
@@ -107,6 +107,7 @@ class MVertex{
// IO routines
void writeMSH(FILE *fp, bool binary=false, bool saveParametric=false,
double scalingFactor=1.0);
+ void writePLY2(FILE *fp);
void writeVRML(FILE *fp, double scalingFactor=1.0);
void writeUNV(FILE *fp, double scalingFactor=1.0);
void writeVTK(FILE *fp, bool binary=false, double scalingFactor=1.0,
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index fa25e3e3625f94867255b15205a526a8c78d934a..fc51040902c5d9e3fd655a882b0caa3008d66348 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -1624,9 +1624,9 @@ void partitionAndRemesh(GFaceCompound *gf)
if(gf->nbSplit > 0) method = MULTILEVEL;
else method = LAPLACIAN;
- multiscalePartition *msp = new multiscalePartition(elements, abs(gf->nbSplit), method);
-
- //gf->partitionFaceCM();
+ int allowType = gf->allowPartition();
+ multiscalePartition *msp = new multiscalePartition(elements, abs(gf->nbSplit),
+ method, allowType);
int NF = msp->getNumberOfParts();
int numv = gf->model()->getMaxElementaryNumber(0) + 1;
@@ -1781,8 +1781,8 @@ void partitionAndRemesh(GFaceCompound *gf)
}
double t3 = Cpu();
- Msg::Info("*** Mesh of surface %d done by assembly remeshed faces (%g s)",
- gf->tag(), t3-t2);
+ Msg::Info("*** Mesh of surface %d done by assembly %d remeshed faces (%g s)",
+ gf->tag(), NF, t3-t2);
Msg::Info("-----------------------------------------------------------");
gf->coherenceNormals();
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index de64974bef8706185b35beb1295b37d5b98977d0..06005de6fad87f41dc1828c1c255b136181a36d4 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -26,7 +26,99 @@
#include "ANN/ANN.h"
#endif
-void voronoiDual(GRegion *gr)
+void printVoronoi(GRegion *gr)
+{
+ 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;
+ 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 voronoi nodes
+ 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);
+ //}//uncomment this
+
+ }
+ fprintf(outfile,"};\n");
+ fclose(outfile);
+
+ //print scalar lines
+ FILE *outfile2;
+ outfile2 = fopen("edges.pos", "w");
+ 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;
+ 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());
+ }
+ fprintf(outfile2,"};\n");
+ fclose(outfile2);
+
+}
+
+void skeletonFromVoronoi(GRegion *gr)
{
std::vector<MTetrahedron*> elements = gr->tetrahedra;
@@ -50,10 +142,10 @@ void voronoiDual(GRegion *gr)
}
printf("Dmax =%g \n", Dmax);
- printf("printing nodes \n");
+ printf("printing skeleton nodes \n");
FILE *outfile;
- outfile = fopen("nodes.pos", "w");
- fprintf(outfile, "View \"voronoi nodes\" {\n");
+ outfile = fopen("skeletonNodes.pos", "w");
+ fprintf(outfile, "View \"skeleton nodes\" {\n");
for(unsigned int i = 0; i < elements.size(); i++){
MTetrahedron *ele = elements[i];
SPoint3 pc = ele->circumcenter();
@@ -62,7 +154,7 @@ void voronoiDual(GRegion *gr)
ele->xyz2uvw(x, uvw);
if(ele->isInside(uvw[0], uvw[1], uvw[2])){
double radius = ele->getCircumRadius();
- if(radius > Dmax/5.) {
+ if(radius > Dmax/10.) {
candidates.insert(pc);
fprintf(outfile,"SP(%g,%g,%g) {%g};\n",
pc.x(), pc.y(), pc.z(), radius);
@@ -72,151 +164,76 @@ void voronoiDual(GRegion *gr)
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);
+ FILE *outfile2;
+ outfile2 = fopen("skeletonEdges.pos", "w");
+ fprintf(outfile2, "View \"skeleton 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;
+ std::set<SPoint3>::iterator itse=seeds.find(endPt);
+ std::set<SPoint3>::iterator its=candidates.find(endPt);
+ if(itse != seeds.end()){
+ printf("found seed =%g %g %g \n", endPt.x(), endPt.y(), endPt.z());
+ seeds.erase(itse);
+ beginPt = *(seeds.begin());
+ std::set<SPoint3>::iterator itsee=candidates.find(beginPt);
+ if (itsee != candidates.end()) 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);
+ 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());
- // }
+
}
@@ -527,7 +544,10 @@ void MeshDelaunayVolume(std::vector<GRegion*> ®ions)
gr->set(faces);
//EMI VORONOI FOR CENRTERLINE
- if (Msg::GetVerbosity() == 20) voronoiDual(gr);
+ if (Msg::GetVerbosity() == 20) {
+ printVoronoi(gr);
+ skeletonFromVoronoi(gr);
+ }
// now do insertion of points
insertVerticesInRegion(gr);
diff --git a/Mesh/multiscalePartition.cpp b/Mesh/multiscalePartition.cpp
index 4d901129d1808cb78352f9f7db7b0988c8a264eb..96751985b04a1b3d7c44a026d25d568571a16c51 100644
--- a/Mesh/multiscalePartition.cpp
+++ b/Mesh/multiscalePartition.cpp
@@ -189,7 +189,7 @@ static int getAspectRatio(std::vector<MElement *> &elements,
//double H = obbox.getMaxSize();
double D = H;
- if (boundaries.size() > 0 ) D = 0.0;
+ if (boundaries.size() > 0 ) D = 10e4;
for (unsigned int i = 0; i < boundaries.size(); i++){
std::set<MVertex*> vb;
std::vector<MEdge> iBound = boundaries[i];
@@ -206,7 +206,7 @@ static int getAspectRatio(std::vector<MElement *> &elements,
bb +=pt;
}
double iD = norm(SVector3(bb.max(), bb.min()));
- D = std::max(D, iD);
+ D = std::min(D, iD);
//SOrientedBoundingBox obboxD = SOrientedBoundingBox::buildOBB(vBounds);
//D = std::max(D, obboxD.getMaxSize());
@@ -289,7 +289,8 @@ static void printLevel(std::vector<MElement *> &elements, int recur, int region)
}
multiscalePartition::multiscalePartition(std::vector<MElement *> &elements,
- int nbParts, typeOfPartition method)
+ int nbParts, typeOfPartition method,
+ int allowPartition)
{
options = CTX::instance()->partitionOptions;
options.num_partitions = nbParts;
@@ -305,6 +306,8 @@ multiscalePartition::multiscalePartition(std::vector<MElement *> &elements,
level->region = 0;
levels.push_back(level);
+ onlyMultilevel = false;
+ if (allowPartition == 2) onlyMultilevel = true;
partition(*level, nbParts, method);
@@ -361,12 +364,18 @@ void multiscalePartition::partition(partitionLevel & level, int nbParts,
nextLevel->recur,nextLevel->region, genus, AR, nbParts);
partition(*nextLevel, nbParts, MULTILEVEL);
}
- else if (genus == 0 && AR > 5 || genus == 0 && NB > 1){
+ else if (genus == 0 && AR > 5 ){//|| genus == 0 && NB > 1){
int nbParts = 2;
- Msg::Info("Mesh partition: level (%d-%d) is ZERO-GENUS (AR=%d NB=%d) ---> LAPLACIAN partition %d parts",
- nextLevel->recur,nextLevel->region, AR, NB, nbParts);
- //partition(*nextLevel, nbParts, MULTILEVEL);
- partition(*nextLevel, nbParts, LAPLACIAN);
+ if(!onlyMultilevel){
+ Msg::Info("Mesh partition: level (%d-%d) is ZERO-GENUS (AR=%d NB=%d) ---> LAPLACIAN partition %d parts",
+ nextLevel->recur,nextLevel->region, AR, NB, nbParts);
+ partition(*nextLevel, nbParts, LAPLACIAN);
+ }
+ else {
+ Msg::Info("Mesh partition: level (%d-%d) is ZERO-GENUS (AR=%d NB=%d) ---> MULTILEVEL partition %d parts",
+ nextLevel->recur,nextLevel->region, AR, NB, nbParts);
+ partition(*nextLevel, nbParts, MULTILEVEL);
+ }
}
else {
Msg::Info("*** Mesh partition: level (%d-%d) is ZERO-GENUS (AR=%d, NB=%d)",
diff --git a/Mesh/multiscalePartition.h b/Mesh/multiscalePartition.h
index 15c4cecfa06a1632f63a622935646cb84cbe0e7e..f1ff46b60461e7d1203a75142dc1911abf2254d5 100644
--- a/Mesh/multiscalePartition.h
+++ b/Mesh/multiscalePartition.h
@@ -29,10 +29,12 @@ class multiscalePartition{
std::vector<partitionLevel*> levels;
void partition(partitionLevel &level, int nbParts, typeOfPartition method);
int totalParts;
+ bool onlyMultilevel;
meshPartitionOptions options;
public:
- multiscalePartition(std::vector<MElement *> &elements, int nbParts, typeOfPartition method);
+ multiscalePartition(std::vector<MElement *> &elements, int nbParts,
+ typeOfPartition method, int allowPartition);
int assembleAllPartitions();
void setNumberOfPartitions(int &nbParts);
int getNumberOfParts(){return totalParts;}
diff --git a/Numeric/Numeric.cpp b/Numeric/Numeric.cpp
index 14fd1cc78045ea3ef20336487e77a38c42203bd1..25ae05b37053ddfaae5ce583f40d4d766d93b2c5 100644
--- a/Numeric/Numeric.cpp
+++ b/Numeric/Numeric.cpp
@@ -476,6 +476,20 @@ double ComputeScalarRep(int numComp, double *val)
return 0.;
}
+void eigenvalue2x2(double mat[2][2], double v[2])
+{
+
+ double a=1;
+ double b=-(mat[0][0]+mat[1][1]);
+ double c= (mat[0][0]*mat[1][1])-(mat[0][1]*mat[1][0]);
+
+ double det = b*b-4.*a*c;
+
+ v[0] = (-b+sqrt(det))/(2*a);
+ v[1] = (-b-sqrt(det))/(2*a);
+
+}
+
void eigenvalue(double mat[3][3], double v[3])
{
// characteristic polynomial of T : find v root of
diff --git a/Numeric/Numeric.h b/Numeric/Numeric.h
index 33643dd6184f094a12bcb3166e82ed6f74d7dc6b..7aa6d64fbb5228cf9ee30c9dd102e74f899ab3ad 100644
--- a/Numeric/Numeric.h
+++ b/Numeric/Numeric.h
@@ -72,6 +72,7 @@ void planarQuad_xyz2xy(double *x, double *y, double *z, double *xn, double *yn);
double computeInnerRadiusForQuad(double *x, double *y, int i);
char float2char(float f);
float char2float(char c);
+void eigenvalue2x2(double mat[2][2], double v[2]);
void eigenvalue(double mat[3][3], double re[3]);
void FindCubicRoots(const double coeff[4], double re[3], double im[3]);
void eigsort(double d[3]);
diff --git a/Parser/Gmsh.tab.cpp b/Parser/Gmsh.tab.cpp
index 7e76eb6bb6e006b91828e9fb16d3090517c6b924..2a6ccc7f83788000c321e6e356e248df38244b65 100644
--- a/Parser/Gmsh.tab.cpp
+++ b/Parser/Gmsh.tab.cpp
@@ -7145,6 +7145,10 @@ yyreduce:
(yyval.i) = 2;
else if(!strcmp((yyvsp[(1) - (1)].c), "Conformal_NoSplit"))
(yyval.i) = -2;
+ else if(!strcmp((yyvsp[(1) - (1)].c), "Harmonic_SplitMetis"))
+ (yyval.i) = 3;
+ else if(!strcmp((yyvsp[(1) - (1)].c), "Conformal_SplitMetis"))
+ (yyval.i) = -3;
Free((yyvsp[(1) - (1)].c));
;}
break;
diff --git a/Parser/Gmsh.y b/Parser/Gmsh.y
index 6a1e97ee8b0469b8890d18a4f5097253c2a37e79..b3b056ad2d507eaba90c7d709b3fdd654eb39bad 100644
--- a/Parser/Gmsh.y
+++ b/Parser/Gmsh.y
@@ -2899,6 +2899,10 @@ CompoundMap :
$$ = 2;
else if(!strcmp($1, "Conformal_NoSplit"))
$$ = -2;
+ else if(!strcmp($1, "Harmonic_SplitMetis"))
+ $$ = 3;
+ else if(!strcmp($1, "Conformal_SplitMetis"))
+ $$ = -3;
Free($1);
}
;
diff --git a/Solver/multiscaleLaplace.cpp b/Solver/multiscaleLaplace.cpp
index d23372a6b722d6cd90f23b69a3b866a30b17f37e..ff10d5525737aea6727f469c3cc24596c7a5dc4e 100644
--- a/Solver/multiscaleLaplace.cpp
+++ b/Solver/multiscaleLaplace.cpp
@@ -706,6 +706,7 @@ static void printLevel_onlysmall(const char* fn,
std::map<MVertex*,SPoint2> *coordinates,
double version,
double tolerance){
+
std::vector<MElement *> small;
double dx[3] = {0,0,0};
int COUNT = 0;
@@ -1101,12 +1102,12 @@ void multiscaleLaplace::parametrize (multiscaleLaplaceLevel & level){
}
//Save multiscale meshes
- std::string name1(level._name+"real.msh");
- std::string name2(level._name+"param.msh");
- std::string name3(level._name+"param_small.msh");
- printLevel (name1.c_str(),level.elements,0,2.2);
- printLevel (name2.c_str(),level.elements,&level.coordinates,2.2);
- printLevel_onlysmall (name3.c_str(),level.elements,&level.coordinates,2.2,1.e-15);
+ // std::string name1(level._name+"real.msh");
+ // std::string name2(level._name+"param.msh");
+ // std::string name3(level._name+"param_small.msh");
+ // printLevel (name1.c_str(),level.elements,0,2.2);
+ // printLevel (name2.c_str(),level.elements,&level.coordinates,2.2);
+ // printLevel_onlysmall (name3.c_str(),level.elements,&level.coordinates,2.2,1.e-15);
//For every small region compute a new parametrization
Msg::Info("Level (%d-%d): %d connected small regions",level.recur, level.region, regions.size());
@@ -1215,22 +1216,22 @@ void multiscaleLaplace::cut (std::vector<MElement *> &elements)
recur_cut_ (1.0, M_PI, 0.0, root,left,right);
connected_left_right(left, right);
- printLevel ("Rootcut-left.msh",left,0,2.2);
- printLevel ("Rootcut-right.msh",right,0,2.2);
- printLevel ("Rootcut-all.msh",elements, 0,2.2);
+ // printLevel ("Rootcut-left.msh",left,0,2.2);
+ // printLevel ("Rootcut-right.msh",right,0,2.2);
+ // printLevel ("Rootcut-all.msh",elements, 0,2.2);
- printLevel ("Rootcut-left-param.msh",left,&root->coordinates,2.2);
- printLevel_onlysmall ("Rootcut-left-param10.msh",left,&root->coordinates,2.2,1.e-10);
- printLevel_onlysmall ("Rootcut-left-param12.msh",left,&root->coordinates,2.2,1.e-12);
- printLevel_onlysmall ("Rootcut-left-param15.msh",left,&root->coordinates,2.2,1.e-15);
+ // printLevel ("Rootcut-left-param.msh",left,&root->coordinates,2.2);
+ // printLevel_onlysmall ("Rootcut-left-param10.msh",left,&root->coordinates,2.2,1.e-10);
+ // printLevel_onlysmall ("Rootcut-left-param12.msh",left,&root->coordinates,2.2,1.e-12);
+ // printLevel_onlysmall ("Rootcut-left-param15.msh",left,&root->coordinates,2.2,1.e-15);
- printLevel ("Rootcut-right-param.msh",right,&root->coordinates,2.2);
- printLevel_onlysmall ("Rootcut-right-param10.msh",right,&root->coordinates,2.2,1.e-10);
- printLevel_onlysmall ("Rootcut-right-param12.msh",right,&root->coordinates,2.2,1.e-12);
- printLevel_onlysmall ("Rootcut-right-param15.msh",right,&root->coordinates,2.2,1.e-15);
+ // printLevel ("Rootcut-right-param.msh",right,&root->coordinates,2.2);
+ // printLevel_onlysmall ("Rootcut-right-param10.msh",right,&root->coordinates,2.2,1.e-10);
+ // printLevel_onlysmall ("Rootcut-right-param12.msh",right,&root->coordinates,2.2,1.e-12);
+ // printLevel_onlysmall ("Rootcut-right-param15.msh",right,&root->coordinates,2.2,1.e-15);
- printLevel_onlysmall ("Rootcut-all-param12.msh",elements,&root->coordinates,2.2,1.e-12);
- printLevel_onlysmall ("Rootcut-all-param15.msh",elements,&root->coordinates,2.2,1.e-15);
+ // printLevel_onlysmall ("Rootcut-all-param12.msh",elements,&root->coordinates,2.2,1.e-12);
+ // printLevel_onlysmall ("Rootcut-all-param15.msh",elements,&root->coordinates,2.2,1.e-15);
if ( elements.size() != left.size()+right.size()) {
Msg::Error("Cutting laplace wrong nb elements (%d) != left + right (%d)", elements.size(), left.size()+right.size());
diff --git a/benchmarks/stl/squirrel.stl b/benchmarks/stl/squirrel.stl
new file mode 100644
index 0000000000000000000000000000000000000000..a238048c4a1648d9dcc25489d5c552bfb22a64a0
Binary files /dev/null and b/benchmarks/stl/squirrel.stl differ