diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index bd79868215947c061d95d7a94dc6385c4f747bda..fa6bb9272c65a055680e8e69b006b0810ab9d639 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -517,6 +517,7 @@ void GFaceCompound::fillNeumannBCS_Plane() const
if (_interior_loops.size()==1 && _type != SQUARE) return;
GModel::current()->setFactory("Gmsh");
+ //GModel::current()->setFactory("OpenCASCADE");
std::vector<std::vector<GFace *> > myFaceLoops;
std::vector<GFace *> myFaces;
for(std::list<std::list<GEdge*> >::const_iterator iloop = _interior_loops.begin();
@@ -2246,7 +2247,7 @@ bool GFaceCompound::checkTopology() const
double D = H;
if (_interior_loops.size() > 0) D = getSizeBB(_U0);
int AR1 = (int) checkAspectRatio();
- int AR2 = (int) (H/D + 0.5);
+ int AR2 = (int) floor(H/D+0.5);
int AR = std::max(AR1, AR2);
if (G != 0 || Nb < 1){
diff --git a/Geo/SVector3.h b/Geo/SVector3.h
index 0cd61b25c0de9f90c39aaea437108937a41cde7e..645058a3849be89edd3e532cfeda283c796c9830 100644
--- a/Geo/SVector3.h
+++ b/Geo/SVector3.h
@@ -142,6 +142,7 @@ inline void buildOrthoBasis_naive(SVector3 &dir, SVector3 &dir1, SVector3 &dir2)
dir2.normalize();
}
+//given a normal, build tangent and binormal
inline void buildOrthoBasis(SVector3 &normal, SVector3 &tangent, SVector3 &binormal)
{
//pick any Unit-Vector that's not parallel to normal:
@@ -164,4 +165,24 @@ inline void buildOrthoBasis(SVector3 &normal, SVector3 &tangent, SVector3 &binor
}
+//given a normal and guess tangent, build binormal
+inline void buildOrthoBasis2(SVector3 &normal, SVector3 &tangent, SVector3 &binormal)
+{
+
+ normal.normalize();
+ tangent.normalize();
+
+ //build a binormal from tangent and normal:
+ binormal = crossprod(tangent, normal);
+ double t1 = binormal.normalize();
+
+ //and correct the tangent from the binormal and the normal.
+ tangent = crossprod(normal, binormal);
+ double t2 = tangent.normalize();
+
+ if (t1 == 0.0 || t2 == 0.0)
+ buildOrthoBasis_naive(normal, tangent, binormal);
+
+}
+
#endif
diff --git a/Mesh/CenterlineField.cpp b/Mesh/CenterlineField.cpp
index 9f2cb2e7943a7d011e7c9d72b6271e23142b534a..ab7111e6faadb90a07fac8414ce9da2cb80e2321 100644
--- a/Mesh/CenterlineField.cpp
+++ b/Mesh/CenterlineField.cpp
@@ -47,7 +47,7 @@ void erase(std::vector<MLine*>& lines, MLine* l) {
double computeLength(std::vector<MLine*> lines){
double length= 0.0;
- for (int i = 0; i< lines.size(); i++){
+ for (unsigned int i = 0; i< lines.size(); i++){
length += lines[i]->getLength();
}
return length;
@@ -358,13 +358,14 @@ Centerline::~Centerline(){
if(nodesR) annDeallocPts(nodesR);
delete[]index;
delete[]dist;
+
}
void Centerline::importFile(std::string fileName){
current = GModel::current();
std::vector<GFace*> currentFaces = current->bindingsGetFaces();
- for (int i = 0; i < currentFaces.size(); i++){
+ for (unsigned int i = 0; i < currentFaces.size(); i++){
GFace *gf = currentFaces[i];
if (gf->geomType() == GEntity::DiscreteSurface){
for(unsigned int j = 0; j < gf->triangles.size(); j++)
@@ -389,7 +390,7 @@ void Centerline::importFile(std::string fileName){
int maxN = 0.0;
std::vector<GEdge*> modEdges = mod->bindingsGetEdges();
- for (int i = 0; i < modEdges.size(); i++){
+ for (unsigned int i = 0; i < modEdges.size(); i++){
GEdge *ge = modEdges[i];
for(unsigned int j = 0; j < ge->lines.size(); j++){
MLine *l = ge->lines[j];
@@ -514,7 +515,7 @@ void Centerline::distanceToSurface(){
//COMPUTE WITH REVERSE ANN TREE (SURFACE POINTS IN TREE)
std::set<MVertex*> allVS;
- for(int j = 0; j < triangles.size(); j++)
+ for(unsigned int j = 0; j < triangles.size(); j++)
for(int k = 0; k<3; k++) allVS.insert(triangles[j]->getVertex(k));
int nbSNodes = allVS.size();
nodesR = annAllocPts(nbSNodes, 3);
@@ -592,12 +593,12 @@ void Centerline::buildKdTree(){
kdtree = new ANNkd_tree(nodes, nbNodes, 3);
- for(unsigned int i = 0; i < nbNodes; ++i){
+ for(int i = 0; i < nbNodes; ++i){
fprintf(f, "SP(%g,%g,%g){%g};\n",
nodes[i][0], nodes[i][1],nodes[i][2],1.0);
}
fprintf(f,"};\n");
- fclose(f);
+ fclose(f);
}
@@ -611,7 +612,6 @@ void Centerline::createSplitCompounds(){
// Remesh new faces (Compound Lines and Compound Surfaces)
Msg::Info("*** Starting parametrize compounds:");
- double t0 = Cpu();
//Parametrize Compound Lines
for (int i=0; i < NE; i++){
@@ -673,7 +673,7 @@ void Centerline::cleanMesh(){
mySplitMesh->addPhysicalEntity(2*NF+1);
current->setPhysicalName("surface mesh", 2, 2*NF+1);
current->add(mySplitMesh);
- for (int i=0; i < discFaces.size(); i++){
+ for (unsigned int i=0; i < discFaces.size(); i++){
GFace *gfc = current->getFaceByTag(NF+i+1);
for(unsigned int j = 0; j < gfc->triangles.size(); ++j){
MTriangle *t = gfc->triangles[j];
@@ -710,7 +710,7 @@ void Centerline::cleanMesh(){
GFace *gf = current->getFaceByTag(i+1);
current->remove(gf);
}
- for (int i=0; i < discFaces.size(); i++){
+ for (unsigned int i=0; i < discFaces.size(); i++){
GFace *gfc = current->getFaceByTag(NF+i+1);
current->remove(gfc);
}
@@ -731,7 +731,6 @@ void Centerline::createFaces(){
for(unsigned int i = 0; i < triangles.size(); ++i)
for(int j = 0; j < 3; j++)
e2e.insert(std::make_pair(triangles[i]->getEdge(j), triangles[i]));
- int iGroup = 0;
while(!e2e.empty()){
std::set<MTriangle*> group;
std::set<MEdge, Less_Edge> touched;
@@ -754,7 +753,6 @@ void Centerline::createFaces(){
Msg::Info("Centerline action (cutMesh) has cut surface mesh in %d faces ", (int)faces.size());
//create discFaces
- int numBef = current->getMaxElementaryNumber(2) + 1;
for(unsigned int i = 0; i < faces.size(); ++i){
int numF = current->getMaxElementaryNumber(2) + 1;
discreteFace *f = new discreteFace(current, numF);
@@ -762,7 +760,7 @@ void Centerline::createFaces(){
discFaces.push_back(f);
std::set<MVertex*> myVertices;
std::vector<MTriangle*> myFace = faces[i];
- for(int j= 0; j< myFace.size(); j++){
+ for(unsigned int j= 0; j< myFace.size(); j++){
MTriangle *t = myFace[j];
f->triangles.push_back(t);
for (int k= 0; k< 3; k++){
@@ -794,7 +792,7 @@ void Centerline::createClosedVolume(){
current->setFactory("Gmsh");
std::vector<std::vector<GFace *> > myFaceLoops;
std::vector<GFace *> myFaces;
- for (int i = 0; i< boundEdges.size(); i++){
+ for (unsigned int i = 0; i< boundEdges.size(); i++){
std::vector<std::vector<GEdge *> > myEdgeLoops;
std::vector<GEdge *> myEdges;
GEdge * gec = current->getEdgeByTag(NE+boundEdges[i]->tag());
@@ -867,7 +865,7 @@ void Centerline::cutMesh(){
int nbSplit = (int)floor(AR / 3. + 0.5);
double li = L/nbSplit;
double lc = 0.0;
- for (int j= 0; j < lines.size(); j++){
+ for (unsigned int j= 0; j < lines.size(); j++){
lc += lines[j]->getLength();
if (lc > li && nbSplit > 1) {
MVertex *v1 = lines[j]->getVertex(0);
@@ -1043,34 +1041,59 @@ void Centerline::operator() (double x, double y, double z, SMetric3 &metr, GEnt
double xyz[3] = {x,y,z};
//take xyz = closest point on boundary in case we are on the planar IN/OUT FACES or in VOLUME
+ bool onTubularSurface = true;
+ double ds = 0.0;
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);
+ onTubularSurface = false;
+ kdtreeR->annkSearch(xyz, 1, index, dist);
+ ds = sqrt(dist[0]);
xyz[0] = nodesR[index[0]][0];
xyz[1] = nodesR[index[0]][1];
xyz[2] = nodesR[index[0]][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;
+ ANNidxArray index2 = new ANNidx[2];
+ ANNdistArray dist2 = new ANNdist[2];
+ kdtree->annkSearch(xyz, 2, index2, dist2);
+ double radMax = sqrt(dist2[0]);
SVector3 p0(nodes[index2[0]][0], nodes[index2[0]][1], nodes[index2[0]][2]);
SVector3 p1(nodes[index2[1]][0], nodes[index2[1]][1], nodes[index2[1]][2]);
- SVector3 dir_t = p1-p0;
- SVector3 dir_n1, dir_n2;
- buildOrthoBasis(dir_t,dir_n1,dir_n2);
-
- double lc_t = 3.*lc;
+
+ //dir_a = direction along the centerline
+ //dir_n = normal direction of the disk
+ //dir_t = tangential direction of the disk
+ SVector3 dir_a = p1-p0;
+ SVector3 dir_n(x-p0.x(), y-p0.y(), z-p0.z());
+ SVector3 dir_t;
+ buildOrthoBasis2(dir_a, dir_n, dir_t);
+
+ double lc = 2*M_PI*radMax/nbPoints;
+ double lc_a = 3.*lc;
+ double lc_n, lc_t;
+ if (onTubularSurface){
+ lc_n = lc_t = lc;
+ }
+ else{
+ double e = radMax/3.;
+ double hn = e/10.;
+ double rm = std::max(radMax-ds, radMax-e);
+ lc_t = 2*M_PI*rm/nbPoints;
+ //lc_n = lc_t = lc;
+ //double ratio = 1.02; //1. + (lc_t-hn)/e;
+ //printf("ratio =%g \n", ratio);
+ //lc_n = ds*(ratio-1) + hn;
+ if (ds < e) lc_n = hn;
+ else lc_n = lc_t;
+ }
+ double lam_a = 1./(lc_a*lc_a);
+ double lam_n = 1./(lc_n*lc_n);
double lam_t = 1./(lc_t*lc_t);
- double lam_n = 1./(lc*lc);
- metr = SMetric3(lam_t,lam_n,lam_n, dir_t, dir_n1, dir_n2);
+
+ metr = SMetric3(lam_a,lam_n,lam_t, dir_a, dir_n, dir_t);
delete[]index2;
delete[]dist2;
diff --git a/Mesh/CenterlineField.h b/Mesh/CenterlineField.h
index f43e6241fd2fc1d57dddf2eb27d821f16dca5609..efba5de307a36dd0d07529ebcb816adf121bbe7a 100644
--- a/Mesh/CenterlineField.h
+++ b/Mesh/CenterlineField.h
@@ -66,6 +66,7 @@ class Centerline : public Field{
bool is_cut;
bool is_closed;
+
//all (unique) lines of centerlines
std::vector<MLine*> lines;
//the stuctured tree of the centerlines