GCC Code Coverage Report

Directory: src/oiseau/
File: src/oiseau/dg/nodal/ref_triangle.cpp
Date: 2025-05-24 01:28:39
Exec Total Coverage
Lines: 143 143 100.0%
Functions: 10 10 100.0%
Branches: 116 214 54.2%
Line Branch Exec Source
1 // Copyright (C) 2025 Tiago V. L. Amorim (@tiagovla)
2 //
3 // This file is part of oiseau (https://github.com/tiagovla/oiseau)
4 //
5 // SPDX-License-Identifier: GPL-3.0-or-later
6
7 #include "oiseau/dg/nodal/ref_triangle.hpp"
8
9 #include <array>
10 #include <cmath>
11 #include <cstddef>
12 #include <numbers>
13 #include <xtensor-blas/xlinalg.hpp>
14 #include <xtensor/core/xmath.hpp>
15 #include <xtensor/core/xoperation.hpp>
16 #include <xtensor/core/xshape.hpp>
17 #include <xtensor/core/xtensor_forward.hpp>
18 #include <xtensor/generators/xbuilder.hpp>
19 #include <xtensor/misc/xmanipulation.hpp>
20 #include <xtensor/views/xslice.hpp>
21 #include <xtensor/views/xview.hpp>
22
23 #include "oiseau/dg/nodal/ref_element.hpp"
24 #include "oiseau/dg/nodal/utils.hpp"
25 #include "oiseau/utils/math.hpp"
26
27 namespace oiseau::dg::nodal {
28
29 16 RefTriangle::RefTriangle(unsigned order) : RefElement(order) {
30 16 this->m_np = ((order + 1) * (order + 2)) / 2;
31 16 this->m_nfp = order + 1;
32
2/4
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 16 times.
✗ Branch 5 not taken.
 16 this->m_r = detail::equilateral_xy_to_rs(detail::generate_triangle_nodes(this->m_order));
33
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 this->m_v = this->vandermonde(this->m_r);
34
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 this->m_gv = this->grad_vandermonde(this->m_r);
35
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 this->m_d = this->grad_operator(this->m_v, this->m_gv);
36 16 }
37
38 241 xt::xarray<double> RefTriangle::basis_function(const xt::xarray<double> &ab, int i, int j) {
39
2/4
✓ Branch 1 taken 241 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 241 times.
✗ Branch 5 not taken.
 241 xt::xarray<double> a = xt::col(ab, 0);
40
2/4
✓ Branch 1 taken 241 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 241 times.
✗ Branch 5 not taken.
 241 xt::xarray<double> b = xt::col(ab, 1);
41
1/2
✓ Branch 1 taken 241 times.
✗ Branch 2 not taken.
 241 auto h1 = oiseau::utils::jacobi_p(i, 0.0, 0.0, a);
42
1/2
✓ Branch 1 taken 241 times.
✗ Branch 2 not taken.
 241 auto h2 = oiseau::utils::jacobi_p(j, 2.0 * i + 1.0, 0.0, b);
43
1/2
✓ Branch 6 taken 241 times.
✗ Branch 7 not taken.
 482 return std::numbers::sqrt2 * h1 * h2 * xt::pow(1 - b, i);
44 241 }
45
46 241 xt::xarray<double> RefTriangle::grad_basis_function(const xt::xarray<double> &ab, int i, int j) {
47
2/4
✓ Branch 1 taken 241 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 241 times.
✗ Branch 5 not taken.
 241 xt::xarray<double> a = xt::col(ab, 0);
48
2/4
✓ Branch 1 taken 241 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 241 times.
✗ Branch 5 not taken.
 241 xt::xarray<double> b = xt::col(ab, 1);
49
50
1/2
✓ Branch 1 taken 241 times.
✗ Branch 2 not taken.
 241 xt::xarray<double> fa = oiseau::utils::jacobi_p(i, 0.0, 0.0, a);
51
1/2
✓ Branch 1 taken 241 times.
✗ Branch 2 not taken.
 241 xt::xarray<double> dfa = oiseau::utils::grad_jacobi_p(i, 0.0, 0.0, a);
52
1/2
✓ Branch 1 taken 241 times.
✗ Branch 2 not taken.
 241 xt::xarray<double> gb = oiseau::utils::jacobi_p(j, 2.0 * i + 1.0, 0.0, b);
53
1/2
✓ Branch 1 taken 241 times.
✗ Branch 2 not taken.
 241 xt::xarray<double> dgb = oiseau::utils::grad_jacobi_p(j, 2.0 * i + 1.0, 0.0, b);
54
55
1/2
✓ Branch 2 taken 241 times.
✗ Branch 3 not taken.
 241 xt::xarray<double> dmodedr = dfa * gb;
56
1/2
✓ Branch 5 taken 241 times.
✗ Branch 6 not taken.
 241 xt::xarray<double> dmodeds = dfa * (gb * (0.5 * (1 + a)));
57
58
1/2
✓ Branch 5 taken 241 times.
✗ Branch 6 not taken.
 241 xt::xarray<double> tmp = dgb * xt::pow(0.5 * (1 - b), i);
59
60
2/2
✓ Branch 0 taken 165 times.
✓ Branch 1 taken 76 times.
 241 if (i > 0) {
61
1/2
✓ Branch 4 taken 165 times.
✗ Branch 5 not taken.
 165 dmodedr *= xt::pow(0.5 * (1 - b), i - 1);
62
1/2
✓ Branch 4 taken 165 times.
✗ Branch 5 not taken.
 165 dmodeds *= xt::pow(0.5 * (1 - b), i - 1);
63
1/2
✓ Branch 6 taken 165 times.
✗ Branch 7 not taken.
 165 tmp += -0.5 * i * gb * xt::pow(0.5 * (1 - b), i - 1);
64 }
65
66
1/2
✓ Branch 3 taken 241 times.
✗ Branch 4 not taken.
 241 dmodeds = dmodeds + fa * tmp;
67
1/2
✓ Branch 2 taken 241 times.
✗ Branch 3 not taken.
 241 dmodedr = pow(2.0, i + 0.5) * dmodedr;
68
1/2
✓ Branch 2 taken 241 times.
✗ Branch 3 not taken.
 241 dmodeds = pow(2.0, i + 0.5) * dmodeds;
69
70
2/4
✓ Branch 2 taken 241 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 241 times.
✗ Branch 6 not taken.
 482 return xt::stack(xt::xtuple(dmodedr, dmodeds), 1);
71 241 }
72
73 16 xt::xarray<double> RefTriangle::vandermonde(const xt::xarray<double> &rs) const {
74
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 const auto ab = detail::rs_to_ab(rs);
75
76 16 const std::size_t n_points = rs.shape()[0];
77 16 const std::size_t n_basis = this->m_np; // (order+1)(order+2)/2
78
79
1/2
✓ Branch 2 taken 16 times.
✗ Branch 3 not taken.
 16 xt::xarray<double> output = xt::zeros<double>({n_points, n_basis});
80
81 16 std::size_t index = 0;
82
2/2
✓ Branch 0 taken 76 times.
✓ Branch 1 taken 16 times.
 92 for (unsigned i = 0; i <= this->m_order; ++i) {
83
2/2
✓ Branch 0 taken 241 times.
✓ Branch 1 taken 76 times.
 317 for (unsigned j = 0; j <= this->m_order - i; ++j, ++index) {
84
3/6
✓ Branch 1 taken 241 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 241 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 241 times.
✗ Branch 8 not taken.
 241 xt::col(output, index) = this->basis_function(ab, i, j);
85 }
86 }
87 16 return output;
88 16 }
89
90 16 xt::xarray<double> RefTriangle::grad_vandermonde(const xt::xarray<double> &rs) const {
91
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 const auto ab = detail::rs_to_ab(rs);
92
93 16 const std::size_t n_points = rs.shape()[0];
94 16 const std::size_t n_basis = this->m_np; // (order+1)(order+2)/2
95 16 const std::size_t dimensions = 2; // dr and ds
96
97
1/2
✓ Branch 2 taken 16 times.
✗ Branch 3 not taken.
 16 xt::xarray<double> output = xt::zeros<double>({n_points, n_basis, dimensions});
98
99 16 std::size_t index = 0;
100
2/2
✓ Branch 0 taken 76 times.
✓ Branch 1 taken 16 times.
 92 for (unsigned i = 0; i <= this->m_order; ++i) {
101
2/2
✓ Branch 0 taken 241 times.
✓ Branch 1 taken 76 times.
 317 for (unsigned j = 0; j <= this->m_order - i; ++j, ++index) {
102
1/2
✓ Branch 1 taken 241 times.
✗ Branch 2 not taken.
 241 auto tmp = this->grad_basis_function(ab, i, j);
103
2/4
✓ Branch 3 taken 241 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 241 times.
✗ Branch 7 not taken.
 241 xt::view(output, xt::all(), index, xt::all()) = tmp;
104 241 }
105 }
106 16 return output;
107 16 }
108
109 16 xt::xarray<double> RefTriangle::grad_operator(const xt::xarray<double> &v,
110 const xt::xarray<double> &gv) const {
111
1/2
✓ Branch 3 taken 16 times.
✗ Branch 4 not taken.
 16 auto gvr = xt::view(gv, xt::all(), xt::all(), 0);
112
1/2
✓ Branch 3 taken 16 times.
✗ Branch 4 not taken.
 16 auto gvs = xt::view(gv, xt::all(), xt::all(), 1);
113
114 16 const auto vt = xt::transpose(v);
115 16 const auto dvt_dr = xt::transpose(gvr);
116 16 const auto dvt_ds = xt::transpose(gvs);
117
118
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 const auto dr = xt::transpose(xt::linalg::solve(vt, dvt_dr));
119
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 const auto ds = xt::transpose(xt::linalg::solve(vt, dvt_ds));
120
121
2/4
✓ Branch 2 taken 16 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 16 times.
✗ Branch 6 not taken.
 32 return xt::stack(xt::xtuple(dr, ds), 2);
122 16 }
123
124 }; // namespace oiseau::dg::nodal
125
126 namespace oiseau::dg::nodal::detail {
127
128 32 xt::xarray<double> rs_to_ab(const xt::xarray<double> &rs) {
129
1/2
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
 32 auto r = xt::col(rs, 0);
130
1/2
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
 32 auto s = xt::col(rs, 1);
131
1/2
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
 32 auto ab = xt::zeros_like(rs);
132
1/2
✓ Branch 2 taken 32 times.
✗ Branch 3 not taken.
 32 xt::xarray<double> denominator = 1.0 - s;
133
1/2
✓ Branch 2 taken 32 times.
✗ Branch 3 not taken.
 32 xt::xarray<bool> is_one = xt::isclose(s, 1.0); // handle s == 1 with tolerance
134
1/2
✓ Branch 6 taken 32 times.
✗ Branch 7 not taken.
 32 xt::xarray<double> a = xt::where(is_one, -1.0, 2.0 * (1.0 + r) / denominator - 1.0);
135
2/4
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 32 times.
✗ Branch 5 not taken.
 32 xt::col(ab, 0) = a;
136
2/4
✓ Branch 1 taken 32 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 32 times.
✗ Branch 5 not taken.
 32 xt::col(ab, 1) = s;
137 32 return ab;
138 32 }
139
140 16 xt::xarray<double> generate_triangle_nodes(unsigned order) {
141 16 constexpr auto alp_opt =
142 std::to_array({0.0000, 0.0000, std::numbers::sqrt2, 0.1001, 0.2751, 0.9800, 1.0999, 1.2832,
143 1.3648, 1.4773, 1.4959, 1.5743, 1.5770, 1.6223, 1.6258});
144
1/2
✓ Branch 0 taken 16 times.
✗ Branch 1 not taken.
 16 double alpha = (order < 16) ? alp_opt[order - 1] : 5.0 / 3.0;
145 16 std::size_t n_p = (order + 1) * (order + 2) / 2;
146
1/2
✓ Branch 2 taken 16 times.
✗ Branch 3 not taken.
 16 xt::xarray<double> l1 = xt::zeros<double>({n_p});
147
1/2
✓ Branch 2 taken 16 times.
✗ Branch 3 not taken.
 16 xt::xarray<double> l3 = xt::zeros<double>({n_p});
148 16 double inv_n = 1.0 / order;
149
150 16 std::size_t idx = 0;
151
2/2
✓ Branch 0 taken 76 times.
✓ Branch 1 taken 16 times.
 92 for (std::size_t i = 0; i < order + 1; ++i) {
152
2/2
✓ Branch 0 taken 241 times.
✓ Branch 1 taken 76 times.
 317 for (std::size_t j = 0; j < order + 1 - i; ++j, ++idx) {
153
1/2
✓ Branch 1 taken 241 times.
✗ Branch 2 not taken.
 241 l1[idx] = i * inv_n;
154
1/2
✓ Branch 1 taken 241 times.
✗ Branch 2 not taken.
 241 l3[idx] = j * inv_n;
155 }
156 }
157
158 16 auto l2 = 1.0 - l1 - l3;
159 16 auto blend1 = 4 * l2 * l3;
160 16 auto blend2 = 4 * l1 * l3;
161 16 auto blend3 = 4 * l1 * l2;
162
163 16 auto e1 = -l2 + l3;
164 16 auto e2 = -l3 + l1;
165 16 auto e3 = -l1 + l2;
166
167
2/4
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 16 times.
✗ Branch 5 not taken.
 16 auto warpf1 = detail::warp_factor(order, e1);
168
2/4
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 16 times.
✗ Branch 5 not taken.
 16 auto warpf2 = detail::warp_factor(order, e2);
169
2/4
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 16 times.
✗ Branch 5 not taken.
 16 auto warpf3 = detail::warp_factor(order, e3);
170
171 16 auto warp1 = blend1 * warpf1 * xt::square(1 + alpha * l1);
172 16 auto warp2 = blend2 * warpf2 * xt::square(1 + alpha * l2);
173 16 auto warp3 = blend3 * warpf3 * xt::square(1 + alpha * l3);
174
175
2/4
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 16 times.
✗ Branch 6 not taken.
 32 xt::xarray<double> out = xt::zeros<double>(xt::dynamic_shape<std::size_t>{n_p, 2});
176
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 auto &&x = xt::col(out, 0);
177
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 auto &&y = xt::col(out, 1);
178
179
1/2
✓ Branch 3 taken 16 times.
✗ Branch 4 not taken.
 16 x = -l2 + l3;
180
1/2
✓ Branch 6 taken 16 times.
✗ Branch 7 not taken.
 16 y = (-l2 - l3 + 2.0 * l1) / std::numbers::sqrt3;
181
182 16 double pi = std::numbers::pi;
183
1/2
✓ Branch 5 taken 16 times.
✗ Branch 6 not taken.
 16 x += warp1 + std::cos(2 * pi / 3) * warp2 + std::cos(4 * pi / 3) * warp3;
184
1/2
✓ Branch 4 taken 16 times.
✗ Branch 5 not taken.
 16 y += std::sin(2 * pi / 3) * warp2 + std::sin(4 * pi / 3) * warp3;
185
186 16 return out;
187 16 }
188
189 16 xt::xarray<double> equilateral_xy_to_rs(const xt::xarray<double> &xy) {
190 16 constexpr double sqrt3 = std::numbers::sqrt3;
191
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 auto x = xt::col(xy, 0);
192
1/2
✓ Branch 1 taken 16 times.
✗ Branch 2 not taken.
 16 auto y = xt::col(xy, 1);
193 16 auto l1 = (std::numbers::sqrt3 * y + 1.0) / 3.0;
194 16 auto l2 = (-3.0 * x - sqrt3 * y + 2.0) / 6.0;
195 16 auto l3 = (3.0 * x - sqrt3 * y + 2.0) / 6.0;
196
197
1/2
✓ Branch 2 taken 16 times.
✗ Branch 3 not taken.
 16 xt::xarray<double> l_coords(xy.shape());
198
2/4
✓ Branch 4 taken 16 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 16 times.
✗ Branch 8 not taken.
 16 xt::col(l_coords, 0) = -l2 + l3 - l1;
199
2/4
✓ Branch 4 taken 16 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 16 times.
✗ Branch 8 not taken.
 16 xt::col(l_coords, 1) = -l2 - (l3 - l1);
200 16 return l_coords;
201 16 }
202
203 48 xt::xarray<double> warp_factor(unsigned order, const xt::xarray<double> &rs) {
204
1/2
✓ Branch 1 taken 48 times.
✗ Branch 2 not taken.
 48 auto lgl_r = oiseau::dg::nodal::utils::jacobi_gl(order, 0.0, 0.0);
205
206
1/2
✓ Branch 2 taken 48 times.
✗ Branch 3 not taken.
 48 xt::xarray<double> r_eq = xt::linspace<double>(-1.0, 1.0, order + 1);
207
2/4
✓ Branch 2 taken 48 times.
✗ Branch 3 not taken.
✓ Branch 5 taken 48 times.
✗ Branch 6 not taken.
 48 xt::xarray<double> v_eq(xt::dynamic_shape<std::size_t>({r_eq.size(), order + 1}));
208
2/2
✓ Branch 0 taken 228 times.
✓ Branch 1 taken 48 times.
 276 for (size_t i = 0; i < order + 1; ++i) {
209
3/6
✓ Branch 1 taken 228 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 228 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 228 times.
✗ Branch 9 not taken.
 228 xt::view(v_eq, xt::all(), i) = oiseau::utils::jacobi_p(i, 0.0, 0.0, r_eq);
210 }
211
212 xt::xarray<double> p_mat =
213
2/4
✓ Branch 2 taken 48 times.
✗ Branch 3 not taken.
✓ Branch 6 taken 48 times.
✗ Branch 7 not taken.
 96 xt::zeros<double>(xt::dynamic_shape<std::size_t>{order + 1, rs.size()});
214
215
2/2
✓ Branch 0 taken 228 times.
✓ Branch 1 taken 48 times.
 276 for (unsigned i = 0; i < order + 1; ++i) {
216
3/6
✓ Branch 1 taken 228 times.
✗ Branch 2 not taken.
✓ Branch 5 taken 228 times.
✗ Branch 6 not taken.
✓ Branch 8 taken 228 times.
✗ Branch 9 not taken.
 228 xt::view(p_mat, i, xt::all()) = oiseau::utils::jacobi_p(i, 0.0, 0.0, rs);
217 }
218
219
1/2
✓ Branch 2 taken 48 times.
✗ Branch 3 not taken.
 48 auto l_mat = xt::linalg::solve(xt::transpose(v_eq), p_mat);
220
1/2
✓ Branch 4 taken 48 times.
✗ Branch 5 not taken.
 48 auto warp = xt::linalg::dot(xt::transpose(l_mat), (-r_eq + lgl_r));
221 48 auto zerof = xt::abs(rs) < 1.0 - 1e-10;
222 48 auto sf = 1.0 - xt::pow((zerof * rs), 2);
223
1/2
✓ Branch 5 taken 48 times.
✗ Branch 6 not taken.
 48 warp = warp / sf + warp * (zerof - 1);
224 48 return warp;
225 48 }
226
227 } // namespace oiseau::dg::nodal::detail
228