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

Improved fast high-order BL mesh curving (better meta-element, geometrical optimization for 2D P2)

parent cc094aa7
Branches
Tags
No related merge requests found
Hide whitespace changes
Inline Side-by-side
contrib/HighOrderMeshOptimizer/MetaEl.cpp
+ 224
− 92
View file @ e339c0c2
...@@ -42,70 +42,136 @@ ...@@ -42,70 +42,136 @@
std::map<int, MetaEl::metaInfoType> MetaEl::_metaInfo; std::map<int, MetaEl::metaInfoType> MetaEl::_metaInfo;
//namespace {
//
//
//void getEdgeIndQuad(int order, std::vector<int> edInd) {
// edInd.clear();
// const int nbEdVert = 4*(order-1);
// edInd.resize(nbEdVert);
// for (int iV = 0; iV < nbEdVert; iV++) edInd[iV] = iV;
//}
//
//
//
//}
MetaEl::metaInfoType::metaInfoType(int type, int order) MetaEl::metaInfoType::metaInfoType(int type, int order)
{ {
static const double TOL = 1e-10;
typedef std::vector<int>::iterator VIntIter; typedef std::vector<int>::iterator VIntIter;
int iBaseFace = 0, iTopFace = 0, nbEdVert = 0; // Get basic info on dimension, faces and vertices
int iBaseFace = 0, iTopFace = 0, nbEdVert = 0, nbPrimTopVert = 0;
int baseFaceType = 0;
switch (type) { switch (type) {
case TYPE_QUA: case TYPE_QUA:
iBaseFace = 0; iTopFace = 2; nbEdVert = 4 + 4*(order-1); dim = 2;
iBaseFace = 0; iTopFace = 2;
nbEdVert = 4 + 4*(order-1); nbPrimTopVert = 2;
baseFaceType = TYPE_LIN;
break; break;
case TYPE_PRI: case TYPE_PRI:
dim = 3;
iBaseFace = 0; iTopFace = 1; nbEdVert = 6 + 9*(order-1); iBaseFace = 0; iTopFace = 1; nbEdVert = 6 + 9*(order-1);
nbPrimTopVert = 3;
baseFaceType = TYPE_TRI;
break; break;
case TYPE_HEX: case TYPE_HEX:
dim = 3;
iBaseFace = 0; iTopFace = 5; nbEdVert = 8 + 12*(order-1); iBaseFace = 0; iTopFace = 5; nbEdVert = 8 + 12*(order-1);
nbPrimTopVert = 4;
baseFaceType = TYPE_QUA;
break; break;
default: default:
Msg::Error("MetaEl not implemented for element of type %d", type); Msg::Error("MetaEl not implemented for element of type %d", type);
dim = 0;
nbVert = 0; nbVert = 0;
return; return;
} }
// Get HO nodal basis // Get HO nodal base
// const int tag = ElementType::getTag(type, order, true); // Get tag for incomplete element type
const int tag = ElementType::getTag(type, order); // Get tag for complete element type const int tag = ElementType::getTag(type, order); // Get tag for complete element type
if (!tag) return; if (!tag) return;
const nodalBasis *basis = BasisFactory::getNodalBasis(tag); const nodalBasis *basis = BasisFactory::getNodalBasis(tag);
nbVert = basis->getNumShapeFunctions(); nbVert = basis->getNumShapeFunctions();
points = basis->points; points = basis->points;
std::set<int> foundVerts;
// Vertices on base and top faces // Get nodal bases (HO and linear) on base face
const int baseFaceTag = ElementType::getTag(baseFaceType, order); // Get tag for base face basis
if (!baseFaceTag) return;
baseFaceBasis = BasisFactory::getNodalBasis(baseFaceTag);
const int nbBaseVert = baseFaceBasis->getNumShapeFunctions();
baseGradShapeFuncVal.resize(nbBaseVert, 3*nbBaseVert);
const int linBaseFaceTag = ElementType::getTag(baseFaceType, 1); // Get tag for base face linear basis
if (!linBaseFaceTag) return;
linBaseFaceBasis = BasisFactory::getNodalBasis(linBaseFaceTag);
const int nbLinBaseShapeFunc = linBaseFaceBasis->getNumShapeFunctions();
baseLinShapeFuncVal.resize(nbBaseVert, nbLinBaseShapeFunc);
// Compute value of nodal bases of base face at base vertices
const fullMatrix<double> basePoints = baseFaceBasis->getReferenceNodes();
for (int iV = 0; iV < baseFaceBasis->getNumShapeFunctions(); iV++) {
const double &u = basePoints(iV, 0), &v = basePoints(iV, 1);
double linShapeFunc[1256];
linBaseFaceBasis->f(u, v, 0., linShapeFunc);
for (int iSF = 0; iSF < nbLinBaseShapeFunc; iSF++)
baseLinShapeFuncVal(iV, iSF) = linShapeFunc[iSF];
double gradShapeFunc[1256][3];
baseFaceBasis->df(u, v, 0., gradShapeFunc);
for (int iSF = 0; iSF < nbBaseVert; iSF++) {
baseGradShapeFuncVal(iV, 3*iSF) = gradShapeFunc[iSF][0];
baseGradShapeFuncVal(iV, 3*iSF+1) = gradShapeFunc[iSF][1];
baseGradShapeFuncVal(iV, 3*iSF+2) = gradShapeFunc[iSF][2];
}
}
// Indices of vertices on base face
std::set<int> foundVerts;
baseInd = basis->getClosure(basis->getClosureId(iBaseFace, 1, 0)); baseInd = basis->getClosure(basis->getClosureId(iBaseFace, 1, 0));
topInd = basis->getClosure(basis->getClosureId(iTopFace, 0, 0));
foundVerts.insert(baseInd.begin(), baseInd.end()); foundVerts.insert(baseInd.begin(), baseInd.end());
// Indices of vertices on top face
topInd = basis->getClosure(basis->getClosureId(iTopFace, 0, 0));
foundVerts.insert(topInd.begin(), topInd.end()); foundVerts.insert(topInd.begin(), topInd.end());
freeTopInd = std::vector<int>(topInd.begin()+nbPrimTopVert, topInd.end());
freeTop2Base.resize(freeTopInd.size());
for(int iVFT = 0; iVFT < freeTopInd.size(); iVFT++) {
const int &indFT = freeTopInd[iVFT];
const double &uTop = points(indFT, 0);
const double &vTop = points(indFT, 1);
for (int iVB = 0; iVB < baseInd.size(); iVB++) { // Find corresponding base vertex
const int &indB = baseInd[iVB];
const double diffU = uTop - points(indB, 0);
const double diffV = (dim == 3) ? vTop - points(indB, 1) : 0.;
if (diffU*diffU+diffV*diffV < TOL) {
freeTop2Base[iVFT] = iVB;
break;
}
}
}
// Vertices on edges (excluding base and top faces) // Indices of vertices on edges (excluding base and top faces)
for (int iV = 0; iV < nbEdVert; iV++) for (int iV = 0; iV < nbEdVert; iV++)
if (foundVerts.find(iV) == foundVerts.end()) { if (foundVerts.find(iV) == foundVerts.end()) {
edgeInd.push_back(iV); edgeInd.push_back(iV);
foundVerts.insert(iV); foundVerts.insert(iV);
} }
// Remaining face and interior vertices // Indices of remaining face and interior vertices
for(int iV = 0; iV < nbVert; iV++) for(int iV = 0; iV < nbVert; iV++) {
if (foundVerts.find(iV) == foundVerts.end()) { if (foundVerts.find(iV) == foundVerts.end()) {
otherInd.push_back(iV); otherInd.push_back(iV);
foundVerts.insert(iV); foundVerts.insert(iV);
const double &uFree = points(iV, 0);
const double &vFree = points(iV, 1);
for (int iVB = 0; iVB < baseInd.size(); iVB++) { // Find corresponding base vertex
const int &indB = baseInd[iVB];
const double diffU = uFree - points(indB, 0);
const double diffV = (dim == 3) ? vFree - points(indB, 1) : 0.;
if (diffU*diffU+diffV*diffV < TOL) {
other2Base.push_back(iVB);
break;
}
}
for (int iVT = 0; iVT < baseInd.size(); iVT++) { // Find corresponding top vertex
const int &indT = topInd[iVT];
const double diffU = uFree - points(indT, 0);
const double diffV = (dim == 3) ? vFree - points(indT, 1) : 0.;
if (diffU*diffU+diffV*diffV < TOL) {
other2Top.push_back(iVT);
break;
}
}
} }
}
} }
...@@ -122,6 +188,44 @@ const MetaEl::metaInfoType &MetaEl::getMetaInfo(int elType, int order) ...@@ -122,6 +188,44 @@ const MetaEl::metaInfoType &MetaEl::getMetaInfo(int elType, int order)
void MetaEl::computeBaseNorm(const SVector3 &metaNorm,
const std::vector<MVertex*> &baseVert,
const std::vector<MVertex*> &topPrimVert,
std::vector<SVector3> &baseNorm)
{
const int nbFaceNodes = baseVert.size(), nbPrimFaceNodes = topPrimVert.size();
// Compute thickness at each primary vertex
std::vector<double> linThick(nbPrimFaceNodes);
for (int iV = 0; iV < nbPrimFaceNodes; iV++)
linThick[iV] = topPrimVert[iV]->distance(baseVert[iV]);
// Compute normal at each base vertex
std::vector<double> thick(nbFaceNodes, 0.);
baseNorm.resize(nbFaceNodes);
for (int iV = 0; iV < nbFaceNodes; iV++) {
for (int iSF = 0; iSF < nbPrimFaceNodes; iSF++)
thick[iV] += linThick[iSF] * _mInfo.baseLinShapeFuncVal(iV, iSF);
double jac[3][2] = {0., 0., 0., 0., 0., 0.};
for (int iSF = 0; iSF < nbFaceNodes; iSF++) {
MVertex *vert = baseVert[iSF];
const double xyz[3] = {vert->x(), vert->y(), vert->z()};
for (int iXYZ = 0; iXYZ < 3; iXYZ++) {
jac[iXYZ][0] += xyz[iXYZ] * _mInfo.baseGradShapeFuncVal(iV, 3*iSF);
jac[iXYZ][1] += xyz[iXYZ] * _mInfo.baseGradShapeFuncVal(iV, 3*iSF+1);
}
}
SVector3 dXYZdU(jac[0][0], jac[1][0], jac[2][0]);
SVector3 dXYZdV = (_mInfo.dim == 2) ? metaNorm :
SVector3(jac[0][1], jac[1][1], jac[2][1]);
baseNorm[iV] = crossprod(dXYZdV, dXYZdU);
baseNorm[iV].normalize();
baseNorm[iV] *= thick[iV];
}
}
MetaEl::MetaEl(int type, int order, const std::vector<MVertex*> &baseVert, MetaEl::MetaEl(int type, int order, const std::vector<MVertex*> &baseVert,
const std::vector<MVertex*> &topPrimVert) : const std::vector<MVertex*> &topPrimVert) :
_mInfo(getMetaInfo(type, order)), _metaEl(0), _metaEl0(0) _mInfo(getMetaInfo(type, order)), _metaEl(0), _metaEl0(0)
...@@ -137,26 +241,58 @@ MetaEl::MetaEl(int type, int order, const std::vector<MVertex*> &baseVert, ...@@ -137,26 +241,58 @@ MetaEl::MetaEl(int type, int order, const std::vector<MVertex*> &baseVert,
// Add copies of vertices in base & top faces (only first-order vertices for top face) // Add copies of vertices in base & top faces (only first-order vertices for top face)
_metaVert.resize(nbVert); _metaVert.resize(nbVert);
for (int iV = 0; iV < baseInd.size(); iV++) for (int iV = 0; iV < baseInd.size(); iV++) {
_metaVert[baseInd[iV]] = new MVertex(*baseVert[iV]); const MVertex* const &bVert = baseVert[iV];
GEntity *geomEnt = bVert->onWhat();
if (geomEnt->dim() == 0)
_metaVert[baseInd[iV]] = new MVertex(*bVert);
else if (geomEnt->dim() == 1) {
double u;
bVert->getParameter(0, u);
_metaVert[baseInd[iV]] = new MEdgeVertex(bVert->x(), bVert->y(),
bVert->z(), geomEnt, u);
}
else if (geomEnt->dim() == 2) {
double u, v;
bVert->getParameter(0, u);
bVert->getParameter(1, v);
_metaVert[baseInd[iV]] = new MEdgeVertex(bVert->x(), bVert->y(),
bVert->z(), geomEnt, u, v);
}
}
for (int iV = 0; iV < topPrimVert.size(); iV++) for (int iV = 0; iV < topPrimVert.size(); iV++)
_metaVert[topInd[iV]] = new MVertex(*topPrimVert[iV]); _metaVert[topInd[iV]] = new MVertex(*topPrimVert[iV]);
// Create first-order meta-element // Create first-order meta-element and normals to base face
int faceType;
SVector3 metaNorm;
switch (type) { switch (type) {
case TYPE_QUA: case TYPE_QUA: {
_metaEl0 = new MQuadrangle(_metaVert); _metaEl0 = new MQuadrangle(_metaVert);
SVector3 v01(_metaVert[0]->point(), _metaVert[1]->point());
SVector3 v03(_metaVert[0]->point(), _metaVert[3]->point());
metaNorm = crossprod(v01, v03);
faceType = TYPE_LIN;
break; break;
case TYPE_PRI: }
case TYPE_PRI: {
_metaEl0 = new MPrism(_metaVert); _metaEl0 = new MPrism(_metaVert);
metaNorm = SVector3(0.);
faceType = TYPE_TRI;
break; break;
case TYPE_HEX: }
case TYPE_HEX: {
_metaEl0 = new MHexahedron(_metaVert); _metaEl0 = new MHexahedron(_metaVert);
metaNorm = SVector3(0.);
faceType = TYPE_QUA;
break; break;
default: }
default: {
Msg::Error("SuperEl not implemented for element of type %d", type); Msg::Error("SuperEl not implemented for element of type %d", type);
return; return;
}
} }
computeBaseNorm(metaNorm, baseVert, topPrimVert, _baseNorm);
// Add HO vertices in top face // Add HO vertices in top face
for (int iV = topPrimVert.size(); iV < topInd.size(); iV++) { for (int iV = topPrimVert.size(); iV < topInd.size(); iV++) {
...@@ -174,24 +310,10 @@ MetaEl::MetaEl(int type, int order, const std::vector<MVertex*> &baseVert, ...@@ -174,24 +310,10 @@ MetaEl::MetaEl(int type, int order, const std::vector<MVertex*> &baseVert,
_metaVert[ind] = new MVertex(p.x(), p.y(), p.z()); _metaVert[ind] = new MVertex(p.x(), p.y(), p.z());
} }
// Get incomplete high-order meta-element basis
const int tag = ElementType::getTag(type, order, true); // Get tag for incomplete element type
if (!tag) return;
const nodalBasis *incBasis = BasisFactory::getNodalBasis(tag);
// Add remaining face and interior points // Add remaining face and interior points
for (int iV = 0; iV < otherInd.size(); iV++) { for (int iV = 0; iV < otherInd.size(); iV++)
const int ind = otherInd[iV]; _metaVert[otherInd[iV]] = new MVertex(0., 0., 0.);
double shapeFunc[1256]; placeOtherNodes();
incBasis->f(points(ind, 0), points(ind, 1), points(ind, 2), shapeFunc);
double x = 0., y = 0., z = 0.;
for (int iSF = 0; iSF < incBasis->getNumShapeFunctions(); iSF++) {
x += shapeFunc[iSF] * _metaVert[iSF]->x();
y += shapeFunc[iSF] * _metaVert[iSF]->y();
z += shapeFunc[iSF] * _metaVert[iSF]->z();
}
_metaVert[ind] = new MVertex(x, y, z);
}
// Create high-order meta-element // Create high-order meta-element
switch (type) { switch (type) {
...@@ -205,37 +327,10 @@ MetaEl::MetaEl(int type, int order, const std::vector<MVertex*> &baseVert, ...@@ -205,37 +327,10 @@ MetaEl::MetaEl(int type, int order, const std::vector<MVertex*> &baseVert,
_metaEl = new MHexahedronN(_metaVert, order); _metaEl = new MHexahedronN(_metaVert, order);
break; break;
} }
// // Scale meta-element if not valid
// // TODO: Scale depending on target Jmin?
// // TODO: Could be improved by using complete meta-element and use optimization
// for (int iter = 0; iter < 10; iter++) {
// if (_metaEl->distoShapeMeasure() > 0) break;
// if (iter == 0)
// Msg::Warning("A meta-element has a negative distortion, trying to scale");
// for (int i = 0; i < topPrimVert.size(); ++i) { // Move top primary vert.
// MVertex *&vb = _metaVert[baseInd[i]], *&vt = _metaVert[topInd[i]];
// SPoint3 pb = vb->point(), pt = vt->point();
// pt += SPoint3(0.25*(pt-pb));
// vt->setXYZ(pt.x(), pt.y(), pt.z());
// }
// for (int i=topPrimVert.size(); i<topInd.size(); ++i) { // Recompute HO vert. in top face
// SPoint3 p;
// const int ind = topInd[i];
// _metaEl0->pnt(points(ind,0), points(ind,1), points(ind,2),p);
// _metaVert[ind]->setXYZ(p.x(), p.y(), p.z());
// }
// for (int i=0; i<otherInd.size(); ++i) { // Recompute vert. not in base & top faces
// SPoint3 p;
// const int ind = otherInd[i];
// _metaEl0->pnt(points(ind,0), points(ind,1), points(ind,2),p);
// _metaVert[ind]->setXYZ(p.x(), p.y(), p.z());
// }
// }
} }
MetaEl::~MetaEl() MetaEl::~MetaEl()
{ {
for (int i = 0; i < _metaVert.size(); i++) delete _metaVert[i]; for (int i = 0; i < _metaVert.size(); i++) delete _metaVert[i];
...@@ -246,27 +341,62 @@ MetaEl::~MetaEl() ...@@ -246,27 +341,62 @@ MetaEl::~MetaEl()
void MetaEl::curveTop(double factor) // Place free interior and face vertices equidistantly between top and base faces
void MetaEl::placeOtherNodes()
{
for (int iVO = 0; iVO < _mInfo.otherInd.size(); iVO++) {
const int &indF = _mInfo.otherInd[iVO];
const int &iVB = _mInfo.other2Base[iVO], indB = _mInfo.baseInd[iVB];
const int &iVT = _mInfo.other2Top[iVO], indT = _mInfo.topInd[iVT];
const int iUVWNorm = _mInfo.dim - 1;
const double fact = 0.5 * (_mInfo.points(indF, iUVWNorm)+1.);
SPoint3 pB = _metaVert[indB]->point(), pT = _metaVert[indT]->point();
const double newX = (1.-fact)*pB.x() + fact*pT.x();
const double newY = (1.-fact)*pB.y() + fact*pT.y();
const double newZ = (1.-fact)*pB.z() + fact*pT.z();
_metaVert[indF]->setXYZ(newX, newY, newZ);
}
}
void MetaEl::setCurvedTop(double factor)
{
// Place HO vertices of top face so that meta-elemnt has almost uniform thickness
for (int iVFT = 0; iVFT < _mInfo.freeTopInd.size(); iVFT++) {
const int &indFT = _mInfo.freeTopInd[iVFT];
const int &iVB = _mInfo.freeTop2Base[iVFT], indB = _mInfo.baseInd[iVB];
const int iUVWNorm = _mInfo.dim - 1;
const SPoint3 pB = _metaVert[indB]->point();
const SVector3 &norm = _baseNorm[iVB];
const double newX = pB.x() + factor*norm.x();
const double newY = pB.y() + factor*norm.y();
const double newZ = pB.z() + factor*norm.z();
_metaVert[indFT]->setXYZ(newX, newY, newZ);
}
// Place the other free nodes equidistantly between base and top faces
placeOtherNodes();
}
void MetaEl::setFlatTop()
{ {
// Get info on meta-element type // Get info on meta-element type
// const int &nbVert = _mInfo.nbVert;
const fullMatrix<double> &points = _mInfo.points; const fullMatrix<double> &points = _mInfo.points;
const std::vector<int> &baseInd = _mInfo.baseInd; const std::vector<int> &freeTopInd = _mInfo.topInd;
const std::vector<int> &topInd = _mInfo.topInd;
// const std::vector<int> &edgeInd = _mInfo.edgeInd;
// const std::vector<int> &otherInd = _mInfo.otherInd;
// Compute displacement of HO vertices in base face // Put top vertices on straight meta-element
for (int iV = 0; iV < baseInd.size(); iV++) { for (int iVFT = 0; iVFT < freeTopInd.size(); iVFT++) {
SPoint3 p0B, p0T; const int &indFT = freeTopInd[iVFT];
const int indB = baseInd[iV], indT = topInd[iV]; SPoint3 p;
_metaEl0->pnt(points(indB, 0), points(indB, 1), points(indB, 2), p0B); _metaEl0->pnt(points(indFT, 0), points(indFT, 1), points(indFT, 2), p);
_metaEl0->pnt(points(indT, 0), points(indT, 1), points(indT, 2), p0T); _metaVert[indFT]->setXYZ(p.x(), p.y(), p.z());
SPoint3 pB = _metaVert[indB]->point();
_metaVert[indT]->x() = p0T.x() + factor * (pB.x()-p0B.x());
_metaVert[indT]->y() = p0T.y() + factor * (pB.y()-p0B.y());
_metaVert[indT]->z() = p0T.z() + factor * (pB.z()-p0B.z());
} }
// Place the other free nodes equidistantly between base and top faces
placeOtherNodes();
} }
...@@ -279,6 +409,7 @@ bool MetaEl::isPointIn(const SPoint3 p) const ...@@ -279,6 +409,7 @@ bool MetaEl::isPointIn(const SPoint3 p) const
} }
bool MetaEl::straightToCurved(double *xyzS, double *xyzC) const bool MetaEl::straightToCurved(double *xyzS, double *xyzC) const
{ {
double uvw[3]; double uvw[3];
...@@ -295,6 +426,7 @@ bool MetaEl::straightToCurved(double *xyzS, double *xyzC) const ...@@ -295,6 +426,7 @@ bool MetaEl::straightToCurved(double *xyzS, double *xyzC) const
} }
std::string MetaEl::printPOS() std::string MetaEl::printPOS()
{ {
std::vector<MVertex*> verts; std::vector<MVertex*> verts;
......
contrib/HighOrderMeshOptimizer/MetaEl.h
+ 14
− 4
View file @ e339c0c2
...@@ -37,9 +37,11 @@ class MetaEl ...@@ -37,9 +37,11 @@ class MetaEl
{ {
public: public:
MetaEl(int type, int order, const std::vector<MVertex*> &baseVert, MetaEl(int type, int order, const std::vector<MVertex*> &baseVert,
const std::vector<MVertex*> &topPrimVert); const std::vector<MVertex*> &topPrimVert);
~MetaEl(); ~MetaEl();
void curveTop(double factor); void placeOtherNodes();
void setCurvedTop(double factor);
void setFlatTop();
bool isOK() const { return _metaEl; } bool isOK() const { return _metaEl; }
bool isPointIn(const SPoint3 p) const; bool isPointIn(const SPoint3 p) const;
bool straightToCurved(double *xyzS, double *xyzC) const; bool straightToCurved(double *xyzS, double *xyzC) const;
...@@ -57,9 +59,12 @@ public: ...@@ -57,9 +59,12 @@ public:
private: private:
struct metaInfoType { struct metaInfoType {
int nbVert; int dim, nbVert;
fullMatrix<double> points; fullMatrix<double> points;
std::vector<int> baseInd, topInd, edgeInd, otherInd; fullMatrix<double> baseLinShapeFuncVal, baseGradShapeFuncVal;
const nodalBasis *baseFaceBasis, *linBaseFaceBasis;
std::vector<int> baseInd, topInd, freeTopInd, edgeInd, otherInd;
std::vector<int> freeTop2Base, other2Base, other2Top;
metaInfoType(int type, int order); metaInfoType(int type, int order);
}; };
static std::map<int, metaInfoType> _metaInfo; static std::map<int, metaInfoType> _metaInfo;
...@@ -67,8 +72,13 @@ private: ...@@ -67,8 +72,13 @@ private:
const metaInfoType &_mInfo; const metaInfoType &_mInfo;
std::vector<MVertex*> _metaVert; std::vector<MVertex*> _metaVert;
MElement *_metaEl, *_metaEl0; MElement *_metaEl, *_metaEl0;
std::vector<SVector3> _baseNorm;
const metaInfoType &getMetaInfo(int elType, int order); const metaInfoType &getMetaInfo(int elType, int order);
void computeBaseNorm(const SVector3 &metaNorm,
const std::vector<MVertex*> &baseVert,
const std::vector<MVertex*> &topPrimVert,
std::vector<SVector3> &baseNorm);
}; };
#endif // _METAEL_H_ #endif // _METAEL_H_
contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp
+ 226
− 30
View file @ e339c0c2
...@@ -49,6 +49,7 @@ ...@@ -49,6 +49,7 @@
#include "BasisFactory.h" #include "BasisFactory.h"
#include "MetaEl.h" #include "MetaEl.h"
#include "qualityMeasuresJacobian.h" #include "qualityMeasuresJacobian.h"
#include "CADDistances.h"
...@@ -673,11 +674,197 @@ void curveElement(const MetaEl &metaElt, std::set<MVertex*> &movedVert, ...@@ -673,11 +674,197 @@ void curveElement(const MetaEl &metaElt, std::set<MVertex*> &movedVert,
double calcCADDistSq2D(GEntity *geomEnt, MElement *bndElt)
{
const int nbVert = bndElt->getNumVertices();
const int nbPrimVert = bndElt->getNumPrimaryVertices();
const GradientBasis *gb = BasisFactory::getGradientBasis(FuncSpaceData(bndElt));
// Coordinates of nodes
fullMatrix<double> nodesXYZ(nbVert, 3);
for (int iV = 0; iV < nbVert; iV++) {
MVertex *vert = bndElt->getVertex(iV);
nodesXYZ(iV, 0) = vert->x();
nodesXYZ(iV, 1) = vert->y();
nodesXYZ(iV, 2) = vert->z();
}
// Compute distance
const GEdge *ge = geomEnt->cast2Edge();
const Range<double> parBounds = ge->parBounds(0);
std::vector<SVector3> tanCAD(nbVert);
for (int iV = 0; iV < nbVert; iV++) {
MVertex *vert = bndElt->getVertex(iV);
double tCAD;
if (vert->onWhat() == geomEnt) // If HO vertex, ...
vert->getParameter(0, tCAD); // ... get stored param. coord. (can be only line).
else { // Otherwise, get param. coord. from CAD.
if (ge->getBeginVertex() &&
(ge->getBeginVertex()->mesh_vertices[0] == vert))
tCAD = parBounds.low();
else if (ge->getEndVertex() &&
(ge->getEndVertex()->mesh_vertices[0] == vert))
tCAD = parBounds.high();
else
tCAD = ge->parFromPoint(vert->point());
}
tanCAD[iV] = ge->firstDer(tCAD); // Compute tangent at vertex
tanCAD[iV].normalize(); // Normalize tangent
}
return taylorDistanceSq1D(gb, nodesXYZ, tanCAD);
}
void optimizeCADDist2DP2(GEntity *geomEnt, std::vector<MVertex*> &baseVert)
{
static const int NPTS = 1000;
MLine3 bndMetaElt(baseVert);
MVertex* &vert = baseVert[2];
const double xS = vert->x(), yS = vert->y();
const GEdge *ge = geomEnt->cast2Edge();
const Range<double> parBounds = ge->parBounds(0);
const double du = (parBounds.high()-parBounds.low())/double(NPTS-1);
double uMin = 0., distSqMin = 1e300;
double uDbg;
vert->getParameter(0, uDbg);
for (double u = parBounds.low(); u < parBounds.high(); u += du) {
vert->setParameter(0, u);
GPoint gp = ge->point(u);
vert->setXYZ(gp.x(), gp.y(), gp.z());
const double distSq = calcCADDistSq2D(geomEnt, &bndMetaElt);
if (distSq < distSqMin) { uMin = u; distSqMin = distSq; }
}
vert->setParameter(0, uMin);
GPoint gp = ge->point(uMin);
vert->setXYZ(gp.x(), gp.y(), gp.z());
}
// void calcCADDistSqAndGradients(int dim, GEntity *geomEnt, MElement *bndElt,
// double &dist, std::vector<double> &gradDist)
// {
// const int nbVert = bndElt->getNumVertices();
// const int nbPrimVert = bndElt->getNumPrimaryVertices();
// const GradientBasis *gb = BasisFactory::getGradientBasis(FuncSpaceData(bndElt));
// // Coordinates of nodes
// fullMatrix<double> nodesXYZ(nbVert, 3);
// for (int iV = 0; iV < nbVert; iV++) {
// MVertex *vert = bndElt->getVertex(iV);
// nodesXYZ(iV, 0) = vert->x();
// nodesXYZ(iV, 1) = vert->y();
// nodesXYZ(iV, 2) = vert->z();
// }
// // Compute distance and gradients
// if (dim == 2) { // 2D
// const GEdge *ge = geomEnt->cast2Edge();
// const Range<double> parBounds = ge->parBounds(0);
// const double eps = 1.e-6 * (parBounds.high()-parBounds.low());
// std::vector<SVector3> tanCAD(nbVert);
// for (int iV = 0; iV < nbVert; iV++) {
// MVertex *vert = bndElt->getVertex(iV);
// double tCAD;
// if (iV >= nbPrimVert) // If HO vertex, ...
// vert->getParameter(0, tCAD); // ... get stored param. coord. (can be only line).
// else { // Otherwise, get param. coord. from CAD.
// if (ge->getBeginVertex() &&
// (ge->getBeginVertex()->mesh_vertices[0] == vert))
// tCAD = parBounds.low();
// else if (ge->getEndVertex() &&
// (ge->getEndVertex()->mesh_vertices[0] == vert))
// tCAD = parBounds.high();
// else
// tCAD = ge->parFromPoint(vert->point());
// }
// tanCAD[iV] = ge->firstDer(tCAD); // Compute tangent at vertex
// tanCAD[iV].normalize(); // Normalize tangent
// }
// dist = taylorDistanceSq1D(gb, nodesXYZ, tanCAD);
// for (int iV = nbPrimVert; iV < nbVert; iV++) {
// const double xS = nodesXYZ(iV, 0), yS = nodesXYZ(iV, 1),
// zS = nodesXYZ(iV, 2); // Save coord. of perturbed node for FD
// const SVector3 tanCADS = tanCAD[iV]; // Save tangent to CAD at perturbed node
// double tCAD;
// bndElt->getVertex(iV)->getParameter(0, tCAD);
// tCAD += eps; // New param. coord. of perturbed node
// GPoint gp = ge->point(tCAD); // New coord. of perturbed node
// nodesXYZ(iV, 0) = gp.x();
// nodesXYZ(iV, 1) = gp.y();
// nodesXYZ(iV, 2) = gp.z();
// tanCAD[iV] = ge->firstDer(tCAD); // New tangent to CAD at perturbed node
// tanCAD[iV].normalize(); // Normalize new tangent
// const double sDistDiff = taylorDistanceSq1D(gb, nodesXYZ, tanCAD); // Compute distance with perturbed node
// gradDist[iV-nbPrimVert] = (sDistDiff-dist) / eps; // Compute gradient
// nodesXYZ(iV, 0) = xS; nodesXYZ(iV, 1) = yS; nodesXYZ(iV, 2) = zS; // Restore coord. of perturbed node
// tanCAD[iV] = tanCADS; // Restore tan. to CAD at perturbed node
// }
// }
// else { // 3D
// const GFace *gf = geomEnt->cast2Face();
// const Range<double> parBounds0 = gf->parBounds(0), parBounds1 = gf->parBounds(1);
// const double eps0 = 1.e-6 * (parBounds0.high()-parBounds0.low());
// const double eps1 = 1.e-6 * (parBounds1.high()-parBounds1.low());
// std::vector<SVector3> normCAD(nbVert);
// for (int iV = 0; iV < nbVert; iV++) {
// MVertex *vert = bndElt->getVertex(iV);
// SPoint2 pCAD;
// if (iV >= nbPrimVert) { // If HO vertex, get parameters
// vert->getParameter(0, pCAD[0]);
// vert->getParameter(1, pCAD[1]);
// }
// else
// reparamMeshVertexOnFace(vert, gf, pCAD); // If not free vertex, reparametrize on surface
// normCAD[iV] = gf->normal(pCAD); // Compute normal at vertex
// normCAD[iV].normalize(); // Normalize normal
// }
// dist = taylorDistanceSq2D(gb, nodesXYZ, normCAD);
// for (int iV = nbPrimVert; iV < nbVert; iV++) {
// MVertex *vert = bndElt->getVertex(iV);
// const double xS = nodesXYZ(iV, 0), yS = nodesXYZ(iV, 1),
// zS = nodesXYZ(iV, 2); // Save coord. of perturbed node for FD
// const SVector3 normCADS = normCAD[iV]; // Save normal to CAD at perturbed node
// SPoint2 pCAD0; // New param. coord. of perturbed node in 1st dir.
// vert->getParameter(0, pCAD0[0]);
// pCAD0 += eps0;
// vert->getParameter(1, pCAD0[1]);
// GPoint gp0 = gf->point(pCAD0); // New coord. of perturbed node in 1st dir.
// nodesXYZ(iV, 0) = gp0.x();
// nodesXYZ(iV, 1) = gp0.y();
// nodesXYZ(iV, 2) = gp0.z();
// normCAD[iV] = gf->normal(pCAD0); // New normal to CAD at perturbed node in 1st dir.
// normCAD[iV].normalize(); // Normalize new normal
// const double sDistDiff0 = taylorDistanceSq2D(gb, nodesXYZ, normCAD); // Compute distance with perturbed node in 1st dir.
// gradDist[2*(iV-nbPrimVert)] = (sDistDiff0-dist) / eps0; // Compute gradient in 1st dir.
// SPoint2 pCAD1; // New param. coord. of perturbed node in 2nd dir.
// vert->getParameter(0, pCAD1[0]);
// vert->getParameter(1, pCAD1[1]);
// pCAD1 += eps1;
// GPoint gp1 = gf->point(pCAD1); // New coord. of perturbed node in 2nd dir.
// nodesXYZ(iV, 0) = gp1.x();
// nodesXYZ(iV, 1) = gp1.y();
// nodesXYZ(iV, 2) = gp1.z();
// normCAD[iV] = gf->normal(pCAD1); // New normal to CAD at perturbed node in 2nd dir.
// normCAD[iV].normalize(); // Normalize new normal
// double sDistDiff1 = taylorDistanceSq2D(gb, nodesXYZ, normCAD); // Compute distance with perturbed node in 2nd dir.
// gradDist[2*(iV-nbPrimVert)+1] = (sDistDiff1-dist) / eps1; // Compute gradient in 2nd dir.
// nodesXYZ(iV, 0) = xS; // Restore coord. of perturbed node
// nodesXYZ(iV, 1) = yS;
// nodesXYZ(iV, 2) = zS;
// normCAD[iV] = normCADS; // Restore tan. to CAD at perturbed node
// }
// }
// }
double curveAndMeasureAboveEl(MetaEl &metaElt, MElement *lastElt, double curveAndMeasureAboveEl(MetaEl &metaElt, MElement *lastElt,
MElement *aboveElt, double deformFact) MElement *aboveElt, double deformFact)
{ {
double minJacDet, maxJacDet; double minJacDet, maxJacDet;
metaElt.curveTop(deformFact); metaElt.setCurvedTop(deformFact);
std::set<MVertex*> movedVertDum; std::set<MVertex*> movedVertDum;
curveElement(metaElt, movedVertDum, lastElt); curveElement(metaElt, movedVertDum, lastElt);
jacobianBasedQuality::minMaxJacobianDeterminant(aboveElt, minJacDet, maxJacDet); jacobianBasedQuality::minMaxJacobianDeterminant(aboveElt, minJacDet, maxJacDet);
...@@ -686,44 +873,53 @@ double curveAndMeasureAboveEl(MetaEl &metaElt, MElement *lastElt, ...@@ -686,44 +873,53 @@ double curveAndMeasureAboveEl(MetaEl &metaElt, MElement *lastElt,
void curveColumn(int metaElType, const std::vector<MVertex*> &baseVert, void curveColumn(GEntity *geomEnt, int metaElType,
std::vector<MVertex*> &baseVert,
const std::vector<MVertex*> &topPrimVert, MElement *aboveElt, const std::vector<MVertex*> &topPrimVert, MElement *aboveElt,
std::vector<MElement*> &blob, std::set<MVertex*> &movedVert, std::vector<MElement*> &blob, std::set<MVertex*> &movedVert,
DbgOutputMeta &dbgOut) DbgOutputMeta &dbgOut)
{ {
static const double MINQUAL = 0.1, TOL = 0.01, MAXITER = 10; static const double MINQUAL = 0.01, TOL = 0.01, MAXITER = 10;
// Order
const int order = blob[0]->getPolynomialOrder();
// If 2D P2, modify base vertices to minimize distance between wall edge and CAD
if ((metaElType == TYPE_QUA) &&
(blob[0]->getPolynomialOrder() == 2)) {
optimizeCADDist2DP2(geomEnt, baseVert);
}
// Create meta-element // Create meta-element
int order = blob[0]->getPolynomialOrder();
MetaEl metaElt(metaElType, order, baseVert, topPrimVert); MetaEl metaElt(metaElType, order, baseVert, topPrimVert);
// Curve top face of meta-element while avoiding breaking the element above // Curve top face of meta-element while avoiding breaking the element above
MElement* &lastElt = blob.back(); // MElement* &lastElt = blob.back();
// std::cout << "DBGTT: lastElt = " << lastElt->getNum() << ", aboveElt " << aboveElt->getNum() << " min. SJ: " << aboveElt->minScaledJacobian(); // double minJacDet, maxJacDet;
double minJacDet, maxJacDet; // double deformMin = 0., qualDeformMin = 1.;
double deformMin = 0., qualDeformMin = 1.; // double deformMax = 1.;
double deformMax = 1.; // double qualDeformMax = curveAndMeasureAboveEl(metaElt, lastElt, aboveElt,
double qualDeformMax = curveAndMeasureAboveEl(metaElt, lastElt, aboveElt, // deformMax);
deformMax); // if (qualDeformMax < MINQUAL) { // Max deformation makes element above invalid
if (qualDeformMax < MINQUAL) { // Max deformation makes element above invalid // for (int iter = 0; iter < MAXITER; iter++) { // Bisection to find max. deformation that makes element above valid
for (int iter = 0; iter < MAXITER; iter++) { // Bisection to find max. deformation that makes element above valid // const double deformMid = 0.5 * (deformMin + deformMax);
const double deformMid = 0.5 * (deformMin + deformMax); // const double qualDeformMid = curveAndMeasureAboveEl(metaElt, lastElt,
const double qualDeformMid = curveAndMeasureAboveEl(metaElt, lastElt, // aboveElt, deformMid);
aboveElt, deformMid); // if (std::abs(deformMax-deformMin) < TOL) break;
if (std::abs(deformMax-deformMin) < TOL) break; // const bool signDeformMax = (qualDeformMax < MINQUAL);
const bool signDeformMax = (qualDeformMax < MINQUAL); // const bool signDeformMid = (qualDeformMid < MINQUAL);
const bool signDeformMid = (qualDeformMid < MINQUAL); // if (signDeformMid == signDeformMax) deformMax = deformMid;
if (signDeformMid == signDeformMax) deformMax = deformMid; // else deformMin = deformMid;
else deformMin = deformMid; // }
} // metaElt.setFlatTop();
} // }
for (int iV = 0; iV < lastElt->getNumVertices(); iV++) // for (int iV = 0; iV < lastElt->getNumVertices(); iV++)
movedVert.insert(lastElt->getVertex(iV)); // movedVert.insert(lastElt->getVertex(iV));
// std::cout << " -> " << aboveElt->minScaledJacobian() << "\n";
dbgOut.addMetaEl(metaElt); dbgOut.addMetaEl(metaElt);
for (int iEl = 0; iEl < blob.size()-1; iEl++) for (int iEl = 0; iEl < blob.size(); iEl++)
// for (int iEl = 0; iEl < blob.size()-1; iEl++)
curveElement(metaElt, movedVert, blob[iEl]); curveElement(metaElt, movedVert, blob[iEl]);
} }
...@@ -754,7 +950,7 @@ void curveMeshFromBnd(MEdgeVecMEltMap &ed2el, MFaceVecMEltMap &face2el, ...@@ -754,7 +950,7 @@ void curveMeshFromBnd(MEdgeVecMEltMap &ed2el, MFaceVecMEltMap &face2el,
std::set<MVertex*> movedVert; std::set<MVertex*> movedVert;
for(std::list<MElement*>::iterator itBE = bndEl.begin(); for(std::list<MElement*>::iterator itBE = bndEl.begin();
itBE != bndEl.end(); itBE++) { // Loop over bnd. elements itBE != bndEl.end(); itBE++) { // Loop over bnd. elements
const int bndType = (*itBE)->getType(); const int bndType = (*itBE)->getType();
int metaElType; int metaElType;
bool foundCol; bool foundCol;
...@@ -790,8 +986,8 @@ void curveMeshFromBnd(MEdgeVecMEltMap &ed2el, MFaceVecMEltMap &face2el, ...@@ -790,8 +986,8 @@ void curveMeshFromBnd(MEdgeVecMEltMap &ed2el, MFaceVecMEltMap &face2el,
DbgOutputCol dbgOutCol; DbgOutputCol dbgOutCol;
dbgOutCol.addBlob(blob); dbgOutCol.addBlob(blob);
dbgOutCol.write("col_KO", (*itBE)->getNum()); dbgOutCol.write("col_KO", (*itBE)->getNum());
curveColumn(metaElType, baseVert, topPrimVert, aboveElt, blob, movedVert, curveColumn(bndEnt, metaElType, baseVert, topPrimVert, aboveElt, blob,
dbgOut); movedVert, dbgOut);
dbgOutCol.write("col_OK", (*itBE)->getNum()); dbgOutCol.write("col_OK", (*itBE)->getNum());
} }
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment