diff --git a/contrib/HighOrderMeshOptimizer/MetaEl.cpp b/contrib/HighOrderMeshOptimizer/MetaEl.cpp
index 17fd8b3e913532cf1d4193a39ef707933a17d1e6..a180ff5e05d2f2d2708ab83233a700c11fcb5ff4 100644
--- a/contrib/HighOrderMeshOptimizer/MetaEl.cpp
+++ b/contrib/HighOrderMeshOptimizer/MetaEl.cpp
@@ -42,70 +42,136 @@
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)
{
+ static const double TOL = 1e-10;
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) {
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;
case TYPE_PRI:
+ dim = 3;
iBaseFace = 0; iTopFace = 1; nbEdVert = 6 + 9*(order-1);
+ nbPrimTopVert = 3;
+ baseFaceType = TYPE_TRI;
break;
case TYPE_HEX:
+ dim = 3;
iBaseFace = 0; iTopFace = 5; nbEdVert = 8 + 12*(order-1);
+ nbPrimTopVert = 4;
+ baseFaceType = TYPE_QUA;
break;
default:
Msg::Error("MetaEl not implemented for element of type %d", type);
+ dim = 0;
nbVert = 0;
return;
}
- // Get HO nodal basis
-// const int tag = ElementType::getTag(type, order, true); // Get tag for incomplete element type
+ // Get HO nodal base
const int tag = ElementType::getTag(type, order); // Get tag for complete element type
if (!tag) return;
const nodalBasis *basis = BasisFactory::getNodalBasis(tag);
nbVert = basis->getNumShapeFunctions();
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));
- topInd = basis->getClosure(basis->getClosureId(iTopFace, 0, 0));
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());
+ 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++)
if (foundVerts.find(iV) == foundVerts.end()) {
edgeInd.push_back(iV);
foundVerts.insert(iV);
}
- // Remaining face and interior vertices
- for(int iV = 0; iV < nbVert; iV++)
+ // Indices of remaining face and interior vertices
+ for(int iV = 0; iV < nbVert; iV++) {
if (foundVerts.find(iV) == foundVerts.end()) {
otherInd.push_back(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)
+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,
const std::vector<MVertex*> &topPrimVert) :
_mInfo(getMetaInfo(type, order)), _metaEl(0), _metaEl0(0)
@@ -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)
_metaVert.resize(nbVert);
- for (int iV = 0; iV < baseInd.size(); iV++)
- _metaVert[baseInd[iV]] = new MVertex(*baseVert[iV]);
+ for (int iV = 0; iV < baseInd.size(); 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++)
_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) {
- case TYPE_QUA:
+ case TYPE_QUA: {
_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;
- case TYPE_PRI:
+ }
+ case TYPE_PRI: {
_metaEl0 = new MPrism(_metaVert);
+ metaNorm = SVector3(0.);
+ faceType = TYPE_TRI;
break;
- case TYPE_HEX:
+ }
+ case TYPE_HEX: {
_metaEl0 = new MHexahedron(_metaVert);
+ metaNorm = SVector3(0.);
+ faceType = TYPE_QUA;
break;
- default:
+ }
+ default: {
Msg::Error("SuperEl not implemented for element of type %d", type);
return;
+ }
}
+ computeBaseNorm(metaNorm, baseVert, topPrimVert, _baseNorm);
// Add HO vertices in top face
for (int iV = topPrimVert.size(); iV < topInd.size(); iV++) {
@@ -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());
}
- // 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
- for (int iV = 0; iV < otherInd.size(); iV++) {
- const int ind = otherInd[iV];
- double shapeFunc[1256];
- 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);
- }
+ for (int iV = 0; iV < otherInd.size(); iV++)
+ _metaVert[otherInd[iV]] = new MVertex(0., 0., 0.);
+ placeOtherNodes();
// Create high-order meta-element
switch (type) {
@@ -205,37 +327,10 @@ MetaEl::MetaEl(int type, int order, const std::vector<MVertex*> &baseVert,
_metaEl = new MHexahedronN(_metaVert, order);
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()
{
for (int i = 0; i < _metaVert.size(); i++) delete _metaVert[i];
@@ -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
- // const int &nbVert = _mInfo.nbVert;
const fullMatrix<double> &points = _mInfo.points;
- const std::vector<int> &baseInd = _mInfo.baseInd;
- const std::vector<int> &topInd = _mInfo.topInd;
- // const std::vector<int> &edgeInd = _mInfo.edgeInd;
- // const std::vector<int> &otherInd = _mInfo.otherInd;
+ const std::vector<int> &freeTopInd = _mInfo.topInd;
- // Compute displacement of HO vertices in base face
- for (int iV = 0; iV < baseInd.size(); iV++) {
- SPoint3 p0B, p0T;
- const int indB = baseInd[iV], indT = topInd[iV];
- _metaEl0->pnt(points(indB, 0), points(indB, 1), points(indB, 2), p0B);
- _metaEl0->pnt(points(indT, 0), points(indT, 1), points(indT, 2), p0T);
- 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());
+ // Put top vertices on straight meta-element
+ for (int iVFT = 0; iVFT < freeTopInd.size(); iVFT++) {
+ const int &indFT = freeTopInd[iVFT];
+ SPoint3 p;
+ _metaEl0->pnt(points(indFT, 0), points(indFT, 1), points(indFT, 2), p);
+ _metaVert[indFT]->setXYZ(p.x(), p.y(), p.z());
}
+
+ // Place the other free nodes equidistantly between base and top faces
+ placeOtherNodes();
}
@@ -279,6 +409,7 @@ bool MetaEl::isPointIn(const SPoint3 p) const
}
+
bool MetaEl::straightToCurved(double *xyzS, double *xyzC) const
{
double uvw[3];
@@ -295,6 +426,7 @@ bool MetaEl::straightToCurved(double *xyzS, double *xyzC) const
}
+
std::string MetaEl::printPOS()
{
std::vector<MVertex*> verts;
diff --git a/contrib/HighOrderMeshOptimizer/MetaEl.h b/contrib/HighOrderMeshOptimizer/MetaEl.h
index e39a2d61c7e6d28eb9235aa631c4a48406a9f556..b4ddd08fc4e1791ac3b36060a4be94625c9d9041 100644
--- a/contrib/HighOrderMeshOptimizer/MetaEl.h
+++ b/contrib/HighOrderMeshOptimizer/MetaEl.h
@@ -37,9 +37,11 @@ class MetaEl
{
public:
MetaEl(int type, int order, const std::vector<MVertex*> &baseVert,
- const std::vector<MVertex*> &topPrimVert);
+ const std::vector<MVertex*> &topPrimVert);
~MetaEl();
- void curveTop(double factor);
+ void placeOtherNodes();
+ void setCurvedTop(double factor);
+ void setFlatTop();
bool isOK() const { return _metaEl; }
bool isPointIn(const SPoint3 p) const;
bool straightToCurved(double *xyzS, double *xyzC) const;
@@ -57,9 +59,12 @@ public:
private:
struct metaInfoType {
- int nbVert;
+ int dim, nbVert;
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);
};
static std::map<int, metaInfoType> _metaInfo;
@@ -67,8 +72,13 @@ private:
const metaInfoType &_mInfo;
std::vector<MVertex*> _metaVert;
MElement *_metaEl, *_metaEl0;
+ std::vector<SVector3> _baseNorm;
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_
diff --git a/contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp b/contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp
index 536c15fc8435b15e3f1fc18d086e93886db4e5f6..76da75f9c8faf5cf52f7e534888bb187964a165d 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomFastCurving.cpp
@@ -49,6 +49,7 @@
#include "BasisFactory.h"
#include "MetaEl.h"
#include "qualityMeasuresJacobian.h"
+#include "CADDistances.h"
@@ -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,
MElement *aboveElt, double deformFact)
{
double minJacDet, maxJacDet;
- metaElt.curveTop(deformFact);
+ metaElt.setCurvedTop(deformFact);
std::set<MVertex*> movedVertDum;
curveElement(metaElt, movedVertDum, lastElt);
jacobianBasedQuality::minMaxJacobianDeterminant(aboveElt, minJacDet, maxJacDet);
@@ -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,
std::vector<MElement*> &blob, std::set<MVertex*> &movedVert,
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
- int order = blob[0]->getPolynomialOrder();
MetaEl metaElt(metaElType, order, baseVert, topPrimVert);
// Curve top face of meta-element while avoiding breaking the element above
- MElement* &lastElt = blob.back();
- // std::cout << "DBGTT: lastElt = " << lastElt->getNum() << ", aboveElt " << aboveElt->getNum() << " min. SJ: " << aboveElt->minScaledJacobian();
- double minJacDet, maxJacDet;
- double deformMin = 0., qualDeformMin = 1.;
- double deformMax = 1.;
- double qualDeformMax = curveAndMeasureAboveEl(metaElt, lastElt, aboveElt,
- deformMax);
- 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
- const double deformMid = 0.5 * (deformMin + deformMax);
- const double qualDeformMid = curveAndMeasureAboveEl(metaElt, lastElt,
- aboveElt, deformMid);
- if (std::abs(deformMax-deformMin) < TOL) break;
- const bool signDeformMax = (qualDeformMax < MINQUAL);
- const bool signDeformMid = (qualDeformMid < MINQUAL);
- if (signDeformMid == signDeformMax) deformMax = deformMid;
- else deformMin = deformMid;
- }
- }
- for (int iV = 0; iV < lastElt->getNumVertices(); iV++)
- movedVert.insert(lastElt->getVertex(iV));
- // std::cout << " -> " << aboveElt->minScaledJacobian() << "\n";
+ // MElement* &lastElt = blob.back();
+ // double minJacDet, maxJacDet;
+ // double deformMin = 0., qualDeformMin = 1.;
+ // double deformMax = 1.;
+ // double qualDeformMax = curveAndMeasureAboveEl(metaElt, lastElt, aboveElt,
+ // deformMax);
+ // 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
+ // const double deformMid = 0.5 * (deformMin + deformMax);
+ // const double qualDeformMid = curveAndMeasureAboveEl(metaElt, lastElt,
+ // aboveElt, deformMid);
+ // if (std::abs(deformMax-deformMin) < TOL) break;
+ // const bool signDeformMax = (qualDeformMax < MINQUAL);
+ // const bool signDeformMid = (qualDeformMid < MINQUAL);
+ // if (signDeformMid == signDeformMax) deformMax = deformMid;
+ // else deformMin = deformMid;
+ // }
+ // metaElt.setFlatTop();
+ // }
+ // for (int iV = 0; iV < lastElt->getNumVertices(); iV++)
+ // movedVert.insert(lastElt->getVertex(iV));
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]);
}
@@ -754,7 +950,7 @@ void curveMeshFromBnd(MEdgeVecMEltMap &ed2el, MFaceVecMEltMap &face2el,
std::set<MVertex*> movedVert;
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();
int metaElType;
bool foundCol;
@@ -790,8 +986,8 @@ void curveMeshFromBnd(MEdgeVecMEltMap &ed2el, MFaceVecMEltMap &face2el,
DbgOutputCol dbgOutCol;
dbgOutCol.addBlob(blob);
dbgOutCol.write("col_KO", (*itBE)->getNum());
- curveColumn(metaElType, baseVert, topPrimVert, aboveElt, blob, movedVert,
- dbgOut);
+ curveColumn(bndEnt, metaElType, baseVert, topPrimVert, aboveElt, blob,
+ movedVert, dbgOut);
dbgOutCol.write("col_OK", (*itBE)->getNum());
}