Fix case when the grids overlap a little

dune/contact/projections/closestpointprojection.hh

@@ -65,7 +65,7 @@ void ClosestPointProjection<BoundaryPatchType>::project(const BoundaryPatchType&
                                                         const BoundaryPatchType& mortar,
                                                         const ctype couplingDist)
 {
-    ctype eps = 1e-5;
+    ctype eps = 1e-11;
 
     // initialize
     this->obsDirections_.resize(nonmortar.numVertices());
@@ -131,19 +131,20 @@ void ClosestPointProjection<BoundaryPatchType>::project(const BoundaryPatchType&
         // if we found a feasible closest point save the distance and
         if (bestTri != -1) {
             this->obstacles_[localIdx] = dist;
-
-            // if obstacle distance is very small, use minus the outer mortar normal as direction
-            if (dist<1e-5) {
+            // if there is overlap, fix orientation of the direction
+            if (dist<this->overlap_ and (projDirection*mortarSegments[bestTri].unitOuterNormal)>eps) {
+                this->obsDirections_[localIdx] = projDirection;
+                this->obsDirections_[localIdx] *= -1;
+            // if obstacle distance is very small, use the inner mortar normal as direction
+            } else if (dist<eps) {
                 this->obsDirections_[localIdx] = mortarSegments[bestTri].unitOuterNormal;
                 this->obsDirections_[localIdx] *= -1;
             } else
                 this->obsDirections_[localIdx] = projDirection;
 
         } else {
-
             this->obstacles_[localIdx] = std::numeric_limits<ctype>::max();
             this->obsDirections_[localIdx] = averagedNonmortarNormals[globalIdx];
-
         }
     }
 }