Merging develop into master

apf/apfMesh.cc

@@ -175,6 +175,11 @@ bool Mesh::canSnap()
   return gmi_can_eval(getModel());
 }
 
+bool Mesh::canGetClosestPoint()
+{
+  return gmi_can_get_closest_point(getModel());
+}
+
 bool Mesh::canGetModelNormal()
 {
   return gmi_has_normal(getModel());
@@ -252,6 +257,16 @@ bool Mesh::isInClosureOf(ModelEntity* g, ModelEntity* target){
   return (res == 1) ? true : false;
 }
 
+bool Mesh::isOnModel(ModelEntity* g, Vector3 p, double scale)
+{
+  Vector3 to;
+  double param[2];
+  gmi_ent* c = (gmi_ent*)g;
+  gmi_closest_point(getModel(), c, &p[0], &to[0], param);
+  double ratio = (to - p).getLength() / scale;
+  return ratio < 0.001;
+}
+
 void Mesh::getPoint(MeshEntity* e, int node, Vector3& p)
 {
   getVector(coordinateField,e,node,p);

apf/apfMesh.h

@@ -285,6 +285,8 @@ class Mesh
     ModelEntity* findModelEntity(int type, int tag);
     /** \brief return true if the geometric model supports snapping */
     bool canSnap();
+    /** \brief return true if the geometric model supports get closest point */
+    bool canGetClosestPoint();
     /** \brief return true if the geometric model supports normal computation */
     bool canGetModelNormal();
     /** \brief evaluate parametric coordinate (p) as a spatial point (x) */
@@ -310,6 +312,8 @@ class Mesh
       Vector3 const& param, Vector3& x);
     /** \brief checks if g is in the closure of the target */
     bool isInClosureOf(ModelEntity* g, ModelEntity* target);
+    /** \brief checks if p is on model g */
+    bool isOnModel(ModelEntity* g, Vector3 p, double scale);
     /** \brief get the distribution of the mesh's coordinate field */
     FieldShape* getShape() const;
     /** \brief get the mesh's coordinate field */

crv/crvBezier.h

@@ -131,7 +131,7 @@ void getGregoryBlendedTransformationCoefficients(int blend, int type,
     apf::NewArray<double>& c);
 
 /** \brief a per entity version of above */
-void snapToInterpolate(apf::Mesh2* m, apf::MeshEntity* e);
+void snapToInterpolate(apf::Mesh2* m, apf::MeshEntity* e, bool isNew = false);
 
 /** \brief compute the matrix to transform between Bezier and Lagrange Points
     \details this is a support function, not actual ever needed.

crv/crvCurveMesh.cc

@@ -44,7 +44,7 @@ void convertInterpolationPoints(apf::Mesh2* m, apf::MeshEntity* e,
   apf::destroyElement(elem);
 }
 
-void snapToInterpolate(apf::Mesh2* m, apf::MeshEntity* e)
+void snapToInterpolate(apf::Mesh2* m, apf::MeshEntity* e, bool isNew)
 {
   PCU_ALWAYS_ASSERT(m->canSnap());
   int type = m->getType(e);
@@ -56,9 +56,18 @@ void snapToInterpolate(apf::Mesh2* m, apf::MeshEntity* e)
     m->setPoint(e,0,pt);
     return;
   }
+  // e is an edge or a face
+  // either way, get a length-scale by computing
+  // the distance b/w first two downward verts
+  apf::MeshEntity* down[12];
+  m->getDownward(e, 0, down);
+  apf::Vector3 p0, p1;
+  m->getPoint(down[0], 0, p0);
+  m->getPoint(down[1], 0, p1);
+  double lengthScale = (p1 - p0).getLength();
   apf::FieldShape * fs = m->getShape();
   int non = fs->countNodesOn(type);
-  apf::Vector3 p, xi, pt(0,0,0);
+  apf::Vector3 p, xi, pt0, pt(0,0,0);
   for(int i = 0; i < non; ++i){
     apf::ModelEntity* g = m->toModel(e);
     fs->getNodeXi(type,i,xi);
@@ -67,7 +76,13 @@ void snapToInterpolate(apf::Mesh2* m, apf::MeshEntity* e)
     else
       transferParametricOnTriSplit(m,e,xi,p);
     m->snapToModel(g,p,pt);
-    m->setPoint(e,i,pt);
+    if (isNew || !m->canGetClosestPoint()) {
+      m->setPoint(e,i,pt);
+      continue;
+    }
+    m->getPoint(e,i,pt0);
+    if (!m->isOnModel(g, pt0, lengthScale))
+      m->setPoint(e,i,pt);
   }
 }
 

crv/crvShapeHandler.cc

@@ -149,7 +149,7 @@ class BezierTransfer : public ma::SolutionTransfer
       for (int i = 0; i < ne; ++i){
         if(shouldSnap && midEdgeVerts[i] &&
             isBoundaryEntity(mesh,midEdgeVerts[i])){
-          snapToInterpolate(mesh,midEdgeVerts[i]);
+          snapToInterpolate(mesh,midEdgeVerts[i],true);
         }
       }
       int np = mesh->getShape()->getEntityShape(parentType)->countNodes();
@@ -173,7 +173,7 @@ class BezierTransfer : public ma::SolutionTransfer
 
           // do snapping here, inside refinement
           if (isBdryEnt && shouldSnap){
-            snapToInterpolate(mesh,newEntities[i]);
+            snapToInterpolate(mesh,newEntities[i],true);
           } else if (useLinear && !isBdryEnt) {
             // boundary entities that don't get snapped should interpolate
             if(childType == apf::Mesh::EDGE){
@@ -538,7 +538,7 @@ class BezierHandler : public ma::ShapeHandler
         if (newType != apf::Mesh::EDGE)
         {
           if (snap && P > 2){
-            snapToInterpolate(mesh,newEntities[i]);
+            snapToInterpolate(mesh,newEntities[i],true);
             convertInterpolationPoints(mesh,newEntities[i],n,ni,bt->coeffs[2]);
           } else {
             for (int j = 0; j < ni; ++j){
@@ -548,7 +548,7 @@ class BezierHandler : public ma::ShapeHandler
           }
         } else {
           if (snap){
-            snapToInterpolate(mesh,newEntities[i]);
+            snapToInterpolate(mesh,newEntities[i],true);
             convertInterpolationPoints(mesh,newEntities[i],n,ni,bt->coeffs[1]);
           } else {
             setLinearEdgePoints(mesh,newEntities[i]);

gmi/gmi.c

@@ -85,6 +85,11 @@ int gmi_can_eval(struct gmi_model* m)
   return m->ops->eval != NULL;
 }
 
+int gmi_can_get_closest_point(struct gmi_model* m)
+{
+  return m->ops->closest_point != NULL;
+}
+
 int gmi_has_normal(struct gmi_model* m)
 {
   return m->ops->normal != NULL;

gmi/gmi.h

@@ -145,6 +145,8 @@ struct gmi_ent* gmi_find(struct gmi_model* m, int dim, int tag);
 struct gmi_set* gmi_adjacent(struct gmi_model* m, struct gmi_ent* e, int dim);
 /** \brief check whether the model implements gmi_eval */
 int gmi_can_eval(struct gmi_model* m);
+/** \brief check whether the model gmi_closest_point */
+int gmi_can_get_closest_point(struct gmi_model* m);
 /** \brief check whether the model implements gmi_normal */
 int gmi_has_normal(struct gmi_model* m);
 /** \brief evaluate the parametric definition of a model boundary entity