diff --git a/contrib/HighOrderMeshOptimizer/CADDistances.cpp b/contrib/HighOrderMeshOptimizer/CADDistances.cpp
index 96b517c97d291f2a5a55bd4ce434f5a29f3507e2..980c5d196e358f3f55669be158cafd2a417a0a8f 100644
--- a/contrib/HighOrderMeshOptimizer/CADDistances.cpp
+++ b/contrib/HighOrderMeshOptimizer/CADDistances.cpp
@@ -353,6 +353,7 @@ double discreteDistance(MLine *l, GEdge *ed, double tol, int meshDiscr, int geom
double t0, t1;
reparamMeshVertexOnEdge(l->getVertex(0), ed, t0);
reparamMeshVertexOnEdge(l->getVertex(1), ed, t1);
+// if (t1 == 0.) t1 = 1.; // DBG: Workaround bug closed curves
parametricLineGEdge l1(ed, t0, t1);
l1.discretize(dpts1, ts1, tol, 5, 45);
}
@@ -370,7 +371,6 @@ double discreteDistance(MLine *l, GEdge *ed, double tol, int meshDiscr, int geom
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);
@@ -432,7 +432,7 @@ double taylorDistanceSq2D(const GradientBasis *gb, const fullMatrix<double> &nod
const std::vector<SVector3> &normCAD)
{
const int nV = nodesXYZ.size1();
- fullMatrix<double> dxyzdX(nV, 3),dxyzdY(nV, 3);
+ fullMatrix<double> dxyzdX(nV, 3), dxyzdY(nV, 3);
gb->getGradientsFromNodes(nodesXYZ, &dxyzdX, &dxyzdY, 0);
double distSq = 0.;
for (int i=0; i<nV; i++) {
@@ -512,8 +512,8 @@ void distanceFromElementsToGeometry(GModel *gm, int dim, std::map<MElement*,doub
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++){
+ 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;
@@ -575,7 +575,6 @@ double distanceToGeometry(GModel *gm, int distType, double tol, int meshDiscr, i
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) {
@@ -585,7 +584,6 @@ double distanceToGeometry(GModel *gm, int distType, double tol, int meshDiscr, i
}
}
}
-// printf("DISTANCE TO GEOMETRY : 1D AND 2D PART %22.15E\n",Obj);
return maxDist;
}
diff --git a/contrib/HighOrderMeshOptimizer/OptHomObjContribCADDist.h b/contrib/HighOrderMeshOptimizer/OptHomObjContribCADDist.h
index fa4aed8d80368f28c5f60eccaa5b8fcf4f9f0d97..23a9f7b52e6b532d0a64cf43459c77d4002ea9d2 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomObjContribCADDist.h
+++ b/contrib/HighOrderMeshOptimizer/OptHomObjContribCADDist.h
@@ -30,7 +30,7 @@ protected:
template<class FuncType>
ObjContribCADDistSq<FuncType>::ObjContribCADDistSq(double weight, double refDist) :
- ObjContrib("CADDist", FuncType::getNamePrefix()+"CADDist"),
+ ObjContrib("ScaledCADDistSq", FuncType::getNamePrefix()+"ScaledCADDistSq"),
_mesh(0), _weight(weight), _refDist(refDist)
{
}
@@ -72,11 +72,11 @@ bool ObjContribCADDistSq<FuncType>::addContrib(double &Obj, alglib::real_1d_arra
const double dFact = _weight * FuncType::computeDiff(f);
for (int i=0; i<nVEl; i++) {
const int iFVi = _mesh->bndEl2FV(iBndEl, i);
- if (iFVi >= 0) { // Skip if not free vertex
- if (bndDim == 1) gradObj[_mesh->indPCFV(iFVi, 0)] += gradF[i] * dFact; // 2D
- else { // 3D
- gradObj[_mesh->indPCFV(iFVi, 0)] += gradF[2*i] * dFact;
- gradObj[_mesh->indPCFV(iFVi, 1)] += gradF[2*i+1] * dFact;
+ if (iFVi >= 0) { // Skip if not free vertex
+ if (bndDim == 1) gradObj[_mesh->indPCFV(iFVi, 0)] += gradF[i] * dFact; // 2D
+ else { // 3D
+ gradObj[_mesh->indPCFV(iFVi, 0)] += gradF[2*i] * dFact; // Deriv. w.r.t. 1st param.coord (edge or face vertex)
+ if (_mesh->nPCFV(iFVi) > 1) gradObj[_mesh->indPCFV(iFVi, 1)] += gradF[2*i+1] * dFact; // Deriv. w.r.t. 2nd param. coord. (only if face vertex)
}
}
}
diff --git a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
index 9f6feb58ec450bc0a016ebc5e1af748f668fc035..51a05488bdd52b8034ac437e7040c75a1c04ce06 100644
--- a/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
+++ b/contrib/HighOrderMeshOptimizer/OptHomRun.cpp
@@ -806,10 +806,10 @@ void HighOrderMeshOptimizerNew(GModel *gm, OptHomParameters &p)
minJacBarFunc.setTarget(p.BARRIER_MIN, 1.);
ObjContribScaledJac<ObjContribFuncBarrierFixMinMovMax> minMaxJacBarFunc(1.);
minMaxJacBarFunc.setTarget(p.BARRIER_MAX, 1.);
- ObjContribCADDistSq<ObjContribFuncSimpleTargetMax> CADDistFunc(p.optCADWeight, p.optCADDistMax);
- CADDistFunc.setTarget(0.);
-// ObjContribCADDistSq<ObjContribFuncBarrierMovMax> CADDistFunc(p.optCADWeight, p.optCADDistMax);
-// CADDistFunc.setTarget(1., 0.);
+// 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.);
diff --git a/contrib/MeshOptimizer/MeshOptCommon.cpp b/contrib/MeshOptimizer/MeshOptCommon.cpp
index d0a1cbfde62ee05f59ad7ebdb1eb273029cc8cf6..4f82cba30d9451363b69f2f1474c2e13b1b9c87c 100644
--- a/contrib/MeshOptimizer/MeshOptCommon.cpp
+++ b/contrib/MeshOptimizer/MeshOptCommon.cpp
@@ -121,8 +121,9 @@ bool MeshOptPatchDef::testElInDist(const SPoint3 &P, double limDist,
const SPoint3 A = faceVert[0]->point();
const SPoint3 B = faceVert[1]->point();
const SPoint3 C = faceVert[2]->point();
- if (faceVert.size() == 3)
+ if (faceVert.size() == 3) {
if (testTriSphereIntersect(A, B, C, P, limDistSq)) return true;
+ }
else {
const SPoint3 D = faceVert[3]->point();
if (testTriSphereIntersect(A, B, C, P, limDistSq) ||
diff --git a/contrib/MeshOptimizer/MeshOptPatch.cpp b/contrib/MeshOptimizer/MeshOptPatch.cpp
index 32a8dfb09fb2195b477b3f38c1c1eb7e479cebd7..1aff778db46a6bbb3f214f178489ddc71aeadd80 100644
--- a/contrib/MeshOptimizer/MeshOptPatch.cpp
+++ b/contrib/MeshOptimizer/MeshOptPatch.cpp
@@ -633,37 +633,62 @@ void Patch::scaledCADDistSqAndGradients(int iBndEl, double &scaledDist,
const int &iVi = iV[i], &iFVi = iFV[i];
MVertex* &vert = _vert[iVi];
SPoint2 pCAD;
- if ((iFVi >= 0) && (vert->onWhat() == gf)) // If free vertex and on surface, ...
- pCAD = SPoint2(_uvw[iFVi].x(), _uvw[iFVi].y()); // ... get stored param. coord.
+ if (iFVi >= 0) { // If free vertex...
+ if (vert->onWhat() == gf) // If on surface, ...
+ pCAD = SPoint2(_uvw[iFVi].x(), _uvw[iFVi].y()); // ... get stored param. coord.
+ else { // Otherwise, reparametrize on surface
+ const GEdge *ge = vert->onWhat()->cast2Edge();
+ pCAD = ge->reparamOnFace(gf, _uvw[iFVi].x(), 1);
+ }
+ }
else
- reparamMeshVertexOnFace(vert, gf, pCAD); // Otherwise, get param. coord. from CAD.
+ reparamMeshVertexOnFace(vert, gf, pCAD); // If not free vertex, reparametrize on surface
normCAD[i] = gf->normal(pCAD); // Compute normal at vertex
normCAD[i].normalize(); // Normalize normal
}
scaledDist = _invRefCADDistSq * taylorDistanceSq2D(gb, nodesXYZ, normCAD);
+// std::cout << "DBGTT: bnd el. " << _bndEl[iBndEl]->getNum() << ": scaledDist = " << scaledDist << "\n";
for (int i=0; i<nV; i++) {
- const int &iFVi = iFV[i];
+ const int &iVi = iV[i], &iFVi = iFV[i];
if (iFVi < 0) continue; // Skip if not free vertex
const double xS = nodesXYZ(i, 0), yS = nodesXYZ(i, 1), zS = nodesXYZ(i, 2); // Save coord. of perturbed node for FD
- const SVector3 tanCADS = normCAD[i]; // Save normal to CAD at perturbed node
- const SPoint2 pCAD0 = SPoint2(_uvw[iFVi].x()+eps0, _uvw[iFVi].y()); // New param. coord. of perturbed node in 1st dir.
- GPoint gp0 = gf->point(pCAD0); // New coord. of perturbed node in 1st dir.
- 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
- 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 = _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.
+ const SVector3 normCADS = normCAD[i]; // Save normal to CAD at perturbed node
+ if (_nPCFV[iFVi] == 2) { // Vertex classified on surface, 2D gradient
+ const SPoint2 pCAD0 = SPoint2(_uvw[iFVi].x()+eps0, _uvw[iFVi].y()); // New param. coord. of perturbed node in 1st dir.
+ GPoint gp0 = gf->point(pCAD0); // New coord. of perturbed node in 1st dir.
+ 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
+ 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 = _invRefCADDistSq *
+ taylorDistanceSq2D(gb, nodesXYZ, normCAD); // Compute distance with perturbed node in 2nd dir.
+ gradScaledDist[2*i+1] = (sDistDiff1-scaledDist) / eps1; // Compute gradient in 2nd dir.
+ }
+ else if (_nPCFV[iFVi] == 1) { // Vertex classified on edge, 1D gradient
+ MVertex* &vert = _vert[iVi];
+ const GEdge *ge = vert->onWhat()->cast2Edge();
+ const Range<double> parBounds = ge->parBounds(0);
+ const double eps = 1.e-6 * (parBounds.high()-parBounds.low());
+ const double tCAD = _uvw[iFVi].x() + eps; // New param. coord. of perturbed node
+ GPoint gp = ge->point(tCAD); // New coord. of perturbed node
+ nodesXYZ(i, 0) = gp.x(); nodesXYZ(i, 1) = gp.y(); nodesXYZ(i, 2) = gp.z();
+ SPoint2 pCAD = gf->parFromPoint(SPoint3(gp.x(), gp.y(), gp.z()), true); // Get param. coord. of perturbed node in face from CAD
+ normCAD[i] = gf->normal(pCAD); // New normal to CAD at perturbed node
+ normCAD[i].normalize(); // Normalize new normal
+ const double sDistDiff = _invRefCADDistSq *
+ taylorDistanceSq2D(gb, nodesXYZ, normCAD); // Compute distance with perturbed node
+ gradScaledDist[2*i] = (sDistDiff-scaledDist) / eps; // Compute gradient
+ }
+ else std::cout << "DBGTT: Inconsistent _nPCFV(iFVi), vert. " << _vert[iVi]->getNum() << "\n";
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
+ normCAD[i] = normCADS; // Restore tan. to CAD at perturbed node
}
}
}
diff --git a/contrib/MeshOptimizer/MeshOptimizer.cpp b/contrib/MeshOptimizer/MeshOptimizer.cpp
index 3365908d13e72f816483ce3abf214adc2ab6b0ae..b8445e941ed76643bc42aec86adffef158f8d2fb 100644
--- a/contrib/MeshOptimizer/MeshOptimizer.cpp
+++ b/contrib/MeshOptimizer/MeshOptimizer.cpp
@@ -590,7 +590,7 @@ void meshOptimizer(GModel *gm, MeshOptParameters &par)
elElMap::iterator bndElIt = el2BndEl.find(el);
if (bndElIt != el2BndEl.end()) {
MElement* &bndEl = bndElIt->second;
- if (par.patchDef->bndElBadness(bndEl, bndEl2Ent[bndEl])) badElts.insert(el);
+ if (par.patchDef->bndElBadness(bndEl, bndEl2Ent[bndEl]) < 0.) badElts.insert(el);
}
}
}