diff --git a/dumux/common/geometry/CMakeLists.txt b/dumux/common/geometry/CMakeLists.txt index 41dbff92bf6e1cb405241110acc2d25bf25f3766..f5e830fc3a5a0b382f1f805664295de1142f487f 100644 --- a/dumux/common/geometry/CMakeLists.txt +++ b/dumux/common/geometry/CMakeLists.txt @@ -10,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/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/test/common/geometry/CMakeLists.txt b/test/common/geometry/CMakeLists.txt index 42222c5fff97304266a86384a1e8436ae4b74b0b..6e57be7f6d0a318670f5bbc84edc0d64000ddcc0 100644 --- a/test/common/geometry/CMakeLists.txt +++ b/test/common/geometry/CMakeLists.txt @@ -7,6 +7,7 @@ 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_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; +}