From 7e3df39b4a0a9f50d6507158ae1eee03cb6e3adb Mon Sep 17 00:00:00 2001
From: Emilie Marchandise <emilie.marchandise@uclouvain.be>
Date: Tue, 16 Apr 2013 09:45:18 +0000
Subject: [PATCH]
---
Common/Octree.cpp | 1 +
Geo/GFaceCompound.cpp | 135 ++++++++++--------
Geo/GFaceCompound.h | 2 +
Geo/GModel.cpp | 11 +-
Geo/GModelIO_GEO.cpp | 12 +-
Mesh/CenterlineField.cpp | 10 +-
Mesh/directions3D.cpp | 33 +++--
Mesh/meshGEdge.cpp | 5 +-
Mesh/meshGFaceBamg.cpp | 23 ---
Mesh/meshGFaceElliptic.cpp | 3 +-
Mesh/meshGRegion.cpp | 8 +-
Mesh/surfaceFiller.cpp | 16 +--
.../centerlines/carotid_centerlines.geo | 8 +-
.../centerlines/cerebral_centerlines.geo | 6 +-
benchmarks/centerlines/cyl_centerlines.geo | 5 +-
15 files changed, 142 insertions(+), 136 deletions(-)
diff --git a/Common/Octree.cpp b/Common/Octree.cpp
index 2751c8924c..783a1839bd 100644
--- a/Common/Octree.cpp
+++ b/Common/Octree.cpp
@@ -93,3 +93,4 @@ void Octree_SearchAll(double *pt, Octree *myOctree, std::vector<void*> *output)
searchAllElements(myOctree->root, pt, myOctree->info, myOctree->function_BB,
myOctree->function_inElement, output);
}
+
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index 40aa04d35e..aa9680097a 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -1196,7 +1196,6 @@ GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
getBoundingEdges();
-
_mapping = HARMONIC;
_type = UNITCIRCLE;
if(toc == RADIAL_BASIS) _mapping = RBF;
@@ -1219,7 +1218,9 @@ GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
fillTris.clear();
#if defined(HAVE_ANN)
+ kdtree = NULL;
uv_kdtree = NULL;
+ uv_nodes = NULL;
nodes = NULL;
index = new ANNidx[2];
dist = new ANNdist[2];
@@ -1273,6 +1274,8 @@ GFaceCompound::GFaceCompound(GModel *m, int tag, std::list<GFace*> &compound,
#if defined(HAVE_ANN)
uv_kdtree = NULL;
+ kdtree = NULL;
+ uv_nodes = NULL;
nodes = NULL;
index = new ANNidx[2];
dist = new ANNdist[2];
@@ -1294,6 +1297,8 @@ GFaceCompound::~GFaceCompound()
delete MONE;
#if defined(HAVE_ANN)
if(uv_kdtree) delete uv_kdtree;
+ if(kdtree) delete kdtree;
+ if(uv_nodes) annDeallocPts(uv_nodes);
if(nodes) annDeallocPts(nodes);
delete[]index;
delete[]dist;
@@ -1930,19 +1935,9 @@ SPoint2 GFaceCompound::parFromPoint(const SPoint3 &p, bool onSurface) const
GPoint GFaceCompound::pointInRemeshedOctree(double par1, double par2) const
{
-
- //if not meshed yet
- if (meshStatistics.status != GFace::DONE || triangles.size()+quadrangles.size() == 0) {
- GPoint gp (3,3,0,this);
- gp.setNoSuccess();
- return gp;
- }
-
+
//create new octree with new mesh elements
if(!octNew){
- //FILE * of = fopen("myOCTREE.pos","w");
- //fprintf(of, "View \"\"{\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]);
@@ -1971,12 +1966,6 @@ GPoint GFaceCompound::pointInRemeshedOctree(double par1, double par2) const
myParamElems.push_back(new MTriangle(news[0],news[1],news[2],e->getNum()));
else if (e->getType() == TYPE_QUA) {
myParamElems.push_back(new MQuadrangle(news[0],news[1],news[2],news[3],e->getNum()));
- // fprintf(of, "SQ(%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g,%g};\n",
- // news[0]->x(), news[0]->y(), news[0]->z(),
- // news[1]->x(), news[1]->y(), news[1]->z(),
- // news[2]->x(), news[2]->y(), news[2]->z(),
- // news[3]->x(), news[3]->y(), news[3]->z(),
- // (double)e->getNum(), (double)e->getNum(), (double)e->getNum(), (double)e->getNum());
}
}
@@ -1984,21 +1973,18 @@ GPoint GFaceCompound::pointInRemeshedOctree(double par1, double par2) const
//build kdtree boundary nodes in parametric space
#if defined(HAVE_ANN)
- nodes = annAllocPts(myBCNodes.size(), 3);
+ uv_nodes = annAllocPts(myBCNodes.size(), 3);
std::set<SPoint2>::iterator itp = myBCNodes.begin();
int ind = 0;
while (itp != myBCNodes.end()){
SPoint2 pt = *itp;
- //fprintf(of, "SP(%g,%g,%g){%g};\n", pt.x(), pt.y(), 0.0, 10000);
- nodes[ind][0] = pt.x();
- nodes[ind][1] = pt.y();
- nodes[ind][2] = 0.0;
+ uv_nodes[ind][0] = pt.x();
+ uv_nodes[ind][1] = pt.y();
+ uv_nodes[ind][2] = 0.0;
itp++; ind++;
}
- uv_kdtree = new ANNkd_tree(nodes, myBCNodes.size(), 3);
+ uv_kdtree = new ANNkd_tree(uv_nodes, myBCNodes.size(), 3);
#endif
- //fprintf(of,"};\n");
- //fclose(of);
}
//now use new octree to find point
@@ -2021,23 +2007,19 @@ GPoint GFaceCompound::pointInRemeshedOctree(double par1, double par2) const
gp.x() = e->interpolate(valX,UV[0],UV[1],UV[2]);
gp.y() = e->interpolate(valY,UV[0],UV[1],UV[2]);
gp.z() = e->interpolate(valZ,UV[0],UV[1],UV[2]);
- //printf("found point in new octree (UV=%g %g) %g %g %g (E=%D)\n",par1, par2, gp.x(), gp.y(), gp.z(), e->getNum());
return gp;
}
//if element not found in new octree find closest point
else{
- //printf("point not found on octnew --> look closest point with kdtree \n");
GPoint gp(50,50,50);
#if defined(HAVE_ANN)
double pt[3] = {par1,par2,0.0};
uv_kdtree->annkSearch(pt, 2, index, dist);
- SPoint3 p1(nodes[index[0]][0], nodes[index[0]][1], nodes[index[0]][2]);
- SPoint3 p2(nodes[index[1]][0], nodes[index[1]][1], nodes[index[1]][2]);
+ SPoint3 p1(uv_nodes[index[0]][0], uv_nodes[index[0]][1], uv_nodes[index[0]][2]);
+ SPoint3 p2(uv_nodes[index[1]][0], uv_nodes[index[1]][1], uv_nodes[index[1]][2]);
SPoint3 pnew; double d;
signedDistancePointLine(p1,p2, SPoint3(par1,par2,0.0), d, pnew);
- //printf("p1=%g %g p2=%g %g \n", p1.x(), p1.y(), p2.x(), p2.y());
- //printf("UV=%g %g UVnew =%g %g \n", par1,par2, pnew.x(), pnew.y());
-
+
double uvw[3]={pnew.x(),pnew.y(), 0.0};
double UV[3];
MElement *e = octNew->find(pnew.x(), pnew.y(), 0.0,-1,true);
@@ -2054,10 +2036,8 @@ GPoint GFaceCompound::pointInRemeshedOctree(double par1, double par2) const
gp.x() = e->interpolate(valX,UV[0],UV[1],UV[2]);
gp.y() = e->interpolate(valY,UV[0],UV[1],UV[2]);
gp.z() = e->interpolate(valZ,UV[0],UV[1],UV[2]);
- //printf("found closest point (UV=%g %g) %g %g %g \n",pnew.x(), pnew.y(), gp.y(), gp.z());
}
else{
- Msg::Error("Point not found in kdtree");
gp.setNoSuccess();
}
#else
@@ -2086,24 +2066,42 @@ GPoint GFaceCompound::point(double par1, double par2) const
double par[2] = {par1,par2};
GFaceCompoundTriangle *lt;
getTriangle(par1, par2, <, U,V);
- if(!lt && _mapping != RBF){
- printf("ARRG POINT NOT FOUND--> should improve octree search \n");
- //exit(1);
- //printf("POINT no success %d tris %d quad \n", triangles.size(), quadrangles.size());
- GPoint gp = pointInRemeshedOctree(par1,par2);
- gp.setNoSuccess();
- return gp;
- }
- else if (!lt && _mapping == RBF){
- if (fabs(par1) > 1 || fabs(par2) > 1){
- GPoint gp (0,0,0,this);
+ if(!lt){
+ //printf("POINT (%g %g) NOT FOUND --> find closest \n", par1,par2);
+ if (meshStatistics.status != GFace::DONE ) {
+ double pt[3] = {par1,par2, 0.0};
+ kdtree->annkSearch(pt, 1, index, dist);
+ lt = &(_gfct[index[0]]);
+
+ SPoint3 pnew_a, pnew_b, pnew_c; double d_a, d_b, d_c;
+ signedDistancePointLine(lt->p1,lt->p2, SPoint3(par1,par2,0.0), d_a, pnew_a);
+ signedDistancePointLine(lt->p1,lt->p3, SPoint3(par1,par2,0.0), d_b, pnew_b);
+ signedDistancePointLine(lt->p2,lt->p3, SPoint3(par1,par2,0.0), d_c, pnew_c);
+ double u,v;
+ if (d_a <= d_b && d_a <= d_c) { u = pnew_a.x(); v= pnew_a.y();}
+ else if (d_b <= d_a && d_b <= d_c) { u = pnew_b.x(); v= pnew_b.y();}
+ else { u = pnew_c.x(); v= pnew_c.y();}
+ double M[2][2],X[2],R[2];
+ const SPoint3 p0 = (lt)->p1;
+ const SPoint3 p1 = (lt)->p2;
+ const SPoint3 p2 = (lt)->p3;
+ M[0][0] = p1.x() - p0.x();
+ M[0][1] = p2.x() - p0.x();
+ M[1][0] = p1.y() - p0.y();
+ M[1][1] = p2.y() - p0.y();
+ R[0] = (u - p0.x());
+ R[1] = (v - p0.y());
+ sys2x2(M, R, X);
+ U = X[0];
+ 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;
}
- double x, y, z;
- SVector3 dXdu, dXdv;
- _rbf->UVStoXYZ(par1, par2,x,y,z, dXdu, dXdv);
- return GPoint(x,y,z);
}
if (lt->gf->geomType() != GEntity::DiscreteSurface){
@@ -2178,13 +2176,16 @@ Pair<SVector3,SVector3> GFaceCompound::firstDer(const SPoint2 ¶m) const
double U,V;
GFaceCompoundTriangle *lt;
getTriangle(param.x(), param.y(), <, U,V);
- if(!lt){
- printf("ARRG FIRSTDER POINT NOT FOUND --> should improve octree search \n");
- exit(1);
- return Pair<SVector3, SVector3>(SVector3(1, 0, 0), SVector3(0, 1, 0));
+ MTriangle *tri;
+ if (lt) tri = lt->tri;
+ else {
+ //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;
}
- MTriangle *tri = lt->tri;
SVector3 dXdu1 = firstDerivatives[tri->getVertex(0)].first();
SVector3 dXdu2 = firstDerivatives[tri->getVertex(1)].first();
SVector3 dXdu3 = firstDerivatives[tri->getVertex(2)].first();
@@ -2370,6 +2371,7 @@ void GFaceCompound::getTriangle(double u, double v,
*lt = (GFaceCompoundTriangle*)Octree_Search(uv, oct);
if(!(*lt)){
+ _u = 0.0; _v = 0.0;
return;
}
@@ -2407,11 +2409,12 @@ void GFaceCompound::buildOct() const
}
}
- // make bounding box larger up to (absolute) geometrical tolerance
- double eps = 0.0; //CTX::instance()->geom.tolerance;
- SPoint3 bbmin = bb.min(), bbmax = bb.max(), bbeps(eps, eps, eps);
- bbmin -= bbeps;
- bbmax += bbeps;
+ //ANN octree
+ std::set<MVertex*> allVS;
+ nodes = annAllocPts(count, 3);
+
+ // make bounding box
+ SPoint3 bbmin = bb.min(), bbmax = bb.max();
double origin[3] = {bbmin.x(), bbmin.y(), bbmin.z()};
double ssize[3] = {bbmax.x() - bbmin.x(),
bbmax.y() - bbmin.y(),
@@ -2486,6 +2489,11 @@ void GFaceCompound::buildOct() const
SVector3 dXdv(dXdxi * inv[0][1] + dXdeta * inv[1][1]);
firstElemDerivatives[(MElement*)t] = Pair<SVector3,SVector3>(dXdu,dXdv);
+ // build ANN kdtree
+ nodes[count][0] = (it0->second.x() + it1->second.x() + it2->second.x())/3.0 ;
+ nodes[count][1] = (it0->second.y() + it1->second.y() + it2->second.y())/3.0 ;
+ nodes[count][2] = 0.0;
+
Octree_Insert(&_gfct[count], oct);
count++;
}
@@ -2493,6 +2501,7 @@ void GFaceCompound::buildOct() const
nbT = count;
Octree_Arrange(oct);
+
//smooth first derivatives at vertices
if(adjv.size() == 0){
std::vector<MTriangle*> allTri;
@@ -2516,6 +2525,9 @@ void GFaceCompound::buildOct() const
firstDerivatives[v] = Pair<SVector3, SVector3>(dXdu, dXdv);
}
+ //build ANN kdtree
+ kdtree = new ANNkd_tree(nodes, count, 3);
+
printStuff();
}
@@ -3028,9 +3040,6 @@ GPoint GFaceCompound::intersectionWithCircle(const SVector3 &n1, const SVector3
ct->p1 * ( 1.-uv[0]-uv[1] ) +
ct->p2 *uv[0] +
ct->p3 *uv[1] ;
- // GPoint ttt = point(pp.x(),pp.y());
- // printf("%g %g %g vs %g %g %g\n",ttt.x(),ttt.y(),ttt.z(),s.x(),s.y(),s.z());
- // printf("%d/%d\n",i,nbT);
return GPoint(s.x(),s.y(),s.z(),this,pp);
}
}
diff --git a/Geo/GFaceCompound.h b/Geo/GFaceCompound.h
index da913da8b4..225886c99e 100644
--- a/Geo/GFaceCompound.h
+++ b/Geo/GFaceCompound.h
@@ -103,6 +103,8 @@ class GFaceCompound : public GFace {
linearSystem <double> *_lsys;
#if defined(HAVE_ANN)
mutable ANNkd_tree *uv_kdtree;
+ mutable ANNkd_tree *kdtree;
+ mutable ANNpointArray uv_nodes;
mutable ANNpointArray nodes;
ANNidxArray index;
ANNdistArray dist;
diff --git a/Geo/GModel.cpp b/Geo/GModel.cpp
index fc52c849b1..d6c4b63659 100644
--- a/Geo/GModel.cpp
+++ b/Geo/GModel.cpp
@@ -2339,9 +2339,6 @@ GEdge* GModel::addCompoundEdge(std::vector<GEdge*> edges, int num){
Curve *c = Create_Curve(num, MSH_SEGM_COMPOUND, 1, NULL, NULL, -1, -1, 0., 1.);
for(unsigned int i= 0; i < edges.size(); i++)
c->compound.push_back(edges[i]->tag());
-
- // Curve *c = Create_Curve(num, MSH_SEGM_DISCRETE, 1,
- // NULL, NULL, -1, -1, 0., 1.);
List_T *points = Tree2List(getGEOInternals()->Points);
GVertex *gvb = gec->getBeginVertex();
GVertex *gve = gec->getEndVertex();
@@ -2363,6 +2360,14 @@ GEdge* GModel::addCompoundEdge(std::vector<GEdge*> edges, int num){
End_Curve(c);
Tree_Add(getGEOInternals()->Curves, &c);
CreateReversedCurve(c);
+
+ // c->Method = gec->meshAttributes.method;
+ // c->nbPointsTransfinite = gec->meshAttributes.nbPointsTransfinite;
+ // c->typeTransfinite = gec->meshAttributes.typeTransfinite ;
+ // c->coeffTransfinite = gec->meshAttributes.coeffTransfinite ;
+ // c->Extrude = gec->meshAttributes.extrude ;
+ // c->ReverseMesh = gec->meshAttributes.reverseMesh ;
+
}
return gec;
diff --git a/Geo/GModelIO_GEO.cpp b/Geo/GModelIO_GEO.cpp
index cb16f2cf6e..f048e1b267 100644
--- a/Geo/GModelIO_GEO.cpp
+++ b/Geo/GModelIO_GEO.cpp
@@ -136,7 +136,7 @@ int GModel::importGEOInternals()
Curve *c;
List_Read(curves, i, &c);
if(c->Num >= 0){
- GEdge *e = getEdgeByTag(c->Num);
+ GEdge *e = getEdgeByTag(c->Num);
if(!e && c->Typ == MSH_SEGM_COMPOUND){
std::vector<GEdge*> comp;
for(unsigned int j = 0; j < c->compound.size(); j++){
@@ -144,6 +144,12 @@ int GModel::importGEOInternals()
if(ge) comp.push_back(ge);
}
e = new GEdgeCompound(this, c->Num, comp);
+ e->meshAttributes.method = c->Method;
+ e->meshAttributes.nbPointsTransfinite = c->nbPointsTransfinite;
+ e->meshAttributes.typeTransfinite = c->typeTransfinite;
+ e->meshAttributes.coeffTransfinite = c->coeffTransfinite;
+ e->meshAttributes.extrude = c->Extrude;
+ e->meshAttributes.reverseMesh = c->ReverseMesh;
add(e);
}
else if(!e && c->beg && c->end){
@@ -156,8 +162,8 @@ int GModel::importGEOInternals()
e = new gmshEdge(this, c, 0, 0);
add(e);
}
- else
- e->resetMeshAttributes();
+
+
if(!c->Visible) e->setVisibility(0);
if(c->Color.type) e->setColor(c->Color.mesh);
if(c->degenerated) {
diff --git a/Mesh/CenterlineField.cpp b/Mesh/CenterlineField.cpp
index 445c732c30..c4bfd84f5b 100644
--- a/Mesh/CenterlineField.cpp
+++ b/Mesh/CenterlineField.cpp
@@ -637,9 +637,11 @@ void Centerline::createSplitCompounds()
int num_gec = NE+i+1;
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;
+ 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;
}
// Parametrize Compound surfaces
@@ -892,7 +894,7 @@ void Centerline::cutMesh()
double AR = L/D;
// printf("*** Centerline branch %d (AR=%.1f) \n", edges[i].tag, AR);
- int nbSplit = (int)floor(AR/2 + 1.9); //AR/2 + 0.9
+ int nbSplit = (int)floor(AR/2 + 0.9); //AR/2 + 0.9
if( nbSplit > 1 ){
//printf("->> cut branch in %d parts \n", nbSplit);
double li = L/nbSplit;
diff --git a/Mesh/directions3D.cpp b/Mesh/directions3D.cpp
index 3356af0dc8..07456ece99 100644
--- a/Mesh/directions3D.cpp
+++ b/Mesh/directions3D.cpp
@@ -879,8 +879,6 @@ void Frame_field::recur_connect_vert(MVertex *v,
if (touched.find(v) != touched.end()) return;
touched.insert(v);
- //std::map<MVertex*, STensor3>::iterator iterv = crossField.find(v);
- //STensor3 cross = iterv->second;
for (std::multimap <MVertex*,MVertex*>::iterator it = v2v.lower_bound(v);
it != v2v.upper_bound(v) ; ++it){
@@ -914,22 +912,27 @@ void Frame_field::recur_connect_vert(MVertex *v,
}
}
- 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));
- //exit(1);
- return;
- }
+ // 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;
//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(Id(j),i);
+ newmat(i,j) = nextCross(i,Id(j)) ;
}
}
- STensor3 newcross;
+ STensor3 newcross(0.0);
newcross.setMat(newmat);
crossField[iter->first] = newcross;
newcross.print("newcross");
@@ -975,6 +978,8 @@ void Frame_field::continuousCrossField(GRegion *gr, GFace *gf){
std::map<MVertex*, STensor3>::iterator iter = crossField.find(beginV);
STensor3 bCross = iter->second;
recur_connect_vert (beginV,bCross,v2v,touched);
+
+ printf("touched =%d vert =%d \n", touched.size(), vertex_to_vertices.size());
}
@@ -988,8 +993,8 @@ void Frame_field::saveCrossField(const std::string& filename, double scale, bool
it != crossField.end(); it++){
SPoint3 p = it->first->point();
STensor3 m = it->second;
- double val1 = eulerAngleFromQtn(cross3D(m).rotationTo(origin))*180./M_PI, val2=val1;
- //double val1=1.0, val2=1.0;
+ //double val1 = eulerAngleFromQtn(cross3D(m).rotationTo(origin))*180./M_PI, val2=val1;
+ double val1=1.0, val2=1.0;
p1 = SPoint3(p.x() + k*m.get_m11(),
p.y() + k*m.get_m21(),
p.z() + k*m.get_m31());
@@ -998,7 +1003,7 @@ void Frame_field::saveCrossField(const std::string& filename, double scale, bool
p.y() - k*m.get_m21(),
p.z() - k*m.get_m31());
if(full) print_segment(p,p1,val1,val2,file);
- //val1=2.0; val2=2.0;
+ val1=2.0; val2=2.0;
p1 = SPoint3(p.x() + k*m.get_m12(),
p.y() + k*m.get_m22(),
p.z() + k*m.get_m32());
@@ -1007,7 +1012,7 @@ void Frame_field::saveCrossField(const std::string& filename, double scale, bool
p.y() - k*m.get_m22(),
p.z() - k*m.get_m32());
if(full) print_segment(p,p1,val1,val2,file);
- //val1=3.0; val2=3.0;
+ val1=3.0; val2=3.0;
p1 = SPoint3(p.x() + k*m.get_m13(),
p.y() + k*m.get_m23(),
p.z() + k*m.get_m33());
diff --git a/Mesh/meshGEdge.cpp b/Mesh/meshGEdge.cpp
index 00a719728d..3bcd86a9d3 100644
--- a/Mesh/meshGEdge.cpp
+++ b/Mesh/meshGEdge.cpp
@@ -379,7 +379,7 @@ void meshGEdge::operator() (GEdge *ge)
CTX::instance()->mesh.lcIntegrationPrecision);
}
else{
- a = Integration(ge, t_begin, t_end, F_Lc, Points,
+ a = Integration(ge, t_begin, t_end, F_Lc, Points,
CTX::instance()->mesh.lcIntegrationPrecision);
}
@@ -391,7 +391,8 @@ void meshGEdge::operator() (GEdge *ge)
}
a = smoothPrimitive(ge, sqrt(CTX::instance()->mesh.smoothRatio), Points);
N = std::max(ge->minimumMeshSegments() + 1, (int)(a + 1.99));
- }
+
+}
// force odd number of points if blossom is used for recombination
if(ge->meshAttributes.method != MESH_TRANSFINITE &&
diff --git a/Mesh/meshGFaceBamg.cpp b/Mesh/meshGFaceBamg.cpp
index 6e4b3fbd11..0a41d9045b 100644
--- a/Mesh/meshGFaceBamg.cpp
+++ b/Mesh/meshGFaceBamg.cpp
@@ -190,29 +190,6 @@ void meshGFaceBamg(GFace *gf){
MElementOctree *_octree = NULL;
if (hasCompounds){
_octree = new MElementOctree(myParamElems);
- // //EMI PRINT TRIS
- // FILE * fi = fopen("TRIS.pos","w");
- // fprintf(fi, "View \"\"{\n");
- // for( int i =0; i< gf->triangles.size(); i++){
- // int nodes [3] = {gf->triangles[i]->getVertex(0)->getIndex(),
- // gf->triangles[i]->getVertex(1)->getIndex(),
- // gf->triangles[i]->getVertex(2)->getIndex()};
- // double u1(bamgVertices[nodes[0]][0]);
- // double u2(bamgVertices[nodes[1]][0]);
- // double u3(bamgVertices[nodes[2]][0]);
- // double v1(bamgVertices[nodes[0]][1]);
- // double v2(bamgVertices[nodes[1]][1]);
- // double v3(bamgVertices[nodes[2]][1]);
- // fprintf(fi, "ST(%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,"
- // "%22.15E,%22.15E){%d,%d,%d};\n",
- // u1, v1, 0.0,
- // u2, v2, 0.0,
- // u3, v3, 0.0,
- // i, i, i);
- // }
- // fprintf(fi,"};\n");
- // fclose(fi);
- // //END EMI PRINT TRIS
}
Mesh2 *refinedBamgMesh = 0;
diff --git a/Mesh/meshGFaceElliptic.cpp b/Mesh/meshGFaceElliptic.cpp
index 8c837dbadb..6571842881 100644
--- a/Mesh/meshGFaceElliptic.cpp
+++ b/Mesh/meshGFaceElliptic.cpp
@@ -611,7 +611,8 @@ static bool computeRingVertices(GFace *gf, Centerline *center,
N = ge1->mesh_vertices.size() + 1;
int N2 = ge2->mesh_vertices.size() + 1;
if (N != N2 || N%2 != 0){
- Msg::Error("You should have an equal pair number of points in centerline field N=%d N2=%d ", N, N2);
+ Msg::Error("You should an even number nbPoints in centerline =%d \n", N);
+ exit(1);
return false;
}
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 9a6a454942..4e78b3f24e 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -635,15 +635,15 @@ void MeshDelaunayVolume(std::vector<GRegion*> ®ions)
insertVerticesInRegion(gr);
}
- // //emi test frame field
+ //emi test frame field
// int NumSmooth = 10;//CTX::instance()->mesh.smoothCrossField
// std::cout << "NumSmooth = " << NumSmooth << std::endl;
// if(NumSmooth && (gr->dim() == 3)){
// double scale = gr->bounds().diag()*1e-2;
// Frame_field::initRegion(gr,NumSmooth);
// Frame_field::saveCrossField("cross0.pos",scale);
- // Frame_field::smoothRegion(gr,NumSmooth);
- // Frame_field::saveCrossField("cross1.pos",scale);
+ // //Frame_field::smoothRegion(gr,NumSmooth);
+ // //Frame_field::saveCrossField("cross1.pos",scale);
// GFace *gf = GModel::current()->getFaceByTag(2);
// Frame_field::continuousCrossField(gr,gf);
// Frame_field::saveCrossField("cross2.pos",scale);
@@ -651,7 +651,7 @@ void MeshDelaunayVolume(std::vector<GRegion*> ®ions)
// Frame_field::init_region(gr);
// Frame_field::clear();
// exit(1);
- // //fin test emi
+ //fin test emi
if (sqr.buildPyramids (gr->model())){
// relocate vertices if pyramids
diff --git a/Mesh/surfaceFiller.cpp b/Mesh/surfaceFiller.cpp
index 2780bd5c1a..f547098ba3 100644
--- a/Mesh/surfaceFiller.cpp
+++ b/Mesh/surfaceFiller.cpp
@@ -192,13 +192,13 @@ bool compute4neighbors (GFace *gf, // the surface
FieldManager *fields = gf->model()->getFields();
if(fields->getBackgroundField() > 0){
Field *f = fields->get(fields->getBackgroundField());
- if (!f->isotropic()){
- (*f)(v_center->x(),v_center->y(),v_center->z(), metricField,gf);
- }
- else {
+ if (!f->isotropic()){
+ (*f)(v_center->x(),v_center->y(),v_center->z(), metricField,gf);
+ }
+ else {
L = (*f)(v_center->x(),v_center->y(),v_center->z(), gf);
metricField = SMetric3(1./(L*L));
- }
+ }
}
@@ -229,8 +229,6 @@ bool compute4neighbors (GFace *gf, // the surface
// normalize vector t1 that is tangent to gf at midpoint
SVector3 t1 = basis_u * cos (quadAngle) + basis_v * sin (quadAngle) ;
t1.normalize();
-
- // printf("%d %g %g %g -- %g %g %g\n",gf->tag(),s1.x(),s1.y(),s1.z(),t1.x(),t1.y(),t1.z());
// compute the second direction t2 and normalize (t1,t2,n) is the tangent frame
SVector3 t2 = crossprod(n,t1);
@@ -239,9 +237,7 @@ bool compute4neighbors (GFace *gf, // the surface
if (DIR == 0 && crossf)fprintf(crossf,"VP(%g,%g,%g) {%g,%g,%g};\n",v_center->x(),v_center->y(),v_center->z(),t2.x(),t2.y(),t2.z());
if (DIR == 0 && crossf)fprintf(crossf,"VP(%g,%g,%g) {%g,%g,%g};\n",v_center->x(),v_center->y(),v_center->z(),-t1.x(),-t1.y(),-t1.z());
if (DIR == 0 && crossf)fprintf(crossf,"VP(%g,%g,%g) {%g,%g,%g};\n",v_center->x(),v_center->y(),v_center->z(),-t2.x(),-t2.y(),-t2.z());
- // if (DIR == 0 && crossf)fprintf(crossf,"VP(%g,%g,%g) {%g,%g,%g};\n",v_center->x(),v_center->y(),v_center->z(),s1.x(),s1.y(),s1.z());
- // if (DIR == 0 && crossf)fprintf(crossf,"VP(%g,%g,%g) {%g,%g,%g};\n",v_center->x(),v_center->y(),v_center->z(),s2.x(),s2.y(),s2.z());
-
+
double size_1 = sqrt(1. / dot(t1,metricField,t1));
double size_2 = sqrt(1. / dot(t2,metricField,t2));
diff --git a/benchmarks/centerlines/carotid_centerlines.geo b/benchmarks/centerlines/carotid_centerlines.geo
index f3fa2d9754..8b4816e93f 100644
--- a/benchmarks/centerlines/carotid_centerlines.geo
+++ b/benchmarks/centerlines/carotid_centerlines.geo
@@ -1,19 +1,19 @@
-Mesh.Algorithm = 8; //(1=MeshAdapt, 2=Automatic, 5=Delaunay, 6=Frontal, 7=bamg, 8=delquad)
+Mesh.Algorithm = 9; //(1=MeshAdapt, 2=Automatic, 5=Delaunay, 6=Frontal, 7=bamg, 8=delquad)
Mesh.Algorithm3D = 1;//(1=tetgen, 4=netgen, 5=FrontalDel, 6=FrontalHex, 7=MMG3D, 9=R-tree
//Mesh.SmoothCrossField = 20;
-Mesh.Smoothing=0;
+//Mesh.Smoothing=0;
Mesh.LcIntegrationPrecision = 1.e-5;
Mesh.RecombineAll = 1;
-//Mesh.Bunin = 100;
+//Mesh.Bunin = 20;
Merge "carotid.stl";
Field[1] = Centerline;
Field[1].FileName = "centerlinesCAROTID.vtk";
-Field[1].nbPoints = 20;
+Field[1].nbPoints = 30;
Field[1].nbElemLayer = 4;
Field[1].hLayer = 0.2; //percent of vessel radius
diff --git a/benchmarks/centerlines/cerebral_centerlines.geo b/benchmarks/centerlines/cerebral_centerlines.geo
index da6689dc35..b1e260d915 100644
--- a/benchmarks/centerlines/cerebral_centerlines.geo
+++ b/benchmarks/centerlines/cerebral_centerlines.geo
@@ -1,4 +1,4 @@
-Mesh.Algorithm = 6; //(1=MeshAdapt, 2=Automatic, 5=Delaunay, 6=Frontal, 7=bamg, 8=delquad)
+Mesh.Algorithm = 8; //(1=MeshAdapt, 2=Automatic, 5=Delaunay, 6=Frontal, 7=bamg, 8=delquad)
//Mesh.Algorithm3D = 7; //(1=tetgen, 4=netgen, 7=mmg3D
Mesh.LcIntegrationPrecision = 1.e-2;
@@ -6,11 +6,11 @@ Mesh.LcIntegrationPrecision = 1.e-2;
Merge "cerebral.stl";
Mesh.RecombineAll = 1;
-Mesh.Bunin = 50;
+//Mesh.Bunin = 50;
Field[1] = Centerline;
Field[1].FileName = "centerlinesCEREBRAL.vtk";
-Field[1].nbPoints = 21;
+Field[1].nbPoints = 20;
Field[1].nbElemLayer = 4;
Field[1].hLayer = 0.2;//percent of vessel radius
diff --git a/benchmarks/centerlines/cyl_centerlines.geo b/benchmarks/centerlines/cyl_centerlines.geo
index 5bc1fc6cbf..02d517ee3c 100644
--- a/benchmarks/centerlines/cyl_centerlines.geo
+++ b/benchmarks/centerlines/cyl_centerlines.geo
@@ -2,8 +2,9 @@ Mesh.Algorithm = 8; //(1=MeshAdapt, 2=Automatic, 5=Delaunay, 6=Frontal, 7=bamg,
Mesh.Algorithm3D = 9; //(1=tetgen, 4=netgen, 7=mmg3D, 9=Rtree)
Mesh.RecombineAll=1;
-Mesh.Recombine3DAll=1;
+//Mesh.Recombine3DAll=1;
Mesh.Smoothing=0;
+Mesh.SmoothCrossField=0;
Mesh.LcIntegrationPrecision = 1.e-3;
@@ -13,7 +14,7 @@ Merge "cylemi.stl";
Field[1] = Centerline;
Field[1].FileName = "centerlinesCYL.vtk";
-Field[1].nbPoints = 35;
+Field[1].nbPoints = 12;
Field[1].closeVolume =1;
Field[1].reMesh =1;
--
GitLab