From 0ef1953f2cd7a4744fe415a456ed8fd5d1debace Mon Sep 17 00:00:00 2001
From: Christophe Geuzaine <cgeuzaine@ulg.ac.be>
Date: Sat, 16 Jun 2012 14:17:30 +0000
Subject: [PATCH] cleanup
---
Geo/MElementCut.cpp | 32 ++-
Geo/MVertexBoundaryLayerData.cpp | 16 +-
Mesh/BDS.cpp | 150 ++++++------
Mesh/CMakeLists.txt | 6 +-
Mesh/CenterlineField.cpp | 325 +++++++++++++-------------
Mesh/HighOrder.cpp | 173 +++++++-------
Mesh/highOrderTools.cpp | 71 +++---
Mesh/meshGFaceBoundaryLayers.cpp | 14 +-
Mesh/meshGRegionDelaunayInsertion.cpp | 258 ++++++++++----------
Mesh/meshMetric.cpp | 95 ++++----
Mesh/qualityMeasures.cpp | 6 +-
11 files changed, 547 insertions(+), 599 deletions(-)
diff --git a/Geo/MElementCut.cpp b/Geo/MElementCut.cpp
index b9ba13a212..865054fd6a 100644
--- a/Geo/MElementCut.cpp
+++ b/Geo/MElementCut.cpp
@@ -186,7 +186,7 @@ void MPolygon::_initVertices()
}
}
if(!found){
- for(k = 0; k < multiEdges.size(); k++)
+ for(k = 0; k < (int)multiEdges.size(); k++)
if(multiEdges[k].isInside(ed.getVertex(0)) &&
multiEdges[k].isInside(ed.getVertex(1))){
found = true; break;
@@ -343,7 +343,6 @@ void MLineChild::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
*npts = 0;
if(_intpt) delete [] _intpt;
if(!_orig) return;
- double jac[3][3];
_intpt = new IntPt[getNGQLPts(pOrder)];
int nptsi;
IntPt *ptsi;
@@ -469,7 +468,6 @@ void MLineBorder::getIntegrationPoints(int pOrder, int *npts, IntPt **pts)
MVertex v1(uvw[1][0], uvw[1][1], uvw[1][2]);
MLine ll(&v0, &v1);
ll.getIntegrationPoints(pOrder, &nptsi, &ptsi);
- double jac[3][3];
for(int ip = 0; ip < nptsi; ip++){
const double u = ptsi[ip].pt[0];
const double v = ptsi[ip].pt[1];
@@ -595,7 +593,7 @@ static void elementSplitMesh(MElement *e, fullMatrix<double> &verticesLs,
std::map<int, int> newElemTags[4],
std::map<int, int> newPhysTags[4],
std::map<int, int> borderElemTags[2],
- std::map<int, int> borderPhysTags[2],
+ std::map<int, int> borderPhysTags[2],
int gLsTag)
{
int elementary = ge->tag();
@@ -623,7 +621,7 @@ static void elementSplitMesh(MElement *e, fullMatrix<double> &verticesLs,
lsTag = -1;
splitElem = true;
break;
- }
+ }
}
// int lsTag = 1; //negative ls
// for(int k = 0; k < e->getNumVertices(); k++){
@@ -744,7 +742,7 @@ static void elementSplitMesh(MElement *e, fullMatrix<double> &verticesLs,
getBorderTag(gLsTag, c, newPhysTags[2][0], borderPhysTags[1]);
if(mf.getNumVertices() == 3){
MFace f1 = MFace(vertexMap[mf.getVertex(0)->getNum()],
- vertexMap[mf.getVertex(1)->getNum()],
+ vertexMap[mf.getVertex(1)->getNum()],
vertexMap[mf.getVertex(2)->getNum()]);
int i = elements[2][reg].size() - 1;
for(; i >= 0; i--){
@@ -753,7 +751,7 @@ static void elementSplitMesh(MElement *e, fullMatrix<double> &verticesLs,
if(f1 == f2) break;
}
if(i < 0){
- MTriangle *tri = new MTriangle(f1.getVertex(0), f1.getVertex(1), f1.getVertex(2));
+ MTriangle *tri = new MTriangle(f1.getVertex(0), f1.getVertex(1), f1.getVertex(2));
elements[2][reg].push_back(tri);
}
else{
@@ -764,7 +762,7 @@ static void elementSplitMesh(MElement *e, fullMatrix<double> &verticesLs,
}
else if(mf.getNumVertices() == 4){
MFace f1 = MFace(vertexMap[mf.getVertex(0)->getNum()],
- vertexMap[mf.getVertex(1)->getNum()],
+ vertexMap[mf.getVertex(1)->getNum()],
vertexMap[mf.getVertex(2)->getNum()],
vertexMap[mf.getVertex(3)->getNum()]);
int i;
@@ -786,7 +784,7 @@ static void elementSplitMesh(MElement *e, fullMatrix<double> &verticesLs,
}
}
if(physTag) assignLsPhysical(GM, reg, 2, physicals, physTag, gLsTag);
- }
+ }
}
}
}
@@ -1480,7 +1478,7 @@ GModel *buildCutMesh(GModel *gm, gLevelset *ls,
MElement *e = gmEntities[i]->getMeshElement(j);
e->setVolumePositive();
elementSplitMesh(e, verticesLs, gmEntities[i], gm, numEle, vertexMap, newParents,
- newDomains, elements, physicals, newElemTags, newPhysTags,
+ newDomains, elements, physicals, newElemTags, newPhysTags,
borderElemTags, borderPhysTags, ls->getTag());
cutGM->getMeshPartitions().insert(e->getPartition());
}
@@ -1516,8 +1514,8 @@ GModel *buildCutMesh(GModel *gm, gLevelset *ls,
}
// Create elementary and physical for non connected border lines
- if(borders[0].size() > nbBorders && gmEntities[i]->dim() == 2 &&
- i == gm->getNumVertices() + gm->getNumEdges() + gm->getNumFaces() - 1){
+ if((int)borders[0].size() > nbBorders && gmEntities[i]->dim() == 2 &&
+ (int)i == gm->getNumVertices() + gm->getNumEdges() + gm->getNumFaces() - 1){
int k = 0;
for (std::map<int, std::vector<MElement*> >::iterator it = elements[1].begin();
it != elements[1].end(); it++){
@@ -1530,12 +1528,12 @@ GModel *buildCutMesh(GModel *gm, gLevelset *ls,
int nLR = lsLineRegs[j];
bool havePhys = physicals[1][nLR].size();
while(1){
- std::vector<MElement*> conLines;
+ std::vector<MElement*> conLines;
conLines.push_back(elements[1][nLR][0]);
elements[1][nLR].erase(elements[1][nLR].begin());
MVertex *v1 = conLines[0]->getVertex(0);
MVertex *v2 = conLines[0]->getVertex(1);
- for(int k = 0; k < elements[1][nLR].size(); ){
+ for(unsigned int k = 0; k < elements[1][nLR].size(); ){
MVertex *va = elements[1][nLR][k]->getVertex(0);
MVertex *vb = elements[1][nLR][k]->getVertex(1);
if(va == v1 || vb == v1 || va == v2 || vb == v2){
@@ -1556,11 +1554,11 @@ GModel *buildCutMesh(GModel *gm, gLevelset *ls,
if(newPhys)
assignLsPhysical(gm, newReg, 1, physicals, newPhys,
lines[lines.size() - 1]->lsTag());
- for(int k = 0; k < conLines.size(); k++)
+ for(unsigned int k = 0; k < conLines.size(); k++)
elements[1][newReg].push_back(conLines[k]);
}
else {
- for(int k = 0; k < conLines.size(); k++)
+ for(unsigned int k = 0; k < conLines.size(); k++)
elements[1][nLR].push_back(conLines[k]);
break;
}
@@ -1595,7 +1593,7 @@ GModel *buildCutMesh(GModel *gm, gLevelset *ls,
}
printf("PHYS\n");
for(int i=0;i<4;i++)
- for(std::map<int, std::map<int, std::string> >::iterator it=physicals[i].begin();it!=physicals[i].end();it++)
+ for(std::map<int, std::map<int, std::string> >::iterator it=physicals[i].begin();it!=physicals[i].end();it++)
for(std::map<int, std::string>::iterator it2 = it->second.begin(); it2!=it->second.end(); it2++)
printf(" dim=%d reg=%d phys=%d \"%s\"\n",i,it->first,it2->first,it2->second.c_str());
printf("\n");
diff --git a/Geo/MVertexBoundaryLayerData.cpp b/Geo/MVertexBoundaryLayerData.cpp
index d709fe1199..2623e34e37 100644
--- a/Geo/MVertexBoundaryLayerData.cpp
+++ b/Geo/MVertexBoundaryLayerData.cpp
@@ -7,28 +7,30 @@
std::vector<MVertex*>* MVertexBoundaryLayerData::getChildren(int i)
{
- if (i < this->children.size() && i >= 0) {
+ if (i < (int)this->children.size() && i >= 0) {
return &(children[i]);
- } else {
+ }
+ else {
return 0;
}
}
-int MVertexBoundaryLayerData::getNumChildren(int i)
+int MVertexBoundaryLayerData::getNumChildren(int i)
{
- if (i < this->children.size() && i >= 0) {
+ if (i < (int)this->children.size() && i >= 0) {
return this->children[i].size();
- } else {
+ }
+ else {
return -1;
}
}
-int MVertexBoundaryLayerData::getNumChildrenFamilies()
+int MVertexBoundaryLayerData::getNumChildrenFamilies()
{
return this->children.size();
}
-void MVertexBoundaryLayerData::addChildrenFamily(std::vector<MVertex*> family)
+void MVertexBoundaryLayerData::addChildrenFamily(std::vector<MVertex*> family)
{
this->children.push_back(family);
}
diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index 75826a78f7..09ae596fb3 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -36,15 +36,15 @@ void outputScalarField(std::list<BDS_Face*> t, const char *iii, int param, GFace
fprintf(f, "SQ(%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){%d,%d,%d,%d};\n",
pts[0]->X, pts[0]->Y, pts[0]->Z,
pts[1]->X, pts[1]->Y, pts[1]->Z,
- pts[2]->X, pts[2]->Y, pts[2]->Z,
- pts[3]->X, pts[3]->Y, pts[3]->Z,
+ pts[2]->X, pts[2]->Y, pts[2]->Z,
+ pts[3]->X, pts[3]->Y, pts[3]->Z,
pts[0]->iD, pts[1]->iD, pts[2]->iD, pts[3]->iD);
}
else{
fprintf(f, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%d,%d,%d};\n",
pts[0]->X, pts[0]->Y, pts[0]->Z,
pts[1]->X, pts[1]->Y, pts[1]->Z,
- pts[2]->X, pts[2]->Y, pts[2]->Z,
+ pts[2]->X, pts[2]->Y, pts[2]->Z,
pts[0]->iD, pts[1]->iD, pts[2]->iD);
}
}
@@ -56,7 +56,7 @@ void outputScalarField(std::list<BDS_Face*> t, const char *iii, int param, GFace
pts[2]->X, pts[2]->Y, pts[2]->Z,
pts[3]->X, pts[3]->Y, pts[3]->Z,
gf->curvatureDiv(SPoint2(pts[0]->u, pts[0]->v)),
- gf->curvatureDiv(SPoint2(pts[1]->u, pts[1]->v)),
+ gf->curvatureDiv(SPoint2(pts[1]->u, pts[1]->v)),
gf->curvatureDiv(SPoint2(pts[2]->u, pts[2]->v)),
gf->curvatureDiv(SPoint2(pts[3]->u, pts[3]->v)));
}
@@ -108,12 +108,14 @@ void normal_triangle(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3,
norme(c);
}
+/*
static double surface_triangle(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3)
{
double c[3];
vector_triangle(p1, p2, p3, c);
return 0.5 * sqrt(c[0] * c[0] + c[1] * c[1] + c[2] * c[2]);
}
+*/
static double surface_triangle_param(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3)
{
@@ -157,8 +159,8 @@ BDS_Point *BDS_Mesh::add_point(int num, double x, double y, double z)
}
BDS_Point *BDS_Mesh::add_point(int num, double u, double v, GFace *gf)
-{
- GPoint gp = gf->point(u*scalingU,v*scalingV);
+{
+ GPoint gp = gf->point(u*scalingU,v*scalingV);
BDS_Point *pp = new BDS_Point(num, gp.x(), gp.y(), gp.z());
pp->u = u;
pp->v = v;
@@ -209,7 +211,7 @@ int Intersect_Edges_2d(double x1, double y1, double x2, double y2,
// double p2[2] = {x2,y2};
// double q1[2] = {x3,y3};
// double q2[2] = {x4,y4};
-
+
// double rp1 = robustPredicates::orient2d(p1, p2, q1);
// double rp2 = robustPredicates::orient2d(p1, p2, q2);
// double rq1 = robustPredicates::orient2d(q1, q2, p1);
@@ -239,7 +241,7 @@ int Intersect_Edges_2d(double x1, double y1, double x2, double y2,
}
BDS_Edge *BDS_Mesh::recover_edge_fast(BDS_Point *p1, BDS_Point *p2){
-
+
std::list<BDS_Face*> ts;
p1->getTriangles(ts);
std::list<BDS_Face*>::iterator it = ts.begin();
@@ -264,7 +266,7 @@ BDS_Edge *BDS_Mesh::recover_edge_fast(BDS_Point *p1, BDS_Point *p2){
}
BDS_Edge *BDS_Mesh::recover_edge(int num1, int num2, bool &_fatal,
- std::set<EdgeToRecover> *e2r,
+ std::set<EdgeToRecover> *e2r,
std::set<EdgeToRecover> *not_recovered)
{
BDS_Edge *e = find_edge(num1, num2);
@@ -274,14 +276,14 @@ BDS_Edge *BDS_Mesh::recover_edge(int num1, int num2, bool &_fatal,
BDS_Point *p1 = find_point(num1);
BDS_Point *p2 = find_point(num2);
-
+
if(!p1 || !p2) {
Msg::Fatal("Could not find points %d or %d in BDS mesh", num1, num2);
return 0;
}
Msg::Debug("edge %d %d has to be recovered", num1, num2);
-
+
int ix = 0;
double x[2];
while(1){
@@ -299,9 +301,9 @@ BDS_Edge *BDS_Mesh::recover_edge(int num1, int num2, bool &_fatal,
p2->u, p2->v,x)){
// intersect
if(e2r && e2r->find(EdgeToRecover(e->p1->iD, e->p2->iD, 0)) != e2r->end()){
- std::set<EdgeToRecover>::iterator itr1 =
- e2r->find(EdgeToRecover(e->p1->iD, e->p2->iD, 0));
- std::set<EdgeToRecover>::iterator itr2 =
+ std::set<EdgeToRecover>::iterator itr1 =
+ e2r->find(EdgeToRecover(e->p1->iD, e->p2->iD, 0));
+ std::set<EdgeToRecover>::iterator itr2 =
e2r->find(EdgeToRecover(num1, num2, 0));
Msg::Debug("edge %d %d on model edge %d cannot be recovered because"
" it intersects %d %d on model edge %d", num1, num2, itr2->ge->tag(),
@@ -312,7 +314,7 @@ BDS_Edge *BDS_Mesh::recover_edge(int num1, int num2, bool &_fatal,
selfIntersection = true;
}
if (e->numfaces() != e->numTriangles()) return 0;
- intersected.push_back(e);
+ intersected.push_back(e);
}
++it;
}
@@ -321,11 +323,11 @@ BDS_Edge *BDS_Mesh::recover_edge(int num1, int num2, bool &_fatal,
// if(ix > 300){
// Msg::Warning("edge %d %d cannot be recovered after %d iterations, trying again",
-// num1, num2, ix);
+// num1, num2, ix);
// ix = 0;
// }
// printf("%d %d\n",intersected.size(),ix);
-
+
if(!intersected.size() || ix > 1000){
BDS_Edge *eee = find_edge(num1, num2);
if(!eee){
@@ -338,9 +340,9 @@ BDS_Edge *BDS_Mesh::recover_edge(int num1, int num2, bool &_fatal,
}
return eee;
}
-
+
int ichoice = ix++ % intersected.size();
- //bool success =
+ //bool success =
swap_edge(intersected[ichoice], BDS_SwapEdgeTestQuality(false, false));
// printf("trying to swop %d %d = %d (%d %d)\n", intersected[ichoice]->p1->iD,
// intersected[ichoice]->p2->iD, success, intersected[ichoice]->deleted,
@@ -436,7 +438,7 @@ BDS_Face *BDS_Mesh::add_triangle(BDS_Edge *e1, BDS_Edge *e2, BDS_Edge *e3)
{
// BDS_Face *tfound = find_triangle(e1, e2, e3);
// if(tfound) return tfound;
-
+
BDS_Face *t = new BDS_Face(e1, e2, e3);
triangles.push_back(t);
return t;
@@ -551,17 +553,17 @@ void BDS_Mesh::cleanup()
it++;
}
}
- {
+ {
std::list<BDS_Edge*>::iterator it = edges.begin();
while(it != edges.end()){
if((*it)->deleted){
delete *it;
it = edges.erase(it);
- }
+ }
else
it++;
- }
- }
+ }
+ }
}
BDS_Mesh::~BDS_Mesh()
@@ -575,9 +577,9 @@ BDS_Mesh::~BDS_Mesh()
bool BDS_Mesh::split_face(BDS_Face *f, BDS_Point *mid)
{
- BDS_Point *p1 = f->e1->commonvertex(f->e2);
- BDS_Point *p2 = f->e3->commonvertex(f->e2);
- BDS_Point *p3 = f->e1->commonvertex(f->e3);
+ BDS_Point *p1 = f->e1->commonvertex(f->e2);
+ BDS_Point *p2 = f->e3->commonvertex(f->e2);
+ BDS_Point *p3 = f->e1->commonvertex(f->e3);
BDS_Edge *p1_mid = new BDS_Edge(p1, mid);
edges.push_back(p1_mid);
BDS_Edge *p2_mid = new BDS_Edge(p2, mid);
@@ -636,7 +638,7 @@ bool BDS_Mesh::split_edge(BDS_Edge *e, BDS_Point *mid)
BDS_Point *p2 = e->p2;
e->oppositeof(op);
-
+
BDS_Point *pts1[4];
e->faces(0)->getNodes(pts1);
@@ -651,7 +653,7 @@ bool BDS_Mesh::split_edge(BDS_Edge *e, BDS_Point *mid)
}
}
- // we should project
+ // we should project
BDS_GeomEntity *g1 = 0, *g2 = 0, *ge = e->g;
BDS_Edge *p1_op1 = find_edge(p1, op[0], e->faces(0));
@@ -726,24 +728,24 @@ bool BDS_Mesh::split_edge(BDS_Edge *e, BDS_Point *mid)
bool BDS_SwapEdgeTestQuality::operator () (BDS_Point *_p1, BDS_Point *_p2,
BDS_Point *_q1, BDS_Point *_q2) const
-{
+{
if(!testSmallTriangles){
double p1 [2] = {_p1->u,_p1->v};
double p2 [2] = {_p2->u,_p2->v};
double op1[2] = {_q1->u,_q1->v};
double op2[2] = {_q2->u,_q2->v};
-
+
double ori_t1 = robustPredicates::orient2d(op1, p1, op2);
double ori_t2 = robustPredicates::orient2d(op1, op2, p2);
-
+
// printf("%d %d %d %d %g %g\n",_p1->iD,_p2->iD,_q1->iD,_q2->iD,ori_t1,ori_t2);
return (ori_t1 * ori_t2 > 0); // the quadrangle was strictly convex !
}
-
- double s1 = fabs(surface_triangle_param(_p1, _p2, _q1));
- double s2 = fabs(surface_triangle_param(_p1, _p2, _q2));
- double s3 = fabs(surface_triangle_param(_p1, _q1, _q2));
- double s4 = fabs(surface_triangle_param(_p2, _q1, _q2));
+
+ double s1 = fabs(surface_triangle_param(_p1, _p2, _q1));
+ double s2 = fabs(surface_triangle_param(_p1, _p2, _q2));
+ double s3 = fabs(surface_triangle_param(_p1, _q1, _q2));
+ double s4 = fabs(surface_triangle_param(_p2, _q1, _q2));
if(fabs(s1 + s2 - s3 - s4) > 1.e-12 * (s1 + s2)) return false;
if(s3 < .02 * (s1 + s2) || s4 < .02 * (s1 + s2))
return false;
@@ -757,9 +759,9 @@ bool BDS_SwapEdgeTestQuality::operator () (BDS_Point *_p1, BDS_Point *_p2, BDS_P
{
if(!testQuality) return true;
double n[3], q[3], on[3], oq[3];
- normal_triangle(_p1, _p2, _p3, n);
- normal_triangle(_q1, _q2, _q3, q);
- normal_triangle(_op1, _op2, _op3, on);
+ normal_triangle(_p1, _p2, _p3, n);
+ normal_triangle(_q1, _q2, _q3, q);
+ normal_triangle(_op1, _op2, _op3, on);
normal_triangle(_oq1, _oq2, _oq3, oq);
double cosnq; prosca(n, q, &cosnq);
@@ -770,21 +772,21 @@ bool BDS_SwapEdgeTestQuality::operator () (BDS_Point *_p1, BDS_Point *_p2, BDS_P
double qb1 = qmTriangle(_op1, _op2, _op3, QMTRI_RHO);
double qb2 = qmTriangle(_oq1, _oq2, _oq3, QMTRI_RHO);
- // we swap for a better configuration
+ // we swap for a better configuration
double mina = std::min(qa1,qa2);
double minb = std::min(qb1,qb2);
// if(cosnq < .3 && cosonq > .5 && minb > .1)
// printf("mina = %g minb = %g cos %g %g\n",mina,minb,cosnq,cosonq);
-
+
if(cosnq < .3 && cosonq > .5 && minb > .1) return true;
-
+
if(minb > mina) return true;
// if(mina > minb && cosnq <= cosonq)return true;
return false;
}
-void swap_config(BDS_Edge *e,
+void swap_config(BDS_Edge *e,
BDS_Point **p11, BDS_Point **p12, BDS_Point **p13,
BDS_Point **p21, BDS_Point **p22, BDS_Point **p23,
BDS_Point **p31, BDS_Point **p32, BDS_Point **p33,
@@ -810,7 +812,7 @@ void swap_config(BDS_Edge *e,
break;
}
}
-
+
if(orientation == 1) {
*p11 = p1;
*p12 = p2;
@@ -865,7 +867,7 @@ bool BDS_Mesh::swap_edge(BDS_Edge *e, const BDS_SwapEdgeTest &theTest)
if(e->deleted)
return false;
-
+
int nbFaces = e->numfaces();
if(nbFaces != 2)
return false;
@@ -910,7 +912,7 @@ bool BDS_Mesh::swap_edge(BDS_Edge *e, const BDS_SwapEdgeTest &theTest)
op[1], op[0], p2))
return false;
}
-
+
if(!theTest(p1, p2, op[0], op[1]))
return false;
@@ -971,7 +973,7 @@ int BDS_Edge::numTriangles() const
// moving one of its vertices
static bool test_move_point_parametric_quad(BDS_Point *p, double u, double v, BDS_Face *t)
-{
+{
BDS_Point *pts[4];
t->getNodes(pts);
@@ -983,9 +985,9 @@ static bool test_move_point_parametric_quad(BDS_Point *p, double u, double v, BD
double ori_init1 = robustPredicates::orient2d(pa, pb, pc);
double ori_init2 = robustPredicates::orient2d(pc, pd, pa);
- if(p == pts[0]){
- pa[0] = u;
- pa[1] = v;
+ if(p == pts[0]){
+ pa[0] = u;
+ pa[1] = v;
}
else if(p == pts[1]){
pb[0] = u;
@@ -1003,7 +1005,7 @@ static bool test_move_point_parametric_quad(BDS_Point *p, double u, double v, BD
Msg::Error("Something wrong in move_point_parametric_quad");
return false;
}
-
+
double ori_final1 = robustPredicates::orient2d(pa, pb, pc);
double ori_final2 = robustPredicates::orient2d(pc, pd, pa);
// allow to move a point when a triangle was flat
@@ -1011,7 +1013,7 @@ static bool test_move_point_parametric_quad(BDS_Point *p, double u, double v, BD
}
static bool test_move_point_parametric_triangle(BDS_Point *p, double u, double v, BDS_Face *t)
-{
+{
if (t->e4)
return test_move_point_parametric_quad(p, u, v, t);
BDS_Point *pts[4];
@@ -1030,9 +1032,9 @@ static bool test_move_point_parametric_triangle(BDS_Point *p, double u, double v
double ori_init = robustPredicates::orient2d(pa, pb, pc);
- if(p == pts[0]){
- pa[0] = u;
- pa[1] = v;
+ if(p == pts[0]){
+ pa[0] = u;
+ pa[1] = v;
}
else if(p == pts[1]){
pb[0] = u;
@@ -1047,7 +1049,7 @@ static bool test_move_point_parametric_triangle(BDS_Point *p, double u, double v
a[0] = pb[0] - pa[0]; a[1] = pb[1] - pa[1];
b[0] = pc[0] - pa[0]; b[1] = pc[1] - pa[1];
-
+
double area_final = fabs(a[0] * b[1] - a[1] * b[0]);
if(area_final < 0.1 * area_init) return false;
double ori_final = robustPredicates::orient2d(pa, pb, pc);
@@ -1081,7 +1083,7 @@ bool BDS_Mesh::collapse_edge_parametric(BDS_Edge *e, BDS_Point *p)
BDS_Point *pt[3][1024];
BDS_GeomEntity *gs[1024];
int ept[2][1024];
- BDS_GeomEntity *egs[1024];
+ BDS_GeomEntity *egs[1024];
int nt = 0;
{
p->getTriangles(t);
@@ -1191,7 +1193,7 @@ bool BDS_Mesh::smooth_point_centroid(BDS_Point *p, GFace *gf, bool test_quality)
while(ited != itede) {
BDS_Edge *e = *ited;
BDS_Point *n = e->othervertex(p);
- double fact = 1.0;
+ double fact = 1.0;
sTot += fact;
U += n->u * fact;
V += n->v * fact;
@@ -1201,7 +1203,7 @@ bool BDS_Mesh::smooth_point_centroid(BDS_Point *p, GFace *gf, bool test_quality)
LC += n->lc() * fact;
++ited;
}
- U /= (sTot);
+ U /= (sTot);
V /= (sTot);
LC /= (sTot);
XX/= (sTot);
@@ -1209,14 +1211,14 @@ bool BDS_Mesh::smooth_point_centroid(BDS_Point *p, GFace *gf, bool test_quality)
ZZ/= (sTot);
GPoint gp;double uv[2];
- if (isSphere){
+ if (isSphere){
gp = gf->closestPoint(SPoint3(XX, YY, ZZ), uv);
U = gp.u();
V = gp.v();
}
else
gp = gf->point(U * scalingU, V * scalingV);
-
+
if (!gp.succeeded()){
return false;
}
@@ -1231,16 +1233,16 @@ bool BDS_Mesh::smooth_point_centroid(BDS_Point *p, GFace *gf, bool test_quality)
double newWorst = 1.0;
double oldWorst = 1.0;
while(it != ite) {
- BDS_Face *t = *it;
+ BDS_Face *t = *it;
BDS_Point *n[4];
t->getNodes(n);
p->u = U;
p->v = V;
- double snew = fabs(surface_triangle_param(n[0], n[1], n[2]));
+ double snew = fabs(surface_triangle_param(n[0], n[1], n[2]));
s1 += snew;
p->u = oldU;
p->v = oldV;
- double sold = fabs(surface_triangle_param(n[0], n[1], n[2]));
+ double sold = fabs(surface_triangle_param(n[0], n[1], n[2]));
s2 += sold;
// printf("%22.15E %22.15E\n", snew, sold);
if(snew < .1 * sold) return false;
@@ -1257,20 +1259,20 @@ bool BDS_Mesh::smooth_point_centroid(BDS_Point *p, GFace *gf, bool test_quality)
p->Z = oldZ;
normal_triangle(n[0], n[1], n[2], norm2);
oldWorst = std::min(oldWorst, qmTriangle(*it, QMTRI_RHO));
- double ps;
+ double ps;
prosca(norm1, norm2, &ps);
if(ps < .5) return false;
}
++it;
}
-
+
// printf("%22.15E %22.15E %22.15E\n",s1,s2,fabs(s2-s1));
if(fabs(s2-s1) > 1.e-14 * (s2 + s1)) return false;
-
+
if(test_quality && newWorst < oldWorst){
return false;
}
-
+
p->u = U;
p->v = V;
p->lc() = LC;
@@ -1281,7 +1283,7 @@ bool BDS_Mesh::smooth_point_centroid(BDS_Point *p, GFace *gf, bool test_quality)
while(eit != p->edges.end()) {
(*eit)->update();
++eit;
- }
+ }
return true;
}
@@ -1291,10 +1293,10 @@ bool BDS_Mesh::smooth_point_parametric(BDS_Point *p, GFace *gf)
if(!p->config_modified)return false;
if(p->g && p->g->classif_degree <= 1)
return false;
-
+
double U = 0;
double V = 0;
- double tot_length = 0;
+ double tot_length = 0;
double LC = 0;
@@ -1311,10 +1313,10 @@ bool BDS_Mesh::smooth_point_parametric(BDS_Point *p, GFace *gf)
U += n[i]->u;
V += n[i]->v;
LC += n[i]->lc();
- tot_length += 1;
+ tot_length += 1;
}
++it;
- }
+ }
U /= tot_length;
V /= tot_length;
LC /= p->edges.size();
diff --git a/Mesh/CMakeLists.txt b/Mesh/CMakeLists.txt
index 8322189c6d..bb03a19a77 100644
--- a/Mesh/CMakeLists.txt
+++ b/Mesh/CMakeLists.txt
@@ -4,10 +4,8 @@
# bugs and problems to <gmsh@geuz.org>.
set(SRC
-CenterlineField.cpp
Generator.cpp
- Field.cpp
- meshGEdge.cpp
+ meshGEdge.cpp
meshGEdgeExtruded.cpp
meshGFace.cpp
meshGFaceTransfinite.cpp meshGFaceExtruded.cpp
@@ -36,6 +34,8 @@ CenterlineField.cpp
Levy3D.cpp
periodical.cpp
directions3D.cpp
+ Field.cpp
+ CenterlineField.cpp
)
file(GLOB HDR RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.h)
diff --git a/Mesh/CenterlineField.cpp b/Mesh/CenterlineField.cpp
index 18d73f6021..2bc4531fb9 100644
--- a/Mesh/CenterlineField.cpp
+++ b/Mesh/CenterlineField.cpp
@@ -13,6 +13,7 @@
#include <map>
#include <set>
#include <list>
+#include <algorithm>
#include <string>
#include <fstream>
#include <sstream>
@@ -37,20 +38,19 @@
#include "Curvature.h"
#include "MElement.h"
#include "Context.h"
+
#if defined(HAVE_ANN)
#include <ANN/ANN.h>
-#include <algorithm>
-void erase(std::vector<MLine*>& lines, MLine* l) {
+static void erase(std::vector<MLine*>& lines, MLine* l)
+{
std::vector<MLine*>::iterator it = std::remove(lines.begin(), lines.end(), l);
lines.erase(it, lines.end());
}
-
-
-double computeLength(std::vector<MLine*> lines){
-
+static double computeLength(std::vector<MLine*> lines)
+{
double length= 0.0;
for (unsigned int i = 0; i< lines.size(); i++){
length += lines[i]->getLength();
@@ -58,8 +58,8 @@ double computeLength(std::vector<MLine*> lines){
return length;
}
-bool isClosed(std::set<MEdge, Less_Edge> &theCut){
-
+static bool isClosed(std::set<MEdge, Less_Edge> &theCut)
+{
std::multiset<MVertex*> boundV;
std::set<MEdge,Less_Edge>::iterator it = theCut.begin();
for (; it!= theCut.end(); it++){
@@ -73,34 +73,10 @@ bool isClosed(std::set<MEdge, Less_Edge> &theCut){
}
}
return true;
-
- // std::list<MEdge> segments;
- // std::map<MVertex*, MEdge> boundv;
- // std::set<MEdge,Less_Edge>::iterator it = theCut.begin();
- // for (; it!= theCut.end(); it++){
- // segments.push_back(*it);
- // }
-
- // // find a lonely MEdge
- // for (std::list<MEdge>::iterator it = segments.begin();
- // it != segments.end(); ++it){
- // MVertex *vL = it->getVertex(0);
- // MVertex *vR = it->getVertex(1);
- // std::map<MVertex*,MEdge>::iterator it1 = boundv.find(vL);
- // if (it1 == boundv.end()) boundv.insert(std::make_pair(vL,*it));
- // else boundv.erase(it1);
- // std::map<MVertex*,MEdge>::iterator it2 = boundv.find(vR);
- // if (it2 == boundv.end()) boundv.insert(std::make_pair(vR,*it));
- // else boundv.erase(it2);
- // }
-
- // if (boundv.size() == 0) return true;
- // else return false;
-
}
-void orderMLines(std::vector<MLine*> &lines, MVertex *vB, MVertex *vE){
-
+static void orderMLines(std::vector<MLine*> &lines, MVertex *vB, MVertex *vE)
+{
std::vector<MLine*> _m;
std::list<MLine*> segments;
std::map<MVertex*, MLine*> boundv;
@@ -198,14 +174,6 @@ void orderMLines(std::vector<MLine*> &lines, MVertex *vB, MVertex *vE){
for (unsigned int i = 0; i < lines.size(); i++)
_orientation[i] = !_orientation[i];
}
-
- // if (junctions.find(lines[0]->getVertex(0)) != junctions.end()) vB = lines[0]->getVertex(0);
- // else if (junctions.find(lines[0]->getVertex(1)) != junctions.end()) vB = lines[0]->getVertex(1);
- // else printf("no vB junc found %d %d \n", lines[0]->getVertex(0)->getNum(), lines[0]->getVertex(1)->getNum());
- // if (junctions.find(lines[lines.size()-1]->getVertex(0)) != junctions.end()) vE = lines[lines.size()-1]->getVertex(0);
- // else if (junctions.find(lines[lines.size()-1]->getVertex(1)) != junctions.end()) vE = lines[lines.size()-1]->getVertex(1);
- // else printf("no vE junc found %d %d \n", lines[0]->getVertex(0)->getNum(), lines[0]->getVertex(1)->getNum());
- // printf("in order vB= %d =%d \n", vB->getNum(), vE->getNum());
}
static void recurConnectByMEdge(const MEdge &e,
@@ -229,13 +197,12 @@ static void recurConnectByMEdge(const MEdge &e,
}
}
-
-void cutTriangle(MTriangle *tri,
- std::map<MEdge,MVertex*,Less_Edge> &cutEdges,
- std::vector<MVertex*> &cutVertices,
- std::vector<MTriangle*> &newTris,
- std::set<MEdge,Less_Edge> &newCut){
-
+static void cutTriangle(MTriangle *tri,
+ std::map<MEdge,MVertex*,Less_Edge> &cutEdges,
+ std::vector<MVertex*> &cutVertices,
+ std::vector<MTriangle*> &newTris,
+ std::set<MEdge,Less_Edge> &newCut)
+{
MVertex *c[3] = {0,0,0};
for (int j=0;j<3;j++){
MEdge ed = tri->getEdge(j);
@@ -319,8 +286,8 @@ void cutTriangle(MTriangle *tri,
}
}
-Centerline::Centerline(std::string fileName): kdtree(0), kdtreeR(0), nodes(0), nodesR(0){
-
+Centerline::Centerline(std::string fileName): kdtree(0), kdtreeR(0), nodes(0), nodesR(0)
+{
recombine = CTX::instance()->mesh.recombineAll;
index = new ANNidx[1];
@@ -339,8 +306,9 @@ Centerline::Centerline(std::string fileName): kdtree(0), kdtreeR(0), nodes(0), n
is_extruded = 0;
}
-Centerline::Centerline(): kdtree(0), kdtreeR(0), nodes(0), nodesR(0){
+Centerline::Centerline(): kdtree(0), kdtreeR(0), nodes(0), nodesR(0)
+{
index = new ANNidx[1];
dist = new ANNdist[1];
@@ -353,19 +321,28 @@ Centerline::Centerline(): kdtree(0), kdtreeR(0), nodes(0), nodesR(0){
is_closed = 0;
is_extruded = 0;
- options["closeVolume"] = new FieldOptionInt(is_closed, "Action: Create In/Outlet planar faces");
- options["extrudeWall"] = new FieldOptionInt(is_extruded, "Action: Extrude wall");
- options["reMesh"] = new FieldOptionInt(is_cut, "Action: Cut the initial mesh in different mesh partitions using the centerlines");
- callbacks["run"] = new FieldCallbackGeneric<Centerline>(this, &Centerline::run, "Run actions (closeVolume, extrudeWall, cutMesh) \n");
-
- options["FileName"] = new FieldOptionString (fileName, "File name for the centerlines", &update_needed);
- options["nbPoints"] = new FieldOptionInt(nbPoints, "Number of mesh elements in a circle");
- options["nbElemLayer"] = new FieldOptionInt(nbElemLayer, "Number of mesh elements the extruded layer");
- options["hLayer"] = new FieldOptionDouble(hLayer, "Thickness (% of radius) of the extruded layer");
-
+ options["closeVolume"] = new FieldOptionInt
+ (is_closed, "Action: Create In/Outlet planar faces");
+ options["extrudeWall"] = new FieldOptionInt
+ (is_extruded, "Action: Extrude wall");
+ options["reMesh"] = new FieldOptionInt
+ (is_cut, "Action: Cut the initial mesh in different mesh partitions using the "
+ "centerlines");
+ callbacks["run"] = new FieldCallbackGeneric<Centerline>
+ (this, &Centerline::run, "Run actions (closeVolume, extrudeWall, cutMesh) \n");
+
+ options["FileName"] = new FieldOptionString
+ (fileName, "File name for the centerlines", &update_needed);
+ options["nbPoints"] = new FieldOptionInt
+ (nbPoints, "Number of mesh elements in a circle");
+ options["nbElemLayer"] = new FieldOptionInt
+ (nbElemLayer, "Number of mesh elements the extruded layer");
+ options["hLayer"] = new FieldOptionDouble
+ (hLayer, "Thickness (% of radius) of the extruded layer");
}
-Centerline::~Centerline(){
+Centerline::~Centerline()
+{
if (mod) delete mod;
if(kdtree) delete kdtree;
if(kdtreeR) delete kdtreeR;
@@ -373,11 +350,10 @@ Centerline::~Centerline(){
if(nodesR) annDeallocPts(nodesR);
delete[]index;
delete[]dist;
-
}
-void Centerline::importFile(std::string fileName){
-
+void Centerline::importFile(std::string fileName)
+{
current = GModel::current();
std::vector<GFace*> currentFaces = current->bindingsGetFaces();
for (unsigned int i = 0; i < currentFaces.size(); i++){
@@ -426,12 +402,10 @@ void Centerline::importFile(std::string fileName){
}
createBranches(maxN);
-
}
-
-void Centerline::createBranches(int maxN){
-
+void Centerline::createBranches(int maxN)
+{
//sort colored lines and create edges
std::vector<std::vector<MLine*> > color_edges;
color_edges.resize(maxN);
@@ -498,15 +472,20 @@ void Centerline::createBranches(int maxN){
}
else itl++;
}
- if (vB == vE) { Msg::Error("Begin and end points branch are the same \n");break;}
+ if (vB == vE) {
+ Msg::Error("Begin and end points branch are the same \n");
+ break;
+ }
orderMLines(myLines, vB, vE);
std::vector<Branch> children;
- Branch newBranch ={ tag++, myLines, computeLength(myLines), vB, vE, children, 1.e6, 0.0};
+ Branch newBranch ={ tag++, myLines, computeLength(myLines), vB, vE,
+ children, 1.e6, 0.0};
edges.push_back(newBranch) ;
}
}
- Msg::Info("Centerline: in/outlets =%d branches =%d ", (int)color_edges.size()+1, (int)edges.size());
+ Msg::Info("Centerline: in/outlets =%d branches =%d ",
+ (int)color_edges.size()+1, (int)edges.size());
//create children
for(unsigned int i = 0; i < edges.size(); ++i) {
@@ -528,7 +507,8 @@ void Centerline::createBranches(int maxN){
}
-void Centerline::distanceToSurface(){
+void Centerline::distanceToSurface()
+{
Msg::Info("Centerline: computing distance to surface mesh ");
//COMPUTE WITH REVERSE ANN TREE (SURFACE POINTS IN TREE)
@@ -561,8 +541,9 @@ void Centerline::distanceToSurface(){
}
}
-void Centerline::computeRadii(){
+void Centerline::computeRadii()
+{
for(unsigned int i = 0; i < edges.size(); ++i) {
std::vector<MLine*> lines = edges[i].lines;
for(unsigned int j = 0; j < lines.size(); ++j) {
@@ -577,8 +558,9 @@ void Centerline::computeRadii(){
}
}
-void Centerline::buildKdTree(){
+void Centerline::buildKdTree()
+{
FILE * f = fopen("myPOINTS.pos","w");
fprintf(f, "View \"\"{\n");
@@ -619,11 +601,10 @@ void Centerline::buildKdTree(){
}
fprintf(f,"};\n");
fclose(f);
-
}
-void Centerline::createSplitCompounds(){
-
+void Centerline::createSplitCompounds()
+{
//number of discrete vertices, edges, faces and regions for the mesh
NV = current->getMaxElementaryNumber(0);
NE = current->getMaxElementaryNumber(1);
@@ -641,7 +622,7 @@ 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);
+ current->addCompoundEdge(e_compound,num_gec);
}
// Parametrize Compound surfaces
@@ -654,19 +635,18 @@ void Centerline::createSplitCompounds(){
Msg::Info("Create Compound Surface (%d) = %d discrete face",
num_gfc, pf->tag());
- GFace *gfc = current->addCompoundFace(f_compound, 7, 0, num_gfc); //1=conf_spectral 4=convex_circle, 7=conf_fe
+ //1=conf_spectral 4=convex_circle, 7=conf_fe
+ GFace *gfc = current->addCompoundFace(f_compound, 7, 0, num_gfc);
gfc->meshAttributes.recombine = recombine;
gfc->addPhysicalEntity(1);
current->setPhysicalName("wall", 2, 1);//tag 1
}
-
}
-void Centerline::cleanMesh(){
-
-
+void Centerline::cleanMesh()
+{
return; // does not work yet !
////////////////////////////////
@@ -740,8 +720,9 @@ void Centerline::cleanMesh(){
current->createTopologyFromMesh();
}
-void Centerline::createFaces(){
+void Centerline::createFaces()
+{
std::vector<std::vector<MTriangle*> > faces;
std::multimap<MEdge, MTriangle*, Less_Edge> e2e;
@@ -767,7 +748,8 @@ void Centerline::createFaces(){
it != touched.end(); ++it)
e2e.erase(*it);
}
- Msg::Info("Centerline: action (cutMesh) has cut surface mesh in %d faces ", (int)faces.size());
+ Msg::Info("Centerline: action (cutMesh) has cut surface mesh in %d faces ",
+ (int)faces.size());
//create discFaces
for(unsigned int i = 0; i < faces.size(); ++i){
@@ -789,11 +771,10 @@ void Centerline::createFaces(){
f->mesh_vertices.insert(f->mesh_vertices.begin(),
myVertices.begin(), myVertices.end());
}
-
}
-void Centerline::createClosedVolume(GEdge *gin, std::vector<GEdge*> boundEdges){
-
+void Centerline::createClosedVolume(GEdge *gin, std::vector<GEdge*> boundEdges)
+{
current->setFactory("Gmsh");
std::vector<std::vector<GFace *> > myFaceLoops;
std::vector<GFace *> myFaces;
@@ -817,7 +798,8 @@ void Centerline::createClosedVolume(GEdge *gin, std::vector<GEdge*> boundEdges){
myFaces.push_back(newFace);
}
- Msg::Info("Centerline: action (closeVolume) has created %d in/out planar faces ", (int)boundEdges.size());
+ Msg::Info("Centerline: action (closeVolume) has created %d in/out planar faces ",
+ (int)boundEdges.size());
for (int i = 0; i < NF; i++){
GFace * gf;
@@ -834,8 +816,8 @@ void Centerline::createClosedVolume(GEdge *gin, std::vector<GEdge*> boundEdges){
}
-void Centerline::extrudeBoundaryLayerWall(GEdge* gin, std::vector<GEdge*> boundEdges){
-
+void Centerline::extrudeBoundaryLayerWall(GEdge* gin, std::vector<GEdge*> boundEdges)
+{
Msg::Info("Centerline: extrude boundary layer wall (%d, %g%%R) ", nbElemLayer, hLayer);
//orient extrude direction outward
@@ -853,19 +835,20 @@ void Centerline::extrudeBoundaryLayerWall(GEdge* gin, std::vector<GEdge*> boundE
int shift = 0;
if(is_cut) shift = NE;
for (int i= 0; i< NF; i++){
- GFace *gfc ;
+ GFace *gfc ;
if (is_cut) gfc = current->getFaceByTag(NF+i+1);
else gfc = current->getFaceByTag(i+1);
current->setFactory("Gmsh");
//view -5 to scale hLayer by radius in BoundaryLayers.cpp
- std::vector<GEntity*> extrudedE = current->extrudeBoundaryLayer(gfc, nbElemLayer, hLayer, dir, -5);
+ std::vector<GEntity*> extrudedE = current->extrudeBoundaryLayer
+ (gfc, nbElemLayer, hLayer, dir, -5);
GFace *eFace = (GFace*) extrudedE[0];
eFace->addPhysicalEntity(5);
current->setPhysicalName("outerWall", 2, 5);//tag 5
GRegion *eRegion = (GRegion*) extrudedE[1];
eRegion->addPhysicalEntity(6);
current->setPhysicalName("wallVolume", 3, 6);//tag 6
- for (int j= 2; j < extrudedE.size(); j++){
+ for (unsigned int j = 2; j < extrudedE.size(); j++){
GFace *elFace = (GFace*) extrudedE[j];
std::list<GEdge*> l_edges = elFace->edges();
for(std::list<GEdge*>::iterator it = l_edges.begin(); it != l_edges.end(); it++){
@@ -881,20 +864,20 @@ void Centerline::extrudeBoundaryLayerWall(GEdge* gin, std::vector<GEdge*> boundE
current->setPhysicalName("outletRings", 2, 8);//tag 8
}
}
- }
+ }
}
}
}
-
-void Centerline::run(){
-
+void Centerline::run()
+{
double t1 = Cpu();
if (update_needed){
std::ifstream input;
input.open(fileName.c_str());
- if(StatFile(fileName)) Msg::Fatal("Centerline file '%s' does not exist", fileName.c_str());
+ if(StatFile(fileName))
+ Msg::Fatal("Centerline file '%s' does not exist", fileName.c_str());
importFile(fileName);
buildKdTree();
update_needed = false;
@@ -909,9 +892,9 @@ void Centerline::run(){
NV = current->getMaxElementaryNumber(0);
NE = current->getMaxElementaryNumber(1);
NF = current->getMaxElementaryNumber(2);
- NR = current->getMaxElementaryNumber(3);
+ NR = current->getMaxElementaryNumber(3);
}
-
+
//identify the boundary edges by looping over all discreteFaces
std::vector<GEdge*> boundEdges;
double dist_inlet = 1.e6;
@@ -920,7 +903,8 @@ void Centerline::run(){
GFace *gf = current->getFaceByTag(i+1);
std::list<GEdge*> l_edges = gf->edges();
for(std::list<GEdge*>::iterator it = l_edges.begin(); it != l_edges.end(); it++){
- std::vector<GEdge*>::iterator ite = std::find(boundEdges.begin(), boundEdges.end(), *it);
+ std::vector<GEdge*>::iterator ite = std::find(boundEdges.begin(),
+ boundEdges.end(), *it);
if (ite != boundEdges.end()) boundEdges.erase(ite);
else boundEdges.push_back(*it);
GVertex *gv = (*it)->getBeginVertex();
@@ -929,7 +913,7 @@ void Centerline::run(){
if(dist < dist_inlet){
dist_inlet = dist;
gin = *it;
- }
+ }
}
}
@@ -938,12 +922,12 @@ void Centerline::run(){
double t2 = Cpu();
Msg::Info("Centerline operators computed in %g (s) ",t2-t1);
-
}
-void Centerline::cutMesh(){
-
- Msg::Info("Centerline: action (cutMesh) splits surface mesh (%d tris) using %s ", triangles.size(), fileName.c_str());
+void Centerline::cutMesh()
+{
+ Msg::Info("Centerline: action (cutMesh) splits surface mesh (%d tris) using %s ",
+ triangles.size(), fileName.c_str());
//splitMesh
for(unsigned int i = 0; i < edges.size(); i++){
@@ -955,8 +939,8 @@ void Centerline::cutMesh(){
//if ( edges[i].children.size()) printf("children (%d) = ", edges[i].children.size());
//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);
if( nbSplit > 1 ){
printf("->> cut branch in %d parts \n", nbSplit);
double li = L/nbSplit;
@@ -969,7 +953,7 @@ void Centerline::cutMesh(){
SVector3 pt(v1->x(), v1->y(), v1->z());
SVector3 dir(v2->x()-v1->x(),v2->y()-v1->y(),v2->z()-v1->z());
std::map<MLine*,double>::iterator itr = radiusl.find(lines[j]);
- bool cutted = cutByDisk(pt, dir, itr->second);
+ cutByDisk(pt, dir, itr->second);
nbSplit--;
lc = 0.0;
}
@@ -985,16 +969,17 @@ void Centerline::cutMesh(){
SVector3 dir(v2->x()-v1->x(),v2->y()-v1->y(),v2->z()-v1->z());
printf("-->> cut branch at bifurcation \n");
std::map<MLine*,double>::iterator itr = radiusl.find(lines[lines.size()-1]);
- bool cutted = cutByDisk(pt, dir, itr->second);
+ //bool cutted =
+ cutByDisk(pt, dir, itr->second);
// if(!cutted){
- // int l = lines.size()-1-lines.size()/(4*nbSplit); //chech this!
- // v1 = lines[l]->getVertex(1);
- // v2 = lines[l]->getVertex(0);
- // pt = SVector3(v1->x(), v1->y(), v1->z());
- // dir = SVector3(v2->x()-v1->x(),v2->y()-v1->y(),v2->z()-v1->z());
- // printf("-->> cut bifurcation NEW \n");
- // itr = radiusl.find(lines[l]);
- // cutted = cutByDisk(pt, dir, itr->second);
+ // int l = lines.size()-1-lines.size()/(4*nbSplit); //chech this!
+ // v1 = lines[l]->getVertex(1);
+ // v2 = lines[l]->getVertex(0);
+ // pt = SVector3(v1->x(), v1->y(), v1->z());
+ // dir = SVector3(v2->x()-v1->x(),v2->y()-v1->y(),v2->z()-v1->z());
+ // printf("-->> cut bifurcation NEW \n");
+ // itr = radiusl.find(lines[l]);
+ // cutted = cutByDisk(pt, dir, itr->second);
// }
}
}
@@ -1012,11 +997,10 @@ void Centerline::cutMesh(){
createSplitCompounds();
Msg::Info("Splitting mesh by centerlines done ");
-
}
-bool Centerline::cutByDisk(SVector3 &PT, SVector3 &NORM, double &maxRad){
-
+bool Centerline::cutByDisk(SVector3 &PT, SVector3 &NORM, double &maxRad)
+{
double a = NORM.x();
double b = NORM.y();
double c = NORM.z();
@@ -1026,7 +1010,7 @@ bool Centerline::cutByDisk(SVector3 &PT, SVector3 &NORM, double &maxRad){
int maxStep = 20;
const double EPS = 0.007;
- //std::set<MEdge,Less_Edge> allEdges;
+ //std::set<MEdge,Less_Edge> allEdges;
std::vector<MEdge> allEdges;
for(unsigned int i = 0; i < triangles.size(); i++){
for ( unsigned int j= 0; j < 3; j++){
@@ -1051,7 +1035,7 @@ bool Centerline::cutByDisk(SVector3 &PT, SVector3 &NORM, double &maxRad){
// for (std::set<MEdge,Less_Edge>::iterator it = allEdges.begin();
// it != allEdges.end() ; ++it){
// MEdge me = *it;
- for (int j=0; j < allEdges.size(); j++){
+ for (unsigned int j = 0; j < allEdges.size(); j++){
MEdge me = allEdges[j];
SVector3 P1(me.getVertex(0)->x(),me.getVertex(0)->y(), me.getVertex(0)->z());
SVector3 P2(me.getVertex(1)->x(),me.getVertex(1)->y(), me.getVertex(1)->z());
@@ -1102,7 +1086,7 @@ bool Centerline::cutByDisk(SVector3 &PT, SVector3 &NORM, double &maxRad){
}
}
-
+
if (step < maxStep){
//printf("cutByDisk OK step =%d \n", step);
return true;
@@ -1114,12 +1098,13 @@ bool Centerline::cutByDisk(SVector3 &PT, SVector3 &NORM, double &maxRad){
}
-double Centerline::operator() (double x, double y, double z, GEntity *ge){
-
+double Centerline::operator() (double x, double y, double z, GEntity *ge)
+{
if (update_needed){
std::ifstream input;
input.open(fileName.c_str());
- if(StatFile(fileName)) Msg::Fatal("Centerline file '%s' does not exist", fileName.c_str());
+ if(StatFile(fileName))
+ Msg::Fatal("Centerline file '%s' does not exist", fileName.c_str());
importFile(fileName);
buildKdTree();
update_needed = false;
@@ -1128,10 +1113,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
- //or in the volume
+ //or in the volume
bool isCompound = false;
if(ge){
- if (ge->dim() == 2 && ge->geomType() == GEntity::CompoundSurface) isCompound = true;
+ if (ge->dim() == 2 && ge->geomType() == GEntity::CompoundSurface) isCompound = true;
std::list<GFace*> cFaces;
if (isCompound) cFaces = ((GFaceCompound*)ge)->getCompounds();
if ( ge->dim() == 3 || (ge->dim() == 2 && ge->geomType() == GEntity::Plane) ||
@@ -1151,7 +1136,7 @@ double Centerline::operator() (double x, double y, double z, GEntity *ge){
//double cmax, cmin;
//SVector3 dirMax,dirMin;
//cmax = ge->curvatures(SPoint2(u, v),&dirMax, &dirMin, &cmax,&cmin);
- //cmax = ge->curvatureMax(SPoint2(u,v));
+ //cmax = ge->curvatureMax(SPoint2(u,v));
//double radC = 1./cmax;
double lc = 2*M_PI*rad/nbPoints;
@@ -1161,13 +1146,13 @@ double Centerline::operator() (double x, double y, double z, GEntity *ge){
}
-
-void Centerline::operator() (double x, double y, double z, SMetric3 &metr, GEntity *ge){
-
+void Centerline::operator() (double x, double y, double z, SMetric3 &metr, GEntity *ge)
+{
if (update_needed){
std::ifstream input;
input.open(fileName.c_str());
- if(StatFile(fileName)) Msg::Fatal("Centerline file '%s' does not exist", fileName.c_str());
+ if(StatFile(fileName))
+ Msg::Fatal("Centerline file '%s' does not exist", fileName.c_str());
importFile(fileName);
buildKdTree();
update_needed = false;
@@ -1175,10 +1160,12 @@ 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
+ //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;
+ bool isCompound = (ge->dim() == 2 && ge->geomType() == GEntity::CompoundSurface) ?
+ true : false;
bool onInOutlets = (ge->geomType() == GEntity::Plane) ? true: false;
std::list<GFace*> cFaces;
if (isCompound) cFaces = ((GFaceCompound*)ge)->getCompounds();
@@ -1201,7 +1188,7 @@ void Centerline::operator() (double x, double y, double z, SMetric3 &metr, GEnt
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]);
-
+
delete[]index2;
delete[]dist2;
@@ -1216,16 +1203,15 @@ void Centerline::operator() (double x, double y, double z, SMetric3 &metr, GEnt
SVector3 dMin, dMax;
curvature.vertexNodalValuesAndDirections(vertices[index[0]],&dMax, &dMin, &cMax, &cMin, 1);
curvature.vertexNodalValues(vertices[index[0]], curv, 0);
- double c0 =std::abs(curv);
- double sign = (curv > 0.0) ? -1.0: 1.0;
+ double sign = (curv > 0.0) ? -1.0: 1.0;
double beta = CTX::instance()->mesh.smoothRatio;
double ratio = 1.1;
double thickness = radMax/3.;
double rho = radMax;
double hwall_n = thickness/20.;
- double hwall_t = 2*M_PI*rho/nbPoints;
+ double hwall_t = 2*M_PI*rho/nbPoints;
double hfar = radMax/5.;
double lc_a = 3.*hwall_t;
@@ -1235,36 +1221,34 @@ void Centerline::operator() (double x, double y, double z, SMetric3 &metr, GEnt
SVector3 dir_a = p1-p0; dir_a.normalize();
SVector3 dir_n(xyz[0]-p0.x(), xyz[1]-p0.y(), xyz[2]-p0.z()); dir_n.normalize();
SVector3 dir_cross = crossprod(dir_a,dir_n); dir_cross.normalize();
- SVector3 dir_t1, dir_t2, dir_a1, dir_a2;
+ SVector3 dir_t1, dir_t2, dir_a1, dir_a2;
buildOrthoBasis2(dir_n, dir_t1, dir_t2);
buildOrthoBasis2(dir_a, dir_a1, dir_a2);
-
+
double ll1 = ds*(ratio-1) + hwall_n;
double lc_n = std::min(ll1,hfar);
- double ll2 = hwall_t *(rho+sign*ds)/rho ; //sign * ds*(ratio-1) + hwall_t;
+ double ll2 = hwall_t *(rho+sign*ds)/rho ; //sign * ds*(ratio-1) + hwall_t;
if (ds > thickness) ll2 = hwall_t *(rho+sign*thickness)/rho ; //lc_a; //hfar
double lc_t = std::min(lc_n*CTX::instance()->mesh.anisoMax, ll2); //std::min(ll2,hfar));
-
+
lc_n = std::max(lc_n, CTX::instance()->mesh.lcMin);
lc_n = std::min(lc_n, CTX::instance()->mesh.lcMax);
lc_t = std::max(lc_t, CTX::instance()->mesh.lcMin);
lc_t = std::min(lc_t, CTX::instance()->mesh.lcMax);
-
+
//curvature metric
SMetric3 curvMetric;
if (ds <= thickness){
- metr = metricBasedOnSurfaceCurvature(dMin, dMax, cMin, cMax, radMax, beta, lc_n, lc_t, hwall_t);
+ metr = metricBasedOnSurfaceCurvature(dMin, dMax, cMin, cMax, radMax,
+ beta, lc_n, lc_t, hwall_t);
}
else if (ds > thickness && onInOutlets){
- curvMetric = buildMetricTangentToCurve(dir_n,lc_n,lc_t);
+ curvMetric = buildMetricTangentToCurve(dir_n,lc_n,lc_t);
}
else if (ds > thickness && !onInOutlets){
- //curvMetric = buildMetricTangentToCurve(dir_n,lc_n,lc_a);
- //metr = SMetric3(1./(lc_a*lc_a), 1./(hfar*hfar), 1./(hfar*hfar), dir_a, dir_a1, dir_a2);
- //metr = intersection_conserveM1(metr,curvMetric);
- curvMetric = buildMetricTangentToCurve(dir_n,lc_n,lc_a);
- metr = SMetric3(1./(lc_a*lc_a), 1./(lc_n*lc_n), 1./(lc_t*lc_t), dir_a, dir_n, dir_cross);
- metr = intersection_conserveM1(metr,curvMetric);
+ curvMetric = buildMetricTangentToCurve(dir_n,lc_n,lc_a);
+ metr = SMetric3(1./(lc_a*lc_a), 1./(lc_n*lc_n), 1./(lc_t*lc_t), dir_a, dir_n, dir_cross);
+ metr = intersection_conserveM1(metr,curvMetric);
}
return;
@@ -1275,11 +1259,13 @@ void Centerline::operator() (double x, double y, double z, SMetric3 &metr, GEnt
}
-SMetric3 Centerline::metricBasedOnSurfaceCurvature(SVector3 dirMin, SVector3 dirMax, double cmin, double cmax, double radMax,
- double beta, double lc_n, double lc_t, double hwall_t){
-
+SMetric3 Centerline::metricBasedOnSurfaceCurvature(SVector3 dirMin, SVector3 dirMax,
+ double cmin, double cmax, double radMax,
+ double beta, double lc_n, double lc_t,
+ double hwall_t)
+{
double lcMin = ((2 * M_PI *radMax) /( 20*nbPoints )); //CTX::instance()->mesh.lcMin;
- double lcMax = ((2 * M_PI *radMax) /( 0.05*nbPoints)); //CTX::instance()->mesh.lcMax; //
+ double lcMax = ((2 * M_PI *radMax) /( 0.05*nbPoints)); //CTX::instance()->mesh.lcMax;
if (cmin == 0) cmin =1.e-12;
if (cmax == 0) cmax =1.e-12;
double lambda1 = ((2 * M_PI) /( fabs(cmin) * nbPoints ) );
@@ -1287,12 +1273,12 @@ SMetric3 Centerline::metricBasedOnSurfaceCurvature(SVector3 dirMin, SVector3 dir
double betaM = 1*beta;
double rhoMin = 1./cmin;
- double rhoMax = 1./cmax;
- //double oneOverD2_min = .5/(lc_t*lc_t) * (1. + sqrt (1. + ( 4.*cmin*cmin*lc_t*lc_t*lc_t*lc_t/ (lc_n*lc_n*beta*beta))));
- //double oneOverD2_max = .5/(lc_t*lc_t) * (1. + sqrt (1. + ( 4.*cmax*cmax*lc_t*lc_t*lc_t*lc_t/ (lc_n*lc_n*beta*beta))));
- double oneOverD2_min = 1./(2.*rhoMin*rhoMin*(betaM*betaM-1))*(sqrt(1+ (4.*rhoMin*rhoMin*(betaM*betaM-1))/(lc_n*lc_n))-1.);
- double oneOverD2_max = 1./(2.*rhoMax*rhoMax*(beta*beta-1))*(sqrt(1+ (4.*rhoMax*rhoMax*(beta*beta-1))/(lc_n*lc_n))-1.);
-
+ double rhoMax = 1./cmax;
+ double oneOverD2_min = 1./(2.*rhoMin*rhoMin*(betaM*betaM-1)) *
+ (sqrt(1+ (4.*rhoMin*rhoMin*(betaM*betaM-1))/(lc_n*lc_n))-1.);
+ double oneOverD2_max = 1./(2.*rhoMax*rhoMax*(beta*beta-1)) *
+ (sqrt(1+ (4.*rhoMax*rhoMax*(beta*beta-1))/(lc_n*lc_n))-1.);
+
lambda1 = sqrt(1./oneOverD2_min);
lambda2 = sqrt(1./oneOverD2_max);
@@ -1305,13 +1291,14 @@ SMetric3 Centerline::metricBasedOnSurfaceCurvature(SVector3 dirMin, SVector3 dir
lambda1 = std::min(lambda1, lcMax);
lambda2 = std::min(lambda2, lcMax);
- SMetric3 curvMetric (1./(lambda1*lambda1),1./(lambda2*lambda2), 1./(lc_n*lc_n), dirMin, dirMax, dirNorm); // 1.e-6
+ SMetric3 curvMetric (1./(lambda1*lambda1),1./(lambda2*lambda2),
+ 1./(lc_n*lc_n), dirMin, dirMax, dirNorm); // 1.e-6
return curvMetric;
}
-void Centerline::printSplit() const{
-
+void Centerline::printSplit() const
+{
FILE * f = fopen("mySPLIT.pos","w");
fprintf(f, "View \"\"{\n");
for(unsigned int i = 0; i < edges.size(); ++i){
diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index cb2d12993e..66a20eefde 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -47,8 +47,8 @@ static bool computeEquidistantParameters1(GEdge *ge, double u0, double uN, int N
GPoint p1 = ge->point(uN);
double du = 1. / (N - 1);
u[0] = u0;
- // printf("starting with %g %g %g\n",p0.x(),p0.y(),u0);
- // printf("ending with %g %g %g\n",p1.x(),p1.y(),uN);
+ // printf("starting with %g %g %g\n",p0.x(),p0.y(),u0);
+ // printf("ending with %g %g %g\n",p1.x(),p1.y(),uN);
for (int i = 1; i < N; i++){
SPoint3 pi (p0.x() + i * du * (p1.x()-p0.x()),
p0.y() + i * du * (p1.y()-p0.y()),
@@ -56,7 +56,7 @@ static bool computeEquidistantParameters1(GEdge *ge, double u0, double uN, int N
double t;
GPoint gp = ge->closestPoint(pi, t);
u[i] = gp.u();
- // printf("going to %g %g u %g\n",pi.x(),pi.y(),gp.u());
+ // printf("going to %g %g u %g\n",pi.x(),pi.y(),gp.u());
}
return true;
}
@@ -120,6 +120,7 @@ static bool computeEquidistantParameters0(GEdge *ge, double u0, double uN, int N
}
return false;
}
+
// 1 = geodesics
static int method_for_computing_intermediary_points = 0;
static bool computeEquidistantParameters(GEdge *ge, double u0, double uN, int N,
@@ -145,20 +146,18 @@ static void myresid(int N, GFace *gf, double *u, double *v, fullVector<double> &
static bool computeEquidistantParameters1(GFace *gf, double u0, double uN,
double v0, double vN, int N,
- double *u, double *v){
- // printf("coucou\n");
+ double *u, double *v)
+{
GPoint p0 = gf->point(u0,v0);
GPoint p1 = gf->point(uN,vN);
double du = 1. / (N - 1);
u[0] = u0;
u[0] = u0;
v[0] = v0;
- // printf("starting with %g %g %g\n",p0.x(),p0.y(),u0);
- // printf("ending with %g %g %g\n",p1.x(),p1.y(),uN);
for (int i = 1; i < N; i++){
- SPoint3 pi (p0.x() + i * du * (p1.x()-p0.x()),
- p0.y() + i * du * (p1.y()-p0.y()),
- p0.z() + i * du * (p1.z()-p0.z()));
+ SPoint3 pi(p0.x() + i * du * (p1.x()-p0.x()),
+ p0.y() + i * du * (p1.y()-p0.y()),
+ p0.z() + i * du * (p1.z()-p0.z()));
SPoint2 t;
GPoint gp = gf->closestPoint(pi, t);
u[i] = gp.u();
@@ -255,7 +254,6 @@ static bool computeEquidistantParameters(GFace *gf, double u0, double uN,
return false;
}
-
static void getEdgeVertices(GEdge *ge, MElement *ele, std::vector<MVertex*> &ve,
edgeContainer &edgeVertices, bool linear,
int nPts = 1)
@@ -321,7 +319,6 @@ static void getEdgeVertices(GEdge *ge, MElement *ele, std::vector<MVertex*> &ve,
int count = u0<u1? j + 1 : nPts + 1 - (j + 1);
GPoint pc = ge->point(US[count]);
v = new MEdgeVertex(pc.x(), pc.y(), pc.z(), ge,US[count]);
- // printf("Edge %d(%g) %d(%g) new vertex %g %g at %g\n",v0->getNum(),u0,v1->getNum(),u1,v->x(),v->y(), US[count]);
}
temp.push_back(v);
// this destroys the ordering of the mesh vertices on the edge
@@ -773,7 +770,8 @@ static void getRegionVertices(GRegion *gr, MElement *incomplete, MElement *ele,
Msg::Error("getRegionVertices not implemented for order %i", nPts+1);
break;
}
- start = ( nPts+2) * ( (nPts+2) + 1) * (2*(nPts+2) + 1) / 6 - (nPts-1) * ( (nPts-1) + 1) * (2*(nPts-1) + 1) / 6;
+ start = ( nPts+2) * ( (nPts+2) + 1) * (2*(nPts+2) + 1) / 6 -
+ (nPts-1) * ( (nPts-1) + 1) * (2*(nPts-1) + 1) / 6;
break;
}
@@ -986,23 +984,26 @@ static void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
p->getVertex(3), p->getVertex(4), p->getVertex(5),
ve[0], ve[1], ve[2], ve[3], ve[4], ve[5], ve[6], ve[7], ve[8],
vf[0], vf[1], vf[2]));
- } else {
-
- //getFaceVertices(gr, p, vf, faceVertices, edgeVertices, linear, nPts);
- //ve.insert(ve.end(), vf.begin(), vf.end());
- //MPrismN incpl(p->getVertex(0), p->getVertex(1), p->getVertex(2), p->getVertex(3),
- // p->getVertex(4), p->getVertex(5), ve, nPts + 1);
- //getRegionVertices(gr, &incpl, p, vr, linear, nPts);
- //if (nPts == 0) {
- // printf("Screwed\n");
- // }
- //ve.insert(ve.end(), vr.begin(), vr.end());
- //MPrism* n = new MPrismN(p->getVertex(0), p->getVertex(1), p->getVertex(2),
- // p->getVertex(3), p->getVertex(4), p->getVertex(5),
- // ve, nPts+1);
- //if (!mappingIsInvertible(n))
- // Msg::Warning("Found invalid curved volume element (# %d in list)", i);
- //prisms2.push_back(n);
+ }
+ else {
+ Msg::Error("PrismN generation not implemented");
+ /*
+ getFaceVertices(gr, p, vf, faceVertices, edgeVertices, linear, nPts);
+ ve.insert(ve.end(), vf.begin(), vf.end());
+ MPrismN incpl(p->getVertex(0), p->getVertex(1), p->getVertex(2), p->getVertex(3),
+ p->getVertex(4), p->getVertex(5), ve, nPts + 1);
+ getRegionVertices(gr, &incpl, p, vr, linear, nPts);
+ if (nPts == 0) {
+ printf("Screwed\n");
+ }
+ ve.insert(ve.end(), vr.begin(), vr.end());
+ MPrism* n = new MPrismN(p->getVertex(0), p->getVertex(1), p->getVertex(2),
+ p->getVertex(3), p->getVertex(4), p->getVertex(5),
+ ve, nPts+1);
+ if (!mappingIsInvertible(n))
+ Msg::Warning("Found invalid curved volume element (# %d in list)", i);
+ prisms2.push_back(n);
+ */
}
}
delete p;
@@ -1030,7 +1031,7 @@ static void setHighOrder(GRegion *gr, edgeContainer &edgeVertices,
}
}
else {
- throw;
+ Msg::Error("PyramidN generation not implemented");
/*
getFaceVertices(gr, p, vf, faceVertices, edgeVertices, linear, nPts);
ve.insert(ve.end(), vf.begin(), vf.end());
@@ -1234,40 +1235,36 @@ void printJacobians(GModel *m, const char *nm)
fclose(f);
}
-void getMeshInfoForHighOrder (GModel *gm,
- int &meshOrder,
- bool &complete,
- bool &CAD){
+void getMeshInfoForHighOrder(GModel *gm, int &meshOrder, bool &complete,
+ bool &CAD)
+{
meshOrder = -1;
CAD = true;
for (GModel::riter itr = gm->firstRegion(); itr != gm->lastRegion(); ++itr) {
if ((*itr)->getNumMeshElements()){
- meshOrder = (*itr)->getMeshElement(0)->getPolynomialOrder();
- complete = (meshOrder <= 2) ? 1 : (*itr)->getMeshElement(0)->getNumVolumeVertices();
+ meshOrder = (*itr)->getMeshElement(0)->getPolynomialOrder();
+ complete = (meshOrder <= 2) ? 1 : (*itr)->getMeshElement(0)->getNumVolumeVertices();
break;
- }
+ }
}
for (GModel::fiter itf = gm->firstFace(); itf != gm->lastFace(); ++itf) {
if ((*itf)->getNumMeshElements()){
if (meshOrder == -1) {
- meshOrder = (*itf)->getMeshElement(0)->getPolynomialOrder();
- complete = (meshOrder <= 2) ? 1 : (*itf)->getMeshElement(0)->getNumFaceVertices();
+ meshOrder = (*itf)->getMeshElement(0)->getPolynomialOrder();
+ complete = (meshOrder <= 2) ? 1 : (*itf)->getMeshElement(0)->getNumFaceVertices();
if ((*itf)->geomType() == GEntity::DiscreteSurface)CAD = false;
break;
}
- }
+ }
}
}
-
-
-
-void ElasticAnalogy ( GModel *m, double threshold, bool onlyVisible) {
-
+void ElasticAnalogy ( GModel *m, double threshold, bool onlyVisible)
+{
bool CAD, complete;
int meshOrder;
- getMeshInfoForHighOrder (m,meshOrder, complete, CAD);
+ getMeshInfoForHighOrder (m,meshOrder, complete, CAD);
highOrderTools hot(m);
// now we smooth mesh the internal vertices of the faces
// we do that model face by model face
@@ -1340,7 +1337,8 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete, bool onlyVisi
int counter = 1;
for(GModel::eiter it = m->firstEdge(); it != m->lastEdge(); ++it) {
- Msg::StatusBar(2, true, "Meshing curves order %d (%i/%i)...", order, counter, m->getNumEdges());
+ Msg::StatusBar(2, true, "Meshing curves order %d (%i/%i)...",
+ order, counter, m->getNumEdges());
counter++;
if (onlyVisible && !(*it)->getVisibility())continue;
setHighOrder(*it, edgeVertices, linear, nPts);
@@ -1348,7 +1346,8 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete, bool onlyVisi
counter = 1;
for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it) {
- Msg::StatusBar(2, true, "Meshing surfaces order %d (%i/%i)...", order, counter, m->getNumFaces());
+ Msg::StatusBar(2, true, "Meshing surfaces order %d (%i/%i)...",
+ order, counter, m->getNumFaces());
counter++;
if (onlyVisible && !(*it)->getVisibility())continue;
setHighOrder(*it, edgeVertices, faceVertices, linear, incomplete, nPts);
@@ -1366,7 +1365,8 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete, bool onlyVisi
counter = 1;
for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it) {
- Msg::StatusBar(2, true, "Meshing volumes order %d (%i/%i)...", order, counter, m->getNumRegions());
+ Msg::StatusBar(2, true, "Meshing volumes order %d (%i/%i)...",
+ order, counter, m->getNumRegions());
counter++;
if (onlyVisible && !(*it)->getVisibility())continue;
setHighOrder(*it, edgeVertices, faceVertices, linear, incomplete, nPts);
@@ -1389,35 +1389,36 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete, bool onlyVisi
// hot.ensureMinimumDistorsion(0.1);
checkHighOrderTriangles("Final surface mesh", m, bad, worst);
}
-
+
// m->writeMSH("CORRECTED.msh");
Msg::StatusBar(2, true, "Done meshing order %d (%g s)", order, t2 - t1);
}
-void computeDistanceFromMeshToGeometry (GModel *m, distanceFromMeshToGeometry_t &dist) {
- for (GModel::eiter itEdge = m->firstEdge(); itEdge != m->lastEdge(); ++itEdge) {
+void computeDistanceFromMeshToGeometry (GModel *m, distanceFromMeshToGeometry_t &dist)
+{
+ for (GModel::eiter itEdge = m->firstEdge(); itEdge != m->lastEdge(); ++itEdge) {
double d2,dmax;
- (*itEdge)->computeDistanceFromMeshToGeometry (d2,dmax);
+ (*itEdge)->computeDistanceFromMeshToGeometry (d2,dmax);
dist.d2[*itEdge] = d2;
dist.d_max[*itEdge] = dmax;
}
for (GModel::fiter itFace = m->firstFace(); itFace != m->lastFace(); ++itFace) {
-
+
}
}
-
-void SetHighOrderComplete (GModel *m, bool onlyVisible){
+void SetHighOrderComplete(GModel *m, bool onlyVisible)
+{
faceContainer faceVertices;
for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it) {
- if (onlyVisible && !(*it)->getVisibility())continue;
+ if (onlyVisible && !(*it)->getVisibility()) continue;
std::vector<MTriangle*> newT;
std::vector<MQuadrangle*> newQ;
- for (int i=0;i<(*it)->triangles.size();i++){
+ for (unsigned int i = 0; i < (*it)->triangles.size(); i++){
MTriangle *t = (*it)->triangles[i];
- std::vector<MVertex*> vf, vt;
+ std::vector<MVertex*> vf, vt;
int nPts = t->getPolynomialOrder() - 1;
MTriangle TEMP (t->getVertex(0), t->getVertex(1), t->getVertex(2));
getFaceVertices (*it, t, t, vf, faceVertices, false, nPts);
@@ -1426,74 +1427,75 @@ void SetHighOrderComplete (GModel *m, bool onlyVisible){
MTriangleN *newTr = new MTriangleN(t->getVertex(0), t->getVertex(1), t->getVertex(2),
vt, nPts + 1);
newT.push_back(newTr);
-
+
delete t;
}
(*it)->triangles = newT;
- for (int i=0;i<(*it)->quadrangles.size();i++){
+ for (unsigned int i = 0; i < (*it)->quadrangles.size(); i++){
MQuadrangle *t = (*it)->quadrangles[i];
- std::vector<MVertex*> vf, vt;
+ std::vector<MVertex*> vf, vt;
int nPts = t->getPolynomialOrder() - 1;
MQuadrangle TEMP (t->getVertex(0), t->getVertex(1), t->getVertex(2), t->getVertex(3));
getFaceVertices (*it, t, &TEMP, vf, faceVertices, false, nPts);
for (int j=4;j<t->getNumVertices();j++)vt.push_back(t->getVertex(j));
vt.insert(vt.end(), vf.begin(), vf.end());
- newQ.push_back(new MQuadrangleN(t->getVertex(0), t->getVertex(1), t->getVertex(2), t->getVertex(3),
+ newQ.push_back(new MQuadrangleN(t->getVertex(0), t->getVertex(1),
+ t->getVertex(2), t->getVertex(3),
vt, nPts + 1));
-
delete t;
}
(*it)->quadrangles = newQ;
-
std::set<MVertex*> newV;
- for (int i=0;i<(*it)->getNumMeshElements();++i){
+ for (unsigned int i = 0; i < (*it)->getNumMeshElements(); ++i){
MElement *e = (*it)->getMeshElement(i);
for (int j=0;j<e->getNumVertices();j++)newV.insert(e->getVertex(j));
}
(*it)->mesh_vertices.clear();
(*it)->mesh_vertices.insert((*it)->mesh_vertices.begin(), newV.begin(), newV.end());
-
-
- }
+ }
}
-
-void SetHighOrderInComplete (GModel *m, bool onlyVisible){
+void SetHighOrderInComplete (GModel *m, bool onlyVisible)
+{
std::set<MVertex*> toDelete;
for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it) {
- if (onlyVisible && !(*it)->getVisibility())continue;
+ if (onlyVisible && !(*it)->getVisibility()) continue;
std::vector<MTriangle*> newT;
- for (int i=0;i<(*it)->triangles.size();i++){
+ for (unsigned int i = 0; i < (*it)->triangles.size(); i++){
MTriangle *t = (*it)->triangles[i];
- std::vector<MVertex*> vt;
+ std::vector<MVertex*> vt;
int order = t->getPolynomialOrder();
- for (int j=3;j<t->getNumVertices()-t->getNumFaceVertices();j++)vt.push_back(t->getVertex(j));
- for (int j=t->getNumVertices()-t->getNumFaceVertices();j < t->getNumVertices();j++)toDelete.insert(t->getVertex(j));
+ for (int j=3;j<t->getNumVertices()-t->getNumFaceVertices();j++)
+ vt.push_back(t->getVertex(j));
+ for (int j = t->getNumVertices()-t->getNumFaceVertices();
+ j < t->getNumVertices(); j++)
+ toDelete.insert(t->getVertex(j));
newT.push_back(new MTriangleN(t->getVertex(0), t->getVertex(1), t->getVertex(2),
vt, order));
-
delete t;
}
(*it)->triangles = newT;
std::vector<MQuadrangle*> newQ;
- for (int i=0;i<(*it)->quadrangles.size();i++){
+ for (unsigned int i = 0; i < (*it)->quadrangles.size(); i++){
MQuadrangle *t = (*it)->quadrangles[i];
- std::vector<MVertex*> vt;
+ std::vector<MVertex*> vt;
int nPts = t->getPolynomialOrder() - 1;
- for (int j=4;j<t->getNumVertices()-t->getNumFaceVertices();j++)vt.push_back(t->getVertex(j));
- newQ.push_back(new MQuadrangleN(t->getVertex(0), t->getVertex(1), t->getVertex(2), t->getVertex(3),
- vt, nPts + 1));
+ for (int j = 4; j < t->getNumVertices()-t->getNumFaceVertices(); j++)
+ vt.push_back(t->getVertex(j));
+ newQ.push_back(new MQuadrangleN(t->getVertex(0), t->getVertex(1),
+ t->getVertex(2), t->getVertex(3),
+ vt, nPts + 1));
delete t;
}
(*it)->quadrangles = newQ;
std::vector<MVertex*> newV;
int numd = 0;
- for (int i=0;i<(*it)->mesh_vertices.size();++i){
+ for (unsigned int i = 0; i < (*it)->mesh_vertices.size(); ++i){
if (toDelete.find((*it)->mesh_vertices[i]) == toDelete.end())
newV.push_back((*it)->mesh_vertices[i]);
else{
@@ -1501,9 +1503,6 @@ void SetHighOrderInComplete (GModel *m, bool onlyVisible){
numd++;
}
}
- printf("%d vertices deleted\n");
(*it)->mesh_vertices = newV;
}
-
}
-
diff --git a/Mesh/highOrderTools.cpp b/Mesh/highOrderTools.cpp
index 71df3f6a7a..2e966677a2 100644
--- a/Mesh/highOrderTools.cpp
+++ b/Mesh/highOrderTools.cpp
@@ -42,7 +42,6 @@ void highOrderTools::moveToStraightSidedLocation(MElement *e) const
MVertex *v = e->getVertex(i);
std::map<MVertex*,SVector3>::const_iterator it = _straightSidedLocation.find(v);
if (it != _straightSidedLocation.end()){
- // printf("move %d : %g %g -> %g %g\n",v->getNum(),v->x() ,v->y() , it->second.x() , it->second.y() );
v->x() = it->second.x();
v->y() = it->second.y();
v->z() = it->second.z();
@@ -50,7 +49,6 @@ void highOrderTools::moveToStraightSidedLocation(MElement *e) const
}
}
-
void highOrderTools::ensureMinimumDistorsion(MElement *e, double threshold)
{
double dist = e->distoShapeMeasure();
@@ -97,7 +95,6 @@ void highOrderTools::ensureMinimumDistorsion(MElement *e, double threshold)
// printf("element done\n");
}
-
static void getDistordedElements(const std::vector<MElement*> &v,
const double & threshold,
std::vector<MElement*> &d,
@@ -145,7 +142,6 @@ static void addOneLayer(const std::vector<MElement*> &v,
}
}
-
double highOrderTools::applySmoothingTo(GFace *gf, double tres, bool mixed)
{
if (!gf){
@@ -159,7 +155,6 @@ double highOrderTools::applySmoothingTo(GFace *gf, double tres, bool mixed)
return applySmoothingTo(v,tres, mixed);
}
-
void highOrderTools::ensureMinimumDistorsion (double threshold)
{
std::vector<MElement*> v;
@@ -179,10 +174,9 @@ void highOrderTools::ensureMinimumDistorsion (double threshold)
ensureMinimumDistorsion(v,threshold);
}
-
/*
- void highOrderTools::applySmoothingTo(GRegion *gr)
- {
+void highOrderTools::applySmoothingTo(GRegion *gr)
+{
std::vector<MElement*> v;
v.insert(v.begin(), gr->tetrahedra.begin(),gr->tetrahedra.end());
v.insert(v.begin(), gr->hexahedra.begin(),gr->hexahedra.end());
@@ -190,14 +184,14 @@ void highOrderTools::ensureMinimumDistorsion (double threshold)
Msg::Info("Smoothing high order mesh : model region %d (%d elements)", gr->tag(),
v.size());
applySmoothingTo(v,gr);
- }
- */
+}
+*/
void highOrderTools::computeMetricInfo(GFace *gf,
- MElement *e,
- fullMatrix<double> &J,
- fullMatrix<double> &JT,
- fullVector<double> &D)
+ MElement *e,
+ fullMatrix<double> &J,
+ fullMatrix<double> &JT,
+ fullVector<double> &D)
{
int nbNodes = e->getNumVertices();
// printf("ELEMENT --\n");
@@ -207,7 +201,7 @@ void highOrderTools::computeMetricInfo(GFace *gf,
// printf("%g %g vs %g %g %g\n",param.x(),param.y(),
// e->getVertex(j)->x(),e->getVertex(j)->y(),e->getVertex(j)->z());
- Pair<SVector3,SVector3> der = gf->firstDer(param);
+ Pair<SVector3,SVector3> der = gf->firstDer(param);
int XJ = j;
int YJ = j + nbNodes;
@@ -229,7 +223,7 @@ void highOrderTools::computeMetricInfo(GFace *gf,
JT(VJ,ZJ) = der.second().z();
SVector3 ss = getSSL(e->getVertex(j));
- GPoint gp = gf->point(param);
+ GPoint gp = gf->point(param);
D(XJ) = (gp.x() - ss.x());
D(YJ) = (gp.y() - ss.y());
D(ZJ) = (gp.z() - ss.z());
@@ -240,7 +234,7 @@ void highOrderTools::computeMetricInfo(GFace *gf,
ss.x(),ss.y(),ss.z(),
e->getVertex(j)->x(),e->getVertex(j)->y(),e->getVertex(j)->z());
*/
-
+
}
}
@@ -339,10 +333,10 @@ void highOrderTools::applySmoothingTo(std::vector<MElement*> & all, GFace *gf)
}
double highOrderTools::smooth_metric_(std::vector<MElement*> & v,
- GFace *gf,
- dofManager<double> &myAssembler,
- std::set<MVertex*> &verticesToMove,
- elasticityTerm &El)
+ GFace *gf,
+ dofManager<double> &myAssembler,
+ std::set<MVertex*> &verticesToMove,
+ elasticityTerm &El)
{
std::set<MVertex*>::iterator it;
double dx = 0.0;
@@ -408,13 +402,13 @@ double highOrderTools::smooth_metric_(std::vector<MElement*> & v,
return dx;
}
-
-highOrderTools::highOrderTools (GModel *gm) : _gm(gm), _dim(2), _tag(111)
+highOrderTools::highOrderTools (GModel *gm) : _gm(gm), _tag(111), _dim(2)
{
computeStraightSidedPositions();
}
-highOrderTools::highOrderTools (GModel *gm, GModel *mesh, int order) : _gm(gm), _dim(2), _tag(111)
+highOrderTools::highOrderTools (GModel *gm, GModel *mesh, int order)
+ : _gm(gm), _tag(111), _dim(2)
{
GeomMeshMatcher::instance()->forceTomatch(gm,mesh,1.e-5);
GeomMeshMatcher::instance()->destroy();
@@ -552,17 +546,18 @@ void highOrderTools::computeStraightSidedPositions ()
}
}
- Msg::Info("highOrderTools has been set up : %d nodes are considered",_straightSidedLocation.size());
+ Msg::Info("highOrderTools has been set up : %d nodes are considered",
+ _straightSidedLocation.size());
}
// apply a displacement that does not create elements that are
// distorted over a value "thres"
double highOrderTools::apply_incremental_displacement (double max_incr,
- std::vector<MElement*> & v,
- bool mixed,
- double thres,
- char *meshName,
- std::vector<MElement*> & disto)
+ std::vector<MElement*> & v,
+ bool mixed,
+ double thres,
+ char *meshName,
+ std::vector<MElement*> & disto)
{
#ifdef HAVE_PETSC
// assume that the mesh is OK, yet already curved
@@ -699,7 +694,7 @@ double highOrderTools::apply_incremental_displacement (double max_incr,
// uncurve elements that are invalid
void highOrderTools::ensureMinimumDistorsion(std::vector<MElement*> &all,
- double threshold)
+ double threshold)
{
for(int tries = 0; tries < 100; tries++){
double minD;
@@ -727,7 +722,8 @@ double highOrderTools::applySmoothingTo (std::vector<MElement*> &all,
// _gm->writeMSH("straightSided.msh");
char sm[] = "sm.msh";
- double percentage_of_what_is_left = apply_incremental_displacement (1., all, mixed, -100000000, sm, all);
+ double percentage_of_what_is_left = apply_incremental_displacement
+ (1., all, mixed, -100000000, sm, all);
// ensureMinimumDistorsion (all, threshold);
return 1.;
@@ -735,9 +731,11 @@ double highOrderTools::applySmoothingTo (std::vector<MElement*> &all,
while(1){
char NN[256];
sprintf(NN,"smoothing-%d.msh",ITER++);
- percentage_of_what_is_left = apply_incremental_displacement (1.,all, mixed, threshold,NN,all);
+ percentage_of_what_is_left = apply_incremental_displacement
+ (1.,all, mixed, threshold,NN,all);
percentage += (1.-percentage) * percentage_of_what_is_left/100.;
- Msg::Info("The smoother was able to do %3d percent of the motion",(int)(percentage*100.));
+ Msg::Info("The smoother was able to do %3d percent of the motion",
+ (int)(percentage*100.));
if (percentage_of_what_is_left == 0.0) break;
else if (percentage_of_what_is_left == 100.)break;
}
@@ -749,9 +747,10 @@ double highOrderTools::applySmoothingTo (std::vector<MElement*> &all,
}
extern void printJacobians(GModel *m, const char *nm);
-void highOrderTools::makePosViewWithJacobians (const char *fn){
+
+void highOrderTools::makePosViewWithJacobians (const char *fn)
+{
printJacobians(_gm,fn);
}
#endif
-
diff --git a/Mesh/meshGFaceBoundaryLayers.cpp b/Mesh/meshGFaceBoundaryLayers.cpp
index cd830de491..722db88afd 100644
--- a/Mesh/meshGFaceBoundaryLayers.cpp
+++ b/Mesh/meshGFaceBoundaryLayers.cpp
@@ -135,7 +135,7 @@ BoundaryLayerColumns* buidAdditionalPoints2D (GFace *gf)
// assume that the initial mesh has been created i.e. that there exist
// triangles inside the domain. Triangles are used to define
// exterior normals
- for (int i=0;i<gf->triangles.size();i++){
+ for (unsigned int i = 0; i < gf->triangles.size(); i++){
SPoint2 p0,p1,p2;
MVertex *v0 = gf->triangles[i]->getVertex(0);
MVertex *v1 = gf->triangles[i]->getVertex(1);
@@ -230,10 +230,9 @@ BoundaryLayerColumns* buidAdditionalPoints2D (GFace *gf)
for (std::multimap<MEdge,SVector3,Less_Edge>::iterator itm =
_columns->_normals.lower_bound(e2);
itm != _columns->_normals.upper_bound(e2); ++itm) N2.push_back(itm->second);
- double LL;
if (N1.size() == N2.size()){
// if (N1.size() > 1)printf("%d sides\n",N1.size());
- for (int SIDE = 0; SIDE < N1.size() ; SIDE++){
+ for (unsigned int SIDE = 0; SIDE < N1.size() ; SIDE++){
// IF THE ANGLE IS GREATER THAN THRESHOLD, ADD DIRECTIONS !!
double angle = computeAngle (gf,e1,e2,N1[SIDE],N2[SIDE]);
// if (N1.size() > 1)printf("angle = %g\n",angle);
@@ -293,7 +292,7 @@ BoundaryLayerColumns* buidAdditionalPoints2D (GFace *gf)
// if (_dirs.size() > 1)printf("%d directions\n",_dirs.size());
// now create the BL points
- for (int DIR=0;DIR<_dirs.size();DIR++){
+ for (unsigned int DIR=0;DIR<_dirs.size();DIR++){
SPoint2 p;
SVector3 n = _dirs[DIR];
@@ -328,10 +327,7 @@ BoundaryLayerColumns* buidAdditionalPoints2D (GFace *gf)
buildMeshMetric(gf, poffset, m, metric);
const double l2 = n.x()*n.x()*metric[0] + 2*n.x()*n.y()*metric[1] + n.y()*n.y()*metric[2] ;
l = 1./sqrt(l2);
- // if (_dirs.size() > 1) printf("l = %g metric = %g %g %g dim %d tag %d \n",l,metric[0],metric[1],metric[2],current->onWhat()->dim(),current->onWhat()->tag());
- // printf("%g %g\n",l,LL);
if (l >= blf->hfar){
- // printf("stopping %g %g\n",l,LL);
break;
}
// printf("%g %g %g \n",current->x(),current->y(),blf->current_distance);
@@ -364,7 +360,7 @@ BoundaryLayerColumns* buidAdditionalPoints2D (GFace *gf)
// printf("pnew %g %g new point %g %g n %g %g\n",pnew.x(),pnew.y(),gp.x(),gp.y(),n.x(),n.y());
_metrics.push_back(m);
// const double l = n[0]*m(0,0) +;
- if (_column.size() > nbCol)break; // FIXME
+ if ((int)_column.size() > nbCol) break; // FIXME
p = pnew;
}
// if (_dirs.size() > 1)printf("adding column with %d nodes\n",_column.size());
@@ -381,7 +377,7 @@ BoundaryLayerColumns* buidAdditionalPoints2D (GFace *gf)
MVertex *v = *it;
for (int i=0;i<_columns->getNbColumns(v);i++){
const BoundaryLayerData &data = _columns->getColumn(v,i);
- for (int j=0;j<data._column.size();j++){
+ for (unsigned int j = 0; j < data._column.size(); j++){
MVertex *blv = data._column[j];
fprintf(f,"SP(%g,%g,%g){%d};\n",blv->x(),blv->y(),blv->z(),v->getNum());
}
diff --git a/Mesh/meshGRegionDelaunayInsertion.cpp b/Mesh/meshGRegionDelaunayInsertion.cpp
index 8aef968963..5aedc64369 100644
--- a/Mesh/meshGRegionDelaunayInsertion.cpp
+++ b/Mesh/meshGRegionDelaunayInsertion.cpp
@@ -22,7 +22,6 @@
int MTet4::radiusNorm = 2;
static double LIMIT_ = 1;
-
static bool isActive(MTet4 *t, double limit_, int &active)
{
if (t->isDeleted()) return false;
@@ -35,23 +34,21 @@ static bool isActive(MTet4 *t, double limit_, int &active)
return false;
}
-
-
int MTet4::inCircumSphere(const double *p) const
{
- double pa[3] = {base->getVertex(0)->x(),
- base->getVertex(0)->y(),
+ double pa[3] = {base->getVertex(0)->x(),
+ base->getVertex(0)->y(),
base->getVertex(0)->z()};
- double pb[3] = {base->getVertex(1)->x(),
+ double pb[3] = {base->getVertex(1)->x(),
base->getVertex(1)->y(),
base->getVertex(1)->z()};
double pc[3] = {base->getVertex(2)->x(),
- base->getVertex(2)->y(),
+ base->getVertex(2)->y(),
base->getVertex(2)->z()};
double pd[3] = {base->getVertex(3)->x(),
base->getVertex(3)->y(),
base->getVertex(3)->z()};
- double result = robustPredicates::insphere(pa, pb, pc, pd, (double*)p) *
+ double result = robustPredicates::insphere(pa, pb, pc, pd, (double*)p) *
robustPredicates::orient3d(pa, pb, pc, pd);
return (result > 0) ? 1 : 0;
}
@@ -130,12 +127,12 @@ static bool isActive(MTet4 *t, double limit_, int &i, std::set<MFace,Less_Face>
return false;
}
-void connectTets(std::list<MTet4*> &l) { connectTets(l.begin(), l.end()); }
+void connectTets(std::list<MTet4*> &l) { connectTets(l.begin(), l.end()); }
void connectTets(std::vector<MTet4*> &l) { connectTets(l.begin(), l.end()); }
-void recurFindCavity(std::list<faceXtet> & shell,
- std::list<MTet4*> & cavity,
- MVertex *v ,
+void recurFindCavity(std::list<faceXtet> & shell,
+ std::list<MTet4*> & cavity,
+ MVertex *v ,
MTet4 *t)
{
// Msg::Info("tet %d %d %d %d",t->tet()->getVertex(0)->getNum(),
@@ -144,7 +141,7 @@ void recurFindCavity(std::list<faceXtet> & shell,
// t->tet()->getVertex(3)->getNum());
// invariant : this one has to be inserted in the cavity
- // consider this tet deleted
+ // consider this tet deleted
// remove its reference to its neighbors
t->setDeleted(true);
// the cavity that has to be removed
@@ -165,7 +162,7 @@ void recurFindCavity(std::list<faceXtet> & shell,
}
}
-bool insertVertex(MVertex *v,
+bool insertVertex(MVertex *v,
MTet4 *t,
MTet4Factory &myFactory,
std::set<MTet4*,compareTet4Ptr> &allTets,
@@ -174,10 +171,10 @@ bool insertVertex(MVertex *v,
std::set<MTet4*,compareTet4Ptr> *activeTets = 0 )
{
std::list<faceXtet> shell;
- std::list<MTet4*> cavity;
+ std::list<MTet4*> cavity;
std::list<MTet4*> new_cavity;
- recurFindCavity(shell, cavity, v, t);
+ recurFindCavity(shell, cavity, v, t);
// check that volume is conserved
double newVolume = 0;
@@ -208,8 +205,8 @@ bool insertVertex(MVertex *v,
vSizesBGM[tr->getVertex(2)->getIndex()] +
vSizesBGM[tr->getVertex(3)->getIndex()]);
double LL = std::min(lc, lcBGM);
-
- MTet4 *t4 = myFactory.Create(tr, vSizes, vSizesBGM);
+
+ MTet4 *t4 = myFactory.Create(tr, vSizes, vSizesBGM);
t4->setOnWhat(t->onWhat());
double d1 = sqrt((it->v[0]->x() - v->x()) * (it->v[0]->x() - v->x()) +
@@ -230,7 +227,7 @@ bool insertVertex(MVertex *v,
// all faces and ensure that the tets are all oriented the same.
new_cavity.push_back(t4);
MTet4 *otherSide = it->t1->getNeigh(it->i1);
-
+
if (otherSide)
new_cavity.push_back(otherSide);
// if (!it->t1->isDeleted())throw;
@@ -239,8 +236,8 @@ bool insertVertex(MVertex *v,
}
// OK, the cavity is star shaped
if (fabs(oldVolume - newVolume) < 1.e-10 * oldVolume &&
- !onePointIsTooClose){
- connectTets(new_cavity.begin(), new_cavity.end());
+ !onePointIsTooClose){
+ connectTets(new_cavity.begin(), new_cavity.end());
allTets.insert(newTets, newTets + shell.size());
if (activeTets){
@@ -260,9 +257,9 @@ bool insertVertex(MVertex *v,
else { // The cavity is NOT star shaped
// printf("hola %d %g %g %22.15E\n",onePointIsTooClose, oldVolume, newVolume, 100.*fabs(oldVolume - newVolume) / oldVolume);
for (unsigned int i = 0; i <shell.size(); i++) myFactory.Free(newTets[i]);
- delete [] newTets;
+ delete [] newTets;
ittet = cavity.begin();
- ittete = cavity.end();
+ ittete = cavity.end();
while(ittet != ittete){
(*ittet)->setDeleted(false);
++ittet;
@@ -281,8 +278,8 @@ static void setLcs(MTetrahedron *t, std::map<MVertex*, double> &vSizes)
double dy = vi->y()-vj->y();
double dz = vi->z()-vj->z();
double l = sqrt(dx * dx + dy * dy + dz * dz);
- std::map<MVertex*,double>::iterator iti = vSizes.find(vi);
- std::map<MVertex*,double>::iterator itj = vSizes.find(vj);
+ std::map<MVertex*,double>::iterator iti = vSizes.find(vi);
+ std::map<MVertex*,double>::iterator itj = vSizes.find(vj);
if (iti == vSizes.end() || iti->second > l) vSizes[vi] = l;
if (itj == vSizes.end() || itj->second > l) vSizes[vj] = l;
}
@@ -301,13 +298,13 @@ static void setLcs(MTetrahedron *t, std::map<MVertex*, double> &vSizes)
// 4th argument will disappear when the reclassification of vertices will be done
bool find_triangle_in_model(GModel *model, MTriangle *tri, GFace **gfound, bool force)
-{
+{
static compareMTriangleLexicographic cmp;
-
+
GModel::fiter fit = model->firstFace() ;
while (fit != model->lastFace()){
- bool found = std::binary_search((*fit)->triangles.begin(),
- (*fit)->triangles.end(),
+ bool found = std::binary_search((*fit)->triangles.begin(),
+ (*fit)->triangles.end(),
tri, cmp);
if(found){
*gfound = *fit;
@@ -331,23 +328,23 @@ GRegion *getRegionFromBoundingFaces(GModel *model,
std::list <GFace *> _faces = (*git)->faces();
// printf("region %d %d faces\n",(*git)->tag(),_faces.size());
// for (std::list<GFace *>::iterator it = _faces.begin(); it != _faces.end(); ++it){
- // printf(" %d",(*it)->tag());
+ // printf(" %d",(*it)->tag());
// }
// printf("\n");
-
+
if(_faces.size() == faces_bound.size()){
bool ok = true;
for (std::list<GFace *>::iterator it = _faces.begin(); it != _faces.end(); ++it){
if(faces_bound.find (*it) == faces_bound.end()) ok = false;
}
if (ok) return *git;
- }
+ }
++git;
}
return 0;
}
-
-void recur_classify(MTet4 *t, std::list<MTet4*> &theRegion,
+
+void recur_classify(MTet4 *t, std::list<MTet4*> &theRegion,
std::set<GFace*> &faces_bound, GRegion *bidon,
GModel *model, const fs_cont &search)
{
@@ -355,7 +352,7 @@ void recur_classify(MTet4 *t, std::list<MTet4*> &theRegion,
if (t->onWhat()) return; // should never return here...
theRegion.push_back(t);
t->setOnWhat(bidon);
-
+
bool FF[4] = {0,0,0,0};
for (int i = 0; i < 4; i++){
@@ -370,14 +367,14 @@ void recur_classify(MTet4 *t, std::list<MTet4*> &theRegion,
if (faces_bound.find(gfound) == faces_bound.end())
faces_bound.insert(gfound);
}
-
+
// MTriangle tri (t->tet()->getVertex (faces[i][0]),
// t->tet()->getVertex (faces[i][1]),
// t->tet()->getVertex (faces[i][2]));
// GFace *gfound;
// if (FF[i] = find_triangle_in_model(model, &tri, &gfound, false)){
// if (faces_bound.find(gfound) == faces_bound.end())
-// faces_bound.insert(gfound);
+// faces_bound.insert(gfound);
// }
}
}
@@ -397,7 +394,7 @@ void adaptMeshGRegion::operator () (GRegion *gr)
allTets.push_back(new MTet4(gr->tetrahedra[i], qm));
}
gr->tetrahedra.clear();
-
+
connectTets(allTets.begin(),allTets.end());
double t1 = Cpu();
@@ -422,7 +419,7 @@ void adaptMeshGRegion::operator () (GRegion *gr)
double low = (double)i / nbRanges;
double high = (double)(i + 1) / nbRanges;
if (qual >= low && qual < high) quality_ranges[i]++;
- }
+ }
}
}
Msg::Info("Adaptation : START with %12.5E QBAD %12.5E QAVG %12.5E",
@@ -432,16 +429,16 @@ void adaptMeshGRegion::operator () (GRegion *gr)
double high = (double)(i + 1) / nbRanges;
Msg::Info("Opti : %3.2f < QUAL < %3.2f : %9d elements ",
low, high, quality_ranges[i]);
- }
- }
-
+ }
+ }
+
double qMin = 0.5;
double sliverLimit = 0.2;
int nbESwap = 0, nbFSwap = 0, nbReloc = 0, nbCollapse = 0;
-
+
while (1){
- std::vector<MTet4*> newTets;
+ std::vector<MTet4*> newTets;
for (CONTAINER::iterator it = allTets.begin(); it!=allTets.end(); ++it){
if (!(*it)->isDeleted()){
for (int i = 0; i < 4; i++){
@@ -453,7 +450,7 @@ void adaptMeshGRegion::operator () (GRegion *gr)
}
}
}
-
+
printf("nbCollapses = %d\n", nbCollapse);
for (CONTAINER::iterator it = allTets.begin(); it!=allTets.end(); ++it){
@@ -469,7 +466,7 @@ void adaptMeshGRegion::operator () (GRegion *gr)
}
}
}
-
+
illegals.clear();
for (int i = 0; i < nbRanges; i++) quality_ranges[i] = 0;
@@ -480,7 +477,7 @@ void adaptMeshGRegion::operator () (GRegion *gr)
for (int i = 0; i < 6; i++){
if (edgeSwap(newTets, *it, i, qm)) {
nbESwap++;
- break;
+ break;
}
}
if (!(*it)->isDeleted()){
@@ -489,13 +486,13 @@ void adaptMeshGRegion::operator () (GRegion *gr)
double low = (double)i / nbRanges;
double high = (double)(i + 1)/ nbRanges;
if (qq >= low && qq < high) quality_ranges[i]++;
- }
+ }
}
}
}
-
+
if (!newTets.size()) break;
-
+
// add all the new tets in the container
for(unsigned int i = 0; i < newTets.size(); i++){
if(!newTets[i]->isDeleted()){
@@ -505,7 +502,7 @@ void adaptMeshGRegion::operator () (GRegion *gr)
delete newTets[i]->tet();
delete newTets[i];
}
- }
+ }
// relocate vertices
for (CONTAINER::iterator it = allTets.begin(); it != allTets.end(); ++it){
@@ -517,7 +514,7 @@ void adaptMeshGRegion::operator () (GRegion *gr)
}
}
}
-
+
double totalVolumeb = 0.0;
double worst = 1.0;
double avg = 0;
@@ -537,11 +534,11 @@ void adaptMeshGRegion::operator () (GRegion *gr)
nbESwap, nbFSwap, nbReloc, totalVolumeb, worst, avg / count, t2 - t1);
break;
}
-
+
int nbSlivers = 0;
for(unsigned int i = 0; i < illegals.size(); i++)
if(!(illegals[i]->isDeleted())) nbSlivers++;
-
+
if (nbSlivers){
Msg::Info("Opti : %d illegal tets are still in the mesh, trying to remove them",
nbSlivers);
@@ -555,8 +552,8 @@ void adaptMeshGRegion::operator () (GRegion *gr)
double high = (double)(i + 1) / nbRanges;
Msg::Info("Opti : %3.2f < QUAL < %3.2f : %9d elements",
low, high, quality_ranges[i]);
- }
-
+ }
+
for (CONTAINER::iterator it = allTets.begin(); it != allTets.end(); ++it){
if (!(*it)->isDeleted()){
gr->tetrahedra.push_back((*it)->tet());
@@ -564,12 +561,12 @@ void adaptMeshGRegion::operator () (GRegion *gr)
}
else{
delete (*it)->tet();
- delete *it;
+ delete *it;
}
}
}
-//template <class CONTAINER, class DATA>
+//template <class CONTAINER, class DATA>
void optimizeMesh(GRegion *gr, const qualityMeasure4Tet &qm)
{
typedef std::list<MTet4 *> CONTAINER ;
@@ -580,7 +577,7 @@ void optimizeMesh(GRegion *gr, const qualityMeasure4Tet &qm)
allTets.push_back(t);
}
gr->tetrahedra.clear();
-
+
connectTets(allTets.begin(), allTets.end());
double t1 = Cpu();
@@ -605,7 +602,7 @@ void optimizeMesh(GRegion *gr, const qualityMeasure4Tet &qm)
double low = (double)i / nbRanges;
double high = (double)(i + 1) / nbRanges;
if (qual >= low && qual < high) quality_ranges[i]++;
- }
+ }
}
}
Msg::Info("Opti : START with %12.5E QBAD %12.5E QAVG %12.5E",
@@ -615,16 +612,16 @@ void optimizeMesh(GRegion *gr, const qualityMeasure4Tet &qm)
double high = (double)(i + 1) / nbRanges;
Msg::Info("Opti : %3.2f < QUAL < %3.2f : %9d elements",
low, high, quality_ranges[i]);
- }
- }
-
+ }
+ }
+
double qMin = 0.5;
double sliverLimit = 0.1;
int nbESwap = 0, nbFSwap = 0, nbReloc = 0;
-
+
while (1){
- std::vector<MTet4*> newTets;
+ std::vector<MTet4*> newTets;
for (CONTAINER::iterator it = allTets.begin(); it != allTets.end(); ++it){
if (!(*it)->isDeleted()){
double qq = (*it)->getQuality();
@@ -638,7 +635,7 @@ void optimizeMesh(GRegion *gr, const qualityMeasure4Tet &qm)
}
}
}
-
+
illegals.clear();
for (int i = 0; i < nbRanges; i++) quality_ranges[i] = 0;
@@ -658,11 +655,11 @@ void optimizeMesh(GRegion *gr, const qualityMeasure4Tet &qm)
double low = (double)i / nbRanges;
double high = (double)(i + 1) / nbRanges;
if (qq >= low && qq < high) quality_ranges[i]++;
- }
+ }
}
}
}
-
+
if (0 && !newTets.size()){
int nbSlivers = 0;
int nbSliversWeCanDoSomething = 0;
@@ -674,11 +671,11 @@ void optimizeMesh(GRegion *gr, const qualityMeasure4Tet &qm)
}
Msg::Info("Opti : %d Sliver Removals", nbSliversWeCanDoSomething);
}
-
+
if (!newTets.size()){
break;
}
-
+
// add all the new tets in the container
for(unsigned int i = 0; i < newTets.size(); i++){
if(!newTets[i]->isDeleted()){
@@ -688,7 +685,7 @@ void optimizeMesh(GRegion *gr, const qualityMeasure4Tet &qm)
delete newTets[i]->tet();
delete newTets[i];
}
- }
+ }
// relocate vertices
for (CONTAINER::iterator it = allTets.begin();it!=allTets.end();++it){
@@ -719,7 +716,7 @@ void optimizeMesh(GRegion *gr, const qualityMeasure4Tet &qm)
Msg::Info("Opti : (%d,%d,%d) = %12.5E QBAD %12.5E QAVG %12.5E (%8.3f sec)",
nbESwap, nbFSwap, nbReloc, totalVolumeb, worst, avg / count, t2 - t1);
}
-
+
if (illegals.size()){
Msg::Info("Opti : %d illegal tets are still in the mesh", illegals.size());
}
@@ -732,7 +729,7 @@ void optimizeMesh(GRegion *gr, const qualityMeasure4Tet &qm)
double high = (double)(i + 1) / nbRanges;
Msg::Info("Opti : %3.2f < QUAL < %3.2f : %9d elements",
low, high, quality_ranges[i]);
- }
+ }
for (CONTAINER::iterator it = allTets.begin(); it != allTets.end(); ++it){
if (!(*it)->isDeleted()){
@@ -741,13 +738,13 @@ void optimizeMesh(GRegion *gr, const qualityMeasure4Tet &qm)
}
else{
delete (*it)->tet();
- delete *it;
+ delete *it;
}
}
}
-static double tetcircumcenter(double a[3], double b[3], double c[3], double d[3],
+static double tetcircumcenter(double a[3], double b[3], double c[3], double d[3],
double circumcenter[3], double *xi, double *eta, double *zeta){
double xba, yba, zba, xca, yca, zca, xda, yda, zda;
double balength, calength, dalength;
@@ -756,7 +753,7 @@ static double tetcircumcenter(double a[3], double b[3], double c[3], double d[3]
double xcrossbc, ycrossbc, zcrossbc;
double denominator;
double xcirca, ycirca, zcirca;
-
+
/* Use coordinates relative to point `a' of the tetrahedron. */
xba = b[0] - a[0];
yba = b[1] - a[1];
@@ -781,14 +778,14 @@ static double tetcircumcenter(double a[3], double b[3], double c[3], double d[3]
xcrossbc = yba * zca - yca * zba;
ycrossbc = zba * xca - zca * xba;
zcrossbc = xba * yca - xca * yba;
-
+
/* Calculate the denominator of the formulae. */
/* Use orient3d() from http://www.cs.cmu.edu/~quake/robust.html */
/* to ensure a correctly signed (and reasonably accurate) result, */
/* avoiding any possibility of division by zero. */
const double xxx = robustPredicates::orient3d(b, c, d, a);
denominator = 0.5 / xxx;
-
+
/* Calculate offset (from `a') of circumcenter. */
xcirca = (balength * xcrosscd + calength * xcrossdb + dalength * xcrossbc) *
denominator;
@@ -799,15 +796,15 @@ static double tetcircumcenter(double a[3], double b[3], double c[3], double d[3]
circumcenter[0] = xcirca + a[0];
circumcenter[1] = ycirca + a[1];
circumcenter[2] = zcirca + a[2];
-
- /*
+
+ /*
printf(" %g %g %g %g\n",
sqrt((a[0]-xcirca)*(a[0]-xcirca)+(a[1]-ycirca)*(a[1]-ycirca)+(a[2]-zcirca)*(a[2]-zcirca)),
sqrt((b[0]-xcirca)*(b[0]-xcirca)+(b[1]-ycirca)*(b[1]-ycirca)+(b[2]-zcirca)*(b[2]-zcirca)),
sqrt((c[0]-xcirca)*(c[0]-xcirca)+(c[1]-ycirca)*(c[1]-ycirca)+(c[2]-zcirca)*(c[2]-zcirca)),
sqrt((d[0]-xcirca)*(d[0]-xcirca)+(d[1]-ycirca)*(d[1]-ycirca)+(d[2]-zcirca)*(d[2]-zcirca)) );
*/
-
+
if (xi != (double *) NULL) {
/* To interpolate a linear function at the circumcenter, define a */
/* coordinate system with a xi-axis directed from `a' to `b', */
@@ -825,7 +822,7 @@ static double tetcircumcenter(double a[3], double b[3], double c[3], double d[3]
}
-void insertVerticesInRegion (GRegion *gr)
+void insertVerticesInRegion (GRegion *gr)
{
//printf("sizeof MTet4 = %d sizeof MTetrahedron %d sizeof(MVertex) %d\n",
// sizeof(MTet4), sizeof(MTetrahedron), sizeof(MVertex));
@@ -841,14 +838,14 @@ void insertVerticesInRegion (GRegion *gr)
std::map<MVertex*, double> vSizesMap;
for(unsigned int i = 0; i < gr->tetrahedra.size(); i++)
setLcs(gr->tetrahedra[i], vSizesMap);
- for(std::map<MVertex*, double>::iterator it = vSizesMap.begin();
+ for(std::map<MVertex*, double>::iterator it = vSizesMap.begin();
it != vSizesMap.end(); ++it){
it->first->setIndex(NUM++);
vSizes.push_back(it->second);
vSizesBGM.push_back(it->second);
}
}
-
+
for(unsigned int i = 0; i < gr->tetrahedra.size(); i++)
allTets.insert(myFactory.Create(gr->tetrahedra[i], vSizes,vSizesBGM));
@@ -867,7 +864,7 @@ void insertVerticesInRegion (GRegion *gr)
std::set<GFace *> faces_bound;
GRegion *bidon = (GRegion*)123;
double _t1 = Cpu();
- Msg::Debug("start with a non classified tet");
+ Msg::Debug("start with a non classified tet");
recur_classify(*it, theRegion, faces_bound, bidon, gr->model(), search);
double _t2 = Cpu();
Msg::Debug("found %d tets with %d faces (%g sec for the classification)",
@@ -875,15 +872,15 @@ void insertVerticesInRegion (GRegion *gr)
GRegion *myGRegion = getRegionFromBoundingFaces(gr->model(), faces_bound);
Msg::Info("a region is found %p",myGRegion);
if(myGRegion) // a geometrical region associated to the list of faces has been found
- for(std::list<MTet4*>::iterator it2 = theRegion.begin();
+ for(std::list<MTet4*>::iterator it2 = theRegion.begin();
it2 != theRegion.end(); ++it2) (*it2)->setOnWhat(myGRegion);
- else // the tets are in the void
+ else // the tets are in the void
for(std::list<MTet4*>::iterator it2 = theRegion.begin();
it2 != theRegion.end(); ++it2)(*it2)->setDeleted(true);
}
}
search.clear();
-
+
for(MTet4Factory::iterator it = allTets.begin(); it!=allTets.end(); ++it){
(*it)->setNeigh(0, 0);
(*it)->setNeigh(1, 0);
@@ -902,7 +899,7 @@ void insertVerticesInRegion (GRegion *gr)
Msg::Error("No tetrahedra in region %d %d", gr->tag(), allTets.size());
break;
}
-
+
MTet4 *worst = *allTets.begin();
if(worst->isDeleted()){
@@ -917,26 +914,26 @@ void insertVerticesInRegion (GRegion *gr)
double center[3];
double uvw[3];
MTetrahedron *base = worst->tet();
- double pa[3] = {base->getVertex(0)->x(),
- base->getVertex(0)->y(),
+ double pa[3] = {base->getVertex(0)->x(),
+ base->getVertex(0)->y(),
base->getVertex(0)->z()};
- double pb[3] = {base->getVertex(1)->x(),
+ double pb[3] = {base->getVertex(1)->x(),
base->getVertex(1)->y(),
base->getVertex(1)->z()};
double pc[3] = {base->getVertex(2)->x(),
- base->getVertex(2)->y(),
+ base->getVertex(2)->y(),
base->getVertex(2)->z()};
double pd[3] = {base->getVertex(3)->x(),
base->getVertex(3)->y(),
base->getVertex(3)->z()};
-
+
// double UU,VV,WW;
tetcircumcenter(pa,pb,pc,pd, center,&uvw[0],&uvw[1],&uvw[2] );
// worst->tet()->xyz2uvw(center, uvw);
// printf("%12.5E %12.5E %12.5E -- %12.5E %12.5E %12.5E \n",
// UU,VV,WW,uvw[0],uvw[1],uvw[2]);
/*
- if (!worst->tet()->isInside(uvw[0], uvw[1], uvw[2]) &&
+ if (!worst->tet()->isInside(uvw[0], uvw[1], uvw[2]) &&
worst->getRadius() > 20){
uvw[0] = uvw[1] = uvw[2] = 0.25;
center[0] = 0.25*(pa[0]+pb[0]+pc[0]+pd[0]);
@@ -945,10 +942,10 @@ void insertVerticesInRegion (GRegion *gr)
}
*/
- if(worst->tet()->isInside(uvw[0], uvw[1], uvw[2])){
+ if(worst->tet()->isInside(uvw[0], uvw[1], uvw[2])){
MVertex *v = new MVertex(center[0], center[1], center[2], worst->onWhat());
v->setIndex(NUM++);
- double lc1 =
+ double lc1 =
(1 - uvw[0] - uvw[1] - uvw[2]) * vSizes[worst->tet()->getVertex(0)->getIndex()] +
uvw[0] * vSizes[worst->tet()->getVertex(1)->getIndex()] +
uvw[1] * vSizes[worst->tet()->getVertex(2)->getIndex()] +
@@ -1003,7 +1000,7 @@ void insertVerticesInRegion (GRegion *gr)
}
}
-
+
while(1){
if(allTets.begin() == allTets.end()) break;
MTet4 *worst = *allTets.begin();
@@ -1012,50 +1009,32 @@ void insertVerticesInRegion (GRegion *gr)
worst->tet() = 0;
}
myFactory.Free(worst);
- allTets.erase(allTets.begin());
+ allTets.erase(allTets.begin());
}
}
-MVertex * optimalPointFrontal (GRegion *gr,
- MTet4 *worst,
- int active_face,
- std::vector<double> &vSizes,
- std::vector<double> &vSizesBGM){
+MVertex * optimalPointFrontal(GRegion *gr,
+ MTet4 *worst,
+ int active_face,
+ std::vector<double> &vSizes,
+ std::vector<double> &vSizesBGM)
+{
double centerTet[3], centerFace[3];
- MTetrahedron *base = worst->tet();
faceXtet fxt ( worst, active_face );
double pa[3] = {fxt.v[0]->x(),fxt.v[0]->y(),fxt.v[0]->z()};
double pb[3] = {fxt.v[1]->x(),fxt.v[1]->y(),fxt.v[1]->z()};
double pc[3] = {fxt.v[2]->x(),fxt.v[2]->y(),fxt.v[2]->z()};
circumCenterXYZ(pa, pb, pc, centerFace);
worst->circumcenter(centerTet);
-
+
SVector3 dir (centerTet[0] - centerFace[0],
centerTet[1] - centerFace[1],
centerTet[2] - centerFace[2]);
-
- const double q = dir.norm();
dir.normalize();
SVector3 rDir (pa[0] - centerFace[0],
pa[1] - centerFace[1],
pa[2] - centerFace[2]);
-
- const double p = 0.5 * rDir.norm();
-
- const double rhoM1 = 0.33333 * (vSizes[fxt.v[0]->getIndex()] +
- vSizes[fxt.v[1]->getIndex()] +
- vSizes[fxt.v[2]->getIndex()] );
- const double rhoM2 = 0.33333 * (vSizesBGM[fxt.v[0]->getIndex()] +
- vSizesBGM[fxt.v[1]->getIndex()] +
- vSizesBGM[fxt.v[2]->getIndex()] );
-
- const double rhoM = std::min(rhoM1, rhoM2);
-
- const double HEIGHT = 1/sqrt(3.);
-
- const double rhoM_hat = std::min(std::max(rhoM, p), (p * p + q * q) / (2 * q));
- const double d = (rhoM_hat + sqrt (rhoM_hat * rhoM_hat - p * p)) * HEIGHT;
// const double a = 2*p *3 /sqrt(3);
const double a = .025;
@@ -1071,7 +1050,6 @@ MVertex * optimalPointFrontal (GRegion *gr,
void bowyerWatsonFrontalLayers(GRegion *gr, bool hex)
{
-
std::vector<double> vSizes;
std::vector<double> vSizesBGM;
MTet4Factory myFactory(1600000);
@@ -1092,14 +1070,14 @@ void bowyerWatsonFrontalLayers(GRegion *gr, bool hex)
std::map<MVertex*, double> vSizesMap;
for(unsigned int i = 0; i < gr->tetrahedra.size(); i++)
setLcs(gr->tetrahedra[i], vSizesMap);
- for(std::map<MVertex*, double>::iterator it = vSizesMap.begin();
+ for(std::map<MVertex*, double>::iterator it = vSizesMap.begin();
it != vSizesMap.end(); ++it){
it->first->setIndex(NUM++);
vSizes.push_back(it->second);
vSizesBGM.push_back(it->second);
}
}
-
+
for(unsigned int i = 0; i < gr->tetrahedra.size(); i++)
allTets.insert(myFactory.Create(gr->tetrahedra[i], vSizes,vSizesBGM));
@@ -1108,14 +1086,14 @@ void bowyerWatsonFrontalLayers(GRegion *gr, bool hex)
fs_cont search;
buildFaceSearchStructure(gr->model(), search);
-
+
for(MTet4Factory::iterator it = allTets.begin(); it != allTets.end(); ++it){
if(!(*it)->onWhat()){
std::list<MTet4*> theRegion;
std::set<GFace *> faces_bound;
GRegion *bidon = (GRegion*)123;
double _t1 = Cpu();
- Msg::Debug("start with a non classified tet");
+ Msg::Debug("start with a non classified tet");
recur_classify(*it, theRegion, faces_bound, bidon, gr->model(), search);
double _t2 = Cpu();
Msg::Debug("found %d tets with %d faces (%g sec for the classification)",
@@ -1123,9 +1101,9 @@ void bowyerWatsonFrontalLayers(GRegion *gr, bool hex)
GRegion *myGRegion = getRegionFromBoundingFaces(gr->model(), faces_bound);
// Msg::Info("a region is found %p",myGRegion);
if(myGRegion) // a geometrical region associated to the list of faces has been found
- for(std::list<MTet4*>::iterator it2 = theRegion.begin();
+ for(std::list<MTet4*>::iterator it2 = theRegion.begin();
it2 != theRegion.end(); ++it2) (*it2)->setOnWhat(myGRegion);
- else // the tets are in the void
+ else // the tets are in the void
for(std::list<MTet4*>::iterator it2 = theRegion.begin();
it2 != theRegion.end(); ++it2)(*it2)->setDeleted(true);
}
@@ -1160,12 +1138,12 @@ void bowyerWatsonFrontalLayers(GRegion *gr, bool hex)
std::set<MTet4*, compareTet4Ptr> activeTetsNotInFront;
while (1){
-
+
if (!activeTets.size())break;
// printf("%d active tets %d tets\n",activeTets.size(),allTets.size());
- std::set<MTet4*,compareTet4Ptr>::iterator WORST_ITER = activeTets.begin();
+ std::set<MTet4*,compareTet4Ptr>::iterator WORST_ITER = activeTets.begin();
MTet4 *worst = (*WORST_ITER);
if(worst->isDeleted()){
activeTets.erase(WORST_ITER);
@@ -1173,27 +1151,27 @@ void bowyerWatsonFrontalLayers(GRegion *gr, bool hex)
}
else {
activeTets.erase(WORST_ITER);
- if (isActive(worst, LIMIT_, active_face,&_front) &&
+ if (isActive(worst, LIMIT_, active_face,&_front) &&
worst->getRadius() > LIMIT_){
// printf("worst = %12.5E\n",worst->getRadius());
-
+
if(ITER++ % 5000 == 0)
Msg::Info("%d points created -- Worst tet radius is %g",
vSizes.size(), worst->getRadius());
-
+
MVertex *v = optimalPointFrontal (gr,worst,active_face,vSizes,vSizesBGM);
v->setIndex(NUM++);
vSizes.push_back(.025);
vSizesBGM.push_back(.025);
if(!worst->inCircumSphere(v) ||
- !insertVertex(v, worst, myFactory, allTets, vSizes,vSizesBGM,&activeTets)){
+ !insertVertex(v, worst, myFactory, allTets, vSizes,vSizesBGM,&activeTets)){
myFactory.changeTetRadius(allTets.begin(), 0.);
if(v) delete v;
}
else{
// printf("yeah ! one new vertex \n");
- v->onWhat()->mesh_vertices.push_back(v);
+ v->onWhat()->mesh_vertices.push_back(v);
// if (ITER == 100)break;
}
}
@@ -1236,8 +1214,8 @@ void bowyerWatsonFrontalLayers(GRegion *gr, bool hex)
worst->tet() = 0;
}
myFactory.Free(worst);
- allTets.erase(allTets.begin());
+ allTets.erase(allTets.begin());
}
MTet4::radiusNorm = 2;
LIMIT_ = 1;
-}
+}
diff --git a/Mesh/meshMetric.cpp b/Mesh/meshMetric.cpp
index cb61ef8ab9..d05d177a4f 100644
--- a/Mesh/meshMetric.cpp
+++ b/Mesh/meshMetric.cpp
@@ -12,7 +12,9 @@
#include "MElementOctree.h"
#include "OS.h"
-static void increaseStencil(MVertex *v, v2t_cont &adj, std::vector<MElement*> <){
+/*
+static void increaseStencil(MVertex *v, v2t_cont &adj, std::vector<MElement*> <)
+{
std::set<MElement*> stencil;
std::set<MVertex*> vs;
stencil.insert(lt.begin(),lt.end());
@@ -28,8 +30,10 @@ static void increaseStencil(MVertex *v, v2t_cont &adj, std::vector<MElement*> &l
lt.clear();
lt.insert(lt.begin(),stencil.begin(),stencil.end());
}
+*/
-meshMetric::meshMetric(GModel *gm){
+meshMetric::meshMetric(GModel *gm)
+{
_dim = gm->getDim();
std::map<MElement*, MElement*> newP;
std::map<MElement*, MElement*> newD;
@@ -57,7 +61,8 @@ meshMetric::meshMetric(GModel *gm){
}
-meshMetric::meshMetric(std::vector<MElement*> elements){
+meshMetric::meshMetric(std::vector<MElement*> elements)
+{
_dim = elements[0]->getDim();
std::map<MElement*, MElement*> newP;
std::map<MElement*, MElement*> newD;
@@ -71,7 +76,9 @@ meshMetric::meshMetric(std::vector<MElement*> elements){
_octree = new MElementOctree(_elements);
}
-void meshMetric::addMetric(int technique, simpleFunction<double> *fct, std::vector<double> parameters){
+void meshMetric::addMetric(int technique, simpleFunction<double> *fct,
+ std::vector<double> parameters)
+{
needMetricUpdate=true;
_fct = fct;
hmin = parameters.size() >=3 ? parameters[1] : CTX::instance()->mesh.lcMin;
@@ -87,8 +94,8 @@ void meshMetric::addMetric(int technique, simpleFunction<double> *fct, std::vect
computeMetric();
}
-void meshMetric::intersectMetrics(){
-
+void meshMetric::intersectMetrics()
+{
if (!setOfMetrics.size()){
std::cout << " meshMetric::intersectMetrics: Can't intersect metrics, no metric registered ! " << std::endl;
return;
@@ -110,8 +117,8 @@ void meshMetric::intersectMetrics(){
}
-void meshMetric::exportInfo(const char * fileendname){
-
+void meshMetric::exportInfo(const char * fileendname)
+{
if (needMetricUpdate) intersectMetrics();
std::stringstream sg,sm,sl,sh;
sg << "meshmetric_gradients_" << fileendname;
@@ -188,18 +195,17 @@ void meshMetric::exportInfo(const char * fileendname){
out_metric.close();
out_ls.close();
out_hess.close();
-
}
-
-meshMetric::~meshMetric(){
+meshMetric::~meshMetric()
+{
if (_octree) delete _octree;
for (unsigned int i=0; i< _elements.size(); i++)
delete _elements[i];
}
-void meshMetric::computeValues(){
-
+void meshMetric::computeValues()
+{
v2t_cont :: iterator it = _adj.begin();
while (it != _adj.end()) {
std::vector<MElement*> lt = it->second;
@@ -209,10 +215,9 @@ void meshMetric::computeValues(){
}
}
-
-
// Determines set of vertices to use for least squares
-std::vector<MVertex*> getLSBlob(unsigned int minNbPt, v2t_cont::iterator it, v2t_cont &adj) {
+std::vector<MVertex*> getLSBlob(unsigned int minNbPt, v2t_cont::iterator it, v2t_cont &adj)
+{
// static const double RADFACT = 3;
//
@@ -248,13 +253,11 @@ std::vector<MVertex*> getLSBlob(unsigned int minNbPt, v2t_cont::iterator it, v2t
}
-
-
// Compute derivatives and second order derivatives using least squares
// 2D LS system: a_i0*x^2+a_i1*x*y+a_i2*y^2+a_i3*x+a_i4*y+a_i5=b_i
// 3D LS system: a_i0*x^2+a_i1*x*y+a_i2*x*z+a_i3*y^2+a_i4*y*z+a_i5*z^2+a_i6*x+a_i7*y+a_i8*z+a_i9=b_i
-void meshMetric::computeHessian(){
-
+void meshMetric::computeHessian()
+{
unsigned int sysDim = (_dim == 2) ? 6 : 10;
unsigned int minNbPtBlob = 3*sysDim;
@@ -301,13 +304,10 @@ void meshMetric::computeHessian(){
dgrads[1][ver] = SVector3(d2udxy,d2udy2,d2udyz);
dgrads[2][ver] = SVector3(d2udxz,d2udyz,d2udz2);
}
-
}
-
-
-void meshMetric::computeMetricLevelSet(){
-
+void meshMetric::computeMetricLevelSet()
+{
int metricNumber = setOfMetrics.size();
for (v2t_cont::iterator it=_adj.begin(); it!=_adj.end(); it++) {
@@ -362,10 +362,8 @@ void meshMetric::computeMetricLevelSet(){
}
-
-
-void meshMetric::computeMetricHessian(){
-
+void meshMetric::computeMetricHessian()
+{
int metricNumber = setOfMetrics.size();
for (v2t_cont::iterator it=_adj.begin(); it!=_adj.end(); it++) {
@@ -403,10 +401,8 @@ void meshMetric::computeMetricHessian(){
}
-
-
-void meshMetric::computeMetricFrey(){
-
+void meshMetric::computeMetricFrey()
+{
int metricNumber = setOfMetrics.size();
for (v2t_cont::iterator it=_adj.begin(); it!=_adj.end(); it++) {
@@ -466,13 +462,10 @@ void meshMetric::computeMetricFrey(){
setOfDetMetric[metricNumber].insert(std::make_pair(ver,sqrt(metric.determinant())));
}
-
}
-
-
-void meshMetric::computeMetricEigenDir(){
-
+void meshMetric::computeMetricEigenDir()
+{
int metricNumber = setOfMetrics.size();
for (v2t_cont::iterator it=_adj.begin(); it!=_adj.end(); it++) {
@@ -559,10 +552,8 @@ void meshMetric::computeMetricEigenDir(){
}
-
-
-void meshMetric::computeMetricIsoLinInterp(){
-
+void meshMetric::computeMetricIsoLinInterp()
+{
int metricNumber = setOfMetrics.size();
for (v2t_cont::iterator it=_adj.begin(); it!=_adj.end(); it++) {
@@ -589,10 +580,8 @@ void meshMetric::computeMetricIsoLinInterp(){
}
-
-
-void meshMetric::computeMetricScaledHessian(){
-
+void meshMetric::computeMetricScaledHessian()
+{
int metricNumber = setOfMetrics.size();
std::list<double> lambda1, lambda2, lambda3;
@@ -658,10 +647,8 @@ void meshMetric::computeMetricScaledHessian(){
}
-
-
-void meshMetric::computeMetric(){
-
+void meshMetric::computeMetric()
+{
//printf("%d elements are considered in the meshMetric \n",(int)_elements.size());
computeValues();
@@ -679,9 +666,8 @@ void meshMetric::computeMetric(){
}
-
-
-double meshMetric::operator() (double x, double y, double z, GEntity *ge) {
+double meshMetric::operator() (double x, double y, double z, GEntity *ge)
+{
if (needMetricUpdate) intersectMetrics();
if (!setOfMetrics.size()){
std::cout << "meshMetric::operator() : No metric defined ! " << std::endl;
@@ -705,7 +691,8 @@ double meshMetric::operator() (double x, double y, double z, GEntity *ge) {
return 1.e22;
}
-void meshMetric::operator() (double x, double y, double z, SMetric3 &metr, GEntity *ge) {
+void meshMetric::operator() (double x, double y, double z, SMetric3 &metr, GEntity *ge)
+{
if (needMetricUpdate) intersectMetrics();
if (!setOfMetrics.size()){
std::cout << "meshMetric::operator() : No metric defined ! " << std::endl;
diff --git a/Mesh/qualityMeasures.cpp b/Mesh/qualityMeasures.cpp
index 4a16e477f1..425fe59d19 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;
@@ -464,11 +464,11 @@ double qmTriangleAngles (MTriangle *e) {
double c;
prosca(v1,v2,&c);
double x = acos(c)-M_PI/3;
- double angle = (x+M_PI/3)/M_PI*180;
+ //double angle = (x+M_PI/3)/M_PI*180;
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);
+ // minAngle = std::min(angle, minAngle);
// printf("Angle %g ", angle);
// printf("Quality %g\n",quality);
}
--
GitLab