From 2b8b05feb2e75ccdd2f8d4e9cd1b4c3b0c305a78 Mon Sep 17 00:00:00 2001
From: Thomas Toulorge <thomas.toulorge@mines-paristech.fr>
Date: Fri, 9 Aug 2013 17:34:36 +0000
Subject: [PATCH] Cleaned up and fixed bugs in high-order fast curving tool
---
.../OptHomFastCurving.cpp | 306 +-----------------
contrib/HighOrderMeshOptimizer/SuperEl.cpp | 32 +-
2 files changed, 17 insertions(+), 321 deletions(-)
diff --git a/contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp b/contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp
index dbcb9d030f..a758a22276 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp
@@ -117,36 +117,6 @@ void exportMeshToDassault(GModel *gm, const std::string &fn, int dim)
-//double distMaxStraight(MElement *el)
-//{
-// const polynomialBasis *lagrange = (polynomialBasis*)el->getFunctionSpace();
-// const polynomialBasis *lagrange1 = (polynomialBasis*)el->getFunctionSpace(1);
-// int nV = lagrange->points.size1();
-// int nV1 = lagrange1->points.size1();
-// SPoint3 sxyz[256];
-// for (int i = 0; i < nV1; ++i) {
-// sxyz[i] = el->getVertex(i)->point();
-// }
-// for (int i = nV1; i < nV; ++i) {
-// double f[256];
-// lagrange1->f(lagrange->points(i, 0), lagrange->points(i, 1),
-// lagrange->points(i, 2), f);
-// for (int j = 0; j < nV1; ++j)
-// sxyz[i] += sxyz[j] * f[j];
-// }
-//
-// double maxdx = 0.0;
-// for (int iV = nV1; iV < nV; iV++) {
-// SVector3 d = el->getVertex(iV)->point()-sxyz[iV];
-// double dx = d.norm();
-// if (dx > maxdx) maxdx = dx;
-// }
-//
-// return maxdx;
-//}
-
-
-
void calcVertex2Elements(int dim, GEntity *entity,
std::map<MVertex*, std::vector<MElement *> > &vertex2elements)
{
@@ -160,100 +130,6 @@ void calcVertex2Elements(int dim, GEntity *entity,
-//void getNormals2D(const std::vector<MVertex*> &edgeVert, SVector3 &n0, SVector3 &n1) {
-//
-//// if (v0->onWhat() != 0) {
-// if (false) {
-// GEdge *ent = edgeVert.front()->onWhat()->cast2Edge();
-// std::map<MVertex*, SVector3, std::less<MVertex*> > &allNormals = ent->getNormals();
-// n0 = allNormals[edgeVert[0]];
-// n1 = allNormals[edgeVert[1]];
-// }
-// else {
-// const int &Nsf = edgeVert.size();
-// const nodalBasis* fs = BasisFactory::getNodalBasis(nodalBasis::getTag(TYPE_LIN,Nsf-1,false));
-// double dsf0[Nsf][3], dsf1[Nsf][3];
-// fs->df(fs->points(0,0),0.,0.,dsf0);
-// fs->df(fs->points(1,0),0.,0.,dsf1);
-// double t0x = 0., t0y = 0., t1x = 0., t1y = 0.;
-// for (int i=0; i<Nsf; ++i) {
-// t0x += dsf0[i][0]*edgeVert[i]->x();
-// t0y += dsf0[i][0]*edgeVert[i]->y();
-// t1x += dsf1[i][0]*edgeVert[i]->x();
-// t1y += dsf1[i][0]*edgeVert[i]->y();
-// }
-// n0[0] = -t0y; n0[1] = t0x; n0[2] = 0.;
-// n0.normalize();
-// n1[0] = -t1y; n1[1] = t1x; n1[2] = 0.;
-// n1.normalize();
-// }
-//
-//}
-
-
-
-//void getBaseVertsAndNormals(std::set<MElement*> badElts, std::map<MElement*,std::vector<MVertex*> > &faces,
-// std::map<MVertex*,SVector3> &normals) {
-//
-// typedef std::map<MVertex*,std::vector<SVector3> > normalMap;
-// normalMap allNormals;
-//
-// for (std::set<MElement*>::iterator itBE = badElts.begin(); itBE != badElts.end(); itBE++) {
-// std::vector<MVertex*> baseVert;
-// for (int i=0; i<(*itBE)->getNumEdges(); i++) {
-// (*itBE)->getEdgeVertices(i,baseVert);
-// if (isEdgeCurved(baseVert)) {
-// const SPoint3 p0 = baseVert[0]->point(), p1 = baseVert[1]->point(), pC = (*itBE)->barycenter(true);
-// SVector3 n(p0.y()-p1.y(), p1.x()-p0.x(), 0.), p0C(p0,pC);
-// if (dot(n,p0C) < 0.) n *= -1.;
-// n.normalize();
-// allNormals[baseVert[0]].push_back(n);
-// allNormals[baseVert[1]].push_back(n);
-// faces[*itBE] = baseVert;
-// break;
-// }
-// }
-// }
-//
-// for(normalMap::iterator itAN = allNormals.begin(); itAN != allNormals.end(); itAN++) {
-// std::vector<SVector3> &normSet = itAN->second;
-//// std::cout << "DBGTT: Considering node " << itAN->first->getNum() << "\n";
-// SVector3 norm(0.);
-// for (std::vector<SVector3>::iterator itN = normSet.begin(); itN != normSet.end(); itN++) {
-//// std::cout << "DBGTT: -> accumulating normal (" << itN->x() << ", "<< itN->y() << ", "<< itN->z() << ")\n";
-// norm += *itN;
-// }
-// norm *= 1./double(normSet.size());
-// normals[itAN->first] = norm;
-// }
-//
-//}
-
-
-
-//bool isEdgeCurved(const std::vector<MVertex*> &edgeVert) {
-//
-// static const double eps = 1.e-8;
-//
-// const double n = edgeVert.size();
-// const MVertex *vS = edgeVert[0];
-// const double vSEx = edgeVert[1]->x()-vS->x(), vSEy = edgeVert[1]->y()-vS->y();
-// const double normSE = vSEx*vSEx+vSEy*vSEy;
-//
-// double distSq = 0.;
-// for (int i=2; i<n; ++i) {
-// const double vSix = edgeVert[i]->x()-vS->x(), vSiy = edgeVert[i]->y()-vS->y();
-// const double fact = (vSix*vSEx+vSiy*vSEy)/normSE;
-// const double dvx = vSix-fact*vSEx, dvy = vSiy-fact*vSEy;
-// distSq += dvx*dvx+dvy*dvy;
-// }
-//
-// return (distSq > eps*normSE);
-//
-//}
-
-
-
// Among edges connected to a given vertex, return the direction of the one that is closest to the given normal
// Return the given normal if no edge is sufficiently close
SVector3 getNormalEdge(MVertex *vert, const SVector3 &n,
@@ -295,96 +171,6 @@ SVector3 getNormalEdge(MVertex *vert, const SVector3 &n,
-//// Detect whether edge/face is curved, and give (non-oriented) normal
-//bool isCurvedAndNormal(int type, int order, const std::vector<MVertex*> &faceVert, SVector3 &normal) {
-//
-// static const double eps = 1.e-8;
-//
-// const nodalBasis *lagrange = BasisFactory::getNodalBasis(ElementType::getTag(type,order,false));
-// const nodalBasis *lagrange1 = BasisFactory::getNodalBasis(ElementType::getTag(type,1,false));
-// int nV = lagrange->points.size1();
-// int nV1 = lagrange1->points.size1();
-// SPoint3 sxyz[256];
-// for (int i = 0; i < nV1; ++i) sxyz[i] = faceVert[i]->point();
-// for (int i = nV1; i < nV; ++i) {
-// double f[256];
-// lagrange1->f(lagrange->points(i, 0), lagrange->points(i, 1), lagrange->points(i, 2), f);
-// for (int j = 0; j < nV1; ++j)
-// sxyz[i] += sxyz[j] * f[j];
-// }
-//
-// double sumDistSq = 0.;
-// for (int iV = nV1; iV < nV; iV++) sumDistSq += SVector3(sxyz[iV],faceVert[iV]->point()).normSq();
-//
-// double scale;
-// if (type == TYPE_LIN) {
-// const SPoint3 &p0 = sxyz[0], &p1 = sxyz[1];
-// normal = SVector3(p0.y()-p1.y(),p1.x()-p0.x(),0.);
-// scale = normal.normSq();
-// }
-// else {
-// const SPoint3 &p0 = sxyz[0], &p1 = sxyz[1], &p2 = sxyz[2];
-// SVector3 p01(p0,p1), p02(p0,p2);
-// normal = crossprod(p01,p02);
-// scale = normal.norm();
-// }
-// normal.normalize();
-//
-//// std::cout << "DBGTT: v0 is " << faceVert[0]->getNum() << ", v1 is " << faceVert[1]->getNum() << ", sumDistSq = " << sumDistSq << ", scale = " << scale << ", test = " << (sumDistSq > eps*double(nV)*scale) << "\n";
-// return (sumDistSq > eps*double(nV)*scale);
-//
-//}
-//
-//
-//
-//void getBaseVertsAndNormals(std::set<MElement*> badElts,
-// const std::map<MVertex*, std::vector<MElement*> > &vertex2elements,
-// std::map<MElement*,int> &baseType,
-// std::map<MElement*,std::vector<MVertex*> > &baseVert,
-// std::map<MElement*,std::vector<SVector3> > &normals) {
-//
-// for (std::set<MElement*>::iterator itBE = badElts.begin(); itBE != badElts.end(); itBE++) {
-// const bool is2D = ((*itBE)->getDim() == 2);
-// const int order = (*itBE)->getPolynomialOrder();
-// const int numFaces = is2D ? (*itBE)->getNumEdges() : (*itBE)->getNumFaces();
-// for (int iFace=0; iFace<numFaces; iFace++) { // Loop on edges/faces
-// std::vector<MVertex*> bv; // Vertices of edge/face
-// int type, numPrimVert;
-// if (is2D) {
-// (*itBE)->getEdgeVertices(iFace,bv);
-// type = TYPE_LIN;
-// numPrimVert = 2;
-// }
-// else {
-// (*itBE)->getFaceVertices(iFace,bv);
-// MFace face = (*itBE)->getFace(iFace);
-// numPrimVert = face.getNumVertices();
-// type = (numPrimVert == 3) ? TYPE_TRI : TYPE_QUA;
-// }
-// SVector3 n; // Normal to straight edge/face
-// if (isCurvedAndNormal(type,order,bv,n)) { // Check whether edge/face is curved and compute straight normal
-// std::cout << "DBGTT: Curved edge/face in el. " << (*itBE)->getNum() << "\n";
-// if (baseVert.find(*itBE) != baseVert.end()) { // If more than 1 curved edge/face in bad el., drop it
-// baseVert.erase(*itBE);
-// normals.erase(*itBE);
-// std::cout << "DBGTT: More than 1 curved edge/face detected in el. " << (*itBE)->getNum() << "\n";
-// break;
-// }
-// baseType[*itBE] = type;
-// baseVert[*itBE] = bv;
-// const SPoint3 p0 = bv[0]->point(), pC = (*itBE)->barycenter(true);
-// SVector3 p0C(p0,pC);
-// if (dot(n,p0C) < 0.) n *= -1.; // Make straight normal in-going
-// for (int iV=0; iV<numPrimVert; iV++) // Compute normals to prim. vert. of edge/face
-// normals[*itBE].push_back(getNormalEdge(bv[iV],n,vertex2elements));
-// }
-// }
-// }
-//
-//}
-
-
-
// Detect whether edge/face is curved, and give normal
bool isCurvedAndNormal(int type, int order, const std::vector<MVertex*> &faceVert,
const SPoint3 pBar, SVector3 &normal, double &maxDist) {
@@ -431,9 +217,9 @@ bool isCurvedAndNormal(int type, int order, const std::vector<MVertex*> &faceVer
maxDist = std::max(maxDist,fabs(normalDisp));
}
- std::cout << "DBGTT: v0 is " << faceVert[0]->getNum() << ", v1 is " << faceVert[1]->getNum()
- << ", v2 is " << faceVert[2]->getNum() << ", maxDist = " << maxDist
- << ", scale = " << scale << ", test = " << (maxDist > eps*scale) << "\n";
+// std::cout << "DBGTT: v0 is " << faceVert[0]->getNum() << ", v1 is " << faceVert[1]->getNum()
+// << ", v2 is " << faceVert[2]->getNum() << ", maxDist = " << maxDist
+// << ", scale = " << scale << ", test = " << (maxDist > eps*scale) << "\n";
return (maxDist > eps*scale);
}
@@ -480,7 +266,7 @@ bool getCurvedFace(MElement* badEl,
double maxDist; // TOFIX: Max. normal. dist. to straight in curved face
const SPoint3 pBar = badEl->barycenter(true); // Bary. of el. to help orienting normal
if (isCurvedAndNormal(type,order,bv,pBar,n,maxDist)) { // Check whether edge/face is curved and compute straight normal
- std::cout << "DBGTT: Curved edge/face in el. " << badEl->getNum() << "\n";
+// std::cout << "DBGTT: Curved edge/face in el. " << badEl->getNum() << "\n";
if (found) { // If more than 1 curved edge/face in bad el., drop it
std::cout << "DBGTT: More than 1 curved edge/face detected in el. " << badEl->getNum() << "\n";
return false;
@@ -563,16 +349,8 @@ void curveMesh(std::map<MVertex*, std::vector<MElement *> > &vertex2elements,
const int nbBadElts = badElements.size();
std::vector<MElement*> badElts;
-// std::vector<int> types;
-// std::vector<std::vector<MVertex*> > faces;
-// std::vector<std::vector<SVector3> > &normals;
-// std::vector<double> &maxDist;
std::vector<SuperEl*> superElts;
badElts.reserve(nbBadElts);
-// types.reserve(nbBadElts);
-// faces.reserve(nbBadElts);
-// normals.reserve(nbBadElts);
-// maxDist.reserve(nbBadElts);
superElts.reserve(nbBadElts);
std::ofstream of("dum.pos");
@@ -585,13 +363,7 @@ void curveMesh(std::map<MVertex*, std::vector<MElement *> > &vertex2elements,
double maxDistBE;
if (getCurvedFace(*itBE,vertex2elements,faceBE,normalsBE,typeBE,maxDistBE)) {
badElts.push_back(*itBE);
-// types.push_back(typeBE);
-// faces.push_back(faceBE);
-// normals.push_back(normalsBE);
-// maxDist.push_back(maxDistBE);
superElts.push_back(new SuperEl(*itBE,maxDistBE*p.distanceFactor,typeBE,faceBE,normalsBE));
-// superElts.push_back(new SuperEl(*itBE,distMaxStraight(*itBE)*p.distanceFactor,typeBE,faceBE,normalsBE));
-// movedVert.insert(faceBE.begin(),faceBE.end()); // Do not touch vert. of curved face
of << superElts.back()->printPOS();
}
}
@@ -601,26 +373,28 @@ void curveMesh(std::map<MVertex*, std::vector<MElement *> > &vertex2elements,
std::set<MVertex*> movedVert;
for (int iBE=0; iBE<badElts.size(); ++iBE) {
-// std::set<MElement*> blob = getSurroundingBlob(badElts[iBE], p.nbLayers, vertex2elements, p.distanceFactor);
std::set<MElement*> blob = getSuperElBlob(badElts[iBE], vertex2elements, superElts[iBE]);
- std::cout << "DBGTT: Blob of bad el. " << badElts[iBE]->getNum() << " contains elts.";
- for (std::set<MElement*>::iterator itE = blob.begin(); itE != blob.end(); ++itE) std::cout << " " << (*itE)->getNum();
- std::cout << "\n";
+// std::cout << "DBGTT: Blob of bad el. " << badElts[iBE]->getNum() << " contains elts.";
+// for (std::set<MElement*>::iterator itE = blob.begin(); itE != blob.end(); ++itE) std::cout << " " << (*itE)->getNum();
+// std::cout << "\n";
makeStraight(badElts[iBE],movedVert); // Make bad. el. straight
- for (std::set<MElement*>::iterator itE = blob.begin(); itE != blob.end(); ++itE)
+ for (std::set<MElement*>::iterator itE = blob.begin(); itE != blob.end(); ++itE) {
+ makeStraight(*itE,movedVert);
for (int i = 0; i < (*itE)->getNumVertices(); ++i) { // For each vert. of each el. in blob
+// for (int i = (*itE)->getNumPrimaryVertices(); i < (*itE)->getNumVertices(); ++i) { // For each vert. of each el. in blob
MVertex* vert = (*itE)->getVertex(i);
if (movedVert.find(vert) == movedVert.end()) { // If vert. not already moved
double xyzS[3] = {vert->x(), vert->y(), vert->z()}, xyzC[3];
if (superElts[iBE]->straightToCurved(xyzS,xyzC)) {
- std::cout << "DBGTT: moving vertex " << vert->getNum() << " from (" << xyzS[0] << "," << xyzS[1] << "," << xyzS[2] << ") to (" << xyzC[0] << "," << xyzC[1] << "," << xyzC[2] << ")\n";
+// std::cout << "DBGTT: moving vertex " << vert->getNum() << " from (" << xyzS[0] << "," << xyzS[1] << "," << xyzS[2] << ") to (" << xyzC[0] << "," << xyzC[1] << "," << xyzC[2] << ")\n";
vert->setXYZ(xyzC[0],xyzC[1],xyzC[2]);
movedVert.insert(vert);
}
- else std::cout << "DBGTT: Failed to move vertex " << vert->getNum() << " with bad. el " << badElts[iBE]->getNum() << "\n";
+// else std::cout << "DBGTT: Failed to move vertex " << vert->getNum() << " with bad. el " << badElts[iBE]->getNum() << "\n";
}
// else std::cout << "DBGTT: Already moved vertex " << vert->getNum() << " with bad. el " << badElts[iBE]->getNum() << "\n";
}
+ }
}
}
@@ -631,60 +405,6 @@ void curveMesh(std::map<MVertex*, std::vector<MElement *> > &vertex2elements,
-//static void optimizeBL(std::map<MVertex*, std::vector<MElement *> > &vertex2elements,
-// std::set<MElement*> &badElements, OptHomParameters &p, int samples)
-//{
-//
-// const int nbBadElts = badElements.size();
-// std::vector<MElement*> &badElts();
-// std::map<MElement*,int> types; // Type (LIN, TRI, QUA) of curved face in bad. el.
-// std::map<MElement*,std::vector<MVertex*> > faces; // Vertices of curved face in bad. el.
-// std::map<MElement*,std::vector<SVector3> > normals; // Normal to each prim. vertex of curved face in bad. el.
-// getBaseVertsAndNormals(badElements,vertex2elements,types,faces,normals);
-//
-// std::ofstream of("dum.pos");
-// of << "View \" \"{\n";
-//
-// for (std::set<MElement*>::const_iterator itBE = badElements.begin(); itBE != badElements.end(); ++itBE) {
-//
-// if (faces.find(*itBE) == faces.end()) continue; // No base edge/face -> no super el. for this bad. el.
-//
-// SuperEl se(*itBE,distMaxStraight(*itBE)*p.distanceFactor,types[*itBE],faces[*itBE],normals[*itBE]);
-//
-// std::set<MElement*> blob = getSurroundingBlob(*itBE, p.nbLayers, vertex2elements, p.distanceFactor);
-// std::set<MVertex*> toMove;
-// for (std::set<MElement*>::iterator itE = blob.begin(); itE != blob.end(); ++itE)
-// for (int i = 0; i < (*itE)->getNumVertices(); ++i) toMove.insert((*itE)->getVertex(i));
-//
-//// std::cout << "DBGTT: Before makeStraight:\n";
-//// se.printCoord();
-// makeStraight(*itBE); // Make bad. el. straight
-//// for (std::set<MElement*>::iterator itE = blob.begin(); itE != blob.end(); ++itE) makeStraight(*itE);
-//// std::cout << "DBGTT: After makeStraight:\n";
-//// se.printCoord();
-//
-// // Apply transformation straight -> curved to all vertices in super el.
-// for (std::set<MVertex*>::iterator itV = toMove.begin(); itV != toMove.end(); ++itV) {
-// double xyzS[3] = {(*itV)->x(), (*itV)->y(), (*itV)->z()}, xyzC[3];
-// if (se.straightToCurved(xyzS,xyzC)) {
-//// std::cout << "DBGTT: moving vertex " << (*itV)->getNum() << " from (" << xyzS[0] << "," << xyzS[1] << "," << xyzS[2] << ") to (" << xyzC[0] << "," << xyzC[1] << "," << xyzC[2] << ")\n";
-// (*itV)->setXYZ(xyzC[0],xyzC[1],xyzC[2]);
-// }
-// else std::cout << "DBGTT: Failed to move vertex " << (*itV)->getNum() << " with bad. el " << (*itBE)->getNum() << "\n";
-// }
-//
-//// std::cout << se.printPOS();
-// of << se.printPOS();
-//
-// }
-//
-// of << "};\n";
-// of.close();
-//
-//}
-
-
-
void HighOrderMeshFastCurving(GModel *gm, FastCurvingParameters &p)
{
diff --git a/contrib/HighOrderMeshOptimizer/SuperEl.cpp b/contrib/HighOrderMeshOptimizer/SuperEl.cpp
index 81e8a00bcf..76a4bece92 100644
--- a/contrib/HighOrderMeshOptimizer/SuperEl.cpp
+++ b/contrib/HighOrderMeshOptimizer/SuperEl.cpp
@@ -146,7 +146,7 @@ void SuperEl::createSuperElPrism(MElement *badEl, double dist, const std::vector
const int Nv = prismBasis->getNumShapeFunctions(); // Number of vertices in HO prism
// Store/create all vertices for super-element
- if (badEl->getTypeForMSH() == MSH_PRI_18) {
+ if (p == 2) {
_superVert.resize(18);
_superVert[0] = v0; // First-order vertices
_superVert[1] = v1;
@@ -167,7 +167,7 @@ void SuperEl::createSuperElPrism(MElement *badEl, double dist, const std::vector
_superEl = new MPrism18(_superVert);
}
else {
- std::cout << "ERROR: Type of prism not supported for bad. el. " << badEl->getNum() << "\n";
+ std::cout << "ERROR: Prismatic superEl. of order " << p << " not supported for bad. el. " << badEl->getNum() << "\n";
_superEl = 0;
_superEl0 = 0;
return;
@@ -212,7 +212,7 @@ void SuperEl::createSuperElHex(MElement *badEl, double dist, const std::vector<M
const int Nv = prismBasis->getNumShapeFunctions(); // Number of vertices in HO hex
// Store/create all vertices for super-element
- if (badEl->getTypeForMSH() == MSH_HEX_27) {
+ if (p == 2) {
_superVert.resize(27);
_superVert[0] = v0; // First-order vertices
_superVert[1] = v1;
@@ -237,7 +237,7 @@ void SuperEl::createSuperElHex(MElement *badEl, double dist, const std::vector<M
_superEl = new MHexahedron27(_superVert);
}
else {
- std::cout << "ERROR: Type of prism not supported for bad. el. " << badEl->getNum() << "\n";
+ std::cout << "ERROR: Hex. superEl. of order " << p << " not supported for bad. el. " << badEl->getNum() << "\n";
_superEl = 0;
_superEl0 = 0;
return;
@@ -277,30 +277,6 @@ bool SuperEl::straightToCurved(double *xyzS, double *xyzC) const {
-//std::string SuperEl::printPOS() {
-//
-// std::vector<MVertex*> verts;
-// _superEl0->getVertices(verts);
-// std::string posStr(_superEl0->getStringForPOS());
-// int n = _superEl0->getNumVertices();
-//
-// std::ostringstream oss;
-//
-// oss << posStr << "(";
-// for(int i = 0; i < n; i++){
-// if(i) oss << ",";
-// oss << _superVert[i]->x() << "," << _superVert[i]->y() << "," << _superVert[i]->z();
-// }
-// oss << "){0";
-// for(int i = 0; i < n; i++) oss << ",0.";
-// oss << "};\n";
-//
-// return oss.str();
-//
-//}
-
-
-
std::string SuperEl::printPOS() {
std::vector<MVertex*> verts;
--
GitLab