diff --git a/test/common/geometry/test_0d3d_intersection.cc b/test/common/geometry/test_0d3d_intersection.cc
index f75a161a42493c66a9c6d1ae83f4b2ae08a633a3..c52c93cfe561d3a0a89ff810aba5ea1913d3224c 100644
--- a/test/common/geometry/test_0d3d_intersection.cc
+++ b/test/common/geometry/test_0d3d_intersection.cc
@@ -2,94 +2,139 @@
#include <iostream>
#include <algorithm>
+#include <functional>
#include <dune/common/exceptions.hh>
#include <dune/common/parallel/mpihelper.hh>
#include <dune/common/fvector.hh>
+#include <dune/geometry/type.hh>
#include <dune/geometry/multilineargeometry.hh>
+#include <dumux/common/math.hh>
#include <dumux/common/geometry/intersectspointgeometry.hh>
#ifndef DOXYGEN
-template<int dimworld = 3, class Geometry>
+template<class Geometry>
bool testIntersection(const Geometry& geo,
- const Dune::FieldVector<double, dimworld>& p,
- bool foundExpected = false)
+ const typename Geometry::GlobalCoordinate& p,
+ bool foundExpected)
{
bool found = Dumux::intersectsPointGeometry(p, geo);
if (!found && foundExpected)
- std::cerr << "Failed detecting intersection with " << p << std::endl;
+ std::cerr << " Failed detecting intersection with " << p << std::endl;
else if (found && foundExpected)
- std::cout << "Found intersection with " << p << std::endl;
+ std::cout << " Found intersection with " << p << std::endl;
else if (found && !foundExpected)
- std::cerr << "Found false positive: intersection with " << p << std::endl;
+ std::cerr << " Found false positive: intersection with " << p << std::endl;
else if (!found && !foundExpected)
- std::cout << "No intersection with " << p << std::endl;
+ 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)
+{
+ 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;
+
+ 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(),
+ [&transform](const auto& p) { return transform(p); });
+ 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));
+
+ // test hexahedron-point intersections
+ 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}});
+ std::transform(cornersHex.begin(), cornersHex.end(), cornersHex.begin(),
+ [&transform](const auto& p) { return transform(p); });
+ auto hexahedron = Geo(Dune::GeometryTypes::hexahedron, cornersHex);
+
+ for (const auto& corner : cornersHex)
+ returns.push_back(testIntersection(hexahedron, corner, true));
+
+ // test pyramid-point intersections
+ 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}});
+ std::transform(cornersPyramid.begin(), cornersPyramid.end(), cornersPyramid.begin(),
+ [&transform](const auto& p) { return transform(p); });
+ auto pyramid = Geo(Dune::GeometryTypes::pyramid, cornersPyramid);
+
+ for (const auto& corner : cornersPyramid)
+ returns.push_back(testIntersection(pyramid, corner, true));
+
+ // test prism-point intersections
+ 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}});
+ std::transform(cornersPrism.begin(), cornersPrism.end(), cornersPrism.begin(),
+ [&transform](const auto& p) { return transform(p); });
+ auto prism = Geo(Dune::GeometryTypes::prism, cornersPrism);
+
+ for (const auto& corner : cornersPrism)
+ returns.push_back(testIntersection(prism, corner, true));
+}
+
+template<class ctype>
+auto createTransformation(const ctype scale,
+ const Dune::FieldVector<ctype, 3>& translate,
+ const Dune::FieldVector<ctype, 3>& rotationAxis,
+ const ctype rotationAngle)
+{
+ std::cout << "\n\n Created transformation with"
+ << " scaling: " << scale
+ << ", translation: " << translate
+ << ", rotationAxis: " << rotationAxis
+ << ", rotationAngle: " << rotationAngle << std::endl;
+ const ctype sinAngle = std::sin(rotationAngle);
+ const ctype cosAngle = std::cos(rotationAngle);
+ return [=](Dune::FieldVector<ctype, 3> p){
+ p *= scale;
+ p += 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
{
// maybe initialize mpi
Dune::MPIHelper::instance(argc, argv);
-
- constexpr int dimworld = 3;
-
+
// collect returns to determine exit code
std::vector<bool> returns;
- for (auto scaling : {1.0, 1e3, 1e12, 1e-12})
+ for (const double scaling : {1.0, 1e3, 1e12, 1e-12})
{
- std::cout << "Test with scaling " << scaling << std::endl;
-
- using Points = std::vector<Dune::FieldVector<double, dimworld>>;
- using Geo = Dune::MultiLinearGeometry<double, 3, dimworld>;
-
- // test tetrahedron-point intersections
- std::cout << "test tetrahedron-point intersections" << std::endl;
-
- auto cornersTetrahedron = Points({{0.0, 0.0, 0.0}, {1.0*scaling, 0.0, 0.0}, {0.0, 1.0*scaling, 0.0}, {0.0, 0.0, 1.0*scaling}});
- auto tetrahedron = Geo(Dune::GeometryTypes::tetrahedron, cornersTetrahedron);
-
- returns.push_back(testIntersection(tetrahedron, {0.0, 0.0, 0.0}, true));
- returns.push_back(testIntersection(tetrahedron, {0.0, 0.0, 0.5*scaling}, true));
- returns.push_back(testIntersection(tetrahedron, {0.25*scaling, 0.25*scaling, 0.5*scaling}, true));
- returns.push_back(testIntersection(tetrahedron, {0.5*scaling, 0.5*scaling, 0.5*scaling}));
- returns.push_back(testIntersection(tetrahedron, {1.01*scaling, 0.0, 0.0}));
- returns.push_back(testIntersection(tetrahedron, {0.5*scaling, 0.0, 0.51*scaling}));
-
- // test hexahedron-point intersections
- std::cout << "test hexahedron-point intersections" << std::endl;
-
- auto cornersHexahedron = Points({{0.0, 0.0, 0.0}, {1.0*scaling, 0.0, 0.0}, {0.0, 1.0*scaling, 0.0}, {1.0*scaling, 1.0*scaling, 0.0}, {0.0, 0.0, 1.0*scaling}, {1.0*scaling, 0.0, 1.0*scaling}, {0.0, 1.0*scaling, 1.0*scaling}, {1.0*scaling, 1.0*scaling, 1.0*scaling}});
- auto hexahedron = Geo(Dune::GeometryTypes::hexahedron, cornersHexahedron);
-
- returns.push_back(testIntersection(hexahedron, {0.0, 0.0, 0.0}, true));
-
- // test pyramid-point intersections
- std::cout << "test pyramid-point intersections" << std::endl;
-
- auto cornersPyramid = Points({{0.0, 0.0, 0.0}, {1.0*scaling, 0.0, 0.0}, {0.0, 1.0*scaling, 0.0}, {1.0*scaling, 1.0*scaling, 0.0}, {0.5*scaling, 0.5*scaling, 1.0*scaling}});
- auto pyramid = Geo(Dune::GeometryTypes::pyramid, cornersPyramid);
-
- returns.push_back(testIntersection(pyramid, {0.0, 0.0, 0.0}, true));
-
- // test prism-point intersections
- std::cout << "test prism-point intersections" << std::endl;
-
- auto cornersPrism = Points({{0.0, 0.0, 0.0}, {1.0*scaling, 0.0, 0.0}, {0.0, 1.0*scaling, 0.0}, {0.0, 0.0, 1.0*scaling}, {1.0*scaling, 0.0, 1.0*scaling}, {0.0, 1.0*scaling, 1.0*scaling}});
- auto prism = Geo(Dune::GeometryTypes::prism, cornersPrism);
-
- returns.push_back(testIntersection(prism, {0.0, 0.0, 0.0}, true));
+ const auto transform = createTransformation(scaling, {0.0, 0.0, 0.0}, {1.0, 0.0, 0.0}, 0.0);
+ runTest<double>(returns, transform);
}
// determine the exit code
- if (std::any_of(returns.begin(), returns.end(), [](bool i){ return !i; }))
+ if (std::any_of(returns.begin(), returns.end(), std::logical_not<bool>{}))
return 1;
- std::cout << "All tests passed!" << std::endl;
+ std::cout << "\n++++++++++++++++++++++\n"
+ << "All tests passed!"
+ << "\n++++++++++++++++++++++" << std::endl;
return 0;
}