diff --git a/dumux/multidomain/embedded/circlepoints.hh b/dumux/multidomain/embedded/circlepoints.hh index 64319d829260a58aed54db7a4e6b9b9512745c1c..ae9aee1a705b5f7964650138e2a6437306af6c58 100644 --- a/dumux/multidomain/embedded/circlepoints.hh +++ b/dumux/multidomain/embedded/circlepoints.hh @@ -28,38 +28,46 @@ #include <vector> #include <cmath> +#include <dumux/common/math.hh> #include <dune/common/exceptions.hh> -namespace Dumux { -namespace EmbeddedCoupling { +namespace Dumux::EmbeddedCoupling { /*! * \ingroup EmbeddedCoupling - * \brief returns a vector of points on a circle + * \param points a vector of points to be filled + * \param sincos vector with [sin(a0), cos(a0), sin(a1), cos(a1), ...] for each circumferential sample point ai * \param center the circle center * \param normal the normal to the circle plane * \param radius the circle radius - * \param numPoints the number of points + * \note This version allows for a more efficient circle point generator + * if the circumferential positions are fixed and can be reused. + * In this case, sine and cosine of the corresponding angles can be precomputed. */ -template<class GlobalPosition> -std::vector<GlobalPosition> circlePoints(const GlobalPosition& center, - const GlobalPosition& normal, - const typename GlobalPosition::value_type radius, - const std::size_t numPoints = 20) +template<class GlobalPosition, class Scalar> +void circlePoints(std::vector<GlobalPosition>& points, + const std::vector<Scalar>& sincos, + const GlobalPosition& center, + const GlobalPosition& normal, + const Scalar radius) { - using std::abs; using std::sin; using std::cos; + 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; - static_assert(GlobalPosition::dimension == 3, "Only implemented for world dimension 3"); - std::vector<GlobalPosition> points(numPoints); + // resize the points vector + const std::size_t numPoints = sincos.size()/2; + points.resize(numPoints); // make sure n is a unit vector auto n = normal; n /= n.two_norm(); - // caculate a vector u perpendicular to n + // calculate a vector u perpendicular to n GlobalPosition u; if (abs(n[0]) < eps && abs(n[1]) < eps) if (abs(n[2]) < eps) @@ -75,26 +83,47 @@ std::vector<GlobalPosition> circlePoints(const GlobalPosition& center, auto tangent = crossProduct(u, n); tangent *= radius/tangent.two_norm(); - // the parameter with an offset - ctype t = 0 + 0.1; // insert the vertices for (std::size_t i = 0; i < numPoints; ++i) { - points[i] = GlobalPosition({u[0]*cos(t) + tangent[0]*sin(t) + center[0], - u[1]*cos(t) + tangent[1]*sin(t) + center[1], - u[2]*cos(t) + tangent[2]*sin(t) + center[2]}); + points[i] = GlobalPosition({u[0]*sincos[2*i+1] + tangent[0]*sincos[2*i] + center[0], + u[1]*sincos[2*i+1] + tangent[1]*sincos[2*i] + center[1], + u[2]*sincos[2*i+1] + tangent[2]*sincos[2*i] + center[2]}); + } +} - t += 2*M_PI/numPoints; +/*! + * \ingroup EmbeddedCoupling + * \brief returns a vector of points on a circle + * \param center the circle center + * \param normal the normal to the circle plane + * \param radius the circle radius + * \param numPoints the number of points + */ +template<class GlobalPosition, class Scalar> +std::vector<GlobalPosition> circlePoints(const GlobalPosition& center, + const GlobalPosition& normal, + const Scalar radius, + const std::size_t numPoints = 20) +{ + std::vector<GlobalPosition> points; - // periodic t - if(t > 2*M_PI) - t -= 2*M_PI; + // precompute the sin/cos + std::vector<Scalar> sincos(2*numPoints); + Scalar t = 0 + 0.1; // add arbitrary offset + for (std::size_t i = 0; i < numPoints; ++i) + { + using std::sin; using std::cos; + sincos[2*i] = sin(t); + sincos[2*i + 1] = cos(t); + t += 2*M_PI/numPoints; + if(t > 2*M_PI) t -= 2*M_PI; } + circlePoints(points, sincos, center, normal, radius); return points; } -} // end namespace EmbeddedCoupling -} // end namespace Dumux +} // end namespace Dumux::EmbeddedCoupling #endif diff --git a/test/common/geometry/CMakeLists.txt b/test/common/geometry/CMakeLists.txt index 90a842d59af9ba98dc4ecc13f2698cf1aa7ffa7c..92ea7081395239030c2a8532113e8f4475afef4e 100644 --- a/test/common/geometry/CMakeLists.txt +++ b/test/common/geometry/CMakeLists.txt @@ -8,6 +8,7 @@ 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_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) dune_symlink_to_source_files(FILES ball.msh) dumux_add_test(SOURCES test_intersectionentityset.cc diff --git a/test/common/geometry/test_circlepoints.cc b/test/common/geometry/test_circlepoints.cc new file mode 100644 index 0000000000000000000000000000000000000000..3b5a02750c9d628731436429f5a596572fdb78cf --- /dev/null +++ b/test/common/geometry/test_circlepoints.cc @@ -0,0 +1,44 @@ +#include <config.h> + +#include <cmath> +#include <vector> +#include <algorithm> +#include <dune/common/exceptions.hh> +#include <dune/common/float_cmp.hh> +#include <dune/common/fvector.hh> + +#include <dumux/multidomain/embedded/circlepoints.hh> + +int main (int argc, char *argv[]) +{ + using namespace Dumux; + + using Point = Dune::FieldVector<double, 3>; + const Point center(0.0); + const Point normal({1.0, 1.0, 0.0}); // doesn't need to be a unit normal vector + const double radius = 3.0; + static constexpr std::size_t numPoints = 100; + + // make a circle and and ellipse with equal major radii -> should yield the same points + const auto circle = EmbeddedCoupling::circlePoints(center, normal, radius, numPoints); + const auto checkRadius = [&](const auto& p){ return Dune::FloatCmp::eq(p.two_norm(), radius, 1e-7); }; + if (!std::all_of(circle.begin(), circle.end(), checkRadius)) + DUNE_THROW(Dune::Exception, "Not all point have the same radius!"); + + Point unitNormal = normal; + unitNormal /= unitNormal.two_norm(); + bool correctOrientation = true; + for (int i = 1; i < circle.size(); ++i) + { + auto n = crossProduct(circle[i]-center, circle[i-1]-center); + n /= n.two_norm(); + correctOrientation &= Dune::FloatCmp::eq(std::abs(n*unitNormal), 1.0, 1e-10); + } + if (!correctOrientation) + DUNE_THROW(Dune::Exception, "Not all points have the correct orientation!"); + + if (circle.size() != numPoints) + DUNE_THROW(Dune::Exception, "Wrong number of points!"); + + return 0; +}