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