diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 8aaef75f2fb9e58a1f81ca068d3d27621b1572d3..201134c40424f416a1ac3d7f0d06e50b2c0034fb 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -15,6 +15,7 @@
#if defined(HAVE_SOLVER)
+#include "Options.h"
#include "MLine.h"
#include "MTriangle.h"
#include "Numeric.h"
@@ -43,6 +44,7 @@
#include "Curvature.h"
#include "MPoint.h"
#include "Numeric.h"
+#include "meshGFace.h"
static void fixEdgeToValue(GEdge *ed, double value, dofManager<double> &myAssembler)
{
@@ -422,6 +424,130 @@ static MTriangle* findInvertedTri(v2t_cont &adjv, MVertex*v0, MVertex*v1, MVerte
return myTri;
+}
+
+void GFaceCompound::orientFillTris(std::list<MTriangle*> loopfillTris) const{
+
+ //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());
+
+}
+
+void GFaceCompound::printFillTris() const{
+
+ 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", 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);
+ }
+ }
+}
+
+void GFaceCompound::fillNeumannBCS_Plane() const
+{
+
+ Msg::Info("Meshing %d interior holes with planes ", _interior_loops.size()-1);
+
+ fillTris.clear();
+ fillNodes.clear();
+
+ if (_interior_loops.size()==1 && _type != SQUARE) return;
+
+ GModel::current()->setFactory("Gmsh");
+ std::vector<std::vector<GFace *> > myFaceLoops;
+ std::vector<GFace *> myFaces;
+ 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<std::vector<GEdge *> > myEdgeLoops;
+ std::vector<GEdge*> myEdges;
+ for (std::list<GEdge*>::iterator itl = loop.begin(); itl != loop.end(); itl++)
+ myEdges.push_back(*itl);
+ myEdgeLoops.push_back(myEdges);
+ GFace *newFace = GModel::current()->addPlanarFace(myEdgeLoops);
+ fillFaces.push_back(newFace);
+ int meshingAlgo = CTX::instance()->mesh.algo2d;
+ opt_mesh_algo2d(0, GMSH_SET, 1.0); //mesh adapt
+ meshGFace mgf;
+ mgf(newFace);
+ opt_mesh_algo2d(0, GMSH_SET, meshingAlgo);
+ for(unsigned int i = 0; i < newFace->triangles.size(); ++i){
+ loopfillTris.push_back(newFace->triangles[i]);
+ fillNodes.insert(newFace->triangles[i]->getVertex(0));
+ fillNodes.insert(newFace->triangles[i]->getVertex(1));
+ fillNodes.insert(newFace->triangles[i]->getVertex(2));
+ }
+ //mod->remove(newFace);
+ }
+ orientFillTris(loopfillTris);
+ }
+
+ printFillTris();
+
}
void GFaceCompound::fillNeumannBCS() const
{
@@ -483,83 +609,13 @@ void GFaceCompound::fillNeumannBCS() const
loopfillTris.push_back(new MTriangle(v4, c, v5));
loopfillTris.push_back(new MTriangle(v0, v3, v1));
}
-
- }
-
- //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;
- }
}
+ orientFillTris(loopfillTris);
}
- 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", 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);
- }
- }
-}
+ printFillTris();
+}
bool GFaceCompound::trivial() const
{
return false;
@@ -650,9 +706,8 @@ bool GFaceCompound::checkOrientation(int iter, bool moveBoundaries) const
if(!oriented && iter < iterMax){
if (moveBoundaries){
if (iter ==0) Msg::Info("--- Flipping : moving boundaries.");
- Msg::Info("--- Moving Boundary - iter %d -",iter);
+ Msg::Debug("--- Moving Boundary - iter %d -",iter);
convexBoundary(nTot/fabs(nTot));
- one2OneMap();
printStuff(iter);
return checkOrientation(iter+1, moveBoundaries);
}
@@ -660,9 +715,8 @@ bool GFaceCompound::checkOrientation(int iter, bool moveBoundaries) const
if (iter ==0) Msg::Info("--- Flipping : applying cavity checks.");
Msg::Debug("--- Cavity Check - iter %d -",iter);
bool success = one2OneMap();
- if (success) return checkOrientation(iter+1, moveBoundaries);
+ if (success) return checkOrientation(iter+1);
}
-
}
else if (iter > 0 && iter < iterMax){
Msg::Info("--- Flipping : no more flips (%d iter)", iter);
@@ -890,20 +944,20 @@ bool GFaceCompound::parametrize() const
if(allNodes.empty()) buildAllNodes();
if (_type != SQUARE){
- printf("ordering vertices \n");
bool success = orderVertices(_U0, _ordered, _coords);
if(!success) {Msg::Error("Could not order vertices on boundary");exit(1);}
}
+ fillNeumannBCS_Plane();
+ //fillNeumannBCS();
+
// Convex parametrization
if (_mapping == CONVEX){
- Msg::Info("Parametrizing surface %d with 'convex map'", tag());
- fillNeumannBCS();
+ Msg::Info("Parametrizing surface %d with 'convex map'", tag());
parametrize(ITERU,CONVEX);
parametrize(ITERV,CONVEX);
if (_type==MEANPLANE){
checkOrientation(0, true);
- // printStuff(22);
// _type = ALREADYFIXED;
// parametrize(ITERU,CONVEX);
// parametrize(ITERV,CONVEX);
@@ -913,15 +967,12 @@ bool GFaceCompound::parametrize() const
// Laplace parametrization
else if (_mapping == HARMONIC){
Msg::Info("Parametrizing surface %d with 'harmonic map'", tag());
- fillNeumannBCS();
parametrize(ITERU,HARMONIC);
parametrize(ITERV,HARMONIC);
if (_type == MEANPLANE) checkOrientation(0, true);
}
// Conformal map parametrization
else if (_mapping == CONFORMAL){
-
- fillNeumannBCS();
std::vector<MVertex *> vert;
bool oriented, hasOverlap;
if (_type == SPECTRAL){
@@ -933,7 +984,7 @@ bool GFaceCompound::parametrize() const
parametrize_conformal(0, NULL, NULL);
}
printStuff(55);
- oriented = checkOrientation(0);
+ oriented = checkOrientation(0);
printStuff(66);
if (!oriented) oriented = checkOrientation(0, true);
printStuff(77);
@@ -970,7 +1021,6 @@ bool GFaceCompound::parametrize() const
_type = UNITCIRCLE;
coordinates.clear();
Octree_Delete(oct);
- fillNeumannBCS();
parametrize(ITERU,CONVEX);
parametrize(ITERV,CONVEX);
checkOrientation(0);
@@ -978,6 +1028,10 @@ bool GFaceCompound::parametrize() const
}
}
+ for (int i= 0; i< fillFaces.size(); i++){
+ GModel::current()->remove(fillFaces[i]);
+ }
+
return paramOK;
}
@@ -1250,6 +1304,7 @@ GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
GFaceCompound::~GFaceCompound()
{
+
if(oct){
Octree_Delete(oct);
delete [] _gfct;
diff --git a/Geo/GFaceCompound.h b/Geo/GFaceCompound.h
index 1a7444e635599d184ef3df3faca44d01bebd72a3..d9fffb7fa7cd0e71242fc91411c4b566ea480217 100644
--- a/Geo/GFaceCompound.h
+++ b/Geo/GFaceCompound.h
@@ -57,7 +57,8 @@ class GFaceCompound : public GFace {
CONVEX_CIRCLE=4,CONVEX_PLANE=5, HARMONIC_SQUARE=6, CONFORMAL_FE=7} typeOfCompound;
typedef enum {HARMONIC=0,CONFORMAL=1, RBF=2, CONVEX=3} typeOfMapping;
typedef enum {UNITCIRCLE, MEANPLANE, SQUARE, ALREADYFIXED,SPECTRAL, FE} typeOfIsomorphism;
- void computeNormals(std::map<MVertex*, SVector3> &normals) const;
+ mutable int nbSplit;
+
protected:
mutable std::set<MVertex *> ov;
mutable GRbf *_rbf;
@@ -77,21 +78,32 @@ class GFaceCompound : public GFace {
mutable std::map<MVertex*, SVector3> _normals;
mutable std::list<MTriangle*> fillTris;
mutable std::set<MVertex*> fillNodes;
+ mutable std::vector<GFace*> fillFaces;
mutable std::vector<MVertex*> _ordered;
mutable std::vector<double> _coords;
mutable std::map<MVertex*, int> _mapV;
linearSystem <double> *_lsys;
void buildOct() const ;
void buildAllNodes() const;
+
+ //different type of parametrizations
void parametrize(iterationStep, typeOfMapping) const;
bool parametrize_conformal(int iter, MVertex *v1, MVertex *v2) const;
bool parametrize_conformal_spectral() const;
- void compute_distance() const;
+
+ //check for parametrizations
bool checkOrientation(int iter, bool moveBoundaries=false) const;
bool checkOverlap(std::vector<MVertex *> &vert) const;
bool one2OneMap() const;
void convexBoundary(double nTot) const;
double checkAspectRatio() const;
+
+ //tools for filling interior holes of surfaces
+ void fillNeumannBCS() const;
+ void fillNeumannBCS_Plane() const;
+ void orientFillTris(std::list<MTriangle*> loopfillTris)const;
+ void printFillTris()const;
+
void computeNormals () const;
void getBoundingEdges();
void getUniqueEdges(std::set<GEdge*> &_unique);
@@ -99,13 +111,12 @@ class GFaceCompound : public GFace {
void getTriangle(double u, double v, GFaceCompoundTriangle **lt,
double &_u, double &_v) const;
virtual double locCurvature(MTriangle *t, double u, double v) const;
- void printStuff(int iNewton=0) const;
- bool trivial() const;
+
double getSizeH() const;
double getSizeBB(const std::list<GEdge* > &elist) const;
- void fillNeumannBCS() const;
- /* double sumAngles(std::vector<MVertex*> ordered) const; */
-
+ bool trivial() const;
+ void printStuff(int iNewton=0) const;
+
public:
GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
std::list<GEdge*> &U0, typeOfCompound typ = HARMONIC_CIRCLE,
@@ -117,11 +128,12 @@ class GFaceCompound : public GFace {
typeOfCompound typ = HARMONIC_CIRCLE,
int allowPartition=1,
linearSystem<double>* lsys =0);
- virtual ~GFaceCompound();
+ ~GFaceCompound();
+
Range<double> parBounds(int i) const
{ return trivial() ? (*(_compound.begin()))->parBounds(i) : Range<double>(-1, 1); }
- virtual GPoint point(double par1, double par2) const;
- typeOfCompound getTypeOfCompound() { return _toc;}
+
+ virtual GPoint point(double par1, double par2) const;
SPoint2 parFromPoint(const SPoint3 &p, bool onSurface=true) const;
virtual Pair<SVector3,SVector3> firstDer(const SPoint2 ¶m) const;
virtual void secondDer(const SPoint2 &, SVector3 *, SVector3 *, SVector3 *) const;
@@ -131,19 +143,24 @@ class GFaceCompound : public GFace {
virtual SPoint2 getCoordinates(MVertex *v) const;
virtual double curvatureMax(const SPoint2 ¶m) const;
virtual double curvatures(const SPoint2 ¶m, SVector3 *dirMax, SVector3 *dirMin,
- double *curvMax, double *curvMin) const;
- virtual int genusGeom () const;
- virtual bool checkTopology() const;
- bool parametrize() const ;
+ double *curvMax, double *curvMin) const;
+ bool parametrize() const;
+ void computeNormals(std::map<MVertex*, SVector3> &normals) const;
void coherenceNormals();
void coherencePatches() const;
+ virtual int genusGeom () const;
+ virtual bool checkTopology() const;
+
virtual std::list<GFace*> getCompounds() const { return _compound; }
- mutable int nbSplit;
+ typeOfCompound getTypeOfCompound() { return _toc;}
int getNbSplit() const { return nbSplit; }
int allowPartition() const{ return _allowPartition; }
void setType(typeOfIsomorphism type){ _type=type;}
- // useful for mesh generators ----------------------------------------
- GPoint intersectionWithCircle (const SVector3 &n1, const SVector3 &n2, const SVector3 &p, const double &d, double uv[2]) const;
+
+ // useful for mesh generators
+ GPoint intersectionWithCircle (const SVector3 &n1, const SVector3 &n2, const SVector3 &p,
+ const double &d, double uv[2]) const;
+
private:
mutable typeOfCompound _toc;
mutable typeOfMapping _mapping;
diff --git a/Geo/GModelFactory.cpp b/Geo/GModelFactory.cpp
index d0037053ed7d8b5e3296fbab8663c30ab4688840..d173108491d3a94c427e4fea3c27f4138d84c65b 100644
--- a/Geo/GModelFactory.cpp
+++ b/Geo/GModelFactory.cpp
@@ -64,7 +64,7 @@ GEdge *GeoFactory::addLine(GModel *gm, GVertex *start, GVertex *end)
GFace *GeoFactory::addPlanarFace(GModel *gm, std::vector< std::vector<GEdge *> > edges)
{
- //create line loops
+ //create line loops
int nLoops = edges.size();
std::vector<EdgeLoop *> vecLoops;
for (int i=0; i< nLoops; i++){
@@ -107,9 +107,7 @@ GFace *GeoFactory::addPlanarFace(GModel *gm, std::vector< std::vector<GEdge *> >
for(int i = 0; i < gec.size(); i++)
c->compound.push_back(gec[i]->tag());
}
- else {
- printf("type not implemented \n"); exit(1);
- }
+
c->Control_Points = List_Create(2, 1, sizeof(Vertex *));
List_Add(c->Control_Points, &vertb);
List_Add(c->Control_Points, &verte);
@@ -121,7 +119,7 @@ GFace *GeoFactory::addPlanarFace(GModel *gm, std::vector< std::vector<GEdge *> >
CreateReversedCurve(c);
List_Delete(temp);
}
- List_Add(temp, &numEdge);
+ List_Add(temp, &numEdge);
}
int num = gm->getMaxElementaryNumber(2) + 1+i;
@@ -144,6 +142,7 @@ GFace *GeoFactory::addPlanarFace(GModel *gm, std::vector< std::vector<GEdge *> >
int numl = vecLoops[i]->Num;
List_Add(temp, &numl);
}
+
setSurfaceGeneratrices(s, temp);
List_Delete(temp);
End_Surface(s);
diff --git a/Geo/Geo.cpp b/Geo/Geo.cpp
index 3f91a450d35a96accbe913599b8275287e9fb7ea..b7d4b7938c9b206889ad18513d126a669ce41da8 100644
--- a/Geo/Geo.cpp
+++ b/Geo/Geo.cpp
@@ -3484,8 +3484,8 @@ void sortEdgesInLoop(int num, List_T *edges)
if((c = FindCurve(j))){
List_Add(temp, &c);
if(c->Typ == MSH_SEGM_DISCRETE){
- Msg::Warning("Aborting line loop sort for discrete edge: hope you know "
- "what you're doing ;-)");
+ Msg::Debug("Aborting line loop sort for discrete edge: hope you know "
+ "what you're doing ;-)");
return;
}
}
diff --git a/Mesh/CenterlineField.cpp b/Mesh/CenterlineField.cpp
index 040bafa1b3b0ac91c79d35c1ff50181e4c922746..b143c86064f1ac8135f0960d063daf504e51f61f 100644
--- a/Mesh/CenterlineField.cpp
+++ b/Mesh/CenterlineField.cpp
@@ -634,9 +634,9 @@ void Centerline::createSplitCompounds(){
int num_gfc = NF+i+1;
Msg::Info("Parametrize Compound Surface (%d) = %d discrete face",
num_gfc, pf->tag());
- GFaceCompound::typeOfCompound typ = GFaceCompound::CONVEX_CIRCLE;
+ //GFaceCompound::typeOfCompound typ = GFaceCompound::CONVEX_CIRCLE;
//GFaceCompound::typeOfCompound typ = GFaceCompound::HARMONIC_PLANE;
- //GFaceCompound::typeOfMapping typ = GFaceCompound::CONFORMAL_SPECTRAL;
+ GFaceCompound::typeOfCompound typ = GFaceCompound::CONFORMAL_SPECTRAL;
GFaceCompound *gfc = new GFaceCompound(current, num_gfc, f_compound, U0,
typ, 0);
gfc->meshAttributes.recombine = recombine;
@@ -1008,11 +1008,10 @@ double Centerline::operator() (double x, double y, double z, GEntity *ge){
double xyz[3] = {x,y,z};
+ //take xyz = closest point on boundary in case we are on the planar in/out faces
bool isCompound = (ge->dim() == 2 && ge->geomType() == GEntity::CompoundSurface) ? true: false;
std::list<GFace*> cFaces;
if (isCompound) cFaces = ((GFaceCompound*)ge)->getCompounds();
-
- //take xyz = closest point on boundary in case we are on the planar in/out faces
if ( ge->dim() == 3 || (ge->dim() == 2 && ge->geomType() == GEntity::Plane) ||
(isCompound && (*cFaces.begin())->geomType() == GEntity::Plane) ){
int num_neighbours = 1;
@@ -1042,13 +1041,24 @@ void Centerline::operator() (double x, double y, double z, SMetric3 &metr, GEnt
update_needed = false;
}
- //double lc = operator()(x,y,z,ge);
- //metr = SMetric3(1./(lc*lc));
+ double xyz[3] = {x,y,z};
+
+ //take xyz = closest point on boundary in case we are on the planar in/out faces
+ bool isCompound = (ge->dim() == 2 && ge->geomType() == GEntity::CompoundSurface) ? true: false;
+ std::list<GFace*> cFaces;
+ if (isCompound) cFaces = ((GFaceCompound*)ge)->getCompounds();
+ if ( ge->dim() == 3 || (ge->dim() == 2 && ge->geomType() == GEntity::Plane) ||
+ (isCompound && (*cFaces.begin())->geomType() == GEntity::Plane) ){
+ int num_neighbours = 1;
+ kdtreeR->annkSearch(xyz, num_neighbours, index, dist);
+ xyz[0] = nodesR[index[0]][0];
+ xyz[1] = nodesR[index[0]][1];
+ xyz[2] = nodesR[index[0]][2];
+ }
- double xyz[3] = {x,y,z };
- ANNidxArray index2 = new ANNidx[2];
- ANNdistArray dist2 = new ANNdist[2];
int num_neighbours = 2;
+ ANNidxArray index2 = new ANNidx[num_neighbours];
+ ANNdistArray dist2 = new ANNdist[num_neighbours];
kdtree->annkSearch(xyz, num_neighbours, index2, dist2);
double d = sqrt(dist2[0]);
double lc = 2*M_PI*d/nbPoints;
@@ -1058,7 +1068,7 @@ void Centerline::operator() (double x, double y, double z, SMetric3 &metr, GEnt
SVector3 dir1, dir2;
buildOrthoBasis(dir0,dir1,dir2);
- double lcA = 4.*lc;
+ double lcA = 3.*lc;
double lam1 = 1./(lcA*lcA);
double lam2 = 1./(lc*lc);
metr = SMetric3(lam1,lam2,lam2, dir0, dir1, dir2);
diff --git a/Mesh/meshMetric.cpp b/Mesh/meshMetric.cpp
index e4720b58bec280e43962085e5e94362606ba79e2..6d31391c1a19d590afbff1ad6dd3e8f9df08e66f 100644
--- a/Mesh/meshMetric.cpp
+++ b/Mesh/meshMetric.cpp
@@ -372,8 +372,8 @@ void meshMetric::computeMetric(){
computeValues();
computeHessian_FE();
- //computeHessian_LS();
-
+ //computeHessian_LS();
+
int metricNumber = setOfMetrics.size();
v2t_cont :: iterator it = _adj.begin();