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

Cleaned up and fixed bugs in high-order fast curving tool

parent ad86e02b
No related branches found
No related tags found
No related merge requests found
Show whitespace changes
Inline Side-by-side
contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp
+ 13
− 293
View file @ 2b8b05fe
......@@ -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,87 +373,35 @@ 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";
}
}
}
}
//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();
//
//}
}
......
contrib/HighOrderMeshOptimizer/SuperEl.cpp
+ 4
− 28
View file @ 2b8b05fe
......@@ -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;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment