From ff469316f8889dfe001a0340ec66ae4c11224c66 Mon Sep 17 00:00:00 2001
From: Thomas Toulorge <thomas.toulorge@mines-paristech.fr>
Date: Fri, 9 May 2014 17:59:50 +0000
Subject: [PATCH] Improved blob creation in OptHOM
---
contrib/HighOrderMeshOptimizer/OptHomRun.cpp | 50 +++++++++++++++-----
1 file changed, 39 insertions(+), 11 deletions(-)
diff --git a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
index 3630ae260b..f99522a5ea 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
@@ -158,13 +158,50 @@ static std::set<MVertex *> getAllBndVertices
return bnd;
}
+// Approximate test of intersection element with circle/sphere by sampling
+static bool testElInDist(const SPoint3 p, double limDist, MElement *el)
+{
+ const double sampleLen = 0.5*limDist; // Distance between sample points
+
+ if (el->getDim() == 2) { // 2D?
+ for (int iEd = 0; iEd < el->getNumEdges(); iEd++) { // Loop over edges of element
+ MEdge ed = el->getEdge(iEd);
+ const int nPts = int(ed.length()/sampleLen)+2; // Nb of sample points based on edge length
+ for (int iPt = 0; iPt < nPts; iPt++) { // Loop over sample points
+ const SPoint3 pt = ed.interpolate(iPt/float(nPts-1));
+ if (p.distance(pt) < limDist) return true;
+ }
+ }
+ }
+ else { // 3D
+ for (int iFace = 0; iFace < el->getNumFaces(); iFace++) { // Loop over faces of element
+ MFace face = el->getFace(iFace);
+ double lMax = 0.; // Max. edge length in face
+ const int nVert = face.getNumVertices();
+ for (int iEd = 0; iEd < nVert; iEd++)
+ lMax = std::max(lMax, face.getEdge(iEd).length());
+ const int nPts = int(lMax/sampleLen)+2; // Nb of sample points based on max. edge length in face
+ for (int iPt0 = 0; iPt0 < nPts; iPt0++) {
+ const double u = iPt0/float(nPts-1);
+ for (int iPt1 = 0; iPt1 < nPts; iPt1++) { // Loop over sample points
+ const double vMax = (nVert == 3) ? 1.-u : 1.;
+ const SPoint3 pt = face.interpolate(u, vMax*iPt1/float(nPts-1));
+ if (p.distance(pt) < limDist) return true;
+ }
+ }
+ }
+ }
+
+ return false;
+}
+
static std::set<MElement*> getSurroundingBlob
(MElement *el, int depth,
const std::map<MVertex*, std::vector<MElement*> > &vertex2elements,
const double distFactor, int forceDepth, bool optPrimSurfMesh)
{
- const SPoint3 p = el->barycenter();
+ const SPoint3 p = el->barycenter(true);
const double dist = el->maxDistToStraight();
const double limDist = ((optPrimSurfMesh && (dist < 1.e-10)) ?
el->getOuterRadius() : dist) * distFactor;
@@ -183,21 +220,12 @@ static std::set<MElement*> getSurroundingBlob
((*it)->getVertex(i))->second;
for (std::vector<MElement*>::const_iterator itN = neighbours.begin();
itN != neighbours.end(); ++itN){
- if ((d < forceDepth) || (p.distance((*itN)->barycenter()) < limDist)){
+ if ((d < forceDepth) || testElInDist(p, limDist, *itN)){
// Assume that if an el is too far, its neighbours are too far as well
if (blob.insert(*itN).second) currentLayer.push_back(*itN);
}
}
}
-// for (int i = 0; i < (*it)->getNumPrimaryVertices(); ++i) {
-// if ((d < forceDepth) || (p.distance((*it)->getVertex(i)->point()) < limDist)) {
-// const std::vector<MElement*> &neighbours = vertex2elements.find((*it)->getVertex(i))->second;
-// for (std::vector<MElement*>::const_iterator itN = neighbours.begin(); itN != neighbours.end(); ++itN) {
-// // Assume that if an el is too far, its neighbours are too far as well
-// if (blob.insert(*itN).second) currentLayer.push_back(*itN);
-// }
-// }
-// }
}
lastLayer = currentLayer;
}
--
GitLab