Commit 5fb834c0 authored by Timo Koch's avatar Timo Koch
Browse files

Merge branch 'feature/1d3d-intersection' into 'master'

Extend intersection tests for 1d3d

Closes #987

See merge request !2453
parents e954822c 8a7ab7e4
......@@ -773,6 +773,12 @@ private:
facetCorners = {{2, 0, 6, 4}, {1, 3, 5, 7}, {0, 1, 4, 5},
{3, 2, 7, 6}, {1, 0, 3, 2}, {4, 5, 6, 7}};
break;
case 6: // prism
facetCorners = {{0, 2, 1}, {3, 4, 5}, {0, 1, 3, 4}, {0, 3, 2, 5}, {1, 2, 4, 5}};
break;
case 5: // pyramid
facetCorners = {{0, 2, 1, 3}, {0, 1, 4}, {1, 3, 4}, {2, 4, 3}, {0, 4, 2}};
break;
case 4: // tetrahedron
facetCorners = {{1, 0, 2}, {0, 1, 3}, {0, 3, 2}, {1, 2, 3}};
break;
......
......@@ -2,127 +2,180 @@
#include <iostream>
#include <algorithm>
#include <initializer_list>
#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/geometry/geometryintersection.hh>
#include <dumux/common/math.hh>
#include "transformation.hh"
#include <dumux/geometry/geometryintersection.hh>
#ifndef DOXYGEN
template<int dimworld = 3>
Dune::MultiLinearGeometry<double, 1, dimworld>
makeLine(std::initializer_list<Dune::FieldVector<double, dimworld>>&& c)
Dune::MultiLinearGeometry<double, 1, 3>
makeLine(std::initializer_list<Dune::FieldVector<double, 3>>&& c)
{
return {Dune::GeometryTypes::line, c};
}
template<int dimworld = 3>
bool testIntersection(const Dune::MultiLinearGeometry<double, dimworld, dimworld>& polyhedron,
const Dune::MultiLinearGeometry<double, 1, dimworld>& line,
bool foundExpected = true)
bool testIntersection(const Dune::MultiLinearGeometry<double, 3, 3>& polyhedron,
const Dune::MultiLinearGeometry<double, 1, 3>& line,
bool foundExpected, bool verbose)
{
using Test = Dumux::GeometryIntersection<Dune::MultiLinearGeometry<double,dimworld,dimworld>,
Dune::MultiLinearGeometry<double,1,dimworld> >;
using Test = Dumux::GeometryIntersection<Dune::MultiLinearGeometry<double, 3, 3>,
Dune::MultiLinearGeometry<double, 1, 3> >;
typename Test::Intersection intersection;
bool found = Test::intersection(polyhedron, line, intersection);
if (!found && foundExpected)
std::cerr << "Failed detecting intersection with " << line.corner(0) << " " << line.corner(1) << std::endl;
else if (found && foundExpected)
std::cout << "Found intersection with " << line.corner(0) << " " << line.corner(1) << std::endl;
{
std::cerr << " Failed detecting intersection with " << line.corner(0) << " " << line.corner(1) << std::endl;
}
else if (found && !foundExpected)
std::cerr << "Found false positive: intersection with " << line.corner(0) << " " << line.corner(1) << std::endl;
else if (!found && !foundExpected)
std::cout << "No intersection with " << line.corner(0) << " " << line.corner(1) << std::endl;
{
std::cerr << " Found false positive: intersection with " << line.corner(0) << " " << line.corner(1) << std::endl;
}
if (verbose)
{
if (found && foundExpected)
std::cout << " Found intersection with " << line.corner(0) << " " << line.corner(1) << std::endl;
else if (!found && !foundExpected)
std::cout << " No intersection with " << line.corner(0) << " " << line.corner(1) << std::endl;
}
return (found == foundExpected);
}
#endif
int main (int argc, char *argv[])
namespace Dumux {
void testEdgesAndDiagonals(std::vector<bool>& returns,
const Dune::MultiLinearGeometry<double, 3, 3>& polyhedron,
bool verbose)
{
// maybe initialize mpi
Dune::MPIHelper::instance(argc, argv);
for (int i = 0; i < polyhedron.corners(); ++i)
{
for (int j = 0; j < polyhedron.corners(); ++j)
{
const auto& ci = polyhedron.corner(i);
const auto& cj = polyhedron.corner(j);
returns.push_back(testIntersection(polyhedron, makeLine({ci, cj}), true, verbose));
returns.push_back(testIntersection(polyhedron, makeLine({ci, cj + (ci-cj)}), true, verbose));
returns.push_back(testIntersection(polyhedron, makeLine({ci - (ci-cj), cj}), true, verbose));
}
}
}
constexpr int dimworld = 3;
constexpr int dim = 3;
template<class Transformation>
void runIntersectionTest(std::vector<bool>& returns, const Transformation& transform, bool verbose)
{
using Points = std::vector<Dune::FieldVector<double, 3>>;
using Geo = Dune::MultiLinearGeometry<double, 3, 3>;
// test cube-line intersections
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}
});
// test tetrahedron-line intersections
if (verbose) std::cout << "\n -- Test tetrahedron-line intersections" << std::endl;
auto cornersTetrahedron = 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(cornersTetrahedron.begin(), cornersTetrahedron.end(), cornersTetrahedron.begin(),
[&](const auto& p) { return transform(p); });
const auto tetrahedron = Geo(Dune::GeometryTypes::tetrahedron, cornersTetrahedron);
// test all edges and diagonals in both orientations
testEdgesAndDiagonals(returns, tetrahedron, verbose);
returns.push_back(testIntersection(tetrahedron, makeLine({{transform({0.25, 0.25, 0.0})}, {transform({0.25, 0.25, 1.0})}}), true, verbose));
returns.push_back(testIntersection(tetrahedron, makeLine({{transform({-1.0, 0.25, 0.5})}, {transform({1.0, 0.25, 0.5})}}), true, verbose));
returns.push_back(testIntersection(tetrahedron, makeLine({{transform({1.0, 1.0, 1.0})}, {transform({-1.0, -1.0, -1.0})}}), true, verbose));
returns.push_back(testIntersection(tetrahedron, makeLine({{transform({1.5, 0.0, 0.5})}, {transform({0.0, 1.5, 0.5})}}), false, verbose));
returns.push_back(testIntersection(tetrahedron, makeLine({{transform({0.0, 0.0, 0.0})}, {transform({0.0, 0.0, -1.0})}}), false, verbose));
returns.push_back(testIntersection(tetrahedron, makeLine({{transform({1.0, 1.0, 0.0})}, {transform({0.0, 0.0, 2.0})}}), false, verbose));
returns.push_back(testIntersection(tetrahedron, makeLine({{transform({1.0, 0.0, 0.1})}, {transform({0.0, 1.0, 0.1})}}), false, verbose));
returns.push_back(testIntersection(tetrahedron, makeLine({{transform({0.0, 0.0, -0.1})}, {transform({1.0, 1.0, -0.1})}}), false, verbose));
// lines with a small offset
returns.push_back(testIntersection(tetrahedron, makeLine({{transform({0.0, -1e-5, 0.0})}, {transform({1.0, 0.0, 0.0})}}), false, verbose));
returns.push_back(testIntersection(tetrahedron, makeLine({{transform({0.0, -1e-5, 0.0})}, {transform({1.0, -1e-5, 0.0})}}), false, verbose));
returns.push_back(testIntersection(tetrahedron, makeLine({{transform({0.0, 0.0, 0.0})}, {transform({1.0, -1e-5, 0.0})}}), false, verbose));
// test hexahedron-line intersections
if (verbose) std::cout << "\n -- Test hexahedron-line intersections" << std::endl;
auto cornersHexahedron = 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(cornersHexahedron.begin(), cornersHexahedron.end(), cornersHexahedron.begin(),
[&](const auto& p) { return transform(p); });
auto hexahedron = Geo(Dune::GeometryTypes::hexahedron, cornersHexahedron);
// test all edges and diagonals in both orientations
testEdgesAndDiagonals(returns, hexahedron, verbose);
returns.push_back(testIntersection(hexahedron, makeLine({{transform({0.0, 0.0, 2.0})}, {transform({1.0, 1.0, 2.0})}}), false, verbose));
returns.push_back(testIntersection(hexahedron, makeLine({{transform({0.0, 0.0, 1.1})}, {transform({1.0, 1.0, 1.1})}}), false, verbose));
returns.push_back(testIntersection(hexahedron, makeLine({{transform({1.1, 1.1, 0.0})}, {transform({1.1, 1.1, 1.0})}}), false, verbose));
returns.push_back(testIntersection(hexahedron, makeLine({{transform({1.1, 0.0, 0.0})}, {transform({1.1, 1.0, 1.0})}}), false, verbose));
returns.push_back(testIntersection(hexahedron, makeLine({{transform({0.0, -0.1, 0.0})}, {transform({1.0, -0.1, 0.0})}}), false, verbose));
// test pyramid-line intersections
if (verbose) std::cout << "\n -- Test pyramid-line 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(),
[&](const auto& p) { return transform(p); });
auto pyramid = Geo(Dune::GeometryTypes::pyramid, cornersPyramid);
// test all edges and diagonals in both orientations
testEdgesAndDiagonals(returns, pyramid, verbose);
returns.push_back(testIntersection(pyramid, makeLine({{transform({0.0, 0.0, 1.0})}, {transform({1.0, 0.0, 1.0})}}), false, verbose));
returns.push_back(testIntersection(pyramid, makeLine({{transform({0.0, 0.0, 1.0})}, {transform({0.0, 1.0, 1.0})}}), false, verbose));
returns.push_back(testIntersection(pyramid, makeLine({{transform({0.0, 0.0, -0.1})}, {transform({1.0, 1.0, -0.1})}}), false, verbose));
returns.push_back(testIntersection(pyramid, makeLine({{transform({0.0, 1.1, 0.0})}, {transform({1.0, 1.1, 0.0})}}), false, verbose));
returns.push_back(testIntersection(pyramid, makeLine({{transform({0.4, 0.0, 1.0})}, {transform({0.4, 1.0, 1.0})}}), false, verbose));
// test prism-line intersections
if (verbose) std::cout << "\n -- Test prism-line 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(),
[&](const auto& p) { return transform(p); });
auto prism = Geo(Dune::GeometryTypes::prism, cornersPrism);
// test all edges and diagonals in both orientations
testEdgesAndDiagonals(returns, prism, verbose);
returns.push_back(testIntersection(prism, makeLine({{transform({1.0, 1.0, 0.0})}, {transform({1.0, 1.0, 1.0})}}), false, verbose));
returns.push_back(testIntersection(prism, makeLine({{transform({2.0, 0.0, 0.5})}, {transform({0.0, 2.0, 0.5})}}), false, verbose));
returns.push_back(testIntersection(prism, makeLine({{transform({1.1, 0.0, 0.0})}, {transform({1.1, 0.0, 1.0})}}), false, verbose));
returns.push_back(testIntersection(prism, makeLine({{transform({-0.1, 0.0, 1.0})}, {transform({-0.1, 1.0, 1.0})}}), false, verbose));
returns.push_back(testIntersection(prism, makeLine({{transform({1.0, 0.0, 1.1})}, {transform({0.0, 1.0, 1.1})}}), false, verbose));
}
} // end namespace Dumux
Dune::MultiLinearGeometry<double, dim, dimworld>
cube(Dune::GeometryTypes::cube(dimworld), cubeCorners);
int main (int argc, char *argv[])
{
using namespace Dumux;
// collect returns to determine exit code
std::vector<bool> returns;
constexpr bool verbose = false;
// the tests
returns.push_back(testIntersection(cube, makeLine({{0.0, 0.0, 0.0}, {1.0, 0.0, 0.0}})));
returns.push_back(testIntersection(cube, makeLine({{0.0, 0.0, 0.0}, {0.0, 1.0, 0.0}})));
returns.push_back(testIntersection(cube, makeLine({{0.0, 0.0, 0.0}, {0.0, 0.0, 1.0}})));
returns.push_back(testIntersection(cube, makeLine({{1.0, 0.0, 0.0}, {0.0, 0.0, 0.0}})));
returns.push_back(testIntersection(cube, makeLine({{1.0, 0.0, 0.0}, {1.0, 0.0, 1.0}})));
returns.push_back(testIntersection(cube, makeLine({{1.0, 0.0, 0.0}, {1.0, 1.0, 0.0}})));
returns.push_back(testIntersection(cube, makeLine({{0.0, 1.0, 0.0}, {0.0, 0.0, 0.0}})));
returns.push_back(testIntersection(cube, makeLine({{0.0, 1.0, 0.0}, {0.0, 1.0, 1.0}})));
returns.push_back(testIntersection(cube, makeLine({{0.0, 1.0, 0.0}, {1.0, 1.0, 0.0}})));
returns.push_back(testIntersection(cube, makeLine({{1.0, 1.0, 0.0}, {1.0, 1.0, 1.0}})));
returns.push_back(testIntersection(cube, makeLine({{1.0, 1.0, 0.0}, {0.0, 1.0, 0.0}})));
returns.push_back(testIntersection(cube, makeLine({{1.0, 1.0, 0.0}, {1.0, 0.0, 0.0}})));
returns.push_back(testIntersection(cube, makeLine({{0.0, 0.0, 1.0}, {1.0, 0.0, 1.0}})));
returns.push_back(testIntersection(cube, makeLine({{0.0, 0.0, 1.0}, {0.0, 1.0, 1.0}})));
returns.push_back(testIntersection(cube, makeLine({{0.0, 0.0, 1.0}, {0.0, 0.0, 0.0}})));
returns.push_back(testIntersection(cube, makeLine({{1.0, 0.0, 1.0}, {0.0, 0.0, 1.0}})));
returns.push_back(testIntersection(cube, makeLine({{1.0, 0.0, 1.0}, {1.0, 0.0, 0.0}})));
returns.push_back(testIntersection(cube, makeLine({{1.0, 0.0, 1.0}, {1.0, 1.0, 1.0}})));
returns.push_back(testIntersection(cube, makeLine({{0.0, 1.0, 1.0}, {0.0, 0.0, 1.0}})));
returns.push_back(testIntersection(cube, makeLine({{0.0, 1.0, 1.0}, {0.0, 1.0, 0.0}})));
returns.push_back(testIntersection(cube, makeLine({{0.0, 1.0, 1.0}, {1.0, 1.0, 1.0}})));
returns.push_back(testIntersection(cube, makeLine({{1.0, 1.0, 1.0}, {1.0, 1.0, 0.0}})));
returns.push_back(testIntersection(cube, makeLine({{1.0, 1.0, 1.0}, {0.0, 1.0, 1.0}})));
returns.push_back(testIntersection(cube, makeLine({{1.0, 1.0, 1.0}, {1.0, 0.0, 1.0}})));
returns.push_back(testIntersection(cube, makeLine({{0.0, 0.0, 0.0}, {1.0, 1.0, 1.0}})));
returns.push_back(testIntersection(cube, makeLine({{0.5, 0.5, 0.5}, {0.5, 0.5, -2.0}})));
returns.push_back(testIntersection(cube, makeLine({{0.5, 0.5, 0.0}, {0.5, 0.5, -2.0}}), false));
returns.push_back(testIntersection(cube, makeLine({{1.0, 1.0, 1.0}, {2.0, 2.0, 2.0}}), false));
// test tetrahedron-line intersections
std::vector<Dune::FieldVector<double, dimworld>> tetCorners({
{0.0, 0.0, 0.0}, {1.0, 0.0, 0.0}, {0.0, 1.0, 0.0}, {0.0, 0.0, 1.0}
});
Dune::MultiLinearGeometry<double, dim, dimworld>
tet(Dune::GeometryTypes::simplex(dimworld), tetCorners);
// the tests
returns.push_back(testIntersection(tet, makeLine({{0.0, 0.0, 0.0}, {1.0, 1.0, 1.0}})));
returns.push_back(testIntersection(tet, makeLine({{0.25, 0.25, 0.0}, {0.25, 0.25, 1.0}})));
returns.push_back(testIntersection(tet, makeLine({{-1.0, 0.25, 0.5}, {1.0, 0.25, 0.5}})));
returns.push_back(testIntersection(tet, makeLine({{1.0, 1.0, 1.0}, {-1.0, -1.0, -1.0}})));
returns.push_back(testIntersection(tet, makeLine({{1.5, 0.0, 0.5}, {0.0, 1.5, 0.5}}), false));
returns.push_back(testIntersection(tet, makeLine({{0.0, 0.0, 0.0}, {0.0, 0.0, -1.0}}), false));
using Vec = Dune::FieldVector<double, 3>;
for (const double scaling : {1.0, 1e3, 1e12, 1e-12})
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, make3DTransformation<double>(scaling, Vec(translation), rotAxis, angle, true), verbose);
// 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;
}
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment