diff --git a/CMakeLists.txt b/CMakeLists.txt
index bc64947300f71bbb2cd343b88bb29c05cd16729b..2e81dd63dd0af83d932bf65f26e4189bbfa1e4e2 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -149,7 +149,7 @@ set(GMSH_API
contrib/HighOrderMeshOptimizer/OptHomRun.h contrib/HighOrderMeshOptimizer/ParamCoord.h
contrib/HighOrderMeshOptimizer/OptHomFastCurving.h contrib/HighOrderMeshOptimizer/SuperEl.h
contrib/HighOrderMeshOptimizer/OptHomIntegralBoundaryDist.h
- contrib/HighOrderMeshOptimizer/OptHomCADDist.h
+ contrib/HighOrderMeshOptimizer/CADDistances.h
contrib/HighOrderMeshOptimizer/OptHomObjContribScaledJac.h
contrib/HighOrderMeshOptimizer/OptHomObjContribMetricMin.h
contrib/HighOrderMeshOptimizer/OptHomObjContribCADDist.h
diff --git a/contrib/HighOrderMeshOptimizer/CADDistances.cpp b/contrib/HighOrderMeshOptimizer/CADDistances.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..bf4b1066057456375e698c73e579e34c3e193697
--- /dev/null
+++ b/contrib/HighOrderMeshOptimizer/CADDistances.cpp
@@ -0,0 +1,588 @@
+// TODO: Header
+
+#include "SPoint3.h"
+#include "SVector3.h"
+#include "nodalBasis.h"
+#include "GEdge.h"
+#include "GFace.h"
+#include "MElement.h"
+#include "MLine.h"
+#include "MTriangle.h"
+#include "BasisFactory.h"
+#include "GModel.h"
+#if defined(HAVE_ANN)
+#include "ANN/ANN.h"
+#endif
+#include "CADDistances.h"
+
+
+namespace {
+
+
+// ---------------------------------- Discretization of curves ----------------------------------
+
+
+class parametricLine
+{
+ void recurDiscretize(const double &t1, const double &t2, const SPoint3 &p1,
+ const SPoint3 &p2, std::vector<SPoint3> &dpts,
+ std::vector<double> &ts, double Prec, int minDepth,
+ int maxDepth, int depth) const;
+public:
+ virtual ~parametricLine() {}
+ virtual SPoint3 operator()(double xi) const = 0;
+ virtual SVector3 derivative(double xi) const = 0;
+ virtual SVector3 secondDerivative(double xi) const = 0;
+ SVector3 curvature(double xi) const;
+ void discretize(std::vector<SPoint3> &dpts,std::vector<double> &ts,
+ double Prec, int minDepth, int maxDepth,
+ double t0 = 0., double t1 = 1.) const;
+};
+
+
+class parametricLineNodalBasis : public parametricLine
+{
+ const nodalBasis &_basis;
+ const std::vector<SPoint3> &_xyz;
+
+public:
+ parametricLineNodalBasis(const nodalBasis &basis, const std::vector<SPoint3> &xyz);
+ virtual SPoint3 operator()(double xi) const;
+ virtual SVector3 derivative(double xi) const;
+ virtual SVector3 secondDerivative(double xi) const;
+};
+
+
+class parametricLineGEdge : public parametricLine
+{
+ const GEdge *_edge;
+ double _t0, _t1;
+
+public:
+ parametricLineGEdge(const GEdge *edge, double t0, double t1);
+ virtual SPoint3 operator()(double xi) const;
+ virtual SVector3 derivative(double xi) const;
+ virtual SVector3 secondDerivative(double xi) const;
+};
+
+
+void parametricLine::recurDiscretize(const double &t1, const double &t2, const SPoint3 &p1,
+ const SPoint3 &p2, std::vector<SPoint3> &dpts,
+ std::vector<double> &ts, double Prec, int minDepth,
+ int maxDepth, int depth) const
+{
+ double t = 0.5 * (t2 + t1);
+ SPoint3 p = (*this)(t);
+ SVector3 dx(p, (p1+p2)*0.5);
+ if (((depth >= minDepth) && (dx.norm() < Prec)) || (depth >= maxDepth)) {
+ dpts.push_back(p);
+ ts.push_back(t);
+ dpts.push_back(p2);
+ ts.push_back(t2);
+ }
+ else {
+ recurDiscretize(t1, t, p1, p, dpts, ts, Prec, minDepth, maxDepth, depth+1);
+ recurDiscretize(t, t2, p, p2, dpts, ts, Prec, minDepth, maxDepth, depth+1);
+ }
+}
+
+
+void parametricLine::discretize(std::vector<SPoint3> &dpts,
+ std::vector<double> &ts,
+ double Prec, int minDepth, int maxDepth, double t0, double t1) const
+{
+ dpts.push_back((*this)(t0));
+ ts.push_back(t0);
+ recurDiscretize(t0, t1, dpts[0], (*this)(t1), dpts, ts, Prec, minDepth, maxDepth, 0);
+}
+
+
+SVector3 parametricLine::curvature(double xi) const
+{
+ SVector3 xp = derivative(xi);
+ SVector3 xpp = secondDerivative(xi);
+ const double nxp = xp.norm(), onxp = 1./nxp;
+ SVector3 c = (onxp*onxp*onxp)*(xpp*nxp-xp*dot(xp,xpp)*onxp);
+ return c;
+}
+
+
+parametricLineNodalBasis::parametricLineNodalBasis(const nodalBasis &basis,
+ const std::vector<SPoint3> &xyz):
+ _basis(basis), _xyz(xyz)
+{
+}
+
+
+SPoint3 parametricLineNodalBasis::operator()(double xi) const
+{
+ double psi[_xyz.size()];
+ SPoint3 p(0., 0., 0.);
+ _basis.f(-1 + 2 * xi, 0., 0., psi);
+ for (size_t j = 0; j < _xyz.size(); ++j) {
+ p[0] += psi[j] * _xyz[j].x();
+ p[1] += psi[j] * _xyz[j].y();
+ p[2] += psi[j] * _xyz[j].z();
+ }
+ return p;
+}
+
+
+SVector3 parametricLineNodalBasis::derivative(double xi) const
+{
+ double dpsi[_xyz.size()][3];
+ SVector3 p(0., 0., 0.);
+ _basis.df(-1 + 2 * xi, 0., 0., dpsi);
+ for (size_t j = 0; j < _xyz.size(); ++j) {
+ p[0] += dpsi[j][0] * _xyz[j].x();
+ p[1] += dpsi[j][0] * _xyz[j].y();
+ p[2] += dpsi[j][0] * _xyz[j].z();
+ }
+ return p;
+}
+
+
+SVector3 parametricLineNodalBasis::secondDerivative(double xi) const
+{
+ double ddpsi[_xyz.size()][3][3];
+ SVector3 p(0, 0, 0);
+ _basis.ddf(-1 + 2 * xi, 0, 0, ddpsi);
+ for (size_t j = 0; j < _xyz.size(); ++j) {
+ p[0] += ddpsi[j][0][0] * _xyz[j].x();
+ p[1] += ddpsi[j][0][0] * _xyz[j].y();
+ p[2] += ddpsi[j][0][0] * _xyz[j].z();
+ }
+ return p;
+}
+
+
+parametricLineGEdge::parametricLineGEdge(const GEdge *edge, double t0, double t1):
+ _edge(edge), _t0(t0), _t1(t1)
+{
+}
+
+
+SPoint3 parametricLineGEdge::operator()(double xi) const
+{
+ GPoint gp = _edge->point(_t0 + (_t1 - _t0) * xi);
+ return SPoint3 (gp.x(), gp.y(), gp.z());
+}
+
+
+SVector3 parametricLineGEdge::derivative(double xi) const
+{
+ return _edge->firstDer(_t0 + (_t1 - _t0) * xi);
+}
+
+
+SVector3 parametricLineGEdge::secondDerivative(double xi) const
+{
+ return _edge->secondDer(_t0 + (_t1 - _t0) * xi);
+}
+
+
+void oversample(std::vector<SPoint3> &s, double tol)
+{
+ std::vector<SPoint3> t;
+ for (unsigned int i=1;i<s.size();i++){
+ SPoint3 p0 = s[i-1];
+ SPoint3 p1 = s[i];
+ double d = p0.distance(p1);
+ int N = (int) (d / tol);
+ t.push_back(p0);
+ for (int j=1;j<N;j++){
+ const double xi = (double) j/ N;
+ t.push_back(p0 + (p1-p0)*xi);
+ }
+ }
+ t.push_back(s[s.size()-1]);
+ s = t;
+}
+
+
+// ---------------------------------- Discrete Fréchet distance ----------------------------------
+
+
+double recurMinMax(int i, int j, fullMatrix<double> &CA,
+ const std::vector<SPoint3> &P, const std::vector<SPoint3> &Q)
+{
+ double CAij;
+ if (CA(i,j) > -1.) {
+ CAij = CA(i,j);
+ }
+ else if (i==0 && j==0) {
+ CA(i,j) = P[0].distance(Q[0]); // update the CA permanent
+ CAij = CA(i,j); // set the current relevant value
+ }
+ else if (i>0 && j==0) {
+ CA(i,j) = std::max(recurMinMax(i-1, 0, CA, P, Q), P[i].distance(Q[1]));
+ CAij = CA(i,j);
+ }
+ else if (i==0 && j>0) {
+ CA(i,j) = std::max(recurMinMax(0, j-1, CA, P, Q), P[0].distance(Q[j]));
+ CAij = CA(i,j);
+ }
+ else if (i>0 && j>0) {
+ double temp = std::min(std::min(recurMinMax(i-1, j, CA, P, Q),
+ recurMinMax(i-1, j-1, CA, P, Q)), recurMinMax(i, j-1, CA, P, Q));
+ CA(i,j) = std::max(temp, P[i].distance(Q[j]));
+ CAij = CA(i,j);
+ }
+ else {
+ CA(i,j) = 1.e22;
+ CAij = CA(i,j);
+ }
+ return CAij;
+}
+
+
+double discreteFrechetDistanceBrute(const std::vector<SPoint3> &P,
+ const std::vector<SPoint3> &Q)
+{
+ const int sP = P.size();
+ const int sQ = Q.size();
+ fullMatrix<double> CA(sP, sQ);
+ CA.setAll(-1.);
+ double cm = recurMinMax(sP-1, sQ-1, CA, P, Q);
+ return cm;
+}
+
+
+// ---------------------------------- Discrete Hausdorff distance ----------------------------------
+
+
+double oneSidedMaxDist(const std::vector<SPoint3> &dpts1, const std::vector<SPoint3> &dpts2,
+ int &iMaxDist, int &jMaxDist)
+{
+ double maxDist = 0.;
+ for (unsigned int i=0; i<dpts1.size(); i++) {
+ double dist = 1.e22;
+ int jDist = 0;
+ for (unsigned int j=0; j<dpts2.size(); j++) {
+ double h = dpts1[i].distance(dpts2[j]);
+ if (h < dist){
+ dist = h;
+ jDist = j;
+ }
+ }
+ if (dist > maxDist) {
+ maxDist = dist;
+ jMaxDist = jDist;
+ iMaxDist = i;
+ }
+ }
+ return maxDist;
+}
+
+
+#if defined(HAVE_ANN)
+double oneSidedMaxDistFast(const std::vector<SPoint3> &dptsA, const ANNpointArray &ptsB,
+ ANNkd_tree *treeB, double tol)
+{
+ ANNidx idx[1];
+ ANNdist distSq[1];
+ ANNdist maxDistSq = 0.;
+ for (unsigned int i=0; i<dptsA.size(); i++) {
+ ANNcoord xyz[3] = {dptsA[i].x(), dptsA[i].y(), dptsA[i].z()};
+ treeB->annkSearch(xyz, 1, idx, distSq);
+ if (distSq[0] > maxDistSq) maxDistSq = distSq[0];
+ }
+ return sqrt(maxDistSq);
+}
+#endif
+
+
+double discreteHausdorffDistanceBrute(const std::vector<SPoint3> &dpts1,
+ const std::vector<SPoint3> &dpts2)
+{
+ int i1, j1, i2, j2;
+ double d1 = oneSidedMaxDist(dpts1, dpts2, i1, j1);
+ double d2 = oneSidedMaxDist(dpts2, dpts1, i2, j2);
+ return (d1 > d2) ? d1 : d2;
+}
+
+
+double discreteHausdorffDistanceFast(const std::vector<SPoint3> &dpts1,
+ const std::vector<SPoint3> &dpts2, double tol)
+{
+#if defined(HAVE_ANN)
+ ANNpointArray pts1 = annAllocPts(dpts1.size(), 3);
+ for (unsigned int k=0; k<dpts1.size(); k++) {
+ pts1[k][0] = dpts1[k].x();
+ pts1[k][1] = dpts1[k].y();
+ pts1[k][2] = dpts1[k].z();
+ }
+ ANNkd_tree *tree1 = new ANNkd_tree(pts1, dpts1.size(), 3);
+
+ ANNpointArray pts2 = annAllocPts(dpts2.size(), 3);
+ for (unsigned int k=0; k<dpts2.size(); k++) {
+ pts2[k][0] = dpts2[k].x();
+ pts2[k][1] = dpts2[k].y();
+ pts2[k][2] = dpts2[k].z();
+ }
+ ANNkd_tree *tree2 = new ANNkd_tree(pts2, dpts2.size(), 3);
+
+ double d1 = oneSidedMaxDistFast(dpts1, pts2, tree2, tol);
+ double d2 = oneSidedMaxDistFast(dpts2, pts1, tree1, tol);
+
+ delete tree1, tree2;
+ annDeallocPts(pts1);
+ annDeallocPts(pts2);
+
+ return (d1 > d2) ? d1 : d2;
+#else
+ Msg::Fatal("Gmsh should be compiled using ANN");
+#endif
+}
+
+
+template<int distDef>
+double discreteDistance(MLine *l, GEdge *ed, double tol, int meshDiscr, int geomDiscr)
+{
+ if (ed->geomType() == GEntity::Line) return 0.;
+
+ std::vector<SPoint3> dpts1, dpts2;
+ std::vector<double> ts1, ts2;
+
+ if (geomDiscr == CADDIST_DECASTELJAU)
+ ed->discretize(tol, dpts1, ts1);
+ else {
+ double t0, t1;
+ reparamMeshVertexOnEdge(l->getVertex(0), ed, t0);
+ reparamMeshVertexOnEdge(l->getVertex(1), ed, t1);
+ parametricLineGEdge l1(ed, t0, t1);
+ l1.discretize(dpts1, ts1, tol, 5, 45);
+ }
+ oversample(dpts1, tol);
+
+ if (meshDiscr == CADDIST_DECASTELJAU)
+ l->discretize(tol, dpts2, ts2);
+ else {
+ const nodalBasis &basis = *(l->getFunctionSpace());
+ std::vector<SPoint3> xyz;
+ const int nV = l->getNumVertices();
+ for (int i=0; i<nV; ++i) xyz.push_back(l->getVertex(i)->point());
+ parametricLineNodalBasis l2(basis, xyz);
+ int minDepth = std::ceil(std::log(nV)/std::log(2));
+ l2.discretize(dpts2, ts2, tol, minDepth, 10*minDepth);
+ }
+ oversample(dpts2, tol);
+// std::cout << "DBGTT: " << dpts1.size() << " points on model, " << dpts2.size() << " points on mesh\n";
+
+ switch (distDef) {
+ case 1: return discreteHausdorffDistanceBrute(dpts1, dpts2);
+ case 2: return discreteHausdorffDistanceFast(dpts1, dpts2, tol);
+ case 3: return discreteFrechetDistanceBrute(dpts1, dpts2);
+ }
+ return -1.;
+}
+
+
+}
+
+
+double discreteHausdorffDistanceBruteEdge(MLine *l, GEdge *ed, double tol,
+ int meshDiscr, int geomDiscr)
+{
+ return discreteDistance<1>(l, ed, tol, meshDiscr, geomDiscr);
+}
+
+
+double discreteHausdorffDistanceFastEdge(MLine *l, GEdge *ed, double tol,
+ int meshDiscr, int geomDiscr)
+{
+ return discreteDistance<2>(l, ed, tol, meshDiscr, geomDiscr);
+}
+
+
+double discreteFrechetDistanceEdge(MLine *l, GEdge *ed, double tol,
+ int meshDiscr, int geomDiscr)
+{
+ return discreteDistance<3>(l, ed, tol, meshDiscr, geomDiscr);
+}
+
+
+// ---------------------------------- Taylor-based distance ----------------------------------
+
+
+double taylorDistanceSq1D(const GradientBasis *gb, const fullMatrix<double> &nodesXYZ,
+ const std::vector<SVector3> &tanCAD)
+{
+ const int nV = nodesXYZ.size1();
+ fullMatrix<double> dxyzdX(nV, 3);
+ gb->getGradientsFromNodes(nodesXYZ, &dxyzdX, 0, 0);
+// const double dx = nodesXYZ(1, 0) - nodesXYZ(0, 0), dy = nodesXYZ(1, 1) - nodesXYZ(0, 1),
+// dz = nodesXYZ(1, 2) - nodesXYZ(0, 2), h = 0.5*sqrt(dx*dx+dy*dy+dz*dz)/double(nV-1);
+ double distSq = 0.;
+ for (int i=0; i<nV; i++) {
+ SVector3 tanMesh(dxyzdX(i, 0), dxyzdX(i, 1), dxyzdX(i, 2));
+ const double h = 0.25*tanMesh.normalize(); // Half of "local edge length"
+ SVector3 diff = (dot(tanCAD[i], tanMesh) > 0) ?
+ tanCAD[i] - tanMesh : tanCAD[i] + tanMesh;
+ distSq += h*h*diff.normSq();
+ }
+ return distSq;
+}
+
+
+double taylorDistanceSq2D(const GradientBasis *gb, const fullMatrix<double> &nodesXYZ,
+ const std::vector<SVector3> &normCAD)
+{
+ const int nV = nodesXYZ.size1();
+ fullMatrix<double> dxyzdX(nV, 3),dxyzdY(nV, 3);
+ gb->getGradientsFromNodes(nodesXYZ, &dxyzdX, &dxyzdY, 0);
+ double distSq = 0.;
+ for (int i=0; i<nV; i++) {
+ const double nz = dxyzdX(i, 0) * dxyzdY(i, 1) - dxyzdX(i, 1) * dxyzdY(i, 0);
+ const double ny = -dxyzdX(i, 0) * dxyzdY(i, 2) + dxyzdX(i, 2) * dxyzdY(i, 0);
+ const double nx = dxyzdX(i, 1) * dxyzdY(i, 2) - dxyzdX(i, 2) * dxyzdY(i, 1);
+ SVector3 normMesh(nx, ny, nz);
+ const double h = 0.25*sqrt(normMesh.normalize()); // Half of sqrt of "local area", to be adjusted w.r.t. el. type?
+ SVector3 diff = (dot(normCAD[i], normMesh) > 0) ?
+ normCAD[i] - normMesh : normCAD[i] + normMesh;
+ distSq += h*h*diff.normSq();
+ }
+ return distSq;
+}
+
+
+double taylorDistanceEdge(MLine *l, GEdge *ge)
+{
+ const int nV = l->getNumVertices();
+ const GradientBasis *gb = BasisFactory::getGradientBasis(FuncSpaceData(l));
+
+ // Coordinates of vertices
+ fullMatrix<double> nodesXYZ(nV, 3);
+ l->getNodesCoord(nodesXYZ);
+
+ // Tangent to CAD at vertices
+ std::vector<SVector3> tanCAD(nV);
+ for (int i=0; i<nV; i++) {
+ double tCAD;
+ reparamMeshVertexOnEdge(l->getVertex(i), ge, tCAD);
+ tanCAD[i] = ge->firstDer(tCAD);
+ tanCAD[i].normalize();
+ }
+
+ // Compute distance
+ return sqrt(taylorDistanceSq1D(gb, nodesXYZ, tanCAD));
+}
+
+
+double taylorDistanceFace(MElement *el, GFace *gf)
+{
+ const int nV = el->getNumVertices();
+ const GradientBasis *gb = BasisFactory::getGradientBasis(FuncSpaceData(el));
+
+ // Coordinates of vertices
+ fullMatrix<double> nodesXYZ(nV, 3);
+ el->getNodesCoord(nodesXYZ);
+
+ // Normal to CAD at vertices
+ std::vector<SVector3> normCAD(nV);
+ for (int i=0; i<nV; i++) {
+ SPoint2 pCAD;
+ reparamMeshVertexOnFace(el->getVertex(i), gf, pCAD);
+ normCAD[i] = gf->normal(pCAD);
+ normCAD[i].normalize();
+ }
+
+ // Compute distance
+ return sqrt(taylorDistanceSq2D(gb, nodesXYZ, normCAD));
+}
+
+
+// ---------------------------------- General functions ----------------------------------
+
+
+void distanceFromElementsToGeometry(GModel *gm, int dim, std::map<MElement*,double> &distances)
+{
+ std::map<MEdge,double,Less_Edge> dist2Edge;
+ for (GModel::eiter it = gm->firstEdge(); it != gm->lastEdge(); ++it) {
+ if ((*it)->geomType() == GEntity::Line) continue;
+ for (unsigned int i = 0; i < (*it)->lines.size(); i++) {
+ double d = taylorDistanceEdge((*it)->lines[i], *it);
+ MEdge e = (*it)->lines[i]->getEdge(0);
+ dist2Edge[e] = d;
+ }
+ }
+
+ std::map<MFace,double,Less_Face> dist2Face;
+ for(GModel::fiter it = gm->firstFace(); it != gm->lastFace(); ++it){
+ if ((*it)->geomType() == GEntity::Plane)continue;
+ for (unsigned int i = 0; i < (*it)->triangles.size(); i++){
+ double d = taylorDistanceFace((*it)->triangles[i], *it);
+ MFace f = (*it)->triangles[i]->getFace(0);
+ dist2Face[f] = d;
+ }
+ }
+
+ std::vector<GEntity*> entities;
+ gm->getEntities(entities);
+ for (int iEnt = 0; iEnt < entities.size(); ++iEnt) {
+ GEntity* &entity = entities[iEnt];
+ if (entity->dim() != dim) continue;
+ for (int iEl = 0; iEl < entity->getNumMeshElements(); iEl++) { // Detect bad elements
+ MElement *element = entity->getMeshElement(iEl);
+ double d = 0.;
+ for (int iEdge = 0; iEdge < element->getNumEdges(); ++iEdge) {
+ MEdge e = element->getEdge(iEdge);
+ std::map<MEdge,double,Less_Edge>::iterator it = dist2Edge.find(e);
+ if (it != dist2Edge.end()) d += it->second;
+ }
+ for (int iFace = 0; iFace < element->getNumFaces(); ++iFace) {
+ MFace f = element->getFace(iFace);
+ std::map<MFace,double,Less_Face>::iterator it = dist2Face.find(f);
+ if (it != dist2Face.end()) d += it->second;
+ }
+ distances[element] = d;
+ }
+ }
+}
+
+
+double distanceToGeometry(GModel *gm, int distType, double tol, int meshDiscr, int geomDiscr)
+{
+ double maxDist = 0.;
+
+ for (GModel::eiter it = gm->firstEdge(); it != gm->lastEdge(); ++it) {
+ if ((*it)->geomType() == GEntity::Line) continue;
+ for (unsigned int i = 0; i < (*it)->lines.size(); i++) {
+ double dist;
+ switch (distType) {
+ case CADDIST_TAYLOR:
+ dist = taylorDistanceEdge((*it)->lines[i], *it);
+ break;
+ case CADDIST_FRECHET:
+ dist = discreteFrechetDistanceEdge((*it)->lines[i], *it,
+ tol, meshDiscr, geomDiscr);
+ break;
+ case CADDIST_HAUSFAST:
+ dist = discreteHausdorffDistanceFastEdge((*it)->lines[i], *it,
+ tol, meshDiscr, geomDiscr);
+ break;
+ case CADDIST_HAUSBRUTE:
+ dist = discreteHausdorffDistanceBruteEdge((*it)->lines[i], *it,
+ tol, meshDiscr, geomDiscr);
+ break;
+ default:
+ Msg::Error("Wrong CAD distance type in distanceToGeometry");
+ break;
+ }
+ maxDist = std::max(dist, maxDist);
+ }
+ }
+// printf("DISTANCE TO GEOMETRY : 1D PART %22.15E\n",Obj);
+
+ if (distType == CADDIST_TAYLOR) {
+ for(GModel::fiter it = gm->firstFace(); it != gm->lastFace(); ++it) {
+ if ((*it)->geomType() == GEntity::Plane) continue;
+ for (unsigned int i = 0; i < (*it)->triangles.size(); i++) {
+ maxDist = std::max(taylorDistanceFace((*it)->triangles[i], *it), maxDist);
+ }
+ }
+ }
+// printf("DISTANCE TO GEOMETRY : 1D AND 2D PART %22.15E\n",Obj);
+
+ return maxDist;
+}
diff --git a/contrib/HighOrderMeshOptimizer/CADDistances.h b/contrib/HighOrderMeshOptimizer/CADDistances.h
new file mode 100644
index 0000000000000000000000000000000000000000..fb35f138b9a0f689f59e80938c753fe08ec3c620
--- /dev/null
+++ b/contrib/HighOrderMeshOptimizer/CADDistances.h
@@ -0,0 +1,43 @@
+// TODO: Header
+
+#ifndef _CADDISTANCES_H_
+#define _CADDISTANCES_H_
+
+#include "fullMatrix.h"
+#include <vector>
+#include <map>
+
+
+class GradientBasis;
+class MLine;
+class MElement;
+class GEdge;
+class GFace;
+class SVector3;
+class GModel;
+
+
+enum { CADDIST_GEN, CADDIST_DECASTELJAU };
+enum { CADDIST_TAYLOR, CADDIST_FRECHET, CADDIST_HAUSFAST, CADDIST_HAUSBRUTE };
+
+
+double discreteFrechetDistanceEdge(MLine *l, GEdge *ed, double tol,
+ int meshDiscr = CADDIST_GEN, int geomDiscr = CADDIST_GEN);
+double discreteHausdorffDistanceBruteEdge(MLine *l, GEdge *ed, double tol,
+ int meshDiscr = CADDIST_GEN, int geomDiscr = CADDIST_GEN);
+double discreteHausdorffDistanceFastEdge(MLine *l, GEdge *ed, double tol,
+ int meshDiscr = CADDIST_GEN, int geomDiscr = CADDIST_GEN);
+double taylorDistanceSq1D(const GradientBasis *gb, const fullMatrix<double> &nodesXYZ,
+ const std::vector<SVector3> &tanCAD);
+double taylorDistanceSq2D(const GradientBasis *gb, const fullMatrix<double> &nodesXYZ,
+ const std::vector<SVector3> &normCAD);
+double taylorDistanceEdge(MLine *l, GEdge *ge);
+double taylorDistanceFace(MElement *el, GFace *gf);
+
+void distanceFromElementsToGeometry(GModel *gm, int dim, std::map<MElement*, double> &distances);
+double distanceToGeometry(GModel *gm, int distType = CADDIST_TAYLOR, double tol = 1e-3,
+ int meshDiscr = CADDIST_DECASTELJAU,
+ int geomDiscr = CADDIST_DECASTELJAU);
+
+
+#endif /* _CADDISTANCES_H_ */
diff --git a/contrib/HighOrderMeshOptimizer/CMakeLists.txt b/contrib/HighOrderMeshOptimizer/CMakeLists.txt
index 29223ca3cd54a09fcdeafb1d3fa363c22560a162..768ccb4b4050497ec186b764e97d0a6d9925f3f8 100644
--- a/contrib/HighOrderMeshOptimizer/CMakeLists.txt
+++ b/contrib/HighOrderMeshOptimizer/CMakeLists.txt
@@ -13,6 +13,7 @@ set(SRC
SuperEl.cpp
OptHomElastic.cpp
OptHomFastCurving.cpp
+ CADDistances.cpp
)
file(GLOB_RECURSE HDR RELATIVE ${CMAKE_CURRENT_SOURCE_DIR} *.hpp)
diff --git a/contrib/HighOrderMeshOptimizer/OptHomCADDist.cpp b/contrib/HighOrderMeshOptimizer/OptHomCADDist.cpp
index 5e3d16b0076d037f17173129a13f6c44cd3d9b16..9cd6685fcfd2e1f8e33625d1066f128b3fd487e6 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomCADDist.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomCADDist.cpp
@@ -175,72 +175,72 @@ double MFaceGFaceDistance(MElement *el, GFace *gf)
}
-void distanceFromElementsToGeometry(GModel *gm, int dim, std::map<MElement*,double> &distances){
-
- std::map<MEdge,double,Less_Edge> dist2Edge;
- for (GModel::eiter it = gm->firstEdge(); it != gm->lastEdge(); ++it){
- if ((*it)->geomType() == GEntity::Line)continue;
- for (unsigned int i=0;i<(*it)->lines.size(); i++){
- double d = MLineGEdgeDistance ( (*it)->lines[i] , *it );
- MEdge e = (*it)->lines[i]->getEdge(0);
- dist2Edge[e] = d;
- }
- }
-
- // printf("DISTANCE TO GEOMETRY : 1D PART %22.15E\n",Obj);
-
- std::map<MFace,double,Less_Face> dist2Face;
- for(GModel::fiter it = gm->firstFace(); it != gm->lastFace(); ++it){
- if ((*it)->geomType() == GEntity::Plane)continue;
- for (unsigned int i=0;i<(*it)->triangles.size(); i++){
- double d = MFaceGFaceDistance ( (*it)->triangles[i] , *it );
- MFace f = (*it)->triangles[i]->getFace(0);
- dist2Face[f] = d;
- }
- }
-
- std::vector<GEntity*> entities;
- gm->getEntities(entities);
- for (int iEnt = 0; iEnt < entities.size(); ++iEnt) {
- GEntity* &entity = entities[iEnt];
- if (entity->dim() != dim) continue;
- for (int iEl = 0; iEl < entity->getNumMeshElements();iEl++) { // Detect bad elements
- MElement *element = entity->getMeshElement(iEl);
- double d = 0.0;
- for (int iEdge = 0; iEdge < element->getNumEdges(); ++iEdge) {
- MEdge e = element->getEdge(iEdge);
- std::map<MEdge,double,Less_Edge>::iterator it = dist2Edge.find(e);
- if(it != dist2Edge.end())d+=it->second;
- }
- for (int iFace = 0; iFace < element->getNumFaces(); ++iFace) {
- MFace f = element->getFace(iFace);
- std::map<MFace,double,Less_Face>::iterator it = dist2Face.find(f);
- if(it != dist2Face.end())d+=it->second;
- }
- distances[element] = d;
- }
- }
-}
-
-
-double distanceToGeometry(GModel *gm)
-{
- double Obj = 0.0;
-
- for (GModel::eiter it = gm->firstEdge(); it != gm->lastEdge(); ++it) {
- if ((*it)->geomType() == GEntity::Line) continue;
- for (unsigned int i=0;i<(*it)->lines.size(); i++)
- Obj = std::max(MLineGEdgeDistance((*it)->lines[i], *it), Obj);
- }
- printf("DISTANCE TO GEOMETRY : 1D PART %22.15E\n",Obj);
-
- for(GModel::fiter it = gm->firstFace(); it != gm->lastFace(); ++it) {
- if ((*it)->geomType() == GEntity::Plane) continue;
- for (unsigned int i=0;i<(*it)->triangles.size(); i++) {
- Obj = std::max(Obj,MFaceGFaceDistance( (*it)->triangles[i] , *it ));
- }
- }
- printf("DISTANCE TO GEOMETRY : 1D AND 2D PART %22.15E\n",Obj);
-
- return Obj;
-}
+//void distanceFromElementsToGeometry(GModel *gm, int dim, std::map<MElement*,double> &distances){
+//
+// std::map<MEdge,double,Less_Edge> dist2Edge;
+// for (GModel::eiter it = gm->firstEdge(); it != gm->lastEdge(); ++it){
+// if ((*it)->geomType() == GEntity::Line)continue;
+// for (unsigned int i=0;i<(*it)->lines.size(); i++){
+// double d = MLineGEdgeDistance ( (*it)->lines[i] , *it );
+// MEdge e = (*it)->lines[i]->getEdge(0);
+// dist2Edge[e] = d;
+// }
+// }
+//
+// // printf("DISTANCE TO GEOMETRY : 1D PART %22.15E\n",Obj);
+//
+// std::map<MFace,double,Less_Face> dist2Face;
+// for(GModel::fiter it = gm->firstFace(); it != gm->lastFace(); ++it){
+// if ((*it)->geomType() == GEntity::Plane)continue;
+// for (unsigned int i=0;i<(*it)->triangles.size(); i++){
+// double d = MFaceGFaceDistance ( (*it)->triangles[i] , *it );
+// MFace f = (*it)->triangles[i]->getFace(0);
+// dist2Face[f] = d;
+// }
+// }
+//
+// std::vector<GEntity*> entities;
+// gm->getEntities(entities);
+// for (int iEnt = 0; iEnt < entities.size(); ++iEnt) {
+// GEntity* &entity = entities[iEnt];
+// if (entity->dim() != dim) continue;
+// for (int iEl = 0; iEl < entity->getNumMeshElements();iEl++) { // Detect bad elements
+// MElement *element = entity->getMeshElement(iEl);
+// double d = 0.0;
+// for (int iEdge = 0; iEdge < element->getNumEdges(); ++iEdge) {
+// MEdge e = element->getEdge(iEdge);
+// std::map<MEdge,double,Less_Edge>::iterator it = dist2Edge.find(e);
+// if(it != dist2Edge.end())d+=it->second;
+// }
+// for (int iFace = 0; iFace < element->getNumFaces(); ++iFace) {
+// MFace f = element->getFace(iFace);
+// std::map<MFace,double,Less_Face>::iterator it = dist2Face.find(f);
+// if(it != dist2Face.end())d+=it->second;
+// }
+// distances[element] = d;
+// }
+// }
+//}
+//
+//
+//double distanceToGeometry(GModel *gm)
+//{
+// double Obj = 0.0;
+//
+// for (GModel::eiter it = gm->firstEdge(); it != gm->lastEdge(); ++it) {
+// if ((*it)->geomType() == GEntity::Line) continue;
+// for (unsigned int i=0;i<(*it)->lines.size(); i++)
+// Obj = std::max(MLineGEdgeDistance((*it)->lines[i], *it), Obj);
+// }
+// printf("DISTANCE TO GEOMETRY : 1D PART %22.15E\n",Obj);
+//
+// for(GModel::fiter it = gm->firstFace(); it != gm->lastFace(); ++it) {
+// if ((*it)->geomType() == GEntity::Plane) continue;
+// for (unsigned int i=0;i<(*it)->triangles.size(); i++) {
+// Obj = std::max(Obj,MFaceGFaceDistance( (*it)->triangles[i] , *it ));
+// }
+// }
+// printf("DISTANCE TO GEOMETRY : 1D AND 2D PART %22.15E\n",Obj);
+//
+// return Obj;
+//}
diff --git a/contrib/HighOrderMeshOptimizer/OptHomCADDist.h b/contrib/HighOrderMeshOptimizer/OptHomCADDist.h
index d91deb6f96d689f90426c076090aa2e2fe962092..2354ffa365437d08ab48d2bed718ddf0a798ef13 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomCADDist.h
+++ b/contrib/HighOrderMeshOptimizer/OptHomCADDist.h
@@ -24,8 +24,8 @@ double distToCAD2D(const GradientBasis *gb, const fullMatrix<double> &nodesXYZ,
double MLineGEdgeDistance (MLine *l, GEdge *ge);
double MFaceGFaceDistance(MElement *el, GFace *gf);
-void distanceFromElementsToGeometry(GModel *gm, int dim, std::map<MElement*, double> &distances);
-double distanceToGeometry(GModel *gm);
+//void distanceFromElementsToGeometry(GModel *gm, int dim, std::map<MElement*, double> &distances);
+//double distanceToGeometry(GModel *gm);
#endif /* OPTHOMCADDIST_H_ */
diff --git a/contrib/HighOrderMeshOptimizer/OptHomObjContribCADDist.h b/contrib/HighOrderMeshOptimizer/OptHomObjContribCADDist.h
index 7c31c4178e78af8372780b2da3bc1c67167d8fe0..fa4aed8d80368f28c5f60eccaa5b8fcf4f9f0d97 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomObjContribCADDist.h
+++ b/contrib/HighOrderMeshOptimizer/OptHomObjContribCADDist.h
@@ -7,11 +7,11 @@
template<class FuncType>
-class ObjContribCADDist : public ObjContrib, public FuncType
+class ObjContribCADDistSq : public ObjContrib, public FuncType
{
public:
- ObjContribCADDist(double weight, double refDist);
- virtual ~ObjContribCADDist() {}
+ ObjContribCADDistSq(double weight, double refDist);
+ virtual ~ObjContribCADDistSq() {}
virtual ObjContrib *copy() const;
virtual void initialize(Patch *mesh);
virtual bool fail() { return false; }
@@ -29,7 +29,7 @@ protected:
template<class FuncType>
-ObjContribCADDist<FuncType>::ObjContribCADDist(double weight, double refDist) :
+ObjContribCADDistSq<FuncType>::ObjContribCADDistSq(double weight, double refDist) :
ObjContrib("CADDist", FuncType::getNamePrefix()+"CADDist"),
_mesh(0), _weight(weight), _refDist(refDist)
{
@@ -37,24 +37,24 @@ ObjContribCADDist<FuncType>::ObjContribCADDist(double weight, double refDist) :
template<class FuncType>
-ObjContrib *ObjContribCADDist<FuncType>::copy() const
+ObjContrib *ObjContribCADDistSq<FuncType>::copy() const
{
- return new ObjContribCADDist<FuncType>(*this);
+ return new ObjContribCADDistSq<FuncType>(*this);
}
template<class FuncType>
-void ObjContribCADDist<FuncType>::initialize(Patch *mesh)
+void ObjContribCADDistSq<FuncType>::initialize(Patch *mesh)
{
_mesh = mesh;
- _mesh->initScaledCADDist(_refDist);
+ _mesh->initScaledCADDistSq(_refDist);
updateMinMax();
FuncType::initialize(_min, _max);
}
template<class FuncType>
-bool ObjContribCADDist<FuncType>::addContrib(double &Obj, alglib::real_1d_array &gradObj)
+bool ObjContribCADDistSq<FuncType>::addContrib(double &Obj, alglib::real_1d_array &gradObj)
{
_min = BIGVAL;
_max = -BIGVAL;
@@ -65,7 +65,7 @@ bool ObjContribCADDist<FuncType>::addContrib(double &Obj, alglib::real_1d_array
const int nVEl = _mesh->nNodBndEl(iBndEl);
double f;
std::vector<double> gradF(nVEl*bndDim);
- _mesh->scaledCADDistAndGradients(iBndEl, f, gradF);
+ _mesh->scaledCADDistSqAndGradients(iBndEl, f, gradF);
_min = std::min(_min, f);
_max = std::max(_max, f);
Obj += FuncType::compute(f) * _weight;
@@ -87,7 +87,7 @@ bool ObjContribCADDist<FuncType>::addContrib(double &Obj, alglib::real_1d_array
template<class FuncType>
-void ObjContribCADDist<FuncType>::updateMinMax()
+void ObjContribCADDistSq<FuncType>::updateMinMax()
{
_min = BIGVAL;
_max = -BIGVAL;
@@ -98,7 +98,7 @@ void ObjContribCADDist<FuncType>::updateMinMax()
const int nVEl = _mesh->nNodBndEl(iBndEl);
double f;
std::vector<double> dumGradF(nVEl*bndDim);
- _mesh->scaledCADDistAndGradients(iBndEl, f, dumGradF);
+ _mesh->scaledCADDistSqAndGradients(iBndEl, f, dumGradF);
_min = std::min(_min, f);
_max = std::max(_max, f);
}
diff --git a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
index 0cd51f319212a06a09d96bc917c769f3ad2b992c..9f6feb58ec450bc0a016ebc5e1af748f668fc035 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
@@ -42,6 +42,7 @@
#include "MHexahedron.h"
#include "MPrism.h"
#include "MLine.h"
+#include "CADDistances.h"
#include "OS.h"
#include <stack>
@@ -586,10 +587,10 @@ static void optimizeOneByOne
#include "OptHomIntegralBoundaryDist.h"
-double ComputeDistanceToGeometry (GModel* gm)
-{
- return distanceToGeometry(gm);
-}
+//double ComputeDistanceToGeometry (GModel* gm)
+//{
+// return distanceToGeometry(gm);
+//}
double ComputeDistanceToGeometry (GEntity *ge , int distanceDefinition, double tolerance)
@@ -691,7 +692,7 @@ void HighOrderMeshOptimizer(GModel *gm, OptHomParameters &p)
#include "GEntity.h"
//#include "MElement.h"
//#include "OptHomRun.h"
-#include "OptHomCADDist.h"
+#include "CADDistances.h"
#include "MeshOptCommon.h"
#include "MeshOptObjContribFunc.h"
#include "MeshOptObjContrib.h"
@@ -756,11 +757,13 @@ double HOPatchDefParameters::bndElBadness(MElement *el, GEntity* gEnt) const
if (optCAD) {
if (el->getType() == TYPE_LIN) { // 2D
if (gEnt->geomType() != GEntity::Line) // Straight geometric line -> no distance
- return MLineGEdgeDistance(static_cast<MLine*>(el), gEnt->cast2Edge());
+ return optCADDistMax -
+ taylorDistanceEdge(static_cast<MLine*>(el), gEnt->cast2Edge());
}
else { // 3D
if (gEnt->geomType() != GEntity::Plane) // Straight geometric plance -> no distance
- return MFaceGFaceDistance(el, gEnt->cast2Face());
+ return optCADDistMax -
+ taylorDistanceFace(el, gEnt->cast2Face());
}
}
return 1.;
@@ -803,15 +806,19 @@ void HighOrderMeshOptimizerNew(GModel *gm, OptHomParameters &p)
minJacBarFunc.setTarget(p.BARRIER_MIN, 1.);
ObjContribScaledJac<ObjContribFuncBarrierFixMinMovMax> minMaxJacBarFunc(1.);
minMaxJacBarFunc.setTarget(p.BARRIER_MAX, 1.);
- ObjContribCADDist<ObjContribFuncSimpleTargetMax> CADDistFunc(p.optCADWeight, p.optCADDistMax);
- CADDistFunc.setTarget(1.);
+ ObjContribCADDistSq<ObjContribFuncSimpleTargetMax> CADDistFunc(p.optCADWeight, p.optCADDistMax);
+ CADDistFunc.setTarget(0.);
+// ObjContribCADDistSq<ObjContribFuncBarrierMovMax> CADDistFunc(p.optCADWeight, p.optCADDistMax);
+// CADDistFunc.setTarget(1., 0.);
+ ObjContribScaledJac<ObjContribFuncBarrierFixMin> minJacFixBarFunc(1.);
+ minJacFixBarFunc.setTarget(p.BARRIER_MIN, 1.);
MeshOptPass minJacPass;
minJacPass.maxParamUpdates = p.optPassMax;
minJacPass.maxOptIter = p.itMax;
minJacPass.contrib.push_back(&nodeDistFunc);
minJacPass.contrib.push_back(&minJacBarFunc);
- if (p.optCAD) minJacPass.contrib.push_back(&CADDistFunc);
+// if (p.optCAD) minJacPass.contrib.push_back(&CADDistFunc);
par.pass.push_back(minJacPass);
if (p.BARRIER_MAX > 0.) {
@@ -820,10 +827,20 @@ void HighOrderMeshOptimizerNew(GModel *gm, OptHomParameters &p)
minMaxJacPass.maxOptIter = p.itMax;
minMaxJacPass.contrib.push_back(&nodeDistFunc);
minMaxJacPass.contrib.push_back(&minMaxJacBarFunc);
- if (p.optCAD) minMaxJacPass.contrib.push_back(&CADDistFunc);
+// if (p.optCAD) minMaxJacPass.contrib.push_back(&CADDistFunc);
par.pass.push_back(minMaxJacPass);
}
+ if (p.optCAD) {
+ MeshOptPass maxCADDistPass;
+ maxCADDistPass.maxParamUpdates = p.optPassMax;
+ maxCADDistPass.maxOptIter = p.itMax;
+ maxCADDistPass.contrib.push_back(&nodeDistFunc);
+ maxCADDistPass.contrib.push_back(&minJacFixBarFunc);
+ maxCADDistPass.contrib.push_back(&CADDistFunc);
+ par.pass.push_back(maxCADDistPass);
+ }
+
meshOptimizer(gm, par);
p.CPU = par.CPU;
diff --git a/contrib/HighOrderMeshOptimizer/OptHomRun.h b/contrib/HighOrderMeshOptimizer/OptHomRun.h
index 56befba458bd5793e29dba930aa6df0bf8ddb661..df3e4601612f30cffa27f3c67e0658a2555c45cd 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomRun.h
+++ b/contrib/HighOrderMeshOptimizer/OptHomRun.h
@@ -77,6 +77,6 @@ void HighOrderMeshOptimizer(GModel *gm, OptHomParameters &p);
void HighOrderMeshOptimizerNew(GModel *gm, OptHomParameters &p);
// distanceDefinition 1) Hausdorff 2) Area/Length 3) Frechet (not done)
double ComputeDistanceToGeometry (GEntity *ge , int distanceDefinition,double tolerance) ;
-double ComputeDistanceToGeometry (GModel*gm);
+//double ComputeDistanceToGeometry (GModel*gm);
#endif
diff --git a/contrib/MeshOptimizer/MeshOptObjContribFunc.cpp b/contrib/MeshOptimizer/MeshOptObjContribFunc.cpp
index cd5757325c612de14504b348183341504b061d07..ca67b9630f997e0a5e540025a48a0d37f730e847 100644
--- a/contrib/MeshOptimizer/MeshOptObjContribFunc.cpp
+++ b/contrib/MeshOptimizer/MeshOptObjContribFunc.cpp
@@ -86,6 +86,13 @@ bool ObjContribFuncBarrierMovMax::stagnated(double vMin, double vMax)
}
+void ObjContribFuncBarrierFixMin::initialize(double vMin, double vMax)
+{
+ const double margin = (vMin == 0.) ? _defaultMargin : MARGINCOEFF*fabs(vMin);
+ _barrier = vMin - margin;
+}
+
+
void ObjContribFuncBarrierFixMinMovMax::initialize(double vMin, double vMax)
{
ObjContribFuncBarrierMovMax::initialize(vMin, vMax);
diff --git a/contrib/MeshOptimizer/MeshOptObjContribFunc.h b/contrib/MeshOptimizer/MeshOptObjContribFunc.h
index ec12b3fde04ff8bf572bb51e208c19364f541c0f..3f31118654f94347c290ef7e2a766c868e4cd4da 100644
--- a/contrib/MeshOptimizer/MeshOptObjContribFunc.h
+++ b/contrib/MeshOptimizer/MeshOptObjContribFunc.h
@@ -13,7 +13,7 @@ protected:
void initialize(double vMin, double vMax) {}
void updateParameters(double vMin, double vMax) {}
bool targetReached(double vMin, double vMax) { return true; }
- bool stagnated(double vMin, double vMax) { return false; };
+ bool stagnated(double vMin, double vMax) { return false; }
double compute(double v) { return v; }
double computeDiff(double v) { return 1.; }
};
@@ -79,6 +79,19 @@ protected:
};
+class ObjContribFuncBarrierFixMin : public ObjContribFuncBarrier
+{
+protected:
+ std::string getNamePrefix() { return "BarrierFixMin"; }
+ void initialize(double vMin, double vMax);
+ void updateParameters(double vMin, double vMax) {}
+ bool targetReached(double vMin, double vMax) { return (vMin >= _barrier); }
+ bool stagnated(double vMin, double vMax) { return false; }
+ inline double compute(double v);
+ inline double computeDiff(double v);
+};
+
+
class ObjContribFuncBarrierFixMinMovMax : public ObjContribFuncBarrierMovMax
{
protected:
@@ -134,6 +147,20 @@ inline double ObjContribFuncBarrierMovMax::computeDiff(double v)
}
+inline double ObjContribFuncBarrierFixMin::compute(double v)
+{
+ if (v > _barrier) return logBarrier(v, _barrier, _opt);
+ else return 1e300;
+}
+
+
+inline double ObjContribFuncBarrierFixMin::computeDiff(double v)
+{
+ if (v > _barrier) return diffLogBarrier(v, _barrier, _opt);
+ else return -1e300;
+}
+
+
inline double ObjContribFuncBarrierFixMinMovMax::compute(double v)
{
double obj;
diff --git a/contrib/MeshOptimizer/MeshOptPatch.cpp b/contrib/MeshOptimizer/MeshOptPatch.cpp
index 0f657ddf023bcc13e72badf983a200173d1d9d4e..32a8dfb09fb2195b477b3f38c1c1eb7e479cebd7 100644
--- a/contrib/MeshOptimizer/MeshOptPatch.cpp
+++ b/contrib/MeshOptimizer/MeshOptPatch.cpp
@@ -36,7 +36,8 @@
#include "BasisFactory.h"
#include "CondNumBasis.h"
#include "OptHomIntegralBoundaryDist.h"
-#include "OptHomCADDist.h"
+//#include "OptHomCADDist.h"
+#include "CADDistances.h"
#include "qualityMeasures.h"
#include "MeshOptPatch.h"
@@ -344,9 +345,9 @@ void Patch::initMetricMin()
}
-void Patch::initScaledCADDist(double refCADDist)
+void Patch::initScaledCADDistSq(double refCADDist)
{
- _invRefCADDist = 1./refCADDist;
+ _invRefCADDistSq = 1./(refCADDist*refCADDist);
}
@@ -565,8 +566,8 @@ bool Patch::bndDistAndGradients(int iEl, double &f, std::vector<double> &gradF,
}
-void Patch::scaledCADDistAndGradients(int iBndEl, double &scaledDist,
- std::vector<double> &gradScaledDist)
+void Patch::scaledCADDistSqAndGradients(int iBndEl, double &scaledDist,
+ std::vector<double> &gradScaledDist)
{
const std::vector<int> &iV = _bndEl2V[iBndEl], &iFV = _bndEl2FV[iBndEl];
const int nV = iV.size();
@@ -604,8 +605,7 @@ void Patch::scaledCADDistAndGradients(int iBndEl, double &scaledDist,
tanCAD[i] = ge->firstDer(tCAD); // Compute tangent at vertex
tanCAD[i].normalize(); // Normalize tangent
}
- const double edLength = _xyz[iV[0]].distance(_xyz[iV[1]]); // Get edge length
- scaledDist = _invRefCADDist * distToCAD1D(gb, nodesXYZ, tanCAD, edLength);
+ scaledDist = _invRefCADDistSq * taylorDistanceSq1D(gb, nodesXYZ, tanCAD);
for (int i=0; i<nV; i++) {
const int &iFVi = iFV[i];
if (iFVi < 0) continue; // Skip if not free vertex
@@ -616,7 +616,8 @@ void Patch::scaledCADDistAndGradients(int iBndEl, double &scaledDist,
nodesXYZ(i, 0) = gp.x(); nodesXYZ(i, 1) = gp.y(); nodesXYZ(i, 2) = gp.z();
tanCAD[i] = ge->firstDer(tCAD); // New tangent to CAD at perturbed node
tanCAD[i].normalize(); // Normalize new tangent
- double sDistDiff = _invRefCADDist * distToCAD1D(gb, nodesXYZ, tanCAD, edLength); // Compute distance with perturbed node
+ const double sDistDiff = _invRefCADDistSq *
+ taylorDistanceSq1D(gb, nodesXYZ, tanCAD); // Compute distance with perturbed node
gradScaledDist[i] = (sDistDiff-scaledDist) / eps; // Compute gradient
nodesXYZ(i, 0) = xS; nodesXYZ(i, 1) = yS; nodesXYZ(i, 2) = zS; // Restore coord. of perturbed node
tanCAD[i] = tanCADS; // Restore tan. to CAD at perturbed node
@@ -639,7 +640,7 @@ void Patch::scaledCADDistAndGradients(int iBndEl, double &scaledDist,
normCAD[i] = gf->normal(pCAD); // Compute normal at vertex
normCAD[i].normalize(); // Normalize normal
}
- scaledDist = _invRefCADDist * distToCAD2D(gb, nodesXYZ, normCAD);
+ scaledDist = _invRefCADDistSq * taylorDistanceSq2D(gb, nodesXYZ, normCAD);
for (int i=0; i<nV; i++) {
const int &iFVi = iFV[i];
if (iFVi < 0) continue; // Skip if not free vertex
@@ -650,14 +651,16 @@ void Patch::scaledCADDistAndGradients(int iBndEl, double &scaledDist,
nodesXYZ(i, 0) = gp0.x(); nodesXYZ(i, 1) = gp0.y(); nodesXYZ(i, 2) = gp0.z();
normCAD[i] = gf->normal(pCAD0); // New normal to CAD at perturbed node in 1st dir.
normCAD[i].normalize(); // Normalize new normal
- double sDistDiff0 = _invRefCADDist * distToCAD2D(gb, nodesXYZ, normCAD); // Compute distance with perturbed node in 1st dir.
+ const double sDistDiff0 = _invRefCADDistSq *
+ taylorDistanceSq2D(gb, nodesXYZ, normCAD); // Compute distance with perturbed node in 1st dir.
gradScaledDist[2*i] = (sDistDiff0-scaledDist) / eps0; // Compute gradient in 1st dir.
const SPoint2 pCAD1 = SPoint2(_uvw[iFVi].x(), _uvw[iFVi].y()+eps1); // New param. coord. of perturbed node in 2nd dir.
GPoint gp1 = gf->point(pCAD1); // New coord. of perturbed node in 2nd dir.
nodesXYZ(i, 0) = gp1.x(); nodesXYZ(i, 1) = gp1.y(); nodesXYZ(i, 2) = gp1.z();
normCAD[i] = gf->normal(pCAD1); // New normal to CAD at perturbed node in 2nd dir.
normCAD[i].normalize(); // Normalize new normal
- double sDistDiff1 = _invRefCADDist * distToCAD2D(gb, nodesXYZ, normCAD); // Compute distance with perturbed node in 2nd dir.
+ double sDistDiff1 = _invRefCADDistSq *
+ taylorDistanceSq2D(gb, nodesXYZ, normCAD); // Compute distance with perturbed node in 2nd dir.
gradScaledDist[2*i+1] = (sDistDiff1-scaledDist) / eps0; // Compute gradient in 2nd dir.
nodesXYZ(i, 0) = xS; nodesXYZ(i, 1) = yS; nodesXYZ(i, 2) = zS; // Restore coord. of perturbed node
normCAD[i] = tanCADS; // Restore tan. to CAD at perturbed node
diff --git a/contrib/MeshOptimizer/MeshOptPatch.h b/contrib/MeshOptimizer/MeshOptPatch.h
index e0d7fb1bc98647a8be321705d12d204a2565ab99..baa3b0b7b2dc4f780a487e80e8db2825c499d1bf 100644
--- a/contrib/MeshOptimizer/MeshOptPatch.h
+++ b/contrib/MeshOptimizer/MeshOptPatch.h
@@ -95,9 +95,9 @@ public:
void initMetricMin();
void metricMinAndGradients(int iEl, std::vector<double> &sJ, std::vector<double> &gSJ);
bool bndDistAndGradients(int iEl, double &f, std::vector<double> &gradF, double eps);
- void initScaledCADDist(double refCADDist);
- void scaledCADDistAndGradients(int iBndEl, double &scaledDist,
- std::vector<double> &gradScaledDist);
+ void initScaledCADDistSq(double refCADDist);
+ void scaledCADDistSqAndGradients(int iBndEl, double &scaledDist,
+ std::vector<double> &gradScaledDist);
// Mesh quality
inline const int &nIJacEl(int iEl) { return _nIJacEl[iEl]; }
@@ -158,7 +158,7 @@ private:
std::vector<MElement*> _bndEl; // Boundary elements
std::vector<std::vector<int> > _bndEl2V, _bndEl2FV; // Vertices & corresponding free vertices on the boundary elements
std::vector<GEntity*> _bndEl2Ent; // Geometric entities corresponding to the boundary elements
- double _invRefCADDist;
+ double _invRefCADDistSq;
void calcNormalEl2D(int iEl, NormalScaling scaling,
fullMatrix<double> &elNorm, bool ideal);
diff --git a/wrappers/gmshpy/gmshMesh.i b/wrappers/gmshpy/gmshMesh.i
index 361cfed3e9c4ef1aa06fb98b0df7cb4bda1579fb..92ae6e7c0aa88b78e633ae513c27bb35d2f6f953 100644
--- a/wrappers/gmshpy/gmshMesh.i
+++ b/wrappers/gmshpy/gmshMesh.i
@@ -19,6 +19,7 @@
#include "OptHomElastic.h"
#include "OptHomFastCurving.h"
#include "MeshQualityOptimizer.h"
+ #include "CADDistances.h"
#endif
#if defined(HAVE_METIS) || defined(HAVE_CHACO)
#include "meshPartition.h"
@@ -59,6 +60,7 @@ namespace std {
%include "OptHomElastic.h"
%include "OptHomFastCurving.h"
%include "MeshQualityOptimizer.h"
+%include "CADDistances.h"
#endif
#if defined(HAVE_METIS) || defined(HAVE_CHACO)
%include "meshPartition.h"