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

Added local surface parametrization for high-order mesh optimization, worked... 

Added local surface parametrization for high-order mesh optimization, worked on one-by-one strategy
parent abf1bf29
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
contrib/HighOrderMeshOptimizer/OptHomMesh.cpp
+ 1
− 1
View file @ 1d3cc55c
......@@ -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;
}
......
contrib/HighOrderMeshOptimizer/OptHomRun.cpp
+ 96
− 95
View file @ 1d3cc55c
......@@ -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
......
contrib/HighOrderMeshOptimizer/ParamCoord.cpp
+ 94
− 9
View file @ 1d3cc55c
......@@ -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);
}
}
contrib/HighOrderMeshOptimizer/ParamCoord.h
+ 32
− 18
View file @ 1d3cc55c
......@@ -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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment