diff --git a/Common/CommandLine.cpp b/Common/CommandLine.cpp index e89b4b883e754ba37fb8a67708d04203540bed95..7dbb2746ecb654510e215afa3db32bdef26ff5bc 100644 --- a/Common/CommandLine.cpp +++ b/Common/CommandLine.cpp @@ -101,7 +101,6 @@ void PrintUsage(const char *name) Msg::Direct(" -fontsize int Specify the font size for the GUI"); Msg::Direct(" -theme string Specify FLTK GUI theme"); Msg::Direct(" -display string Specify display"); - Msg::Direct(" -showCompounds Shows the underlying surfaces/edges/mesh of compounds"); #endif Msg::Direct("Other options:"); Msg::Direct(" - Parse input files, then exit"); @@ -749,7 +748,7 @@ void GetOptions(int argc, char *argv[]) Msg::Fatal("Missing argument"); } else if(!strcmp(argv[i] + 1, "showCompounds")) { - CTX::instance()->showCompounds = 1; + CTX::instance()->geom.hideCompounds = 0; i++; } #endif diff --git a/Common/Context.cpp b/Common/Context.cpp index bb864a3ea7612c330087fffcf7f5579b56f4f080..cd8d6386d17dd6a565b97f8ca5660ea9793246f5 100644 --- a/Common/Context.cpp +++ b/Common/Context.cpp @@ -41,7 +41,7 @@ CTX::CTX() bgmFileName = ""; createAppendMeshStatReport = 0; lc = 1.; - min[0] = min[1] = min[2] = max[2] = 0.; + min[0] = min[1] = min[2] = max[2] = 0.; max[0] = max[1] = 1.; // for nice view when adding point in new model cg[0] = cg[1] = cg[2] = 0.; polygonOffset = 0; @@ -54,7 +54,7 @@ CTX::CTX() post.draw = 1; post.combineTime = 0; // try to combineTime views at startup lock = 0; // very primitive locking - + #if defined(HAVE_FLTK) glFontEnum = FL_HELVETICA; #else @@ -62,7 +62,6 @@ CTX::CTX() #endif forcedBBox = 0; hideUnselected = 0; - showCompounds = 0; numWindows = numTiles = 1; deltaFontSize = 0; recentFiles.resize(5); diff --git a/Common/Context.h b/Common/Context.h index 6e8f6a362cf749744a9fc6d40201e3e67ef7bdff..2eec2b6416112bbed17fa17af1fd54bcc3996bda 100644 --- a/Common/Context.h +++ b/Common/Context.h @@ -53,6 +53,7 @@ struct contextGeometryOptions { int occSewFaces, occConnectFaces; int copyMeshingMethod, exactExtrusion; int matchGeomAndMesh; + int hideCompounds; }; class CTX { @@ -190,8 +191,6 @@ class CTX { int printing; // hide all unselected entities? int hideUnselected; - // hide underlying curves and surfaces of compounds (makes work a lot easier) - int showCompounds; // temporary storage of rotation, translation, scale (until the GUI // is ready) double tmpRotation[3], tmpTranslation[3], tmpScale[3], tmpQuaternion[4]; diff --git a/Common/DefaultOptions.h b/Common/DefaultOptions.h index bcca4c18c8978ee2b232fa8a1d1313efe0783eb4..9c762e393ed4524a6e051cbde66e6b89b64ca295 100644 --- a/Common/DefaultOptions.h +++ b/Common/DefaultOptions.h @@ -653,6 +653,8 @@ StringXNumber GeometryOptions_Number[] = { { F|O, "ExtrudeSplinePoints" , opt_geometry_extrude_spline_points, 5. , "Number of control points for splines created during extrusion" }, + { F|O, "HideCompounds" , opt_geometry_hide_compounds, 1. , + "Hide entities that make up compound entities?" }, { F|O, "HighlightOrphans" , opt_geometry_highlight_orphans, 0. , "Highlight orphan entities (lines connected to a single surface, etc.)?" }, diff --git a/Common/Options.cpp b/Common/Options.cpp index 320783ca5fcd49849b47f78f9818c5dab907897e..db3ab7fff4dfdc9c4a2545b0375f6ec8a3d6cd98 100644 --- a/Common/Options.cpp +++ b/Common/Options.cpp @@ -3646,6 +3646,21 @@ double opt_geometry_auto_coherence(OPT_ARGS_NUM) return CTX::instance()->geom.autoCoherence; } +double opt_geometry_hide_compounds(OPT_ARGS_NUM) +{ + if(action & GMSH_SET){ + if(CTX::instance()->geom.hideCompounds != val) + CTX::instance()->mesh.changed |= (ENT_LINE | ENT_SURFACE); + CTX::instance()->geom.hideCompounds = (int)val; + } +#if defined(HAVE_FLTK) + // if(FlGui::available() && (action & GMSH_GUI)) + // FlGui::instance()->options->geo.butt[10]->value + // (CTX::instance()->geom.hideCompounds); +#endif + return CTX::instance()->geom.hideCompounds; +} + double opt_geometry_highlight_orphans(OPT_ARGS_NUM) { if(action & GMSH_SET) diff --git a/Common/Options.h b/Common/Options.h index b2c1d7f88d97e28986e9d6f51987c728eac2a70e..82e81bc682b5748264a8a3a9ad7a201661bff857 100644 --- a/Common/Options.h +++ b/Common/Options.h @@ -284,6 +284,7 @@ double opt_geometry_offset0(OPT_ARGS_NUM); double opt_geometry_offset1(OPT_ARGS_NUM); double opt_geometry_offset2(OPT_ARGS_NUM); double opt_geometry_auto_coherence(OPT_ARGS_NUM); +double opt_geometry_hide_compounds(OPT_ARGS_NUM); double opt_geometry_highlight_orphans(OPT_ARGS_NUM); double opt_geometry_tolerance(OPT_ARGS_NUM); double opt_geometry_normals(OPT_ARGS_NUM); diff --git a/Fltk/onelabWindow.cpp b/Fltk/onelabWindow.cpp index 33baf52cc3ffc93b3fe97530550df69ee79435e7..4aa18a04651ded41d5dcb46ab60da19ac98211dc 100644 --- a/Fltk/onelabWindow.cpp +++ b/Fltk/onelabWindow.cpp @@ -533,11 +533,11 @@ static void updateOnelabGraphs() drawContext::global()->draw(); } -static void importPhysicalGroups(GModel *m) +static void importPhysicalGroups(onelab::client *c, GModel *m) { std::map<int, std::vector<GEntity*> > groups[4]; m->getPhysicalGroups(groups); - + // create "read-only" onelab groups } static void runGmshClient(const std::string &action) @@ -565,7 +565,7 @@ static void runGmshClient(const std::string &action) // the model name has changed modelName = GModel::current()->getName(); geometry_reload_cb(0, 0); - importPhysicalGroups(GModel::current()); + importPhysicalGroups(c, GModel::current()); } } else if(action == "compute"){ @@ -576,7 +576,7 @@ static void runGmshClient(const std::string &action) // changed modelName = GModel::current()->getName(); geometry_reload_cb(0, 0); - importPhysicalGroups(GModel::current()); + importPhysicalGroups(c, GModel::current()); if(FlGui::instance()->onelab->meshAuto()){ mesh_3d_cb(0, 0); CreateOutputFile(mshFileName, CTX::instance()->mesh.fileFormat); diff --git a/Geo/GEdge.cpp b/Geo/GEdge.cpp index 6443d386119173f00b701cd1244bb89fb6584078..df1b99bf8f1893b3eb1889a6f5362bd422093149 100644 --- a/Geo/GEdge.cpp +++ b/Geo/GEdge.cpp @@ -79,9 +79,9 @@ MElement *const *GEdge::getStartElementType(int type) const } MElement *GEdge::getMeshElement(unsigned int index) -{ +{ if(index < lines.size()) - return lines[index]; + return lines[index]; return 0; } @@ -140,7 +140,7 @@ SOrientedBoundingBox GEdge::getOBB() SPoint3 pt2(getEndVertex()->x(), getEndVertex()->y(), getEndVertex()->z()); vertices.push_back(pt1); vertices.push_back(pt2); - } + } else if(geomType() != DiscreteCurve && geomType() != BoundaryLayerCurve){ Range<double> tr = this->parBounds(0); // N can be choosen arbitrarily, but 10 points seems reasonable @@ -151,7 +151,7 @@ SOrientedBoundingBox GEdge::getOBB() SPoint3 pt(p.x(), p.y(), p.z()); vertices.push_back(pt); } - } + } else { SPoint3 dummy(0, 0, 0); vertices.push_back(dummy); @@ -163,15 +163,16 @@ SOrientedBoundingBox GEdge::getOBB() void GEdge::setVisibility(char val, bool recursive) { - if (getCompound() && !CTX::instance()->showCompounds) { - // use visibility info of compound edge if this edge belongs to it + if (getCompound() && CTX::instance()->geom.hideCompounds) { + // use visibility info of compound edge if this edge belongs to it GEntity::setVisibility(0); if(v0) v0->setVisibility(0); if(v1) v1->setVisibility(0); bool val2 = getCompound()->getVisibility(); if(getCompound()->v0) getCompound()->v0->setVisibility(val2); if(getCompound()->v1) getCompound()->v1->setVisibility(val2); - } else { + } + else { GEntity::setVisibility(val); if(recursive){ if(v0) v0->setVisibility(val); @@ -192,7 +193,7 @@ void GEdge::writeGEO(FILE *fp) if(!getBeginVertex() || !getEndVertex() || geomType() == DiscreteCurve) return; if(geomType() == Line){ - fprintf(fp, "Line(%d) = {%d, %d};\n", + fprintf(fp, "Line(%d) = {%d, %d};\n", tag(), getBeginVertex()->tag(), getEndVertex()->tag()); } else{ @@ -205,7 +206,7 @@ void GEdge::writeGEO(FILE *fp) for(int i = 1; i < N; i++){ double u = umin + (double)i / N * (umax - umin); GPoint p = point(u); - fprintf(fp, "Point(p%d + %d) = {%.16g, %.16g, %.16g};\n", + fprintf(fp, "Point(p%d + %d) = {%.16g, %.16g, %.16g};\n", tag(), i, p.x(), p.y(), p.z()); } fprintf(fp, "Spline(%d) = {%d", tag(), getBeginVertex()->tag()); @@ -215,7 +216,7 @@ void GEdge::writeGEO(FILE *fp) } if(meshAttributes.Method == MESH_TRANSFINITE){ - fprintf(fp, "Transfinite Line {%d} = %d", + fprintf(fp, "Transfinite Line {%d} = %d", tag() * (meshAttributes.typeTransfinite > 0 ? 1 : -1), meshAttributes.nbPointsTransfinite); if(meshAttributes.typeTransfinite){ @@ -267,11 +268,11 @@ double GEdge::curvature(double par) const { SVector3 d1 = firstDer(par); SVector3 d2 = secondDer(par); - + double one_over_norm = 1. / norm(d1); SVector3 cross_prod = crossprod(d1,d2); - + return ( norm(cross_prod) * one_over_norm * one_over_norm * one_over_norm ); } @@ -292,35 +293,35 @@ double GEdge::length(const double &u0, const double &u1, const int nbQuadPoints) /* consider a curve x(t) and a point y - - use a golden section algorithm that minimizes - min_t \|x(t)-y\| + use a golden section algorithm that minimizes + + min_t \|x(t)-y\| */ const double GOLDEN = (1. + sqrt(5.)) / 2.; const double GOLDEN2 = 2 - GOLDEN; - + // x1 and x3 are the current bounds; the minimum is between them. // x2 is the center point, which is closer to x1 than to x3 -double goldenSectionSearch(const GEdge *ge, const SPoint3 &q, double x1, +double goldenSectionSearch(const GEdge *ge, const SPoint3 &q, double x1, double x2, double x3, double tau) -{ +{ // Create a new possible center in the area between x2 and x3, closer to x2 double x4 = x2 + GOLDEN2 * (x3 - x2); - + // Evaluate termination criterion if (fabs(x3 - x1) < tau * (fabs(x2) + fabs(x4))) return (x3 + x1) / 2; - + const SVector3 dp4 = q - ge->position(x4); const SVector3 dp2 = q - ge->position(x2); - + const double d4 = dp4.norm(); const double d2 = dp2.norm(); - + if (d4 < d2) return goldenSectionSearch(ge, q, x2, x4, x3, tau); else @@ -330,13 +331,13 @@ double goldenSectionSearch(const GEdge *ge, const SPoint3 &q, double x1, GPoint GEdge::closestPoint(const SPoint3 &q, double &t) const { // printf("looking for closest point in curve %d to point %g %g\n",tag(),q.x(),q.y()); - + const int nbSamples = 100; - + double tolerance = 1.e-12; Range<double> interval = parBounds(0); - + double tMin = std::min(interval.high(), interval.low()); double tMax = std::max(interval.high(), interval.low()); @@ -352,7 +353,7 @@ GPoint GEdge::closestPoint(const SPoint3 &q, double &t) const DMIN = D; } } - + // printf("parameter %g as an initial guess (dist = %g)\n",topt,DMIN); if (topt == tMin) @@ -366,7 +367,7 @@ GPoint GEdge::closestPoint(const SPoint3 &q, double &t) const const double D = dp.norm(); // printf("after golden section parameter %g (dist = %g)\n",t,D); - + return point(t); } @@ -377,7 +378,7 @@ double GEdge::parFromPoint(const SPoint3 &P) const return t; } -bool GEdge::XYZToU(const double X, const double Y, const double Z, +bool GEdge::XYZToU(const double X, const double Y, const double Z, double &u, const double relax) const { const int MaxIter = 25; @@ -391,12 +392,12 @@ bool GEdge::XYZToU(const double X, const double Y, const double Z, double uMax = uu.high(); SVector3 Q(X, Y, Z), P; - + double init[NumInitGuess]; - - for (int i = 0; i < NumInitGuess; i++) + + for (int i = 0; i < NumInitGuess; i++) init[i] = uMin + (uMax - uMin) / (NumInitGuess - 1) * i; - + for(int i = 0; i < NumInitGuess; i++){ u = init[i]; double uNew = u; @@ -405,30 +406,30 @@ bool GEdge::XYZToU(const double X, const double Y, const double Z, SVector3 dPQ = P - Q; err2 = dPQ.norm(); - - if (err2 < 1.e-8 * CTX::instance()->lc) return true; - + + if (err2 < 1.e-8 * CTX::instance()->lc) return true; + while(iter++ < MaxIter && err2 > 1e-8 * CTX::instance()->lc) { SVector3 der = firstDer(u); uNew = u - relax * dot(dPQ,der) / dot(der,der); uNew = std::min(uMax,std::max(uMin,uNew)); P = position(uNew); - + dPQ = P - Q; err2 = dPQ.norm(); err = fabs(uNew - u); u = uNew; - } - + } + if (err2 < 1e-8 * CTX::instance()->lc) return true; } - + if(relax > 1.e-2) { - // Msg::Info("point %g %g %g on edge %d : Relaxation factor = %g", + // Msg::Info("point %g %g %g on edge %d : Relaxation factor = %g", // Q.x(), Q.y(), Q.z(), 0.75 * relax); return XYZToU(Q.x(), Q.y(), Q.z(), u, 0.75 * relax); } - + // Msg::Error("Could not converge reparametrisation of point (%e,%e,%e) on edge %d", // Q.x(), Q.y(), Q.z(), tag()); return false; @@ -446,5 +447,5 @@ void GEdge::replaceEndingPoints(GVertex *replOfv0, GVertex *replOfv1) v1->delEdge(this); replOfv1->addEdge(this); v1 = replOfv1; - } + } } diff --git a/Geo/GFace.cpp b/Geo/GFace.cpp index dcf554ec36dfe648fed1bf641fec63dbec6f8674..5046574b1176a815f59301f8ec12c7d7c1fc21ba 100644 --- a/Geo/GFace.cpp +++ b/Geo/GFace.cpp @@ -263,7 +263,7 @@ std::list<GVertex*> GFace::vertices() const void GFace::setVisibility(char val, bool recursive) { - if (getCompound() && !CTX::instance()->showCompounds) { + if (getCompound() && CTX::instance()->geom.hideCompounds) { GEntity::setVisibility(0); for (std::list<GEdge*>::iterator it = l_edges.begin(); it != l_edges.end(); ++it) (*it)->setVisibility(0, true); @@ -276,7 +276,8 @@ void GFace::setVisibility(char val, bool recursive) else (*it)->setVisibility(val2, true); } - } else { + } + else { GEntity::setVisibility(val); if(recursive){ for (std::list<GEdge*>::iterator it = l_edges.begin(); it != l_edges.end(); ++it) @@ -955,7 +956,7 @@ GPoint GFace::closestPoint(const SPoint3 &queryPoint, const double initialGuess[ // Creating the optimisation problem // alglib::ae_int_t dim = 2; - alglib::ae_int_t corr = 2; // Num of corrections in the scheme in [3,7] + alglib::ae_int_t corr = 2; // Num of corrections in the scheme in [3,7] alglib::minlbfgsstate state; alglib::real_1d_array x; // = "[0.5,0.5]"; std::vector<double> initial_conditions(dim); diff --git a/Geo/GModelIO_Geo.cpp b/Geo/GModelIO_Geo.cpp index 982a19425ba45d42287d1ff1b58b83ba503129ae..5a8aee90237d2082653df0f026becfe262cc5543 100644 --- a/Geo/GModelIO_Geo.cpp +++ b/Geo/GModelIO_Geo.cpp @@ -145,9 +145,9 @@ int GModel::importGEOInternals() } e = new GEdgeCompound(this, c->Num, comp); add(e); - if (!CTX::instance()->showCompounds) + if(CTX::instance()->geom.hideCompounds) for (std::vector<GEdge*>::iterator it = comp.begin(); it != comp.end(); ++it) - (*it)->setVisibility(0,true); + (*it)->setVisibility(0, true); } else if(!e && c->beg && c->end){ e = new gmshEdge(this, c, @@ -202,9 +202,9 @@ int GModel::importGEOInternals() f->meshAttributes.Method = s->Method; f->meshAttributes.extrude = s->Extrude; add(f); - if (!CTX::instance()->showCompounds) + if (CTX::instance()->geom.hideCompounds) for (std::list<GFace*>::iterator it = comp.begin(); it != comp.end(); ++it) - (*it)->setVisibility(0,true); + (*it)->setVisibility(0, true); if(s->EmbeddedCurves){ for(int i = 0; i < List_Nbr(s->EmbeddedCurves); i++){ Curve *c; diff --git a/Geo/GModelVertexArrays.cpp b/Geo/GModelVertexArrays.cpp index 4902edaa623faf3e0c5909302e7e7b923b576337..1c452dba3b261ae1397c2d895db5d4b8c9fdde54 100644 --- a/Geo/GModelVertexArrays.cpp +++ b/Geo/GModelVertexArrays.cpp @@ -7,8 +7,6 @@ #include "GmshMessage.h" #include "GmshDefines.h" #include "GModel.h" -#include "GFaceCompound.h" -#include "GEdgeCompound.h" #include "MLine.h" #include "MTriangle.h" #include "MQuadrangle.h" @@ -62,7 +60,7 @@ static unsigned int getColorByElement(MElement *ele) else if(CTX::instance()->mesh.colorCarousel == 3){ // by partition return CTX::instance()->color.mesh.carousel[abs(ele->getPartition() % 20)]; } - else{ + else{ // by elementary or physical entity (this is not perfect (since // e.g. a triangle can have no vertices categorized on a surface), // but it's the best we can do "fast" since we don't store the @@ -78,9 +76,9 @@ static unsigned int getColorByElement(MElement *ele) static double evalClipPlane(int clip, double x, double y, double z) { - return CTX::instance()->clipPlane[clip][0] * x + - CTX::instance()->clipPlane[clip][1] * y + - CTX::instance()->clipPlane[clip][2] * z + + return CTX::instance()->clipPlane[clip][0] * x + + CTX::instance()->clipPlane[clip][1] * y + + CTX::instance()->clipPlane[clip][2] * z + CTX::instance()->clipPlane[clip][3]; } @@ -109,12 +107,12 @@ bool isElementVisible(MElement *ele) q = ele->etaShapeMeasure(); else q = ele->gammaShapeMeasure(); - if(q < CTX::instance()->mesh.qualityInf || + if(q < CTX::instance()->mesh.qualityInf || q > CTX::instance()->mesh.qualitySup) return false; } if(CTX::instance()->mesh.radiusSup) { double r = ele->maxEdge(); - if(r < CTX::instance()->mesh.radiusInf || + if(r < CTX::instance()->mesh.radiusInf || r > CTX::instance()->mesh.radiusSup) return false; } if(CTX::instance()->clipWholeElements){ @@ -125,7 +123,7 @@ bool isElementVisible(MElement *ele) } else{ double d = intersectClipPlane(clip, ele); - if(ele->getDim() == 3 && + if(ele->getDim() == 3 && CTX::instance()->clipOnlyDrawIntersectingVolume && d){ hidden = true; break; @@ -189,7 +187,7 @@ static void addElementsInArrays(GEntity *e, std::vector<T*> &elements, MElement *ele = elements[i]; if(!isElementVisible(ele) || ele->getDim() < 1) continue; - + unsigned int c = getColorByElement(ele); unsigned int col[4] = {c, c, c, c}; @@ -253,15 +251,8 @@ class initMeshGEdge { public: void operator () (GEdge *e) { - if(!e->getVisibility()) { - if(e->getCompound()) { - if(!e->getCompound()->getVisibility()) return; - } - else - return; - } - e->deleteVertexArrays(); + if(!e->getVisibility()) return; e->setAllElementsVisible(CTX::instance()->mesh.lines && areAllElementsVisible(e->lines)); @@ -290,7 +281,7 @@ class initMeshGFace { { int num = 0; if(CTX::instance()->mesh.surfacesEdges){ - num += (3 * f->triangles.size() + 4 * f->quadrangles.size() + + num += (3 * f->triangles.size() + 4 * f->quadrangles.size() + 4 * f->polygons.size()) / 2; if(CTX::instance()->mesh.explode != 1.) num *= 2; if(_curved) num *= 2; @@ -310,16 +301,10 @@ class initMeshGFace { public: void operator () (GFace *f) { - if(!f->getVisibility()) { - if(f->getCompound()) { - if(!f->getCompound()->getVisibility()) return; - } else - return; - } - f->deleteVertexArrays(); + if(!f->getVisibility()) return; f->setAllElementsVisible - (CTX::instance()->mesh.triangles && areAllElementsVisible(f->triangles) && + (CTX::instance()->mesh.triangles && areAllElementsVisible(f->triangles) && CTX::instance()->mesh.quadrangles && areAllElementsVisible(f->quadrangles)); bool edg = CTX::instance()->mesh.surfacesEdges; @@ -343,7 +328,7 @@ class initMeshGRegion { bool _curved; int _estimateIfClipped(int num) { - if(CTX::instance()->clipWholeElements && + if(CTX::instance()->clipWholeElements && CTX::instance()->clipOnlyDrawIntersectingVolume){ for(int clip = 0; clip < 6; clip++){ if(CTX::instance()->mesh.clip & (1 << clip)) @@ -385,10 +370,9 @@ class initMeshGRegion { } public: void operator () (GRegion *r) - { - if(!r->getVisibility()) return; - + { r->deleteVertexArrays(); + if(!r->getVisibility()) return; r->setAllElementsVisible (CTX::instance()->mesh.tetrahedra && areAllElementsVisible(r->tetrahedra) && CTX::instance()->mesh.hexahedra && areAllElementsVisible(r->hexahedra) && @@ -398,7 +382,7 @@ class initMeshGRegion { bool edg = CTX::instance()->mesh.volumesEdges; bool fac = CTX::instance()->mesh.volumesFaces; if(edg || fac){ - _curved = (areSomeElementsCurved(r->tetrahedra) || + _curved = (areSomeElementsCurved(r->tetrahedra) || areSomeElementsCurved(r->hexahedra) || areSomeElementsCurved(r->prisms) || areSomeElementsCurved(r->pyramids)); @@ -425,7 +409,7 @@ void GModel::fillVertexArrays() if(status >= 1 && CTX::instance()->mesh.changed & ENT_LINE) std::for_each(firstEdge(), lastEdge(), initMeshGEdge()); - + if(status >= 2 && CTX::instance()->mesh.changed & ENT_SURFACE){ if(normals) delete normals; normals = new smooth_normals(CTX::instance()->mesh.angleSmoothNormals); @@ -433,7 +417,7 @@ void GModel::fillVertexArrays() std::for_each(firstFace(), lastFace(), initSmoothNormalsGFace()); std::for_each(firstFace(), lastFace(), initMeshGFace()); } - + if(status >= 3 && CTX::instance()->mesh.changed & ENT_VOLUME) std::for_each(firstRegion(), lastRegion(), initMeshGRegion()); } diff --git a/Mesh/meshGFace.cpp b/Mesh/meshGFace.cpp index 125529c0556b5fa08e3163a122f94712884ee7bb..8a4f82a191a85a5959ad93e103f7c1ed8e7f1f08 100644 --- a/Mesh/meshGFace.cpp +++ b/Mesh/meshGFace.cpp @@ -101,7 +101,7 @@ static void copyMesh(GFace *source, GFace *target) } if (count < 2) return; - + const double t1u = param_target[0].x(), t1v = param_target[0].y(); const double t2u = param_target[1].x(), t2v = param_target[1].y(); const double s1u = param_source[0].x(), s1v = param_source[0].y(); @@ -133,7 +133,7 @@ static void copyMesh(GFace *source, GFace *target) for (unsigned i = 0; i < source->triangles.size(); i++){ MVertex *vt[3]; for (int j = 0; j < 3; j++){ - MVertex *vs = source->triangles[i]->getVertex(j); + MVertex *vs = source->triangles[i]->getVertex(j); vt[j] = vs2vt[vs]; if (!vt[j]){ SPoint2 p; @@ -202,7 +202,7 @@ static bool noSeam(GFace *gf) } static void remeshUnrecoveredEdges(std::map<MVertex*, BDS_Point*> &recoverMapInv, - std::set<EdgeToRecover> &edgesNotRecovered, + std::set<EdgeToRecover> &edgesNotRecovered, std::list<GFace*> &facesToRemesh) { facesToRemesh.clear(); @@ -218,7 +218,7 @@ static void remeshUnrecoveredEdges(std::map<MVertex*, BDS_Point*> &recoverMapInv dem(*it); } } - + // add a new point in the middle of the intersecting segment int p1 = itr->p1; int p2 = itr->p2; @@ -256,13 +256,13 @@ static void remeshUnrecoveredEdges(std::map<MVertex*, BDS_Point*> &recoverMapInv 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(); if(v2->onWhat() == g1 && v2->onWhat() == g2) 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()); if(lc2 == -1) @@ -300,10 +300,10 @@ static bool algoDelaunay2D(GFace *gf) return false; if(CTX::instance()->mesh.algo2d == ALGO_2D_DELAUNAY || - CTX::instance()->mesh.algo2d == ALGO_2D_BAMG || - CTX::instance()->mesh.algo2d == ALGO_2D_FRONTAL || - CTX::instance()->mesh.algo2d == ALGO_2D_FRONTAL_QUAD || - CTX::instance()->mesh.algo2d == ALGO_2D_BAMG) + CTX::instance()->mesh.algo2d == ALGO_2D_BAMG || + CTX::instance()->mesh.algo2d == ALGO_2D_FRONTAL || + CTX::instance()->mesh.algo2d == ALGO_2D_FRONTAL_QUAD || + CTX::instance()->mesh.algo2d == ALGO_2D_BAMG) return true; if(CTX::instance()->mesh.algo2d == ALGO_2D_AUTO && gf->geomType() == GEntity::Plane) @@ -312,7 +312,7 @@ static bool algoDelaunay2D(GFace *gf) return false; } -static void computeElementShapes(GFace *gf, double &worst, double &avg, +static void computeElementShapes(GFace *gf, double &worst, double &avg, double &best, int &nT, int &greaterThan) { worst = 1.e22; @@ -331,7 +331,7 @@ static void computeElementShapes(GFace *gf, double &worst, double &avg, avg /= nT; } -static bool recoverEdge(BDS_Mesh *m, GEdge *ge, +static bool recoverEdge(BDS_Mesh *m, GEdge *ge, std::map<MVertex*, BDS_Point*> &recoverMapInv, std::set<EdgeToRecover> *e2r, std::set<EdgeToRecover> *notRecovered, int pass) @@ -341,7 +341,7 @@ static bool recoverEdge(BDS_Mesh *m, GEdge *ge, 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++){ @@ -352,14 +352,14 @@ static bool recoverEdge(BDS_Mesh *m, GEdge *ge, if(itpstart != recoverMapInv.end() && itpend != recoverMapInv.end()){ BDS_Point *pstart = itpstart->second; BDS_Point *pend = itpend->second; - if(pass == 1) + 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 { if (_fatallyFailed) Msg::Error("Unable to recover an edge %g %g && %g %g (%d/%d)", - vstart->x(), vstart->y(), vend->x(), vend->y(), i, + vstart->x(), vstart->y(), vend->x(), vend->y(), i, ge->mesh_vertices.size()); return !_fatallyFailed; } @@ -369,7 +369,7 @@ static bool recoverEdge(BDS_Mesh *m, GEdge *ge, if(pass == 2 && ge->getBeginVertex()){ MVertex *vstart = *(ge->getBeginVertex()->mesh_vertices.begin()); - MVertex *vend = *(ge->getEndVertex()->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()){ @@ -454,7 +454,7 @@ void filterOverlappingElements(int dim, std::vector<MElement*> &e, for (int j=0;j<el->getNumVertices();++j){ MVertex *v = el->getVertex(j); std::vector<MElement *> inters = octree.findAll(v->x(),v->y(),v->z(),dim); - std::vector<MElement *> inters2; + std::vector<MElement *> inters2; for (int k=0;k<inters.size();k++){ bool found = false; for (int l=0;l<inters[k]->getNumVertices();l++){ @@ -468,7 +468,7 @@ void filterOverlappingElements(int dim, std::vector<MElement*> &e, } if (intersection){ printf("intersection found\n"); - einter.push_back(el); + einter.push_back(el); } else { eout.push_back(el); @@ -508,14 +508,14 @@ void modifyInitialMeshForTakingIntoAccountBoundaryLayers(GFace *gf) for (int l=0;l < N ;++l){ MVertex *v11,*v12,*v21,*v22; v21 = c1._column[l]; - v22 = c2._column[l]; + v22 = c2._column[l]; if (l == 0){ v11 = v1; v12 = v2; } else { v11 = c1._column[l-1]; - v12 = c2._column[l-1]; + v12 = c2._column[l-1]; } MEdge dv2 (v21,v22); @@ -539,7 +539,7 @@ void modifyInitialMeshForTakingIntoAccountBoundaryLayers(GFace *gf) 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); @@ -547,14 +547,14 @@ void modifyInitialMeshForTakingIntoAccountBoundaryLayers(GFace *gf) for (int l=0;l < N ;++l){ MVertex *v11,*v12,*v21,*v22; v21 = c1._column[l]; - v22 = c2._column[l]; + v22 = c2._column[l]; if (l == 0){ v11 = v; v12 = v; } else { v11 = c1._column[l-1]; - v12 = c2._column[l-1]; + v12 = c2._column[l-1]; } if (v11 != v12){ blQuads.push_back(new MQuadrangle(v11,v12,v22,v21)); @@ -577,27 +577,27 @@ void modifyInitialMeshForTakingIntoAccountBoundaryLayers(GFace *gf) fprintf(ff2,"};\n"); fclose(ff2); - + std::vector<MElement*> els,newels,oldels; for (int i=0;i<blQuads.size();i++)els.push_back(blQuads[i]); filterOverlappingElements (2,els,newels,oldels); blQuads.clear(); for (int i=0;i<newels.size();i++)blQuads.push_back((MQuadrangle*)newels[i]); for (int i=0;i<oldels.size();i++)delete oldels[i]; - - for (int i=0;i<blQuads.size();i++){ + + for (int i=0;i<blQuads.size();i++){ addOrRemove(blQuads[i]->getVertex(0),blQuads[i]->getVertex(1),bedges); addOrRemove(blQuads[i]->getVertex(1),blQuads[i]->getVertex(2),bedges); addOrRemove(blQuads[i]->getVertex(2),blQuads[i]->getVertex(3),bedges); addOrRemove(blQuads[i]->getVertex(3),blQuads[i]->getVertex(0),bedges); - for (int j=0;j<4;j++) + for (int j=0;j<4;j++) if(blQuads[i]->getVertex(j)->onWhat() == gf)verts.insert(blQuads[i]->getVertex(j)); } - for (int i=0;i<blTris.size();i++){ + for (int i=0;i<blTris.size();i++){ addOrRemove(blTris[i]->getVertex(0),blTris[i]->getVertex(1),bedges); addOrRemove(blTris[i]->getVertex(1),blTris[i]->getVertex(2),bedges); addOrRemove(blTris[i]->getVertex(2),blTris[i]->getVertex(0),bedges); - for (int j=0;j<3;j++) + for (int j=0;j<3;j++) if(blTris[i]->getVertex(j)->onWhat() == gf)verts.insert(blTris[i]->getVertex(j)); } @@ -621,8 +621,8 @@ void modifyInitialMeshForTakingIntoAccountBoundaryLayers(GFace *gf) deMeshGFace kil_; kil_(gf); - meshGenerator(gf, 0, 0, true , true, &hop); - + meshGenerator(gf, 0, 0, true , true, &hop); + gf->quadrangles = blQuads; gf->triangles.insert(gf->triangles.begin(),blTris.begin(),blTris.end()); gf->mesh_vertices.insert(gf->mesh_vertices.begin(),verts.begin(),verts.end()); @@ -631,7 +631,7 @@ void modifyInitialMeshForTakingIntoAccountBoundaryLayers(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 meshGenerator(GFace *gf, int RECUR_ITER, bool repairSelfIntersecting1dMesh, bool onlyInitialMesh, bool debug, @@ -667,7 +667,7 @@ bool meshGenerator(GFace *gf, int RECUR_ITER, } else Msg::Info("Degenerated mesh on edge %d", (*ite)->tag()); - ++ite; + ++ite; } std::list<GEdge*> emb_edges = gf->embeddedEdges(); @@ -675,7 +675,7 @@ bool meshGenerator(GFace *gf, int RECUR_ITER, 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(), @@ -683,7 +683,7 @@ bool meshGenerator(GFace *gf, int RECUR_ITER, } ++ite; } - + // add embedded vertices std::list<GVertex*> emb_vertx = gf->embeddedVertices(); std::list<GVertex*>::iterator itvx = emb_vertx.begin(); @@ -692,8 +692,8 @@ bool meshGenerator(GFace *gf, int RECUR_ITER, (*itvx)->mesh_vertices.end()); ++itvx; } - - // add additional vertices + + // add additional vertices all_vertices.insert(gf->additionalVertices.begin(), gf->additionalVertices.end()); @@ -715,7 +715,7 @@ bool meshGenerator(GFace *gf, int RECUR_ITER, std::vector<BDS_Point*> points(all_vertices.size()); SBoundingBox3d bbox; int count = 0; - for(std::set<MVertex*>::iterator it = all_vertices.begin(); + for(std::set<MVertex*>::iterator it = all_vertices.begin(); it != all_vertices.end(); it++){ MVertex *here = *it; GEntity *ge = here->onWhat(); @@ -745,9 +745,9 @@ bool meshGenerator(GFace *gf, int RECUR_ITER, DocRecord doc(points.size() + 4); { for(unsigned int i = 0; i < points.size(); i++){ - double XX = CTX::instance()->mesh.randFactor * LC2D * (double)rand() / + double XX = CTX::instance()->mesh.randFactor * LC2D * (double)rand() / (double)RAND_MAX; - double YY = CTX::instance()->mesh.randFactor * LC2D * (double)rand() / + double YY = CTX::instance()->mesh.randFactor * LC2D * (double)rand() / (double)RAND_MAX; doc.points[i].where.h = points[i]->u + XX; doc.points[i].where.v = points[i]->v + YY; @@ -755,7 +755,7 @@ bool meshGenerator(GFace *gf, int RECUR_ITER, doc.points[i].adjacent = NULL; } - + // increase the size of the bounding box bbox *= 2.5; @@ -775,7 +775,7 @@ bool meshGenerator(GFace *gf, int RECUR_ITER, doc.points[points.size() + ip].adjacent = 0; doc.points[points.size() + ip].data = pp; } - + // Use "fast" inhouse recursive algo to generate the triangulation. // At this stage the triangulation is not what we need // -) It does not necessary recover the boundaries @@ -784,7 +784,7 @@ bool meshGenerator(GFace *gf, int RECUR_ITER, Msg::Debug("Meshing of the convex hull (%d points)", points.size()); doc.MakeMeshWithPoints(); Msg::Debug("Meshing of the convex hull (%d points) done", points.size()); - + for(int i = 0; i < doc.numTriangles; i++) { BDS_Point *p1 = (BDS_Point*)doc.points[doc.triangles[i].a].data; BDS_Point *p2 = (BDS_Point*)doc.points[doc.triangles[i].b].data; @@ -818,8 +818,8 @@ bool meshGenerator(GFace *gf, int RECUR_ITER, recoverEdge(m, *ite, recoverMapInv, &edgesToRecover, &edgesNotRecovered, 1); ++ite; } - - + + // effectively recover the medge ite = edges.begin(); while(ite != edges.end()){ @@ -832,7 +832,7 @@ bool meshGenerator(GFace *gf, int RECUR_ITER, } ++ite; } - + Msg::Debug("Recovering %d mesh Edges (%d not recovered)", edgesToRecover.size(), edgesNotRecovered.size()); @@ -843,18 +843,18 @@ bool meshGenerator(GFace *gf, int RECUR_ITER, sstream << " " << itr->ge->tag(); Msg::Warning(":-( There are %d intersections in the 1D mesh (curves%s)", edgesNotRecovered.size(), sstream.str().c_str()); - if (repairSelfIntersecting1dMesh) + if (repairSelfIntersecting1dMesh) Msg::Warning("8-| Gmsh splits those edges and tries again"); - + if(debug){ char name[245]; - sprintf(name, "surface%d-not_yet_recovered-real-%d.msh", gf->tag(), + sprintf(name, "surface%d-not_yet_recovered-real-%d.msh", gf->tag(), RECUR_ITER); gf->model()->writeMSH(name); } - + std::list<GFace *> facesToRemesh; - if(repairSelfIntersecting1dMesh) + if(repairSelfIntersecting1dMesh) remeshUnrecoveredEdges(recoverMapInv, edgesNotRecovered, facesToRemesh); else{ std::set<EdgeToRecover>::iterator itr = edgesNotRecovered.begin(); @@ -884,10 +884,10 @@ bool meshGenerator(GFace *gf, int RECUR_ITER, if(RECUR_ITER > 0) Msg::Warning(":-) Gmsh was able to recover all edges after %d iterations", - RECUR_ITER); - + RECUR_ITER); + Msg::Debug("Boundary Edges recovered for surface %d", gf->tag()); - + // look for a triangle that has a negative node and recursively // tag all exterior triangles { @@ -901,7 +901,7 @@ bool meshGenerator(GFace *gf, int RECUR_ITER, BDS_Face *t = *itt; BDS_Point *n[4]; t->getNodes(n); - if (n[0]->iD < 0 || n[1]->iD < 0 || + if (n[0]->iD < 0 || n[1]->iD < 0 || n[2]->iD < 0 ) { recur_tag(t, &CLASS_EXTERIOR); break; @@ -909,7 +909,7 @@ bool meshGenerator(GFace *gf, int RECUR_ITER, ++itt; } } - + // now find an edge that has belongs to one of the exterior // triangles { @@ -954,7 +954,7 @@ bool meshGenerator(GFace *gf, int RECUR_ITER, ++ite; } } - + ite = emb_edges.begin(); while(ite != emb_edges.end()){ if(!(*ite)->isMeshDegenerated()) @@ -962,9 +962,9 @@ bool meshGenerator(GFace *gf, int RECUR_ITER, ++ite; } - // compute characteristic lengths at vertices + // compute characteristic lengths at vertices if (!onlyInitialMesh){ - Msg::Debug("Computing mesh size field at mesh vertices %d", + Msg::Debug("Computing mesh size field at mesh vertices %d", edgesToRecover.size()); for(int i = 0; i < doc.numPoints; i++){ BDS_Point *pp = (BDS_Point*)doc.points[i].data; @@ -981,7 +981,7 @@ bool meshGenerator(GFace *gf, int RECUR_ITER, pp->lcBGM() = BGM_MeshSize(ge, u, 0, here->x(), here->y(), here->z()); } else - pp->lcBGM() = MAX_LC; + pp->lcBGM() = MAX_LC; pp->lc() = pp->lcBGM(); } } @@ -1004,11 +1004,11 @@ bool meshGenerator(GFace *gf, int RECUR_ITER, if(e->numfaces() == 0) m->del_edge(e); else{ - if(!e->g) + if(!e->g) e->g = &CLASS_F; - if(!e->p1->g || e->p1->g->classif_degree > e->g->classif_degree) + 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) + if(!e->p2->g || e->p2->g->classif_degree > e->g->classif_degree) e->p2->g = e->g; } ++ite; @@ -1027,7 +1027,7 @@ bool meshGenerator(GFace *gf, int RECUR_ITER, sprintf(name, "surface%d-recovered-param.pos", gf->tag()); outputScalarField(m->triangles, name, 1); } - + if(1){ std::list<BDS_Face*>::iterator itt = m->triangles.begin(); while (itt != m->triangles.end()){ @@ -1050,7 +1050,7 @@ bool meshGenerator(GFace *gf, int RECUR_ITER, ++itt; } } - + if (Msg::GetVerbosity() == 10){ GEdge *ge = new discreteEdge(gf->model(), 1000, 0, 0); MElementOctree octree(gf->model()); @@ -1078,7 +1078,7 @@ bool meshGenerator(GFace *gf, int RECUR_ITER, // printf("coucou here !!!\n"); // backgroundMesh::unset(); // buildBackGroundMesh (gf); - // } + // } refineMeshBDS(gf, *m, CTX::instance()->mesh.refineSteps, true, &recoverMapInv); optimizeMeshBDS(gf, *m, 2); @@ -1087,7 +1087,7 @@ bool meshGenerator(GFace *gf, int RECUR_ITER, optimizeMeshBDS(gf, *m, 2); // if(CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine || 1) { // backgroundMesh::unset(); - // } + // } } /* computeMeshSizeFieldAccuracy(gf, *m, gf->meshStatistics.efficiency_index, @@ -1103,7 +1103,7 @@ bool meshGenerator(GFace *gf, int RECUR_ITER, // BOUNDARY LAYER if (!onlyInitialMesh) modifyInitialMeshForTakingIntoAccountBoundaryLayers(gf); - + // the delaunay algo is based directly on internal gmsh structures // BDS mesh is passed in order not to recompute local coordinates of // vertices @@ -1131,12 +1131,12 @@ bool meshGenerator(GFace *gf, int RECUR_ITER, outputScalarField(m->triangles, name,1); } if(CTX::instance()->mesh.remove4triangles) - removeFourTrianglesNodes(gf,false); + removeFourTrianglesNodes(gf,false); // delete the mesh delete m; - if((CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine) && + if((CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine) && !CTX::instance()->mesh.optimizeLloyd && !onlyInitialMesh) recombineIntoQuadsIterative(gf); @@ -1197,21 +1197,21 @@ static bool buildConsecutiveListOfVertices(GFace *gf, GEdgeLoop &gel, GEdgeLoop::iter it = gel.begin(); - if(MYDEBUG) - printf("face %d with %d edges case %d\n", gf->tag(), + if(MYDEBUG) + 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 ; 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()); - + 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)); @@ -1226,7 +1226,7 @@ static bool buildConsecutiveListOfVertices(GFace *gf, GEdgeLoop &gel, 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> > + 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); @@ -1292,7 +1292,7 @@ static bool buildConsecutiveListOfVertices(GFace *gf, GEdgeLoop &gel, 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", + printf("%g %g dist_reversed = %12.5E\n", first_coord_reversed.x(), first_coord_reversed.y(), d_reversed); if(d_reversed < tol){ coords.clear(); @@ -1339,7 +1339,7 @@ static bool buildConsecutiveListOfVertices(GFace *gf, GEdgeLoop &gel, } return false; } - + std::vector<MVertex*> edgeLoop; if(found._sign == 1){ edgeLoop.push_back(found.ge->getBeginVertex()->mesh_vertices[0]); @@ -1351,11 +1351,11 @@ static bool buildConsecutiveListOfVertices(GFace *gf, GEdgeLoop &gel, for(int i = found.ge->mesh_vertices.size() - 1; i >= 0; i--) edgeLoop.push_back(found.ge->mesh_vertices[i]); } - + if(MYDEBUG) 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++){ MVertex *here = edgeLoop[i]; @@ -1375,7 +1375,7 @@ static bool buildConsecutiveListOfVertices(GFace *gf, GEdgeLoop &gel, } else pp->lcBGM() = MAX_LC; - + pp->lc() = pp->lcBGM(); m->add_geom (ge->tag(), ge->dim()); BDS_GeomEntity *g = m->get_geom(ge->tag(), ge->dim()); @@ -1393,7 +1393,7 @@ static bool buildConsecutiveListOfVertices(GFace *gf, GEdgeLoop &gel, if(MYDEBUG) printf("last coord %g %g\n", last_coord.x(), last_coord.y()); result.insert(result.end(), edgeLoop_BDS.begin(), edgeLoop_BDS.end()); } - + // It has worked, so we add all the points to the recover map recoverMap.insert(recoverMapLocal.begin(), recoverMapLocal.end()); @@ -1422,20 +1422,20 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true) SBoundingBox3d bbox; int nbPointsTotal = 0; { - for(std::list<GEdgeLoop>::iterator it = gf->edgeLoops.begin(); + for(std::list<GEdgeLoop>::iterator it = gf->edgeLoops.begin(); 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++){ - if(buildConsecutiveListOfVertices(gf, *it, edgeLoop_BDS, bbox, m, - recoverMap, nbPointsLocal, + if(buildConsecutiveListOfVertices(gf, *it, edgeLoop_BDS, bbox, m, + recoverMap, nbPointsLocal, nbPointsTotal, fact[i] * LC2D)){ ok = true; break; } - if(buildConsecutiveListOfVertices(gf, *it, edgeLoop_BDS, bbox, m, + if(buildConsecutiveListOfVertices(gf, *it, edgeLoop_BDS, bbox, m, recoverMap, nbPointsLocal, nbPointsTotal, fact[i] * LC2D, true)){ ok = true; @@ -1463,9 +1463,9 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true) 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 XX = CTX::instance()->mesh.randFactor * LC2D * (double)rand() / (double)RAND_MAX; - double YY = CTX::instance()->mesh.randFactor * LC2D * (double)rand() / + 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; @@ -1483,7 +1483,7 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true) bb[0] = new MVertex(bbox.min().x(), bbox.min().y(), 0, 0, -1); 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); + bb[3] = new MVertex(bbox.max().x(), bbox.max().y(), 0, 0, -4); 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); @@ -1495,7 +1495,7 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true) doc.points[nbPointsTotal+ip].data = pp; } for(int ip = 0; ip < 4; ip++) delete bb[ip]; - + // Use "fast" inhouse recursive algo to generate the triangulation // At this stage the triangulation is not what we need // -) It does not necessary recover the boundaries @@ -1503,7 +1503,7 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true) // loop is the outer one) Msg::Debug("Meshing of the convex hull (%d points)", nbPointsTotal); doc.MakeMeshWithPoints(); - + for(int i = 0; i < doc.numTriangles; i++){ BDS_Point *p1 = (BDS_Point*)doc.points[doc.triangles[i].a].data; BDS_Point *p2 = (BDS_Point*)doc.points[doc.triangles[i].b].data; @@ -1556,7 +1556,7 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true) BDS_Face *t = *itt; BDS_Point *n[4]; t->getNodes(n); - if (n[0]->iD < 0 || n[1]->iD < 0 || + if (n[0]->iD < 0 || n[1]->iD < 0 || n[2]->iD < 0 ) { recur_tag(t, &CLASS_EXTERIOR); break; @@ -1564,7 +1564,7 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true) ++itt; } } - + // now find an edge that has belongs to one of the exterior // triangles { @@ -1601,9 +1601,9 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true) ++itt; } } - + m->cleanup(); - + { std::list<BDS_Edge*>::iterator ite = m->edges.begin(); while (ite != m->edges.end()){ @@ -1642,7 +1642,7 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true) // printf("coucou here !!!\n"); // backgroundMesh::unset(); // buildBackGroundMesh (gf); - // } + // } refineMeshBDS(gf, *m, CTX::instance()->mesh.refineSteps, true); optimizeMeshBDS(gf, *m, 2); refineMeshBDS(gf, *m, -CTX::instance()->mesh.refineSteps, false); @@ -1657,9 +1657,9 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true) gf->meshStatistics.status = GFace::DONE; // if(CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine || 1) { // backgroundMesh::unset(); - // } + // } } - + // fill the small gmsh structures { std::set<BDS_Point*, PointLessThan>::iterator itp = m->points.begin(); @@ -1674,7 +1674,7 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true) ++itp; } } - + { std::list<BDS_Face*>::iterator itt = m->triangles.begin(); while (itt != m->triangles.end()){ @@ -1700,7 +1700,7 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true) ++itt; } } - + if(debug){ char name[245]; sprintf(name, "surface%d-final-real.pos", gf->tag()); @@ -1708,7 +1708,7 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true) sprintf(name, "surface%d-final-param.pos", gf->tag()); outputScalarField(m->triangles, name, 1); } - + if(algoDelaunay2D(gf)){ if(CTX::instance()->mesh.algo2d == ALGO_2D_FRONTAL) bowyerWatsonFrontal(gf); @@ -1717,18 +1717,18 @@ static bool meshGeneratorPeriodic(GFace *gf, bool debug = true) else if(CTX::instance()->mesh.algo2d == ALGO_2D_DELAUNAY || CTX::instance()->mesh.algo2d == ALGO_2D_AUTO) bowyerWatson(gf); - else + else meshGFaceBamg(gf); laplaceSmoothing(gf,CTX::instance()->mesh.nbSmoothing); } - - // delete the mesh + + // delete the mesh delete m; - if((CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine) && + if((CTX::instance()->mesh.recombineAll || gf->meshAttributes.recombine) && !CTX::instance()->mesh.optimizeLloyd) recombineIntoQuads(gf); - + computeElementShapes(gf, gf->meshStatistics.worst_element_shape, gf->meshStatistics.average_element_shape, gf->meshStatistics.best_element_shape, @@ -1746,7 +1746,7 @@ void deMeshGFace::operator() (GFace *gf) } // for debugging, change value from -1 to -100; -int debugSurface = -1; //-1; +int debugSurface = -1; //-1; void meshGFace::operator() (GFace *gf) { @@ -1775,7 +1775,7 @@ void meshGFace::operator() (GFace *gf) gf->meshStatistics.status = GFace::PENDING; return; } - Msg::Info("Meshing face %d (%s) as a copy of %d", gf->tag(), + Msg::Info("Meshing face %d (%s) as a copy of %d", gf->tag(), gf->getTypeString().c_str(), gf->meshMaster()); copyMesh(gff, gf); gf->meshStatistics.status = GFace::DONE; @@ -1812,7 +1812,7 @@ void meshGFace::operator() (GFace *gf) Msg::Debug("Generating the mesh"); if ((gf->getNativeType() != GEntity::AcisModel || (!gf->periodic(0) && !gf->periodic(1))) && - (noSeam(gf) || gf->getNativeType() == GEntity::GmshModel || + (noSeam(gf) || gf->getNativeType() == GEntity::GmshModel || gf->edgeLoops.empty())){ meshGenerator(gf, 0, repairSelfIntersecting1dMesh, onlyInitialMesh, debugSurface >= 0 || debugSurface == -100); @@ -1822,7 +1822,7 @@ void meshGFace::operator() (GFace *gf) (gf, debugSurface >= 0 || debugSurface == -100)) Msg::Error("Impossible to mesh periodic face %d", gf->tag()); } - + Msg::Debug("Type %d %d triangles generated, %d internal vertices", gf->geomType(), gf->triangles.size(), gf->mesh_vertices.size()); @@ -1834,7 +1834,7 @@ void meshGFace::operator() (GFace *gf) twoPassesMesh--; if (backgroundMesh::current()){ backgroundMesh::unset(); - } + } if (CTX::instance()->mesh.saveAll){ backgroundMesh::set(gf); char name[256]; @@ -1849,22 +1849,22 @@ void meshGFace::operator() (GFace *gf) bool checkMeshCompound(GFaceCompound *gf, std::list<GEdge*> &edges) { bool isMeshed = false; -#if defined(HAVE_SOLVER) +#if defined(HAVE_SOLVER) bool correctTopo = gf->checkTopology(); if (!correctTopo && gf->allowPartition()){ partitionAndRemesh((GFaceCompound*) gf); isMeshed = true; return isMeshed; } - + bool correctParam = gf->parametrize(); - + if (!correctParam && gf->allowPartition()){ partitionAndRemesh((GFaceCompound*) gf); isMeshed = true; return isMeshed; } - + //Replace edges by their compounds std::set<GEdge*> mySet; std::list<GEdge*>::iterator it = edges.begin(); @@ -1872,7 +1872,7 @@ bool checkMeshCompound(GFaceCompound *gf, std::list<GEdge*> &edges) if((*it)->getCompound()){ mySet.insert((*it)->getCompound()); } - else{ + else{ mySet.insert(*it); } ++it; @@ -1898,9 +1898,9 @@ void partitionAndRemesh(GFaceCompound *gf) typeOfPartition method; if(gf->nbSplit > 0) method = MULTILEVEL; else method = LAPLACIAN; - + int allowType = gf->allowPartition(); - multiscalePartition *msp = new multiscalePartition(elements, abs(gf->nbSplit), + multiscalePartition *msp = new multiscalePartition(elements, abs(gf->nbSplit), method, allowType); int NF = msp->getNumberOfParts(); @@ -1908,15 +1908,15 @@ void partitionAndRemesh(GFaceCompound *gf) int nume = gf->model()->getMaxElementaryNumber(1) + 1; int numf = gf->model()->getMaxElementaryNumber(2) + 1; std::vector<discreteFace*> pFaces; - createPartitionFaces(gf->model(), elements, NF, pFaces); - + createPartitionFaces(gf->model(), elements, NF, pFaces); + gf->model()->createTopologyFromFaces(pFaces); - + double tmult = Cpu(); - Msg::Info("Multiscale Partition SUCCESSFULLY PERFORMED : %d parts (%g s)", + Msg::Info("Multiscale Partition SUCCESSFULLY PERFORMED : %d parts (%g s)", NF, tmult - tbegin); gf->model()->writeMSH("multiscalePARTS.msh", 2.2, false, true); - + // Remesh new faces (Compound Lines and Compound Surfaces) Msg::Info("*** Starting parametrize compounds:"); double t0 = Cpu(); @@ -1926,9 +1926,9 @@ void partitionAndRemesh(GFaceCompound *gf) for (int i=0; i < NE; i++){ std::vector<GEdge*>e_compound; GEdge *pe = gf->model()->getEdgeByTag(nume+i);//partition edge - e_compound.push_back(pe); + e_compound.push_back(pe); int num_gec = nume + NE + i ; - Msg::Info("Parametrize Compound Line (%d) = %d discrete edge", + Msg::Info("Parametrize Compound Line (%d) = %d discrete edge", num_gec, pe->tag()); GEdgeCompound *gec = new GEdgeCompound(gf->model(), num_gec, e_compound); gf->model()->add(gec); @@ -1937,16 +1937,16 @@ void partitionAndRemesh(GFaceCompound *gf) } // Parametrize Compound surfaces - std::set<MVertex*> allNod; + std::set<MVertex*> allNod; std::list<GEdge*> U0; for (int i=0; i < NF; i++){ std::list<GFace*> f_compound; - GFace *pf = gf->model()->getFaceByTag(numf+i);//partition face + GFace *pf = gf->model()->getFaceByTag(numf+i);//partition face int num_gfc = numf + NF + i ; - f_compound.push_back(pf); + f_compound.push_back(pf); Msg::Info("Parametrize Compound Surface (%d) = %d discrete face", num_gfc, pf->tag()); - + GFaceCompound *gfc = new GFaceCompound(gf->model(), num_gfc, f_compound, U0, gf->getTypeOfMapping()); @@ -1958,7 +1958,7 @@ void partitionAndRemesh(GFaceCompound *gf) double t1 = Cpu(); Msg::Info("*** Parametrize compounds done (%g s)", t1-t0); - + Msg::Info("*** Starting meshing 1D edges ...:"); for (int i = 0; i < NE; i++){ GEdge *gec = gf->model()->getEdgeByTag(nume + NE + i); @@ -1981,21 +1981,21 @@ void partitionAndRemesh(GFaceCompound *gf) MTriangle *t = gfc->triangles[j]; std::vector<MVertex *> v(3); for(int k = 0; k < 3; k++){ - v[k] = t->getVertex(k); + v[k] = t->getVertex(k); allNod.insert(v[k]); } gf->triangles.push_back(new MTriangle(v[0], v[1], v[2])); - } + } for(unsigned int j = 0; j < gfc->quadrangles.size(); ++j){ MQuadrangle *q = gfc->quadrangles[j]; std::vector<MVertex *> v(4); for(int k = 0; k < 4; k++){ - v[k] = q->getVertex(k); + v[k] = q->getVertex(k); allNod.insert(v[k]); } gf->quadrangles.push_back(new MQuadrangle(v[0], v[1], v[2], v[3])); - } - + } + } // Removing discrete Vertices - Edges - Faces @@ -2007,13 +2007,13 @@ void partitionAndRemesh(GFaceCompound *gf) for (int i=0; i < NE; i++){ GEdge *gec = gf->model()->getEdgeByTag(nume+NE+i); GEdge *pe = gf->model()->getEdgeByTag(nume+i); - gf->model()->remove(pe); + gf->model()->remove(pe); gf->model()->remove(gec); } for (int i=0; i < NF; i++){ GFace *gfc = gf->model()->getFaceByTag(numf+NF+i); GFace *pf = gf->model()->getFaceByTag(numf+i); - gf->model()->remove(pf); + gf->model()->remove(pf); gf->model()->remove(gfc); } @@ -2066,9 +2066,9 @@ void partitionAndRemesh(GFaceCompound *gf) Msg::Info("*** Mesh of surface %d done by assembly %d remeshed faces (%g s)", gf->tag(), NF, t3-t2); Msg::Info("-----------------------------------------------------------"); - + gf->coherenceNormals(); - gf->meshStatistics.status = GFace::DONE; + gf->meshStatistics.status = GFace::DONE; #endif } @@ -2090,7 +2090,7 @@ void orientMeshGFace::operator()(GFace *gf) // In old versions we checked the orientation by comparing the // orientation of a line element on the boundary w.r.t. its // connected surface element. This is probably better than what - // follows, but + // follows, but // * it failed when the 3D Delaunay changes the 1D mesh (since we // don't recover it yet) // * it failed with OpenCASCADE geometries, where surface orientions @@ -2107,7 +2107,7 @@ void orientMeshGFace::operator()(GFace *gf) SPoint2 param; if(v->onWhat() == gf && v->getParameter(0, param[0]) && v->getParameter(1, param[1])){ - SVector3 nf = gf->normal(param); + SVector3 nf = gf->normal(param); SVector3 ne = e->getFace(0).normal(); if(dot(ne, nf) < 0){ Msg::Debug("Reverting orientation of mesh in face %d", gf->tag()); @@ -2118,7 +2118,7 @@ void orientMeshGFace::operator()(GFace *gf) } } } - + // if we could not find such an element, just try to evaluate the // normal at the barycenter of an element on the surface for(unsigned int i = 0; i < gf->getNumMeshElements(); i++){ @@ -2135,7 +2135,7 @@ void orientMeshGFace::operator()(GFace *gf) } if(ok){ param *= 1. / e->getNumVertices(); - SVector3 nf = gf->normal(param); + SVector3 nf = gf->normal(param); SVector3 ne = e->getFace(0).normal(); if(dot(ne, nf) < 0){ Msg::Debug("Reverting orientation of mesh in face %d", gf->tag()); diff --git a/benchmarks/stl/mobilette.geo b/benchmarks/stl/mobilette.geo index fda2e1b8306970477c5bad0cc9f3348f351e20d5..1a24ac3ad017865d94b88fa37a031fce85d6b824 100644 --- a/benchmarks/stl/mobilette.geo +++ b/benchmarks/stl/mobilette.geo @@ -1,10 +1,11 @@ -Mesh.Algorithm = 8; //(1=MeshAdapt, 2=Automatic, 5=Delaunay, 6=Frontal, 7=bamg) +Mesh.Algorithm = 2; //(1=MeshAdapt, 2=Automatic, 5=Delaunay, 6=Frontal, 7=bamg) Mesh.CharacteristicLengthMin=1.5; Mesh.CharacteristicLengthMax=2.5; Mesh.RemeshAlgorithm=1; -Mesh.RemeshParametrization=1;//(0) harmonic (1) conformal +Mesh.RemeshParametrization=1;//(0) harmonic (1) conformal Mesh.RecombinationAlgorithm=1; -Mesh.RecombineAll=1; +//Mesh.RecombineAll=1; + // merge reclassified STL Merge "mobilette-class.msh"; @@ -27,8 +28,5 @@ EndFor Surface Loop(1) = {s : s + #ss[]-1}; Volume(1) = {1}; -Mesh.RecombineAll=1; -Mesh.RecombinationAlgorithm=1; - -//Physical Surface(1) = {s : s + #ss[]-1}; -//Physical Volume(1) = 1; +Physical Surface(1) = {s : s + #ss[]-1}; +Physical Volume(1) = 1;