diff --git a/Fltk/statisticsWindow.cpp b/Fltk/statisticsWindow.cpp
index beb42ddf23d994cef53e03ecc1a89fbe207f4c2b..e7506606ad7595ddc17724b59af9c05218d5f112 100644
--- a/Fltk/statisticsWindow.cpp
+++ b/Fltk/statisticsWindow.cpp
@@ -201,183 +201,184 @@ void statisticsWindow::compute(bool elementQuality)
{
#if 0
{
- // MINIMUM ANGLES
- double minAngle = 120.0;
+ // MINIMUM MAXIMUM ANGLES
+ double minAngle = 1.0; //M_PI;
double meanAngle = 0.0;
int count = 0;
std::vector<GEntity*> entities;
GModel::current()->getEntities(entities);
std::map<int, std::vector<double> > d;
for(unsigned int i = 0; i < entities.size(); i++){
- if(entities[i]->dim() < 2) continue;
- for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){
- MElement *e = entities[i]->getMeshElement(j);
- double angle = e->angleShapeMeasure();
- minAngle = std::min(minAngle, angle);
- meanAngle += angle;
- count++;
+ if(entities[i]->dim() == 3) {// continue;//<3
+ for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){
+ MElement *e = entities[i]->getMeshElement(j);
+ double angle = e->angleShapeMeasure();
+ minAngle = std::min(minAngle, angle);
+ meanAngle += angle;
+ count++;
+ }
}
}
meanAngle = meanAngle / count;
- printf("Angles = min=%g av=%g \n", minAngle, meanAngle);
+ printf("Angles min =%g av=%g nbhex=%d\n", minAngle, meanAngle, count);
}
- {
- // MESH DEGREE VERTICES
- std::vector<GEntity*> entities;
- std::set<MEdge, Less_Edge> edges;
- GModel::current()->getEntities(entities);
- std::map<MVertex*, int > vert2Deg;
- for(unsigned int i = 0; i < entities.size(); i++){
- if(entities[i]->dim() < 2 ) continue;
- // if(entities[i]->tag() < 100) continue;
- for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){
- MElement *e = entities[i]->getMeshElement(j);
- for(unsigned int k = 0; k < e->getNumEdges(); k++){
- edges.insert(e->getEdge(k));
- }
- for(unsigned int k = 0; k < e->getNumVertices(); k++){
- MVertex *v = e->getVertex(k);
- if (v->onWhat()->dim() < 2) continue;
- std::map<MVertex*, int >::iterator it = vert2Deg.find(v);
- if (it == vert2Deg.end()) {
- vert2Deg.insert(std::make_pair(v,1));
- }
- else{
- int nbE = it->second+1;
- it->second = nbE;
- }
- }
- }
- }
- int dMin = 10;
- int dMax = 0;
- int d4 = 0;
- int nbElems = vert2Deg.size();
- std::map<MVertex*, int >::const_iterator itmap = vert2Deg.begin();
- for(; itmap !=vert2Deg.end(); itmap++){
- MVertex *v = itmap->first;
- int nbE = itmap->second;
- dMin = std::min(nbE, dMin);
- dMax = std::max(nbE, dMax);
- if (nbE == 4) d4 += 1;
- }
- if (nbElems > 0)
- printf("Stats degree vertices: dMin=%d , dMax=%d, d4=%g \n",
- dMin, dMax, (double)d4/nbElems);
- FieldManager *fields = GModel::current()->getFields();
- Field *f = fields->get(fields->background_field);
- int nbEdges = edges.size();
- printf("nb edges =%d \n", nbEdges);
- if(system("rm qualEdges.txt"));
- FILE *fp = fopen("qualEdges.txt", "w");
- std::vector<int> qualE;
- int nbS = 50;
- qualE.resize(nbS);
- if(fields.getBackgroundField() > 0){
- std::set<MEdge, Less_Edge>::iterator it = edges.begin();
- double sum = 0;
- for (; it !=edges.end();++it){
- MVertex *v0 = it->getVertex(0);
- MVertex *v1 = it->getVertex(1);
- double l = 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()));
- double lf = (*f)(0.5*(v0->x()+v1->x()), 0.5*(v0->y()+v1->y()),
- 0.5*(v0->z()+v1->z()),v0->onWhat());
- double el = l/lf;
- int index = (int) ceil(el*nbS*0.5);
- qualE[index]+= 1;
- double e = (l>lf) ? lf/l : l/lf;
- sum += e - 1.0;
- }
- double tau = exp ((1./edges.size()) * sum);
- printf("N edges = %d tau = %g\n",(int)edges.size(),tau);
-
- double ibegin = 2./(2*nbS);
- double inext = 2./nbS;
- for (int i= 0; i< qualE.size(); i++){
- fprintf(fp, "0 0 0 0 %g 0 0 %g \n", ibegin+i*inext , (double)qualE[i]/nbEdges);
- }
-
- }
- fclose(fp);
- }
- {
- std::vector<GEntity*> entities;
- std::set<MEdge, Less_Edge> edges;
- GModel::current()->getEntities(entities);
- std::map<MVertex*, int > vert2Deg;
- for(unsigned int i = 0; i < entities.size(); i++){
- if(entities[i]->dim() < 2 ) continue;
- if(entities[i]->tag() != 10) continue;
- for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){
- MElement *e = entities[i]->getMeshElement(j);
- for(unsigned int k = 0; k < e->getNumEdges(); k++){
- edges.insert(e->getEdge(k));
- }
- for(unsigned int k = 0; k < e->getNumVertices(); k++){
- MVertex *v = e->getVertex(k);
- if (v->onWhat()->dim() < 2) continue;
- std::map<MVertex*, int >::iterator it = vert2Deg.find(v);
- if (it == vert2Deg.end()){
- vert2Deg.insert(std::make_pair(v,1));
- }
- else{
- int nbE = it->second+1;
- it->second = nbE;
- }
- }
- }
- }
- int dMin = 10;
- int dMax = 0;
- int d4 = 0;
- int nbElems = vert2Deg.size();
- std::map<MVertex*, int >::const_iterator itmap = vert2Deg.begin();
- for(; itmap !=vert2Deg.end(); itmap++){
- MVertex *v = itmap->first;
- int nbE = itmap->second;
- dMin = std::min(nbE, dMin);
- dMax = std::max(nbE, dMax);
- if (nbE == 4) d4 += 1;
- }
- if (nbElems > 0) printf("Stats degree vertices: dMin=%d , dMax=%d, d4=%g \n",
- dMin, dMax, (double)d4/nbElems);
- FieldManager *fields = GModel::current()->getFields();
- Field *f = fields->get(fields.getBackgroundField());
- int nbEdges = edges.size();
- if(system("rm qualEdges.txt"));
- FILE *fp = fopen("qualEdges.txt", "w");
- std::vector<int> qualE;
- int nbS = 50;
- qualE.resize(nbS);
- if(fields.getBackgroundField() > 0){
- std::set<MEdge, Less_Edge>::iterator it = edges.begin();
- double sum = 0;
- for (; it !=edges.end();++it){
- MVertex *v0 = it->getVertex(0);
- MVertex *v1 = it->getVertex(1);
- double l = 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()));
- double lf = (*f)(0.5*(v0->x()+v1->x()), 0.5*(v0->y()+v1->y()),
- 0.5*(v0->z()+v1->z()),v0->onWhat());
- double el = l/lf;
- int index = (int) ceil(el*nbS*0.5);
- qualE[index]+= 1;
- double e = (l>lf) ? lf/l : l/lf;
- sum += e - 1.0;
- }
- double tau = exp ((1./edges.size()) * sum);
- double ibegin = 2./(2*nbS);
- double inext = 2./nbS;
- for (int i= 0; i< qualE.size(); i++){
- fprintf(fp, "0 0 0 0 %g 0 0 %g \n", ibegin+i*inext , (double)qualE[i]/nbEdges);
- }
- }
- fclose(fp);
- }
+ // {
+ // // MESH DEGREE VERTICES
+ // std::vector<GEntity*> entities;
+ // std::set<MEdge, Less_Edge> edges;
+ // GModel::current()->getEntities(entities);
+ // std::map<MVertex*, int > vert2Deg;
+ // for(unsigned int i = 0; i < entities.size(); i++){
+ // if(entities[i]->dim() < 2 ) continue;
+ // // if(entities[i]->tag() < 100) continue;
+ // for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){
+ // MElement *e = entities[i]->getMeshElement(j);
+ // for(unsigned int k = 0; k < e->getNumEdges(); k++){
+ // edges.insert(e->getEdge(k));
+ // }
+ // for(unsigned int k = 0; k < e->getNumVertices(); k++){
+ // MVertex *v = e->getVertex(k);
+ // if (v->onWhat()->dim() < 2) continue;
+ // std::map<MVertex*, int >::iterator it = vert2Deg.find(v);
+ // if (it == vert2Deg.end()) {
+ // vert2Deg.insert(std::make_pair(v,1));
+ // }
+ // else{
+ // int nbE = it->second+1;
+ // it->second = nbE;
+ // }
+ // }
+ // }
+ // }
+ // int dMin = 10;
+ // int dMax = 0;
+ // int d4 = 0;
+ // int nbElems = vert2Deg.size();
+ // std::map<MVertex*, int >::const_iterator itmap = vert2Deg.begin();
+ // for(; itmap !=vert2Deg.end(); itmap++){
+ // MVertex *v = itmap->first;
+ // int nbE = itmap->second;
+ // dMin = std::min(nbE, dMin);
+ // dMax = std::max(nbE, dMax);
+ // if (nbE == 4) d4 += 1;
+ // }
+ // if (nbElems > 0)
+ // printf("Stats degree vertices: dMin=%d , dMax=%d, d4=%g \n",
+ // dMin, dMax, (double)d4/nbElems);
+ // FieldManager *fields = GModel::current()->getFields();
+ // Field *f = fields->get(fields->background_field);
+ // int nbEdges = edges.size();
+ // printf("nb edges =%d \n", nbEdges);
+ // if(system("rm qualEdges.txt"));
+ // FILE *fp = fopen("qualEdges.txt", "w");
+ // std::vector<int> qualE;
+ // int nbS = 50;
+ // qualE.resize(nbS);
+ // if(fields.getBackgroundField() > 0){
+ // std::set<MEdge, Less_Edge>::iterator it = edges.begin();
+ // double sum = 0;
+ // for (; it !=edges.end();++it){
+ // MVertex *v0 = it->getVertex(0);
+ // MVertex *v1 = it->getVertex(1);
+ // double l = 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()));
+ // double lf = (*f)(0.5*(v0->x()+v1->x()), 0.5*(v0->y()+v1->y()),
+ // 0.5*(v0->z()+v1->z()),v0->onWhat());
+ // double el = l/lf;
+ // int index = (int) ceil(el*nbS*0.5);
+ // qualE[index]+= 1;
+ // double e = (l>lf) ? lf/l : l/lf;
+ // sum += e - 1.0;
+ // }
+ // double tau = exp ((1./edges.size()) * sum);
+ // printf("N edges = %d tau = %g\n",(int)edges.size(),tau);
+
+ // double ibegin = 2./(2*nbS);
+ // double inext = 2./nbS;
+ // for (int i= 0; i< qualE.size(); i++){
+ // fprintf(fp, "0 0 0 0 %g 0 0 %g \n", ibegin+i*inext , (double)qualE[i]/nbEdges);
+ // }
+
+ // }
+ // fclose(fp);
+ // }
+ // {
+ // std::vector<GEntity*> entities;
+ // std::set<MEdge, Less_Edge> edges;
+ // GModel::current()->getEntities(entities);
+ // std::map<MVertex*, int > vert2Deg;
+ // for(unsigned int i = 0; i < entities.size(); i++){
+ // if(entities[i]->dim() < 2 ) continue;
+ // if(entities[i]->tag() != 10) continue;
+ // for(unsigned int j = 0; j < entities[i]->getNumMeshElements(); j++){
+ // MElement *e = entities[i]->getMeshElement(j);
+ // for(unsigned int k = 0; k < e->getNumEdges(); k++){
+ // edges.insert(e->getEdge(k));
+ // }
+ // for(unsigned int k = 0; k < e->getNumVertices(); k++){
+ // MVertex *v = e->getVertex(k);
+ // if (v->onWhat()->dim() < 2) continue;
+ // std::map<MVertex*, int >::iterator it = vert2Deg.find(v);
+ // if (it == vert2Deg.end()){
+ // vert2Deg.insert(std::make_pair(v,1));
+ // }
+ // else{
+ // int nbE = it->second+1;
+ // it->second = nbE;
+ // }
+ // }
+ // }
+ // }
+ // int dMin = 10;
+ // int dMax = 0;
+ // int d4 = 0;
+ // int nbElems = vert2Deg.size();
+ // std::map<MVertex*, int >::const_iterator itmap = vert2Deg.begin();
+ // for(; itmap !=vert2Deg.end(); itmap++){
+ // MVertex *v = itmap->first;
+ // int nbE = itmap->second;
+ // dMin = std::min(nbE, dMin);
+ // dMax = std::max(nbE, dMax);
+ // if (nbE == 4) d4 += 1;
+ // }
+ // if (nbElems > 0) printf("Stats degree vertices: dMin=%d , dMax=%d, d4=%g \n",
+ // dMin, dMax, (double)d4/nbElems);
+ // FieldManager *fields = GModel::current()->getFields();
+ // Field *f = fields->get(fields.getBackgroundField());
+ // int nbEdges = edges.size();
+ // if(system("rm qualEdges.txt"));
+ // FILE *fp = fopen("qualEdges.txt", "w");
+ // std::vector<int> qualE;
+ // int nbS = 50;
+ // qualE.resize(nbS);
+ // if(fields.getBackgroundField() > 0){
+ // std::set<MEdge, Less_Edge>::iterator it = edges.begin();
+ // double sum = 0;
+ // for (; it !=edges.end();++it){
+ // MVertex *v0 = it->getVertex(0);
+ // MVertex *v1 = it->getVertex(1);
+ // double l = 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()));
+ // double lf = (*f)(0.5*(v0->x()+v1->x()), 0.5*(v0->y()+v1->y()),
+ // 0.5*(v0->z()+v1->z()),v0->onWhat());
+ // double el = l/lf;
+ // int index = (int) ceil(el*nbS*0.5);
+ // qualE[index]+= 1;
+ // double e = (l>lf) ? lf/l : l/lf;
+ // sum += e - 1.0;
+ // }
+ // double tau = exp ((1./edges.size()) * sum);
+ // double ibegin = 2./(2*nbS);
+ // double inext = 2./nbS;
+ // for (int i= 0; i< qualE.size(); i++){
+ // fprintf(fp, "0 0 0 0 %g 0 0 %g \n", ibegin+i*inext , (double)qualE[i]/nbEdges);
+ // }
+ // }
+ // fclose(fp);
+ // }
#endif
int num = 0;
diff --git a/Geo/GFaceCompound.cpp b/Geo/GFaceCompound.cpp
index b22b554f8ad00aa462f57fbd79403939d11ccde8..bd1b92b082bf304b55be2446d1f8a53d93ddaf4f 100644
--- a/Geo/GFaceCompound.cpp
+++ b/Geo/GFaceCompound.cpp
@@ -175,13 +175,18 @@ static bool computeCGKernelPolygon(std::map<MVertex*,SPoint3> &coordinates,
int nbPts = cavV.size();
fullMatrix<double> u(100,2);
- int i = 0;
+ int ipt = 0;
for(std::vector<MVertex*>::iterator it = cavV.begin(); it != cavV.end(); it++){
SPoint3 vsp = coordinates[*it];
- u(i,0) = vsp[0];
- u(i,1) = vsp[1];
- i++;
+ u(ipt,0) = vsp[0];
+ u(ipt,1) = vsp[1];
+ ucg += u(ipt,0);
+ vcg += u(ipt,1);
+ ipt++;
}
+ ucg /= ipt;
+ vcg /= ipt;
+
double eps = -5.e-7;
int N = nbPts;
@@ -238,6 +243,8 @@ static bool computeCGKernelPolygon(std::map<MVertex*,SPoint3> &coordinates,
int nbFinal = setP.size();
if(nbFinal > 0){
+ ucg = 0.0;
+ vcg = 0.0;
for(std::set<int>::iterator it =setP.begin(); it != setP.end(); it++){
ucg += u(*it,0);
vcg += u(*it,1);
@@ -679,11 +686,14 @@ bool GFaceCompound::checkOrientation(int iter, bool moveBoundaries) const
else if (!moveBoundaries){
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);
+ oriented = one2OneMap();
+ printStuff(iter);
+ iter++;
+ if (!oriented) return checkOrientation(iter);
}
}
- else if (iter > 0 && iter < iterMax){
+
+ if (iter > 0 && iter < iterMax){
Msg::Info("--- Flipping : no more flips (%d iter)", iter);
}
@@ -857,11 +867,11 @@ bool GFaceCompound::one2OneMap() const
std::vector<MVertex*> cavV;
myPolygon(vTri, cavV);
bool success = computeCGKernelPolygon(coordinates, cavV, u_cg, v_cg);
- if (success){
+ //if (success){ //if not succes compute with CG polygon
nbRepair++;
SPoint3 p_cg(u_cg,v_cg,0.0);
coordinates[v] = p_cg;
- }
+ //}
}
}
if (nbRepair == 0) return false;
@@ -918,24 +928,24 @@ bool GFaceCompound::parametrize() const
// Conformal map parametrization
else if (_mapping == CONFORMAL){
std::vector<MVertex *> vert;
- bool oriented;
+ bool oriented, overlap;
if (_type == SPECTRAL){
Msg::Info("Parametrizing surface %d with 'spectral conformal map'", tag());
- parametrize_conformal_spectral();
+ overlap = parametrize_conformal_spectral();
}
- else if (_type == FE){
+ else {
Msg::Info("Parametrizing surface %d with 'FE conformal map'", tag());
- parametrize_conformal(0, NULL, NULL);
+ overlap = parametrize_conformal(0, NULL, NULL);
}
- printStuff(55);
- oriented = checkOrientation(0, true);
- printStuff(77);
- if (_type==SPECTRAL && (!oriented || checkOverlap(vert)) ){
+ //printStuff(55);
+ oriented = checkOrientation(0);
+ //printStuff(77);
+ if (_type==SPECTRAL && (!oriented || overlap) ){
Msg::Warning("!!! parametrization switched to 'FE conformal' map");
- parametrize_conformal(0, NULL, NULL);
- oriented = checkOrientation(0, true);
+ overlap = parametrize_conformal(0, NULL, NULL);
+ oriented = checkOrientation(0);
}
- if (!oriented || checkOverlap(vert)){
+ if (!oriented || overlap){
Msg::Warning("$$$ parametrization switched to 'convex' map");
_type = UNITCIRCLE;
parametrize(ITERU,CONVEX);
@@ -973,6 +983,7 @@ bool GFaceCompound::parametrize() const
if (_mapping != RBF){
if (!checkOrientation(0)){
+ printStuff(22);
Msg::Info("### parametrization switched to 'convex map' onto circle");
printStuff(33);
_type = UNITCIRCLE;
diff --git a/Geo/MHexahedron.cpp b/Geo/MHexahedron.cpp
index 615a23d2091af589489823d8cc117d243b794eb0..070e21af38ddcd60ba398ac7b556dad4aa201885 100644
--- a/Geo/MHexahedron.cpp
+++ b/Geo/MHexahedron.cpp
@@ -8,6 +8,7 @@
#include "Context.h"
#include "polynomialBasis.h"
#include "MQuadrangle.h"
+#include "qualityMeasures.h"
int MHexahedron::getVolumeSign()
{
double mat[3][3];
@@ -32,6 +33,39 @@ void MHexahedron::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
*pts = getGQHPts(pOrder);
}
+double MHexahedron::angleShapeMeasure()
+{
+
+#if defined(HAVE_MESH)
+ double angleMax = 0.0;
+ double angleMin = M_PI;
+ double zeta = 0.0;
+ for (int i=0; i<getNumFaces(); i++){
+ std::vector<MVertex*> vv;
+ vv.push_back(getFace(i).getVertex(0));
+ vv.push_back(getFace(i).getVertex(1));
+ vv.push_back(getFace(i).getVertex(2));
+ vv.push_back(getFace(i).getVertex(3));
+ // MVertex *v0 = new MVertex(0, 0, 0); vv.push_back(v0);
+ // MVertex *v1 = new MVertex(1., 0, 0);vv.push_back(v1);
+ // MVertex *v2 = new MVertex(2., 1., 0);vv.push_back(v2);
+ // MVertex *v3 = new MVertex(1, 1., 0);vv.push_back(v3);
+ for (int j=0; j<4; j++){
+ SVector3 a(vv[(j+2)%4]->x()-vv[(j+1)%4]->x(),vv[(j+2)%4]->y()-vv[(j+1)%4]->y(),vv[(j+2)%4]->z()-vv[(j+1)%4]->z() );
+ SVector3 b(vv[(j+1)%4]->x()-vv[(j)%4]->x(), vv[(j+1)%4]->y()-vv[(j)%4]->y(), vv[(j+1)%4]->z()-vv[(j)%4]->z() );
+ double angle = acos( dot(a,b)/(norm(a)*norm(b))); //*180/M_PI;
+ angleMax = std::max(angleMax, angle);
+ angleMin = std::min(angleMin, angle);
+ }
+ //printf("angle max =%g min =%g \n", angleMax*180/M_PI, angleMin*180/M_PI);
+ }
+ zeta = 1.-std::max((angleMax-0.5*M_PI)/(0.5*M_PI),(0.5*M_PI-angleMin)/(0.5*M_PI));
+ return zeta;
+#else
+ return 1.;
+#endif
+
+}
double MHexahedron::getInnerRadius()
{
//Only for vertically aligned elements (not inclined)
diff --git a/Geo/MHexahedron.h b/Geo/MHexahedron.h
index c50db1af48e4b0837405436d52a7cdd82be06eff..4cadf664f2197904838c92b82453a9291a436012 100644
--- a/Geo/MHexahedron.h
+++ b/Geo/MHexahedron.h
@@ -89,6 +89,7 @@ class MHexahedron : public MElement {
_v[faces_hexa(num, 3)]);
}
virtual double getInnerRadius();
+ virtual double angleShapeMeasure();
virtual void getFaceInfo (const MFace & face, int &ithFace, int &sign, int &rot)const;
virtual int getNumFacesRep(){ return 12; }
virtual void getFaceRep(int num, double *x, double *y, double *z, SVector3 *n)
diff --git a/Mesh/CenterlineField.cpp b/Mesh/CenterlineField.cpp
index 519667629a6e71cf828325dbfeb039c15ee3ad8e..d6bc9feeabab966a440abc4248800bfb453e4e1f 100644
--- a/Mesh/CenterlineField.cpp
+++ b/Mesh/CenterlineField.cpp
@@ -869,7 +869,7 @@ void Centerline::cutMesh()
//for (int k= 0; k< edges[i].children.size() ; k++) printf("%d ", edges[i].children[k].tag);
//printf("\n");
- int nbSplit = (int)floor(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;
@@ -1135,7 +1135,7 @@ void Centerline::operator() (double x, double y, double z, SMetric3 &metr, GEnt
curvature.vertexNodalValues(vertices[index[0]], curv, 0);
double sign = (curv > 0.0) ? -1.0: 1.0;
- double beta = CTX::instance()->mesh.smoothRatio;
+ double beta = CTX::instance()->mesh.smoothRatio; //beta = 1.25 better !
double ratio = 1.1;
double thickness = radMax/3.;
diff --git a/Mesh/meshGFaceElliptic.cpp b/Mesh/meshGFaceElliptic.cpp
index 43b8d5b4fcafbffea12aa33a276a10f104f29f64..6ff27d4f17d2dae212f069d2fae8e43942d8f6ca 100644
--- a/Mesh/meshGFaceElliptic.cpp
+++ b/Mesh/meshGFaceElliptic.cpp
@@ -41,7 +41,7 @@ static void printQuads(GFace *gf, fullMatrix<SPoint2> uv, std::vector<MQuadrangl
FILE *f = fopen(name,"w");
fprintf(f,"View \"%s\" {\n",name);
- for (int i = 0; i < uv.size1(); i++)
+ for (int i = 1; i < uv.size1()-1; i++)
for (int j = 0; j < uv.size2(); j++)
fprintf(f,"SP(%g,%g,%g) {%d};\n", uv(i,j).x(), uv(i,j).y(), 0.0, i);
@@ -85,11 +85,17 @@ static void printParamGrid(GFace *gf, std::vector<MVertex*> vert1, std::vector<M
FILE *f = fopen(name,"w");
fprintf(f,"View \"%s\" {\n",name);
- for (unsigned int i = 0; i < p1.size(); i++)
- fprintf(f,"SP(%g,%g,%g) {%d};\n", p1[i].x(), p1[i].y(), 0.0, i);
+ // for (unsigned int i = 0; i < p1.size(); i++)
+ // fprintf(f,"SP(%g,%g,%g) {%d};\n", p1[i].x(), p1[i].y(), 0.0, i);
+ // for (unsigned int j = 0; j < p2.size(); j++)
+ // fprintf(f,"SP(%g,%g,%g) {%d};\n", p2[j].x(), p2[j].y(), 0.0, 100+j);
- for (unsigned int j = 0; j < p2.size(); j++)
- fprintf(f,"SP(%g,%g,%g) {%d};\n", p2[j].x(), p2[j].y(), 0.0, 100+j);
+ for (unsigned int i = 0; i < p1.size()-1; i++)
+ fprintf(f,"SL(%g,%g,%g,%g,%g,%g) {%d,%d};\n", p1[i].x(), p1[i].y(), 0.0, p1[i+1].x(), p1[i+1].y(), 0.0, 1, 1);
+ for (unsigned int i = 0; i < p2.size()-1; i++)
+ fprintf(f,"SL(%g,%g,%g,%g,%g,%g) {%d,%d};\n", p2[i].x(), p2[i].y(), 0.0, p2[i+1].x(), p2[i+1].y(), 0.0, 1, 1);
+ fprintf(f,"SL(%g,%g,%g,%g,%g,%g) {%d,%d};\n", p1[p1.size()-1].x(), p1[ p1.size()-1].y(), 0.0, p1[0].x(), p1[0].y(), 0.0, 1, 1);
+ fprintf(f,"SL(%g,%g,%g,%g,%g,%g) {%d,%d};\n", p2[p2.size()-1].x(), p2[ p2.size()-1].y(), 0.0, p2[0].x(), p2[0].y(), 0.0, 1, 1);
fprintf(f,"};\n");
fclose(f);
@@ -235,7 +241,7 @@ static void transfiniteSmoother(GFace* gf,
int N = uv.size2();
int jStart = isPeriodic ? 0 : 1;
- int numSmooth = 100;
+ int numSmooth = 150;
fullMatrix<SPoint2> uvold = uv;
for(int k = 0; k < numSmooth; k++){
double norm = 0.0;
@@ -750,7 +756,6 @@ bool createRegularTwoCircleGridPeriodic (Centerline *center, GFace *gf)
updateFaceQuads(gf, quads, newv);
//exit(1);
-
//printParamGrid(gf, vert1, vert2, e00,e22,e02,e02,e02,e02, quads);
return true;
diff --git a/Mesh/qualityMeasures.cpp b/Mesh/qualityMeasures.cpp
index 425fe59d194110f04065d9b90fac9ed0c7f053a1..ca48affd43236c4f973ef25f4cd4c89edc79c1b6 100644
--- a/Mesh/qualityMeasures.cpp
+++ b/Mesh/qualityMeasures.cpp
@@ -429,7 +429,7 @@ double qmTriangleAngles (MTriangle *e) {
rot[2][0]= 0; rot[2][1]=0; rot[2][2]=1;
double tmp[3][3];
- // double minAngle = 120.0;
+ //double minAngle = 120.0;
for (int i = 0; i < e->getNumPrimaryVertices(); i++) {
const double u = i == 1 ? 1 : 0;
const double v = i == 2 ? 1 : 0;
@@ -468,13 +468,12 @@ double qmTriangleAngles (MTriangle *e) {
double quality = (atan(a*(x+M_PI/9)) + atan(a*(M_PI/9-x)))/den;
worst_quality = std::min(worst_quality, quality);
- // minAngle = std::min(angle, minAngle);
- // printf("Angle %g ", angle);
+ //minAngle = std::min(angle, minAngle);
+ //printf("Angle %g ", angle);
// printf("Quality %g\n",quality);
}
- // printf("MinAngle %g ", minAngle);
- // printf("\n");
- // return minAngle;
+ //printf("MinAngle %g \n", minAngle);
+ //return minAngle;
return worst_quality;
}
@@ -530,12 +529,13 @@ double qmQuadrangleAngles (MQuadrangle *e) {
double c;
prosca(v1,v2,&c);
- // printf("Youhou %g %g\n",c,acos(c)*180/M_PI);
double x = fabs(acos(c))-M_PI/2;
+ //double angle = fabs(acos(c))*180/M_PI;
double quality = (atan(a*(x+M_PI/4)) + atan(a*(2*M_PI/4 - (x+M_PI/4))))/den;
worst_quality = std::min(worst_quality, quality);
}
-
+
return worst_quality;
+
}
diff --git a/benchmarks/centerlines/aneurysm_centerlines.geo b/benchmarks/centerlines/aneurysm_centerlines.geo
index 02831d26ce8272000a13f7777efe56a23af1cda6..ba4987ebd668e066d7c5b549ffee3186ca71343e 100644
--- a/benchmarks/centerlines/aneurysm_centerlines.geo
+++ b/benchmarks/centerlines/aneurysm_centerlines.geo
@@ -4,7 +4,7 @@ Mesh.Algorithm3D = 7; //(1=tetgen, 4=netgen, 7=mmg3D
Mesh.LcIntegrationPrecision = 1.e-2;
Mesh.RecombineAll = 1;
-Mesh.Bunin = 170;
+//Mesh.Bunin = 120;
Merge "aneu_ext.stl";
diff --git a/benchmarks/centerlines/aorta_centerlines.geo b/benchmarks/centerlines/aorta_centerlines.geo
index 52c55b59f3cacbc843bb985a39f64e4d7c2a6df8..a97090ea88dd0265631478c192a2294fffdd7b2f 100644
--- a/benchmarks/centerlines/aorta_centerlines.geo
+++ b/benchmarks/centerlines/aorta_centerlines.geo
@@ -1,24 +1,29 @@
Mesh.Algorithm = 8; //(1=MeshAdapt, 2=Automatic, 5=Delaunay, 6=Frontal, 7=bamg, 8=delquad)
-Mesh.Algorithm3D = 7; //(1=tetgen, 4=netgen, 5=FrontalDel, 6=FrontalHex, 7=MMG3D, 9=R-tree
+Mesh.Algorithm3D = 8; //(1=tetgen, 4=netgen, 5=FrontalDel, 6=FrontalHex, 7=MMG3D, 9=R-tree
Mesh.LcIntegrationPrecision = 1.e-3;
Mesh.RecombineAll = 1;
-Mesh.Bunin = 150;
+Mesh.Bunin = 60;
Merge "aorta2.stl";
Field[1] = Centerline;
Field[1].FileName = "centerlinesAORTA.vtk";
-Field[1].nbPoints = 20; //33;
+Field[1].nbPoints = 25; //33;
Field[1].nbElemLayer = 4;
Field[1].hLayer = 0.2;//percent of vessel radius
Field[1].closeVolume =1;
-//Field[1].extrudeWall =1;
+Field[1].extrudeWall =1;
Field[1].reMesh =1;
Field[1].run;
+
Background Field = 1;
+
+
+
+
diff --git a/benchmarks/centerlines/bypass_centerlines.geo b/benchmarks/centerlines/bypass_centerlines.geo
index a21967ad288997a459a81fdac6dd44fb233576de..916b8a5aaa5ecd0779f1e5e66309b0fa3ae287f8 100644
--- a/benchmarks/centerlines/bypass_centerlines.geo
+++ b/benchmarks/centerlines/bypass_centerlines.geo
@@ -1,4 +1,4 @@
-Mesh.Algorithm = 8; //(1=MeshAdapt, 2=Automatic, 5=Delaunay, 6=Frontal, 7=bamg, 8=delquad)
+Mesh.Algorithm = 7; //(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;
@@ -10,7 +10,7 @@ Merge "bypass.stl";
Field[1] = Centerline;
Field[1].FileName = "centerlinesBYPASS.msh";
-Field[1].nbPoints = 25;
+Field[1].nbPoints = 15;
Field[1].nbElemLayer = 4;
Field[1].hLayer = 0.2;//percent of vessel radius
diff --git a/benchmarks/centerlines/carotid_centerlines.geo b/benchmarks/centerlines/carotid_centerlines.geo
index 3f87cc1f7ef8dedb272c47a3a89b3b28b10d1a0d..d0a2d072441ef07dcab9fa617cecc23af324e1a2 100644
--- a/benchmarks/centerlines/carotid_centerlines.geo
+++ b/benchmarks/centerlines/carotid_centerlines.geo
@@ -4,7 +4,7 @@ Mesh.Algorithm3D = 7;
Mesh.LcIntegrationPrecision = 1.e-2;
Mesh.RecombineAll = 1;
-Mesh.Bunin = 130;
+Mesh.Bunin = 100;
Merge "carotid.stl";
@@ -16,7 +16,7 @@ Field[1].nbElemLayer = 4;
Field[1].hLayer = 0.2; //percent of vessel radius
Field[1].closeVolume =1;
-//Field[1].extrudeWall =1;
+Field[1].extrudeWall =1;
Field[1].reMesh =1;
Field[1].run;
diff --git a/benchmarks/centerlines/lung_centerlines.geo b/benchmarks/centerlines/lung_centerlines.geo
index 59ddbc2844a779b6c343eb29027e9b6f0d07c1ba..c1fc3361159b21081b1f54886dd229c73ad7686a 100644
--- a/benchmarks/centerlines/lung_centerlines.geo
+++ b/benchmarks/centerlines/lung_centerlines.geo
@@ -1,7 +1,8 @@
-Mesh.Algorithm = 5; //(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 = 1; //(1=tetgen, 4=netgen, 7=mmg3D
-//Mesh.RecombineAll = 1;
+Mesh.RecombineAll = 1;
+Mesh.Bunin = 100;
Mesh.LcIntegrationPrecision = 1.e-2;
Merge "lung.stl";