diff --git a/Mesh/meshGFaceBDS.cpp b/Mesh/meshGFaceBDS.cpp
index 1d7710a6e5cc50a63482f0812847e8a8427d103a..66f0079f3a51bd12a60fca141ae9d2bcf2a5493a 100644
--- a/Mesh/meshGFaceBDS.cpp
+++ b/Mesh/meshGFaceBDS.cpp
@@ -32,21 +32,21 @@ double computeEdgeLinearLength(BDS_Point *p1, BDS_Point *p2)
return l;
}
-extern double lengthInfniteNorm(const double p[2], const double q[2],
+extern double lengthInfniteNorm(const double p[2], const double q[2],
const double quadAngle);
inline double computeEdgeLinearLength(BDS_Point *p1, BDS_Point *p2, GFace *f,
double SCALINGU, double SCALINGV)
{
// FIXME SUPER HACK
-
+
// if (CTX::instance()->mesh.recombineAll || f->meshAttributes.recombine || 1){
// double quadAngle = backgroundMesh::current()->getAngle (0.5 * (p1->u + p2->u) * SCALINGU, 0.5 * (p1->v + p2->v) * SCALINGV,0);
// const double a [2] = {p1->u,p1->v};
// const double b [2] = {p2->u,p2->v};
// return lengthInfniteNorm(a,b, quadAngle);
// }
-
+
GPoint GP = f->point(SPoint2(0.5 * (p1->u + p2->u) * SCALINGU,
0.5 * (p1->v + p2->v) * SCALINGV));
@@ -153,7 +153,7 @@ inline double computeEdgeMiddleCoord(BDS_Point *p1, BDS_Point *p2, GFace *f,
return 0.25 * (3 * l2 - l1) / l2;
}
-inline double computeEdgeLinearLength(BDS_Edge *e, GFace *f,
+inline double computeEdgeLinearLength(BDS_Edge *e, GFace *f,
double SCALINGU, double SCALINGV)
{
// FIXME !!!
@@ -165,7 +165,7 @@ inline double computeEdgeLinearLength(BDS_Edge *e, GFace *f,
double NewGetLc(BDS_Point *p)
{
- return Extend1dMeshIn2dSurfaces() ?
+ return Extend1dMeshIn2dSurfaces() ?
std::min(p->lc(), p->lcBGM()) : p->lcBGM();
}
@@ -173,7 +173,7 @@ static double correctLC_(BDS_Point *p1,BDS_Point *p2, GFace *f,
double SCALINGU, double SCALINGV)
{
double l1 = NewGetLc(p1);
- double l2 = NewGetLc(p2);
+ double l2 = NewGetLc(p2);
double l = 0.5 * (l1 + l2);
if(CTX::instance()->mesh.lcFromCurvature){
@@ -193,18 +193,18 @@ static double correctLC_(BDS_Point *p1,BDS_Point *p2, GFace *f,
double NewGetLc(BDS_Edge *e, GFace *f, double SCALINGU, double SCALINGV)
{
double linearLength = computeEdgeLinearLength(e, f, SCALINGU, SCALINGV);
- double l = correctLC_ (e->p1,e->p2,f, SCALINGU, SCALINGV);
+ double l = correctLC_ (e->p1,e->p2,f, SCALINGU, SCALINGV);
return linearLength / l;
}
double NewGetLc(BDS_Point *p1, BDS_Point *p2, GFace *f, double su, double sv)
{
double linearLength = computeEdgeLinearLength(p1, p2, f, su, sv);
- double l = correctLC_ (p1,p2,f, su, sv);
+ double l = correctLC_ (p1,p2,f, su, sv);
return linearLength / l;
}
-void computeMeshSizeFieldAccuracy(GFace *gf, BDS_Mesh &m, double &avg,
+void computeMeshSizeFieldAccuracy(GFace *gf, BDS_Mesh &m, double &avg,
double &max_e, double &min_e, int &nE, int &GS)
{
std::list<BDS_Edge*>::iterator it = m.edges.begin();
@@ -243,18 +243,18 @@ bool edgeSwapTestAngle(BDS_Edge *e, double min_cos)
normal_triangle(n1[0], n1[1], n1[2], norm1);
normal_triangle(n2[0], n2[1], n2[2], norm2);
double cosa;prosca (norm1, norm2, &cosa);
- return cosa > min_cos;
+ return cosa > min_cos;
}
bool evalSwapForOptimize(BDS_Edge *e, GFace *gf, BDS_Mesh &m)
{
- if(e->numfaces() != 2) return false;
+ if(e->numfaces() != 2) return false;
BDS_Point *p11, *p12, *p13;
BDS_Point *p21, *p22, *p23;
BDS_Point *p31, *p32, *p33;
BDS_Point *p41, *p42, *p43;
- swap_config(e, &p11, &p12, &p13, &p21, &p22, &p23, &p31, &p32, &p33, &p41,
+ swap_config(e, &p11, &p12, &p13, &p21, &p22, &p23, &p31, &p32, &p33, &p41,
&p42, &p43);
// First, evaluate what we gain in element quality if the
@@ -269,7 +269,7 @@ bool evalSwapForOptimize(BDS_Edge *e, GFace *gf, BDS_Mesh &m)
bool qualShouldSwap = qb > 2*qa;
bool qualCouldSwap = !(qb < qa * .5) && qb > .05;
- // then evaluate if swap produces smoother surfaces
+ // then evaluate if swap produces smoother surfaces
double norm11[3];
double norm12[3];
double norm21[3];
@@ -281,8 +281,8 @@ bool evalSwapForOptimize(BDS_Edge *e, GFace *gf, BDS_Mesh &m)
double cosa; prosca(norm11, norm12, &cosa);
double cosb; prosca(norm21, norm22, &cosb);
// double smoothIndicator = cosb - cosa;
- // bool smoothShouldSwap = (cosa < 0.1 && cosb > 0.3);
- // bool smoothCouldSwap = !(cosb < cosa * .5);
+ // bool smoothShouldSwap = (cosa < 0.1 && cosb > 0.3);
+ // bool smoothCouldSwap = !(cosb < cosa * .5);
double la = computeEdgeLinearLength(p11, p12);
double la_ = computeEdgeLinearLength(p11, p12, gf, m.scalingU, m.scalingV);
@@ -294,8 +294,8 @@ bool evalSwapForOptimize(BDS_Edge *e, GFace *gf, BDS_Mesh &m)
double distanceIndicator = LA - LB;
bool distanceShouldSwap = (LB < .5 * LA) && LA > 1.e-2;
- bool distanceCouldSwap = !(LB > 2 * LA) || LB < 1.e-2;
-
+ bool distanceCouldSwap = !(LB > 2 * LA) || LB < 1.e-2;
+
if (20 * qa < qb) return true;
// if (LB > .025 && distanceIndicator > 0 && qb > .25) return true;
// if (LB > .05 && distanceIndicator > 0 && qb > .15) return true;
@@ -314,13 +314,13 @@ bool evalSwapForOptimize(BDS_Edge *e, GFace *gf, BDS_Mesh &m)
// if (distanceCouldSwap && qualShouldSwap) return true;
if (cosa < 0 && cosb > 0 && qb > 0.0)
return true;
- return false;
+ return false;
}
bool edgeSwapTestDelaunay(BDS_Edge *e, GFace *gf)
{
BDS_Point *op[2];
-
+
if(!e->p1->config_modified && !e->p2->config_modified) return false;
if(e->numfaces() != 2) return false;
@@ -333,24 +333,24 @@ bool edgeSwapTestDelaunay(BDS_Edge *e, GFace *gf)
double op2x[3] = {op[1]->X, op[1]->Y, op[1]->Z};
double fourth[3];
fourthPoint(p1x, p2x, op1x, fourth);
- double result = robustPredicates::insphere(p1x, p2x, op1x, fourth, op2x) *
- robustPredicates::orient3d(p1x, p2x, op1x, fourth);
+ double result = robustPredicates::insphere(p1x, p2x, op1x, fourth, op2x) *
+ robustPredicates::orient3d(p1x, p2x, op1x, fourth);
return result > 0.;
}
bool edgeSwapTestHighOrder(BDS_Edge *e,GFace *gf)
{
- // must evaluate the swap with the perspectve of
- // the generation of 2 high order elements
+ // must evaluate the swap with the perspectve of
+ // the generation of 2 high order elements
// The rationale is to consider the edges as
- // exactly matching curves and surfaces
+ // exactly matching curves and surfaces
return false;
}
bool edgeSwapTestDelaunayAniso(BDS_Edge *e, GFace *gf, std::set<swapquad> &configs)
{
BDS_Point *op[2];
-
+
if(!e->p1->config_modified && ! e->p2->config_modified) return false;
if(e->numfaces() != 2) return false;
@@ -360,10 +360,10 @@ bool edgeSwapTestDelaunayAniso(BDS_Edge *e, GFace *gf, std::set<swapquad> &confi
swapquad sq(e->p1->iD, e->p2->iD, op[0]->iD, op[1]->iD);
if (configs.find(sq) != configs.end()) return false;
configs.insert(sq);
-
+
double edgeCenter[2] = {0.5 * (e->p1->u + e->p2->u),
0.5 * (e->p1->v + e->p2->v)};
-
+
double p1[2] = {e->p1->u, e->p1->v};
double p2[2] = {e->p2->u, e->p2->v};
double p3[2] = {op[0]->u, op[0]->v};
@@ -372,15 +372,15 @@ bool edgeSwapTestDelaunayAniso(BDS_Edge *e, GFace *gf, std::set<swapquad> &confi
buildMetric(gf, edgeCenter, metric);
if (!inCircumCircleAniso(gf, p1, p2, p3, p4, metric)){
return false;
- }
+ }
return true;
}
bool evalSwap(BDS_Edge *e, double &qa, double &qb)
{
BDS_Point *op[2];
-
- if(e->numfaces() != 2) return false;
+
+ if(e->numfaces() != 2) return false;
e->oppositeof (op);
double qa1 = qmTriangle(e->p1, e->p2, op[0], QMTRI_RHO);
double qa2 = qmTriangle(e->p1, e->p2, op[1], QMTRI_RHO);
@@ -394,18 +394,18 @@ bool evalSwap(BDS_Edge *e, double &qa, double &qb)
int edgeSwapTestQuality(BDS_Edge *e, double fact=1.1, bool force=false)
{
BDS_Point *op[2];
-
+
if (!force)
if(!e->p1->config_modified && ! e->p2->config_modified) return 0;
-
+
if(e->numfaces() != 2) return 0;
-
+
e->oppositeof (op);
if (!force)
if (!edgeSwapTestAngle(e, cos(CTX::instance()->mesh.allowSwapEdgeAngle * M_PI / 180.)))
return -1;
-
+
double qa1 = qmTriangle(e->p1, e->p2, op[0], QMTRI_RHO);
double qa2 = qmTriangle(e->p1, e->p2, op[1], QMTRI_RHO);
double qb1 = qmTriangle(e->p1, op[0], op[1], QMTRI_RHO);
@@ -431,14 +431,14 @@ void swapEdgePass(GFace *gf, BDS_Mesh &m, int &nb_swap)
double qual = (CTX::instance()->mesh.algo2d == ALGO_2D_MESHADAPT_OLD) ? 1 : 5;
int result = edgeSwapTestQuality(*it, qual);
if (CTX::instance()->mesh.algo2d == ALGO_2D_MESHADAPT_OLD){
- if (m.swap_edge(*it, BDS_SwapEdgeTestQuality(true))) nb_swap++;
+ if (m.swap_edge(*it, BDS_SwapEdgeTestQuality(true))) nb_swap++;
}
else if (result >= 0 && edgeSwapTestDelaunay(*it, gf)){
if (m.swap_edge(*it, BDS_SwapEdgeTestQuality(false))) nb_swap++;
}
}
++it;
- }
+ }
}
void delaunayizeBDS(GFace *gf, BDS_Mesh &m, int &nb_swap)
@@ -453,7 +453,7 @@ void delaunayizeBDS(GFace *gf, BDS_Mesh &m, int &nb_swap)
while (1){
if (NN2++ >= NN1)break;
if (!(*it)->deleted){
- if (edgeSwapTestDelaunayAniso(*it, gf, configs)){
+ if (edgeSwapTestDelaunayAniso(*it, gf, configs)){
if (m.swap_edge(*it , BDS_SwapEdgeTestQuality(false))){
NSW++;
}
@@ -466,25 +466,18 @@ void delaunayizeBDS(GFace *gf, BDS_Mesh &m, int &nb_swap)
}
}
-// A test for spheres
-static void midpointsphere(GFace *gf,
- double u1,
- double v1,
- double u2,
- double v2,
- double &u12,
- double &v12,
- SPoint3 ¢er,
- double r)
+/*
+static void midpointsphere(GFace *gf, double u1, double v1, double u2, double v2,
+ double &u12, double &v12, SPoint3 ¢er, double r)
{
GPoint p1 = gf->point(u1, v1);
GPoint p2 = gf->point(u2, v2);
-
+
SVector3 DIR ((p1.x()+p2.x())/2.0 - center.x(),
(p1.y()+p2.y())/2.0 - center.y(),
(p1.z()+p2.z())/2.0 - center.z());
DIR.normalize();
-
+
double X,Y,Z;
double guess [2] = {0.5 * (u1 + u2), 0.5 * (v1 + v2)};
@@ -492,7 +485,6 @@ static void midpointsphere(GFace *gf,
v12 = guess[1];
if ( (v1 > 0.7*M_PI/2 && v2 > 0.7*M_PI/2) ||
(v1 < -0.7*M_PI/2 && v2 < -0.7*M_PI/2) ){
- // printf("coucou\n");
X = center.x() + DIR.x() * r;
Y = center.y() + DIR.y() * r;
Z = center.z() + DIR.z() * r;
@@ -509,14 +501,11 @@ static void midpointsphere(GFace *gf,
u12 = proj.u();
v12 = proj.v();
}
- // printf("%g %g %g -- %g\n",center.x(),center.y(),center.z(),r);
- // printf("%g %g -- %g %g -- %g %g -- %g %g\n",
- // u1, v1, u2, v2, u12, v12, 0.5 * (u1 + u2), 0.5 * (v1 + v2));
- return;
}
+*/
void splitEdgePass(GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split)
-{
+{
std::list<BDS_Edge*>::iterator it = m.edges.begin();
std::vector<std::pair<double, BDS_Edge*> > edges;
@@ -542,12 +531,12 @@ void splitEdgePass(GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split)
V = coord * e->p1->v + (1 - coord) * e->p2->v;
GPoint gpp = gf->point(m.scalingU*U,m.scalingV*V);
- if (gpp.succeeded()){
+ if (gpp.succeeded()){
mid = m.add_point(++m.MAXPOINTNUMBER, gpp.x(),gpp.y(),gpp.z());
mid->u = U;
mid->v = V;
if (backgroundMesh::current() && 0){
- mid->lc() = mid->lcBGM() =
+ mid->lc() = mid->lcBGM() =
backgroundMesh::current()->operator()
((coord * e->p1->u + (1 - coord) * e->p2->u)*m.scalingU,
(coord * e->p1->v + (1 - coord) * e->p2->v)*m.scalingV,
@@ -599,7 +588,7 @@ void collapseEdgePass(GFace *gf, BDS_Mesh &m, double MINE_, int MAXNP, int &nb_c
}
}
-void collapseEdgePassUnSorted(GFace *gf, BDS_Mesh &m, double MINE_, int MAXNP,
+void collapseEdgePassUnSorted(GFace *gf, BDS_Mesh &m, double MINE_, int MAXNP,
int &nb_collaps)
{
int NN1 = m.edges.size();
@@ -608,7 +597,7 @@ void collapseEdgePassUnSorted(GFace *gf, BDS_Mesh &m, double MINE_, int MAXNP,
while (1){
if(NN2++ >= NN1) break;
-
+
if(!(*it)->deleted){
double lone = NewGetLc(*it, gf, m.scalingU, m.scalingV);
if(!(*it)->deleted && (*it)->numfaces() == 2 && lone < MINE_){
@@ -630,14 +619,14 @@ void smoothVertexPass(GFace *gf, BDS_Mesh &m, int &nb_smooth, bool q)
// FIXME SUPER HACK
// return;
std::set<BDS_Point*,PointLessThan>::iterator itp = m.points.begin();
- while(itp != m.points.end()){
+ while(itp != m.points.end()){
if(m.smooth_point_centroid(*itp, gf,q))
nb_smooth ++;
++itp;
}
}
-void refineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
+void refineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
const bool computeNodalSizeField,
std::map<MVertex*, BDS_Point*> *recoverMapInv)
{
@@ -673,7 +662,7 @@ void refineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
int ne = 0;
while(it != ite){
double l = (*it)->length();
- if ((*it)->g && (*it)->g->classif_degree == 1){
+ if ((*it)->g && (*it)->g->classif_degree == 1){
L = ne ? std::max(L, l) : l;
ne++;
}
@@ -694,7 +683,7 @@ void refineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
const double MINE_ = 0.6666, MAXE_ = 1.4;
while (1){
-
+
// we count the number of local mesh modifs.
int nb_split = 0;
int nb_smooth = 0;
@@ -712,7 +701,7 @@ void refineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
if (!(*it)->deleted){
(*it)->p1->config_modified = false;
(*it)->p2->config_modified = false;
- double lone = NewGetLc(*it, gf, m.scalingU, m.scalingV);
+ double lone = NewGetLc(*it, gf, m.scalingU, m.scalingV);
maxL = std::max(maxL, lone);
minL = std::min(minL, lone);
}
@@ -754,7 +743,7 @@ void refineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
t_sw += t7 - t6;
t_col += t4 - t3;
t_sm += t6 - t5;
- m.cleanup();
+ m.cleanup();
IT++;
Msg::Debug(" iter %3d minL %8.3f/%8.3f maxL %8.3f/%8.3f : "
@@ -762,7 +751,7 @@ void refineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
IT, minL, minE, maxL, maxE, nb_split, nb_swap, nb_collaps, nb_smooth);
if (nb_split == 0 && nb_collaps == 0) break;
- }
+ }
double t_total = t_spl + t_sw + t_col + t_sm;
if(!t_total) t_total = 1.e-6;
@@ -792,8 +781,8 @@ void invalidEdgesPeriodic(BDS_Mesh &m, std::map<BDS_Point*, MVertex*> *recoverMa
BDS_Edge *e = *it;
if (!e->deleted && e->numfaces() == 1){
std::map<BDS_Point*, MVertex*>::iterator itp1 = recoverMap->find(e->p1);
- std::map<BDS_Point*, MVertex*>::iterator itp2 = recoverMap->find(e->p2);
- if (itp1 != recoverMap->end() && itp2 != recoverMap->end() &&
+ std::map<BDS_Point*, MVertex*>::iterator itp2 = recoverMap->find(e->p2);
+ if (itp1 != recoverMap->end() && itp2 != recoverMap->end() &&
itp1->second == itp2->second){
degenerated.insert(itp1->second);
}
@@ -808,19 +797,19 @@ void invalidEdgesPeriodic(BDS_Mesh &m, std::map<BDS_Point*, MVertex*> *recoverMa
BDS_Edge *e = *it;
if (!e->deleted && e->numfaces() == 2){
std::map<BDS_Point*, MVertex*>::iterator itp1 = recoverMap->find(e->p1);
- std::map<BDS_Point*, MVertex*>::iterator itp2 = recoverMap->find(e->p2);
- if (itp1 != recoverMap->end() &&
- itp2 != recoverMap->end() &&
+ std::map<BDS_Point*, MVertex*>::iterator itp2 = recoverMap->find(e->p2);
+ if (itp1 != recoverMap->end() &&
+ itp2 != recoverMap->end() &&
itp1->second == itp2->second) toSplit.insert(e);
else if (itp1 != recoverMap->end() && itp2 == recoverMap->end()){
std::pair<MVertex*, BDS_Point*> a(itp1->second, e->p2);
- std::map<std::pair<MVertex*, BDS_Point*>, BDS_Edge*>::iterator itf =
+ std::map<std::pair<MVertex*, BDS_Point*>, BDS_Edge*>::iterator itf =
touchPeriodic.find(a);
if (itf == touchPeriodic.end()) touchPeriodic[a] = e;
else if (degenerated.find(itp1->second) == degenerated.end()){
toSplit.insert(e); toSplit.insert(itf->second);
}
- }
+ }
else if (itp1 == recoverMap->end() && itp2 != recoverMap->end()){
std::pair<MVertex*, BDS_Point*> a(itp2->second, e->p1);
std::map<std::pair<MVertex*, BDS_Point*>, BDS_Edge*>::iterator itf =
@@ -829,7 +818,7 @@ void invalidEdgesPeriodic(BDS_Mesh &m, std::map<BDS_Point*, MVertex*> *recoverMa
else if (degenerated.find(itp2->second) == degenerated.end()){
toSplit.insert(e); toSplit.insert(itf->second);
}
- }
+ }
}
++it;
}
@@ -842,7 +831,7 @@ void invalidEdgesPeriodic(BDS_Mesh &m, std::map<BDS_Point*, MVertex*> *recoverMa
// if not do not allow p1 v and p2 v, split them
// if p1 p2 exists and it is internal, split it
-int solveInvalidPeriodic(GFace *gf, BDS_Mesh &m,
+int solveInvalidPeriodic(GFace *gf, BDS_Mesh &m,
std::map<BDS_Point*, MVertex*> *recoverMap)
{
std::set<BDS_Edge*> toSplit;
@@ -856,7 +845,7 @@ int solveInvalidPeriodic(GFace *gf, BDS_Mesh &m,
BDS_Point *mid ;
mid = m.add_point(++m.MAXPOINTNUMBER,
coord * e->p1->u + (1 - coord) * e->p2->u,
- coord * e->p1->v + (1 - coord) * e->p2->v, gf);
+ coord * e->p1->v + (1 - coord) * e->p2->v, gf);
mid->lcBGM() = BGM_MeshSize(gf,
(coord * e->p1->u + (1 - coord) * e->p2->u) * m.scalingU,
(coord * e->p1->v + (1 - coord) * e->p2->v) * m.scalingV,
@@ -871,7 +860,7 @@ int solveInvalidPeriodic(GFace *gf, BDS_Mesh &m,
return toSplit.size();
}
-void optimizeMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
+void optimizeMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
std::map<BDS_Point*,MVertex*> *recoverMap=0)
{
int nb_swap;
@@ -885,19 +874,19 @@ void optimizeMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
std::list<BDS_Edge*>::iterator it = m.edges.begin();
while (1){
if (NN2++ >= NN1)break;
- if (evalSwapForOptimize(*it, gf, m))
+ if (evalSwapForOptimize(*it, gf, m))
m.swap_edge(*it, BDS_SwapEdgeTestQuality(false));
++it;
}
- m.cleanup();
+ m.cleanup();
int nb_smooth;
smoothVertexPass(gf, m, nb_smooth, true);
}
}
-
+
if (recoverMap){
while(solveInvalidPeriodic(gf, m, recoverMap)){
- }
+ }
}
}
@@ -931,9 +920,9 @@ BDS_Mesh *gmsh2BDS(std::list<GFace*> &l)
reparamMeshVertexOnFace(e->getVertex(j), gf, param);
p[j]->u = param[0];
p[j]->v = param[1];
- m->add_geom(e->getVertex(j)->onWhat()->tag(),
+ m->add_geom(e->getVertex(j)->onWhat()->tag(),
e->getVertex(j)->onWhat()->dim());
- BDS_GeomEntity *g = m->get_geom(e->getVertex(j)->onWhat()->tag(),
+ BDS_GeomEntity *g = m->get_geom(e->getVertex(j)->onWhat()->tag(),
e->getVertex(j)->onWhat()->dim());
p[j]->g = g;
}
@@ -955,9 +944,9 @@ void collapseSmallEdges(GModel &gm)
}
BDS_Mesh *pm = gmsh2BDS(faces);
outputScalarField(pm->triangles, "all.pos", 0);
-
+
for (GModel::eiter eit = gm.firstEdge(); eit != gm.lastEdge(); eit++){
- }
+ }
delete pm;
}
diff --git a/Mesh/meshGFaceBamg.cpp b/Mesh/meshGFaceBamg.cpp
index 5c0498505e0bbf668a50c4b3d65663814cfad009..deb1c1fc13996b4e3d34eb67faa7c6d88271e15e 100644
--- a/Mesh/meshGFaceBamg.cpp
+++ b/Mesh/meshGFaceBamg.cpp
@@ -32,7 +32,7 @@ Mesh2 *Bamg(Mesh2 *Thh, double * args,double *mm11,double *mm12,double *mm22, bo
static void computeMeshMetricsForBamg(GFace *gf, int numV,
- Vertex2 *bamgVertices,
+ Vertex2 *bamgVertices,
double *mm11, double *mm12, double *mm22)
{
// char name[245];
@@ -48,7 +48,7 @@ static void computeMeshMetricsForBamg(GFace *gf, int numV,
double v = bamgVertices[i][1];
GPoint gp = gf->point(SPoint2(u, v));
SMetric3 m = BGM_MeshMetric(gf, u, v, gp.x(), gp.y(), gp.z());
-
+
// compute the derivatives of the parametrization
Pair<SVector3, SVector3> der = gf->firstDer(SPoint2(u, v));
@@ -65,7 +65,7 @@ static void computeMeshMetricsForBamg(GFace *gf, int numV,
bamg::Metric M1(R(0,0),R(1,0),R(1,1));
mm11[i] = M1.a11;
mm12[i] = M1.a21;
- mm22[i] = M1.a22;
+ mm22[i] = M1.a22;
}
}
@@ -84,7 +84,7 @@ void meshGFaceBamg(GFace *gf){
GEdge *gec = (GEdge*)(*it)->getCompound();
mySet.insert(gec);
}
- else{
+ else{
mySet.insert(*it);
}
++it;
@@ -103,17 +103,17 @@ void meshGFaceBamg(GFace *gf){
std::set<MVertex*> all;
std::map<int,MVertex*> recover;
for (unsigned int i = 0; i < gf->triangles.size(); i++){
- for (unsigned int j = 0; j < 3; j++)
+ for (unsigned int j = 0; j < 3; j++)
all.insert(gf->triangles[i]->getVertex(j));
}
- Vertex2 *bamgVertices = new Vertex2[all.size()];
+ Vertex2 *bamgVertices = new Vertex2[all.size()];
int index = 0;
for(std::set<MVertex*>::iterator it = all.begin(); it!=all.end(); ++it){
- if ((*it)->onWhat()->dim() <= 1){
+ if ((*it)->onWhat()->dim() <= 1){
//for(std::set<MVertex*>::iterator it = bcVertex.begin(); it!=bcVertex.end(); ++it){
SPoint2 p;
- bool success = reparamMeshVertexOnFace(*it, gf, p);
+ reparamMeshVertexOnFace(*it, gf, p);
bamgVertices[index][0] = p.x();
bamgVertices[index][1] = p.y();
bamgVertices[index].lab = index;
@@ -127,18 +127,18 @@ void meshGFaceBamg(GFace *gf){
//FIXME : SEAMS should have to be taken into account here !!!
if ((*it)->onWhat()->dim() >= 2){
SPoint2 p;
- bool success = reparamMeshVertexOnFace(*it, gf, p);
+ reparamMeshVertexOnFace(*it, gf, p);
bamgVertices[index][0] = p.x();
bamgVertices[index][1] = p.y();
recover[index] = *it;
(*it)->setIndex(index++);
}
}
-
+
std::vector<MElement*> myParamElems;
std::vector<MVertex*> newVert;
Triangle2 *bamgTriangles = new Triangle2[gf->triangles.size()];
- for (unsigned int i = 0; i < gf->triangles.size(); i++){
+ for (unsigned int i = 0; i < gf->triangles.size(); i++){
int nodes [3] = {gf->triangles[i]->getVertex(0)->getIndex(),
gf->triangles[i]->getVertex(1)->getIndex(),
gf->triangles[i]->getVertex(2)->getIndex()};
@@ -149,8 +149,8 @@ void meshGFaceBamg(GFace *gf){
double v2(bamgVertices[nodes[1]][1]);
double v3(bamgVertices[nodes[2]][1]);
if (hasCompounds){
- MVertex *vv1 = new MVertex(u1,v1,0.0);
- MVertex *vv2 = new MVertex(u2,v2,0.0);
+ MVertex *vv1 = new MVertex(u1,v1,0.0);
+ MVertex *vv2 = new MVertex(u2,v2,0.0);
MVertex *vv3 = new MVertex(u3,v3,0.0);
newVert.push_back(vv1);
newVert.push_back(vv2);
@@ -166,7 +166,7 @@ void meshGFaceBamg(GFace *gf){
}
bamgTriangles[i].init(bamgVertices, nodes, gf->tag());
}
-
+
int numEdges = 0;
for (std::list<GEdge*>::iterator it = edges.begin(); it != edges.end(); ++it){
numEdges += (*it)->lines.size();
@@ -177,12 +177,13 @@ void meshGFaceBamg(GFace *gf){
for (std::list<GEdge*>::iterator it = edges.begin(); it != edges.end(); ++it){
for (unsigned int i = 0; i < (*it)->lines.size(); ++i){
int nodes [2] = {(*it)->lines[i]->getVertex(0)->getIndex(),
- (*it)->lines[i]->getVertex(1)->getIndex()};
+ (*it)->lines[i]->getVertex(1)->getIndex()};
bamgBoundary[count].init(bamgVertices, nodes, (*it)->tag());
- bamgBoundary[count++].lab = count;
+ bamgBoundary[count].lab = count;
+ count++;
}
}
-
+
Mesh2 *bamgMesh = new Mesh2 (all.size(), gf->triangles.size(), numEdges,
bamgVertices, bamgTriangles, bamgBoundary);
@@ -213,13 +214,13 @@ void meshGFaceBamg(GFace *gf){
// fclose(fi);
// //END EMI PRINT TRIS
}
-
+
Mesh2 *refinedBamgMesh = 0;
int iterMax = 11;
for (int k= 0; k < iterMax; k++){
-
+
int nbVert = bamgMesh->nv;
-
+
double *mm11 = new double[nbVert];
double *mm12 = new double[nbVert];
double *mm22 = new double[nbVert];
@@ -228,8 +229,8 @@ void meshGFaceBamg(GFace *gf){
args[16] = CTX::instance()->mesh.anisoMax;
args[ 7] = CTX::instance()->mesh.smoothRatio;
//args[ 21] = 90.0;//cutoffrad = 90 degree
- computeMeshMetricsForBamg (gf, nbVert, bamgMesh->vertices, mm11,mm12,mm22);
-
+ computeMeshMetricsForBamg (gf, nbVert, bamgMesh->vertices, mm11,mm12,mm22);
+
try{
refinedBamgMesh = Bamg(bamgMesh, args, mm11, mm12, mm22, false);
Msg::Info("bamg succeeded %d vertices %d triangles",
@@ -242,7 +243,7 @@ void meshGFaceBamg(GFace *gf){
delete [] mm11;
delete [] mm12;
delete [] mm22;
-
+
int nT = bamgMesh->nt;
int nTnow = refinedBamgMesh->nt;
@@ -250,7 +251,7 @@ void meshGFaceBamg(GFace *gf){
bamgMesh = refinedBamgMesh;
if (fabs((double)(nTnow - nT)) < 0.01 * nT) break;
}
-
+
std::map<int,MVertex*> yetAnother;
for (int i = 0; i < refinedBamgMesh->nv; i++){
Vertex2 &v = refinedBamgMesh->vertices[i];
@@ -260,12 +261,12 @@ void meshGFaceBamg(GFace *gf){
// }
//If point not found because compound edges have been remeshed and boundary triangles have changed
//then we call our new octree
- if ( !gp.succeeded() && hasCompounds){
+ if ( !gp.succeeded() && hasCompounds){
double uvw[3] = {v[0],v[1], 0.0};
double UV[3];
double initialTol = MElement::getTolerance();
- MElement::setTolerance(1.e-2);
- MElement *e = _octree->find(v[0],v[1], 0.0, -1);
+ MElement::setTolerance(1.e-2);
+ MElement *e = _octree->find(v[0],v[1], 0.0, -1);
MElement::setTolerance(initialTol);
if (e){
e->xyz2uvw(uvw,UV);
@@ -291,7 +292,7 @@ void meshGFaceBamg(GFace *gf){
}
}
- for (unsigned int i = 0; i < gf->triangles.size(); i++){
+ for (unsigned int i = 0; i < gf->triangles.size(); i++){
delete gf->triangles[i];
}
gf->triangles.clear();
@@ -302,9 +303,9 @@ void meshGFaceBamg(GFace *gf){
Vertex2 &v3 = t[2];
gf->triangles.push_back(new MTriangle(yetAnother[(*refinedBamgMesh)(v1)],
yetAnother[(*refinedBamgMesh)(v2)],
- yetAnother[(*refinedBamgMesh)(v3)]));
+ yetAnother[(*refinedBamgMesh)(v3)]));
}
-
+
//delete pointers
if (refinedBamgMesh) delete refinedBamgMesh;
if ( _octree) delete _octree;
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index 42ce669dadd20e381139eefadbf488b8ec9fb5ea..888d1183559de13ef29bad58ea1ac1bd6897bd1b 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -27,8 +27,9 @@ double LIMIT_ = 0.5 * sqrt(2.0) * 1;
int MTri3::radiusNorm = 2;
-// This function controls the frontal algorithm
-static bool isBoundary(MTri3 *t, double limit_, int &active){
+/*
+static bool isBoundary(MTri3 *t, double limit_, int &active)
+{
if (t->isDeleted()) return false;
for (active = 0; active < 3; active++){
MTri3 *neigh = t->getNeigh(active);
@@ -36,6 +37,7 @@ static bool isBoundary(MTri3 *t, double limit_, int &active){
}
return false;
}
+*/
static bool isActive(MTri3 *t, double limit_, int &active)
{
@@ -66,7 +68,6 @@ static bool isActive(MTri3 *t, double limit_, int &i, std::set<MEdge,Less_Edge>
return false;
}
-
static void updateActiveEdges(MTri3 *t, double limit_, std::set<MEdge,Less_Edge> &front)
{
if (t->isDeleted()) return;
@@ -81,7 +82,6 @@ static void updateActiveEdges(MTri3 *t, double limit_, std::set<MEdge,Less_Edge>
}
}
-
bool circumCenterMetricInTriangle(MTriangle *base, const double *metric,
const std::vector<double> &Us,
const std::vector<double> &Vs)
@@ -259,7 +259,9 @@ int inCircumCircleAniso(GFace *gf, MTriangle *base,
return d2 < Radius2;
}
-MTri3::MTri3(MTriangle *t, double lc, SMetric3 *metric, const std::vector<double> *Us, const std::vector<double> *Vs, GFace *gf)
+MTri3::MTri3(MTriangle *t, double lc, SMetric3 *metric,
+ const std::vector<double> *Us,
+ const std::vector<double> *Vs, GFace *gf)
: deleted(false), base(t)
{
neigh[0] = neigh[1] = neigh[2] = 0;
@@ -290,7 +292,8 @@ MTri3::MTri3(MTriangle *t, double lc, SMetric3 *metric, const std::vector<double
double midpoint[2] = {(p1[0]+p2[0]+p3[0])/3.0,(p1[1]+p2[1]+p3[1])/3.0};
- double quadAngle = backgroundMesh::current() ? backgroundMesh::current()->getAngle (midpoint[0],midpoint[1],0) : 0.0;
+ double quadAngle = backgroundMesh::current() ?
+ backgroundMesh::current()->getAngle (midpoint[0],midpoint[1],0) : 0.0;
double x0 = p1[0] * cos(quadAngle) + p1[1] * sin(quadAngle);
double y0 = -p1[0] * sin(quadAngle) + p1[1] * cos(quadAngle);
@@ -317,7 +320,6 @@ MTri3::MTri3(MTriangle *t, double lc, SMetric3 *metric, const std::vector<double
}
}
-
int MTri3::inCircumCircle(const double *p) const
{
double pa[3] =
@@ -540,7 +542,6 @@ bool insertVertex(bool force, GFace *gf, MVertex *v, double *param , MTri3 *t,
}
else{
recurFindCavityAniso(gf, shell, cavity, metric, param, t, Us, Vs);
- // printf("RECURSIVELY FIND A CAVITY %d %d \n",shell.size(),cavity.size());
}
// check that volume is conserved
@@ -579,13 +580,15 @@ bool insertVertex(bool force, GFace *gf, MVertex *v, double *param , MTri3 *t,
double d2 = sqrt((it->v[1]->x() - v->x()) * (it->v[1]->x() - v->x()) +
(it->v[1]->y() - v->y()) * (it->v[1]->y() - v->y()) +
(it->v[1]->z() - v->z()) * (it->v[1]->z() - v->z()));
- const double MID[3] = {0.5*(it->v[0]->x()+it->v[1]->x()),0.5*(it->v[0]->y()+it->v[1]->y()),0.5*(it->v[0]->z()+it->v[1]->z())};
- double d3 = sqrt((MID[0] - v->x()) * (MID[0] - v->x()) + (MID[1] - v->y()) * (MID[1] - v->y()) + (MID[2] - v->z()) * (MID[2] - v->z()));
+ const double MID[3] = {0.5*(it->v[0]->x()+it->v[1]->x()),
+ 0.5*(it->v[0]->y()+it->v[1]->y()),
+ 0.5*(it->v[0]->z()+it->v[1]->z())};
+ double d3 = sqrt((MID[0] - v->x()) * (MID[0] - v->x()) +
+ (MID[1] - v->y()) * (MID[1] - v->y()) +
+ (MID[2] - v->z()) * (MID[2] - v->z()));
if ((d1 < LL * .25 || d2 < LL * .25 || d3 < LL * .25) && !force) {
- // printf("TOO CLOSE %g %g %g %g %g %g\n",d1,d2,d3,LL,lc,lcBGM);
onePointIsTooClose = true;
}
- // if (t4->getRadius () < LIMIT_ / 2) onePointIsTooClose = true;
newTris[k++] = t4;
// all new triangles are pushed front in order to be able to
@@ -604,7 +607,6 @@ bool insertVertex(bool force, GFace *gf, MVertex *v, double *param , MTri3 *t,
!onePointIsTooClose){
connectTris(new_cavity.begin(), new_cavity.end());
allTets.insert(newTris, newTris + shell.size());
- // printf("shell.size() = %d triangles.size() = %d \n",shell.size(),allTets.size());
if (activeTets){
for (std::list<MTri3*>::iterator i = new_cavity.begin(); i != new_cavity.end(); ++i){
int active_edge;
@@ -617,10 +619,9 @@ bool insertVertex(bool force, GFace *gf, MVertex *v, double *param , MTri3 *t,
delete [] newTris;
return true;
}
+
// The cavity is NOT star shaped
else{
- // printf("cavity not star shaped %22.15E vs %22.15E\n",oldVolume,newVolume);
- // printf("shell.size() = %d triangles.size() = %d \n",shell.size(),allTets.size());
for (unsigned int i = 0; i < shell.size(); i++) delete newTris[i];
delete [] newTris;
ittet = cavity.begin();
@@ -687,26 +688,32 @@ static MTri3* search4Triangle (MTri3 *t, double pt[2],
SPoint3 q1(pt[0],pt[1],0);
int ITER = 0;
while (1){
- // printf("%d %d %d\n",t->tri()->getVertex(0)->getIndex(),t->tri()->getVertex(1)->getIndex() ,t->tri()->getVertex(2)->getIndex());
- SPoint3 q2((Us[t->tri()->getVertex(0)->getIndex()] + Us[t->tri()->getVertex(1)->getIndex()] + Us[t->tri()->getVertex(2)->getIndex()])/3.0,
- (Vs[t->tri()->getVertex(0)->getIndex()] + Vs[t->tri()->getVertex(1)->getIndex()] + Vs[t->tri()->getVertex(2)->getIndex()])/3.0,
+ SPoint3 q2((Us[t->tri()->getVertex(0)->getIndex()] +
+ Us[t->tri()->getVertex(1)->getIndex()] +
+ Us[t->tri()->getVertex(2)->getIndex()]) / 3.0,
+ (Vs[t->tri()->getVertex(0)->getIndex()] +
+ Vs[t->tri()->getVertex(1)->getIndex()] +
+ Vs[t->tri()->getVertex(2)->getIndex()]) / 3.0,
0);
int i;
- for (i=0;i<3;i++){
- SPoint3 p1 ( Us[t->tri()->getVertex(i == 0 ? 2 : i-1)->getIndex()], Vs[t->tri()->getVertex(i == 0 ? 2 : i-1)->getIndex()], 0);
- SPoint3 p2 ( Us[t->tri()->getVertex(i)->getIndex()], Vs[t->tri()->getVertex(i)->getIndex()], 0);
+ for (i = 0; i < 3; i++){
+ SPoint3 p1 (Us[t->tri()->getVertex(i == 0 ? 2 : i-1)->getIndex()],
+ Vs[t->tri()->getVertex(i == 0 ? 2 : i-1)->getIndex()], 0);
+ SPoint3 p2 (Us[t->tri()->getVertex(i)->getIndex()],
+ Vs[t->tri()->getVertex(i)->getIndex()], 0);
double xcc[2];
- if (intersection_segments(p1,p2,q1,q2,xcc))break;
+ if (intersection_segments(p1, p2, q1, q2, xcc)) break;
}
- if (i>=3)break;
- t = t->getNeigh(i);
- if (!t)break;
+ if (i >= 3) break;
+ t = t->getNeigh(i);
+ if (!t) break;
bool inside = invMapUV(t->tri(), pt, Us, Vs, uv, 1.e-8);
- if (inside) {return t;}
- if (ITER++ > AllTris.size())break;
+ if (inside) return t;
+ if (ITER++ > (int)AllTris.size()) break;
}
return 0;
- for(std::set<MTri3*,compareTri3Ptr>::iterator itx = AllTris.begin(); itx != AllTris.end();++itx){
+ for(std::set<MTri3*,compareTri3Ptr>::iterator itx = AllTris.begin();
+ itx != AllTris.end();++itx){
if (!(*itx)->isDeleted()){
inside = invMapUV((*itx)->tri(), pt, Us, Vs, uv, 1.e-8);
if (inside){
@@ -761,15 +768,15 @@ static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it
}
}
else {
- ptin = search4Triangle (worst, center, Us, Vs, AllTris,uv);
- // if (!ptin){
- // printf("strange : %g %g seems to be out of the domain for face %d\n",center[0],center[1],gf->tag());
- // }
- if (ptin)inside = true;
+ ptin = search4Triangle (worst, center, Us, Vs, AllTris,uv);
+ // if (!ptin){
+ // printf("strange : %g %g seems to be out of the domain for face %d\n",
+ // center[0],center[1],gf->tag());
+ // }
+ if (ptin) inside = true;
}
- // if (!ptin)ptin = worst;
- if ( inside) {
+ if (inside) {
// we use here local coordinates as real coordinates
// x,y and z will be computed hereafter
// Msg::Info("Point is inside");
@@ -799,10 +806,8 @@ static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it
(false, gf, v, center, worst, AllTris,ActiveTris, vSizes, vSizesBGM,vMetricsBGM,
Us, Vs, metric) ) {
- MTriangle *base = worst->tri();
- Msg::Debug("Point %g %g cannot be inserted because %d", center[0], center[1], p.succeeded() );
-
- // Msg::Info("Point %g %g cannot be inserted because %d", center[0], center[1], p.succeeded() );
+ Msg::Debug("Point %g %g cannot be inserted because %d",
+ center[0], center[1], p.succeeded() );
AllTris.erase(it);
worst->forceRadius(-1);
@@ -817,8 +822,8 @@ static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it
}
}
else {
+ /*
MTriangle *base = worst->tri();
- /*
Msg::Info("Point %g %g is outside (%g %g , %g %g , %g %g)",
center[0], center[1],
Us[base->getVertex(0)->getIndex()],
@@ -835,7 +840,6 @@ static bool insertAPoint(GFace *gf, std::set<MTri3*,compareTri3Ptr>::iterator it
}
}
-
void bowyerWatson(GFace *gf, int MAXPNT)
{
std::set<MTri3*,compareTri3Ptr> AllTris;
@@ -873,7 +877,8 @@ void bowyerWatson(GFace *gf, int MAXPNT)
Msg::Debug("%7d points created -- Worst tri radius is %8.3f",
vSizes.size(), worst->getRadius());
double center[2],metric[3],r2;
- if (worst->getRadius() < /*1.333333/(sqrt(3.0))*/0.5 * sqrt(2.0) || vSizes.size() > MAXPNT) break;
+ if (worst->getRadius() < /*1.333333/(sqrt(3.0))*/0.5 * sqrt(2.0) ||
+ (int)vSizes.size() > MAXPNT) break;
circUV(worst->tri(), Us, Vs, center, gf);
MTriangle *base = worst->tri();
double pa[2] = {(Us[base->getVertex(0)->getIndex()] +
@@ -903,11 +908,9 @@ void bowyerWatson(GFace *gf, int MAXPNT)
}
/*
- Let's try a frontal delaunay approach now that the delaunay mesher is stable
-
- We use the approach of Rebay (JCP 1993) that we extend to anisotropic metrics
-
- The point isertion is done on the Voronoi Edges
+ Let's try a frontal delaunay approach now that the delaunay mesher is stable.
+ We use the approach of Rebay (JCP 1993) that we extend to anisotropic metrics.
+ The point isertion is done on the Voronoi Edges.
*/
double lengthInfniteNorm(const double p[2], const double q[2],
@@ -925,7 +928,8 @@ double lengthInfniteNorm(const double p[2], const double q[2],
void circumCenterInfinite (MTriangle *base, double quadAngle,
const std::vector<double> &Us,
- const std::vector<double> &Vs, double *x) {
+ const std::vector<double> &Vs, double *x)
+{
double pa[2] = {Us[base->getVertex(0)->getIndex()],
Vs[base->getVertex(0)->getIndex()]};
double pb[2] = {Us[base->getVertex(1)->getIndex()],
@@ -946,7 +950,6 @@ void circumCenterInfinite (MTriangle *base, double quadAngle,
x[1] = 0.5 * (ymax-ymin);
}
-
static double lengthMetric(const double p[2], const double q[2],
const double metric[3])
{
@@ -987,7 +990,8 @@ double optimalPointFrontal (GFace *gf,
std::vector<double> &vSizes,
std::vector<double> &vSizesBGM,
double newPoint[2],
- double metric[3]){
+ double metric[3])
+{
double center[2],r2;
MTriangle *base = worst->tri();
circUV(base, Us, Vs, center, gf);
@@ -1002,7 +1006,6 @@ double optimalPointFrontal (GFace *gf,
// compute the middle point of the edge
int ip1 = active_edge - 1 < 0 ? 2 : active_edge - 1;
int ip2 = active_edge;
- int ip3 = (active_edge+1)%3;
double P[2] = {Us[base->getVertex(ip1)->getIndex()],
Vs[base->getVertex(ip1)->getIndex()]};
@@ -1037,7 +1040,8 @@ double optimalPointFrontal (GFace *gf,
// double d = (rhoM_hat + sqrt (rhoM_hat * rhoM_hat - p * p)) / RATIO;
//const double d = 100./RATIO;
- // printf("(%g %g) (%g %g) %g %g %g %g %g %g\n",P[0],P[1],Q[0],Q[1],RATIO,p,q,rhoM,rhoM_hat,d);
+ // printf("(%g %g) (%g %g) %g %g %g %g %g %g\n",
+ // P[0],P[1],Q[0],Q[1],RATIO,p,q,rhoM,rhoM_hat,d);
// printf("size %12.5E\n",vSizesBGM[base->getVertex(ip1)->getIndex()]);
const double rhoM1 = 0.5*
(vSizes[base->getVertex(ip1)->getIndex()] +
@@ -1087,9 +1091,11 @@ void optimalPointFrontalB (GFace *gf,
std::vector<double> &vSizes,
std::vector<double> &vSizesBGM,
double newPoint[2],
- double metric[3]){
+ double metric[3])
+{
// as a starting point, let us use the "fast algo"
- double d = optimalPointFrontal (gf,worst,active_edge,Us,Vs,vSizes,vSizesBGM,newPoint,metric);
+ double d = optimalPointFrontal (gf,worst,active_edge,Us,Vs,vSizes,
+ vSizesBGM,newPoint,metric);
int ip1 = active_edge - 1 < 0 ? 2 : active_edge - 1;
int ip2 = active_edge;
MVertex *v1 = worst->tri()->getVertex(ip1);
@@ -1134,7 +1140,6 @@ void optimalPointFrontalB (GFace *gf,
else {
Msg::Debug("--- Non optimal point found -----------");
}
- // fclose(f);
}
void bowyerWatsonFrontal(GFace *gf)
@@ -1160,11 +1165,10 @@ void bowyerWatsonFrontal(GFace *gf)
else if ((*it)->getRadius() < LIMIT_)break;
}
-
// insert points
while (1){
-
- /* if(ITER % 100== 0){
+
+ /* if(ITER % 100== 0){
char name[245];
sprintf(name,"delfr2d%d-ITER%d.pos",gf->tag(),ITER);
_printTris (name, AllTris, Us,Vs,true);
@@ -1197,11 +1201,11 @@ void bowyerWatsonFrontal(GFace *gf)
*/
}
-// char name[245];
-// sprintf(name,"frontal%d-real.pos", gf->tag());
-// _printTris (name, AllTris, Us, Vs,false);
-// sprintf(name,"frontal%d-param.pos", gf->tag());
-// _printTris (name, AllTris, Us, Vs,true);
+ // char name[245];
+ // sprintf(name,"frontal%d-real.pos", gf->tag());
+ // _printTris (name, AllTris, Us, Vs,false);
+ // sprintf(name,"frontal%d-param.pos", gf->tag());
+ // _printTris (name, AllTris, Us, Vs,true);
transferDataStructure(gf, AllTris, Us, Vs);
// in case of boundary layer meshing
#if defined(HAVE_ANN)
@@ -1217,7 +1221,6 @@ void bowyerWatsonFrontal(GFace *gf)
#endif
}
-
void optimalPointFrontalQuad (GFace *gf,
MTri3* worst,
int active_edge,
@@ -1226,7 +1229,8 @@ void optimalPointFrontalQuad (GFace *gf,
std::vector<double> &vSizes,
std::vector<double> &vSizesBGM,
double newPoint[2],
- double metric[3]){
+ double metric[3])
+{
MTriangle *base = worst->tri();
int ip1 = active_edge - 1 < 0 ? 2 : active_edge - 1;
int ip2 = active_edge;
@@ -1290,7 +1294,6 @@ void optimalPointFrontalQuad (GFace *gf,
npy = temp;
}
-
newPoint[0] = midpoint[0] + cos(quadAngle) * npx - sin(quadAngle) * npy;
newPoint[1] = midpoint[1] + sin(quadAngle) * npx + cos(quadAngle) * npy;
@@ -1301,7 +1304,6 @@ void optimalPointFrontalQuad (GFace *gf,
}
}
-
void optimalPointFrontalQuadB (GFace *gf,
MTri3* worst,
int active_edge,
@@ -1310,9 +1312,8 @@ void optimalPointFrontalQuadB (GFace *gf,
std::vector<double> &vSizes,
std::vector<double> &vSizesBGM,
double newPoint[2],
- double metric[3]){
-
-
+ double metric[3])
+{
optimalPointFrontalQuad (gf,worst,active_edge,Us,Vs,vSizes,vSizesBGM,newPoint,metric);
return;
MTriangle *base = worst->tri();
@@ -1325,8 +1326,6 @@ void optimalPointFrontalQuadB (GFace *gf,
Vs[base->getVertex(ip1)->getIndex()]};
double Q[2] = {Us[base->getVertex(ip2)->getIndex()],
Vs[base->getVertex(ip2)->getIndex()]};
- double O[2] = {Us[base->getVertex(ip3)->getIndex()],
- Vs[base->getVertex(ip3)->getIndex()]};
double midpoint[2] = {0.5 * (P[0] + Q[0]), 0.5 * (P[1] + Q[1])};
// compute background mesh data
@@ -1334,7 +1333,7 @@ void optimalPointFrontalQuadB (GFace *gf,
// double quadAngle = 0;
double center[2];
circumCenterInfinite (base, quadAngle,Us,Vs,center);
-
+
// rotate the points with respect to the angle
double XP1 = 0.5*(Q[0] - P[0]);
double YP1 = 0.5*(Q[1] - P[1]);
@@ -1357,15 +1356,15 @@ void optimalPointFrontalQuadB (GFace *gf,
(Q - P) . (x,y) = 0 --> (-x-tx,2x-y-ty) . (x,y) = 0 --> -x^2 -tx^2 + 2xy - y^2 -t y^2 = 0
--> t (x^2 + y^2) = -x^2 + 2xy - y^2 --> t = - (x-y)^2 / (x^2 + y^2) = -1 + 2 xy /(x^2+y^2)
- The L2 distance between Q and the line is
-
- h^2 = (-x - tx)^2 + (2x-y-ty)^2 = x^2 (1+t)^2 +
-
+ The L2 distance between Q and the line is
+
+ h^2 = (-x - tx)^2 + (2x-y-ty)^2 = x^2 (1+t)^2 +
+
*/
if (fabs(xp)>=fabs(yp)){
- t = -1. + 2*fabs(xp*yp)/(xp*xp + yp*yp);
+ t = -1. + 2*fabs(xp*yp)/(xp*xp + yp*yp);
factor = fabs(xp) / sqrt(xp*xp+yp*yp);
}
else {
@@ -1386,10 +1385,10 @@ void optimalPointFrontalQuadB (GFace *gf,
N2.normalize();
if (dot (N2,O3-PMID) < 0) N2 = N2*(-1.0);
- const double rhoM1 = 0.5 *
+ const double rhoM1 = 0.5 *
(vSizes[base->getVertex(ip1)->getIndex()] +
vSizes[base->getVertex(ip2)->getIndex()] ) ;// * RATIO;
- const double rhoM2 = 0.5 *
+ const double rhoM2 = 0.5 *
(vSizesBGM[base->getVertex(ip1)->getIndex()] +
vSizesBGM[base->getVertex(ip2)->getIndex()] ) ;// * RATIO;
const double a = Extend1dMeshIn2dSurfaces() ? std::min(rhoM1, rhoM2) : rhoM2;
@@ -1417,7 +1416,7 @@ void optimalPointFrontalQuadB (GFace *gf,
else {
Msg::Warning("--- NEVER GO THERE -----------");
}
- }
+ }
}
}
double uvt[3] = {newPoint[0],newPoint[1],0.0};
@@ -1431,22 +1430,21 @@ void optimalPointFrontalQuadB (GFace *gf,
}
else {
newPoint[0] = -10000;
- newPoint[1] = 100000;
+ newPoint[1] = 100000;
Msg::Warning("--- Non optimal point found -----------");
}
// buildMetric(gf, newPoint, metric);
}
-
-void buildBackGroundMesh (GFace *gf) {
-
+void buildBackGroundMesh (GFace *gf)
+{
// printf("build bak mesh\n");
quadsToTriangles(gf, 100000);
if (!backgroundMesh::current()) {
std::vector<MTriangle*> TR;
// std::vector<MQuadrangle*> QR;
- for(int i=0;i<gf->triangles.size();i++){
+ for(unsigned int i = 0; i < gf->triangles.size(); i++){
TR.push_back(new MTriangle(gf->triangles[i]->getVertex(0),
gf->triangles[i]->getVertex(1),
gf->triangles[i]->getVertex(2)));
@@ -1540,7 +1538,7 @@ void bowyerWatsonFrontalLayers(GFace *gf, bool quad)
if (!ActiveTris.size())break;
- /* if (1 || gf->tag() == 1900){
+ /* if (1 || gf->tag() == 1900){
char name[245];
sprintf(name,"x_GFace_%d_Layer_%d.pos",gf->tag(),ITER);
_printTris (name, AllTris, Us,Vs,true);
@@ -1550,30 +1548,33 @@ void bowyerWatsonFrontalLayers(GFace *gf, bool quad)
MTri3 *worst = (*WORST_ITER);
ActiveTris.erase(WORST_ITER);
- if (!worst->isDeleted() && (ITERATION > max_layers ? isActive(worst, LIMIT_, active_edge) : isActive(worst, LIMIT_, active_edge,&_front) ) &&
+ if (!worst->isDeleted() &&
+ (ITERATION > max_layers ? isActive(worst, LIMIT_, active_edge) :
+ isActive(worst, LIMIT_, active_edge,&_front) ) &&
worst->getRadius() > LIMIT_){
- // for (active_edge = 0 ; active_edge < 0 ; active_edge ++){
- // if (active_edges[active_edge])break;
- // }
- // Msg::Info("%7d points created -- Worst tri infinite radius is %8.3f -- front size %6d",
- // vSizes.size(), worst->getRadius(),_front.size());
+ // for (active_edge = 0 ; active_edge < 0 ; active_edge ++){
+ // if (active_edges[active_edge])break;
+ // }
+ // Msg::Info("%7d points created -- Worst tri infinite radius is %8.3f -- "
+ // "front size %6d", vSizes.size(), worst->getRadius(),_front.size());
if(ITER++ % 5000 == 0)
- Msg::Debug("%7d points created -- Worst tri infinite radius is %8.3f -- front size %6d",
- vSizes.size(), worst->getRadius(),_front.size());
+ Msg::Debug("%7d points created -- Worst tri infinite radius is %8.3f -- "
+ "front size %6d", vSizes.size(), worst->getRadius(),_front.size());
// compute the middle point of the edge
double newPoint[2],metric[3]={1,0,1};
- if (quad)optimalPointFrontalQuadB (gf,worst,active_edge,Us,Vs,vSizes,vSizesBGM,newPoint,metric);
- else optimalPointFrontalB (gf,worst,active_edge,Us,Vs,vSizes,vSizesBGM,newPoint,metric);
-
+ if (quad) optimalPointFrontalQuadB (gf,worst,active_edge,Us,Vs,vSizes,
+ vSizesBGM,newPoint,metric);
+ else optimalPointFrontalB (gf,worst,active_edge,Us,Vs,vSizes,
+ vSizesBGM,newPoint,metric);
// printf("start INSERT A POINT %g %g \n",newPoint[0],newPoint[1]);
insertAPoint(gf, AllTris.end(), newPoint, 0, Us, Vs, vSizes,
vSizesBGM, vMetricsBGM, AllTris, &ActiveTris, worst);
- // else if (!worst->isDeleted() && worst->getRadius() > LIMIT_){
- // ActiveTrisNotInFront.insert(worst);
- // }
- // printf("-----------------> size %d\n",AllTris.size());
+ // else if (!worst->isDeleted() && worst->getRadius() > LIMIT_){
+ // ActiveTrisNotInFront.insert(worst);
+ // }
+ // printf("-----------------> size %d\n",AllTris.size());
/*
if(ITER % 1== 0){
@@ -1588,28 +1589,27 @@ void bowyerWatsonFrontalLayers(GFace *gf, bool quad)
}
}
_front.clear();
- it = ActiveTrisNotInFront.begin();
- //it = AllTris.begin();
- for ( ; it!=ActiveTrisNotInFront.end();++it){
- // for ( ; it!=AllTris.end();++it){
+ it = ActiveTrisNotInFront.begin();
+ //it = AllTris.begin();
+ for ( ; it!=ActiveTrisNotInFront.end();++it){
+ // for ( ; it!=AllTris.end();++it){
if((*it)->getRadius() > LIMIT_ && isActive(*it,LIMIT_,active_edge)){
ActiveTris.insert(*it);
updateActiveEdges(*it, LIMIT_, _front);
}
}
- // Msg::Info("%d active tris %d front edges %d not in front",ActiveTris.size(),_front.size(),ActiveTrisNotInFront.size());
- if (!ActiveTris.size())break;
+ // Msg::Info("%d active tris %d front edges %d not in front",
+ // ActiveTris.size(),_front.size(),ActiveTrisNotInFront.size());
+ if (!ActiveTris.size()) break;
}
- char name[245];
- // sprintf(name,"frontal%d-real.pos", gf->tag());
- // _printTris (name, AllTris, Us, Vs,false);
- // sprintf(name,"frontal%d-param.pos", gf->tag());
- // _printTris (name, AllTris, Us, Vs,true);
+ // char name[245];
+ // sprintf(name,"frontal%d-real.pos", gf->tag());
+ // _printTris (name, AllTris, Us, Vs,false);
+ // sprintf(name,"frontal%d-param.pos", gf->tag());
+ // _printTris (name, AllTris, Us, Vs,true);
transferDataStructure(gf, AllTris, Us, Vs);
MTri3::radiusNorm = 2;
LIMIT_ = 0.5 * sqrt(2.0) * 1;
backgroundMesh::unset();
}
-
-
diff --git a/Mesh/meshGFaceLloyd.cpp b/Mesh/meshGFaceLloyd.cpp
index d31f9f321ec0707b8013ac48a81c30ebc8da45c7..df06bbeee90592a7a1c2cac3b3eb4c9c54e1d3df 100644
--- a/Mesh/meshGFaceLloyd.cpp
+++ b/Mesh/meshGFaceLloyd.cpp
@@ -541,7 +541,7 @@ void smoothing::optimize_face(GFace* gf){
gf->additionalVertices.end());
// clear the list of additional vertices
gf->additionalVertices.clear();
-
+
free(initial_conditions);
free(variables_scales);
}
@@ -1145,7 +1145,7 @@ void lpcvt::clip_cells(DocRecord& triangulator,GFace* gf){
}
void lpcvt::clear(){
- int i;
+ unsigned int i;
for(i=0;i<fifo.size();i++){
fifo.pop();
}
@@ -1192,7 +1192,7 @@ void lpcvt::print_voronoi1(){
}
void lpcvt::print_voronoi2(){
- int i;
+ unsigned int i;
int j;
int num;
SPoint2 p1,p2;
@@ -1322,7 +1322,7 @@ void lpcvt::write(DocRecord& triangulator,GFace* gf,int p){
}
void lpcvt::eval(DocRecord& triangulator,std::vector<SVector3>& gradients,double& energy,int p){
- int i;
+ unsigned int i;
int index;
int index1;
int index2;
@@ -2067,7 +2067,7 @@ segment segment_list::get_segment(int index){
}
bool segment_list::add_segment(int index1,int index2,int reference){
- int i;
+ unsigned int i;
for(i=0;i<segments.size();i++){
if(segments[i].equal(index1,index2)) return 0;
}
diff --git a/Mesh/meshGFaceOptimize.cpp b/Mesh/meshGFaceOptimize.cpp
index 87b4ad663b4eb71b58d04aadf70ad27a12b6ee13..0af7a0d05f52eb9d5cc77961e77c38a3d05e7e09 100644
--- a/Mesh/meshGFaceOptimize.cpp
+++ b/Mesh/meshGFaceOptimize.cpp
@@ -53,7 +53,8 @@ static int diff2 (MElement *q1, MElement *q2)
return 0;
}
-static int SANITY_(GFace *gf,MQuadrangle *q1, MQuadrangle *q2, int count){
+static int SANITY_(GFace *gf,MQuadrangle *q1, MQuadrangle *q2, int count)
+{
for(unsigned int i = 0; i < gf->quadrangles.size(); i++){
MQuadrangle *e1 = gf->quadrangles[i];
MQuadrangle *e2 = q1;
@@ -458,10 +459,7 @@ static int _removeTwoQuadsNodes(GFace *gf)
std::set<MElement*> touched;
std::set<MVertex*> vtouched;
while (it != adj.end()) {
- bool skip=false;
- double surfaceRef=0;
if(it->second.size()==2 && it->first->onWhat()->dim() == 2) {
- const std::vector<MElement*> < = it->second;
MElement *q1 = it->second[0];
MElement *q2 = it->second[1];
if (q1->getNumVertices() == 4 &&
@@ -543,8 +541,8 @@ static bool _isItAGoodIdeaToCollapseThatVertex (GFace *gf,
MVertex *v1,
MVertex *v2,
MVertex *v,
- double FACT){
-
+ double FACT)
+{
double surface_old = surfaceFaceUV(q,gf);
double surface_new = 0;
// double worst_quality_old = q->etaShapeMeasure();
@@ -604,16 +602,14 @@ static bool _isItAGoodIdeaToCollapseThatVertex (GFace *gf,
}
return true;
-
}
-
static bool _isItAGoodIdeaToMoveThatVertex (GFace *gf,
const std::vector<MElement*> &e1,
MVertex *v1,
const SPoint2 &before,
- const SPoint2 &after){
-
+ const SPoint2 &after)
+{
double surface_old = 0;
double surface_new = 0;
@@ -663,8 +659,8 @@ static int _quadWithOneVertexOnBoundary (GFace *gf,
MVertex *v1,
MVertex *v2,
MVertex *v3,
- MVertex *v4){
- //return 0;
+ MVertex *v4)
+{
if (v1->onWhat()->dim() < 2 &&
v2->onWhat()->dim() == 2 &&
v3->onWhat()->dim() == 2 &&
@@ -710,21 +706,12 @@ static int _quadWithOneVertexOnBoundary (GFace *gf,
return 0;
}
-static int _countCommon(std::vector<MElement*> &a, std::vector<MElement*> &b) {
- int count = 0;
- for (unsigned int i=0;i<a.size();i++){
- for (unsigned int j=0;j<b.size();j++){
- if (a[i]==b[j])count++;
- }
- }
- return count;
-}
-
-//see paper from Bunin
-//Guy Bunin (2006) Non-Local Topological Cleanup ,15th International Meshing Roundtable.
-//this is our interpretation of the algorithm.
+// see paper from Bunin Guy Bunin (2006) Non-Local Topological Cleanup ,15th
+// International Meshing Roundtable. This is our interpretation of the
+// algorithm.
static std::vector<MVertex*> closestPathBetweenTwoDefects (v2t_cont &adj,
- v2t_cont :: iterator it) {
+ v2t_cont :: iterator it)
+{
// get neighboring vertices
std::stack<std::vector<MVertex*> > paths;
@@ -738,8 +725,8 @@ static std::vector<MVertex*> closestPathBetweenTwoDefects (v2t_cont &adj,
paths.pop();
it = adj.find (path[path.size() - 1]);
std::set<MVertex *> neigh;
- for (int i = 0; i < it->second.size(); i++){
- for (int j=0;j<4;j++){
+ for (unsigned int i = 0; i < it->second.size(); i++){
+ for (int j = 0; j < 4; j++){
MEdge e = it->second[i]->getEdge(j);
if (e.getVertex(0) == it->first) neigh.insert(e.getVertex(1));
else if (e.getVertex(1) == it->first) neigh.insert(e.getVertex(0));
@@ -778,7 +765,7 @@ static std::vector<MVertex*> closestPathBetweenTwoDefects (v2t_cont &adj,
static MVertex * createNewVertex (GFace *gf, SPoint2 p){
GPoint gp = gf->point(p);
if (!gp.succeeded()) {
- //printf("**** ARRG new vertex not created \n");
+ //printf("**** ARRG new vertex not created \n");
return NULL;
}
return new MFaceVertex(gp.x(),gp.y(),gp.z(),gf,gp.u(),gp.v());
@@ -820,10 +807,11 @@ void createRegularMesh (GFace *gf,
std::vector<MVertex*> &e34, int sign34,
MVertex *v4, SPoint2 &c4,
std::vector<MVertex*> &e41, int sign41,
- std::vector<MQuadrangle*> &q){
+ std::vector<MQuadrangle*> &q)
+{
int N = e12.size();
int M = e23.size();
-
+
//create points using transfinite interpolation
std::vector<std::vector<MVertex*> > tab(M+2);
@@ -848,7 +836,7 @@ void createRegularMesh (GFace *gf,
double v = (double) (j+1)/ ((double)(M+1));
MVertex *v12 = (sign12 >0) ? e12[i] : e12 [N-1-i];
MVertex *v23 = (sign23 >0) ? e23[j] : e23 [M-1-j];
- MVertex *v34 = (sign34 <0) ? e34[i] : e34 [N-1-i];
+ MVertex *v34 = (sign34 <0) ? e34[i] : e34 [N-1-i];
MVertex *v41 = (sign41 <0) ? e41[j] : e41 [M-1-j];
SPoint2 p12; reparamMeshVertexOnFace(v12, gf, p12);
SPoint2 p23; reparamMeshVertexOnFace(v23, gf, p23);
@@ -882,9 +870,9 @@ void createRegularMesh (GFace *gf,
void updateQuadCavity (GFace *gf,
v2t_cont &adj,
std::vector<MElement*> &oldq,
- std::vector<MQuadrangle*> &newq){
-
- for (int i=0;i<oldq.size();i++){
+ std::vector<MQuadrangle*> &newq)
+{
+ for (unsigned int i = 0; i < oldq.size(); i++){
gf->quadrangles.erase(std::remove(gf->quadrangles.begin(),
gf->quadrangles.end(),oldq[i]));
for (int j=0;j<4;j++){
@@ -895,7 +883,7 @@ void updateQuadCavity (GFace *gf,
}
// delete oldq[i];
}
- for (int i=0;i<newq.size();i++){
+ for (unsigned int i = 0; i < newq.size(); i++){
gf->quadrangles.push_back(newq[i]);
for (int j=0;j<4;j++){
MVertex *v = newq[i]->getVertex(j);
@@ -913,7 +901,6 @@ void updateQuadCavity (GFace *gf,
}
}
-
struct quadBlob {
GFace *gf;
int maxCavitySize;
@@ -926,19 +913,22 @@ struct quadBlob {
static fullMatrix<double> M3 ;
static fullMatrix<double> M5 ;
- quadBlob(v2t_cont &a, std::vector<MVertex*> & path, GFace *g, int maxC) : adj(a),gf(g),maxCavitySize(maxC),iBlob(0)
+ quadBlob(v2t_cont &a, std::vector<MVertex*> & path, GFace *g, int maxC)
+ : gf(g),maxCavitySize(maxC),iBlob(0),adj(a)
{
expand_blob(path);
}
- inline bool inBlob(MElement* e) const {
+ inline bool inBlob(MElement* e) const
+ {
return std::find(quads.begin(), quads.end(), e) != quads.end();
}
- inline bool inBoundary(MVertex *v, std::vector<MVertex*> &b) const {
+ inline bool inBoundary(MVertex *v, std::vector<MVertex*> &b) const
+ {
return std::find(b.begin(), b.end(), v) != b.end();
}
- inline void setBlobNumber(int number) {iBlob = number;}
-
- static void computeMatrices(){
+ inline void setBlobNumber(int number) { iBlob = number; }
+ static void computeMatrices()
+ {
if (matricesDone)return;
matricesDone = true;
M3.resize(6,6);
@@ -960,33 +950,33 @@ struct quadBlob {
}
M5.invertInPlace();
}
-
- void expand_blob (std::vector<MVertex*> & path) {
+ void expand_blob (std::vector<MVertex*> & path)
+ {
std::set<MElement*> e;
e.insert(quads.begin(),quads.end());
- for (int i=0;i<path.size();i++){
+ for (unsigned int i = 0; i < path.size(); i++){
v2t_cont :: const_iterator it = adj.find(path[i]);
e.insert(it->second.begin(),it->second.end());
}
quads.clear();
quads.insert(quads.begin(),e.begin(),e.end());
std::set<MVertex*> all;
- for (int i=0;i<quads.size();i++)
- for (int j=0;j<4;j++)
+ for (unsigned int i = 0; i < quads.size(); i++)
+ for (int j = 0; j < 4; j++)
all.insert(quads[i]->getVertex(j));
nodes.clear();
nodes.insert(nodes.begin(),all.begin(),all.end());
bnodes.clear();
}
-
- void buildBoundaryNodes(){
+ void buildBoundaryNodes()
+ {
bnodes.clear();
- for (int i=0;i<nodes.size();i++){
+ for (unsigned int i = 0; i < nodes.size(); i++){
v2t_cont :: const_iterator it = adj.find(nodes[i]);
if (it->first->onWhat()->dim() < 2) bnodes.push_back(nodes[i]);
else {
bool found = false;
- for (int j=0;j<it->second.size();j++){
+ for (unsigned int j = 0; j < it->second.size(); j++){
if (!inBlob(it->second[j])){
found = true;
}
@@ -997,28 +987,29 @@ struct quadBlob {
}
}
}
-
- int topologicalAngle(MVertex*v) const {
+ int topologicalAngle(MVertex*v) const
+ {
int typical = 0;
v2t_cont :: const_iterator it = adj.find(v);
int outside = 0;
- for (int j=0;j<it->second.size();j++)
- if (!inBlob(it->second[j]))outside++;
+ for (unsigned int j = 0; j < it->second.size(); j++)
+ if (!inBlob(it->second[j])) outside++;
- if (v->onWhat()->dim() == 1)typical = 2;
- else if (v->onWhat()->dim() == 0)typical = 3;
+ if (v->onWhat()->dim() == 1) typical = 2;
+ else if (v->onWhat()->dim() == 0) typical = 3;
return outside - 2 + typical;
}
- int construct_meshable_blob (int iter) {
+ int construct_meshable_blob (int iter)
+ {
int blobsize = quads.size();
//printf("*** starting with a blob of size %d (ITER=%d)\n",blobsize, iter);
while(1){
- if(quads.size() > maxCavitySize)return 0;
- if(quads.size() >= gf->quadrangles.size())return 0;
+ if((int)quads.size() > maxCavitySize) return 0;
+ if(quads.size() >= gf->quadrangles.size()) return 0;
buildBoundaryNodes();
int topo_convex_region = 1;
- for (int i=0; i<bnodes.size();i++){
+ for (unsigned int i = 0; i < bnodes.size(); i++){
MVertex *v = bnodes[i];
// do not do anything in boundary layers !!!
if (v->onWhat()->dim() == 2){
@@ -1033,19 +1024,18 @@ struct quadBlob {
}
}
if (topo_convex_region){
- if (meshable(iter))return 1;
+ if (meshable(iter)) return 1;
else expand_blob(bnodes);
}
- if (blobsize == quads.size())return 0;
+ if (blobsize == (int)quads.size()) return 0;
blobsize = quads.size();
}// while 1
}
-
- bool orderBNodes () {
-
+ bool orderBNodes ()
+ {
MVertex *v = 0;
std::vector<MVertex*> temp;
- for (int i=0;i<bnodes.size();i++){
+ for (unsigned int i = 0; i < bnodes.size(); i++){
if (topologicalAngle(bnodes[i]) > 0){
v = bnodes[i];
break;
@@ -1058,12 +1048,12 @@ struct quadBlob {
v2t_cont :: const_iterator it = adj.find(v);
int ss = temp.size();
std::vector<MElement*> elems = it->second;
-
+
//EMI: try to orient BNodes the same as quad orientations
- for (int j=0;j<elems.size();j++){
+ for (unsigned int j = 0; j < elems.size(); j++){
bool found =false;
if(inBlob(elems[j])){
- for (int i=0;i<4;i++){
+ for (int i = 0; i < 4; i++){
MVertex *v0 = elems[j]->getVertex(i%4);
MVertex *v1 = elems[j]->getVertex((i+1)%4);
if (v0 == v && inBoundary(v1,bnodes) && !inBoundary(v1,temp)) {
@@ -1076,7 +1066,7 @@ struct quadBlob {
}
if (found) break;
}
-
+
//JF this does not orient quads
// for (int j=0;j<elems.size();j++){
// bool found = false;
@@ -1104,11 +1094,11 @@ struct quadBlob {
// if (found)break;
// }
- if (temp.size() == ss)return false;
- if (temp.size() == bnodes.size())break;
+ if ((int)temp.size() == ss) return false;
+ if (temp.size() == bnodes.size()) break;
}
- //change orientation
+ //change orientation
// bool oriented = false;
// MVertex *vB1 = temp[0];
// MVertex *vB2 = temp[1];
@@ -1121,29 +1111,28 @@ struct quadBlob {
// }
// if (oriented) break;
// }
- // std::vector<MVertex*> temp2;
- // temp2.push_back(temp[0]);
- // if (!oriented) {
- // std::reverse(temp.begin(), temp.end());
- // for (int i = 0; i< temp.size()-1; i++) temp2.push_back(temp[i]);
- // temp = temp2;
- // }
-
+ // std::vector<MVertex*> temp2;
+ // temp2.push_back(temp[0]);
+ // if (!oriented) {
+ // std::reverse(temp.begin(), temp.end());
+ // for (int i = 0; i< temp.size()-1; i++) temp2.push_back(temp[i]);
+ // temp = temp2;
+ // }
+
bnodes = temp;
return true;
}
-
- void printBlob(int iter, int method){
-
+ void printBlob(int iter, int method)
+ {
if (!CTX::instance()->mesh.saveAll) return;
char name[234];
sprintf(name,"blob%d_%d_%d.pos", iBlob, iter, method);
FILE *f = fopen(name,"w");
fprintf(f,"View \"%s\" {\n",name);
- for (int i=0;i<quads.size();i++){
+ for (unsigned int i = 0; i < quads.size(); i++){
quads[i]->writePOS(f,true,false,false,false,false,false);
}
- for (int i=0;i<bnodes.size();i++){
+ for (unsigned int i = 0; i < bnodes.size(); i++){
fprintf(f,"SP(%g,%g,%g) {%d};\n",
bnodes[i]->x(),
bnodes[i]->y(),
@@ -1151,16 +1140,13 @@ struct quadBlob {
}
fprintf(f,"};\n");
fclose(f);
-
}
-
void smooth (int);
-
- bool meshable (int iter) {
-
+ bool meshable (int iter)
+ {
int ncorners = 0;
MVertex *corners[5];
- for (int i=0;i<bnodes.size();i++){
+ for (unsigned int i = 0; i < bnodes.size(); i++){
if (topologicalAngle(bnodes[i]) > 0) ncorners ++;
if (ncorners > 5)return false;
}
@@ -1173,7 +1159,7 @@ struct quadBlob {
if (!orderBNodes () )return false;
int side = -1;
int count[5] = {0,0,0,0,0};
- for (int i=0;i<bnodes.size();i++){
+ for (unsigned int i = 0; i < bnodes.size(); i++){
if (topologicalAngle(bnodes[i]) > 0){
side++;
corners[side] = (bnodes[i]);
@@ -1182,7 +1168,7 @@ struct quadBlob {
}
if (ncorners == 3){
-
+
fullVector<double> rhs(6);
fullVector<double> result(6);
rhs(3) = count[0]+1.;
@@ -1240,7 +1226,7 @@ struct quadBlob {
e012_20, -1,
q);
- createRegularMesh (gf,
+ createRegularMesh (gf,
v012,p012,
e012_01, +1,
v01,p01,
@@ -1278,7 +1264,7 @@ struct quadBlob {
for (int i=0;i<a2-1;i++)e1_2.push_back(bnodes[i+1 + a1]);
for (int i=0;i<a1-1;i++)e2_3.push_back(bnodes[i+1 + a1 + a2]);
for (int i=0;i<a2-1;i++)e3_0.push_back(bnodes[i+1 + a1 + a2 + a1]);
-
+
std::vector<MQuadrangle*> q;
createRegularMesh (gf,
v0,p0,
@@ -1452,7 +1438,7 @@ fullMatrix<double> quadBlob::M5;
static int _defectsRemovalBunin(GFace *gf, int maxCavitySize)
{
-
+
if (maxCavitySize == 0)return 0;
v2t_cont adj;
std::vector<MElement*> c;
@@ -1472,7 +1458,7 @@ static int _defectsRemovalBunin(GFace *gf, int maxCavitySize)
double t2 = Cpu();
double C = (t2-t1);
double A=0,B=0;
- for (int i=0;i<defects.size();i++){
+ for (unsigned int i = 0; i < defects.size(); i++){
it = adj.find(defects[i]);
if (it->first->onWhat()->dim() == 2 && it->second.size() != 4 && it->second.size() != 0) {
double t3 = Cpu();
@@ -1504,7 +1490,6 @@ static int _defectsRemovalBunin(GFace *gf, int maxCavitySize)
// quad
static int _splitFlatQuads(GFace *gf, double minQual)
{
-
v2t_cont adj;
buildVertexToElement(gf->triangles,adj);
buildVertexToElement(gf->quadrangles,adj);
@@ -1513,8 +1498,6 @@ static int _splitFlatQuads(GFace *gf, double minQual)
std::vector<MQuadrangle*> added_quadrangles;
std::set<MElement*> deleted_quadrangles;
while (it != adj.end()) {
- bool skip=false;
- double surfaceRef=0;
// split that guy if too bad
if(it->second.size()==1 && it->first->onWhat()->dim() == 1) {
const std::vector<MElement*> < = it->second;
@@ -1567,7 +1550,7 @@ static int _splitFlatQuads(GFace *gf, double minQual)
}
++it;
}
- for (int i=0;i<gf->quadrangles.size();i++){
+ for (unsigned int i = 0; i < gf->quadrangles.size(); i++){
if (deleted_quadrangles.find(gf->quadrangles[i]) == deleted_quadrangles.end()){
added_quadrangles.push_back(gf->quadrangles[i]);
}
@@ -1609,16 +1592,21 @@ static int _removeDiamonds(GFace *gf)
touched.find(v2) == touched.end() &&
touched.find(v3) == touched.end() &&
touched.find(v4) == touched.end() ) {
- if (_quadWithOneVertexOnBoundary (gf,touched,diamonds,deleted,it2->second,it4->second,it1->second,q,v1,v2,v3,v4)){}
- else if (_quadWithOneVertexOnBoundary (gf,touched,diamonds,deleted,it3->second,it1->second,it2->second,q,v2,v3,v4,v1)){}
- else if (_quadWithOneVertexOnBoundary (gf,touched,diamonds,deleted,it4->second,it2->second,it3->second,q,v3,v4,v1,v2)){}
- else if (_quadWithOneVertexOnBoundary (gf,touched,diamonds,deleted,it1->second,it3->second,it4->second,q,v4,v1,v2,v3)){}
+ if (_quadWithOneVertexOnBoundary (gf,touched,diamonds,deleted,it2->second,
+ it4->second,it1->second,q,v1,v2,v3,v4)){}
+ else if (_quadWithOneVertexOnBoundary (gf,touched,diamonds,deleted,it3->second,
+ it1->second,it2->second,q,v2,v3,v4,v1)){}
+ else if (_quadWithOneVertexOnBoundary (gf,touched,diamonds,deleted,it4->second,
+ it2->second,it3->second,q,v3,v4,v1,v2)){}
+ else if (_quadWithOneVertexOnBoundary (gf,touched,diamonds,deleted,it1->second,
+ it3->second,it4->second,q,v4,v1,v2,v3)){}
else if (v1->onWhat()->dim() == 2 &&
v2->onWhat()->dim() == 2 &&
v3->onWhat()->dim() == 2 &&
v4->onWhat()->dim() == 2 &&
it1->second.size() ==3 && it3->second.size() == 3 &&
- _isItAGoodIdeaToCollapseThatVertex (gf, it1->second, it3->second, q, v1, v3, v3,10.)){
+ _isItAGoodIdeaToCollapseThatVertex (gf, it1->second, it3->second,
+ q, v1, v3, v3,10.)){
touched.insert(v1);
touched.insert(v2);
touched.insert(v3);
@@ -1637,7 +1625,8 @@ static int _removeDiamonds(GFace *gf)
v1->onWhat()->dim() == 2 &&
v3->onWhat()->dim() == 2 &&
it1->second.size() ==3 && it3->second.size() == 3 &&
- _isItAGoodIdeaToCollapseThatVertex (gf, it1->second, it3->second, q, v1, v3, v1,10.)){
+ _isItAGoodIdeaToCollapseThatVertex (gf, it1->second, it3->second,
+ q, v1, v3, v1,10.)){
touched.insert(v1);
touched.insert(v2);
touched.insert(v3);
@@ -1656,7 +1645,8 @@ static int _removeDiamonds(GFace *gf)
v2->onWhat()->dim() == 2 &&
v4->onWhat()->dim() == 2 &&
it2->second.size() == 3 && it4->second.size() == 3 &&
- _isItAGoodIdeaToCollapseThatVertex (gf, it2->second, it4->second, q, v2, v4, v4,10.)){
+ _isItAGoodIdeaToCollapseThatVertex (gf, it2->second, it4->second,
+ q, v2, v4, v4,10.)){
touched.insert(v1);
touched.insert(v2);
touched.insert(v3);
@@ -1675,7 +1665,8 @@ static int _removeDiamonds(GFace *gf)
v2->onWhat()->dim() == 2 &&
v4->onWhat()->dim() == 2 &&
it2->second.size() == 3 && it4->second.size() == 3 &&
- _isItAGoodIdeaToCollapseThatVertex (gf, it2->second, it4->second, q, v2, v4, v2,10.)){
+ _isItAGoodIdeaToCollapseThatVertex (gf, it2->second, it4->second,
+ q, v2, v4, v2,10.)){
touched.insert(v1);
touched.insert(v2);
touched.insert(v3);
@@ -1714,7 +1705,8 @@ static int _removeDiamonds(GFace *gf)
return diamonds.size();
}
-int removeDiamonds(GFace *gf){
+int removeDiamonds(GFace *gf)
+{
int nbRemove = 0;
while(1){
int x = _removeDiamonds(gf);
@@ -1725,7 +1717,8 @@ int removeDiamonds(GFace *gf){
return nbRemove;
}
-int splitFlatQuads(GFace *gf, double q){
+int splitFlatQuads(GFace *gf, double q)
+{
int nbRemove = 0;
while(1){
int x = _splitFlatQuads(gf,q);
@@ -1748,14 +1741,17 @@ struct opti_data_vertex_relocation {
MVertex *v;
GFace *gf;
double minq,eps;
- opti_data_vertex_relocation (const std::vector<MElement*> & _e, MVertex * _v, GFace *_gf) : e(_e), v(_v), gf(_gf)
+ opti_data_vertex_relocation (const std::vector<MElement*> & _e,
+ MVertex * _v, GFace *_gf)
+ : e(_e), v(_v), gf(_gf)
{
minq = -2;
eps = minq / (1. - minq);
}
- double f() const {
+ double f() const
+ {
double val = 0.0;
- for (int i=0;i < e.size();i++){
+ for (unsigned int i = 0; i < e.size(); i++){
const double q = (e[i]->getNumVertices() == 4) ? e[i]->etaShapeMeasure() :
e[i]->gammaShapeMeasure();
const double dd = (1 + eps) * q - eps;
@@ -1765,7 +1761,8 @@ struct opti_data_vertex_relocation {
}
return val;
}
- void df(const double &F, const double &U, const double &V, double &dfdu, double &dfdv) {
+ void df(const double &F, const double &U, const double &V, double &dfdu, double &dfdv)
+ {
const double eps = 1.e-6;
set_(U+eps,V);
const double FU = f();
@@ -1775,7 +1772,8 @@ struct opti_data_vertex_relocation {
dfdu = -(F-FU)/eps;
dfdv = -(F-FV)/eps;
}
- void set_(double U, double V){
+ void set_(double U, double V)
+ {
GPoint gp = gf->point(U,V);
if (!gp.succeeded()) Msg::Error("point not OK \n");
v->x() = gp.x();
@@ -1786,7 +1784,9 @@ struct opti_data_vertex_relocation {
}
};
-void bfgs_callback_vertex_relocation (const alglib::real_1d_array& x, double& func, alglib::real_1d_array& grad, void* ptr) {
+void bfgs_callback_vertex_relocation (const alglib::real_1d_array& x,
+ double& func, alglib::real_1d_array& grad, void* ptr)
+{
opti_data_vertex_relocation* w = static_cast<opti_data_vertex_relocation*>(ptr);
w->set_(x[0],x[1]);
func = w->f();
@@ -1827,12 +1827,12 @@ static void _relocateVertexOpti(GFace *gf, MVertex *ver,
minlbfgsresults(state,x,rep);
double minq = 0;
- for (int i=0;i < data.e.size();i++){
+ for (unsigned int i = 0; i < data.e.size(); i++){
const double q = (data.e[i]->getNumVertices() == 4) ? data.e[i]->etaShapeMeasure() :
data.e[i]->gammaShapeMeasure();
minq = std::min(q,minq);
}
- if (minq < 0.01)data.set_(U,V);
+ if (minq < 0.01) data.set_(U,V);
else printf("OKIDOKI\n");
// if (rep.terminationtype != 4){
// data.set_(U,V);
@@ -1908,9 +1908,10 @@ static void _relocateVertex(GFace *gf, MVertex *ver,
}
}
-void quadBlob::smooth (int niter) {
- for (int i=0;i<niter;i++){
- for (int j=0;j<nodes.size();j++){
+void quadBlob::smooth (int niter)
+{
+ for (int i = 0; i < niter; i++){
+ for (unsigned int j = 0; j < nodes.size(); j++){
v2t_cont::iterator it = adj.find(nodes[j]);
_relocateVertex(gf, it->first, it->second);
}
@@ -1980,7 +1981,7 @@ int _edgeSwapQuadsForBetterQuality(GFace *gf)
return 0;
}
- MQuadrangle *q1A, *q2A, *q1B, *q2B;
+ MQuadrangle *q1A, *q2A, *q1B, *q2B;
if (rot1 > 0.0 && rot2 > 0.0){
q1A = new MQuadrangle (v11,v22,v2,v12);
q2A = new MQuadrangle (v22,v11,v1,v21);
@@ -1996,7 +1997,7 @@ int _edgeSwapQuadsForBetterQuality(GFace *gf)
double worst_quality_A = std::min(q1A->etaShapeMeasure(),q2A->etaShapeMeasure());
double worst_quality_B = std::min(q1B->etaShapeMeasure(),q2B->etaShapeMeasure());
- bool c1,c2,ca1,ca2,cb1,cb2;
+ bool ca1,ca2,cb1,cb2;
double old_surface = surfaceFaceUV(e1,gf) + surfaceFaceUV(e2,gf) ;
double new_surface_A = surfaceFaceUV(q1A,gf,&ca1) + surfaceFaceUV(q2A,gf,&ca2) ;
double new_surface_B = surfaceFaceUV(q1B,gf,&cb1) + surfaceFaceUV(q2B,gf,&cb2) ;
@@ -2052,7 +2053,8 @@ int _edgeSwapQuadsForBetterQuality(GFace *gf)
return COUNT;
}
-static int edgeSwapQuadsForBetterQuality ( GFace *gf ) {
+static int edgeSwapQuadsForBetterQuality (GFace *gf)
+{
int COUNT = 0;
while(1){
int k = _edgeSwapQuadsForBetterQuality (gf);
@@ -2063,12 +2065,13 @@ static int edgeSwapQuadsForBetterQuality ( GFace *gf ) {
return COUNT;
}
-static std::vector<MVertex*> computeBoundingPoints (const std::vector<MElement*> <,MVertex *v = 0){
-
+static std::vector<MVertex*> computeBoundingPoints (const std::vector<MElement*> <,
+ MVertex *v = 0)
+{
// get boundary edges
std::set<MEdge,Less_Edge> edges;
- for (int i=0;i<lt.size();i++){
- for (int j=0;j<lt[i]->getNumEdges();j++){
+ for (unsigned int i = 0; i < lt.size(); i++){
+ for (int j = 0; j < lt[i]->getNumEdges(); j++){
std::set<MEdge,Less_Edge>::iterator it = edges.find(lt[i]->getEdge(j));
if (it == edges.end())
edges.insert(lt[i]->getEdge(j));
@@ -2127,13 +2130,11 @@ static std::vector<MVertex*> computeBoundingPoints (const std::vector<MElement*>
return result;
}
-
-int postProcessExtraEdges (GFace *gf, std::vector<std::pair<MElement*,MElement*> > &toProcess){
-
+int postProcessExtraEdges (GFace *gf, std::vector<std::pair<MElement*,MElement*> > &toProcess)
+{
// some pairs of elements that should be recombined are not neighbor elements (see our paper)
// so we recombine them in another way (adding a new node)
-
Msg::Debug("%d extra edges to be processed",(int)toProcess.size());
// return 1;
@@ -2144,7 +2145,7 @@ int postProcessExtraEdges (GFace *gf, std::vector<std::pair<MElement*,MElement*>
std::set<MElement*>deleted;
- for (int i=0; i<toProcess.size() ; i++){
+ for (unsigned int i = 0; i < toProcess.size(); i++){
MElement *t1 = toProcess[i].first;
MElement *t2 = toProcess[i].second;
MVertex *common = 0;
@@ -2177,7 +2178,7 @@ int postProcessExtraEdges (GFace *gf, std::vector<std::pair<MElement*,MElement*>
// create a new vertex
else {
SPoint2 p;
- bool success = reparamMeshVertexOnFace(it->first,gf, p);
+ reparamMeshVertexOnFace(it->first,gf, p);
MFaceVertex *newVertex = new MFaceVertex (it->first->x(),
it->first->y(),
it->first->z(),
@@ -2231,10 +2232,10 @@ int postProcessExtraEdges (GFace *gf, std::vector<std::pair<MElement*,MElement*>
std::vector<MElement*> newAdj;
newAdj.push_back(q1);
newAdj.push_back(q2);
- for (int k=0;k<it->second.size();k++){
+ for (unsigned int k = 0; k < it->second.size(); k++){
if (it->second[k]->getNumVertices() == 4){
newAdj.push_back(it->second[k]);
- for (int vv=0;vv<4;vv++){
+ for (int vv = 0; vv < 4; vv++){
if (it->first == it->second[k]->getVertex(vv))
it->second[k]->setVertex(vv,newVertex);
}
@@ -2451,8 +2452,6 @@ int edgeSwapPass(GFace *gf, std::set<MTri3*, compareTri3Ptr> &allTris,
return nbSwapTot;
}
-
-
inline double computeEdgeAdimLength(MVertex *v1, MVertex *v2, GFace *f,
const std::vector<double> &Us,
const std::vector<double> &Vs,
@@ -2690,7 +2689,8 @@ static int _recombineIntoQuads(GFace *gf, int recur_level, bool cubicGraph = 1)
it->second.first->getNumVertices() == 3){
for (int i=0;i<2;i++){
MVertex *v = it->first.getVertex(i);
- std::map<MVertex*,std::pair<MElement*,MElement*> > :: iterator itv = makeGraphPeriodic.find(v);
+ std::map<MVertex*,std::pair<MElement*,MElement*> > :: iterator itv =
+ makeGraphPeriodic.find(v);
if (itv == makeGraphPeriodic.end()){
makeGraphPeriodic[v] = std::make_pair(it->second.first,(MElement*)0);
}
@@ -2725,7 +2725,7 @@ static int _recombineIntoQuads(GFace *gf, int recur_level, bool cubicGraph = 1)
int *elist = new int [2*ecount];
int *elen = new int [ecount];
- for (int i=0;i<pairs.size();++i){
+ for (unsigned int i = 0; i < pairs.size(); ++i){
elist[2*i] = t2n[pairs[i].t1];
elist[2*i+1] = t2n[pairs[i].t2];
//elen [i] = (int) pairs[i].angle;
@@ -2733,9 +2733,9 @@ static int _recombineIntoQuads(GFace *gf, int recur_level, bool cubicGraph = 1)
double angle = atan2(pairs[i].n1->y()-pairs[i].n2->y(),
pairs[i].n1->x()-pairs[i].n2->x());
- double x = .5*(pairs[i].n1->x()+pairs[i].n2->x());
- double y = .5*(pairs[i].n1->y()+pairs[i].n2->y());
- double opti = atan2(y,x);
+ //double x = .5*(pairs[i].n1->x()+pairs[i].n2->x());
+ //double y = .5*(pairs[i].n1->y()+pairs[i].n2->y());
+ //double opti = atan2(y,x);
//angle -= opti;
while (angle < 0 || angle > M_PI/2){
if (angle < 0) angle += M_PI/2;
@@ -2811,8 +2811,9 @@ static int _recombineIntoQuads(GFace *gf, int recur_level, bool cubicGraph = 1)
}
}
if (removed.size() == 0) return _recombineIntoQuads(gf, 11);
- for (int i=0;i<gf->triangles.size();++i){
- if (removed.find(gf->triangles[i]) == removed.end())triangles2.push_back(gf->triangles[i]);
+ for (unsigned int i = 0; i < gf->triangles.size(); ++i){
+ if (removed.find(gf->triangles[i]) ==
+ removed.end())triangles2.push_back(gf->triangles[i]);
else delete gf->triangles[i];
}
gf->triangles=triangles2;
@@ -2841,7 +2842,6 @@ static int _recombineIntoQuads(GFace *gf, int recur_level, bool cubicGraph = 1)
MElement *t2 = n2t[i2];
touched.insert(t1);
touched.insert(t2);
- int orientation = 0;
MVertex *other = 0;
for(int i = 0; i < 3; i++) {
if (t1->getVertex(0) != t2->getVertex(i) &&
@@ -2871,15 +2871,17 @@ static int _recombineIntoQuads(GFace *gf, int recur_level, bool cubicGraph = 1)
free(elist);
pairs.clear();
if (recur_level == 0)
- Msg::Debug("Perfect Match Succeeded in Quadrangulation (%g sec)",matzeit);
+ Msg::Debug("Perfect Match Succeeded in Quadrangulation (%g sec)", matzeit);
else
- Msg::Info(" :-) Perfect Match Succeeded in Quadrangulation after Splits (%g sec)",matzeit);
+ Msg::Info(" :-) Perfect Match Succeeded in Quadrangulation after Splits (%g sec)",
+ matzeit);
}
}
#else
- Msg::Warning("Gmsh should be compiled with the Blossom IV code and CONCORDE in order to allow the Blossom optimization");
+ Msg::Warning("Gmsh should be compiled with the Blossom IV code and CONCORDE "
+ "in order to allow the Blossom optimization");
#endif
}
@@ -2971,31 +2973,43 @@ void recombineIntoQuads(GFace *gf,
double t6 = Cpu();
int maxCavitySize = CTX::instance()->mesh.bunin;
optistatus[4] = _defectsRemovalBunin(gf,maxCavitySize);
- if(optistatus[4] && saveAll){ sprintf(NAME,"iter%dB.msh",COUNT++); gf->model()->writeMSH(NAME); }
+ if(optistatus[4] && saveAll){
+ sprintf(NAME,"iter%dB.msh",COUNT++); gf->model()->writeMSH(NAME);
+ }
double t7 = Cpu();
double t1 = Cpu();
optistatus[0] = splitFlatQuads(gf, .3);
- if(optistatus[0] && saveAll){ sprintf(NAME,"iter%dSF.msh",COUNT++); gf->model()->writeMSH(NAME);}
+ if(optistatus[0] && saveAll){
+ sprintf(NAME,"iter%dSF.msh",COUNT++); gf->model()->writeMSH(NAME);
+ }
double t2 = Cpu();
optistatus[1] = removeTwoQuadsNodes(gf);
- if(optistatus[1] && saveAll){ sprintf(NAME,"iter%dTQ.msh",COUNT++); gf->model()->writeMSH(NAME);}
+ if(optistatus[1] && saveAll){
+ sprintf(NAME,"iter%dTQ.msh",COUNT++); gf->model()->writeMSH(NAME);
+ }
double t3 = Cpu();
optistatus[2] = removeDiamonds(gf);
- if(optistatus[2] && saveAll){ sprintf(NAME,"iter%dD.msh",COUNT++); gf->model()->writeMSH(NAME); }
+ if(optistatus[2] && saveAll){
+ sprintf(NAME,"iter%dD.msh",COUNT++); gf->model()->writeMSH(NAME);
+ }
+
double t4 = Cpu();
laplaceSmoothing(gf,CTX::instance()->mesh.nbSmoothing);
double t5 = Cpu();
optistatus[3] = edgeSwapQuadsForBetterQuality(gf);
- if(optistatus[3] && saveAll){ sprintf(NAME,"iter%dS.msh",COUNT++); gf->model()->writeMSH(NAME); }
+ if(optistatus[3] && saveAll){
+ sprintf(NAME,"iter%dS.msh",COUNT++); gf->model()->writeMSH(NAME);
+ }
tFQ += (t2-t1);
tTQ += (t3-t2);
tDI += (t4-t3);
tLA += (t5-t4);
tSW += (t6-t5);
tBU += (t7-t6);
- if (!(optistatus[0]+optistatus[1]+optistatus[2]+optistatus[3]+optistatus[4])) break;
+ if (!(optistatus[0]+optistatus[1]+optistatus[2]+optistatus[3]+optistatus[4]))
+ break;
bool theSame = true;
for (int i=0;i<5;i++){
if (optistatus[i] != oldoptistatus[i])theSame = false;
@@ -3004,9 +3018,10 @@ void recombineIntoQuads(GFace *gf,
if(theSame)break;
if (ITER++ >= 20)break;
}
- Msg::Debug("Opti times FQ(%4.3f) tQ(%4.3f) DI(%4.3f) LA(%4.3f) SW(%4.3f) BUN(%4.3f)",tFQ,tTQ,tDI,tLA,tSW,tBU);
+ Msg::Debug("Opti times FQ(%4.3f) tQ(%4.3f) DI(%4.3f) LA(%4.3f) "
+ "SW(%4.3f) BUN(%4.3f)",tFQ,tTQ,tDI,tLA,tSW,tBU);
}
-
+
if(haveParam) laplaceSmoothing(gf,CTX::instance()->mesh.nbSmoothing);
}
@@ -3029,9 +3044,8 @@ void recombineIntoQuads(GFace *gf,
void quadsToTriangles(GFace *gf, double minqual)
{
-
std::vector<MQuadrangle*> qds;
- for (int i=0;i<gf->quadrangles.size();i++){
+ for (unsigned int i = 0; i < gf->quadrangles.size(); i++){
MQuadrangle *q = gf->quadrangles[i];
if (q->etaShapeMeasure() < minqual){
MTriangle *t11 = new MTriangle (q->getVertex(0),q->getVertex(1),q->getVertex(2));
@@ -3185,24 +3199,21 @@ void Temporary::read_data(std::string file_name)
int i,j,number;
double x,y,z;
MElement*element;
- PView*view;
PViewData*data;
PView::readMSH(file_name,-1);
- std::vector<PView*> list = PView::list;
- data = list[0]->getData();
+ data = PView::list[0]->getData();
for(i = 0;i < data->getNumEntities(0);i++)
{
if(data->skipEntity(0,i)) continue;
- for(j = 0;j < data->getNumElements(0,i);j++)
- {
- if(data->skipElement(0,i,j)) continue;
- element = data->getElement(0,i,j);
- number = element->getNum();
- data->getValue(0,i,j,0,x);
- data->getValue(0,i,j,1,y);
- data->getValue(0,i,j,2,z);
- gradients[number] = SVector3(x,y,z);
- }
+ for(j = 0;j < data->getNumElements(0,i);j++){
+ if(data->skipElement(0,i,j)) continue;
+ element = data->getElement(0,i,j);
+ number = element->getNum();
+ data->getValue(0,i,j,0,x);
+ data->getValue(0,i,j,1,y);
+ data->getValue(0,i,j,2,z);
+ gradients[number] = SVector3(x,y,z);
+ }
}
#endif
}
@@ -3220,6 +3231,6 @@ void Temporary::quadrilaterize(std::string file_name,double _w1,double _w2,doubl
for(iterator = model->firstFace();iterator != model->lastFace();iterator++)
{
face = *iterator;
- _recombineIntoQuads(face,1,1);
+ _recombineIntoQuads(face,1,1);
}
}
diff --git a/Mesh/meshGRegion.cpp b/Mesh/meshGRegion.cpp
index 64d68cb87b1a7812e3ea4cecf99bf2af09604368..baa11d7055aa242c51e6a2ab2157624492c62814 100644
--- a/Mesh/meshGRegion.cpp
+++ b/Mesh/meshGRegion.cpp
@@ -53,9 +53,10 @@ void printVoronoi(GRegion *gr, std::set<SPoint3> &candidates)
itmap->second = allTets;
}
}
- for (int j=0; j<4; j++){
+ for (int j = 0; j < 4; j++){
MFace f = ele->getFace(j);
- std::map<MFace, std::vector<MTetrahedron*>, Less_Face >::iterator itmap = face2Tet.find(f);
+ std::map<MFace, std::vector<MTetrahedron*>, Less_Face >::iterator itmap =
+ face2Tet.find(f);
if (itmap == face2Tet.end()){
std::vector<MTetrahedron*> oneTet;
oneTet.push_back(ele);
diff --git a/Mesh/meshGRegionMMG3D.cpp b/Mesh/meshGRegionMMG3D.cpp
index ae0f8bec465fd69805cddb376f96f0abf30969a4..643ee81d298c5ae8672c4d9ef7979daf16ebbfd5 100644
--- a/Mesh/meshGRegionMMG3D.cpp
+++ b/Mesh/meshGRegionMMG3D.cpp
@@ -27,7 +27,7 @@ static void MMG2gmsh(GRegion *gr, MMG_pMesh mmg, std::map<int,MVertex*> &mmg2gms
{
std::map<int,MVertex*> kToMVertex;
for (int k=1;k<= mmg->np ; k++){
- MMG_pPoint ppt = &mmg->point[k];
+ MMG_pPoint ppt = &mmg->point[k];
if (ppt->tag & M_UNUSED) continue;
std::map<int,MVertex*>::iterator it = mmg2gmsh.find(ppt->ref);
if (it == mmg2gmsh.end()){
@@ -36,10 +36,10 @@ static void MMG2gmsh(GRegion *gr, MMG_pMesh mmg, std::map<int,MVertex*> &mmg2gms
kToMVertex[k] = v;
}
else kToMVertex[k] = it->second;
- }
-
- for (int k=1; k<=mmg->ne; k++) {
- MMG_pTetra ptetra = &mmg->tetra[k];
+ }
+
+ for (int k=1; k<=mmg->ne; k++) {
+ MMG_pTetra ptetra = &mmg->tetra[k];
if ( ptetra->v[0] ){
MVertex *v1 = kToMVertex[ptetra->v[0]];
MVertex *v2 = kToMVertex[ptetra->v[1]];
@@ -54,26 +54,26 @@ static void MMG2gmsh(GRegion *gr, MMG_pMesh mmg, std::map<int,MVertex*> &mmg2gms
}
}
-static void gmsh2MMG(GRegion *gr, MMG_pMesh mmg, MMG_pSol sol,
+static void gmsh2MMG(GRegion *gr, MMG_pMesh mmg, MMG_pSol sol,
std::map<int,MVertex*> &mmg2gmsh)
{
mmg->ne = gr->tetrahedra.size();
std::set<MVertex*> allVertices;
- for (int i=0;i< gr->tetrahedra.size() ; i++){
+ for (unsigned int i = 0; i < gr->tetrahedra.size(); i++){
allVertices.insert(gr->tetrahedra[i]->getVertex(0));
allVertices.insert(gr->tetrahedra[i]->getVertex(1));
allVertices.insert(gr->tetrahedra[i]->getVertex(2));
allVertices.insert(gr->tetrahedra[i]->getVertex(3));
}
mmg->np = sol->np = allVertices.size();
-
+
std::list<GFace*> f = gr->faces();
mmg->nt = 0;
for (std::list<GFace*>::iterator it = f.begin(); it != f.end() ; ++it){
mmg->nt += (*it)->triangles.size();
}
-
+
mmg->npmax = sol->npmax = 1000000;
mmg->ntmax = 700000;
mmg->nemax = 7000000;
@@ -85,15 +85,15 @@ static void gmsh2MMG(GRegion *gr, MMG_pMesh mmg, MMG_pSol sol,
mmg->adja = (int*)calloc(4*mmg->nemax+5,sizeof(int));
sol->offset = 6;
- sol->met = (double*)calloc(sol->npmax+1,sol->offset*sizeof(double));
- sol->metold = (double*)calloc(sol->npmax+1,sol->offset*sizeof(double));
+ sol->met = (double*)calloc(sol->npmax+1,sol->offset*sizeof(double));
+ sol->metold = (double*)calloc(sol->npmax+1,sol->offset*sizeof(double));
std::map<MVertex*,std::pair<double,int> > LCS;
for (std::list<GFace*>::iterator it = f.begin(); it != f.end() ; ++it){
- for (int i=0;i<(*it)->triangles.size();i++){
+ for (unsigned int i = 0; i < (*it)->triangles.size(); i++){
MTriangle *t = (*it)->triangles[i];
double L = t->maxEdge();
- for (int k=0;k<3;k++){
+ for (int k = 0; k < 3; k++){
MVertex *v = t->getVertex(k);
std::map<MVertex*,std::pair<double,int> >::iterator itv = LCS.find(v);
if (itv != LCS.end()){
@@ -108,13 +108,13 @@ static void gmsh2MMG(GRegion *gr, MMG_pMesh mmg, MMG_pSol sol,
}
//printf("%d vertices %d on faces\n", (int) allVertices.size(), (int) LCS.size());
-
+
int k=1;
int count = 1;//sol->offset;
std::map<int,int> gmsh2mmg_num;
for (std::set<MVertex*>::iterator it = allVertices.begin();
it != allVertices.end(); ++it){
- MMG_pPoint ppt = &mmg->point[k];
+ MMG_pPoint ppt = &mmg->point[k];
ppt->c[0] = (*it)->x();
ppt->c[1] = (*it)->y();
@@ -123,7 +123,7 @@ static void gmsh2MMG(GRegion *gr, MMG_pMesh mmg, MMG_pSol sol,
gmsh2mmg_num[(*it)->getNum()] = k;
MVertex *v = *it;
- double U=0,V=0;
+ double U = 0, V = 0;
if (v->onWhat()->dim() == 1){
v->getParameter(0,U);
}
@@ -131,8 +131,8 @@ static void gmsh2MMG(GRegion *gr, MMG_pMesh mmg, MMG_pSol sol,
v->getParameter(0,U);
v->getParameter(1,V);
}
- //double lc = BGM_MeshSize(v->onWhat(), U,V,v->x(), v->y(), v->z());
- SMetric3 m = BGM_MeshMetric(v->onWhat(), U,V,v->x(), v->y(), v->z());
+ //double lc = BGM_MeshSize(v->onWhat(), U,V,v->x(), v->y(), v->z());
+ SMetric3 m = BGM_MeshMetric(v->onWhat(), U,V,v->x(), v->y(), v->z());
std::map<MVertex*,std::pair<double,int> >::iterator itv = LCS.find(v);
if (itv != LCS.end()){
@@ -140,7 +140,7 @@ static void gmsh2MMG(GRegion *gr, MMG_pMesh mmg, MMG_pSol sol,
//if (CTX::instance()->mesh.lcExtendFromBoundary){
double LL = itv->second.first/itv->second.second;
SMetric3 l4(1./(LL*LL));
- SMetric3 MM = intersection_conserve_mostaniso (l4, m);
+ SMetric3 MM = intersection_conserve_mostaniso (l4, m);
m = MM;
//lc = std::min(LL,lc);
// }
@@ -153,16 +153,16 @@ static void gmsh2MMG(GRegion *gr, MMG_pMesh mmg, MMG_pSol sol,
sol->met[count++] = m(2,1);
sol->met[count++] = m(2,2);
// printf("%g %g %g %g %g %g\n",m(0,0),m(0,1),m(0,2),m(1,1),m(1,2),m(2,2));
-
+
// for (int i=0; i<sol->offset; i++) {
// sol->met[isol + i] = lc;
// printf("sol[%d] = %12.5E\n",isol + i,lc);
// }
k++;
}
-
- for (k=1; k<=mmg->ne; k++) {
- MMG_pTetra ptetra = &mmg->tetra[k];
+
+ for (k = 1; k <= mmg->ne; k++) {
+ MMG_pTetra ptetra = &mmg->tetra[k];
ptetra->v[0] = gmsh2mmg_num[gr->tetrahedra[k-1]->getVertex(0)->getNum()];
ptetra->v[1] = gmsh2mmg_num[gr->tetrahedra[k-1]->getVertex(1)->getNum()];
ptetra->v[2] = gmsh2mmg_num[gr->tetrahedra[k-1]->getVertex(2)->getNum()];
@@ -172,27 +172,28 @@ static void gmsh2MMG(GRegion *gr, MMG_pMesh mmg, MMG_pSol sol,
k = 1;
for (std::list<GFace*>::iterator it = f.begin(); it != f.end() ; ++it){
- for (int i=0;i<(*it)->triangles.size();i++){
- MMG_pTria ptriangle = &mmg->tria[k];
+ for (unsigned int i = 0; i < (*it)->triangles.size(); i++){
+ MMG_pTria ptriangle = &mmg->tria[k];
ptriangle->v[0] = gmsh2mmg_num[(*it)->triangles[i]->getVertex(0)->getNum()];
ptriangle->v[1] = gmsh2mmg_num[(*it)->triangles[i]->getVertex(1)->getNum()];
ptriangle->v[2] = gmsh2mmg_num[(*it)->triangles[i]->getVertex(2)->getNum()];
ptriangle->ref = (*it)->tag();
k++;
}
- }
+ }
//mmg->disp = 0;
-
+
}
-static void updateSizes(GRegion *gr, MMG_pMesh mmg, MMG_pSol sol, std::map<int,MVertex*> &mmg2gmsh)
+static void updateSizes(GRegion *gr, MMG_pMesh mmg, MMG_pSol sol,
+ std::map<int,MVertex*> &mmg2gmsh)
{
std::list<GFace*> f = gr->faces();
-
+
std::map<MVertex*,std::pair<double,int> > LCS;
// if (CTX::instance()->mesh.lcExtendFromBoundary){
for (std::list<GFace*>::iterator it = f.begin(); it != f.end() ; ++it){
- for (int i=0;i<(*it)->triangles.size();i++){
+ for (unsigned int i = 0; i < (*it)->triangles.size(); i++){
MTriangle *t = (*it)->triangles[i];
double L = t->maxEdge();
for (int k=0;k<3;k++){
@@ -209,24 +210,24 @@ static void updateSizes(GRegion *gr, MMG_pMesh mmg, MMG_pSol sol, std::map<int,M
}
}
// }
-
+
int count = 1;
for (int k=1 ; k<=mmg->np; k++){
- MMG_pPoint ppt = &mmg->point[k];
+ MMG_pPoint ppt = &mmg->point[k];
if (ppt->tag & M_UNUSED) continue;
-
- SMetric3 m = BGM_MeshMetric(gr, 0,0,ppt->c[0],ppt->c[1],ppt->c[2]);
-
+
+ SMetric3 m = BGM_MeshMetric(gr, 0,0,ppt->c[0],ppt->c[1],ppt->c[2]);
+
std::map<int,MVertex*>::iterator it = mmg2gmsh.find(k);
-
+
if (it != mmg2gmsh.end() && CTX::instance()->mesh.lcExtendFromBoundary){
std::map<MVertex*,std::pair<double,int> >::iterator itv = LCS.find(it->second);
if (itv != LCS.end()){
double LL = itv->second.first/itv->second.second;
SMetric3 l4(1./(LL*LL));
// printf("adding a size %g\n",LL);
- SMetric3 MM = intersection_conserve_mostaniso (l4, m);
+ SMetric3 MM = intersection_conserve_mostaniso (l4, m);
m = MM;
}
}
@@ -235,7 +236,7 @@ static void updateSizes(GRegion *gr, MMG_pMesh mmg, MMG_pSol sol, std::map<int,M
m(1,1) += 1.e-12;
m(2,2) += 1.e-12;
}
-
+
sol->met[count++] = m(0,0);
sol->met[count++] = m(1,0);
sol->met[count++] = m(2,0);
@@ -244,7 +245,7 @@ static void updateSizes(GRegion *gr, MMG_pMesh mmg, MMG_pSol sol, std::map<int,M
sol->met[count++] = m(2,2);
}
free(sol->metold);
- sol->metold = (double*)calloc(sol->npmax+1,sol->offset*sizeof(double));
+ sol->metold = (double*)calloc(sol->npmax+1,sol->offset*sizeof(double));
}
static void freeMMG(MMG_pMesh mmgMesh, MMG_pSol mmgSol)
@@ -255,7 +256,7 @@ static void freeMMG(MMG_pMesh mmgMesh, MMG_pSol mmgSol)
free(mmgMesh->tria);
free(mmgMesh->tetra);
free(mmgMesh);
- //if ( mmgSol->npfixe ){
+ //if ( mmgSol->npfixe ){
free(mmgSol->met);
free(mmgSol->metold);
//}
@@ -264,38 +265,40 @@ static void freeMMG(MMG_pMesh mmgMesh, MMG_pSol mmgSol)
void refineMeshMMG(GRegion *gr)
{
- MMG_pMesh mmg = (MMG_pMesh)calloc(1,sizeof(MMG_Mesh));
- MMG_pSol sol = (MMG_pSol)calloc(1,sizeof(MMG_Sol));
+ MMG_pMesh mmg = (MMG_pMesh)calloc(1,sizeof(MMG_Mesh));
+ MMG_pSol sol = (MMG_pSol)calloc(1,sizeof(MMG_Sol));
std::map<int,MVertex*> mmg2gmsh;
gmsh2MMG (gr, mmg, sol,mmg2gmsh);
-
+
int iterMax = 11;
for (int ITER=0;ITER<iterMax;ITER++){
- int nT = mmg->ne;
+ int nT = mmg->ne;
int verb_mmg = (Msg::GetVerbosity() > 9) ? -1 : 0;
int opt[9] = {1,0,64,0,0,0, verb_mmg , 0,0};
Msg::Debug("-------- GMSH LAUNCHES MMG3D ---------------");
- mmg3d::MMG_mmg3dlib(opt,mmg,sol);
+ mmg3d::MMG_mmg3dlib(opt,mmg,sol);
Msg::Debug("-------- MG3D TERMINATED -------------------");
Msg::Info("MMG3D succeeded (ITER=%d) %d vertices %d tetrahedra",
ITER, mmg->np, mmg->ne);
// Here we should interact with BGM
updateSizes(gr,mmg, sol,mmg2gmsh);
- int nTnow = mmg->ne;
+ int nTnow = mmg->ne;
if (fabs((double)(nTnow - nT)) < 0.05 * nT) break;
- }
+ }
- //char test[] = "test.mesh";
+ //char test[] = "test.mesh";
//MMG_saveMesh(mmg, test);
gr->deleteVertexArrays();
- for (int i=0;i<gr->tetrahedra.size();++i)delete gr->tetrahedra[i];
+ for (unsigned int i = 0; i < gr->tetrahedra.size();++i)
+ delete gr->tetrahedra[i];
gr->tetrahedra.clear();
- for (int i=0;i<gr->mesh_vertices.size();++i)delete gr->mesh_vertices[i];
+ for (unsigned int i = 0; i < gr->mesh_vertices.size(); ++i)
+ delete gr->mesh_vertices[i];
gr->mesh_vertices.clear();
-
+
MMG2gmsh(gr, mmg, mmg2gmsh);
freeMMG(mmg, sol);
}
diff --git a/contrib/bamg/Mesh2d.hpp b/contrib/bamg/Mesh2d.hpp
index 4c54887ceae1b681c29dd92d291761e2b70705a8..9606a379e899ca0f45edace77b693f763cd42e6a 100644
--- a/contrib/bamg/Mesh2d.hpp
+++ b/contrib/bamg/Mesh2d.hpp
@@ -1,11 +1,11 @@
-// la regle de programmation 3
-#include "assertion.hpp"
+// la regle de programmation 3
+#include "assertion.hpp"
#include <cstdlib>
using namespace std;
// definition R
#include "Label.hpp"
#include "R2.hpp"
-
+
class Vertex2 : public R2,public Label {
friend ostream& operator <<(ostream& f, const Vertex2 & v )
@@ -17,7 +17,7 @@ class Vertex2 : public R2,public Label {
// Vertex2(R2 P,int r=0): R2(P),Label(r){}
private: // pas de copie pour ne pas perdre l'adresse
Vertex2(const Vertex2 &);
- void operator=(const Vertex2 &);
+ void operator=(const Vertex2 &);
};
@@ -27,7 +27,7 @@ class Vertex2 : public R2,public Label {
class Triangle2: public Label {
- // variable prive
+ // variable prive
Vertex2 *vertices[3]; // an array of 3 pointer to vertex
public:
R area;
@@ -38,20 +38,20 @@ class Triangle2: public Label {
Vertex2 * & operator()(int i) {
ASSERTION(i>=0 && i <3);
return vertices[i];} // to see traingle as a array of vertex
- void init(Vertex2 * v0,int * iv,int r)
+ void init(Vertex2 * v0,int * iv,int r)
{ vertices[0]=v0+iv[0];
vertices[1]=v0+iv[1];
- vertices[2]=v0+iv[2];
+ vertices[2]=v0+iv[2];
R2 AB(*vertices[0],*vertices[1]);
R2 AC(*vertices[0],*vertices[2]);
area = (AB^AC)*0.5;
lab=r;
- ASSERTION(area>0);
+ ASSERTION(area>0);
}
R2 Edge(int i) const {ASSERTION(i>=0 && i <3);
return R2(*vertices[(i+1)%3],*vertices[(i+2)%3]);}// opposite edge vertex i
R2 H(int i) const { ASSERTION(i>=0 && i <3);
- R2 E=Edge(i);return E.perp()/(2*area);} // heigth
+ R2 E=Edge(i);return E.perp()/(2*area);} // heigth
void Gradlambda(R2 * GradL) const
{
@@ -63,30 +63,30 @@ class Triangle2: public Label {
R lenEdge(int i) const {ASSERTION(i>=0 && i <3);
R2 E=Edge(i);return sqrt((E,E));}
- // Transformation: $\hat{K} \mapsto K$
+ // Transformation: $\hat{K} \mapsto K$
R2 operator()(const R2 & Phat) const {
const R2 &A =*vertices[0];
const R2 &B =*vertices[1];
const R2 &C =*vertices[2];
- return (1-Phat.x- Phat.y)* A + Phat.x *B +Phat.y*C ;}
-
+ return (1-Phat.x- Phat.y)* A + Phat.x *B +Phat.y*C ;}
+
private:
Triangle2(const Triangle2 &); // pas de construction par copie
- void operator=(const Triangle2 &);// pas affectation par copy
+ void operator=(const Triangle2 &);// pas affectation par copy
// Ajoute FH dernier cours ----
public:
R mesure() const {return area;}
- static const int nv=3;
+ static const int nv=3;
};
-// la classe pour le maillage du bord, ici des Segment:
+// la classe pour le maillage du bord, ici des Segment:
class Seg: public Label {
public:
- static const int nv=2; // the nomber of vertices
+ static const int nv=2; // the nomber of vertices
private:
Vertex2 *vertices[nv]; // an array of 2 pointer to vertex
public:
@@ -99,35 +99,35 @@ public:
return *vertices[i];} // to see the segment as a array of vertex
void init(Vertex2 * v0,int * iv,int r) // the true constructor
- {
+ {
vertices[0]=v0+iv[0];
vertices[1]=v0+iv[1];
R2 AB(*vertices[0],*vertices[1]);
l= AB.norme();
lab=r;
- ASSERTION(l>0);
+ ASSERTION(l>0);
}
-
- // Transformation: $[0,1] \mapsto K$
+
+ // Transformation: $[0,1] \mapsto K$
R2 operator()(const R & Phat) const {
const R2 &A =*vertices[0];
const R2 &B =*vertices[1];
- return (1-Phat)* A + Phat *B ;}
+ return (1-Phat)* A + Phat *B ;}
R mesure() const {return l;}
private:
Seg(const Seg &); // pas de construction par copie
- void operator=(const Seg &);// pas affectation par copy
+ void operator=(const Seg &);// pas affectation par copy
};
-class Mesh2
-{
+class Mesh2
+{
public:
- typedef Triangle2 Element;
+ typedef Triangle2 Element;
typedef Seg BorderElement;
typedef R2 Rd;
typedef Vertex2 Vertex;
@@ -141,7 +141,7 @@ public:
Seg & be(int i) const {return borderelements[CheckBE(i)];} // boundary element ..
Mesh2(const char * filename); // read on a file
Mesh2(int nnv,int nnt,int nnbe, Vertex2 *v,Triangle2 *t, Seg *b_e)
- : nv(nnv),nt(nnt),nbe(nnbe),vertices(v),triangles(t),borderelements(b_e),area(0.),peri(0.)
+ : nv(nnv),nt(nnt),nbe(nnbe),area(0.),peri(0.),vertices(v),triangles(t),borderelements(b_e)
{
for(int i=0;i<nt;++i) area += triangles[i].mesure();
for(int i=0;i<nbe;++i) peri += borderelements[i].mesure();
@@ -153,18 +153,18 @@ public:
int operator()(const Seg & v) const {return CheckBE(&v - borderelements);}
int operator()(const Seg * v) const{return CheckBE(v - borderelements);}
int operator()(int it,int j) const {return (*this)(triangles[it][j]);}// Nu vertex j of triangle it
- // to check the bound
- int CheckV(int i) const { ASSERTION(i>=0 && i < nv); return i;}
+ // to check the bound
+ int CheckV(int i) const { ASSERTION(i>=0 && i < nv); return i;}
int CheckT(int i) const { ASSERTION(i>=0 && i < nt); return i;}
int CheckBE(int i) const { ASSERTION(i>=0 && i < nbe); return i;}
- ~Mesh2() {
- delete [] vertices;
+ ~Mesh2() {
+ delete [] vertices;
delete [] triangles;
delete [] borderelements;}
private:
Mesh2(const Mesh2 &); // pas de construction par copie
- void operator=(const Mesh2 &);// pas affectation par copy
+ void operator=(const Mesh2 &);// pas affectation par copy
};
-
+
diff --git a/contrib/mmg3d/build/sources/libmmg3d.h b/contrib/mmg3d/build/sources/libmmg3d.h
index 22c002499ccaac300baea921c11557f883ce538e..4f1afacec9c2ef3675bffb9501d7d2e506fa6782 100644
--- a/contrib/mmg3d/build/sources/libmmg3d.h
+++ b/contrib/mmg3d/build/sources/libmmg3d.h
@@ -3,32 +3,32 @@ Logiciel initial: MMG3D Version 4.0
Co-auteurs : Cecile Dobrzynski et Pascal Frey.
Propriétaires :IPB - UPMC -INRIA.
-Copyright © 2004-2005-2006-2007-2008-2009-2010-2011,
+Copyright © 2004-2005-2006-2007-2008-2009-2010-2011,
diffusé sous les termes et conditions de la licence publique générale de GNU
-Version 3 ou toute version ultérieure.
+Version 3 ou toute version ultérieure.
Ce fichier est une partie de MMG3D.
MMG3D est un logiciel libre ; vous pouvez le redistribuer et/ou le modifier
suivant les termes de la licence publique générale de GNU
Version 3 ou toute version ultérieure.
-MMG3D est distribué dans l'espoir qu'il sera utile, mais SANS
-AUCUNE GARANTIE ; sans même garantie de valeur marchande.
+MMG3D est distribué dans l'espoir qu'il sera utile, mais SANS
+AUCUNE GARANTIE ; sans même garantie de valeur marchande.
Voir la licence publique générale de GNU pour plus de détails.
-MMG3D est diffusé en espérant qu’il sera utile,
-mais SANS AUCUNE GARANTIE, ni explicite ni implicite,
-y compris les garanties de commercialisation ou
-d’adaptation dans un but spécifique.
+MMG3D est diffusé en espérant qu’il sera utile,
+mais SANS AUCUNE GARANTIE, ni explicite ni implicite,
+y compris les garanties de commercialisation ou
+d’adaptation dans un but spécifique.
Reportez-vous à la licence publique générale de GNU pour plus de détails.
-Vous devez avoir reçu une copie de la licence publique générale de GNU
-en même temps que ce document.
+Vous devez avoir reçu une copie de la licence publique générale de GNU
+en même temps que ce document.
Si ce n’est pas le cas, aller voir <http://www.gnu.org/licenses/>.
-/****************************************************************************
+****************************************************************************
Initial software: MMG3D Version 4.0
Co-authors: Cecile Dobrzynski et Pascal Frey.
Owners: IPB - UPMC -INRIA.
-Copyright © 2004-2005-2006-2007-2008-2009-2010-2011,
-spread under the terms and conditions of the license GNU General Public License
+Copyright © 2004-2005-2006-2007-2008-2009-2010-2011,
+spread under the terms and conditions of the license GNU General Public License
as published Version 3, or (at your option) any later version.
This file is part of MMG3D
@@ -41,7 +41,7 @@ but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
-along with MMG3D. If not, see <http://www.gnu.org/licenses/>.
+along with MMG3D. If not, see <http://www.gnu.org/licenses/>.
****************************************************************************/
typedef struct {
double c[3];
@@ -57,7 +57,7 @@ typedef struct {
int mark;
double qual;
int ref,bdryref[4];
- unsigned char flag,edge,tabedg;
+ unsigned char flag,edge,tabedg;
unsigned char bdryinfo[6];
} MMG_Tetra;
typedef MMG_Tetra * MMG_pTetra;
@@ -71,7 +71,7 @@ typedef MMG_Tria * MMG_pTria;
typedef struct {
int np,ver;
double *mv,*cold;
- short *alpha;
+ short *alpha;
} MMG_Displ;
typedef MMG_Displ * MMG_pDispl;
@@ -81,7 +81,7 @@ typedef struct {
short imprim,option,memory,rn,rn2;
int bucksiz;
double delta,dt;
- double min[3],max[3];
+ double min[3],max[3];
} MMG_Info;
@@ -104,7 +104,7 @@ typedef MMG_Mesh * MMG_pMesh;
typedef struct {
int np,npfixe,npmax,ver;
double *met,hmin,hmax;
- char *name;
+ char *name;
double *metold;
unsigned char offset;
} MMG_Sol;