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

Added use of boundary layer mesh data structures for high-order fast curving... 

Added use of boundary layer mesh data structures for high-order fast curving (not robust). Added in SetOrderN cleaning up of pointers to destroyed elements in BL mesh data structures.
parent 13c41c3b
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Geo/boundaryLayersData.h
+ 4
− 0
View file @ 4c84d2b7
...@@ -72,6 +72,10 @@ public: ...@@ -72,6 +72,10 @@ public:
_elemColumns.clear(); _elemColumns.clear();
_fans.clear(); _fans.clear();
} }
void clearElementData () {
_toFirst.clear();
_elemColumns.clear();
}
iter begin() { return _data.begin(); } iter begin() { return _data.begin(); }
iter end() { return _data.end(); } iter end() { return _data.end(); }
......
Mesh/HighOrder.cpp
+ 15
− 0
View file @ 4c84d2b7
...@@ -22,6 +22,7 @@ ...@@ -22,6 +22,7 @@
#include "fullMatrix.h" #include "fullMatrix.h"
#include "BasisFactory.h" #include "BasisFactory.h"
#include "MVertexRTree.h" #include "MVertexRTree.h"
#include "OptHomFastCurving.h"
#include <sstream> #include <sstream>
...@@ -1634,6 +1635,7 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete, bool onlyVisi ...@@ -1634,6 +1635,7 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete, bool onlyVisi
Msg::ProgressMeter(++counter, nTot, false, msg); Msg::ProgressMeter(++counter, nTot, false, msg);
if (onlyVisible && !(*it)->getVisibility()) continue; if (onlyVisible && !(*it)->getVisibility()) continue;
setHighOrder(*it, newHOVert[*it], edgeVertices, faceVertices, linear, incomplete, nPts); setHighOrder(*it, newHOVert[*it], edgeVertices, faceVertices, linear, incomplete, nPts);
if ((*it)->getColumns() != 0) (*it)->getColumns()->clearElementData(); // Clear obsolete element data from BL data structures
} }
for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it) { for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it) {
...@@ -1641,6 +1643,7 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete, bool onlyVisi ...@@ -1641,6 +1643,7 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete, bool onlyVisi
Msg::ProgressMeter(++counter, nTot, false, msg); Msg::ProgressMeter(++counter, nTot, false, msg);
if (onlyVisible && !(*it)->getVisibility())continue; if (onlyVisible && !(*it)->getVisibility())continue;
setHighOrder(*it, newHOVert[*it], edgeVertices, faceVertices, linear, incomplete, nPts); setHighOrder(*it, newHOVert[*it], edgeVertices, faceVertices, linear, incomplete, nPts);
if ((*it)->getColumns() != 0) (*it)->getColumns()->clearElementData(); // Clear obsolete element data from BL data structures
} }
// Update all high order vertices // Update all high order vertices
...@@ -1651,6 +1654,18 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete, bool onlyVisi ...@@ -1651,6 +1654,18 @@ void SetOrderN(GModel *m, int order, bool linear, bool incomplete, bool onlyVisi
for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it) for(GModel::riter it = m->firstRegion(); it != m->lastRegion(); ++it)
updateHighOrderVertices(*it, newHOVert[*it], onlyVisible); updateHighOrderVertices(*it, newHOVert[*it], onlyVisible);
// 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; }
if (p.dim == 0)
for(GModel::fiter it = m->firstFace(); it != m->lastFace(); ++it)
if ((*it)->getNumMeshElements() > 0) { p.dim = 2; break; }
if (p.dim > 0)
HighOrderMeshFastCurving(GModel::current(), p);
updatePeriodicEdgesAndFaces(m); updatePeriodicEdgesAndFaces(m);
double t2 = Cpu(); double t2 = Cpu();
......
contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp
+ 181
− 126
View file @ 4c84d2b7
...@@ -209,25 +209,31 @@ void getOppositeEdgeTri(MElement *el, const MEdge &elBaseEd, MEdge &elTopEd, ...@@ -209,25 +209,31 @@ 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, std::vector<MElement*> &blob, MEdge &elBaseEd, const MEdge &topEd,
MElement* &aboveElt) std::vector<MElement*> &blob, 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 (edLenMin > edLenMax*p.maxRho) break; if (stopGeom) {
const double dp = dot(elBaseEd.normal(), elTopEd.normal()); if (edLenMin > edLenMax*p.maxRho) break;
if (std::abs(dp) < maxDP) break; const double dp = dot(elBaseEd.normal(), elTopEd.normal());
if (std::abs(dp) < maxDP) break;
}
blob.push_back(el); blob.push_back(el);
elBaseEd = elTopEd; elBaseEd = elTopEd;
...@@ -237,12 +243,15 @@ void getColumnQuad(MEdgeVecMEltMap &ed2el, const FastCurvingParameters &p, ...@@ -237,12 +243,15 @@ void getColumnQuad(MEdgeVecMEltMap &ed2el, const FastCurvingParameters &p,
void getColumnTri(MEdgeVecMEltMap &ed2el, const FastCurvingParameters &p, void getColumnTri(MEdgeVecMEltMap &ed2el, const FastCurvingParameters &p,
MEdge &elBaseEd, std::vector<MElement*> &blob, MEdge &elBaseEd, const MEdge &topEd,
MElement* &aboveElt) std::vector<MElement*> &blob, 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++) {
...@@ -250,13 +259,17 @@ void getColumnTri(MEdgeVecMEltMap &ed2el, const FastCurvingParameters &p, ...@@ -250,13 +259,17 @@ 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);
const SVector3 normBase = elBaseEd.normal(); SVector3 normBase, normTop0;
const SVector3 normTop0 = elTopEd0.normal(); if (stopGeom) {
if (std::abs(dot(normBase, normTop0)) < maxDPIn) break; normBase = elBaseEd.normal();
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];
...@@ -265,14 +278,19 @@ void getColumnTri(MEdgeVecMEltMap &ed2el, const FastCurvingParameters &p, ...@@ -265,14 +278,19 @@ 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);
const SVector3 normTop1 = elTopEd1.normal(); SVector3 normTop1;
if (std::abs(dot(normTop0, normTop1)) < maxDPIn) break; if (stopGeom) {
normTop1 = elTopEd1.normal();
if (std::abs(dot(normTop0, normTop1)) < maxDPIn) break;
}
// Check stop criteria // Check stop criteria
const double edLenMin = std::min(edLenMin0, edLenMin1); if (stopGeom) {
const double edLenMax = std::max(edLenMax0, edLenMax1); const double edLenMin = std::min(edLenMin0, edLenMin1);
if (edLenMin > edLenMax*p.maxRho) break; const double edLenMax = std::max(edLenMax0, edLenMax1);
if (std::abs(dot(normBase, normTop1)) < maxDP) break; 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 // Add elements to blob and pass top edge to next layer
blob.push_back(el0); blob.push_back(el0);
...@@ -283,8 +301,9 @@ void getColumnTri(MEdgeVecMEltMap &ed2el, const FastCurvingParameters &p, ...@@ -283,8 +301,9 @@ void getColumnTri(MEdgeVecMEltMap &ed2el, const FastCurvingParameters &p,
bool getColumn2D(MEdgeVecMEltMap &ed2el, const FastCurvingParameters &p, bool getColumn2D(MEdgeVecMEltMap &ed2el, BoundaryLayerColumns *blCols,
const MEdge &baseEd, std::vector<MVertex*> &baseVert, const FastCurvingParameters &p, const MEdge &baseEd,
std::vector<MVertex*> &baseVert,
std::vector<MVertex*> &topPrimVert, std::vector<MVertex*> &topPrimVert,
std::vector<MElement*> &blob, MElement* &aboveElt) std::vector<MElement*> &blob, MElement* &aboveElt)
{ {
...@@ -296,10 +315,23 @@ bool getColumn2D(MEdgeVecMEltMap &ed2el, const FastCurvingParameters &p, ...@@ -296,10 +315,23 @@ bool getColumn2D(MEdgeVecMEltMap &ed2el, const FastCurvingParameters &p,
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, blob, aboveElt); getColumnTri(ed2el, p, elBaseEd, topEd, blob, aboveElt);
else getColumnQuad(ed2el, p, elBaseEd, blob, aboveElt); else getColumnQuad(ed2el, p, elBaseEd, topEd, blob, aboveElt);
topPrimVert.resize(2); topPrimVert.resize(2);
topPrimVert[0] = elBaseEd.getVertex(0); topPrimVert[0] = elBaseEd.getVertex(0);
...@@ -873,8 +905,8 @@ double curveAndMeasureAboveEl(MetaEl &metaElt, MElement *lastElt, ...@@ -873,8 +905,8 @@ double curveAndMeasureAboveEl(MetaEl &metaElt, MElement *lastElt,
void curveColumn(GEntity *geomEnt, int metaElType, void curveColumn(const FastCurvingParameters &p, GEntity *geomEnt,
std::vector<MVertex*> &baseVert, int metaElType, std::vector<MVertex*> &baseVert,
const std::vector<MVertex*> &topPrimVert, MElement *aboveElt, const std::vector<MVertex*> &topPrimVert, MElement *aboveElt,
std::vector<MElement*> &blob, std::set<MVertex*> &movedVert, std::vector<MElement*> &blob, std::set<MVertex*> &movedVert,
DbgOutputMeta &dbgOut) DbgOutputMeta &dbgOut)
...@@ -884,49 +916,97 @@ void curveColumn(GEntity *geomEnt, int metaElType, ...@@ -884,49 +916,97 @@ void curveColumn(GEntity *geomEnt, int metaElType,
// Order // Order
const int order = blob[0]->getPolynomialOrder(); const int order = blob[0]->getPolynomialOrder();
// If 2D P2, modify base vertices to minimize distance between wall edge and CAD // If 2D P2 and allowed, modify base vertices to minimize distance between wall edge and CAD
if ((metaElType == TYPE_QUA) && if ((p.optimizeGeometry) && (metaElType == TYPE_QUA) && (order == 2))
(blob[0]->getPolynomialOrder() == 2)) {
optimizeCADDist2DP2(geomEnt, baseVert); optimizeCADDist2DP2(geomEnt, baseVert);
}
// Create meta-element // Create meta-element
MetaEl metaElt(metaElType, order, baseVert, topPrimVert); MetaEl metaElt(metaElType, order, baseVert, topPrimVert);
// Curve top face of meta-element while avoiding breaking the element above // If allowed, curve top face of meta-element while avoiding breaking the element above
// MElement* &lastElt = blob.back(); if (p.curveOuterBL) {
// double minJacDet, maxJacDet; MElement* &lastElt = blob.back();
// double deformMin = 0., qualDeformMin = 1.; double minJacDet, maxJacDet;
// double deformMax = 1.; double deformMin = 0., qualDeformMin = 1.;
// double qualDeformMax = curveAndMeasureAboveEl(metaElt, lastElt, aboveElt, double deformMax = 1.;
// deformMax); double qualDeformMax = curveAndMeasureAboveEl(metaElt, lastElt, aboveElt,
// if (qualDeformMax < MINQUAL) { // Max deformation makes element above invalid deformMax);
// for (int iter = 0; iter < MAXITER; iter++) { // Bisection to find max. deformation that makes element above valid if (qualDeformMax < MINQUAL) { // Max deformation makes element above invalid
// const double deformMid = 0.5 * (deformMin + deformMax); for (int iter = 0; iter < MAXITER; iter++) { // Bisection to find max. deformation that makes element above valid
// const double qualDeformMid = curveAndMeasureAboveEl(metaElt, lastElt, const double deformMid = 0.5 * (deformMin + deformMax);
// aboveElt, deformMid); const double qualDeformMid = curveAndMeasureAboveEl(metaElt, lastElt,
// if (std::abs(deformMax-deformMin) < TOL) break; aboveElt, deformMid);
// const bool signDeformMax = (qualDeformMax < MINQUAL); if (std::abs(deformMax-deformMin) < TOL) break;
// const bool signDeformMid = (qualDeformMid < MINQUAL); const bool signDeformMax = (qualDeformMax < MINQUAL);
// if (signDeformMid == signDeformMax) deformMax = deformMid; const bool signDeformMid = (qualDeformMid < MINQUAL);
// else deformMin = deformMid; if (signDeformMid == signDeformMax) deformMax = deformMid;
// } else deformMin = deformMid;
// metaElt.setFlatTop(); }
// } metaElt.setFlatTop();
// for (int iV = 0; iV < lastElt->getNumVertices(); iV++) }
// movedVert.insert(lastElt->getVertex(iV)); for (int iV = 0; iV < lastElt->getNumVertices(); iV++)
movedVert.insert(lastElt->getVertex(iV));
}
// Curve elements
dbgOut.addMetaEl(metaElt); dbgOut.addMetaEl(metaElt);
for (int iEl = 0; iEl < blob.size(); iEl++) for (int iEl = 0; iEl < blob.size(); iEl++)
// for (int iEl = 0; iEl < blob.size()-1; iEl++)
curveElement(metaElt, movedVert, blob[iEl]); curveElement(metaElt, movedVert, blob[iEl]);
} }
void curveMeshFromBndElt(MEdgeVecMEltMap &ed2el, MFaceVecMEltMap &face2el,
GEntity *bndEnt, BoundaryLayerColumns *blCols,
MElement *bndElt, std::set<MVertex*> movedVert,
const FastCurvingParameters &p,
DbgOutputMeta &dbgOut)
{
const int bndType = bndElt->getType();
int metaElType;
bool foundCol;
std::vector<MVertex*> baseVert, topPrimVert;
std::vector<MElement*> blob;
MElement *aboveElt = 0;
if (bndType == TYPE_LIN) { // 1D boundary
MVertex *vb0 = bndElt->getVertex(0);
MVertex *vb1 = bndElt->getVertex(1);
metaElType = TYPE_QUA;
MEdge baseEd(vb0, vb1);
foundCol = getColumn2D(ed2el, blCols, p, baseEd, baseVert,
topPrimVert, blob, aboveElt);
}
else { // 2D boundary
MVertex *vb0 = bndElt->getVertex(0);
MVertex *vb1 = bndElt->getVertex(1);
MVertex *vb2 = bndElt->getVertex(2);
MVertex *vb3;
if (bndType == TYPE_QUA) {
vb3 = bndElt->getVertex(3);
metaElType = TYPE_HEX;
}
else {
vb3 = 0;
metaElType = TYPE_PRI;
}
MFace baseFace(vb0, vb1, vb2, vb3);
foundCol = getColumn3D(face2el, p, baseFace, baseVert, topPrimVert,
blob, aboveElt);
}
if (!foundCol) return; // Skip bnd. el. if top vertices not found
DbgOutputCol dbgOutCol;
dbgOutCol.addBlob(blob);
dbgOutCol.write("col_KO", bndElt->getNum());
curveColumn(p, bndEnt, metaElType, baseVert, topPrimVert, aboveElt, blob,
movedVert, dbgOut);
dbgOutCol.write("col_OK", bndElt->getNum());
}
void curveMeshFromBnd(MEdgeVecMEltMap &ed2el, MFaceVecMEltMap &face2el, void curveMeshFromBnd(MEdgeVecMEltMap &ed2el, MFaceVecMEltMap &face2el,
GEntity *bndEnt, const FastCurvingParameters &p) GEntity *bndEnt, BoundaryLayerColumns *blCols,
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;
...@@ -943,54 +1023,16 @@ void curveMeshFromBnd(MEdgeVecMEltMap &ed2el, MFaceVecMEltMap &face2el, ...@@ -943,54 +1023,16 @@ void curveMeshFromBnd(MEdgeVecMEltMap &ed2el, MFaceVecMEltMap &face2el,
insertIfCurved(gFace->quadrangles[i], bndEl); insertIfCurved(gFace->quadrangles[i], bndEl);
} }
else else
Msg::Error("Cannot treat model entity %i of dim %i", bndEnt->tag(), Msg::Error("Cannot process model entity %i of dim %i", bndEnt->tag(),
bndEnt->dim()); bndEnt->dim());
// Loop over boundary elements to curve them by columns
DbgOutputMeta dbgOut; DbgOutputMeta dbgOut;
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
const int bndType = (*itBE)->getType(); curveMeshFromBndElt(ed2el, face2el, bndEnt, blCols,
int metaElType; *itBE, movedVert, p, dbgOut);
bool foundCol;
std::vector<MVertex*> baseVert, topPrimVert;
std::vector<MElement*> blob;
MElement *aboveElt = 0;
if (bndType == TYPE_LIN) { // 1D boundary
MVertex *vb0 = (*itBE)->getVertex(0);
MVertex *vb1 = (*itBE)->getVertex(1);
metaElType = TYPE_QUA;
MEdge baseEd(vb0, vb1);
foundCol = getColumn2D(ed2el, p, baseEd, baseVert,
topPrimVert, blob, aboveElt);
}
else { // 2D boundary
MVertex *vb0 = (*itBE)->getVertex(0);
MVertex *vb1 = (*itBE)->getVertex(1);
MVertex *vb2 = (*itBE)->getVertex(2);
MVertex *vb3;
if (bndType == TYPE_QUA) {
vb3 = (*itBE)->getVertex(3);
metaElType = TYPE_HEX;
}
else {
vb3 = 0;
metaElType = TYPE_PRI;
}
MFace baseFace(vb0, vb1, vb2, vb3);
foundCol = getColumn3D(face2el, p, baseFace, baseVert, topPrimVert,
blob, aboveElt);
}
if (!foundCol) continue; // Skip bnd. el. if top vertices not found
DbgOutputCol dbgOutCol;
dbgOutCol.addBlob(blob);
dbgOutCol.write("col_KO", (*itBE)->getNum());
curveColumn(bndEnt, metaElType, baseVert, topPrimVert, aboveElt, blob,
movedVert, dbgOut);
dbgOutCol.write("col_OK", (*itBE)->getNum());
}
dbgOut.write("meta-elements", bndEnt->tag()); dbgOut.write("meta-elements", bndEnt->tag());
} }
...@@ -1004,39 +1046,52 @@ void curveMeshFromBnd(MEdgeVecMEltMap &ed2el, MFaceVecMEltMap &face2el, ...@@ -1004,39 +1046,52 @@ void curveMeshFromBnd(MEdgeVecMEltMap &ed2el, MFaceVecMEltMap &face2el,
void HighOrderMeshFastCurving(GModel *gm, FastCurvingParameters &p) void HighOrderMeshFastCurving(GModel *gm, FastCurvingParameters &p)
{ {
double t1 = Cpu(); double t1 = Cpu();
Msg::StatusBar(true, "Curving high order boundary layer mesh...");
// Retrieve geometric entities
std::vector<GEntity*> allGEnt;
gm->getEntities(allGEnt);
// Curve mesh for non-straight boundary entities
for (int iEnt = 0; iEnt < allGEnt.size(); ++iEnt) {
// Retrive entity
GEntity* &gEnt = allGEnt[iEnt];
if (gEnt->dim() != p.dim) continue;
// Compute edge/face -> elt. connectivity
Msg::Info("Computing connectivity for entity %i...", gEnt->tag());
MEdgeVecMEltMap ed2el;
MFaceVecMEltMap face2el;
if (p.dim == 2) calcEdge2Elements(allGEnt[iEnt], ed2el);
else calcFace2Elements(allGEnt[iEnt], face2el);
// Retrieve boundary entities
std::vector<GEntity*> bndEnts;
if (p.dim == 2) {
std::list<GEdge*> gEds = gEnt->edges();
bndEnts = std::vector<GEntity*>(gEds.begin(), gEds.end());
}
else {
std::list<GFace*> gFaces = gEnt->faces();
bndEnts = std::vector<GEntity*>(gFaces.begin(), gFaces.end());
}
Msg::StatusBar(true, "Optimizing high order mesh..."); // Retrieve boudary layer data
std::vector<GEntity*> allEntities; BoundaryLayerColumns *blCols = (p.dim == 2) ?
gm->getEntities(allEntities); gEnt->cast2Face()->getColumns() :
gEnt->cast2Region()->getColumns();
// Compute vert. -> elt. connectivity
Msg::Info("Computing connectivity..."); // Curve mesh from each boundary entity
MEdgeVecMEltMap ed2el; for (int iBndEnt = 0; iBndEnt < bndEnts.size(); iBndEnt++) {
MFaceVecMEltMap face2el; GEntity* &bndEnt = bndEnts[iBndEnt];
for (int iEnt = 0; iEnt < allEntities.size(); ++iEnt) if (p.onlyVisible && !bndEnt->getVisibility()) continue;
if (p.dim == 2) calcEdge2Elements(allEntities[iEnt], ed2el); const GEntity::GeomType bndType = bndEnt->geomType();
else calcFace2Elements(allEntities[iEnt], face2el); if ((bndType == GEntity::Line) || (bndType == GEntity::Plane)) continue;
Msg::Info("Curving elements for boundary entity %d...", bndEnt->tag());
// Build multimap of each geometric entity to its boundaries curveMeshFromBnd(ed2el, face2el, bndEnt, blCols, p);
std::multimap<GEntity*,GEntity*> entities; }
for (int iEnt = 0; iEnt < allEntities.size(); ++iEnt) {
GEntity *dummy = 0;
GEntity* &entity = allEntities[iEnt];
if (entity->dim() == p.dim-1 &&
(!p.onlyVisible || entity->getVisibility()) &&
(entity->geomType() != GEntity::Plane)) // Consider non-plane boundary entities
entities.insert(std::pair<GEntity*,GEntity*>(dummy, entity));
}
// Loop over geometric entities
for (std::multimap<GEntity*,GEntity*>::iterator itBE = entities.begin();
itBE != entities.end(); itBE++) {
GEntity *domEnt = itBE->first, *bndEnt = itBE->second;
Msg::Info("Curving elements for boundary entity %d...", bndEnt->tag());
curveMeshFromBnd(ed2el, face2el, bndEnt, p);
} }
double t2 = Cpu(); double t2 = Cpu();
Msg::StatusBar(true, "Done curving high order boundary layer mesh (%g s)", t2-t1);
Msg::StatusBar(true, "Done curving high order mesh (%g s)", t2-t1);
} }
contrib/HighOrderMeshOptimizer/OptHomFastCurving.h
+ 5
− 1
View file @ 4c84d2b7
...@@ -35,13 +35,17 @@ class GModel; ...@@ -35,13 +35,17 @@ 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 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
double maxRho; // Maximum ratio min/max edge/face size for 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 double maxAngle; // Maximum angle between layers of elements to curve in BL
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), maxNumLayers(100), maxRho(0.5), dim(3), onlyVisible(true), useBLData(false), optimizeGeometry(false),
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.
Please register or to comment