Fixed bugs in treatment of boundary elements in mesh optimizer. Fixed...

Fixed bugs in treatment of boundary elements in mesh optimizer. Fixed optimization of 2D CAD distance. Changed CAD distance contribution (now moving barrier).

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);
@@ -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;
 }

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)
 {
 }
@@ -75,8 +75,8 @@ bool ObjContribCADDistSq<FuncType>::addContrib(double &Obj, alglib::real_1d_arra
       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;
+          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)
         }
       }
     }

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.);
 

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) ||

contrib/MeshOptimizer/MeshOptPatch.cpp

@@ -633,19 +633,27 @@ 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, ...
+      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 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();
@@ -661,9 +669,26 @@ void Patch::scaledCADDistSqAndGradients(int iBndEl, double &scaledDist,
         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.
+        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
     }
   }
 }

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);
           }
         }
       }