GCC Code Coverage Report

Directory: ../../../builds/dumux-repositories/
File: /builds/dumux-repositories/dumux/dumux/multidomain/embedded/localrefinementquadrature.hh
Date: 2024-09-21 20:52:54
Exec Total Coverage
Lines: 59 61 96.7%
Functions: 6 11 54.5%
Branches: 88 98 89.8%
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 A quadrature rule using local refinement to approximate partitioned elements
11 */
12
13 #ifndef DUMUX_MULTIDOMAIN_EMBEDDED_LOCAL_REFINEMENT_QUADRATURE_HH
14 #define DUMUX_MULTIDOMAIN_EMBEDDED_LOCAL_REFINEMENT_QUADRATURE_HH
15
16 #include <array>
17 #include <vector>
18 #include <utility>
19 #include <algorithm>
20
21 #include <dune/common/exceptions.hh>
22
23 namespace Dumux {
24
25 /*!
26 * \ingroup EmbeddedCoupling
27 * \brief A quadrature rule using local refinement to approximate partitioned elements
28 *
29 * An indicator function provides a partition of the given geometry. The function returns the partition id
30 * for a given position. The algorithm then virtually refines the geometry using local refinemt to
31 * approximate the partition boundaries. One integration point per partition is added at the centroid
32 * of the partition.
33 *
34 * This can be used, e.g. for low-order integration of the coupling term
35 * in the projection scheme over the interface facets.
36 * Each facet may be coupled to multiple 1D elements.
37 * See Koch (2021) https://arxiv.org/abs/2106.06358 which also describes this algorithm.
38 */
39 template<class Geometry, class IndicatorFunction>
40
1/4
✓ Branch 0 taken 1512 times.
✗ Branch 1 not taken.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
 3024 class LocalRefinementSimplexQuadrature
41 {
42 using IndicatorTriple = std::array<std::size_t, 3>;
43 using GlobalPosition = typename Geometry::GlobalCoordinate;
44 using LocalPosition = typename Geometry::LocalCoordinate;
45 using Corners = std::array<GlobalPosition, 3>;
46
47 class QuadPoint
48 {
49 LocalPosition localPos_;
50 double weight_;
51 std::size_t indicator_;
52
53 public:
54 2252 QuadPoint(LocalPosition&& localPos, double weight, std::size_t i)
55 2252 : localPos_(weight*std::move(localPos))
56 , weight_(weight)
57 2252 , indicator_(i)
58 {}
59
60 54361 void add(LocalPosition&& localPos, double weight)
61 {
62 108722 localPos_ += weight*localPos;
63 54361 weight_ += weight;
64 54361 }
65
66 void finalize()
67 {
68 2252 localPos_ /= weight_;
69 }
70
71 4504 const LocalPosition& position() const { return localPos_; }
72 ✗ double weight() const { return weight_; }
73 ✗ std::size_t indicator() const { return indicator_; }
74 };
75
76 class Triangle
77 {
78 Corners corners_;
79 public:
80 Triangle(const Corners& c) : corners_(c) {}
81 44204 Triangle(Corners&& c) : corners_(std::move(c)) {}
82 526584 const GlobalPosition& operator[](std::size_t i) const { return corners_[i]; }
83
84 21990 GlobalPosition center() const
85 {
86 21990 GlobalPosition center(0.0);
87
2/2
✓ Branch 0 taken 65970 times.
✓ Branch 1 taken 21990 times.
 131940 for (const auto& c : corners_)
88 65970 center += c;
89 21990 center /= corners_.size();
90 21990 return center;
91 }
92
93 58125 auto volume() const
94 {
95 174375 const auto ab = corners_[1] - corners_[0];
96 174375 const auto ac = corners_[2] - corners_[0];
97 58125 return crossProduct(ab, ac).two_norm()*0.5;
98 }
99 };
100
101 public:
102 1512 LocalRefinementSimplexQuadrature(const Geometry& geo, std::size_t maxLevel, const IndicatorFunction& i)
103 : ind_(i)
104 , maxLevel_(maxLevel)
105
0/2
✗ Branch 1 not taken.
✗ Branch 2 not taken.
 1512 , geometry_(geo)
106 {
107 static constexpr int dim = Geometry::mydimension;
108 static_assert(dim == 2, "Only triangles are supported so far");
109
110 1512 if (geo.corners() != (dim+1))
111 DUNE_THROW(Dune::InvalidStateException, "Only simplex geometries are allowed");
112
113 // evaluate indicator
114 1512 const auto tri = Triangle{{ geo.corner(0), geo.corner(1), geo.corner(2) }};
115 3024 const auto triple = IndicatorTriple{{ ind_(tri[0]), ind_(tri[1]), ind_(tri[2]) }};
116 1512 volume_ = tri.volume();
117
118 // add new sub-qps recursively by adaptive refinement
119
1/2
✓ Branch 1 taken 1512 times.
✗ Branch 2 not taken.
 1512 addSimplex_(tri, 0, triple, triple);
120
121
4/4
✓ Branch 0 taken 2252 times.
✓ Branch 1 taken 1512 times.
✓ Branch 2 taken 2252 times.
✓ Branch 3 taken 1512 times.
 9040 for (auto& qp : qps_)
122 2252 qp.finalize();
123 1512 }
124
125 1512 auto begin() const { return qps_.begin(); }
126 1512 auto end() const { return qps_.end(); }
127 auto size() const { return qps_.size(); }
128
129 private:
130 44204 void addSimplex_(const Triangle& tri, std::size_t level, const IndicatorTriple& triple, const IndicatorTriple& childTriple)
131 {
132 // if all indicator are the same just add the triangle
133
10/10
✓ Branch 0 taken 21990 times.
✓ Branch 1 taken 22214 times.
✓ Branch 2 taken 21990 times.
✓ Branch 3 taken 22214 times.
✓ Branch 4 taken 21990 times.
✓ Branch 5 taken 22214 times.
✓ Branch 6 taken 21990 times.
✓ Branch 7 taken 22214 times.
✓ Branch 8 taken 21990 times.
✓ Branch 9 taken 22214 times.
 221020 if (std::all_of(childTriple.begin(), childTriple.end(), [a0=childTriple[0]](auto a){ return (a == a0); }))
134 {
135 // the volumes need to sum up to 0.5 (triangle reference element volume)
136
26/28
✓ Branch 0 taken 35 times.
✓ Branch 1 taken 280 times.
✓ Branch 2 taken 35 times.
✓ Branch 3 taken 280 times.
✓ Branch 4 taken 124 times.
✓ Branch 5 taken 156 times.
✓ Branch 6 taken 124 times.
✓ Branch 7 taken 156 times.
✓ Branch 8 taken 62 times.
✓ Branch 9 taken 94 times.
✓ Branch 10 taken 62 times.
✓ Branch 11 taken 94 times.
✓ Branch 12 taken 94 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 94 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 461 times.
✓ Branch 17 taken 1337 times.
✓ Branch 18 taken 461 times.
✓ Branch 19 taken 1337 times.
✓ Branch 20 taken 10260 times.
✓ Branch 21 taken 8839 times.
✓ Branch 22 taken 10260 times.
✓ Branch 23 taken 8839 times.
✓ Branch 24 taken 9069 times.
✓ Branch 25 taken 419 times.
✓ Branch 26 taken 9069 times.
✓ Branch 27 taken 419 times.
 127370 auto it = std::find_if(qps_.begin(), qps_.end(), [&](const auto& qp){ return (qp.indicator() == childTriple[0]); });
137
6/6
✓ Branch 0 taken 20105 times.
✓ Branch 1 taken 1885 times.
✓ Branch 2 taken 20105 times.
✓ Branch 3 taken 1885 times.
✓ Branch 4 taken 20105 times.
✓ Branch 5 taken 1885 times.
 65970 if (it != qps_.end())
138 20105 it->add(geometry_.local(tri.center()), 0.5*tri.volume()/volume_);
139 else
140 3770 qps_.emplace_back(geometry_.local(tri.center()), 0.5*tri.volume()/volume_, childTriple[0]);
141 }
142
143 // if they are different but the maximum level is reached, split volume in three
144
2/2
✓ Branch 0 taken 11541 times.
✓ Branch 1 taken 10673 times.
 22214 else if (level == maxLevel_)
145 {
146
2/2
✓ Branch 0 taken 34623 times.
✓ Branch 1 taken 11541 times.
 46164 for (int i = 0; i < 3; ++i)
147 {
148 // the volumes need to sum up to 0.5 (triangle reference element volume)
149
26/28
✓ Branch 0 taken 50 times.
✓ Branch 1 taken 567 times.
✓ Branch 2 taken 50 times.
✓ Branch 3 taken 567 times.
✓ Branch 4 taken 219 times.
✓ Branch 5 taken 348 times.
✓ Branch 6 taken 219 times.
✓ Branch 7 taken 348 times.
✓ Branch 8 taken 135 times.
✓ Branch 9 taken 213 times.
✓ Branch 10 taken 135 times.
✓ Branch 11 taken 213 times.
✓ Branch 12 taken 213 times.
✗ Branch 13 not taken.
✓ Branch 14 taken 213 times.
✗ Branch 15 not taken.
✓ Branch 16 taken 915 times.
✓ Branch 17 taken 2558 times.
✓ Branch 18 taken 915 times.
✓ Branch 19 taken 2558 times.
✓ Branch 20 taken 16632 times.
✓ Branch 21 taken 16088 times.
✓ Branch 22 taken 16632 times.
✓ Branch 23 taken 16088 times.
✓ Branch 24 taken 16092 times.
✓ Branch 25 taken 321 times.
✓ Branch 26 taken 16092 times.
✓ Branch 27 taken 321 times.
 210315 auto it = std::find_if(qps_.begin(), qps_.end(), [&](const auto& qp){ return (qp.indicator() == childTriple[i]); });
150
6/6
✓ Branch 0 taken 34256 times.
✓ Branch 1 taken 367 times.
✓ Branch 2 taken 34256 times.
✓ Branch 3 taken 367 times.
✓ Branch 4 taken 34256 times.
✓ Branch 5 taken 367 times.
 103869 if (it != qps_.end())
151 34256 it->add(geometry_.local(tri[i]), 0.5*tri.volume()/volume_/3.0);
152 else
153 734 qps_.emplace_back(geometry_.local(tri[i]), 0.5*tri.volume()/volume_/3.0, childTriple[i]);
154 }
155 }
156
157 // otherwise refine and go recursively over all children
158 else
159 {
160 10673 const auto children = refine(tri);
161
2/2
✓ Branch 3 taken 42692 times.
✓ Branch 4 taken 10673 times.
 74711 for (const auto& c : children)
162 {
163 // make sure to recompute indicator every time since new indices might come up
164 85384 const auto grandChildTriple = IndicatorTriple{{ ind_(c[0]), ind_(c[1]), ind_(c[2]) }};
165 42692 addSimplex_(c, level+1, triple, grandChildTriple);
166 }
167 }
168 44204 }
169
170 10673 std::array<Triangle, 4> refine(const Triangle& tri) const
171 {
172 42692 const auto p01 = 0.5*(tri[0] + tri[1]);
173 42692 const auto p12 = 0.5*(tri[1] + tri[2]);
174 42692 const auto p20 = 0.5*(tri[2] + tri[0]);
175
176 return {{
177 21346 {{ tri[0], p01, p20 }},
178 21346 {{ tri[1], p01, p12 }},
179 21346 {{ tri[2], p12, p20 }},
180 {{ p12, p01, p20 }}
181 74711 }};
182 }
183
184 const IndicatorFunction &ind_;
185 std::size_t maxLevel_;
186 const Geometry& geometry_;
187 double volume_;
188
189 std::vector<QuadPoint> qps_;
190 };
191
192 } // end namespace Dumux
193
194 #endif
195