Skip to content
Snippets Groups Projects
Commit 07e53574 authored by Thomas Toulorge's avatar Thomas Toulorge
Browse files

Removed use of boundary layer mesh data structures for fast high-order curving (unreliable)

parent 4c84d2b7
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Mesh/HighOrder.cpp
+ 0
− 1
View file @ 07e53574
...@@ -1656,7 +1656,6 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete, bool onlyVisi ...@@ -1656,7 +1656,6 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete, bool onlyVisi
// Determine mesh dimension and curve BL elements // Determine mesh dimension and curve BL elements
FastCurvingParameters p; FastCurvingParameters p;
p.useBLData = true;
p.dim = 0; p.dim = 0;
for (GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it) for (GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it)
if ((*it)->getNumMeshElements() > 0) { p.dim = 3; break; } if ((*it)->getNumMeshElements() > 0) { p.dim = 3; break; }
......
contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp
+ 28
− 66
View file @ 07e53574
...@@ -209,31 +209,25 @@ void getOppositeEdgeTri(MElement *el, const MEdge &elBaseEd, MEdge &elTopEd, ...@@ -209,31 +209,25 @@ void getOppositeEdgeTri(MElement *el, const MEdge &elBaseEd, MEdge &elTopEd,
void getColumnQuad(MEdgeVecMEltMap &ed2el, const FastCurvingParameters &p, void getColumnQuad(MEdgeVecMEltMap &ed2el, const FastCurvingParameters &p,
MEdge &elBaseEd, const MEdge &topEd, MEdge &elBaseEd, std::vector<MElement*> &blob,
std::vector<MElement*> &blob, MElement* &aboveElt) MElement* &aboveElt)
{ {
const double maxDP = std::cos(p.maxAngle); 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; MElement *el = 0;
for (int iLayer = 0; iLayer < p.maxNumLayers; iLayer++) { for (int iLayer = 0; iLayer < p.maxNumLayers; iLayer++) {
std::vector<MElement*> newElts = ed2el[elBaseEd]; std::vector<MElement*> newElts = ed2el[elBaseEd];
el = (newElts[0] == el) ? newElts[1] : newElts[0]; el = (newElts[0] == el) ? newElts[1] : newElts[0];
aboveElt = el; aboveElt = el;
if ((!stopGeom) && (elBaseEd == topEd)) break;
if (el->getType() != TYPE_QUA) break; if (el->getType() != TYPE_QUA) break;
MEdge elTopEd; MEdge elTopEd;
double edLenMin, edLenMax; double edLenMin, edLenMax;
getOppositeEdgeQuad(el, elBaseEd, elTopEd, edLenMin, edLenMax); getOppositeEdgeQuad(el, elBaseEd, elTopEd, edLenMin, edLenMax);
if (stopGeom) { if (edLenMin > edLenMax*p.maxRho) break;
if (edLenMin > edLenMax*p.maxRho) break; const double dp = dot(elBaseEd.normal(), elTopEd.normal());
const double dp = dot(elBaseEd.normal(), elTopEd.normal()); if (std::abs(dp) < maxDP) break;
if (std::abs(dp) < maxDP) break;
}
blob.push_back(el); blob.push_back(el);
elBaseEd = elTopEd; elBaseEd = elTopEd;
...@@ -243,15 +237,12 @@ void getColumnQuad(MEdgeVecMEltMap &ed2el, const FastCurvingParameters &p, ...@@ -243,15 +237,12 @@ void getColumnQuad(MEdgeVecMEltMap &ed2el, const FastCurvingParameters &p,
void getColumnTri(MEdgeVecMEltMap &ed2el, const FastCurvingParameters &p, void getColumnTri(MEdgeVecMEltMap &ed2el, const FastCurvingParameters &p,
MEdge &elBaseEd, const MEdge &topEd, MEdge &elBaseEd, std::vector<MElement*> &blob,
std::vector<MElement*> &blob, MElement* &aboveElt) MElement* &aboveElt)
{ {
const double maxDP = std::cos(p.maxAngle); const double maxDP = std::cos(p.maxAngle);
const double maxDPIn = std::cos(p.maxAngleInner); 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; MElement *el0 = 0, *el1 = 0;
for (int iLayer = 0; iLayer < p.maxNumLayers; iLayer++) { for (int iLayer = 0; iLayer < p.maxNumLayers; iLayer++) {
...@@ -259,17 +250,13 @@ void getColumnTri(MEdgeVecMEltMap &ed2el, const FastCurvingParameters &p, ...@@ -259,17 +250,13 @@ void getColumnTri(MEdgeVecMEltMap &ed2el, const FastCurvingParameters &p,
std::vector<MElement*> newElts0 = ed2el[elBaseEd]; std::vector<MElement*> newElts0 = ed2el[elBaseEd];
el0 = (newElts0[0] == el1) ? newElts0[1] : newElts0[0]; el0 = (newElts0[0] == el1) ? newElts0[1] : newElts0[0];
aboveElt = el0; aboveElt = el0;
if ((!stopGeom) && (elBaseEd == topEd)) break;
if (el0->getType() != TYPE_TRI) break; if (el0->getType() != TYPE_TRI) break;
MEdge elTopEd0; MEdge elTopEd0;
double edLenMin0, edLenMax0; double edLenMin0, edLenMax0;
getOppositeEdgeTri(el0, elBaseEd, elTopEd0, edLenMin0, edLenMax0); getOppositeEdgeTri(el0, elBaseEd, elTopEd0, edLenMin0, edLenMax0);
SVector3 normBase, normTop0; const SVector3 normBase = elBaseEd.normal();
if (stopGeom) { const SVector3 normTop0 = elTopEd0.normal();
normBase = elBaseEd.normal(); if (std::abs(dot(normBase, normTop0)) < maxDPIn) break;
normTop0 = elTopEd0.normal();
if (std::abs(dot(normBase, normTop0)) < maxDPIn) break;
}
// Get second element in layer // Get second element in layer
std::vector<MElement*> newElts1 = ed2el[elTopEd0]; std::vector<MElement*> newElts1 = ed2el[elTopEd0];
...@@ -278,19 +265,14 @@ void getColumnTri(MEdgeVecMEltMap &ed2el, const FastCurvingParameters &p, ...@@ -278,19 +265,14 @@ void getColumnTri(MEdgeVecMEltMap &ed2el, const FastCurvingParameters &p,
MEdge elTopEd1; MEdge elTopEd1;
double edLenMin1, edLenMax1; double edLenMin1, edLenMax1;
getOppositeEdgeTri(el1, elTopEd0, elTopEd1, edLenMin1, edLenMax1); getOppositeEdgeTri(el1, elTopEd0, elTopEd1, edLenMin1, edLenMax1);
SVector3 normTop1; const SVector3 normTop1 = elTopEd1.normal();
if (stopGeom) { if (std::abs(dot(normTop0, normTop1)) < maxDPIn) break;
normTop1 = elTopEd1.normal();
if (std::abs(dot(normTop0, normTop1)) < maxDPIn) break;
}
// Check stop criteria // Check stop criteria
if (stopGeom) { const double edLenMin = std::min(edLenMin0, edLenMin1);
const double edLenMin = std::min(edLenMin0, edLenMin1); const double edLenMax = std::max(edLenMax0, edLenMax1);
const double edLenMax = std::max(edLenMax0, edLenMax1); if (edLenMin > edLenMax*p.maxRho) break;
if (edLenMin > edLenMax*p.maxRho) break; if (std::abs(dot(normBase, normTop1)) < maxDP) break;
if (std::abs(dot(normBase, normTop1)) < maxDP) break;
}
// Add elements to blob and pass top edge to next layer // Add elements to blob and pass top edge to next layer
blob.push_back(el0); blob.push_back(el0);
...@@ -301,9 +283,8 @@ void getColumnTri(MEdgeVecMEltMap &ed2el, const FastCurvingParameters &p, ...@@ -301,9 +283,8 @@ void getColumnTri(MEdgeVecMEltMap &ed2el, const FastCurvingParameters &p,
bool getColumn2D(MEdgeVecMEltMap &ed2el, BoundaryLayerColumns *blCols, bool getColumn2D(MEdgeVecMEltMap &ed2el, const FastCurvingParameters &p,
const FastCurvingParameters &p, const MEdge &baseEd, const MEdge &baseEd, std::vector<MVertex*> &baseVert,
std::vector<MVertex*> &baseVert,
std::vector<MVertex*> &topPrimVert, std::vector<MVertex*> &topPrimVert,
std::vector<MElement*> &blob, MElement* &aboveElt) std::vector<MElement*> &blob, MElement* &aboveElt)
{ {
...@@ -315,23 +296,10 @@ bool getColumn2D(MEdgeVecMEltMap &ed2el, BoundaryLayerColumns *blCols, ...@@ -315,23 +296,10 @@ bool getColumn2D(MEdgeVecMEltMap &ed2el, BoundaryLayerColumns *blCols,
el->getEdgeVertices(iFirstElEd, baseVert); el->getEdgeVertices(iFirstElEd, baseVert);
MEdge elBaseEd(baseVert[0], baseVert[1]); 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 // Sweep column upwards by choosing largest edges in each element
if (el->getType() == TYPE_TRI) if (el->getType() == TYPE_TRI)
getColumnTri(ed2el, p, elBaseEd, topEd, blob, aboveElt); getColumnTri(ed2el, p, elBaseEd, blob, aboveElt);
else getColumnQuad(ed2el, p, elBaseEd, topEd, blob, aboveElt); else getColumnQuad(ed2el, p, elBaseEd, blob, aboveElt);
topPrimVert.resize(2); topPrimVert.resize(2);
topPrimVert[0] = elBaseEd.getVertex(0); topPrimVert[0] = elBaseEd.getVertex(0);
...@@ -957,8 +925,8 @@ void curveColumn(const FastCurvingParameters &p, GEntity *geomEnt, ...@@ -957,8 +925,8 @@ void curveColumn(const FastCurvingParameters &p, GEntity *geomEnt,
void curveMeshFromBndElt(MEdgeVecMEltMap &ed2el, MFaceVecMEltMap &face2el, void curveMeshFromBndElt(MEdgeVecMEltMap &ed2el, MFaceVecMEltMap &face2el,
GEntity *bndEnt, BoundaryLayerColumns *blCols, GEntity *bndEnt, MElement *bndElt,
MElement *bndElt, std::set<MVertex*> movedVert, std::set<MVertex*> movedVert,
const FastCurvingParameters &p, const FastCurvingParameters &p,
DbgOutputMeta &dbgOut) DbgOutputMeta &dbgOut)
{ {
...@@ -973,7 +941,7 @@ void curveMeshFromBndElt(MEdgeVecMEltMap &ed2el, MFaceVecMEltMap &face2el, ...@@ -973,7 +941,7 @@ void curveMeshFromBndElt(MEdgeVecMEltMap &ed2el, MFaceVecMEltMap &face2el,
MVertex *vb1 = bndElt->getVertex(1); MVertex *vb1 = bndElt->getVertex(1);
metaElType = TYPE_QUA; metaElType = TYPE_QUA;
MEdge baseEd(vb0, vb1); MEdge baseEd(vb0, vb1);
foundCol = getColumn2D(ed2el, blCols, p, baseEd, baseVert, foundCol = getColumn2D(ed2el, p, baseEd, baseVert,
topPrimVert, blob, aboveElt); topPrimVert, blob, aboveElt);
} }
else { // 2D boundary else { // 2D boundary
...@@ -1005,8 +973,7 @@ void curveMeshFromBndElt(MEdgeVecMEltMap &ed2el, MFaceVecMEltMap &face2el, ...@@ -1005,8 +973,7 @@ void curveMeshFromBndElt(MEdgeVecMEltMap &ed2el, MFaceVecMEltMap &face2el,
void curveMeshFromBnd(MEdgeVecMEltMap &ed2el, MFaceVecMEltMap &face2el, void curveMeshFromBnd(MEdgeVecMEltMap &ed2el, MFaceVecMEltMap &face2el,
GEntity *bndEnt, BoundaryLayerColumns *blCols, GEntity *bndEnt, const FastCurvingParameters &p)
const FastCurvingParameters &p)
{ {
// Build list of bnd. elements to consider // Build list of bnd. elements to consider
std::list<MElement*> bndEl; std::list<MElement*> bndEl;
...@@ -1031,8 +998,7 @@ void curveMeshFromBnd(MEdgeVecMEltMap &ed2el, MFaceVecMEltMap &face2el, ...@@ -1031,8 +998,7 @@ void curveMeshFromBnd(MEdgeVecMEltMap &ed2el, MFaceVecMEltMap &face2el,
std::set<MVertex*> movedVert; std::set<MVertex*> movedVert;
for(std::list<MElement*>::iterator itBE = bndEl.begin(); for(std::list<MElement*>::iterator itBE = bndEl.begin();
itBE != bndEl.end(); itBE++) // Loop over bnd. elements itBE != bndEl.end(); itBE++) // Loop over bnd. elements
curveMeshFromBndElt(ed2el, face2el, bndEnt, blCols, curveMeshFromBndElt(ed2el, face2el, bndEnt, *itBE, movedVert, p, dbgOut);
*itBE, movedVert, p, dbgOut);
dbgOut.write("meta-elements", bndEnt->tag()); dbgOut.write("meta-elements", bndEnt->tag());
} }
...@@ -1076,19 +1042,15 @@ void HighOrderMeshFastCurving(GModel *gm, FastCurvingParameters &p) ...@@ -1076,19 +1042,15 @@ void HighOrderMeshFastCurving(GModel *gm, FastCurvingParameters &p)
bndEnts = std::vector<GEntity*>(gFaces.begin(), gFaces.end()); 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 // Curve mesh from each boundary entity
for (int iBndEnt = 0; iBndEnt < bndEnts.size(); iBndEnt++) { for (int iBndEnt = 0; iBndEnt < bndEnts.size(); iBndEnt++) {
GEntity* &bndEnt = bndEnts[iBndEnt]; GEntity* &bndEnt = bndEnts[iBndEnt];
if (p.onlyVisible && !bndEnt->getVisibility()) continue; if (p.onlyVisible && !bndEnt->getVisibility()) continue;
const GEntity::GeomType bndType = bndEnt->geomType(); const GEntity::GeomType bndType = bndEnt->geomType();
if ((bndType == GEntity::Line) || (bndType == GEntity::Plane)) continue; if ((bndType == GEntity::Line) || (bndType == GEntity::Plane)) continue;
Msg::Info("Curving elements for boundary entity %d...", bndEnt->tag()); Msg::Info("Curving elements in entity %d for boundary entity %d...",
curveMeshFromBnd(ed2el, face2el, bndEnt, blCols, p); gEnt->tag(), bndEnt->tag());
curveMeshFromBnd(ed2el, face2el, bndEnt, p);
} }
} }
......
contrib/HighOrderMeshOptimizer/OptHomFastCurving.h
+ 1
− 2
View file @ 07e53574
...@@ -35,7 +35,6 @@ class GModel; ...@@ -35,7 +35,6 @@ class GModel;
struct FastCurvingParameters { struct FastCurvingParameters {
int dim ; // Spatial dimension of the mesh to curve int dim ; // Spatial dimension of the mesh to curve
bool onlyVisible ; // Apply curving to visible entities ONLY? 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 optimizeGeometry; // Optimize boundary edges/faces to fit geometry?
bool curveOuterBL; // Curve also the outer surface of the boundary layer? bool curveOuterBL; // Curve also the outer surface of the boundary layer?
int maxNumLayers; // Maximum number of layers of elements to curve in BL int maxNumLayers; // Maximum number of layers of elements to curve in BL
...@@ -44,7 +43,7 @@ struct FastCurvingParameters { ...@@ -44,7 +43,7 @@ struct FastCurvingParameters {
double maxAngleInner; // Maximum angle between edges/faces within layers of triangles/tets to curve in BL double maxAngleInner; // Maximum angle between edges/faces within layers of triangles/tets to curve in BL
FastCurvingParameters () : FastCurvingParameters () :
dim(3), onlyVisible(true), useBLData(false), optimizeGeometry(false), dim(3), onlyVisible(true), optimizeGeometry(false),
curveOuterBL(false), maxNumLayers(100), maxRho(0.5), curveOuterBL(false), maxNumLayers(100), maxRho(0.5),
maxAngle(3.1415927*10./180.), maxAngleInner(3.1415927*30./180.) maxAngle(3.1415927*10./180.), maxAngleInner(3.1415927*30./180.)
{ {
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment