From 7bf9e339d6f680d7aaf9fddc812a6dc969e56d80 Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Sat, 4 Apr 2020 19:22:33 +0200
Subject: [PATCH] [test][0d3d] Run intersection test with various
transformations
---
.../common/geometry/test_0d3d_intersection.cc | 98 +++++++++++++------
1 file changed, 69 insertions(+), 29 deletions(-)
diff --git a/test/common/geometry/test_0d3d_intersection.cc b/test/common/geometry/test_0d3d_intersection.cc
index 5bff2e0698..16c4899888 100644
--- a/test/common/geometry/test_0d3d_intersection.cc
+++ b/test/common/geometry/test_0d3d_intersection.cc
@@ -3,7 +3,9 @@
#include <iostream>
#include <algorithm>
#include <functional>
+#include <type_traits>
+#include <dune/common/classname.hh>
#include <dune/common/exceptions.hh>
#include <dune/common/fvector.hh>
#include <dune/geometry/type.hh>
@@ -16,28 +18,42 @@
template<class Geometry>
bool testIntersection(const Geometry& geo,
const typename Geometry::GlobalCoordinate& p,
- bool foundExpected)
+ bool foundExpected, bool verbose)
{
bool found = Dumux::intersectsPointGeometry(p, geo);
if (!found && foundExpected)
- std::cerr << " Failed detecting intersection with " << p << std::endl;
- else if (found && foundExpected)
- std::cout << " Found intersection with " << p << std::endl;
+ {
+ std::cerr << " Failed detecting intersection of " << geo.type();
+ for (int i = 0; i < geo.corners(); ++i)
+ std::cerr << " (" << geo.corner(i) << ")";
+ std::cerr << " with point: " << p << std::endl;
+ }
else if (found && !foundExpected)
- std::cerr << " Found false positive: intersection with " << p << std::endl;
- else if (!found && !foundExpected)
- std::cout << " No intersection with " << p << std::endl;
+ {
+ std::cerr << " Found false positive: intersection of " << geo.type();
+ for (int i = 0; i < geo.corners(); ++i)
+ std::cerr << " (" << geo.corner(i) << ")";
+ std::cerr << " with point: " << p << std::endl;
+ }
+ if (verbose)
+ {
+ if (found && foundExpected)
+ std::cout << " Found intersection with " << p << std::endl;
+ else if (!found && !foundExpected)
+ std::cout << " No intersection with " << p << std::endl;
+ }
return (found == foundExpected);
}
-template<class ctype, class Transformation>
-void runTest(std::vector<bool>& returns, const Transformation& transform)
+template<class Transformation>
+void runIntersectionTest(std::vector<bool>& returns, const Transformation& transform, bool verbose)
{
+ using ctype = typename std::decay_t<decltype(transform({0.0}))>::value_type;
using Geo = Dune::MultiLinearGeometry<ctype, 3, 3>;
using Points = std::vector<typename Geo::GlobalCoordinate>;
// test tetrahedron-point intersections
- std::cout << "\n -- Test tetrahedron-point intersections" << std::endl;
+ if (verbose) std::cout << "\n -- Test tetrahedron-point intersections" << std::endl;
auto cornersTet = Points({{0.0, 0.0, 0.0}, {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}});
std::transform(cornersTet.begin(), cornersTet.end(), cornersTet.begin(),
@@ -45,15 +61,15 @@ void runTest(std::vector<bool>& returns, const Transformation& transform)
const auto tetrahedron = Geo(Dune::GeometryTypes::tetrahedron, cornersTet);
for (const auto& corner : cornersTet)
- returns.push_back(testIntersection(tetrahedron, corner, true));
- returns.push_back(testIntersection(tetrahedron, transform({0.0, 0.0, 0.5}), true));
- returns.push_back(testIntersection(tetrahedron, transform({0.25, 0.25, 0.5}), true));
- returns.push_back(testIntersection(tetrahedron, transform({0.5, 0.5, 0.5}), false));
- returns.push_back(testIntersection(tetrahedron, transform({1.01, 0.0, 0.0}), false));
- returns.push_back(testIntersection(tetrahedron, transform({0.5, 0.0, 0.51}), false));
+ returns.push_back(testIntersection(tetrahedron, corner, true, verbose));
+ returns.push_back(testIntersection(tetrahedron, transform({0.0, 0.0, 0.5}), true, verbose));
+ returns.push_back(testIntersection(tetrahedron, transform({0.25, 0.25, 0.5}), true, verbose));
+ returns.push_back(testIntersection(tetrahedron, transform({0.5, 0.5, 0.5}), false, verbose));
+ returns.push_back(testIntersection(tetrahedron, transform({1.01, 0.0, 0.0}), false, verbose));
+ returns.push_back(testIntersection(tetrahedron, transform({0.5, 0.0, 0.51}), false, verbose));
// test hexahedron-point intersections
- std::cout << "\n -- Test hexahedron-point intersections" << std::endl;
+ if (verbose) std::cout << "\n -- Test hexahedron-point intersections" << std::endl;
auto cornersHex = Points({{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}});
@@ -62,10 +78,15 @@ void runTest(std::vector<bool>& returns, const Transformation& transform)
auto hexahedron = Geo(Dune::GeometryTypes::hexahedron, cornersHex);
for (const auto& corner : cornersHex)
- returns.push_back(testIntersection(hexahedron, corner, true));
+ returns.push_back(testIntersection(hexahedron, corner, true, verbose));
+ returns.push_back(testIntersection(hexahedron, transform({0.5, 0.5, 0.5}), true, verbose));
+ returns.push_back(testIntersection(hexahedron, transform({0.01, 0.01, 0.0001}), true, verbose));
+ returns.push_back(testIntersection(hexahedron, transform({1.01, 0.5, 0.5}), false, verbose));
+ returns.push_back(testIntersection(hexahedron, transform({2.0, 2.0, 2.0}), false, verbose));
+ returns.push_back(testIntersection(hexahedron, transform({-0.5, -0.0, -0.51}), false, verbose));
// test pyramid-point intersections
- std::cout << "\n -- Test pyramid-point intersections" << std::endl;
+ if (verbose) std::cout << "\n -- Test pyramid-point intersections" << std::endl;
auto cornersPyramid = Points({{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.5, 0.5, 1.0}});
@@ -74,10 +95,18 @@ void runTest(std::vector<bool>& returns, const Transformation& transform)
auto pyramid = Geo(Dune::GeometryTypes::pyramid, cornersPyramid);
for (const auto& corner : cornersPyramid)
- returns.push_back(testIntersection(pyramid, corner, true));
+ returns.push_back(testIntersection(pyramid, corner, true, verbose));
+ returns.push_back(testIntersection(pyramid, transform({0.5, 0.5, 0.0}), true, verbose));
+ returns.push_back(testIntersection(pyramid, transform({0.5, 0.5, 0.7}), true, verbose));
+ returns.push_back(testIntersection(pyramid, transform({0.5, 0.5, -0.0001}), false, verbose));
+ returns.push_back(testIntersection(pyramid, transform({0.25, 0.25, 0.5}), true, verbose));
+ returns.push_back(testIntersection(pyramid, transform({0.25, 0.75, 0.5}), true, verbose));
+ returns.push_back(testIntersection(pyramid, transform({0.75, 0.75, 0.5}), true, verbose));
+ returns.push_back(testIntersection(pyramid, transform({0.25, 0.25, 0.5001}), false, verbose));
+ returns.push_back(testIntersection(pyramid, transform({1.0, 1.0, 0.0001}), false, verbose));
// test prism-point intersections
- std::cout << "\n -- Test prism-point intersections" << std::endl;
+ if (verbose) std::cout << "\n -- Test prism-point intersections" << std::endl;
auto cornersPrism = Points({{0.0, 0.0, 0.0}, {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0},
{0.0, 0.0, 1.0}, {1.0, 0.0, 1.0}, {0.0, 1.0, 1.0}});
@@ -86,7 +115,13 @@ void runTest(std::vector<bool>& returns, const Transformation& transform)
auto prism = Geo(Dune::GeometryTypes::prism, cornersPrism);
for (const auto& corner : cornersPrism)
- returns.push_back(testIntersection(prism, corner, true));
+ returns.push_back(testIntersection(prism, corner, true, verbose));
+ returns.push_back(testIntersection(prism, transform({0.25, 0.0, 0.25}), true, verbose));
+ returns.push_back(testIntersection(prism, transform({0.0, 0.25, 0.25}), true, verbose));
+ returns.push_back(testIntersection(prism, transform({0.25, 0.25, 1.0}), true, verbose));
+ returns.push_back(testIntersection(prism, transform({0.25, 0.25, 0.0}), true, verbose));
+ returns.push_back(testIntersection(prism, transform({0.25, 0.25, -0.0001}), false, verbose));
+ returns.push_back(testIntersection(prism, transform({0.25, 0.25, 1.0001}), false, verbose));
}
template<class ctype>
@@ -95,8 +130,9 @@ auto createTransformation(const ctype scale,
const Dune::FieldVector<ctype, 3>& rotationAxis,
const ctype rotationAngle)
{
- std::cout << "\n\n Created transformation with"
- << " scaling: " << scale
+ std::cout << "Intersection test with transformation:"
+ << " ctype: " << Dune::className<ctype>()
+ << ", scaling: " << scale
<< ", translation: " << translate
<< ", rotationAxis: " << rotationAxis
<< ", rotationAngle: " << rotationAngle << std::endl;
@@ -104,25 +140,29 @@ auto createTransformation(const ctype scale,
const ctype cosAngle = std::cos(rotationAngle);
return [=](Dune::FieldVector<ctype, 3> p){
p *= scale;
- p += translate;
+ p.axpy(scale, translate);
auto tp = p;
tp *= cosAngle;
tp.axpy(sinAngle, Dumux::crossProduct(rotationAxis, p));
return tp.axpy((1.0-cosAngle)*(rotationAxis*p), rotationAxis);
};
}
+
#endif
int main (int argc, char *argv[]) try
{
// collect returns to determine exit code
std::vector<bool> returns;
+ constexpr bool verbose = false;
+ using Vec = Dune::FieldVector<double, 3>;
for (const double scaling : {1.0, 1e3, 1e12, 1e-12})
- {
- const auto transform = createTransformation(scaling, {0.0, 0.0, 0.0}, {1.0, 0.0, 0.0}, 0.0);
- runTest<double>(returns, transform);
- }
+ for (const double translation : {0.0, 1.0})
+ for (const double angle : {0.0, 0.2*M_PI, 0.5*M_PI, 0.567576567*M_PI, M_PI})
+ for (const auto& rotAxis : {Vec(std::sqrt(3.0)/3.0), Vec({std::sqrt(2.0)/2.0, std::sqrt(2.0)/2.0, 0.0})})
+ runIntersectionTest(returns,
+ createTransformation<double>(scaling, Vec(translation), rotAxis, angle), verbose);
// determine the exit code
if (std::any_of(returns.begin(), returns.end(), std::logical_not<bool>{}))
--
GitLab