From cec0a50c71507de8441a2f8b978b8acebd5591ca Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Fri, 12 Jun 2020 14:49:55 +0200
Subject: [PATCH 1/3] [md][embedded][geometry] Add cylinder integration tool

---
 .../embedded/cylinderintegration.hh           | 468 ++++++++++++++++++
 1 file changed, 468 insertions(+)
 create mode 100644 dumux/multidomain/embedded/cylinderintegration.hh

diff --git a/dumux/multidomain/embedded/cylinderintegration.hh b/dumux/multidomain/embedded/cylinderintegration.hh
new file mode 100644
index 0000000000..e54fd5515e
--- /dev/null
+++ b/dumux/multidomain/embedded/cylinderintegration.hh
@@ -0,0 +1,468 @@
+// -*- 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
+ * \ingroup EmbeddedCoupling
+ * \brief Integration over cylindrical and elliptic cylindrical domains
+ * Lowest order integration formulas that are mostly useful to evenly distribute
+ * mass in a cylindrical domain.
+ */
+
+#ifndef DUMUX_MULTIDOMAIN_EMBEDDED_CYLINDERINTEGRATOR_HH
+#define DUMUX_MULTIDOMAIN_EMBEDDED_CYLINDERINTEGRATOR_HH
+
+#include <cmath>
+#include <algorithm>
+#include <numeric>
+
+#include <dune/common/fvector.hh>
+
+#include <dumux/common/math.hh>
+#include <dumux/multidomain/embedded/circlepoints.hh>
+
+namespace Dumux::EmbeddedCoupling {
+
+/*!
+ * \ingroup EmbeddedCoupling
+ * \brief Helper class to integrate over a cylinder domain
+ * \note This is mostly useful if the integral is known and the
+ *       integral mass has to be locally distributed to some non-matching domain
+ * \note The algorithm creates (almost) evenly spaced area elements on a circle
+ *       and extrudes this pattern into a cylinder. Each sub-cell is an integration
+ *       point with constant ansatz on the sub-cell (lowest order integration formula).
+ *       Hence to improve the quality of the integral, increase the sample size.
+ * \note The volume integral of a constant function is exact
+ * See Beckers & Beckers (2012) doi:10.1016/j.comgeo.2012.01.011 for circle distribution
+ * and T. Koch PhD thesis (2020), Section 7.2.4.
+ */
+template<class Scalar>
+class CylinderIntegration
+{
+    using GlobalPosition = Dune::FieldVector<Scalar, 3>;
+
+public:
+    /*!
+     * \brief Constructor
+     * \param rStep characteristic relative integration length (real number between 0 and 1)
+     * \note half the characteristic length means 2*2*2=8 times more integration points
+     */
+    CylinderIntegration(const Scalar rStep)
+    {
+        using std::ceil; using std::clamp;
+        rSamples_ = std::max<std::size_t>(1, ceil(1.0/clamp(rStep, 1.0/std::numeric_limits<std::size_t>::max(), 1.0)));
+        initialize_();
+    }
+
+    /*!
+     * \brief Constructor
+     * \param rStep characteristic relative integration length in r-direction (real number between 0 and 1)
+     * \param zSamples characteristic relative integration length in z-direction (real number between 0 and 1)
+     * \note Use this constructor to achieve a non-balanced (away from 1) aspect ratio between r and z-direction
+     */
+    CylinderIntegration(const Scalar rStep, const Scalar zStep)
+    : zStepFixed_(true)
+    {
+        using std::ceil; using std::clamp;
+        rSamples_ = std::max<std::size_t>(1, ceil(1.0/clamp(rStep, 1.0/std::numeric_limits<std::size_t>::max(), 1.0)));
+        zSamples_ = std::max<std::size_t>(1, ceil(1.0/clamp(zStep, 1.0/std::numeric_limits<std::size_t>::max(), 1.0)));
+        initialize_();
+    }
+
+    /*!
+     * \brief Set the geometry of the cylinder
+     * \param bottomCenter bottom center position
+     * \param topCenter top center position
+     * \param radius cylinder radius
+     * \param verbosity the verbosity level (default: 0 -> no terminal output)
+     */
+    void setGeometry(const GlobalPosition& bottomCenter,
+                     const GlobalPosition& topCenter,
+                     const Scalar radius,
+                     int verbosity = 0)
+    {
+        const auto zAxis = (topCenter-bottomCenter);
+        const auto height = zAxis.two_norm();
+        const auto zAxisUnitVector = zAxis/height;
+
+        // compute step size in r-direction
+        const auto rStep = radius/Scalar(rSamples_);
+
+        // compute zSamples from r samples if not specified by the user
+        using std::ceil;
+        if (!zStepFixed_)
+            zSamples_ = std::max<std::size_t>(1, ceil(height/rStep));
+        const auto zStep = height/Scalar(zSamples_);
+
+        // compute offsets for index calculation
+        auto kOffset = numRingSamples_;
+        std::partial_sum(kOffset.begin(), kOffset.end(), kOffset.begin());
+
+        // compute total number of samples
+        numPoints_ = zSamples_*circleSamples_;
+        if (verbosity > 0)
+            std::cout << "CylinderIntegration -- number of integration points: " << numPoints_ << std::endl;
+
+        // resize integration points and integration elements
+        integrationElement_.resize(numPoints_);
+        points_.resize(numPoints_);
+
+        // reserve enough memory for the points on each ring
+        std::vector<GlobalPosition> ringPoints;
+        ringPoints.reserve(numRingSamples_.back());
+
+        // compute integration points and integration elements
+        for (std::size_t i = 0; i < zSamples_; ++i)
+        {
+            // for each cylinder slice i
+            auto sliceCenter = bottomCenter;
+            sliceCenter.axpy((Scalar(i)+0.5)*zStep, zAxisUnitVector);
+
+            // generate circle points for each ring j
+            for (std::size_t j = 0; j < rSamples_; ++j)
+            {
+                const auto r = (Scalar(j)+0.5)*rStep;
+                circlePoints(ringPoints, sincos_[j], sliceCenter, zAxisUnitVector, r);
+                const auto ringOffest = j > 0 ? kOffset[j-1] : 0;
+                const auto ringSamples = numRingSamples_[j];
+                // volume of each element in ring
+                // total ring volume is given by M_PI*rStep*rStep*zStep*((j+1)^2 - j^2)
+                const auto integrationElement = M_PI*rStep*rStep*zStep*(1.0 + 2.0*Scalar(j))/ringSamples;
+                for (std::size_t k = 0; k < ringSamples; ++k)
+                {
+                    const std::size_t idx = k + ringOffest + circleSamples_*i;
+                    points_[idx] = ringPoints[k];
+                    integrationElement_[idx] = integrationElement;
+                }
+            }
+        }
+    }
+
+    //! The integration element of the ith integration point
+    Scalar integrationElement(std::size_t i) const
+    { return integrationElement_[i]; }
+
+    //! The ith integration point
+    const GlobalPosition& integrationPoint(std::size_t i) const
+    { return points_[i]; }
+
+    //! The number of integration points
+    std::size_t size() const
+    { return numPoints_; }
+
+private:
+    void initialize_()
+    {
+        // precompute number of cells in ring and total circle samples
+        circleSamples_ = 0;
+        numRingSamples_.resize(rSamples_);
+        for (int j = 0; j < rSamples_; ++j)
+        {
+            // number of cells in ring j
+            using std::floor;
+            numRingSamples_[j] = floor(2*M_PI*(Scalar(j)+0.5));
+            circleSamples_ += numRingSamples_[j];
+        }
+
+        // optimization because calling too many sin/cos can be really expensive
+        sincos_.resize(rSamples_);
+        for (int j = 0; j < rSamples_; ++j)
+        {
+            const auto numPoints = numRingSamples_[j];
+            sincos_[j].resize(2*numPoints);
+            Scalar t = 0 + 0.1;
+            for (std::size_t i = 0; i < numPoints; ++i)
+            {
+                using std::sin; using std::cos;
+                sincos_[j][2*i] = sin(t);
+                sincos_[j][2*i + 1] = cos(t);
+                t += 2*M_PI/numPoints;
+                if (t > 2*M_PI) t -= 2*M_PI;
+            }
+        }
+    }
+
+    bool zStepFixed_ = false;
+    std::size_t zSamples_, rSamples_, numPoints_, circleSamples_;
+    std::vector<Scalar> integrationElement_;
+    std::vector<GlobalPosition> points_;
+    std::vector<std::size_t> numRingSamples_;
+    std::vector<std::vector<Scalar>> sincos_;
+};
+
+namespace Impl {
+//! check if a point is in an ellipse
+template<class GlobalPosition>
+inline bool pointInEllipse(const GlobalPosition& p,
+                           const GlobalPosition& center,
+                           const GlobalPosition& firstAxis,
+                           const GlobalPosition& secondAxis,
+                           const GlobalPosition& normal,
+                           const typename GlobalPosition::value_type a,
+                           const typename GlobalPosition::value_type b)
+{
+    const auto d = p-center;
+    // early return if point is not on the ellipse plane
+    if (d*normal > 1e-7*a)
+        return false;
+
+    // check if it's actually within the ellipse
+    const auto da = (d*firstAxis);
+    const auto db = (d*secondAxis);
+
+    return (da*da/(a*a) + db*db/(b*b) < 1.0);
+}
+} // end namespace Impl
+
+/*!
+ * \ingroup EmbeddedCoupling
+ * \brief Helper class to integrate over an elliptic cylinder domain
+ * \note The cylinder caps do not have to be orthogonal to the centerline axis but top and bottom cap are parallel planes
+ * \note This is mostly useful if the integral is known and the
+ *       integral mass has to be locally distributed to some non-matching domain
+ * \note The algorithm is based on evenly spaced area elements with rejection sampling.
+ *       To improve the quality of the integral, increase the sample size.
+ * \note The volume integral of a constant function is exact
+ * See Koch et al. (2020) https://doi.org/10.1016/j.jcp.2020.109369
+ */
+template<class Scalar>
+class EllipticCylinderIntegration
+{
+    using GlobalPosition = Dune::FieldVector<Scalar, 3>;
+
+public:
+    /*!
+     * \brief Constructor
+     * \param aStep characteristic relative integration length (real number between 0 and 1)
+     * \note half the characteristic length means 2*2*2=8 times more integration points
+     */
+    explicit EllipticCylinderIntegration(const Scalar aStep)
+    : aCharLength_(aStep)
+    {}
+
+    /*!
+     * \brief Set the geometry of elliptic cylinder
+     * \param bottomCenter bottom center position
+     * \param topCenter top center position
+     * \param firstAxis first ellipse axis
+     * \param secondAxis second ellipse axis
+     * \param verbosity the verbosity level (default: 0 -> no terminal output)
+     */
+    void setGeometry(const GlobalPosition& bottomCenter,
+                     const GlobalPosition& topCenter,
+                     GlobalPosition firstAxis,
+                     GlobalPosition secondAxis,
+                     int verbosity = 0)
+    {
+        auto zAxis = (topCenter-bottomCenter);
+        const auto length = zAxis.two_norm();
+        zAxis /= length;
+        const auto a = firstAxis.two_norm();
+        const auto b = secondAxis.two_norm();
+        firstAxis /= a;
+        secondAxis /= b;
+        const auto normal = crossProduct(firstAxis, secondAxis);
+        const auto height = (zAxis*normal)*length;
+
+        const auto aStep = aCharLength_*a;
+        using std::floor;
+        const std::size_t aSamples = std::max<std::size_t>(floor(2*a/aStep), 1);
+        const std::size_t bSamples = std::max<std::size_t>(floor(2*b/aStep), 1);
+        const std::size_t zSamples = std::max<std::size_t>(floor(height/aStep), 1);
+        const auto bStep = 2*b / Scalar(bSamples);
+        const auto zStep = length / Scalar(zSamples);
+        const auto zStepFactor = length/height;
+
+        integrationElement_ = aStep*bStep*zStep/zStepFactor;
+        points_.clear();
+        points_.reserve(aSamples*bSamples*zSamples);
+
+        // walk over lattice grid and reject points outside the ellipse
+        auto startZ = bottomCenter; startZ.axpy(0.5*zStep, zAxis);
+        auto startAB = startZ;
+        startAB.axpy(-a + 0.5*aStep, firstAxis);
+        startAB.axpy(-b + 0.5*bStep, secondAxis);
+        for (std::size_t as = 0; as < aSamples; ++as)
+        {
+            auto posA = startAB;
+            posA.axpy(as*aStep, firstAxis);
+            for (std::size_t bs = 0; bs < bSamples; ++bs)
+            {
+                auto pos = posA;
+                pos.axpy(bs*bStep, secondAxis);
+                if (Impl::pointInEllipse(pos, startZ, firstAxis, secondAxis, normal, a, b))
+                    points_.emplace_back(std::move(pos));
+            }
+        }
+
+        const auto abPoints = points_.size();
+        numPoints_ = abPoints*zSamples;
+        points_.resize(numPoints_);
+
+        for (std::size_t zs = 1; zs < zSamples; ++zs)
+        {
+            auto zOffset = zAxis; zOffset *= zs*zStep;
+            std::transform(points_.begin(), points_.begin() + abPoints, points_.begin() + zs*abPoints,
+                           [&zOffset](const auto& p){ return p + zOffset; });
+        }
+
+        // error correction to obtain exact integral of constant functions
+        const auto meanLocalError = (height*a*b*M_PI - integrationElement_*numPoints_)/Scalar(numPoints_);
+        integrationElement_ += meanLocalError;
+
+        if (verbosity > 0)
+        {
+            std::cout << "EllipticCylinderIntegration -- number of integration points: " << numPoints_ <<"\n";
+            if (verbosity > 1)
+                std::cout << "  -- a: " << a << ", b: " << b << ", h: " << height << "\n"
+                          << "  -- aSamples: " << aSamples << ", bSamples: " << bSamples << ", zSamples: " << zSamples << "\n"
+                          << "  -- volume: " << integrationElement_*numPoints_ << " (expected " << height*a*b*M_PI << ")\n"
+                          << "  -- corrected a volume error of: " << meanLocalError/integrationElement_*100 << "%\n";
+        }
+    }
+
+    //! obtain ith integration element
+    Scalar integrationElement(std::size_t i) const
+    { return integrationElement_; }
+
+    //! obtain ith integration point
+    const GlobalPosition& integrationPoint(std::size_t i) const
+    { return points_[i]; }
+
+    //! number of samples points
+    std::size_t size() const
+    { return numPoints_; }
+
+private:
+    const Scalar aCharLength_;
+    std::size_t numPoints_;
+    Scalar integrationElement_;
+    std::vector<GlobalPosition> points_;
+};
+
+/*!
+ * \ingroup EmbeddedCoupling
+ * \brief Helper class to integrate over an elliptic domain
+ * \note This is mostly useful if the integral is known and the
+ *       integral mass has to be locally distributed to some non-matching domain
+ * \note The algorithm is based on evenly spaced area elements with rejection sampling.
+ *       To improve the quality of the integral, increase the sample size.
+ * \note The area integral of a constant function is exact
+ * See Koch et al. (2020) https://doi.org/10.1016/j.jcp.2020.109369
+ */
+template<class Scalar>
+class EllipseIntegration
+{
+    using GlobalPosition = Dune::FieldVector<Scalar, 3>;
+
+public:
+    /*!
+     * \brief Constructor
+     * \param aStep characteristic relative integration length (real number between 0 and 1)
+     * \note half the characteristic length means 2*2=4 times more integration points
+     */
+    explicit EllipseIntegration(const Scalar aStep)
+    : aCharLength_(aStep)
+    {}
+
+    /*!
+     * \brief set geometry of an ellipse
+     * \param center the center position
+     * \param firstAxis first ellipse axis
+     * \param secondAxis second ellipse axis
+     * \param verbosity the verbosity level (default: 0 -> no terminal output)
+     */
+    void setGeometry(const GlobalPosition& center,
+                     GlobalPosition firstAxis,
+                     GlobalPosition secondAxis,
+                     int verbosity = 0)
+    {
+        const auto a = firstAxis.two_norm();
+        const auto b = secondAxis.two_norm();
+        firstAxis /= a;
+        secondAxis /= b;
+        const auto normal = crossProduct(firstAxis, secondAxis);
+
+        const auto aStep = aCharLength_*a;
+        using std::floor;
+        const std::size_t aSamples = std::max<std::size_t>(floor(2*a/aStep), 1);
+        const std::size_t bSamples = std::max<std::size_t>(floor(2*b/aStep), 1);
+        const auto bStep = 2*b / Scalar(bSamples);
+
+        integrationElement_ = aStep*bStep;
+        points_.reserve(aSamples*bSamples);
+
+        // walk over lattice grid an reject point outside the ellipse
+        auto startZ = center;
+        auto startAB = startZ;
+        startAB.axpy(-a + 0.5*aStep, firstAxis);
+        startAB.axpy(-b + 0.5*bStep, secondAxis);
+        for (std::size_t as = 0; as < aSamples; ++as)
+        {
+            auto posA = startAB;
+            posA.axpy(as*aStep, firstAxis);
+            for (std::size_t bs = 0; bs < bSamples; ++bs)
+            {
+                auto pos = posA;
+                pos.axpy(bs*bStep, secondAxis);
+                if (Impl::pointInEllipse(pos, startZ, firstAxis, secondAxis, normal, a, b))
+                    points_.emplace_back(std::move(pos));
+            }
+        }
+
+        // store number of sample points
+        const auto abPoints = points_.size();
+        numPoints_ = abPoints;
+
+        // error correction to obtain exact integral of constant functions
+        const auto meanLocalError = (a*b*M_PI - integrationElement_*numPoints_)/Scalar(numPoints_);
+        integrationElement_ += meanLocalError;
+
+        if (verbosity > 0)
+        {
+            std::cout << "EllipseIntegration -- number of integration points: " << numPoints_ << "\n";
+            if (verbosity > 1)
+                std::cout << "  -- a: " << a << ", b: " << b << "\n"
+                          << "  -- aSamples: " << aSamples << ", bSamples: " << bSamples << "\n"
+                          << "  -- area: " << integrationElement_*numPoints_ << " (expected " << a*b*M_PI << ")\n"
+                          << "  -- corrected an area error of: " << meanLocalError/integrationElement_*100 << "%\n";
+        }
+    }
+
+    //! obtain ith integration element
+    Scalar integrationElement(std::size_t i) const
+    { return integrationElement_; }
+
+    //! obtain ith integration point
+    const GlobalPosition& integrationPoint(std::size_t i) const
+    { return points_[i]; }
+
+    //! number of integration points
+    std::size_t size() const
+    { return numPoints_; }
+
+private:
+    const Scalar aCharLength_;
+    std::size_t numPoints_;
+    Scalar integrationElement_;
+    std::vector<GlobalPosition> points_;
+};
+
+} // end namespace Dumux::EmbeddedCoupling
+
+#endif
-- 
GitLab


From d729771661219a50660d23d66d7f6195a0541b76 Mon Sep 17 00:00:00 2001
From: Martin Schneider <martin.schneider@iws.uni-stuttgart.de>
Date: Tue, 16 Jun 2020 15:35:23 +0200
Subject: [PATCH 2/3] [md][embedded][geometry] ellipse point calculation in
 separate function

---
 .../embedded/cylinderintegration.hh           | 184 +++++++++---------
 1 file changed, 97 insertions(+), 87 deletions(-)

diff --git a/dumux/multidomain/embedded/cylinderintegration.hh b/dumux/multidomain/embedded/cylinderintegration.hh
index e54fd5515e..dbcb901a4f 100644
--- a/dumux/multidomain/embedded/cylinderintegration.hh
+++ b/dumux/multidomain/embedded/cylinderintegration.hh
@@ -30,6 +30,9 @@
 #include <cmath>
 #include <algorithm>
 #include <numeric>
+#include <iterator>
+#include <utility>
+#include <tuple>
 
 #include <dune/common/fvector.hh>
 
@@ -227,6 +230,49 @@ inline bool pointInEllipse(const GlobalPosition& p,
 
     return (da*da/(a*a) + db*db/(b*b) < 1.0);
 }
+
+//! construct evenly distributed integration points on an ellipse
+template<class GlobalPosition>
+inline std::pair<std::vector<GlobalPosition>, typename GlobalPosition::value_type>
+ellipseIntegrationPoints(const GlobalPosition& center,
+                         const GlobalPosition& firstUnitAxis,
+                         const GlobalPosition& secondUnitAxis,
+                         typename GlobalPosition::value_type a,
+                         typename GlobalPosition::value_type b,
+                         const GlobalPosition& normal,
+                         typename GlobalPosition::value_type characteristicLength)
+{
+    // choose step size in a-direction
+    const auto aStep = characteristicLength;
+
+    using std::floor;
+    const std::size_t aSamples = std::max<std::size_t>(floor(2*a/aStep), 1);
+    const std::size_t bSamples = std::max<std::size_t>(floor(2*b/aStep), 1);
+    const auto bStep = 2*b / bSamples;
+
+    // reserve upper limit estimate for memory needed
+    std::vector<GlobalPosition> points;
+    points.reserve(aSamples*bSamples);
+
+    // walk over lattice grid and reject points outside the ellipse
+    auto startAB = center;
+    startAB.axpy(-a + 0.5*aStep, firstUnitAxis);
+    startAB.axpy(-b + 0.5*bStep, secondUnitAxis);
+    for (std::size_t as = 0; as < aSamples; ++as)
+    {
+        auto posA = startAB;
+        posA.axpy(as*aStep, firstUnitAxis);
+        for (std::size_t bs = 0; bs < bSamples; ++bs)
+        {
+            auto pos = posA;
+            pos.axpy(bs*bStep, secondUnitAxis);
+            if (pointInEllipse(pos, center, firstUnitAxis, secondUnitAxis, normal, a, b))
+                points.emplace_back(std::move(pos));
+        }
+    }
+    // return points and integration element
+    return {points, aStep*bStep};
+}
 } // end namespace Impl
 
 /*!
@@ -248,80 +294,67 @@ class EllipticCylinderIntegration
 public:
     /*!
      * \brief Constructor
-     * \param aStep characteristic relative integration length (real number between 0 and 1)
+     * \param relCharLength characteristic relative integration length (real number between 0 and 1)
      * \note half the characteristic length means 2*2*2=8 times more integration points
      */
-    explicit EllipticCylinderIntegration(const Scalar aStep)
-    : aCharLength_(aStep)
+    explicit EllipticCylinderIntegration(const Scalar relCharLength)
+    : relCharLength_(relCharLength)
     {}
 
     /*!
      * \brief Set the geometry of elliptic cylinder
      * \param bottomCenter bottom center position
      * \param topCenter top center position
-     * \param firstAxis first ellipse axis
-     * \param secondAxis second ellipse axis
+     * \param firstAxis first ellipse axis (length corresponding to axis length)
+     * \param secondAxis second ellipse axis (length corresponding to axis length)
      * \param verbosity the verbosity level (default: 0 -> no terminal output)
      */
     void setGeometry(const GlobalPosition& bottomCenter,
                      const GlobalPosition& topCenter,
-                     GlobalPosition firstAxis,
-                     GlobalPosition secondAxis,
+                     const GlobalPosition& firstAxis,
+                     const GlobalPosition& secondAxis,
                      int verbosity = 0)
     {
+        const auto a = firstAxis.two_norm();
+        const auto b = secondAxis.two_norm();
+        const auto characteristicLength = relCharLength_*a;
+
+        auto firstUnitAxis = firstAxis; firstUnitAxis /= a;
+        auto secondUnitAxis = secondAxis; secondUnitAxis /= b;
+        const auto normal = crossProduct(firstUnitAxis, secondUnitAxis);
+
         auto zAxis = (topCenter-bottomCenter);
         const auto length = zAxis.two_norm();
         zAxis /= length;
-        const auto a = firstAxis.two_norm();
-        const auto b = secondAxis.two_norm();
-        firstAxis /= a;
-        secondAxis /= b;
-        const auto normal = crossProduct(firstAxis, secondAxis);
         const auto height = (zAxis*normal)*length;
-
-        const auto aStep = aCharLength_*a;
         using std::floor;
-        const std::size_t aSamples = std::max<std::size_t>(floor(2*a/aStep), 1);
-        const std::size_t bSamples = std::max<std::size_t>(floor(2*b/aStep), 1);
-        const std::size_t zSamples = std::max<std::size_t>(floor(height/aStep), 1);
-        const auto bStep = 2*b / Scalar(bSamples);
+        const std::size_t zSamples = std::max<std::size_t>(floor(height/characteristicLength), 1);
         const auto zStep = length / Scalar(zSamples);
         const auto zStepFactor = length/height;
 
-        integrationElement_ = aStep*bStep*zStep/zStepFactor;
-        points_.clear();
-        points_.reserve(aSamples*bSamples*zSamples);
-
         // walk over lattice grid and reject points outside the ellipse
-        auto startZ = bottomCenter; startZ.axpy(0.5*zStep, zAxis);
-        auto startAB = startZ;
-        startAB.axpy(-a + 0.5*aStep, firstAxis);
-        startAB.axpy(-b + 0.5*bStep, secondAxis);
-        for (std::size_t as = 0; as < aSamples; ++as)
-        {
-            auto posA = startAB;
-            posA.axpy(as*aStep, firstAxis);
-            for (std::size_t bs = 0; bs < bSamples; ++bs)
-            {
-                auto pos = posA;
-                pos.axpy(bs*bStep, secondAxis);
-                if (Impl::pointInEllipse(pos, startZ, firstAxis, secondAxis, normal, a, b))
-                    points_.emplace_back(std::move(pos));
-            }
-        }
+        auto startZ = bottomCenter;
+        startZ.axpy(0.5*zStep, zAxis);
+        auto [ellipseIPs, ellipseIntegrationElement]
+            = Impl::ellipseIntegrationPoints(startZ, firstUnitAxis, secondUnitAxis, a, b, normal, characteristicLength);
 
-        const auto abPoints = points_.size();
+        // compute total number of points
+        const auto abPoints = ellipseIPs.size();
         numPoints_ = abPoints*zSamples;
-        points_.resize(numPoints_);
 
+        // first move in the first layer of ellipse integration points
+        points_.clear();
+        points_.reserve(numPoints_);
+        std::move(ellipseIPs.begin(), ellipseIPs.end(), std::back_inserter(points_));
+
+        // then extrude along the z axis
+        auto zStepVector = zAxis; zStepVector *= zStep;
         for (std::size_t zs = 1; zs < zSamples; ++zs)
-        {
-            auto zOffset = zAxis; zOffset *= zs*zStep;
-            std::transform(points_.begin(), points_.begin() + abPoints, points_.begin() + zs*abPoints,
-                           [&zOffset](const auto& p){ return p + zOffset; });
-        }
+            std::transform(points_.end() - abPoints, points_.end(), std::back_inserter(points_),
+                           [&](const auto& p){ return p + zStepVector; });
 
-        // error correction to obtain exact integral of constant functions
+        // computation and error correction for integral element to obtain exact integral of constant functions
+        integrationElement_ = ellipseIntegrationElement*zStep/zStepFactor;
         const auto meanLocalError = (height*a*b*M_PI - integrationElement_*numPoints_)/Scalar(numPoints_);
         integrationElement_ += meanLocalError;
 
@@ -330,7 +363,6 @@ public:
             std::cout << "EllipticCylinderIntegration -- number of integration points: " << numPoints_ <<"\n";
             if (verbosity > 1)
                 std::cout << "  -- a: " << a << ", b: " << b << ", h: " << height << "\n"
-                          << "  -- aSamples: " << aSamples << ", bSamples: " << bSamples << ", zSamples: " << zSamples << "\n"
                           << "  -- volume: " << integrationElement_*numPoints_ << " (expected " << height*a*b*M_PI << ")\n"
                           << "  -- corrected a volume error of: " << meanLocalError/integrationElement_*100 << "%\n";
         }
@@ -349,7 +381,7 @@ public:
     { return numPoints_; }
 
 private:
-    const Scalar aCharLength_;
+    const Scalar relCharLength_;
     std::size_t numPoints_;
     Scalar integrationElement_;
     std::vector<GlobalPosition> points_;
@@ -373,63 +405,42 @@ class EllipseIntegration
 public:
     /*!
      * \brief Constructor
-     * \param aStep characteristic relative integration length (real number between 0 and 1)
+     * \param relCharLength characteristic relative integration length (real number between 0 and 1)
      * \note half the characteristic length means 2*2=4 times more integration points
      */
-    explicit EllipseIntegration(const Scalar aStep)
-    : aCharLength_(aStep)
+    explicit EllipseIntegration(const Scalar relCharLength)
+    : relCharLength_(relCharLength)
     {}
 
     /*!
      * \brief set geometry of an ellipse
      * \param center the center position
-     * \param firstAxis first ellipse axis
-     * \param secondAxis second ellipse axis
+     * \param firstAxis first ellipse axis (length corresponding to axis length)
+     * \param secondAxis second ellipse axis (length corresponding to axis length)
      * \param verbosity the verbosity level (default: 0 -> no terminal output)
      */
     void setGeometry(const GlobalPosition& center,
-                     GlobalPosition firstAxis,
-                     GlobalPosition secondAxis,
+                     const GlobalPosition& firstAxis,
+                     const GlobalPosition& secondAxis,
                      int verbosity = 0)
     {
         const auto a = firstAxis.two_norm();
         const auto b = secondAxis.two_norm();
-        firstAxis /= a;
-        secondAxis /= b;
-        const auto normal = crossProduct(firstAxis, secondAxis);
+        const auto characteristicLength = relCharLength_*a;
 
-        const auto aStep = aCharLength_*a;
-        using std::floor;
-        const std::size_t aSamples = std::max<std::size_t>(floor(2*a/aStep), 1);
-        const std::size_t bSamples = std::max<std::size_t>(floor(2*b/aStep), 1);
-        const auto bStep = 2*b / Scalar(bSamples);
-
-        integrationElement_ = aStep*bStep;
-        points_.reserve(aSamples*bSamples);
-
-        // walk over lattice grid an reject point outside the ellipse
-        auto startZ = center;
-        auto startAB = startZ;
-        startAB.axpy(-a + 0.5*aStep, firstAxis);
-        startAB.axpy(-b + 0.5*bStep, secondAxis);
-        for (std::size_t as = 0; as < aSamples; ++as)
-        {
-            auto posA = startAB;
-            posA.axpy(as*aStep, firstAxis);
-            for (std::size_t bs = 0; bs < bSamples; ++bs)
-            {
-                auto pos = posA;
-                pos.axpy(bs*bStep, secondAxis);
-                if (Impl::pointInEllipse(pos, startZ, firstAxis, secondAxis, normal, a, b))
-                    points_.emplace_back(std::move(pos));
-            }
-        }
+        auto firstUnitAxis = firstAxis; firstUnitAxis /= a;
+        auto secondUnitAxis = secondAxis; secondUnitAxis /= b;
+        const auto normal = crossProduct(firstUnitAxis, secondUnitAxis);
+
+        // generate integration points
+        std::tie(points_, integrationElement_)
+            = Impl::ellipseIntegrationPoints(center, firstUnitAxis, secondUnitAxis, a, b, normal, characteristicLength);
 
         // store number of sample points
         const auto abPoints = points_.size();
         numPoints_ = abPoints;
 
-        // error correction to obtain exact integral of constant functions
+        // computation and error correction for integral element to obtain exact integral of constant functions
         const auto meanLocalError = (a*b*M_PI - integrationElement_*numPoints_)/Scalar(numPoints_);
         integrationElement_ += meanLocalError;
 
@@ -438,7 +449,6 @@ public:
             std::cout << "EllipseIntegration -- number of integration points: " << numPoints_ << "\n";
             if (verbosity > 1)
                 std::cout << "  -- a: " << a << ", b: " << b << "\n"
-                          << "  -- aSamples: " << aSamples << ", bSamples: " << bSamples << "\n"
                           << "  -- area: " << integrationElement_*numPoints_ << " (expected " << a*b*M_PI << ")\n"
                           << "  -- corrected an area error of: " << meanLocalError/integrationElement_*100 << "%\n";
         }
@@ -457,7 +467,7 @@ public:
     { return numPoints_; }
 
 private:
-    const Scalar aCharLength_;
+    const Scalar relCharLength_;
     std::size_t numPoints_;
     Scalar integrationElement_;
     std::vector<GlobalPosition> points_;
-- 
GitLab


From ec7e69538d0944d16b87537caff6c15b3cca8b16 Mon Sep 17 00:00:00 2001
From: Timo Koch <timo.koch@iws.uni-stuttgart.de>
Date: Fri, 12 Jun 2020 14:50:23 +0200
Subject: [PATCH 3/3] [test] Add test for cylinder integration

---
 test/common/geometry/CMakeLists.txt           |   1 +
 .../geometry/test_cylinderintegration.cc      | 168 ++++++++++++++++++
 2 files changed, 169 insertions(+)
 create mode 100644 test/common/geometry/test_cylinderintegration.cc

diff --git a/test/common/geometry/CMakeLists.txt b/test/common/geometry/CMakeLists.txt
index 92ea708139..f1a6d093ce 100644
--- a/test/common/geometry/CMakeLists.txt
+++ b/test/common/geometry/CMakeLists.txt
@@ -9,6 +9,7 @@ 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)
+dumux_add_test(SOURCES test_cylinderintegration.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_cylinderintegration.cc b/test/common/geometry/test_cylinderintegration.cc
new file mode 100644
index 0000000000..f756c24dbb
--- /dev/null
+++ b/test/common/geometry/test_cylinderintegration.cc
@@ -0,0 +1,168 @@
+#include <config.h>
+
+#include <iostream>
+#include <iomanip>
+#include <cmath>
+
+#include <dune/common/float_cmp.hh>
+#include <dune/common/fvector.hh>
+#include <dune/common/exceptions.hh>
+
+#include <dumux/common/math.hh>
+#include <dumux/multidomain/embedded/cylinderintegration.hh>
+#include "transformation.hh"
+
+namespace Dumux {
+
+template<class Scalar>
+void checkVolumeIntegral(const Scalar volIntegralComputed, const Scalar volIntegralExact)
+{
+    std::cout << std::setprecision(7) << std::scientific
+              << "Volume integral: " << volIntegralComputed << " (computed), "
+                                     << volIntegralExact << " (exact)" << "\n";
+    if (Dune::FloatCmp::ne(volIntegralComputed, volIntegralExact, 1e-7))
+        DUNE_THROW(Dune::Exception, "Volume integral is not exact!");
+}
+
+template<class Scalar>
+void checkAreaIntegral(const Scalar areaIntegralComputed, const Scalar areaIntegralExact)
+{
+    std::cout << std::setprecision(7) << std::scientific
+              << "Area integral: " << areaIntegralComputed << " (computed), "
+                                   << areaIntegralExact << " (exact)" << "\n";
+    if (Dune::FloatCmp::ne(areaIntegralComputed, areaIntegralExact, 1e-7))
+        DUNE_THROW(Dune::Exception, "Area integral is not exact!");
+}
+
+template<class GlobalPosition>
+void checkCentroid(const GlobalPosition& centroidComputed, const GlobalPosition& centroidExact)
+{
+    std::cout << std::setprecision(7) << std::scientific
+              << "Centroid: " << centroidComputed << " (computed), "
+                              << centroidExact << " (exact)" << "\n";
+    if (Dune::FloatCmp::ne(centroidComputed, centroidExact, 1e-7))
+        DUNE_THROW(Dune::Exception, "Centroid is wrong!");
+}
+
+template<class TestData>
+void checkCylinderIntegration(const TestData& data, const double characteristicLength)
+{
+    using Point = typename TestData::Point;
+
+    std::cout << std::endl;
+    EmbeddedCoupling::CylinderIntegration<double> cylinderIntegration(characteristicLength);
+    cylinderIntegration.setGeometry(data.bottomCenter, data.topCenter, data.radius, /*verbosity=*/1);
+    double volumeIntegral = 0.0;
+    Point centroid({0.0, 0.0, 0.0});
+    for (std::size_t i = 0; i < cylinderIntegration.size(); ++i)
+    {
+        volumeIntegral += cylinderIntegration.integrationElement(i);
+        centroid.axpy(cylinderIntegration.integrationElement(i), cylinderIntegration.integrationPoint(i));
+    }
+    centroid /= volumeIntegral;
+
+    checkVolumeIntegral(volumeIntegral, data.exactVolumeIntegral);
+    checkCentroid(centroid, data.exactCentroid);
+}
+
+template<class TestData, class EllipseData>
+void checkEllipticCylinderIntegration(const TestData& data, const EllipseData& ell, const double characteristicLength)
+{
+    using Point = typename TestData::Point;
+
+    std::cout << std::endl;
+    EmbeddedCoupling::EllipticCylinderIntegration<double> ellCylIntegration(characteristicLength);
+    ellCylIntegration.setGeometry(data.bottomCenter, data.topCenter, ell.firstAxis, ell.secondAxis, /*verbosity=*/1);
+    double volumeIntegral = 0.0;
+    Point centroid({0.0, 0.0, 0.0});
+    for (std::size_t i = 0; i < ellCylIntegration.size(); ++i)
+    {
+        volumeIntegral += ellCylIntegration.integrationElement(i);
+        centroid.axpy(ellCylIntegration.integrationElement(i), ellCylIntegration.integrationPoint(i));
+    }
+    centroid /= volumeIntegral;
+
+    checkVolumeIntegral(volumeIntegral, data.exactVolumeIntegral);
+    checkCentroid(centroid, data.exactCentroid);
+
+    // check top cap ellipse
+    std::cout << std::endl;
+    EmbeddedCoupling::EllipseIntegration<double> ellIntegration(characteristicLength);
+    ellIntegration.setGeometry(data.topCenter, ell.firstAxis, ell.secondAxis, /*verbosity=*/1);
+    double areaIntegral = 0.0;
+    Point topCentroid({0.0, 0.0, 0.0});
+    for (std::size_t i = 0; i < ellIntegration.size(); ++i)
+    {
+        areaIntegral += ellIntegration.integrationElement(i);
+        topCentroid.axpy(ellIntegration.integrationElement(i), ellIntegration.integrationPoint(i));
+    }
+    topCentroid /= areaIntegral;
+
+    checkAreaIntegral(areaIntegral, ell.firstAxis.two_norm()*ell.secondAxis.two_norm()*M_PI);
+    checkCentroid(topCentroid, data.topCenter);
+}
+
+} // end namespace Dumux
+
+int main(int argc, char* argv[])
+{
+    using namespace Dumux;
+    using Point = Dune::FieldVector<double, 3>;
+
+    struct TestData {
+        using Point = Dune::FieldVector<double, 3>;
+        Point bottomCenter, topCenter;
+        double radius;
+        double exactVolumeIntegral;
+        Point exactCentroid;
+    };
+
+    TestData data;
+    data.bottomCenter = Point({1.2, 2.3, 4.5});
+    data.topCenter = Point({8.2, 4.3, 6.5});
+    data.radius = 4.0;
+    data.exactVolumeIntegral = M_PI*data.radius*data.radius*(data.topCenter-data.bottomCenter).two_norm();
+    data.exactCentroid = data.bottomCenter; data.exactCentroid.axpy(0.5, data.topCenter-data.bottomCenter);
+
+    /////////////////////////////////////////////////////////////////////
+    // Determine cylinder volume and centroid by integration
+    ////////////////////////////////////////////////////////////////////
+    checkCylinderIntegration(data, 0.05);
+
+    /////////////////////////////////////////////////////////////////////
+    // Determine elliptic cylinder volume and centroid by integration
+    // Construct an perfect cylinder (a == b == radius and orthogonal caps)
+    ////////////////////////////////////////////////////////////////////
+    struct EllipseData { Point firstAxis, secondAxis; };
+    auto topNormal = data.topCenter-data.bottomCenter; topNormal /= topNormal.two_norm();
+    auto firstAxis = crossProduct(topNormal, Point({1.0, 0.0, 0.0}));
+    firstAxis /= firstAxis.two_norm();
+    auto secondAxis = crossProduct(topNormal, firstAxis);
+    // scale both axis with the same radius -> circle
+    firstAxis *= data.radius; secondAxis *= data.radius;
+    checkEllipticCylinderIntegration(data, EllipseData{firstAxis, secondAxis}, 0.05);
+
+    /////////////////////////////////////////////////////////////////////
+    // Determine elliptic cylinder volume and centroid by integration
+    // Construct an elliptic cylinder (a == 4*b and orthogonal caps)
+    ////////////////////////////////////////////////////////////////////
+
+    firstAxis *= 2; secondAxis *= 0.5;
+    checkEllipticCylinderIntegration(data, EllipseData{firstAxis, secondAxis}, 0.025);
+
+    /////////////////////////////////////////////////////////////////////
+    // Determine elliptic cylinder volume and centroid by integration
+    // Construct an elliptic cylinder (a == 2*b and tilted caps)
+    ////////////////////////////////////////////////////////////////////
+
+    // rotate about the first major axis
+    const auto rotationAngle = M_PI*0.25;
+    auto rotationAxis = firstAxis; rotationAxis /= rotationAxis.two_norm();
+    const auto rotate = make3DTransformation(1.0, Point(0.0), rotationAxis, rotationAngle, /*verbose*/false);
+    secondAxis = rotate(secondAxis);
+    // scale second axis so that the volume stays the same as for the orthonal cap case
+    secondAxis /= std::cos(rotationAngle);
+    checkEllipticCylinderIntegration(data, EllipseData{firstAxis, secondAxis}, 0.025);
+
+    return 0;
+}
-- 
GitLab