diff --git a/.clang-format b/.clang-format
index 27259a73fb47150addb12dcd8a4a1569009bafaa..8fad88790abb39639111cde25eec3134aa771de5 100644
--- a/.clang-format
+++ b/.clang-format
@@ -40,7 +40,6 @@ BreakBeforeBraces: Custom
BreakBeforeInheritanceComma: true
BreakBeforeTernaryOperators: false
BreakConstructorInitializersBeforeComma: false
-BreakConstructorInitializers: BeforeComma
BreakAfterJavaFieldAnnotations: false
BreakStringLiterals: true
ColumnLimit: 80
@@ -54,11 +53,11 @@ DerivePointerAlignment: false
DisableFormat: false
ExperimentalAutoDetectBinPacking: false
FixNamespaceComments: true
-ForEachMacros:
+ForEachMacros:
- foreach
- Q_FOREACH
- BOOST_FOREACH
-IncludeCategories:
+IncludeCategories:
- Regex: '^"(llvm|llvm-c|clang|clang-c)/'
Priority: 2
- Regex: '^(<|"(gtest|gmock|isl|json)/)'
@@ -105,4 +104,3 @@ Standard: Cpp03
TabWidth: 4
UseTab: Never
...
-
diff --git a/Geo/OCCEdge.cpp b/Geo/OCCEdge.cpp
index 3ecd71268239475e1d8f379edb139fa3d1571e47..a451333f1ea0e1a83a899878323bc194aacfcca6 100644
--- a/Geo/OCCEdge.cpp
+++ b/Geo/OCCEdge.cpp
@@ -278,8 +278,8 @@ int OCCEdge::minimumMeshSegments() const
int np;
// if it is a seam, then return 1
- if (l_faces.size() == 1 && isSeam (l_faces[0]))return 1;
-
+ if(l_faces.size() == 1 && isSeam(l_faces[0])) return 1;
+
if(geomType() == Line)
np = GEdge::minimumMeshSegments();
else
diff --git a/Geo/OCCFace.cpp b/Geo/OCCFace.cpp
index a9d7bb5006eba7b8009fc829c64c729a4cc58273..a10dd2b6ebb972656264185676edc11f54270d89 100644
--- a/Geo/OCCFace.cpp
+++ b/Geo/OCCFace.cpp
@@ -50,10 +50,8 @@ OCCFace::OCCFace(GModel *m, TopoDS_Face _s, int num)
setup();
if(model()->getOCCInternals())
model()->getOCCInternals()->bind(s, num);
- // TEST
- if (tag() == 10){
- writeBREP("s10.brep");
- }
+
+ // if(tag() == 10) writeBREP("s10.brep");
}
OCCFace::~OCCFace()
diff --git a/Geo/OCCRegion.cpp b/Geo/OCCRegion.cpp
index f90de851122341e16ea5876415497a3b613daf99..c2b6a02e8c8e8f57c00cd6323a4817317e0b4e79 100644
--- a/Geo/OCCRegion.cpp
+++ b/Geo/OCCRegion.cpp
@@ -29,10 +29,8 @@ OCCRegion::OCCRegion(GModel *m, TopoDS_Solid _s, int num)
setup();
if(model()->getOCCInternals())
model()->getOCCInternals()->bind(s, num);
- // if (tag() == 1){
- // char name[256];
- // sprintf(name,"v%d.brep",tag());
- // writeBREP(name);
+
+ // if(tag() == 1) writeBREP("v1.brep");
}
OCCRegion::~OCCRegion()
@@ -121,7 +119,8 @@ GEntity::GeomType OCCRegion::geomType() const
return Volume;
}
-void OCCRegion::writeBREP (const char *filename){
+void OCCRegion::writeBREP(const char *filename)
+{
BRep_Builder b;
TopoDS_Compound c;
b.MakeCompound(c);
@@ -129,6 +128,4 @@ void OCCRegion::writeBREP (const char *filename){
BRepTools::Write(c, filename);
}
-
-
#endif
diff --git a/Geo/OCCRegion.h b/Geo/OCCRegion.h
index 050e0c1c3b7a24cf52e8dbe6362f81f328028f74..99e3bd8a4f3e4b26b5c75bada03cc7987604f267 100644
--- a/Geo/OCCRegion.h
+++ b/Geo/OCCRegion.h
@@ -24,8 +24,8 @@ class OCCRegion : public GRegion {
virtual GeomType geomType() const;
ModelType getNativeType() const { return OpenCascadeModel; }
void * getNativePtr() const { return (void*)&s; }
- TopoDS_Solid getTopoDS_Shape() {return s;}
- void writeBREP (const char *filename);
+ TopoDS_Solid getTopoDS_Shape() { return s; }
+ void writeBREP(const char *filename);
};
#endif
diff --git a/Mesh/BDS.cpp b/Mesh/BDS.cpp
index 4d5cc13dde3d13a2bde5a23a8ede900a741336e6..12263e2f19dbb6031280f46fd602c73f97ecd446 100644
--- a/Mesh/BDS.cpp
+++ b/Mesh/BDS.cpp
@@ -17,61 +17,63 @@
#include "meshGFaceDelaunayInsertion.h"
#include "qualityMeasures.h"
-double _COS_N ( BDS_Point *_p1, BDS_Point *_p2, BDS_Point *_p3 , GFace *gf){
-
+double _COS_N(BDS_Point *_p1, BDS_Point *_p2, BDS_Point *_p3, GFace *gf)
+{
double n[3];
normal_triangle(_p1, _p2, _p3, n);
- SVector3 N1 = gf->normal (SPoint2(_p1->u,_p1->v));
- SVector3 N2 = gf->normal (SPoint2(_p2->u,_p2->v));
- SVector3 N3 = gf->normal (SPoint2(_p3->u,_p3->v));
- SVector3 N = N1+N2+N3;
+ SVector3 N1 = gf->normal(SPoint2(_p1->u, _p1->v));
+ SVector3 N2 = gf->normal(SPoint2(_p2->u, _p2->v));
+ SVector3 N3 = gf->normal(SPoint2(_p3->u, _p3->v));
+ SVector3 N = N1 + N2 + N3;
N.normalize();
- return N.x()*n[0]+N.y()*n[1]+N.z()*n[2];
+ return N.x() * n[0] + N.y() * n[1] + N.z() * n[2];
}
-
-double BDS_Face_Validity (GFace *gf, BDS_Face *f){
+double BDS_Face_Validity(GFace *gf, BDS_Face *f)
+{
// if (gf->model()->// return 1;
BDS_Point *pts[4];
f->getNodes(pts);
- if (pts[0]->degenerated+pts[1]->degenerated+pts[2]->degenerated < 2){
+ if(pts[0]->degenerated + pts[1]->degenerated + pts[2]->degenerated < 2) {
double qa1 = qmTriangle::gamma(pts[0], pts[1], pts[2]);
- return qa1*_COS_N (pts[0],pts[1],pts[2],gf);
+ return qa1 * _COS_N(pts[0], pts[1], pts[2], gf);
}
return 1.0;
}
-void outputScalarField(std::vector<BDS_Face*> &t, const char *iii, int param, GFace *gf)
+void outputScalarField(std::vector<BDS_Face *> &t, const char *iii, int param,
+ GFace *gf)
{
- if (gf){
-
- FILE* view_c = Fopen("param_mesh_as_it_is_in_3D.pos","w");
- if(!view_c){
+ if(gf) {
+ FILE *view_c = Fopen("param_mesh_as_it_is_in_3D.pos", "w");
+ if(!view_c) {
Msg::Error("Could not open file param_mesh_as_it_is_in_3D.pos");
return;
}
- fprintf(view_c,"View \"paramC\"{\n");
- std::vector<BDS_Face*>::iterator tit = t.begin();
- std::vector<BDS_Face*>::iterator tite = t.end();
+ fprintf(view_c, "View \"paramC\"{\n");
+ std::vector<BDS_Face *>::iterator tit = t.begin();
+ std::vector<BDS_Face *>::iterator tite = t.end();
while(tit != tite) {
BDS_Point *pts[4];
- if (!(*tit)->deleted){
- (*tit)->getNodes(pts);
- for (int j=0;j<3;j++){
- double U1 = pts[j]->degenerated ? pts[(j+1)%3]->u : pts[j]->u;
- double U2 = pts[(j+1)%3]->degenerated ? pts[j]->u : pts[(j+1)%3]->u ;
- SPoint2 p1(U1,pts[j]->v);
- SPoint2 p2(U2,pts[(j+1)%3]->v);
- SPoint2 prev = p1;
- for (int k=1;k<30;k++){
- double t = (double)k/29;
- SPoint2 p = p1*(1.-t) + p2*t;
- GPoint pa = gf->point(p.x(),p.y());
- GPoint pb = gf->point(prev.x(),prev.y());
- fprintf(view_c,"SL(%g,%g,%g,%g,%g,%g){1,1,1};\n",pa.x(),pa.y(),pa.z(),pb.x(),pb.y(),pb.z());
- prev = p;
- }
- }
+ if(!(*tit)->deleted) {
+ (*tit)->getNodes(pts);
+ for(int j = 0; j < 3; j++) {
+ double U1 = pts[j]->degenerated ? pts[(j + 1) % 3]->u : pts[j]->u;
+ double U2 =
+ pts[(j + 1) % 3]->degenerated ? pts[j]->u : pts[(j + 1) % 3]->u;
+ SPoint2 p1(U1, pts[j]->v);
+ SPoint2 p2(U2, pts[(j + 1) % 3]->v);
+ SPoint2 prev = p1;
+ for(int k = 1; k < 30; k++) {
+ double t = (double)k / 29;
+ SPoint2 p = p1 * (1. - t) + p2 * t;
+ GPoint pa = gf->point(p.x(), p.y());
+ GPoint pb = gf->point(prev.x(), prev.y());
+ fprintf(view_c, "SL(%g,%g,%g,%g,%g,%g){1,1,1};\n", pa.x(), pa.y(),
+ pa.z(), pb.x(), pb.y(), pb.z());
+ prev = p;
+ }
+ }
}
++tit;
}
@@ -91,53 +93,52 @@ void outputScalarField(std::vector<BDS_Face*> &t, const char *iii, int param, GF
BDS_Point *pts[4];
if(!(*tit)->deleted) {
(*tit)->getNodes(pts);
- //double v = BDS_Face_Validity (gf, *tit);
- if(!param && gf){
- if (pts[0]->degenerated+pts[1]->degenerated+pts[2]->degenerated < 2)
- {
- fprintf(f, "ST(%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E){%g,%g,%g};\n",
- pts[0]->X, pts[0]->Y, pts[0]->Z,
- pts[1]->X, pts[1]->Y, pts[1]->Z,
- pts[2]->X, pts[2]->Y, pts[2]->Z,(double)pts[0]->iD, (double)pts[1]->iD, (double)pts[2]->iD);
- }
+ // double v = BDS_Face_Validity (gf, *tit);
+ if(!param && gf) {
+ if(pts[0]->degenerated + pts[1]->degenerated + pts[2]->degenerated <
+ 2) {
+ fprintf(f,
+ "ST(%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%"
+ "22.15E,%22.15E){%g,%g,%g};\n",
+ pts[0]->X, pts[0]->Y, pts[0]->Z, pts[1]->X, pts[1]->Y,
+ pts[1]->Z, pts[2]->X, pts[2]->Y, pts[2]->Z,
+ (double)pts[0]->iD, (double)pts[1]->iD, (double)pts[2]->iD);
+ }
}
- if(param && gf){
- if (pts[0]->degenerated+pts[1]->degenerated+pts[2]->degenerated > 1)
- {
- }
- else if (pts[0]->degenerated)
- fprintf(f, "SQ(%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g,%g};\n",
- pts[1]->u, pts[1]->v, 0.0,
- pts[2]->u, pts[2]->v, 0.0,
- pts[2]->u, pts[0]->v, 0.0,
- pts[1]->u, pts[0]->v, 0.0,
- // v,v,v,v);
- (double)pts[1]->iD, (double)pts[2]->iD, (double)pts[0]->iD, (double)pts[0]->iD);
- else if (pts[1]->degenerated)
- fprintf(f, "SQ(%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g,%g};\n",
- pts[2]->u, pts[2]->v, 0.0,
- pts[0]->u, pts[0]->v, 0.0,
- pts[0]->u, pts[1]->v, 0.0,
- pts[2]->u, pts[1]->v, 0.0,
- //v,v,v,v);
- (double)pts[2]->iD, (double)pts[0]->iD, (double)pts[1]->iD, (double)pts[1]->iD);
- else if (pts[2]->degenerated)
- fprintf(f, "SQ(%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g,%g};\n",
- pts[0]->u, pts[0]->v, 0.0,
- pts[1]->u, pts[1]->v, 0.0,
- pts[1]->u, pts[2]->v, 0.0,
- pts[0]->u, pts[2]->v, 0.0,
- //v,v,v,v);
- (double)pts[0]->iD, (double)pts[1]->iD, (double)pts[2]->iD, (double)pts[2]->iD);
- else{
- fprintf(f, "ST(%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E){%g,%g,%g};\n",
- pts[0]->u, pts[0]->v, 0.0,
- pts[1]->u, pts[1]->v, 0.0,
- pts[2]->u, pts[2]->v, 0.0,
- // v,v,v);
- (double)pts[0]->iD, (double)pts[1]->iD, (double)pts[2]->iD);
-
- }
+ if(param && gf) {
+ if(pts[0]->degenerated + pts[1]->degenerated + pts[2]->degenerated >
+ 1) {
+ }
+ else if(pts[0]->degenerated)
+ fprintf(f, "SQ(%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g,%g};\n",
+ pts[1]->u, pts[1]->v, 0.0, pts[2]->u, pts[2]->v, 0.0,
+ pts[2]->u, pts[0]->v, 0.0, pts[1]->u, pts[0]->v, 0.0,
+ // v,v,v,v);
+ (double)pts[1]->iD, (double)pts[2]->iD, (double)pts[0]->iD,
+ (double)pts[0]->iD);
+ else if(pts[1]->degenerated)
+ fprintf(f, "SQ(%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g,%g};\n",
+ pts[2]->u, pts[2]->v, 0.0, pts[0]->u, pts[0]->v, 0.0,
+ pts[0]->u, pts[1]->v, 0.0, pts[2]->u, pts[1]->v, 0.0,
+ // v,v,v,v);
+ (double)pts[2]->iD, (double)pts[0]->iD, (double)pts[1]->iD,
+ (double)pts[1]->iD);
+ else if(pts[2]->degenerated)
+ fprintf(f, "SQ(%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){%g,%g,%g,%g};\n",
+ pts[0]->u, pts[0]->v, 0.0, pts[1]->u, pts[1]->v, 0.0,
+ pts[1]->u, pts[2]->v, 0.0, pts[0]->u, pts[2]->v, 0.0,
+ // v,v,v,v);
+ (double)pts[0]->iD, (double)pts[1]->iD, (double)pts[2]->iD,
+ (double)pts[2]->iD);
+ else {
+ fprintf(f,
+ "ST(%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%22.15E,%"
+ "22.15E,%22.15E){%g,%g,%g};\n",
+ pts[0]->u, pts[0]->v, 0.0, pts[1]->u, pts[1]->v, 0.0,
+ pts[2]->u, pts[2]->v, 0.0,
+ // v,v,v);
+ (double)pts[0]->iD, (double)pts[1]->iD, (double)pts[2]->iD);
+ }
}
}
++tit;
@@ -146,11 +147,8 @@ void outputScalarField(std::vector<BDS_Face*> &t, const char *iii, int param, GF
fclose(f);
}
-
BDS_Vector::BDS_Vector(const BDS_Point &p2, const BDS_Point &p1)
- : x(p2.X - p1.X)
- , y(p2.Y - p1.Y)
- , z(p2.Z - p1.Z)
+ : x(p2.X - p1.X), y(p2.Y - p1.Y), z(p2.Z - p1.Z)
{
}
@@ -178,7 +176,6 @@ void normal_triangle(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3, double c[3])
norme(c);
}
-
double surface_triangle_param(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3)
{
// FIXME
@@ -186,19 +183,19 @@ double surface_triangle_param(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3)
// SEEMS TO BE THE CASE WITH OCC
double c;
- if (p1->degenerated+p2->degenerated+p3->degenerated > 1)
- c = 0;//vector_triangle_parametric(p1, p2, p3, c);
- else if (p1->degenerated){
+ if(p1->degenerated + p2->degenerated + p3->degenerated > 1)
+ c = 0; // vector_triangle_parametric(p1, p2, p3, c);
+ else if(p1->degenerated) {
double du = fabs(p3->u - p2->u);
- c = 2*fabs (0.5*(p3->v+p2->v) - p1->v) * du;
+ c = 2 * fabs(0.5 * (p3->v + p2->v) - p1->v) * du;
}
- else if (p2->degenerated){
+ else if(p2->degenerated) {
double du = fabs(p3->u - p1->u);
- c = 2*fabs (0.5*(p3->v+p1->v) - p2->v) * du;
+ c = 2 * fabs(0.5 * (p3->v + p1->v) - p2->v) * du;
}
- else if (p3->degenerated){
+ else if(p3->degenerated) {
double du = fabs(p2->u - p1->u);
- c = 2*fabs (0.5*(p2->v+p1->v) - p3->v) * du;
+ c = 2 * fabs(0.5 * (p2->v + p1->v) - p3->v) * du;
}
else
c = vector_triangle_parametric(p1, p2, p3);
@@ -321,9 +318,9 @@ BDS_Edge *BDS_Mesh::recover_edge_fast(BDS_Point *p1, BDS_Point *p2)
BDS_Face *f = e->otherFace(t);
if(!f->e4) {
BDS_Point *p2b = f->oppositeVertex(e);
- if (p2 == p2b){
- if (swap_edge(e, BDS_SwapEdgeTestRecover(), true)){
- return find_edge (p1,p2->iD);
+ if(p2 == p2b) {
+ if(swap_edge(e, BDS_SwapEdgeTestRecover(), true)) {
+ return find_edge(p1, p2->iD);
}
}
}
@@ -423,8 +420,7 @@ BDS_Edge *BDS_Mesh::recover_edge(int num1, int num2, bool &_fatal,
std::vector<int>::size_type ichoice = 0;
bool success = false;
- while (!success && ichoice < intersected.size())
- {
+ while(!success && ichoice < intersected.size()) {
success = swap_edge(intersected[ichoice++], BDS_SwapEdgeTestRecover());
}
@@ -465,15 +461,21 @@ BDS_Face *BDS_Mesh::find_triangle(BDS_Edge *e1, BDS_Edge *e2, BDS_Edge *e3)
{
for(int i = 0; i < e1->numfaces(); i++) {
BDS_Face *t = e1->faces(i);
- if(is_equivalent(e1, e2, e3, t->e1, t->e2, t->e3)) { return t; }
+ if(is_equivalent(e1, e2, e3, t->e1, t->e2, t->e3)) {
+ return t;
+ }
}
for(int i = 0; i < e2->numfaces(); i++) {
BDS_Face *t = e2->faces(i);
- if(is_equivalent(e1, e2, e3, t->e1, t->e2, t->e3)) { return t; }
+ if(is_equivalent(e1, e2, e3, t->e1, t->e2, t->e3)) {
+ return t;
+ }
}
for(int i = 0; i < e3->numfaces(); i++) {
BDS_Face *t = e3->faces(i);
- if(is_equivalent(e1, e2, e3, t->e1, t->e2, t->e3)) { return t; }
+ if(is_equivalent(e1, e2, e3, t->e1, t->e2, t->e3)) {
+ return t;
+ }
}
return 0;
}
@@ -540,9 +542,8 @@ void BDS_Mesh::add_geom(int p1, int p2)
geom.insert(new BDS_GeomEntity(p1, p2));
}
-void BDS_Edge::computeNeighborhood(BDS_Point *pts1[4],
- BDS_Point *pts2[4],
- BDS_Point * oface[2]) const
+void BDS_Edge::computeNeighborhood(BDS_Point *pts1[4], BDS_Point *pts2[4],
+ BDS_Point *oface[2]) const
{
oface[0] = oface[1] = 0;
if(faces(0)) {
@@ -567,7 +568,6 @@ void BDS_Edge::computeNeighborhood(BDS_Point *pts1[4],
void BDS_Edge::oppositeof(BDS_Point *oface[2]) const
{
-
// FIXME !!!!!!!!!!!!
oface[0] = oface[1] = 0;
@@ -703,14 +703,15 @@ bool BDS_Mesh::split_edge(BDS_Edge *e, BDS_Point *mid)
// printf("coucou %d %d %d %d\n",p1->iD,p2->iD,op[0]->iD,op[1]->iD);
double ori0 = fabs(surface_triangle_param(p2, p1, op[0])) +
- fabs(surface_triangle_param(p2, p1, op[1]));
+ fabs(surface_triangle_param(p2, p1, op[1]));
double ori1 = fabs(surface_triangle_param(mid, p1, op[1]));
double ori2 = fabs(surface_triangle_param(mid, op[1], p2));
double ori3 = fabs(surface_triangle_param(mid, p2, op[0]));
double ori4 = fabs(surface_triangle_param(mid, op[0], p1));
double eps = 1.e-2;
- if (ori1< eps*ori0 || ori2 < eps*ori0 || ori3 < eps*ori0 || ori4 < eps*ori0){
+ if(ori1 < eps * ori0 || ori2 < eps * ori0 || ori3 < eps * ori0 ||
+ ori4 < eps * ori0) {
// printf("%g %g %g %g %g\n",ori0,ori1,ori2,ori3,ori4);
// return false;
}
@@ -796,7 +797,6 @@ bool BDS_Mesh::split_edge(BDS_Edge *e, BDS_Point *mid)
bool BDS_SwapEdgeTestRecover::operator()(BDS_Point *_p1, BDS_Point *_p2,
BDS_Point *_q1, BDS_Point *_q2) const
{
-
double p1[2] = {_p1->u, _p1->v};
double p2[2] = {_p2->u, _p2->v};
double op1[2] = {_q1->u, _q1->v};
@@ -818,9 +818,6 @@ bool BDS_SwapEdgeTestRecover::operator()(BDS_Point *_p1, BDS_Point *_p2,
return true;
}
-
-
-
// This function does actually the swap without taking into account
// the feasability of the operation. Those conditions have to be
// taken into account before doing the edge swap
@@ -828,17 +825,15 @@ bool BDS_SwapEdgeTestRecover::operator()(BDS_Point *_p1, BDS_Point *_p2,
bool BDS_SwapEdgeTestQuality::operator()(BDS_Point *_p1, BDS_Point *_p2,
BDS_Point *_q1, BDS_Point *_q2) const
{
-
if(!testSmallTriangles) return true;
double s1 = fabs(surface_triangle_param(_p1, _p2, _q1));
double s2 = fabs(surface_triangle_param(_p1, _p2, _q2));
double s3 = fabs(surface_triangle_param(_p1, _q1, _q2));
double s4 = fabs(surface_triangle_param(_p2, _q1, _q2));
- if(fabs(s1 + s2 - s3 - s4) > 1.e-12*(s3 + s4)) {
+ if(fabs(s1 + s2 - s3 - s4) > 1.e-12 * (s3 + s4)) {
return false;
}
- if(s3 < .02 * (s1 + s2) || s4 < .02 * (s1 + s2))
- return false;
+ if(s3 < .02 * (s1 + s2) || s4 < .02 * (s1 + s2)) return false;
/*
if(!testSmallTriangles) {
@@ -903,7 +898,6 @@ bool BDS_SwapEdgeTestQuality::operator()(BDS_Point *_p1, BDS_Point *_p2,
}
*/
-
bool BDS_SwapEdgeTestQuality::operator()(BDS_Point *_p1, BDS_Point *_p2,
BDS_Point *_p3, BDS_Point *_q1,
BDS_Point *_q2, BDS_Point *_q3,
@@ -929,57 +923,53 @@ bool BDS_SwapEdgeTestQuality::operator()(BDS_Point *_p1, BDS_Point *_p2,
// if(cosnq < 0.3 && cosonq > 0.5 && minb > 0.1) return true;
return minb > mina;
-
}
-
-bool BDS_SwapEdgeTestNormals::operator () (BDS_Point *_p1, BDS_Point *_p2,
- BDS_Point *_q1, BDS_Point *_q2) const
+bool BDS_SwapEdgeTestNormals::operator()(BDS_Point *_p1, BDS_Point *_p2,
+ BDS_Point *_q1, BDS_Point *_q2) const
{
double s1 = fabs(surface_triangle_param(_p1, _p2, _q1));
double s2 = fabs(surface_triangle_param(_p1, _p2, _q2));
double s3 = fabs(surface_triangle_param(_p1, _q1, _q2));
double s4 = fabs(surface_triangle_param(_p2, _q1, _q2));
- if(fabs(s1 + s2 - s3 - s4) > 1.e-12*(s3 + s4)) {
+ if(fabs(s1 + s2 - s3 - s4) > 1.e-12 * (s3 + s4)) {
return false;
}
- if(s3 < .02 * (s1 + s2) || s4 < .02 * (s1 + s2))
- return false;
+ if(s3 < .02 * (s1 + s2) || s4 < .02 * (s1 + s2)) return false;
return true;
}
-bool BDS_SwapEdgeTestNormals::operator () (BDS_Point *_p1, BDS_Point *_p2, BDS_Point *_p3,
- BDS_Point *_q1, BDS_Point *_q2, BDS_Point *_q3,
- BDS_Point *_op1, BDS_Point *_op2, BDS_Point *_op3,
- BDS_Point *_oq1, BDS_Point *_oq2, BDS_Point *_oq3) const
+bool BDS_SwapEdgeTestNormals::operator()(BDS_Point *_p1, BDS_Point *_p2,
+ BDS_Point *_p3, BDS_Point *_q1,
+ BDS_Point *_q2, BDS_Point *_q3,
+ BDS_Point *_op1, BDS_Point *_op2,
+ BDS_Point *_op3, BDS_Point *_oq1,
+ BDS_Point *_oq2, BDS_Point *_oq3) const
{
-
double qa1 = qmTriangle::gamma(_p1, _p2, _p3);
double qa2 = qmTriangle::gamma(_q1, _q2, _q3);
double qb1 = qmTriangle::gamma(_op1, _op2, _op3);
double qb2 = qmTriangle::gamma(_oq1, _oq2, _oq3);
- double mina = std::min(qa1,qa2);
- double minb = std::min(qb1,qb2);
-
- if(minb > 5*mina) return true;
+ double mina = std::min(qa1, qa2);
+ double minb = std::min(qb1, qb2);
- double OLD = std::min(_ori*qa1*_COS_N ( _p1, _p2, _p3 , gf),
- _ori*qa2*_COS_N ( _q1, _q2, _q3 , gf));
+ if(minb > 5 * mina) return true;
- double NEW = std::min(_ori*qb1*_COS_N ( _op1, _op2, _op3 , gf),
- _ori*qb2*_COS_N ( _oq1, _oq2, _oq3 , gf));
+ double OLD = std::min(_ori * qa1 * _COS_N(_p1, _p2, _p3, gf),
+ _ori * qa2 * _COS_N(_q1, _q2, _q3, gf));
+ double NEW = std::min(_ori * qb1 * _COS_N(_op1, _op2, _op3, gf),
+ _ori * qb2 * _COS_N(_oq1, _oq2, _oq3, gf));
// printf("%g %g\n",OLD, NEW);
- if (OLD < 0.2 && OLD < NEW)
- return true;
+ if(OLD < 0.2 && OLD < NEW) return true;
return false;
-
}
-bool BDS_Mesh::swap_edge(BDS_Edge *e, const BDS_SwapEdgeTest &theTest, bool force)
+bool BDS_Mesh::swap_edge(BDS_Edge *e, const BDS_SwapEdgeTest &theTest,
+ bool force)
{
/*
p1
@@ -1000,8 +990,7 @@ bool BDS_Mesh::swap_edge(BDS_Edge *e, const BDS_SwapEdgeTest &theTest, bool forc
BDS_Point *p1 = e->p1;
BDS_Point *p2 = e->p2;
- if(e->deleted)
- return false;
+ if(e->deleted) return false;
int nbFaces = e->numfaces();
if(nbFaces != 2) return false;
@@ -1010,36 +999,33 @@ bool BDS_Mesh::swap_edge(BDS_Edge *e, const BDS_SwapEdgeTest &theTest, bool forc
const int CHECK1 = -1, CHECK2 = -1;
-
BDS_Point *op[2];
BDS_Point *pts1[4];
BDS_Point *pts2[4];
e->computeNeighborhood(pts1, pts2, op);
// if (op[0] == op[1])return false;
- if (p1->iD == CHECK1 && p2->iD == CHECK2){
- printf("-------------- e %d %d --> %d %d\n",p1->iD,p2->iD, op[0]->iD, op[1]->iD);
- printf("-------------- %d %d %d\n",pts1[0]->iD,pts1[1]->iD,pts1[2]->iD);
- printf("-------------- %d %d %d\n",pts2[0]->iD,pts2[1]->iD,pts2[2]->iD);
+ if(p1->iD == CHECK1 && p2->iD == CHECK2) {
+ printf("-------------- e %d %d --> %d %d\n", p1->iD, p2->iD, op[0]->iD,
+ op[1]->iD);
+ printf("-------------- %d %d %d\n", pts1[0]->iD, pts1[1]->iD, pts1[2]->iD);
+ printf("-------------- %d %d %d\n", pts2[0]->iD, pts2[1]->iD, pts2[2]->iD);
}
-
- if (! force &&
- !p1->config_modified &&
- !p2->config_modified &&
- !op[0]->config_modified &&
- !op[1]->config_modified )return false;
+ if(!force && !p1->config_modified && !p2->config_modified &&
+ !op[0]->config_modified && !op[1]->config_modified)
+ return false;
// AVOID CREATING POINTS WITH 2 NEIGHBORING TRIANGLES
// std::vector<BDS_Face*> f1 = p1->getTriangles();
// std::vector<BDS_Face*> f2 = p2->getTriangles();
- // if (p1->g && p1->g->classif_degree == 2 && p1->edges.size() <= 4)return false;
- // if (p2->g && p2->g->classif_degree == 2 && p2->edges.size() <= 4)return false;
- // if (p1->g && p1->g->classif_degree == 1 && p1->edges.size() <= 3)return false;
- // if (p2->g && p2->g->classif_degree == 1 && p2->edges.size() <= 3)return false;
+ // if (p1->g && p1->g->classif_degree == 2 && p1->edges.size() <= 4)return
+ // false; if (p2->g && p2->g->classif_degree == 2 && p2->edges.size() <=
+ // 4)return false; if (p1->g && p1->g->classif_degree == 1 &&
+ // p1->edges.size() <= 3)return false; if (p2->g && p2->g->classif_degree ==
+ // 1 && p2->edges.size() <= 3)return false;
// || p2->edges.size() == 4)return false;
- if (p1->iD == CHECK1 && p2->iD == CHECK2)
- printf("topology OK \n");
+ if(p1->iD == CHECK1 && p2->iD == CHECK2) printf("topology OK \n");
BDS_GeomEntity *g1 = 0, *g2 = 0, *ge = e->g;
@@ -1064,13 +1050,9 @@ bool BDS_Mesh::swap_edge(BDS_Edge *e, const BDS_SwapEdgeTest &theTest, bool forc
return false;
}
- if(!theTest(p1, p2, op[0], op[1]))
- return false;
-
- if (p1->iD == CHECK1 && p2->iD == CHECK2)
- printf("TEST2 OK\n");
-
+ if(!theTest(p1, p2, op[0], op[1])) return false;
+ if(p1->iD == CHECK1 && p2->iD == CHECK2) printf("TEST2 OK\n");
BDS_Edge *p1_op1 = find_edge(p1, op[0], e->faces(0));
BDS_Edge *op1_p2 = find_edge(op[0], p2, e->faces(0));
@@ -1217,7 +1199,6 @@ static bool test_move_point_parametric_triangle(BDS_Point *p, double u,
return ori_init * ori_final > 0;
}
-
/*
VERTICES C AND D LOSE ONE EDGE !!
AS IF TWO EDGE SWAPS WERE DONE
@@ -1234,36 +1215,37 @@ static bool test_move_point_parametric_triangle(BDS_Point *p, double u,
x-------x
*/
-
-
bool BDS_Mesh::collapse_edge_parametric(BDS_Edge *e, BDS_Point *p, bool force)
{
-
- if(!force && e->numfaces() != 2)
- return false;
- if(!force && p->g && p->g->classif_degree == 0)
- return false;
- // not really ok but 'til now this is the best choice not to do collapses on model edges
- if(!force && p->g && p->g->classif_degree == 1)
- return false;
- if(!force && e->g && p->g) {
- if(e->g->classif_degree == 2 && p->g != e->g)
- return false;
+ if(!force && e->numfaces() != 2) return false;
+ if(!force && p->g && p->g->classif_degree == 0) return false;
+ // not really ok but 'til now this is the best choice not to do collapses on
+ // model edges
+ if(!force && p->g && p->g->classif_degree == 1) return false;
+ if(!force && e->g && p->g) {
+ if(e->g->classif_degree == 2 && p->g != e->g) return false;
}
- if (e->numfaces() == 2){
+ if(e->numfaces() == 2) {
BDS_Point *oface[2];
e->oppositeof(oface);
- if (!force && oface[0]->g && oface[0]->g->classif_degree == 2 && oface[0]->edges.size() <= 4)return false;
- if (!force && oface[1]->g && oface[1]->g->classif_degree == 2 && oface[1]->edges.size() <= 4)return false;
- if (!force && oface[0]->g && oface[0]->g->classif_degree == 1 && oface[0]->edges.size() <= 3)return false;
- if (!force && oface[1]->g && oface[1]->g->classif_degree == 1 && oface[1]->edges.size() <= 3)return false;
+ if(!force && oface[0]->g && oface[0]->g->classif_degree == 2 &&
+ oface[0]->edges.size() <= 4)
+ return false;
+ if(!force && oface[1]->g && oface[1]->g->classif_degree == 2 &&
+ oface[1]->edges.size() <= 4)
+ return false;
+ if(!force && oface[0]->g && oface[0]->g->classif_degree == 1 &&
+ oface[0]->edges.size() <= 3)
+ return false;
+ if(!force && oface[1]->g && oface[1]->g->classif_degree == 1 &&
+ oface[1]->edges.size() <= 3)
+ return false;
}
- std::vector<BDS_Face*> t = p->getTriangles();
+ std::vector<BDS_Face *> t = p->getTriangles();
BDS_Point *o = e->othervertex(p);
-
BDS_Point *pt[3][1024];
BDS_GeomEntity *gs[1024];
int ept[2][1024];
@@ -1285,20 +1267,22 @@ bool BDS_Mesh::collapse_edge_parametric(BDS_Edge *e, BDS_Point *p, bool force)
pt[0][nt] = (pts[0] == p) ? o : pts[0];
pt[1][nt] = (pts[1] == p) ? o : pts[1];
pt[2][nt] = (pts[2] == p) ? o : pts[2];
- double snew = std::abs(surface_triangle_param(pt[0][nt], pt[1][nt], pt[2][nt]));
- if (!force && snew < .02 *sold) {
- // printf("argh\n");
- return false;
- }
- area_new += snew;
+ double snew =
+ std::abs(surface_triangle_param(pt[0][nt], pt[1][nt], pt[2][nt]));
+ if(!force && snew < .02 * sold) {
+ // printf("argh\n");
+ return false;
+ }
+ area_new += snew;
++nt;
}
++it;
}
}
- if (!force && nt == 2) return false;
+ if(!force && nt == 2) return false;
- if (!force && fabs(area_old - area_new) > 1.e-12 * (area_old + area_new))return false;
+ if(!force && fabs(area_old - area_new) > 1.e-12 * (area_old + area_new))
+ return false;
{
std::vector<BDS_Face *>::iterator it = t.begin();
@@ -1508,67 +1492,64 @@ bool BDS_Mesh::smooth_point_centroid(BDS_Point *p, GFace *gf, bool test_quality)
return true;
}
-
-
-
-
#else
// Tutte's simple smoothing
// other implementations are coming
bool BDS_Mesh::smooth_point_centroid(BDS_Point *p, GFace *gf, bool hard)
{
- //return false;
+ // return false;
if(p->degenerated) return false;
if(!p->config_modified) return false;
- int CHECK =-1;
+ int CHECK = -1;
if(p->g && p->g->classif_degree <= 1) return false;
if(p->g && p->g->classif_tag < 0) {
p->config_modified = true;
return true;
}
- std::vector<BDS_Edge*>::const_iterator eit = p->edges.begin();
+ std::vector<BDS_Edge *>::const_iterator eit = p->edges.begin();
std::vector<BDS_Edge *>::const_iterator itede = p->edges.end();
while(eit != itede) {
if((*eit)->numfaces() == 1) return false;
eit++;
}
- if (p->iD == CHECK)printf("point %d connected to ",CHECK);
+ if(p->iD == CHECK) printf("point %d connected to ", CHECK);
- double XX=0,YY=0,ZZ=0;
+ double XX = 0, YY = 0, ZZ = 0;
double U = 0;
double V = 0;
double LC = 0;
double oldU = p->u;
double oldV = p->v;
- std::vector<BDS_Face*> ts = p->getTriangles();
+ std::vector<BDS_Face *> ts = p->getTriangles();
double sTot = p->edges.size();
double fact = 0.0;
double ENERGY = 0.0;
eit = p->edges.begin();
- if (eit == itede) {
- Msg::Debug("Hidden bug ... I should have deleted a point but I still do not know why it segfault when I do it :-) ");
+ if(eit == itede) {
+ Msg::Debug("Hidden bug ... I should have deleted a point but I still do "
+ "not know why it segfault when I do it :-) ");
return false;
}
while(eit != itede) {
- BDS_Edge *e = *eit;
+ BDS_Edge *e = *eit;
BDS_Point *n = e->othervertex(p);
double NU = n->degenerated ? p->u : n->u;
- if (p->iD == CHECK)printf("%d ",n->iD);
- double du = sqrt((NU-oldU)*(NU-oldU)+(n->v-oldV)*(n->v-oldV));
- if (du == 0){
+ if(p->iD == CHECK) printf("%d ", n->iD);
+ double du = sqrt((NU - oldU) * (NU - oldU) + (n->v - oldV) * (n->v - oldV));
+ if(du == 0) {
return false;
}
double length = e->length();
- ENERGY += length*length;
- double factloc = du/ length;
- U += NU * factloc ;
- V += n->v * factloc;
+ ENERGY += length * length;
+ double factloc = du / length;
+ U += NU * factloc;
+ V += n->v * factloc;
fact += factloc;
XX += n->X;
YY += n->Y;
@@ -1576,16 +1557,15 @@ bool BDS_Mesh::smooth_point_centroid(BDS_Point *p, GFace *gf, bool hard)
LC += n->lc();
++eit;
}
- if (p->iD == CHECK)printf("\n");
+ if(p->iD == CHECK) printf("\n");
// printf("%g\n",fact);
// sTot *= fact;
U /= (fact);
V /= (fact);
LC /= (sTot);
- XX/= (sTot);
- YY/= (sTot);
- ZZ/= (sTot);
-
+ XX /= (sTot);
+ YY /= (sTot);
+ ZZ /= (sTot);
GPoint gp;
@@ -1599,49 +1579,50 @@ bool BDS_Mesh::smooth_point_centroid(BDS_Point *p, GFace *gf, bool hard)
}
}
else*/
- gp = gf->point(U , V );
+ gp = gf->point(U, V);
- if (!gp.succeeded()){
+ if(!gp.succeeded()) {
return false;
}
eit = p->edges.begin();
double ENERGY_NEW = 0;
while(eit != itede) {
- BDS_Edge *e = *eit;
+ BDS_Edge *e = *eit;
BDS_Point *n = e->othervertex(p);
- double l2 = (gp.x()-n->X)*(gp.x()-n->X)+
- (gp.y()-n->Y)*(gp.y()-n->Y)+
- (gp.z()-n->Z)*(gp.z()-n->Z);
+ double l2 = (gp.x() - n->X) * (gp.x() - n->X) +
+ (gp.y() - n->Y) * (gp.y() - n->Y) +
+ (gp.z() - n->Z) * (gp.z() - n->Z);
ENERGY_NEW += l2;
++eit;
}
// simple strategy has failed to reduce energy
- if (p->iD == CHECK)
- printf("SIMPLE CENTROID SCHEME %g %g\n",ENERGY_NEW,ENERGY);
+ if(p->iD == CHECK)
+ printf("SIMPLE CENTROID SCHEME %g %g\n", ENERGY_NEW, ENERGY);
- if (ENERGY_NEW > ENERGY/* || hard*/){
- double uv[2] = {U,V};
- gp = gf->closestPoint(SPoint3(XX, YY, ZZ),uv);
+ if(ENERGY_NEW > ENERGY /* || hard*/) {
+ double uv[2] = {U, V};
+ gp = gf->closestPoint(SPoint3(XX, YY, ZZ), uv);
U = gp.u();
V = gp.v();
// return false;
eit = p->edges.begin();
ENERGY_NEW = 0;
while(eit != itede) {
- BDS_Edge *e = *eit;
+ BDS_Edge *e = *eit;
BDS_Point *n = e->othervertex(p);
- double l2 = (gp.x()-n->X)*(gp.x()-n->X)+
- (gp.y()-n->Y)*(gp.y()-n->Y)+
- (gp.z()-n->Z)*(gp.z()-n->Z);
+ double l2 = (gp.x() - n->X) * (gp.x() - n->X) +
+ (gp.y() - n->Y) * (gp.y() - n->Y) +
+ (gp.z() - n->Z) * (gp.z() - n->Z);
ENERGY_NEW += l2;
++eit;
}
- if (p->iD == CHECK)
- printf("PROJECTION : %g %g\n",ENERGY_NEW,ENERGY);
- if (ENERGY_NEW > ENERGY){
- Msg::Debug("Impossible to move vertex %d using simple strategies... leaving it there",p->iD);
+ if(p->iD == CHECK) printf("PROJECTION : %g %g\n", ENERGY_NEW, ENERGY);
+ if(ENERGY_NEW > ENERGY) {
+ Msg::Debug("Impossible to move vertex %d using simple strategies... "
+ "leaving it there",
+ p->iD);
return false;
}
}
@@ -1653,8 +1634,8 @@ bool BDS_Mesh::smooth_point_centroid(BDS_Point *p, GFace *gf, bool hard)
double s1 = 0, s2 = 0;
double newWorst = 1.0, oldWorst = 1.0;
- //double OLD = 1, NEW = 1;
- std::vector<BDS_Face*>::const_iterator it = ts.begin();
+ // double OLD = 1, NEW = 1;
+ std::vector<BDS_Face *>::const_iterator it = ts.begin();
while(it != ts.end()) {
BDS_Face *t = *it;
BDS_Point *n[4];
@@ -1668,7 +1649,7 @@ bool BDS_Mesh::smooth_point_centroid(BDS_Point *p, GFace *gf, bool hard)
p->v = oldV;
double sold = std::abs(surface_triangle_param(n[0], n[1], n[2]));
s2 += sold;
- if (snew < 0.02 * sold)return false;
+ if(snew < 0.02 * sold) return false;
p->X = gp.x();
p->Y = gp.y();
p->Z = gp.z();
@@ -1682,13 +1663,13 @@ bool BDS_Mesh::smooth_point_centroid(BDS_Point *p, GFace *gf, bool hard)
}
// if (p->edges.size() == 3)printf("3 -> %22.15E\n",fabs(s2-s1));
- if(fabs(s2-s1) > 1.e-12 * (s2 + s1)) {
+ if(fabs(s2 - s1) > 1.e-12 * (s2 + s1)) {
// if (p->iD == CHECK)
// printf("PARAMETRIC TRIANGLES OVERLAP %22.15E!!\n",fabs(s2-s1));
// else
return false;
}
- if (newWorst < oldWorst) return false;
+ if(newWorst < oldWorst) return false;
// if (p->edges.size() == 3)printf("OK \n");
// if (OLD < 0 && NEW > OLD){
@@ -1696,7 +1677,6 @@ bool BDS_Mesh::smooth_point_centroid(BDS_Point *p, GFace *gf, bool hard)
// }
// if (NEW < 0)return false;
-
p->u = U;
p->v = V;
p->lc() = LC;
@@ -1713,10 +1693,8 @@ bool BDS_Mesh::smooth_point_centroid(BDS_Point *p, GFace *gf, bool hard)
#endif
-
bool BDS_Mesh::smooth_point_parametric(BDS_Point *const point, GFace *const gf)
{
-
if(!point->config_modified) return false;
if(point->g && point->g->classif_degree <= 1) return false;
@@ -1747,7 +1725,7 @@ bool BDS_Mesh::smooth_point_parametric(BDS_Point *const point, GFace *const gf)
it = triangles.begin();
while(it != triangles.end()) {
BDS_Face *t = *it;
- if(!test_move_point_parametric_triangle(point, U, V, t)){
+ if(!test_move_point_parametric_triangle(point, U, V, t)) {
// printf("coucou %g %g -> %g %g\n", point->u, point->v,U,V);
return false;
}
diff --git a/Mesh/BDS.h b/Mesh/BDS.h
index 65266d02fb4bb3f1b5b367548c2ec68bdbd6e486..4925b7e839aa2faff0e6bf28fb421ecf2bf52e51 100644
--- a/Mesh/BDS.h
+++ b/Mesh/BDS.h
@@ -33,11 +33,7 @@ public:
int classif_tag;
int classif_degree;
- BDS_GeomEntity(int a, int b)
- : classif_tag(a)
- , classif_degree(b)
- {
- }
+ BDS_GeomEntity(int a, int b) : classif_tag(a), classif_degree(b) {}
~BDS_GeomEntity() {}
@@ -123,9 +119,7 @@ public:
}
BDS_Vector(const BDS_Point &p2, const BDS_Point &p1);
BDS_Vector(const double X = 0., const double Y = 0., const double Z = 0.)
- : x(X)
- , y(Y)
- , z(Z)
+ : x(X), y(Y), z(Z)
{
}
static double t;
@@ -160,10 +154,13 @@ public:
{
edges.erase(std::remove(edges.begin(), edges.end(), e), edges.end());
}
- std::vector<BDS_Face*> getTriangles() const;
- BDS_Point(int id, double x=0, double y=0, double z=0)
+ std::vector<BDS_Face *> getTriangles() const;
+ BDS_Point(int id, double x = 0, double y = 0, double z = 0)
: _lcBGM(1.e22), _lcPTS(1.e22), X(x), Y(y), Z(z), u(0), v(0),
- config_modified(true), degenerated(false), /*_degeneratedTo(NULL),*/ iD(id), g(0) {}
+ config_modified(true), degenerated(false),
+ /*_degeneratedTo(NULL),*/ iD(id), g(0)
+ {
+ }
};
class BDS_Edge {
@@ -171,9 +168,7 @@ class BDS_Edge {
std::vector<BDS_Face *> _faces;
public:
- BDS_Edge(BDS_Point *A, BDS_Point *B)
- : deleted(false)
- , g(0)
+ BDS_Edge(BDS_Point *A, BDS_Point *B) : deleted(false), g(0)
{
if(*A < *B) {
p1 = A;
@@ -230,10 +225,9 @@ public:
std::bind2nd(std::equal_to<BDS_Face *>(), t)),
_faces.end());
}
- void oppositeof(BDS_Point * oface[2]) const;
- void computeNeighborhood(BDS_Point * oface[2],
- BDS_Point *t1[4],
- BDS_Point *t2[4] ) const;
+ void oppositeof(BDS_Point *oface[2]) const;
+ void computeNeighborhood(BDS_Point *oface[2], BDS_Point *t1[4],
+ BDS_Point *t2[4]) const;
void update()
{
_length = std::sqrt((p1->X - p2->X) * (p1->X - p2->X) +
@@ -250,12 +244,7 @@ public:
class BDS_Face {
public:
BDS_Face(BDS_Edge *A, BDS_Edge *B, BDS_Edge *C, BDS_Edge *D = 0)
- : deleted(false)
- , e1(A)
- , e2(B)
- , e3(C)
- , e4(D)
- , g(0)
+ : deleted(false), e1(A), e2(B), e3(C), e4(D), g(0)
{
e1->addface(this);
e2->addface(this);
@@ -292,13 +281,13 @@ public:
}
inline void getNodes(BDS_Point *_n[4]) const
{
- if (!e4){
+ if(!e4) {
_n[0] = e1->commonvertex(e3);
_n[1] = e1->commonvertex(e2);
_n[2] = e2->commonvertex(e3);
_n[3] = 0;
}
- else{
+ else {
_n[0] = e1->commonvertex(e4);
_n[1] = e1->commonvertex(e2);
_n[2] = e2->commonvertex(e3);
@@ -360,9 +349,7 @@ public:
class BDS_SwapEdgeTestRecover : public BDS_SwapEdgeTest {
public:
- BDS_SwapEdgeTestRecover()
- {
- }
+ BDS_SwapEdgeTestRecover() {}
virtual bool operator()(BDS_Point *p1, BDS_Point *p2, BDS_Point *q1,
BDS_Point *q2) const;
virtual bool operator()(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3,
@@ -372,14 +359,12 @@ public:
virtual ~BDS_SwapEdgeTestRecover() {}
};
-
class BDS_SwapEdgeTestQuality : public BDS_SwapEdgeTest {
bool testQuality, testSmallTriangles;
public:
BDS_SwapEdgeTestQuality(bool a, bool b = true)
- : testQuality(a)
- , testSmallTriangles(b)
+ : testQuality(a), testSmallTriangles(b)
{
}
virtual bool operator()(BDS_Point *p1, BDS_Point *p2, BDS_Point *q1,
@@ -391,27 +376,24 @@ public:
virtual ~BDS_SwapEdgeTestQuality() {}
};
-class BDS_SwapEdgeTestNormals : public BDS_SwapEdgeTest
-{
+class BDS_SwapEdgeTestNormals : public BDS_SwapEdgeTest {
GFace *gf;
double _ori;
- public:
- BDS_SwapEdgeTestNormals(GFace *_gf, double ori): gf(_gf), _ori(ori) {}
- virtual bool operator() (BDS_Point *p1, BDS_Point *p2,
- BDS_Point *q1, BDS_Point *q2) const;
- virtual bool operator() (BDS_Point *p1, BDS_Point *p2, BDS_Point *p3,
- BDS_Point *q1, BDS_Point *q2, BDS_Point *q3,
- BDS_Point *op1, BDS_Point *op2, BDS_Point *op3,
- BDS_Point *oq1, BDS_Point *oq2, BDS_Point *oq3) const ;
-};
+public:
+ BDS_SwapEdgeTestNormals(GFace *_gf, double ori) : gf(_gf), _ori(ori) {}
+ virtual bool operator()(BDS_Point *p1, BDS_Point *p2, BDS_Point *q1,
+ BDS_Point *q2) const;
+ virtual bool operator()(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3,
+ BDS_Point *q1, BDS_Point *q2, BDS_Point *q3,
+ BDS_Point *op1, BDS_Point *op2, BDS_Point *op3,
+ BDS_Point *oq1, BDS_Point *oq2, BDS_Point *oq3) const;
+};
-struct EdgeToRecover
-{
- int p1,p2;
+struct EdgeToRecover {
+ int p1, p2;
GEdge *ge;
- EdgeToRecover(int _p1, int _p2, GEdge *_ge)
- : ge(_ge)
+ EdgeToRecover(int _p1, int _p2, GEdge *_ge) : ge(_ge)
{
if(_p1 < _p2) {
p1 = _p1;
@@ -435,10 +417,7 @@ class BDS_Mesh {
public:
int MAXPOINTNUMBER;
double Min[3], Max[3], LC;
- BDS_Mesh(int _MAXX = 0)
- : MAXPOINTNUMBER(_MAXX)
- {
- }
+ BDS_Mesh(int _MAXX = 0) : MAXPOINTNUMBER(_MAXX) {}
virtual ~BDS_Mesh();
BDS_Mesh(const BDS_Mesh &other);
std::set<BDS_GeomEntity *, GeomLessThan> geom;
@@ -477,10 +456,12 @@ public:
BDS_Edge *recover_edge_fast(BDS_Point *p1, BDS_Point *p2);
/// Can invalidate the iterators for \p edge
- bool swap_edge(BDS_Edge*, const BDS_SwapEdgeTest &theTest, bool force = false);
- bool collapse_edge_parametric(BDS_Edge*, BDS_Point*, bool = false);
- bool smooth_point_parametric(BDS_Point * const p, GFace * const gf);
- bool smooth_point_centroid(BDS_Point *p, GFace *gf, bool test_quality=false);
+ bool swap_edge(BDS_Edge *, const BDS_SwapEdgeTest &theTest,
+ bool force = false);
+ bool collapse_edge_parametric(BDS_Edge *, BDS_Point *, bool = false);
+ bool smooth_point_parametric(BDS_Point *const p, GFace *const gf);
+ bool smooth_point_centroid(BDS_Point *p, GFace *gf,
+ bool test_quality = false);
bool split_edge(BDS_Edge *, BDS_Point *);
bool edge_constraint(BDS_Point *p1, BDS_Point *p2);
// Global operators
@@ -488,11 +469,11 @@ public:
};
void normal_triangle(BDS_Point *p1, BDS_Point *p2, BDS_Point *p3, double c[3]);
-void outputScalarField(std::vector<BDS_Face*> &t, const char *fn, int param, GFace *gf=0);
+void outputScalarField(std::vector<BDS_Face *> &t, const char *fn, int param,
+ GFace *gf = 0);
void recur_tag(BDS_Face *t, BDS_GeomEntity *g);
-int Intersect_Edges_2d(double x1, double y1, double x2, double y2,
- double x3, double y3, double x4, double y4,
- double x[2]);
-double BDS_Face_Validity (GFace *gf, BDS_Face *f);
+int Intersect_Edges_2d(double x1, double y1, double x2, double y2, double x3,
+ double y3, double x4, double y4, double x[2]);
+double BDS_Face_Validity(GFace *gf, BDS_Face *f);
#endif
diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp
index f0d2f6fd218ede050b5895948a0f7c1928f74082..8575e2f7ac2458bf3270f3f4c9a0ebd65bcd27ae 100644
--- a/Mesh/meshGFace.cpp
+++ b/Mesh/meshGFace.cpp
@@ -44,108 +44,116 @@
// define this to use the old initial delaunay
#define OLD_CODE_DELAUNAY 1
-void copyMesh(GEdge*, GEdge*, int);
+void copyMesh(GEdge *, GEdge *, int);
-void derivativeP2 (double u, double v, double dfdu[6], double dfdv[6]){
+void derivativeP2(double u, double v, double dfdu[6], double dfdv[6])
+{
/* const double sf[6][6] = { { 1, -3, -3, 4, 2, 2},
{ 0, -1, 0, 0, 2, 0},
{ 0, 0, -1, 0, 0, 2},
{ 0, 4, 0, -4, -4, 0},
{ 0, 0, 0, 4, 0, 0},
{ 0, 0, 4, -4, 0, -4} };
- */
- dfdu[0] = -3 + 4*v + 4*u;
- dfdu[1] = -1 + 4*u;
- dfdu[2] = 0 ;
- dfdu[3] = 4 -4*v - 8*u;
- dfdu[4] = 4*v ;
- dfdu[5] = -4*v ;
-
- dfdv[0] = -3 + 4*u + 4*v;
- dfdv[1] = 0 ;
- dfdv[2] = -1 + 4*v;
- dfdv[3] = -4*u ;
- dfdv[4] = 4*u ;
- dfdv[5] = 4 -4*u - 8*v;
+ */
+ dfdu[0] = -3 + 4 * v + 4 * u;
+ dfdu[1] = -1 + 4 * u;
+ dfdu[2] = 0;
+ dfdu[3] = 4 - 4 * v - 8 * u;
+ dfdu[4] = 4 * v;
+ dfdu[5] = -4 * v;
+
+ dfdv[0] = -3 + 4 * u + 4 * v;
+ dfdv[1] = 0;
+ dfdv[2] = -1 + 4 * v;
+ dfdv[3] = -4 * u;
+ dfdv[4] = 4 * u;
+ dfdv[5] = 4 - 4 * u - 8 * v;
}
-void jac_corners_p2 (double *xa, double *xb, double *xc, double *xab, double *xbc, double *xca, double J[6])
+void jac_corners_p2(double *xa, double *xb, double *xc, double *xab,
+ double *xbc, double *xca, double J[6])
{
double *x[6] = {xa, xb, xc, xab, xbc, xca};
- double nodes[6][2] = { {0, 0}, {1, 0}, {0, 1}, {0.5, 0}, {0.5, 0.5}, {0, 0.5} };
- for (int i=0; i<6; i++){
+ double nodes[6][2] = {{0, 0}, {1, 0}, {0, 1}, {0.5, 0}, {0.5, 0.5}, {0, 0.5}};
+ for(int i = 0; i < 6; i++) {
double u = nodes[i][0];
double v = nodes[i][1];
- double dfdu[6],dfdv[6];
- derivativeP2 (u, v, dfdu, dfdv);
- double dxdu = 0, dxdv = 0, dydu=0, dydv=0;
- for (int j=0; j<6; j++){
+ double dfdu[6], dfdv[6];
+ derivativeP2(u, v, dfdu, dfdv);
+ double dxdu = 0, dxdv = 0, dydu = 0, dydv = 0;
+ for(int j = 0; j < 6; j++) {
dxdu += x[j][0] * dfdu[j];
dxdv += x[j][0] * dfdv[j];
dydu += x[j][1] * dfdu[j];
dydv += x[j][1] * dfdv[j];
}
- J[i] = dxdu * dydv - dydu * dxdv;
+ J[i] = dxdu * dydv - dydu * dxdv;
}
}
-double validity_p2triangle_formula (double *xa, double *xb, double *xc, double *xab, double *xbc, double *xca)
+double validity_p2triangle_formula(double *xa, double *xb, double *xc,
+ double *xab, double *xbc, double *xca)
{
- //return 1;
+ // return 1;
// return orient2d(xa,xb,xc);
// double EPS = 1.e-3;
double J[6];
- jac_corners_p2 (xa, xb, xc, xab, xbc, xca, J);
- double bez[6] = {J[0], J[1], J[2], 2*J[3] - 0.5*(J[0]+J[1]), 2*J[4] - 0.5*(J[1]+J[2]), 2*J[5] - 0.5*(J[0]+J[2])};
+ jac_corners_p2(xa, xb, xc, xab, xbc, xca, J);
+ double bez[6] = {J[0],
+ J[1],
+ J[2],
+ 2 * J[3] - 0.5 * (J[0] + J[1]),
+ 2 * J[4] - 0.5 * (J[1] + J[2]),
+ 2 * J[5] - 0.5 * (J[0] + J[2])};
double _MIN = 1.e22;
- for (int i=0; i<6; i++)_MIN = bez[i] < _MIN ? bez[i] : _MIN;
+ for(int i = 0; i < 6; i++) _MIN = bez[i] < _MIN ? bez[i] : _MIN;
return _MIN;
}
-
-bool pointInsideParametricDomain (std::vector<SPoint2> &bnd, SPoint2 &p, SPoint2 &out, int &N){
+bool pointInsideParametricDomain(std::vector<SPoint2> &bnd, SPoint2 &p,
+ SPoint2 &out, int &N)
+{
int count = 0;
- for (size_t i=0;i<bnd.size(); i+=2){
+ for(size_t i = 0; i < bnd.size(); i += 2) {
SPoint2 p1 = bnd[i];
- SPoint2 p2 = bnd[i+1];
- double a = robustPredicates::orient2d(p1,p2,p);
- double b = robustPredicates::orient2d(p1,p2,out);
- if (a*b < 0){
- a = robustPredicates::orient2d(p,out,p1);
- b = robustPredicates::orient2d(p,out,p2);
- if (a*b < 0)count ++;
+ SPoint2 p2 = bnd[i + 1];
+ double a = robustPredicates::orient2d(p1, p2, p);
+ double b = robustPredicates::orient2d(p1, p2, out);
+ if(a * b < 0) {
+ a = robustPredicates::orient2d(p, out, p1);
+ b = robustPredicates::orient2d(p, out, p2);
+ if(a * b < 0) count++;
}
}
N = count;
- if (count % 2 == 0)return false;
+ if(count % 2 == 0) return false;
return true;
}
-
void trueBoundary(const char *iii, GFace *gf, std::vector<SPoint2> &bnd)
{
// FILE* view_t = Fopen(iii,"w");
// fprintf(view_t,"View \"True Boundary\"{\n");
std::vector<GEdge *> edg = gf->edges();
- std::set<GEdge*> edges (edg.begin(), edg.end());
-
- for (std::set<GEdge*>::iterator it = edges.begin() ; it != edges.end() ; ++it){
+ std::set<GEdge *> edges(edg.begin(), edg.end());
+
+ for(std::set<GEdge *>::iterator it = edges.begin(); it != edges.end(); ++it) {
GEdge *ge = *it;
Range<double> r = ge->parBounds(0);
- SPoint2 p[300];
+ SPoint2 p[300];
int NITER = ge->isSeam(gf) ? 2 : 1;
- for (int i=0;i<NITER;i++){
+ for(int i = 0; i < NITER; i++) {
int count = NITER == 2 ? 300 : 300;
- for (int k=0;k<count;k++){
- double t = (double)k/(count-1);
- double xi = r.low() + (r.high() - r.low()) * t;
- p[k] = ge->reparamOnFace(gf, xi, i);
- if (k > 0){
- // fprintf(view_t,"SL(%g,%g,%g,%g,%g,%g){1,1};\n",p[k-1].x(),p[k-1].y(),0.0,
- // p[k].x(),p[k].y(),0.0);
- bnd.push_back(p[k-1]);
- bnd.push_back(p[k]);
- }
+ for(int k = 0; k < count; k++) {
+ double t = (double)k / (count - 1);
+ double xi = r.low() + (r.high() - r.low()) * t;
+ p[k] = ge->reparamOnFace(gf, xi, i);
+ if(k > 0) {
+ // fprintf(view_t,"SL(%g,%g,%g,%g,%g,%g){1,1};\n",p[k-1].x(),p[k-1].y(),0.0,
+ // p[k].x(),p[k].y(),0.0);
+ bnd.push_back(p[k - 1]);
+ bnd.push_back(p[k]);
+ }
}
}
}
@@ -153,7 +161,6 @@ void trueBoundary(const char *iii, GFace *gf, std::vector<SPoint2> &bnd)
// fclose(view_t);
}
-
static void computeElementShapes(GFace *gf, double &worst, double &avg,
double &best, int &nT, int &greaterThan)
{
@@ -162,50 +169,51 @@ static void computeElementShapes(GFace *gf, double &worst, double &avg,
best = 0.0;
nT = 0;
greaterThan = 0;
- for(unsigned int i = 0; i < gf->triangles.size(); i++){
+ for(unsigned int i = 0; i < gf->triangles.size(); i++) {
double q = qmTriangle::gamma(gf->triangles[i]);
if(q > .9) greaterThan++;
avg += q;
worst = std::min(worst, q);
- best = std::max(best, q);
+ best = std::max(best, q);
nT++;
}
avg /= nT;
}
-class quadMeshRemoveHalfOfOneDMesh
-{
+class quadMeshRemoveHalfOfOneDMesh {
GFace *_gf;
- public:
- std::map<GEdge*,std::vector<MLine*> > _backup;
- std::map<MEdge, MVertex*,Less_Edge> _middle;
+
+public:
+ std::map<GEdge *, std::vector<MLine *> > _backup;
+ std::map<MEdge, MVertex *, Less_Edge> _middle;
// remove one point every two and remember middle points
- quadMeshRemoveHalfOfOneDMesh(GFace* gf) : _gf(gf)
+ quadMeshRemoveHalfOfOneDMesh(GFace *gf) : _gf(gf)
{
// only do it if a recombination has to be done
if((CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine) &&
- CTX::instance()->mesh.algoRecombine == 2){
- // printf("GFace %d removing half of the points in the 1D mesh\n",gf->tag());
- std::vector<GEdge*> const& edges = gf->edges();
- std::vector<GEdge*>::const_iterator ite = edges.begin();
- while(ite != edges.end()){
- if(!(*ite)->isMeshDegenerated()){
- std::vector<MLine*> temp;
+ CTX::instance()->mesh.algoRecombine == 2) {
+ // printf("GFace %d removing half of the points in the 1D
+ // mesh\n",gf->tag());
+ std::vector<GEdge *> const &edges = gf->edges();
+ std::vector<GEdge *>::const_iterator ite = edges.begin();
+ while(ite != edges.end()) {
+ if(!(*ite)->isMeshDegenerated()) {
+ std::vector<MLine *> temp;
(*ite)->mesh_vertices.clear();
- for(unsigned int i = 0; i< (*ite)->lines.size(); i+=2){
- if(i+1 >= (*ite)->lines.size()){
+ for(unsigned int i = 0; i < (*ite)->lines.size(); i += 2) {
+ if(i + 1 >= (*ite)->lines.size()) {
Msg::Error("1D mesh cannot be divided by 2");
break;
}
MVertex *v1 = (*ite)->lines[i]->getVertex(0);
MVertex *v2 = (*ite)->lines[i]->getVertex(1);
- MVertex *v3 = (*ite)->lines[i+1]->getVertex(1);
- v2->x() = 0.5*(v1->x()+v3->x());
- v2->y() = 0.5*(v1->y()+v3->y());
- v2->z() = 0.5*(v1->z()+v3->z());
- temp.push_back(new MLine(v1,v3));
+ MVertex *v3 = (*ite)->lines[i + 1]->getVertex(1);
+ v2->x() = 0.5 * (v1->x() + v3->x());
+ v2->y() = 0.5 * (v1->y() + v3->y());
+ v2->z() = 0.5 * (v1->z() + v3->z());
+ temp.push_back(new MLine(v1, v3));
if(v1->onWhat() == *ite) (*ite)->mesh_vertices.push_back(v1);
- _middle[MEdge(v1,v3)] = v2;
+ _middle[MEdge(v1, v3)] = v2;
}
_backup[*ite] = (*ite)->lines;
// printf("line %d goes from %d to %d\n",
@@ -214,35 +222,35 @@ class quadMeshRemoveHalfOfOneDMesh
}
++ite;
}
- CTX::instance()->mesh.lcFactor *=2.0;
+ CTX::instance()->mesh.lcFactor *= 2.0;
}
}
void subdivide()
{
- std::vector<MQuadrangle*> qnew;
+ std::vector<MQuadrangle *> qnew;
- std::map<MEdge, MVertex*, Less_Edge> eds;
+ std::map<MEdge, MVertex *, Less_Edge> eds;
- for(unsigned int i = 0; i < _gf->triangles.size(); i++){
+ for(unsigned int i = 0; i < _gf->triangles.size(); i++) {
MVertex *v[3];
SPoint2 m[3];
- for(int j = 0; j < 3; j++){
+ for(int j = 0; j < 3; j++) {
MEdge E = _gf->triangles[i]->getEdge(j);
SPoint2 p1, p2;
- reparamMeshEdgeOnFace(E.getVertex(0),E.getVertex(1),_gf,p1,p2);
+ reparamMeshEdgeOnFace(E.getVertex(0), E.getVertex(1), _gf, p1, p2);
std::map<MEdge, MVertex *, Less_Edge>::iterator it = _middle.find(E);
std::map<MEdge, MVertex *, Less_Edge>::iterator it2 = eds.find(E);
m[j] = p1;
- if(it == _middle.end() && it2 == eds.end()){
- GPoint gp = _gf->point((p1+p2)*0.5);
- double XX = 0.5*(E.getVertex(0)->x()+E.getVertex(1)->x());
- double YY = 0.5*(E.getVertex(0)->y()+E.getVertex(1)->y());
- double ZZ = 0.5*(E.getVertex(0)->z()+E.getVertex(1)->z());
- v[j] = new MFaceVertex(XX,YY,ZZ,_gf,gp.u(),gp.v());
+ if(it == _middle.end() && it2 == eds.end()) {
+ GPoint gp = _gf->point((p1 + p2) * 0.5);
+ double XX = 0.5 * (E.getVertex(0)->x() + E.getVertex(1)->x());
+ double YY = 0.5 * (E.getVertex(0)->y() + E.getVertex(1)->y());
+ double ZZ = 0.5 * (E.getVertex(0)->z() + E.getVertex(1)->z());
+ v[j] = new MFaceVertex(XX, YY, ZZ, _gf, gp.u(), gp.v());
_gf->mesh_vertices.push_back(v[j]);
eds[E] = v[j];
}
- else if(it == _middle.end()){
+ else if(it == _middle.end()) {
v[j] = it2->second;
}
else {
@@ -250,38 +258,41 @@ class quadMeshRemoveHalfOfOneDMesh
v[j]->onWhat()->mesh_vertices.push_back(v[j]);
}
}
- GPoint gp = _gf->point((m[0]+m[1]+m[2])*(1./3.));
- double XX = (v[0]->x()+v[1]->x()+v[2]->x())*(1./3.);
- double YY = (v[0]->y()+v[1]->y()+v[2]->y())*(1./3.);
- double ZZ = (v[0]->z()+v[1]->z()+v[2]->z())*(1./3.);
- MFaceVertex *vmid = new MFaceVertex(XX,YY,ZZ,_gf,gp.u(),gp.v());
+ GPoint gp = _gf->point((m[0] + m[1] + m[2]) * (1. / 3.));
+ double XX = (v[0]->x() + v[1]->x() + v[2]->x()) * (1. / 3.);
+ double YY = (v[0]->y() + v[1]->y() + v[2]->y()) * (1. / 3.);
+ double ZZ = (v[0]->z() + v[1]->z() + v[2]->z()) * (1. / 3.);
+ MFaceVertex *vmid = new MFaceVertex(XX, YY, ZZ, _gf, gp.u(), gp.v());
_gf->mesh_vertices.push_back(vmid);
- qnew.push_back(new MQuadrangle(_gf->triangles[i]->getVertex(0),v[0],vmid,v[2]));
- qnew.push_back(new MQuadrangle(_gf->triangles[i]->getVertex(1),v[1],vmid,v[0]));
- qnew.push_back(new MQuadrangle(_gf->triangles[i]->getVertex(2),v[2],vmid,v[1]));
+ qnew.push_back(
+ new MQuadrangle(_gf->triangles[i]->getVertex(0), v[0], vmid, v[2]));
+ qnew.push_back(
+ new MQuadrangle(_gf->triangles[i]->getVertex(1), v[1], vmid, v[0]));
+ qnew.push_back(
+ new MQuadrangle(_gf->triangles[i]->getVertex(2), v[2], vmid, v[1]));
delete _gf->triangles[i];
}
_gf->triangles.clear();
- for(unsigned int i = 0; i < _gf->quadrangles.size(); i++){
+ for(unsigned int i = 0; i < _gf->quadrangles.size(); i++) {
MVertex *v[4];
SPoint2 m[4];
- for(int j = 0; j < 4; j++){
+ for(int j = 0; j < 4; j++) {
MEdge E = _gf->quadrangles[i]->getEdge(j);
SPoint2 p1, p2;
- reparamMeshEdgeOnFace(E.getVertex(0),E.getVertex(1),_gf,p1,p2);
+ reparamMeshEdgeOnFace(E.getVertex(0), E.getVertex(1), _gf, p1, p2);
std::map<MEdge, MVertex *, Less_Edge>::iterator it = _middle.find(E);
std::map<MEdge, MVertex *, Less_Edge>::iterator it2 = eds.find(E);
m[j] = p1;
- if(it == _middle.end() && it2 == eds.end()){
- GPoint gp = _gf->point((p1+p2)*0.5);
- double XX = 0.5*(E.getVertex(0)->x()+E.getVertex(1)->x());
- double YY = 0.5*(E.getVertex(0)->y()+E.getVertex(1)->y());
- double ZZ = 0.5*(E.getVertex(0)->z()+E.getVertex(1)->z());
- v[j] = new MFaceVertex(XX,YY,ZZ,_gf,gp.u(),gp.v());
+ if(it == _middle.end() && it2 == eds.end()) {
+ GPoint gp = _gf->point((p1 + p2) * 0.5);
+ double XX = 0.5 * (E.getVertex(0)->x() + E.getVertex(1)->x());
+ double YY = 0.5 * (E.getVertex(0)->y() + E.getVertex(1)->y());
+ double ZZ = 0.5 * (E.getVertex(0)->z() + E.getVertex(1)->z());
+ v[j] = new MFaceVertex(XX, YY, ZZ, _gf, gp.u(), gp.v());
_gf->mesh_vertices.push_back(v[j]);
eds[E] = v[j];
}
- else if(it == _middle.end()){
+ else if(it == _middle.end()) {
v[j] = it2->second;
}
else {
@@ -289,36 +300,41 @@ class quadMeshRemoveHalfOfOneDMesh
v[j]->onWhat()->mesh_vertices.push_back(v[j]);
}
}
- GPoint gp = _gf->point((m[0]+m[1]+m[2]+m[3])*0.25);
- // FIXME : NOT EXACTLY CORRECT, BUT THAT'S THE PLACE WE WANT THE POINT TO RESIDE
- double XX = 0.25*(v[0]->x()+v[1]->x()+v[2]->x()+v[3]->x());
- double YY = 0.25*(v[0]->y()+v[1]->y()+v[2]->y()+v[3]->y());
- double ZZ = 0.25*(v[0]->z()+v[1]->z()+v[2]->z()+v[3]->z());
- MVertex *vmid = new MFaceVertex(XX,YY,ZZ,_gf,gp.u(),gp.v());
+ GPoint gp = _gf->point((m[0] + m[1] + m[2] + m[3]) * 0.25);
+ // FIXME : NOT EXACTLY CORRECT, BUT THAT'S THE PLACE WE WANT THE POINT TO
+ // RESIDE
+ double XX = 0.25 * (v[0]->x() + v[1]->x() + v[2]->x() + v[3]->x());
+ double YY = 0.25 * (v[0]->y() + v[1]->y() + v[2]->y() + v[3]->y());
+ double ZZ = 0.25 * (v[0]->z() + v[1]->z() + v[2]->z() + v[3]->z());
+ MVertex *vmid = new MFaceVertex(XX, YY, ZZ, _gf, gp.u(), gp.v());
_gf->mesh_vertices.push_back(vmid);
- qnew.push_back(new MQuadrangle(_gf->quadrangles[i]->getVertex(0),v[0],vmid,v[3]));
- qnew.push_back(new MQuadrangle(_gf->quadrangles[i]->getVertex(1),v[1],vmid,v[0]));
- qnew.push_back(new MQuadrangle(_gf->quadrangles[i]->getVertex(2),v[2],vmid,v[1]));
- qnew.push_back(new MQuadrangle(_gf->quadrangles[i]->getVertex(3),v[3],vmid,v[2]));
+ qnew.push_back(
+ new MQuadrangle(_gf->quadrangles[i]->getVertex(0), v[0], vmid, v[3]));
+ qnew.push_back(
+ new MQuadrangle(_gf->quadrangles[i]->getVertex(1), v[1], vmid, v[0]));
+ qnew.push_back(
+ new MQuadrangle(_gf->quadrangles[i]->getVertex(2), v[2], vmid, v[1]));
+ qnew.push_back(
+ new MQuadrangle(_gf->quadrangles[i]->getVertex(3), v[3], vmid, v[2]));
delete _gf->quadrangles[i];
}
_gf->quadrangles = qnew;
- // printf("%d triangles %d quads\n",_gf->triangles.size(),_gf->quadrangles.size());
+ // printf("%d triangles %d
+ // quads\n",_gf->triangles.size(),_gf->quadrangles.size());
}
void finish()
{
if((CTX::instance()->mesh.recombineAll || _gf->meshAttributes.recombine) &&
- CTX::instance()->mesh.algoRecombine == 2){
+ CTX::instance()->mesh.algoRecombine == 2) {
// recombine the elements on the half mesh
- CTX::instance()->mesh.lcFactor /=2.0;
- recombineIntoQuads(_gf,true,true,.1,true);
+ CTX::instance()->mesh.lcFactor /= 2.0;
+ recombineIntoQuads(_gf, true, true, .1, true);
// Msg::Info("subdividing");
subdivide();
// _gf->model()->writeMSH("hop2.msh");
restore();
- recombineIntoQuads(_gf,true,true,1.e-3,false);
- computeElementShapes(_gf,
- _gf->meshStatistics.worst_element_shape,
+ recombineIntoQuads(_gf, true, true, 1.e-3, false);
+ computeElementShapes(_gf, _gf->meshStatistics.worst_element_shape,
_gf->meshStatistics.average_element_shape,
_gf->meshStatistics.best_element_shape,
_gf->meshStatistics.nbTriangle,
@@ -327,11 +343,11 @@ class quadMeshRemoveHalfOfOneDMesh
}
void restore()
{
- std::vector<GEdge*> const& edges = _gf->edges();
- std::vector<GEdge*>::const_iterator ite = edges.begin();
- while(ite != edges.end()){
- for(unsigned int i = 0; i < (*ite)->lines.size(); i++){
- delete (*ite)->lines[i];
+ std::vector<GEdge *> const &edges = _gf->edges();
+ std::vector<GEdge *>::const_iterator ite = edges.begin();
+ while(ite != edges.end()) {
+ for(unsigned int i = 0; i < (*ite)->lines.size(); i++) {
+ delete(*ite)->lines[i];
}
(*ite)->lines = _backup[*ite];
++ite;
@@ -341,45 +357,47 @@ class quadMeshRemoveHalfOfOneDMesh
static void copyMesh(GFace *source, GFace *target)
{
- std::map<MVertex*, MVertex*> vs2vt;
+ std::map<MVertex *, MVertex *> vs2vt;
// add principal vertex pairs
- std::vector<GVertex*> const& s_vtcs = source->vertices();
- std::vector<GVertex*> const& t_vtcs = target->vertices();
+ std::vector<GVertex *> const &s_vtcs = source->vertices();
+ std::vector<GVertex *> const &t_vtcs = target->vertices();
if(s_vtcs.size() != t_vtcs.size()) {
Msg::Info("Periodicity imposed on topologically incompatible surfaces"
- "(%d vs %d bounding vertices)", s_vtcs.size(), t_vtcs.size());
+ "(%d vs %d bounding vertices)",
+ s_vtcs.size(), t_vtcs.size());
}
- std::set<GVertex*> checkVtcs(s_vtcs.begin(), s_vtcs.end());
+ std::set<GVertex *> checkVtcs(s_vtcs.begin(), s_vtcs.end());
- for(std::vector<GVertex*>::const_iterator tvIter = t_vtcs.begin();
+ for(std::vector<GVertex *>::const_iterator tvIter = t_vtcs.begin();
tvIter != t_vtcs.end(); ++tvIter) {
-
GVertex *gvt = *tvIter;
- std::map<GVertex*, GVertex*>::iterator gvsIter = target->vertexCounterparts.find(gvt);
+ std::map<GVertex *, GVertex *>::iterator gvsIter =
+ target->vertexCounterparts.find(gvt);
if(gvsIter == target->vertexCounterparts.end()) {
Msg::Error("Periodic meshing of surface %d with surface %d: "
"point %d has no periodic counterpart",
target->tag(), source->tag(), gvt->tag());
}
- else{
+ else {
GVertex *gvs = gvsIter->second;
if(checkVtcs.find(gvs) == checkVtcs.end()) {
if(gvs)
- Msg::Error("Periodic meshing of surface %d with surface %d: "
- "point %d has periodic counterpart %d outside of source surface",
- target->tag(), source->tag(), gvt->tag(), gvs->tag());
+ Msg::Error(
+ "Periodic meshing of surface %d with surface %d: "
+ "point %d has periodic counterpart %d outside of source surface",
+ target->tag(), source->tag(), gvt->tag(), gvs->tag());
else
Msg::Error("Periodic meshing of surface %d with surface %d: "
"point %d has no periodic counterpart",
target->tag(), source->tag(), gvt->tag());
}
- if(gvs){
+ if(gvs) {
MVertex *vs = gvs->mesh_vertices[0];
MVertex *vt = gvt->mesh_vertices[0];
vs2vt[vs] = vt;
@@ -390,43 +408,44 @@ static void copyMesh(GFace *source, GFace *target)
// add corresponding edge nodes assuming edges were correctly meshed already
- std::vector<GEdge*> s_edges = source->edges();
- std::vector<GEdge*> t_edges = target->edges();
+ std::vector<GEdge *> s_edges = source->edges();
+ std::vector<GEdge *> t_edges = target->edges();
- std::set<GEdge*> checkEdges;
- checkEdges.insert(s_edges.begin(),s_edges.end());
+ std::set<GEdge *> checkEdges;
+ checkEdges.insert(s_edges.begin(), s_edges.end());
- for(std::vector<GEdge*>::iterator te_iter = t_edges.begin();
+ for(std::vector<GEdge *>::iterator te_iter = t_edges.begin();
te_iter != t_edges.end(); ++te_iter) {
+ GEdge *get = *te_iter;
- GEdge* get = *te_iter;
-
- std::map<GEdge*, std::pair<GEdge*, int> >::iterator gesIter =
+ std::map<GEdge *, std::pair<GEdge *, int> >::iterator gesIter =
target->edgeCounterparts.find(get);
if(gesIter == target->edgeCounterparts.end()) {
Msg::Error("Periodic meshing of surface %d with surface %d: "
"curve %d has no periodic counterpart",
target->tag(), source->tag(), get->tag());
}
- else{
- GEdge* ges = gesIter->second.first;
+ else {
+ GEdge *ges = gesIter->second.first;
if(checkEdges.find(ges) == checkEdges.end()) {
- Msg::Error("Periodic meshing of surface %d with surface %d: "
- "curve %d has periodic counterpart %d outside of get surface",
- target->tag(), source->tag(), get->tag(), ges->tag());
+ Msg::Error(
+ "Periodic meshing of surface %d with surface %d: "
+ "curve %d has periodic counterpart %d outside of get surface",
+ target->tag(), source->tag(), get->tag(), ges->tag());
}
if(get->mesh_vertices.size() != ges->mesh_vertices.size()) {
Msg::Error("Periodic meshing of surface %d with surface %d: "
"curve %d has %d vertices, whereas correspondant %d has %d",
- target->tag(), source->tag(),
- get->tag(), get->mesh_vertices.size(),
- ges->tag(), ges->mesh_vertices.size());
+ target->tag(), source->tag(), get->tag(),
+ get->mesh_vertices.size(), ges->tag(),
+ ges->mesh_vertices.size());
}
int orientation = gesIter->second.second;
- int is = orientation == 1 ? 0 : get->mesh_vertices.size()-1;
- for(unsigned it=0;it<get->mesh_vertices.size();it++,is+=orientation) {
- MVertex* vs = ges->mesh_vertices[is];
- MVertex* vt = get->mesh_vertices[it];
+ int is = orientation == 1 ? 0 : get->mesh_vertices.size() - 1;
+ for(unsigned it = 0; it < get->mesh_vertices.size();
+ it++, is += orientation) {
+ MVertex *vs = ges->mesh_vertices[is];
+ MVertex *vt = get->mesh_vertices[it];
vs2vt[vs] = vt;
target->correspondingVertices[vt] = vs;
}
@@ -434,9 +453,9 @@ static void copyMesh(GFace *source, GFace *target)
}
// now transform
- std::vector<double>& tfo = target->affineTransform;
+ std::vector<double> &tfo = target->affineTransform;
- for(unsigned int i = 0; i < source->mesh_vertices.size(); i++){
+ for(unsigned int i = 0; i < source->mesh_vertices.size(); i++) {
MVertex *vs = source->mesh_vertices[i];
SPoint2 XXX;
@@ -444,28 +463,29 @@ static void copyMesh(GFace *source, GFace *target)
double res[4] = {0., 0., 0., 0.};
int idx = 0;
for(int i = 0; i < 4; i++)
- for(int j = 0; j < 4; j++)
- res[i] += tfo[idx++] * ps[j];
+ for(int j = 0; j < 4; j++) res[i] += tfo[idx++] * ps[j];
SPoint3 tp(res[0], res[1], res[2]);
XXX = target->parFromPoint(tp);
GPoint gp = target->point(XXX);
- MVertex *vt = new MFaceVertex(gp.x(), gp.y(), gp.z(), target, gp.u(), gp.v());
+ MVertex *vt =
+ new MFaceVertex(gp.x(), gp.y(), gp.z(), target, gp.u(), gp.v());
target->mesh_vertices.push_back(vt);
target->correspondingVertices[vt] = vs;
vs2vt[vs] = vt;
}
- for(unsigned i = 0; i < source->triangles.size(); i++){
+ for(unsigned i = 0; i < source->triangles.size(); i++) {
MVertex *vt[3];
- for(int j = 0; j < 3; j++){
+ for(int j = 0; j < 3; j++) {
MVertex *vs = source->triangles[i]->getVertex(j);
vt[j] = vs2vt[vs];
}
- if(!vt[0] || !vt[1] ||!vt[2]){
- Msg::Error("Problem in mesh copying procedure %p %p %p %d %d %d",
- vt[0], vt[1], vt[2], source->triangles[i]->getVertex(0)->onWhat()->dim(),
+ if(!vt[0] || !vt[1] || !vt[2]) {
+ Msg::Error("Problem in mesh copying procedure %p %p %p %d %d %d", vt[0],
+ vt[1], vt[2],
+ source->triangles[i]->getVertex(0)->onWhat()->dim(),
source->triangles[i]->getVertex(1)->onWhat()->dim(),
source->triangles[i]->getVertex(2)->onWhat()->dim());
return;
@@ -473,12 +493,12 @@ static void copyMesh(GFace *source, GFace *target)
target->triangles.push_back(new MTriangle(vt[0], vt[1], vt[2]));
}
- for(unsigned i = 0; i < source->quadrangles.size(); i++){
+ for(unsigned i = 0; i < source->quadrangles.size(); i++) {
MVertex *v1 = vs2vt[source->quadrangles[i]->getVertex(0)];
MVertex *v2 = vs2vt[source->quadrangles[i]->getVertex(1)];
MVertex *v3 = vs2vt[source->quadrangles[i]->getVertex(2)];
MVertex *v4 = vs2vt[source->quadrangles[i]->getVertex(3)];
- if(!v1 || !v2 || !v3 || !v4){
+ if(!v1 || !v2 || !v3 || !v4) {
Msg::Error("Problem in mesh copying procedure %p %p %p %p %d %d %d %d",
v1, v2, v3, v4,
source->quadrangles[i]->getVertex(0)->onWhat()->dim(),
@@ -496,30 +516,32 @@ void fourthPoint(double *p1, double *p2, double *p3, double *p4)
circumCenterXYZ(p1, p2, p3, c);
double vx[3] = {p2[0] - p1[0], p2[1] - p1[1], p2[2] - p1[2]};
double vy[3] = {p3[0] - p1[0], p3[1] - p1[1], p3[2] - p1[2]};
- double vz[3]; prodve(vx, vy, vz);
+ double vz[3];
+ prodve(vx, vy, vz);
norme(vz);
- double R = sqrt((p1[0] - c[0]) * (p1[0] - c[0]) +
- (p1[1] - c[1]) * (p1[1] - c[1]) +
- (p1[2] - c[2]) * (p1[2] - c[2]));
+ double R =
+ sqrt((p1[0] - c[0]) * (p1[0] - c[0]) + (p1[1] - c[1]) * (p1[1] - c[1]) +
+ (p1[2] - c[2]) * (p1[2] - c[2]));
p4[0] = c[0] + R * vz[0];
p4[1] = c[1] + R * vz[1];
p4[2] = c[2] + R * vz[2];
}
-
-static void remeshUnrecoveredEdges(std::map<MVertex*, BDS_Point*> &recoverMapInv,
- std::set<EdgeToRecover> &edgesNotRecovered,
- std::list<GFace*> &facesToRemesh)
+static void
+remeshUnrecoveredEdges(std::map<MVertex *, BDS_Point *> &recoverMapInv,
+ std::set<EdgeToRecover> &edgesNotRecovered,
+ std::list<GFace *> &facesToRemesh)
{
facesToRemesh.clear();
deMeshGFace dem;
std::set<EdgeToRecover>::iterator itr = edgesNotRecovered.begin();
- for(; itr != edgesNotRecovered.end(); ++itr){
- std::vector<GFace*> l_faces = itr->ge->faces();
+ for(; itr != edgesNotRecovered.end(); ++itr) {
+ std::vector<GFace *> l_faces = itr->ge->faces();
// un-mesh model faces adjacent to the model edge
- for(std::vector<GFace*>::iterator it = l_faces.begin(); it != l_faces.end(); ++it){
- if((*it)->triangles.size() || (*it)->quadrangles.size()){
+ for(std::vector<GFace *>::iterator it = l_faces.begin();
+ it != l_faces.end(); ++it) {
+ if((*it)->triangles.size() || (*it)->quadrangles.size()) {
facesToRemesh.push_back(*it);
dem(*it);
}
@@ -533,41 +555,48 @@ static void remeshUnrecoveredEdges(std::map<MVertex*, BDS_Point*> &recoverMapInv
GVertex *g2 = itr->ge->getEndVertex();
Range<double> bb = itr->ge->parBounds(0);
- std::vector<MLine*> newLines;
+ std::vector<MLine *> newLines;
- for(int i = 0; i < N; i++){
+ for(int i = 0; i < N; i++) {
MVertex *v1 = itr->ge->lines[i]->getVertex(0);
MVertex *v2 = itr->ge->lines[i]->getVertex(1);
- std::map<MVertex*, BDS_Point*>::iterator itp1 = recoverMapInv.find(v1);
- std::map<MVertex*, BDS_Point*>::iterator itp2 = recoverMapInv.find(v2);
- if(itp1 != recoverMapInv.end() && itp2 != recoverMapInv.end()){
+ std::map<MVertex *, BDS_Point *>::iterator itp1 = recoverMapInv.find(v1);
+ std::map<MVertex *, BDS_Point *>::iterator itp2 = recoverMapInv.find(v2);
+ if(itp1 != recoverMapInv.end() && itp2 != recoverMapInv.end()) {
BDS_Point *pp1 = itp1->second;
BDS_Point *pp2 = itp2->second;
- if((pp1->iD == p1 && pp2->iD == p2) || (pp1->iD == p2 && pp2->iD == p1)){
+ if((pp1->iD == p1 && pp2->iD == p2) ||
+ (pp1->iD == p2 && pp2->iD == p1)) {
double t1;
double lc1 = -1;
- if(v1->onWhat() == g1) t1 = bb.low();
- else if(v1->onWhat() == g2) t1 = bb.high();
+ if(v1->onWhat() == g1)
+ t1 = bb.low();
+ else if(v1->onWhat() == g2)
+ t1 = bb.high();
else {
- MEdgeVertex *ev1 = (MEdgeVertex*)v1;
+ MEdgeVertex *ev1 = (MEdgeVertex *)v1;
lc1 = ev1->getLc();
v1->getParameter(0, t1);
}
double t2;
double lc2 = -1;
- if(v2->onWhat() == g1) t2 = bb.low();
- else if(v2->onWhat() == g2) t2 = bb.high();
+ if(v2->onWhat() == g1)
+ t2 = bb.low();
+ else if(v2->onWhat() == g2)
+ t2 = bb.high();
else {
- MEdgeVertex *ev2 = (MEdgeVertex*)v2;
+ MEdgeVertex *ev2 = (MEdgeVertex *)v2;
lc2 = ev2->getLc();
v2->getParameter(0, t2);
}
// periodic
if(v1->onWhat() == g1 && v1->onWhat() == g2)
- t1 = fabs(t2-bb.low()) < fabs(t2-bb.high()) ? bb.low() : bb.high();
+ t1 =
+ fabs(t2 - bb.low()) < fabs(t2 - bb.high()) ? bb.low() : bb.high();
if(v2->onWhat() == g1 && v2->onWhat() == g2)
- t2 = fabs(t1-bb.low()) < fabs(t1-bb.high()) ? bb.low() : bb.high();
+ t2 =
+ fabs(t1 - bb.low()) < fabs(t1 - bb.high()) ? bb.low() : bb.high();
if(lc1 == -1)
lc1 = BGM_MeshSize(v1->onWhat(), 0, 0, v1->x(), v1->y(), v1->z());
@@ -577,12 +606,13 @@ static void remeshUnrecoveredEdges(std::map<MVertex*, BDS_Point*> &recoverMapInv
double t = 0.5 * (t2 + t1);
double lc = 0.5 * (lc1 + lc2);
GPoint V = itr->ge->point(t);
- MEdgeVertex * newv = new MEdgeVertex(V.x(), V.y(), V.z(), itr->ge, t, 0, lc);
+ MEdgeVertex *newv =
+ new MEdgeVertex(V.x(), V.y(), V.z(), itr->ge, t, 0, lc);
newLines.push_back(new MLine(v1, newv));
newLines.push_back(new MLine(newv, v2));
delete itr->ge->lines[i];
}
- else{
+ else {
newLines.push_back(itr->ge->lines[i]);
}
}
@@ -594,7 +624,7 @@ static void remeshUnrecoveredEdges(std::map<MVertex*, BDS_Point*> &recoverMapInv
itr->ge->lines = newLines;
itr->ge->mesh_vertices.clear();
N = itr->ge->lines.size();
- for(int i = 1; i < N; i++){
+ for(int i = 1; i < N; i++) {
itr->ge->mesh_vertices.push_back(itr->ge->lines[i]->getVertex(0));
}
}
@@ -602,7 +632,6 @@ static void remeshUnrecoveredEdges(std::map<MVertex*, BDS_Point*> &recoverMapInv
static bool algoDelaunay2D(GFace *gf)
{
-
if(gf->getMeshingAlgo() == ALGO_2D_DELAUNAY ||
gf->getMeshingAlgo() == ALGO_2D_BAMG ||
gf->getMeshingAlgo() == ALGO_2D_FRONTAL ||
@@ -619,59 +648,66 @@ static bool algoDelaunay2D(GFace *gf)
}
static bool recoverEdge(BDS_Mesh *m, GEdge *ge,
- std::map<MVertex*, BDS_Point*> &recoverMapInv,
+ std::map<MVertex *, BDS_Point *> &recoverMapInv,
std::set<EdgeToRecover> *e2r,
std::set<EdgeToRecover> *notRecovered, int pass)
{
BDS_GeomEntity *g = 0;
- if(pass == 2){
+ if(pass == 2) {
m->add_geom(ge->tag(), 1);
g = m->get_geom(ge->tag(), 1);
}
bool _fatallyFailed;
- for(unsigned int i = 0; i < ge->lines.size(); i++){
+ for(unsigned int i = 0; i < ge->lines.size(); i++) {
MVertex *vstart = ge->lines[i]->getVertex(0);
MVertex *vend = ge->lines[i]->getVertex(1);
- std::map<MVertex*, BDS_Point*>::iterator itpstart = recoverMapInv.find(vstart);
- std::map<MVertex*, BDS_Point*>::iterator itpend = recoverMapInv.find(vend);
- if(itpstart != recoverMapInv.end() && itpend != recoverMapInv.end()){
+ std::map<MVertex *, BDS_Point *>::iterator itpstart =
+ recoverMapInv.find(vstart);
+ std::map<MVertex *, BDS_Point *>::iterator itpend =
+ recoverMapInv.find(vend);
+ if(itpstart != recoverMapInv.end() && itpend != recoverMapInv.end()) {
BDS_Point *pstart = itpstart->second;
BDS_Point *pend = itpend->second;
if(pass == 1)
e2r->insert(EdgeToRecover(pstart->iD, pend->iD, ge));
- else{
- BDS_Edge *e = m->recover_edge(pstart->iD, pend->iD, _fatallyFailed, e2r, notRecovered);
- if(e) e->g = g;
+ else {
+ BDS_Edge *e = m->recover_edge(pstart->iD, pend->iD, _fatallyFailed, e2r,
+ notRecovered);
+ if(e)
+ e->g = g;
else {
- if(_fatallyFailed){
- Msg::Error("Unable to recover the edge %d (%d/%d) on GEdge %d (on GFace %d)",
- ge->lines[i]->getNum(), i+1, ge->lines.size(), ge->tag(),
- ge->faces().back()->tag());
- //outputScalarField(m->triangles, "wrongmesh.pos", 0);
- //outputScalarField(m->triangles, "wrongparam.pos", 1);
- }
+ if(_fatallyFailed) {
+ Msg::Error(
+ "Unable to recover the edge %d (%d/%d) on GEdge %d (on GFace %d)",
+ ge->lines[i]->getNum(), i + 1, ge->lines.size(), ge->tag(),
+ ge->faces().back()->tag());
+ // outputScalarField(m->triangles, "wrongmesh.pos", 0);
+ // outputScalarField(m->triangles, "wrongparam.pos", 1);
+ }
return !_fatallyFailed;
}
}
}
}
- if(pass == 2 && ge->getBeginVertex()){
+ if(pass == 2 && ge->getBeginVertex()) {
MVertex *vstart = *(ge->getBeginVertex()->mesh_vertices.begin());
MVertex *vend = *(ge->getEndVertex()->mesh_vertices.begin());
- std::map<MVertex*, BDS_Point*>::iterator itpstart = recoverMapInv.find(vstart);
- std::map<MVertex*, BDS_Point*>::iterator itpend = recoverMapInv.find(vend);
- if(itpstart != recoverMapInv.end() && itpend != recoverMapInv.end()){
+ std::map<MVertex *, BDS_Point *>::iterator itpstart =
+ recoverMapInv.find(vstart);
+ std::map<MVertex *, BDS_Point *>::iterator itpend =
+ recoverMapInv.find(vend);
+ if(itpstart != recoverMapInv.end() && itpend != recoverMapInv.end()) {
BDS_Point *pstart = itpstart->second;
BDS_Point *pend = itpend->second;
- if(!pstart->g){
+ if(!pstart->g) {
m->add_geom(pstart->iD, 0);
BDS_GeomEntity *g0 = m->get_geom(pstart->iD, 0);
pstart->g = g0;
}
- if(!pend->g){
+ if(!pend->g) {
m->add_geom(pend->iD, 0);
BDS_GeomEntity *g0 = m->get_geom(pend->iD, 0);
pend->g = g0;
@@ -683,15 +719,14 @@ static bool recoverEdge(BDS_Mesh *m, GEdge *ge,
}
void BDS2GMSH(BDS_Mesh *m, GFace *gf,
- std::map<BDS_Point*, MVertex*, PointLessThan> &recoverMap)
+ std::map<BDS_Point *, MVertex *, PointLessThan> &recoverMap)
{
{
- std::set<BDS_Point*,PointLessThan>::iterator itp = m->points.begin();
- while(itp != m->points.end()){
+ std::set<BDS_Point *, PointLessThan>::iterator itp = m->points.begin();
+ while(itp != m->points.end()) {
BDS_Point *p = *itp;
- if(recoverMap.find(p) == recoverMap.end()){
- MVertex *v = new MFaceVertex
- (p->X, p->Y, p->Z, gf, p->u, p->v);
+ if(recoverMap.find(p) == recoverMap.end()) {
+ MVertex *v = new MFaceVertex(p->X, p->Y, p->Z, gf, p->u, p->v);
recoverMap[p] = v;
gf->mesh_vertices.push_back(v);
}
@@ -699,22 +734,22 @@ void BDS2GMSH(BDS_Mesh *m, GFace *gf,
}
}
{
- std::vector<BDS_Face*>::iterator itt = m->triangles.begin();
- while(itt != m->triangles.end()){
+ std::vector<BDS_Face *>::iterator itt = m->triangles.begin();
+ while(itt != m->triangles.end()) {
BDS_Face *t = *itt;
- if(!t->deleted){
+ if(!t->deleted) {
BDS_Point *n[4];
t->getNodes(n);
MVertex *v1 = recoverMap[n[0]];
MVertex *v2 = recoverMap[n[1]];
MVertex *v3 = recoverMap[n[2]];
- if(!n[3]){
+ if(!n[3]) {
// when a singular point is present, degenerated triangles may be
// created, for example on a sphere that contains one pole
if(v1 != v2 && v1 != v3 && v2 != v3)
gf->triangles.push_back(new MTriangle(v1, v2, v3));
}
- else{
+ else {
MVertex *v4 = recoverMap[n[3]];
gf->quadrangles.push_back(new MQuadrangle(v1, v2, v3, v4));
}
@@ -724,195 +759,204 @@ void BDS2GMSH(BDS_Mesh *m, GFace *gf,
}
}
-static void addOrRemove(MVertex *v1, MVertex *v2, std::set<MEdge,Less_Edge> & bedges,
- std::set<MEdge,Less_Edge> &removed)
+static void addOrRemove(MVertex *v1, MVertex *v2,
+ std::set<MEdge, Less_Edge> &bedges,
+ std::set<MEdge, Less_Edge> &removed)
{
- MEdge e(v1,v2);
- if(removed.find(e) != removed.end())return;
- std::set<MEdge,Less_Edge>::iterator it = bedges.find(e);
- if(it == bedges.end())bedges.insert(e);
+ MEdge e(v1, v2);
+ if(removed.find(e) != removed.end()) return;
+ std::set<MEdge, Less_Edge>::iterator it = bedges.find(e);
+ if(it == bedges.end())
+ bedges.insert(e);
else {
bedges.erase(it);
removed.insert(e);
}
}
-static void modifyInitialMeshForBoundaryLayers(GFace *gf,
- std::vector<MQuadrangle*> &blQuads,
- std::vector<MTriangle*> &blTris,
- std::set<MVertex*> &verts,
- bool debug)
+static void modifyInitialMeshForBoundaryLayers(
+ GFace *gf, std::vector<MQuadrangle *> &blQuads,
+ std::vector<MTriangle *> &blTris, std::set<MVertex *> &verts, bool debug)
{
if(!buildAdditionalPoints2D(gf)) return;
- BoundaryLayerColumns* _columns = gf->getColumns();
+ BoundaryLayerColumns *_columns = gf->getColumns();
- std::set<MEdge,Less_Edge> bedges;
- std::set<MEdge,Less_Edge> removed;
+ std::set<MEdge, Less_Edge> bedges;
+ std::set<MEdge, Less_Edge> removed;
- std::vector<GEdge*> edges = gf->edges();
- std::vector<GEdge*> const& embedded_edges = gf->embeddedEdges();
+ std::vector<GEdge *> edges = gf->edges();
+ std::vector<GEdge *> const &embedded_edges = gf->embeddedEdges();
edges.insert(edges.begin(), embedded_edges.begin(), embedded_edges.end());
- std::vector<GEdge*>::iterator ite = edges.begin();
+ std::vector<GEdge *>::iterator ite = edges.begin();
FILE *ff2 = 0;
- if(debug) ff2 = Fopen("tato.pos","w");
- if(ff2) fprintf(ff2,"View \" \"{\n");
+ if(debug) ff2 = Fopen("tato.pos", "w");
+ if(ff2) fprintf(ff2, "View \" \"{\n");
- std::vector<MLine*> _lines;
+ std::vector<MLine *> _lines;
- while(ite != edges.end()){
- for(unsigned int i = 0; i< (*ite)->lines.size(); i++){
+ while(ite != edges.end()) {
+ for(unsigned int i = 0; i < (*ite)->lines.size(); i++) {
_lines.push_back((*ite)->lines[i]);
MVertex *v1 = (*ite)->lines[i]->getVertex(0);
MVertex *v2 = (*ite)->lines[i]->getVertex(1);
- MEdge dv(v1,v2);
- addOrRemove(v1,v2,bedges, removed);
- for(unsigned int SIDE = 0 ; SIDE < _columns->_normals.count(dv); SIDE ++){
- std::vector<MElement*> myCol;
+ MEdge dv(v1, v2);
+ addOrRemove(v1, v2, bedges, removed);
+ for(unsigned int SIDE = 0; SIDE < _columns->_normals.count(dv); SIDE++) {
+ std::vector<MElement *> myCol;
edgeColumn ec = _columns->getColumns(v1, v2, SIDE);
- const BoundaryLayerData & c1 = ec._c1;
- const BoundaryLayerData & c2 = ec._c2;
+ const BoundaryLayerData &c1 = ec._c1;
+ const BoundaryLayerData &c2 = ec._c2;
int N = std::min(c1._column.size(), c2._column.size());
if(!N) continue;
- for(int l = 0; l < N; ++l){
+ for(int l = 0; l < N; ++l) {
MVertex *v11, *v12, *v21, *v22;
v21 = c1._column[l];
v22 = c2._column[l];
- if(l == 0){
+ if(l == 0) {
v11 = v1;
v12 = v2;
}
else {
- v11 = c1._column[l-1];
- v12 = c2._column[l-1];
+ v11 = c1._column[l - 1];
+ v12 = c2._column[l - 1];
}
MQuadrangle *qq = new MQuadrangle(v11, v21, v22, v12);
- qq->setPartition(l+1);
+ qq->setPartition(l + 1);
myCol.push_back(qq);
blQuads.push_back(qq);
if(ff2)
- fprintf(ff2,"SQ(%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){%d,%d,%d,%d};\n",
- v11->x(),v11->y(),v11->z(),
- v12->x(),v12->y(),v12->z(),
- v22->x(),v22->y(),v22->z(),
- v21->x(),v21->y(),v21->z(),l+1,l+1,l+1,l+1);
+ fprintf(ff2,
+ "SQ(%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){%d,%d,%d,%d};\n",
+ v11->x(), v11->y(), v11->z(), v12->x(), v12->y(), v12->z(),
+ v22->x(), v22->y(), v22->z(), v21->x(), v21->y(), v21->z(),
+ l + 1, l + 1, l + 1, l + 1);
}
if(c1._column.size() != c2._column.size()) {
- MVertex *v11,*v12,*v;
- v11 = c1._column[N-1];
- v12 = c2._column[N-1];
- v = c1._column.size() > c2._column.size() ?
- c1._column[N] : c2._column[N];
+ MVertex *v11, *v12, *v;
+ v11 = c1._column[N - 1];
+ v12 = c2._column[N - 1];
+ v = c1._column.size() > c2._column.size() ? c1._column[N] :
+ c2._column[N];
MTriangle *qq = new MTriangle(v11, v12, v);
- qq->setPartition(N+1);
+ qq->setPartition(N + 1);
myCol.push_back(qq);
blTris.push_back(qq);
if(ff2)
- fprintf(ff2,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%d,%d,%d};\n",
- v11->x(),v11->y(),v11->z(),
- v12->x(),v12->y(),v12->z(),
- v->x(),v->y(),v->z(),N+1,N+1,N+1);
+ fprintf(ff2, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%d,%d,%d};\n",
+ v11->x(), v11->y(), v11->z(), v12->x(), v12->y(), v12->z(),
+ v->x(), v->y(), v->z(), N + 1, N + 1, N + 1);
}
// int M = std::max(c1._column.size(),c2._column.size());
- for(unsigned int l = 0; l <myCol.size(); l++)
+ for(unsigned int l = 0; l < myCol.size(); l++)
_columns->_toFirst[myCol[l]] = myCol[0];
_columns->_elemColumns[myCol[0]] = myCol;
}
- }
+ }
++ite;
}
for(BoundaryLayerColumns::iterf itf = _columns->beginf();
- itf != _columns->endf() ; ++itf){
+ itf != _columns->endf(); ++itf) {
MVertex *v = itf->first;
int nbCol = _columns->getNbColumns(v);
- for(int i = 0; i < nbCol - 1; i++){
- const BoundaryLayerData & c1 = _columns->getColumn(v,i);
- const BoundaryLayerData & c2 = _columns->getColumn(v,i+1);
- int N = std::min(c1._column.size(),c2._column.size());
+ for(int i = 0; i < nbCol - 1; i++) {
+ const BoundaryLayerData &c1 = _columns->getColumn(v, i);
+ const BoundaryLayerData &c2 = _columns->getColumn(v, i + 1);
+ int N = std::min(c1._column.size(), c2._column.size());
// printf("%d %d\n",c1._column.size(),c2._column.size());
- std::vector<MElement*> myCol;
- for(int l=0;l < N ;++l){
- MVertex *v11,*v12,*v21,*v22;
+ std::vector<MElement *> myCol;
+ for(int l = 0; l < N; ++l) {
+ MVertex *v11, *v12, *v21, *v22;
v21 = c1._column[l];
v22 = c2._column[l];
- if(l == 0){
+ if(l == 0) {
v11 = v;
v12 = v;
- }
+ }
else {
- v11 = c1._column[l-1];
- v12 = c2._column[l-1];
+ v11 = c1._column[l - 1];
+ v12 = c2._column[l - 1];
}
- if(v11 != v12){
- MQuadrangle *qq = new MQuadrangle(v11,v12,v22,v21);
- qq->setPartition(l+1);
+ if(v11 != v12) {
+ MQuadrangle *qq = new MQuadrangle(v11, v12, v22, v21);
+ qq->setPartition(l + 1);
myCol.push_back(qq);
blQuads.push_back(qq);
if(ff2)
- fprintf(ff2,"SQ(%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){%d,%d,%d,%d};\n",
- v11->x(),v11->y(),v11->z(),
- v12->x(),v12->y(),v12->z(),
- v22->x(),v22->y(),v22->z(),
- v21->x(),v21->y(),v21->z(),l+1,l+1,l+1,l+1);
+ fprintf(ff2,
+ "SQ(%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g){%d,%d,%d,%d};\n",
+ v11->x(), v11->y(), v11->z(), v12->x(), v12->y(), v12->z(),
+ v22->x(), v22->y(), v22->z(), v21->x(), v21->y(), v21->z(),
+ l + 1, l + 1, l + 1, l + 1);
}
else {
- MTriangle *qq = new MTriangle(v,v22,v21);
- qq->setPartition(l+1);
+ MTriangle *qq = new MTriangle(v, v22, v21);
+ qq->setPartition(l + 1);
myCol.push_back(qq);
blTris.push_back(qq);
if(ff2)
- fprintf(ff2,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%d,%d,%d};\n",
- v->x(),v->y(),v->z(),
- v22->x(),v22->y(),v22->z(),
- v21->x(),v21->y(),v21->z(),l+1,l+1,l+1);
+ fprintf(ff2, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g){%d,%d,%d};\n", v->x(),
+ v->y(), v->z(), v22->x(), v22->y(), v22->z(), v21->x(),
+ v21->y(), v21->z(), l + 1, l + 1, l + 1);
}
}
- for(unsigned int l=0;l<myCol.size();l++)_columns->_toFirst[myCol[l]] = myCol[0];
+ for(unsigned int l = 0; l < myCol.size(); l++)
+ _columns->_toFirst[myCol[l]] = myCol[0];
_columns->_elemColumns[myCol[0]] = myCol;
}
}
- if(ff2){
- fprintf(ff2,"};\n");
+ if(ff2) {
+ fprintf(ff2, "};\n");
fclose(ff2);
}
- filterOverlappingElements(_lines, blTris,blQuads,_columns->_elemColumns,_columns->_toFirst);
+ filterOverlappingElements(_lines, blTris, blQuads, _columns->_elemColumns,
+ _columns->_toFirst);
- for(unsigned int i = 0; i < blQuads.size();i++){
- addOrRemove(blQuads[i]->getVertex(0),blQuads[i]->getVertex(1),bedges, removed);
- addOrRemove(blQuads[i]->getVertex(1),blQuads[i]->getVertex(2),bedges, removed);
- addOrRemove(blQuads[i]->getVertex(2),blQuads[i]->getVertex(3),bedges, removed);
- addOrRemove(blQuads[i]->getVertex(3),blQuads[i]->getVertex(0),bedges, removed);
+ for(unsigned int i = 0; i < blQuads.size(); i++) {
+ addOrRemove(blQuads[i]->getVertex(0), blQuads[i]->getVertex(1), bedges,
+ removed);
+ addOrRemove(blQuads[i]->getVertex(1), blQuads[i]->getVertex(2), bedges,
+ removed);
+ addOrRemove(blQuads[i]->getVertex(2), blQuads[i]->getVertex(3), bedges,
+ removed);
+ addOrRemove(blQuads[i]->getVertex(3), blQuads[i]->getVertex(0), bedges,
+ removed);
for(int j = 0; j < 4; j++)
- if(blQuads[i]->getVertex(j)->onWhat() == gf)verts.insert(blQuads[i]->getVertex(j));
- }
- for(unsigned int i = 0; i < blTris.size(); i++){
- addOrRemove(blTris[i]->getVertex(0),blTris[i]->getVertex(1),bedges, removed);
- addOrRemove(blTris[i]->getVertex(1),blTris[i]->getVertex(2),bedges, removed);
- addOrRemove(blTris[i]->getVertex(2),blTris[i]->getVertex(0),bedges, removed);
+ if(blQuads[i]->getVertex(j)->onWhat() == gf)
+ verts.insert(blQuads[i]->getVertex(j));
+ }
+ for(unsigned int i = 0; i < blTris.size(); i++) {
+ addOrRemove(blTris[i]->getVertex(0), blTris[i]->getVertex(1), bedges,
+ removed);
+ addOrRemove(blTris[i]->getVertex(1), blTris[i]->getVertex(2), bedges,
+ removed);
+ addOrRemove(blTris[i]->getVertex(2), blTris[i]->getVertex(0), bedges,
+ removed);
for(int j = 0; j < 3; j++)
- if(blTris[i]->getVertex(j)->onWhat() == gf)verts.insert(blTris[i]->getVertex(j));
+ if(blTris[i]->getVertex(j)->onWhat() == gf)
+ verts.insert(blTris[i]->getVertex(j));
}
- discreteEdge ne(gf->model(), 444444,0,
- (*edges.begin())->getEndVertex());
- std::vector<GEdge*> hop;
- std::set<MEdge,Less_Edge>::iterator it = bedges.begin();
+ discreteEdge ne(gf->model(), 444444, 0, (*edges.begin())->getEndVertex());
+ std::vector<GEdge *> hop;
+ std::set<MEdge, Less_Edge>::iterator it = bedges.begin();
FILE *ff = 0;
- if(debug) ff = Fopen("toto.pos","w");
+ if(debug) ff = Fopen("toto.pos", "w");
if(ff) fprintf(ff, "View \" \"{\n");
- for(; it != bedges.end(); ++it){
- ne.lines.push_back(new MLine(it->getVertex(0),it->getVertex(1)));
+ for(; it != bedges.end(); ++it) {
+ ne.lines.push_back(new MLine(it->getVertex(0), it->getVertex(1)));
if(ff)
- fprintf(ff, "SL(%g,%g,%g,%g,%g,%g){1,1};\n",
- it->getVertex(0)->x(),it->getVertex(0)->y(),it->getVertex(0)->z(),
- it->getVertex(1)->x(),it->getVertex(1)->y(),it->getVertex(1)->z());
+ fprintf(ff, "SL(%g,%g,%g,%g,%g,%g){1,1};\n", it->getVertex(0)->x(),
+ it->getVertex(0)->y(), it->getVertex(0)->z(),
+ it->getVertex(1)->x(), it->getVertex(1)->y(),
+ it->getVertex(1)->z());
}
- if(ff){
+ if(ff) {
fprintf(ff, "};\n");
fclose(ff);
}
@@ -920,48 +964,49 @@ static void modifyInitialMeshForBoundaryLayers(GFace *gf,
deMeshGFace kil_;
kil_(gf);
- meshGenerator(gf, 0, 0, true , false, &hop);
+ meshGenerator(gf, 0, 0, true, false, &hop);
}
-static bool improved_translate(GFace* gf,MVertex* vertex,SVector3& v1,SVector3& v2)
+static bool improved_translate(GFace *gf, MVertex *vertex, SVector3 &v1,
+ SVector3 &v2)
{
- double x,y;
+ double x, y;
double angle;
SPoint2 point;
- SVector3 s1,s2;
+ SVector3 s1, s2;
SVector3 normal;
- SVector3 basis_u,basis_v;
- Pair<SVector3,SVector3> derivatives;
+ SVector3 basis_u, basis_v;
+ Pair<SVector3, SVector3> derivatives;
reparamMeshVertexOnFace(vertex, gf, point);
x = point.x();
y = point.y();
- angle = backgroundMesh::current()->getAngle(x,y,0.0);
+ angle = backgroundMesh::current()->getAngle(x, y, 0.0);
derivatives = gf->firstDer(point);
s1 = derivatives.first();
s2 = derivatives.second();
- normal = crossprod(s1,s2);
+ normal = crossprod(s1, s2);
basis_u = s1;
basis_u.normalize();
- basis_v = crossprod(normal,basis_u);
+ basis_v = crossprod(normal, basis_u);
basis_v.normalize();
- v1 = basis_u*cos(angle) + basis_v*sin(angle);
- v2 = crossprod(v1,normal);
+ v1 = basis_u * cos(angle) + basis_v * sin(angle);
+ v2 = crossprod(v1, normal);
v2.normalize();
return 1;
}
-static void directions_storage(GFace* gf)
+static void directions_storage(GFace *gf)
{
- std::set<MVertex*> vertices;
- for(unsigned int i = 0; i < gf->getNumMeshElements(); i++){
- MElement* element = gf->getMeshElement(i);
- for(std::size_t j = 0; j < element->getNumVertices(); j++){
+ std::set<MVertex *> vertices;
+ for(unsigned int i = 0; i < gf->getNumMeshElements(); i++) {
+ MElement *element = gf->getMeshElement(i);
+ for(std::size_t j = 0; j < element->getNumVertices(); j++) {
MVertex *vertex = element->getVertex(j);
vertices.insert(vertex);
}
@@ -974,18 +1019,19 @@ static void directions_storage(GFace* gf)
gf->storage3.clear();
gf->storage4.clear();
- for(std::set<MVertex*>::iterator it = vertices.begin();
- it != vertices.end(); it++){
+ for(std::set<MVertex *>::iterator it = vertices.begin(); it != vertices.end();
+ it++) {
SPoint2 point;
SVector3 v1;
SVector3 v2;
- bool ok = improved_translate(gf,*it,v1,v2);
- if(ok){
- gf->storage1.push_back(SPoint3((*it)->x(),(*it)->y(),(*it)->z()));
+ bool ok = improved_translate(gf, *it, v1, v2);
+ if(ok) {
+ gf->storage1.push_back(SPoint3((*it)->x(), (*it)->y(), (*it)->z()));
gf->storage2.push_back(v1);
gf->storage3.push_back(v2);
- reparamMeshVertexOnFace(*it,gf,point);
- gf->storage4.push_back(backgroundMesh::current()->operator()(point.x(),point.y(),0.0));
+ reparamMeshVertexOnFace(*it, gf, point);
+ gf->storage4.push_back(
+ backgroundMesh::current()->operator()(point.x(), point.y(), 0.0));
}
}
@@ -994,42 +1040,42 @@ static void directions_storage(GFace* gf)
// Builds An initial triangular mesh that respects the boundaries of
// the domain, including embedded points and surfaces
-bool meshGenerator(GFace *gf, int RECUR_ITER,
- bool repairSelfIntersecting1dMesh,
- bool onlyInitialMesh,
- bool debug,
- std::vector<GEdge*> *replacement_edges)
+bool meshGenerator(GFace *gf, int RECUR_ITER, bool repairSelfIntersecting1dMesh,
+ bool onlyInitialMesh, bool debug,
+ std::vector<GEdge *> *replacement_edges)
{
// onlyInitialMesh=true;
BDS_GeomEntity CLASS_F(1, 2);
BDS_GeomEntity CLASS_EXTERIOR(1, 3);
- std::map<BDS_Point*, MVertex*,PointLessThan> recoverMap;
- std::map<MVertex*, BDS_Point*> recoverMapInv;
- std::vector<GEdge*> edges = replacement_edges ? *replacement_edges : gf->edges();
+ std::map<BDS_Point *, MVertex *, PointLessThan> recoverMap;
+ std::map<MVertex *, BDS_Point *> recoverMapInv;
+ std::vector<GEdge *> edges =
+ replacement_edges ? *replacement_edges : gf->edges();
// build a set with all points of the boundaries
- std::set<MVertex*, MVertexLessThanNum> all_vertices, boundary;
- std::vector<GEdge*>::const_iterator ite = edges.begin();
+ std::set<MVertex *, MVertexLessThanNum> all_vertices, boundary;
+ std::vector<GEdge *>::const_iterator ite = edges.begin();
FILE *fdeb = NULL;
- if(debug && RECUR_ITER == 0){
+ if(debug && RECUR_ITER == 0) {
char name[245];
sprintf(name, "surface%d-boundary-real.pos", gf->tag());
- fdeb = fopen (name,"w");
- fprintf(fdeb,"View \"\"{\n");
+ fdeb = fopen(name, "w");
+ fprintf(fdeb, "View \"\"{\n");
}
- while(ite != edges.end()){
+ while(ite != edges.end()) {
if((*ite)->isSeam(gf)) return false;
- if(!(*ite)->isMeshDegenerated()){
- for(unsigned int i = 0; i< (*ite)->lines.size(); i++){
+ if(!(*ite)->isMeshDegenerated()) {
+ for(unsigned int i = 0; i < (*ite)->lines.size(); i++) {
MVertex *v1 = (*ite)->lines[i]->getVertex(0);
MVertex *v2 = (*ite)->lines[i]->getVertex(1);
- if (fdeb){
- fprintf(fdeb,"SL(%g,%g,%g,%g,%g,%g){%d,%d};\n",
- v1->x(),v1->y(),v1->z(),v2->x(),v2->y(),v2->z(),(*ite)->tag(),(*ite)->tag());
- }
+ if(fdeb) {
+ fprintf(fdeb, "SL(%g,%g,%g,%g,%g,%g){%d,%d};\n", v1->x(), v1->y(),
+ v1->z(), v2->x(), v2->y(), v2->z(), (*ite)->tag(),
+ (*ite)->tag());
+ }
all_vertices.insert(v1);
all_vertices.insert(v2);
@@ -1048,23 +1094,25 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
++ite;
}
- if (fdeb){
- fprintf(fdeb,"};\n");
+ if(fdeb) {
+ fprintf(fdeb, "};\n");
fclose(fdeb);
}
- if(boundary.size()){
- Msg::Error("The 1D mesh seems not to be forming a closed loop (%d boundary points are considered once)",boundary.size());
+ if(boundary.size()) {
+ Msg::Error("The 1D mesh seems not to be forming a closed loop (%d boundary "
+ "points are considered once)",
+ boundary.size());
gf->meshStatistics.status = GFace::FAILED;
return false;
}
- std::vector<GEdge*> const& emb_edges = gf->embeddedEdges();
+ std::vector<GEdge *> const &emb_edges = gf->embeddedEdges();
ite = emb_edges.begin();
- while(ite != emb_edges.end()){
- if(!(*ite)->isMeshDegenerated()){
+ while(ite != emb_edges.end()) {
+ if(!(*ite)->isMeshDegenerated()) {
all_vertices.insert((*ite)->mesh_vertices.begin(),
- (*ite)->mesh_vertices.end() );
+ (*ite)->mesh_vertices.end());
all_vertices.insert((*ite)->getBeginVertex()->mesh_vertices.begin(),
(*ite)->getBeginVertex()->mesh_vertices.end());
all_vertices.insert((*ite)->getEndVertex()->mesh_vertices.begin(),
@@ -1074,9 +1122,9 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
}
// add embedded vertices
- std::list<GVertex*> emb_vertx = gf->embeddedVertices();
- std::list<GVertex*>::iterator itvx = emb_vertx.begin();
- while(itvx != emb_vertx.end()){
+ std::list<GVertex *> emb_vertx = gf->embeddedVertices();
+ std::list<GVertex *>::iterator itvx = emb_vertx.begin();
+ while(itvx != emb_vertx.end()) {
all_vertices.insert((*itvx)->mesh_vertices.begin(),
(*itvx)->mesh_vertices.end());
++itvx;
@@ -1086,19 +1134,19 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
all_vertices.insert(gf->additionalVertices.begin(),
gf->additionalVertices.end());
-
- if(all_vertices.size() < 3){
+ if(all_vertices.size() < 3) {
Msg::Warning("Mesh Generation of Model Face %d Skipped: "
"Only %d mesh vertices on the contours",
gf->tag(), all_vertices.size());
gf->meshStatistics.status = GFace::DONE;
return true;
}
- if(all_vertices.size() == 3){
+ if(all_vertices.size() == 3) {
MVertex *vv[3];
int i = 0;
- for(std::set<MVertex*, MVertexLessThanNum>::iterator it = all_vertices.begin();
- it != all_vertices.end(); it++){
+ for(std::set<MVertex *, MVertexLessThanNum>::iterator it =
+ all_vertices.begin();
+ it != all_vertices.end(); it++) {
vv[i++] = *it;
}
gf->triangles.push_back(new MTriangle(vv[0], vv[1], vv[2]));
@@ -1110,11 +1158,12 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
// meshing procedure
BDS_Mesh *m = new BDS_Mesh;
- std::vector<BDS_Point*> points(all_vertices.size());
+ std::vector<BDS_Point *> points(all_vertices.size());
SBoundingBox3d bbox;
int count = 0;
- for(std::set<MVertex*, MVertexLessThanNum>::iterator it = all_vertices.begin();
- it != all_vertices.end(); it++){
+ for(std::set<MVertex *, MVertexLessThanNum>::iterator it =
+ all_vertices.begin();
+ it != all_vertices.end(); it++) {
MVertex *here = *it;
GEntity *ge = here->onWhat();
SPoint2 param;
@@ -1141,17 +1190,16 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
SVector3 dd(bbox.max(), bbox.min());
double LC2D = norm(dd);
DocRecord doc(points.size() + 4);
- for(unsigned int i = 0; i < points.size(); i++){
+ for(unsigned int i = 0; i < points.size(); i++) {
double XX = CTX::instance()->mesh.randFactor * LC2D * (double)rand() /
- (double)RAND_MAX;
+ (double)RAND_MAX;
double YY = CTX::instance()->mesh.randFactor * LC2D * (double)rand() /
- (double)RAND_MAX;
+ (double)RAND_MAX;
// printf("%22.15E %22.15E \n",XX,YY);
doc.points[i].where.h = points[i]->u + XX;
doc.points[i].where.v = points[i]->v + YY;
doc.points[i].data = points[i];
doc.points[i].adjacent = NULL;
-
}
// increase the size of the bounding box
@@ -1163,7 +1211,7 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
{bbox.min().x(), bbox.max().y()},
{bbox.max().x(), bbox.min().y()},
{bbox.max().x(), bbox.max().y()}};
- for(int ip = 0; ip < 4; ip++){
+ for(int ip = 0; ip < 4; ip++) {
BDS_Point *pp = m->add_point(-ip - 1, bb[ip][0], bb[ip][1], gf);
m->add_geom(gf->tag(), 2);
BDS_GeomEntity *g = m->get_geom(gf->tag(), 2);
@@ -1180,55 +1228,53 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
// -) It contains triangles outside the domain (the first edge
// loop is the outer one)
Msg::Debug("Meshing of the convex hull (%d points)", points.size());
- try{
+ try {
doc.MakeMeshWithPoints();
- }
- catch(const char *err){
+ } catch(const char *err) {
Msg::Error("%s", err);
}
Msg::Debug("Meshing of the convex hull (%d points) done", points.size());
-
for(int i = 0; i < doc.numTriangles; i++) {
int a = doc.triangles[i].a;
int b = doc.triangles[i].b;
int c = doc.triangles[i].c;
int n = doc.numPoints;
- if(a < 0 || a >= n || b < 0 || b >= n || c < 0 || c >= n){
+ if(a < 0 || a >= n || b < 0 || b >= n || c < 0 || c >= n) {
Msg::Warning("Skipping bad triangle %d", i);
continue;
}
- BDS_Point *p1 = (BDS_Point*)doc.points[doc.triangles[i].a].data;
- BDS_Point *p2 = (BDS_Point*)doc.points[doc.triangles[i].b].data;
- BDS_Point *p3 = (BDS_Point*)doc.points[doc.triangles[i].c].data;
+ BDS_Point *p1 = (BDS_Point *)doc.points[doc.triangles[i].a].data;
+ BDS_Point *p2 = (BDS_Point *)doc.points[doc.triangles[i].b].data;
+ BDS_Point *p3 = (BDS_Point *)doc.points[doc.triangles[i].c].data;
m->add_triangle(p1->iD, p2->iD, p3->iD);
}
}
#else
{
- std::vector<MVertex*> v;
- std::vector<MTriangle*> result;
+ std::vector<MVertex *> v;
+ std::vector<MTriangle *> result;
v.insert(v.end(), all_vertices.begin(), all_vertices.end());
- std::map<MVertex*,SPoint3> pos;
+ std::map<MVertex *, SPoint3> pos;
for(unsigned int i = 0; i < v.size(); i++) {
MVertex *v0 = v[i];
- BDS_Point *p0 = recoverMapInv[v0];
- pos[v0] = SPoint3(v0->x(),v0->y(),v0->z());
- v0->setXYZ(p0->u,p0->v,0.0);
+ BDS_Point *p0 = recoverMapInv[v0];
+ pos[v0] = SPoint3(v0->x(), v0->y(), v0->z());
+ v0->setXYZ(p0->u, p0->v, 0.0);
}
delaunayMeshIn2D(v, result, 0);
- for(unsigned int i = 0; i < v.size()-4; i++) {
+ for(unsigned int i = 0; i < v.size() - 4; i++) {
MVertex *v0 = v[i];
SPoint3 pp = pos[v0];
- v0->setXYZ(pp.x(),pp.y(),pp.z());
+ v0->setXYZ(pp.x(), pp.y(), pp.z());
}
// add the four corners
- for(int ip = 0; ip < 4; ip++){
- MVertex *vv = v[v.size()-ip-1];
- BDS_Point *pp = m->add_point(-ip - 1, vv->x(),vv->y(), gf);
+ for(int ip = 0; ip < 4; ip++) {
+ MVertex *vv = v[v.size() - ip - 1];
+ BDS_Point *pp = m->add_point(-ip - 1, vv->x(), vv->y(), gf);
m->add_geom(gf->tag(), 2);
recoverMapInv[vv] = pp;
BDS_GeomEntity *g = m->get_geom(gf->tag(), 2);
@@ -1239,20 +1285,20 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
MVertex *v0 = result[i]->getVertex(0);
MVertex *v1 = result[i]->getVertex(1);
MVertex *v2 = result[i]->getVertex(2);
- BDS_Point *p0 = recoverMapInv[v0];
- BDS_Point *p1 = recoverMapInv[v1];
- BDS_Point *p2 = recoverMapInv[v2];
+ BDS_Point *p0 = recoverMapInv[v0];
+ BDS_Point *p1 = recoverMapInv[v1];
+ BDS_Point *p2 = recoverMapInv[v2];
m->add_triangle(p0->iD, p1->iD, p2->iD);
}
}
#endif
- if(debug && RECUR_ITER == 0){
+ if(debug && RECUR_ITER == 0) {
char name[245];
sprintf(name, "surface%d-initial-real.pos", gf->tag());
- outputScalarField(m->triangles, name, 0,gf);
+ outputScalarField(m->triangles, name, 0, gf);
sprintf(name, "surface%d-initial-param.pos", gf->tag());
- outputScalarField(m->triangles, name, 1,gf);
+ outputScalarField(m->triangles, name, 1, gf);
}
// Recover the boundary edges and compute characteristic lenghts using mesh
@@ -1262,24 +1308,26 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
std::set<EdgeToRecover> edgesToRecover;
std::set<EdgeToRecover> edgesNotRecovered;
ite = edges.begin();
- while(ite != edges.end()){
+ while(ite != edges.end()) {
if(!(*ite)->isMeshDegenerated())
- recoverEdge(m, *ite, recoverMapInv, &edgesToRecover, &edgesNotRecovered, 1);
+ recoverEdge(m, *ite, recoverMapInv, &edgesToRecover, &edgesNotRecovered,
+ 1);
++ite;
}
ite = emb_edges.begin();
- while(ite != emb_edges.end()){
+ while(ite != emb_edges.end()) {
if(!(*ite)->isMeshDegenerated())
- recoverEdge(m, *ite, recoverMapInv, &edgesToRecover, &edgesNotRecovered, 1);
+ recoverEdge(m, *ite, recoverMapInv, &edgesToRecover, &edgesNotRecovered,
+ 1);
++ite;
}
-
// effectively recover the medge
ite = edges.begin();
- while(ite != edges.end()){
- if(!(*ite)->isMeshDegenerated()){
- if(!recoverEdge(m, *ite, recoverMapInv, &edgesToRecover, &edgesNotRecovered, 2)){
+ while(ite != edges.end()) {
+ if(!(*ite)->isMeshDegenerated()) {
+ if(!recoverEdge(m, *ite, recoverMapInv, &edgesToRecover,
+ &edgesNotRecovered, 2)) {
delete m;
gf->meshStatistics.status = GFace::FAILED;
return false;
@@ -1288,10 +1336,10 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
++ite;
}
- Msg::Debug("Recovering %d mesh Edges (%d not recovered)", edgesToRecover.size(),
- edgesNotRecovered.size());
+ Msg::Debug("Recovering %d mesh Edges (%d not recovered)",
+ edgesToRecover.size(), edgesNotRecovered.size());
- if(edgesNotRecovered.size()){
+ if(edgesNotRecovered.size()) {
std::ostringstream sstream;
for(std::set<EdgeToRecover>::iterator itr = edgesNotRecovered.begin();
itr != edgesNotRecovered.end(); ++itr)
@@ -1301,7 +1349,7 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
if(repairSelfIntersecting1dMesh)
Msg::Warning("8-| Gmsh splits those edges and tries again");
- if(debug){
+ if(debug) {
char name[245];
sprintf(name, "surface%d-not_yet_recovered-real-%d.msh", gf->tag(),
RECUR_ITER);
@@ -1311,11 +1359,11 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
std::list<GFace *> facesToRemesh;
if(repairSelfIntersecting1dMesh)
remeshUnrecoveredEdges(recoverMapInv, edgesNotRecovered, facesToRemesh);
- else{
+ else {
std::set<EdgeToRecover>::iterator itr = edgesNotRecovered.begin();
- //int *_error = new int[3 * edgesNotRecovered.size()];
+ // int *_error = new int[3 * edgesNotRecovered.size()];
int I = 0;
- for(; itr != edgesNotRecovered.end(); ++itr){
+ for(; itr != edgesNotRecovered.end(); ++itr) {
int p1 = itr->p1;
int p2 = itr->p2;
int tag = itr->ge->tag();
@@ -1325,15 +1373,14 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
//_error[3 * I + 2] = tag;
I++;
}
- //throw _error;
+ // throw _error;
}
// delete the mesh
delete m;
if(RECUR_ITER < 10 && facesToRemesh.size() == 0)
- return meshGenerator
- (gf, RECUR_ITER + 1, repairSelfIntersecting1dMesh, onlyInitialMesh,
- debug, replacement_edges);
+ return meshGenerator(gf, RECUR_ITER + 1, repairSelfIntersecting1dMesh,
+ onlyInitialMesh, debug, replacement_edges);
return false;
}
@@ -1346,18 +1393,17 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
// look for a triangle that has a negative node and recursively tag all
// exterior triangles
{
- std::vector<BDS_Face*>::iterator itt = m->triangles.begin();
- while(itt != m->triangles.end()){
+ std::vector<BDS_Face *>::iterator itt = m->triangles.begin();
+ while(itt != m->triangles.end()) {
(*itt)->g = 0;
++itt;
}
itt = m->triangles.begin();
- while(itt != m->triangles.end()){
+ while(itt != m->triangles.end()) {
BDS_Face *t = *itt;
BDS_Point *n[4];
t->getNodes(n);
- if(n[0]->iD < 0 || n[1]->iD < 0 ||
- n[2]->iD < 0 ) {
+ if(n[0]->iD < 0 || n[1]->iD < 0 || n[2]->iD < 0) {
recur_tag(t, &CLASS_EXTERIOR);
break;
}
@@ -1367,40 +1413,40 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
// now find an edge that has belongs to one of the exterior triangles
{
- std::vector<BDS_Edge*>::iterator ite = m->edges.begin();
- while(ite != m->edges.end()){
+ std::vector<BDS_Edge *>::iterator ite = m->edges.begin();
+ while(ite != m->edges.end()) {
BDS_Edge *e = *ite;
- if(e->g && e->numfaces() == 2){
- if(e->faces(0)->g == &CLASS_EXTERIOR){
+ if(e->g && e->numfaces() == 2) {
+ if(e->faces(0)->g == &CLASS_EXTERIOR) {
recur_tag(e->faces(1), &CLASS_F);
break;
}
- else if(e->faces(1)->g == &CLASS_EXTERIOR){
+ else if(e->faces(1)->g == &CLASS_EXTERIOR) {
recur_tag(e->faces(0), &CLASS_F);
break;
}
}
++ite;
}
- std::vector<BDS_Face*>::iterator itt = m->triangles.begin();
- while(itt != m->triangles.end()){
+ std::vector<BDS_Face *>::iterator itt = m->triangles.begin();
+ while(itt != m->triangles.end()) {
if((*itt)->g == &CLASS_EXTERIOR) (*itt)->g = 0;
++itt;
}
}
{
- std::vector<BDS_Edge*>::iterator ite = m->edges.begin();
- while(ite != m->edges.end()){
+ std::vector<BDS_Edge *>::iterator ite = m->edges.begin();
+ while(ite != m->edges.end()) {
BDS_Edge *e = *ite;
- if(e->g && e->numfaces() == 2){
+ if(e->g && e->numfaces() == 2) {
BDS_Point *oface[2];
e->oppositeof(oface);
- if(oface[0]->iD < 0){
+ if(oface[0]->iD < 0) {
recur_tag(e->faces(1), &CLASS_F);
break;
}
- else if(oface[1]->iD < 0){
+ else if(oface[1]->iD < 0) {
recur_tag(e->faces(0), &CLASS_F);
break;
}
@@ -1410,43 +1456,45 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
}
ite = emb_edges.begin();
- while(ite != emb_edges.end()){
+ while(ite != emb_edges.end()) {
if(!(*ite)->isMeshDegenerated())
- recoverEdge(m, *ite, recoverMapInv, &edgesToRecover, &edgesNotRecovered, 2);
+ recoverEdge(m, *ite, recoverMapInv, &edgesToRecover, &edgesNotRecovered,
+ 2);
++ite;
}
// compute characteristic lengths at vertices
- if(CTX::instance()->mesh.algo2d != ALGO_2D_BAMG && !onlyInitialMesh){
- Msg::Debug("Computing mesh size field at mesh vertices %d",
- edgesToRecover.size());
- std::set<BDS_Point*, PointLessThan>::iterator it = m->points.begin();
- for(; it != m->points.end();++it){
- // for(int i = 0; i < doc.numPoints; i++){
- // BDS_Point *pp = (BDS_Point*)doc.points[i].data;
- BDS_Point *pp = *it;
- std::map<BDS_Point*, MVertex*,PointLessThan>::iterator itv = recoverMap.find(pp);
- if(itv != recoverMap.end()){
- MVertex *here = itv->second;
- GEntity *ge = here->onWhat();
- if(ge->dim() == 0){
- pp->lcBGM() = BGM_MeshSize(ge, 0, 0, here->x(), here->y(), here->z());
- }
- else if(ge->dim() == 1){
- double u;
- here->getParameter(0, u);
- pp->lcBGM() = BGM_MeshSize(ge, u, 0, here->x(), here->y(), here->z());
- }
- else
- pp->lcBGM() = MAX_LC;
- pp->lc() = pp->lcBGM();
+ if(CTX::instance()->mesh.algo2d != ALGO_2D_BAMG && !onlyInitialMesh) {
+ Msg::Debug("Computing mesh size field at mesh vertices %d",
+ edgesToRecover.size());
+ std::set<BDS_Point *, PointLessThan>::iterator it = m->points.begin();
+ for(; it != m->points.end(); ++it) {
+ // for(int i = 0; i < doc.numPoints; i++){
+ // BDS_Point *pp = (BDS_Point*)doc.points[i].data;
+ BDS_Point *pp = *it;
+ std::map<BDS_Point *, MVertex *, PointLessThan>::iterator itv =
+ recoverMap.find(pp);
+ if(itv != recoverMap.end()) {
+ MVertex *here = itv->second;
+ GEntity *ge = here->onWhat();
+ if(ge->dim() == 0) {
+ pp->lcBGM() = BGM_MeshSize(ge, 0, 0, here->x(), here->y(), here->z());
+ }
+ else if(ge->dim() == 1) {
+ double u;
+ here->getParameter(0, u);
+ pp->lcBGM() = BGM_MeshSize(ge, u, 0, here->x(), here->y(), here->z());
}
+ else
+ pp->lcBGM() = MAX_LC;
+ pp->lc() = pp->lcBGM();
}
+ }
}
// delete useless stuff
- std::vector<BDS_Face*>::iterator itt = m->triangles.begin();
- while(itt != m->triangles.end()){
+ std::vector<BDS_Face *>::iterator itt = m->triangles.begin();
+ while(itt != m->triangles.end()) {
BDS_Face *t = *itt;
if(!t->g) m->del_face(t);
++itt;
@@ -1454,14 +1502,13 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
m->cleanup();
{
- std::vector<BDS_Edge*>::iterator ite = m->edges.begin();
- while(ite != m->edges.end()){
+ std::vector<BDS_Edge *>::iterator ite = m->edges.begin();
+ while(ite != m->edges.end()) {
BDS_Edge *e = *ite;
if(e->numfaces() == 0)
m->del_edge(e);
- else{
- if(!e->g)
- e->g = &CLASS_F;
+ else {
+ if(!e->g) e->g = &CLASS_F;
if(!e->p1->g || e->p1->g->classif_degree > e->g->classif_degree)
e->p1->g = e->g;
if(!e->p2->g || e->p2->g->classif_degree > e->g->classif_degree)
@@ -1476,29 +1523,29 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
m->del_point(m->find_point(-3));
m->del_point(m->find_point(-4));
- if(debug){
+ if(debug) {
char name[245];
sprintf(name, "surface%d-recovered-real.pos", gf->tag());
- outputScalarField(m->triangles, name, 0,gf);
+ outputScalarField(m->triangles, name, 0, gf);
sprintf(name, "surface%d-recovered-param.pos", gf->tag());
- outputScalarField(m->triangles, name, 1,gf);
+ outputScalarField(m->triangles, name, 1, gf);
}
- if(1){
- std::vector<BDS_Face*>::iterator itt = m->triangles.begin();
- while(itt != m->triangles.end()){
+ if(1) {
+ std::vector<BDS_Face *>::iterator itt = m->triangles.begin();
+ while(itt != m->triangles.end()) {
BDS_Face *t = *itt;
- if(!t->deleted){
+ if(!t->deleted) {
BDS_Point *n[4];
t->getNodes(n);
MVertex *v1 = recoverMap[n[0]];
MVertex *v2 = recoverMap[n[1]];
MVertex *v3 = recoverMap[n[2]];
- if(!n[3]){
+ if(!n[3]) {
if(v1 != v2 && v1 != v3 && v2 != v3)
gf->triangles.push_back(new MTriangle(v1, v2, v3));
}
- else{
+ else {
MVertex *v4 = recoverMap[n[3]];
gf->quadrangles.push_back(new MQuadrangle(v1, v2, v3, v4));
}
@@ -1512,26 +1559,28 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
Msg::Debug("Delaunizing the initial mesh");
delaunayizeBDS(gf, *m, nb_swap);
}
- //gf->triangles.clear();
- //gf->quadrangles.clear();
+ // gf->triangles.clear();
+ // gf->quadrangles.clear();
// only delete the mesh data stored in the base GFace class
gf->GFace::deleteMesh();
Msg::Debug("Starting to add internal points");
// start mesh generation
- if(!algoDelaunay2D(gf) && !onlyInitialMesh){
- // if(CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine || 1) {
+ if(!algoDelaunay2D(gf) && !onlyInitialMesh) {
+ // if(CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine ||
+ // 1) {
// backgroundMesh::unset();
// buildBackGroundMesh(gf);
// }
refineMeshBDS(gf, *m, CTX::instance()->mesh.refineSteps, true,
- &recoverMapInv,NULL);
+ &recoverMapInv, NULL);
optimizeMeshBDS(gf, *m, 2);
refineMeshBDS(gf, *m, CTX::instance()->mesh.refineSteps, false,
- &recoverMapInv,NULL);
+ &recoverMapInv, NULL);
optimizeMeshBDS(gf, *m, 2);
- // if(CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine || 1) {
+ // if(CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine ||
+ // 1) {
// backgroundMesh::unset();
// }
}
@@ -1542,7 +1591,8 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
gf->meshStatistics.smallest_edge_length,
gf->meshStatistics.nbEdge,
gf->meshStatistics.nbGoodLength);
- printf("=== Efficiency index is tau=%g\n", gf->meshStatistics.efficiency_index);
+ printf("=== Efficiency index is tau=%g\n",
+ gf->meshStatistics.efficiency_index);
*/
gf->meshStatistics.status = GFace::DONE;
@@ -1550,9 +1600,9 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
// fill the small gmsh structures
BDS2GMSH(m, gf, recoverMap);
- std::vector<MQuadrangle*> blQuads;
- std::vector<MTriangle*> blTris;
- std::set<MVertex*> verts;
+ std::vector<MQuadrangle *> blQuads;
+ std::vector<MTriangle *> blTris;
+ std::set<MVertex *> verts;
bool infty = false;
if(gf->getMeshingAlgo() == ALGO_2D_FRONTAL_QUAD ||
@@ -1568,32 +1618,33 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
// the delaunay algo is based directly on internal gmsh structures BDS mesh is
// passed in order not to recompute local coordinates of vertices
- if(algoDelaunay2D(gf) && !onlyInitialMesh){
- if(gf->getMeshingAlgo() == ALGO_2D_FRONTAL){
+ if(algoDelaunay2D(gf) && !onlyInitialMesh) {
+ if(gf->getMeshingAlgo() == ALGO_2D_FRONTAL) {
bowyerWatsonFrontal(gf);
}
- else if(gf->getMeshingAlgo() == ALGO_2D_FRONTAL_QUAD){
- bowyerWatsonFrontalLayers(gf,true);
+ else if(gf->getMeshingAlgo() == ALGO_2D_FRONTAL_QUAD) {
+ bowyerWatsonFrontalLayers(gf, true);
}
- else if(gf->getMeshingAlgo() == ALGO_2D_PACK_PRLGRMS){
+ else if(gf->getMeshingAlgo() == ALGO_2D_PACK_PRLGRMS) {
bowyerWatsonParallelograms(gf);
}
- else if(gf->getMeshingAlgo() == ALGO_2D_PACK_PRLGRMS_CSTR){
- bowyerWatsonParallelogramsConstrained(gf,gf->constr_vertices);
+ else if(gf->getMeshingAlgo() == ALGO_2D_PACK_PRLGRMS_CSTR) {
+ bowyerWatsonParallelogramsConstrained(gf, gf->constr_vertices);
}
else if(gf->getMeshingAlgo() == ALGO_2D_DELAUNAY ||
- gf->getMeshingAlgo() == ALGO_2D_AUTO){
+ gf->getMeshingAlgo() == ALGO_2D_AUTO) {
bowyerWatson(gf);
}
else {
- bowyerWatson(gf,15000);
+ bowyerWatson(gf, 15000);
meshGFaceBamg(gf);
}
- if(!infty || !(CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine))
+ if(!infty ||
+ !(CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine))
laplaceSmoothing(gf, CTX::instance()->mesh.nbSmoothing, infty);
}
- if(debug){
+ if(debug) {
char name[256];
// sprintf(name, "trueBoundary%d.pos", gf->tag());
// std::vector<SPoint2> bnd;
@@ -1601,14 +1652,16 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
sprintf(name, "real%d.pos", gf->tag());
outputScalarField(m->triangles, name, 0, gf);
sprintf(name, "param%d.pos", gf->tag());
- outputScalarField(m->triangles, name,1,gf);
+ outputScalarField(m->triangles, name, 1, gf);
}
delete m;
- gf->quadrangles.insert(gf->quadrangles.begin(),blQuads.begin(),blQuads.end());
- gf->triangles.insert(gf->triangles.begin(),blTris.begin(),blTris.end());
- gf->mesh_vertices.insert(gf->mesh_vertices.begin(),verts.begin(),verts.end());
+ gf->quadrangles.insert(gf->quadrangles.begin(), blQuads.begin(),
+ blQuads.end());
+ gf->triangles.insert(gf->triangles.begin(), blTris.begin(), blTris.end());
+ gf->mesh_vertices.insert(gf->mesh_vertices.begin(), verts.begin(),
+ verts.end());
splitElementsInBoundaryLayerIfNeeded(gf);
@@ -1623,11 +1676,11 @@ bool meshGenerator(GFace *gf, int RECUR_ITER,
gf->meshStatistics.nbGoodQuality);
gf->mesh_vertices.insert(gf->mesh_vertices.end(),
- gf->additionalVertices.begin(),
- gf->additionalVertices.end());
+ gf->additionalVertices.begin(),
+ gf->additionalVertices.end());
gf->additionalVertices.clear();
- if(CTX::instance()->mesh.algo3d==ALGO_3D_RTREE){
+ if(CTX::instance()->mesh.algo3d == ALGO_3D_RTREE) {
directions_storage(gf);
}
@@ -1656,24 +1709,22 @@ static void printMesh1d(int iEdge, int seam, std::vector<SPoint2> &m)
}
*/
-static bool buildConsecutiveListOfVertices(GFace *gf, GEdgeLoop &gel,
- std::vector<BDS_Point*> &result,
- SBoundingBox3d &bbox, BDS_Mesh *m,
- std::map<BDS_Point*, MVertex*,PointLessThan>
- &recoverMap,
- int &count, int countTot, double tol,
- bool seam_the_first = false)
+static bool buildConsecutiveListOfVertices(
+ GFace *gf, GEdgeLoop &gel, std::vector<BDS_Point *> &result,
+ SBoundingBox3d &bbox, BDS_Mesh *m,
+ std::map<BDS_Point *, MVertex *, PointLessThan> &recoverMap, int &count,
+ int countTot, double tol, bool seam_the_first = false)
{
// for each edge, we build a list of points that are the mapping of the edge
// points on the face for seams, we build the list for every side for closed
// loops, we build it on both senses
- std::map<GEntity*, std::vector<SPoint2> > meshes;
- std::map<GEntity*, std::vector<SPoint2> > meshes_seam;
+ std::map<GEntity *, std::vector<SPoint2> > meshes;
+ std::map<GEntity *, std::vector<SPoint2> > meshes_seam;
const int MYDEBUG = false;
- std::map<BDS_Point*, MVertex*,PointLessThan> recoverMapLocal;
+ std::map<BDS_Point *, MVertex *, PointLessThan> recoverMapLocal;
result.clear();
count = 0;
@@ -1684,21 +1735,21 @@ static bool buildConsecutiveListOfVertices(GFace *gf, GEdgeLoop &gel,
printf("face %d with %d edges case %d\n", gf->tag(),
(int)gf->edges().size(), seam_the_first);
- while(it != gel.end()){
- GEdgeSigned ges = *it ;
+ while(it != gel.end()) {
+ GEdgeSigned ges = *it;
std::vector<SPoint2> mesh1d;
std::vector<SPoint2> mesh1d_seam;
bool seam = ges.ge->isSeam(gf);
- //if(seam) printf("face %d has seam %d\n", gf->tag(), ges.ge->tag());
+ // if(seam) printf("face %d has seam %d\n", gf->tag(), ges.ge->tag());
Range<double> range = ges.ge->parBounds(0);
MVertex *here = ges.ge->getBeginVertex()->mesh_vertices[0];
mesh1d.push_back(ges.ge->reparamOnFace(gf, range.low(), 1));
if(seam) mesh1d_seam.push_back(ges.ge->reparamOnFace(gf, range.low(), -1));
- for(unsigned int i = 0; i < ges.ge->mesh_vertices.size(); i++){
+ for(unsigned int i = 0; i < ges.ge->mesh_vertices.size(); i++) {
double u;
here = ges.ge->mesh_vertices[i];
here->getParameter(0, u);
@@ -1708,9 +1759,10 @@ static bool buildConsecutiveListOfVertices(GFace *gf, GEdgeLoop &gel,
here = ges.ge->getEndVertex()->mesh_vertices[0];
mesh1d.push_back(ges.ge->reparamOnFace(gf, range.high(), 1));
if(seam) mesh1d_seam.push_back(ges.ge->reparamOnFace(gf, range.high(), -1));
- meshes.insert(std::pair<GEntity*,std::vector<SPoint2> >(ges.ge, mesh1d));
- if(seam) meshes_seam.insert(std::pair<GEntity*,std::vector<SPoint2> >
- (ges.ge, mesh1d_seam));
+ meshes.insert(std::pair<GEntity *, std::vector<SPoint2> >(ges.ge, mesh1d));
+ if(seam)
+ meshes_seam.insert(
+ std::pair<GEntity *, std::vector<SPoint2> >(ges.ge, mesh1d_seam));
// printMesh1d(ges.ge->tag(), seam, mesh1d);
// if(seam) printMesh1d (ges.ge->tag(), seam, mesh1d_seam);
it++;
@@ -1722,13 +1774,12 @@ static bool buildConsecutiveListOfVertices(GFace *gf, GEdgeLoop &gel,
SPoint2 last_coord(0, 0);
int counter = 0;
- while(unordered.size()){
- if(MYDEBUG)
- printf("unordered.size() = %d\n", (int)unordered.size());
+ while(unordered.size()) {
+ if(MYDEBUG) printf("unordered.size() = %d\n", (int)unordered.size());
std::list<GEdgeSigned>::iterator it = unordered.begin();
- std::vector<SPoint2> coords;
+ std::vector<SPoint2> coords;
- while(it != unordered.end()){
+ while(it != unordered.end()) {
std::vector<SPoint2> mesh1d;
std::vector<SPoint2> mesh1d_seam;
std::vector<SPoint2> mesh1d_reversed;
@@ -1736,18 +1787,23 @@ static bool buildConsecutiveListOfVertices(GFace *gf, GEdgeLoop &gel,
GEdge *ge = (*it).ge;
bool seam = ge->isSeam(gf);
mesh1d = meshes[ge];
- if(seam){ mesh1d_seam = meshes_seam[ge]; }
- mesh1d_reversed.insert(mesh1d_reversed.begin(), mesh1d.rbegin(), mesh1d.rend());
- if(seam) mesh1d_seam_reversed.insert(mesh1d_seam_reversed.begin(),
- mesh1d_seam.rbegin(),mesh1d_seam.rend());
- if(!counter){
+ if(seam) {
+ mesh1d_seam = meshes_seam[ge];
+ }
+ mesh1d_reversed.insert(mesh1d_reversed.begin(), mesh1d.rbegin(),
+ mesh1d.rend());
+ if(seam)
+ mesh1d_seam_reversed.insert(mesh1d_seam_reversed.begin(),
+ mesh1d_seam.rbegin(), mesh1d_seam.rend());
+ if(!counter) {
counter++;
- if(seam && seam_the_first){
+ if(seam && seam_the_first) {
coords = ((*it)._sign == 1) ? mesh1d_seam : mesh1d_seam_reversed;
found = (*it);
- Msg::Info("This test case would have failed in previous Gmsh versions ;-)");
+ Msg::Info(
+ "This test case would have failed in previous Gmsh versions ;-)");
}
- else{
+ else {
coords = ((*it)._sign == 1) ? mesh1d : mesh1d_reversed;
found = (*it);
}
@@ -1756,48 +1812,47 @@ static bool buildConsecutiveListOfVertices(GFace *gf, GEdgeLoop &gel,
unordered.erase(it);
break;
}
- else{
- if(MYDEBUG)
- printf("Followed by edge = %d\n", (*it).ge->tag());
+ else {
+ if(MYDEBUG) printf("Followed by edge = %d\n", (*it).ge->tag());
SPoint2 first_coord = mesh1d[0];
double d = -1, d_reversed = -1, d_seam = -1, d_seam_reversed = -1;
d = dist2(last_coord, first_coord);
if(MYDEBUG)
printf("%g %g dist = %12.5E\n", first_coord.x(), first_coord.y(), d);
- if(d < tol){
+ if(d < tol) {
coords.clear();
coords = mesh1d;
- found = GEdgeSigned(1,ge);
+ found = GEdgeSigned(1, ge);
unordered.erase(it);
goto Finalize;
}
SPoint2 first_coord_reversed = mesh1d_reversed[0];
d_reversed = dist2(last_coord, first_coord_reversed);
if(MYDEBUG)
- printf("%g %g dist_reversed = %12.5E\n",
- first_coord_reversed.x(), first_coord_reversed.y(), d_reversed);
- if(d_reversed < tol){
+ printf("%g %g dist_reversed = %12.5E\n", first_coord_reversed.x(),
+ first_coord_reversed.y(), d_reversed);
+ if(d_reversed < tol) {
coords.clear();
coords = mesh1d_reversed;
- found = (GEdgeSigned(-1,ge));
+ found = (GEdgeSigned(-1, ge));
unordered.erase(it);
goto Finalize;
}
- if(seam){
+ if(seam) {
SPoint2 first_coord_seam = mesh1d_seam[0];
SPoint2 first_coord_seam_reversed = mesh1d_seam_reversed[0];
- d_seam = dist2(last_coord,first_coord_seam);
+ d_seam = dist2(last_coord, first_coord_seam);
if(MYDEBUG) printf("dist_seam = %12.5E\n", d_seam);
- if(d_seam < tol){
+ if(d_seam < tol) {
coords.clear();
coords = mesh1d_seam;
- found = (GEdgeSigned(1,ge));
+ found = (GEdgeSigned(1, ge));
unordered.erase(it);
goto Finalize;
}
d_seam_reversed = dist2(last_coord, first_coord_seam_reversed);
if(MYDEBUG) printf("dist_seam_reversed = %12.5E\n", d_seam_reversed);
- if(d_seam_reversed < tol){
+ if(d_seam_reversed < tol) {
coords.clear();
coords = mesh1d_seam_reversed;
found = GEdgeSigned(-1, ge);
@@ -1811,27 +1866,28 @@ static bool buildConsecutiveListOfVertices(GFace *gf, GEdgeLoop &gel,
}
Finalize:
if(MYDEBUG) printf("Finalize, found %d points\n", (int)coords.size());
- if(coords.size() == 0){
+ if(coords.size() == 0) {
// It has not worked : either tolerance is wrong or the first seam edge
// has to be taken with the other parametric coordinates (because it is
// only present once in the closure of the domain).
- for(std::map<BDS_Point*, MVertex*,PointLessThan>::iterator it = recoverMapLocal.begin();
- it != recoverMapLocal.end(); ++it){
+ for(std::map<BDS_Point *, MVertex *, PointLessThan>::iterator it =
+ recoverMapLocal.begin();
+ it != recoverMapLocal.end(); ++it) {
m->del_point(it->first);
}
return false;
}
- std::vector<MVertex*> edgeLoop;
- if(found._sign == 1){
- if(found.ge->getBeginVertex()){
+ std::vector<MVertex *> edgeLoop;
+ if(found._sign == 1) {
+ if(found.ge->getBeginVertex()) {
edgeLoop.push_back(found.ge->getBeginVertex()->mesh_vertices[0]);
for(unsigned int i = 0; i < found.ge->mesh_vertices.size(); i++)
edgeLoop.push_back(found.ge->mesh_vertices[i]);
}
}
- else{
- if(found.ge->getEndVertex()){
+ else {
+ if(found.ge->getEndVertex()) {
edgeLoop.push_back(found.ge->getEndVertex()->mesh_vertices[0]);
for(int i = found.ge->mesh_vertices.size() - 1; i >= 0; i--)
edgeLoop.push_back(found.ge->mesh_vertices[i]);
@@ -1839,35 +1895,36 @@ static bool buildConsecutiveListOfVertices(GFace *gf, GEdgeLoop &gel,
}
if(MYDEBUG)
- printf("edge %d size %d size %d\n",
- found.ge->tag(), (int)edgeLoop.size(), (int)coords.size());
+ printf("edge %d size %d size %d\n", found.ge->tag(), (int)edgeLoop.size(),
+ (int)coords.size());
- std::vector<BDS_Point*> edgeLoop_BDS;
- for(unsigned int i = 0; i < edgeLoop.size(); i++){
+ std::vector<BDS_Point *> edgeLoop_BDS;
+ for(unsigned int i = 0; i < edgeLoop.size(); i++) {
MVertex *here = edgeLoop[i];
GEntity *ge = here->onWhat();
BDS_Point *pp = 0;
- if(ge->dim() == 0){
+ if(ge->dim() == 0) {
// Point might already be part of other loop
- for(std::map<BDS_Point*, MVertex*, PointLessThan>::iterator it = recoverMap.begin();
- it != recoverMap.end(); ++it){
- if(it->second == here){
+ for(std::map<BDS_Point *, MVertex *, PointLessThan>::iterator it =
+ recoverMap.begin();
+ it != recoverMap.end(); ++it) {
+ if(it->second == here) {
pp = it->first;
break;
}
}
}
- if(pp == 0){
+ if(pp == 0) {
double U, V;
SPoint2 param = coords[i];
U = param.x();
V = param.y();
pp = m->add_point(count + countTot, U, V, gf);
- if(ge->dim() == 0){
+ if(ge->dim() == 0) {
pp->lcBGM() = BGM_MeshSize(ge, 0, 0, here->x(), here->y(), here->z());
}
- else if(ge->dim() == 1){
+ else if(ge->dim() == 1) {
double u;
here->getParameter(0, u);
pp->lcBGM() = BGM_MeshSize(ge, u, 0, here->x(), here->y(), here->z());
@@ -1880,8 +1937,8 @@ static bool buildConsecutiveListOfVertices(GFace *gf, GEdgeLoop &gel,
BDS_GeomEntity *g = m->get_geom(ge->tag(), ge->dim());
pp->g = g;
if(MYDEBUG)
- printf("point %3d (%8.5f %8.5f : %8.5f %8.5f) (%2d,%2d)\n",
- count, pp->u, pp->v, param.x(), param.y(), pp->g->classif_tag,
+ printf("point %3d (%8.5f %8.5f : %8.5f %8.5f) (%2d,%2d)\n", count,
+ pp->u, pp->v, param.x(), param.y(), pp->g->classif_tag,
pp->g->classif_degree);
bbox += SPoint3(U, V, 0);
}
@@ -1902,15 +1959,13 @@ static bool buildConsecutiveListOfVertices(GFace *gf, GEdgeLoop &gel,
static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
{
-
- std::map<BDS_Point*, MVertex*, PointLessThan> recoverMap;
-
+ std::map<BDS_Point *, MVertex *, PointLessThan> recoverMap;
+
char name[245];
sprintf(name, "trueBoundary%d.pos", gf->tag());
std::vector<SPoint2> true_boundary;
trueBoundary(name, gf, true_boundary);
-
Range<double> rangeU = gf->parBounds(0);
Range<double> rangeV = gf->parBounds(1);
@@ -1923,31 +1978,31 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
// procedure
BDS_Mesh *m = new BDS_Mesh;
- std::vector<std::vector<BDS_Point*> > edgeLoops_BDS;
+ std::vector<std::vector<BDS_Point *> > edgeLoops_BDS;
SBoundingBox3d bbox;
int nbPointsTotal = 0;
{
for(std::list<GEdgeLoop>::iterator it = gf->edgeLoops.begin();
- it != gf->edgeLoops.end(); it++){
- std::vector<BDS_Point* > edgeLoop_BDS;
+ it != gf->edgeLoops.end(); it++) {
+ std::vector<BDS_Point *> edgeLoop_BDS;
int nbPointsLocal;
const double fact[4] = {1.e-12, 1.e-7, 1.e-5, 1.e-3};
bool ok = false;
- for(int i = 0; i < 4; i++){
+ for(int i = 0; i < 4; i++) {
if(buildConsecutiveListOfVertices(gf, *it, edgeLoop_BDS, bbox, m,
recoverMap, nbPointsLocal,
- nbPointsTotal, fact[i] * LC2D)){
+ nbPointsTotal, fact[i] * LC2D)) {
ok = true;
break;
}
- if(buildConsecutiveListOfVertices(gf, *it, edgeLoop_BDS, bbox, m,
- recoverMap, nbPointsLocal,
- nbPointsTotal, fact[i] * LC2D, true)){
+ if(buildConsecutiveListOfVertices(
+ gf, *it, edgeLoop_BDS, bbox, m, recoverMap, nbPointsLocal,
+ nbPointsTotal, fact[i] * LC2D, true)) {
ok = true;
break;
}
}
- if(!ok){
+ if(!ok) {
gf->meshStatistics.status = GFace::FAILED;
Msg::Error("The 1D mesh seems not to be forming a closed loop");
delete m;
@@ -1958,7 +2013,7 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
}
}
- if(nbPointsTotal < 3){
+ if(nbPointsTotal < 3) {
Msg::Warning("Mesh Generation of Model Face %d Skipped: "
"Only %d Mesh Vertices on The Contours",
gf->tag(), nbPointsTotal);
@@ -1966,11 +2021,12 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
delete m;
return true;
}
- if(nbPointsTotal == 3){
+ if(nbPointsTotal == 3) {
MVertex *vv[3];
int i = 0;
- for(std::map<BDS_Point*, MVertex*, PointLessThan>::iterator it = recoverMap.begin();
- it != recoverMap.end(); it++){
+ for(std::map<BDS_Point *, MVertex *, PointLessThan>::iterator it =
+ recoverMap.begin();
+ it != recoverMap.end(); it++) {
vv[i++] = it->second;
}
gf->triangles.push_back(new MTriangle(vv[0], vv[1], vv[2]));
@@ -1979,64 +2035,65 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
return true;
}
-
// Use a divide & conquer type algorithm to create a triangulation. We add to
// the triangulation a box with 4 points that encloses the domain.
std::vector<int> edgesEmbedded;
-#if 1 //OLD_CODE_DELAUNAY
+#if 1 // OLD_CODE_DELAUNAY
{
int count = 0;
/////////////////////////////////////////////////////////////////
// Embedded Vertices
// add embedded vertices
- std::list<GVertex*> emb_vertx = gf->embeddedVertices();
- std::list<GVertex*>::iterator itvx = emb_vertx.begin();
+ std::list<GVertex *> emb_vertx = gf->embeddedVertices();
+ std::list<GVertex *>::iterator itvx = emb_vertx.begin();
- std::map<MVertex*, std::set<BDS_Point*> > invertedRecoverMap;
- for (std::map<BDS_Point*, MVertex*, PointLessThan>::iterator it = recoverMap.begin();
- it != recoverMap.end(); it++){
+ std::map<MVertex *, std::set<BDS_Point *> > invertedRecoverMap;
+ for(std::map<BDS_Point *, MVertex *, PointLessThan>::iterator it =
+ recoverMap.begin();
+ it != recoverMap.end(); it++) {
invertedRecoverMap[it->second].insert(it->first);
}
int pNum = m->MAXPOINTNUMBER;
- nbPointsTotal += emb_vertx.size();
+ nbPointsTotal += emb_vertx.size();
{
- std::vector<GEdge*> const& emb_edges = gf->embeddedEdges();
- std::vector<GEdge*>::const_iterator ite = emb_edges.begin();
- std::set<MVertex*> vs;
- while(ite != emb_edges.end()){
- for(unsigned int i = 0; i< (*ite)->lines.size(); i++){
- for(unsigned int j = 0; j< 2; j++){
- MVertex *v = (*ite)->lines[i]->getVertex(j);
- if (invertedRecoverMap.find(v) == invertedRecoverMap.end() && vs.find(v) == vs.end()){
- vs.insert(v);
- }
- }
- }
- ++ite;
+ std::vector<GEdge *> const &emb_edges = gf->embeddedEdges();
+ std::vector<GEdge *>::const_iterator ite = emb_edges.begin();
+ std::set<MVertex *> vs;
+ while(ite != emb_edges.end()) {
+ for(unsigned int i = 0; i < (*ite)->lines.size(); i++) {
+ for(unsigned int j = 0; j < 2; j++) {
+ MVertex *v = (*ite)->lines[i]->getVertex(j);
+ if(invertedRecoverMap.find(v) == invertedRecoverMap.end() &&
+ vs.find(v) == vs.end()) {
+ vs.insert(v);
+ }
+ }
+ }
+ ++ite;
}
- nbPointsTotal += vs.size();
+ nbPointsTotal += vs.size();
}
- DocRecord doc(nbPointsTotal + 4 );
+ DocRecord doc(nbPointsTotal + 4);
- while(itvx != emb_vertx.end()){
+ while(itvx != emb_vertx.end()) {
MVertex *v = (*itvx)->mesh_vertices[0];
- double uv[2]={0,0};
- GPoint gp = gf->closestPoint (SPoint3(v->x(),v->y(),v->z()),uv);
+ double uv[2] = {0, 0};
+ GPoint gp = gf->closestPoint(SPoint3(v->x(), v->y(), v->z()), uv);
BDS_Point *pp = m->add_point(++pNum, gp.u(), gp.v(), gf);
- m->add_geom(-(*itvx)->tag(),0);
- pp->g = m->get_geom(-(*itvx)->tag(),0);
- pp->lcBGM() = BGM_MeshSize(*itvx, 0, 0, v->x(),v->y(),v->z());
+ m->add_geom(-(*itvx)->tag(), 0);
+ pp->g = m->get_geom(-(*itvx)->tag(), 0);
+ pp->lcBGM() = BGM_MeshSize(*itvx, 0, 0, v->x(), v->y(), v->z());
pp->lc() = pp->lcBGM();
// printf("%g\n",pp->lc());
recoverMap[pp] = v;
double XX = CTX::instance()->mesh.randFactor * LC2D * (double)rand() /
- (double)RAND_MAX;
+ (double)RAND_MAX;
double YY = CTX::instance()->mesh.randFactor * LC2D * (double)rand() /
- (double)RAND_MAX;
+ (double)RAND_MAX;
doc.points[count].where.h = pp->u + XX;
doc.points[count].where.v = pp->v + YY;
doc.points[count].adjacent = NULL;
@@ -2046,104 +2103,109 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
}
// nbPointsTotal += count;
- std::vector<GEdge*> const& emb_edges = gf->embeddedEdges();
- std::vector<GEdge*>::const_iterator ite = emb_edges.begin();
- std::set<MVertex*> vs;
+ std::vector<GEdge *> const &emb_edges = gf->embeddedEdges();
+ std::vector<GEdge *>::const_iterator ite = emb_edges.begin();
+ std::set<MVertex *> vs;
std::map<MVertex *, BDS_Point *> facile;
- while(ite != emb_edges.end()){
- m->add_geom(-(*ite)->tag(),1);
- for(unsigned int i = 0; i< (*ite)->lines.size(); i++){
- for(unsigned int j = 0; j< 2; j++){
- MVertex *v = (*ite)->lines[i]->getVertex(j);
- BDS_Point *pp = 0;
- const std::map<MVertex*, std::set<BDS_Point*> >::iterator it = invertedRecoverMap.find(v);
- if (it != invertedRecoverMap.end())
- {
- if (it->second.size() > 1) {
- const GEdge *edge = (*ite);
- const Range<double> parBounds = edge->parBounds(0);
- GPoint firstPoint = edge->point(parBounds.low());
- GPoint lastPoint = edge->point(parBounds.high());
- double param;
- if (v->point().distance(SPoint3(firstPoint.x(), firstPoint.y(), firstPoint.z()))
- < v->point().distance(SPoint3(lastPoint.x(), lastPoint.y(), lastPoint.z()))) {
- // Vertex lies on first point of edge
- param = parBounds.low();
- }
- else {
- // Vertex lies on last point of edge
- param = parBounds.high();
- }
- SPoint2 pointOnSurface = edge->reparamOnFace(gf, param, 1);
-
- const std::set<BDS_Point*> &possiblePoints = it->second;
- for (std::set<BDS_Point*>::iterator pntIt = possiblePoints.begin(); pntIt != possiblePoints.end();
- ++pntIt) {
- if (pointOnSurface.distance(SPoint2((*pntIt)->u, (*pntIt)->v)) < 1e-10) {
- pp = (*pntIt);
- break;
+ while(ite != emb_edges.end()) {
+ m->add_geom(-(*ite)->tag(), 1);
+ for(unsigned int i = 0; i < (*ite)->lines.size(); i++) {
+ for(unsigned int j = 0; j < 2; j++) {
+ MVertex *v = (*ite)->lines[i]->getVertex(j);
+ BDS_Point *pp = 0;
+ const std::map<MVertex *, std::set<BDS_Point *> >::iterator it =
+ invertedRecoverMap.find(v);
+ if(it != invertedRecoverMap.end()) {
+ if(it->second.size() > 1) {
+ const GEdge *edge = (*ite);
+ const Range<double> parBounds = edge->parBounds(0);
+ GPoint firstPoint = edge->point(parBounds.low());
+ GPoint lastPoint = edge->point(parBounds.high());
+ double param;
+ if(v->point().distance(
+ SPoint3(firstPoint.x(), firstPoint.y(), firstPoint.z())) <
+ v->point().distance(
+ SPoint3(lastPoint.x(), lastPoint.y(), lastPoint.z()))) {
+ // Vertex lies on first point of edge
+ param = parBounds.low();
+ }
+ else {
+ // Vertex lies on last point of edge
+ param = parBounds.high();
+ }
+ SPoint2 pointOnSurface = edge->reparamOnFace(gf, param, 1);
+
+ const std::set<BDS_Point *> &possiblePoints = it->second;
+ for(std::set<BDS_Point *>::iterator pntIt =
+ possiblePoints.begin();
+ pntIt != possiblePoints.end(); ++pntIt) {
+ if(pointOnSurface.distance(SPoint2((*pntIt)->u, (*pntIt)->v)) <
+ 1e-10) {
+ pp = (*pntIt);
+ break;
+ }
+ }
+ if(pp == 0) {
+ Msg::Error("Embedded edge vertex %d is on the seam edge of "
+ "surface %d and no appropriate point could be "
+ "found!\n",
+ v->getNum(), gf->tag());
+ }
+ }
+ else {
+ pp = *(it->second.begin());
+ }
+ facile[v] = pp;
+ }
+ if(pp == 0 && vs.find(v) == vs.end()) {
+ vs.insert(v);
+ double uv[2] = {0, 0};
+ GPoint gp = gf->closestPoint(SPoint3(v->x(), v->y(), v->z()), uv);
+ BDS_Point *pp = m->add_point(++pNum, gp.u(), gp.v(), gf);
+ pp->g = m->get_geom(-(*ite)->tag(), 1);
+ if(v->onWhat()->dim() == 0)
+ pp->lcBGM() =
+ BGM_MeshSize(v->onWhat(), 0, 0, v->x(), v->y(), v->z());
+ else {
+ double uu;
+ v->getParameter(0, uu);
+ pp->lcBGM() = BGM_MeshSize(*ite, uu, 0, v->x(), v->y(), v->z());
+ }
+ pp->lc() = pp->lcBGM();
+ // printf("%g\n",pp->lc());
+ recoverMap[pp] = v;
+ facile[v] = pp;
+ double XX = CTX::instance()->mesh.randFactor * LC2D *
+ (double)rand() / (double)RAND_MAX;
+ double YY = CTX::instance()->mesh.randFactor * LC2D *
+ (double)rand() / (double)RAND_MAX;
+ doc.points[count].where.h = pp->u + XX;
+ doc.points[count].where.v = pp->v + YY;
+ doc.points[count].adjacent = NULL;
+ doc.points[count].data = pp;
+ count++;
}
}
- if (pp == 0){
- Msg::Error("Embedded edge vertex %d is on the seam edge of surface %d and no appropriate point could be found!\n",
- v->getNum(), gf->tag());
- }
- }
- else {
- pp = *(it->second.begin());
- }
- facile[v] = pp;
- }
- if (pp == 0 && vs.find(v) == vs.end()){
- vs.insert(v);
- double uv[2]={0,0};
- GPoint gp = gf->closestPoint (SPoint3(v->x(),v->y(),v->z()),uv);
- BDS_Point *pp = m->add_point(++pNum, gp.u(), gp.v(), gf);
- pp->g = m->get_geom(-(*ite)->tag(),1);
- if (v->onWhat()->dim() == 0)
- pp->lcBGM() = BGM_MeshSize(v->onWhat(), 0, 0, v->x(),v->y(),v->z());
- else {
- double uu;
- v->getParameter(0,uu);
- pp->lcBGM() = BGM_MeshSize(*ite, uu, 0, v->x(),v->y(),v->z());
- }
- pp->lc() = pp->lcBGM();
- // printf("%g\n",pp->lc());
- recoverMap[pp] = v;
- facile[v] = pp;
- double XX = CTX::instance()->mesh.randFactor * LC2D * (double)rand() /
- (double)RAND_MAX;
- double YY = CTX::instance()->mesh.randFactor * LC2D * (double)rand() /
- (double)RAND_MAX;
- doc.points[count].where.h = pp->u + XX;
- doc.points[count].where.v = pp->v + YY;
- doc.points[count].adjacent = NULL;
- doc.points[count].data = pp;
- count++;
- }
- }
}
- for(unsigned int i = 0; i< (*ite)->lines.size(); i++){
- BDS_Point *p0 = facile[(*ite)->lines[i]->getVertex(0)];
- BDS_Point *p1 = facile[(*ite)->lines[i]->getVertex(1)];
- edgesEmbedded.push_back(p0->iD);
- edgesEmbedded.push_back(p1->iD);
+ for(unsigned int i = 0; i < (*ite)->lines.size(); i++) {
+ BDS_Point *p0 = facile[(*ite)->lines[i]->getVertex(0)];
+ BDS_Point *p1 = facile[(*ite)->lines[i]->getVertex(1)];
+ edgesEmbedded.push_back(p0->iD);
+ edgesEmbedded.push_back(p1->iD);
}
++ite;
}
/////////////////////////////////////////////////////////////////
-
-
- for(unsigned int i = 0; i < edgeLoops_BDS.size(); i++){
- std::vector<BDS_Point*> &edgeLoop_BDS = edgeLoops_BDS[i];
- for(unsigned int j = 0; j < edgeLoop_BDS.size(); j++){
+ for(unsigned int i = 0; i < edgeLoops_BDS.size(); i++) {
+ std::vector<BDS_Point *> &edgeLoop_BDS = edgeLoops_BDS[i];
+ for(unsigned int j = 0; j < edgeLoop_BDS.size(); j++) {
BDS_Point *pp = edgeLoop_BDS[j];
double XX = CTX::instance()->mesh.randFactor * LC2D * (double)rand() /
- (double)RAND_MAX;
+ (double)RAND_MAX;
double YY = CTX::instance()->mesh.randFactor * LC2D * (double)rand() /
- (double)RAND_MAX;
+ (double)RAND_MAX;
doc.points[count].where.h = pp->u + XX;
doc.points[count].where.v = pp->v + YY;
doc.points[count].adjacent = NULL;
@@ -2156,7 +2218,7 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
// the domain, use negative number to distinguish those fake
// vertices
- if(du / dv < 1200 && dv / du < 1200){
+ if(du / dv < 1200 && dv / du < 1200) {
// FIX A BUG HERE IF THE SIZE OF THE BOX IS ZERO
bbox.makeCube();
}
@@ -2167,15 +2229,15 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
bb[1] = new MVertex(bbox.min().x(), bbox.max().y(), 0, 0, -2);
bb[2] = new MVertex(bbox.max().x(), bbox.min().y(), 0, 0, -3);
bb[3] = new MVertex(bbox.max().x(), bbox.max().y(), 0, 0, -4);
- for(int ip = 0; ip < 4; ip++){
+ for(int ip = 0; ip < 4; ip++) {
BDS_Point *pp = m->add_point(-ip - 1, bb[ip]->x(), bb[ip]->y(), gf);
m->add_geom(gf->tag(), 2);
BDS_GeomEntity *g = m->get_geom(gf->tag(), 2);
pp->g = g;
- doc.points[nbPointsTotal+ip].where.h = bb[ip]->x();
- doc.points[nbPointsTotal+ip].where.v = bb[ip]->y();
- doc.points[nbPointsTotal+ip].adjacent = 0;
- doc.points[nbPointsTotal+ip].data = pp;
+ doc.points[nbPointsTotal + ip].where.h = bb[ip]->x();
+ doc.points[nbPointsTotal + ip].where.v = bb[ip]->y();
+ doc.points[nbPointsTotal + ip].adjacent = 0;
+ doc.points[nbPointsTotal + ip].data = pp;
}
for(int ip = 0; ip < 4; ip++) delete bb[ip];
@@ -2186,25 +2248,24 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
// loop is the outer one)
Msg::Debug("Meshing of the convex hull (%d points)", nbPointsTotal);
- try{
+ try {
doc.MakeMeshWithPoints();
- }
- catch(const char *err){
+ } catch(const char *err) {
Msg::Error("%s", err);
}
- for(int i = 0; i < doc.numTriangles; i++){
+ for(int i = 0; i < doc.numTriangles; i++) {
int a = doc.triangles[i].a;
int b = doc.triangles[i].b;
int c = doc.triangles[i].c;
int n = doc.numPoints;
- if(a < 0 || a >= n || b < 0 || b >= n || c < 0 || c >= n){
+ if(a < 0 || a >= n || b < 0 || b >= n || c < 0 || c >= n) {
Msg::Warning("Skipping bad triangle %d", i);
continue;
}
- BDS_Point *p1 = (BDS_Point*)doc.points[doc.triangles[i].a].data;
- BDS_Point *p2 = (BDS_Point*)doc.points[doc.triangles[i].b].data;
- BDS_Point *p3 = (BDS_Point*)doc.points[doc.triangles[i].c].data;
+ BDS_Point *p1 = (BDS_Point *)doc.points[doc.triangles[i].a].data;
+ BDS_Point *p2 = (BDS_Point *)doc.points[doc.triangles[i].b].data;
+ BDS_Point *p3 = (BDS_Point *)doc.points[doc.triangles[i].c].data;
m->add_triangle(p1->iD, p2->iD, p3->iD);
}
}
@@ -2212,11 +2273,11 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
{
/// FIXME FOR PERIODIC : some MVertices should be duplicated...
/// Still to be done...
- std::vector<MVertex*> v;
- std::map<MVertex*, BDS_Point*> recoverMapInv;
- for(unsigned int i = 0; i < edgeLoops_BDS.size(); i++){
- std::vector<BDS_Point*> &edgeLoop_BDS = edgeLoops_BDS[i];
- for(unsigned int j = 0; j < edgeLoop_BDS.size(); j++){
+ std::vector<MVertex *> v;
+ std::map<MVertex *, BDS_Point *> recoverMapInv;
+ for(unsigned int i = 0; i < edgeLoops_BDS.size(); i++) {
+ std::vector<BDS_Point *> &edgeLoop_BDS = edgeLoops_BDS[i];
+ for(unsigned int j = 0; j < edgeLoop_BDS.size(); j++) {
BDS_Point *pp = edgeLoop_BDS[j];
v.push_back(recoverMap[pp]);
recoverMapInv[recoverMap[pp]] = pp;
@@ -2224,26 +2285,26 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
}
// printf("coucou2 %d verices\n",v.size());
- std::map<MVertex*,SPoint3> pos;
+ std::map<MVertex *, SPoint3> pos;
for(unsigned int i = 0; i < v.size(); i++) {
MVertex *v0 = v[i];
- BDS_Point *p0 = recoverMapInv[v0];
- pos[v0] = SPoint3(v0->x(),v0->y(),v0->z());
- v0->setXYZ(p0->u,p0->v,0.0);
+ BDS_Point *p0 = recoverMapInv[v0];
+ pos[v0] = SPoint3(v0->x(), v0->y(), v0->z());
+ v0->setXYZ(p0->u, p0->v, 0.0);
}
- std::vector<MTriangle*> result;
+ std::vector<MTriangle *> result;
delaunayMeshIn2D(v, result, 0);
- for(unsigned int i = 0; i < v.size()-4; i++) {
+ for(unsigned int i = 0; i < v.size() - 4; i++) {
MVertex *v0 = v[i];
SPoint3 pp = pos[v0];
- v0->setXYZ(pp.x(),pp.y(),pp.z());
+ v0->setXYZ(pp.x(), pp.y(), pp.z());
}
// add the four corners
- for(int ip = 0; ip < 4; ip++){
- MVertex *vv = v[v.size()-ip-1];
- BDS_Point *pp = m->add_point(-ip - 1, vv->x(),vv->y(), gf);
+ for(int ip = 0; ip < 4; ip++) {
+ MVertex *vv = v[v.size() - ip - 1];
+ BDS_Point *pp = m->add_point(-ip - 1, vv->x(), vv->y(), gf);
m->add_geom(gf->tag(), 2);
recoverMapInv[vv] = pp;
BDS_GeomEntity *g = m->get_geom(gf->tag(), 2);
@@ -2254,77 +2315,76 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
MVertex *v0 = result[i]->getVertex(0);
MVertex *v1 = result[i]->getVertex(1);
MVertex *v2 = result[i]->getVertex(2);
- BDS_Point *p0 = recoverMapInv[v0];
- BDS_Point *p1 = recoverMapInv[v1];
- BDS_Point *p2 = recoverMapInv[v2];
+ BDS_Point *p0 = recoverMapInv[v0];
+ BDS_Point *p1 = recoverMapInv[v1];
+ BDS_Point *p2 = recoverMapInv[v2];
m->add_triangle(p0->iD, p1->iD, p2->iD);
}
}
#endif
-
// Recover the boundary edges and compute characteristic lenghts using mesh
// edge spacing
BDS_GeomEntity CLASS_F(1, 2);
BDS_GeomEntity CLASS_E(1, 1);
BDS_GeomEntity CLASS_EXTERIOR(3, 2);
- if(debug){
+ if(debug) {
char name[245];
sprintf(name, "surface%d-initial-real.pos", gf->tag());
- outputScalarField(m->triangles, name, 0,gf);
+ outputScalarField(m->triangles, name, 0, gf);
sprintf(name, "surface%d-initial-param.pos", gf->tag());
- outputScalarField(m->triangles, name, 1,gf);
+ outputScalarField(m->triangles, name, 1, gf);
}
bool _fatallyFailed;
- for(unsigned int i = 0; i < edgesEmbedded.size()/2; i++){
- BDS_Edge * e = m->recover_edge
- (edgesEmbedded[2*i], edgesEmbedded[2*i+1], _fatallyFailed);
- if(!e){
- Msg::Error("Impossible to recover the edge %d %d",
- edgesEmbedded[2*i], edgesEmbedded[2*i+1]);
+ for(unsigned int i = 0; i < edgesEmbedded.size() / 2; i++) {
+ BDS_Edge *e = m->recover_edge(edgesEmbedded[2 * i],
+ edgesEmbedded[2 * i + 1], _fatallyFailed);
+ if(!e) {
+ Msg::Error("Impossible to recover the edge %d %d", edgesEmbedded[2 * i],
+ edgesEmbedded[2 * i + 1]);
gf->meshStatistics.status = GFace::FAILED;
delete m;
return false;
}
- else e->g = &CLASS_E;
+ else
+ e->g = &CLASS_E;
}
- for(unsigned int i = 0; i < edgeLoops_BDS.size(); i++){
- std::vector<BDS_Point*> &edgeLoop_BDS = edgeLoops_BDS[i];
- for(unsigned int j = 0; j < edgeLoop_BDS.size(); j++){
- BDS_Edge * e = m->recover_edge
- (edgeLoop_BDS[j]->iD, edgeLoop_BDS[(j + 1) % edgeLoop_BDS.size()]->iD, _fatallyFailed);
- if(!e){
+ for(unsigned int i = 0; i < edgeLoops_BDS.size(); i++) {
+ std::vector<BDS_Point *> &edgeLoop_BDS = edgeLoops_BDS[i];
+ for(unsigned int j = 0; j < edgeLoop_BDS.size(); j++) {
+ BDS_Edge *e = m->recover_edge(
+ edgeLoop_BDS[j]->iD, edgeLoop_BDS[(j + 1) % edgeLoop_BDS.size()]->iD,
+ _fatallyFailed);
+ if(!e) {
Msg::Error("Impossible to recover the edge %d %d", edgeLoop_BDS[j]->iD,
edgeLoop_BDS[(j + 1) % edgeLoop_BDS.size()]->iD);
gf->meshStatistics.status = GFace::FAILED;
delete m;
return false;
}
- else e->g = &CLASS_E;
+ else
+ e->g = &CLASS_E;
}
}
-
-
// look for a triangle that has a negative node and recursively tag all
// exterior triangles
{
- std::vector<BDS_Face*>::iterator itt = m->triangles.begin();
- while(itt != m->triangles.end()){
+ std::vector<BDS_Face *>::iterator itt = m->triangles.begin();
+ while(itt != m->triangles.end()) {
(*itt)->g = 0;
++itt;
}
itt = m->triangles.begin();
- while(itt != m->triangles.end()){
+ while(itt != m->triangles.end()) {
BDS_Face *t = *itt;
BDS_Point *n[4];
t->getNodes(n);
- if(n[0]->iD < 0 || n[1]->iD < 0 ||
- n[2]->iD < 0 ) {
+ if(n[0]->iD < 0 || n[1]->iD < 0 || n[2]->iD < 0) {
recur_tag(t, &CLASS_EXTERIOR);
break;
}
@@ -2334,23 +2394,23 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
// now find an edge that has belongs to one of the exterior triangles
{
- std::vector<BDS_Edge*>::const_iterator ite = m->edges.begin();
- while(ite != m->edges.end()){
+ std::vector<BDS_Edge *>::const_iterator ite = m->edges.begin();
+ while(ite != m->edges.end()) {
BDS_Edge *edge = *ite;
- if(edge->g && edge->numfaces() == 2){
- if(edge->faces(0)->g == &CLASS_EXTERIOR){
+ if(edge->g && edge->numfaces() == 2) {
+ if(edge->faces(0)->g == &CLASS_EXTERIOR) {
recur_tag(edge->faces(1), &CLASS_F);
break;
}
- else if(edge->faces(1)->g == &CLASS_EXTERIOR){
+ else if(edge->faces(1)->g == &CLASS_EXTERIOR) {
recur_tag(edge->faces(0), &CLASS_F);
break;
}
}
++ite;
}
- std::vector<BDS_Face*>::iterator itt = m->triangles.begin();
- while(itt != m->triangles.end()){
+ std::vector<BDS_Face *>::iterator itt = m->triangles.begin();
+ while(itt != m->triangles.end()) {
if((*itt)->g == &CLASS_EXTERIOR) (*itt)->g = 0;
++itt;
}
@@ -2358,11 +2418,11 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
// delete useless stuff
{
- std::vector<BDS_Face*>::iterator itt = m->triangles.begin();
- while(itt != m->triangles.end()){
+ std::vector<BDS_Face *>::iterator itt = m->triangles.begin();
+ while(itt != m->triangles.end()) {
BDS_Face *t = *itt;
- if(!t->g){
- m->del_face (t);
+ if(!t->g) {
+ m->del_face(t);
}
++itt;
}
@@ -2371,17 +2431,18 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
m->cleanup();
{
- std::vector<BDS_Edge*>::const_iterator ite = m->edges.begin();
- while(ite != m->edges.end()){
+ std::vector<BDS_Edge *>::const_iterator ite = m->edges.begin();
+ while(ite != m->edges.end()) {
BDS_Edge *edge = *ite;
if(edge->numfaces() == 0)
m->del_edge(edge);
- else{
- if(!edge->g)
- edge->g = &CLASS_F;
- if(!edge->p1->g || edge->p1->g->classif_degree > edge->g->classif_degree)
+ else {
+ if(!edge->g) edge->g = &CLASS_F;
+ if(!edge->p1->g ||
+ edge->p1->g->classif_degree > edge->g->classif_degree)
edge->p1->g = edge->g;
- if(!edge->p2->g || edge->p2->g->classif_degree > edge->g->classif_degree)
+ if(!edge->p2->g ||
+ edge->p2->g->classif_degree > edge->g->classif_degree)
edge->p2->g = edge->g;
}
++ite;
@@ -2393,19 +2454,17 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
m->del_point(m->find_point(-3));
m->del_point(m->find_point(-4));
- if(debug){
+ if(debug) {
char name[245];
sprintf(name, "surface%d-recovered-real.pos", gf->tag());
- outputScalarField(m->triangles, name, 0,gf);
+ outputScalarField(m->triangles, name, 0, gf);
sprintf(name, "surface%d-recovered-param.pos", gf->tag());
- outputScalarField(m->triangles, name, 1,gf);
+ outputScalarField(m->triangles, name, 1, gf);
}
-
// start mesh generation for periodic face
- if(!algoDelaunay2D(gf)){
-
+ if(!algoDelaunay2D(gf)) {
// optimizeMeshBDS(gf, *m, 2, &recoverMap); // fix periodic shit
// if(debug){
// char name[245];
@@ -2415,17 +2474,18 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
// outputScalarField(m->triangles, name, 1);
// }
-
- // need for a BGM for cross field
- // if(CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine || 1) {
+ // need for a BGM for cross field
+ // if(CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine
+ // || 1) {
// printf("coucou here !!!\n");
// backgroundMesh::unset();
// buildBackGroundMesh(gf);
// }
- modifyInitialMeshToRemoveDegeneracies(gf, *m,&recoverMap);
+ modifyInitialMeshToRemoveDegeneracies(gf, *m, &recoverMap);
- refineMeshBDS(gf, *m, CTX::instance()->mesh.refineSteps, true,NULL,&recoverMap, &true_boundary);
- if(debug){
+ refineMeshBDS(gf, *m, CTX::instance()->mesh.refineSteps, true, NULL,
+ &recoverMap, &true_boundary);
+ if(debug) {
char name[245];
sprintf(name, "surface%d-phase1-real.pos", gf->tag());
outputScalarField(m->triangles, name, 0, gf);
@@ -2433,14 +2493,15 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
outputScalarField(m->triangles, name, 1, gf);
}
optimizeMeshBDS(gf, *m, 2);
- if(debug){
+ if(debug) {
char name[245];
sprintf(name, "surface%d-phase2-real.pos", gf->tag());
outputScalarField(m->triangles, name, 0, gf);
sprintf(name, "surface%d-phase2-param.pos", gf->tag());
outputScalarField(m->triangles, name, 1, gf);
}
- refineMeshBDS(gf, *m, -CTX::instance()->mesh.refineSteps, false,NULL,&recoverMap, &true_boundary);
+ refineMeshBDS(gf, *m, -CTX::instance()->mesh.refineSteps, false, NULL,
+ &recoverMap, &true_boundary);
// if(debug){
// char name[245];
// sprintf(name, "surface%d-phase3-real.pos", gf->tag());
@@ -2461,14 +2522,14 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
gf->meshStatistics.nbGoodLength);*/
gf->meshStatistics.status = GFace::DONE;
-
- if(debug){
+ if(debug) {
char name[245];
sprintf(name, "surface%d-just-real.pos", gf->tag());
outputScalarField(m->triangles, name, 0, gf);
}
- // if(CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine || 1) {
+ // if(CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine ||
+ // 1) {
// backgroundMesh::unset();
// }
}
@@ -2476,55 +2537,56 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
// This is a structure that we need only for periodic cases. We will duplicate
// the vertices (MVertex) that are on seams
- std::map<MVertex*, MVertex*> equivalence;
- std::map<MVertex*, SPoint2> parametricCoordinates;
- if(algoDelaunay2D(gf)){
- std::map<MVertex*, BDS_Point*> invertMap;
- std::map<BDS_Point*, MVertex*, PointLessThan>::iterator it = recoverMap.begin();
- while(it != recoverMap.end()){
+ std::map<MVertex *, MVertex *> equivalence;
+ std::map<MVertex *, SPoint2> parametricCoordinates;
+ if(algoDelaunay2D(gf)) {
+ std::map<MVertex *, BDS_Point *> invertMap;
+ std::map<BDS_Point *, MVertex *, PointLessThan>::iterator it =
+ recoverMap.begin();
+ while(it != recoverMap.end()) {
// we have twice vertex MVertex with 2 different coordinates
MVertex *mv1 = it->second;
BDS_Point *bds = it->first;
- std::map<MVertex*, BDS_Point*>::iterator invIt = invertMap.find(mv1);
- if(invIt != invertMap.end()){
+ std::map<MVertex *, BDS_Point *>::iterator invIt = invertMap.find(mv1);
+ if(invIt != invertMap.end()) {
// create a new "fake" vertex that will be destroyed afterwards
MVertex *mv2 = 0;
if(mv1->onWhat()->dim() == 1) {
double t;
- mv1->getParameter(0,t);
- mv2 = new MEdgeVertex(mv1->x(),mv1->y(),mv1->z(),mv1->onWhat(), t, 0,
- ((MEdgeVertex*)mv1)->getLc());
+ mv1->getParameter(0, t);
+ mv2 = new MEdgeVertex(mv1->x(), mv1->y(), mv1->z(), mv1->onWhat(), t,
+ 0, ((MEdgeVertex *)mv1)->getLc());
}
else if(mv1->onWhat()->dim() == 0) {
- mv2 = new MVertex(mv1->x(),mv1->y(),mv1->z(),mv1->onWhat());
+ mv2 = new MVertex(mv1->x(), mv1->y(), mv1->z(), mv1->onWhat());
}
else
Msg::Error("Could not reconstruct seam");
- if(mv2){
+ if(mv2) {
it->second = mv2;
equivalence[mv2] = mv1;
- parametricCoordinates[mv2] = SPoint2(bds->u,bds->v);
+ parametricCoordinates[mv2] = SPoint2(bds->u, bds->v);
invertMap[mv2] = bds;
}
}
else {
- parametricCoordinates[mv1] = SPoint2(bds->u,bds->v);
+ parametricCoordinates[mv1] = SPoint2(bds->u, bds->v);
invertMap[mv1] = bds;
}
++it;
}
// recoverMap.insert(new_relations.begin(), new_relations.end());
}
- // Msg::Info("%d points that are duplicated for Delaunay meshing", equivalence.size());
+ // Msg::Info("%d points that are duplicated for Delaunay meshing",
+ // equivalence.size());
// fill the small gmsh structures
{
- std::set<BDS_Point*, PointLessThan>::iterator itp = m->points.begin();
- while(itp != m->points.end()){
+ std::set<BDS_Point *, PointLessThan>::iterator itp = m->points.begin();
+ while(itp != m->points.end()) {
BDS_Point *p = *itp;
- if(recoverMap.find(p) == recoverMap.end()){
- MVertex *v = new MFaceVertex
- (p->X, p->Y, p->Z, gf, p->u, p->v);
+ if(recoverMap.find(p) == recoverMap.end()) {
+ MVertex *v = new MFaceVertex(p->X, p->Y, p->Z, gf, p->u, p->v);
recoverMap[p] = v;
gf->mesh_vertices.push_back(v);
}
@@ -2532,27 +2594,27 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
}
}
- std::map<MTriangle*, BDS_Face*> invert_map;
+ std::map<MTriangle *, BDS_Face *> invert_map;
{
- std::vector<BDS_Face*>::iterator itt = m->triangles.begin();
- while(itt != m->triangles.end()){
+ std::vector<BDS_Face *>::iterator itt = m->triangles.begin();
+ while(itt != m->triangles.end()) {
BDS_Face *t = *itt;
- if(!t->deleted){
+ if(!t->deleted) {
BDS_Point *n[4];
t->getNodes(n);
MVertex *v1 = recoverMap[n[0]];
MVertex *v2 = recoverMap[n[1]];
MVertex *v3 = recoverMap[n[2]];
- if(!n[3]){
+ if(!n[3]) {
// when a singular point is present, degenerated triangles may be
// created, for example on a sphere that contains one pole
- if(v1 != v2 && v1 != v3 && v2 != v3){
+ if(v1 != v2 && v1 != v3 && v2 != v3) {
// we are in the periodic case. if we aim at using delaunay mesh
// generation in thoses cases, we should double some of the vertices
gf->triangles.push_back(new MTriangle(v1, v2, v3));
}
}
- else{
+ else {
MVertex *v4 = recoverMap[n[3]];
gf->quadrangles.push_back(new MQuadrangle(v1, v2, v3, v4));
}
@@ -2561,7 +2623,7 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
}
}
- if(debug){
+ if(debug) {
char name[245];
sprintf(name, "surface%d-final-real.pos", gf->tag());
outputScalarField(m->triangles, name, 0, gf);
@@ -2575,52 +2637,54 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
gf->getMeshingAlgo() == ALGO_2D_PACK_PRLGRMS_CSTR)
infty = true;
- if(infty)
- buildBackGroundMesh(gf, &equivalence, ¶metricCoordinates);
+ if(infty) buildBackGroundMesh(gf, &equivalence, ¶metricCoordinates);
// boundary layer
- std::vector<MQuadrangle*> blQuads;
- std::vector<MTriangle*> blTris;
- std::set<MVertex*> verts;
+ std::vector<MQuadrangle *> blQuads;
+ std::vector<MTriangle *> blTris;
+ std::set<MVertex *> verts;
modifyInitialMeshForBoundaryLayers(gf, blQuads, blTris, verts, debug);
- gf->quadrangles.insert(gf->quadrangles.begin(), blQuads.begin(), blQuads.end());
+ gf->quadrangles.insert(gf->quadrangles.begin(), blQuads.begin(),
+ blQuads.end());
gf->triangles.insert(gf->triangles.begin(), blTris.begin(), blTris.end());
- gf->mesh_vertices.insert(gf->mesh_vertices.begin(), verts.begin(), verts.end());
+ gf->mesh_vertices.insert(gf->mesh_vertices.begin(), verts.begin(),
+ verts.end());
- if(algoDelaunay2D(gf)){
+ if(algoDelaunay2D(gf)) {
if(gf->getMeshingAlgo() == ALGO_2D_FRONTAL)
- bowyerWatsonFrontal(gf, &equivalence, ¶metricCoordinates, &true_boundary);
+ bowyerWatsonFrontal(gf, &equivalence, ¶metricCoordinates,
+ &true_boundary);
else if(gf->getMeshingAlgo() == ALGO_2D_FRONTAL_QUAD)
- bowyerWatsonFrontalLayers(gf,true, &equivalence, ¶metricCoordinates);
+ bowyerWatsonFrontalLayers(gf, true, &equivalence, ¶metricCoordinates);
else if(gf->getMeshingAlgo() == ALGO_2D_PACK_PRLGRMS)
- bowyerWatsonParallelograms(gf,&equivalence, ¶metricCoordinates);
+ bowyerWatsonParallelograms(gf, &equivalence, ¶metricCoordinates);
else if(gf->getMeshingAlgo() == ALGO_2D_PACK_PRLGRMS_CSTR)
- bowyerWatsonParallelogramsConstrained(gf, gf->constr_vertices, &equivalence,
- ¶metricCoordinates);
+ bowyerWatsonParallelogramsConstrained(
+ gf, gf->constr_vertices, &equivalence, ¶metricCoordinates);
else if(gf->getMeshingAlgo() == ALGO_2D_DELAUNAY ||
gf->getMeshingAlgo() == ALGO_2D_AUTO)
bowyerWatson(gf, 1000000000, &equivalence, ¶metricCoordinates);
else
meshGFaceBamg(gf);
- if(!infty || !(CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine))
+ if(!infty ||
+ !(CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine))
laplaceSmoothing(gf, CTX::instance()->mesh.nbSmoothing, infty);
}
- if(debug){
- char name[256];
- sprintf(name, "real%d.pos", gf->tag());
- outputScalarField(m->triangles, name, 0, gf);
- sprintf(name, "param%d.pos", gf->tag());
- outputScalarField(m->triangles, name,1, gf);
- }
-
+ if(debug) {
+ char name[256];
+ sprintf(name, "real%d.pos", gf->tag());
+ outputScalarField(m->triangles, name, 0, gf);
+ sprintf(name, "param%d.pos", gf->tag());
+ outputScalarField(m->triangles, name, 1, gf);
+ }
// delete the mesh
delete m;
if((CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine) &&
CTX::instance()->mesh.algoRecombine != 2)
- recombineIntoQuads(gf,true,false);
+ recombineIntoQuads(gf, true, false);
computeElementShapes(gf, gf->meshStatistics.worst_element_shape,
gf->meshStatistics.average_element_shape,
@@ -2633,7 +2697,8 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true)
void deMeshGFace::operator()(GFace *gf)
{
- if(gf->geomType() == GEntity::DiscreteSurface && !CTX::instance()->meshDiscrete)
+ if(gf->geomType() == GEntity::DiscreteSurface &&
+ !CTX::instance()->meshDiscrete)
return;
gf->deleteMesh();
gf->meshStatistics.status = GFace::PENDING;
@@ -2648,7 +2713,7 @@ void meshGFace::operator()(GFace *gf, bool print)
{
gf->model()->setCurrentMeshEntity(gf);
- if(debugSurface >= 0 && gf->tag() != debugSurface){
+ if(debugSurface >= 0 && gf->tag() != debugSurface) {
gf->meshStatistics.status = GFace::DONE;
return;
}
@@ -2666,10 +2731,10 @@ void meshGFace::operator()(GFace *gf, bool print)
// because meshGenerator never called
if(MeshTransfiniteSurface(gf)) return;
if(MeshExtrudedSurface(gf)) return;
- if(gf->meshMaster() != gf){
- GFace *gff = dynamic_cast<GFace*>(gf->meshMaster());
- if(gff){
- if(gff->meshStatistics.status != GFace::DONE){
+ if(gf->meshMaster() != gf) {
+ GFace *gff = dynamic_cast<GFace *>(gf->meshMaster());
+ if(gff) {
+ if(gff->meshStatistics.status != GFace::DONE) {
gf->meshStatistics.status = GFace::PENDING;
return;
}
@@ -2685,24 +2750,25 @@ void meshGFace::operator()(GFace *gf, bool print)
const char *algo = "Unknown";
- switch(gf->getMeshingAlgo()){
- case ALGO_2D_MESHADAPT : algo = "MeshAdapt"; break;
- case ALGO_2D_FRONTAL : algo = "Frontal"; break;
- case ALGO_2D_FRONTAL_QUAD : algo = "Frontal Quad"; break;
- case ALGO_2D_DELAUNAY : algo = "Delaunay"; break;
- case ALGO_2D_BAMG : algo = "Bamg"; break;
- case ALGO_2D_PACK_PRLGRMS : algo = "Square Packing"; break;
- case ALGO_2D_AUTO :
+ switch(gf->getMeshingAlgo()) {
+ case ALGO_2D_MESHADAPT: algo = "MeshAdapt"; break;
+ case ALGO_2D_FRONTAL: algo = "Frontal"; break;
+ case ALGO_2D_FRONTAL_QUAD: algo = "Frontal Quad"; break;
+ case ALGO_2D_DELAUNAY: algo = "Delaunay"; break;
+ case ALGO_2D_BAMG: algo = "Bamg"; break;
+ case ALGO_2D_PACK_PRLGRMS: algo = "Square Packing"; break;
+ case ALGO_2D_AUTO:
algo = (gf->geomType() == GEntity::Plane) ? "Delaunay" : "MeshAdapt";
break;
}
- if(!algoDelaunay2D(gf)){
+ if(!algoDelaunay2D(gf)) {
algo = "MeshAdapt";
}
if(print)
- Msg::Info("Meshing surface %d (%s, %s)", gf->tag(), gf->getTypeString().c_str(), algo);
+ Msg::Info("Meshing surface %d (%s, %s)", gf->tag(),
+ gf->getTypeString().c_str(), algo);
// compute loops on the fly (indices indicate start and end points of a loop;
// loops are not yet oriented)
@@ -2713,16 +2779,16 @@ void meshGFace::operator()(GFace *gf, bool print)
quadMeshRemoveHalfOfOneDMesh halfmesh(gf);
bool singularEdges = false;
- std::vector<GEdge*>::const_iterator ite = gf->edges().begin();
- while(ite != gf->edges().end()){
- if((*ite)->isSeam(gf)) singularEdges= true;
- if((*ite)->isMeshDegenerated())singularEdges = true;
+ std::vector<GEdge *>::const_iterator ite = gf->edges().begin();
+ while(ite != gf->edges().end()) {
+ if((*ite)->isSeam(gf)) singularEdges = true;
+ if((*ite)->isMeshDegenerated()) singularEdges = true;
ite++;
}
-
- if (gf->getNativeType() != GEntity::GmshModel && (gf->periodic(0) || gf->periodic(1) || singularEdges)){
- if(!meshGeneratorPeriodic
- (gf, debugSurface >= 0 || debugSurface == -100))
+
+ if(gf->getNativeType() != GEntity::GmshModel &&
+ (gf->periodic(0) || gf->periodic(1) || singularEdges)) {
+ if(!meshGeneratorPeriodic(gf, debugSurface >= 0 || debugSurface == -100))
Msg::Error("Impossible to mesh periodic face %d", gf->tag());
}
else {
@@ -2735,7 +2801,7 @@ void meshGFace::operator()(GFace *gf, bool print)
halfmesh.finish();
- if(gf->getNumMeshElements() == 0){
+ if(gf->getNumMeshElements() == 0) {
Msg::Warning("Surface %d consists of no elements", gf->tag());
}
}
@@ -2777,23 +2843,25 @@ static bool getGFaceNormalFromBary(GFace *gf, MElement *el, SVector3 &nf)
}
static void getGFaceOrientation(GFace *gf, BoundaryLayerColumns *blc,
- bool existBL, bool fromVert,
- int &orientNonBL, int &orientBL)
+ bool existBL, bool fromVert, int &orientNonBL,
+ int &orientBL)
{
for(unsigned int iEl = 0; iEl < gf->getNumMeshElements(); iEl++) {
MElement *e = gf->getMeshElement(iEl);
- const bool isBLEl = existBL &&
- (blc->_toFirst.find(e) != blc->_toFirst.end());
+ const bool isBLEl =
+ existBL && (blc->_toFirst.find(e) != blc->_toFirst.end());
SVector3 nf;
// Check only if orientation of BL/non-BL el. not already known
if((!isBLEl && orientNonBL == 0) || (isBLEl && orientBL == 0)) {
const bool found = fromVert ? getGFaceNormalFromVert(gf, e, nf) :
- getGFaceNormalFromBary(gf, e, nf);
+ getGFaceNormalFromBary(gf, e, nf);
if(found) {
SVector3 ne = e->getFace(0).normal();
const int orient = (dot(ne, nf) > 0.) ? 1 : -1;
- if(isBLEl) orientBL = orient;
- else orientNonBL = orient;
+ if(isBLEl)
+ orientBL = orient;
+ else
+ orientNonBL = orient;
}
}
// Stop when orientation found for non-BL and BL el.
@@ -2809,7 +2877,7 @@ void orientMeshGFace::operator()(GFace *gf)
gf->model()->setCurrentMeshEntity(gf);
if(gf->geomType() == GEntity::DiscreteSurface ||
- gf->geomType() == GEntity::BoundaryLayerSurface){
+ gf->geomType() == GEntity::BoundaryLayerSurface) {
// don't do anything
}
else {
@@ -2855,14 +2923,15 @@ void orientMeshGFace::operator()(GFace *gf)
if(orientNonBL == -1) e->reverse();
}
else // If el. in BL
- // ... reverse if needed
- if(orientBL == -1) e->reverse();
+ // ... reverse if needed
+ if(orientBL == -1)
+ e->reverse();
}
}
else // If no BL, reverse all elements if needed
if(orientNonBL == -1)
- for(unsigned int iEl = 0; iEl < gf->getNumMeshElements(); iEl++)
- gf->getMeshElement(iEl)->reverse();
+ for(unsigned int iEl = 0; iEl < gf->getNumMeshElements(); iEl++)
+ gf->getMeshElement(iEl)->reverse();
}
// Apply user-specified mesh orientation constraints
diff --git a/Mesh/meshGFace.h b/Mesh/meshGFace.h
index 395df17fd5a99eb7ed66ac155acd6509a058212b..43a81b4a15cc235762d9c9744c45d04c457ccdcc 100644
--- a/Mesh/meshGFace.h
+++ b/Mesh/meshGFace.h
@@ -25,8 +25,7 @@ class meshGFace {
public:
meshGFace(bool r = true)
- : repairSelfIntersecting1dMesh(r)
- , onlyInitialMesh(false)
+ : repairSelfIntersecting1dMesh(r), onlyInitialMesh(false)
{
}
void operator()(GFace *, bool print = true);
@@ -61,7 +60,9 @@ bool checkMeshCompound(GFaceCompound *gf, std::list<GEdge *> &edges);
bool meshGenerator(GFace *gf, int RECUR_ITER, bool repairSelfIntersecting1dMesh,
bool onlyInitialMesh, bool debug = true,
std::vector<GEdge *> *replacement_edges = 0);
-bool pointInsideParametricDomain (std::vector<SPoint2> &bnd, SPoint2 &p, SPoint2 &out, int &N);
-double validity_p2triangle_formula (double *xa, double *xb, double *xc, double *xab, double *xbc, double *xca);
+bool pointInsideParametricDomain(std::vector<SPoint2> &bnd, SPoint2 &p,
+ SPoint2 &out, int &N);
+double validity_p2triangle_formula(double *xa, double *xb, double *xc,
+ double *xab, double *xbc, double *xca);
#endif
diff --git a/Mesh/meshGFaceBDS.cpp b/Mesh/meshGFaceBDS.cpp
index 5ff290e3ea716d4a039013dac7ab0d45c50dee01..61d009ffbf7ae0d7d19f2e26452ffca61f1e198b 100644
--- a/Mesh/meshGFaceBDS.cpp
+++ b/Mesh/meshGFaceBDS.cpp
@@ -24,26 +24,28 @@
#include "Field.h"
#include "OS.h"
-
-static void getDegeneratedVertices (BDS_Mesh &m,
- std::map<BDS_Point*, MVertex*,PointLessThan> *recoverMap,
- std::set<MVertex *, MVertexLessThanNum> °enerated,
- std::vector<BDS_Edge*> °enerated_edges){
+static void getDegeneratedVertices(
+ BDS_Mesh &m, std::map<BDS_Point *, MVertex *, PointLessThan> *recoverMap,
+ std::set<MVertex *, MVertexLessThanNum> °enerated,
+ std::vector<BDS_Edge *> °enerated_edges)
+{
degenerated.clear();
- std::vector<BDS_Edge*>::iterator it = m.edges.begin();
+ std::vector<BDS_Edge *>::iterator it = m.edges.begin();
- while (it != m.edges.end()){
+ while(it != m.edges.end()) {
BDS_Edge *e = *it;
- if (!e->deleted && e->numfaces() == 1){
- std::map<BDS_Point*, MVertex*,PointLessThan>::iterator itp1 = recoverMap->find(e->p1);
- std::map<BDS_Point*, MVertex*,PointLessThan>::iterator itp2 = recoverMap->find(e->p2);
- if (itp1 != recoverMap->end() && itp2 != recoverMap->end() &&
- itp1->second == itp2->second){
+ if(!e->deleted && e->numfaces() == 1) {
+ std::map<BDS_Point *, MVertex *, PointLessThan>::iterator itp1 =
+ recoverMap->find(e->p1);
+ std::map<BDS_Point *, MVertex *, PointLessThan>::iterator itp2 =
+ recoverMap->find(e->p2);
+ if(itp1 != recoverMap->end() && itp2 != recoverMap->end() &&
+ itp1->second == itp2->second) {
degenerated.insert(itp1->second);
- // itp1->first->_degeneratedTo = itp2->first;
- // itp2->first->_degeneratedTo = itp1->first;
- degenerated_edges.push_back(e);
+ // itp1->first->_degeneratedTo = itp2->first;
+ // itp2->first->_degeneratedTo = itp1->first;
+ degenerated_edges.push_back(e);
}
}
++it;
@@ -60,11 +62,9 @@ double computeEdgeLinearLength(BDS_Point *p1, BDS_Point *p2)
inline double computeEdgeLinearLength(BDS_Point *p1, BDS_Point *p2, GFace *f)
{
- GPoint GP = f->point(SPoint2(0.5 * (p1->u + p2->u),
- 0.5 * (p1->v + p2->v)));
+ GPoint GP = f->point(SPoint2(0.5 * (p1->u + p2->u), 0.5 * (p1->v + p2->v)));
- if (!GP.succeeded())
- return computeEdgeLinearLength(p1,p2);
+ if(!GP.succeeded()) return computeEdgeLinearLength(p1, p2);
const double dx1 = p1->X - GP.x();
const double dy1 = p1->Y - GP.y();
@@ -93,11 +93,11 @@ double NewGetLc(BDS_Point *point)
point->lcBGM();
}
-static double correctLC_(BDS_Point *p1,BDS_Point *p2, GFace *f)
+static double correctLC_(BDS_Point *p1, BDS_Point *p2, GFace *f)
{
double l1 = NewGetLc(p1);
double l2 = NewGetLc(p2);
- double l = .5*(l1+l2);
+ double l = .5 * (l1 + l2);
const double coord = 0.5;
double U = coord * p1->u + (1 - coord) * p2->u;
@@ -107,13 +107,13 @@ static double correctLC_(BDS_Point *p1,BDS_Point *p2, GFace *f)
double lmid = BGM_MeshSize(f, U, V, gpp.x(), gpp.y(), gpp.z());
l = std::min(l, lmid);
- if(CTX::instance()->mesh.lcFromCurvature){
+ if(CTX::instance()->mesh.lcFromCurvature) {
double l3 = l;
double const lcmin = std::min(std::min(l1, l2), l3);
- l1 = std::min(lcmin*1.2,l1);
- l2 = std::min(lcmin*1.2,l2);
- l3 = std::min(lcmin*1.2,l3);
- l = (l1+l2+l3)/3.0;
+ l1 = std::min(lcmin * 1.2, l1);
+ l2 = std::min(lcmin * 1.2, l2);
+ l3 = std::min(lcmin * 1.2, l3);
+ l = (l1 + l2 + l3) / 3.0;
}
return l;
}
@@ -127,25 +127,25 @@ double NewGetLc(BDS_Edge *const edge, GFace *const face)
double NewGetLc(BDS_Point *p1, BDS_Point *p2, GFace *f)
{
double linearLength = computeEdgeLinearLength(p1, p2, f);
- double l = correctLC_ (p1,p2,f);
+ double l = correctLC_(p1, p2, f);
return linearLength / l;
}
void computeMeshSizeFieldAccuracy(GFace *gf, BDS_Mesh &m, double &avg,
- double &max_e, double &min_e, int &nE, int &GS)
+ double &max_e, double &min_e, int &nE,
+ int &GS)
{
-
- std::vector<BDS_Edge*>::const_iterator it = m.edges.begin();
+ std::vector<BDS_Edge *>::const_iterator it = m.edges.begin();
avg = 0.0;
min_e = 1.e22;
max_e = 0;
nE = 0;
GS = 0;
- while (it != m.edges.end()){
- if (!(*it)->deleted){
+ while(it != m.edges.end()) {
+ if(!(*it)->deleted) {
double const lone = NewGetLc(*it, gf);
- if (lone > 1.0 / std::sqrt(2.0) && lone < std::sqrt(2.0)) GS++;
- avg += lone >1 ? (1. / lone) - 1. : lone - 1.;
+ if(lone > 1.0 / std::sqrt(2.0) && lone < std::sqrt(2.0)) GS++;
+ avg += lone > 1 ? (1. / lone) - 1. : lone - 1.;
max_e = std::max(max_e, lone);
min_e = std::min(min_e, lone);
nE++;
@@ -173,19 +173,19 @@ static bool edgeSwapTestAngle(BDS_Edge *e, double min_cos)
return prosca(norm1, norm2) > min_cos;
}
-
-static bool edgeSwapTestDelaunayAniso(BDS_Edge *e, GFace *gf, std::set<swapquad> &configs)
+static 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->p1->config_modified && !e->p2->config_modified) return false;
if(e->numfaces() != 2) return false;
e->oppositeof(op);
swapquad sq(e->p1->iD, e->p2->iD, op[0]->iD, op[1]->iD);
- if (configs.find(sq) != configs.end()) return false;
+ if(configs.find(sq) != configs.end()) return false;
configs.insert(sq);
double edgeCenter[2] = {0.5 * (e->p1->u + e->p2->u),
@@ -197,7 +197,7 @@ static bool edgeSwapTestDelaunayAniso(BDS_Edge *e, GFace *gf, std::set<swapquad>
double p4[2] = {op[1]->u, op[1]->v};
double metric[3];
buildMetric(gf, edgeCenter, metric);
- if (!inCircumCircleAniso(gf, p1, p2, p3, p4, metric)){
+ if(!inCircumCircleAniso(gf, p1, p2, p3, p4, metric)) {
return false;
}
return true;
@@ -207,9 +207,10 @@ static int edgeSwapTest(GFace *gf, BDS_Edge *e)
{
BDS_Point *op[2];
- const double THRESH = cos(CTX::instance()->mesh.allowSwapEdgeAngle * M_PI / 180.);
+ const double THRESH =
+ cos(CTX::instance()->mesh.allowSwapEdgeAngle * M_PI / 180.);
- e->oppositeof (op);
+ e->oppositeof(op);
double qa1 = qmTriangle::gamma(e->p1, e->p2, op[0]);
double qa2 = qmTriangle::gamma(e->p1, e->p2, op[1]);
@@ -220,28 +221,33 @@ static int edgeSwapTest(GFace *gf, BDS_Edge *e)
// if(qb > 15*qa) return 1;
- if (!edgeSwapTestAngle(e, THRESH))
- return -1;
+ if(!edgeSwapTestAngle(e, THRESH)) return -1;
- if(qb > qa) return 1;
- else return -1;
+ if(qb > qa)
+ return 1;
+ else
+ return -1;
}
-static void swapEdgePass(GFace *gf, BDS_Mesh &m, int &nb_swap, int FINALIZE = 0, double orientation = 1.0)
+static void swapEdgePass(GFace *gf, BDS_Mesh &m, int &nb_swap, int FINALIZE = 0,
+ double orientation = 1.0)
{
BDS_SwapEdgeTest *qual;
- if (FINALIZE && gf->getNativeType() != GEntity::GmshModel) qual = new BDS_SwapEdgeTestNormals(gf, orientation);
- else qual = new BDS_SwapEdgeTestQuality (true,true);
+ if(FINALIZE && gf->getNativeType() != GEntity::GmshModel)
+ qual = new BDS_SwapEdgeTestNormals(gf, orientation);
+ else
+ qual = new BDS_SwapEdgeTestQuality(true, true);
// qual = new BDS_SwapEdgeTestQuality (true,true);
typedef std::vector<BDS_Edge *>::size_type size_type;
size_type origSize = m.edges.size();
- for(size_type index = 0; index < 2*origSize && index < m.edges.size(); ++index) {
- if(!m.edges.at(index)->deleted && m.edges.at(index)->numfaces() == 2) {
+ for(size_type index = 0; index < 2 * origSize && index < m.edges.size();
+ ++index) {
+ if(!m.edges.at(index)->deleted && m.edges.at(index)->numfaces() == 2) {
int const result = FINALIZE ? 1 : edgeSwapTest(gf, m.edges.at(index));
if(result >= 0) {
- if(m.swap_edge(m.edges.at(index), *qual)) {
- ++nb_swap;
- }
+ if(m.swap_edge(m.edges.at(index), *qual)) {
+ ++nb_swap;
+ }
}
}
}
@@ -256,13 +262,14 @@ void delaunayizeBDS(GFace *gf, BDS_Mesh &m, int &nb_swap)
nb_swap = 0;
std::set<swapquad> configs;
while(1) {
-
size_type NSW = 0;
for(size_type index = 0; index < m.edges.size(); ++index) {
if(!m.edges.at(index)->deleted) {
if(edgeSwapTestDelaunayAniso(m.edges.at(index), gf, configs)) {
- if(m.swap_edge(m.edges.at(index), BDS_SwapEdgeTestQuality(false))) { ++NSW; }
+ if(m.swap_edge(m.edges.at(index), BDS_SwapEdgeTestQuality(false))) {
+ ++NSW;
+ }
}
}
}
@@ -271,14 +278,15 @@ void delaunayizeBDS(GFace *gf, BDS_Mesh &m, int &nb_swap)
}
}
-bool edges_sort(std::pair<double, BDS_Edge*> a, std::pair<double, BDS_Edge*> b)
+bool edges_sort(std::pair<double, BDS_Edge *> a,
+ std::pair<double, BDS_Edge *> b)
{
// don't compare pointers: it leads to non-deterministic behavior
// if (a.first == b.first){
// return ((*a.second) < (*b.second));
// }
- if (std::abs(a.first - b.first) < 1e-10){
- if (a.second->p1->iD == b.second->p1->iD)
+ if(std::abs(a.first - b.first) < 1e-10) {
+ if(a.second->p1->iD == b.second->p1->iD)
return (a.second->p2->iD < b.second->p2->iD);
else
return (a.second->p1->iD < b.second->p1->iD);
@@ -287,7 +295,7 @@ bool edges_sort(std::pair<double, BDS_Edge*> a, std::pair<double, BDS_Edge*> b)
return (a.first < b.first);
}
-static bool middlePoint (GFace *gf, BDS_Edge *e, double &u, double &v)
+static bool middlePoint(GFace *gf, BDS_Edge *e, double &u, double &v)
{
// try that
@@ -303,117 +311,126 @@ static bool middlePoint (GFace *gf, BDS_Edge *e, double &u, double &v)
double Z2 = e->p2->Z;
int iter = 0;
- while (1){
- u = 0.5*(u1+u2);
- v = 0.5*(v1+v2);
- GPoint gpp = gf->point(u,v);
+ while(1) {
+ u = 0.5 * (u1 + u2);
+ v = 0.5 * (v1 + v2);
+ GPoint gpp = gf->point(u, v);
double X = gpp.x();
double Y = gpp.y();
double Z = gpp.z();
- double l1 = std::sqrt((X-X1)*(X-X1)+(Y-Y1)*(Y-Y1)+(Z-Z1)*(Z-Z1));
- double l2 = std::sqrt((X-X2)*(X-X2)+(Y-Y2)*(Y-Y2)+(Z-Z2)*(Z-Z2));
+ double l1 = std::sqrt((X - X1) * (X - X1) + (Y - Y1) * (Y - Y1) +
+ (Z - Z1) * (Z - Z1));
+ double l2 = std::sqrt((X - X2) * (X - X2) + (Y - Y2) * (Y - Y2) +
+ (Z - Z2) * (Z - Z2));
// 1 ------ p -- 2
- if (l1 > 1.2*l2){
+ if(l1 > 1.2 * l2) {
// printf("1 %g 2 %g \n",l1,l2);
- u2 = u; v2 = v;
+ u2 = u;
+ v2 = v;
}
- else if (l2 > 1.2*l1){
+ else if(l2 > 1.2 * l1) {
// printf("1 %g 2 %g \n",l1,l2);
- u1 = u; v1 = v;
+ u1 = u;
+ v1 = v;
}
- else break;
- if (iter++ > 10){
- u = 0.5*(e->p1->u+e->p2->u);
- v = 0.5*(e->p1->v+e->p2->v);
+ else
+ break;
+ if(iter++ > 10) {
+ u = 0.5 * (e->p1->u + e->p2->u);
+ v = 0.5 * (e->p1->v + e->p2->v);
return false;
}
}
return true;
}
-
-static void splitEdgePass(GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split, std::vector<SPoint2> *true_boundary)
+static void splitEdgePass(GFace *gf, BDS_Mesh &m, double MAXE_, int &nb_split,
+ std::vector<SPoint2> *true_boundary)
{
- std::vector<std::pair<double, BDS_Edge*> > edges;
+ std::vector<std::pair<double, BDS_Edge *> > edges;
- SPoint2 out (10,10);
+ SPoint2 out(10, 10);
- std::vector<BDS_Edge*>::const_iterator it = m.edges.begin();
- while (it != m.edges.end()){
- if(!(*it)->deleted && (*it)->numfaces() == 2 && (*it)->g->classif_degree == 2){
+ std::vector<BDS_Edge *>::const_iterator it = m.edges.begin();
+ while(it != m.edges.end()) {
+ if(!(*it)->deleted && (*it)->numfaces() == 2 &&
+ (*it)->g->classif_degree == 2) {
double lone = NewGetLc(*it, gf);
- if(lone > MAXE_)edges.push_back(std::make_pair(-lone, *it));
+ if(lone > MAXE_) edges.push_back(std::make_pair(-lone, *it));
}
++it;
}
std::sort(edges.begin(), edges.end(), edges_sort);
- std::vector<BDS_Point*> mids(edges.size());
+ std::vector<BDS_Point *> mids(edges.size());
bool faceDiscrete = gf->geomType() == GEntity::DiscreteSurface;
- for (unsigned int i = 0; i < edges.size(); ++i){
+ for(unsigned int i = 0; i < edges.size(); ++i) {
BDS_Edge *e = edges[i].second;
BDS_Point *mid = NULL;
- if (!e->deleted){
+ if(!e->deleted) {
double U1 = e->p1->u;
double U2 = e->p2->u;
- if (e->p1->degenerated)U1 = U2;
- if (e->p2->degenerated)U2 = U1;
- double U = 0.5*(U1+U2);
- double V = 0.5*(e->p1->v+e->p2->v);
- if (faceDiscrete)middlePoint (gf, e, U, V);
+ if(e->p1->degenerated) U1 = U2;
+ if(e->p2->degenerated) U2 = U1;
+ double U = 0.5 * (U1 + U2);
+ double V = 0.5 * (e->p1->v + e->p2->v);
+ if(faceDiscrete) middlePoint(gf, e, U, V);
- GPoint gpp = gf->point(U,V);
+ GPoint gpp = gf->point(U, V);
bool inside = true;
- if (true_boundary){
- SPoint2 pp(U,V);
- int N;
- if (!pointInsideParametricDomain (*true_boundary, pp, out, N)){
- inside = false;
- // FILE *f = fopen("TOTO.pos","a");
- // fprintf(f,"SP(%g,%g,0){%d};\n",pp.x(),pp.y(),N);
- // fclose(f);
- }
+ if(true_boundary) {
+ SPoint2 pp(U, V);
+ int N;
+ if(!pointInsideParametricDomain(*true_boundary, pp, out, N)) {
+ inside = false;
+ // FILE *f = fopen("TOTO.pos","a");
+ // fprintf(f,"SP(%g,%g,0){%d};\n",pp.x(),pp.y(),N);
+ // fclose(f);
+ }
}
- if (inside && gpp.succeeded()){
- mid = m.add_point(++m.MAXPOINTNUMBER, gpp.x(),gpp.y(),gpp.z());
- mid->u = U;
- mid->v = V;
- mid->lc() = 0.5 * (e->p1->lc() + e->p2->lc());
- mid->lcBGM() = BGM_MeshSize (gf,U,V, mid->X,mid->Y,mid->Z);
+ if(inside && gpp.succeeded()) {
+ mid = m.add_point(++m.MAXPOINTNUMBER, gpp.x(), gpp.y(), gpp.z());
+ mid->u = U;
+ mid->v = V;
+ mid->lc() = 0.5 * (e->p1->lc() + e->p2->lc());
+ mid->lcBGM() = BGM_MeshSize(gf, U, V, mid->X, mid->Y, mid->Z);
}
}
mids[i] = mid;
}
- for (unsigned int i = 0; i < edges.size(); ++i){
+ for(unsigned int i = 0; i < edges.size(); ++i) {
BDS_Edge *e = edges[i].second;
- if (!e->deleted){
+ if(!e->deleted) {
BDS_Point *mid = mids[i];
- if (mid){
- if(!m.split_edge(e, mid)) m.del_point(mid);
- else nb_split++;
+ if(mid) {
+ if(!m.split_edge(e, mid))
+ m.del_point(mid);
+ else
+ nb_split++;
}
}
}
}
-double getMaxLcWhenCollapsingEdge(GFace *gf, BDS_Mesh &m, BDS_Edge *e, BDS_Point *p)
+double getMaxLcWhenCollapsingEdge(GFace *gf, BDS_Mesh &m, BDS_Edge *e,
+ BDS_Point *p)
{
BDS_Point *o = e->othervertex(p);
double maxLc = 0.0;
- std::vector<BDS_Edge*> edges(p->edges);
- std::vector<BDS_Edge*>::iterator eit = edges.begin();
- while (eit != edges.end()) {
+ std::vector<BDS_Edge *> edges(p->edges);
+ std::vector<BDS_Edge *>::iterator eit = edges.begin();
+ while(eit != edges.end()) {
BDS_Point *newP1 = 0, *newP2 = 0;
- if ((*eit)->p1 == p){
+ if((*eit)->p1 == p) {
newP1 = o;
newP2 = (*eit)->p2;
}
- else if ((*eit)->p2 == p){
+ else if((*eit)->p2 == p) {
newP1 = (*eit)->p1;
newP2 = o;
}
@@ -455,54 +472,56 @@ double getMaxLcWhenCollapsingEdge(GFace *gf, BDS_Mesh &m, BDS_Edge *e, BDS_Point
return false;
}
*/
-void collapseEdgePass(GFace *gf, BDS_Mesh &m, double MINE_, int MAXNP, int &nb_collaps)
+void collapseEdgePass(GFace *gf, BDS_Mesh &m, double MINE_, int MAXNP,
+ int &nb_collaps)
{
// return;
- std::vector<BDS_Edge*>::const_iterator it = m.edges.begin();
- std::vector<std::pair<double, BDS_Edge*> > edges;
+ std::vector<BDS_Edge *>::const_iterator it = m.edges.begin();
+ std::vector<std::pair<double, BDS_Edge *> > edges;
- while (it != m.edges.end()){
- if(!(*it)->deleted && (*it)->numfaces() == 2 && (*it)->g->classif_degree == 2){
+ while(it != m.edges.end()) {
+ if(!(*it)->deleted && (*it)->numfaces() == 2 &&
+ (*it)->g->classif_degree == 2) {
double lone = NewGetLc(*it, gf);
- if(lone < MINE_) edges.push_back (std::make_pair(lone, *it));
+ if(lone < MINE_) edges.push_back(std::make_pair(lone, *it));
}
++it;
}
std::sort(edges.begin(), edges.end(), edges_sort);
- for (unsigned int i = 0; i < edges.size(); i++){
+ for(unsigned int i = 0; i < edges.size(); i++) {
BDS_Edge *e = edges[i].second;
- if(!e->deleted){
+ if(!e->deleted) {
double lone1 = 0.;
bool collapseP1Allowed = false;
// lone1 = getMaxLcWhenCollapsingEdge(gf, m, e, e->p1);
- collapseP1Allowed = true;//std::abs(lone1-1.0) < std::abs(edges[i].first - 1.0);
+ collapseP1Allowed =
+ true; // std::abs(lone1-1.0) < std::abs(edges[i].first - 1.0);
double lone2 = 0.;
bool collapseP2Allowed = false;
- if (e->p2->iD > MAXNP){
- // lone2 = getMaxLcWhenCollapsingEdge(gf, m, e, e->p2);
- collapseP2Allowed = true;//std::abs(lone2-1.0) < std::abs(edges[i].first - 1.0);
+ if(e->p2->iD > MAXNP) {
+ // lone2 = getMaxLcWhenCollapsingEdge(gf, m, e, e->p2);
+ collapseP2Allowed =
+ true; // std::abs(lone2-1.0) < std::abs(edges[i].first - 1.0);
}
BDS_Point *p = 0;
- if (collapseP1Allowed && collapseP2Allowed){
- if (std::abs(lone1 - lone2) < 1e-12)
+ if(collapseP1Allowed && collapseP2Allowed) {
+ if(std::abs(lone1 - lone2) < 1e-12)
p = e->p1->iD < e->p2->iD ? e->p1 : e->p2;
else
p = std::abs(lone1 - 1.0) < std::abs(lone2 - 1.0) ? e->p1 : e->p2;
}
- else if (collapseP1Allowed && !collapseP2Allowed)
+ else if(collapseP1Allowed && !collapseP2Allowed)
p = e->p1;
- else if (collapseP1Allowed && !collapseP2Allowed)
+ else if(collapseP1Allowed && !collapseP2Allowed)
p = e->p2;
bool res = false;
- if(p)
- res = m.collapse_edge_parametric(e, p);
- if(res)
- nb_collaps++;
+ if(p) res = m.collapse_edge_parametric(e, p);
+ if(res) nb_collaps++;
}
}
}
@@ -510,15 +529,15 @@ void collapseEdgePass(GFace *gf, BDS_Mesh &m, double MINE_, int MAXNP, int &nb_c
void smoothVertexPass(GFace *gf, BDS_Mesh &m, int &nb_smooth, bool q)
{
// FIXME SUPER HACK
- std::set<BDS_Point*,PointLessThan>::iterator itp = m.points.begin();
- while(itp != m.points.end()){
- if(m.smooth_point_centroid(*itp, gf,q))
- nb_smooth ++;
+ std::set<BDS_Point *, PointLessThan>::iterator itp = m.points.begin();
+ while(itp != m.points.end()) {
+ if(m.smooth_point_centroid(*itp, gf, q)) nb_smooth++;
++itp;
}
}
-void CHECK_STRANGE (const char *c,BDS_Mesh &m){
+void CHECK_STRANGE(const char *c, BDS_Mesh &m)
+{
return;
#if 0
int strange = 0;
@@ -559,115 +578,124 @@ void CHECK_STRANGE (const char *c,BDS_Mesh &m){
#endif
}
-void computeNodalSizes(GFace *gf, BDS_Mesh &m,
- std::map<BDS_Point*, MVertex*,PointLessThan> *recoverMap){
-
- std::set<BDS_Point*,PointLessThan>::iterator itp = m.points.begin();
- while (itp != m.points.end()){
- std::vector<BDS_Edge*>::iterator it = (*itp)->edges.begin();
- std::vector<BDS_Edge*>::iterator ite = (*itp)->edges.end();
+void computeNodalSizes(
+ GFace *gf, BDS_Mesh &m,
+ std::map<BDS_Point *, MVertex *, PointLessThan> *recoverMap)
+{
+ std::set<BDS_Point *, PointLessThan>::iterator itp = m.points.begin();
+ while(itp != m.points.end()) {
+ std::vector<BDS_Edge *>::iterator it = (*itp)->edges.begin();
+ std::vector<BDS_Edge *>::iterator ite = (*itp)->edges.end();
double L = 0;
int ne = 0;
- while(it != ite){
+ while(it != ite) {
double l = (*it)->length();
- if ((*it)->g && (*it)->g->classif_degree == 1){
- L = ne ? std::max(L, l) : l;
- ne++;
- }
+ if((*it)->g && (*it)->g->classif_degree == 1) {
+ L = ne ? std::max(L, l) : l;
+ ne++;
+ }
++it;
}
- if ((*itp)->g && (*itp)->g->classif_tag > 0){
- if (!ne) L = MAX_LC;
- if(CTX::instance()->mesh.lcFromPoints)
- (*itp)->lc() = L;
+ if((*itp)->g && (*itp)->g->classif_tag > 0) {
+ if(!ne) L = MAX_LC;
+ if(CTX::instance()->mesh.lcFromPoints) (*itp)->lc() = L;
(*itp)->lcBGM() = L;
}
++itp;
}
return;
{
- std::set<BDS_Point*,PointLessThan>::iterator itp = m.points.begin();
- while (itp != m.points.end()){
+ std::set<BDS_Point *, PointLessThan>::iterator itp = m.points.begin();
+ while(itp != m.points.end()) {
(*itp)->lc() = MAX_LC;
++itp;
}
- std::vector<BDS_Edge*>::const_iterator it = m.edges.begin();
- while (it != m.edges.end()){
+ std::vector<BDS_Edge *>::const_iterator it = m.edges.begin();
+ while(it != m.edges.end()) {
bool degenerated = false;
- if (recoverMap){
- std::map<BDS_Point*, MVertex*,PointLessThan>::iterator itp1 = recoverMap->find((*it)->p1);
- std::map<BDS_Point*, MVertex*,PointLessThan>::iterator itp2 = recoverMap->find((*it)->p2);
- if (itp1 != recoverMap->end() && itp2 != recoverMap->end() &&
- itp1->second == itp2->second){
- degenerated = true;
- }
+ if(recoverMap) {
+ std::map<BDS_Point *, MVertex *, PointLessThan>::iterator itp1 =
+ recoverMap->find((*it)->p1);
+ std::map<BDS_Point *, MVertex *, PointLessThan>::iterator itp2 =
+ recoverMap->find((*it)->p2);
+ if(itp1 != recoverMap->end() && itp2 != recoverMap->end() &&
+ itp1->second == itp2->second) {
+ degenerated = true;
+ }
}
- if (!degenerated && (*it)->g && (*it)->g->classif_degree == 1){
- double L = (*it)->length();
- if ((*it)->p1->lc() == MAX_LC)(*it)->p1->lc() = L;
- else (*it)->p1->lc() = std::max (L, (*it)->p1->lc());
- if ((*it)->p2->lc() == MAX_LC)(*it)->p2->lc() = L;
- else (*it)->p2->lc() = std::max (L, (*it)->p2->lc());
+ if(!degenerated && (*it)->g && (*it)->g->classif_degree == 1) {
+ double L = (*it)->length();
+ if((*it)->p1->lc() == MAX_LC)
+ (*it)->p1->lc() = L;
+ else
+ (*it)->p1->lc() = std::max(L, (*it)->p1->lc());
+ if((*it)->p2->lc() == MAX_LC)
+ (*it)->p2->lc() = L;
+ else
+ (*it)->p2->lc() = std::max(L, (*it)->p2->lc());
}
++it;
}
itp = m.points.begin();
- std::vector<BDS_Point*> pts;
- while (itp != m.points.end()){
- if ((*itp)->lc() == MAX_LC)pts.push_back(*itp);
+ std::vector<BDS_Point *> pts;
+ while(itp != m.points.end()) {
+ if((*itp)->lc() == MAX_LC) pts.push_back(*itp);
++itp;
}
int iter = 0;
- while (1){
+ while(1) {
bool allTouched = true;
- for (size_t i=0; i< pts.size(); i++){
- std::vector<BDS_Edge*>::const_iterator it = pts[i]->edges.begin();
- std::vector<BDS_Edge*>::const_iterator ite = pts[i]->edges.end();
- while(it != ite){
- BDS_Point *p = (*it)->othervertex (pts[i]);
- if (p->lc() != MAX_LC){
- if (pts[i]->lc() == MAX_LC)pts[i]->lc() = p->lc();
- else pts[i]->lc() = 0.5*(p->lc()+pts[i]->lc());
- }
- ++it;
- }
- if (pts[i]->lc() == MAX_LC) allTouched = false;
+ for(size_t i = 0; i < pts.size(); i++) {
+ std::vector<BDS_Edge *>::const_iterator it = pts[i]->edges.begin();
+ std::vector<BDS_Edge *>::const_iterator ite = pts[i]->edges.end();
+ while(it != ite) {
+ BDS_Point *p = (*it)->othervertex(pts[i]);
+ if(p->lc() != MAX_LC) {
+ if(pts[i]->lc() == MAX_LC)
+ pts[i]->lc() = p->lc();
+ else
+ pts[i]->lc() = 0.5 * (p->lc() + pts[i]->lc());
+ }
+ ++it;
+ }
+ if(pts[i]->lc() == MAX_LC) allTouched = false;
}
- if (iter++ >= 5 || allTouched)break;
+ if(iter++ >= 5 || allTouched) break;
}
}
}
-void modifyInitialMeshToRemoveDegeneracies(GFace *gf, BDS_Mesh &m,
- std::map<BDS_Point*, MVertex*,PointLessThan> *recoverMap)
+void modifyInitialMeshToRemoveDegeneracies(
+ GFace *gf, BDS_Mesh &m,
+ std::map<BDS_Point *, MVertex *, PointLessThan> *recoverMap)
{
std::set<MVertex *, MVertexLessThanNum> degenerated;
- std::vector<BDS_Edge*> degenerated_edges;
- getDegeneratedVertices (m,recoverMap,degenerated,degenerated_edges);
+ std::vector<BDS_Edge *> degenerated_edges;
+ getDegeneratedVertices(m, recoverMap, degenerated, degenerated_edges);
// printf("%d degenerated vertices\n",degenerated.size());
- for (std::map<BDS_Point*, MVertex*,PointLessThan>::iterator it = recoverMap->begin();
- it != recoverMap->end(); ++it){
- std::set<MVertex *, MVertexLessThanNum>::iterator it2 = degenerated.find(it->second);
- if (it2 != degenerated.end()){
+ for(std::map<BDS_Point *, MVertex *, PointLessThan>::iterator it =
+ recoverMap->begin();
+ it != recoverMap->end(); ++it) {
+ std::set<MVertex *, MVertexLessThanNum>::iterator it2 =
+ degenerated.find(it->second);
+ if(it2 != degenerated.end()) {
it->first->degenerated = true;
-
}
}
- for (size_t i = 0;i<degenerated_edges.size();i++){
- m.collapse_edge_parametric(degenerated_edges[i], degenerated_edges[i]->p1, true);
+ for(size_t i = 0; i < degenerated_edges.size(); i++) {
+ m.collapse_edge_parametric(degenerated_edges[i], degenerated_edges[i]->p1,
+ true);
recoverMap->erase(degenerated_edges[i]->p1);
}
}
-
-
void refineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
const bool computeNodalSizeField,
- std::map<MVertex*, BDS_Point*> *recoverMapInv,
- std::map<BDS_Point*, MVertex*,PointLessThan> *recoverMap,
- std::vector<SPoint2> *true_boundary)
+ std::map<MVertex *, BDS_Point *> *recoverMapInv,
+ std::map<BDS_Point *, MVertex *, PointLessThan> *recoverMap,
+ std::vector<SPoint2> *true_boundary)
{
// return;
int IT = 0;
@@ -675,13 +703,13 @@ void refineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
// classify correctly the embedded vertices use a negative model
// face number to avoid mesh motion
- if(recoverMapInv){
- std::list<GVertex*> emb_vertx = gf->embeddedVertices();
- std::list<GVertex*>::iterator itvx = emb_vertx.begin();
- while(itvx != emb_vertx.end()){
+ if(recoverMapInv) {
+ std::list<GVertex *> emb_vertx = gf->embeddedVertices();
+ std::list<GVertex *>::iterator itvx = emb_vertx.begin();
+ while(itvx != emb_vertx.end()) {
MVertex *v = *((*itvx)->mesh_vertices.begin());
- std::map<MVertex*, BDS_Point*>::iterator itp = recoverMapInv->find(v);
- if(itp != recoverMapInv->end()){
+ std::map<MVertex *, BDS_Point *>::iterator itp = recoverMapInv->find(v);
+ if(itp != recoverMapInv->end()) {
BDS_Point *p = itp->second;
m.add_geom(-1, 2);
p->g = m.get_geom(-1, 2);
@@ -692,16 +720,15 @@ void refineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
// if (recoverMap)outputScalarField(m.triangles, "init.pos", 1, gf);
-
// If asked, compute nodal size field using 1D Mesh
- if (computeNodalSizeField)computeNodalSizes(gf, m,recoverMap);
+ if(computeNodalSizeField) computeNodalSizes(gf, m, recoverMap);
- double t_spl = 0, t_sw = 0,t_col = 0,t_sm = 0;
+ double t_spl = 0, t_sw = 0, t_col = 0, t_sm = 0;
const double MINE_ = 0.7, MAXE_ = 1.4;
- while (1){
- CHECK_STRANGE ("start",m);
+ while(1) {
+ CHECK_STRANGE("start", m);
// we count the number of local mesh modifs.
int nb_split = 0;
@@ -713,11 +740,11 @@ void refineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
double minL = 1.e22, maxL = 0;
int NN1 = m.edges.size();
int NN2 = 0;
- std::vector<BDS_Edge*>::iterator it = m.edges.begin();
+ std::vector<BDS_Edge *>::iterator it = m.edges.begin();
- while (1){
- if (NN2++ >= NN1)break;
- if (!(*it)->deleted){
+ while(1) {
+ if(NN2++ >= NN1) break;
+ if(!(*it)->deleted) {
(*it)->p1->config_modified = false;
(*it)->p2->config_modified = false;
double lone = NewGetLc(*it, gf);
@@ -727,33 +754,33 @@ void refineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
++it;
}
- if ((minL > MINE_ && maxL < MAXE_) || IT > (abs(NIT))) break;
+ if((minL > MINE_ && maxL < MAXE_) || IT > (abs(NIT))) break;
double maxE = MAXE_;
double minE = MINE_;
double t1 = Cpu();
// outputScalarField(m.triangles, "0.pos", 1, gf);
splitEdgePass(gf, m, maxE, nb_split, true_boundary);
// outputScalarField(m.triangles, "a.pos", 1, gf);
- CHECK_STRANGE ("split",m);
+ CHECK_STRANGE("split", m);
double t2 = Cpu();
swapEdgePass(gf, m, nb_swap);
- CHECK_STRANGE ("swap",m);
+ CHECK_STRANGE("swap", m);
// outputScalarField(m.triangles, "b.pos", 1, gf);
double t3 = Cpu();
collapseEdgePass(gf, m, minE, MAXNP, nb_collaps);
- CHECK_STRANGE ("collapse",m);
+ CHECK_STRANGE("collapse", m);
// outputScalarField(m.triangles, "c.pos", 1, gf);
double t4 = Cpu();
- swapEdgePass ( gf, m, nb_swap);
- CHECK_STRANGE ("swap",m);
+ swapEdgePass(gf, m, nb_swap);
+ CHECK_STRANGE("swap", m);
double t5 = Cpu();
smoothVertexPass(gf, m, nb_smooth, false);
double t6 = Cpu();
- CHECK_STRANGE ("smmooth",m);
- swapEdgePass ( gf, m, nb_swap);
- CHECK_STRANGE ("swap",m);
+ CHECK_STRANGE("smmooth", m);
+ swapEdgePass(gf, m, nb_swap);
+ CHECK_STRANGE("swap", m);
double t7 = Cpu();
// char name[256]; sprintf(name,"iter%d.pos",IT);
// outputScalarField(m.triangles, name, 1, gf);
@@ -761,11 +788,11 @@ void refineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
// }
// clean up the mesh
t_spl += t2 - t1;
- t_sw += t3 - t2;
- t_sw += t5 - t4;
- t_sw += t7 - t6;
+ t_sw += t3 - t2;
+ t_sw += t5 - t4;
+ t_sw += t7 - t6;
t_col += t4 - t3;
- t_sm += t6 - t5;
+ t_sm += t6 - t5;
m.cleanup();
// char nn[256];
// sprintf(nn,"ITER%d.pos",IT);
@@ -775,10 +802,10 @@ void refineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
IT++;
Msg::Debug(" iter %3d minL %8.3f/%8.3f maxL %8.3f/%8.3f : "
"%6d splits, %6d swaps, %6d collapses, %6d moves",
- IT, minL, minE, maxL, maxE, nb_split, nb_swap, nb_collaps, nb_smooth);
+ IT, minL, minE, maxL, maxE, nb_split, nb_swap, nb_collaps,
+ nb_smooth);
// getchar();
- if (nb_split == 0 && nb_collaps == 0) break;
-
+ if(nb_split == 0 && nb_collaps == 0) break;
}
// outputScalarField(m.triangles, "before.pos", 1, gf);
@@ -786,56 +813,60 @@ void refineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
int bad = 0;
int invalid = 0;
- if (gf->getNativeType() != GEntity::GmshModel){
- for (size_t i=0;i<m.triangles.size();i++){
+ if(gf->getNativeType() != GEntity::GmshModel) {
+ for(size_t i = 0; i < m.triangles.size(); i++) {
BDS_Point *pts[4];
m.triangles[i]->getNodes(pts);
- double val = BDS_Face_Validity (gf, m.triangles[i]);
+ double val = BDS_Face_Validity(gf, m.triangles[i]);
invalid += val < 0 ? 1 : 0;
}
- double orientation = invalid > (int)m.triangles.size()/2 ? -1.0 : 1.0;
+ double orientation = invalid > (int)m.triangles.size() / 2 ? -1.0 : 1.0;
- while (ITER++ <10){
+ while(ITER++ < 10) {
// printf("START\n");
bad = 0;
invalid = 0;
- for (size_t i=0;i<m.triangles.size();i++){
- if (!m.triangles[i]->deleted){
- BDS_Point *pts[4];
- m.triangles[i]->getNodes(pts);
- pts[0]->config_modified = false;
- pts[1]->config_modified = false;
- pts[2]->config_modified = false;
- }
+ for(size_t i = 0; i < m.triangles.size(); i++) {
+ if(!m.triangles[i]->deleted) {
+ BDS_Point *pts[4];
+ m.triangles[i]->getNodes(pts);
+ pts[0]->config_modified = false;
+ pts[1]->config_modified = false;
+ pts[2]->config_modified = false;
+ }
}
- for (size_t i=0;i<m.triangles.size();i++){
- BDS_Point *pts[4];
- m.triangles[i]->getNodes(pts);
- // if (pts[0]->degenerated + pts[1]->degenerated + pts[2]->degenerated < 2){
- double val = orientation * BDS_Face_Validity (gf, m.triangles[i]);
- if (val <= 0.2){
- if (!m.triangles[i]->deleted && val <= 0)invalid++;
- pts[0]->config_modified = true;
- pts[1]->config_modified = true;
- pts[2]->config_modified = true;
- bad ++;
- if (val < 0) invalid++;
- }
- // }
+ for(size_t i = 0; i < m.triangles.size(); i++) {
+ BDS_Point *pts[4];
+ m.triangles[i]->getNodes(pts);
+ // if (pts[0]->degenerated + pts[1]->degenerated +
+ // pts[2]->degenerated < 2){
+ double val = orientation * BDS_Face_Validity(gf, m.triangles[i]);
+ if(val <= 0.2) {
+ if(!m.triangles[i]->deleted && val <= 0) invalid++;
+ pts[0]->config_modified = true;
+ pts[1]->config_modified = true;
+ pts[2]->config_modified = true;
+ bad++;
+ if(val < 0) invalid++;
+ }
+ // }
}
- if (bad != 0){
- int nb_swap = 0;
- int nb_smooth = 0;
- swapEdgePass ( gf, m, nb_swap, 1, orientation);
- smoothVertexPass(gf, m, nb_smooth, true);
- if ((nb_swap == 0 && nb_smooth == 0) || ITER == 10){
- if (invalid && !computeNodalSizeField)Msg::Warning("Meshing surface %d : %d elements remain invalid\n", gf->tag(),invalid);
- break;
- }
+ if(bad != 0) {
+ int nb_swap = 0;
+ int nb_smooth = 0;
+ swapEdgePass(gf, m, nb_swap, 1, orientation);
+ smoothVertexPass(gf, m, nb_smooth, true);
+ if((nb_swap == 0 && nb_smooth == 0) || ITER == 10) {
+ if(invalid && !computeNodalSizeField)
+ Msg::Warning("Meshing surface %d : %d elements remain invalid\n",
+ gf->tag(), invalid);
+ break;
+ }
}
else {
- // Msg::Info("Meshing surface %d : all elements are oriented properly\n", gf->tag());
- break;
+ // Msg::Info("Meshing surface %d : all elements are oriented
+ // properly\n", gf->tag());
+ break;
}
}
}
@@ -863,146 +894,158 @@ void refineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
Msg::Debug(" ---------------------------------------");
Msg::Debug(" CPU Report ");
Msg::Debug(" ---------------------------------------");
- Msg::Debug(" CPU SWAP %12.5E sec (%3.0f %%)", t_sw,100 * t_sw / t_total);
- Msg::Debug(" CPU SPLIT %12.5E sec (%3.0f %%) ", t_spl,100 * t_spl / t_total);
- Msg::Debug(" CPU COLLAPS %12.5E sec (%3.0f %%) ", t_col,100 * t_col / t_total);
- Msg::Debug(" CPU SMOOTH %12.5E sec (%3.0f %%) ", t_sm,100 * t_sm / t_total);
+ Msg::Debug(" CPU SWAP %12.5E sec (%3.0f %%)", t_sw, 100 * t_sw / t_total);
+ Msg::Debug(" CPU SPLIT %12.5E sec (%3.0f %%) ", t_spl,
+ 100 * t_spl / t_total);
+ Msg::Debug(" CPU COLLAPS %12.5E sec (%3.0f %%) ", t_col,
+ 100 * t_col / t_total);
+ Msg::Debug(" CPU SMOOTH %12.5E sec (%3.0f %%) ", t_sm, 100 * t_sm / t_total);
Msg::Debug(" ---------------------------------------");
- Msg::Debug(" CPU TOTAL %12.5E sec ",t_total);
+ Msg::Debug(" CPU TOTAL %12.5E sec ", t_total);
Msg::Debug(" ---------------------------------------");
}
-bool degeneratedTriangle (BDS_Face *f, std::map<BDS_Point*, MVertex*,PointLessThan> *recoverMap){
-
- MVertex *v[3] = {0,0,0};
+bool degeneratedTriangle(
+ BDS_Face *f, std::map<BDS_Point *, MVertex *, PointLessThan> *recoverMap)
+{
+ MVertex *v[3] = {0, 0, 0};
BDS_Point *n[4];
f->getNodes(n);
- for (int i=0;i<3;i++){
- std::map<BDS_Point*, MVertex*,PointLessThan>::iterator itp1 = recoverMap->find(n[i]);
- if (itp1 != recoverMap->end())v[i] = itp1->second;
+ for(int i = 0; i < 3; i++) {
+ std::map<BDS_Point *, MVertex *, PointLessThan>::iterator itp1 =
+ recoverMap->find(n[i]);
+ if(itp1 != recoverMap->end()) v[i] = itp1->second;
}
- if (v[0] && (v[0] == v[1] || v[0] == v[2]))return true;
- if (v[1] && (v[1] == v[2] || v[0] == v[1]))return true;
- if (v[2] && (v[1] == v[2] || v[0] == v[2]))return true;
+ if(v[0] && (v[0] == v[1] || v[0] == v[2])) return true;
+ if(v[1] && (v[1] == v[2] || v[0] == v[1])) return true;
+ if(v[2] && (v[1] == v[2] || v[0] == v[2])) return true;
return false;
}
-
-void invalidEdgesPeriodic(BDS_Mesh &m,
- std::map<BDS_Point*, MVertex*,PointLessThan> *recoverMap,
- std::set<BDS_Edge*, EdgeLessThan> &toSplit)
+void invalidEdgesPeriodic(
+ BDS_Mesh &m, std::map<BDS_Point *, MVertex *, PointLessThan> *recoverMap,
+ std::set<BDS_Edge *, EdgeLessThan> &toSplit)
{
-
// first look for degenerated vertices
std::set<MVertex *, MVertexLessThanNum> degenerated;
- std::vector<BDS_Edge*> degenerated_edges;
- getDegeneratedVertices (m,recoverMap,degenerated,degenerated_edges);
- Msg::Debug("Found %d degenerated vertices",degenerated.size());
+ std::vector<BDS_Edge *> degenerated_edges;
+ getDegeneratedVertices(m, recoverMap, degenerated, degenerated_edges);
+ Msg::Debug("Found %d degenerated vertices", degenerated.size());
toSplit.clear();
- std::vector<BDS_Edge*>::const_iterator it = m.edges.begin();
- std::map<std::pair<MVertex*, BDS_Point*>, BDS_Edge *> touchPeriodic;
- while (it != m.edges.end()){
+ std::vector<BDS_Edge *>::const_iterator it = m.edges.begin();
+ std::map<std::pair<MVertex *, BDS_Point *>, BDS_Edge *> touchPeriodic;
+ while(it != m.edges.end()) {
BDS_Edge *e = *it;
- if (!e->deleted && e->numfaces() == 2){
- std::map<BDS_Point*, MVertex*,PointLessThan>::iterator itp1 = recoverMap->find(e->p1);
- std::map<BDS_Point*, MVertex*,PointLessThan>::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 =
+ if(!e->deleted && e->numfaces() == 2) {
+ std::map<BDS_Point *, MVertex *, PointLessThan>::iterator itp1 =
+ recoverMap->find(e->p1);
+ std::map<BDS_Point *, MVertex *, PointLessThan>::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 =
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);
+ 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 =
+ 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 =
touchPeriodic.find(a);
- if (itf == touchPeriodic.end()) touchPeriodic[a] = e;
- else if (degenerated.find(itp2->second) == degenerated.end()){
- toSplit.insert(e); toSplit.insert(itf->second);
+ if(itf == touchPeriodic.end())
+ touchPeriodic[a] = e;
+ else if(degenerated.find(itp2->second) == degenerated.end()) {
+ toSplit.insert(e);
+ toSplit.insert(itf->second);
}
}
}
++it;
}
- //return;
- Msg::Debug("%d degenerated edges should be split",toSplit.size());
+ // return;
+ Msg::Debug("%d degenerated edges should be split", toSplit.size());
size_t aa = toSplit.size();
// same edges should be split
{
- std::set<BDS_Point*, PointLessThan>::iterator it = m.points.begin();
- while (it != m.points.end()){
- std::map<BDS_Point*, MVertex*,PointLessThan>::iterator itp2 = recoverMap->find(*it);
- if (itp2 == recoverMap->end()){
- std::vector<BDS_Edge*>::const_iterator ite = (*it)->edges.begin();
- std::map<MVertex*,BDS_Edge*> mp;
- while (ite != (*it)->edges.end()){
- if (!(*ite)->deleted && (*ite)->numfaces() == 2){
- if (!degeneratedTriangle ((*ite)->faces(0),recoverMap) &&
- !degeneratedTriangle ((*ite)->faces(1),recoverMap)){
- BDS_Point *other = (*it) == (*ite)->p1 ? (*ite)->p2 : (*ite)->p1;
- std::map<BDS_Point*, MVertex*,PointLessThan>::iterator itp1 = recoverMap->find(other);
- if (itp1 != recoverMap->end()){
- if (mp.find(itp1->second) != mp.end()){
- toSplit.insert(*ite); toSplit.insert(mp[itp1->second]);
- }
- mp[itp1->second] = *ite;
- }
- }
- }
- ++ite;
- }
+ std::set<BDS_Point *, PointLessThan>::iterator it = m.points.begin();
+ while(it != m.points.end()) {
+ std::map<BDS_Point *, MVertex *, PointLessThan>::iterator itp2 =
+ recoverMap->find(*it);
+ if(itp2 == recoverMap->end()) {
+ std::vector<BDS_Edge *>::const_iterator ite = (*it)->edges.begin();
+ std::map<MVertex *, BDS_Edge *> mp;
+ while(ite != (*it)->edges.end()) {
+ if(!(*ite)->deleted && (*ite)->numfaces() == 2) {
+ if(!degeneratedTriangle((*ite)->faces(0), recoverMap) &&
+ !degeneratedTriangle((*ite)->faces(1), recoverMap)) {
+ BDS_Point *other = (*it) == (*ite)->p1 ? (*ite)->p2 : (*ite)->p1;
+ std::map<BDS_Point *, MVertex *, PointLessThan>::iterator itp1 =
+ recoverMap->find(other);
+ if(itp1 != recoverMap->end()) {
+ if(mp.find(itp1->second) != mp.end()) {
+ toSplit.insert(*ite);
+ toSplit.insert(mp[itp1->second]);
+ }
+ mp[itp1->second] = *ite;
+ }
+ }
+ }
+ ++ite;
+ }
}
++it;
}
}
- Msg::Debug("%d similar edges should be split",toSplit.size()-aa);
+ Msg::Debug("%d similar edges should be split", toSplit.size() - aa);
// printf("%d splits\n",toSplit.size());
-
-
}
// consider p1 and p2, 2 vertices that are different in the parametric
// plane BUT are the same in the real plane
// consider a vertex v that is internal
-// if p1 is the start and the end of a degenerated edge, then allow edges p1 v and p2 v
-// if not do not allow p1 v and p2 v, split them
-// if p1 p2 exists and it is internal, split it
+// if p1 is the start and the end of a degenerated edge, then allow edges p1 v
+// and p2 v 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,
- std::map<BDS_Point*, MVertex*,PointLessThan> *recoverMap)
+int solveInvalidPeriodic(
+ GFace *gf, BDS_Mesh &m,
+ std::map<BDS_Point *, MVertex *, PointLessThan> *recoverMap)
{
- if (!recoverMap) return 0;
+ if(!recoverMap) return 0;
- BDS_GeomEntity *g = m.get_geom(gf->tag(),2);
+ BDS_GeomEntity *g = m.get_geom(gf->tag(), 2);
// return 0;
- Msg::Debug("Solving invalid edges/triangles in initial meshes of periodic face %d",gf->tag());
+ Msg::Debug(
+ "Solving invalid edges/triangles in initial meshes of periodic face %d",
+ gf->tag());
- std::set<BDS_Edge*, EdgeLessThan> toSplit;
+ std::set<BDS_Edge *, EdgeLessThan> toSplit;
invalidEdgesPeriodic(m, recoverMap, toSplit);
- std::set<BDS_Edge*>::iterator ite = toSplit.begin();
+ std::set<BDS_Edge *>::iterator ite = toSplit.begin();
- Msg::Debug("%d edges of the initial mesh are going to be split",toSplit.size());
+ Msg::Debug("%d edges of the initial mesh are going to be split",
+ toSplit.size());
- for (;ite !=toSplit.end(); ++ite){
+ for(; ite != toSplit.end(); ++ite) {
BDS_Edge *e = *ite;
- if (!e->deleted && e->numfaces() == 2){
+ if(!e->deleted && e->numfaces() == 2) {
const double coord = 0.5;
- BDS_Point *mid ;
+ 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);
- mid->lcBGM() = BGM_MeshSize(gf,
- (coord * e->p1->u + (1 - coord) * e->p2->u),
- (coord * e->p1->v + (1 - coord) * e->p2->v),
- mid->X, mid->Y, mid->Z);
+ mid->lcBGM() = BGM_MeshSize(
+ gf, (coord * e->p1->u + (1 - coord) * e->p2->u),
+ (coord * e->p1->v + (1 - coord) * e->p2->v), mid->X, mid->Y, mid->Z);
mid->lc() = 0.5 * (e->p1->lc() + e->p2->lc());
mid->u = coord * e->p1->u + (1 - coord) * e->p2->u;
mid->v = coord * e->p1->v + (1 - coord) * e->p2->v;
@@ -1017,9 +1060,9 @@ int solveInvalidPeriodic(GFace *gf, BDS_Mesh &m,
return toSplit.size();
}
-
-void optimizeMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
- std::map<BDS_Point*,MVertex*,PointLessThan> *recoverMap=0)
+void optimizeMeshBDS(
+ GFace *gf, BDS_Mesh &m, const int NIT,
+ std::map<BDS_Point *, MVertex *, PointLessThan> *recoverMap = 0)
{
return;
}
@@ -1028,19 +1071,19 @@ void optimizeMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
// build the BDS from a list of GFace
// This is a TRUE copy
-BDS_Mesh *gmsh2BDS(std::list<GFace*> &l)
+BDS_Mesh *gmsh2BDS(std::list<GFace *> &l)
{
BDS_Mesh *m = new BDS_Mesh;
- for (std::list<GFace*>::iterator it = l.begin(); it != l.end(); ++it){
+ for(std::list<GFace *>::iterator it = l.begin(); it != l.end(); ++it) {
GFace *gf = *it;
m->add_geom(gf->tag(), 2);
BDS_GeomEntity *g2 = m->get_geom(gf->tag(), 2);
- for (unsigned int i = 0; i < gf->triangles.size(); i++){
+ for(unsigned int i = 0; i < gf->triangles.size(); i++) {
MTriangle *e = gf->triangles[i];
BDS_Point *p[3];
- for (int j = 0; j < 3; j++){
+ for(int j = 0; j < 3; j++) {
p[j] = m->find_point(e->getVertex(j)->getNum());
- if (!p[j]) {
+ if(!p[j]) {
p[j] = m->add_point(e->getVertex(j)->getNum(), e->getVertex(j)->x(),
e->getVertex(j)->y(), e->getVertex(j)->z());
SPoint2 param;
diff --git a/Mesh/meshGFaceBDS.h b/Mesh/meshGFaceBDS.h
index 9927f218a43ba7c24feeb0c2dfc0c08d83a23c62..80e8a0a3350e62d416133ee1050b9a6e2487e582 100644
--- a/Mesh/meshGFaceBDS.h
+++ b/Mesh/meshGFaceBDS.h
@@ -17,20 +17,24 @@ class BDS_Point;
class MVertex;
class BDS_Mesh;
-void computeMeshSizeFieldAccuracy(GFace *gf, BDS_Mesh &m, double &avg,
- double &max_e, double &min_e, int &nE, int &GS);
-void refineMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
- const bool computeNodalSizeField,
- std::map<MVertex*, BDS_Point*> *recoverMapInv=0,
- std::map<BDS_Point*, MVertex*,PointLessThan> *recoverMap=0, std::vector<SPoint2> *true_boundary=0);
-void optimizeMeshBDS(GFace *gf, BDS_Mesh &m, const int NIT,
- std::map<BDS_Point*, MVertex*,PointLessThan> *recoverMap=0);
+void computeMeshSizeFieldAccuracy(GFace *gf, BDS_Mesh &m, double &avg,
+ double &max_e, double &min_e, int &nE,
+ int &GS);
+void refineMeshBDS(
+ GFace *gf, BDS_Mesh &m, const int NIT, const bool computeNodalSizeField,
+ std::map<MVertex *, BDS_Point *> *recoverMapInv = 0,
+ std::map<BDS_Point *, MVertex *, PointLessThan> *recoverMap = 0,
+ std::vector<SPoint2> *true_boundary = 0);
+void optimizeMeshBDS(
+ GFace *gf, BDS_Mesh &m, const int NIT,
+ std::map<BDS_Point *, MVertex *, PointLessThan> *recoverMap = 0);
void delaunayizeBDS(GFace *gf, BDS_Mesh &m, int &nb_swap);
void collapseSmallEdges(GModel &gm);
-BDS_Mesh *gmsh2BDS(std::list<GFace*> &l);
+BDS_Mesh *gmsh2BDS(std::list<GFace *> &l);
double computeEdgeLinearLength(BDS_Point *p1, BDS_Point *p2);
void smoothVertexPass(GFace *gf, BDS_Mesh &m, int &nb_smooth, bool q);
-void modifyInitialMeshToRemoveDegeneracies(GFace *gf, BDS_Mesh &m,
- std::map<BDS_Point*, MVertex*,PointLessThan> *recoverMap);
+void modifyInitialMeshToRemoveDegeneracies(
+ GFace *gf, BDS_Mesh &m,
+ std::map<BDS_Point *, MVertex *, PointLessThan> *recoverMap);
#endif
diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp
index 9999177376d8b20c02afbe28823c820ed9fd20b9..77a16aed8cf2610190f31446fe035fb3eca7da1f 100644
--- a/Mesh/meshGFaceDelaunayInsertion.cpp
+++ b/Mesh/meshGFaceDelaunayInsertion.cpp
@@ -56,16 +56,15 @@ void _printTris(char *name, ITERATOR it, ITERATOR end, bidimMeshData *data)
fprintf(ff, "View\"test\"{\n");
while(it != end) {
MTri3 *worst = *it;
- if (!worst->isDeleted()){
- if (data)
- fprintf(ff,"ST(%g,%g,%g,%g,%g,%g,%g,%g,%g) {%d,%d,%d};\n",
+ if(!worst->isDeleted()) {
+ if(data)
+ fprintf(ff, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g) {%d,%d,%d};\n",
data->Us[data->getIndex((worst)->tri()->getVertex(0))],
data->Vs[data->getIndex((worst)->tri()->getVertex(0))], 0.0,
data->Us[data->getIndex((worst)->tri()->getVertex(1))],
data->Vs[data->getIndex((worst)->tri()->getVertex(1))], 0.0,
data->Us[data->getIndex((worst)->tri()->getVertex(2))],
- data->Vs[data->getIndex((worst)->tri()->getVertex(2))],
- 0.0,
+ data->Vs[data->getIndex((worst)->tri()->getVertex(2))], 0.0,
(worst)->tri()->getVertex(0)->getNum(),
(worst)->tri()->getVertex(1)->getNum(),
(worst)->tri()->getVertex(2)->getNum());
@@ -114,7 +113,9 @@ static bool isActive(MTri3 *t, double limit_, int &i,
int ip1 = i - 1 < 0 ? 2 : i - 1;
int ip2 = i;
MEdge me(t->tri()->getVertex(ip1), t->tri()->getVertex(ip2));
- if(front->find(me) != front->end()) { return true; }
+ if(front->find(me) != front->end()) {
+ return true;
+ }
}
}
return false;
@@ -344,8 +345,7 @@ int inCircumCircleAniso(GFace *gf, MTriangle *base, const double *uv,
MTri3::MTri3(MTriangle *t, double lc, SMetric3 *metric, bidimMeshData *data,
GFace *gf)
- : deleted(false)
- , base(t)
+ : deleted(false), base(t)
{
neigh[0] = neigh[1] = neigh[2] = 0;
double center[3];
@@ -357,7 +357,9 @@ MTri3::MTri3(MTriangle *t, double lc, SMetric3 *metric, bidimMeshData *data,
base->getVertex(2)->z()};
if(!metric) {
- if(radiusNorm == 3) { circum_radius = 1. / base->gammaShapeMeasure(); }
+ if(radiusNorm == 3) {
+ circum_radius = 1. / base->gammaShapeMeasure();
+ }
else if(radiusNorm == 2) {
circumCenterXYZ(pa, pb, pc, center);
const double dx = base->getVertex(0)->x() - center[0];
@@ -453,7 +455,9 @@ void connectTris(Iterator beg, Iterator end, std::vector<edgeXface> &conn)
while(beg != end) {
if(!(*beg)->isDeleted()) {
- for(int j = 0; j < 3; j++) { conn.push_back(edgeXface(*beg, j)); }
+ for(int j = 0; j < 3; j++) {
+ conn.push_back(edgeXface(*beg, j));
+ }
}
++beg;
}
@@ -807,19 +811,20 @@ int insertVertexB(std::list<edgeXface> &shell, std::list<MTri3 *> &cavity,
}
// double din = t->getInnerRadius();
- double d1 = distance(v0,v);
- double d2 = distance(v1,v);
- double d3 = distance(v0,v1);
- SVector3 v0v1 (v1->x()-v0->x(),v1->y()-v0->y(),v1->z()-v0->z());
- SVector3 v0v (v->x()-v0->x(),v->y()-v0->y(),v->z()-v0->z());
- SVector3 pv = crossprod(v0v1,v0v);
+ double d1 = distance(v0, v);
+ double d2 = distance(v1, v);
+ double d3 = distance(v0, v1);
+ SVector3 v0v1(v1->x() - v0->x(), v1->y() - v0->y(), v1->z() - v0->z());
+ SVector3 v0v(v->x() - v0->x(), v->y() - v0->y(), v->z() - v0->z());
+ SVector3 pv = crossprod(v0v1, v0v);
double d4 = pv.norm() / d3;
// avoid angles that are too obtuse
double cosv = ((d1 * d1 + d2 * d2 - d3 * d3) / (2. * d1 * d2));
- if ((d1 < LL * .5 || d2 < LL * .5 || d4 < LL * .10 || cosv < -.9999) && !force) {
+ if((d1 < LL * .5 || d2 < LL * .5 || d4 < LL * .10 || cosv < -.9999) &&
+ !force) {
onePointIsTooClose = true;
- //printf("%12.5E %12.5E %12.5E %12.5E \n",d1,d2,LL,cosv);
+ // printf("%12.5E %12.5E %12.5E %12.5E \n",d1,d2,LL,cosv);
}
newTris[k++] = t4;
@@ -1005,7 +1010,9 @@ static MTri3 *search4Triangle(MTri3 *t, double pt[2], bidimMeshData &data,
itx != AllTris.end(); ++itx) {
if(!(*itx)->isDeleted()) {
inside = invMapUV((*itx)->tri(), pt, data, uv, 1.e-8);
- if(inside) { return *itx; }
+ if(inside) {
+ return *itx;
+ }
}
}
printf("argh %g %g!!!!\n", pt[0], pt[1]);
@@ -1421,10 +1428,9 @@ bool optimalPointFrontalB(GFace *gf, MTri3 *worst, int active_edge,
return true;
}
-void bowyerWatsonFrontal(GFace *gf,
- std::map<MVertex* , MVertex*>* equivalence,
- std::map<MVertex*, SPoint2> * parametricCoordinates,
- std::vector<SPoint2> *true_boundary)
+void bowyerWatsonFrontal(GFace *gf, std::map<MVertex *, MVertex *> *equivalence,
+ std::map<MVertex *, SPoint2> *parametricCoordinates,
+ std::vector<SPoint2> *true_boundary)
{
std::set<MTri3 *, compareTri3Ptr> AllTris;
std::set<MTri3 *, compareTri3Ptr> ActiveTris;
@@ -1452,10 +1458,10 @@ void bowyerWatsonFrontal(GFace *gf,
break;
}
- Range<double> RU = gf->parBounds (0);
- Range<double> RV = gf->parBounds (1);
- SPoint2 FAR (2*RU.high(),2*RV.high());
-
+ Range<double> RU = gf->parBounds(0);
+ Range<double> RV = gf->parBounds(1);
+ SPoint2 FAR(2 * RU.high(), 2 * RV.high());
+
// insert points
int ITERATION = 0;
while(1) {
@@ -1481,28 +1487,30 @@ void bowyerWatsonFrontal(GFace *gf,
Msg::Debug("%7d points created -- Worst tri radius is %8.3f",
gf->mesh_vertices.size(), worst->getRadius());
double newPoint[2], metric[3];
- //optimalPointFrontal (gf,worst,active_edge,Us,Vs,vSizes,vSizesBGM,newPoint,metric);
- if (optimalPointFrontalB (gf,worst,active_edge,DATA,newPoint,metric)){
- // printf("iteration %d passes first round\n",ITERATION);
- SPoint2 NP(newPoint[0],newPoint[1]);
- int nnnn;
- if (!true_boundary ||
- pointInsideParametricDomain (*true_boundary, NP, FAR, nnnn))
- insertAPoint(gf, AllTris.end(), newPoint, metric, DATA, AllTris, &ActiveTris, worst, NULL, testStarShapeness);
+ // optimalPointFrontal
+ // (gf,worst,active_edge,Us,Vs,vSizes,vSizesBGM,newPoint,metric);
+ if(optimalPointFrontalB(gf, worst, active_edge, DATA, newPoint, metric)) {
+ // printf("iteration %d passes first round\n",ITERATION);
+ SPoint2 NP(newPoint[0], newPoint[1]);
+ int nnnn;
+ if(!true_boundary ||
+ pointInsideParametricDomain(*true_boundary, NP, FAR, nnnn))
+ insertAPoint(gf, AllTris.end(), newPoint, metric, DATA, AllTris,
+ &ActiveTris, worst, NULL, testStarShapeness);
}
}
}
// nbSwaps = edgeSwapPass(gf, AllTris, SWCR_QUAL, DATA);
char name[245];
- sprintf(name,"delFrontal_GFace_%d.pos",gf->tag());
- _printTris (name, AllTris.begin(), AllTris.end(), &DATA);
- sprintf(name,"delFrontal_GFace_%d_Real.pos",gf->tag());
- _printTris (name, AllTris.begin(), AllTris.end(), NULL);
- // sprintf(name,"delFrontal_GFace_%d_Layer_Real%d.pos",gf->tag(),ITERATION);
- // _printTris (name, AllTris.begin(), AllTris.end(),NULL);
- // sprintf(name,"delFrontal_GFace_%d_Layer_%d_Active.pos",gf->tag(),ITERATION);
- // _printTris (name, ActiveTris.begin(), ActiveTris.end(), &DATA);
+ sprintf(name, "delFrontal_GFace_%d.pos", gf->tag());
+ _printTris(name, AllTris.begin(), AllTris.end(), &DATA);
+ sprintf(name, "delFrontal_GFace_%d_Real.pos", gf->tag());
+ _printTris(name, AllTris.begin(), AllTris.end(), NULL);
+ // sprintf(name,"delFrontal_GFace_%d_Layer_Real%d.pos",gf->tag(),ITERATION);
+ // _printTris (name, AllTris.begin(), AllTris.end(),NULL);
+ // sprintf(name,"delFrontal_GFace_%d_Layer_%d_Active.pos",gf->tag(),ITERATION);
+ // _printTris (name, ActiveTris.begin(), ActiveTris.end(), &DATA);
transferDataStructure(gf, AllTris, DATA);
// removeThreeTrianglesNodes(gf);
@@ -2012,7 +2020,9 @@ MTri3 *getTriToBreak(MVertex *v, std::vector<MTri3 *> &t, int &NB_GLOBAL_SEARCH,
// last inserted is used as starting point
// we know it is not deleted
unsigned int k = t.size() - 1;
- while(t[k]->isDeleted()) { k--; }
+ while(t[k]->isDeleted()) {
+ k--;
+ }
MTri3 *start = t[k];
start = search4Triangle(start, v, (int)t.size(), ITER);
if(start) return start;
@@ -2116,7 +2126,9 @@ void delaunayMeshIn2D(std::vector<MVertex *> &v,
for(unsigned int k = 0; k < std::min(cavity.size(), shell.size()); k++) {
cavity[k]->setDeleted(false);
- for(unsigned int l = 0; l < 3; l++) { cavity[k]->setNeigh(l, 0); }
+ for(unsigned int l = 0; l < 3; l++) {
+ cavity[k]->setNeigh(l, 0);
+ }
}
T = Cpu();
connectTris(extended_cavity.begin(), extended_cavity.end(), conn);
diff --git a/Mesh/meshGFaceDelaunayInsertion.h b/Mesh/meshGFaceDelaunayInsertion.h
index 5c916c18b925b7944de0dc09babfe932665d05bc..0327ba36db1c1afad6ba0950be1386409c862c62 100644
--- a/Mesh/meshGFaceDelaunayInsertion.h
+++ b/Mesh/meshGFaceDelaunayInsertion.h
@@ -65,8 +65,7 @@ struct bidimMeshData {
}
bidimMeshData(std::map<MVertex *, MVertex *> *e = 0,
std::map<MVertex *, SPoint2> *p = 0)
- : equivalence(e)
- , parametricCoordinates(p)
+ : equivalence(e), parametricCoordinates(p)
{
}
};
@@ -153,21 +152,22 @@ public:
}
};
-void connectTriangles(std::list<MTri3*> &);
-void connectTriangles(std::vector<MTri3*> &);
-void connectTriangles(std::set<MTri3*, compareTri3Ptr> &AllTris);
-void bowyerWatson(GFace *gf, int MAXPNT= 1000000000,
- std::map<MVertex*, MVertex*>* equivalence= 0,
- std::map<MVertex*, SPoint2> * parametricCoordinates = 0);
-void bowyerWatsonFrontal(GFace *gf,
- std::map<MVertex*, MVertex*>* equivalence= 0,
- std::map<MVertex*, SPoint2> * parametricCoordinates = 0, std::vector<SPoint2> *true_boundary = 0);
-void bowyerWatsonFrontalLayers(GFace *gf, bool quad,
- std::map<MVertex*, MVertex*>* equivalence= 0,
- std::map<MVertex*, SPoint2> * parametricCoordinates = 0);
-void bowyerWatsonParallelograms(GFace *gf,
- std::map<MVertex*, MVertex*>* equivalence= 0,
- std::map<MVertex*, SPoint2> * parametricCoordinates = 0);
+void connectTriangles(std::list<MTri3 *> &);
+void connectTriangles(std::vector<MTri3 *> &);
+void connectTriangles(std::set<MTri3 *, compareTri3Ptr> &AllTris);
+void bowyerWatson(GFace *gf, int MAXPNT = 1000000000,
+ std::map<MVertex *, MVertex *> *equivalence = 0,
+ std::map<MVertex *, SPoint2> *parametricCoordinates = 0);
+void bowyerWatsonFrontal(
+ GFace *gf, std::map<MVertex *, MVertex *> *equivalence = 0,
+ std::map<MVertex *, SPoint2> *parametricCoordinates = 0,
+ std::vector<SPoint2> *true_boundary = 0);
+void bowyerWatsonFrontalLayers(
+ GFace *gf, bool quad, std::map<MVertex *, MVertex *> *equivalence = 0,
+ std::map<MVertex *, SPoint2> *parametricCoordinates = 0);
+void bowyerWatsonParallelograms(
+ GFace *gf, std::map<MVertex *, MVertex *> *equivalence = 0,
+ std::map<MVertex *, SPoint2> *parametricCoordinates = 0);
void bowyerWatsonParallelogramsConstrained(
GFace *gf, const std::set<MVertex *> &constr_vertices,
std::map<MVertex *, MVertex *> *equivalence = 0,
@@ -186,10 +186,7 @@ struct edgeXface {
MTri3 *t1;
int i1;
int ori;
- edgeXface(MTri3 *_t, int iFac)
- : t1(_t)
- , i1(iFac)
- , ori(1)
+ edgeXface(MTri3 *_t, int iFac) : t1(_t), i1(iFac), ori(1)
{
v[0] = t1->tri()->getVertex(iFac == 0 ? 2 : iFac - 1);
v[1] = t1->tri()->getVertex(iFac);