From 7cdba047dc627c41115b10abac09e89b3d354b2a Mon Sep 17 00:00:00 2001
From: Thomas Toulorge <thomas.toulorge@mines-paristech.fr>
Date: Thu, 9 Mar 2017 15:08:48 +0000
Subject: [PATCH] Improved boundary layer detection for fast high-order curving
---
Fltk/highOrderToolsWindow.cpp | 3 +-
.../OptHomFastCurving.cpp | 107 +++++++++++++-----
.../OptHomFastCurving.h | 13 ++-
3 files changed, 86 insertions(+), 37 deletions(-)
diff --git a/Fltk/highOrderToolsWindow.cpp b/Fltk/highOrderToolsWindow.cpp
index 9e035fb266..3a49732803 100644
--- a/Fltk/highOrderToolsWindow.cpp
+++ b/Fltk/highOrderToolsWindow.cpp
@@ -124,7 +124,7 @@ static void chooseopti_cb(Fl_Widget *w, void *data)
o->choice[0]->deactivate();
o->choice[3]->deactivate();
o->value[1]->deactivate();
- o->value[2]->activate();
+ o->value[2]->deactivate();
o->value[3]->deactivate();
o->value[4]->deactivate();
o->value[5]->deactivate();
@@ -197,7 +197,6 @@ static void highordertools_runopti_cb(Fl_Widget *w, void *data)
FastCurvingParameters p;
p.onlyVisible = onlyVisible;
p.dim = dim;
- p.maxNumLayers = (int) o->value[2]->value();
HighOrderMeshFastCurving(GModel::current(), p);
break;
}
diff --git a/contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp b/contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp
index 4d057c9d00..afbe2d0f4b 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp
@@ -27,27 +27,28 @@
//
// Contributors: Thomas Toulorge, Jonathan Lambrechts
-#include <stdio.h>
+#include <cstdio>
#include <sstream>
#include <fstream>
-#include <iterator>
-#include <string.h>
-#include "GmshConfig.h"
+#include <string>
+#include <cmath>
#include "OptHomFastCurving.h"
+#include "GmshConfig.h"
#include "GModel.h"
#include "Gmsh.h"
+#include "MLine.h"
#include "MTriangle.h"
#include "MQuadrangle.h"
#include "MTetrahedron.h"
#include "MHexahedron.h"
#include "MPrism.h"
-#include "MLine.h"
+#include "MEdge.h"
+#include "MFace.h"
#include "OS.h"
-#include <stack>
#include "SVector3.h"
#include "BasisFactory.h"
-#include "Field.h"
#include "MetaEl.h"
+#include "qualityMeasuresJacobian.h"
@@ -205,15 +206,72 @@ void getOppositeEdge(MElement *el, const MEdge &elBaseEd, MEdge &elTopEd,
}
+void getColumnQuad(MEdgeVecMEltMap &ed2el, const FastCurvingParameters &p,
+ MEdge &elBaseEd, std::vector<MElement*> &blob)
+{
+ const double maxDP = std::cos(p.maxAngle);
+
+ MElement *el = 0;
+
+ for (int iLayer = 0; iLayer < p.maxNumLayers; iLayer++) {
+ std::vector<MElement*> newElts = ed2el[elBaseEd];
+ el = (newElts[0] == el) ? newElts[1] : newElts[0];
+ if (el->getType() != TYPE_QUA) break;
+ MEdge elTopEd;
+ double edLenMin, edLenMax;
+ getOppositeEdge(el, elBaseEd, elTopEd, edLenMin, edLenMax);
+
+ if (edLenMin > edLenMax*p.maxRho) break;
+ const double dp = dot(elBaseEd.normal(), elTopEd.normal());
+ if (std::abs(dp) < maxDP) break;
+
+ blob.push_back(el);
+ elBaseEd = elTopEd;
+ }
+}
-bool getColumn2D(MEdgeVecMEltMap &ed2el,
- int maxNumLayers, const MEdge &baseEd,
- std::vector<MVertex*> &baseVert,
+
+void getColumnTri(MEdgeVecMEltMap &ed2el, const FastCurvingParameters &p,
+ MEdge &elBaseEd, std::vector<MElement*> &blob)
+{
+ const double maxDP = std::cos(p.maxAngle);
+
+ MElement *el0 = 0, *el1 = 0;
+
+ for (int iLayer = 0; iLayer < p.maxNumLayers; iLayer++) {
+ std::vector<MElement*> newElts0 = ed2el[elBaseEd];
+ el0 = (newElts0[0] == el1) ? newElts0[1] : newElts0[0];
+ if (el0->getType() != TYPE_TRI) break;
+ MEdge elMidEd;
+ double edLenMin0, edLenMax0;
+ getOppositeEdge(el0, elBaseEd, elMidEd, edLenMin0, edLenMax0);
+
+ std::vector<MElement*> newElts1 = ed2el[elMidEd];
+ el1 = (newElts1[0] == el0) ? newElts1[1] : newElts1[0];
+ if (el1->getType() != TYPE_TRI) break;
+ MEdge elTopEd;
+ double edLenMin1, edLenMax1;
+ getOppositeEdge(el1, elMidEd, elTopEd, edLenMin1, edLenMax1);
+
+ const double edLenMin = std::min(edLenMin0, edLenMin1);
+ const double edLenMax = std::max(edLenMax0, edLenMax1);
+ if (edLenMin > edLenMax*p.maxRho) break;
+ const double dp = dot(elBaseEd.normal(), elTopEd.normal());
+ if (std::abs(dp) < maxDP) break;
+
+ blob.push_back(el0);
+ blob.push_back(el1);
+ elBaseEd = elTopEd;
+ }
+}
+
+
+
+bool getColumn2D(MEdgeVecMEltMap &ed2el, const FastCurvingParameters &p,
+ const MEdge &baseEd, std::vector<MVertex*> &baseVert,
std::vector<MVertex*> &topPrimVert,
std::vector<MElement*> &blob)
{
- static const double tol = 2.;
-
// Get first element and base vertices
std::vector<MElement*> firstElts = ed2el[baseEd];
if (firstElts.size() != 1) {
@@ -225,18 +283,10 @@ bool getColumn2D(MEdgeVecMEltMap &ed2el,
const int iFirstElEd = getElementEdge(baseEd, el);
el->getEdgeVertices(iFirstElEd, baseVert);
MEdge elBaseEd(baseVert[0], baseVert[1]);
- MEdge elTopEd;
// Sweep column upwards by choosing largest edges in each element
- for (int iLayer = 0; iLayer < maxNumLayers; iLayer++) {
- double edLenMin, edLenMax;
- getOppositeEdge(el, elBaseEd, elTopEd, edLenMin, edLenMax);
- if (edLenMax < edLenMin*tol) break;
- blob.push_back(el);
- std::vector<MElement*> newElts = ed2el[elTopEd];
- el = (newElts[0] == el) ? newElts[1] : newElts[0];
- elBaseEd = elTopEd;
- }
+ if (el->getType() == TYPE_TRI) getColumnTri(ed2el, p, elBaseEd, blob);
+ else getColumnQuad(ed2el, p, elBaseEd, blob);
// TODO: improve by taking all vertices?
topPrimVert.resize(2);
@@ -248,9 +298,8 @@ bool getColumn2D(MEdgeVecMEltMap &ed2el,
// TODO: Implement 3D
-bool getColumn3D(MFaceVecMEltMap &face2el,
- int maxNumLayers, const MFace &baseFace,
- std::vector<MVertex*> &baseVert,
+bool getColumn3D(MFaceVecMEltMap &face2el, const FastCurvingParameters &p,
+ const MFace &baseFace, std::vector<MVertex*> &baseVert,
std::vector<MVertex*> &topPrimVert,
std::vector<MElement*> &blob)
{
@@ -303,7 +352,7 @@ private:
void curveMeshFromBnd(MEdgeVecMEltMap &ed2el, MFaceVecMEltMap &face2el,
- GEntity *bndEnt, FastCurvingParameters &p)
+ GEntity *bndEnt, const FastCurvingParameters &p)
{
// Build list of bnd. elements to consider
std::list<MElement*> bndEl;
@@ -338,8 +387,7 @@ void curveMeshFromBnd(MEdgeVecMEltMap &ed2el, MFaceVecMEltMap &face2el,
MVertex *vb1 = (*itBE)->getVertex(1);
metaElType = TYPE_QUA;
MEdge baseEd(vb0, vb1);
- foundCol = getColumn2D(ed2el, p.maxNumLayers, baseEd, baseVert,
- topPrimVert, blob);
+ foundCol = getColumn2D(ed2el, p, baseEd, baseVert, topPrimVert, blob);
}
else { // 2D boundary
MVertex *vb0 = (*itBE)->getVertex(0);
@@ -355,8 +403,7 @@ void curveMeshFromBnd(MEdgeVecMEltMap &ed2el, MFaceVecMEltMap &face2el,
metaElType = TYPE_PRI;
}
MFace baseFace(vb0, vb1, vb2, vb3);
- foundCol = getColumn3D(face2el, p.maxNumLayers, baseFace, baseVert,
- topPrimVert, blob);
+ foundCol = getColumn3D(face2el, p, baseFace, baseVert, topPrimVert, blob);
}
if (!foundCol) continue; // Skip bnd. el. if top vertices not found
int order = blob[0]->getPolynomialOrder();
diff --git a/contrib/HighOrderMeshOptimizer/OptHomFastCurving.h b/contrib/HighOrderMeshOptimizer/OptHomFastCurving.h
index 7fd62f20be..f34571fc6f 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomFastCurving.h
+++ b/contrib/HighOrderMeshOptimizer/OptHomFastCurving.h
@@ -33,12 +33,15 @@
class GModel;
struct FastCurvingParameters {
- int dim ; // Which dimension to curve
- bool onlyVisible ; // Apply curving to visible entities ONLY
- int maxNumLayers; // Maximum number of element layers to consider when trying to detect BL
+ int dim ; // Spatial dimension of the mesh to curve
+ bool onlyVisible ; // Apply curving to visible entities ONLY?
+ int maxNumLayers; // Maximum number of layers of elements to curve in BL
+ double maxRho; // Maximum ratio min/max edge/face size for elements to curve in BL
+ double maxAngle; // Maximum angle between layers of elements to curve in BL
- FastCurvingParameters ()
- : dim(3) , onlyVisible(true), maxNumLayers(6)
+ FastCurvingParameters () :
+ dim(3) , onlyVisible(true), maxNumLayers(100), maxRho(0.5),
+ maxAngle(3.1415926535897932*10./180.)
{
}
};
--
GitLab