diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 72b4d8865f88f418f832746ca29ebe555061f3ee..46c86ddce26fd03e480af12b592fe1f1f177125d 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -1929,6 +1929,7 @@ GPoint GFaceCompound::pointInRemeshedOctree(double par1, double par2) const
//create new octree with new mesh elements
if(!octNew){
+ printf("create new octrre \n");
std::vector<MElement *> myElems;
for (unsigned int i = 0; i < triangles.size(); i++) myElems.push_back(triangles[i]);
for (unsigned int i = 0; i < quadrangles.size(); i++) myElems.push_back(quadrangles[i]);
@@ -1960,10 +1961,13 @@ GPoint GFaceCompound::pointInRemeshedOctree(double par1, double par2) const
}
}
+ printf("size octrre new = %d \n", myParamElems.size());
octNew = new MElementOctree(myParamElems);
//build kdtree boundary nodes in parametric space
- uv_nodes = annAllocPts(myBCNodes.size(), 3);
+ printf("build bc kdtree \n");
+ int nbBCNodes = myBCNodes.size();
+ uv_nodes = annAllocPts(nbBCNodes, 3);
std::set<SPoint2>::iterator itp = myBCNodes.begin();
int ind = 0;
while (itp != myBCNodes.end()){
@@ -1973,16 +1977,13 @@ GPoint GFaceCompound::pointInRemeshedOctree(double par1, double par2) const
uv_nodes[ind][2] = 0.0;
itp++; ind++;
}
- uv_kdtree = new ANNkd_tree(uv_nodes, myBCNodes.size(), 3);
+ uv_kdtree = new ANNkd_tree(uv_nodes, nbBCNodes, 3);
}
//now use new octree to find point
double uvw[3]={par1,par2, 0.0};
double UV[3];
- double initialTol = MElement::getTolerance();
- MElement::setTolerance(1.e-2);
MElement *e = octNew->find(par1,par2, 0.0,-1,false);
- MElement::setTolerance(initialTol);
if (e){
e->xyz2uvw(uvw,UV);
double valX[8], valY[8], valZ[8];
@@ -2000,6 +2001,7 @@ GPoint GFaceCompound::pointInRemeshedOctree(double par1, double par2) const
}
//if element not found in new octree find closest point
else{
+ printf("not found in new octree \n");
GPoint gp(50,50,50);
double pt[3] = {par1,par2,0.0};
uv_kdtree->annkSearch(pt, 2, index, dist);
@@ -2053,7 +2055,14 @@ GPoint GFaceCompound::point(double par1, double par2) const
getTriangle(par1, par2, <, U,V);
if(!lt){
//printf("POINT (%g %g) NOT FOUND --> find closest \n", par1,par2);
- if (meshStatistics.status != GFace::DONE ) {
+ if (meshStatistics.status == GFace::DONE && triangles.size()+quadrangles.size() != 0) {
+ //printf("look in remeshed octree \n");
+ GPoint gp = pointInRemeshedOctree(par1,par2);
+ gp.setNoSuccess();
+ return gp;
+ }
+ else{
+ //printf("look in kdtree \n");
double pt[3] = {par1,par2, 0.0};
kdtree->annkSearch(pt, 1, index, dist);
lt = &(_gfct[index[0]]);
@@ -2081,12 +2090,6 @@ GPoint GFaceCompound::point(double par1, double par2) const
V = X[1];
//printf("found closest point (%g %g) U V =%g %g \n", u,v, U,V);
}
- else{
- //printf("look in remeshed octree \n");
- GPoint gp = pointInRemeshedOctree(par1,par2);
- gp.setNoSuccess();
- return gp;
- }
}
if (lt->gf->geomType() != GEntity::DiscreteSurface){
@@ -2164,8 +2167,8 @@ Pair<SVector3,SVector3> GFaceCompound::firstDer(const SPoint2 ¶m) const
MTriangle *tri;
if (lt) tri = lt->tri;
else {
- //printf("FIRSTDER POINT NOT FOUND --> kdtree \n");
- //printf("uv=%g %g \n", param.x(), param.y());
+ printf("FIRSTDER POINT NOT FOUND --> kdtree \n");
+ printf("uv=%g %g \n", param.x(), param.y());
double pt[3] = {param.x(), param.y(), 0.0};
kdtree->annkSearch(pt, 1, index, dist);
tri = (&_gfct[index[0]])->tri;
diff --git a/Mesh/CenterlineField.cpp b/Mesh/CenterlineField.cpp
index c4bfd84f5b6e6859ebcdf380d78fb4d7c757a3c2..e0d6e32684fbf185237ef35532948d6ba7cde35e 100644
--- a/Mesh/CenterlineField.cpp
+++ b/Mesh/CenterlineField.cpp
@@ -638,10 +638,12 @@ void Centerline::createSplitCompounds()
Msg::Info("Create Compound Line (%d) = %d discrete edge",
num_gec, pe->tag());
GEdge *gec = current->addCompoundEdge(e_compound,num_gec);
- gec->meshAttributes.method = MESH_TRANSFINITE;
- gec->meshAttributes.nbPointsTransfinite = nbPoints+1;
- gec->meshAttributes.typeTransfinite = 0;
- gec->meshAttributes.coeffTransfinite = 1.0;
+ if (CTX::instance()->mesh.algo2d != ALGO_2D_BAMG){
+ gec->meshAttributes.method = MESH_TRANSFINITE;
+ gec->meshAttributes.nbPointsTransfinite = nbPoints+1;
+ gec->meshAttributes.typeTransfinite = 0;
+ gec->meshAttributes.coeffTransfinite = 1.0;
+ }
}
// Parametrize Compound surfaces
diff --git a/Mesh/directions3D.cpp b/Mesh/directions3D.cpp
index 07456ece995d186965fbd767f7cfcd44bb856c11..f2041fedcc992f01bf8299e33672f2e59c3a7d8e 100644
--- a/Mesh/directions3D.cpp
+++ b/Mesh/directions3D.cpp
@@ -872,73 +872,82 @@ void Frame_field::save(const std::vector<std::pair<SPoint3, STensor3> > data,
file.close();
}
-void Frame_field::recur_connect_vert(MVertex *v,
+void Frame_field::recur_connect_vert(FILE *fi, int count,
+ MVertex *v,
STensor3 &cross,
std::multimap<MVertex*,MVertex*> &v2v,
std::set<MVertex*> &touched){
if (touched.find(v) != touched.end()) return;
touched.insert(v);
+
+ count++;
for (std::multimap <MVertex*,MVertex*>::iterator it = v2v.lower_bound(v);
it != v2v.upper_bound(v) ; ++it){
+
MVertex *nextV = it->second;
-
- //change orientation
- ///------------------
- printf("*********************** Changing orientation \n");
-
- //compute dot product (N0,R0,A0)^T dot (Ni,Ri,Ai)
- //where N,R,A are the normal, radial and axial direction s
+ if (touched.find(nextV) == touched.end()){
+
+ //compute dot product (N0,R0,A0) dot (Ni,Ri,Ai)^T
+ //where N,R,A are the 3 directions
std::map<MVertex*, STensor3>::iterator iter = crossField.find(nextV);
STensor3 nextCross = iter->second;
- STensor3 crossT = cross.transpose();
- STensor3 prod = crossT.operator*=(nextCross);
- cross.print("cross ");
- nextCross.print("nextCross ");
- prod.print("product");
+ STensor3 nextCrossT = nextCross.transpose();
+ STensor3 prod = cross.operator*=(nextCrossT);
fullMatrix<double> mat(3,3); prod.getMat(mat);
-
+
//find biggest dot product
fullVector<int> Id(3);
- for (int i = 0; i < 3; i++){
+ Id(0) = Id(1) = Id(2) = 0;
+ for (int j = 0; j < 3; j++){
double maxVal = 0.0;
- for (int j = 0; j < 3; j++){
+ for (int i = 0; i < 3; i++){
double val = fabs(mat(i,j));
if( val > maxVal ){
maxVal = val;
- Id(i) = j;
+ Id(j) = i;
}
}
}
- // if (Id(0) +Id(1)+ Id(2) != 3 || (Id(0) == 1 && Id(1)==1 && Id(2)==1)) {
- // std::cout << "This should not happen: sum should be 0+1+2" << std::endl;
- // printf("Id =%d %d %d \n", Id(0), Id(1), Id(2));
- // Id(0) = 0; Id(1) = 1; Id(2) = 2;
- // //break;
- // //exit(1);
- // return;
- // }
-
- Id(0) = 0; Id(1) = 1; Id(2) = 2;
-
+ //check
+ if (Id(0) +Id(1)+ Id(2) != 3 || (Id(0) == 1 && Id(1)==1 && Id(2)==1)) {
+ std::cout << "This should not happen: sum should be 0+1+2" << std::endl;
+ printf("Id =%d %d %d \n", Id(0), Id(1), Id(2));
+ return;
+ }
+
//create new cross
- printf("Id =%d %d %d \n", Id(0), Id(1), Id(2));
fullMatrix<double> newmat(3,3);
for (int i = 0; i < 3; i++){
for (int j = 0; j < 3; j++){
- newmat(i,j) = nextCross(i,Id(j)) ;
+ newmat(i,j) = nextCross(Id(i),j) ;
}
}
STensor3 newcross(0.0);
newcross.setMat(newmat);
crossField[iter->first] = newcross;
- newcross.print("newcross");
+
+ //print info
+ printf("************** COUNT = %d \n", count);
+ if (Id(0) == 0 && Id(1) == 1 && Id(2) == 2)
+ printf("orient OK Id=%d %d %d\n", Id(0), Id(1), Id(2));
+ else{
+ printf("change orientation Id=%d %d %d \n", Id(0), Id(1), Id(2));
+ cross.print("cross ");
+ nextCross.print("nextCross ");
+ prod.print("product");
+ newcross.print("newcross");
+ }
+
+ fprintf(fi,"SP(%g,%g,%g) {%g};\n",nextV->x(),nextV->y(),nextV->z(), (double)count);
//continue recursion
- recur_connect_vert (nextV, newcross, v2v,touched);
+ recur_connect_vert (fi, count, nextV, newcross, v2v,touched);
+ }
+
}
}
@@ -977,8 +986,16 @@ void Frame_field::continuousCrossField(GRegion *gr, GFace *gf){
std::set<MVertex*> touched;
std::map<MVertex*, STensor3>::iterator iter = crossField.find(beginV);
STensor3 bCross = iter->second;
- recur_connect_vert (beginV,bCross,v2v,touched);
+
+ FILE *fi = fopen ("cross_recur.pos","w");
+ fprintf(fi,"View \"\"{\n");
+ fprintf(fi,"SP(%g,%g,%g) {%g};\n",beginV->x(),beginV->y(),beginV->z(), 0.0);
+ int count = 0;
+
+ recur_connect_vert (fi, count, beginV,bCross,v2v,touched);
+ fprintf(fi,"};\n");
+ fclose (fi);
printf("touched =%d vert =%d \n", touched.size(), vertex_to_vertices.size());
}
diff --git a/Mesh/directions3D.h b/Mesh/directions3D.h
index ce418ef83456a78b7ee3090dbb09199ad51d4394..a1f67511520df980b0fc3cec65bfd2016df3ab85 100644
--- a/Mesh/directions3D.h
+++ b/Mesh/directions3D.h
@@ -68,7 +68,7 @@ class Frame_field{
const std::string& filename);
static void saveCrossField(const std::string& filename, double scale, bool full=true);
static void continuousCrossField(GRegion *gr, GFace *gf);
- static void recur_connect_vert(MVertex *v,STensor3 &cross, std::multimap<MVertex*,MVertex*> &v2v, std::set<MVertex*> &touched);
+ static void recur_connect_vert(FILE*fi, int count, MVertex *v,STensor3 &cross, std::multimap<MVertex*,MVertex*> &v2v, std::set<MVertex*> &touched);
static void save_energy(GRegion* gr, const std::string& filename);
static void save_dist(const std::string& filename);
static void checkAnnData(GEntity* ge, const std::string& filename);
diff --git a/benchmarks/centerlines/cyl_centerlines.geo b/benchmarks/centerlines/cyl_centerlines.geo
index 02d517ee3c49f7ccfaf96770e1c9bc68cf805cdd..dbd4470cbad950118341aa77a0e316ce15fc1943 100644
--- a/benchmarks/centerlines/cyl_centerlines.geo
+++ b/benchmarks/centerlines/cyl_centerlines.geo
@@ -1,10 +1,10 @@
-Mesh.Algorithm = 8; //(1=MeshAdapt, 2=Automatic, 5=Delaunay, 6=Frontal, 7=bamg, 8=delquad, 9=packing)
-Mesh.Algorithm3D = 9; //(1=tetgen, 4=netgen, 7=mmg3D, 9=Rtree)
+Mesh.Algorithm = 7; //(1=MeshAdapt, 2=Automatic, 5=Delaunay, 6=Frontal, 7=bamg, 8=delquad, 9=packing)
+Mesh.Algorithm3D = 7; //(1=tetgen, 4=netgen, 7=mmg3D, 9=Rtree)
-Mesh.RecombineAll=1;
+//Mesh.RecombineAll=1;
//Mesh.Recombine3DAll=1;
-Mesh.Smoothing=0;
-Mesh.SmoothCrossField=0;
+//Mesh.Smoothing=0;
+//Mesh.SmoothCrossField=0;
Mesh.LcIntegrationPrecision = 1.e-3;
@@ -14,7 +14,7 @@ Merge "cylemi.stl";
Field[1] = Centerline;
Field[1].FileName = "centerlinesCYL.vtk";
-Field[1].nbPoints = 12;
+Field[1].nbPoints = 16;
Field[1].closeVolume =1;
Field[1].reMesh =1;