From 1d3cc55c15dba6b127ebe613a144a64d0757fec4 Mon Sep 17 00:00:00 2001
From: Thomas Toulorge <thomas.toulorge@mines-paristech.fr>
Date: Fri, 18 Oct 2013 17:39:40 +0000
Subject: [PATCH] Added local surface parametrization for high-order mesh
optimization, worked on one-by-one strategy
---
contrib/HighOrderMeshOptimizer/OptHomMesh.cpp | 2 +-
contrib/HighOrderMeshOptimizer/OptHomRun.cpp | 191 +++++++++---------
contrib/HighOrderMeshOptimizer/ParamCoord.cpp | 103 +++++++++-
contrib/HighOrderMeshOptimizer/ParamCoord.h | 50 +++--
4 files changed, 223 insertions(+), 123 deletions(-)
diff --git a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
index 4b81c780dd..e24731670d 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
@@ -73,7 +73,7 @@ Mesh::Mesh(const std::map<MElement*,GEntity*> &element2entity,
if (vDim == 3) param = new ParamCoordPhys3D();
else if (hasParam) param = new ParamCoordParent(vert);
else {
- if (vDim == 2) param = new ParamCoordPhys2D(); //todo: make 2d local surf. param
+ if (vDim == 2) param = new ParamCoordLocalSurf(vert);
else param = new ParamCoordLocalLine(vert);
nonGeoMove = true;
}
diff --git a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
index 8c0b20135c..7a013dec96 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
@@ -499,111 +499,112 @@ static void optimizeOneByOne
for (int iterBlob=0; iterBlob<p.maxAdaptBlob; iterBlob++) {
-// OptHOM *opt;
-//
-// // First step: small blob with unsafe optimization (only 1st-order
-// // bnd. vertices fixed)
-// std::set<MElement*> toOptimizePrim = getSurroundingBlob
-// (worstEl, nbLayers, vertex2elements, distanceFactor, 0);
-// std::set<MVertex*> toFix1 = getPrimBndVertices(toOptimizePrim, vertex2elements);
-// std::set<MElement*> toOptimize1;
-// std::set_difference(toOptimizePrim.begin(),toOptimizePrim.end(),
-// badElts.begin(),badElts.end(), // Do not optimize badElts
-// std::inserter(toOptimize1, toOptimize1.end()));
-// Msg::Info("Optimizing primary blob %i (max. %i remaining) composed of"
-// " %4d elements", iBadEl, badElts.size(), toOptimize1.size());
-// fflush(stdout);
-// opt = new OptHOM(toOptimize1, toFix1, p.fixBndNodes);
-// //std::ostringstream ossI1;
-// //ossI1 << "initial_primary_" << iter << ".msh";
-// //opt->mesh.writeMSH(ossI1.str().c_str());
-// success = opt->optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN, p.BARRIER_MAX,
-// false, samples, p.itMax, p.optPassMax);
-//
-// // Second step: add external layer, check it and optimize it safely (all
-// // bnd. vertices fixed) if new broken element
-// if (success > 0) {
-// opt->mesh.updateGEntityPositions();
-// std::set<MElement*> layer = addBlobLayer(toOptimizePrim, vertex2elements);
-// if (detectNewBrokenElement(layer, badElts, p)) {
-// delete opt;
-// std::set<MVertex*> toFix2 = getAllBndVertices(toOptimizePrim, vertex2elements);
-// std::set<MElement*> toOptimize2;
-// std::set_difference(toOptimizePrim.begin(),toOptimizePrim.end(),
-// badElts.begin(),badElts.end(),
-// std::inserter(toOptimize2, toOptimize2.end()));
-// Msg::Info("Optimizing corrective blob %i (max. %i remaining) "
-// "composed of %4d elements", iBadEl, badElts.size(),
-// toOptimize2.size());
-// fflush(stdout);
-// opt = new OptHOM(toOptimize2, toFix2, p.fixBndNodes);
-// //std::ostringstream ossI1;
-// //ossI1 << "initial_corrective_" << iter << ".msh";
-// //opt->mesh.writeMSH(ossI1.str().c_str());
-// success = opt->optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN,
-// p.BARRIER_MAX, false, samples, p.itMax, p.optPassMax);
-// if (success >= 0 && p.BARRIER_MIN_METRIC > 0) {
-// Msg::Info("Jacobian optimization succeed, starting svd optimization");
-// success = opt->optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN_METRIC,
-// p.BARRIER_MAX, true, samples, p.itMax, p.optPassMax);
-// }
-// }
-// else {
-// Msg::Info("Primary blob %i did not break new elements, "
-// "no correction needed", iBadEl);
-// fflush(stdout);
-// }
-// }
+ OptHOM *opt;
+ // First step: small blob with unsafe optimization (only 1st-order
+ // bnd. vertices fixed)
std::set<MElement*> toOptimizePrim = getSurroundingBlob
- (worstEl, nbLayers, vertex2elements, distanceFactor, 1, p.optPrimSurfMesh);
-// std::set<MElement*> layer = addBlobLayer(toOptimizePrim, vertex2elements);
- std::set<MVertex*> toFix = getAllBndVertices(toOptimizePrim, vertex2elements);
- std::set<MElement*> toOptimize;
+ (worstEl, nbLayers, vertex2elements, distanceFactor, 0, p.optPrimSurfMesh);
+ std::set<MVertex*> toFix1 = getPrimBndVertices(toOptimizePrim, vertex2elements);
+ std::set<MElement*> toOptimize1;
std::set_difference(toOptimizePrim.begin(),toOptimizePrim.end(),
- badElts.begin(),badElts.end(),
- std::inserter(toOptimize, toOptimize.end()));
- Msg::Info("Optimizing blob %i (max. %i remaining) "
- "composed of %4d elements", iBadEl, badElts.size(),
- toOptimize.size());
+ badElts.begin(),badElts.end(), // Do not optimize badElts
+ std::inserter(toOptimize1, toOptimize1.end()));
+ Msg::Info("Optimizing primary blob %i (max. %i remaining) composed of"
+ " %4d elements", iBadEl, badElts.size(), toOptimize1.size());
fflush(stdout);
- OptHOM *opt = new OptHOM(element2entity, toOptimize, toFix, p.fixBndNodes);
+ opt = new OptHOM(element2entity, toOptimize1, toFix1, p.fixBndNodes);
//std::ostringstream ossI1;
- //ossI1 << "initial_corrective_" << iter << ".msh";
+ //ossI1 << "initial_primary_" << iter << ".msh";
//opt->mesh.writeMSH(ossI1.str().c_str());
- success = opt->optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN,
- p.BARRIER_MAX, false, samples, p.itMax, p.optPassMax);
- if (success >= 0 && p.BARRIER_MIN_METRIC > 0) {
- Msg::Info("Jacobian optimization succeed, starting svd optimization");
- success = opt->optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN_METRIC,
- p.BARRIER_MAX, true, samples, p.itMax, p.optPassMax);
- }
+ success = opt->optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN, p.BARRIER_MAX,
+ false, samples, p.itMax, p.optPassMax);
- // Measure min and max Jac., update mesh
- if ((success > 0) || (iterBlob == p.maxAdaptBlob-1)) {
- double minJac, maxJac, distMaxBND, distAvgBND;
- opt->recalcJacDist();
- opt->getJacDist(minJac, maxJac, distMaxBND, distAvgBND);
- p.minJac = std::min(p.minJac,minJac);
- p.maxJac = std::max(p.maxJac,maxJac);
+ // Second step: add external layer, check it and optimize it safely (all
+ // bnd. vertices fixed) if new broken element
+ if (success > 0) {
opt->mesh.updateGEntityPositions();
+ std::set<MElement*> layer = addBlobLayer(toOptimizePrim, vertex2elements);
+ if (detectNewBrokenElement(layer, badElts, p)) {
+ delete opt;
+ std::set<MVertex*> toFix2 = getAllBndVertices(toOptimizePrim, vertex2elements);
+ std::set<MElement*> toOptimize2;
+ std::set_difference(toOptimizePrim.begin(),toOptimizePrim.end(),
+ badElts.begin(),badElts.end(),
+ std::inserter(toOptimize2, toOptimize2.end()));
+ Msg::Info("Optimizing corrective blob %i (max. %i remaining) "
+ "composed of %4d elements", iBadEl, badElts.size(),
+ toOptimize2.size());
+ fflush(stdout);
+ opt = new OptHOM(element2entity, toOptimize2, toFix2, p.fixBndNodes);
+ //std::ostringstream ossI1;
+ //ossI1 << "initial_corrective_" << iter << ".msh";
+ //opt->mesh.writeMSH(ossI1.str().c_str());
+ success = opt->optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN,
+ p.BARRIER_MAX, false, samples, p.itMax, p.optPassMax);
+ if (success >= 0 && p.BARRIER_MIN_METRIC > 0) {
+ Msg::Info("Jacobian optimization succeed, starting svd optimization");
+ success = opt->optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN_METRIC,
+ p.BARRIER_MAX, true, samples, p.itMax, p.optPassMax);
+ }
+ }
+ else {
+ Msg::Info("Primary blob %i did not break new elements, "
+ "no correction needed", iBadEl);
+ fflush(stdout);
+ }
}
- //std::ostringstream ossI2;
- //ossI2 << "final_ITER_" << iter << ".msh";
- //temp.mesh.writeMSH(ossI2.str().c_str());
- if (success <= 0) {
- distanceFactor *= p.adaptBlobDistFact;
- nbLayers *= p.adaptBlobLayerFact;
- Msg::Info("Blob %i failed (adapt #%i), adapting with increased size",
- iBadEl, iterBlob);
-// if (iterBlob == p.maxAdaptBlob-1) {
- std::ostringstream ossI2;
- ossI2 << "final_" << iBadEl << ".msh";
- opt->mesh.writeMSH(ossI2.str().c_str());
-// }
- }
- else break;
+// std::set<MElement*> toOptimizePrim = getSurroundingBlob
+// (worstEl, nbLayers, vertex2elements, distanceFactor, 1, p.optPrimSurfMesh);
+//// std::set<MElement*> layer = addBlobLayer(toOptimizePrim, vertex2elements);
+// std::set<MVertex*> toFix = getAllBndVertices(toOptimizePrim, vertex2elements);
+// std::set<MElement*> toOptimize;
+// std::set_difference(toOptimizePrim.begin(),toOptimizePrim.end(),
+// badElts.begin(),badElts.end(),
+// std::inserter(toOptimize, toOptimize.end()));
+// Msg::Info("Optimizing blob %i (max. %i remaining) "
+// "composed of %4d elements", iBadEl, badElts.size(),
+// toOptimize.size());
+// fflush(stdout);
+// OptHOM *opt = new OptHOM(element2entity, toOptimize, toFix, p.fixBndNodes);
+// //std::ostringstream ossI1;
+// //ossI1 << "initial_corrective_" << iter << ".msh";
+// //opt->mesh.writeMSH(ossI1.str().c_str());
+// success = opt->optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN,
+// p.BARRIER_MAX, false, samples, p.itMax, p.optPassMax);
+// if (success >= 0 && p.BARRIER_MIN_METRIC > 0) {
+// Msg::Info("Jacobian optimization succeed, starting svd optimization");
+// success = opt->optimize(p.weightFixed, p.weightFree, p.BARRIER_MIN_METRIC,
+// p.BARRIER_MAX, true, samples, p.itMax, p.optPassMax);
+// }
+//
+// // Measure min and max Jac., update mesh
+// if ((success > 0) || (iterBlob == p.maxAdaptBlob-1)) {
+// double minJac, maxJac, distMaxBND, distAvgBND;
+// opt->recalcJacDist();
+// opt->getJacDist(minJac, maxJac, distMaxBND, distAvgBND);
+// p.minJac = std::min(p.minJac,minJac);
+// p.maxJac = std::max(p.maxJac,maxJac);
+// opt->mesh.updateGEntityPositions();
+// }
+//
+// //std::ostringstream ossI2;
+// //ossI2 << "final_ITER_" << iter << ".msh";
+// //temp.mesh.writeMSH(ossI2.str().c_str());
+// if (success <= 0) {
+// distanceFactor *= p.adaptBlobDistFact;
+// nbLayers *= p.adaptBlobLayerFact;
+// Msg::Info("Blob %i failed (adapt #%i), adapting with increased size",
+// iBadEl, iterBlob);
+//// if (iterBlob == p.maxAdaptBlob-1) {
+// std::ostringstream ossI2;
+// ossI2 << "final_" << iBadEl << ".msh";
+// opt->mesh.writeMSH(ossI2.str().c_str());
+//// }
+// }
+// else break;
+
}
//#pragma omp critical
diff --git a/contrib/HighOrderMeshOptimizer/ParamCoord.cpp b/contrib/HighOrderMeshOptimizer/ParamCoord.cpp
index 1361610e1d..97ce05ee4a 100644
--- a/contrib/HighOrderMeshOptimizer/ParamCoord.cpp
+++ b/contrib/HighOrderMeshOptimizer/ParamCoord.cpp
@@ -117,21 +117,106 @@ void ParamCoordParent::gXyz2gUvw(const SPoint3 &uvw,
}
-ParamCoordLocalLine::ParamCoordLocalLine(MVertex* v) : dir(0.)
+namespace {
+
+SVector3 getLineElTangent(MElement *el, int iNode) {
+
+ double gsf[1256][3], u, v, w;
+ el->getNode(iNode,u,v,w);
+// el->getGradShapeFunctions(u,v,w,gsf);
+ el->getGradShapeFunctions(u,v,w,gsf,1);
+
+ SVector3 dxyzdu(0.);
+// int nSF = el->getNumShapeFunctions()();
+ int nSF = el->getNumPrimaryVertices();
+ for (int j=0; j<nSF; j++) {
+ const SPoint3 p = el->getVertex(j)->point();
+ dxyzdu(0) += gsf[j][0]*p.x();
+ dxyzdu(1) += gsf[j][0]*p.y();
+ dxyzdu(2) += gsf[j][0]*p.z();
+ }
+ dxyzdu.normalize();
+
+ return dxyzdu;
+
+}
+
+SVector3 getSurfElNormal(MElement *el, int iNode) {
+
+ double gsf[1256][3], u, v, w;
+ el->getNode(iNode,u,v,w);
+// el->getGradShapeFunctions(u,v,w,gsf);
+ el->getGradShapeFunctions(u,v,w,gsf,1);
+
+ SVector3 dxyzdu(0.), dxyzdv(0.);
+// int nSF = el->getNumShapeFunctions()();
+ int nSF = el->getNumPrimaryVertices();
+ for (int j=0; j<nSF; j++) {
+ const SPoint3 p = el->getVertex(j)->point();
+ dxyzdu(0) += gsf[j][0]*p.x();
+ dxyzdu(1) += gsf[j][0]*p.y();
+ dxyzdu(2) += gsf[j][0]*p.z();
+ dxyzdv(0) += gsf[j][1]*p.x();
+ dxyzdv(1) += gsf[j][1]*p.y();
+ dxyzdv(2) += gsf[j][1]*p.z();
+ }
+
+ SVector3 normal = crossprod(dxyzdu,dxyzdv);
+ normal.normalize();
+ return normal;
+
+}
+
+}
+
+ParamCoordLocalLine::ParamCoordLocalLine(MVertex* v) :
+ dir(0.), x0(v->x()), y0(v->y()), z0(v->z())
{
- GEdge *ge = static_cast<GEdge*>(v->onWhat());
+ GEntity *ge = v->onWhat();
+ const unsigned nEl = ge->getNumMeshElements();
- for (std::vector<MLine*>::iterator it = ge->lines.begin(); it != ge->lines.end(); it++) {
+ for (unsigned iEl = 0; iEl < nEl; iEl++) {
+ MElement *el = ge->getMeshElement(iEl);
std::vector<MVertex*> lVerts;
- (*it)->getVertices(lVerts);
- if (std::find(lVerts.begin(),lVerts.end(),v) == lVerts.end()) {
- SVector3 tan(lVerts[0]->point(),lVerts[1]->point());
- tan.normalize();
- dir += tan;
+ el->getVertices(lVerts);
+ std::vector<MVertex*>::iterator itV = std::find(lVerts.begin(),lVerts.end(),v);
+ if (itV != lVerts.end()) {
+ const int iNode = std::distance(lVerts.begin(),itV);
+ dir += getLineElTangent(el,iNode);
}
}
-
dir.normalize();
}
+
+ParamCoordLocalSurf::ParamCoordLocalSurf(MVertex* v) : x0(v->x()), y0(v->y()), z0(v->z())
+{
+
+ GEntity *ge = v->onWhat();
+ const unsigned nEl = ge->getNumMeshElements();
+
+ SVector3 n(0.);
+ for (unsigned iEl = 0; iEl < nEl; iEl++) {
+ MElement *el = ge->getMeshElement(iEl);
+ std::vector<MVertex*> lVerts;
+ el->getVertices(lVerts);
+ std::vector<MVertex*>::iterator itV = std::find(lVerts.begin(),lVerts.end(),v);
+ if (itV != lVerts.end()) {
+ const int iNode = std::distance(lVerts.begin(),itV);
+ n += getSurfElNormal(el,iNode);
+ }
+ }
+ n.normalize();
+
+ if (fabs(fabs(dot(n,SVector3(1.,0.,0.)))-1.) < 1.e-10) { // If normal is x-axis, take y- and z- axis as dir.
+ dir0 = SVector3(0.,1.,0.);
+ dir1 = SVector3(0.,0.,1.);
+ }
+ else {
+ dir0 = SVector3(1.-n.x()*n.x(),-n.x()*n.y(),-n.x()*n.z()); // 1st dir. = (normalized) proj. of e_x in plane
+ dir0.normalize();
+ dir1 = crossprod(dir0,n);
+ }
+
+}
diff --git a/contrib/HighOrderMeshOptimizer/ParamCoord.h b/contrib/HighOrderMeshOptimizer/ParamCoord.h
index 3c509cb83e..a0bb4f597e 100644
--- a/contrib/HighOrderMeshOptimizer/ParamCoord.h
+++ b/contrib/HighOrderMeshOptimizer/ParamCoord.h
@@ -65,23 +65,6 @@ public:
}
};
-class ParamCoordPhys2D : public ParamCoord
-{
-public:
- SPoint3 getUvw(MVertex* v) { return v->point(); }
- SPoint3 uvw2Xyz(const SPoint3 &uvw) { return uvw; }
- void gXyz2gUvw(const SPoint3 &uvw, const SPoint3 &gXyz, SPoint3 &gUvw) { gUvw = gXyz; }
- void gXyz2gUvw(const SPoint3 &uvw, const std::vector<SPoint3> &gXyz, std::vector<SPoint3> &gUvw)
- {
- std::vector<SPoint3>::iterator itUvw=gUvw.begin();
- for (std::vector<SPoint3>::const_iterator itXyz=gXyz.begin(); itXyz != gXyz.end();
- itXyz++) {
- *itUvw = *itXyz;
- itUvw++;
- }
- }
-};
-
class ParamCoordParent : public ParamCoord
{
public:
@@ -102,7 +85,9 @@ class ParamCoordLocalLine : public ParamCoord
public:
ParamCoordLocalLine(MVertex* v);
SPoint3 getUvw(MVertex* v) { return SPoint3(0.,0.,0.); }
- SPoint3 uvw2Xyz(const SPoint3 &uvw) { return SPoint3(uvw[0]*dir[0],uvw[0]*dir[1],uvw[0]*dir[2]); }
+ SPoint3 uvw2Xyz(const SPoint3 &uvw) {
+ return SPoint3(x0+uvw[0]*dir[0],y0+uvw[0]*dir[1],z0+uvw[0]*dir[2]);
+ }
void gXyz2gUvw(const SPoint3 &uvw, const SPoint3 &gXyz, SPoint3 &gUvw) {
gUvw[0] = gXyz.x()*dir[0] + gXyz.y()*dir[1] + gXyz.z()*dir[2];
}
@@ -115,7 +100,36 @@ public:
}
}
protected:
+ double x0, y0, z0;
SVector3 dir;
};
+class ParamCoordLocalSurf : public ParamCoord
+{
+public:
+ ParamCoordLocalSurf(MVertex* v);
+ SPoint3 getUvw(MVertex* v) { return SPoint3(0.,0.,0.); }
+ SPoint3 uvw2Xyz(const SPoint3 &uvw) {
+ return SPoint3(x0+uvw[0]*dir0[0]+uvw[1]*dir1[0],
+ y0+uvw[0]*dir0[1]+uvw[1]*dir1[1],
+ z0+uvw[0]*dir0[2]+uvw[1]*dir1[2]);
+ }
+ void gXyz2gUvw(const SPoint3 &uvw, const SPoint3 &gXyz, SPoint3 &gUvw) {
+ gUvw[0] = gXyz.x()*dir0[0] + gXyz.y()*dir0[1] + gXyz.z()*dir0[2];
+ gUvw[1] = gXyz.x()*dir1[0] + gXyz.y()*dir1[1] + gXyz.z()*dir1[2];
+ }
+ void gXyz2gUvw(const SPoint3 &uvw, const std::vector<SPoint3> &gXyz, std::vector<SPoint3> &gUvw) {
+ std::vector<SPoint3>::iterator itUvw = gUvw.begin();
+ for (std::vector<SPoint3>::const_iterator itXyz=gXyz.begin();
+ itXyz != gXyz.end(); itXyz++) {
+ (*itUvw)[0] = itXyz->x()*dir0[0] + itXyz->y()*dir0[1] + itXyz->z()*dir0[2];
+ (*itUvw)[1] = itXyz->x()*dir1[0] + itXyz->y()*dir1[1] + itXyz->z()*dir1[2];
+ itUvw++;
+ }
+ }
+protected:
+ double x0, y0, z0;
+ SVector3 dir0, dir1;
+};
+
#endif
--
GitLab