diff --git a/dumux/geometry/geometryintersection.hh b/dumux/geometry/geometryintersection.hh index 5522ff36e45f2f01de23618d68df9c42f1274de2..a768b79dd999713eaa34929dcff35f16b7e7e19f 100644 --- a/dumux/geometry/geometryintersection.hh +++ b/dumux/geometry/geometryintersection.hh @@ -1316,15 +1316,80 @@ public: /*! * \brief Colliding two segments * \param geo1/geo2 The geometries to intersect - * \param intersection Container to store corners of intersection segment - * \note this overload is used when point-like intersections are seeked - * \todo implement this query + * \param intersection The intersection point + * \note This overload is used when point-like intersections are seeked */ template<class P = Policy, std::enable_if_t<P::dimIntersection == 0, int> = 0> static bool intersection(const Geometry1& geo1, const Geometry2& geo2, Intersection& intersection) { static_assert(int(dimworld) == int(Geometry2::coorddimension), "Can only collide geometries of same coordinate dimension"); - DUNE_THROW(Dune::NotImplemented, "segment-segment intersection detection for point-like intersections"); + + const auto v1 = geo1.corner(1) - geo1.corner(0); + const auto v2 = geo2.corner(1) - geo2.corner(0); + const auto ac = geo2.corner(0) - geo1.corner(0); + + const auto v1Norm = v1.two_norm(); + const auto eps = eps_*v1Norm; + const auto eps2 = eps*eps; + const auto eps3 = eps2*eps; + + const auto n = crossProduct(v1, v2); + + // first check if segments are parallel + using std::abs; + if ( n.two_norm2() < eps3*v1Norm ) + { + // check if they lie on the same line + const auto acNorm = ac.two_norm(); + if (abs(ac*v1) < acNorm*v1Norm) + return false; + + // they lie on the same line, + // if so, point intersection can only be one of the corners + const auto sp = v1*v2; + if ( sp < 0.0 ) + { + if ( acNorm < eps ) + { intersection = geo2.corner(0); return true; } + + if ( (geo2.corner(1) - geo1.corner(1)).two_norm2() < eps2 ) + { intersection = geo2.corner(1); return true; } + } + else + { + if ( (geo2.corner(1) - geo1.corner(0)).two_norm2() < eps2 ) + { intersection = geo2.corner(1); return true; } + + if ( (geo2.corner(0) - geo1.corner(1)).two_norm2() < eps2 ) + { intersection = geo2.corner(0); return true; } + } + + // no intersection + return false; + } + + // in-plane normal vector on v1 + auto v1Normal = crossProduct(v1, n); + v1Normal /= v1Normal.two_norm(); + + // intersection point on v2 in local coords + const auto t2 = - 1.0*(ac*v1Normal) / (v2*v1Normal); + + // check if the local coords are valid + if (t2 < -1.0*eps_ || t2 > 1.0 + eps_) + return false; + + // compute global coordinates + auto isPoint = geo2.global(t2); + + // check if point lies on v1 + if (intersectsPointGeometry(isPoint, geo1)) + { + intersection = std::move(isPoint); + return true; + } + + return false; } /*!