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;