From 1ed40a57a7536e01c20f2d9ed19996474d99866c Mon Sep 17 00:00:00 2001 From: ws Date: Mon, 27 Aug 2018 16:56:42 +0200 Subject: [PATCH 1/3] Improvement in 2D frontal meshing: in case of a circular surface, we might have all elements with a too large radius. The frontal algorithm will then start from the elements next to the boundary to see whether these could be refined by adding a point. If only 2 elements are then in the cavity, this point is only allowed to add when the quality improves. However, in case of a circular plane surface with a uniform coarseness, it might happen that all elements along the boundary are part of cavities of only 2 elements and adding points is worsen the element quality and thus this node is not added. The radius of the element is then set to -1. If no nodes are added then, no element at the inside will be a candidate for refinement as they are all connected to elements with a too high radius or too low (-1). To avoid this, we should set the radius of the elements which were not refined due to worsen element quality to a value which indicates that the element is okay (so not -1, but a value between 0 and LIMIT_). Now, the internal triangles are candidates for refinement as well, and the mesh looks much more uniform, as it should be. --- Mesh/meshGFaceDelaunayInsertion.cpp | 27 ++++++++++++++++++++++++++- 1 file changed, 26 insertions(+), 1 deletion(-) diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp index 00893c820..d756fd8b9 100644 --- a/Mesh/meshGFaceDelaunayInsertion.cpp +++ b/Mesh/meshGFaceDelaunayInsertion.cpp @@ -1019,6 +1019,7 @@ int insertVertexB(std::list &shell, std::list &cavity, // DT_INSERT_VERTEX += t2-t1; if(std::abs(oldVolume - newVolume) > EPS * oldVolume) return -3; if(onePointIsTooClose) return -4; + if(cavity.size() == 2) return -6; return -5; } } @@ -1247,9 +1248,33 @@ static bool insertAPoint(GFace *gf, Msg::Debug("Point %g %g cannot be inserted because it is out of the " "parametric domain)", center[0], center[1]); + if (result == -6) + Msg::Debug("Point %g %g cannot be inserted because it is not improving " + "the quality of the 2 elements", + center[0], center[1]); AllTris.erase(it); - worst->forceRadius(-1); + + // As in case of result = -6 it might happen that inside we would like to refine the mesh, + // we should not force the radius to -1, but pretend that the radius is good enough (less than LIMIT_ but higher than 0). + if (result == -6){ + Msg::Info("Forcing radius to %g\n", 0.5 * LIMIT_); + worst->forceRadius(0.5 * LIMIT_); + + // Go over its neighbours to check whether they should become active + for (size_t index = 0; index < 3; index++){ + const auto pNeighbour = worst->getNeigh(index); + if (pNeighbour == 0) continue; + int active_edge; + if (isActive(pNeighbour, LIMIT_, active_edge) && pNeighbour->getRadius() > LIMIT_){ + if ((*ActiveTris).find(pNeighbour) == (*ActiveTris).end()) + (*ActiveTris).insert(pNeighbour); + } + } + } + else + worst->forceRadius(-1); + AllTris.insert(worst); delete v; for(std::list::iterator itc = cavity.begin(); -- GitLab From 3bc16d57c7e4dc6862a0e427a75b88373f5f11e7 Mon Sep 17 00:00:00 2001 From: ws Date: Mon, 27 Aug 2018 17:27:31 +0200 Subject: [PATCH 2/3] Do not use "auto" --- Mesh/meshGFaceDelaunayInsertion.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp index d756fd8b9..a4fb98086 100644 --- a/Mesh/meshGFaceDelaunayInsertion.cpp +++ b/Mesh/meshGFaceDelaunayInsertion.cpp @@ -1263,7 +1263,7 @@ static bool insertAPoint(GFace *gf, // Go over its neighbours to check whether they should become active for (size_t index = 0; index < 3; index++){ - const auto pNeighbour = worst->getNeigh(index); + MTri3 *pNeighbour = worst->getNeigh(index); // TODO C++11 use auto if (pNeighbour == 0) continue; int active_edge; if (isActive(pNeighbour, LIMIT_, active_edge) && pNeighbour->getRadius() > LIMIT_){ -- GitLab From f3ba150c89c833c21d06794ffb18abd45ab16c4e Mon Sep 17 00:00:00 2001 From: Christophe Geuzaine Date: Tue, 28 Aug 2018 09:31:08 +0200 Subject: [PATCH 3/3] fix crash --- .gitlab-ci.yml | 1 + Mesh/meshGFaceDelaunayInsertion.cpp | 37 ++++++++++++++++------------- 2 files changed, 21 insertions(+), 17 deletions(-) diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml index 3aa470f8e..e86ebb510 100644 --- a/.gitlab-ci.yml +++ b/.gitlab-ci.yml @@ -63,6 +63,7 @@ windows64_ci: - cd build - bash -c "/usr/bin/cmake -DCMAKE_PREFIX_PATH='/usr/local/opencascade;/usr/local;/usr/x86_64-w64-mingw32/sys-root/mingw' -DCMAKE_C_COMPILER=/usr/bin/x86_64-w64-mingw32-gcc.exe -DCMAKE_CXX_COMPILER=/usr/bin/x86_64-w64-mingw32-g++.exe -DCMAKE_Fortran_COMPILER=/usr/bin/x86_64-w64-mingw32-gfortran.exe -DCMAKE_RC_COMPILER=/usr/bin/x86_64-w64-mingw32-windres.exe -DPETSC_ARCH=complex_mumps_seq -DPETSC_DIR=/home/Administrateur/petsc -DSLEPC_DIR=/home/Administrateur/slepc -DENABLE_OS_SPECIFIC_INSTALL=1 ${EXTRA_OPTION} .." - bash -c "/usr/bin/make -j 8" + - bash -c "ctest -j 4 --output-on-failure" tags: - windows64 - shared diff --git a/Mesh/meshGFaceDelaunayInsertion.cpp b/Mesh/meshGFaceDelaunayInsertion.cpp index a4fb98086..ef0b3c55e 100644 --- a/Mesh/meshGFaceDelaunayInsertion.cpp +++ b/Mesh/meshGFaceDelaunayInsertion.cpp @@ -130,7 +130,7 @@ void _printTris(char *name, ITERATOR it, ITERATOR end, bidimMeshData *data, GFac (worst)->tri()->getVertex(0)->getNum(), (worst)->tri()->getVertex(0)->getNum(), (worst)->tri()->getVertex(1)->getNum(), - (worst)->tri()->getVertex(2)->getNum()); + (worst)->tri()->getVertex(2)->getNum()); } else if (!deg[0] && deg[1] && !deg[2]){ fprintf(ff, "SQ(%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g) {%d,%d,%d,%d};\n", @@ -138,7 +138,7 @@ void _printTris(char *name, ITERATOR it, ITERATOR end, bidimMeshData *data, GFac (worst)->tri()->getVertex(1)->getNum(), (worst)->tri()->getVertex(1)->getNum(), (worst)->tri()->getVertex(2)->getNum(), - (worst)->tri()->getVertex(0)->getNum()); + (worst)->tri()->getVertex(0)->getNum()); } else if (!deg[0] && !deg[1] && deg[2]){ fprintf(ff, "SQ(%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g,%g) {%d,%d,%d,%d};\n", @@ -146,7 +146,7 @@ void _printTris(char *name, ITERATOR it, ITERATOR end, bidimMeshData *data, GFac (worst)->tri()->getVertex(2)->getNum(), (worst)->tri()->getVertex(2)->getNum(), (worst)->tri()->getVertex(0)->getNum(), - (worst)->tri()->getVertex(1)->getNum()); + (worst)->tri()->getVertex(1)->getNum()); } else if (!deg[0] && !deg[1] && !deg[2]){ fprintf(ff, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g) {%d,%d,%d};\n", @@ -155,7 +155,7 @@ void _printTris(char *name, ITERATOR it, ITERATOR end, bidimMeshData *data, GFac (worst)->tri()->getVertex(1)->getNum(), (worst)->tri()->getVertex(2)->getNum()); } - } + } else { fprintf(ff, "ST(%g,%g,%g,%g,%g,%g,%g,%g,%g) {%d,%d,%d};\n", u1,v1,0.,u2,v2,0.,u3,v3,0., @@ -1255,20 +1255,23 @@ static bool insertAPoint(GFace *gf, AllTris.erase(it); - // As in case of result = -6 it might happen that inside we would like to refine the mesh, - // we should not force the radius to -1, but pretend that the radius is good enough (less than LIMIT_ but higher than 0). + // As in case of result = -6 it might happen that inside we would like to + // refine the mesh, we should not force the radius to -1, but pretend that + // the radius is good enough (less than LIMIT_ but higher than 0). if (result == -6){ - Msg::Info("Forcing radius to %g\n", 0.5 * LIMIT_); + Msg::Debug("Forcing radius to %g", 0.5 * LIMIT_); worst->forceRadius(0.5 * LIMIT_); - - // Go over its neighbours to check whether they should become active - for (size_t index = 0; index < 3; index++){ - MTri3 *pNeighbour = worst->getNeigh(index); // TODO C++11 use auto - if (pNeighbour == 0) continue; - int active_edge; - if (isActive(pNeighbour, LIMIT_, active_edge) && pNeighbour->getRadius() > LIMIT_){ - if ((*ActiveTris).find(pNeighbour) == (*ActiveTris).end()) - (*ActiveTris).insert(pNeighbour); + if(ActiveTris){ + // Go over its neighbours to check whether they should become active + for (size_t index = 0; index < 3; index++){ + MTri3 *pNeighbour = worst->getNeigh(index); // TODO C++11 use auto + if (pNeighbour == 0) continue; + int active_edge; + if(isActive(pNeighbour, LIMIT_, active_edge) && + pNeighbour->getRadius() > LIMIT_){ + if((*ActiveTris).find(pNeighbour) == (*ActiveTris).end()) + (*ActiveTris).insert(pNeighbour); + } } } } @@ -1597,7 +1600,7 @@ void bowyerWatsonFrontal(GFace *gf, std::map *equivalence, std::set degenerated; getDegeneratedVertices (gf,degenerated); - + buildMeshGenerationDataStructures(gf, AllTris, DATA); // delaunise the initial mesh -- GitLab