| 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 Geometry | ||
| 10 | * \author Timo Koch | ||
| 11 | * \brief A quadrature based on refinement | ||
| 12 | */ | ||
| 13 | |||
| 14 | #ifndef DUMUX_GEOMETRY_REFINEMENT_QUADRATURERULE_HH | ||
| 15 | #define DUMUX_GEOMETRY_REFINEMENT_QUADRATURERULE_HH | ||
| 16 | |||
| 17 | #include <dune/geometry/quadraturerules.hh> | ||
| 18 | #include <dune/geometry/virtualrefinement.hh> | ||
| 19 | #include <dumux/common/math.hh> | ||
| 20 | |||
| 21 | namespace Dumux { | ||
| 22 | |||
| 23 | /*! | ||
| 24 | * \ingroup Geometry | ||
| 25 | * \brief A "quadrature" based on virtual refinement | ||
| 26 | */ | ||
| 27 | template<typename ct, int mydim> | ||
| 28 | class RefinementQuadratureRule final : public Dune::QuadratureRule<ct, mydim> | ||
| 29 | { | ||
| 30 | public: | ||
| 31 | //! The space dimension | ||
| 32 | static constexpr int dim = mydim; | ||
| 33 | |||
| 34 | //! The highest quadrature order available | ||
| 35 | static constexpr int highest_order = 1; | ||
| 36 | |||
| 37 |
1/4✓ Branch 0 taken 1 times.
✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
2 | ~RefinementQuadratureRule() final {} |
| 38 | |||
| 39 | 1 | RefinementQuadratureRule(Dune::GeometryType type, int levels) | |
| 40 |
1/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
|
1 | : Dune::QuadratureRule<ct, dim>(type, highest_order) |
| 41 | { | ||
| 42 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | auto& refinement = Dune::buildRefinement<dim, ct>(type, type); |
| 43 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | const auto tag = Dune::refinementLevels(levels); |
| 44 | |||
| 45 | // get the vertices | ||
| 46 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | const auto numVertices = refinement.nVertices(tag); |
| 47 |
2/6✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
|
2 | std::vector<Dune::FieldVector<ct, dim>> points; points.reserve(numVertices); |
| 48 |
2/6✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
2 | auto vIt = refinement.vBegin(tag); |
| 49 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
2 | const auto vItEnd = refinement.vEnd(tag); |
| 50 |
3/4✓ Branch 1 taken 46 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 45 times.
✓ Branch 4 taken 1 times.
|
46 | for (; vIt != vItEnd; ++vIt) |
| 51 |
3/8✓ Branch 1 taken 45 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 45 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 45 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
45 | points.emplace_back(vIt.coords()); |
| 52 | |||
| 53 | // go over all elements and get center and volume | ||
| 54 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | const auto numElements = refinement.nElements(tag); |
| 55 |
1/2✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
|
1 | this->reserve(numElements); |
| 56 |
3/8✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 1 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 1 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
3 | auto eIt = refinement.eBegin(tag); |
| 57 |
2/4✓ Branch 1 taken 1 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 1 times.
✗ Branch 4 not taken.
|
3 | const auto eItEnd = refinement.eEnd(tag); |
| 58 |
3/4✓ Branch 1 taken 65 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 64 times.
✓ Branch 4 taken 1 times.
|
65 | for (; eIt != eItEnd; ++eIt) |
| 59 | { | ||
| 60 | // quadrature point in the centroid of the sub-element and weight is its volume | ||
| 61 |
3/8✓ Branch 1 taken 64 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 64 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 64 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
64 | const auto weight = computeVolume_(type, points, eIt.vertexIndices()); |
| 62 |
3/8✓ Branch 1 taken 64 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 64 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 64 times.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
64 | this->emplace_back(eIt.coords(), weight); |
| 63 | } | ||
| 64 | 1 | } | |
| 65 | |||
| 66 | private: | ||
| 67 | template<class VertexIndices> | ||
| 68 | 64 | ct computeVolume_(Dune::GeometryType t, const std::vector<Dune::FieldVector<ct, dim>>& points, const VertexIndices& indices) const | |
| 69 | { | ||
| 70 |
1/2✓ Branch 0 taken 64 times.
✗ Branch 1 not taken.
|
64 | if (t != Dune::GeometryTypes::simplex(2)) |
| 71 | ✗ | DUNE_THROW(Dune::NotImplemented, "Only implemented for 2d simplices"); | |
| 72 | |||
| 73 | 192 | const auto ab = points[indices[1]] - points[indices[0]]; | |
| 74 | 128 | const auto ac = points[indices[2]] - points[indices[0]]; | |
| 75 | using std::abs; | ||
| 76 | 192 | return abs(crossProduct(ab, ac))*0.5; | |
| 77 | } | ||
| 78 | }; | ||
| 79 | |||
| 80 | } // end namespace Dumux | ||
| 81 | |||
| 82 | #endif | ||
| 83 |