From 3cfe8b5bc00dbe9eb6e83a1bb791abcdc413a683 Mon Sep 17 00:00:00 2001
From: Thomas Toulorge <thomas.toulorge@mines-paristech.fr>
Date: Mon, 20 Mar 2017 10:38:20 +0000
Subject: [PATCH] Improved robustness of fast curving of high-order boundary
layer mesh
---
contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp | 11 ++++++++++-
1 file changed, 10 insertions(+), 1 deletion(-)
diff --git a/contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp b/contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp
index 9366e9a960..d620a3e1c6 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp
@@ -218,6 +218,7 @@ void getColumnQuad(MEdgeVecMEltMap &ed2el, const FastCurvingParameters &p,
for (int iLayer = 0; iLayer < p.maxNumLayers; iLayer++) {
std::vector<MElement*> newElts = ed2el[elBaseEd];
+ if ((iLayer > 0) && (newElts.size() < 2)) { aboveElt = 0; break; }
el = (newElts[0] == el) ? newElts[1] : newElts[0];
aboveElt = el;
if (el->getType() != TYPE_QUA) break;
@@ -248,6 +249,7 @@ void getColumnTri(MEdgeVecMEltMap &ed2el, const FastCurvingParameters &p,
for (int iLayer = 0; iLayer < p.maxNumLayers; iLayer++) {
// Get first element in layer
std::vector<MElement*> newElts0 = ed2el[elBaseEd];
+ if ((iLayer > 0) && (newElts0.size() < 2)) { aboveElt = 0; break; }
el0 = (newElts0[0] == el1) ? newElts0[1] : newElts0[0];
aboveElt = el0;
if (el0->getType() != TYPE_TRI) break;
@@ -260,6 +262,7 @@ void getColumnTri(MEdgeVecMEltMap &ed2el, const FastCurvingParameters &p,
// Get second element in layer
std::vector<MElement*> newElts1 = ed2el[elTopEd0];
+ if (newElts1.size() < 2) { aboveElt = 0; break; }
el1 = (newElts1[0] == el0) ? newElts1[1] : newElts1[0];
if (el1->getType() != TYPE_TRI) break;
MEdge elTopEd1;
@@ -449,6 +452,7 @@ void getColumnTet(MFaceVecMEltMap &face2el,
for (int iLayer = 0; iLayer < p.maxNumLayers; iLayer++) {
// Get first element in layer
std::vector<MElement*> newElts0 = face2el[elBaseFace];
+ if ((iLayer > 0) && (newElts0.size() < 2)) { aboveElt = 0; break; }
el0 = (newElts0[0] == el2) ? newElts0[1] : newElts0[0];
aboveElt = el0;
if (el0->getType() != TYPE_TET) break;
@@ -461,6 +465,7 @@ void getColumnTet(MFaceVecMEltMap &face2el,
// Get second element in layer
std::vector<MElement*> newElts1 = face2el[elTopFace0];
+ if (newElts1.size() < 2) { aboveElt = 0; break; }
el1 = (newElts1[0] == el0) ? newElts1[1] : newElts1[0];
if (el1->getType() != TYPE_TET) break;
MFace elTopFace1;
@@ -471,6 +476,7 @@ void getColumnTet(MFaceVecMEltMap &face2el,
// Get third element in layer
std::vector<MElement*> newElts2 = face2el[elTopFace1];
+ if (newElts2.size() < 2) { aboveElt = 0; break; }
el2 = (newElts2[0] == el1) ? newElts2[1] : newElts2[0];
if (el2->getType() != TYPE_TET) break;
MFace elTopFace2;
@@ -507,6 +513,7 @@ void getColumnPrismHex(int elType, MFaceVecMEltMap &face2el,
for (int iLayer = 0; iLayer < p.maxNumLayers; iLayer++) {
std::vector<MElement*> newElts = face2el[elBaseFace];
+ if ((iLayer > 0) && (newElts.size() < 2)) { aboveElt = 0; break; }
el = (newElts[0] == el) ? newElts[1] : newElts[0];
aboveElt = el;
if (el->getType() != elType) break;
@@ -863,10 +870,11 @@ void optimizeCADDist2DP2(GEntity *geomEnt, std::vector<MVertex*> &baseVert)
double curveAndMeasureAboveEl(MetaEl &metaElt, MElement *lastElt,
MElement *aboveElt, double deformFact)
{
- double minJacDet, maxJacDet;
metaElt.setCurvedTop(deformFact);
std::set<MVertex*> movedVertDum;
curveElement(metaElt, movedVertDum, lastElt);
+ if (aboveElt == 0) return 1.;
+ double minJacDet, maxJacDet;
jacobianBasedQuality::minMaxJacobianDeterminant(aboveElt, minJacDet, maxJacDet);
return minJacDet/maxJacDet;
}
@@ -965,6 +973,7 @@ void curveMeshFromBndElt(MEdgeVecMEltMap &ed2el, MFaceVecMEltMap &face2el,
DbgOutputCol dbgOutCol;
dbgOutCol.addBlob(blob);
dbgOutCol.write("col_KO", bndElt->getNum());
+ if (aboveElt == 0) std::cout << "DBGTT: aboveElt = 0 for bnd. elt. " << bndElt->getNum() << std::endl;
curveColumn(p, bndEnt, metaElType, baseVert, topPrimVert, aboveElt, blob,
movedVert, dbgOut);
dbgOutCol.write("col_OK", bndElt->getNum());
--
GitLab