diff --git a/Mesh/HighOrder.cpp b/Mesh/HighOrder.cpp
index a138c4bfd92ae2ebd5396a71cf3fa0f6b30148e3..614bc87e3e4a3f3da91e1f07aeaf56725d3d3222 100644
--- a/Mesh/HighOrder.cpp
+++ b/Mesh/HighOrder.cpp
@@ -1656,7 +1656,6 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete, bool onlyVisi
// Determine mesh dimension and curve BL elements
FastCurvingParameters p;
- p.useBLData = true;
p.dim = 0;
for (GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it)
if ((*it)->getNumMeshElements() > 0) { p.dim = 3; break; }
diff --git a/contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp b/contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp
index 1a7078580baf396787c513d7445d54a57f9dc288..9366e9a960ba75c672b143e03b356880193473af 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp
@@ -209,31 +209,25 @@ void getOppositeEdgeTri(MElement *el, const MEdge &elBaseEd, MEdge &elTopEd,
void getColumnQuad(MEdgeVecMEltMap &ed2el, const FastCurvingParameters &p,
- MEdge &elBaseEd, const MEdge &topEd,
- std::vector<MElement*> &blob, MElement* &aboveElt)
+ MEdge &elBaseEd, std::vector<MElement*> &blob,
+ MElement* &aboveElt)
{
const double maxDP = std::cos(p.maxAngle);
- // Enable geometric stop criteria if top edge not known (i.e. if no BL data)
- bool stopGeom = (topEd.getVertex(0) == 0) || (topEd.getVertex(1) == 0);
-
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];
aboveElt = el;
- if ((!stopGeom) && (elBaseEd == topEd)) break;
if (el->getType() != TYPE_QUA) break;
MEdge elTopEd;
double edLenMin, edLenMax;
getOppositeEdgeQuad(el, elBaseEd, elTopEd, edLenMin, edLenMax);
- if (stopGeom) {
- if (edLenMin > edLenMax*p.maxRho) break;
- const double dp = dot(elBaseEd.normal(), elTopEd.normal());
- if (std::abs(dp) < maxDP) break;
- }
+ 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;
@@ -243,15 +237,12 @@ void getColumnQuad(MEdgeVecMEltMap &ed2el, const FastCurvingParameters &p,
void getColumnTri(MEdgeVecMEltMap &ed2el, const FastCurvingParameters &p,
- MEdge &elBaseEd, const MEdge &topEd,
- std::vector<MElement*> &blob, MElement* &aboveElt)
+ MEdge &elBaseEd, std::vector<MElement*> &blob,
+ MElement* &aboveElt)
{
const double maxDP = std::cos(p.maxAngle);
const double maxDPIn = std::cos(p.maxAngleInner);
- // Enable geometric stop criteria if top edge not known (i.e. if no BL data)
- bool stopGeom = (topEd.getVertex(0) == 0) || (topEd.getVertex(1) == 0);
-
MElement *el0 = 0, *el1 = 0;
for (int iLayer = 0; iLayer < p.maxNumLayers; iLayer++) {
@@ -259,17 +250,13 @@ void getColumnTri(MEdgeVecMEltMap &ed2el, const FastCurvingParameters &p,
std::vector<MElement*> newElts0 = ed2el[elBaseEd];
el0 = (newElts0[0] == el1) ? newElts0[1] : newElts0[0];
aboveElt = el0;
- if ((!stopGeom) && (elBaseEd == topEd)) break;
if (el0->getType() != TYPE_TRI) break;
MEdge elTopEd0;
double edLenMin0, edLenMax0;
getOppositeEdgeTri(el0, elBaseEd, elTopEd0, edLenMin0, edLenMax0);
- SVector3 normBase, normTop0;
- if (stopGeom) {
- normBase = elBaseEd.normal();
- normTop0 = elTopEd0.normal();
- if (std::abs(dot(normBase, normTop0)) < maxDPIn) break;
- }
+ const SVector3 normBase = elBaseEd.normal();
+ const SVector3 normTop0 = elTopEd0.normal();
+ if (std::abs(dot(normBase, normTop0)) < maxDPIn) break;
// Get second element in layer
std::vector<MElement*> newElts1 = ed2el[elTopEd0];
@@ -278,19 +265,14 @@ void getColumnTri(MEdgeVecMEltMap &ed2el, const FastCurvingParameters &p,
MEdge elTopEd1;
double edLenMin1, edLenMax1;
getOppositeEdgeTri(el1, elTopEd0, elTopEd1, edLenMin1, edLenMax1);
- SVector3 normTop1;
- if (stopGeom) {
- normTop1 = elTopEd1.normal();
- if (std::abs(dot(normTop0, normTop1)) < maxDPIn) break;
- }
+ const SVector3 normTop1 = elTopEd1.normal();
+ if (std::abs(dot(normTop0, normTop1)) < maxDPIn) break;
// Check stop criteria
- if (stopGeom) {
- const double edLenMin = std::min(edLenMin0, edLenMin1);
- const double edLenMax = std::max(edLenMax0, edLenMax1);
- if (edLenMin > edLenMax*p.maxRho) break;
- if (std::abs(dot(normBase, normTop1)) < maxDP) break;
- }
+ const double edLenMin = std::min(edLenMin0, edLenMin1);
+ const double edLenMax = std::max(edLenMax0, edLenMax1);
+ if (edLenMin > edLenMax*p.maxRho) break;
+ if (std::abs(dot(normBase, normTop1)) < maxDP) break;
// Add elements to blob and pass top edge to next layer
blob.push_back(el0);
@@ -301,9 +283,8 @@ void getColumnTri(MEdgeVecMEltMap &ed2el, const FastCurvingParameters &p,
-bool getColumn2D(MEdgeVecMEltMap &ed2el, BoundaryLayerColumns *blCols,
- const FastCurvingParameters &p, const MEdge &baseEd,
- std::vector<MVertex*> &baseVert,
+bool getColumn2D(MEdgeVecMEltMap &ed2el, const FastCurvingParameters &p,
+ const MEdge &baseEd, std::vector<MVertex*> &baseVert,
std::vector<MVertex*> &topPrimVert,
std::vector<MElement*> &blob, MElement* &aboveElt)
{
@@ -315,23 +296,10 @@ bool getColumn2D(MEdgeVecMEltMap &ed2el, BoundaryLayerColumns *blCols,
el->getEdgeVertices(iFirstElEd, baseVert);
MEdge elBaseEd(baseVert[0], baseVert[1]);
- // If enabled and boundary layer data exists, retrieve top edge of column
- MEdge topEd;
- if (p.useBLData && (blCols != 0) && (blCols->size() > 0)) {
- const BoundaryLayerData* blData[2] = {0, 0};
- blData[0] = &(blCols->getColumn(baseVert[0], elBaseEd));
- blData[1] = &(blCols->getColumn(baseVert[1], elBaseEd));
- if ((blData[0] != 0) && (blData[1] != 0)) {
- MVertex *topPrimVert0 = blData[0]->_column.back();
- MVertex *topPrimVert1 = blData[1]->_column.back();
- topEd = MEdge(topPrimVert0, topPrimVert1);
- }
- }
-
// Sweep column upwards by choosing largest edges in each element
if (el->getType() == TYPE_TRI)
- getColumnTri(ed2el, p, elBaseEd, topEd, blob, aboveElt);
- else getColumnQuad(ed2el, p, elBaseEd, topEd, blob, aboveElt);
+ getColumnTri(ed2el, p, elBaseEd, blob, aboveElt);
+ else getColumnQuad(ed2el, p, elBaseEd, blob, aboveElt);
topPrimVert.resize(2);
topPrimVert[0] = elBaseEd.getVertex(0);
@@ -957,8 +925,8 @@ void curveColumn(const FastCurvingParameters &p, GEntity *geomEnt,
void curveMeshFromBndElt(MEdgeVecMEltMap &ed2el, MFaceVecMEltMap &face2el,
- GEntity *bndEnt, BoundaryLayerColumns *blCols,
- MElement *bndElt, std::set<MVertex*> movedVert,
+ GEntity *bndEnt, MElement *bndElt,
+ std::set<MVertex*> movedVert,
const FastCurvingParameters &p,
DbgOutputMeta &dbgOut)
{
@@ -973,7 +941,7 @@ void curveMeshFromBndElt(MEdgeVecMEltMap &ed2el, MFaceVecMEltMap &face2el,
MVertex *vb1 = bndElt->getVertex(1);
metaElType = TYPE_QUA;
MEdge baseEd(vb0, vb1);
- foundCol = getColumn2D(ed2el, blCols, p, baseEd, baseVert,
+ foundCol = getColumn2D(ed2el, p, baseEd, baseVert,
topPrimVert, blob, aboveElt);
}
else { // 2D boundary
@@ -1005,8 +973,7 @@ void curveMeshFromBndElt(MEdgeVecMEltMap &ed2el, MFaceVecMEltMap &face2el,
void curveMeshFromBnd(MEdgeVecMEltMap &ed2el, MFaceVecMEltMap &face2el,
- GEntity *bndEnt, BoundaryLayerColumns *blCols,
- const FastCurvingParameters &p)
+ GEntity *bndEnt, const FastCurvingParameters &p)
{
// Build list of bnd. elements to consider
std::list<MElement*> bndEl;
@@ -1031,8 +998,7 @@ void curveMeshFromBnd(MEdgeVecMEltMap &ed2el, MFaceVecMEltMap &face2el,
std::set<MVertex*> movedVert;
for(std::list<MElement*>::iterator itBE = bndEl.begin();
itBE != bndEl.end(); itBE++) // Loop over bnd. elements
- curveMeshFromBndElt(ed2el, face2el, bndEnt, blCols,
- *itBE, movedVert, p, dbgOut);
+ curveMeshFromBndElt(ed2el, face2el, bndEnt, *itBE, movedVert, p, dbgOut);
dbgOut.write("meta-elements", bndEnt->tag());
}
@@ -1076,19 +1042,15 @@ void HighOrderMeshFastCurving(GModel *gm, FastCurvingParameters &p)
bndEnts = std::vector<GEntity*>(gFaces.begin(), gFaces.end());
}
- // Retrieve boudary layer data
- BoundaryLayerColumns *blCols = (p.dim == 2) ?
- gEnt->cast2Face()->getColumns() :
- gEnt->cast2Region()->getColumns();
-
// Curve mesh from each boundary entity
for (int iBndEnt = 0; iBndEnt < bndEnts.size(); iBndEnt++) {
GEntity* &bndEnt = bndEnts[iBndEnt];
if (p.onlyVisible && !bndEnt->getVisibility()) continue;
const GEntity::GeomType bndType = bndEnt->geomType();
if ((bndType == GEntity::Line) || (bndType == GEntity::Plane)) continue;
- Msg::Info("Curving elements for boundary entity %d...", bndEnt->tag());
- curveMeshFromBnd(ed2el, face2el, bndEnt, blCols, p);
+ Msg::Info("Curving elements in entity %d for boundary entity %d...",
+ gEnt->tag(), bndEnt->tag());
+ curveMeshFromBnd(ed2el, face2el, bndEnt, p);
}
}
diff --git a/contrib/HighOrderMeshOptimizer/OptHomFastCurving.h b/contrib/HighOrderMeshOptimizer/OptHomFastCurving.h
index 940d10577ab8624fdde63c4e3dd57682fe69fe13..c87a4a4fbf19cc1e18e6093b567d1f2957034a6b 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomFastCurving.h
+++ b/contrib/HighOrderMeshOptimizer/OptHomFastCurving.h
@@ -35,7 +35,6 @@ class GModel;
struct FastCurvingParameters {
int dim ; // Spatial dimension of the mesh to curve
bool onlyVisible ; // Apply curving to visible entities ONLY?
- bool useBLData; // Use boundary layer data structures created by Gmsh?
bool optimizeGeometry; // Optimize boundary edges/faces to fit geometry?
bool curveOuterBL; // Curve also the outer surface of the boundary layer?
int maxNumLayers; // Maximum number of layers of elements to curve in BL
@@ -44,7 +43,7 @@ struct FastCurvingParameters {
double maxAngleInner; // Maximum angle between edges/faces within layers of triangles/tets to curve in BL
FastCurvingParameters () :
- dim(3), onlyVisible(true), useBLData(false), optimizeGeometry(false),
+ dim(3), onlyVisible(true), optimizeGeometry(false),
curveOuterBL(false), maxNumLayers(100), maxRho(0.5),
maxAngle(3.1415927*10./180.), maxAngleInner(3.1415927*30./180.)
{