diff --git a/dumux/geometry/geometryintersection.hh b/dumux/geometry/geometryintersection.hh
index 5522ff36e45f2f01de23618d68df9c42f1274de2..f8ce13bec151b7cc52a9e45b9ebdf77689886f46 100644
--- a/dumux/geometry/geometryintersection.hh
+++ b/dumux/geometry/geometryintersection.hh
@@ -830,6 +830,190 @@ public:
}
};
+/*!
+ * \ingroup Geometry
+ * \brief A class for polygon--polygon intersections in 3d space
+ */
+template <class Geometry1, class Geometry2, class Policy>
+class GeometryIntersection<Geometry1, Geometry2, Policy, 3, 2, 2>
+{
+ enum { dimworld = 3 };
+ enum { dim1 = 2 };
+ enum { dim2 = 2 };
+
+public:
+ using ctype = typename Policy::ctype;
+ using Point = typename Policy::Point;
+ using Intersection = typename Policy::Intersection;
+
+private:
+ static constexpr ctype eps_ = 1.5e-7; // base epsilon for floating point comparisons
+ using ReferenceElementsGeo1 = typename Dune::ReferenceElements<ctype, dim1>;
+ using ReferenceElementsGeo2 = typename Dune::ReferenceElements<ctype, dim2>;
+
+public:
+ /*!
+ * \brief Colliding two polygons
+ * \note First we find the vertex candidates for the intersection region as follows:
+ * Add vertices of first polygon that are inside the second polygon
+ * Add intersections of polygon edges
+ * Remove duplicate points from the list
+ * Compute the convex hull polygon
+ * \param geo1/geo2 The geometries to intersect
+ * \param intersection Container to store the corner points of the polygon (as convex hull)
+ * \note This overload is used when polygon like intersections are seeked
+ */
+ template<class P = Policy, std::enable_if_t<P::dimIntersection == 2, 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");
+
+ // if the geometries are not parallel, there cannot be a polygon intersection
+ const auto a1 = geo1.corner(1) - geo1.corner(0);
+ const auto b1 = geo1.corner(2) - geo1.corner(0);
+ const auto n1 = crossProduct(a1, b1);
+
+ const auto a2 = geo2.corner(1) - geo2.corner(0);
+ const auto b2 = geo2.corner(2) - geo2.corner(0);
+ const auto n2 = crossProduct(a2, b2);
+
+ using std::max;
+ using std::sqrt;
+ const auto a1Norm2 = a1.two_norm2();
+ const auto b1Norm2 = b1.two_norm2();
+ const auto maxNorm2 = max(a1Norm2, b1Norm2);
+ const auto eps2 = maxNorm2*eps_;
+ const auto eps = sqrt(maxNorm2)*eps_;
+ const auto eps3 = eps*eps2;
+
+ // early exit if polygons are not parallel
+ using std::abs;
+ if (crossProduct(n1, n2).two_norm2() > eps2*maxNorm2)
+ return false;
+
+ // they are parallel, verify that they are on the same plane
+ auto d1 = geo2.corner(0) - geo1.corner(0);
+ if (d1.two_norm2() < eps2)
+ d1 = geo2.corner(1) - geo1.corner(0);
+
+ if (abs(d1*n2) > eps3)
+ return false;
+
+ // the candidate intersection points
+ std::vector<Point> points; points.reserve(8);
+
+ // add polygon1 corners that are inside polygon2
+ for (int i = 0; i < geo1.corners(); ++i)
+ if (intersectsPointGeometry(geo1.corner(i), geo2))
+ points.emplace_back(geo1.corner(i));
+
+ const auto numPoints1 = points.size();
+ const bool resultIsGeo1 = numPoints1 == geo1.corners();
+ if (!resultIsGeo1)
+ {
+ // add polygon2 corners that are inside polygon1
+ for (int i = 0; i < geo2.corners(); ++i)
+ if (intersectsPointGeometry(geo2.corner(i), geo1))
+ points.emplace_back(geo2.corner(i));
+
+ const bool resultIsGeo2 = (points.size() - numPoints1) == geo2.corners();
+ if (!resultIsGeo2)
+ {
+ const auto referenceElement1 = ReferenceElementsGeo1::general(geo1.type());
+ const auto referenceElement2 = ReferenceElementsGeo2::general(geo2.type());
+
+ // add intersections of edges
+ using SegGeometry = Dune::MultiLinearGeometry<ctype, 1, dimworld>;
+ using PointPolicy = IntersectionPolicy::PointPolicy<ctype, dimworld>;
+ for (int i = 0; i < referenceElement1.size(dim1-1); ++i)
+ {
+ const auto localEdgeGeom1 = referenceElement1.template geometry<dim1-1>(i);
+ const auto edge1 = SegGeometry(
+ Dune::GeometryTypes::line,
+ std::vector<Point>({
+ geo1.global(localEdgeGeom1.corner(0)),
+ geo1.global(localEdgeGeom1.corner(1))
+ })
+ );
+
+ for (int j = 0; j < referenceElement2.size(dim2-1); ++j)
+ {
+ const auto localEdgeGeom2 = referenceElement2.template geometry<dim2-1>(j);
+ const auto edge2 = SegGeometry(
+ Dune::GeometryTypes::line,
+ std::vector<Point>({
+ geo2.global(localEdgeGeom2.corner(0)),
+ geo2.global(localEdgeGeom2.corner(1))
+ })
+ );
+
+ using EdgeTest = GeometryIntersection<SegGeometry, SegGeometry, PointPolicy>;
+ typename EdgeTest::Intersection edgeIntersection;
+ if (EdgeTest::intersection(edge1, edge2, edgeIntersection))
+ points.emplace_back(edgeIntersection);
+ }
+ }
+ }
+ }
+
+ if (points.empty())
+ return false;
+
+ // remove duplicates
+ std::sort(points.begin(), points.end(), [eps] (const auto& a, const auto& b) -> bool
+ {
+ using std::abs;
+ return (abs(a[0]-b[0]) > eps ? a[0] < b[0]
+ : (abs(a[1]-b[1]) > eps ? a[1] < b[1]
+ : (a[2] < b[2])));
+ });
+
+ const auto squaredEps = eps*eps;
+ points.erase(std::unique(
+ points.begin(), points.end(),
+ [squaredEps] (const auto& a, const auto&b) { return (b-a).two_norm2() < squaredEps; }),
+ points.end()
+ );
+
+ // return false if we don't have at least three unique points
+ if (points.size() < 3)
+ return false;
+
+ // intersection polygon is convex hull of above points
+ intersection = grahamConvexHull<2>(points);
+ assert(!intersection.empty());
+ return true;
+ }
+
+ /*!
+ * \brief Colliding two polygons
+ * \param geo1/geo2 The geometries to intersect
+ * \param intersection Container to store the corners of intersection segment
+ * \note this overload is used when segment-like intersections are seeked
+ * \todo implement this query
+ */
+ template<class P = Policy, std::enable_if_t<P::dimIntersection == 1, 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, "Polygon-polygon intersection detection for segment-like intersections");
+ }
+
+ /*!
+ * \brief Colliding two polygons
+ * \param geo1/geo2 The geometries to intersect
+ * \param intersection The intersection point
+ * \note this overload is used when point-like intersections are seeked
+ * \todo implement this query
+ */
+ 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, "Polygon-polygon intersection detection for touching points");
+ }
+};
+
/*!
* \ingroup Geometry
* \brief A class for polyhedron--polygon intersection in 3d space
@@ -951,18 +1135,20 @@ public:
// remove duplicates
const auto eps = (geo1.corner(0) - geo1.corner(1)).two_norm()*eps_;
- std::sort(points.begin(), points.end(), [&eps](const auto& a, const auto& b) -> bool
+ const auto notEqual = [eps] (auto a, auto b) { using std::abs; return abs(b-a) > eps; };
+ std::sort(points.begin(), points.end(), [notEqual](const auto& a, const auto& b) -> bool
{
- using std::abs;
- return (abs(a[0]-b[0]) > eps ? a[0] < b[0] : (abs(a[1]-b[1]) > eps ? a[1] < b[1] : (a[2] < b[2])));
+ return (notEqual(a[0], b[0]) ? a[0] < b[0]
+ : (notEqual(a[1], b[1]) ? a[1] < b[1]
+ : (a[2] < b[2])));
});
- auto removeIt = std::unique(points.begin(), points.end(), [&eps](const auto& a, const auto&b)
- {
- return (b-a).two_norm() < eps;
- });
-
- points.erase(removeIt, points.end());
+ const auto squaredEps = eps*eps;
+ points.erase(std::unique(
+ points.begin(), points.end(),
+ [squaredEps] (const auto& a, const auto&b) { return (b-a).two_norm2() < squaredEps; }),
+ points.end()
+ );
// return false if we don't have more than three unique points
if (points.size() < 3) return false;
@@ -1316,15 +1502,69 @@ 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 v1Norm2 = v1.two_norm2();
+ const auto eps2 = eps_*v1Norm2;
+
+ const auto n = crossProduct(v1, v2);
+
+ // first check if segments are parallel
+ using std::abs;
+ if ( n.two_norm2() < eps2*v1Norm2 )
+ {
+ // check if they lie on the same line
+ if (crossProduct(v1, ac).two_norm2() > eps2)
+ 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 ( ac.two_norm2() < eps2 )
+ { 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
+ const auto v1Normal = crossProduct(v1, n);
+
+ // 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;
+
+ if (auto ip = geo2.global(t2); intersectsPointGeometry(ip, geo1))
+ { intersection = std::move(ip); return true; }
+
+ return false;
}
/*!
diff --git a/test/geometry/test_1d1d_intersection.cc b/test/geometry/test_1d1d_intersection.cc
index 0c16920668ecec0ba5e870b39d525ae399c7903e..ef5590f642f22dec0eeee9eeea24d1a749b5aaed 100644
--- a/test/geometry/test_1d1d_intersection.cc
+++ b/test/geometry/test_1d1d_intersection.cc
@@ -196,8 +196,8 @@ int main (int argc, char *argv[])
// test for dimWorld = 2
testPointIntersections<2>(returns);
- // TODO: implement and test for dimWorld = 3
- // testPointIntersections<3>(returns);
+ // test for dimWorld = 3
+ testPointIntersections<3>(returns);
// test segment intersection for dimWorld = 2
testSegmentIntersections<2>(returns);
diff --git a/test/geometry/test_2d2d_intersection.cc b/test/geometry/test_2d2d_intersection.cc
index a9951bc455045a6a62f2f072d35c9f6d3d23530d..426b591fe73a8d074ce535cab1b93e64572193c0 100644
--- a/test/geometry/test_2d2d_intersection.cc
+++ b/test/geometry/test_2d2d_intersection.cc
@@ -128,14 +128,68 @@ void testPolygonIntersections(std::vector<bool>& returns)
}
}
+void testParallelPolygons(std::vector<bool>& returns)
+{
+ using Point = Dune::FieldVector<double, 3>;
+
+ for (auto scaling : {1.0, 1e3, 1e12, 1e-12})
+ {
+ const double unit = 1.0*scaling;
+ const double offUnit = (1.0 + 1e-6)*scaling;
+
+ std::cout << "Test with scaling " << scaling << std::endl;
+ const auto tria1 = makeTriangle( Point{{0.0, 0.0, unit}},
+ Point{{unit, 0.0, unit}},
+ Point{{unit, unit, unit}} );
+ const auto tria2 = makeTriangle( Point{{0.0, 0.0, offUnit}},
+ Point{{0.0, unit, offUnit}},
+ Point{{unit, 0.0, offUnit}} );
+ returns.push_back(testPolygonIntersection<3>(tria1, tria2, false));
+ }
+
+ std::cout << std::endl;
+}
+
+void testNonParallelPolygons(std::vector<bool>& returns)
+{
+ using Point = Dune::FieldVector<double, 3>;
+
+ for (auto scaling : {1.0, 1e3, 1e12, 1e-12})
+ {
+ const double unit = 1.0*scaling;
+ const double offUnit = (1.0 + 1e-6)*scaling;
+
+ std::cout << "Test with scaling " << scaling << std::endl;
+ const auto tria1 = makeTriangle( Point{{0.0, 0.0, unit}},
+ Point{{unit, 0.0, unit}},
+ Point{{unit, unit, unit}} );
+ const auto tria2 = makeTriangle( Point{{0.0, 0.0, unit}},
+ Point{{0.0, unit, unit}},
+ Point{{unit, 0.0, offUnit}} );
+ returns.push_back(testPolygonIntersection<3>(tria1, tria2, false));
+ }
+
+ std::cout << std::endl;
+}
+
#endif
int main(int argc, char* argv[])
{
std::vector<bool> returns;
+
+ std::cout << "Testing intersections in 2d space" << std::endl;
testPolygonIntersections<2>(returns);
- // TODO: implement and test intersections in 3d
+ std::cout << "Testing intersecions in 3d space" << std::endl;
+ testPolygonIntersections<3>(returns);
+
+ std::cout << "Testing parallel polygons in 3d space" << std::endl;
+ testParallelPolygons(returns);
+
+ std::cout << "Testing non-parallel polygons in 3d space" << std::endl;
+ testNonParallelPolygons(returns);
+
// TODO: implement and test point and segment intersections
// determine the exit code
diff --git a/test/geometry/test_2d3d_intersection.cc b/test/geometry/test_2d3d_intersection.cc
index 7bd05ed525a24eb25b838711fd915ab388b4016c..bf06825c0d4b01dafcdafc10d1cbd12158902c31 100644
--- a/test/geometry/test_2d3d_intersection.cc
+++ b/test/geometry/test_2d3d_intersection.cc
@@ -48,157 +48,201 @@ void testSegTriangle(const Dune::FieldVector<double, dimworld>& a,
}
}
-int main(int argc, char* argv[])
+template<typename ctype>
+auto makeUnitCube(ctype scale = 1.0)
{
- constexpr int dimworld = 3;
- using Point = Dune::FieldVector<double, dimworld>;
-
- Point a{0.0, 0.0, 0.1};
- Point b{1.0, 0.0, 0.0};
- Point c{0.0, 1.0, 0.0};
-
- Point p0{0.5, 0.5, -1.0};
- Point q0{0.5, 0.5, 1.0};
-
- testSegTriangle(a, b, c, p0, q0);
- testSegTriangle(a, b, c, q0, p0);
-
- Point p1 = a;
- Point q1 = b;
+ using CubeGeometry = Dune::MultiLinearGeometry<ctype, 3, 3>;
+ using GlobalPosition = typename CubeGeometry::GlobalCoordinate;
+
+ const auto unit = scale*1.0;
+ return CubeGeometry(
+ Dune::GeometryTypes::cube(3),
+ std::vector<GlobalPosition>({
+ {0.0, 0.0, 0.0}, {unit, 0.0, 0.0}, {0.0, unit, 0.0}, {unit, unit, 0.0},
+ {0.0, 0.0, unit}, {unit, 0.0, unit}, {0.0, unit, unit}, {unit, unit, unit}
+ })
+ );
+}
- testSegTriangle(a, b, c, p1, q1, false);
- testSegTriangle(a, b, c, q1, p1, false);
+template<typename ctype>
+auto makeHorizontalQuad(ctype scale = 1.0, ctype zPositionUnscaled = 1.0)
+{
+ using QuadGeometry = Dune::MultiLinearGeometry<ctype, 2, 3>;
+ using GlobalPosition = typename QuadGeometry::GlobalCoordinate;
+
+ zPositionUnscaled *= scale;
+ const auto unit = scale*1.0;
+ return QuadGeometry(
+ Dune::GeometryTypes::cube(2),
+ std::vector<GlobalPosition>({
+ {0.0, 0.0, zPositionUnscaled}, {unit, 0.0, zPositionUnscaled},
+ {0.0, unit, zPositionUnscaled}, {unit, unit, zPositionUnscaled}
+ })
+ );
+}
- Point p2{0.0, 0.0, 0.0};
- Point q2{0.0, 0.0, 0.2};
+template<typename ctype>
+auto makeInclinedTriangle(ctype scale = 1.0)
+{
+ using TriaGeometry = Dune::MultiLinearGeometry<ctype, 2, 3>;
+ using GlobalPosition = typename TriaGeometry::GlobalCoordinate;
+ return TriaGeometry(
+ Dune::GeometryTypes::simplex(2),
+ std::vector<GlobalPosition>({
+ {-0.1*scale, -0.1*scale, 0.3*scale},
+ {1.1*scale, -0.1*scale, 0.3*scale},
+ {0.5*scale, 2.0*scale, 0.8*scale}
+ })
+ );
+}
- testSegTriangle(a, b, c, p2, q2);
+int main(int argc, char* argv[])
+{
+ constexpr int dimworld = 3;
+ using Point = Dune::FieldVector<double, dimworld>;
- using Geometry3D = Dune::MultiLinearGeometry<double, 3, dimworld>;
- using Geometry2D = Dune::MultiLinearGeometry<double, 2, dimworld>;
- using Test = Dumux::GeometryIntersection<Geometry3D, Geometry2D>;
- typename Test::Intersection intersectionPolygon;
+ for (auto scale : {1.0, 1e3, 1e12, 1e-12})
+ {
+ std::cout << "Testing scale " << scale << std::endl;
- std::vector<Dune::FieldVector<double, dimworld>> cubeCorners({
- {0.0, 0.0, 0.0}, {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {1.0, 1.0, 0.0},
- {0.0, 0.0, 1.0}, {1.0, 0.0, 1.0}, {0.0, 1.0, 1.0}, {1.0, 1.0, 1.0}
- });
+ const Point a{0.0*scale, 0.0*scale, 0.1*scale};
+ const Point b{1.0*scale, 0.0*scale, 0.0*scale};
+ const Point c{0.0*scale, 1.0*scale, 0.0*scale};
- std::vector<Dune::FieldVector<double, dimworld>> quadCorners({
- {0.0, 0.0, 0.5}, {1.0, 0.0, 0.5}, {0.0, 1.0, 0.5}, {1.0, 1.0, 0.5}
- });
+ const Point p0{0.5*scale, 0.5*scale, -1.0*scale};
+ const Point q0{0.5*scale, 0.5*scale, 1.0*scale};
- std::vector<Dune::FieldVector<double, dimworld>> triCorners({
- {-0.1, -0.1, 0.3}, {1.1, -0.1, 0.3}, {0.5, 2.0, 0.8}
- });
+ testSegTriangle(a, b, c, p0, q0);
+ testSegTriangle(a, b, c, q0, p0);
- Geometry3D cube(Dune::GeometryTypes::cube(dimworld), cubeCorners);
- Geometry2D quad(Dune::GeometryTypes::cube(dimworld-1), quadCorners);
- Geometry2D tri(Dune::GeometryTypes::simplex(dimworld-1), triCorners);
+ const Point p1 = a;
+ const Point q1 = b;
- if (Test::intersection(cube, quad, intersectionPolygon))
- {
- const auto triangulation = Dumux::triangulate<2, dimworld>(intersectionPolygon);
- Dumux::writeVTKPolyDataTriangle(triangulation, "quad_intersections");
- if (triangulation.size() != 4)
- DUNE_THROW(Dune::InvalidStateException, "Found " << triangulation.size() << " instead of 4 intersections!");
- }
- else
- DUNE_THROW(Dune::InvalidStateException, "No intersections found!");
+ testSegTriangle(a, b, c, p1, q1, false);
+ testSegTriangle(a, b, c, q1, p1, false);
- if (Test::intersection(cube, tri, intersectionPolygon))
- {
- const auto triangulation = Dumux::triangulate<2, dimworld>(intersectionPolygon);
- Dumux::writeVTKPolyDataTriangle(triangulation, "tri_intersections");
- if (triangulation.size() != 6)
- DUNE_THROW(Dune::InvalidStateException, "Found " << triangulation.size() << " instead of 6 intersections!");
- }
- else
- DUNE_THROW(Dune::InvalidStateException, "No intersections found!");
-
- // Test a tricky situation intersecting a hexagon with two quadrilaterals.
- // One face of the hexagon is in-plane with the two quadrilaterals, which
- // together enclose the entire hexagon face. Thus, the sum of the areas of
- // the two intersections must be equal to that of the respective hexagon face.
- // The separating edge between the two quadrilaterals passes very close by one
- // of the hex face corners. This situation was taken from an actual mortar
- // application which exposed a bug in the 1d-2d algorithm in 3d space.
- std::vector<Dune::FieldVector<double, dimworld>> hexCorners({
- {0.8619712261141845, 1.206565697102448, 0.4152254463380627},
- {0.8715547951379876, 1.175481516085758, 0.2599429674402531},
- {0.7046736652858085, 1.209207158782034, 0.4447521639960824},
- {0.7307137896935333, 1.171606500732276, 0.2574814772189353},
- {0.8926414294658755, 1.0, 0.4012781039132332},
- {0.8316315415883838, 1.0, 0.2581023823497354},
- {0.7349272235105835, 1.0, 0.4613391568883559},
- {0.7479805294570721, 1.0, 0.3155931875126768}
- });
-
- std::vector<Dune::FieldVector<double, dimworld>> quad1Corners({
- {1.000000000000000000000000000000, 1.0, 1.000000000000000000000000000000},
- {1.000000000000000000000000000000, 1.0, 0.327819111366782434124900191819},
- {0.303996171086607036571081152942, 1.0, 1.000000000000000000000000000000},
- {0.297170743844278550938042826601, 1.0, 0.293625720570556747457402479995}
- });
-
- std::vector<Dune::FieldVector<double, dimworld>> quad2Corners({
- {1.000000000000000000000000000000, 1.0, 0.327819111366782434124900191819},
- {1.000000000000000000000000000000, 1.0, 0.000000000000000000000000000000},
- {0.297170743844278550938042826601, 1.0, 0.293625720570556747457402479995},
- {0.325413399309280815252520824288, 1.0, 0.000000000000000000000000000000}
- });
-
- Geometry3D hex(Dune::GeometryTypes::cube(dimworld), hexCorners);
- Geometry2D quad1(Dune::GeometryTypes::cube(dimworld-1), quad1Corners);
- Geometry2D quad2(Dune::GeometryTypes::cube(dimworld-1), quad2Corners);
-
- // lambda to compute the area of a triangulated intersection
- auto computeArea = [] (const auto& triangulation)
- {
- double a = 0.0;
- for (const auto& t : triangulation)
- a += Geometry2D(Dune::GeometryTypes::simplex(dimworld-1),
- std::vector<Point>(t.begin(), t.end())).volume();
- return a;
- };
-
- double area = 0.0;
- if (Test::intersection(hex, quad1, intersectionPolygon))
- {
- const auto triangulation = Dumux::triangulate<2, dimworld>(intersectionPolygon);
- Dumux::writeVTKPolyDataTriangle(triangulation, "quad1_intersections");
- if (triangulation.size() != 5)
- DUNE_THROW(Dune::InvalidStateException, "Found " << triangulation.size() << " instead of 5 intersections!");
+ const Point p2{0.0*scale, 0.0*scale, 0.0*scale};
+ const Point q2{0.0*scale, 0.0*scale, 0.2*scale};
- area += computeArea(triangulation);
- }
+ testSegTriangle(a, b, c, p2, q2);
- if (Test::intersection(hex, quad2, intersectionPolygon))
- {
- const auto triangulation = Dumux::triangulate<2, dimworld>(intersectionPolygon);
- Dumux::writeVTKPolyDataTriangle(triangulation, "quad2_intersections");
- if (triangulation.size() != 1)
- DUNE_THROW(Dune::InvalidStateException, "Found " << triangulation.size() << " instead of 1 intersections!");
+ const auto cube = makeUnitCube(scale);
+ const auto quad = makeHorizontalQuad(scale, 0.5);
+ using Geometry3D = decltype(cube);
+ using Geometry2D = decltype(quad);
+ using Test = Dumux::GeometryIntersection<Geometry3D, Geometry2D>;
+ typename Test::Intersection intersectionPolygon;
- area += computeArea(triangulation);
+ if (Test::intersection(cube, quad, intersectionPolygon))
+ {
+ const auto triangulation = Dumux::triangulate<2, dimworld>(intersectionPolygon);
+ Dumux::writeVTKPolyDataTriangle(triangulation, "quad_intersections");
+ if (triangulation.size() != 4)
+ DUNE_THROW(Dune::InvalidStateException, "Found " << triangulation.size() << " instead of 4 intersections!");
+ }
+ else
+ DUNE_THROW(Dune::InvalidStateException, "No intersections found!");
+
+ // quad that is slightly above the cube
+ const auto close_quad = makeHorizontalQuad(scale, 1.0+1e-6);
+ if (Test::intersection(cube, close_quad, intersectionPolygon))
+ DUNE_THROW(Dune::InvalidStateException, "Found unexpected intersection with close quad");
+
+ const auto tri = makeInclinedTriangle(scale);
+ if (Test::intersection(cube, tri, intersectionPolygon))
+ {
+ const auto triangulation = Dumux::triangulate<2, dimworld>(intersectionPolygon);
+ Dumux::writeVTKPolyDataTriangle(triangulation, "tri_intersections");
+ if (triangulation.size() != 6)
+ DUNE_THROW(Dune::InvalidStateException, "Found " << triangulation.size() << " instead of 6 intersections!");
+ }
+ else
+ DUNE_THROW(Dune::InvalidStateException, "No intersections found!");
+
+ // Test a tricky situation intersecting a hexagon with two quadrilaterals.
+ // One face of the hexagon is in-plane with the two quadrilaterals, which
+ // together enclose the entire hexagon face. Thus, the sum of the areas of
+ // the two intersections must be equal to that of the respective hexagon face.
+ // The separating edge between the two quadrilaterals passes very close by one
+ // of the hex face corners. This situation was taken from an actual mortar
+ // application which exposed a bug in the 1d-2d algorithm in 3d space.
+ std::vector<Dune::FieldVector<double, dimworld>> hexCorners({
+ {0.8619712261141845*scale, 1.206565697102448*scale, 0.4152254463380627*scale},
+ {0.8715547951379876*scale, 1.175481516085758*scale, 0.2599429674402531*scale},
+ {0.7046736652858085*scale, 1.209207158782034*scale, 0.4447521639960824*scale},
+ {0.7307137896935333*scale, 1.171606500732276*scale, 0.2574814772189353*scale},
+ {0.8926414294658755*scale, 1.0*scale, 0.4012781039132332*scale},
+ {0.8316315415883838*scale, 1.0*scale, 0.2581023823497354*scale},
+ {0.7349272235105835*scale, 1.0*scale, 0.4613391568883559*scale},
+ {0.7479805294570721*scale, 1.0*scale, 0.3155931875126768*scale}
+ });
+
+ std::vector<Dune::FieldVector<double, dimworld>> quad1Corners({
+ {1.000000000000000000000000000000*scale, 1.0*scale, 1.000000000000000000000000000000*scale},
+ {1.000000000000000000000000000000*scale, 1.0*scale, 0.327819111366782434124900191819*scale},
+ {0.303996171086607036571081152942*scale, 1.0*scale, 1.000000000000000000000000000000*scale},
+ {0.297170743844278550938042826601*scale, 1.0*scale, 0.293625720570556747457402479995*scale}
+ });
+
+ std::vector<Dune::FieldVector<double, dimworld>> quad2Corners({
+ {1.000000000000000000000000000000*scale, 1.0*scale, 0.327819111366782434124900191819*scale},
+ {1.000000000000000000000000000000*scale, 1.0*scale, 0.000000000000000000000000000000*scale},
+ {0.297170743844278550938042826601*scale, 1.0*scale, 0.293625720570556747457402479995*scale},
+ {0.325413399309280815252520824288*scale, 1.0*scale, 0.000000000000000000000000000000*scale}
+ });
+
+ Geometry3D hex(Dune::GeometryTypes::cube(dimworld), hexCorners);
+ Geometry2D quad1(Dune::GeometryTypes::cube(dimworld-1), quad1Corners);
+ Geometry2D quad2(Dune::GeometryTypes::cube(dimworld-1), quad2Corners);
+
+ // lambda to compute the area of a triangulated intersection
+ auto computeArea = [] (const auto& triangulation)
+ {
+ double a = 0.0;
+ for (const auto& t : triangulation)
+ a += Geometry2D(Dune::GeometryTypes::simplex(dimworld-1),
+ std::vector<Point>(t.begin(), t.end())).volume();
+ return a;
+ };
+
+ double area = 0.0;
+ if (Test::intersection(hex, quad1, intersectionPolygon))
+ {
+ const auto triangulation = Dumux::triangulate<2, dimworld>(intersectionPolygon);
+ Dumux::writeVTKPolyDataTriangle(triangulation, "quad1_intersections");
+ if (triangulation.size() != 5)
+ DUNE_THROW(Dune::InvalidStateException, "Found " << triangulation.size() << " instead of 5 intersections!");
+
+ area += computeArea(triangulation);
+ }
+
+ if (Test::intersection(hex, quad2, intersectionPolygon))
+ {
+ const auto triangulation = Dumux::triangulate<2, dimworld>(intersectionPolygon);
+ Dumux::writeVTKPolyDataTriangle(triangulation, "quad2_intersections");
+ if (triangulation.size() != 1)
+ DUNE_THROW(Dune::InvalidStateException, "Found " << triangulation.size() << " instead of 1 intersections!");
+
+ area += computeArea(triangulation);
+ }
+
+ // compute area of the intersecting hexagon face
+ std::vector<Dune::FieldVector<double, dimworld>> hexFaceCorners({
+ {0.8926414294658755*scale, 1.0*scale, 0.4012781039132332*scale},
+ {0.8316315415883838*scale, 1.0*scale, 0.2581023823497354*scale},
+ {0.7349272235105835*scale, 1.0*scale, 0.4613391568883559*scale},
+ {0.7479805294570721*scale, 1.0*scale, 0.3155931875126768*scale}
+ });
+ Geometry2D hexFace(Dune::GeometryTypes::cube(dimworld-1), hexFaceCorners);
+ const auto hexFaceArea = hexFace.volume();
+
+ using std::abs;
+ const auto areaMismatch = abs(area-hexFaceArea);
+ std::cout << "Area mismatch is: " << areaMismatch << std::endl;
+ if (areaMismatch > hexFaceArea*1e-7)
+ DUNE_THROW(Dune::InvalidStateException, "Intersection area mismatch too large!");
}
- // compute area of the intersecting hexagon face
- std::vector<Dune::FieldVector<double, dimworld>> hexFaceCorners({
- {0.8926414294658755, 1.0, 0.4012781039132332},
- {0.8316315415883838, 1.0, 0.2581023823497354},
- {0.7349272235105835, 1.0, 0.4613391568883559},
- {0.7479805294570721, 1.0, 0.3155931875126768}
- });
- Geometry2D hexFace(Dune::GeometryTypes::cube(dimworld-1), hexFaceCorners);
- const auto hexFaceArea = hexFace.volume();
-
- using std::abs;
- const auto areaMismatch = abs(area-hexFaceArea);
- std::cout << "Area mismatch is: " << areaMismatch << std::endl;
- if (areaMismatch > 1e-14)
- DUNE_THROW(Dune::InvalidStateException, "Intersection area mismatch too large!");
-
std::cout << "All tests passed!" << std::endl;
return 0;