diff --git a/dumux/common/geometry/CMakeLists.txt b/dumux/common/geometry/CMakeLists.txt
index c18c814cf536b887bf2345cec4220294a1367453..f5e830fc3a5a0b382f1f805664295de1142f487f 100644
--- a/dumux/common/geometry/CMakeLists.txt
+++ b/dumux/common/geometry/CMakeLists.txt
@@ -1,6 +1,7 @@
install(FILES
boundingboxtree.hh
diameter.hh
+distance.hh
geometricentityset.hh
geometryintersection.hh
grahamconvexhull.hh
@@ -9,6 +10,7 @@ intersectionentityset.hh
intersectspointgeometry.hh
intersectspointsimplex.hh
makegeometry.hh
+normal.hh
refinementquadraturerule.hh
triangulation.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/common/geometry)
diff --git a/dumux/common/geometry/distance.hh b/dumux/common/geometry/distance.hh
new file mode 100644
index 0000000000000000000000000000000000000000..47fccf32ea9057aed73b1a2b5ce77c33783d20b5
--- /dev/null
+++ b/dumux/common/geometry/distance.hh
@@ -0,0 +1,181 @@
+/*****************************************************************************
+ * 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 Geometry
+ * \brief Helper functions for distance queries
+ */
+#ifndef DUMUX_GEOMETRY_DISTANCE_HH
+#define DUMUX_GEOMETRY_DISTANCE_HH
+
+#include <dune/common/fvector.hh>
+#include <dune/geometry/quadraturerules.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup Geometry
+ * \brief Compute the average distance from a point to a geometry by integration
+ */
+template<class Geometry>
+inline typename Geometry::ctype
+averageDistancePointGeometry(const typename Geometry::GlobalCoordinate& p,
+ const Geometry& geometry,
+ std::size_t integrationOrder = 2)
+{
+ typename Geometry::ctype avgDist = 0.0;
+ const auto& quad = Dune::QuadratureRules<typename Geometry::ctype, Geometry::mydimension>::rule(geometry.type(), integrationOrder);
+ for (const auto& qp : quad)
+ avgDist += (geometry.global(qp.position())-p).two_norm()*qp.weight()*geometry.integrationElement(qp.position());
+ return avgDist/geometry.volume();
+}
+
+/*!
+ * \ingroup Geometry
+ * \brief Compute the distance from a point to a line through the points a and b
+ */
+template<class Point>
+inline typename Point::value_type
+distancePointLine(const Point& p, const Point& a, const Point& b)
+{
+ const auto ab = b - a;
+ const auto t = (p - a)*ab/ab.two_norm2();
+ auto proj = a;
+ proj.axpy(t, ab);
+ return (proj - p).two_norm();
+}
+
+/*!
+ * \ingroup Geometry
+ * \brief Compute the distance from a point to a line given by a geometry with two corners
+ * \note We currently lack the a representation of a line geometry. This convenience function
+ * assumes a segment geometry (with two corners) is passed which represents a line geometry.
+ */
+template<class Geometry>
+inline typename Geometry::ctype
+distancePointLine(const typename Geometry::GlobalCoordinate& p, const Geometry& geometry)
+{
+ static_assert(Geometry::mydimension == 1, "Geometry has to be a line");
+ const auto& a = geometry.corner(0);
+ const auto& b = geometry.corner(1);
+ return distancePointLine(p, a, b);
+}
+
+/*!
+ * \ingroup Geometry
+ * \brief Compute the distance from a point to the segment connecting the points a and b
+ */
+template<class Point>
+inline typename Point::value_type
+distancePointSegment(const Point& p, const Point& a, const Point& b)
+{
+ const auto ab = b - a;
+ auto t = (p - a)*ab;
+
+ if (t <= 0.0)
+ return (a - p).two_norm();
+
+ const auto lengthSq = ab.two_norm2();
+ if (t >= lengthSq)
+ return (b - p).two_norm();
+
+ auto proj = a;
+ proj.axpy(t/lengthSq, ab);
+ return (proj - p).two_norm();
+}
+
+/*!
+ * \ingroup Geometry
+ * \brief Compute the distance from a point to a given segment geometry
+ */
+template<class Geometry>
+inline typename Geometry::ctype
+distancePointSegment(const typename Geometry::GlobalCoordinate& p, const Geometry& geometry)
+{
+ static_assert(Geometry::mydimension == 1, "Geometry has to be a segment");
+ const auto& a = geometry.corner(0);
+ const auto& b = geometry.corner(1);
+ return distancePointSegment(p, a, b);
+}
+
+/*!
+ * \ingroup Geometry
+ * \brief Compute the shortest distance between two points
+ */
+template<class ctype, int dimWorld>
+inline ctype distance(const Dune::FieldVector<ctype, dimWorld>& a,
+ const Dune::FieldVector<ctype, dimWorld>& b)
+{ return (a-b).two_norm(); }
+
+
+
+namespace Detail {
+
+// helper struct to compute distance between two geometries, specialized below
+template<class Geo1, class Geo2,
+ int dimWorld = Geo1::coorddimension,
+ int dim1 = Geo1::mydimension, int dim2 = Geo2::mydimension>
+struct GeometryDistance
+{
+ static_assert(Geo1::coorddimension == Geo2::coorddimension, "Geometries have to have the same coordinate dimensions");
+ static auto distance(const Geo1& geo1, const Geo2& geo2)
+ {
+ DUNE_THROW(Dune::NotImplemented, "Geometry distance computation not implemented for dimworld = "
+ << dimWorld << ", dim1 = " << dim1 << ", dim2 = " << dim2);
+ }
+};
+
+// distance point-point
+template<class Geo1, class Geo2, int dimWorld>
+struct GeometryDistance<Geo1, Geo2, dimWorld, 0, 0>
+{
+ static_assert(Geo1::coorddimension == Geo2::coorddimension, "Geometries have to have the same coordinate dimensions");
+ static auto distance(const Geo1& geo1, const Geo2& geo2)
+ { return Dumux::distance(geo1.corner(0), geo2.corner(0)); }
+};
+
+// distance segment-point
+template<class Geo1, class Geo2, int dimWorld>
+struct GeometryDistance<Geo1, Geo2, dimWorld, 1, 0>
+{
+ static_assert(Geo1::coorddimension == Geo2::coorddimension, "Geometries have to have the same coordinate dimensions");
+ static auto distance(const Geo1& geo1, const Geo2& geo2)
+ { return distancePointSegment(geo2.corner(0), geo1); }
+};
+
+// distance point-segment
+template<class Geo1, class Geo2, int dimWorld>
+struct GeometryDistance<Geo1, Geo2, dimWorld, 0, 1>
+{
+ static_assert(Geo1::coorddimension == Geo2::coorddimension, "Geometries have to have the same coordinate dimensions");
+ static auto distance(const Geo1& geo1, const Geo2& geo2)
+ { return distancePointSegment(geo1.corner(0), geo2); }
+};
+
+} // end namespace Detail
+
+/*!
+ * \ingroup Geometry
+ * \brief Compute the shortest distance between two geometries
+ */
+template<class Geo1, class Geo2>
+inline auto distance(const Geo1& geo1, const Geo2& geo2)
+{ return Detail::GeometryDistance<Geo1, Geo2>::distance(geo1, geo2); }
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/common/geometry/normal.hh b/dumux/common/geometry/normal.hh
new file mode 100644
index 0000000000000000000000000000000000000000..a6e95ecb170b8c75a6df877f8667742f85585ee7
--- /dev/null
+++ b/dumux/common/geometry/normal.hh
@@ -0,0 +1,77 @@
+/*****************************************************************************
+ * 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 Geometry
+ * \brief Helper functions for normals
+ */
+#ifndef DUMUX_GEOMETRY_NORMAL_HH
+#define DUMUX_GEOMETRY_NORMAL_HH
+
+#include <algorithm>
+#include <dune/common/float_cmp.hh>
+
+namespace Dumux {
+
+/*!
+ * \ingroup Geometry
+ * \brief Create a vector normal to the given one (v is expected to be non-zero)
+ * \note This returns some orthogonal vector with arbitrary length
+ */
+template<class Vector>
+inline Vector normal(const Vector& v)
+{
+ static_assert(Vector::size() > 1, "normal expects a coordinate dimension > 1");
+
+ if constexpr (Vector::size() == 2)
+ return Vector({-v[1], v[0]});
+
+ const auto it = std::find_if(v.begin(), v.end(), [](const auto& x) { return Dune::FloatCmp::ne(x, 0.0); });
+ const auto index = std::distance(v.begin(), it);
+ if (index != Vector::size()-1)
+ {
+ Vector normal(0.0);
+ normal[index] = -v[index+1];
+ normal[index+1] = v[index];
+ return normal;
+ }
+ else
+ {
+ Vector normal(0.0);
+ normal[index-1] = -v[index];
+ normal[index] = v[index-1];
+ return normal;
+
+ }
+}
+
+/*!
+ * \ingroup Geometry
+ * \brief Create a vector normal to the given one (v is expected to be non-zero)
+ * \note This returns some orthogonal vector with unit length
+ */
+template<class Vector>
+inline Vector unitNormal(const Vector& v)
+{
+ auto normal = Dumux::normal(v);
+ normal /= normal.two_norm();
+ return normal;
+}
+
+} // end namespace Dumux
+
+#endif
diff --git a/dumux/multidomain/embedded/circlepoints.hh b/dumux/multidomain/embedded/circlepoints.hh
index ae9aee1a705b5f7964650138e2a6437306af6c58..6145915bb08de6c6c0246cc95eda90b8ca9a6091 100644
--- a/dumux/multidomain/embedded/circlepoints.hh
+++ b/dumux/multidomain/embedded/circlepoints.hh
@@ -28,9 +28,11 @@
#include <vector>
#include <cmath>
-#include <dumux/common/math.hh>
#include <dune/common/exceptions.hh>
+#include <dumux/common/math.hh>
+#include <dumux/common/geometry/normal.hh>
+
namespace Dumux::EmbeddedCoupling {
/*!
@@ -54,11 +56,6 @@ void circlePoints(std::vector<GlobalPosition>& points,
assert(sincos.size() % 2 == 0 && "Sample angles have to be pairs of sin/cos so size needs to be even.");
static_assert(GlobalPosition::dimension == 3, "Only implemented for world dimension 3");
- using std::abs;
- using ctype = typename GlobalPosition::value_type;
-
- constexpr ctype eps = 1.5e-7;
-
// resize the points vector
const std::size_t numPoints = sincos.size()/2;
points.resize(numPoints);
@@ -68,16 +65,7 @@ void circlePoints(std::vector<GlobalPosition>& points,
n /= n.two_norm();
// calculate a vector u perpendicular to n
- GlobalPosition u;
- if (abs(n[0]) < eps && abs(n[1]) < eps)
- if (abs(n[2]) < eps)
- DUNE_THROW(Dune::MathError, "The normal vector has to be non-zero!");
- else
- u = {0, 1, 0};
- else
- u = {-n[1], n[0], 0};
-
- u *= radius/u.two_norm();
+ auto u = unitNormal(n); u*= radius;
// the circle parameterization is p(t) = r*cos(t)*u + r*sin(t)*(n x u) + c
auto tangent = crossProduct(u, n);
diff --git a/dumux/multidomain/embedded/couplingmanager1d3d.hh b/dumux/multidomain/embedded/couplingmanager1d3d.hh
index 0b9c09038ccae294cb8eee4ce11746cde9868b9f..53e647e594d54014599d6d30fa6a6e1c62b7ed31 100644
--- a/dumux/multidomain/embedded/couplingmanager1d3d.hh
+++ b/dumux/multidomain/embedded/couplingmanager1d3d.hh
@@ -33,6 +33,7 @@
#include <dumux/common/properties.hh>
#include <dumux/common/indextraits.hh>
+#include <dumux/common/geometry/distance.hh>
#include <dumux/multidomain/embedded/pointsourcedata.hh>
#include <dumux/multidomain/embedded/integrationpointsource.hh>
#include <dumux/multidomain/embedded/couplingmanagerbase.hh>
@@ -1040,7 +1041,7 @@ public:
this->pointSourceData().emplace_back(std::move(psData));
// compute average distance to bulk cell
- this->averageDistanceToBulkCell().push_back(this->computeDistance(outside.geometry(), globalPos));
+ this->averageDistanceToBulkCell().push_back(averageDistancePointGeometry(globalPos, outside.geometry()));
// export the lowdim coupling stencil
// we insert all vertices / elements and make it unique later
diff --git a/dumux/multidomain/embedded/couplingmanagerbase.hh b/dumux/multidomain/embedded/couplingmanagerbase.hh
index b130b9edd1021e8fbe7d2675fa4e68c3b90dbafb..7e0bc1621e0b82637b51466673d4f64ede749f16 100644
--- a/dumux/multidomain/embedded/couplingmanagerbase.hh
+++ b/dumux/multidomain/embedded/couplingmanagerbase.hh
@@ -36,6 +36,7 @@
#include <dune/geometry/quadraturerules.hh>
#include <dumux/common/properties.hh>
+#include <dumux/common/geometry/distance.hh>
#include <dumux/common/geometry/intersectingentities.hh>
#include <dumux/discretization/method.hh>
#include <dumux/multidomain/couplingmanager.hh>
@@ -318,7 +319,7 @@ public:
this->pointSourceData().emplace_back(std::move(psData));
// compute average distance to bulk cell
- averageDistanceToBulkCell_.push_back(computeDistance(outside.geometry(), globalPos));
+ averageDistanceToBulkCell_.push_back(averageDistancePointGeometry(globalPos, outside.geometry()));
// export the lowdim coupling stencil
// we insert all vertices / elements and make it unique later
@@ -477,16 +478,6 @@ protected:
glue_->build(bulkGridGeometry.boundingBoxTree(), lowDimGridGeometry.boundingBoxTree());
}
- template<class Geometry, class GlobalPosition>
- Scalar computeDistance(const Geometry& geometry, const GlobalPosition& p) const
- {
- Scalar avgDist = 0.0;
- const auto& quad = Dune::QuadratureRules<Scalar, bulkDim>::rule(geometry.type(), 5);
- for (auto&& qp : quad)
- avgDist += (geometry.global(qp.position())-p).two_norm()*qp.weight();
- return avgDist;
- }
-
//! Return reference to point source data vector member
std::vector<PointSourceData>& pointSourceData()
{ return pointSourceData_; }
diff --git a/test/common/geometry/CMakeLists.txt b/test/common/geometry/CMakeLists.txt
index f1a6d093ce7bb5caa35dd1dfdf4a9921e674c87c..6e57be7f6d0a318670f5bbc84edc0d64000ddcc0 100644
--- a/test/common/geometry/CMakeLists.txt
+++ b/test/common/geometry/CMakeLists.txt
@@ -6,6 +6,8 @@ dumux_add_test(SOURCES test_1d3d_intersection.cc LABELS unit)
dumux_add_test(SOURCES test_1d2d_intersection.cc LABELS unit)
dumux_add_test(SOURCES test_2d2d_intersection.cc LABELS unit)
dumux_add_test(SOURCES test_2d3d_intersection.cc LABELS unit)
+dumux_add_test(SOURCES test_distance.cc LABELS unit)
+dumux_add_test(SOURCES test_normal.cc LABELS unit)
dumux_add_test(SOURCES test_graham_convex_hull.cc LABELS unit)
dumux_add_test(SOURCES test_intersectingentity_cartesiangrid.cc LABELS unit)
dumux_add_test(SOURCES test_circlepoints.cc LABELS unit)
diff --git a/test/common/geometry/test_distance.cc b/test/common/geometry/test_distance.cc
new file mode 100644
index 0000000000000000000000000000000000000000..3fbed2c96d991718eaf916e83a39bc3f6a2d1557
--- /dev/null
+++ b/test/common/geometry/test_distance.cc
@@ -0,0 +1,164 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * 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
+ * \brief Test for distance computations.
+ */
+#include <config.h>
+
+#include <vector>
+#include <random>
+#include <cmath>
+
+#include <dune/common/fvector.hh>
+#include <dune/common/float_cmp.hh>
+#include <dune/common/exceptions.hh>
+
+#include <dune/geometry/affinegeometry.hh>
+#include <dune/geometry/multilineargeometry.hh>
+
+#include <dumux/common/geometry/distance.hh>
+#include <dumux/common/geometry/normal.hh>
+
+// helper function to make point geometry from field vector
+template<class Point, int dimWorld>
+Point makePoint(const Dune::FieldVector<typename Point::ctype, dimWorld>& p)
+{
+ using GlobalPosition = Dune::FieldVector<typename Point::ctype, dimWorld>;
+ using Corners = std::vector<GlobalPosition>;
+ return { Dune::GeometryTypes::vertex, Corners{{p}} };
+}
+
+// helper function to make a segment geometry from field vectors
+template<class Segment, int dimWorld>
+Segment makeSegment(const Dune::FieldVector<typename Segment::ctype, dimWorld>& p1,
+ const Dune::FieldVector<typename Segment::ctype, dimWorld>& p2)
+{
+ using GlobalPosition = Dune::FieldVector<typename Segment::ctype, dimWorld>;
+ using Corners = std::vector<GlobalPosition>;
+ return { Dune::GeometryTypes::line, Corners{{p1, p2}} };
+}
+
+// sample a point on a sphere with the given radius
+template<class Point>
+Point samplePointOnSphere(typename Point::ctype radius)
+{
+ static std::default_random_engine generator(std::random_device{}());
+ static constexpr int dimWorld = Point::coorddimension;
+ static_assert(dimWorld > 1, "This only works in 2d or 3d");
+
+ // sample a point from random coordinates
+ std::uniform_real_distribution<double> distro(-radius, radius);
+
+ auto p = dimWorld == 2 ? makePoint<Point, dimWorld>({distro(generator), distro(generator)})
+ : makePoint<Point, dimWorld>({distro(generator), distro(generator), distro(generator)});
+
+ // project it onto the sphere
+ auto pos = p.corner(0);
+ pos *= radius/pos.two_norm();
+ p = makePoint<Point>(pos);
+
+ return p;
+}
+
+// checks the result of a distance computation
+template<class ctype>
+void checkGeometryDistance(ctype expected, ctype computed, const std::string& geometryPairName)
+{
+ if (Dune::FloatCmp::ne(expected, computed))
+ DUNE_THROW(Dune::InvalidStateException, "Unexpected " << geometryPairName << " distance ");
+}
+
+// checks the distances between various points with points/segments/lines
+template<class Point, class Segment>
+void runTests()
+{
+ using namespace Dumux;
+
+ using ctype = typename Point::ctype;
+ using GlobalPosition = typename Point::GlobalCoordinate;
+
+ const auto origin = makePoint<Point>(GlobalPosition(0.0));
+ for (ctype scale : {1e-12, 1.0, 1e12})
+ {
+ // check distances for some random points in a known distance
+ for (int i = 0; i < 20; ++i)
+ {
+ const auto p = samplePointOnSphere<Point>(scale);
+
+ // test point-point distance
+ checkGeometryDistance(scale, distance(origin, p), "point-point");
+
+ // test point-segment distance (where projection is on the segment)
+ auto n = normal(p.corner(0));
+ n *= scale/n.two_norm();
+
+ auto segment = makeSegment<Segment>(origin.corner(0) + n, p.corner(0) + n);
+ checkGeometryDistance(scale, distance(segment, p), "point-segment");
+
+ // test point-segment distance (where projection is NOT on segment)
+ auto v = p.corner(0);
+
+ using std::sqrt;
+ auto segment2 = makeSegment<Segment>(segment.corner(0) - v, segment.corner(1) - v);
+ checkGeometryDistance(sqrt(2.0)*scale, distance(segment2, p), "segment-segment");
+
+ v *= 2.0;
+ auto segment3 = makeSegment<Segment>(segment.corner(0) + v, segment.corner(1) + v);
+ checkGeometryDistance(sqrt(2.0)*scale, distance(segment2, p), "segment-segment");
+
+ // for lines, we should always get a distance of scale
+ checkGeometryDistance(scale, distancePointLine(p.corner(0), segment.corner(0), segment.corner(1)), "segment-segment");
+ checkGeometryDistance(scale, distancePointLine(p.corner(0), segment2.corner(0), segment2.corner(1)), "segment-segment");
+ checkGeometryDistance(scale, distancePointLine(p.corner(0), segment3.corner(0), segment3.corner(1)), "segment-segment");
+ }
+ }
+}
+
+int main(int argc, char** argv) try
+{
+ // affine geometries (2d)
+ runTests< Dune::AffineGeometry<double, 0, 2>, Dune::AffineGeometry<double, 1, 2> >();
+
+ // affine geometries (3d)
+ runTests< Dune::AffineGeometry<double, 0, 3>, Dune::AffineGeometry<double, 1, 3> >();
+
+ // multilinear geometries (2d)
+ runTests< Dune::MultiLinearGeometry<double, 0, 2>, Dune::MultiLinearGeometry<double, 1, 2> >();
+
+ // multilinear geometries (3d)
+ runTests< Dune::MultiLinearGeometry<double, 0, 3>, Dune::MultiLinearGeometry<double, 1, 3> >();
+
+ // mixed types (2d)
+ runTests< Dune::AffineGeometry<double, 0, 2>, Dune::MultiLinearGeometry<double, 1, 2> >();
+ runTests< Dune::MultiLinearGeometry<double, 0, 2>, Dune::AffineGeometry<double, 1, 2> >();
+
+ // mixed types (3d)
+ runTests< Dune::AffineGeometry<double, 0, 3>, Dune::MultiLinearGeometry<double, 1, 3> >();
+ runTests< Dune::MultiLinearGeometry<double, 0, 3>, Dune::AffineGeometry<double, 1, 3> >();
+
+ return 0;
+}
+// //////////////////////////////////
+// Error handler
+// /////////////////////////////////
+catch (const Dune::Exception& e) {
+ std::cout << e << std::endl;
+ return 1;
+}
diff --git a/test/common/geometry/test_normal.cc b/test/common/geometry/test_normal.cc
new file mode 100644
index 0000000000000000000000000000000000000000..2a95def1b644225b2042a9a4e3a31b8b6efd397a
--- /dev/null
+++ b/test/common/geometry/test_normal.cc
@@ -0,0 +1,94 @@
+// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
+// vi: set et ts=4 sw=4 sts=4:
+/*****************************************************************************
+ * 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
+ * \brief Test for normal computation
+ */
+#include <config.h>
+
+#include <random>
+#include <algorithm>
+
+#include <dune/common/exceptions.hh>
+#include <dune/common/fvector.hh>
+
+#include <dumux/common/geometry/normal.hh>
+
+template<class P>
+void testOrthogonality(const P& normal, const P& vector, double scale, const std::string& caseName = "")
+{
+ const auto eps = Dune::FloatCmp::DefaultEpsilon<double, Dune::FloatCmp::absolute>::value()*scale*100;
+ if (Dune::FloatCmp::ne<double, Dune::FloatCmp::absolute>(normal*vector, 0.0, eps))
+ DUNE_THROW(Dune::Exception, "Normal not orthogonal, n: " << normal << ", p: " << vector
+ << ", dim: " << P::size() << ", sp: " << normal*vector
+ << ", scale: " << scale << ", epsilon: " << eps << ", case: " << caseName);
+}
+
+template<int dim, std::size_t repetitions = 100000>
+void testNormal()
+{
+ std::default_random_engine gen(std::random_device{}());
+ std::uniform_real_distribution<double> rnd(-1.0, 1.0), rndIdx(0.0, dim-1);
+ Dune::FieldVector<double, dim> p(0.0);
+
+ for (double scale : {1.0, 1e3, 1e-12, 1e12})
+ {
+ for (std::size_t i = 0; i < repetitions; ++i)
+ {
+ std::generate(p.begin(), p.end(), [&]{ return rnd(gen)*scale; });
+ const auto n = Dumux::normal(p);
+ testOrthogonality(n, p, scale, "rnd");
+ }
+
+ // one coordinate != 0
+ for (std::size_t i = 0; i < repetitions; ++i)
+ {
+ p = 0.0;
+ const auto randomIndex = (std::size_t) std::round(rndIdx(gen));
+ p[randomIndex] = scale;
+ testOrthogonality(Dumux::normal(p), p, scale, "rnd-one-nonzero");
+ p[randomIndex] = scale*1e-20;
+ testOrthogonality(Dumux::normal(p), p, scale, "rnd-one-small");
+ }
+
+ // two coordinates != 0
+ for (std::size_t i = 0; i < repetitions; ++i)
+ {
+ p = 0.0;
+ const auto randomIndex0 = (std::size_t) std::round(rndIdx(gen));
+ const auto randomIndex1 = (std::size_t) std::round(rndIdx(gen));
+ p[randomIndex0] = scale;
+ p[randomIndex1] = scale;
+ testOrthogonality(Dumux::normal(p), p, scale, "rnd-two-nonzero");
+ p[randomIndex0] = scale*1e-20;
+ p[randomIndex1] = scale*1e-20;
+ testOrthogonality(Dumux::normal(p), p, scale, "rnd-two-small");
+ }
+ }
+}
+
+int main(int argc, char** argv)
+{
+ testNormal<2>();
+ testNormal<3>();
+ testNormal<4>();
+ testNormal<10>();
+
+ return 0;
+}