| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- | ||
| 2 | // vi: set et ts=4 sw=4 sts=4: | ||
| 3 | // | ||
| 4 | // SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder | ||
| 5 | // SPDX-License-Identifier: GPL-3.0-or-later | ||
| 6 | // | ||
| 7 | /*! | ||
| 8 | * \file | ||
| 9 | * \ingroup EmbeddedCoupling | ||
| 10 | * \brief Integration over cylindrical and elliptic cylindrical domains | ||
| 11 | * Lowest order integration formulas that are mostly useful to evenly distribute | ||
| 12 | * mass in a cylindrical domain. | ||
| 13 | */ | ||
| 14 | |||
| 15 | #ifndef DUMUX_MULTIDOMAIN_EMBEDDED_CYLINDERINTEGRATOR_HH | ||
| 16 | #define DUMUX_MULTIDOMAIN_EMBEDDED_CYLINDERINTEGRATOR_HH | ||
| 17 | |||
| 18 | #include <cmath> | ||
| 19 | #include <algorithm> | ||
| 20 | #include <numeric> | ||
| 21 | #include <iterator> | ||
| 22 | #include <utility> | ||
| 23 | #include <tuple> | ||
| 24 | #include <limits> | ||
| 25 | |||
| 26 | #include <dune/common/fvector.hh> | ||
| 27 | |||
| 28 | #include <dumux/common/math.hh> | ||
| 29 | #include <dumux/multidomain/embedded/circlepoints.hh> | ||
| 30 | |||
| 31 | namespace Dumux::EmbeddedCoupling { | ||
| 32 | |||
| 33 | /*! | ||
| 34 | * \ingroup EmbeddedCoupling | ||
| 35 | * \brief Helper class to integrate over a cylinder domain | ||
| 36 | * \note This is mostly useful if the integral is known and the | ||
| 37 | * integral mass has to be locally distributed to some non-matching domain | ||
| 38 | * \note The algorithm creates (almost) evenly spaced area elements on a circle | ||
| 39 | * and extrudes this pattern into a cylinder. Each sub-cell is an integration | ||
| 40 | * point with constant ansatz on the sub-cell (lowest order integration formula). | ||
| 41 | * Hence to improve the quality of the integral, increase the sample size. | ||
| 42 | * \note The volume integral of a constant function is exact | ||
| 43 | * See Beckers & Beckers (2012) doi:10.1016/j.comgeo.2012.01.011 for circle distribution | ||
| 44 | * and T. Koch PhD thesis (2020), Section 7.2.4. | ||
| 45 | */ | ||
| 46 | template<class Scalar> | ||
| 47 | class CylinderIntegration | ||
| 48 | { | ||
| 49 | using GlobalPosition = Dune::FieldVector<Scalar, 3>; | ||
| 50 | |||
| 51 | public: | ||
| 52 | /*! | ||
| 53 | * \brief Constructor | ||
| 54 | * \param rStep characteristic relative integration length (positive real number between 0 and 1) | ||
| 55 | * \note half the characteristic length means 2*2*2=8 times more integration points | ||
| 56 | */ | ||
| 57 | 1 | CylinderIntegration(const Scalar rStep) | |
| 58 | 4 | { | |
| 59 | 1 | rSamples_ = characteristicLengthToNumSamples_(rStep); | |
| 60 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | initialize_(); |
| 61 | 1 | } | |
| 62 | |||
| 63 | /*! | ||
| 64 | * \brief Constructor | ||
| 65 | * \param rStep characteristic relative integration length in r-direction (positive real number between 0 and 1) | ||
| 66 | * \param zStep characteristic relative integration length in z-direction (positive real number between 0 and 1) | ||
| 67 | * \note Use this constructor to achieve a non-balanced (away from 1) aspect ratio between r and z-direction | ||
| 68 | */ | ||
| 69 | 6 | CylinderIntegration(const Scalar rStep, const Scalar zStep) | |
| 70 |
0/6✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
24 | : zStepFixed_(true) |
| 71 | { | ||
| 72 | using std::ceil; using std::clamp; | ||
| 73 | 6 | rSamples_ = characteristicLengthToNumSamples_(rStep); | |
| 74 | 6 | zSamples_ = characteristicLengthToNumSamples_(zStep); | |
| 75 |
1/2✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
|
6 | initialize_(); |
| 76 | 6 | } | |
| 77 | |||
| 78 | /*! | ||
| 79 | * \brief Set the geometry of the cylinder | ||
| 80 | * \param bottomCenter bottom center position | ||
| 81 | * \param topCenter top center position | ||
| 82 | * \param radius cylinder radius | ||
| 83 | * \param verbosity the verbosity level (default: 0 -> no terminal output) | ||
| 84 | */ | ||
| 85 | 1425 | void setGeometry(const GlobalPosition& bottomCenter, | |
| 86 | const GlobalPosition& topCenter, | ||
| 87 | const Scalar radius, | ||
| 88 | int verbosity = 0) | ||
| 89 | { | ||
| 90 | 1425 | const auto zAxis = (topCenter-bottomCenter); | |
| 91 | 1425 | const auto height = zAxis.two_norm(); | |
| 92 | 1425 | const auto zAxisUnitVector = zAxis/height; | |
| 93 | |||
| 94 | // compute step size in r-direction | ||
| 95 | 1425 | const auto rStep = radius/Scalar(rSamples_); | |
| 96 | |||
| 97 | // compute zSamples from r samples if not specified by the user | ||
| 98 | using std::ceil; | ||
| 99 |
2/2✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1424 times.
|
1425 | if (!zStepFixed_) |
| 100 |
1/2✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
|
2 | zSamples_ = std::max<std::size_t>(1, ceil(height/rStep)); |
| 101 | 1425 | const auto zStep = height/Scalar(zSamples_); | |
| 102 | |||
| 103 | // compute offsets for index calculation | ||
| 104 |
0/2✗ Branch 1 not taken.
✗ Branch 2 not taken.
|
1425 | auto kOffset = numRingSamples_; |
| 105 |
4/8✓ Branch 0 taken 1425 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1425 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1425 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 1425 times.
✗ Branch 7 not taken.
|
5700 | std::partial_sum(kOffset.begin(), kOffset.end(), kOffset.begin()); |
| 106 | |||
| 107 | // compute total number of samples | ||
| 108 | 1425 | numPoints_ = zSamples_*circleSamples_; | |
| 109 |
2/2✓ Branch 0 taken 1 times.
✓ Branch 1 taken 1424 times.
|
1425 | if (verbosity > 0) |
| 110 |
2/4✓ Branch 2 taken 1 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 1 times.
✗ Branch 6 not taken.
|
2 | std::cout << "CylinderIntegration -- number of integration points: " << numPoints_ << std::endl; |
| 111 | |||
| 112 | // resize integration points and integration elements | ||
| 113 |
1/2✓ Branch 1 taken 1425 times.
✗ Branch 2 not taken.
|
1425 | integrationElement_.resize(numPoints_); |
| 114 |
1/2✓ Branch 1 taken 1425 times.
✗ Branch 2 not taken.
|
1425 | points_.resize(numPoints_); |
| 115 | |||
| 116 | // reserve enough memory for the points on each ring | ||
| 117 |
3/8✓ Branch 1 taken 1425 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1425 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 1425 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
4275 | std::vector<GlobalPosition> ringPoints; |
| 118 |
2/4✓ Branch 1 taken 1425 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1425 times.
✗ Branch 5 not taken.
|
2850 | ringPoints.reserve(numRingSamples_.back()); |
| 119 | |||
| 120 | // compute integration points and integration elements | ||
| 121 |
2/2✓ Branch 0 taken 1462 times.
✓ Branch 1 taken 1425 times.
|
2887 | for (std::size_t i = 0; i < zSamples_; ++i) |
| 122 | { | ||
| 123 | // for each cylinder slice i | ||
| 124 | 1462 | auto sliceCenter = bottomCenter; | |
| 125 | 1462 | sliceCenter.axpy((Scalar(i)+0.5)*zStep, zAxisUnitVector); | |
| 126 | |||
| 127 | // generate circle points for each ring j | ||
| 128 |
2/2✓ Branch 0 taken 23550 times.
✓ Branch 1 taken 1462 times.
|
25012 | for (std::size_t j = 0; j < rSamples_; ++j) |
| 129 | { | ||
| 130 | 23550 | const auto r = (Scalar(j)+0.5)*rStep; | |
| 131 |
2/4✓ Branch 1 taken 23550 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 23550 times.
✗ Branch 5 not taken.
|
47100 | circlePoints(ringPoints, sincos_[j], sliceCenter, zAxisUnitVector, r); |
| 132 |
2/2✓ Branch 0 taken 22088 times.
✓ Branch 1 taken 1462 times.
|
23550 | const auto ringOffest = j > 0 ? kOffset[j-1] : 0; |
| 133 | 23550 | const auto ringSamples = numRingSamples_[j]; | |
| 134 | // volume of each element in ring | ||
| 135 | // total ring volume is given by M_PI*rStep*rStep*zStep*((j+1)^2 - j^2) | ||
| 136 | 23550 | const auto integrationElement = M_PI*rStep*rStep*zStep*(1.0 + 2.0*Scalar(j))/ringSamples; | |
| 137 |
2/2✓ Branch 0 taken 3437684 times.
✓ Branch 1 taken 23550 times.
|
3461234 | for (std::size_t k = 0; k < ringSamples; ++k) |
| 138 | { | ||
| 139 | 3437684 | const std::size_t idx = k + ringOffest + circleSamples_*i; | |
| 140 | 6875368 | points_[idx] = ringPoints[k]; | |
| 141 | 6875368 | integrationElement_[idx] = integrationElement; | |
| 142 | } | ||
| 143 | } | ||
| 144 | } | ||
| 145 | 1425 | } | |
| 146 | |||
| 147 | //! The integration element of the ith integration point | ||
| 148 | Scalar integrationElement(std::size_t i) const | ||
| 149 |
4/4✓ Branch 0 taken 12325157 times.
✓ Branch 1 taken 3792 times.
✓ Branch 2 taken 12325157 times.
✓ Branch 3 taken 3792 times.
|
24799942 | { return integrationElement_[i]; } |
| 150 | |||
| 151 | //! The ith integration point | ||
| 152 | const GlobalPosition& integrationPoint(std::size_t i) const | ||
| 153 | 24753418 | { return points_[i]; } | |
| 154 | |||
| 155 | //! The number of integration points | ||
| 156 | ✗ | std::size_t size() const | |
| 157 | ✗ | { return numPoints_; } | |
| 158 | |||
| 159 | private: | ||
| 160 | 7 | void initialize_() | |
| 161 | { | ||
| 162 | // precompute number of cells in ring and total circle samples | ||
| 163 | 7 | circleSamples_ = 0; | |
| 164 | 7 | numRingSamples_.resize(rSamples_); | |
| 165 |
2/2✓ Branch 1 taken 530 times.
✓ Branch 2 taken 7 times.
|
537 | for (int j = 0; j < rSamples_; ++j) |
| 166 | { | ||
| 167 | // number of cells in ring j | ||
| 168 | using std::floor; | ||
| 169 | 530 | numRingSamples_[j] = floor(2*M_PI*(Scalar(j)+0.5)); | |
| 170 | 530 | circleSamples_ += numRingSamples_[j]; | |
| 171 | } | ||
| 172 | |||
| 173 | // optimization because calling too many sin/cos can be really expensive | ||
| 174 | 7 | sincos_.resize(rSamples_); | |
| 175 |
2/2✓ Branch 1 taken 530 times.
✓ Branch 2 taken 7 times.
|
537 | for (int j = 0; j < rSamples_; ++j) |
| 176 | { | ||
| 177 | 530 | const auto numPoints = numRingSamples_[j]; | |
| 178 | 530 | sincos_[j].resize(2*numPoints); | |
| 179 | 530 | Scalar t = 0 + 0.1; | |
| 180 |
2/2✓ Branch 1 taken 158380 times.
✓ Branch 2 taken 530 times.
|
158910 | for (std::size_t i = 0; i < numPoints; ++i) |
| 181 | { | ||
| 182 | using std::sin; using std::cos; | ||
| 183 |
4/4✓ Branch 0 taken 530 times.
✓ Branch 1 taken 157850 times.
✓ Branch 2 taken 530 times.
✓ Branch 3 taken 157850 times.
|
316760 | sincos_[j][2*i] = sin(t); |
| 184 |
2/2✓ Branch 0 taken 530 times.
✓ Branch 1 taken 157850 times.
|
158380 | sincos_[j][2*i + 1] = cos(t); |
| 185 | 158380 | t += 2*M_PI/numPoints; | |
| 186 |
2/2✓ Branch 0 taken 530 times.
✓ Branch 1 taken 157850 times.
|
158380 | if (t > 2*M_PI) t -= 2*M_PI; |
| 187 | } | ||
| 188 | } | ||
| 189 | 7 | } | |
| 190 | |||
| 191 | // invert characteristic length (between 0 and 1) | ||
| 192 | // and make sure there is no integer wrap around for small ratios | ||
| 193 | ✗ | std::size_t characteristicLengthToNumSamples_(const Scalar cl) const | |
| 194 | { | ||
| 195 | using std::clamp; using std::min; using std::max; using std::ceil; using std::nexttoward; | ||
| 196 | ✗ | const Scalar clampedCl = clamp(cl, std::numeric_limits<Scalar>::min(), 1.0); | |
| 197 | ✗ | const Scalar floatNumSamples = ceil(1.0/clampedCl); | |
| 198 | ✗ | const Scalar largestValidFloat = nexttoward(std::numeric_limits<std::size_t>::max(), 0.0); | |
| 199 | ✗ | return static_cast<std::size_t>(min(floatNumSamples, largestValidFloat)); | |
| 200 | } | ||
| 201 | |||
| 202 | bool zStepFixed_ = false; | ||
| 203 | std::size_t zSamples_, rSamples_, numPoints_, circleSamples_; | ||
| 204 | std::vector<Scalar> integrationElement_; | ||
| 205 | std::vector<GlobalPosition> points_; | ||
| 206 | std::vector<std::size_t> numRingSamples_; | ||
| 207 | std::vector<std::vector<Scalar>> sincos_; | ||
| 208 | }; | ||
| 209 | |||
| 210 | namespace Detail { | ||
| 211 | //! check if a point is in an ellipse | ||
| 212 | template<class GlobalPosition> | ||
| 213 | 10880 | inline bool pointInEllipse(const GlobalPosition& p, | |
| 214 | const GlobalPosition& center, | ||
| 215 | const GlobalPosition& firstAxis, | ||
| 216 | const GlobalPosition& secondAxis, | ||
| 217 | const GlobalPosition& normal, | ||
| 218 | const typename GlobalPosition::value_type a, | ||
| 219 | const typename GlobalPosition::value_type b) | ||
| 220 | { | ||
| 221 | 10880 | const auto d = p-center; | |
| 222 | // early return if point is not on the ellipse plane | ||
| 223 |
1/2✓ Branch 0 taken 10880 times.
✗ Branch 1 not taken.
|
10880 | if (d*normal > 1e-7*a) |
| 224 | return false; | ||
| 225 | |||
| 226 | // check if it's actually within the ellipse | ||
| 227 | const auto da = (d*firstAxis); | ||
| 228 | 10880 | const auto db = (d*secondAxis); | |
| 229 | |||
| 230 | 10880 | return (da*da/(a*a) + db*db/(b*b) < 1.0); | |
| 231 | } | ||
| 232 | |||
| 233 | //! construct evenly distributed integration points on an ellipse | ||
| 234 | template<class GlobalPosition> | ||
| 235 | inline std::pair<std::vector<GlobalPosition>, typename GlobalPosition::value_type> | ||
| 236 | 6 | ellipseIntegrationPoints(const GlobalPosition& center, | |
| 237 | const GlobalPosition& firstUnitAxis, | ||
| 238 | const GlobalPosition& secondUnitAxis, | ||
| 239 | typename GlobalPosition::value_type a, | ||
| 240 | typename GlobalPosition::value_type b, | ||
| 241 | const GlobalPosition& normal, | ||
| 242 | typename GlobalPosition::value_type characteristicLength) | ||
| 243 | { | ||
| 244 | // choose step size in a-direction | ||
| 245 | 6 | const auto aStep = characteristicLength; | |
| 246 | |||
| 247 | using std::floor; | ||
| 248 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
|
6 | const std::size_t aSamples = std::max<std::size_t>(floor(2*a/aStep), 1); |
| 249 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 6 times.
|
6 | const std::size_t bSamples = std::max<std::size_t>(floor(2*b/aStep), 1); |
| 250 | 6 | const auto bStep = 2*b / bSamples; | |
| 251 | |||
| 252 | // reserve upper limit estimate for memory needed | ||
| 253 |
1/4✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
6 | std::vector<GlobalPosition> points; |
| 254 |
1/2✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
|
6 | points.reserve(aSamples*bSamples); |
| 255 | |||
| 256 | // walk over lattice grid and reject points outside the ellipse | ||
| 257 | 6 | auto startAB = center; | |
| 258 | 6 | startAB.axpy(-a + 0.5*aStep, firstUnitAxis); | |
| 259 | 6 | startAB.axpy(-b + 0.5*bStep, secondUnitAxis); | |
| 260 |
2/2✓ Branch 0 taken 400 times.
✓ Branch 1 taken 6 times.
|
406 | for (std::size_t as = 0; as < aSamples; ++as) |
| 261 | { | ||
| 262 | 400 | auto posA = startAB; | |
| 263 | 400 | posA.axpy(as*aStep, firstUnitAxis); | |
| 264 |
2/2✓ Branch 0 taken 10880 times.
✓ Branch 1 taken 400 times.
|
11280 | for (std::size_t bs = 0; bs < bSamples; ++bs) |
| 265 | { | ||
| 266 | 10880 | auto pos = posA; | |
| 267 | 10880 | pos.axpy(bs*bStep, secondUnitAxis); | |
| 268 |
2/2✓ Branch 1 taken 8568 times.
✓ Branch 2 taken 2312 times.
|
10880 | if (pointInEllipse(pos, center, firstUnitAxis, secondUnitAxis, normal, a, b)) |
| 269 |
1/2✓ Branch 1 taken 8568 times.
✗ Branch 2 not taken.
|
8568 | points.emplace_back(std::move(pos)); |
| 270 | } | ||
| 271 | } | ||
| 272 | // return points and integration element | ||
| 273 |
2/4✓ Branch 1 taken 6 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 6 times.
✗ Branch 4 not taken.
|
12 | return {points, aStep*bStep}; |
| 274 | } | ||
| 275 | } // end namespace Detail | ||
| 276 | |||
| 277 | /*! | ||
| 278 | * \ingroup EmbeddedCoupling | ||
| 279 | * \brief Helper class to integrate over an elliptic cylinder domain | ||
| 280 | * \note The cylinder caps do not have to be orthogonal to the centerline axis but top and bottom cap are parallel planes | ||
| 281 | * \note This is mostly useful if the integral is known and the | ||
| 282 | * integral mass has to be locally distributed to some non-matching domain | ||
| 283 | * \note The algorithm is based on evenly spaced area elements with rejection sampling. | ||
| 284 | * To improve the quality of the integral, increase the sample size. | ||
| 285 | * \note The volume integral of a constant function is exact | ||
| 286 | * See Koch et al. (2020) https://doi.org/10.1016/j.jcp.2020.109369 | ||
| 287 | */ | ||
| 288 | template<class Scalar> | ||
| 289 |
1/4✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
6 | class EllipticCylinderIntegration |
| 290 | { | ||
| 291 | using GlobalPosition = Dune::FieldVector<Scalar, 3>; | ||
| 292 | |||
| 293 | public: | ||
| 294 | /*! | ||
| 295 | * \brief Constructor | ||
| 296 | * \param relCharLength characteristic relative integration length (real number between 0 and 1) | ||
| 297 | * \note half the characteristic length means 2*2*2=8 times more integration points | ||
| 298 | */ | ||
| 299 | 3 | explicit EllipticCylinderIntegration(const Scalar relCharLength) | |
| 300 |
2/4✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
|
6 | : relCharLength_(relCharLength) |
| 301 | {} | ||
| 302 | |||
| 303 | /*! | ||
| 304 | * \brief Set the geometry of elliptic cylinder | ||
| 305 | * \param bottomCenter bottom center position | ||
| 306 | * \param topCenter top center position | ||
| 307 | * \param firstAxis first ellipse axis (length corresponding to axis length) | ||
| 308 | * \param secondAxis second ellipse axis (length corresponding to axis length) | ||
| 309 | * \param verbosity the verbosity level (default: 0 -> no terminal output) | ||
| 310 | */ | ||
| 311 | 3 | void setGeometry(const GlobalPosition& bottomCenter, | |
| 312 | const GlobalPosition& topCenter, | ||
| 313 | const GlobalPosition& firstAxis, | ||
| 314 | const GlobalPosition& secondAxis, | ||
| 315 | int verbosity = 0) | ||
| 316 | { | ||
| 317 | 3 | const auto a = firstAxis.two_norm(); | |
| 318 | 3 | const auto b = secondAxis.two_norm(); | |
| 319 | 3 | const auto characteristicLength = relCharLength_*a; | |
| 320 | |||
| 321 | 3 | auto firstUnitAxis = firstAxis; firstUnitAxis /= a; | |
| 322 | 3 | auto secondUnitAxis = secondAxis; secondUnitAxis /= b; | |
| 323 | 3 | const auto normal = crossProduct(firstUnitAxis, secondUnitAxis); | |
| 324 | |||
| 325 | 3 | auto zAxis = (topCenter-bottomCenter); | |
| 326 | 3 | const auto length = zAxis.two_norm(); | |
| 327 | 3 | zAxis /= length; | |
| 328 | 3 | const auto height = (zAxis*normal)*length; | |
| 329 | using std::floor; | ||
| 330 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
|
3 | const std::size_t zSamples = std::max<std::size_t>(floor(height/characteristicLength), 1); |
| 331 | 3 | const auto zStep = length / Scalar(zSamples); | |
| 332 | 3 | const auto zStepFactor = length/height; | |
| 333 | |||
| 334 | // walk over lattice grid and reject points outside the ellipse | ||
| 335 | 3 | auto startZ = bottomCenter; | |
| 336 | 3 | startZ.axpy(0.5*zStep, zAxis); | |
| 337 |
3/8✗ Branch 1 not taken.
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 3 times.
✓ Branch 5 taken 3 times.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
6 | auto [ellipseIPs, ellipseIntegrationElement] |
| 338 | = Detail::ellipseIntegrationPoints(startZ, firstUnitAxis, secondUnitAxis, a, b, normal, characteristicLength); | ||
| 339 | |||
| 340 | // compute total number of points | ||
| 341 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
|
3 | const auto abPoints = ellipseIPs.size(); |
| 342 | 3 | numPoints_ = abPoints*zSamples; | |
| 343 | |||
| 344 | // first move in the first layer of ellipse integration points | ||
| 345 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
|
3 | points_.clear(); |
| 346 |
1/2✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
|
3 | points_.reserve(numPoints_); |
| 347 | 12 | std::move(ellipseIPs.begin(), ellipseIPs.end(), std::back_inserter(points_)); | |
| 348 | |||
| 349 | // then extrude along the z axis | ||
| 350 | 3 | auto zStepVector = zAxis; zStepVector *= zStep; | |
| 351 |
2/2✓ Branch 0 taken 97 times.
✓ Branch 1 taken 3 times.
|
100 | for (std::size_t zs = 1; zs < zSamples; ++zs) |
| 352 |
5/10✓ Branch 1 taken 97 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 97 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 97 times.
✗ Branch 8 not taken.
✓ Branch 10 taken 97 times.
✗ Branch 11 not taken.
✓ Branch 13 taken 97 times.
✗ Branch 14 not taken.
|
485 | std::transform(points_.end() - abPoints, points_.end(), std::back_inserter(points_), |
| 353 | 269640 | [&](const auto& p){ return p + zStepVector; }); | |
| 354 | |||
| 355 | // computation and error correction for integral element to obtain exact integral of constant functions | ||
| 356 | 3 | integrationElement_ = ellipseIntegrationElement*zStep/zStepFactor; | |
| 357 | 3 | const auto meanLocalError = (height*a*b*M_PI - integrationElement_*numPoints_)/Scalar(numPoints_); | |
| 358 | 3 | integrationElement_ += meanLocalError; | |
| 359 | |||
| 360 |
1/2✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
|
3 | if (verbosity > 0) |
| 361 | { | ||
| 362 |
2/4✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 3 times.
✗ Branch 6 not taken.
|
6 | std::cout << "EllipticCylinderIntegration -- number of integration points: " << numPoints_ <<"\n"; |
| 363 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
|
3 | if (verbosity > 1) |
| 364 | ✗ | std::cout << " -- a: " << a << ", b: " << b << ", h: " << height << "\n" | |
| 365 | ✗ | << " -- volume: " << integrationElement_*numPoints_ << " (expected " << height*a*b*M_PI << ")\n" | |
| 366 | ✗ | << " -- corrected a volume error of: " << meanLocalError/integrationElement_*100 << "%\n"; | |
| 367 | } | ||
| 368 | 3 | } | |
| 369 | |||
| 370 | //! obtain ith integration element | ||
| 371 | ✗ | Scalar integrationElement(std::size_t i) const | |
| 372 | ✗ | { return integrationElement_; } | |
| 373 | |||
| 374 | //! obtain ith integration point | ||
| 375 | const GlobalPosition& integrationPoint(std::size_t i) const | ||
| 376 | 278208 | { return points_[i]; } | |
| 377 | |||
| 378 | //! number of samples points | ||
| 379 | ✗ | std::size_t size() const | |
| 380 | ✗ | { return numPoints_; } | |
| 381 | |||
| 382 | private: | ||
| 383 | const Scalar relCharLength_; | ||
| 384 | std::size_t numPoints_; | ||
| 385 | Scalar integrationElement_; | ||
| 386 | std::vector<GlobalPosition> points_; | ||
| 387 | }; | ||
| 388 | |||
| 389 | /*! | ||
| 390 | * \ingroup EmbeddedCoupling | ||
| 391 | * \brief Helper class to integrate over an elliptic domain | ||
| 392 | * \note This is mostly useful if the integral is known and the | ||
| 393 | * integral mass has to be locally distributed to some non-matching domain | ||
| 394 | * \note The algorithm is based on evenly spaced area elements with rejection sampling. | ||
| 395 | * To improve the quality of the integral, increase the sample size. | ||
| 396 | * \note The area integral of a constant function is exact | ||
| 397 | * See Koch et al. (2020) https://doi.org/10.1016/j.jcp.2020.109369 | ||
| 398 | */ | ||
| 399 | template<class Scalar> | ||
| 400 |
2/6✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 3 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
|
6 | class EllipseIntegration |
| 401 | { | ||
| 402 | using GlobalPosition = Dune::FieldVector<Scalar, 3>; | ||
| 403 | |||
| 404 | public: | ||
| 405 | /*! | ||
| 406 | * \brief Constructor | ||
| 407 | * \param relCharLength characteristic relative integration length (real number between 0 and 1) | ||
| 408 | * \note half the characteristic length means 2*2=4 times more integration points | ||
| 409 | */ | ||
| 410 | 3 | explicit EllipseIntegration(const Scalar relCharLength) | |
| 411 |
2/4✓ Branch 1 taken 3 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 3 times.
✗ Branch 5 not taken.
|
6 | : relCharLength_(relCharLength) |
| 412 | {} | ||
| 413 | |||
| 414 | /*! | ||
| 415 | * \brief set geometry of an ellipse | ||
| 416 | * \param center the center position | ||
| 417 | * \param firstAxis first ellipse axis (length corresponding to axis length) | ||
| 418 | * \param secondAxis second ellipse axis (length corresponding to axis length) | ||
| 419 | * \param verbosity the verbosity level (default: 0 -> no terminal output) | ||
| 420 | */ | ||
| 421 | 3 | void setGeometry(const GlobalPosition& center, | |
| 422 | const GlobalPosition& firstAxis, | ||
| 423 | const GlobalPosition& secondAxis, | ||
| 424 | int verbosity = 0) | ||
| 425 | { | ||
| 426 | 3 | const auto a = firstAxis.two_norm(); | |
| 427 | 3 | const auto b = secondAxis.two_norm(); | |
| 428 | 3 | const auto characteristicLength = relCharLength_*a; | |
| 429 | |||
| 430 | 3 | auto firstUnitAxis = firstAxis; firstUnitAxis /= a; | |
| 431 | 3 | auto secondUnitAxis = secondAxis; secondUnitAxis /= b; | |
| 432 | 3 | const auto normal = crossProduct(firstUnitAxis, secondUnitAxis); | |
| 433 | |||
| 434 | // generate integration points | ||
| 435 | 3 | std::tie(points_, integrationElement_) | |
| 436 |
1/2✗ Branch 3 not taken.
✓ Branch 4 taken 3 times.
|
3 | = Detail::ellipseIntegrationPoints(center, firstUnitAxis, secondUnitAxis, a, b, normal, characteristicLength); |
| 437 | |||
| 438 | // store number of sample points | ||
| 439 |
1/2✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
|
3 | const auto abPoints = points_.size(); |
| 440 | 3 | numPoints_ = abPoints; | |
| 441 | |||
| 442 | // computation and error correction for integral element to obtain exact integral of constant functions | ||
| 443 | 3 | const auto meanLocalError = (a*b*M_PI - integrationElement_*numPoints_)/Scalar(numPoints_); | |
| 444 | 3 | integrationElement_ += meanLocalError; | |
| 445 | |||
| 446 |
1/2✓ Branch 0 taken 3 times.
✗ Branch 1 not taken.
|
3 | if (verbosity > 0) |
| 447 | { | ||
| 448 | 6 | std::cout << "EllipseIntegration -- number of integration points: " << numPoints_ << "\n"; | |
| 449 |
1/2✗ Branch 0 not taken.
✓ Branch 1 taken 3 times.
|
3 | if (verbosity > 1) |
| 450 | ✗ | std::cout << " -- a: " << a << ", b: " << b << "\n" | |
| 451 | ✗ | << " -- area: " << integrationElement_*numPoints_ << " (expected " << a*b*M_PI << ")\n" | |
| 452 | ✗ | << " -- corrected an area error of: " << meanLocalError/integrationElement_*100 << "%\n"; | |
| 453 | } | ||
| 454 | 3 | } | |
| 455 | |||
| 456 | //! obtain ith integration element | ||
| 457 | ✗ | Scalar integrationElement(std::size_t i) const | |
| 458 | ✗ | { return integrationElement_; } | |
| 459 | |||
| 460 | //! obtain ith integration point | ||
| 461 | const GlobalPosition& integrationPoint(std::size_t i) const | ||
| 462 | 8568 | { return points_[i]; } | |
| 463 | |||
| 464 | //! number of integration points | ||
| 465 | ✗ | std::size_t size() const | |
| 466 | ✗ | { return numPoints_; } | |
| 467 | |||
| 468 | private: | ||
| 469 | const Scalar relCharLength_; | ||
| 470 | std::size_t numPoints_; | ||
| 471 | Scalar integrationElement_; | ||
| 472 | std::vector<GlobalPosition> points_; | ||
| 473 | }; | ||
| 474 | |||
| 475 | } // end namespace Dumux::EmbeddedCoupling | ||
| 476 | |||
| 477 | #endif | ||
| 478 |