diff --git a/test/common/geometry/test_0d2d_intersection.cc b/test/common/geometry/test_0d2d_intersection.cc
index 25e5308fc1edeaaa5115340f69b9578bc3553f69..8c521e60065fd2036c7bafd27de3baff616343a9 100644
--- a/test/common/geometry/test_0d2d_intersection.cc
+++ b/test/common/geometry/test_0d2d_intersection.cc
@@ -12,204 +12,107 @@
#include <dune/geometry/multilineargeometry.hh>
#include <dumux/common/math.hh>
-#include <dumux/common/geometry/intersectspointgeometry.hh>
-#ifndef DOXYGEN
-template<class Geometry>
-bool testIntersection(const Geometry& geo,
- const typename Geometry::GlobalCoordinate& p,
- bool foundExpected, bool verbose)
-{
- bool found = Dumux::intersectsPointGeometry(p, geo);
- if (!found && foundExpected)
- {
- 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 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);
-}
+#include "transformation.hh"
+#include "test_intersection.hh"
+
+namespace Dumux {
template<int dimworld, 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, 2, dimworld>;
- using Points = std::vector<typename Geo::GlobalCoordinate>;
+ using Points2D = std::vector<Dune::FieldVector<ctype, 2>>;
+ using PointsDimWorld = std::vector<Dune::FieldVector<ctype, dimworld>>;
// test triangle-point intersections
if (verbose) std::cout << "\n -- Test triangle-point intersections" << std::endl;
- if constexpr (dimworld == 2)
- {
- auto cornersTri = Points ({{0.0, 0.0}, {1.0, 0.0}, {0.0, 1.0}});
- std::transform(cornersTri.begin(), cornersTri.end(), cornersTri.begin(),
- [&transform](const auto& p) { return transform(p); });
- auto triangle = Geo(Dune::GeometryTypes::triangle, cornersTri);
-
- for (const auto& corner : cornersTri)
- returns.push_back(testIntersection(triangle, corner, true, verbose));
-
- returns.push_back(testIntersection(triangle, transform({0.25, 0.25}), true, verbose));
- returns.push_back(testIntersection(triangle, transform({0.5, 0.5}), true, verbose));
- returns.push_back(testIntersection(triangle, transform({0.0, 0.5}), true, verbose));
- returns.push_back(testIntersection(triangle, transform({1.01, 0.0}), false, verbose));
- returns.push_back(testIntersection(triangle, transform({0.5, 0.51}), false, verbose));
- returns.push_back(testIntersection(triangle, transform({0.0, -0.01}), false, verbose));
- }
+ auto cornersTri2D = Points2D ({{0.0, 0.0}, {1.0, 0.0}, {0.0, 1.0}});
+ auto cornersTri = PointsDimWorld(3);
+ std::transform(cornersTri2D.begin(), cornersTri2D.end(), cornersTri.begin(),
+ [&](const auto& p) { return transform(p); });
+ auto triangle = Geo(Dune::GeometryTypes::triangle, cornersTri);
- if constexpr (dimworld == 3)
- {
- auto cornersTri = Points ({{0.0, 0.0, 2.0}, {1.0, 0.0, 2.0}, {0.0, 1.0, 2.0}});
- std::transform(cornersTri.begin(), cornersTri.end(), cornersTri.begin(),
- [&transform](const auto& p) { return transform(p); });
- auto triangle = Geo(Dune::GeometryTypes::triangle, cornersTri);
-
- for (const auto& corner : cornersTri)
- returns.push_back(testIntersection(triangle, corner, true, verbose));
-
- returns.push_back(testIntersection(triangle, transform({0.25, 0.25, 2.0}), true, verbose));
- returns.push_back(testIntersection(triangle, transform({0.5, 0.5, 2.0}), true, verbose));
- returns.push_back(testIntersection(triangle, transform({0.0, 0.5, 2.0}), true, verbose));
- returns.push_back(testIntersection(triangle, transform({1.01, 0.0, 2.0}), false, verbose));
- returns.push_back(testIntersection(triangle, transform({0.5, 0.51, 2.0}), false, verbose));
- returns.push_back(testIntersection(triangle, transform({0.0, -0.01, 2.0}), false, verbose));
- }
+ for (const auto& corner : cornersTri)
+ returns.push_back(testIntersection(triangle, corner, true, verbose));
+
+ returns.push_back(testIntersection(triangle, transform({0.25, 0.25}), true, verbose));
+ returns.push_back(testIntersection(triangle, transform({0.5, 0.5}), true, verbose));
+ returns.push_back(testIntersection(triangle, transform({0.0, 0.5}), true, verbose));
+ returns.push_back(testIntersection(triangle, transform({1.01, 0.0}), false, verbose));
+ returns.push_back(testIntersection(triangle, transform({0.5, 0.51}), false, verbose));
+ returns.push_back(testIntersection(triangle, transform({0.0, -0.01}), false, verbose));
// test quadrilateral-point intersections
if (verbose) std::cout << "\n -- Test quadrilateral-point intersections" << std::endl;
- if constexpr (dimworld == 2)
- {
- auto cornersQuad = Points ({{0.0, 0.0}, {1.0, 0.0}, {0.0, 1.0}, {1.0, 1.0}});
- std::transform(cornersQuad.begin(), cornersQuad.end(), cornersQuad.begin(),
- [&transform](const auto& p) { return transform(p); });
- auto quadrilateral = Geo(Dune::GeometryTypes::quadrilateral, cornersQuad);
-
- for (const auto& corner : cornersQuad)
- returns.push_back(testIntersection(quadrilateral, corner, true, verbose));
-
- returns.push_back(testIntersection(quadrilateral, transform({0.5, 0.5}), true, verbose));
- returns.push_back(testIntersection(quadrilateral, transform({0.5, 0.0}), true, verbose));
- returns.push_back(testIntersection(quadrilateral, transform({0.5, 1.0}), true, verbose));
- returns.push_back(testIntersection(quadrilateral, transform({1.01, 1.0}), false, verbose));
- returns.push_back(testIntersection(quadrilateral, transform({0.5, 1.01}), false, verbose));
- returns.push_back(testIntersection(quadrilateral, transform({0.0, -0.01}), false, verbose));
- }
+ auto cornersQuad2D = Points2D ({{0.0, 0.0}, {1.0, 0.0}, {0.0, 1.0}, {1.0, 1.0}});
+ auto cornersQuad = PointsDimWorld(4);
+ std::transform(cornersQuad2D.begin(), cornersQuad2D.end(), cornersQuad.begin(),
+ [&](const auto& p) { return transform(p); });
+ auto quadrilateral = Geo(Dune::GeometryTypes::quadrilateral, cornersQuad);
- if constexpr (dimworld == 3)
- {
- auto cornersQuad = Points ({{0.0, 0.0, 2.0}, {1.0, 0.0, 2.0}, {0.0, 1.0, 2.0}, {1.0, 1.0, 2.0}});
- std::transform(cornersQuad.begin(), cornersQuad.end(), cornersQuad.begin(),
- [&transform](const auto& p) { return transform(p); });
- auto quadrilateral = Geo(Dune::GeometryTypes::quadrilateral, cornersQuad);
-
- for (const auto& corner : cornersQuad)
- returns.push_back(testIntersection(quadrilateral, corner, true, verbose));
-
- returns.push_back(testIntersection(quadrilateral, transform({0.5, 0.5, 2.0}), true, verbose));
- returns.push_back(testIntersection(quadrilateral, transform({0.5, 0.0, 2.0}), true, verbose));
- returns.push_back(testIntersection(quadrilateral, transform({0.5, 1.0, 2.0}), true, verbose));
- returns.push_back(testIntersection(quadrilateral, transform({1.01, 1.0, 2.0}), false, verbose));
- returns.push_back(testIntersection(quadrilateral, transform({0.5, 1.01, 2.0}), false, verbose));
- returns.push_back(testIntersection(quadrilateral, transform({0.0, -0.01, 2.0}), false, verbose));
- }
-}
+ for (const auto& corner : cornersQuad)
+ returns.push_back(testIntersection(quadrilateral, corner, true, verbose));
+
+ returns.push_back(testIntersection(quadrilateral, transform({0.5, 0.5}), true, verbose));
+ returns.push_back(testIntersection(quadrilateral, transform({0.5, 0.0}), true, verbose));
+ returns.push_back(testIntersection(quadrilateral, transform({0.5, 1.0}), true, verbose));
+ returns.push_back(testIntersection(quadrilateral, transform({1.01, 1.0}), false, verbose));
+ returns.push_back(testIntersection(quadrilateral, transform({0.5, 1.01}), false, verbose));
+ returns.push_back(testIntersection(quadrilateral, transform({0.0, -0.01}), false, verbose));
-template<class ctype>
-auto create2DTransformation(const ctype scale,
- const Dune::FieldVector<ctype, 2>& translate,
- const ctype rotationAngle)
-{
- std::cout << "Intersection test with transformation:"
- << " ctype: " << Dune::className<ctype>()
- << ", scaling: " << scale
- << ", translation: " << translate
- << ", rotationAngle: " << rotationAngle << std::endl;
- // Rotation of a vector in two dimensions
- // See rotation matrix (2d) at https://en.wikipedia.org/wiki/Rotation_matrix
- const ctype sinAngle = std::sin(rotationAngle);
- const ctype cosAngle = std::cos(rotationAngle);
- return [=](Dune::FieldVector<ctype, 2> p){
- p *= scale;
- p.axpy(scale, translate);
- auto tp = p;
- tp[0] = p[0]*cosAngle-p[1]*sinAngle;
- tp[1] = p[0]*sinAngle+p[1]*cosAngle;
- return tp;
- };
}
-template<class ctype>
-auto create3DTransformation(const ctype scale,
- const Dune::FieldVector<ctype, 3>& translate,
- const Dune::FieldVector<ctype, 3>& rotationAxis,
- const ctype rotationAngle)
+template<class Transformation>
+void run3DIntersectionTest(std::vector<bool>& returns, const Transformation& transform, bool verbose)
{
- std::cout << "Intersection test with transformation:"
- << " ctype: " << Dune::className<ctype>()
- << ", scaling: " << scale
- << ", translation: " << translate
- << ", rotationAxis: " << rotationAxis
- << ", rotationAngle: " << rotationAngle << std::endl;
- // Rotation of a vector in three dimensions
- // See Rodrigues' rotation formular at https://en.wikipedia.org/wiki/Rodrigues%27_rotation_formula
- const ctype sinAngle = std::sin(rotationAngle);
- const ctype cosAngle = std::cos(rotationAngle);
- return [=](Dune::FieldVector<ctype, 3> p){
- p *= scale;
- 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);
+ // augment 2d vector by a third coordinate
+ using Point = std::decay_t<decltype(transform({0.0, 0.0, 0.0}))>;
+ const auto transform2D3D = [&](Dune::FieldVector<typename Point::value_type, 2> p2D){
+ return transform({p2D[0], p2D[1], 2.0});
};
+
+ runIntersectionTest<3>(returns, transform2D3D, verbose);
}
-#endif
+template<class Transformation>
+void run2DIntersectionTest(std::vector<bool>& returns, const Transformation& transform, bool verbose)
+{ runIntersectionTest<2>(returns, transform, verbose); }
+
+} // end namespace Dumux
int main (int argc, char *argv[]) try
{
+ using namespace Dumux;
+
// collect returns to determine exit code
std::vector<bool> returns;
constexpr bool verbose = false;
{
- constexpr int dimworld = 2;
- using Vec = Dune::FieldVector<double, dimworld>;
+ using Vec = Dune::FieldVector<double, 2>;
std::cout << "Test for point intersection with 2d geometries in 2d corddim" << std::endl;
for (const double scaling : {1.0, 1e3, 1e12, 1e-12})
for (const auto& 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})
- runIntersectionTest<dimworld>(returns, create2DTransformation<double>(scaling, Vec(translation), angle), verbose);
+ run2DIntersectionTest(returns,
+ make2DTransformation<double>(scaling, Vec(translation), angle, true), verbose);
}
{
- constexpr int dimworld = 3;
- using Vec = Dune::FieldVector<double, dimworld>;
+ using Vec = Dune::FieldVector<double, 3>;
std::cout << "Test for point intersection with 2d geometries in 3d corddim" << std::endl;
for (const double scaling : {1.0, 1e3, 1e12, 1e-12})
for (const auto& 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<dimworld>(returns, create3DTransformation<double>(scaling, Vec(translation), rotAxis, angle), verbose);
+ run3DIntersectionTest(returns,
+ make3DTransformation<double>(scaling, Vec(translation), rotAxis, angle, true), verbose);
}
// determine the exit code
diff --git a/test/common/geometry/test_0d3d_intersection.cc b/test/common/geometry/test_0d3d_intersection.cc
index 16c48998884fa06518e290ec3e4f0529f345fc52..505de29f9f3f885e337fd3ed9ca6044fd01a0502 100644
--- a/test/common/geometry/test_0d3d_intersection.cc
+++ b/test/common/geometry/test_0d3d_intersection.cc
@@ -12,38 +12,11 @@
#include <dune/geometry/multilineargeometry.hh>
#include <dumux/common/math.hh>
-#include <dumux/common/geometry/intersectspointgeometry.hh>
-#ifndef DOXYGEN
-template<class Geometry>
-bool testIntersection(const Geometry& geo,
- const typename Geometry::GlobalCoordinate& p,
- bool foundExpected, bool verbose)
-{
- bool found = Dumux::intersectsPointGeometry(p, geo);
- if (!found && foundExpected)
- {
- 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 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);
-}
+#include "transformation.hh"
+#include "test_intersection.hh"
+
+namespace Dumux {
template<class Transformation>
void runIntersectionTest(std::vector<bool>& returns, const Transformation& transform, bool verbose)
@@ -57,7 +30,7 @@ void runIntersectionTest(std::vector<bool>& returns, const Transformation& trans
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& p) { return transform(p); });
const auto tetrahedron = Geo(Dune::GeometryTypes::tetrahedron, cornersTet);
for (const auto& corner : cornersTet)
@@ -74,7 +47,7 @@ void runIntersectionTest(std::vector<bool>& returns, const Transformation& trans
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); });
+ [&](const auto& p) { return transform(p); });
auto hexahedron = Geo(Dune::GeometryTypes::hexahedron, cornersHex);
for (const auto& corner : cornersHex)
@@ -91,7 +64,7 @@ void runIntersectionTest(std::vector<bool>& returns, const Transformation& trans
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); });
+ [&](const auto& p) { return transform(p); });
auto pyramid = Geo(Dune::GeometryTypes::pyramid, cornersPyramid);
for (const auto& corner : cornersPyramid)
@@ -111,7 +84,7 @@ void runIntersectionTest(std::vector<bool>& returns, const Transformation& trans
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); });
+ [&](const auto& p) { return transform(p); });
auto prism = Geo(Dune::GeometryTypes::prism, cornersPrism);
for (const auto& corner : cornersPrism)
@@ -124,34 +97,12 @@ void runIntersectionTest(std::vector<bool>& returns, const Transformation& trans
returns.push_back(testIntersection(prism, transform({0.25, 0.25, 1.0001}), false, verbose));
}
-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 << "Intersection test with transformation:"
- << " ctype: " << Dune::className<ctype>()
- << ", 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.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
+} // end namespace Dumux
int main (int argc, char *argv[]) try
{
+ using namespace Dumux;
+
// collect returns to determine exit code
std::vector<bool> returns;
constexpr bool verbose = false;
@@ -162,7 +113,7 @@ int main (int argc, char *argv[]) try
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);
+ make3DTransformation<double>(scaling, Vec(translation), rotAxis, angle, true), verbose);
// determine the exit code
if (std::any_of(returns.begin(), returns.end(), std::logical_not<bool>{}))
diff --git a/test/common/geometry/test_intersection.hh b/test/common/geometry/test_intersection.hh
new file mode 100644
index 0000000000000000000000000000000000000000..98fc4a29c7e9cac7c93adc9bdc24d963fdb28f6d
--- /dev/null
+++ b/test/common/geometry/test_intersection.hh
@@ -0,0 +1,66 @@
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup Common
+ * \brief Intersection test helpers
+ */
+#ifndef DUMUX_TEST_INTERSECTION_HH
+#define DUMUX_TEST_INTERSECTION_HH
+
+#include <dumux/common/geometry/intersectspointgeometry.hh>
+
+namespace Dumux {
+
+/*!
+ * \file
+ * \ingroup Common
+ * \brief Test intersection overload for point and geometry
+ */
+template<class Geometry>
+bool testIntersection(const Geometry& geo,
+ const typename Geometry::GlobalCoordinate& p,
+ bool foundExpected, bool verbose)
+{
+ const bool found = intersectsPointGeometry(p, geo);
+ if (!found && foundExpected)
+ {
+ 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 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);
+}
+
+} // end namespace Dumux
+
+#endif
diff --git a/test/common/geometry/transformation.hh b/test/common/geometry/transformation.hh
new file mode 100644
index 0000000000000000000000000000000000000000..a72957687fef34de229bbdbd9f16841f07a48647
--- /dev/null
+++ b/test/common/geometry/transformation.hh
@@ -0,0 +1,123 @@
+/*****************************************************************************
+ * See the file COPYING for full copying permissions. *
+ * *
+ * This program is free software: you can redistribute it and/or modify *
+ * it under the terms of the GNU General Public License as published by *
+ * the Free Software Foundation, either version 3 of the License, or *
+ * (at your option) any later version. *
+ * *
+ * This program is distributed in the hope that it will be useful, *
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of *
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
+ * GNU General Public License for more details. *
+ * *
+ * You should have received a copy of the GNU General Public License *
+ * along with this program. If not, see <http://www.gnu.org/licenses/>. *
+ *****************************************************************************/
+/*!
+ * \file
+ * \ingroup Common
+ * \brief Create a transformation function (translation + rotation)
+ */
+#ifndef DUMUX_TEST_GEOMETRY_TRANSFORMATION_HH
+#define DUMUX_TEST_GEOMETRY_TRANSFORMATION_HH
+
+#include <iostream>
+
+#include <dune/common/classname.hh>
+#include <dune/common/fvector.hh>
+
+#include <dumux/common/math.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup Common
+ * \brief Create a transformation function (translation) in 1D
+ */
+template<class ctype>
+auto make1DTransformation(const ctype scale,
+ const Dune::FieldVector<ctype, 1>& translate,
+ bool verbose = true)
+{
+ if (verbose)
+ std::cout << "Created 1D transformation with"
+ << " ctype: " << Dune::className<ctype>()
+ << ", scaling: " << scale
+ << ", translation: " << translate << std::endl;
+
+ return [=](Dune::FieldVector<ctype, 1> p){
+ p *= scale;
+ p.axpy(scale, translate);
+ auto tp = p;
+ return p;
+ };
+}
+
+/*!
+ * \ingroup Common
+ * \brief Create a transformation function (translation + rotation) in 2D
+ * \note https://en.wikipedia.org/wiki/Rotation_matrix
+ */
+template<class ctype>
+auto make2DTransformation(const ctype scale,
+ const Dune::FieldVector<ctype, 2>& translate,
+ const ctype rotationAngle,
+ bool verbose = true)
+{
+ if (verbose)
+ std::cout << "Created 2D transformation with"
+ << " ctype: " << Dune::className<ctype>()
+ << ", scaling: " << scale
+ << ", translation: " << translate
+ << ", rotationAngle: " << rotationAngle << std::endl;
+
+ using std::sin; using std::cos;
+ const ctype sinAngle = sin(rotationAngle);
+ const ctype cosAngle = cos(rotationAngle);
+ return [=](Dune::FieldVector<ctype, 2> p){
+ p *= scale;
+ p.axpy(scale, translate);
+ auto tp = p;
+ tp[0] = p[0]*cosAngle-p[1]*sinAngle;
+ tp[1] = p[0]*sinAngle+p[1]*cosAngle;
+ return tp;
+ };
+}
+
+/*!
+ * \ingroup Common
+ * \brief Create a transformation function (translation + rotation) in 3D
+ * \note https://en.wikipedia.org/wiki/Rodrigues%27_rotation_formula
+ */
+template<class ctype>
+auto make3DTransformation(const ctype scale,
+ const Dune::FieldVector<ctype, 3>& translate,
+ const Dune::FieldVector<ctype, 3>& rotationAxis,
+ const ctype rotationAngle,
+ bool verbose = true)
+{
+ if (verbose)
+ std::cout << "Created transformation with"
+ << " ctype: " << Dune::className<ctype>()
+ << ", scaling: " << scale
+ << ", translation: " << translate
+ << ", rotationAxis: " << rotationAxis
+ << ", rotationAngle: " << rotationAngle << std::endl;
+
+ using std::sin; using std::cos;
+ const ctype sinAngle = sin(rotationAngle);
+ const ctype cosAngle = cos(rotationAngle);
+ return [=](Dune::FieldVector<ctype, 3> p){
+ p *= scale;
+ 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);
+ };
+}
+
+} // end namespace Dumux
+
+#endif